Doxygen Source Code Documentation
SUMA_BrainWrap.c File Reference
#include "SUMA_suma.h"#include "../thd_brainormalize.h"Go to the source code of this file.
Defines | |
| #define | MSK_DBG 0 |
| #define | Meth1 |
Functions | |
| int | SUMA_DidUserQuit (void) |
| float | SUMA_LoadPrepInVol (SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_SurfaceObject **SOhull) |
| int | SUMA_Find_IminImax (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, int ni, float *MinMax, float *MinMax_dist, float *MinMax_over, float *MinMax_over_dist, float *Means, float *undershish, float *overshish, int *dvecind_under, int *dvecind_over, int ShishMax) |
| int | SUMA_Find_IminImax_Avg (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, int ni, float *MinMax, float *MinMax_dist, float *MinMax_over, float *MinMax_over_dist, float *Means, float *undershish, float *overshish, int *dvecind_under, int *dvecind_over, int ShishMax) |
| int | SUMA_SkullMask (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs) |
| creates a crude mask of the brain's limit based on the surface of the skull mask is crude, for purposes of keeping BrainWrap from leaking out when there are very strong shading artifacts | |
| int | SUMA_StretchToFitLeCerveau (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, SUMA_COMM_STRUCT *cs) |
| byte * | SUMA_FindVoxelsInSurface_SLOW (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole) |
| First pass at finding voxels inside a surface. Slow baby, slow. The first approach traces a ray parallel to the X axis from each voxel and counts the number of intersections, in the direction of the ray, with the surface. Odd = in, Even = out. That is Meth1. The other method, projects all triangles to the YZ plane and then finds out if a voxel's projection lies inside the projected triangles. You still have to do a test to count only the triangles that are intersected in the direction of the ray In the end, it turns out to be slower than Meth1. Bad. | |
| short * | SUMA_SurfGridIntersect (SUMA_SurfaceObject *SO, float *NodeIJKlist, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole, THD_3dim_dataset *fillmaskset) |
| short * | SUMA_FindVoxelsInSurface (SUMA_SurfaceObject *SO, SUMA_VOLPAR *VolPar, int *N_inp, int fillhole, THD_3dim_dataset *fillmaskset) |
| finds voxels inside a closed surface. Surface's normals are to point outwards. | |
| float * | SUMA_Suggest_Touchup_Grad (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch) |
| A function to recommend node movement towards a better future. FAILS, mesh density too low. Need 3d edge mask.... | |
| float * | SUMA_Suggest_Touchup (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs, int *N_touch) |
| A function to recommend node movement towards a better future. | |
| int | SUMA_Reposition_Touchup (SUMA_SurfaceObject *SO, SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt, float limtouch, SUMA_COMM_STRUCT *cs) |
| A function to move node surfaces towards a better future. | |
| EDIT_options * | SUMA_BlankAfniEditOptions (void) |
| Creates an empty AFNI EDIT_options structure. All fields are initialized to do nothing. The fields initialized are based on those used in EDIT_one_dataset in edt_onedset.c Free returned pointer with SUMA_free. | |
Variables | |
| SUMA_CommonFields * | SUMAg_CF |
| SUMA_DO * | SUMAg_DOv |
| SUMA_SurfaceViewer * | SUMAg_SVv |
| int | SUMAg_N_SVv |
| int | SUMAg_N_DOv |
| int | InteractiveQuit |
Define Documentation
|
|
Definition at line 1062 of file SUMA_BrainWrap.c. |
|
|
Definition at line 1061 of file SUMA_BrainWrap.c. |
Function Documentation
|
|
Creates an empty AFNI EDIT_options structure. All fields are initialized to do nothing. The fields initialized are based on those used in EDIT_one_dataset in edt_onedset.c Free returned pointer with SUMA_free.
Definition at line 1943 of file SUMA_BrainWrap.c. References EDIT_options::abss, EDIT_options::blur, EDIT_options::clip_bot, EDIT_options::clip_top, EDIT_options::clip_unscaled, EDIT_options::clust_rmm, EDIT_options::clust_vmul, EDIT_options::dilate, EDIT_options::do_zvol, ECFLAG_SAME, EDIT_options::edit_clust, EDIT_options::erode_pv, EDIT_options::fake_dxyz, FCFLAG_NONE, EDIT_options::filter_opt, EDIT_options::filter_rmm, EDIT_options::iv_fim, EDIT_options::iv_thr, EDIT_options::mult, EDIT_options::noneg, EDIT_options::scale, SUMA_calloc, SUMA_ENTRY, SUMA_RETURN, EDIT_options::thrblur, EDIT_options::thresh, EDIT_options::thrfilter_opt, EDIT_options::thrfilter_rmm, EDIT_options::thtoin, and EDIT_options::verbose.
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 }
|
|
|
Definition at line 29 of file SUMA_BrainWrap.c.
00030 {
00031 return(InteractiveQuit);
00032 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 285 of file SUMA_BrainWrap.c. References SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, DSET_NX, DSET_NY, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, THD_ivec3::ijk, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, lmax(), lmin(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_ENTRY, SUMA_MAX_PAIR, SUMA_MIN_PAIR, SUMA_NORM_VEC, SUMA_RETURN, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, THD_3dmm_to_3dind(), THD_dicomm_to_3dmm(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::tm, SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp, and THD_fvec3::xyz. Referenced by SUMA_StretchToFitLeCerveau(), SUMA_Suggest_Touchup(), and SUMA_Suggest_Touchup_Grad().
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 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 427 of file SUMA_BrainWrap.c. References SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, DSET_NX, DSET_NY, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, THD_ivec3::ijk, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, lmax(), lmin(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_ENTRY, SUMA_Get_NodeIncident(), SUMA_MAX_PAIR, SUMA_MIN_PAIR, SUMA_NORM_VEC, SUMA_RETURN, SUMA_SL_Err, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, THD_3dmm_to_3dind(), THD_dicomm_to_3dmm(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::tm, SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp, and THD_fvec3::xyz.
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 }
|
|
||||||||||||||||||||||||
|
finds voxels inside a closed surface. Surface's normals are to point outwards.
Definition at line 1441 of file SUMA_BrainWrap.c. References i, LocalHead, n1, n2, SUMA_SurfaceObject::N_Node, SUMA_SurfaceObject::NodeList, SUMA_Boolean, SUMA_ENTRY, SUMA_etime(), SUMA_free, SUMA_Free_MT_intersect_triangle(), SUMA_malloc, SUMA_MAX_VEC, SUMA_MIN_MAX_SUM_VECMAT_COL, SUMA_MIN_VEC, SUMA_RETURN, SUMA_SL_Crit, SUMA_SurfGridIntersect(), and SUMA_vec_dicomm_to_3dfind(). Referenced by main().
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 }
|
|
||||||||||||||||||||
|
First pass at finding voxels inside a surface. Slow baby, slow. The first approach traces a ray parallel to the X axis from each voxel and counts the number of intersections, in the direction of the ray, with the surface. Odd = in, Even = out. That is Meth1. The other method, projects all triangles to the YZ plane and then finds out if a voxel's projection lies inside the projected triangles. You still have to do a test to count only the triangles that are intersected in the direction of the ray In the end, it turns out to be slower than Meth1. Bad.
Definition at line 1075 of file SUMA_BrainWrap.c. References SUMA_SurfaceObject::FaceSetDim, SUMA_SurfaceObject::FaceSetList, i, LocalHead, n1, n2, SUMA_SurfaceObject::N_FaceSet, SUMA_SurfaceObject::N_Node, SUMA_MT_INTERSECT_TRIANGLE::N_poshits, SUMA_SurfaceObject::NodeList, SUMA_VOLPAR::nx, SUMA_VOLPAR::ny, SUMA_VOLPAR::nz, SUMA_ABS, SUMA_Boolean, SUMA_ENTRY, SUMA_etime(), SUMA_free, SUMA_Free_MT_intersect_triangle(), SUMA_isinpoly(), SUMA_malloc, SUMA_MAX_VEC, SUMA_MIN_MAX_SUM_VECMAT_COL, SUMA_MIN_VEC, SUMA_MT_intersect_triangle(), SUMA_RETURN, SUMA_SIGN, SUMA_SL_Crit, and SUMA_vec_dicomm_to_3dfind().
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 }
|
|
||||||||||||
|
Definition at line 43 of file SUMA_BrainWrap.c. References SUMA_GENERIC_PROG_OPTIONS_STRUCT::cog, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_DX, DSET_DY, DSET_DZ, DSET_NVALS, DSET_NVOX, DSET_NX, DSET_NY, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, EDIT_coerce_scale_type(), SUMA_SurfaceObject::EL, SUMA_SurfaceObject::FaceNormList, SUMA_SurfaceObject::FaceSetList, free, SUMA_SurfaceObject::glar_FaceNormList, SUMA_SurfaceObject::glar_FaceSetList, SUMA_SurfaceObject::glar_NodeList, SUMA_SurfaceObject::glar_NodeNormList, i, THD_ivec3::ijk, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, SUMA_ISINSPHERE::IsIn, SUMA_GENERIC_PROG_OPTIONS_STRUCT::k98mask, SUMA_GENERIC_PROG_OPTIONS_STRUCT::k98maskcnt, SUMA_GENERIC_PROG_OPTIONS_STRUCT::Kill98, LocalHead, SUMA_SurfaceObject::N_FaceSet, SUMA_SurfaceObject::N_Node, SUMA_ISINSPHERE::nIsIn, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_GENERIC_PROG_OPTIONS_STRUCT::nvox, SUMA_GENERIC_PROG_OPTIONS_STRUCT::r, SUMA_Boolean, SUMA_calloc, SUMA_dPercRange(), SUMA_ENTRY, SUMA_free, SUMA_free_Edge_List(), SUMA_Free_IsInSphere(), SUMA_Free_Surface_Object(), SUMA_isinsphere(), SUMA_isSkin(), SUMA_Make_Edge_List_eng(), SUMA_MakeConsistent(), SUMA_malloc, SUMA_MIN_PAIR, SUMA_OrientTriangles(), SUMA_Patch2Surf(), SUMA_qhull_wrap(), SUMA_realloc, SUMA_RECOMPUTE_NORMALS, SUMA_RETURN, SUMA_SL_Crit, SUMA_SL_Err, SUMA_SL_Note, SUMA_SL_Warn, SUMA_SurfaceMetrics(), SUMA_z_doubqsort(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::t, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t98, THD_3dind_to_3dmm(), THD_3dmm_to_dicomm(), THD_BN_XCM, THD_BN_YCM, THD_BN_ZCM, SUMA_GENERIC_PROG_OPTIONS_STRUCT::tm, SUMA_GENERIC_PROG_OPTIONS_STRUCT::VolCM, and THD_fvec3::xyz.
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 }
|
|
||||||||||||||||||||
|
A function to move node surfaces towards a better future. This function used to be a macro gone too large. Definition at line 1789 of file SUMA_BrainWrap.c. References a, SUMA_SurfaceObject::Center, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, SUMA_COMM_STRUCT::kth, LocalHead, SUMA_SurfaceObject::N_Node, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_COMM_STRUCT::Send, SUMA_GENERIC_PROG_OPTIONS_STRUCT::Stop, SUMA_Boolean, SUMA_calloc, SUMA_Chung_Smooth(), SUMA_Chung_Smooth_Weights(), SUMA_COLUMN_MAJOR, SUMA_ENTRY, SUMA_free, SUMA_free2D(), SUMA_IS_EYE_ZONE, SUMA_LH, SUMA_malloc, SUMA_MIN_PAIR, SUMA_NODE_XYZ, SUMA_POINT_AT_DISTANCE, SUMA_RECOMPUTE_NORMALS, SUMA_RETURN, SUMA_ROW_MAJOR, SUMA_SendToSuma(), SUMA_SL_Err, SUMA_SL_Warn, SUMA_Suggest_Touchup(), SUMA_Taubin_Smooth(), and SUMA_UNIT_VEC.
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 }
|
|
||||||||||||||||
|
creates a crude mask of the brain's limit based on the surface of the skull mask is crude, for purposes of keeping BrainWrap from leaking out when there are very strong shading artifacts
Definition at line 572 of file SUMA_BrainWrap.c. References SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, DSET_NX, DSET_NY, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, SUMA_COMM_STRUCT::kth, LocalHead, SUMA_SurfaceObject::N_Node, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_COMM_STRUCT::Send, SUMA_GENERIC_PROG_OPTIONS_STRUCT::send_hull, SUMA_Boolean, SUMA_ENTRY, SUMA_free, SUMA_malloc, SUMA_MAX_NEG_SHISH_JUMP, SUMA_MAX_PAIR, SUMA_NODE_XYZ, SUMA_ONE_SHISH_PLEASE, SUMA_RECOMPUTE_NORMALS, SUMA_RETURN, SUMA_ROW_MAJOR, SUMA_SendToSuma(), SUMA_SL_Crit, SUMA_SL_Warn, SUMA_Taubin_Smooth(), and SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp.
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 }
|
|
||||||||||||||||
|
Main function to stretch a sphere to fit the brain. A modification of BET's algorithm. Definition at line 675 of file SUMA_BrainWrap.c. References a, SUMA_GENERIC_PROG_OPTIONS_STRUCT::bot_lztclip, cbuf, SUMA_SurfaceObject::Center, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dbg_eyenodes, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, SUMA_GENERIC_PROG_OPTIONS_STRUCT::DemoPause, E, SUMA_GENERIC_PROG_OPTIONS_STRUCT::ExpFrac, InteractiveQuit, SUMA_COMM_STRUCT::kth, l, LocalHead, SUMA_GENERIC_PROG_OPTIONS_STRUCT::N_it, SUMA_SurfaceObject::N_Node, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NNsmooth, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NoEyes, r, SUMA_COMM_STRUCT::Send, SUMA_GENERIC_PROG_OPTIONS_STRUCT::smootheach, SUMA_GENERIC_PROG_OPTIONS_STRUCT::Stop, SUMA_GENERIC_PROG_OPTIONS_STRUCT::su1, SUMA_3dSS_DEMO_PAUSE, SUMA_3dSS_INTERACTIVE, SUMA_ABS, SUMA_Boolean, SUMA_DOTP_VEC, SUMA_ENTRY, SUMA_Find_IminImax(), SUMA_free, SUMA_IS_EYE_ZONE, SUMA_LH, SUMA_malloc, SUMA_MAX_PAIR, SUMA_MAX_SHISH_JUMP, SUMA_MEAN_NEIGHB_COORD, SUMA_MEAN_SEGMENT_LENGTH, SUMA_Mesh_Area(), SUMA_NODE_XYZ, SUMA_NORM_VEC, SUMA_PAUSE_PROMPT, SUMA_ReadCharStdin(), SUMA_RECOMPUTE_NORMALS, SUMA_RETURN, SUMA_S_Err, SUMA_SCALE_VEC, SUMA_SendToSuma(), SUMA_SL_Crit, SUMA_SL_Err, SUMA_SL_Warn, SUMA_SUB_VEC, SUMA_Suggest_Touchup(), SUMA_Suggest_Touchup_Grad(), SUMA_WRAP_BRAIN_SMOOTH_NN, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t98, SUMA_GENERIC_PROG_OPTIONS_STRUCT::tm, SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp, SUMA_GENERIC_PROG_OPTIONS_STRUCT::UseExpansion, SUMA_GENERIC_PROG_OPTIONS_STRUCT::var_lzt, and SUMA_GENERIC_PROG_OPTIONS_STRUCT::ztv.
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 }
|
|
||||||||||||||||||||||||
|
A function to recommend node movement towards a better future.
Definition at line 1672 of file SUMA_BrainWrap.c. References a, SUMA_SurfaceObject::Center, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, LocalHead, SUMA_SurfaceObject::N_Node, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_GENERIC_PROG_OPTIONS_STRUCT::r, SUMA_GENERIC_PROG_OPTIONS_STRUCT::Stop, SUMA_Boolean, SUMA_calloc, SUMA_ENTRY, SUMA_Find_IminImax(), SUMA_free, SUMA_malloc, SUMA_MAX_PAIR, SUMA_RETURN, SUMA_SL_Crit, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, and SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp. Referenced by SUMA_Reposition_Touchup(), and SUMA_StretchToFitLeCerveau().
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 }
|
|
||||||||||||||||||||||||
|
A function to recommend node movement towards a better future. FAILS, mesh density too low. Need 3d edge mask....
Definition at line 1522 of file SUMA_BrainWrap.c. References a, SUMA_SurfaceObject::Center, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d1, SUMA_GENERIC_PROG_OPTIONS_STRUCT::d4, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, LocalHead, SUMA_SurfaceObject::N_Node, SUMA_GENERIC_PROG_OPTIONS_STRUCT::NodeDbg, SUMA_SurfaceObject::NodeList, SUMA_SurfaceObject::NodeNormList, SUMA_GENERIC_PROG_OPTIONS_STRUCT::nvox, SUMA_VOLPAR::nx, SUMA_VOLPAR::ny, SUMA_GENERIC_PROG_OPTIONS_STRUCT::r, SUMA_1D_2_3D_index, SUMA_Boolean, SUMA_calloc, SUMA_ENTRY, SUMA_Find_IminImax(), SUMA_free, SUMA_malloc, SUMA_MAX_PAIR, SUMA_NORM_VEC, SUMA_RETURN, SUMA_SL_Crit, SUMA_SL_Warn, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t2, SUMA_GENERIC_PROG_OPTIONS_STRUCT::t98, SUMA_GENERIC_PROG_OPTIONS_STRUCT::tm, SUMA_GENERIC_PROG_OPTIONS_STRUCT::travstp, and SUMA_SurfaceObject::VolPar. Referenced by SUMA_StretchToFitLeCerveau().
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 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 1222 of file SUMA_BrainWrap.c. References DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, SUMA_VOLPAR::dx, SUMA_VOLPAR::dy, SUMA_VOLPAR::dz, EDIT_coerce_scale_type(), SUMA_SurfaceObject::FaceSetDim, SUMA_SurfaceObject::FaceSetList, SUMA_VOLPAR::Hand, LocalHead, n1, n2, SUMA_SurfaceObject::N_FaceSet, SUMA_SurfaceObject::N_Node, SUMA_MT_INTERSECT_TRIANGLE::N_poshits, SUMA_SurfaceObject::NodeDim, SUMA_VOLPAR::nx, SUMA_VOLPAR::ny, nz, SUMA_VOLPAR::nz, p, SUMA_3D_2_1D_index, SUMA_Boolean, SUMA_calloc, SUMA_DIST_FROM_PLANE, SUMA_ENTRY, SUMA_FillToVoxelMask(), SUMA_free, SUMA_Free_MT_intersect_triangle(), SUMA_IN_TRIBOX_INSIDE, SUMA_IN_TRIBOX_OUTSIDE, SUMA_INSIDE_SURFACE, SUMA_INTERSECTS_TRIANGLE_INSIDE, SUMA_INTERSECTS_TRIANGLE_OUTSIDE, SUMA_IS_NEG, SUMA_isVoxelIntersect_Triangle(), SUMA_LH, SUMA_malloc, SUMA_MIN_MAX_SUM_VECMAT_COL, SUMA_MT_intersect_triangle(), SUMA_ON_NODE, SUMA_realloc, SUMA_RETURN, SUMA_SL_Crit, SUMA_SL_Err, SUMA_SL_Warn, SUMA_TRIANGLE_BOUNDING_BOX, SUMA_UNIT_VEC, SUMA_VoxelsInBox(), and THD_mask_fillin_once(). Referenced by SUMA_FindVoxelsInSurface().
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 }
|
Variable Documentation
|
|
Definition at line 27 of file SUMA_BrainWrap.c. Referenced by SUMA_StretchToFitLeCerveau(). |
|
|
Global pointer to structure containing info common to all viewers Definition at line 20 of file SUMA_BrainWrap.c. |
|
|
Global pointer to Displayable Object structure vector Definition at line 21 of file SUMA_BrainWrap.c. |
|
|
Number of DOs stored in DOv Definition at line 24 of file SUMA_BrainWrap.c. |
|
|
Number of SVs stored in SVv Definition at line 23 of file SUMA_BrainWrap.c. |
|
|
Global pointer to the vector containing the various Surface Viewer Structures Definition at line 22 of file SUMA_BrainWrap.c. |
