/* New routine to carry out thinning
13jul99 - bug - when pt 1882 has 4 connections, two to pt 126
   when remove get segementation fauld
15jun99
  REvised routine removeloop:
   1) The was error in old routine - sometimes cut main connection
   2) Now cuts the longest of the 3 connections
NEW VERSION:
  On 12mar98 - added procedure  ("removeloop") to remove small loops:
    point 1 connected to 2 which is connected to 1 -
      Arbitrarily cut 1 to 2 connection.
   Copy procedures in newmaid, but modify to use Standard Template Library, etc.
  Major Changes:
   1) Use MAXCHAR levels rather then 7
   2) Set value that are going to be revoved equal to -MARCHAR rather then +MAXCHAR so can use MAXCHAR as a real value
   3) Set removed value = REMOVE (- =100), rather then 0
*/

#include "drawwindow.h"
//#include "screen.h"
#include "math.h"
#include "trace.h"
#include "bone.h"
#include "pdbtolev.h"
#include "sphere.h"

#include "string.h"
#include "stdlib.h" 

#define REMOVE  -100
float fminden[10]= {0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6}; // Also used in extendfit.c++
extern struct maiddatafile maidfile;
extern struct spheredata sp[NSPHERES];
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version
extern int verbose; //Defined in trace.c++  if = 1,print some stuff;  = 2 print a lot of stuff

void testbadcon(struct bonedata& bone1);

#if 1  // Moved to bone.c++
void initializemap(struct griddata& den1)
{
    int i;
    for(i=0;i<3;i++)
           den1.griddim[i]=den1.amax[i]-den1.amin[i]+2;/* the +2 is needed for iden array*/
    fprintf(fpout,"Allocating den array dimensions %d %d %d\n",
                           den1.griddim[0],den1.griddim[1],den1.griddim[2]);
    scmat3(den1.den,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
    float alpha,beta,gam;
    alpha=den1.cell[4]*PI/180.;
    beta=den1.cell[5]*PI/180.;
    gam=den1.cell[6]*PI/180.0; 
    den1.cosal=cos(alpha);
    den1.sinal=sin(alpha);
    den1.cosbe=cos(beta);
    den1.sinbe=sin(beta);
    den1.cosgam=cos(gam);
    den1.singam=sin(gam);
    den1.cvy=(den1.cosal-den1.cosbe*den1.cosgam)/den1.singam;
    den1.cvz=sqrt(den1.singam*den1.singam-
            (den1.cosal*den1.cosal+den1.cosbe*den1.cosbe
             -2.0*den1.cosal*den1.cosbe*den1.cosgam))/den1.singam;
    for(i=0;i<=2;i++){
       den1.delx[i]=den1.cell[i+1]/den1.na[i];
       den1.ximin[i]=den1.amin[i]*den1.delx[i]; /* used in mcv the minimum x[i] value*/
       den1.ximax[i]=den1.amax[i]*den1.delx[i]; /* max real space value value*/
    }
fprintf(fpout,"min-max: %5.3f %5.3f;  %5.3f %5.3f;  %5.3f %5.3f\n",
  den1.ximin[0],den1.ximax[0],den1.ximin[1],den1.ximax[1],den1.ximin[2],den1.ximax[2]);
printf("   delx[0,1,2] = %5.3f %5.3f %5.3f\n",den1.delx[0],den1.delx[1],den1.delx[2]);
#if 0  // moved to routine initializemapscreen
    for(i=0;i<=2;i++) // Use default contour radius
          den1.crange[i]=sc.contourradius/den1.delx[i];
    den1.contonoff =1;
    den1.selectcontonoff = -1;
    for(i=0;i<3;i++)
       den1.color[i]=bluev[i];
    den1.conlevel = sc.conlevel*maidfile.denscale;  // Set contourleve to default value
    // Allocate arrays for contour routines
    den1.ntri = NTRI/(den1.delx[0]*den1.delx[1]*den1.delx[2]); // Scale ntri relative del = 1
fprintf(fpout,"Allocating arrays for contour routines  ntriangles = %d\n",den1.ntri);
    imat2(den1.tri,den1.ntri,9);
fprintf(fpout,"   Done\n");
#endif
}
#endif

#if 1  // moved to bone.c++ so can be used for batch mode
void readbinarymap(struct griddata& den1)
{
   FILE *fps;
   long i,j,k;

   fprintf(fpout,"Reading the binary density file %s\n",den1.mapname);
   fps=fopen(den1.mapname,"rb");
   fread(&den1.na[0],sizeof(den1.na[0]),3,fps);
   fread(&den1.amin[0],sizeof(den1.amin[0]),3,fps);
   fread(&den1.amax[0],sizeof(den1.amax[0]),3,fps);
   fread(&den1.cell[1],sizeof(den1.cell[1]),6,fps);
   initializemap(den1);
#if 0
        for(i=0;i<3;i++)
           den1.griddim[i]=den1.amax[i]-den1.amin[i]+2;/* the +2 is needed for iden array*/
        fprintf(fpout,"Allocating den array dimensions %d %d %d\n",
                           den1.griddim[0],den1.griddim[1],den1.griddim[2]);
        scmat3(den1.den,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
#endif
   for(k=0;k<=den1.amax[2]-den1.amin[2];k++){
       fprintf(fpout,"Reading z level %d\n",k); 
       for(j=0;j<=den1.amax[1]-den1.amin[1];j++) 
	    for(i=0;i<=den1.amax[0]-den1.amin[0];i++) {
                fread(&den1.den[i][j][k],sizeof(den1.den[i][j][k]),1,fps);
            }
   }
   fclose(fps);
}
#endif  


void bsort(struct griddata& den1,struct bonedata& bone1,std::vector<int> sortv[],
      int bonemin,int sortnum[],int *maxsize)
// in each sort level, points run from 0,1,2, .... <sortnum
{
   int iden,i,j,k;
   for(i=bone1.cmin[0];i<=bone1.cmax[0];i++)
     for(j=bone1.cmin[1];j<=bone1.cmax[1];j++)
        for(k=bone1.cmin[2];k<=bone1.cmax[2];k++){
           if( (iden = (int)den1.den[i][j][k])>=bonemin){
                sortv[iden].push_back(d3tod1(bone1.cmax,i,j,k));// add the 1 dimensional position value to bin sortv 
           }
        }
   *maxsize= 0;
   for(i=bonemin;i<=MAXCHAR;i++){
      sortnum[i]=sortv[i].size();
      if(sortnum[i]> *maxsize)
          *maxsize = sortnum[i];
   }
#if 0
  for(i=bonemin;i<=MAXCHAR;i++)
    fprintf(fpout,"i = %d:  size = %d\n",i,sortv[i].size());
  int totalcap=0; int totalsize= 0;
  for(i=0;i<=MAXCHAR;i++){
     totalsize = totalsize+sortv[i].size();
     totalcap = totalcap+sortv[i].capacity();
  }
  fprintf(fpout," totalcap = %d   total entries = %d\n", totalcap,totalsize);
#endif
}

void setcubeval(int ilev,int ivec[3],int point[30],
             signed char ***iden)
/*
  ivec[0,1,2] is grid position of center point  
sets values of 3x3 cube; numbering is
    according to fig. 1 in paper by Greer.  value = 0 if less then ilev
     =2 if equal to ilev and 3 if greater then ilev;  
     =1 for points that are to be removed at end of pass
*/
{
    int i,ii,jj,kk,j,k,count;
    count=0; /*point #, 1-27*/
    for(k= -1;k<=1;k++)
       for(j= -1;j<=1;j++)
          for(i= -1;i<=1;i++){
             count=count+1;/*runs from 1 to 24 as in paper by Greer*/
             ii=ivec[0]+i; jj=ivec[1]+j; kk=ivec[2]+k; // Must be in map
             if((int)iden[ii][jj][kk]== REMOVE) /*if removed or too low density*/
                   point[count]=0;
//             else if( (int)iden[ii][jj][kk]<=ilev )// I think == ilev works just as well since have already reomved lower values
             else if( (int)iden[ii][jj][kk]==ilev )// IMPORTANT == ilev better than <=
                   point[count]=2;
             else if( (int)iden[ii][jj][kk]== -MAXCHAR ) 
                   point[count]=1; /*point that is going to be removed*/
             else point[count]=3; /* >ilev*/
          }          
} 


void modify(int point[]) /*make bruit force correct for case of 2,4,6,8 postive,
   also make 5 positive*/
{
   if((point[2]==1)&&(point[4]==1)&&(point[6]==1)&&(point[8]==1) )
      point[5]=1;
   else 
   if((point[2]==1)&&(point[10]==1)&&(point[20]==1)&&(point[12]==1) )
      point[11]=1;
   else 
   if((point[4]==1)&&(point[10]==1)&&(point[22]==1)&&(point[16]==1) )
      point[13]=1;
   else 
   if((point[22]==1)&&(point[20]==1)&&(point[24]==1)&&(point[26]==1) )
      point[23]=1;
   else 
   if((point[24]==1)&&(point[12]==1)&&(point[6]==1)&&(point[18]==1) )
      point[15]=1;
   else 
   if((point[18]==1)&&(point[8]==1)&&(point[16]==1)&&(point[26]==1) )
      point[17]=1;
}

 
void convert(int point[],int x[4][12])
{
   int k;
   for(k=1;k<=27;k++){ /*convert point to 0's and 1's*/
       if(point[k]>1) point[k]=1; /* = 2 or 3*/
          else point[k]=0;
   }   
   modify(point); /*correction to Toriwaki algorithm*/
   x[1][2]=point[1];
   x[1][3]=point[2];
   x[1][4]=point[3];
   x[2][2]=point[4];
   x[2][3]=point[5];
   x[2][4]=point[6];
   x[3][2]=point[7];
   x[3][3]=point[8];
   x[3][4]=point[9];

   x[1][1]=point[10];
   x[1][0]=point[11];
   x[1][5]=point[12];
   x[2][1]=point[13];
   x[2][0]=point[14];
   x[2][5]=point[15];
   x[3][1]=point[16];
   x[3][0]=point[17];
   x[3][5]=point[18];

   x[1][8]=point[19];
   x[1][7]=point[20];
   x[1][6]=point[21];
   x[2][8]=point[22];
   x[2][7]=point[23];
   x[2][6]=point[24];
   x[3][8]=point[25];
   x[3][7]=point[26];
   x[3][6]=point[27];
}


int connected(int x[4][12])
{
   int i,j,h,k,ksum,hsum,term2,term3,xc[4][12],con;
     
   for(i=1;i<=3;i++) /*find complement 1-x*/
     for(j=0;j<=8;j++)
       xc[i][j]=1-x[i][j];
   for(h=1;h<=3;h++)
      xc[h][9]=xc[h][1];
   term2=0; 
   for(h=1;h<=3;h=h+2) {
      ksum=0;
      for(k=1;k<=7;k=k+2)
         ksum=ksum+xc[h][k]*xc[2][k];
      term2=term2+xc[h][0]*(1-ksum);
   }
   term3=0;
   for(k=1;k<=7;k=k+2){
      hsum=0;
      for(h=1;h<=3;h=h+2)
         hsum=hsum+xc[h][0]*xc[h][k]*xc[h][k+1]*xc[h][k+2];
      term3=term3+xc[2][k]*(1-xc[2][k+1]*xc[2][k+2]*(1-hsum));
   }
   con=2-term2-term3;
   return con;
}

int tip(int point[])
{
   int i,test;
   test=0;
   for(i=1;i<=27;i++)
     if(i!=14)
       if(point[i]>1){ /* i.e = to or greater than ilev*/
         test=test+1;
         if(test>1) {
             test=1;
             goto loop;
         }
       }
   test=0;
   loop:;
   return test;
}
/* I don't know which of the two "hole" procedures below is better,
    results are similar for both of them. Choose second one because this alloww
    a plane to thin to a net???
 They give identical results on test with albp
*/
#if 0
int hole(int point[])  /*return 0 if hole formed, 1 if not*/  
{
   int i,j,k,result;                           
   if( (point[11]==0)||(point[13]==0)||(point[15]==0)||(point[17]==0) ){
      if( (point[5]==0)||(point[11]==0)||(point[17]==0)||(point[23]==0) ){
        if( (point[5]==0)||(point[13]==0)||(point[15]==0)||(point[23]==0) )
           result= 1;
      }
   }
   else 
      result= 0;
   return result;
}
#endif
#if 1
int hole(int point[])
 /*return 0 if hole formed, 1 if not;
       hole is define by one 6-connected point must be 0"*/  
{
   int result;                           
   if( (point[11]==0)||(point[13]==0)||(point[15]==0)||(point[17]==0)
      ||(point[5]==0)||(point[23]==0) )
           result= 1;
   else 
      result= 0;
   return result;
}
#endif

void thin(int maxsize,struct bonedata& bone1,signed char ***iden,
     int bonemin,std::vector<int> sortv[128],int sortnum[])
/* Similar to thin in old version ("newmaid")
    New aspects:
     1) Use MAXCHAR levels rather then just the 7 in old routine
     2) Use sortv in place of pict  - acts the same - after sorting
          treat like a regular array - reseting value and deleting etc.\  
        This means that must use e.g.  i = 0,  i<sortnum, rather than .begin()
          and .end()
     3) This way, can copy thin.c exactly with
            sortnum  replacing count[j]
         and since sortv[ilev][j] is 1D pos - must convert to 3d
*/
{
    int ilev,ivec[3],point[30],rcount,tor[4][12];
    int jj,j,nrem,totrem,removeany;
    int  *remove;  // allocate = maxsize below
//    vector<int>:: iterator jj;
    
    fprintf(fpout,"Thining  \n");
    remove =new int[maxsize];
    do{/*repeat entire search through density levels until none removed*/
     removeany=0;
fprintf(fpout,"Search for removable points\n");
     for(ilev = bonemin;ilev<=MAXCHAR;ilev++){
//fprintf(fpout,"Working on level= %d\n",ilev);
       do {/*repeat until no more removed for this ilevel*/
         rcount=0; /*counts the number removed*/
         for(jj=0;jj<sortnum[ilev];jj++){ // loop through all points in level # ilev
            d1tod3(bone1.cmax,sortv[ilev][jj],ivec); // convert from 1d to 3d vec[i,j,k]
            setcubeval(ilev,ivec,point,iden);
            if(hole(point)==1){ /*indicates does not form hole*/
               if(tip(point)==1) {/*indicates not a single point or tip*/
                  convert(point,tor); /*convert to numbering of Toriwaki (fig 3a)
                           NOTE THIS INCLUDES THE NEW MODIFICATION*/    
                  if(connected(tor)==1){ /*indicates still connected after delete*/
                     /*assign  point to class that will be removed*/
                     rcount=rcount+1;
                     iden[ivec[0]][ivec[1]][ivec[2]]= -MAXCHAR; /*den level of removable pt*/
                     remove[rcount]=jj;
                  } /*connected*/
             }/*tip*/
            } //hole
          } /*jj*/
        /* This uses efficient algorithm for removing points:  
           1) Finds block of nrem removable points
           2) Then sifts points above block down by nrem until reach next removable pt.
               and then go back to (1)
        */
          if(rcount>0){ /* remove the "removable" points*/ 
            removeany=1;
            totrem=0; /* counts total number removed*/        
            j=remove[1]; /*# of first removable point*/
            remove[rcount+1]=0; /*necessary to reset this for while below*/
            nrem=0;
          while(1){ /*goto statements used to exit loop*/
            do {
              nrem=nrem+1;/* running count of removed points*/
              totrem=totrem+1;
              d1tod3(bone1.cmax,sortv[ilev][j],ivec); // convert from 1d to 3d vec[i,j,k]
              iden[ivec[0]][ivec[1]][ivec[2]]= REMOVE; // removed points = REMOVE
              j=j+1;
              if(j==sortnum[ilev]) goto endj;
            } while( (j==remove[nrem+1])&& (totrem<rcount) );
            while(j!=remove[nrem+1]){ /*shift down by nrem*/
               sortv[ilev][j-nrem]=sortv[ilev][j];
               j=j+1;
               if(j==sortnum[ilev]) goto endj;
            }
          }
          endj:;  
          sortnum[ilev]=sortnum[ilev]-nrem; 
        } 
      } while(rcount>0); /*repeat for ilevel until no more removed*/
    } /*ilev*/
 } while(removeany==1);
}

void checkneigh(int pos,struct bonedata& bone1,
     signed char ***iden,int *toomany,int *maxpos)
// Check neighbors, record ones that are bone pts
// maxpos is the maxiumum value of pos of the connected point, used below
//      to make temporary 1 D array 
{
   int i,j,k,ivec[3];
   struct bonept pt;

   d1tod3(bone1.cmax,pos,ivec); // convert from 1d to 3d vec[i,j,k]
   pt.pos = pos;
   if(pt.pos > *maxpos)
     *maxpos = pt.pos;

if(pos == 38561){
 d1tod3(bone1.cmax,pos,ivec);
 fprintf(fpout,"\n");
}
   pt.num= 0;
   for(k= -1;k<=1;k++)
     for(j= -1;j<=1;j++)
       for(i= -1;i<=1;i++)
         if( (i!=0)||(j!=0)||(k!=0) ) // not center pt
            if( (int)iden[ivec[0]+i][ivec[1]+j][ivec[2]+k] > REMOVE){
               pt.con[ pt.num ]  = d3tod1(bone1.cmax,ivec[0]+i,ivec[1]+j,ivec[2]+k);// con run from 0,1,...pt.num
               pt.num=pt.num+1;
               if(pt.num == MAXCON){
                  *toomany =  *toomany+1;
                  goto quitcheck;
               }
            }
   quitcheck:;
   bone1.ptv.push_back(pt); // copy struct pt to bone1.pt  
}              
            
     
void findconnect(struct bonedata& bone1,signed char ***iden,
       int bonemin,std::vector<int> sortv[128],int sortnum[],int *maxpos)
// Go through all remaining points - find the connect data
//bone.pt  is defined by vector<struct bonept> pt 
// (.pos = 1D pos, .num = number connected (maxiumu of MAXCON); 
//     .con =label of connected pts
{
    int ilev,i,toomany;

    toomany = 0; // records number of pts that have more then MAXCON connections
    for(ilev = bonemin;ilev<=MAXCHAR;ilev++){
       for(i = 0;i<sortnum[ilev];i++){           
          checkneigh(sortv[ilev][i],bone1,iden,&toomany,maxpos);// maxpos already pointer
       }
    }
    fprintf(fpout," Number of bone pts with more than %d connections = %d\n",
                MAXCON,toomany);
}

void  defineiden(struct bonedata& bone1,signed char ***iden,
      struct griddata& den1,int bonemin)
// Make values that surround boundar = REMOVE
{
   int i,j,k;

   for(i=bone1.cmin[0];i<=bone1.cmax[0];i++)
      for(j=bone1.cmin[1];j<=bone1.cmax[1];j++)
        for(k=bone1.cmin[2];k<=bone1.cmax[2];k++){
          if((int)den1.den[i][j][k]>=bonemin)
             iden[i][j][k]=den1.den[i][j][k];
          else
             iden[i][j][k] = REMOVE;
       }
   i=bone1.cmin[0]-1; 
      for(j=bone1.cmin[1]-1;j<=bone1.cmax[1]+1;j++)
        for(k=bone1.cmin[2]-1;k<=bone1.cmax[2]+1;k++){
             iden[i][j][k] = REMOVE;
       }
   i=bone1.cmax[0]+1; 
      for(j=bone1.cmin[1]-1;j<=bone1.cmax[1]+1;j++)
        for(k=bone1.cmin[2]-1;k<=bone1.cmax[2]+1;k++){
             iden[i][j][k] = REMOVE;
       }
   j=bone1.cmin[1]-1; 
     for(i=bone1.cmin[0]-1;i<=bone1.cmax[0]+1;i++)
        for(k=bone1.cmin[2]-1;k<=bone1.cmax[2]+1;k++){
             iden[i][j][k] = REMOVE;
       }
   j=bone1.cmax[1]+1; 
     for(i=bone1.cmin[0]-1;i<=bone1.cmax[0]+1;i++)
        for(k=bone1.cmin[2]-1;k<=bone1.cmax[2]+1;k++){
             iden[i][j][k] = REMOVE;
       }
   k=bone1.cmin[2]-1; 
   for(i=bone1.cmin[0]-1;i<=bone1.cmax[0]+1;i++)
      for(j=bone1.cmin[1]-1;j<=bone1.cmax[1]+1;j++){
             iden[i][j][k] = REMOVE;
       }
   k=bone1.cmax[2]+1; 
   for(i=bone1.cmin[0]-1;i<=bone1.cmax[0]+1;i++)
      for(j=bone1.cmin[1]-1;j<=bone1.cmax[1]+1;j++){
             iden[i][j][k] = REMOVE;
      }
}

void convtonum(int maxpos,struct bonedata& bone1) // convert connecter point in terms of 1 d pos to number of bone pt
// Converts the bone.ptv.con[i] from the 1D pos to the bone pt number
{
   int i,j;
   int totnum = bone1.ptv.end()-bone1.ptv.begin();
   int *mat = new int[maxpos+5]; // allocate 1 D matrix
   for(i = 0;i<totnum;i++) // mat  convert 1dpos to bone pt num;
      mat[ bone1.ptv[i].pos] = i;
   for(i = 0;i<totnum;i++) // mat  convert 1dpos to bone pt num;
     for(j=0;j<bone1.ptv[i].num;j++) // the number of connections
       bone1.ptv[i].con[j] = mat[ bone1.ptv[i].con[j] ]; //Redefine con to be bone number os positon pt 
   delete []mat;  // free memory   
}

void cutcon(struct bonedata& bone1,int pt1,int pt2)
/* Cut the connection between pts pt1 and pt2
*/
{
     int ii,i2;

     for(ii=0;ii<bone1.ptv[pt1].num;ii++){ // Remove connect from cut1
             if(bone1.ptv[pt1].con[ii]== pt2){
                    for(i2=ii+1;i2<bone1.ptv[pt1].num;i2++) // Redefine the other connections
                           bone1.ptv[pt1].con[i2-1]=bone1.ptv[pt1].con[i2];
                    bone1.ptv[pt1].num=bone1.ptv[pt1].num-1;// Reduce number of connections by 1
                    goto firstcut;
             }
     }//ii
     firstcut:;
     for(ii=0;ii<bone1.ptv[pt2].num;ii++){ // Remove connect from cut2
              if(bone1.ptv[pt2].con[ii]== pt1){
                  for(i2=ii+1;i2<bone1.ptv[pt2].num;i2++) // Redefine the other connections
                           bone1.ptv[pt2].con[i2-1]=bone1.ptv[pt2].con[i2];
                  bone1.ptv[pt2].num=bone1.ptv[pt2].num-1;// Reduce number of connections by 1
                  goto quitrout;
               }
     } //ii
     quitrout:;
}

void addcon(struct bonedata& bone1,int pt1,int pt2)
/* Add connection between pts pt1 and pt2
*/
{
   int icon;

#if 0
  if(pt2== 126)
    printf("Connect in addcon\n");
#endif
   icon = bone1.ptv[pt1].num;
   bone1.ptv[pt1].con[icon] = pt2;
   if(icon< MAXCON-1)
     bone1.ptv[pt1].num= icon+1;
   icon = bone1.ptv[pt2].num;
   bone1.ptv[pt2].con[icon] = pt1;
   if(icon< MAXCON-1)
      bone1.ptv[pt2].num= icon+1;
}

 
void makecut(struct griddata& den1,struct bonedata& bone1,
    int *icut,int pt[],int nloop)
/* Given a loop (either nloop =  3 or 4) find the branch to cut and make cut
*/
{
    int jj,i,j,cut1,cut2;
    float vec[6][3],dist,maxdist;

#if 0
for(i=1;i<=nloop;i++)
  if(pt[i]== 6689)
    printf("CUT  CUT  CUT\n");
#endif
                    *icut = *icut+1;
                    for(i=1;i<=nloop;i++){ // See if any pt has >3 connections - cut branch to it
                       if(bone1.ptv[pt[i]].num>3){
                         cut1=i;
                         jj= i-1;
                         if(jj <1) 
                             jj=nloop;                        
                         if(bone1.ptv[pt[jj]].num<=3){
                               cut2= jj;
                               goto cut;
                         }
                         jj= i+1;
                         if(jj >nloop) 
                             jj=1;                        
                         if(bone1.ptv[pt[jj]].num<=3){
                               cut2= jj;
                               goto cut;
                             }
#if 0
                         for(j=1;j<=nloop;j++){ // do not want second point to have >3 connections - because it is probably part of another back to back loop
                           if(j!=i){
                             if(bone1.ptv[pt[j]].num<=3){
                               cut2= j;
                               goto cut;
                             }
                           }
                         } //j
#endif
                         if(cut1<nloop) // If reach, all nloop  point with more than 3 connections
                             cut2=nloop; // arbitrary assignment
                         else
                             cut2=1;
                         goto cut;
                      } // if i.num>3
                    }
                    for(i=1;i<=nloop;i++) //If reach here - no branch with 4 connectin -  Find longest branch and cut it
                       d1tow(den1,bone1.cmax,bone1.ptv[pt[i]].pos,vec[i]);
                    maxdist = 0.0;
                    for(i=1;i<=nloop;i++){
                       j= 1+ i;
                       if(j>nloop)
                          j=1;
                       dist = distv(vec[i],vec[j]);
                       if(dist>maxdist){
                            maxdist = dist;
                            cut1= i; 
                            cut2=j;
                       }
                    } // i
                    cut:;
                    cutcon(bone1,pt[cut1],pt[cut2]);
#if 0
                    for(ii=0;ii<bone1.ptv[pt[cut1]].num;ii++){ // Remove connect from cut1
                       if(bone1.ptv[pt[cut1]].con[ii]== pt[cut2]){
                         for(i2=ii+1;i2<bone1.ptv[pt[cut1]].num;i2++) // Redefine the other connections
                           bone1.ptv[pt[cut1]].con[i2-1]=bone1.ptv[pt[cut1]].con[i2];
                         bone1.ptv[pt[cut1]].num=bone1.ptv[pt[cut1]].num-1;// Reduce number of connections by 1
                         goto firstcut;
                       }
                    }//ii
                    firstcut:;
                    for(ii=0;ii<bone1.ptv[pt[cut2]].num;ii++){ // Remove connect from cut2
                       if(bone1.ptv[pt[cut2]].con[ii]== pt[cut1]){
                         for(i2=ii+1;i2<bone1.ptv[pt[cut2]].num;i2++) // Redefine the other connections
                           bone1.ptv[pt[cut2]].con[i2-1]=bone1.ptv[pt[cut2]].con[i2];
                         bone1.ptv[pt[cut2]].num=bone1.ptv[pt[cut2]].num-1;// Reduce number of connections by 1
                         goto secondcut;
                       }
                    } //ii
                    secondcut:;
#endif
}


void removeloop(struct griddata& den1,struct bonedata& bone1)
/*
NEW 2jul99, try to remove points with 4 connections
  Look for neighbors with 2 con and 3 con. cut connection to
   the 3 con point and attach it to the 2 con pt.
 If all 4 pts have only 2 con, connect the 2 closest to each and
   cut the branch from one to 4con pt
NEW 15jun99 - Cut out used[] - cause problem if two branches separted by 1 point
   Also - cut whichever bond is is the longest
Procedure to remove short loops:  1 connected to 2 connected to 3  connected to 1
  Arbitrarily cut 1 to 3 connection
Assumle 1 of the 3 pts has connections:
Procedure:
  set  used[totpts] = 0 for all pts - after pt is tested - set to 1
  Search through all points:
    If point i has 3 (used[i] = 0) connections
       follow out all connection - steps  - see if returns
  Set all points that come up in test of connection to pt i = used = 1
*/
{
   int icut,totnum,pt[6];
   int i,j,i4,r4,i2,i3,pt3[6],pt2[6],ipt,min[3];
   int ii,jj,icon[10],ncon,bcon;//7aug01 - modification to removing 4 connnect. pt
   float dist,mindist,vec2[6][3],vec3[6][3];

   totnum = bone1.ptv.end()-bone1.ptv.begin();
printf("REMOVING LOOPS\n");
   icut = 0;
   for(pt[1] = 0;pt[1]<totnum;pt[1]++){ // ipt = initial pt, see if loop back to it
//if(pt[1]==2101)
//   printf("\n");
      repeatloop:; // repeat if > 3 connections and find connection
      if(bone1.ptv[pt[1]].num >=3){ // Assume 1 of the 3 points has  >=3 connections
        for(icon[1]=0;icon[1]< bone1.ptv[pt[1]].num;icon[1]++){ // try all points connected to pt1
           pt[2] = bone1.ptv[pt[1]].con[icon[1]]; // the number of the next connected pt
            for(icon[2]=0;icon[2]< bone1.ptv[pt[2]].num;icon[2]++){ // all points connected to pt2
              pt[3] = bone1.ptv[pt[2]].con[icon[2]]; // the number of the next connected pt
              if(pt[3]!=pt[1]){ // do not go back
                for(icon[3]=0;icon[3]< bone1.ptv[pt[3]].num;icon[3]++){ // all points connected to pt3
                  pt[4] = bone1.ptv[pt[3]].con[icon[3]]; // the number of the next connected pt
                  if(pt[4] == pt[1]){ // then connected to pt1 (loop) - remove this connection
                    makecut(den1,bone1,&icut, pt,3); // cut loop of 3 pts
                    goto foundconnection;  // go to next pt1                      
                  } // pt4 == pt1
#if 1
                  if(pt[4]!=pt[2]){ // do not go back
                    for(icon[4]=0;icon[4]< bone1.ptv[pt[4]].num;icon[4]++){ // all points connected to pt3
                      pt[5] = bone1.ptv[pt[4]].con[icon[4]]; // the number of the next connected pt
                      if(pt[5] == pt[1]){ // then connected to pt1 (loop) - remove this connection
                         makecut(den1,bone1,&icut, pt,4); // cut loop of 4 pts
                         goto foundconnection;  // go to next pt1 
                      } //pt[5] == pt[1]
                    } // icon[4] 
                  } // if p4 != pt2
#endif                   
                } // icon3
              } // if p3 != pt1
             }// icon2
          } // icon1
        } // it pt1.num >=3
      foundconnection:;
    }// pt1 
   fprintf(fpout," REMOVED  %d short loops\n",icut);
#if 1
   i4 = r4 =  0;
   printf("Removing points with 4 connections\n");
   for(pt[1] = 0;pt[1]<totnum;pt[1]++){ //Try to remove point with for connections
if(pt[1]==9423)
 printf("\n");
     if(bone1.ptv[pt[1]].num ==4){ // Assume 1 of the 3 points has  >=3 connections
       i2 = i3 =  0;
       i4 = i4 +1;
       for (i=0;i<4;i++)
          icon[i]=bone1.ptv[pt[1]].con[i];
       for(i=0;i<4;i++){ //7aug01 - if 4 connected pt in 3 loop - just remove that connection
         ipt =  bone1.ptv[pt[1]].con[i];
         ncon=bone1.ptv[ipt].num;
         for(ii=0;ii<ncon;ii++){ //check if pt with 4 connections is connected pt con. to ipt
           bcon=bone1.ptv[ipt].con[ii];
           for (jj=0;jj<4;jj++){
              if(bcon==icon[jj]){
                 cutcon(bone1,pt[1],bcon); // cut connection from pt[1] to double con. pt
                 if(verbose>=1)
                    printf("Cutting connect from 4 connect pt. in a 3 pt loop\n");
                 goto nextpt;
              }
           }
         }
         if(bone1.ptv[ipt].num==3){ // if 3 con, store in pt3
            i3 = i3+1;
            pt3[i3] = ipt;
         }
         else if(bone1.ptv[ipt].num==2){ // if 2 con, store in pt2
            i2 = i2+1;
            pt2[i2] = ipt;
         }
       }// i
printf( " %d pts with 2 con, %d with 3\n",i2,i3);
       mindist = 1000;
       if( (i2!=0)&&(i3!=0) ){ //make new con -  find dist for pt2 set to pt3 set with shortest distance
         r4=r4+1; // counts the removed 4 connections
         for(i=1;i<=i2;i++){
            d1tow(den1,bone1.cmax,bone1.ptv[pt2[i]].pos,vec2[i]);// convert pos to world nvec[3]
            for(j=1;j<=i3;j++){
              d1tow(den1,bone1.cmax,bone1.ptv[pt3[j]].pos,vec3[j]);// convert pos to world nvec[3]
              dist = distv(vec3[j],vec2[i]);
              if(dist<mindist){ // store i and j min in min
                 min[1]=i; 
                 min[2]=j;
                 mindist = dist;
              } 
            }//j
         }// i
         cutcon(bone1,pt[1],pt3[ min[2] ]); // cut connection from pt[1] to 3 con pt
         addcon(bone1,pt3[ min[2] ],pt2[ min[1] ]); // add connectin for 3 con pt to 2 con pt
      } //  (i2!=0)&&(i3!=0) 
      else if(i2>=2){ //  Connect 2 2con pts and cut one connect to 4con
         r4=r4+1; // counts the removed 4 connections
         for(i=1;i<=i2;i++){
            d1tow(den1,bone1.cmax,bone1.ptv[pt2[i]].pos,vec2[i]);// convert pos to world nvec[3]
            for(j=i+1;j<=i2;j++){
              d1tow(den1,bone1.cmax,bone1.ptv[pt2[j]].pos,vec3[j]);// convert pos to world nvec[3]
              dist = distv(vec3[j],vec2[i]);
              if(dist<mindist){ // store i and j min in min
                 min[1]=i; 
                 min[2]=j;
                 mindist = dist;
              } 
            }//j
         }// i
         cutcon(bone1,pt2[ min[2] ],pt[1]); // cut connection from pt[1] to 3 con pt
         addcon(bone1,pt2[ min[2] ],pt2[ min[1] ]); // add connectin for  
      } //i2>=2  
    }// if = 4
    nextpt:;
  }//pt[1]
  printf(" %d points with 4 connections, %d of these were removed\n",i4,r4);
#endif
}                   


void findbone(struct griddata& den1,struct bonedata& bone1)
/*  Corresponds to old thinroutines
 den1 = struct map
 For complete map  cmin = 0;   cmax = den1.amax[i]-den1.amin[i])
 bone.minden = minimum value of density that is possibly considered for bone
 iden is initially identical to den for density>minden;and  = Remove = -100 for points < minden  
       As remove points, set iden = -100;
// maxpos is the maxiumum value of pos of the connected point, used below
//      to make temporary 1 D array 
*/
{
   int i,bonemin,sortnum[128],maxsize,maxpos;//MAXCHAR = 127
   signed char ***iden;
   std::vector<int> sortv[128];// den values range from - 128 to +127 

fprintf(fpout," Finding bones for density file %s using minden %f\n",
     den1.mapname,bone1.minden);
 
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension

//   for(i= 0;i<128;i++)
//     sortv[i].erase(sortv[i].begin(),sortv[i].end());
   bone1.ptv.erase(bone1.ptv.begin(),bone1.ptv.end());// delete all previous points in vector ptv
//   vector<int> sortv[128];// den values range from - 128 to +127 
   scmat3(iden,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
   for(i=0;i<3;i++){ // use complete map - stay one point inside boundary so that neight witin map
      bone1.cmin[i] = 1;
      bone1.cmax[i] = den1.amax[i]-den1.amin[i]-1;
   }


   bonemin = bone1.minden*maidfile.denscale;  // lower limit  used to make bone (round down) 
fprintf(fpout,"bonemin = %d\n",bonemin);
   defineiden(bone1,iden,den1,bonemin); // set point on boundary = REMOVE
   bsort(den1,bone1,sortv,bonemin,sortnum,&maxsize);  // sort density into bins; sort[i] is vector of 1 ..... max
             // of integers = 1D postition
   thin(maxsize,bone1,iden,bonemin,sortv,sortnum);
   maxpos = 0;  // used as indicator of maximum value of 1d positions of bone point
   findconnect(bone1,iden,bonemin,sortv,sortnum,&maxpos);
fprintf(fpout,"maxpos = %d\n",maxpos);
   convtonum(maxpos,bone1); // convert connecter point in terms of 1 d pos to number of bone pt
#if 1
   removeloop(den1,bone1);
   removeloop(den1,bone1); // run it again 
#endif
#if 0 // for debugging
  printf("Running testbadcon\n");
  testbadcon(bone1);
#endif
//fprintf(fpout,"begin end = %d %d\n",bone1.ptv.begin(),bone1.ptv.end());
writebonefile(den1,bone1);
   freescmat3(iden,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// free memory array 3 dim array (mat)
readbonefile(bone1);
}



void readbonefile(struct bonedata& bone1)
{
   int i,j,totnum;
   FILE *fp;
   char comment[100]; 
   struct bonept  bpt;

#if 1
if(verbose>=2)
 fprintf(fpout,"Reading bone file %s\n",bone1.bonename);
   bone1.ptv.erase(bone1.ptv.begin(),bone1.ptv.end());// delete all previous points in vector ptv
   fp = fopen(bone1.bonename,"r");
   fscanf(fp,"%s%f",comment,&bone1.minden);
   fscanf(fp,"%s%f%f%f",comment,
      &bone1.mcolor[0],&bone1.mcolor[1],&bone1.mcolor[2]);
   fscanf(fp,"%s%f%f%f",comment,
     &bone1.scolor[0],&bone1.scolor[1],&bone1.scolor[2]);
   fscanf(fp,"%s%d%d%d%d%d%d",comment,
     &bone1.cmin[0],&bone1.cmin[1],&bone1.cmin[2],&bone1.cmax[0],&bone1.cmax[1],&bone1.cmax[2]); 
   fscanf(fp,"%s%d",comment,&totnum);
   for(i=0;i<totnum;++i){
     fscanf(fp,"%d%d", &bpt.pos,&bpt.num);
     for(j=0;j<bpt.num;j++)
       fscanf(fp,"%d",&bpt.con[j]);
     bone1.ptv.push_back(bpt);
   }
   bone1.iexist =2; // Indicates core read and ready to draw
   fclose(fp);
#endif
}



void writebonefile(struct griddata& den1,struct bonedata& bone1)
{
  int i,j,totnum; 
  FILE *fp;
  char ch[20]; 

   strcpy(bone1.bonename,den1.mapname);
   strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
   sprintf(ch,"%1.1f",bone1.minden);
   strcat(bone1.bonename,ch);
   fp=fopen(bone1.bonename,"w");
   fprintf(fp,"Minden= %f\n",bone1.minden);
   fprintf(fp,"Maincolor= %f %f %f\n",
        bone1.mcolor[0],bone1.mcolor[1],bone1.mcolor[2]);
   fprintf(fp,"Sidecolor= %f %f %f\n",
        bone1.scolor[0],bone1.scolor[1],bone1.scolor[2]);
   fprintf(fp,"cmin,cmax= %d %d %d %d %d %d\n",
     bone1.cmin[0],bone1.cmin[1],bone1.cmin[2],bone1.cmax[0],bone1.cmax[1],bone1.cmax[2]); 
   totnum = bone1.ptv.end()-bone1.ptv.begin();
   fprintf(fp,"totnum= %d\n",totnum);  
   for(i=0;i<totnum;++i){
     fprintf(fp,"%d %d", bone1.ptv[i].pos,bone1.ptv[i].num);
     for(j=0;j<bone1.ptv[i].num;j++)
       fprintf(fp," %d",bone1.ptv[i].con[j]);

     fprintf(fp,"\n");
   }
   fclose(fp);
fprintf(fpout,"Writing bone file %s\n",bone1.bonename);
}

void writebonetopdb(struct griddata& den1,struct geometry& geo,struct bonedata& bone1)
/* Writes the bone points inside spheres in spherefile as C atoms
    in pdb
*/
{
    int i,ipt,totnum,ivec[3],count;
    float fpos[3],sprad2[NSPHERES];
    float sigma,bfact;
    char lresnumch[NCHAR],atomnumch[NCHAR],atomname[NCHAR],resname[NCHAR],seg[NCHAR];
    FILE *fpout;

    
    recallspheresbat(geo);
    for(i=1;i<=maidfile.nsphere;i++) // find square of radius
            sprad2[i]=sqr(sp[i].radius);
    totnum = bone1.ptv.end()-bone1.ptv.begin();
    count=0;
    sigma = 1.0;
    bfact = 20.0;
    strcpy(seg,"");
    fpout = fopen("bone.pdb","w");
    for(ipt =1;ipt<totnum;ipt++){
      d1tod3(bone1.cmax,bone1.ptv[ipt].pos,ivec);
      imcv(den1,ivec,fpos);// convert int grid ivec to float world co-ord fvec
      if(checksphere(fpos,sprad2)==1){ // if in one of the sphere
           count = count+1;
           itoa(count,atomnumch);
           strcpy(atomname,"C");
           strcpy(resname,"BON");
           strcpy(lresnumch,atomnumch);
           fprintf(fpout,"ATOM  %5s  %-3s %3s  %4s     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
                     atomnumch,atomname,resname,lresnumch,
                     fpos[0],fpos[1],fpos[2],sigma,bfact,seg);
       }// in sphere
    }// ipt
    fclose(fpout);
}    
     
void testbadcon(struct bonedata& bone1)
/* Test to see if two connected points are identical
This is only used for debugging to check bone file
*/
{
   int i,j,k,icon[10],totnum,ncon; 
   totnum = bone1.ptv.end()-bone1.ptv.begin();
   for(i = 0;i<totnum;i++){ // ipt = initial pt, see if loop back to it
      ncon = bone1.ptv[i].num;
      if(ncon >=3){ // check to see if 2 connections identical
        for(j=0;j<ncon;j++){ // try all points connected to pt1
          icon[j] = bone1.ptv[i].con[j]; // the pt num of the connect point
        }
        for (j=0;j<ncon;j++)
           for(k=0;k<ncon;k++)
             if( (k!=j)&&(icon[j]==icon[k]) )
               printf("ERROR: Bone pt %d:  icon[%d]=%d icon[%d]=%d\n",
                        i,j,icon[j],k,icon[k]);
     }//ncon>=3
   }//i
}
