00001 #include "SUMA_suma.h"
00002 #include "../thd_brainormalize.h"
00003 
00004 #undef STAND_ALONE
00005 
00006 #if defined SUMA_BrainWrap_STANDALONE
00007 #define STAND_ALONE 
00008 #endif
00009 
00010 #ifdef STAND_ALONE
00011 /* these global variables must be declared even if they will not be used by this main */
00012 SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/
00013 SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00014                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00015 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00016 SUMA_DO *SUMAg_DOv = NULL;   /*!< Global pointer to Displayable Object structure vector*/
00017 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00018 SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */
00019 #else
00020 extern SUMA_CommonFields *SUMAg_CF;
00021 extern SUMA_DO *SUMAg_DOv;
00022 extern SUMA_SurfaceViewer *SUMAg_SVv;
00023 extern int SUMAg_N_SVv; 
00024 extern int SUMAg_N_DOv;  
00025 #endif
00026 
00027 static int InteractiveQuit;
00028 
00029 int SUMA_DidUserQuit(void)
00030 {
00031    return(InteractiveQuit);
00032 }
00033 
00034 /* 
00035    volume returned is negative if failed to read volume. 
00036    Volume units are in mm3
00037    \brief Volume preparations for later stripping:
00038    -1 Run SpatNorm (Not done here at the moment)
00039    -2 Get percentile ranges and sphere center a la BET 
00040       (Fast robust automated brain extraction, Stephen M. Smith, HBM 2002 v 17:3 pp 143-155)
00041    -3 Get convex hull of outer surface and shrink to set limit on expansion
00042 */         
00043 float SUMA_LoadPrepInVol (SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_SurfaceObject **SOhull)
00044 {
00045    static char FuncName[]={"SUMA_LoadPrepInVol"};
00046    int orient, i, nxx, nxy, cnt, *isort = NULL,iPercRangeVal[2], tind, znxy, *ijk=NULL, nf = 0, trouble= 0, Force = 0;
00047    THD_fvec3 ncoord, ndicom;
00048    THD_ivec3 nind3;
00049    float vol=-1, mass=-1, volmass = -1, *CoordList = NULL;
00050    double PercRange[2], PercRangeVal[2], *dvec_sort=NULL, *dvec=NULL;
00051    SUMA_ISINSPHERE IsIn;
00052    SUMA_SurfaceObject *SOhullp = NULL;
00053    SUMA_Boolean LocalHead = NOPE;
00054    
00055    SUMA_ENTRY;
00056    
00057    if (Opt->debug > 2) LocalHead = YUP;
00058    
00059    
00060    Opt->nvox = DSET_NVOX( Opt->in_vol );
00061    if (DSET_NVALS( Opt->in_vol) != 1) {
00062       SUMA_SL_Err("Input volume can only have one sub-brick in it.\nUse [.] selectors to choose sub-brick needed.");
00063       SUMA_RETURN(vol);
00064    }
00065    
00066    Opt->VolCM[0] = THD_BN_XCM; Opt->VolCM[1] = THD_BN_YCM; Opt->VolCM[2] = THD_BN_ZCM;
00067    
00068    Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
00069    if (!Opt->dvec) {
00070       SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
00071       SUMA_RETURN(NOPE);
00072    }
00073    EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
00074                            DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) ,      /* input  */
00075                            MRI_double               , Opt->dvec  ) ;   /* output */
00076    
00077    /* estimate the t2, t98 and tm parameters, should do them regionally by Bob's octant method */
00078    if (Opt->t2 < 0) {
00079       PercRange[0] = 2; PercRange[1] = 98;
00080       dvec_sort = SUMA_dPercRange (Opt->dvec, NULL, Opt->nvox, PercRange, PercRangeVal, iPercRangeVal);
00081       if (!dvec_sort) {
00082          SUMA_SL_Err("Failed to find range");
00083          SUMA_RETURN(NOPE);
00084       }
00085       Opt->t2 = PercRangeVal[0]; Opt->t98 = PercRangeVal[1];
00086       if (!Opt->t98) {
00087          SUMA_SL_Err("98 percentile is 0!\nWhat kind of dataset is this?");
00088          SUMA_RETURN(NOPE);
00089       }
00090       
00091       Opt->t = Opt->t2 + 0.1 * (Opt->t98 - Opt->t2);
00092    
00093       Opt->tm = -1;
00094       if (LocalHead) fprintf (SUMA_STDERR,"%s: Wild shot tm = %f\n", FuncName, dvec_sort[(iPercRangeVal[0]+iPercRangeVal[1])/2]);
00095       if (dvec_sort) SUMA_free(dvec_sort); dvec_sort = NULL;
00096    } else {
00097       /* user specified */
00098    }
00099 
00100    /* estimate the volume */
00101    CoordList = (float *)SUMA_calloc(3 * Opt->nvox, sizeof(float) );
00102    if (!CoordList) {
00103       SUMA_SL_Err("Failed to allocate CoordList");
00104       SUMA_RETURN(vol);
00105    }
00106    Opt->cog[0] = Opt->cog[1] = Opt->cog[2] = 0;
00107    nxx = DSET_NX(Opt->in_vol);
00108    nxy = DSET_NX(Opt->in_vol) * DSET_NY(Opt->in_vol);
00109    vol = 0.0; volmass = 0;
00110    for (i=0; i<Opt->nvox; ++i) { 
00111       /* coord of point */
00112          nind3.ijk[2] = (int)(i/nxy); znxy = nind3.ijk[2] * nxy;
00113          nind3.ijk[1] = ( i - znxy ) / nxx;
00114          nind3.ijk[0] = i - nind3.ijk[1] * nxx - znxy;
00115          ncoord = THD_3dind_to_3dmm(Opt->in_vol, nind3);
00116          ndicom = THD_3dmm_to_dicomm(Opt->in_vol,ncoord);
00117          CoordList[3*i] = ndicom.xyz[0];
00118          CoordList[3*i+1] = ndicom.xyz[1];
00119          CoordList[3*i+2] = ndicom.xyz[2];
00120       if (Opt->dvec[i] > Opt->t) {
00121          mass = SUMA_MIN_PAIR(Opt->t98, Opt->dvec[i]);
00122          Opt->cog[0] += mass * ndicom.xyz[0];
00123          Opt->cog[1] += mass * ndicom.xyz[1];
00124          Opt->cog[2] += mass * ndicom.xyz[2];
00125          volmass += mass;
00126          ++vol;
00127       } 
00128    }
00129    Opt->cog[0] = Opt->cog[0]/volmass; Opt->cog[1] = Opt->cog[1]/volmass; Opt->cog[2] = Opt->cog[2]/volmass; 
00130    if (LocalHead) {
00131       fprintf (SUMA_STDERR,"%s: COG %f %f %f\n",FuncName, Opt->cog[0], Opt->cog[1], Opt->cog[2]); 
00132    }
00133    vol *= fabs(DSET_DX(Opt->in_vol) * DSET_DY(Opt->in_vol) * DSET_DZ(Opt->in_vol) );
00134    
00135    /* find the radius */
00136    Opt->r = pow(vol*3.0/(3.14159*4.0), 1/3.0);
00137    if (LocalHead) {
00138          fprintf (SUMA_STDERR,"%s: Volume %f, radius %f\n", FuncName, vol, Opt->r);
00139    }
00140    
00141    /* form a vector of values inside the sphere */
00142    if (Opt->tm < 0) {
00143       /* create a mask of the voxels inside the sphere of radius r */
00144       IsIn =  SUMA_isinsphere (CoordList, Opt->nvox, Opt->cog , Opt->r , 1 );
00145       dvec = (double *)SUMA_malloc(Opt->nvox * sizeof(double));
00146       if (!dvec) { SUMA_SL_Crit("Failed to allocate for dvec"); SUMA_RETURN(-1); }
00147       cnt = 0;
00148       for (i=0; i<IsIn.nIsIn; ++i) { 
00149          dvec[cnt] = Opt->dvec[IsIn.IsIn[i]]; ++cnt;
00150       }
00151       /* now sort dvec */
00152       isort = SUMA_z_doubqsort(dvec, cnt);
00153       Opt->tm = dvec[cnt/2];
00154       if (LocalHead) {
00155          fprintf (SUMA_STDERR,"%s: Real tm %f\nt2 = %f, t = %f, t98 = %f\n", FuncName, Opt->tm, Opt->t2, Opt->t, Opt->t98);
00156       }
00157       if (dvec) SUMA_free(dvec); dvec = NULL;
00158       if (isort) SUMA_free(isort); isort = NULL;
00159       SUMA_Free_IsInSphere(&IsIn);
00160       
00161    }
00162    
00163    /* kill values over 98 % to help not catch the eyes*/
00164    if (Opt->Kill98) {
00165       int max_mask, cnt;
00166       if (LocalHead) SUMA_SL_Note("Killing values > 98");
00167       max_mask = (int)(Opt->nvox*0.1);
00168       Opt->k98mask = (int *)SUMA_calloc(max_mask, sizeof(int));
00169       if (!Opt->k98mask) {
00170          SUMA_SL_Crit("Failed to allocate");
00171          SUMA_RETURN(-1);
00172       }
00173       cnt = 0;
00174       for (i=0; i<Opt->nvox; ++i) {
00175          if (Opt->dvec[i] >= Opt->t98) {
00176             Opt->dvec[i] = 0;
00177             if (cnt == max_mask) {
00178                SUMA_SL_Warn("Too many high values in dset!");
00179                max_mask *= 2;
00180                Opt->k98mask = (int *)SUMA_realloc(Opt->k98mask, max_mask * sizeof(int));
00181                if (!Opt->k98mask) { SUMA_SL_Crit("Failed to allocate"); SUMA_RETURN(-1); }
00182             } 
00183             Opt->k98mask[cnt] = i; ++cnt;
00184          }  
00185       }
00186       Opt->k98maskcnt = cnt;
00187    }
00188    
00189    /* use CoordList to build a convex hull */
00190    if (*SOhull) {
00191       byte * isskin = NULL;
00192       int N_skin;
00193       if (LocalHead) SUMA_SL_Note("Computing Convex Hull ..");
00194       isskin = SUMA_isSkin(Opt->in_vol, Opt->dvec, Opt->t2, &N_skin);
00195       if (!N_skin || !isskin) {
00196          SUMA_SL_Err("Failed to find skin!\nIgnoring skull limit.");
00197       } else {
00198          cnt = 0;
00199          for (i=0; i<Opt->nvox; ++i) { 
00200             if (isskin[i] > Opt->t2) {
00201                CoordList[3*cnt] = CoordList[3*i] ;
00202                CoordList[3*cnt+1] = CoordList[3*i+1] ;
00203                CoordList[3*cnt+2] = CoordList[3*i+2] ;
00204                ++cnt;
00205             }
00206          }
00207          if (! (nf = SUMA_qhull_wrap(cnt, CoordList, &ijk, 1)) ) {
00208             fprintf(SUMA_STDERR,"Error %s:\nFailed in SUMA_qhull_wrap\n", FuncName);
00209             *SOhull = NULL;
00210          } else {
00211             fprintf(SUMA_STDERR,"%s:\nForming hull.\n", FuncName);
00212             *SOhull = SUMA_Patch2Surf(CoordList, cnt, ijk, nf, 3);
00213             SOhullp = *SOhull;
00214             #if 0
00215             if (LocalHead) fprintf(SUMA_STDERR,"%s: Making hull consistently orientated\n", FuncName);
00216             SOhullp->EL = SUMA_Make_Edge_List_eng (SOhullp->FaceSetList, SOhullp->N_FaceSet, SOhullp->N_Node, SOhullp->NodeList, 1, NULL);
00217             if (!SUMA_MakeConsistent (SOhullp->FaceSetList, SOhullp->N_FaceSet, SOhullp->EL, 0, &trouble)) {
00218                SUMA_SL_Err("Failed in SUMA_MakeConsistent");
00219             }
00220             if (LocalHead) fprintf(SUMA_STDERR,"%s: Checking orientation.\n", FuncName);
00221             SUMA_SL_Warn("Stuff shaky here, Attend to it.");
00222             Force = 1; 
00223             orient = SUMA_OrientTriangles (SOhullp->NodeList, SOhullp->N_Node, SOhullp->FaceSetList, SOhullp->N_FaceSet, 1, Force);
00224             if (orient < 0 || Force) { /* flipping was done, dump the edge list */ if (SOhullp->EL) SUMA_free_Edge_List(SOhullp->EL); SOhullp->EL = NULL; }
00225             if (!orient) { fprintf(SUMA_STDERR,"Error %s:\nFailed in SUMA_OrientTriangles\n", FuncName); }
00226             if (LocalHead) {
00227                if (orient < 0) { SUMA_SL_Note("Hull was reoriented"); }
00228                else { SUMA_SL_Note("Hull was properly oriented"); }
00229             }
00230             #endif 
00231                
00232             if (LocalHead) fprintf(SUMA_STDERR,"%s: Refining hull surface.\n", FuncName);
00233             if (LocalHead) fprintf(SUMA_STDERR,"%s: %d nodes, %d triangles.\n", FuncName, SOhullp->N_Node, SOhullp->N_FaceSet);
00234             /* SUMA_Print_Surface_Object(SOhullp, stderr); */
00235             /* if (!SUMA_Subdivide_Mesh(&(SOhullp->NodeList), &(SOhullp->N_Node), &(SOhullp->FaceSetList), &(SOhullp->N_FaceSet), 50)) { */
00236             if (0) {
00237                SUMA_SL_Err("Failed to subdivide mesh");
00238                SUMA_Free_Surface_Object (SOhullp); 
00239                *SOhull = NULL;
00240             } else {
00241                /* SUMA_Print_Surface_Object(SOhullp, stderr); */
00242                if (SOhullp->NodeNormList) SUMA_free(SOhullp->NodeNormList); SOhullp->NodeNormList = NULL;
00243                if (SOhullp->FaceNormList) SUMA_free(SOhullp->FaceNormList); SOhullp->FaceNormList = NULL;
00244                if (SOhullp->glar_NodeList) SUMA_free(SOhullp->glar_NodeList); SOhullp->glar_NodeList = NULL;
00245                if (SOhullp->glar_FaceSetList) SUMA_free(SOhullp->glar_FaceSetList); SOhullp->glar_FaceSetList = NULL;
00246                if (SOhullp->glar_FaceNormList) SUMA_free(SOhullp->glar_FaceNormList); SOhullp->glar_FaceNormList = NULL;
00247                if (SOhullp->glar_NodeNormList) SUMA_free(SOhullp->glar_NodeNormList); SOhullp->glar_NodeNormList = NULL;
00248                if (LocalHead) fprintf(SUMA_STDERR,"%s: %d nodes, %d triangles.\n", FuncName, SOhullp->N_Node, SOhullp->N_FaceSet);
00249                SUMA_RECOMPUTE_NORMALS((SOhullp));
00250                if (!SUMA_SurfaceMetrics(SOhullp, "EdgeList,MemberFace,PolyArea", NULL)) {
00251                   fprintf(SUMA_STDERR,"Error %s: Error in SUMA_SurfaceMetrics\n", FuncName);
00252                   SUMA_Free_Surface_Object (SOhullp); /* that takes care of freeing leftovers in MF */
00253                   *SOhull = NULL;
00254                }
00255             }
00256 
00257          }
00258       
00259          if (ijk) free(ijk); ijk=NULL;
00260          if (isskin) SUMA_free(isskin); isskin = NULL;
00261       }
00262    }
00263    
00264    if (CoordList) SUMA_free(CoordList); CoordList = NULL;
00265    SUMA_RETURN(vol);
00266 }
00267 
00268 
00269 /*!
00270    \param MinMax (float) 2x1: MinMax[0] minimum under node
00271                               MinMax_dist[1] maximum under node
00272    \param MinMax_dist (float) 2x1: MinMax_dist[0] distance of minimum under node
00273                                    MinMax_dist[1] distance of maximum under node
00274    \param MinMax_over (float) 2x1: MinMax_over[0] minimum over node
00275                                    MinMax_over[1] maximum over node
00276                                    (set this one to NULL if you do not want to search over the node.
00277                                    In that case all the _over (or above like Means[2]) values will
00278                                    be undetermined.
00279    \param MinMax_over_dist (float) 2x1: MinMax_over_dist[0] distance of minimum over node
00280                                         MinMax_over_dist[1] distance of maximum over node
00281    \param Means (float *) 3x1 :  Means[0] value at node, 
00282                                  Means[1] mean value below node (only for vals > Opt->t),
00283                                  Means[2] mean value above node  (only for vals > Opt->t)
00284 */
00285 int SUMA_Find_IminImax (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, 
00286                         int ni, 
00287                         float *MinMax, float *MinMax_dist , 
00288                         float *MinMax_over, float *MinMax_over_dist,
00289                         float *Means, 
00290                         float *undershish, float *overshish, 
00291                         int *dvecind_under, int *dvecind_over, int ShishMax)
00292 {
00293    static char FuncName[]={"SUMA_Find_IminImax"};
00294    float d1, d2, travdir[3], d1t, d2t, lmin, lmax, lmin_dist, lmax_dist;
00295    float t2 = Opt->t2, tm = Opt->tm, t = Opt->t; 
00296    THD_fvec3 ncoord, ndicom;
00297    THD_ivec3 nind3;
00298    int nind, istep, istepmax, istep2max, nxx, nxy, nMeans[3], stopint;
00299    SUMA_ENTRY;
00300    
00301    d1 = Opt->d1; d2 = d1/2.0; 
00302    
00303    Means[0] = Means[1] = Means[2] = 0.0;
00304    nMeans[0] = nMeans[1] = nMeans[2] = 0;
00305    lmin = Opt->tm;
00306    lmin_dist = 0.0;
00307    lmax_dist = 0.0;
00308    lmax = Opt->t;
00309    nxx = DSET_NX(Opt->in_vol);
00310    nxy = DSET_NX(Opt->in_vol) * DSET_NY(Opt->in_vol);
00311    istep = 0; istepmax = (int)(ceil((double)d1/Opt->travstp)); istep2max = (int)(ceil((double)d2/Opt->travstp));
00312    stopint = 0;
00313    while (istep <= istepmax) {
00314       /* calculate offset */
00315       travdir[0] = - istep * Opt->travstp * SO->NodeNormList[3*ni]; travdir[1] = -istep * Opt->travstp * SO->NodeNormList[3*ni+1]; travdir[2] = -istep * Opt->travstp * SO->NodeNormList[3*ni+2]; 
00316       
00317       /* get 1d coord of point */
00318       ndicom.xyz[0] = SO->NodeList[3*ni] + travdir[0] ; ndicom.xyz[1] = SO->NodeList[3*ni+1]+ travdir[1]; ndicom.xyz[2] = SO->NodeList[3*ni+2]+ travdir[2]; 
00319       ncoord = THD_dicomm_to_3dmm(Opt->in_vol, ndicom);
00320       nind3 = THD_3dmm_to_3dind(Opt->in_vol, ncoord);
00321       nind = nind3.ijk[0] + nind3.ijk[1] * nxx + nind3.ijk[2] * nxy;
00322       #if 0
00323       if (ni == Opt->NodeDbg) {
00324          fprintf(SUMA_STDERR, "%s: Node %d\n"
00325                               " nind3 = [%d %d %d] voxVal = %.3f\n", 
00326                               FuncName, Opt->NodeDbg, nind3.ijk[0], nind3.ijk[1], nind3.ijk[2],
00327                               Opt->dvec[nind]);
00328       }
00329       #endif
00330       if (undershish && istep < ShishMax) undershish[istep] = Opt->dvec[nind];   /* store values under node */
00331       
00332       /* track the brain mean, only if the value is above brain non-brain threshold
00333       stop inegrating as soon as you hit nonbrain */
00334       if (Opt->dvec[nind] > Opt->t && !stopint) { Means[1] += Opt->dvec[nind]; ++ nMeans[1]; }
00335       else stopint = 1;
00336       if (dvecind_under) dvecind_under[istep] = nind;
00337         
00338       /* find local min */ 
00339       /* lmin = SUMA_MIN_PAIR(lmin, Opt->dvec[nind]); */
00340       if (lmin > Opt->dvec[nind]) { lmin = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmin_dist); }
00341       
00342       if (istep <= istep2max) {
00343          if (lmax < Opt->dvec[nind]) { lmax = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmax_dist); }  
00344       }
00345       
00346       ++istep;
00347    }
00348    
00349    if (istep < ShishMax) { undershish[istep] = -1; if (dvecind_under) dvecind_under[istep] = -1; }/* crown the shish */
00350    
00351    Means[1] /= (float)nMeans[1];
00352    MinMax[0] = SUMA_MAX_PAIR(t2, lmin);
00353    MinMax_dist[0] = lmin_dist;  
00354    MinMax[1] = SUMA_MIN_PAIR(tm, lmax); 
00355    MinMax_dist[1] = lmax_dist;
00356    
00357    /* search for a minimum overhead */
00358    if (MinMax_over) {
00359       lmin = Opt->tm;
00360       lmin_dist = 0.0;
00361       lmax_dist = 0.0;
00362       lmax = Opt->t;
00363       istep = 0; istepmax = (int)(ceil((double)Opt->d4/Opt->travstp));
00364       stopint = 0;
00365       while (istep <= istepmax) {
00366          /* calculate offset */
00367          travdir[0] =  istep * Opt->travstp * SO->NodeNormList[3*ni]; travdir[1] = istep * Opt->travstp * SO->NodeNormList[3*ni+1]; travdir[2] = istep * Opt->travstp * SO->NodeNormList[3*ni+2]; 
00368 
00369          /* get 1d coord of point */
00370          ndicom.xyz[0] = SO->NodeList[3*ni] + travdir[0] ; ndicom.xyz[1] = SO->NodeList[3*ni+1]+ travdir[1]; ndicom.xyz[2] = SO->NodeList[3*ni+2]+ travdir[2]; 
00371          ncoord = THD_dicomm_to_3dmm(Opt->in_vol, ndicom);
00372          nind3 = THD_3dmm_to_3dind(Opt->in_vol, ncoord);
00373          nind = nind3.ijk[0] + nind3.ijk[1] * nxx + nind3.ijk[2] * nxy;
00374          #if 0
00375          if (ni == Opt->NodeDbg) {
00376             fprintf(SUMA_STDERR, "%s: Node %d\n"
00377                                  " nind3 = [%d %d %d] voxVal = %.3f\n", 
00378                                  FuncName, Opt->NodeDbg, nind3.ijk[0], nind3.ijk[1], nind3.ijk[2],
00379                                  Opt->dvec[nind]);
00380          }
00381          #endif
00382          if (!istep) { /* get the value at the node */
00383             Means[0] = Opt->dvec[nind];
00384          }
00385          
00386          if (overshish && istep < ShishMax) overshish[istep] = Opt->dvec[nind];   /* store values under node */
00387          
00388          if (Opt->dvec[nind] > Opt->t && !stopint) { Means[2] += Opt->dvec[nind]; ++ nMeans[2]; } 
00389          else stopint = 1;
00390          if (dvecind_over) dvecind_over[istep] = nind;
00391          
00392          /* find local min and max*/ 
00393          if (lmin > Opt->dvec[nind]) { lmin = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmin_dist); }
00394          if (lmax < Opt->dvec[nind]) { lmax = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmax_dist); }  
00395          
00396          ++istep;
00397       }
00398       
00399       if (istep < ShishMax) { overshish[istep] = -1; if (dvecind_over) dvecind_over[istep] = -1; }/* crown the shish */
00400 
00401       Means[2] /= (float)nMeans[2];
00402       MinMax_over[0] = lmin; 
00403       MinMax_over_dist[0] = lmin_dist;  
00404       MinMax_over[1] = lmax;
00405       MinMax_over_dist[1] = lmax_dist;
00406    }
00407    
00408    SUMA_RETURN(YUP);
00409 }
00410 
00411 /*!
00412    \param MinMax (float) 2x1: MinMax[0] minimum under node
00413                               MinMax_dist[1] maximum under node
00414    \param MinMax_dist (float) 2x1: MinMax_dist[0] distance of minimum under node
00415                                    MinMax_dist[1] distance of maximum under node
00416    \param MinMax_over (float) 2x1: MinMax_over[0] minimum over node
00417                                    MinMax_over[1] maximum over node
00418                                    (set this one to NULL if you do not want to search over the node.
00419                                    In that case all the _over (or above like Means[2]) values will
00420                                    be undetermined.
00421    \param MinMax_over_dist (float) 2x1: MinMax_over_dist[0] distance of minimum over node
00422                                         MinMax_over_dist[1] distance of maximum over node
00423    \param Means (float *) 3x1 :  Means[0] value at node, 
00424                                  Means[1] mean value below node (only for vals > Opt->t),
00425                                  Means[2] mean value above node  (only for vals > Opt->t)
00426 */
00427 int SUMA_Find_IminImax_Avg (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, 
00428                         int ni, 
00429                         float *MinMax, float *MinMax_dist , 
00430                         float *MinMax_over, float *MinMax_over_dist,
00431                         float *Means, 
00432                         float *undershish, float *overshish, 
00433                         int *dvecind_under, int *dvecind_over, int ShishMax)
00434 {
00435    static char FuncName[]={"SUMA_Find_IminImax_Avg"};
00436    float d1, d2, travdir[3], d1t, d2t, lmin, lmax, lmin_dist, lmax_dist;
00437    float t2 = Opt->t2, tm = Opt->tm, t = Opt->t; 
00438    THD_fvec3 ncoord, ndicom;
00439    THD_ivec3 nind3;
00440    int nind, istep, istepmax, istep2max, nxx, nxy, nMeans[3], stopint;
00441    int N_vtnmax = 500, N_vtn, vtn[N_vtnmax]; /* vector of triangles incident triangles to node ni */
00442    
00443    SUMA_ENTRY;
00444    
00445    d1 = Opt->d1; d2 = d1/2.0; 
00446    
00447    vtn[N_vtnmax-1] = -1; /* flag to check for overruns */
00448    Means[0] = Means[1] = Means[2] = 0.0;
00449    nMeans[0] = nMeans[1] = nMeans[2] = 0;
00450    lmin = Opt->tm;
00451    lmin_dist = 0.0;
00452    lmax_dist = 0.0;
00453    lmax = Opt->t;
00454    nxx = DSET_NX(Opt->in_vol);
00455    nxy = DSET_NX(Opt->in_vol) * DSET_NY(Opt->in_vol);
00456    istep = 0; istepmax = (int)(ceil((double)d1/Opt->travstp)); istep2max = (int)(ceil((double)d2/Opt->travstp));
00457    stopint = 0;
00458    while (istep <= istepmax) {
00459       /* calculate offset */
00460       travdir[0] = - istep * Opt->travstp * SO->NodeNormList[3*ni]; travdir[1] = -istep * Opt->travstp * SO->NodeNormList[3*ni+1]; travdir[2] = -istep * Opt->travstp * SO->NodeNormList[3*ni+2]; 
00461       
00462       /* find the set of triangles incident to node ni */
00463       if (!SUMA_Get_NodeIncident(ni, SO, vtn, &N_vtn)) {
00464           SUMA_SL_Err("Failed to find incident triangles.\nDecidement, ca va tres mal.\n");
00465           SUMA_RETURN(NOPE);
00466       }
00467       if (vtn[N_vtnmax-1] != -1) {
00468          SUMA_SL_Err("Way too many incident triangles. Memory corruption likely.");
00469          SUMA_RETURN(NOPE);
00470       }
00471       /* for each triangle, find the voxels that it intersects */
00472       
00473       /* get 1d coord of point */
00474       ndicom.xyz[0] = SO->NodeList[3*ni] + travdir[0] ; ndicom.xyz[1] = SO->NodeList[3*ni+1]+ travdir[1]; ndicom.xyz[2] = SO->NodeList[3*ni+2]+ travdir[2]; 
00475       ncoord = THD_dicomm_to_3dmm(Opt->in_vol, ndicom);
00476       nind3 = THD_3dmm_to_3dind(Opt->in_vol, ncoord);
00477       nind = nind3.ijk[0] + nind3.ijk[1] * nxx + nind3.ijk[2] * nxy;
00478       #if 0
00479       if (ni == Opt->NodeDbg) {
00480          fprintf(SUMA_STDERR, "%s: Node %d\n"
00481                               " nind3 = [%d %d %d] voxVal = %.3f\n", 
00482                               FuncName, Opt->NodeDbg, nind3.ijk[0], nind3.ijk[1], nind3.ijk[2],
00483                               Opt->dvec[nind]);
00484       }
00485       #endif
00486       if (undershish && istep < ShishMax) undershish[istep] = Opt->dvec[nind];   /* store values under node */
00487       
00488       /* track the brain mean, only if the value is above brain non-brain threshold
00489       stop inegrating as soon as you hit nonbrain */
00490       if (Opt->dvec[nind] > Opt->t && !stopint) { Means[1] += Opt->dvec[nind]; ++ nMeans[1]; }
00491       else stopint = 1;
00492       if (dvecind_under) dvecind_under[istep] = nind;
00493         
00494       /* find local min */ 
00495       /* lmin = SUMA_MIN_PAIR(lmin, Opt->dvec[nind]); */
00496       if (lmin > Opt->dvec[nind]) { lmin = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmin_dist); }
00497       
00498       if (istep <= istep2max) {
00499          if (lmax < Opt->dvec[nind]) { lmax = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmax_dist); }  
00500       }
00501       
00502       ++istep;
00503    }
00504    
00505    if (istep < ShishMax) { undershish[istep] = -1; if (dvecind_under) dvecind_under[istep] = -1; }/* crown the shish */
00506    
00507    Means[1] /= (float)nMeans[1];
00508    MinMax[0] = SUMA_MAX_PAIR(t2, lmin);
00509    MinMax_dist[0] = lmin_dist;  
00510    MinMax[1] = SUMA_MIN_PAIR(tm, lmax); 
00511    MinMax_dist[1] = lmax_dist;
00512    
00513    /* search for a minimum overhead */
00514    if (MinMax_over) {
00515       lmin = Opt->tm;
00516       lmin_dist = 0.0;
00517       lmax_dist = 0.0;
00518       lmax = Opt->t;
00519       istep = 0; istepmax = (int)(ceil((double)Opt->d4/Opt->travstp));
00520       stopint = 0;
00521       while (istep <= istepmax) {
00522          /* calculate offset */
00523          travdir[0] =  istep * Opt->travstp * SO->NodeNormList[3*ni]; travdir[1] = istep * Opt->travstp * SO->NodeNormList[3*ni+1]; travdir[2] = istep * Opt->travstp * SO->NodeNormList[3*ni+2]; 
00524 
00525          /* get 1d coord of point */
00526          ndicom.xyz[0] = SO->NodeList[3*ni] + travdir[0] ; ndicom.xyz[1] = SO->NodeList[3*ni+1]+ travdir[1]; ndicom.xyz[2] = SO->NodeList[3*ni+2]+ travdir[2]; 
00527          ncoord = THD_dicomm_to_3dmm(Opt->in_vol, ndicom);
00528          nind3 = THD_3dmm_to_3dind(Opt->in_vol, ncoord);
00529          nind = nind3.ijk[0] + nind3.ijk[1] * nxx + nind3.ijk[2] * nxy;
00530          #if 0
00531          if (ni == Opt->NodeDbg) {
00532             fprintf(SUMA_STDERR, "%s: Node %d\n"
00533                                  " nind3 = [%d %d %d] voxVal = %.3f\n", 
00534                                  FuncName, Opt->NodeDbg, nind3.ijk[0], nind3.ijk[1], nind3.ijk[2],
00535                                  Opt->dvec[nind]);
00536          }
00537          #endif
00538          if (!istep) { /* get the value at the node */
00539             Means[0] = Opt->dvec[nind];
00540          }
00541          
00542          if (overshish && istep < ShishMax) overshish[istep] = Opt->dvec[nind];   /* store values under node */
00543          
00544          if (Opt->dvec[nind] > Opt->t && !stopint) { Means[2] += Opt->dvec[nind]; ++ nMeans[2]; } 
00545          else stopint = 1;
00546          if (dvecind_over) dvecind_over[istep] = nind;
00547          
00548          /* find local min and max*/ 
00549          if (lmin > Opt->dvec[nind]) { lmin = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmin_dist); }
00550          if (lmax < Opt->dvec[nind]) { lmax = Opt->dvec[nind]; SUMA_NORM_VEC(travdir, 3, lmax_dist); }  
00551          
00552          ++istep;
00553       }
00554       
00555       if (istep < ShishMax) { overshish[istep] = -1; if (dvecind_over) dvecind_over[istep] = -1; }/* crown the shish */
00556 
00557       Means[2] /= (float)nMeans[2];
00558       MinMax_over[0] = lmin; 
00559       MinMax_over_dist[0] = lmin_dist;  
00560       MinMax_over[1] = lmax;
00561       MinMax_over_dist[1] = lmax_dist;
00562    }
00563    
00564    SUMA_RETURN(YUP);
00565 }
00566 
00567 /*! 
00568    \brief creates a crude mask of the brain's limit based on the surface of the skull
00569    mask is crude, for purposes of keeping BrainWrap from leaking out when there are 
00570    very strong shading artifacts
00571 */
00572 int SUMA_SkullMask (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs)
00573 {
00574    static char FuncName[]={"SUMA_SkullMask"};
00575    int it=0, in=0, Npass, Done, ShishMax, i_diffmax, kth_buf, N_it;
00576    float *n=NULL, *undershish = NULL, diffmax;
00577    float *tmpptr, *NewCoord = NULL, f3, f4, *dsmooth = NULL;
00578    float *refNodeList = NULL, sksearch = 15;
00579    float U[3], Un;
00580    int  nx, nxy, n_stp;
00581    SUMA_Boolean DoDbg=NOPE, Send_buf;
00582    SUMA_Boolean LocalHead = NOPE;
00583    
00584    SUMA_ENTRY;
00585    
00586    if (Opt->debug > 2) LocalHead = YUP;
00587 
00588    nx = DSET_NX(Opt->in_vol);
00589    nxy = nx * DSET_NY(Opt->in_vol);   
00590    
00591    kth_buf = cs->kth;
00592    Send_buf = cs->Send;
00593    if (!Opt->send_hull) {
00594       cs->Send = NOPE;
00595    }   
00596 
00597    NewCoord = (float *)SUMA_malloc(3*SO->N_Node*sizeof(float));
00598    if (!NewCoord) {
00599       SUMA_SL_Crit("Could not allocate for temp vector.");
00600       SUMA_RETURN(NOPE);
00601    }
00602    
00603    ShishMax = (int)(SUMA_MAX_PAIR(Opt->d1, Opt->d4) / Opt->travstp * 1.2);
00604    undershish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
00605    if (!undershish )  {
00606       SUMA_SL_Crit("Failed to allocate");
00607       SUMA_RETURN(NOPE);
00608    }
00609    
00610    
00611    /* search for a distance for the biggest negative gradient over a distance of 2 cm */
00612    sksearch = 15; /* maximum search into skull */
00613    n_stp = (int)(ceil((double)sksearch/Opt->travstp));
00614    Npass = 0;  Done = 0; N_it = 1;
00615    do {
00616       if (Opt->debug > 1) fprintf(SUMA_STDERR,"%s: Gradient Pass #%d\n", FuncName, Npass);
00617       for (it=0; it < N_it; ++it) {
00618          for (in=0; in<SO->N_Node; ++in) {
00619             n = &(SO->NodeList[3*in]);
00620 
00621             /* using normals for direction, NORMALS SHOULD POINT OUTWARDS...*/
00622             if (n[2] > -45) { /* Do not touch nodes on bottom */
00623                U[0] = -SO->NodeNormList[3*in]; U[1] = -SO->NodeNormList[3*in+1]; U[2] = -SO->NodeNormList[3*in+2]; 
00624                SUMA_ONE_SHISH_PLEASE(n, U, Opt->in_vol, Opt->dvec, Opt->travstp, n_stp, nx, nxy, undershish, ShishMax);
00625                SUMA_MAX_NEG_SHISH_JUMP(undershish, diffmax, i_diffmax);
00626             } else {
00627                i_diffmax = 0;
00628             }
00629             NewCoord[3*in+0] = n[0] + i_diffmax * Opt->travstp * U[0]; 
00630             NewCoord[3*in+1] = n[1] + i_diffmax * Opt->travstp * U[1]; 
00631             NewCoord[3*in+2] = n[2] + i_diffmax * Opt->travstp * U[2]; 
00632             
00633             DoDbg = NOPE;
00634          } /* loop over number of nodes */
00635    
00636             
00637          /* swap vectors */
00638          tmpptr = SO->NodeList;
00639          SO->NodeList = NewCoord; 
00640          NewCoord = tmpptr;
00641          tmpptr = NULL;
00642          
00643          cs->kth = 1;  /*make sure all gets sent at this stage */
00644          dsmooth = SUMA_Taubin_Smooth( SO, NULL,
00645                                     0.6307, -.6732, SO->NodeList,
00646                                     8, 3, SUMA_ROW_MAJOR, dsmooth, cs, NULL);    
00647          memcpy((void*)SO->NodeList, (void *)dsmooth, SO->N_Node * 3 * sizeof(float));
00648          cs->kth = kth_buf; 
00649          
00650          /* recalculate surface normals */
00651          SUMA_RECOMPUTE_NORMALS(SO); 
00652          if (cs->Send) {
00653             if (!SUMA_SendToSuma (SO, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
00654                SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
00655             }
00656          }
00657 
00658       } /* loop over number of iterations */
00659       Done = 1;
00660    } while (!Done);         
00661    
00662    cs->Send = Send_buf;
00663 
00664    if (undershish) SUMA_free(undershish); undershish = NULL;
00665    if (refNodeList) SUMA_free(refNodeList); refNodeList = NULL;
00666    if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL;
00667    if (NewCoord) SUMA_free(NewCoord); NewCoord = NULL;
00668    SUMA_RETURN(YUP);
00669 }
00670 
00671 /*! 
00672    Main function to stretch a sphere to fit the brain. 
00673    A modification of BET's algorithm.
00674 */
00675 int SUMA_StretchToFitLeCerveau (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs)
00676 {
00677    static char FuncName[]={"SUMA_StretchToFitLeCerveau"};
00678    int it=0,it0 = 0, in=0, ii, Npass, Done, ShishMax, i_diffmax, kth_buf, nit, Stage;
00679    float mnc[3], s[3], sn[3], st[3], dp, threshclip, lztfac,
00680          nsn, rmin, rmax, E, F, su1, su2, su3, su4,  
00681          r, u1[3], u2[3], u3[3], u[3],
00682          tm=0.0, t2 = 0.0, t98 = 0.0, Imin, Imax, l=0.0, MinMax[2], 
00683          MinMax_dist[2],MinMax_over[2], MinMax_over_dist[2], Means[3],
00684          *n=NULL, tb = 0.0, tbe = 0.0, t = 0.0, Imin_over=0.0, Imin_dist_over=0.0,
00685          *overshish = NULL, *undershish = NULL, diffmax, *touchup=NULL, *a, P2[2][3], *norm;
00686    float *tmpptr, *NewCoord = NULL, f3, f4, lZt, *dsmooth = NULL;
00687    float *refNodeList = NULL, *OrigNodeList=NULL;
00688    double MaxExp;
00689    int keepgoing=0, N_troub, past_N_troub;
00690    double pastarea, curarea, darea;
00691    char cbuf;
00692    FILE *OutNodeFile = NULL;
00693    SUMA_Boolean DoDbg=NOPE;
00694    SUMA_Boolean LocalHead = NOPE;
00695    
00696    SUMA_ENTRY;
00697    
00698    InteractiveQuit = 0;
00699    
00700    if (Opt->debug > 2) LocalHead = YUP;
00701 
00702    kth_buf = cs->kth;
00703    
00704    t2 = Opt->t2; t98 = Opt->t98; tm = Opt->tm; t = Opt->t;
00705    threshclip = ((Opt->t98 - Opt->t2)/5.0);
00706    
00707    if (LocalHead) fprintf(SUMA_STDERR,"%s: t2 = %f, t = %f, tm = %f, t98 = %f\n", FuncName, t2, t, tm, t98);
00708    
00709    
00710    NewCoord = (float *)SUMA_malloc(3*SO->N_Node*sizeof(float));
00711    if (!NewCoord) {
00712       SUMA_SL_Crit("Could not allocate for temp vector.");
00713       SUMA_RETURN(NOPE);
00714    }
00715    rmin = 3.33; rmax = 10.0;  E = (1/rmin + 1/rmax)/2; F = 6/(1/rmin - 1/rmax); 
00716    su1 = Opt->su1; 
00717    
00718    if (Opt->NodeDbg >= 0) {
00719       char fname[100];
00720       sprintf(fname,"dbg_%d.1D", Opt->NodeDbg);
00721       OutNodeFile = fopen(fname,"w");
00722       if (!OutNodeFile) {
00723          SUMA_S_Err("Failed to open debug file");
00724          SUMA_RETURN(NOPE);
00725       }
00726       fprintf(OutNodeFile, "#0 Node id (in)\n"
00727                            "#1 Iteration (it)\n"
00728                            "#2 Seg length (l)\n"
00729                            "#3-5 Imin, Imax, tb\n"
00730                            "     MinMaxdist, 2 col\n"
00731                            "     MinMaxover, 2 col\n"
00732                            "     MinMaxover_dist, 2 col\n"
00733                            "     Means (at, under, over), 3 col\n"
00734                            "#6-8 Node coord (n)\n"
00735                            "#9 Norm . Disp (dp)\n"
00736                            "#10-12 Mean Neighb Coord (mnc)\n"
00737                            "#13-15 Normal direction (  )\n"
00738                            "#16-18 Disp (s)\n"
00739                            "#19-21 DispT (st)\n"
00740                            "#22-24 DispN (sn)\n"
00741                            "#25-27 constants (su1 su2 su3)\n"
00742                            "#28-30 update vector (u)\n"
00743                            "#31 t2\n"
00744                            "#32 t\n"
00745                            "#33 tm\n"
00746                            "#34 t98\n"
00747                            "#\n"
00748                            "#t2 = %f, t = %f, tm = %f, t98 = %f\n", t2, t, tm, t98);
00749    }
00750    
00751    ShishMax = (int)(SUMA_MAX_PAIR(Opt->d1, Opt->d4) / Opt->travstp * 1.2);
00752    overshish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
00753    undershish = (float *)SUMA_malloc(sizeof(float)*ShishMax); 
00754    OrigNodeList = (float*)SUMA_malloc(sizeof(float)*3*SO->N_Node);
00755    if (!OrigNodeList || !overshish || !undershish) {
00756       SUMA_SL_Crit("Failed to allocate");
00757       SUMA_RETURN(NOPE);
00758    }
00759    memcpy((void*)OrigNodeList, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float));
00760    
00761    pastarea = 0.0; curarea = 0.0; darea = 0.0, past_N_troub = 0;
00762    Npass = 0;   Done = 0; nit = Opt->N_it; it0 = 0;
00763    do {
00764       Stage = 0; /* reset stage */
00765       if (Opt->debug > 1) fprintf(SUMA_STDERR,"%s: Pass #%d\n", FuncName, Npass);
00766       MaxExp = 0.0; /* maximum expansion of any node */
00767       for (it=it0; it < nit; ++it) {
00768          SUMA_MEAN_SEGMENT_LENGTH(SO, l);
00769          if (LocalHead && !(it % 50)) fprintf (SUMA_STDERR,"%s: Iteration %d, l = %f . SO->Center = %f, %f, %f...\n", 
00770             FuncName, it, l, SO->Center[0], SO->Center[1], SO->Center[2]);
00771          if (Opt->var_lzt) {
00772             if (it <= Opt->N_it) lztfac = SUMA_MAX_PAIR(0.2, (1.2* (float)it / (float)Opt->N_it));  /* make things expand quickly in the beginning, to help escape large sections of csf close to outer surface, then tighten grip towards the end*/
00773             else lztfac = 1.2;
00774          } else lztfac = 1.0;
00775          MaxExp = 0.0; /* maximum expansion of any node */
00776          for (in=0; in<SO->N_Node; ++in) {
00777             n = &(SO->NodeList[3*in]);
00778                SUMA_MEAN_NEIGHB_COORD(SO, in, mnc);
00779                s[0] = mnc[0] - n[0]; s[1] = mnc[1] - n[1]; s[2] = mnc[2] - n[2];
00780                SUMA_DOTP_VEC(s, &(SO->NodeNormList[3*in]), dp, 3, float, float);  /* dp is the dot product of s with the normal at in */ 
00781                SUMA_SCALE_VEC(&(SO->NodeNormList[3*in]), sn, dp, 3, float, float); /* sn is the component of s along normal */
00782                SUMA_NORM_VEC(sn, 3, nsn); /* norm of sn */
00783                SUMA_SUB_VEC(s, sn, st, 3, float, float, float);  /* st is the tangential component of s along normal */
00784                SUMA_SCALE_VEC(st, u1, su1, 3, float, float); /* u1 = st * su1 */
00785 
00786                r = ( l * l ) / (2.0 * nsn);
00787                su2 = ( 1 + tanh(F*(1/r - E)))/2.0;
00788                SUMA_SCALE_VEC(sn, u2, su2, 3, float, float); /* u2 = sn * su2 */
00789 
00790                SUMA_Find_IminImax(SO, Opt, in, MinMax, MinMax_dist, MinMax_over, MinMax_over_dist, Means, undershish, overshish, NULL, NULL, ShishMax);  
00791                Imin = MinMax[0]; Imax = MinMax[1];
00792 
00793                if (n[2] - SO->Center[2] < -25) { lZt = SUMA_MAX_PAIR (Opt->bot_lztclip, (Opt->ztv[in] * lztfac)); } /* clip lZt at no less than 0.5 for lowest sections, 
00794                                                                                                             otherwise you might go down to the dumps */
00795                else { lZt = Opt->ztv[in] * lztfac;  }
00796 
00797                if (lZt > 1) lZt = 0.999;
00798 
00799                tb = (Imax - t2) * lZt + t2; 
00800                f3 = 2.0 * (Imin - tb) / (Imax - t2 );
00801 
00802                su3 = Opt->ExpFrac * l * f3 ;
00803 
00804                if (Opt->UseExpansion) {
00805                   tbe = (MinMax_over[1] - t2) * lZt + t2; 
00806                   f4 = 2.0 * (MinMax_over[0] - tbe) / (MinMax_over[1] - t2 );
00807                   su4 = Opt->ExpFrac * l * f4 ; if (su4 < 0) su4 = 0; /* too tight expanding otherwise */
00808                } else su4 = 0;
00809 
00810                /* work the eyes */
00811                if (Opt->NoEyes) {
00812                   if (it > 0.1 * Opt->N_it) { 
00813                      if (SUMA_IS_EYE_ZONE(n,SO->Center)) { /* (n[1] - SO->Center[1]) < -10 && (n[2] - SO->Center[2]) < 0 area with the eyes */
00814                         /* if the mean over is larger than the mean below, go easy, likely to hit the eyes */
00815                         if (!Opt->Stop[in]) {
00816                            if (Opt->dbg_eyenodes) { 
00817                               int kk;
00818                               fprintf (Opt->dbg_eyenodes,"%d\t %f %.3f ", in, Means[2], Means[2]/Means[1]); 
00819                               for (kk=0; kk<ShishMax; ++kk) if (overshish[kk] >= 0) fprintf (Opt->dbg_eyenodes,"%.3f ", overshish[kk]);
00820                               fprintf (Opt->dbg_eyenodes,"\n");
00821                            }
00822                            if (Means[2] > 1.2*Means[1]) {
00823                               SUMA_MAX_SHISH_JUMP(overshish, diffmax, i_diffmax, threshclip);
00824                               Opt->Stop[in] = overshish[i_diffmax];
00825                               if (0 && LocalHead) fprintf (SUMA_STDERR, "%s: Stopping node %d, at values > %f\n", FuncName, in, Opt->Stop[in]);
00826                            } else if (Means[0] >= Opt->t98) {
00827                               Opt->Stop[in] = Opt->t98;
00828                               if (0 && LocalHead) fprintf (SUMA_STDERR, "%s: Stopping node %d, at values > %f\n", FuncName, in, Opt->Stop[in]);
00829                            }
00830                         }
00831                      }
00832                   }
00833                   if (Opt->Stop[in]) {
00834                      if (Means[0] >= Opt->Stop[in]) { su3 = 0; su4 = 0; }
00835                   }
00836                }
00837                /* freezing: a node that has been frozen, freezing may be used during final touchup*/
00838                if (Opt->Stop[in] < 0) { 
00839                   su3 = 0; su4 = 0;   
00840                   if (in == Opt->NodeDbg) fprintf(SUMA_STDERR,"%s:\nNode %d has been frozen\n", FuncName, in);
00841                }
00842                
00843                u[0] = su1 * st[0] + su2 * sn[0] + ( su3 + su4 ) * SO->NodeNormList[3*in]   ;
00844                u[1] = su1 * st[1] + su2 * sn[1] + ( su3 + su4 ) * SO->NodeNormList[3*in+1] ;
00845                u[2] = su1 * st[2] + su2 * sn[2] + ( su3 + su4 ) * SO->NodeNormList[3*in+2] ; 
00846                if (0 && (in == Opt->NodeDbg)) { 
00847                   fprintf(SUMA_STDERR, "%s: Node %d, iter %d, l %.6f\n", FuncName, in, it, l);
00848                   fprintf(SUMA_STDERR, " MinMaxTb = [%.6f %.6f %.6f]\n", Imin, Imax, tb);
00849                   fprintf(SUMA_STDERR, " MinMaxdist=[%.6f %.6f ]\n", MinMax_dist[0], MinMax_dist[1]);
00850                   fprintf(SUMA_STDERR, " MinMaxover     =[%.6f %.6f ]\n", MinMax_over[0], MinMax_over[1]);
00851                   fprintf(SUMA_STDERR, " MinMaxover_dist=[%.6f %.6f ]\n", MinMax_over_dist[0], MinMax_over_dist[1]); 
00852                   fprintf(SUMA_STDERR, " Means (at, under, over) = [%.6f %.6f %.6f]\n", Means[0], Means[1], Means[2]); 
00853                   fprintf(SUMA_STDERR, " Coord(n )= [%.6f %.6f %.6f], dp = %f\n",   n[0],  n[1],  n[2], dp); 
00854                   fprintf(SUMA_STDERR, " Neigh(mnc)=[%.6f %.6f %.6f]\n",   mnc[0],  mnc[1],  mnc[2]); 
00855                   fprintf(SUMA_STDERR, " Norm (  )= [%.6f %.6f %.6f]\n",   SO->NodeNormList[3*in], SO->NodeNormList[3*in+1],SO->NodeNormList[3*in+2]); 
00856                   fprintf(SUMA_STDERR, " Disp (s )= [%.6f %.6f %.6f]\n",   s[0],  s[1],  s[2]); 
00857                   fprintf(SUMA_STDERR, "      (st)= [%.6f %.6f %.6f]\n",  st[0], st[1], st[2]); 
00858                   fprintf(SUMA_STDERR, "      (sn)= [%.6f %.6f %.6f]\n",  sn[0], sn[1], sn[2]); 
00859                   fprintf(SUMA_STDERR, "       su = [%.6f %.6f %.6f]\n",  su1, su2, su3); 
00860                   fprintf(SUMA_STDERR, "        u = [%.6f %.6f %.6f]\n",  u[0], u[1], u[2]); 
00861                   if (OutNodeFile) {
00862                      fprintf(OutNodeFile, " %d %d  %.6f"
00863                                           " %.6f %.6f %.6f"
00864                                           " %.6f %.6f"
00865                                           " %.6f %.6f"
00866                                           " %.6f %.6f"
00867                                           " %.6f %.6f %.6f"
00868                                           " %.6f %.6f %.6f %.6f"
00869                                           " %.6f %.6f %.6f"
00870                                           " %.6f %.6f %.6f"
00871                                           " %.6f %.6f %.6f"
00872                                           " %.6f %.6f %.6f"
00873                                           " %.6f %.6f %.6f"
00874                                           " %.6f %.6f %.6f"
00875                                           " %.6f %.6f %.6f"
00876                                           " %.6f %.6f %.6f %.6f"
00877                                           "\n"
00878                                           , in, it, l
00879                                           , Imin, Imax, tb
00880                                           , MinMax_dist[0], MinMax_dist[1]
00881                                           , MinMax_over[0], MinMax_over[1]
00882                                           , MinMax_over_dist[0], MinMax_over_dist[1]
00883                                           , Means[0], Means[1], Means[2]
00884                                           ,   n[0],  n[1],  n[2], dp
00885                                           ,   mnc[0],  mnc[1],  mnc[2]
00886                                           ,SO->NodeNormList[3*in], SO->NodeNormList[3*in+1],SO->NodeNormList[3*in+2]
00887                                           ,s[0],  s[1],  s[2]
00888                                           ,  st[0], st[1], st[2]
00889                                           ,  sn[0], sn[1], sn[2]
00890                                           ,  su1, su2, su3
00891                                           ,  u[0], u[1], u[2]
00892                                           , t2, t, tm, t98);
00893                   }
00894 
00895                }
00896                NewCoord[3*in] = n[0] + u[0]; NewCoord[3*in+1] = n[1] + u[1]; NewCoord[3*in+2] = n[2] + u[2];       
00897                MaxExp = SUMA_MAX_PAIR (MaxExp, ( sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]) ) );
00898                DoDbg = NOPE;
00899             
00900          } /* loop over number of nodes */
00901          /* swap vectors */
00902             tmpptr = SO->NodeList;
00903             SO->NodeList = NewCoord; 
00904             NewCoord = tmpptr;
00905             tmpptr = NULL;
00906 
00907          if (Opt->NNsmooth) {
00908             if (!(it % Opt->smootheach) && it > 0 && it < 0.75*Opt->N_it) {  /* cannot filter too close to convergence */
00909                SUMA_WRAP_BRAIN_SMOOTH_NN(Opt->NNsmooth, dsmooth, refNodeList);
00910             }
00911          }
00912 
00913          /* recalculate surface normals */
00914          SUMA_RECOMPUTE_NORMALS(SO); 
00915          if (cs->Send) {
00916             if (!SUMA_SendToSuma (SO, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
00917                SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
00918             }
00919          }
00920 
00921       } /* loop over number of iterations */
00922       ++Stage;
00923       if (Opt->DemoPause == SUMA_3dSS_INTERACTIVE) {
00924          fprintf (SUMA_STDERR,"3dSkullStrip Interactive: \n"
00925                               "Increasing number of iterations to reach minimal expansion.\n"
00926                               "Do you want to (C)ontinue, (P)ass or (S)ave this? [C] ");
00927          cbuf = SUMA_ReadCharStdin ('c', 0,"csp");
00928          fprintf (SUMA_STDERR,"%c\n", cbuf);
00929          switch (cbuf) {
00930             case 's':   
00931                fprintf (SUMA_STDERR,"Will check surface, save mask and exit.\n");
00932                InteractiveQuit = 1;
00933                goto CLEAN_RETURN;
00934                break;
00935             case 'p':
00936                fprintf (SUMA_STDERR,"Passing this stage \n");
00937                goto CLEAN_RETURN;
00938                break;
00939             case 'c':
00940                fprintf (SUMA_STDERR,"Continuing with stage.\n");
00941                break;
00942          }                 
00943       }
00944       if (Stage == 1) { /* if the surface is still growing, keep going */
00945          if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("About to be in stage 1"); }
00946          if (LocalHead) fprintf (SUMA_STDERR,"%s: \n In stage 1\n", FuncName);
00947          if (MaxExp > 0.5) {
00948             /* Now, if you still have expansion, continue */
00949             if (!pastarea) { /* first time around, calculate area */
00950                pastarea = SUMA_Mesh_Area(SO, NULL, -1);
00951                if (LocalHead) fprintf (SUMA_STDERR,"%s: \n Stage1: pastarea = %f\n", FuncName, pastarea);
00952                keepgoing = 1;
00953             }else {
00954                curarea = SUMA_Mesh_Area(SO, NULL, -1);
00955                darea = ( curarea - pastarea ) / pastarea; 
00956                if (SUMA_ABS(darea) > 0.1) {
00957                   keepgoing = 1;
00958                } else { 
00959                   keepgoing = 0;
00960                }
00961                pastarea = curarea;
00962             }
00963             if (keepgoing) {
00964                it0 = nit; nit = nit + (int)(Opt->N_it/2.5);
00965                if (LocalHead) fprintf (SUMA_STDERR,"%s: \n Stage1: MaxExp = %f, darea = %f, going for more...\n", FuncName, MaxExp, darea);
00966                Done = 0;
00967             } else {
00968                if (LocalHead) fprintf (SUMA_STDERR,"%s: \n Stage1: satiated, small area differential.\n", FuncName);
00969                ++Stage;
00970             }
00971          } else {
00972             if (LocalHead) fprintf (SUMA_STDERR,"%s: \n Stage1: satiated, low MaxExp\n", FuncName);
00973             ++Stage;
00974          }
00975       }
00976       if (Stage == 2) {
00977          if (0) {
00978             if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("About to be in stage 2"); }
00979             touchup = SUMA_Suggest_Touchup_Grad(SO, Opt, 4.0, cs, &N_troub);
00980             if (!touchup) SUMA_SL_Warn("Failed in SUMA_Suggest_Touchup");
00981          } else {
00982             SUMA_LH( "Bypasing, SUMA_Suggest_Touchup_Grad is bad,"
00983                      " need serious 3D edge detection. "
00984                      " And cone shaped search zone for the gradient crossing  ... ");
00985             N_troub = 0;
00986          }
00987          
00988          if (!N_troub) { 
00989             /* No touchup, done */
00990             ++Stage;
00991          } else {
00992             if (touchup) SUMA_free(touchup); touchup = NULL;
00993             ++Stage;  
00994          }
00995       }
00996       if (Stage == 3) {
00997          if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("About to be in stage 3"); }
00998          /* see if there is room for extension */
00999          touchup = SUMA_Suggest_Touchup(SO, Opt, 4.0, cs, &N_troub);
01000          if (!touchup) SUMA_SL_Warn("Failed in SUMA_Suggest_Touchup");
01001          if (!N_troub) { 
01002             /* No touchup, done */
01003             Done = 1;
01004          } else {
01005             /* reduce tightness where expansion is needed */
01006             if (LocalHead) {
01007                fprintf(SUMA_STDERR,"%s:\n reducing tightness, applying touchup\n", FuncName);
01008             }
01009             for (in=0; in<SO->N_Node; ++in) {
01010                if (touchup[in]) {
01011                   if (Opt->NodeDbg == in) fprintf(SUMA_STDERR,"%s: Acting on node %d, touchup[in] %f, past shrink_fac = %f\n", 
01012                                                          FuncName, in, touchup[in], Opt->ztv[in]);                  
01013                   if (touchup[in] < 1) Opt->ztv[in] *= 0.8;
01014                   else if (touchup[in] < 2) Opt->ztv[in] *= 0.7;
01015                   else if (touchup[in] < 3) Opt->ztv[in] *= 0.6;
01016                   else if (touchup[in] < 4) Opt->ztv[in] *= 0.5;
01017                   else Opt->ztv[in] *= 0.4;
01018                }
01019             }
01020             it0 = nit; nit = nit + (int)(Opt->N_it/2.5);
01021             if (LocalHead) fprintf (SUMA_STDERR,"%s: \n Stage3: %d troubled nodes, going for more...\n", FuncName, N_troub);
01022             Done = 0;
01023             
01024             if (!past_N_troub) { past_N_troub = 0; }
01025             else {
01026                float dtroub;
01027                dtroub = (float)(past_N_troub - N_troub)/(float)past_N_troub;
01028                if (SUMA_ABS(dtroub) > 0.2) {
01029                   Done = 0; /* continue */
01030                } else {
01031                   Done = 1; /* done deal */
01032                }
01033                past_N_troub = N_troub;
01034             }
01035          }
01036          
01037          if (touchup) SUMA_free(touchup); touchup = NULL;
01038       }
01039       if (Stage > 3) {
01040          SUMA_SL_Err("Stage number invalid");
01041          Done = 1;
01042       }
01043       /* put a limit */
01044       if ( (nit > 3 * Opt->N_it) && !Done) {
01045          SUMA_LH("Funding limit reached. Abandoning improvements");
01046          Done = 1;
01047       }
01048    } while (!Done);
01049       
01050    CLEAN_RETURN:
01051    if (undershish) SUMA_free(undershish); undershish = NULL;
01052    if (overshish) SUMA_free(overshish); overshish = NULL;
01053    if (OrigNodeList) SUMA_free(OrigNodeList); OrigNodeList = NULL;
01054    if (refNodeList) SUMA_free(refNodeList); refNodeList = NULL;
01055    if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL;
01056    if (OutNodeFile) fclose(OutNodeFile); OutNodeFile = NULL;
01057    if (NewCoord) SUMA_free(NewCoord); NewCoord = NULL;
01058    SUMA_RETURN(YUP);
01059 }
01060 
01061 #define MSK_DBG 0
01062 #define Meth1   /* SLOW METHOD, undefine it to have an even slower method! */
01063 /*!
01064    \brief First pass at finding voxels inside a surface. Slow baby, slow.
01065    The first approach traces a ray parallel to the X axis from each voxel 
01066    and counts the number of intersections, in the direction of the ray, with
01067    the surface. Odd = in, Even = out. That is Meth1.
01068    The other method, projects all triangles to the YZ plane and then finds out 
01069    if a voxel's projection lies inside the projected triangles. You still have to
01070    do a test to count only the triangles that are intersected in the direction of the ray
01071    In the end, it turns out to be slower than Meth1. Bad
01072    
01073    \sa SUMA_FindVoxelsInSurface
01074 */
01075 byte *SUMA_FindVoxelsInSurface_SLOW (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole) 
01076 {
01077    static char FuncName[]={"SUMA_FindVoxelsInSurface_SLOW"};
01078    byte *isin = NULL, *tmpin = NULL;
01079    int i, N_in, j, k , n, khits, dims[2], N_hits, iii, jjj, niii, ncul;
01080    float *tmpXYZ=NULL, Center[3], MaxDims[3], MinDims[3], aMaxDims, aMinDims, delta_t;
01081    float hdim[3], t0, t1, t2, SOCenter[0], p0[3], p1[3];
01082    SUMA_MT_INTERSECT_TRIANGLE *mti = NULL; 
01083    struct  timeval tti;
01084    int nfound = 0, N_poshits = 0;
01085    int cnt = 0, sgn, sgnref, n1, n2, n3;
01086    float ivec[3] = { 1, 0, 0}, dp = 0.0, cx = 0.0;
01087    SUMA_Boolean LocalHead = NOPE;
01088    #if MSK_DBG
01089       FILE *fdb=fopen("SUMA_FindVoxelsInSurface.1D","w");
01090    #endif
01091    
01092    SUMA_ENTRY;
01093    
01094    SUMA_etime (&tti, 0);
01095    if (LocalHead) {
01096       fprintf(SUMA_STDERR,"%s: Patience...\n", FuncName);
01097    }
01098    N_in = *N_inp = 0;
01099    
01100    /* copy surface coordinates, we're going to ijk land */
01101    tmpXYZ = (float *)SUMA_malloc(SO->N_Node * 3 * sizeof(float));
01102    isin = (byte *)SUMA_malloc(VolPar->nx*VolPar->ny*VolPar->nz * sizeof(byte));
01103    if (!tmpXYZ || !isin) {
01104       SUMA_SL_Crit("Faile to allocate");
01105       SUMA_RETURN(NULL);
01106    }
01107    
01108    memcpy ((void*)tmpXYZ, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float));
01109    
01110    /* transform the surface's coordinates from RAI to 3dfind */
01111    SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, VolPar);
01112    
01113    /* find the new center and bounding box for the surface in the new coordinate system*/
01114    SUMA_MIN_MAX_SUM_VECMAT_COL (tmpXYZ, SO->N_Node, 3, MinDims, MaxDims, SOCenter); 
01115    SOCenter[0] /= SO->N_Node;  
01116    SOCenter[1] /= SO->N_Node;  
01117    SOCenter[2] /= SO->N_Node;  
01118    SUMA_MIN_VEC (MinDims, 3, aMinDims );   
01119    SUMA_MAX_VEC (MaxDims, 3, aMaxDims);    
01120    /* Box half Dim */
01121    for (i=0; i<3; ++i) { hdim[i] = ( MaxDims[i] - MinDims[i]) / 2.0; }
01122    /* Box Center */
01123    for (i=0; i<3; ++i) { Center[i] = MinDims[i] + hdim[i]; }
01124    if (LocalHead) fprintf (SUMA_STDERR,"Center in IJK is: %f %f %f\n"
01125                                        "HlfDim in IJK is: %f %f %f\n"
01126                                        , Center[0], Center[1], Center[2],
01127                                          hdim[0], hdim[1], hdim[2]);
01128    
01129    n = 0; N_in = 0;
01130    for (k=0; k < VolPar->nz; ++k) {
01131       for (j=0; j < VolPar->ny; ++j) {
01132          for (i=0; i < VolPar->nx; ++i) {
01133             isin[n] = 0; /* initialize to out */
01134             /* is voxel center inside bounding box ? */
01135             t0 = hdim[0] - SUMA_ABS(i - Center[0]);   
01136             if (t0 >= 0) {
01137                t1 = hdim[1] - SUMA_ABS(j - Center[1]);
01138                if (t1 >= 0) {
01139                   t2 = hdim[2] - SUMA_ABS(k - Center[2]); 
01140                   if (t2 >= 0)  {
01141                      isin[n] = 1; ++N_in;
01142                   }
01143                }   
01144             }
01145             if (isin[n]) { /* inside bounding box, but is it inside the surface ? */
01146                p0[0] = i; p1[0] = i+1000; p0[1] = p1[1] = j; p0[2] = p1[2] = k; 
01147                #ifdef Meth1
01148                   /* works, but slow as a turtle */
01149                   mti = SUMA_MT_intersect_triangle(p0, p1, tmpXYZ, SO->N_Node, SO->FaceSetList, SO->N_FaceSet, mti);
01150                   if (!(mti->N_poshits % 2)) { /* number of positive direction hits is a multiple of 2 */
01151                      isin[n] = 0; --N_in;
01152                   } 
01153                #else 
01154                   dims[0] = 1; dims[1] = 2; /* rays are along x vector, we will look for intersections with projections on plane y z */
01155                   tmpin = SUMA_isinpoly(p0, tmpXYZ, SO->FaceSetList, SO->N_FaceSet, SO->FaceSetDim, dims, &N_hits, tmpin, NULL);
01156                   N_poshits = 0; sgnref = 0;
01157                   nfound = 0; cnt = 0;
01158                   while (nfound < N_hits) {
01159                      if (tmpin[cnt]) {
01160                         ++nfound; 
01161                         /* face centroid x coordinate, this contraption works as long as dims is [1 2] and as long as
01162                         the triangles making up the mesh are comparable in size to voxel resolution. Otherwise, it is 
01163                         not a particularly pleasing approximation. The proper way would be to do a ray/triangle intersection
01164                         for those hit triangles (use SUMA_MT_isIntersect_Triangle and check for differences in signed distance). 
01165                         But that ends up like the method above with no speed up to write home about NOT WORTH IT*/
01166                         n1 =  SO->FaceSetList[3*cnt]; n2 = SO->FaceSetList[3*cnt+1]; n3 = SO->FaceSetList[3*cnt+2];   \
01167                         cx = (tmpXYZ[3*n1]   + tmpXYZ[3*n2]   + tmpXYZ[3*n3]  )/3; \
01168                         sgn = SUMA_SIGN(cx - i);
01169                         if (nfound == 1) {
01170                            sgnref = sgn;
01171                         }
01172                         /* fprintf(SUMA_STDERR, "%s: nfound %d: dp %f: Ref %d, sgn %d\n", FuncName, nfound, dp, sgnref, sgn); */
01173                         if (sgnref == sgn) ++N_poshits;
01174                      }
01175                      ++cnt;      
01176                   }
01177                   if (N_poshits % 2 == 0) { /* dude's outside */ 
01178                      isin[n] = 0; --N_in;
01179                      /* if (LocalHead) fprintf (SUMA_STDERR,"%d %d %d %d hits\n",i, j, k, N_hits);  */
01180                   }
01181                #endif
01182             }
01183             #if MSK_DBG
01184                if (isin[n]) fprintf(fdb, "%d %d %d %d\n", i, j, k, isin[n]);
01185             #endif
01186             ++n;
01187          }
01188       }               
01189    }         
01190 
01191    *N_inp = N_in;
01192    #if MSK_DBG
01193       fclose(fdb); fdb = NULL;
01194    #endif
01195    
01196    delta_t = SUMA_etime (&tti, 1); 
01197    if (LocalHead) {
01198       fprintf(SUMA_STDERR,"%s: Execution time %f seconds.\n", FuncName, delta_t);
01199    }
01200    
01201    SUMA_free(tmpXYZ); tmpXYZ = NULL;
01202    if (mti) mti = SUMA_Free_MT_intersect_triangle(mti);
01203    if (tmpin) SUMA_free(tmpin); tmpin = NULL;
01204    SUMA_RETURN(isin);
01205 }
01206 
01207 
01208 
01209 /*!
01210    \param NodeIJKlist (float *) the equivalent of SO->NodeList only in i,j,k indices into VolPar's grid.
01211                            In the future, might want to allow for this param to be NULL and have it
01212                            be generated internally from SO->NodeList and VolPar.
01213    \return isin (short *)  1: voxel contains node
01214                            2: voxel is in box containing triangle
01215                            3: voxel is in box containing triangle and is on the side opposite to the normal 
01216                            See SUMA_SURF_GRID_INTERSECT_OPTIONS for the various possible values
01217 
01218    - The first part of this function is the source for function SUMA_GetVoxelsIntersectingTriangle
01219    If you find bugs here, fix them there too. 
01220 
01221 */
01222 short *SUMA_SurfGridIntersect (SUMA_SurfaceObject *SO, float *NodeIJKlist, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole, THD_3dim_dataset *fillmaskset)
01223 {
01224    static char FuncName[]={"SUMA_SurfGridIntersect"};
01225    short *isin=NULL;
01226    byte *ijkmask=NULL, *inmask = NULL, *ijkout = NULL;
01227    float *p1, *p2, *p3, min_v[3], max_v[3], p[3], dist;
01228    float MaxDims[3], MinDims[3], SOCenter[3], dxyz[3];
01229    int nn, nijk, nx, ny, nz, nxy, nxyz, nf, n1, n2, n3, nn3, *voxelsijk=NULL, N_alloc, en;
01230    int N_inbox, nt, nt3, ijkseed = -1, N_in, N_realloc;
01231    byte *fillmaskvec=NULL;
01232    SUMA_Boolean LocalHead = NOPE;
01233    
01234    SUMA_ENTRY;
01235    
01236    if (SO->FaceSetDim != 3 || SO->NodeDim != 3) {
01237       SUMA_SL_Err("SO->FaceSetDim != 3 || SO->NodeDim != 3"); 
01238       SUMA_RETURN(NULL);
01239    }
01240    
01241    nx = VolPar->nx; ny = VolPar->ny; nz = VolPar->nz; nxy = nx * ny; nxyz = nx * ny * nz;
01242    
01243 
01244    isin = (short *)SUMA_calloc(nxyz, sizeof(short));
01245    if (!isin) {
01246       SUMA_SL_Crit("Failed to allocate");
01247       SUMA_RETURN(NULL);
01248    }
01249    
01250    /* mark the node voxels */
01251    for (nn=0; nn<SO->N_Node; ++nn) {
01252       /* find the ijk of each node: */
01253       nn3 = 3*nn;
01254       nijk = SUMA_3D_2_1D_index((int)NodeIJKlist[nn3], (int)NodeIJKlist[nn3+1] , (int)NodeIJKlist[nn3+2], nx , nxy); 
01255       if (nijk < nxyz) { if (!isin[nijk]) { isin[nijk] = SUMA_ON_NODE; ++(*N_inp); } }
01256    }
01257    
01258    
01259    /* cycle through all triangles and find voxels that intersect them */
01260    N_alloc = 2000; /* expected maximum number of voxels in triangle's bounding box */
01261    N_realloc = 0;
01262    voxelsijk = (int *)SUMA_malloc(sizeof(int)*N_alloc*3);
01263    if (!voxelsijk) { SUMA_SL_Crit("Failed to Allocate!"); SUMA_RETURN(NULL);  }   
01264    dxyz[0] = VolPar->dx; dxyz[1] = VolPar->dy; dxyz[2] = VolPar->dz;
01265    for (nf=0; nf<SO->N_FaceSet; ++nf) {
01266       n1 = SO->FaceSetList[SO->FaceSetDim*nf]; n2 = SO->FaceSetList[SO->FaceSetDim*nf+1]; n3 = SO->FaceSetList[SO->FaceSetDim*nf+2];
01267       /* find the bounding box of the triangle */
01268       p1 = &(NodeIJKlist[3*n1]); p2 = &(NodeIJKlist[3*n2]); p3 = &(NodeIJKlist[3*n3]); 
01269       SUMA_TRIANGLE_BOUNDING_BOX(p1, p2, p3, min_v, max_v);
01270       
01271       /* quick check of preallocate size of voxelsijk */
01272       en =((int)(max_v[0] - min_v[0] + 2) * (int)(max_v[1] - min_v[1] + 2) * (int)(max_v[2] - min_v[2] + 2)); 
01273       if ( en > N_alloc) {
01274          ++N_realloc; if (N_realloc > 5) { SUMA_SL_Warn("Reallocating, increase limit to improve speed.\nEither triangles too large or grid too small"); }
01275          N_alloc = 2*en;
01276          voxelsijk = (int *)SUMA_realloc(voxelsijk, 3*N_alloc*sizeof(int));
01277          if (!voxelsijk) { SUMA_SL_Crit("Failed to Allocate!"); SUMA_RETURN(NULL); }
01278       } 
01279       /* find the list of voxels inhabiting this box */
01280       N_inbox = 0;
01281       if (!SUMA_VoxelsInBox(voxelsijk, &N_inbox, min_v, max_v)) {
01282          SUMA_SL_Err("Unexpected error!"); SUMA_RETURN(NULL); 
01283       }
01284       if (!N_inbox) { SUMA_SL_Err("Unexpected error, no voxels in box!"); SUMA_RETURN(NULL);  }
01285       
01286       /* mark these voxels as inside the business */
01287       for (nt=0; nt < N_inbox; ++nt) {
01288          nt3 = 3*nt;
01289          if (voxelsijk[nt3] < nx &&  voxelsijk[nt3+1] < ny &&  voxelsijk[nt3+2] < nz) {
01290             nijk = SUMA_3D_2_1D_index(voxelsijk[nt3], voxelsijk[nt3+1], voxelsijk[nt3+2], nx , nxy);  
01291             if (!isin[nijk]) { 
01292                /* what side of the plane is this voxel on ? */
01293                p[0] = (float)voxelsijk[nt3]; p[1] = (float)voxelsijk[nt3+1]; p[2] = (float)voxelsijk[nt3+2]; 
01294                SUMA_DIST_FROM_PLANE(p1, p2, p3, p, dist);
01295                
01296                if (dist) {
01297                   if (SUMA_IS_NEG(VolPar->Hand * dist)) isin[nijk] = SUMA_IN_TRIBOX_INSIDE;  /* ZSS Added handedness factor. who would have thought? Be damned 3D coord systems! */
01298                   else isin[nijk] = SUMA_IN_TRIBOX_OUTSIDE; 
01299                }
01300                
01301                if (1) { /* does this triangle actually intersect this voxel ?*/
01302                   if (SUMA_isVoxelIntersect_Triangle (p, dxyz, p1, p2, p3)) {
01303                      if (isin[nijk] == SUMA_IN_TRIBOX_INSIDE) isin[nijk] = SUMA_INTERSECTS_TRIANGLE_INSIDE;
01304                      else isin[nijk] = SUMA_INTERSECTS_TRIANGLE_OUTSIDE;
01305                   } 
01306                }
01307                ++(*N_inp);    
01308             }
01309          }
01310       }
01311    }
01312    
01313    /* Now fill in inside of the sphere */
01314       /* create a mask vector*/
01315       ijkmask = (byte *)SUMA_calloc(nxyz, sizeof(byte));
01316       ijkout = (byte *)SUMA_calloc(nxyz, sizeof(byte));
01317       if (!ijkmask) {
01318          SUMA_SL_Crit("Failed to allocate");
01319          SUMA_RETURN(NULL);
01320       }
01321       for (nt=0; nt<nxyz; ++nt) {
01322          if (isin[nt]) {
01323             ijkmask[nt] = 1;
01324             if (isin[nt] == SUMA_IN_TRIBOX_OUTSIDE) ijkout[nt] = 1;
01325          }
01326       }
01327       
01328       /* seed: find the center the surface in the index coordinate system*/
01329       SUMA_MIN_MAX_SUM_VECMAT_COL (NodeIJKlist, SO->N_Node, 3, MinDims, MaxDims, SOCenter); 
01330       SOCenter[0] /= SO->N_Node;  SOCenter[1] /= SO->N_Node;   SOCenter[2] /= SO->N_Node;
01331       {
01332          float u[3], un, p0[3], p1[3];
01333          int Found = 0, cnt;
01334          SUMA_MT_INTERSECT_TRIANGLE *mti = NULL; 
01335 
01336          /* Ray from a node on the surface to the center */
01337          p0[0] = NodeIJKlist[0]; p1[0] = SOCenter[0]; 
01338          p0[1] = NodeIJKlist[1]; p1[1] = SOCenter[1]; 
01339          p0[2] = NodeIJKlist[2]; p1[2] = SOCenter[2]; 
01340          SUMA_UNIT_VEC(p0, p1, u, un);
01341          /* travel along that ray until you find a point inside the surface AND not on the mask */
01342          Found = 0; cnt = 1;
01343          while (!Found && cnt <= un) {
01344             p1[0] = p0[0] + cnt * u[0];
01345             p1[1] = p0[1] + cnt * u[1];
01346             p1[2] = p0[2] + cnt * u[2];
01347             if (LocalHead) {
01348                fprintf(SUMA_STDERR,"%s:\nTrying seed ijk is %d %d %d\n", FuncName, (int)p1[0], (int)p1[1], (int)p1[2]); 
01349             }
01350             ijkseed = SUMA_3D_2_1D_index(p1[0], p1[1], p1[2], nx , nxy);
01351             mti = SUMA_MT_intersect_triangle(p1, SOCenter, NodeIJKlist, SO->N_Node, SO->FaceSetList, SO->N_FaceSet, mti);
01352             if (!(mti->N_poshits % 2)) { /* number of positive direction hits is a multiple of 2 */
01353                /* seed is outside */
01354                SUMA_LH("Seed outside");
01355             } else {
01356                SUMA_LH("Seed inside");
01357                /* seed is inside, is it on the mask ? */
01358                if (!ijkmask[ijkseed]) { SUMA_LH("Seed Accepted");Found = YUP; }
01359                else SUMA_LH("Seed on mask");
01360             }
01361             ++cnt;   
01362          }
01363          if (!Found) {
01364             SUMA_SL_Err("Failed to find seed!");
01365             if (isin) SUMA_free(isin); isin = NULL;
01366             goto CLEAN_EXIT;
01367          }
01368          if (mti) mti = SUMA_Free_MT_intersect_triangle(mti);
01369       }
01370       
01371       /* hide the mask values that are outside the surface */
01372       for (nt=0; nt<nxyz; ++nt) { if (ijkout[nt]) isin[nt] = 0; }
01373       
01374       if (fillmaskset) {
01375                   fillmaskvec = (byte *)SUMA_malloc(nx*ny*nz*sizeof(byte));
01376                   if (!fillmaskvec) { SUMA_SL_Crit("Failed to allocate fillmaskvec"); }
01377                   /* is this value 0 in the fillmaskset */
01378                   EDIT_coerce_scale_type( nx*ny*nz , DSET_BRICK_FACTOR(fillmaskset,0) ,
01379                            DSET_BRICK_TYPE(fillmaskset,0), DSET_ARRAY(fillmaskset, 0) ,      /* input  */
01380                            MRI_byte               , fillmaskvec  ) ;   /* output */
01381       }
01382       
01383       /* get the mask */
01384       if (ijkseed < 0) {
01385          SUMA_SL_Err("Should not be here!");
01386          if (isin) SUMA_free(isin); isin = NULL;
01387          goto CLEAN_EXIT;
01388       }
01389       inmask = SUMA_FillToVoxelMask(ijkmask, ijkseed, nx, ny, nz, &N_in, NULL); 
01390       if (!inmask) {
01391          SUMA_SL_Err("Failed to FillToVoxelMask!");
01392       } else {
01393          if (fillhole) { /* not a good idea to keep in SUMA_FillToVoxelMask */
01394             byte *inmask_prefill=NULL;
01395             
01396             /* keep track of original inmask */
01397             inmask_prefill = (byte *)SUMA_calloc(nx * ny * nz , sizeof(byte));
01398             memcpy ((void*)inmask_prefill, (void *)inmask, nx * ny * nz * sizeof(byte));
01399 
01400             if (LocalHead) fprintf(SUMA_STDERR,"%s:\n filling small holes %d...\n", FuncName, fillhole);
01401             /* fill in any isolated holes in mask */
01402             (void) THD_mask_fillin_once( nx,ny,nz , inmask , fillhole ) ;  /* thd_automask.c */
01403             /* remove filled holes that overlap with 0 values in mask */
01404             if (fillmaskvec) {
01405                for (nt=0; nt<nxyz; ++nt) {
01406                   if (inmask_prefill[nt] == 0 && inmask[nt] == 1 && fillmaskvec[nt] == 0) { /* if a voxel was not in mask before filling, then was filled as a hole but it is zero in the original dset, then fill it not */
01407                      inmask[nt] = 0;
01408                      if (LocalHead) fprintf(SUMA_STDERR,"%s: Shutting fill at %d\n", FuncName, nt);
01409                   }
01410                }
01411             }
01412             SUMA_free(inmask_prefill); inmask_prefill = NULL;
01413          }
01414          /* flag it */
01415          for (nt=0; nt<nxyz; ++nt) {
01416             if (inmask[nt] && !isin[nt]) {
01417                /* it would be nice to fill this value only if the intensity in the spat normed volume is not 0 
01418                The problem is that you have to transform coordinates back from this dataset back to the spatnormed set*/
01419                   isin[nt] = SUMA_INSIDE_SURFACE;
01420             }
01421          }
01422       }
01423       /* now put back the values outside the surface but protect values in mask */
01424       for (nt=0; nt<nxyz; ++nt) { if (ijkout[nt] && !inmask[nt]) isin[nt] = SUMA_IN_TRIBOX_OUTSIDE; }
01425       
01426    CLEAN_EXIT:   
01427    if (ijkout) SUMA_free(ijkout); ijkout = NULL;
01428    if (ijkmask) SUMA_free(ijkmask); ijkmask = NULL;
01429    if (inmask) SUMA_free(inmask); inmask = NULL;
01430    if (voxelsijk) SUMA_free(voxelsijk); voxelsijk = NULL;
01431    SUMA_RETURN(isin);
01432    
01433 }               
01434 
01435 /*!
01436    \brief finds voxels inside a closed surface. Surface's normals are to point outwards.
01437    \sa   SUMA_SurfGridIntersect   
01438    \sa   SUMA_SURF_GRID_INTERSECT_OPTIONS for the various possible values
01439  
01440 */   
01441 short *SUMA_FindVoxelsInSurface (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole, THD_3dim_dataset *fillmaskset) 
01442 {
01443    static char FuncName[]={"SUMA_FindVoxelsInSurface"};
01444    short *isin = NULL;
01445    byte *tmpin = NULL;
01446    int i, N_in, j, k , n, khits, dims[2], N_hits, iii, jjj, niii, ncul;
01447    float *tmpXYZ=NULL, Center[3], MaxDims[3], MinDims[3], aMaxDims, aMinDims, delta_t;
01448    float hdim[3], t0, t1, t2, SOCenter[0], p0[3], p1[3];
01449    SUMA_MT_INTERSECT_TRIANGLE *mti = NULL; 
01450    struct  timeval tti;
01451    int nfound = 0, N_poshits = 0;
01452    int cnt = 0, sgn, sgnref, n1, n2, n3;
01453    float ivec[3] = { 1, 0, 0}, dp = 0.0, cx = 0.0;
01454    SUMA_Boolean LocalHead = NOPE;
01455    
01456 
01457    #if MSK_DBG
01458       FILE *fdb=fopen("SUMA_FindVoxelsInSurface.1D","w");
01459    #endif
01460    
01461    SUMA_ENTRY;
01462    
01463    SUMA_etime (&tti, 0);
01464    if (LocalHead) {
01465       fprintf(SUMA_STDERR,"%s: Patience...\n", FuncName);
01466    }
01467    N_in = *N_inp = 0;
01468    
01469    /* copy surface coordinates, we're going to ijk land */
01470    tmpXYZ = (float *)SUMA_malloc(SO->N_Node * 3 * sizeof(float));
01471    if (!tmpXYZ) {
01472       SUMA_SL_Crit("Faile to allocate");
01473       SUMA_RETURN(NULL);
01474    }
01475    
01476    memcpy ((void*)tmpXYZ, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float));
01477    
01478    /* transform the surface's coordinates from RAI to 3dfind */
01479    SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, VolPar);
01480    
01481    /* find the new center and bounding box for the surface in the new coordinate system*/
01482    SUMA_MIN_MAX_SUM_VECMAT_COL (tmpXYZ, SO->N_Node, 3, MinDims, MaxDims, SOCenter); 
01483    SOCenter[0] /= SO->N_Node;  
01484    SOCenter[1] /= SO->N_Node;  
01485    SOCenter[2] /= SO->N_Node;  
01486    SUMA_MIN_VEC (MinDims, 3, aMinDims );   
01487    SUMA_MAX_VEC (MaxDims, 3, aMaxDims);    
01488    /* Box half Dim */
01489    for (i=0; i<3; ++i) { hdim[i] = ( MaxDims[i] - MinDims[i]) / 2.0; }
01490    /* Box Center */
01491    for (i=0; i<3; ++i) { Center[i] = MinDims[i] + hdim[i]; }
01492    if (LocalHead) fprintf (SUMA_STDERR,"Center in IJK is: %f %f %f\n"
01493                                        "HlfDim in IJK is: %f %f %f\n"
01494                                        , Center[0], Center[1], Center[2],
01495                                          hdim[0], hdim[1], hdim[2]);
01496    
01497    /* Find the voxels that intersect the surface */
01498    isin = SUMA_SurfGridIntersect (SO, tmpXYZ, VolPar, &N_in, fillhole, fillmaskset);               
01499 
01500    *N_inp = N_in;
01501    #if MSK_DBG
01502       fclose(fdb); fdb = NULL;
01503    #endif
01504    
01505    delta_t = SUMA_etime (&tti, 1); 
01506    if (LocalHead) {
01507       fprintf(SUMA_STDERR,"%s: Execution time %f seconds.\n", FuncName, delta_t);
01508    }
01509    
01510    SUMA_free(tmpXYZ); tmpXYZ = NULL;
01511    if (mti) mti = SUMA_Free_MT_intersect_triangle(mti);
01512    if (tmpin) SUMA_free(tmpin); tmpin = NULL;
01513    
01514    SUMA_RETURN(isin);
01515 }
01516 
01517 /*!
01518    \brief A function to recommend node movement towards a better future.
01519    FAILS, mesh density too low. Need 3d edge mask....
01520    \sa SUMA_Reposition_Touchup
01521 */
01522 float *SUMA_Suggest_Touchup_Grad(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch)
01523 {
01524    static char FuncName[]={"SUMA_Suggest_Touchup_Grad"};
01525    int in, N_troub = 0, cond1=0, cond2=0, cond3 = 0, cond4 = 0, nn, ShishMax, *dvecind_under, *dvecind_over, ni, nj, nk, nx, ny, nxy;   
01526    float MinMax_over[2], MinMax[2], MinMax_dist[2], MinMax_over_dist[2], Means[3], tb; 
01527    float U[3], Un, *a, P2[2][3], *norm, shft, *b, gradient_promise; 
01528    float *touchup=NULL, targ[3]; 
01529    float *overshish=NULL, *undershish=NULL, *gradovershish=NULL;
01530    static FILE *GradOut = NULL;
01531    SUMA_Boolean LocalHead = YUP;
01532    
01533    SUMA_ENTRY;
01534       
01535    SUMA_SL_Warn("Don't call me, good for nothing, mesh density too low to work well ...\n");
01536    SUMA_RETURN(NULL);
01537       
01538       if (!GradOut) GradOut = fopen("GradOut.1D", "wa");
01539 
01540       if (Opt->debug > 2) LocalHead = YUP;
01541       *N_touch = 0;
01542       nx = SO->VolPar->nx; ny = SO->VolPar->ny; nxy = nx * ny;
01543       
01544       ShishMax = (int)(SUMA_MAX_PAIR(Opt->d1, Opt->d4) / Opt->travstp * 1.2);
01545       overshish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01546       undershish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01547       gradovershish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01548       dvecind_under = (int *)SUMA_malloc(sizeof(int)*ShishMax);
01549       dvecind_over = (int *)SUMA_malloc(sizeof(int)*ShishMax);
01550       touchup = (float *)SUMA_calloc(SO->N_Node, sizeof(float));   
01551       if (!touchup) { SUMA_SL_Crit("Failed to allocate"); SUMA_RETURN(NULL); }  
01552       for (in=0; in<SO->N_Node; ++in) {   
01553          SUMA_Find_IminImax(SO, Opt, in,  MinMax, MinMax_dist, MinMax_over, MinMax_over_dist, 
01554                            Means, undershish, overshish, dvecind_under, dvecind_over, ShishMax); 
01555          /* Shift the node outwards:   
01556             0- the minimum location of the minimum over the node is not 0
01557             AND
01558             1- if the minimum over the node is closer than the minimum below the node  
01559                or the internal minimum is more than 1mm away  and the internal min is masked to 0 
01560                or the minimum below the node just sucks, > tb and the minimum over the node is < tb
01561                The second part of the condition is meant to better model a finger of cortex that has enough csf on the inner surface 
01562                to have it be set to 0 by SpatNorm this dark csf would keep the brain surface from reaching the outer surface 
01563             AND   
01564             2- if the maximum over the node is much higher than the maximum below the node (fat on skull), 
01565                it must be beyond the location of the minimum over the node   
01566             AND   
01567             3- The new position of the node is not much larger than the current value
01568             AND
01569             4- There are no large negative gradients being crossed on the way out.  
01570          */ 
01571          tb = (MinMax[1] - Opt->t2) * 0.5 + Opt->t2; 
01572          cond1 = (    (MinMax_over_dist[0] < MinMax_dist[0]) 
01573                      || (MinMax_dist[0] > 1 && MinMax[0] >= MinMax_over[0] ) 
01574                      || (MinMax[0] > tb && MinMax_over[0] < MinMax[0]) ); 
01575          if (MinMax_over[1] > MinMax[1] && MinMax_over[1] > 0.9 * Opt->t98 && (MinMax_over_dist[1] < MinMax_over_dist[0]) ) cond2 = 0;  
01576          else cond2 = 1; 
01577          if (MinMax_over[0] > 1.2 * Means[0]) cond3 = 0;
01578          else cond3 = 1; 
01579          if (Opt->NodeDbg == in) {   
01580                a = &(SO->NodeList[3*in]);   
01581                fprintf(SUMA_STDERR, "%s: Debug during touchup for node %d:\n"   
01582                                     " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   
01583                                     " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   
01584                                     " Val at node %.1f, mean below %.1f, mean above %.1f\n"  
01585                                     " zalt= %f, r/2 = %f\n"    
01586                                     " Conditions: %f %d %d %d\n",    
01587                        FuncName, in, 
01588                        MinMax[0], MinMax[1], MinMax_dist[0], MinMax_dist[1], 
01589                        MinMax_over[0], MinMax_over[1], MinMax_over_dist[0], MinMax_over_dist[1], 
01590                        Means[0], Means[1], Means[2],   
01591                        a[2] - SO->Center[2], Opt->r/2, 
01592                        MinMax_over_dist[0] , cond1 , cond2, cond3); 
01593          }  
01594          /* Now search for the first location that has a negative gradient of 1/3 tm*/
01595          nn = 1;
01596          while(nn<ShishMax) {
01597             if (overshish[nn] >= 0) {
01598                gradovershish[nn-1] = overshish[nn] - overshish[nn-1];
01599                if (gradovershish[nn-1] < - Opt->tm/3.0) {
01600                   targ[0] = (nn - 1) * Opt->travstp * SO->NodeNormList[3*in]; 
01601                   targ[1] = (nn - 1) * Opt->travstp * SO->NodeNormList[3*in+1];
01602                   targ[2] = (nn - 1) * Opt->travstp * SO->NodeNormList[3*in+2];
01603                   if (dvecind_over[nn] < 0 || dvecind_over[nn] > Opt->nvox) {
01604                      fprintf(SUMA_STDERR,"Error %s: Bad value for  dvecind_over[%d]=%d\n", FuncName, nn, dvecind_over[in]);
01605                   } else {
01606                      if (MinMax_over_dist[0] && cond1 && cond2 && cond3 && Opt->dvec[dvecind_over[nn]]) { 
01607                         while (nn < ShishMax && dvecind_over[nn] >=0) {
01608                            Opt->dvec[dvecind_over[nn]] = 0;
01609                            if (LocalHead) { 
01610                               SUMA_1D_2_3D_index(dvecind_over[nn], ni, nj, nk, nx, nxy);
01611                               fprintf(SUMA_STDERR,"%s: Zeroing voxel %d %d %d, grad %f, from node %d\n", FuncName, ni, nj, nk, gradovershish[nn-1], in);
01612                               fprintf(GradOut,"%d %d %d\n",ni, nj, nk); 
01613                            }
01614                            ++nn;
01615                         }
01616                      }
01617                   }
01618                }
01619                ++nn;
01620             }else {
01621                gradovershish[nn-1] = 0;
01622                nn = ShishMax;
01623             }
01624          }
01625          /* STUFF BELOW IS USELESS */
01626          
01627          SUMA_NORM_VEC(targ, 3, gradient_promise);
01628          
01629          if (MinMax_over_dist[0] && cond1 && cond2 && cond3) {   
01630                a = &(SO->NodeList[3*in]);   
01631                /* reposition nodes */  
01632                if (0 && LocalHead) { fprintf(SUMA_STDERR, "%s: Suggest repair for node %d (better min above):\n"   
01633                                     " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   
01634                                     " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   
01635                                     " Val at node %.1f, mean below %.1f, mean above %.1f\n"  
01636                                     " zalt= %f, r/2 = %f\n",    
01637                        FuncName, in, 
01638                        MinMax[0], MinMax[1], MinMax_dist[0], MinMax_dist[1], 
01639                        MinMax_over[0], MinMax_over[1], MinMax_over_dist[0], MinMax_over_dist[1], 
01640                        Means[0], Means[1], Means[2],   
01641                        a[2] - SO->Center[2], Opt->r/2); } 
01642                {  /* Used to insist on value and top half if (MinMax_over_dist[0] && (a[2] - SO->Center[2] > 0))  */
01643                                              /* to avoid going into eye balls */
01644                                              /* Then added !SUMA_IS_EYE_ZONE(a,SO->Center) to avoid eyes */  
01645                                              /* but that is no longer necessary with fancier expansion conditions. */
01646                   if (gradient_promise < MinMax_over_dist[0]) {
01647                      touchup[in] = gradient_promise;
01648                   } else touchup[in] = MinMax_over_dist[0]; 
01649                }  
01650             ++N_troub;   
01651          }  
01652          /* nodes in the front of the brain that are sitting in the midst of a very bright area 
01653          should be moved further down until they hit the brain vs non brain area*/ 
01654       }  
01655       
01656       if (dvecind_under) SUMA_free(dvecind_under); dvecind_under= NULL;
01657       if (dvecind_over) SUMA_free(dvecind_over); dvecind_over= NULL;
01658       if (overshish) SUMA_free(overshish); overshish = NULL; 
01659       if (gradovershish) SUMA_free(gradovershish); gradovershish = NULL; 
01660       if (undershish) SUMA_free(undershish); undershish = NULL; 
01661    
01662    *N_touch = N_troub;
01663    if (GradOut) fclose(GradOut); GradOut = NULL;
01664    
01665    SUMA_RETURN(touchup);
01666 }
01667    
01668 /*!
01669    \brief A function to recommend node movement towards a better future.
01670    \sa SUMA_Reposition_Touchup
01671 */
01672 float *SUMA_Suggest_Touchup(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch)
01673 {
01674    static char FuncName[]={"SUMA_Suggest_Touchup"};
01675    int in, N_troub = 0, cond1=0, cond2=0, cond3 = 0, cond4 = 0, nn, ShishMax;   
01676    float MinMax_over[2], MinMax[2], MinMax_dist[2], MinMax_over_dist[2], Means[3], tb; 
01677    float U[3], Un, *a, P2[2][3], *norm, shft, *b; 
01678    float *touchup=NULL; 
01679    float *overshish=NULL, *undershish=NULL, *gradovershish=NULL;
01680    SUMA_Boolean LocalHead = NOPE;
01681    
01682    SUMA_ENTRY;
01683    
01684       if (Opt->debug > 2) LocalHead = YUP;
01685       *N_touch = 0;
01686       
01687       ShishMax = (int)(SUMA_MAX_PAIR(Opt->d1, Opt->d4) / Opt->travstp * 1.2);
01688       overshish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01689       undershish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01690       gradovershish = (float *)SUMA_malloc(sizeof(float)*ShishMax);
01691       touchup = (float *)SUMA_calloc(SO->N_Node, sizeof(float));   
01692       if (!touchup) { SUMA_SL_Crit("Failed to allocate"); SUMA_RETURN(NULL); }  
01693       for (in=0; in<SO->N_Node; ++in) {   
01694          SUMA_Find_IminImax(SO, Opt, in,  MinMax, MinMax_dist, MinMax_over, MinMax_over_dist, Means, undershish, overshish, NULL, NULL, ShishMax); 
01695          /* Shift the node outwards:   
01696             0- the minimum location of the minimum over the node is not 0
01697             AND
01698             1- if the minimum over the node is closer than the minimum below the node  or the internal minimum is more than 1mm away 
01699                and the internal min is masked to 0 or the minimum below the node just sucks, > tb and the minimum over the node is < tb
01700                The second part of the condition is meant to better model a finger of cortex that has enough csf on the inner surface 
01701                to have it be set to 0 by SpatNorm this dark csf would keep the brain surface from reaching the outer surface 
01702             AND   
01703             2- if the maximum over the node is much higher than the maximum below the node (fat on skull), 
01704                it must be beyond the location of the minimum over the node   
01705             AND   
01706             3- The new position of the node is not much larger than the current value
01707             AND
01708             4- The mean above the node is not more than 1.2 * the mean below the node.  
01709          */ 
01710          tb = (MinMax[1] - Opt->t2) * 0.5 + Opt->t2; 
01711          cond1 = (    (MinMax_over_dist[0] < MinMax_dist[0]) || (MinMax_dist[0] > 1 && MinMax[0] >= MinMax_over[0] ) 
01712                      || (MinMax[0] > tb && MinMax_over[0] < MinMax[0]) ); 
01713          /* statement below used to be: if (MinMax_over[1] > MinMax[1] && MinMax_over[1] > 0.9 * Opt->t98 && (MinMax_over_dist[1] < MinMax_over_dist[0]) ) cond2 = 0;  */
01714          if (MinMax_over[1] > 1.2 * MinMax[1] && (MinMax_over_dist[1] < MinMax_over_dist[0]) ) cond2 = 0;  
01715          else cond2 = 1; 
01716          if (MinMax_over[0] > 1.2 * Means[0]) cond3 = 0;
01717          else cond3 = 1;
01718          if (Means[2] > Means[1]) cond4 = 0;  /* March 21 */
01719          else cond4 = 1; 
01720          if (Opt->NodeDbg == in) {   
01721                a = &(SO->NodeList[3*in]);   
01722                fprintf(SUMA_STDERR, "%s: Debug during touchup for node %d:\n"   
01723                                     " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   
01724                                     " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   
01725                                     " Val at node %.1f, mean below %.1f, mean above %.1f\n"  
01726                                     " zalt= %f, r/2 = %f\n"    
01727                                     " Conditions: %f %d %d %d %d\n",    
01728                        FuncName, in, 
01729                        MinMax[0], MinMax[1], MinMax_dist[0], MinMax_dist[1], 
01730                        MinMax_over[0], MinMax_over[1], MinMax_over_dist[0], MinMax_over_dist[1], 
01731                        Means[0], Means[1], Means[2],   
01732                        a[2] - SO->Center[2], Opt->r/2, 
01733                        MinMax_over_dist[0] , cond1 , cond2, cond3, cond4); 
01734          }  
01735          /* what is happening with the gradients NOT USED YET*/
01736          for (nn=1; nn<ShishMax; ++nn) {
01737             if (overshish[nn] >= 0) {
01738                gradovershish[nn-1] = overshish[nn] - overshish[nn-1];
01739                   
01740             }else {
01741                gradovershish[nn-1] = 0;
01742             }
01743          }
01744          
01745          if (MinMax_over_dist[0] && cond1 && cond2 && cond3) {   
01746                a = &(SO->NodeList[3*in]);   
01747                if (cond4) {
01748                   /* reposition nodes */  
01749                   if (Opt->NodeDbg == in) { fprintf(SUMA_STDERR, "%s: Suggest repair for node %d (better min above):\n"   
01750                                        " MinMax     =[%.1f,%.1f] MinMax_dist     = [%.2f,%.2f] \n"   
01751                                        " MinMax_over=[%.1f,%.1f] MinMax_over_dist= [%.2f,%.2f] \n"   
01752                                        " Val at node %.1f, mean below %.1f, mean above %.1f\n"  
01753                                        " zalt= %f, r/2 = %f\n",    
01754                           FuncName, in, 
01755                           MinMax[0], MinMax[1], MinMax_dist[0], MinMax_dist[1], 
01756                           MinMax_over[0], MinMax_over[1], MinMax_over_dist[0], MinMax_over_dist[1], 
01757                           Means[0], Means[1], Means[2],   
01758                           a[2] - SO->Center[2], Opt->r/2); } 
01759                   {  /* Used to insist on value and top half if (MinMax_over_dist[0] && (a[2] - SO->Center[2] > 0))  */
01760                                                 /* to avoid going into eye balls */
01761                                                 /* Then added !SUMA_IS_EYE_ZONE(a,SO->Center) to avoid eyes */  
01762                                                 /* but that is no longer necessary with fancier expansion conditions. */
01763                      touchup[in] = MinMax_over_dist[0]; 
01764                   }  
01765                ++N_troub; 
01766             } else { /* freeze that node, an attempt to control leak into lardo */
01767                Opt->Stop[in] = -1.0;
01768                if (Opt->NodeDbg == in) { fprintf(SUMA_STDERR, "%s: Freezing node %d\n", FuncName, in); }
01769             }  
01770          }  
01771          /* nodes in the front of the brain that are sitting in the midst of a very bright area 
01772          should be moved further down until they hit the brain vs non brain area*/ 
01773       }  
01774       
01775       
01776       if (overshish) SUMA_free(overshish); overshish = NULL; 
01777       if (gradovershish) SUMA_free(gradovershish); gradovershish = NULL; 
01778       if (undershish) SUMA_free(undershish); undershish = NULL; 
01779    
01780    *N_touch = N_troub;
01781    
01782    SUMA_RETURN(touchup);
01783 }
01784 /*!
01785    \brief A function to move node surfaces towards a better future.
01786    
01787    This function used to be a macro gone too large. 
01788 */
01789 int SUMA_Reposition_Touchup(SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs) 
01790 {
01791    static char FuncName[]={"SUMA_Reposition_Touchup"};
01792    byte *fmask=NULL;
01793    int in, N_troub = 0, cond1=0, cond2=0, cond3 = 0, cond4 = 0, nn;   
01794    float MinMax_over[2], MinMax[2], MinMax_dist[2], MinMax_over_dist[2], Means[3], tb; 
01795    float U[3], Un, *a, P2[2][3], *norm, shft, *b; 
01796    float *touchup=NULL, **wgt=NULL, *dsmooth=NULL; 
01797    int nstp, stillmoving, kth_buf, N_Touch;
01798    float stp, *csmooth=NULL, *shftvec=NULL, *targetloc=NULL;
01799    SUMA_Boolean Send_buf;
01800    SUMA_Boolean LocalHead = NOPE;
01801    
01802    SUMA_ENTRY;
01803    
01804       if (Opt->debug > 2) LocalHead = YUP;
01805       
01806       touchup = SUMA_Suggest_Touchup(SO, Opt, limtouch, cs, &N_troub);
01807       if (!touchup) {
01808          SUMA_SL_Err("Failed in SUMA_Suggest_Touchup");
01809          SUMA_RETURN(NOPE);
01810       }
01811       if (!N_troub) {
01812          SUMA_LH("Nothing to do, no trouble nodes.");
01813          SUMA_RETURN(YUP);
01814       }
01815       
01816       if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* %d troubled nodes found\n", FuncName, N_troub); 
01817       if (0){ /* March 22: Smoothing is not such a hot idea because it might cause certain nodes to depart into void and 
01818          then beyond, if you do a 2 pass touchup. Do not do it anymore */
01819           
01820          /* fix the big bumps: Bad for brain surface that has lots of bump on top, like 201, 
01821          that kind of smoothing is better performed on an inner model of the skull .... Someday ....*/
01822          if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* Now filtering...\n", FuncName); 
01823          wgt = SUMA_Chung_Smooth_Weights(SO);   
01824          if (!wgt) { 
01825             SUMA_SL_Err("Failed to compute weights.\n"); 
01826             exit(1); 
01827          }  
01828          dsmooth = SUMA_Chung_Smooth (SO, wgt, 200, 20, touchup, 1, SUMA_COLUMN_MAJOR, NULL, cs);   
01829          if (LocalHead) fprintf (SUMA_STDERR,"%s: ********************* filtering done.\n", FuncName);   
01830       } else { 
01831          dsmooth = touchup; /* carefull at the free business */
01832       }  
01833       /* add the changes */   
01834       #if 1 /* sudden jump, Fast and furious but that causes intersections */ 
01835       for (in=0; in<SO->N_Node; ++in) {   
01836          a = &(SO->NodeList[3*in]);   
01837          if (Opt->Stop[in] >= 0) {  
01838             if (1 || !SUMA_IS_EYE_ZONE(a,SO->Center)) { 
01839                if (a[2] - SO->Center[2] > 10 )  shft = touchup[in]; /* no smooth baby for top where bumpy sulci can occur */ 
01840                else shft = dsmooth[in];   
01841                if (shft) { 
01842                   a = &(SO->NodeList[3*in]);   
01843                   norm = &(SO->NodeNormList[3*in]);  
01844                   SUMA_POINT_AT_DISTANCE(norm, a, SUMA_MIN_PAIR(shft, limtouch), P2);   
01845                   SO->NodeList[3*in] = P2[0][0]; SO->NodeList[3*in+1] = P2[0][1]; SO->NodeList[3*in+2] = P2[0][2];   
01846                   if (0 && LocalHead) fprintf(SUMA_STDERR,"%s: Acting on node %d because it is in top half, boosting by %f\n", 
01847                            FuncName, in, SUMA_MIN_PAIR(shft, limtouch));   
01848                }
01849             }
01850          } else {
01851             if (in == Opt->NodeDbg) fprintf(SUMA_STDERR,"%s:\nNode %d has been frozen before, no cigar.\n", FuncName, in);
01852          }
01853       }  
01854       #else /* smooth operator, slow as hell, does not reach far enough */
01855       /* first figure out where each node is to go: */
01856       shftvec = (float *)SUMA_calloc( SO->N_Node, sizeof(float));
01857       targetloc = (float *)SUMA_malloc(3 * SO->N_Node * sizeof(float));
01858       for (in=0; in<SO->N_Node; ++in) {  
01859          shftvec[in] = 0.0;
01860          a = &(SO->NodeList[3*in]);   
01861          if (1 || !SUMA_IS_EYE_ZONE(a,SO->Center)) { 
01862             if (a[2] - SO->Center[2] > 10 )  shft = touchup[in]; /* no smooth baby for top where bumpy sulci can occur */ 
01863             else shft = dsmooth[in];   
01864             shft = SUMA_MIN_PAIR(shft, limtouch);
01865             if (shft) { 
01866                shftvec[in] = shft;
01867                norm = &(SO->NodeNormList[3*in]); 
01868                targetloc[3*in  ] = shft*norm[0]+a[0];
01869                targetloc[3*in+1] = shft*norm[1]+a[1];
01870                targetloc[3*in+2] = shft*norm[2]+a[2]; 
01871             }
01872          }
01873       }
01874       /* now move slowly towards target, only move those nodes that are to be shifted */
01875       fmask = (byte *)SUMA_calloc(SO->N_Node , sizeof(byte));
01876       stp = 1; /* incremental step in mm */
01877       nstp = 0;
01878       do {
01879          stillmoving = 0;
01880          for (in=0; in<SO->N_Node; ++in) {  
01881             fmask[in] = 0;
01882             a = &(SO->NodeList[3*in]); 
01883             b = &(targetloc[3*in]);
01884             if (shftvec[in]) {
01885                SUMA_UNIT_VEC( a, b, U, Un);
01886                if (Un < stp)  {
01887                   shft = Un;
01888                   shftvec[in] = 0.0; /* no more shifts */
01889                } else {
01890                   shft = stp;
01891                   stillmoving = 1;
01892                }
01893                /* fprintf(SUMA_STDERR,"%s:\n node %d shft=%f\nU=[%f %f %f]\n", FuncName, in, shft, U[0], U[1], U[2]); */
01894                if (shft) { 
01895                   SUMA_POINT_AT_DISTANCE(U, a, shft, P2);   
01896                   SO->NodeList[3*in] = P2[0][0]; SO->NodeList[3*in+1] = P2[0][1]; SO->NodeList[3*in+2] = P2[0][2];
01897                   fmask[in] = 1; /* filter the one that moves */
01898                }         
01899             }
01900          }
01901             Send_buf = cs->Send;
01902             cs->Send = NOPE;
01903             kth_buf = cs->kth;
01904             cs->kth = 1;  /*make sure all gets sent at this stage */
01905          if (1) {/* filter the surface */
01906             fprintf (SUMA_STDERR,"%s: Touchup smoothing.\n", FuncName);
01907             csmooth = SUMA_Taubin_Smooth( SO, NULL,
01908                                           0.6307, -.6732, SO->NodeList,
01909                                           20, 3, SUMA_ROW_MAJOR, csmooth, cs, fmask);    
01910             memcpy((void*)SO->NodeList, (void *)csmooth, SO->N_Node * 3 * sizeof(float));
01911          }      
01912             cs->Send = Send_buf;
01913             SUMA_RECOMPUTE_NORMALS(SO);
01914             if (cs->Send) {
01915                if (!SUMA_SendToSuma (SO, cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
01916                   SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
01917                }
01918             }
01919             cs->kth = kth_buf;
01920          if (LocalHead) fprintf (SUMA_STDERR,"%s: Touchup pass %d complete.\n", FuncName, nstp);  
01921          ++nstp;
01922       } while (stillmoving && nstp < 60);
01923       if (LocalHead) fprintf (SUMA_STDERR,"%s: Stopped growth stillmoving = %d, nstp = %d\n", FuncName, stillmoving, nstp);
01924       #endif 
01925       if (fmask) SUMA_free(fmask); fmask = NULL; 
01926       if (shftvec) SUMA_free(shftvec); shftvec = NULL;
01927       if (targetloc) SUMA_free(targetloc); targetloc = NULL;
01928       if (touchup == dsmooth) dsmooth = NULL;
01929       if (touchup) SUMA_free(touchup); touchup = NULL;   
01930       if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL;   
01931       if (csmooth) SUMA_free(csmooth); csmooth = NULL;   
01932       if (wgt) SUMA_free2D ((char **)wgt, SO->N_Node); wgt = NULL;   
01933    
01934    SUMA_RETURN(YUP);
01935 }
01936 
01937 /*!
01938    \brief Creates an empty AFNI EDIT_options structure.
01939     All fields are initialized to do nothing.
01940     The fields initialized are based on those used in EDIT_one_dataset in edt_onedset.c
01941    Free returned pointer with SUMA_free
01942 */
01943 EDIT_options *SUMA_BlankAfniEditOptions(void)
01944 {
01945    static char FuncName[]={"SUMA_BlankAfniEditOptions"};
01946    EDIT_options *edopt = NULL;
01947 
01948    SUMA_ENTRY;
01949 
01950    edopt = (EDIT_options *)SUMA_calloc(1, sizeof(EDIT_options));
01951    /* set all fields to do nothing */
01952    edopt->thtoin = 0; 
01953    edopt->noneg = 0; 
01954    edopt->abss = 0; 
01955    edopt->clip_bot = 0; 
01956    edopt->clip_top = 0; 
01957    edopt->thresh = 0.0; 
01958    edopt->edit_clust = ECFLAG_SAME - 1;
01959    edopt->clust_rmm = 0.0;
01960         edopt->clust_vmul = 0.0;
01961    edopt->erode_pv  = 0.0;
01962    edopt->dilate = 0;
01963    edopt->filter_opt = FCFLAG_NONE;
01964    edopt->filter_rmm = 0.0;
01965    edopt->thrfilter_opt = FCFLAG_NONE;
01966    edopt->thrfilter_rmm = 0.0;
01967    edopt->blur = 0.0; 
01968    edopt->thrblur = 0.0; 
01969    edopt->scale = 0.0; 
01970    edopt->mult = 0.0; 
01971    edopt->do_zvol = 0; 
01972    edopt->iv_fim = -1;
01973    edopt->iv_thr = -1;
01974    edopt->verbose = 0;
01975    edopt->fake_dxyz = 0;
01976    edopt->clip_unscaled = 0; 
01977 
01978    SUMA_RETURN(edopt);
01979 }