// contour.c++
#include "drawwindow.h"  // also calls global.h
#include "screen.h"
#include "contour.h"
#include "math.h"
#include "trace.h"
#include "peptide.h"

void findtriangles(struct griddata& den1,int cm,int ii[3]);
void simplecont(struct griddata& den1,int imin[],int imax[]);
void interdata(); /*generates lookup table used in interpolating points*/
void interpolate(int cm,int pt1,int pt2,int xpt[3]); /*arbitrary interpolate to 1/2 (=50/100)*/
void findcubeval(struct griddata& den1,int i1,int j,int k,int i[9][2]);
void findcubeclass(int i[9][2],int *cm);
void readcubeclass(char classfile[100]);

int idata[9][9][5],iclass[256][5][30],vert[256][5][9];

extern struct geometry geo;
extern Widget toplevel;


void drawsurface(struct griddata den1)  // simpleest version - just draw all tianglles
{
    int i,iv,j;
    float vertx[3],ix[3];

   glInitNames(); /* reset pick buffer */
//   glPushName(-1); // put -1 on the top of the buffer
   glPushName(0); // put -1 on the top of the buffer
   glLoadName(3); /* put 3 on top of buffer, indicates contour */
    glColor3fv(den1.color);  //This the simpler, vector version
//    glLineWidth(geo.linewidth);
    glLineWidth(1.0);  // keep contour thin
    for(i=1;i<=den1.tri[0][0];i++){
      glBegin(GL_LINE_LOOP);
      for(iv=0;iv<=2;iv++){ /*the vertices*/
         for(j=0;j<=2;j++) /*x,y,z; values go from 0 to 100*/
             ix[j]=den1.tri[i][3*iv+j]/100.0;
         mcv(den1,ix,vertx);
         glVertex3fv(vertx);
      }
      glEnd();
   }// triangles
}      


void contourmap(struct griddata& den1,float selectpos[]) 
//Finds contour window about float point selectpos  - modified from setcontour of old maid  - 
// Thus finds contours for map window - does not find pdb or bone window
{  
   int i, imin[3],imax[3],cpos[3];
printf("selectpos = %f %f %f\n",selectpos[0],selectpos[1],selectpos[2]);
printf("Finding contour of map\n");   
   xtoi(den1,selectpos,cpos); /*convert to i,j,k co-ordinates (0<i<imax[0], etc.)*/
   for(i=0;i<=2;i++){ /*set range used by contouring*/
      imin[i]=cpos[i]-den1.crange[i];
      imax[i]=cpos[i]+den1.crange[i];
      if(imin[i]<0) imin[i]=0;
      if(imax[i]>(den1.amax[i]-den1.amin[i])) 
           imax[i]=den1.amax[i]-den1.amin[i];
   }
   simplecont(den1,imin,imax); // simple version - just find trinagles
}

void simplecont(struct griddata& den1,int imin[],int imax[])
// New simple version  - just finds list of triangles
// See routines in contour.cgeneral to see details, etc.
{
  int ii[3],i1,j,k,i[9][2],cm;

  interdata(); /*lookup table used in interpolation*/
  den1.tri[0][0]=0; // set # of triangles = 0
  for(k=imin[2];k<=imax[2];k++) {/*zplane*/
    for(j=imin[1];j<=imax[1];j++)
      for(i1=imin[0];i1<=imax[0];i1++) {
        findcubeval(den1,i1,j,k,i);/*assign magnitide (i[0]) and gsurface # (i[1])*/
        findcubeclass(i,&cm);
        ii[0]=i1;ii[1]=j;ii[2]=k;
        findtriangles(den1,cm,ii);
        if(den1.tri[0][0]> (den1.ntri -6) ){
            errormessage(toplevel,"Too many triangles in contour.c++, Increase NTRI and recompile or decrease radius");
            goto endcont;
        }
     }/*i1*/
   } /*k*/
   endcont:;
printf("Number of triangles in contour = %d\n",den1.tri[0][0]);
}

void findtriangles(struct griddata& den1,int cm,int ii[3])
/* finds list of triangles that make up surface
   tri[0][0] = # of triangle
   tri[i][0,1,2] = vertex of point 1; tri[i][3,4,5], and tri[i][6,7,8] other vertices
*/
{
   int ntri,count,nter,is,jj,jv,jx,pt1,pt2,xpt[10][3],pointnum;
//   int is,ntri,nter,count,pt1,pt2,xpt[10][3],jj,jx,jv,pointnum;

   if(iclass[cm][0][0]>0){ // there is a surface in this cube
     for(is=1;is<=iclass[cm][0][0];is++){ /*loop through surfaces*/ 
        ntri=iclass[cm][is][0];/*#of triangles in surface*/
        nter=iclass[cm][is][1];/* number of intercepts in surface*/
        count=2;
        for(jj=1;jj<=nter;jj++) {/*find position of intercepts*/
              pt1=iclass[cm][is][count];
              pt2=iclass[cm][is][count+1];
              count=count+2;
              interpolate(cm,pt1,pt2,xpt[jj]); /*finds the x,y,z of point jj*/
        }
        for(jj=1;jj<=ntri;jj++){  /*find co-ord. of triangle jj*/
           den1.tri[0][0]=den1.tri[0][0]+1;
           for(jv=0;jv<=2;jv++) {/*find vertex jv*/
               pointnum=iclass[cm][is][count];
               count=count+1;
               for(jx=0;jx<=2;jx++) { /*find x,y,z  mult. origin by 100, divide by 100 when convert ot float */          
                 den1.tri[ den1.tri[0][0] ][3*jv+jx]=ii[jx]*100+xpt[pointnum][jx];// add cube origin xpt
               }
           } /*jv*/
        } //jj
     }/*is*/
   }// iclass
}      


void readcubeclass(char classfile[100])
/* New version, instead of reading file at run time, it is "included" during
    compilation, this means that cubeclass file does not need to be in 
    same file as maps, etc.
  Only thing left to do is define special case for cm=255
*/
{
  int cm,ii;
#include "inccubeclass"  /* this file sets value of arrays class and vert */
  cm=255;/*special case*/
  iclass[cm][0][0]=1; /*1 surface, containing all vertexes*/
  vert[cm][1][0]=8;
   for(ii=1;ii<=8;ii++)
         vert[cm][1][ii]=ii;
}

void findcubeclass(int i[9][2],int *cm)
{
    int ii;
    *cm=0; /*find number of set cm=(0,255) */
    for(ii=0;ii<8;ii++) {
       *cm= *cm+i[8-ii][0];
       if(ii!=7)
       *cm=2*(*cm);
     }
}


void findcubeval(struct griddata& den1,int i1,int j,int k,int i[9][2])
{ 
   /*aassig magnitude values to each vertex*/
         if((int)den1.den[i1][j][k]>den1.conlevel) i[1][0]=1;else i[1][0]=0;
         if((int)den1.den[i1+1][j][k]>den1.conlevel) i[2][0]=1;else i[2][0]=0;
         if((int)den1.den[i1+1][j+1][k]>den1.conlevel) i[3][0]=1;else i[3][0]=0;
         if((int)den1.den[i1][j+1][k]>den1.conlevel) i[4][0]=1;else i[4][0]=0;
         if((int)den1.den[i1][j][k+1]>den1.conlevel) i[5][0]=1;else i[5][0]=0;
         if((int)den1.den[i1+1][j][k+1]>den1.conlevel) i[6][0]=1;else i[6][0]=0;
         if((int)den1.den[i1+1][j+1][k+1]>den1.conlevel) i[7][0]=1;else i[7][0]=0;
         if((int)den1.den[i1][j+1][k+1]>den1.conlevel) i[8][0]=1;else i[8][0]=0;
}

void interpolate(int cm,int pt1,int pt2,int xpt[3]) /*arbitrary interpolate to 1/2 (=50/100)*/
{
    xpt[ idata[pt1][pt2][1] ]=idata[pt1][pt2][3];/*these are the two*/
    xpt[ idata[pt1][pt2][2] ]=idata[pt1][pt2][4];/*  axis that don't inter*/
    xpt[ idata[pt1][pt2][0] ]= 50;/*in general case need interpolation here*/
} 

void interdata() /*generates lookup table used in interpolating points*/
/* array idata of form: (i,j are vertex numbers)
   idata[i][j][0]= axis (0,1,or 2) that interpolate on
   idata[i][j][1,2]=axis don't itermpolate and
     idata[i][j][1 or 2]= idata[i][j][3 or 4] = 0 or 100 (0 or 1)
*/
{
   int i,j,k;
   idata[1][2][0]=0;idata[1][2][1]=1;idata[1][2][2]=2;idata[1][2][3]=0;idata[1][2][4]=0; 
   idata[2][3][0]=1; idata[2][3][1]=0;idata[2][3][2]=2;idata[2][3][3]=100; idata[2][3][4]=0;
   idata[3][4][0]=0;idata[3][4][1]=1;idata[3][4][2]=2;idata[3][4][3]=100;idata[3][4][4]=0; 
   idata[1][4][0]=1; idata[1][4][1]=0;idata[1][4][2]=2;idata[1][4][3]=0; idata[1][4][4]=0;
   idata[1][5][0]=2;idata[1][5][1]=0;idata[1][5][2]=1;idata[1][5][3]=0;idata[1][5][4]=0; 
   idata[2][6][0]=2;idata[2][6][1]=0;idata[2][6][2]=1;idata[2][6][3]=100;idata[2][6][4]=0; 
   idata[3][7][0]=2;idata[3][7][1]=0;idata[3][7][2]=1;idata[3][7][3]=100;idata[3][7][4]=100; 
   idata[4][8][0]=2;idata[4][8][1]=0;idata[4][8][2]=1;idata[4][8][3]=0;idata[4][8][4]=100; 
   idata[5][6][0]=0;idata[5][6][1]=1;idata[5][6][2]=2; idata[5][6][3]=0;idata[5][6][4]=100;
   idata[6][7][0]=1;idata[6][7][1]=0;idata[6][7][2]=2; idata[6][7][3]=100;idata[6][7][4]=100;
   idata[7][8][0]=0;idata[7][8][1]=1;idata[7][8][2]=2;idata[7][8][3]=100;idata[7][8][4]=100; 
   idata[5][8][0]=1;idata[5][8][1]=0;idata[5][8][2]=2;idata[5][8][3]=0;idata[5][8][4]=100;
   for(i=1;i<=8;i++) /*make e.g. (1,2) = (2,1)*/
      for(j=i+1;j<=8;j++)
         for(k=0;k<=4;k++) 
            idata[j][i][k]=idata[i][j][k];
}      
