00001 /*Program to take the input surface and inflate it by changing each node by a point x mm from it along its normal */
00002 
00003 /* Header FILES */
00004    
00005 #include "SUMA_suma.h"
00006 
00007 SUMA_SurfaceViewer *SUMAg_cSV; /*!< Global pointer to current Surface Viewer structure*/
00008 SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00009                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00010 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00011 SUMA_DO *SUMAg_DOv;     /*!< Global pointer to Displayable Object structure vector*/
00012 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00013 SUMA_CommonFields *SUMAg_CF; /*!< Global pointer to structure containing info common to all viewers */
00014 
00015 /* function declarations */
00016 float maximum(int N, float *inarray);
00017 float minimum(int N, float *inarray);
00018 void cmp_surf_usage ();
00019 
00020 /*
00021 File : SUMA_inflate_compare.c
00022 authors : Shruti Japee and Ziad Saad
00023 Date : Thu August 1 10:24:01 EST 2002
00024 Purpose : To inflate a given surface by x mm and then compare original to inflated.
00025 */
00026 
00027 int main (int argc,char *argv[])
00028 {
00029   static char FuncName[]={"SUMA_inflate_compare"}; 
00030   SUMA_SurfSpecFile Spec;
00031   char *specfilename = NULL;
00032   int kar, id;
00033   SUMA_Boolean brk;
00034   SUMA_Boolean SurfIn = NOPE;
00035 
00036   SUMA_SurfaceObject *Surf1 = NULL, *Surf2=NULL;
00037   char *Surf1_FileName = NULL;
00038   char *Surf2_FileName = NULL;
00039   char *Vol1Parent_FileName=NULL;
00040   char *Vol2Parent_FileName=NULL;
00041   
00042   /* for SureFit surface */
00043   SUMA_SFname *Surf1_SFName = NULL, *Surf2_SFName = NULL;
00044    
00045   /* other variables */
00046   int i;
00047   int num_nodes1;
00048   int num_nodes2;  
00049   float P0[3];
00050   float delta_t; 
00051   float P1[3];
00052   float N0[3];
00053   float maxdistance, mindistance;
00054   float *distance;
00055   float Points[2][3], p1[3], p2[3];
00056   SUMA_COLOR_MAP *MyColMap;
00057   SUMA_SCALE_TO_MAP_OPT *MyOpt;
00058   SUMA_COLOR_SCALED_VECT * MySV;
00059   SUMA_MT_INTERSECT_TRIANGLE *triangle;
00060   SUMA_SURF_NORM SN1;
00061   SUMA_SURF_NORM SN2;
00062   SUMA_SurfaceObject *SO1, *SO2;
00063   struct timeval tt; 
00064   FILE *colorfile;
00065   FILE *distancefile;
00066   char colorfilename[1000];
00067   char distancefilename[1000];
00068   char *tag1 = NULL;
00069   char *tag2 = NULL;
00070   char *state1 = NULL;
00071   char *state2 = NULL;
00072   char *hemi = NULL;
00073   float B_dim[3];
00074   SUMA_ISINBOX isin;
00075   SUMA_PATCH *Patch=NULL;
00076   SUMA_Boolean TryFull = NOPE, FullOnly;
00077   SUMA_MEMBER_FACE_SETS *Memb = NULL;
00078   int *FaceSet_tmp;
00079   int N_FaceSet_tmp;
00080   char *fout = NULL;
00081   int inflation = 10; /* inflation in mm */
00082   
00083   /* allocate space for CommonFields structure */
00084   SUMAg_CF = SUMA_Create_CommonFields ();
00085   if (SUMAg_CF == NULL) {
00086     fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Create_CommonFields\n", FuncName);
00087     exit(1);
00088   }
00089   
00090   if (argc < 5) {
00091     cmp_surf_usage();
00092     exit (1);
00093   }
00094   
00095   /* read in the surfaces */
00096   kar = 1;
00097   brk = NOPE;
00098   SurfIn = NOPE;
00099   FullOnly = YUP;
00100   while (kar < argc) {
00101     /* loop accross command line options */
00102     if ((strcmp(argv[kar], "-h") == 0) || (strcmp(argv[kar], "-help") == 0)) {
00103       cmp_surf_usage ();
00104       exit (1);
00105     }
00106     if (!brk && (strcmp(argv[kar], "-hemi")) == 0) {
00107       kar ++;
00108       if (kar >= argc) {
00109         fprintf (SUMA_STDERR, "need argument after -hemi");
00110         exit (1);
00111          }
00112       hemi = argv[kar];
00113       brk = YUP;
00114     }
00115     if (!brk && (strcmp(argv[kar], "-prefix")) == 0) {
00116       kar ++;
00117       if (kar >= argc) {
00118         fprintf (SUMA_STDERR, "need argument after -prefix");
00119         exit (1);
00120       }
00121       fout = argv[kar];
00122       brk = YUP;
00123     }
00124     if (!brk && (strcmp(argv[kar], "-spec")) == 0) {
00125       kar ++;
00126       if (kar >= argc) {
00127         fprintf (SUMA_STDERR, "need argument after -spec ");
00128         exit (1);
00129       }
00130       specfilename = argv[kar];
00131       brk = YUP;
00132     }
00133     if (!brk && (strcmp(argv[kar], "-box")) == 0) {
00134       kar ++;
00135       if (kar+2 >= argc) {
00136         fprintf (SUMA_STDERR, "need 3 arguments after -box");
00137         exit (1);
00138       }
00139       B_dim[0] = atof(argv[kar]); kar ++;
00140       B_dim[1] = atof(argv[kar]); kar ++;
00141       B_dim[2] = atof(argv[kar]);
00142       
00143       FullOnly = NOPE;
00144       brk = YUP;
00145     }
00146     if (!brk) {
00147       fprintf (SUMA_STDERR,"Error %s: Option %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00148       exit (1);
00149     } 
00150     else {      
00151       brk = NOPE;
00152       kar ++;
00153     }
00154   }/* loop across command line options */
00155   
00156   
00157   if (specfilename == NULL) {
00158     fprintf (SUMA_STDERR,"Error %s: No spec filename specified.\n", FuncName);
00159     exit(1);
00160   }
00161   if (!SUMA_Read_SpecFile (specfilename, &Spec)) {
00162     fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
00163     exit(1);
00164   }     
00165 
00166   /**** loading the first surface *****/
00167   if (SUMA_iswordin(Spec.SurfaceType[0], "FreeSurfer") == 1) {
00168     Surf1_FileName = Spec.SurfaceFile[0];
00169     Surf1 = SUMA_Load_Surface_Object(Surf1_FileName, SUMA_FREE_SURFER, SUMA_ASCII, Vol1Parent_FileName);
00170     tag1 =  "FS";
00171   }
00172   else 
00173     if (SUMA_iswordin(Spec.SurfaceType[0], "SureFit") == 1) {
00174       Surf1_SFName = SUMA_malloc(sizeof(SUMA_SFname));
00175       strcpy(Surf1_SFName->name_coord,Spec.CoordFile[0]);
00176       strcpy(Surf1_SFName->name_topo, Spec.TopoFile[0]);
00177       strcpy(Surf1_SFName->name_param, Spec.SureFitVolParam[0]);
00178       Surf1 = SUMA_Load_Surface_Object(Surf1_SFName, SUMA_SUREFIT, SUMA_ASCII,Vol1Parent_FileName);
00179       tag1 =  "SF";
00180     }
00181   state1 = Spec.State[0];
00182   
00183   /**** loading the second surface *****/
00184   if (SUMA_iswordin(Spec.SurfaceType[0], "FreeSurfer") == 1) {
00185     Surf2_FileName = Spec.SurfaceFile[0];
00186     Surf2 = SUMA_Load_Surface_Object(Surf1_FileName, SUMA_FREE_SURFER, SUMA_ASCII, Vol1Parent_FileName);
00187     tag2 =  "FS";
00188   }
00189   else 
00190     if (SUMA_iswordin(Spec.SurfaceType[0], "SureFit") == 1) {
00191       Surf2_SFName = SUMA_malloc(sizeof(SUMA_SFname));
00192       strcpy(Surf2_SFName->name_coord,Spec.CoordFile[0]);
00193       strcpy(Surf2_SFName->name_topo, Spec.TopoFile[0]);
00194       strcpy(Surf2_SFName->name_param, Spec.SureFitVolParam[0]);
00195       Surf2 = SUMA_Load_Surface_Object(Surf1_SFName, SUMA_SUREFIT, SUMA_ASCII,Vol1Parent_FileName);
00196       tag2 =  "SF";
00197     }
00198   state2 = Spec.State[0];
00199   
00200   if (fout == NULL) { /* form default name */
00201     sprintf(distancefilename, "dist_%s_%s_%s_%dmm.txt",hemi,tag1,state1,inflation);
00202     sprintf(colorfilename, "dist_%s_%s_%s_%dmm.col",hemi,tag1,state1,inflation);
00203   } 
00204   else {
00205     sprintf (colorfilename, "%s.col", fout);
00206     sprintf (distancefilename, "%s.txt", fout);
00207   }
00208   
00209   if (SUMA_filexists(colorfilename) || SUMA_filexists(distancefilename)) {
00210     fprintf (SUMA_STDERR,"Error %s: One or both of output files %s, %s exists.\nWill not overwrite.\n", \
00211              FuncName, distancefilename, colorfilename);
00212     exit(1);
00213   }
00214   
00215   /***********************************************************************************************/
00216 
00217   SO1 = Surf1;
00218   SO2 = Surf2;
00219  
00220   //SUMA_Print_Surface_Object ( SO1, NULL);
00221   //SUMA_Print_Surface_Object ( SO2, NULL);
00222 
00223   num_nodes1 = SO1->N_Node;
00224   num_nodes2 = SO2->N_Node;
00225   
00226   fprintf(SUMA_STDERR, "Number of nodes in surface 1: %d \n", num_nodes1);
00227   fprintf(SUMA_STDERR, "Number of nodes in surface 2: %d \n", num_nodes2);
00228   SN1 = SUMA_SurfNorm(SO1->NodeList,  SO1->N_Node, SO1->FaceSetList, SO1->N_FaceSet);
00229   
00230   /* Take each node in the SO1-> Nodelist and its corresponding SN1->NodeNormList.  This is the normalized normal vector to the node on the first surface.
00231      So each node is P0.  P1 is computed as some point along the normal vector to that node.  Lets say P1 is P0 + 10 mm along the normal. 
00232   
00233    Write out P1 as a surface */
00234    /* for each node on the first surface do the following */
00235   SUMA_etime (&tt, 0);
00236   for (i = 0; i < num_nodes1; i++) {
00237     id = SO1->NodeDim * i;
00238          P0[0] = SO1->NodeList[id];
00239     P0[1] = SO1->NodeList[id+1];
00240     P0[2] = SO1->NodeList[id+2];
00241     N0[0] = SN1.NodeNormList[id];
00242     N0[1] = SN1.NodeNormList[id+1];
00243     N0[2] = SN1.NodeNormList[id+2];
00244     SUMA_POINT_AT_DISTANCE(N0, P0, 10, Points);
00245     P1[0] = Points[0][0];
00246     P1[1] = Points[0][1];
00247     P1[2] = Points[0][2];
00248     //    fprintf(SUMA_STDERR,"P0 values: %f\t%f\t%f, P1 values: %f\t%f\t%f\n", P0[0],P0[1],P0[2],P1[0],P1[1],P1[2]);
00249 
00250    /* Saving the new surfaces */
00251     SO2->NodeList[id] = P1[0];
00252     SO2->NodeList[id+1] = P1[1];
00253     SO2->NodeList[id+2] = P1[2];
00254   }
00255 
00256  /* Now use SO1 and SO2 as the two surfaces and compute the distance */
00257 
00258   SN1 = SUMA_SurfNorm(SO1->NodeList,  SO1->N_Node, SO1->FaceSetList, SO1->N_FaceSet);
00259   SN2 = SUMA_SurfNorm(SO2->NodeList,  SO2->N_Node, SO2->FaceSetList, SO2->N_FaceSet);
00260   if (!FullOnly) {
00261     /* get the Node member structure */
00262     fprintf(SUMA_STDOUT, "%s: Computing MemberFaceSets... \n", FuncName);
00263     Memb = SUMA_MemberFaceSets (SO2->N_Node, SO2->FaceSetList, SO2->N_FaceSet, 3, SO2->idcode_str);
00264     if (Memb->NodeMemberOfFaceSet == NULL) {
00265       fprintf(SUMA_STDERR, "Error %s: Failed in SUMA_MemberFaceSets. \n", FuncName);
00266       exit(1);
00267     }
00268   }
00269   
00270   /* Take each node in the SO1-> Nodelist and its corresponding SN1->NodeNormList.  This is the normalized normal vector to the node on the first surface.
00271      So each node is P0.  P1 is computed as some point along the normal vector to that node.  Lets say P1 is P0 + 10 mm along the normal.  */
00272   
00273   /* for each node on the first surface do the following */
00274   
00275   SUMA_etime (&tt, 0);
00276   distance = SUMA_malloc(num_nodes1*sizeof(float));
00277   for (i = 0; i < num_nodes1; i++) {
00278     //for (i = 0; i < 10; i++) {
00279     id = SO1->NodeDim * i;
00280          P0[0] = SO1->NodeList[id];
00281     P0[1] = SO1->NodeList[id+1];
00282     P0[2] = SO1->NodeList[id+2];
00283     N0[0] = SN1.NodeNormList[id];
00284     N0[1] = SN1.NodeNormList[id+1];
00285     N0[2] = SN1.NodeNormList[id+2];
00286     SUMA_POINT_AT_DISTANCE(N0, P0, 1000, Points);
00287     P1[0] = Points[0][0];
00288     P1[1] = Points[0][1];
00289     P1[2] = Points[0][2];
00290     if (!FullOnly) { /* trying to speed up intersection computations by restricting it to nodes within a box */
00291       TryFull = NOPE;
00292       /* search for nodes on surface 2 within xx mm of P0 */
00293       isin = SUMA_isinbox (SO2->NodeList, SO2->N_Node, P0, B_dim, 0);
00294       if (isin.nIsIn) {
00295         /* find the patch of surface 2 that is formed by those intersection nodes */
00296         Patch = SUMA_getPatch (isin.IsIn, isin.nIsIn, SO2->FaceSetList, SO2->N_FaceSet, Memb, 1);
00297         if (Patch == NULL) {
00298           fprintf(SUMA_STDERR, "Error %s: Null returned from SUMA_getPatch.\n", FuncName);
00299           exit (1);
00300         }
00301         
00302         /* Perform the intersection based on that patch using Shruti's version */
00303         FaceSet_tmp = Patch->FaceSetList;
00304         N_FaceSet_tmp = Patch->N_FaceSet;
00305       }else {
00306         fprintf (SUMA_STDOUT, "%s: No nodes in box about node %d. Trying for full surface intersection.\n", FuncName, i);
00307         TryFull = YUP; /* flag to send it to full intersection */
00308       }
00309     } 
00310     
00311     if (FullOnly || TryFull) {
00312       Patch = NULL;
00313       FaceSet_tmp = SO2->FaceSetList;
00314       N_FaceSet_tmp = SO2->N_FaceSet;
00315     }
00316     
00317     /*now try with the segment from Points[0] to Points[1] returned above. */
00318     p1[0] = Points[0][0]; 
00319     p1[1] = Points[0][1]; 
00320     p1[2] = Points[0][2];
00321     p2[0] = Points[1][0]; 
00322     p2[1] = Points[1][1]; 
00323     p2[2] = Points[1][2];
00324     triangle = SUMA_MT_intersect_triangle(p1,p2, SO2->NodeList, SO2->N_Node, SO2->FaceSetList, SO2->N_FaceSet, NULL);
00325     //SUMA_Show_MT_intersect_triangle(triangle, NULL);
00326     if (triangle->N_hits ==0) {
00327       distance[i] = -1;
00328       // fprintf(SUMA_STDERR, "Could not find hit for node %d in either direction.\n", i);
00329     }
00330     else {
00331       distance[i] = sqrtf(pow(triangle->P[0]-P0[0],2)+pow(triangle->P[1]-P0[1],2)+pow(triangle->P[2]-P0[2],2));
00332       // fprintf(SUMA_STDERR, "distance is : %f\n", distance[i]);
00333     
00334     }
00335          
00336     SUMA_Free_MT_intersect_triangle(triangle); 
00337     if (Patch) SUMA_freePatch(Patch);
00338         
00339     /* calculating the time elapsed and remaining */
00340     if (!(i%100)) {
00341       delta_t = SUMA_etime(&tt, 1);
00342       fprintf (SUMA_STDERR, " [%d]/[%d] %.2f/100%% completed. Dt = %.2f min done of %.2f min total\r" ,  i, num_nodes1, (float)i / num_nodes1 * 100, delta_t/60, delta_t/i * num_nodes1/60);
00343     }      
00344   
00345   }
00346 
00347   /* write out the distance file */
00348   if((distancefile = fopen(distancefilename, "w"))==NULL) {
00349     fprintf(SUMA_STDERR, "Could not open file distance.txt.\n");
00350     exit(1);
00351   }
00352   else {  
00353     for (i=0; i < num_nodes1; ++i) {
00354       fprintf (distancefile,"%d\t%f\n", i, distance[i]);
00355     }
00356     fclose (distancefile);
00357   }
00358   
00359   /* output this distance as a color file */
00360   MyColMap = SUMA_GetStandardMap(SUMA_CMAP_MATLAB_DEF_BYR64);
00361   MyOpt = SUMA_ScaleToMapOptInit();
00362   MySV = SUMA_Create_ColorScaledVect(num_nodes1);
00363   mindistance = minimum(num_nodes1, distance);
00364   maxdistance = maximum(num_nodes1, distance);
00365   SUMA_ScaleToMap(distance,num_nodes1,mindistance, maxdistance, MyColMap,MyOpt,MySV);
00366 
00367   
00368   /* write out the distance color file */
00369   if((colorfile = fopen(colorfilename, "w"))==NULL) {
00370     fprintf(SUMA_STDERR, "Could not open file distance.col.\n");
00371     exit(1);
00372   }
00373   else {
00374     for (i=0; i < num_nodes1; ++i) {
00375       fprintf (colorfile,"%d\t%f\t%f\t%f\n", i, MySV->cV[3*i  ], MySV->cV[3*i+1], MySV->cV[3*i+2]);
00376     }
00377     fclose (colorfile);
00378   }
00379   
00380   if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
00381   
00382   return 1;
00383 }
00384 
00385 
00386 /*************************** FUNCTION DEFINITIONS **********************************/
00387 
00388 void cmp_surf_usage ()
00389 {
00390   printf ("\nUsage:  SUMA_inflate_compare \n\t-spec <Spec file>\n\n");
00391   printf ("\n\t-spec <Spec file>: File containing surface specification. This file is typically \n");
00392   printf ("\t                   generated by @SUMA_Make_Spec_FS (for FreeSurfer surfaces) or \n");
00393   printf ("\t                   @SUMA_Make_Spec_SF (for SureFit surfaces). The Spec file should \n");
00394   printf ("\t                   be located in the directory containing the surfaces.\n");
00395   printf ("\n\t-hemi <left or right>: specify the hemisphere being processed \n");
00396   printf ("\n\t[-prefix <fileprefix>]: Prefix for distance and node color output files.\n");
00397   printf ("\t                 This option is optional. Existing file will not be overwritten.\n");
00398   printf ("\n\t-box <wX wY wZ>: restrict intersection computations for nodes \n");
00399   printf ("\t        contained in a box of w* dimensions. This might speed things\n");
00400   printf ("\t        up sometimes if the box dimension needs not be large.\n");
00401   printf ("\t        This option is pretty much useless.\n"); 
00402   /*
00403     printf ("\n\t[-dev]: This option will give access to options that are not well polished for consumption.\n");
00404     printf ("\n\t        \n");
00405   */
00406   printf ("\n\n\tFor more help: http://afni.nimh.nih.gov/ssc/ziad/SUMA/SUMA_doc.htm\n");
00407   printf ("\n\n\tIf you can't get help here, please get help somewhere.\n");
00408   SUMA_Version(NULL);
00409   
00410   printf ("\n\t Shruti Japee LBC/NIMH/NIH shruti@codon.nih.gov Ziad Saad SSSC/NIMH/NIH ziad@nih.gov \n\n");
00411   exit (0);
00412 }
00413 
00414 
00415 
00416 float minimum (int N, float *inarray)
00417 {
00418   float min = inarray[0];
00419   int i;
00420   for (i = 1; i < N; i++)
00421     {
00422       if (inarray[i] < min)
00423         min = inarray[i];
00424     }
00425   return min;
00426 }
00427 /* determine the max in a vector of floats */
00428 float maximum (int N, float *inarray)
00429 {
00430   float max = inarray[0];
00431   int i;
00432   for (i = 1; i < N; i++)
00433     {
00434       if (inarray[i] > max)
00435         max = inarray[i];
00436     }
00437   return max;
00438 }