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
00012 SUMA_SurfaceViewer *SUMAg_cSV = NULL;
00013 SUMA_SurfaceViewer *SUMAg_SVv = NULL;
00014
00015 int SUMAg_N_SVv = 0;
00016 SUMA_DO *SUMAg_DOv = NULL;
00017 int SUMAg_N_DOv = 0;
00018 SUMA_CommonFields *SUMAg_CF = NULL;
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
00036
00037
00038
00039
00040
00041
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) ,
00075 MRI_double , Opt->dvec ) ;
00076
00077
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
00098 }
00099
00100
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
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
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
00142 if (Opt->tm < 0) {
00143
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
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
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
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) { 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
00235
00236 if (0) {
00237 SUMA_SL_Err("Failed to subdivide mesh");
00238 SUMA_Free_Surface_Object (SOhullp);
00239 *SOhull = NULL;
00240 } else {
00241
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);
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
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
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
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
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];
00331
00332
00333
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
00339
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; }
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
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
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
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) {
00383 Means[0] = Opt->dvec[nind];
00384 }
00385
00386 if (overshish && istep < ShishMax) overshish[istep] = Opt->dvec[nind];
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
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; }
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
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
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];
00442
00443 SUMA_ENTRY;
00444
00445 d1 = Opt->d1; d2 = d1/2.0;
00446
00447 vtn[N_vtnmax-1] = -1;
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
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
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
00472
00473
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];
00487
00488
00489
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
00495
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; }
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
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
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
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) {
00539 Means[0] = Opt->dvec[nind];
00540 }
00541
00542 if (overshish && istep < ShishMax) overshish[istep] = Opt->dvec[nind];
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
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; }
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
00569
00570
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
00612 sksearch = 15;
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
00622 if (n[2] > -45) {
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 }
00635
00636
00637
00638 tmpptr = SO->NodeList;
00639 SO->NodeList = NewCoord;
00640 NewCoord = tmpptr;
00641 tmpptr = NULL;
00642
00643 cs->kth = 1;
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
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 }
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
00673
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;
00765 if (Opt->debug > 1) fprintf(SUMA_STDERR,"%s: Pass #%d\n", FuncName, Npass);
00766 MaxExp = 0.0;
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));
00773 else lztfac = 1.2;
00774 } else lztfac = 1.0;
00775 MaxExp = 0.0;
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);
00781 SUMA_SCALE_VEC(&(SO->NodeNormList[3*in]), sn, dp, 3, float, float);
00782 SUMA_NORM_VEC(sn, 3, nsn);
00783 SUMA_SUB_VEC(s, sn, st, 3, float, float, float);
00784 SUMA_SCALE_VEC(st, u1, su1, 3, float, float);
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);
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)); }
00794
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;
00808 } else su4 = 0;
00809
00810
00811 if (Opt->NoEyes) {
00812 if (it > 0.1 * Opt->N_it) {
00813 if (SUMA_IS_EYE_ZONE(n,SO->Center)) {
00814
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
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 }
00901
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) {
00909 SUMA_WRAP_BRAIN_SMOOTH_NN(Opt->NNsmooth, dsmooth, refNodeList);
00910 }
00911 }
00912
00913
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 }
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) {
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
00949 if (!pastarea) {
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
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
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
01003 Done = 1;
01004 } else {
01005
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;
01030 } else {
01031 Done = 1;
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
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
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
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
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
01111 SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, VolPar);
01112
01113
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
01121 for (i=0; i<3; ++i) { hdim[i] = ( MaxDims[i] - MinDims[i]) / 2.0; }
01122
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;
01134
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]) {
01146 p0[0] = i; p1[0] = i+1000; p0[1] = p1[1] = j; p0[2] = p1[2] = k;
01147 #ifdef Meth1
01148
01149 mti = SUMA_MT_intersect_triangle(p0, p1, tmpXYZ, SO->N_Node, SO->FaceSetList, SO->N_FaceSet, mti);
01150 if (!(mti->N_poshits % 2)) {
01151 isin[n] = 0; --N_in;
01152 }
01153 #else
01154 dims[0] = 1; dims[1] = 2;
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
01162
01163
01164
01165
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
01173 if (sgnref == sgn) ++N_poshits;
01174 }
01175 ++cnt;
01176 }
01177 if (N_poshits % 2 == 0) {
01178 isin[n] = 0; --N_in;
01179
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
01211
01212
01213
01214
01215
01216
01217
01218
01219
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
01251 for (nn=0; nn<SO->N_Node; ++nn) {
01252
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
01260 N_alloc = 2000;
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
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
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
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
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
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;
01298 else isin[nijk] = SUMA_IN_TRIBOX_OUTSIDE;
01299 }
01300
01301 if (1) {
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
01314
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
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
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
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)) {
01353
01354 SUMA_LH("Seed outside");
01355 } else {
01356 SUMA_LH("Seed inside");
01357
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
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
01378 EDIT_coerce_scale_type( nx*ny*nz , DSET_BRICK_FACTOR(fillmaskset,0) ,
01379 DSET_BRICK_TYPE(fillmaskset,0), DSET_ARRAY(fillmaskset, 0) ,
01380 MRI_byte , fillmaskvec ) ;
01381 }
01382
01383
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) {
01394 byte *inmask_prefill=NULL;
01395
01396
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
01402 (void) THD_mask_fillin_once( nx,ny,nz , inmask , fillhole ) ;
01403
01404 if (fillmaskvec) {
01405 for (nt=0; nt<nxyz; ++nt) {
01406 if (inmask_prefill[nt] == 0 && inmask[nt] == 1 && fillmaskvec[nt] == 0) {
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
01415 for (nt=0; nt<nxyz; ++nt) {
01416 if (inmask[nt] && !isin[nt]) {
01417
01418
01419 isin[nt] = SUMA_INSIDE_SURFACE;
01420 }
01421 }
01422 }
01423
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
01437
01438
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
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
01479 SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, VolPar);
01480
01481
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
01489 for (i=0; i<3; ++i) { hdim[i] = ( MaxDims[i] - MinDims[i]) / 2.0; }
01490
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
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
01519
01520
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
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
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
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
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
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 {
01643
01644
01645
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
01653
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
01670
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
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
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
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;
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
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
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 {
01760
01761
01762
01763 touchup[in] = MinMax_over_dist[0];
01764 }
01765 ++N_troub;
01766 } else {
01767 Opt->Stop[in] = -1.0;
01768 if (Opt->NodeDbg == in) { fprintf(SUMA_STDERR, "%s: Freezing node %d\n", FuncName, in); }
01769 }
01770 }
01771
01772
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
01786
01787
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){
01818
01819
01820
01821
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;
01832 }
01833
01834 #if 1
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];
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
01855
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];
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
01875 fmask = (byte *)SUMA_calloc(SO->N_Node , sizeof(byte));
01876 stp = 1;
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;
01889 } else {
01890 shft = stp;
01891 stillmoving = 1;
01892 }
01893
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;
01898 }
01899 }
01900 }
01901 Send_buf = cs->Send;
01902 cs->Send = NOPE;
01903 kth_buf = cs->kth;
01904 cs->kth = 1;
01905 if (1) {
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
01939
01940
01941
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
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 }