00001 #ifdef HAVE_CONFIG_H
00002 #include <config.h>
00003 #endif
00004 
00005 #include "SUMA_suma.h"
00006 
00007 #define SURFPATCH_MAX_SURF 1  /*!< Maximum number of input surfaces */
00008 
00009 #ifdef SUMA_SampBias_STAND_ALONE
00010 #define STAND_ALONE
00011 #endif
00012 
00013 #ifdef STAND_ALONE
00014 /* these global variables must be declared even if they will not be used by this main */
00015 SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/
00016 SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00017                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00018 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00019 SUMA_DO *SUMAg_DOv = NULL;   /*!< Global pointer to Displayable Object structure vector*/
00020 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00021 SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */
00022 #else
00023 extern SUMA_CommonFields *SUMAg_CF;
00024 extern SUMA_DO *SUMAg_DOv;
00025 extern SUMA_SurfaceViewer *SUMAg_SVv;
00026 extern int SUMAg_N_SVv; 
00027 extern int SUMAg_N_DOv;  
00028 #endif
00029 SUMA_Boolean LocalHead = NOPE;
00030 
00031 
00032 void usage_SUMA_SampBias ()
00033 {
00034    static char FuncName[]={"usage_SUMA_SampBias"};
00035    char * s = NULL;
00036    s = SUMA_help_basics();
00037    printf ( "\nUsage:\n"
00038             "  SampBias -spec SPECFILE -surf SURFNAME -plimit limit -dlimit limit -out FILE\n"
00039             "\n"
00040             "  Mandatory parameters:\n"
00041             "     -spec SpecFile: Spec file containing input surfaces.\n"
00042             "     -surf SURFNAME: Name of input surface \n"
00043             "     -plimit limit: maximum length of path along surface in mm.\n"
00044             "                    default is 50 mm\n"
00045             "     -dlimit limit: maximum length of euclidean distance in mm.\n"
00046             "                    default is 1000 mm\n"
00047             "     -out FILE: output dataset\n"
00048             "\n"
00049             "\n"
00050             "%s"
00051             "\n",s); SUMA_free(s); s = NULL;
00052    s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
00053    printf(" blame Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov     \n");
00054    exit (0);
00055 }
00056     
00057 typedef struct {
00058    SUMA_SO_File_Type iType;
00059    char *sv_name;
00060    char *surf_names[SURFPATCH_MAX_SURF];
00061    int N_surf;
00062    char *spec_file;
00063    float plimit;
00064    float dlimit;
00065    char *outfile;
00066    char *histnote;
00067 } SUMA_KUBATEST_OPTIONS;
00068 
00069 
00070 /*!
00071    \brief parse the arguments for SurfSmooth program
00072    
00073    \param argv (char *)
00074    \param argc (int)
00075    \return Opt (SUMA_GETPATCH_OPTIONS *) options structure.
00076                To free it, use 
00077                SUMA_free(Opt->outfile);
00078                SUMA_free(Opt->histnote); 
00079                SUMA_free(Opt);
00080 */
00081 SUMA_KUBATEST_OPTIONS *SUMA_SampBias_ParseInput (char *argv[], int argc, SUMA_KUBATEST_OPTIONS* Opt)
00082 {
00083    static char FuncName[]={"SUMA_SampBias_ParseInput"}; 
00084    int kar, i, ind;
00085    char *outprefix;
00086    SUMA_Boolean brk = NOPE;
00087    SUMA_Boolean LocalHead = NOPE;
00088 
00089    SUMA_ENTRY;
00090 
00091 
00092    kar = 1;
00093    Opt->iType = SUMA_FT_NOT_SPECIFIED;
00094    Opt->sv_name = NULL;
00095    Opt->spec_file = NULL;
00096    Opt->N_surf = -1;
00097    for (i=0; i<SURFPATCH_MAX_SURF; ++i) 
00098       Opt->surf_names[i] = NULL;
00099    brk = NOPE;
00100    Opt->plimit = 50;
00101    Opt->dlimit = 1000;
00102    Opt->outfile = NULL;
00103 
00104    while (kar < argc) 
00105    { /* loop accross command ine options */
00106       /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
00107       if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0)
00108       {
00109          usage_SUMA_SampBias();
00110          exit (0);
00111       }
00112 
00113       SUMA_SKIP_COMMON_OPTIONS(brk, kar);
00114 
00115       if (!brk && (strcmp(argv[kar], "-spec") == 0)) 
00116       {
00117          kar ++;
00118          if (kar >= argc)  
00119          {
00120             fprintf (SUMA_STDERR, "need argument after -spec \n");
00121             exit (1);
00122          }
00123          Opt->spec_file = argv[kar];
00124          brk = YUP;
00125       }
00126       if (!brk && (strcmp(argv[kar], "-surf") == 0))
00127       {
00128          if (kar + 1>= argc)  
00129          {
00130             fprintf (SUMA_STDERR, "need argument after -surf SURF_NAME \n");
00131             exit (1);
00132          }
00133          /*ind = argv[kar][6] - 'A';
00134          if (ind < 0 || ind >= SURFPATCH_MAX_SURF) 
00135          {
00136             fprintf (SUMA_STDERR, "you can use only one surface\n");
00137             exit (1);
00138          }
00139          */
00140          kar ++;
00141          Opt->surf_names[0] = argv[kar];
00142          Opt->N_surf = 1;
00143          brk = YUP;
00144       }
00145 
00146       if (!brk && (strcmp(argv[kar], "-plimit") == 0)) 
00147       {
00148          kar++;
00149          if (kar >= argc)  
00150          {
00151             fprintf (SUMA_STDERR, "need argument after -plimit \n");
00152             exit (1);
00153          }
00154          Opt->plimit = atof(argv[kar]);
00155          brk = YUP;
00156       }
00157      
00158       if (!brk && (strcmp(argv[kar], "-dlimit") == 0)) 
00159       {
00160          kar++;
00161          if (kar >= argc)  
00162          {
00163             fprintf (SUMA_STDERR, "need argument after -dlimit \n");
00164             exit (1);
00165          }
00166          Opt->dlimit = atof(argv[kar]);
00167          brk = YUP;
00168       }
00169      
00170       if (!brk && (strcmp(argv[kar], "-out") == 0)) 
00171       {
00172          kar++;
00173          if (kar >= argc)  
00174          {
00175             fprintf (SUMA_STDERR, "need argument after -out \n");
00176             exit (1);
00177          }
00178          Opt->outfile = argv[kar];
00179          brk = YUP;
00180       }
00181      
00182      if (!brk) 
00183       {
00184          fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00185          exit (1);
00186       } else 
00187       {
00188          brk = NOPE;
00189          kar ++;
00190       }
00191    }
00192 
00193    /* sanity checks */
00194    if (Opt->N_surf != 1) 
00195    {
00196       SUMA_SL_Err("No surface specified.");
00197       exit(1);
00198    }
00199    if (Opt->spec_file == NULL) 
00200    {
00201       SUMA_SL_Err("No spec file specified.");
00202       exit(1);
00203    }
00204    if (Opt->outfile == NULL) 
00205    {
00206       SUMA_SL_Err("No outfile specified.");
00207       exit(1);
00208    }
00209    else if (SUMA_filexists(Opt->outfile))
00210    {
00211       SUMA_SL_Err("that outfile already exists");
00212       exit(1);
00213    }
00214    
00215    Opt->histnote = SUMA_HistString (NULL, argc, argv, NULL);
00216    
00217    SUMA_RETURN (Opt);
00218 }
00219 
00220 void calcWithOffsets(SUMA_SurfaceObject *SO, SUMA_KUBATEST_OPTIONS* Opt)
00221 {
00222    static char FuncName[]={"calcWithOffsets"};
00223    /* initialize OffS */
00224    SUMA_GET_OFFSET_STRUCT *OffS = SUMA_Initialize_getoffsets (SO->N_Node);
00225    struct timeval start_time, start_time_all;
00226    float etime_GetOffset, etime_GetOffset_all;
00227    int i = 0;
00228    FILE* outFile = NULL;
00229    SUMA_Boolean write = YUP;
00230    
00231    SUMA_ENTRY;
00232 
00233    outFile = fopen(Opt->outfile, "w");
00234    fprintf(outFile,  "#Col. 0 Node index\n"
00235                      "#Col. 1 Node for which the ratio in 4 is the largest. (Companion of Node in Col.0)\n"
00236                      "#Col. 2 distance in 3D\n"
00237                      "#Col. 3 shortest surface path\n"
00238                      "#Col. 4 Ratio of path/distance\n");
00239    if (Opt->histnote) {
00240       fprintf(outFile,  "#History:%s\n", Opt->histnote);
00241    }
00242       
00243    SUMA_etime(&start_time_all,0);
00244    for (i=0; i < SO->N_Node; ++i)
00245 /*   i=45552; */
00246    {
00247       float pathD = 0;/*shortest distance along surface to node ii*/
00248       float geomD = 0;/*geometric distance to node ii*/
00249       float ratio = 1;
00250       int j = 0;
00251       int ii = 0;
00252 
00253       /* show me the offset from node 0 */
00254       if (LocalHead) fprintf(SUMA_STDERR,"%s: Calculating offsets from node %d\n",FuncName, i);
00255       if (i == 0) {
00256          SUMA_etime(&start_time,0);
00257       }
00258       SUMA_getoffsets2 (i, SO, Opt->plimit, OffS, NULL, 0);
00259       if (i == 99) {
00260          etime_GetOffset = SUMA_etime(&start_time,1);
00261          fprintf(SUMA_STDERR, "%s: Search to %f mm took %f seconds for %d nodes.\n"
00262                   "Projected completion time: %f minutes\n",
00263                   FuncName, Opt->plimit, etime_GetOffset, i+1,
00264                   etime_GetOffset * SO->N_Node / 60.0 / (i+1));
00265       }
00266       /*find smallest ratio*/
00267       for (j=0; j < OffS->N_Nodes; j++)
00268       {
00269          if( i!=j && OffS->LayerVect[j] >= 0)
00270          {
00271             float x1 = SO->NodeList[i*3+0];
00272             float x2 = SO->NodeList[j*3+0];
00273             float y1 = SO->NodeList[i*3+1];
00274             float y2 = SO->NodeList[j*3+1];
00275             float z1 = SO->NodeList[i*3+2];
00276             float z2 = SO->NodeList[j*3+2];
00277             float dx = x1 - x2;
00278             float dy = y1 - y2;
00279             float dz = z1 - z2;
00280             float d1 = OffS->OffVect[j];
00281             float d2 = sqrt(dx*dx + dy*dy + dz*dz);
00282             float r = d1 / d2;
00283             /*printf("i=%i j=%i offvect=%f dist=%f r=%f layervect=%i\n", i,j,d1,d2,r, OffS->LayerVect[j]); */
00284             if ( d2 < Opt->dlimit && d1 < Opt->plimit && r > ratio )
00285             {
00286                ratio = r;
00287                ii = j;
00288                pathD = d1;
00289                geomD = d2;
00290             }
00291          }
00292       }
00293       if (write)
00294          fprintf(outFile, "%i\t%i\t%f\t%f\t%f\n", i, ii, geomD, pathD, ratio);
00295 
00296       if (LocalHead) 
00297          fprintf(SUMA_STDERR,"%s: Recycling OffS\n", FuncName);
00298       SUMA_Recycle_getoffsets (OffS);
00299       if (LocalHead)
00300          fprintf(SUMA_STDERR,"%s: Done.\n", FuncName);
00301    }
00302    if (write)
00303       fclose(outFile);
00304 
00305    etime_GetOffset_all = SUMA_etime(&start_time_all,1);
00306    fprintf(SUMA_STDERR, "%s: Done.\nSearch to %f mm took %f minutes for %d nodes.\n" ,
00307                             FuncName, Opt->plimit, etime_GetOffset_all / 60.0 , SO->N_Node);
00308    SUMA_Free_getoffsets(OffS);
00309    SUMA_RETURNe;
00310 }
00311 
00312 #ifdef SUMA_SampBias_STAND_ALONE
00313 int main (int argc,char *argv[])
00314 {/* Main */    
00315    static char FuncName[]={"iotest"};
00316    int SO_read = -1;
00317    int i;
00318    SUMA_SurfaceObject *SO = NULL;
00319    SUMA_SurfSpecFile Spec;
00320    void *SO_name = NULL;
00321    SUMA_Boolean LocalHead = NOPE;
00322    SUMA_KUBATEST_OPTIONS Opt;
00323    
00324    SUMA_mainENTRY;
00325    
00326    SUMA_STANDALONE_INIT;
00327    
00328    /* Allocate space for DO structure */
00329    SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
00330    
00331    if (argc < 4)
00332    {
00333       usage_SUMA_SampBias();
00334       exit (1);
00335    }
00336    
00337    SUMA_SampBias_ParseInput (argv, argc, &Opt);
00338    
00339    /* read all surfaces */
00340    if (!SUMA_Read_SpecFile (Opt.spec_file, &Spec))
00341    {
00342       fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
00343       exit(1);
00344    }
00345    SO_read = SUMA_spec_select_surfs(&Spec, Opt.surf_names, SURFPATCH_MAX_SURF, 0);
00346    if ( SO_read != Opt.N_surf )
00347    {
00348        if (SO_read >=0 )
00349           fprintf(SUMA_STDERR,"Error %s:\nFound %d surfaces, expected %d.\n", FuncName,  SO_read, Opt.N_surf);
00350        exit(1);
00351    }
00352    /* now read into SUMAg_DOv */
00353    if (!SUMA_LoadSpec_eng(&Spec, SUMAg_DOv, &SUMAg_N_DOv, Opt.sv_name, 0, SUMAg_CF->DsetList) )
00354    {
00355       fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_LoadSpec_eng\n", FuncName);
00356       exit(1);
00357    }
00358    
00359    /* now identify surface needed */
00360    SO = SUMA_find_named_SOp_inDOv(Opt.surf_names[0], SUMAg_DOv, SUMAg_N_DOv);
00361    if (!SO)
00362    {
00363       fprintf (SUMA_STDERR,"Error %s:\n"
00364                            "Failed to find surface %s\n"
00365                            "in spec file. Use full name.\n",
00366                            FuncName, Opt.surf_names[0]);
00367       exit(1);
00368    }
00369    calcWithOffsets(SO, &Opt);
00370       
00371    SUMA_LH("clean up");
00372    if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) {
00373       SUMA_SL_Err("DO Cleanup Failed!");
00374    }
00375 
00376    if (!SUMA_Free_CommonFields(SUMAg_CF)) {SUMA_SL_Err("SUMAg_CF Cleanup Failed!");}
00377    
00378    if (Opt.outfile) SUMA_free(Opt.outfile);
00379    if (Opt.histnote) SUMA_free(Opt.histnote);
00380    
00381    SUMA_RETURN(0);
00382 }
00383 #endif