00001 #include "SUMA_suma.h"
00002 #include "../thd_brainormalize.h"
00003 #include "../rickr/r_new_resam_dset.h"
00004
00005 #undef STAND_ALONE
00006
00007 #if defined SUMA_BrainWrap_STANDALONE
00008 #define STAND_ALONE
00009 #endif
00010
00011 #ifdef STAND_ALONE
00012
00013 SUMA_SurfaceViewer *SUMAg_cSV = NULL;
00014 SUMA_SurfaceViewer *SUMAg_SVv = NULL;
00015
00016 int SUMAg_N_SVv = 0;
00017 SUMA_DO *SUMAg_DOv = NULL;
00018 int SUMAg_N_DOv = 0;
00019 SUMA_CommonFields *SUMAg_CF = NULL;
00020 #else
00021 extern SUMA_CommonFields *SUMAg_CF;
00022 extern SUMA_DO *SUMAg_DOv;
00023 extern SUMA_SurfaceViewer *SUMAg_SVv;
00024 extern int SUMAg_N_SVv;
00025 extern int SUMAg_N_DOv;
00026 #endif
00027
00028 #ifdef SUMA_BrainWrap_STANDALONE
00029 void usage_SUMA_BrainWrap (SUMA_GENERIC_ARGV_PARSE *ps)
00030 {
00031 static char FuncName[]={"usage_SUMA_BrainWrap"};
00032 char * s = NULL, *sio=NULL, *st = NULL, *sts = NULL;
00033 int i;
00034 s = SUMA_help_basics();
00035 sio = SUMA_help_IO_Args(ps);
00036 sts = SUMA_help_talk(ps);
00037 printf ( "\n"
00038 "Usage: A program to extract the brain from surrounding.\n"
00039 " tissue from MRI T1-weighted images. The largely automated\n"
00040 " process consists of three steps:\n"
00041 " 1- Preprocessing of volume to remove gross spatial image \n"
00042 " non-uniformity artifacts and reposition the brain in\n"
00043 " a reasonable manner for convenience.\n"
00044 " 2- Expand a spherical surface iteratively until it envelopes\n"
00045 " the brain. This is a modified version of the BET algorithm:\n"
00046 " Fast robust automated brain extraction, \n"
00047 " by Stephen M. Smith, HBM 2002 v 17:3 pp 143-155\n"
00048 " Modifications include the use of:\n"
00049 " . outer brain surface\n"
00050 " . expansion driven by data inside and outside the surface\n"
00051 " . avoidance of eyes and ventricles\n"
00052 " . a set of operations to avoid the clipping of certain brain\n"
00053 " areas and reduce leakage into the skull in heavily shaded\n"
00054 " data\n"
00055 " . two additional processing stages to ensure convergence and\n"
00056 " reduction of clipped areas.\n"
00057 " 3- The creation of various masks and surfaces modeling brain\n"
00058 " and portions of the skull\n"
00059 "\n"
00060 " 3dSkullStrip < -input VOL >\n"
00061 " [< -o_TYPE PREFIX >] [< -prefix Vol_Prefix >] \n"
00062 " [< -spatnorm >] [< -no_spatnorm >] [< -write_spatnorm >]\n"
00063 " [< -niter N_ITER >] [< -ld LD >] \n"
00064 " [< -shrink_fac SF >] [< -var_shrink_fac >] \n"
00065 " [< -no_var_shrink_fac >] [< shrink_fac_bot_lim SFBL >]\n"
00066 " [< -pushout >] [< -no_pushout >] [< -exp_frac FRAC]\n"
00067 " [< -touchup >] [< -no_touchup >]\n"
00068 " [< -fill_hole R >] [< -NN_smooth NN_SM >]\n"
00069 " [< -smooth_final SM >] [< -avoid_vent >] [< -no_avoid_vent >]\n"
00070 " [< -use_skull >] [< -no_use_skull >] \n"
00071 " [< -avoid_eyes >] [< -no_avoid_eyes >] \n"
00072 " [< -perc_int PERC_INT >] \n"
00073 " [< -max_inter_iter MII >] [-mask_vol]\n"
00074 " [< -debug DBG >] [< -node_dbg NODE_DBG >]\n"
00075 " [< -demo_pause >]\n"
00076 "\n"
00077 " NOTE: Program is in Beta mode, please report bugs and strange failures\n"
00078 " to ziad@nih.gov\n"
00079 "\n"
00080 " Mandatory parameters:\n"
00081 " -input VOL: Input AFNI (or AFNI readable) volume.\n"
00082 " \n"
00083 "\n"
00084 " Optional Parameters:\n"
00085 " -o_TYPE PREFIX: prefix of output surface.\n"
00086 " where TYPE specifies the format of the surface\n"
00087 " and PREFIX is, well, the prefix.\n"
00088 " TYPE is one of: fs, 1d (or vec), sf, ply.\n"
00089 " More on that below.\n"
00090 " -prefix VOL_PREFIX: prefix of output volume.\n"
00091 " If not specified, the prefix is the same\n"
00092 " as the one used with -o_TYPE.\n"
00093 " The output volume is skull stripped version\n"
00094 " of the input volume. In the earlier version\n"
00095 " of the program, a mask volume was written out.\n"
00096 " You can still get that mask volume instead of the\n"
00097 " skull-stripped volume with the option -mask_vol . \n"
00098 " -mask_vol: Output a mask volume instead of a skull-stripped\n"
00099 " volume.\n"
00100 " The mask volume containes:\n"
00101 " 0 : outside the brain\n"
00102 " 1 : voxels near brain surface on the outside.\n"
00103 " 2 : voxels containing a node of the brain surface.\n"
00104 " 3 : voxels near brain surface on the inside.\n"
00105 " 4 : voxels inside the brain.\n"
00106 " -spat_norm: (Default) Perform spatial normalization first.\n"
00107 " This is a necessary step unless the volume has\n"
00108 " been 'spatnormed' already.\n"
00109 " -no_spatnorm: Do not perform spatial normalization.\n"
00110 " Use this option only when the volume \n"
00111 " has been run through the 'spatnorm' process\n"
00112 " -spatnorm_dxyz DXYZ: Use DXY for the spatial resolution of the\n"
00113 " spatially normalized volume. The default \n"
00114 " is the lowest of all three dimensions.\n"
00115 " For human brains, use DXYZ of 1.0, for\n"
00116 " primate brain, use the default setting.\n"
00117 " -write_spatnorm: Write the 'spatnormed' volume to disk.\n"
00118 " -niter N_ITER: Number of iterations. Default is 250\n"
00119 " For denser meshes, you need more iterations\n"
00120 " N_ITER of 750 works for LD of 50.\n"
00121 " -ld LD: Parameter to control the density of the surface.\n"
00122 " Default is 20. See CreateIcosahedron -help\n"
00123 " for details on this option.\n"
00124 " -shrink_fac SF: Parameter controlling the brain vs non-brain\n"
00125 " intensity threshold (tb). Default is 0.6.\n"
00126 " tb = (Imax - t2) SF + t2 \n"
00127 " where t2 is the 2 percentile value and Imax is the local\n"
00128 " maximum, limited to the median intensity value.\n"
00129 " For more information on tb, t2, etc. read the BET paper\n"
00130 " mentioned above. Note that in 3dSkullStrip, SF can vary across \n"
00131 " iterations and might be automatically clipped in certain areas.\n"
00132 " SF can vary between 0 and 1.\n"
00133 " 0: Intensities < median inensity are considered non-brain\n"
00134 " 1: Intensities < t2 are considered non-brain\n"
00135 " -var_shrink_fac: Vary the shrink factor with the number of\n"
00136 " iterations. This reduces the likelihood of a surface\n"
00137 " getting stuck on large pools of CSF before reaching\n"
00138 " the outer surface of the brain. (Default)\n"
00139 " -no_var_shrink_fac: Do not use var_shrink_fac.\n"
00140 " -shrink_fac_bot_lim SFBL: Do not allow the varying SF to go\n"
00141 " below SFBL . Default 0.65. \n"
00142 " This option helps reduce potential for leakage below \n"
00143 " the cerebellum.\n"
00144 " -pushout: Consider values above each node in addition to values\n"
00145 " below the node when deciding on expansion. (Default)\n"
00146 " -no_pushout: Do not use -pushout.\n"
00147 " -exp_frac FRAC: Speed of expansion (see BET paper). Default is 0.1.\n"
00148 " -touchup: Perform touchup operations at end to include\n"
00149 " areas not covered by surface expansion. (Default)\n"
00150 " -no_touchup: Do not use -touchup\n"
00151 " -fill_hole R: Fill small holes that can result from small surface\n"
00152 " intersections caused by the touchup operation.\n"
00153 " R is the maximum number of pixels on the side of a hole\n"
00154 " that can be filled. Big holes are not filled.\n"
00155 " If you use -touchup, the default R is 10. Otherwise \n"
00156 " the default is 0.\n"
00157 " This is a less than elegant solution to the small\n"
00158 " intersections which are usually eliminated\n"
00159 " automatically. \n"
00160 " -NN_smooth NN_SM: Perform Nearest Neighbor coordinate interpolation\n"
00161 " every few iterations. Default is 72\n"
00162 " -smooth_final SM: Perform final surface smoothing after all iterations.\n"
00163 " Default is 20 smoothing iterations.\n"
00164 " Smoothing is done using Taubin's method, \n"
00165 " see SurfSmooth -help for detail.\n"
00166 " -avoid_vent: avoid ventricles. Default.\n"
00167 " -no_avoid_vent: Do not use -avoid_vent.\n"
00168 " -avoid_eyes: avoid eyes. Default\n"
00169 " -no_avoid_eyes: Do not use -avoid_eyes.\n"
00170 " -use_skull: Use outer skull to limit expansion of surface into\n"
00171 " the skull due to very strong shading artifacts.\n"
00172 " This option is buggy at the moment, use it only \n"
00173 " if you have leakage into skull.\n"
00174 " -no_use_skull: Do not use -use_skull (Default).\n"
00175 " -send_no_skull: Do not send the skull surface to SUMA if you are\n"
00176 " using -talk_suma\n"
00177 " -perc_int PERC_INT: Percentage of segments allowed to intersect\n"
00178 " surface. Ideally this should be 0 (Default). \n"
00179 " However, few surfaces might have small stubborn\n"
00180 " intersections that produce a few holes.\n"
00181 " PERC_INT should be a small number, typically\n"
00182 " between 0 and 0.1\n"
00183 " -max_inter_iter N_II: Number of iteration to remove intersection\n"
00184 " problems. With each iteration, the program\n"
00185 " automatically increases the amount of smoothing\n"
00186 " to get rid of intersections. Default is 4\n"
00187 " -blur_fwhm FWHM: Blur dset after spatial normalization.\n"
00188 " Recommended when you have lots of CSF in brain\n"
00189 " and when you have protruding gyri (finger like)\n"
00190 " Recommended value is 2..4. \n"
00191 " -interactive: Make the program stop at various stages in the \n"
00192 " segmentation process for a prompt from the user\n"
00193 " to continue or skip that stage of processing.\n"
00194 " This option is best used in conjunction with options\n"
00195 " -talk_suma and -feed_afni\n"
00196 " -demo_pause: Pause at various step in the process to facilitate\n"
00197 " interactive demo while 3dSkullStrip is communicating\n"
00198 " with AFNI and SUMA. See 'Eye Candy' mode below and\n"
00199 " -talk_suma option. \n"
00200 "\n"
00201 "%s"
00202 " -visual: Equivalent to using -talk_suma -feed_afni -send_kth 5\n"
00203 "\n"
00204 " -debug DBG: debug levels of 0 (default), 1, 2, 3.\n"
00205 " This is no Rick Reynolds debug, which is oft nicer\n"
00206 " than the results, but it will do.\n"
00207 " -node_dbg NODE_DBG: Output lots of parameters for node\n"
00208 " NODE_DBG for each iteration.\n"
00209 "\n"
00210 " Tips:\n"
00211 " I ran the program with the default parameters on 200+ datasets.\n"
00212 " The results were quite good in all but a couple of instances, here\n"
00213 " are some tips on fixing trouble spots:\n"
00214 "\n"
00215 " Clipping in frontal areas, close to the eye balls:\n"
00216 " + Try -no_avoid_eyes option\n"
00217 " Clipping in general:\n"
00218 " + Use lower -shrink_fac, start with 0.5 then 0.4\n"
00219 " Some lobules are not included:\n"
00220 " + Use a denser mesh (like -ld 50) and increase iterations \n"
00221 " (-niter 750). The program will take much longer to run in that case.\n"
00222 " + Instead of using denser meshes, you could try blurring the data \n"
00223 " before skull stripping. Something like 3dmerge -1blur_fwhm 2 did\n"
00224 " wonders for some of my data with the default options of 3dSkullStrip\n"
00225 " Blurring is a lot faster than increasing mesh density.\n"
00226 " + Use also a smaller -shrink_fac is you have lots of CSF between gyri.\n"
00227 "\n"
00228 " Eye Candy Mode: (previous restrictions removed)\n"
00229 " You can run BrainWarp and have it send successive iterations\n"
00230 " to SUMA and AFNI. This is very helpful in following the\n"
00231 " progression of the algorithm and determining the source\n"
00232 " of trouble, if any.\n"
00233 " Example:\n"
00234 " afni -niml -yesplugouts &\n"
00235 " suma -niml &\n"
00236 " 3dSkullStrip -input Anat+orig -o_ply anat_brain -talk_suma -feed_afni -send_kth 5\n"
00237 "\n"
00238
00239
00240
00241
00242
00243
00244
00245
00246 "%s"
00247 "\n", sio, s);
00248 SUMA_free(s); s = NULL; SUMA_free(st); st = NULL; SUMA_free(sio); sio = NULL; SUMA_free(sts); sts = NULL;
00249 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
00250 printf(" Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \n");
00251 exit (0);
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_BrainWrap_ParseInput (char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
00264 {
00265 static char FuncName[]={"SUMA_BrainWrap_ParseInput"};
00266 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
00267 int kar, i, ind, exists;
00268 char *outname, cview[10];
00269 SUMA_Boolean brk = NOPE;
00270 SUMA_Boolean LocalHead = NOPE;
00271
00272 SUMA_ENTRY;
00273
00274 Opt = SUMA_Alloc_Generic_Prog_Options_Struct();
00275
00276 kar = 1;
00277 Opt->SpatNormDxyz = 0.0;
00278 Opt->spec_file = NULL;
00279 Opt->out_vol_prefix = NULL;
00280 Opt->out_prefix = NULL;
00281 Opt->sv_name = NULL;
00282 Opt->N_surf = -1;
00283 Opt->in_name = NULL;
00284 Opt->cmask = NULL;
00285 Opt->MaskMode = 0;
00286 for (i=0; i<SUMA_GENERIC_PROG_MAX_SURF; ++i) { Opt->surf_names[i] = NULL; }
00287 outname = NULL;
00288 Opt->in_vol = NULL;
00289 Opt->nvox = -1;
00290 Opt->ninmask = -1;
00291 Opt->mcdatav = NULL;
00292 Opt->debug = 0;
00293 Opt->v0 = 0.0;
00294 Opt->v1 = -1.0;
00295 Opt->dvec = NULL;
00296 Opt->SurfFileType = SUMA_PLY;
00297 Opt->SurfFileFormat = SUMA_ASCII;
00298 Opt->xform = SUMA_ISO_XFORM_MASK;
00299 Opt->obj_type = -1;
00300 Opt->obj_type_res = -1;
00301 Opt->XYZ = NULL;
00302 Opt->in_1D = NULL;
00303 Opt->N_XYZ = 0;
00304 Opt->Zt = 0.6;
00305 Opt->ExpFrac = 0.1;
00306 Opt->N_it = 250;
00307 Opt->Icold = 20;
00308 Opt->NodeDbg = -1;
00309 Opt->t2 = Opt->t98 = Opt->t = Opt->tm = -1;
00310 Opt->r = 0;
00311 Opt->d1 = 20;
00312 Opt->su1 = 1;
00313 Opt->UseNew = 1.0;
00314 Opt->d4 = 15;
00315 Opt->ztv = NULL;
00316 Opt->Kill98 = 0;
00317 Opt->NoEyes = 1;
00318 Opt->NNsmooth = 72;
00319 Opt->smootheach = 50;
00320 Opt->avoid_vent = 1;
00321 Opt->smooth_end = 20;
00322 Opt->k98mask = NULL;
00323 Opt->k98maskcnt = 0;
00324 Opt->dbg_eyenodes = NULL;
00325 Opt->travstp = 1.0;
00326 Opt->Stop = NULL;
00327 Opt->MaxIntIter = 4;
00328 Opt->UseExpansion = 1;
00329 Opt->PercInt = 0;
00330 Opt->UseSkull = 0;
00331 Opt->send_hull = 1;
00332 Opt->bot_lztclip = 0.65;
00333 Opt->var_lzt = 1.0;
00334
00335 Opt->DemoPause = SUMA_3dSS_NO_PAUSE;
00336 Opt->DoSpatNorm = 1;
00337 Opt->WriteSpatNorm = 0;
00338 Opt->fillhole = -1;
00339 Opt->iset = NULL;
00340 Opt->SpatShift[0] = Opt->SpatShift[1] = Opt->SpatShift[2] = 0.0;
00341 Opt->OrigSpatNormedSet = NULL;
00342 Opt->in_edvol = NULL;
00343 Opt->blur_fwhm = 0.0;
00344 brk = NOPE;
00345 while (kar < argc) {
00346
00347 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
00348 usage_SUMA_BrainWrap(ps);
00349 exit (0);
00350 }
00351
00352 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
00353
00354 if (!brk && (strcmp(argv[kar], "-pushout") == 0)) {
00355 Opt->UseExpansion = 1;
00356 brk = YUP;
00357 }
00358
00359 if (!brk && (strcmp(argv[kar], "-visual") == 0)) {
00360 ps->cs->talk_suma = 1;
00361 ps->cs->kth = 5;
00362 ps->cs->Feed2Afni = 1;
00363 brk = YUP;
00364 }
00365
00366 if (!brk && (strcmp(argv[kar], "-no_pushout") == 0)) {
00367 Opt->UseExpansion = 0;
00368 brk = YUP;
00369 }
00370
00371 if (!brk && (strcmp(argv[kar], "-spatnorm") == 0)) {
00372 Opt->DoSpatNorm = 1;
00373 brk = YUP;
00374 }
00375
00376 if (!brk && (strcmp(argv[kar], "-no_spatnorm") == 0)) {
00377 Opt->DoSpatNorm = 0;
00378 brk = YUP;
00379 }
00380
00381 if (!brk && (strcmp(argv[kar], "-write_spatnorm") == 0)) {
00382 Opt->WriteSpatNorm = 1;
00383 brk = YUP;
00384 }
00385
00386 if (!brk && (strcmp(argv[kar], "-use_skull") == 0)) {
00387 Opt->UseSkull = 1;
00388 brk = YUP;
00389 }
00390
00391 if (!brk && (strcmp(argv[kar], "-no_use_skull") == 0)) {
00392 Opt->UseSkull = 0;
00393 brk = YUP;
00394 }
00395
00396 if (!brk && (strcmp(argv[kar], "-send_no_skull") == 0)) {
00397 Opt->send_hull = 0;
00398 brk = YUP;
00399 }
00400
00401
00402 if (!brk && (strcmp(argv[kar], "-var_shrink_fac") == 0)) {
00403 Opt->var_lzt = 1;
00404 brk = YUP;
00405 }
00406
00407 if (!brk && (strcmp(argv[kar], "-no_var_shrink_fac") == 0)) {
00408 Opt->var_lzt = 0;
00409 brk = YUP;
00410 }
00411
00412
00413 if (!brk && (strcmp(argv[kar], "-fill_hole") == 0)) {
00414 kar ++;
00415 if (kar >= argc) {
00416 fprintf (SUMA_STDERR, "need argument after -fill_hole \n");
00417 exit (1);
00418 }
00419 Opt->fillhole = atoi(argv[kar]);
00420 if ( (Opt->fillhole < 0 || Opt->fillhole > 50) ) {
00421 fprintf (SUMA_STDERR, "parameter after -fill_hole (%d) should be >= 0 and <= 50 \n", Opt->fillhole);
00422 exit (1);
00423 }
00424 brk = YUP;
00425 }
00426
00427 if (!brk && (strcmp(argv[kar], "-perc_int") == 0)) {
00428 kar ++;
00429 if (kar >= argc) {
00430 fprintf (SUMA_STDERR, "need argument after -perc_int \n");
00431 exit (1);
00432 }
00433 Opt->PercInt = atof(argv[kar]);
00434 if ( (Opt->PercInt < 0 && Opt->PercInt != -1) || Opt->PercInt > 10) {
00435 fprintf (SUMA_STDERR, "parameter after -perc_int should be -1 or between 0 and 10 (have %f) \n", Opt->PercInt);
00436 exit (1);
00437 }
00438 brk = YUP;
00439 }
00440
00441 if (!brk && (strcmp(argv[kar], "-spatnorm_dxyz") == 0)) {
00442 kar ++;
00443 if (kar >= argc) {
00444 fprintf (SUMA_STDERR, "need argument after -spatnorm_dxyz \n");
00445 exit (1);
00446 }
00447 Opt->SpatNormDxyz = atof(argv[kar]);
00448 if ( Opt->SpatNormDxyz < 0 || Opt->SpatNormDxyz > 10) {
00449 fprintf (SUMA_STDERR, "parameter after -spatnorm_dxyz should be between 0 and 10 (have %f) \n",
00450 Opt->SpatNormDxyz);
00451 exit (1);
00452 }
00453 brk = YUP;
00454 }
00455
00456 if (!brk && (strcmp(argv[kar], "-blur_fwhm") == 0)) {
00457 kar ++;
00458 if (kar >= argc) {
00459 fprintf (SUMA_STDERR, "need argument after -blur_fwhm \n");
00460 exit (1);
00461 }
00462 Opt->blur_fwhm = atof(argv[kar]);
00463 if ( Opt->PercInt < 0 || Opt->PercInt > 50) {
00464 fprintf (SUMA_STDERR, "parameter after -blur_fwhm should be between 0 and 50 (have %f) \n", Opt->blur_fwhm);
00465 exit (1);
00466 }
00467 brk = YUP;
00468 }
00469
00470 if (!brk && (strcmp(argv[kar], "-input_1D") == 0)) {
00471 kar ++;
00472 if (kar >= argc) {
00473 fprintf (SUMA_STDERR, "need argument after -input_1D \n");
00474 exit (1);
00475 }
00476 Opt->in_1D = argv[kar];
00477 brk = YUP;
00478 }
00479
00480 if (!brk && (strcmp(argv[kar], "-avoid_vent") == 0)) {
00481 Opt->avoid_vent = 1;
00482 brk = YUP;
00483 }
00484
00485 if (!brk && (strcmp(argv[kar], "-no_avoid_vent") == 0)) {
00486 Opt->avoid_vent = 0;
00487 brk = YUP;
00488 }
00489
00490 if (!brk && (strcmp(argv[kar], "-demo_pause") == 0)) {
00491 Opt->DemoPause = SUMA_3dSS_DEMO_PAUSE;
00492 brk = YUP;
00493 }
00494
00495 if (!brk && (strcmp(argv[kar], "-interactive") == 0)) {
00496 Opt->DemoPause = SUMA_3dSS_INTERACTIVE;
00497 brk = YUP;
00498 }
00499
00500 if (!brk && (strcmp(argv[kar], "-d1") == 0)) {
00501 kar ++;
00502 if (kar >= argc) {
00503 fprintf (SUMA_STDERR, "need argument after -d1 \n");
00504 exit (1);
00505 }
00506 Opt->d1 = atof(argv[kar]);
00507 brk = YUP;
00508 }
00509
00510 if (!brk && (strcmp(argv[kar], "-max_inter_iter") == 0)) {
00511 kar ++;
00512 if (kar >= argc) {
00513 fprintf (SUMA_STDERR, "need argument after -max_inter_iter \n");
00514 exit (1);
00515 }
00516 Opt->MaxIntIter = atoi(argv[kar]);
00517 if (Opt->MaxIntIter < 0 || Opt->MaxIntIter > 10) {
00518 fprintf (SUMA_STDERR, "-max_inter_iter parameter should be between 0 and 10 \n");
00519 exit (1);
00520 }
00521 brk = YUP;
00522 }
00523 if (!brk && (strcmp(argv[kar], "-smooth_final") == 0)) {
00524 kar ++;
00525 if (kar >= argc) {
00526 fprintf (SUMA_STDERR, "need argument after -smooth_final \n");
00527 exit (1);
00528 }
00529 Opt->smooth_end = atoi(argv[kar]);
00530 if (Opt->smooth_end < 0 || Opt->smooth_end > 100 || Opt->smooth_end % 2) {
00531 fprintf (SUMA_STDERR, "irrational or exhuberant values for smooth_final\n Must use even number between 0 to 100 (%d entered)\n", Opt->smooth_end);
00532 exit (1);
00533 }
00534 brk = YUP;
00535 }
00536 if (!brk && (strcmp(argv[kar], "-NNsmooth") == 0)) {
00537 kar ++;
00538 if (kar >= argc) {
00539 fprintf (SUMA_STDERR, "need argument after -NNsmooth \n");
00540 exit (1);
00541 }
00542 Opt->NNsmooth = atoi(argv[kar]);
00543 if (Opt->NNsmooth % 2) {
00544 fprintf (SUMA_STDERR, "Number of smoothing steps must be an even positive number (not %d) \n", Opt->NNsmooth);
00545 exit (1);
00546 }
00547 brk = YUP;
00548 }
00549
00550 if (!brk && ((strcmp(argv[kar], "-shrink_fac") == 0) || (strcmp(argv[kar], "-bf") == 0))) {
00551 kar ++;
00552 if (kar >= argc) {
00553 fprintf (SUMA_STDERR, "need argument after -shrink_fac || -bf \n");
00554 exit (1);
00555 }
00556 Opt->Zt = atof(argv[kar]);
00557 brk = YUP;
00558 }
00559
00560 if (!brk && (strcmp(argv[kar], "-shrink_fac_bot_lim") == 0)) {
00561 kar ++;
00562 if (kar >= argc) {
00563 fprintf (SUMA_STDERR, "need argument after -shrink_fac_bot_lim \n");
00564 exit (1);
00565 }
00566 Opt->bot_lztclip = atof(argv[kar]);
00567 brk = YUP;
00568 }
00569
00570 if (!brk && (strcmp(argv[kar], "-t2") == 0)) {
00571 kar ++;
00572 if (kar >= argc) {
00573 fprintf (SUMA_STDERR, "need argument after -t2 \n");
00574 exit (1);
00575 }
00576 Opt->t2 = atof(argv[kar]);
00577 brk = YUP;
00578 }
00579 if (!brk && (strcmp(argv[kar], "-t98") == 0)) {
00580 kar ++;
00581 if (kar >= argc) {
00582 fprintf (SUMA_STDERR, "need argument after -t98 \n");
00583 exit (1);
00584 }
00585 Opt->t98 = atof(argv[kar]);
00586 brk = YUP;
00587 }
00588 if (!brk && (strcmp(argv[kar], "-t") == 0)) {
00589 kar ++;
00590 if (kar >= argc) {
00591 fprintf (SUMA_STDERR, "need argument after -t \n");
00592 exit (1);
00593 }
00594 Opt->t = atof(argv[kar]);
00595 brk = YUP;
00596 }
00597 if (!brk && (strcmp(argv[kar], "-tm") == 0)) {
00598 kar ++;
00599 if (kar >= argc) {
00600 fprintf (SUMA_STDERR, "need argument after -tm \n");
00601 exit (1);
00602 }
00603 Opt->tm = atof(argv[kar]);
00604 brk = YUP;
00605 }
00606
00607 if (!brk && (strcmp(argv[kar], "-niter") == 0)) {
00608 kar ++;
00609 if (kar >= argc) {
00610 fprintf (SUMA_STDERR, "need argument after -niter \n");
00611 exit (1);
00612 }
00613 Opt->N_it = atoi(argv[kar]);
00614 if (Opt->N_it < 0 || Opt->N_it > 100000) {
00615 fprintf (SUMA_STDERR, "%d is a bad iteration number.\n", Opt->N_it);
00616 exit (1);
00617 }
00618 brk = YUP;
00619 }
00620
00621 if (!brk && (strcmp(argv[kar], "-exp_frac") == 0)) {
00622 kar ++;
00623 if (kar >= argc) {
00624 fprintf (SUMA_STDERR, "need argument after -exp_frac \n");
00625 exit (1);
00626 }
00627 Opt->ExpFrac = atof(argv[kar]);
00628 if (Opt->ExpFrac < -1 || Opt->ExpFrac > 1) {
00629 fprintf (SUMA_STDERR, "%f is a bad expantion fraction.\n", Opt->ExpFrac);
00630 exit (1);
00631 }
00632 brk = YUP;
00633 }
00634 if (!brk && (strcmp(argv[kar], "-node_dbg") == 0)) {
00635 kar ++;
00636 if (kar >= argc) {
00637 fprintf (SUMA_STDERR, "need argument after -node_dbg \n");
00638 exit (1);
00639 }
00640 Opt->NodeDbg = atoi(argv[kar]);
00641 if (Opt->NodeDbg < 0) {
00642 fprintf (SUMA_STDERR, "%d is a bad node number.\n", Opt->NodeDbg);
00643 exit (1);
00644 }
00645 brk = YUP;
00646 }
00647
00648 if (!brk && (strcmp(argv[kar], "-sm_fac") == 0)) {
00649 kar ++;
00650 if (kar >= argc) {
00651 fprintf (SUMA_STDERR, "need argument after -sm_fac \n");
00652 exit (1);
00653 }
00654 Opt->su1 = atof(argv[kar]);
00655 if (Opt->su1 < 0 || Opt->su1 > 1.5) {
00656 fprintf (SUMA_STDERR, "%f is a bad smoothness factor (acceptable range is 0 to 1\nbut values up to 1.5 are allowed for amusement).\n", Opt->su1);
00657 exit (1);
00658 }
00659 brk = YUP;
00660 }
00661
00662 if (!brk && (strcmp(argv[kar], "-ld") == 0)) {
00663 kar ++;
00664 if (kar >= argc) {
00665 fprintf (SUMA_STDERR, "need argument after -ld \n");
00666 exit (1);
00667 }
00668 Opt->Icold = atoi(argv[kar]);
00669 if (Opt->Icold < 0 || Opt->Icold > 300) {
00670 fprintf (SUMA_STDERR, "%d is a bad iteration number.\n", Opt->Icold);
00671 exit (1);
00672 }
00673 brk = YUP;
00674 }
00675
00676 if (!brk && (strcmp(argv[kar], "-debug") == 0)) {
00677 kar ++;
00678 if (kar >= argc) {
00679 fprintf (SUMA_STDERR, "need argument after -debug \n");
00680 exit (1);
00681 }
00682 Opt->debug = atoi(argv[kar]);
00683 if (Opt->debug > 2) { LocalHead = YUP; }
00684 brk = YUP;
00685 }
00686
00687
00688 if (!brk && (strcmp(argv[kar], "-input") == 0)) {
00689 kar ++;
00690 if (kar >= argc) {
00691 fprintf (SUMA_STDERR, "need argument after -input \n");
00692 exit (1);
00693 }
00694 Opt->in_name = SUMA_copy_string(argv[kar]);
00695 brk = YUP;
00696 }
00697
00698 if (!brk && (strcmp(argv[kar], "-prefix") == 0)) {
00699 kar ++;
00700 if (kar >= argc) {
00701 fprintf (SUMA_STDERR, "need argument after -prefix \n");
00702 exit (1);
00703 }
00704 Opt->out_vol_prefix = SUMA_copy_string(argv[kar]);
00705 brk = YUP;
00706 }
00707
00708 if (!brk && ( (strcmp(argv[kar], "-mask_vol") == 0) ) ) {
00709 Opt->MaskMode = 1;
00710 brk = YUP;
00711 }
00712
00713 if (!brk && ( (strcmp(argv[kar], "-touchup") == 0) ) ) {
00714 Opt->UseNew = 1.0;
00715 brk = YUP;
00716 }
00717
00718 if (!brk && ( (strcmp(argv[kar], "-no_touchup") == 0) ) ) {
00719 Opt->UseNew = 0.0;
00720 brk = YUP;
00721 }
00722
00723 if (!brk && (strcmp(argv[kar], "-k98") == 0)) {
00724 SUMA_SL_Warn("Bad option, causes trouble (big clipping) in other parts of brain.");
00725 Opt->Kill98 = 1;
00726 brk = YUP;
00727 }
00728
00729 if (!brk && ( (strcmp(argv[kar], "-noeyes") == 0) || (strcmp(argv[kar], "-no_eyes") == 0) || (strcmp(argv[kar], "-avoid_eyes") == 0)) ) {
00730 Opt->NoEyes = 1;
00731 brk = YUP;
00732 }
00733
00734 if (!brk && ( (strcmp(argv[kar], "-no_noeyes") == 0) || (strcmp(argv[kar], "-no_no_eyes") == 0) || (strcmp(argv[kar], "-no_avoid_eyes") == 0)) ) {
00735 Opt->NoEyes = 0;
00736 brk = YUP;
00737 }
00738
00739 if (!brk && !ps->arg_checked[kar]) {
00740 fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00741 exit (1);
00742 } else {
00743 brk = NOPE;
00744 kar ++;
00745 }
00746 }
00747
00748 if (Opt->fillhole < 0) {
00749 if (Opt->UseExpansion) {
00750 if (Opt->debug) {
00751 SUMA_SL_Note("Setting fill_hole to 10");
00752 }
00753 Opt->fillhole = 10;
00754 } else Opt->fillhole = 0;
00755 }
00756
00757
00758 if (ps->o_N_surfnames) {
00759 Opt->out_prefix = SUMA_copy_string(ps->o_surfnames[0]);
00760 Opt->SurfFileType = ps->o_FT[0];
00761 Opt->SurfFileFormat = ps->o_FF[0];
00762 }
00763
00764 if (!Opt->in_name) {
00765 fprintf (SUMA_STDERR,"Error %s:\n-input must be used.\n", FuncName);
00766 exit(1);
00767 }
00768
00769 if (!SUMA_AfniView (Opt->in_name, cview)) {
00770 fprintf (SUMA_STDERR,"Error %s:\nCould not guess view from input dset %s\n", FuncName, Opt->in_name);
00771 exit(1);
00772 }
00773
00774 if (!Opt->out_prefix) Opt->out_prefix = SUMA_copy_string("skull_strip_out");
00775 if (!Opt->out_vol_prefix) {
00776 if (!Opt->out_prefix) Opt->out_vol_prefix = SUMA_AfniPrefix("skull_strip_out", NULL, NULL, &exists);
00777 else Opt->out_vol_prefix = SUMA_AfniPrefix(Opt->out_prefix, NULL, NULL, &exists);
00778 } else {
00779 char *stmp = Opt->out_vol_prefix;
00780 Opt->out_vol_prefix = SUMA_AfniPrefix(stmp, NULL, NULL, &exists);
00781 SUMA_free(stmp); stmp = NULL;
00782 }
00783 if (SUMA_AfniExistsView(exists, cview)) {
00784 fprintf (SUMA_STDERR,"Error %s:\nOutput dset %s%s exists, will not overwrite\n", FuncName, Opt->out_vol_prefix, cview);
00785 exit(1);
00786 }
00787
00788
00789 if (Opt->t2 >= 0 || Opt->t98 >= 0 || Opt->tm >= 0 || Opt->t >= 0 ) {
00790 if (!(Opt->t2 >= 0 && Opt->t98 > 0 && Opt->tm > 0 && Opt->t > 0)){
00791 fprintf (SUMA_STDERR,"Error %s: \nAll or none of the t parameters are to be specified on command line:\nt2 %f t %f tm %f t98 %f\n",
00792 FuncName, Opt->t2, Opt->t, Opt->tm, Opt->t98);
00793 exit(1);
00794 }
00795 }
00796 SUMA_RETURN(Opt);
00797 }
00798
00799
00800
00801
00802 int main (int argc,char *argv[])
00803 {
00804 static char FuncName[]={"3dSkullStrip"};
00805 int i, N_in = 0, i3, kth_buf, hull_ld;
00806 int ii,jj,kk,ll,ijk , nx,ny,nz , nn, nint = 0 , nseg;
00807 void *SO_name=NULL, *SO_name_hull=NULL;
00808 float vol, *isin_float=NULL, pint, *dsmooth = NULL, XYZrai_shift[3];
00809 SUMA_SurfaceObject *SO = NULL, *SOhull=NULL;
00810 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
00811 char stmp[200], stmphull[200], *hullprefix=NULL, *prefix=NULL, *spatprefix=NULL, cbuf;
00812 SUMA_Boolean exists = NOPE;
00813 SUMA_GENERIC_ARGV_PARSE *ps=NULL;
00814 short *isin = NULL;
00815 SUMA_FORM_AFNI_DSET_STRUCT *OptDs = NULL;
00816 THD_3dim_dataset *dset = NULL;
00817 THD_ivec3 orixyz , nxyz ;
00818 THD_fvec3 dxyz , orgxyz , fv2, originRAIfv;
00819 THD_3dim_dataset *oset = NULL;
00820 MRI_IMAGE *imin=NULL, *imout=NULL, *imout_orig=NULL , *imout_edge=NULL ;
00821
00822 SUMA_Boolean LocalHead = NOPE;
00823
00824 SUMA_mainENTRY;
00825 SUMA_STANDALONE_INIT;
00826
00827
00828
00829 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
00830 ps = SUMA_Parse_IO_Args(argc, argv, "-o;-talk;");
00831
00832 if (argc < 2) {
00833 usage_SUMA_BrainWrap(ps);
00834 exit (1);
00835 }
00836
00837 Opt = SUMA_BrainWrap_ParseInput (argv, argc, ps);
00838
00839 if (Opt->debug > 2) LocalHead = YUP;
00840
00841 SO_name = SUMA_Prefix2SurfaceName(Opt->out_prefix, NULL, NULL, Opt->SurfFileType, &exists);
00842 if (exists && strcmp(Opt->out_prefix,"skull_strip_out")) {
00843 fprintf(SUMA_STDERR,"Error %s:\nOutput file(s) %s* on disk.\nWill not overwrite.\n", FuncName, Opt->out_prefix);
00844 exit(1);
00845 }
00846 hullprefix = SUMA_append_string(Opt->out_prefix,"_hull");
00847 SO_name_hull = SUMA_Prefix2SurfaceName(hullprefix, NULL, NULL, Opt->SurfFileType, &exists);
00848 if (exists) {
00849 fprintf(SUMA_STDERR,"Error %s:\nOutput file(s) %s* on disk.\nWill not overwrite.\n", FuncName, hullprefix);
00850 exit(1);
00851 }
00852 kth_buf = ps->cs->kth;
00853
00854
00855 #if 0
00856 SUMA_SL_Note("Debugging for eye nodes");
00857 Opt->dbg_eyenodes = fopen("eyenodes.1D.dset", "w");
00858 #endif
00859
00860
00861 if (Opt->DoSpatNorm) {
00862 if (Opt->debug) SUMA_SL_Note("Loading dset, performing Spatial Normalization");
00863
00864 Opt->iset = THD_open_dataset( Opt->in_name );
00865 if( !ISVALID_DSET(Opt->iset) ){
00866 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->in_name) ;
00867 exit(1);
00868 }
00869 Opt->iset_hand = SUMA_THD_handedness( Opt->iset );
00870 if (LocalHead) fprintf(SUMA_STDERR,"%s: Handedness of orig dset %d\n", FuncName, Opt->iset_hand);
00871
00872
00873 imin = THD_median_brick( Opt->iset ) ;
00874 if( imin == NULL ){
00875 fprintf(stderr,"**ERROR: can't load dataset %s\n",Opt->in_name) ;
00876 exit(1);
00877 }
00878
00879 if (Opt->SpatNormDxyz) {
00880 if (Opt->debug) SUMA_S_Note("Overriding default resampling");
00881 mri_brainormalize_initialize(Opt->SpatNormDxyz, Opt->SpatNormDxyz, Opt->SpatNormDxyz);
00882 } else {
00883 mri_brainormalize_initialize(Opt->iset->daxes->xxdel, Opt->iset->daxes->yydel, Opt->iset->daxes->zzdel);
00884 }
00885 imin->dx = fabs(Opt->iset->daxes->xxdel) ;
00886 imin->dy = fabs(Opt->iset->daxes->yydel) ;
00887 imin->dz = fabs(Opt->iset->daxes->zzdel) ;
00888
00889 DSET_unload( Opt->iset ) ;
00890
00891
00892 if( DSET_BRICK_TYPE(Opt->iset,0) == MRI_short ||
00893 DSET_BRICK_TYPE(Opt->iset,0) == MRI_byte ){
00894
00895 imout = mri_to_short(1.0,imin) ;
00896 mri_free(imin) ; imin = imout ;
00897 }
00898
00899
00900 mri_brainormalize_verbose( Opt->debug ) ;
00901 if (Opt->fillhole) {
00902 imout = mri_brainormalize( imin , Opt->iset->daxes->xxorient,
00903 Opt->iset->daxes->yyorient,
00904 Opt->iset->daxes->zzorient, &imout_orig, NULL) ;
00905 } else {
00906 imout = mri_brainormalize( imin , Opt->iset->daxes->xxorient,
00907 Opt->iset->daxes->yyorient,
00908 Opt->iset->daxes->zzorient, NULL, NULL) ;
00909 }
00910 mri_free( imin ) ;
00911
00912 if( imout == NULL ){
00913 fprintf(stderr,"**ERROR: normalization fails!?\n"); exit(1);
00914 }
00915
00916 if (imout_orig) {
00917 if (Opt->debug) SUMA_SL_Note("Creating an output dataset in original grid...");
00918
00919 LOAD_FVEC3(originRAIfv , Opt->iset->daxes->xxorg , Opt->iset->daxes->yyorg , Opt->iset->daxes->zzorg) ;
00920 originRAIfv = THD_3dmm_to_dicomm( Opt->iset , originRAIfv ) ;
00921
00922 LOAD_FVEC3(fv2 , Opt->iset->daxes->xxorg + (Opt->iset->daxes->nxx-1)*Opt->iset->daxes->xxdel ,
00923 Opt->iset->daxes->yyorg + (Opt->iset->daxes->nyy-1)*Opt->iset->daxes->yydel ,
00924 Opt->iset->daxes->zzorg + (Opt->iset->daxes->nzz-1)*Opt->iset->daxes->zzdel ) ;
00925 fv2 = THD_3dmm_to_dicomm( Opt->iset , fv2 ) ;
00926
00927 if( originRAIfv.xyz[0] > fv2.xyz[0] ) { float tf; tf = originRAIfv.xyz[0]; originRAIfv.xyz[0] = fv2.xyz[0]; fv2.xyz[0] = tf; }
00928 if( originRAIfv.xyz[1] > fv2.xyz[1] ) { float tf; tf = originRAIfv.xyz[1]; originRAIfv.xyz[1] = fv2.xyz[1]; fv2.xyz[1] = tf; }
00929 if( originRAIfv.xyz[2] > fv2.xyz[2] ) { float tf; tf = originRAIfv.xyz[2]; originRAIfv.xyz[2] = fv2.xyz[2]; fv2.xyz[2] = tf; }
00930
00931 if (LocalHead) {
00932 fprintf(stderr,"++3dSpatNorm (ZSS): RAI origin info: %f %f %f\n", originRAIfv.xyz[0], originRAIfv.xyz[1], originRAIfv.xyz[2]);
00933 }
00934
00935 Opt->OrigSpatNormedSet = EDIT_empty_copy( NULL ) ;
00936 tross_Copy_History( Opt->iset , Opt->OrigSpatNormedSet ) ;
00937 tross_Make_History( "3dSkullStrip" , argc,argv , Opt->OrigSpatNormedSet ) ;
00938
00939 LOAD_IVEC3( nxyz , imout_orig->nx , imout_orig->ny , imout_orig->nz ) ;
00940 LOAD_FVEC3( dxyz , imout_orig->dx , imout_orig->dy , imout_orig->dz ) ;
00941 LOAD_FVEC3( orgxyz , originRAIfv.xyz[0] , originRAIfv.xyz[1] , originRAIfv.xyz[2] ) ;
00942 LOAD_IVEC3( orixyz , ORI_R2L_TYPE , ORI_A2P_TYPE , ORI_I2S_TYPE ) ;
00943
00944 prefix = SUMA_AfniPrefix(Opt->in_name, NULL, NULL, NULL);
00945 if (!prefix) { SUMA_SL_Err("Bad prefix!!!"); exit(1); }
00946 spatprefix = SUMA_append_string(prefix, "_SpatNorm_OrigSpace");
00947 EDIT_dset_items( Opt->OrigSpatNormedSet ,
00948 ADN_prefix , spatprefix ,
00949 ADN_datum_all , imout_orig->kind ,
00950 ADN_nxyz , nxyz ,
00951 ADN_xyzdel , dxyz ,
00952 ADN_xyzorg , orgxyz ,
00953 ADN_xyzorient , orixyz ,
00954 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00955 ADN_view_type , VIEW_ORIGINAL_TYPE ,
00956 ADN_type , HEAD_ANAT_TYPE ,
00957 ADN_func_type , ANAT_BUCK_TYPE ,
00958 ADN_none ) ;
00959
00960 EDIT_substitute_brick( Opt->OrigSpatNormedSet , 0 , imout_orig->kind , mri_data_pointer(imout_orig) ) ;
00961
00962 oset = r_new_resam_dset ( Opt->OrigSpatNormedSet, Opt->iset, 0, 0, 0, NULL, MRI_NN, NULL);
00963 if (!oset) {
00964 fprintf(stderr,"**ERROR: Failed to reslice!?\n"); exit(1);
00965 }
00966 DSET_delete(Opt->OrigSpatNormedSet); Opt->OrigSpatNormedSet = oset; oset = NULL;
00967
00968 if (Opt->WriteSpatNorm) {
00969 SUMA_LH("Writing SpatNormed dset in original space");
00970 DSET_write(Opt->OrigSpatNormedSet) ;
00971 }
00972
00973 if (prefix) SUMA_free(prefix); prefix = NULL;
00974 if (spatprefix) SUMA_free(spatprefix); spatprefix = NULL;
00975 }
00976
00977 if (imout_edge) {
00978 SUMA_SL_Note("Creating an output edge dataset ...");
00979 oset = EDIT_empty_copy( NULL ) ;
00980 tross_Copy_History( Opt->iset , oset ) ;
00981 tross_Make_History( "3dSkullStrip" , argc,argv , oset ) ;
00982
00983 LOAD_IVEC3( nxyz , imout_edge->nx , imout_edge->ny , imout_edge->nz ) ;
00984 LOAD_FVEC3( dxyz , imout_edge->dx , imout_edge->dy , imout_edge->dz ) ;
00985 LOAD_FVEC3( orgxyz , imout_edge->xo , imout_edge->yo , imout_edge->zo ) ;
00986 LOAD_IVEC3( orixyz , ORI_R2L_TYPE , ORI_A2P_TYPE , ORI_I2S_TYPE ) ;
00987
00988 prefix = SUMA_AfniPrefix(Opt->in_name, NULL, NULL, NULL);
00989 if (!prefix) { SUMA_SL_Err("Bad prefix!!!"); exit(1); }
00990 spatprefix = SUMA_append_string(prefix, "_EdgeSpatNorm");
00991 EDIT_dset_items( oset ,
00992 ADN_prefix , spatprefix ,
00993 ADN_datum_all , imout_edge->kind ,
00994 ADN_nxyz , nxyz ,
00995 ADN_xyzdel , dxyz ,
00996 ADN_xyzorg , orgxyz ,
00997 ADN_xyzorient , orixyz ,
00998 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00999 ADN_view_type , VIEW_ORIGINAL_TYPE ,
01000 ADN_type , HEAD_ANAT_TYPE ,
01001 ADN_func_type , ANAT_BUCK_TYPE ,
01002 ADN_none ) ;
01003
01004 EDIT_substitute_brick( oset , 0 , imout_edge->kind , mri_data_pointer(imout_edge) ) ;
01005 if (Opt->WriteSpatNorm) {
01006 SUMA_LH("Writing Edge dset");
01007 DSET_write(oset) ;
01008 }
01009 Opt->in_edvol = oset; oset = NULL;
01010
01011 if (prefix) SUMA_free(prefix); prefix = NULL;
01012 if (spatprefix) SUMA_free(spatprefix); spatprefix = NULL;
01013 fprintf(SUMA_STDERR,"%s: -spatnorm: Expecting %d voxels in in_vol dset (%d %d %d)\n",
01014 FuncName, DSET_NVOX( Opt->in_vol ), DSET_NX( Opt->in_vol ), DSET_NY( Opt->in_vol ), DSET_NZ( Opt->in_vol ));
01015 }
01016
01017 oset = EDIT_empty_copy( NULL ) ;
01018
01019 { char idstr[500], *nid=NULL; sprintf(idstr,"%s_Spatnorm", Opt->iset->idcode.str);
01020 SUMA_NEW_ID(nid, idstr); strncpy(oset->idcode.str, nid, IDCODE_LEN); SUMA_free(nid); nid = NULL;}
01021 tross_Copy_History( Opt->iset , oset ) ;
01022 tross_Make_History( "3dSkullStrip" , argc,argv , oset ) ;
01023
01024 LOAD_IVEC3( nxyz , imout->nx , imout->ny , imout->nz ) ;
01025 LOAD_FVEC3( dxyz , imout->dx , imout->dy , imout->dz ) ;
01026 LOAD_FVEC3( orgxyz , imout->xo , imout->yo , imout->zo ) ;
01027 LOAD_IVEC3( orixyz , ORI_R2L_TYPE , ORI_A2P_TYPE , ORI_I2S_TYPE ) ;
01028
01029 prefix = SUMA_AfniPrefix(Opt->in_name, NULL, NULL, NULL);
01030 if (!prefix) { SUMA_SL_Err("Bad prefix!!!"); exit(1); }
01031 spatprefix = SUMA_append_string(prefix, "_SpatNorm");
01032 EDIT_dset_items( oset ,
01033 ADN_prefix , spatprefix ,
01034 ADN_datum_all , imout->kind ,
01035 ADN_nxyz , nxyz ,
01036 ADN_xyzdel , dxyz ,
01037 ADN_xyzorg , orgxyz ,
01038 ADN_xyzorient , orixyz ,
01039 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
01040 ADN_view_type , VIEW_ORIGINAL_TYPE ,
01041 ADN_type , HEAD_ANAT_TYPE ,
01042 ADN_func_type , ANAT_BUCK_TYPE ,
01043 ADN_none ) ;
01044
01045 EDIT_substitute_brick( oset , 0 , imout->kind , mri_data_pointer(imout) ) ;
01046 if (Opt->WriteSpatNorm) {
01047 SUMA_LH("Writing SpatNormed dset");
01048 DSET_write(oset) ;
01049 }
01050 Opt->in_vol = oset; oset = NULL;
01051
01052 if (prefix) SUMA_free(prefix); prefix = NULL;
01053 if (spatprefix) SUMA_free(spatprefix); spatprefix = NULL;
01054 if( Opt->debug ) fprintf(SUMA_STDERR,"%s: -spatnorm: Expecting %d voxels in in_vol dset (%d %d %d)\n",
01055 FuncName, DSET_NVOX( Opt->in_vol ), DSET_NX( Opt->in_vol ), DSET_NY( Opt->in_vol ), DSET_NZ( Opt->in_vol ));
01056
01057 } else {
01058
01059 SUMA_SL_Note("Loading dset, performing no Spatial Normalization");
01060 Opt->in_vol = THD_open_dataset( Opt->in_name );
01061 if( !ISVALID_DSET(Opt->in_vol) ){
01062 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->in_name) ;
01063 exit(1);
01064 }
01065 if( Opt->debug ) fprintf(SUMA_STDERR,"%s: -no_spatnorm: Expecting %d voxels in in_vol dset (%d %d %d)\n",
01066 FuncName, DSET_NVOX( Opt->in_vol ), DSET_NX( Opt->in_vol ), DSET_NY( Opt->in_vol ), DSET_NZ( Opt->in_vol ));
01067
01068 DSET_load(Opt->in_vol);
01069 if (Opt->fillhole) Opt->OrigSpatNormedSet = Opt->in_vol;
01070
01071 mri_brainormalize_initialize(Opt->in_vol->daxes->xxdel, Opt->in_vol->daxes->yydel, Opt->in_vol->daxes->zzdel);
01072 if (DSET_NX( Opt->in_vol) != THD_BN_nx() || DSET_NY( Opt->in_vol) != THD_BN_ny() || DSET_NZ( Opt->in_vol) != THD_BN_nz() ) {
01073 fprintf(SUMA_STDERR,"Error %s:\n SpatNormed Dset must be %d x %d x %d\n", FuncName, THD_BN_nx(), THD_BN_ny(), THD_BN_nz() );
01074 exit(1);
01075 }
01076 Opt->iset_hand = SUMA_THD_handedness( Opt->in_vol );
01077 if (LocalHead) fprintf(SUMA_STDERR,"%s: Handedness of orig dset %d\n", FuncName, Opt->iset_hand);
01078
01079 }
01080
01081
01082 if (!ISVALID_DSET(Opt->in_vol)) {
01083 if (!Opt->in_name) {
01084 SUMA_SL_Err("NULL input volume.");
01085 exit(1);
01086 } else {
01087 SUMA_SL_Err("invalid volume.");
01088 exit(1);
01089 }
01090 } else if ( DSET_BRICK_TYPE(Opt->in_vol, 0) == MRI_complex) {
01091 SUMA_SL_Err("Can't do complex data.");
01092 exit(1);
01093 }
01094
01095 if (Opt->blur_fwhm) {
01096 SUMA_SL_Note("Blurring...");
01097 EDIT_blur_volume( DSET_NX(Opt->in_vol), DSET_NY(Opt->in_vol), DSET_NZ(Opt->in_vol),
01098 SUMA_ABS((DSET_DX(Opt->in_vol))), SUMA_ABS((DSET_DY(Opt->in_vol))), SUMA_ABS((DSET_DZ(Opt->in_vol))),
01099 DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol,0), 0.42466090*Opt->blur_fwhm) ;
01100 }
01101
01102 if (Opt->UseSkull) { SOhull = SUMA_Alloc_SurfObject_Struct(1); }
01103 else SOhull = NULL;
01104 sprintf(stmphull,"OuterHull");
01105
01106 vol = SUMA_LoadPrepInVol (Opt, &SOhull);
01107 if ( vol <= 0 ) {
01108 SUMA_S_Err("Failed to load/prep volume");
01109 exit(1);
01110 } else {
01111 if (LocalHead) fprintf (SUMA_STDERR,"%s: Got me a volume of %f mm3\n", FuncName, vol);
01112 }
01113
01114
01115
01116 do {
01117
01118 sprintf(stmp,"icobaby_ld%d", Opt->Icold);
01119 SO = SUMA_CreateIcosahedron (Opt->r/2.0, Opt->Icold, Opt->cog, "n", 1);
01120 if (!SO) {
01121 SUMA_S_Err("Failed to create Icosahedron");
01122 exit(1);
01123 }
01124
01125 if (Opt->Stop) SUMA_free(Opt->Stop); Opt->Stop = NULL;
01126 if (!Opt->Stop) {
01127 Opt->Stop = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
01128 if (!Opt->Stop) {
01129 SUMA_SL_Crit("Failed to allocate");
01130 exit(1);
01131 }
01132 }
01133
01134 if (Opt->avoid_vent) {
01135 float U[3], Un, *a, P2[2][3];
01136 for (i=0; i<SO->N_Node; ++i) {
01137
01138 a = &(SO->NodeList[3*i]);
01139 if (a[2] - SO->Center[2] > Opt->r/3 || a[1] - SO->Center[1] > Opt->r/2) {
01140 SUMA_UNIT_VEC(SO->Center, a, U, Un);
01141 SUMA_POINT_AT_DISTANCE_NORM(U, SO->Center, (Un+1.1*Opt->d1), P2);
01142 SO->NodeList[3*i] = P2[0][0]; SO->NodeList[3*i+1] = P2[0][1]; SO->NodeList[3*i+2] = P2[0][2];
01143 }
01144 }
01145 }
01146
01147 Opt->ztv = (float *)SUMA_malloc(sizeof(float)*SO->N_Node);
01148 if (!Opt->ztv) {
01149 SUMA_SL_Crit("Failed to allocate");
01150 exit(1);
01151 }
01152 for (i=0; i<SO->N_Node; ++i) Opt->ztv[i] = Opt->Zt;
01153
01154
01155 SO->VolPar = SUMA_VolParFromDset (Opt->in_vol);
01156 SO->SUMA_VolPar_Aligned = YUP;
01157 if (0){
01158 SUMA_VTI *vti=NULL;
01159 int iii, n1, n2, n3;
01160 float *tmpXYZ=NULL;
01161 SUMA_S_Note("Test Section");
01162
01163 tmpXYZ = (float *)SUMA_malloc(SO->N_Node * 3 * sizeof(float));
01164 if (!tmpXYZ) {
01165 SUMA_SL_Crit("Faile to allocate");
01166 exit(1);
01167 }
01168 memcpy ((void*)tmpXYZ, (void *)SO->NodeList, SO->N_Node * 3 * sizeof(float));
01169
01170
01171 SUMA_vec_dicomm_to_3dfind (tmpXYZ, SO->N_Node, SO->VolPar);
01172 vti = SUMA_CreateVTI(SO->N_FaceSet, NULL);
01173 for (iii=0; iii<SO->N_FaceSet; ++iii) vti->TriIndex[iii] = iii;
01174 n1 = SO->FaceSetList[5*3]; n2 = SO->FaceSetList[5*3+1]; n3 = SO->FaceSetList[5*3+2];
01175 fprintf(SUMA_STDERR, "%s: Triangle 5 has nodes: \n"
01176 " %f %f %f\n"
01177 " %f %f %f\n"
01178 " %f %f %f\n", FuncName,
01179 SO->NodeList[n1*3], SO->NodeList[n1*3+1], SO->NodeList[n1*3+2],
01180 SO->NodeList[n2*3], SO->NodeList[n2*3+1], SO->NodeList[n2*3+2],
01181 SO->NodeList[n3*3], SO->NodeList[n3*3+1], SO->NodeList[n3*3+2]);
01182
01183 vti = SUMA_GetVoxelsIntersectingTriangle(SO, SO->VolPar, tmpXYZ, vti);
01184 vti = SUMA_FreeVTI(vti);
01185 }
01186 if (!SO->State) {SO->State = SUMA_copy_string("3dSkullStrip"); }
01187 if (!SO->Group) {SO->Group = SUMA_copy_string("3dSkullStrip"); }
01188 if (!SO->Label) {SO->Label = SUMA_copy_string(stmp); }
01189
01190
01191 if (SO->Label) { if (SO->idcode_str) SUMA_free(SO->idcode_str); SO->idcode_str = NULL; SUMA_NEW_ID(SO->idcode_str, SO->Label); }
01192
01193 if (!nint && SOhull) {
01194 SOhull->VolPar = SUMA_VolParFromDset (Opt->in_vol);
01195 SOhull->SUMA_VolPar_Aligned = YUP;
01196
01197 if (!SOhull->State) {SOhull->State = SUMA_copy_string("3dSkullStrip"); }
01198 if (!SOhull->Group) {SOhull->Group = SUMA_copy_string("3dSkullStrip"); }
01199 if (!SOhull->Label) {SOhull->Label = SUMA_copy_string(stmphull); }
01200 if (SOhull->Label) { if (SOhull->idcode_str) SUMA_free(SOhull->idcode_str); SOhull->idcode_str = NULL; SUMA_NEW_ID(SOhull->idcode_str,SOhull->Label); }
01201 }
01202
01203
01204 if (ps->cs->talk_suma) {
01205 ps->cs->istream = SUMA_BRAINWRAP_LINE;
01206 ps->cs->afni_istream = SUMA_AFNI_STREAM_INDEX2;
01207 ps->cs->kth = 1;
01208 if (!SUMA_SendToSuma (NULL, ps->cs, NULL, SUMA_NO_DSET_TYPE, 0)) {
01209 SUMA_SL_Err("Failed to initialize SUMA_SendToSuma");
01210 ps->cs->Send = NOPE;
01211 ps->cs->afni_Send = NOPE;
01212 ps->cs->talk_suma = NOPE;
01213 } else if (!SUMA_SendToAfni (ps->cs, NULL, 0)) {
01214 SUMA_SL_Err("Failed to initialize SUMA_SendToAfni");
01215 ps->cs->afni_Send = NOPE;
01216 ps->cs->Send = NOPE;
01217 } else {
01218 if (!nint && SOhull) {
01219 if (Opt->send_hull) {
01220 SUMA_LH("Sending Hull");
01221 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("Sending HULL next"); }
01222 SUMA_SendSumaNewSurface(SOhull, ps->cs);
01223 }
01224 }
01225 SUMA_LH("Sending Ico");
01226 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("Sending Ico next"); }
01227 SUMA_SendSumaNewSurface(SO, ps->cs);
01228
01229 }
01230 ps->cs->kth = kth_buf;
01231 }
01232
01233 if (!nint && Opt->UseSkull) {
01234
01235 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("Shrinking skull hull next"); }
01236 SUMA_SkullMask (SOhull, Opt, ps->cs);
01237
01238 fprintf (SUMA_STDERR,"%s: Locating voxels on skull boundary ...\n", FuncName);
01239 isin = SUMA_FindVoxelsInSurface (SOhull, SO->VolPar, &N_in, 0, NULL);
01240 for (i=0; i<SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz; ++i) { if (isin[i] <= SUMA_ON_NODE) Opt->dvec[i] = 0; }
01241 #if 0
01242 SUMA_SL_Note("Writing hull mask");
01243 {
01244 FILE *fout=fopen("hullmask.1D","w");
01245 int ii, jj, kk;
01246 for (i=0; i<SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz; ++i) {
01247 SUMA_1D_2_3D_index(i, ii, jj, kk, SO->VolPar->nx, SO->VolPar->nx*SO->VolPar->ny);
01248 fprintf(fout,"%d %d %d %d\n",ii, jj, kk, isin[i]);
01249 }
01250 fclose(fout);
01251 }
01252 #endif
01253 if (isin) SUMA_free(isin); isin = NULL;
01254 }
01255
01256
01257 if (Opt->DoSpatNorm && ps->cs->afni_Send) {
01258 SUMA_SL_Note("Sending spatnormed volume to AFNI");
01259 if (!SUMA_SendToAfni(ps->cs, Opt->in_vol, 1)) {
01260 SUMA_SL_Err("Failed to send volume to AFNI");
01261 ps->cs->afni_Send = NOPE;
01262 }
01263 }
01264
01265
01266 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("Brain expansion next"); }
01267 SUMA_StretchToFitLeCerveau (SO, Opt, ps->cs);
01268
01269
01270 if (Opt->PercInt >= 0) {
01271 if (Opt->debug) fprintf(SUMA_STDERR,"%s: Checking for self intersection...\n", FuncName);
01272 nseg = 30 * Opt->Icold * Opt->Icold;
01273 nint = SUMA_isSelfIntersect(SO, (int)(Opt->PercInt * nseg / 100.0));
01274 if (nint < 0) {
01275 SUMA_SL_Err("Error in SUMA_isSelfIntersect. Ignoring self intersection test.");
01276 nint = 0;
01277 }
01278
01279 pint = (float)nint / (float)nseg * 100;
01280 if (LocalHead) fprintf (SUMA_STDERR,"%s: at least %d (%f%%) of segments intersect the surface.\nStopping criterion is %f%%\n", FuncName, nint, pint, Opt->PercInt);
01281 if (pint > Opt->PercInt) {
01282 static int SelfIntRep = 0;
01283 int smstep;
01284 if (SelfIntRep < Opt->MaxIntIter) {
01285
01286 if (!Opt->avoid_vent && !SelfIntRep) {
01287 Opt->avoid_vent = 1;
01288 } else {
01289 smstep = 12 * ( Opt->Icold / 25 ) * ( Opt->Icold / 25 );
01290 if ( smstep < 12) smstep = 12;
01291 else if ( smstep % 2 ) ++smstep;
01292 #if 0
01293 if (!(SelfIntRep % 2)) {
01294 Opt->avoid_vent = 0;
01295 Opt->NNsmooth += smstep;
01296 } else {
01297 Opt->avoid_vent = 1;
01298 }
01299 #else
01300
01301 Opt->NNsmooth += smstep;
01302 #endif
01303 ++SelfIntRep;
01304 }
01305 fprintf(SUMA_STDERR,"Warning %s:****************\n Surface self intersecting! trying again:\n smoothing of %d, avoid_vent %d\n", FuncName, Opt->NNsmooth, (int)Opt->avoid_vent);
01306 SUMA_Free_Surface_Object(SO); SO = NULL;
01307 } else {
01308 fprintf(SUMA_STDERR,"Warning %s: Stubborn intersection remaining at smoothing of %d. Ignoring it.", FuncName, Opt->NNsmooth);
01309 nint = 0;
01310 }
01311 } else {
01312 if (Opt->debug) {
01313 if (nint) fprintf(SUMA_STDERR,"%s: Number of intersection below criterion.\n", FuncName);
01314 else fprintf(SUMA_STDERR,"%s: No intersections found.\n", FuncName);
01315 }
01316 nint = 0;
01317 }
01318 } else {
01319 fprintf(SUMA_STDERR,"%s: Self intersection check turned off.\n", FuncName);
01320 nint = 0;
01321 }
01322 } while (nint != 0);
01323 if (Opt->debug) fprintf(SUMA_STDERR,"%s: Final smoothing of %d\n", FuncName, Opt->NNsmooth);
01324 if (SUMA_DidUserQuit()) {
01325 if (Opt->debug) fprintf(SUMA_STDERR,"%s: straight to end per user demand (%d)...\n", FuncName, SUMA_DidUserQuit());
01326 goto FINISH_UP;
01327 }
01328
01329
01330 if (Opt->UseNew) {
01331 double mval = 255;
01332 if (Opt->debug) fprintf (SUMA_STDERR,"%s: Touchup correction, pass 1 ...\n", FuncName);
01333 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("touchup correction next"); }
01334 if (Opt->DemoPause == SUMA_3dSS_INTERACTIVE) {
01335 fprintf (SUMA_STDERR,"3dSkullStrip Interactive: \n"
01336 "Touchup, pass 1.\n"
01337 "Do you want to (C)ontinue, (P)ass or (S)ave this? [C] ");
01338 cbuf = SUMA_ReadCharStdin ('c', 0,"csp");
01339 switch (cbuf) {
01340 case 's':
01341 fprintf (SUMA_STDERR,"Saving mask as is.\n");
01342 goto FINISH_UP;
01343 break;
01344 case 'p':
01345 fprintf (SUMA_STDERR,"Passing this stage \n");
01346 goto BEAUTY;
01347 break;
01348 case 'c':
01349 fprintf (SUMA_STDERR,"Continuing with stage.\n");
01350 break;
01351 }
01352 }
01353
01354 if (mval < Opt->t98) {
01355 SUMA_SL_Warn("Looks like some values in dset might be larger than 255 !");
01356 mval = Opt->t98+10;
01357 }
01358 if (Opt->k98maskcnt && Opt->k98mask) { for (ii=0; ii<Opt->k98maskcnt; ++ii) Opt->dvec[Opt->k98mask[ii]] = mval; }
01359
01360 SUMA_Reposition_Touchup(SO, Opt, 6, ps->cs);
01361
01362 if (LocalHead) fprintf (SUMA_STDERR,"%s: Touchup correction Done.\n", FuncName);
01363 }
01364 BEAUTY:
01365
01366 if (Opt->smooth_end) {
01367 ps->cs->kth = 1;
01368 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("beauty treatment smoothing next"); }
01369 if (Opt->debug) fprintf (SUMA_STDERR,"%s: The beauty treatment smoothing.\n", FuncName);
01370 if (Opt->DemoPause == SUMA_3dSS_INTERACTIVE) {
01371 fprintf (SUMA_STDERR,"3dSkullStrip Interactive: \n"
01372 "Beauty treatment smoothing.\n"
01373 "Do you want to (C)ontinue, (P)ass or (S)ave this? [C] ");
01374 cbuf = SUMA_ReadCharStdin ('c', 0,"csp");
01375 switch (cbuf) {
01376 case 's':
01377 fprintf (SUMA_STDERR,"Saving mask as is.\n");
01378 goto FINISH_UP;
01379 break;
01380 case 'p':
01381 fprintf (SUMA_STDERR,"Passing this stage \n");
01382 goto TOUCHUP_2;
01383 break;
01384 case 'c':
01385 fprintf (SUMA_STDERR,"Continuing with stage.\n");
01386 break;
01387 }
01388 }
01389 dsmooth = SUMA_Taubin_Smooth( SO, NULL,
01390 0.6307, -.6732, SO->NodeList,
01391 Opt->smooth_end, 3, SUMA_ROW_MAJOR, dsmooth, ps->cs, NULL);
01392 memcpy((void*)SO->NodeList, (void *)dsmooth, SO->N_Node * 3 * sizeof(float));
01393 SUMA_RECOMPUTE_NORMALS(SO);
01394 if (ps->cs->Send) {
01395 if (!SUMA_SendToSuma (SO, ps->cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
01396 SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
01397 }
01398 }
01399 ps->cs->kth = kth_buf;
01400 }
01401 TOUCHUP_2:
01402
01403 if (Opt->UseNew) {
01404 ps->cs->kth = 1;
01405 if (Opt->DemoPause == SUMA_3dSS_DEMO_PAUSE) { SUMA_PAUSE_PROMPT("touchup correction 2 next"); }
01406 if (Opt->debug) fprintf (SUMA_STDERR,"%s: Final touchup correction ...\n", FuncName);
01407 if (Opt->DemoPause == SUMA_3dSS_INTERACTIVE) {
01408 fprintf (SUMA_STDERR,"3dSkullStrip Interactive: \n"
01409 "Touchup, pass 2.\n"
01410 "Do you want to (C)ontinue, (P)ass or (S)ave this? [C] ");
01411 cbuf = SUMA_ReadCharStdin ('c', 0,"csp");
01412 fprintf (SUMA_STDERR,"%c\n", cbuf);
01413 switch (cbuf) {
01414 case 's':
01415 fprintf (SUMA_STDERR,"Saving mask as is.\n");
01416 goto FINISH_UP;
01417 break;
01418 case 'p':
01419 fprintf (SUMA_STDERR,"Passing this stage \n");
01420 goto FINISH_UP;
01421 break;
01422 case 'c':
01423 fprintf (SUMA_STDERR,"Continuing with stage.\n");
01424 break;
01425 }
01426 }
01427
01428 SUMA_Reposition_Touchup(SO, Opt, 2, ps->cs);
01429
01430 if (LocalHead) fprintf (SUMA_STDERR,"%s: Final touchup correction Done.\n", FuncName);
01431 ps->cs->kth = kth_buf;
01432 }
01433
01434 FINISH_UP:
01435
01436 ps->cs->kth = 1;
01437 if (ps->cs->Send) {
01438 if (!SUMA_SendToSuma (SO, ps->cs, (void *)SO->NodeList, SUMA_NODE_XYZ, 1)) {
01439 SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCommunication halted.");
01440 }
01441 }
01442 ps->cs->kth = kth_buf;
01443
01444 if (Opt->DoSpatNorm) {
01445 float ispat, jspat, kspat, iorig, jorig, korig , xrai_orig, yrai_orig, zrai_orig;
01446
01447 SUMA_LH("Changing coords of output surface");
01448
01449
01450 ispat = 0; jspat = 0; kspat = 0;
01451 brainnormalize_coord( ispat , jspat, kspat,
01452 &iorig, &jorig, &korig , Opt->iset,
01453 &xrai_orig, &yrai_orig, &zrai_orig);
01454 Opt->SpatShift[0] = xrai_orig - THD_BN_XORG; Opt->SpatShift[1] = yrai_orig - THD_BN_YORG; Opt->SpatShift[2] = zrai_orig - THD_BN_ZORG;
01455
01456
01457 for (i=0; i<SO->N_Node; ++i) {
01458 i3 = 3*i;
01459 SO->NodeList[i3 + 0] += Opt->SpatShift[0];
01460 SO->NodeList[i3 + 1] += Opt->SpatShift[1];
01461 SO->NodeList[i3 + 2] += Opt->SpatShift[2];
01462 }
01463
01464 if (SOhull) {
01465 for (i=0; i<SOhull->N_Node; ++i) {
01466 i3 = 3*i;
01467 SOhull->NodeList[i3 + 0] += Opt->SpatShift[0];
01468 SOhull->NodeList[i3 + 1] += Opt->SpatShift[1];
01469 SOhull->NodeList[i3 + 2] += Opt->SpatShift[2];
01470 }
01471 }
01472
01473
01474 if (LocalHead) fprintf(SUMA_STDERR,"%s: \nPre: %d %d %d\n", FuncName, SO->VolPar->nx, SO->VolPar->ny, SO->VolPar->nz);
01475 SUMA_Free_VolPar(SO->VolPar); SO->VolPar = NULL;
01476 SO->VolPar = SUMA_VolParFromDset(Opt->iset);
01477 if (LocalHead) fprintf(SUMA_STDERR,"%s: \nPost: %d %d %d\n", FuncName, SO->VolPar->nx, SO->VolPar->ny, SO->VolPar->nz);
01478 if (SOhull) {
01479 SUMA_Free_VolPar(SOhull->VolPar); SOhull->VolPar = NULL;
01480 SOhull->VolPar = SUMA_VolParFromDset (Opt->iset);
01481 }
01482 }
01483
01484
01485 if (strcmp(Opt->out_prefix,"skull_strip_out")) {
01486 fprintf (SUMA_STDERR,"%s: Writing surface ...\n", FuncName);
01487 if (!SUMA_Save_Surface_Object (SO_name, SO, Opt->SurfFileType, Opt->SurfFileFormat, NULL)) {
01488 fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
01489 exit (1);
01490 }
01491 }
01492
01493 if (Opt->UseSkull && SOhull) {
01494 fprintf (SUMA_STDERR,"%s: Writing skull surface ...\n", FuncName);
01495 if (!SUMA_Save_Surface_Object (SO_name_hull, SOhull, Opt->SurfFileType, Opt->SurfFileFormat, NULL)) {
01496 fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
01497 exit (1);
01498 }
01499 }
01500
01501
01502 OptDs = SUMA_New_FormAfniDset_Opt();
01503 if (Opt->out_vol_prefix) {
01504 SUMA_FileName NewName = SUMA_StripPath(Opt->out_vol_prefix);
01505 OptDs->prefix = NewName.FileName; NewName.FileName = NULL;
01506 OptDs->prefix_path = NewName.Path; NewName.Path = NULL;
01507 } else {
01508 OptDs->prefix = SUMA_copy_string("3dSkullStrip");
01509 OptDs->prefix_path = SUMA_copy_string("./");
01510 }
01511 if (Opt->DoSpatNorm) {
01512 OptDs->mset = Opt->iset;
01513 } else {
01514 OptDs->mset = Opt->in_vol;
01515 }
01516
01517
01518 if (Opt->debug) fprintf (SUMA_STDERR,"%s: Locating voxels inside surface ...\n", FuncName);
01519 isin = SUMA_FindVoxelsInSurface (SO, SO->VolPar, &N_in, Opt->fillhole, Opt->OrigSpatNormedSet);
01520 isin_float = (float *)SUMA_malloc(sizeof(float) * SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz);
01521 if (!isin_float) {
01522 SUMA_SL_Crit("Failed to allocate");
01523 exit(1);
01524 }
01525
01526
01527 if (Opt->dvec) { SUMA_free(Opt->dvec); Opt->dvec = NULL; }
01528 Opt->nvox = DSET_NVOX( Opt->OrigSpatNormedSet );
01529 Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
01530 if (!Opt->dvec) {
01531 SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
01532 SUMA_RETURN(NOPE);
01533 }
01534 EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->OrigSpatNormedSet,0) ,
01535 DSET_BRICK_TYPE(Opt->OrigSpatNormedSet,0), DSET_ARRAY(Opt->OrigSpatNormedSet, 0) ,
01536 MRI_double , Opt->dvec ) ;
01537
01538 if (!Opt->MaskMode) {
01539 SUMA_LH("Creating skull-stripped volume");
01540 for (i=0; i<SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz; ++i) {
01541
01542 if (isin[i] >= SUMA_ON_NODE) isin_float[i] = (float)Opt->dvec[i];
01543 else isin_float[i] = 0.0;
01544 }
01545 } else {
01546 SUMA_LH("Creating mask volume");
01547 for (i=0; i<SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz; ++i) {
01548 isin_float[i] = (float)isin[i];
01549 }
01550 }
01551 if (isin) SUMA_free(isin); isin = NULL;
01552
01553 if (Opt->debug) fprintf (SUMA_STDERR,"%s: Writing masked volume ...\n", FuncName);
01554 OptDs->full_list = 1;
01555 OptDs->dval = 1;
01556 dset = SUMA_FormAfnidset (NULL, isin_float, SO->VolPar->nx*SO->VolPar->ny*SO->VolPar->nz, OptDs);
01557
01558
01559 if (1){
01560 EDIT_options *edopt = SUMA_BlankAfniEditOptions();
01561 if (Opt->debug) fprintf (SUMA_STDERR,"%s: Applying a bit of erosion and dilatation \n", FuncName);
01562 edopt->edit_clust = ECFLAG_SAME;
01563 edopt->clust_rmm = SUMA_MAX_PAIR(SUMA_ABS((DSET_DX(OptDs->mset))), SUMA_ABS((DSET_DY(OptDs->mset))));
01564 edopt->clust_rmm = SUMA_MAX_PAIR(SUMA_ABS((DSET_DZ(OptDs->mset))), edopt->clust_rmm)*1.01;
01565 edopt->clust_vmul = 1000*SUMA_ABS((DSET_DX(OptDs->mset)))*SUMA_ABS((DSET_DY(OptDs->mset)))*SUMA_ABS((DSET_DZ(OptDs->mset)));
01566 edopt->erode_pv = 75.0 / 100.0;
01567 edopt->dilate = 1;
01568 EDIT_one_dataset( dset , edopt);
01569 SUMA_free(edopt);
01570 }
01571
01572 if (!dset) {
01573 SUMA_SL_Err("Failed to create output dataset!");
01574 } else {
01575 tross_Make_History( FuncName , argc,argv , dset ) ;
01576 DSET_write(dset) ;
01577 }
01578
01579
01580 if (ps->cs->Send && !ps->cs->GoneBad) {
01581
01582 if (!SUMA_SendToSuma (SO, ps->cs, NULL, SUMA_NODE_XYZ, 2)) {
01583 SUMA_SL_Warn("Failed in SUMA_SendToSuma\nCleanup failed");
01584 }
01585 }
01586
01587 if (dsmooth) SUMA_free(dsmooth); dsmooth = NULL;
01588 if (OptDs) { OptDs->mset = NULL; OptDs = SUMA_Free_FormAfniDset_Opt(OptDs); }
01589 if (dset) { DSET_delete(dset); dset = NULL; }
01590 if (isin) { SUMA_free(isin); isin = NULL; }
01591 if (isin_float) { SUMA_free(isin_float); isin_float = NULL; }
01592 if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
01593 if (Opt) Opt = SUMA_Free_Generic_Prog_Options_Struct(Opt);
01594 if (hullprefix) SUMA_free(hullprefix); hullprefix = NULL;
01595 if (SO_name_hull) SUMA_free(SO_name_hull); SO_name_hull = NULL;
01596 if (SO_name) SUMA_free(SO_name); SO_name = NULL;
01597 if (SO) SUMA_Free_Surface_Object(SO); SO = NULL;
01598 if (SOhull) SUMA_Free_Surface_Object(SOhull); SOhull = NULL;
01599 if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) {
01600 SUMA_SL_Err("DO Cleanup Failed!");
01601 }
01602 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
01603 exit(0);
01604 }
01605 #endif
