00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "SUMA_suma.h"
00012 #include "MarchingCubes/MarchingCubes.h"
00013
00014 #undef STAND_ALONE
00015
00016 #if defined SUMA_IsoSurface_STANDALONE
00017 #define STAND_ALONE
00018 #endif
00019
00020 #ifdef STAND_ALONE
00021
00022 SUMA_SurfaceViewer *SUMAg_cSV = NULL;
00023 SUMA_SurfaceViewer *SUMAg_SVv = NULL;
00024
00025 int SUMAg_N_SVv = 0;
00026 SUMA_DO *SUMAg_DOv = NULL;
00027 int SUMAg_N_DOv = 0;
00028 SUMA_CommonFields *SUMAg_CF = NULL;
00029 #else
00030 extern SUMA_CommonFields *SUMAg_CF;
00031 extern SUMA_DO *SUMAg_DOv;
00032 extern SUMA_SurfaceViewer *SUMAg_SVv;
00033 extern int SUMAg_N_SVv;
00034 extern int SUMAg_N_DOv;
00035 #endif
00036
00037
00038 #define SUMA_SET_MC_DATA(mcb, val, i, j, k) { \
00039 mcb->data[ i + j*mcb->size_x + k*mcb->size_x*mcb->size_y] = val; \
00040 }
00041
00042
00043
00044
00045
00046
00047 SUMA_SurfaceObject *SUMA_MarchingCubesSurface(SUMA_GENERIC_PROG_OPTIONS_STRUCT * Opt)
00048 {
00049 static char FuncName[]={"SUMA_MarchingCubesSurface"};
00050 SUMA_SurfaceObject *SO=NULL;
00051 int nxx, nyy, nzz, cnt, i, j, k, *FaceSetList=NULL;
00052 float *NodeList=NULL;
00053 SUMA_NEW_SO_OPT *nsoopt = NULL;
00054 THD_fvec3 fv, iv;
00055 MCB *mcp ;
00056
00057 SUMA_ENTRY;
00058
00059 if (Opt->obj_type < 0) {
00060 nxx = DSET_NX(Opt->in_vol);
00061 nyy = DSET_NY(Opt->in_vol);
00062 nzz = DSET_NZ(Opt->in_vol);
00063
00064 if (Opt->debug) {
00065 fprintf(SUMA_STDERR,"%s:\nNxx=%d\tNyy=%d\tNzz=%d\n", FuncName, nxx, nyy, nzz);
00066 }
00067
00068 mcp = MarchingCubes(-1, -1, -1);
00069 set_resolution( mcp, nxx, nyy, nzz ) ;
00070 init_all(mcp) ;
00071 if (Opt->debug) fprintf(SUMA_STDERR,"%s:\nSetting data...\n", FuncName);
00072 cnt = 0;
00073 for( k = 0 ; k < mcp->size_z ; k++ ) {
00074 for( j = 0 ; j < mcp->size_y ; j++ ) {
00075 for( i = 0 ; i < mcp->size_x ; i++ ) {
00076 SUMA_SET_MC_DATA ( mcp, Opt->mcdatav[cnt], i, j, k);
00077 ++cnt;
00078 }
00079 }
00080 }
00081
00082 } else {
00083
00084 nxx = nyy = nzz = Opt->obj_type_res;
00085 mcp = MarchingCubes(-1, -1, -1);
00086 set_resolution( mcp, nxx, nyy, nzz) ;
00087 init_all(mcp) ;
00088 compute_data( *mcp , Opt->obj_type) ;
00089 }
00090
00091
00092 if (Opt->debug) fprintf(SUMA_STDERR,"%s:\nrunning MarchingCubes...\n", FuncName);
00093 run(mcp) ;
00094 clean_temps(mcp) ;
00095
00096 if (Opt->debug > 1) {
00097 fprintf(SUMA_STDERR,"%s:\nwriting out NodeList and FaceSetList...\n", FuncName);
00098 write1Dmcb(mcp);
00099 }
00100
00101 if (Opt->debug) {
00102 fprintf(SUMA_STDERR,"%s:\nNow creating SO...\n", FuncName);
00103 }
00104
00105 NodeList = (float *)SUMA_malloc(sizeof(float)*3*mcp->nverts);
00106 FaceSetList = (int *)SUMA_malloc(sizeof(int)*3*mcp->ntrigs);
00107 if (!NodeList || !FaceSetList) {
00108 SUMA_SL_Crit("Failed to allocate!");
00109 SUMA_RETURN(SO);
00110 }
00111
00112 if (Opt->obj_type < 0) {
00113 if (Opt->debug) {
00114 fprintf(SUMA_STDERR,"%s:\nCopying vertices, changing to DICOM \nOrig:(%f %f %f) \nD:(%f %f %f)...\n",
00115 FuncName, DSET_XORG(Opt->in_vol), DSET_YORG(Opt->in_vol), DSET_ZORG(Opt->in_vol),
00116 DSET_DX(Opt->in_vol), DSET_DY(Opt->in_vol), DSET_DZ(Opt->in_vol));
00117 }
00118 for ( i = 0; i < mcp->nverts; i++ ) {
00119 j = 3*i;
00120 fv.xyz[0] = DSET_XORG(Opt->in_vol) + mcp->vertices[i].x * DSET_DX(Opt->in_vol);
00121 fv.xyz[1] = DSET_YORG(Opt->in_vol) + mcp->vertices[i].y * DSET_DY(Opt->in_vol);
00122 fv.xyz[2] = DSET_ZORG(Opt->in_vol) + mcp->vertices[i].z * DSET_DZ(Opt->in_vol);
00123
00124 iv = SUMA_THD_3dmm_to_dicomm( Opt->in_vol->daxes->xxorient, Opt->in_vol->daxes->yyorient, Opt->in_vol->daxes->zzorient, fv );
00125 NodeList[j ] = iv.xyz[0];
00126 NodeList[j+1] = iv.xyz[1];
00127 NodeList[j+2] = iv.xyz[2];
00128 }
00129 for ( i = 0; i < mcp->ntrigs; i++ ) {
00130 j = 3*i;
00131 FaceSetList[j ] = mcp->triangles[i].v3;
00132 FaceSetList[j+1] = mcp->triangles[i].v2;
00133 FaceSetList[j+2] = mcp->triangles[i].v1;
00134 }
00135 } else {
00136
00137 for ( i = 0; i < mcp->nverts; i++ ) {
00138 j = 3*i;
00139 NodeList[j ] = mcp->vertices[i].x;
00140 NodeList[j+1] = mcp->vertices[i].y;
00141 NodeList[j+2] = mcp->vertices[i].z;
00142 }
00143 for ( i = 0; i < mcp->ntrigs; i++ ) {
00144 j = 3*i;
00145 FaceSetList[j ] = mcp->triangles[i].v3;
00146 FaceSetList[j+1] = mcp->triangles[i].v2;
00147 FaceSetList[j+2] = mcp->triangles[i].v1;
00148 }
00149 }
00150
00151
00152 nsoopt = SUMA_NewNewSOOpt();
00153 SO = SUMA_NewSO(&NodeList, mcp->nverts, &FaceSetList, mcp->ntrigs, nsoopt);
00154
00155 if (Opt->debug) {
00156 fprintf(SUMA_STDERR,"%s:\nCleaning mcp...\n", FuncName);
00157 }
00158 clean_all(mcp) ;
00159 free(mcp);
00160 nsoopt=SUMA_FreeNewSOOpt(nsoopt);
00161 SUMA_RETURN(SO);
00162 }
00163
00164
00165
00166 SUMA_Boolean SUMA_Get_isosurface_datasets (SUMA_GENERIC_PROG_OPTIONS_STRUCT * Opt)
00167 {
00168 static char FuncName[]={"SUMA_Get_isosurface_datasets"};
00169 int i;
00170
00171 SUMA_ENTRY;
00172
00173 Opt->in_vol = THD_open_dataset( Opt->in_name );
00174 if (!ISVALID_DSET(Opt->in_vol)) {
00175 if (!Opt->in_name) {
00176 SUMA_SL_Err("NULL input volume.");
00177 SUMA_RETURN(NOPE);
00178 } else {
00179 SUMA_SL_Err("invalid volume.");
00180 SUMA_RETURN(NOPE);
00181 }
00182 } else if ( DSET_BRICK_TYPE(Opt->in_vol, 0) == MRI_complex) {
00183 SUMA_SL_Err("Can't do complex data.");
00184 SUMA_RETURN(NOPE);
00185 }
00186
00187 Opt->nvox = DSET_NVOX( Opt->in_vol );
00188 if (DSET_NVALS( Opt->in_vol) != 1) {
00189 SUMA_SL_Err("Input volume can only have one sub-brick in it.\nUse [.] selectors to choose sub-brick needed.");
00190 SUMA_RETURN(NOPE);
00191 }
00192
00193
00194 Opt->mcdatav = (double *)SUMA_malloc(sizeof(double)*Opt->nvox);
00195 if (!Opt->mcdatav) {
00196 SUMA_SL_Crit("Failed to allocate for maskv");
00197 SUMA_RETURN(NOPE);
00198 }
00199 if (Opt->xform == SUMA_ISO_XFORM_MASK) {
00200 switch (Opt->MaskMode) {
00201 case SUMA_ISO_CMASK:
00202 if (Opt->cmask) {
00203
00204 int clen = strlen( Opt->cmask );
00205 char * cmd;
00206 byte *bmask;
00207
00208
00209 cmd = (char *)malloc((clen + 1) * sizeof(char));
00210 strcpy( cmd, Opt->cmask);
00211
00212 bmask = EDT_calcmask( cmd, &Opt->ninmask );
00213
00214 free( cmd );
00215
00216 if ( bmask == NULL ) {
00217 SUMA_SL_Err("Failed to compute mask from -cmask option");
00218 SUMA_free(Opt->mcdatav); Opt->mcdatav=NULL;
00219 SUMA_RETURN(NOPE);
00220 }
00221 if ( Opt->ninmask != Opt->nvox ) {
00222 SUMA_SL_Err("Input and cmask datasets do not have "
00223 "the same dimensions\n" );
00224 SUMA_free(Opt->mcdatav); Opt->mcdatav=NULL;
00225 SUMA_RETURN(NOPE);
00226 }
00227 Opt->ninmask = THD_countmask( Opt->ninmask, bmask );
00228 for (i=0; i<Opt->nvox; ++i) if (bmask[i]) Opt->mcdatav[i] = (double)bmask[i]; else Opt->mcdatav[i] = -1;
00229 free(bmask);bmask=NULL;
00230 } else {
00231 SUMA_SL_Err("NULL cmask"); SUMA_RETURN(NOPE);
00232 }
00233 break;
00234 case SUMA_ISO_VAL:
00235 case SUMA_ISO_RANGE:
00236
00237 DSET_load(Opt->in_vol);
00238 Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
00239 if (!Opt->dvec) {
00240 SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
00241 SUMA_RETURN(NOPE);
00242 }
00243 EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
00244 DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) ,
00245 MRI_double , Opt->dvec ) ;
00246
00247 PURGE_DSET(Opt->in_vol);
00248
00249 Opt->ninmask = 0;
00250 if (Opt->MaskMode == SUMA_ISO_VAL) {
00251 for (i=0; i<Opt->nvox; ++i) {
00252 if (Opt->dvec[i] == Opt->v0) {
00253 Opt->mcdatav[i] = 1; ++ Opt->ninmask;
00254 } else Opt->mcdatav[i] = -1;
00255 }
00256 } else if (Opt->MaskMode == SUMA_ISO_RANGE) {
00257 for (i=0; i<Opt->nvox; ++i) {
00258 if (Opt->dvec[i] >= Opt->v0 && Opt->dvec[i] < Opt->v1) {
00259 Opt->mcdatav[i] = 1; ++ Opt->ninmask;
00260 } else Opt->mcdatav[i] = -1;
00261 }
00262 } else {
00263 SUMA_SL_Err("Bad Miracle.");
00264 SUMA_RETURN(NOPE);
00265 }
00266 SUMA_free(Opt->dvec); Opt->dvec = NULL;
00267 break;
00268 default:
00269 SUMA_SL_Err("Unexpected value of MaskMode");
00270 SUMA_RETURN(NOPE);
00271 break;
00272 }
00273 } else if (Opt->xform == SUMA_ISO_XFORM_SHIFT) {
00274
00275 DSET_load(Opt->in_vol);
00276 Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
00277 if (!Opt->dvec) {
00278 SUMA_SL_Crit("Failed to allocate for dvec.\nOh misery.");
00279 SUMA_RETURN(NOPE);
00280 }
00281 EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
00282 DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) ,
00283 MRI_double , Opt->dvec ) ;
00284
00285 PURGE_DSET(Opt->in_vol);
00286 Opt->ninmask = Opt->nvox;
00287 for (i=0; i<Opt->nvox; ++i) {
00288 Opt->mcdatav[i] = Opt->dvec[i] - Opt->v0;
00289 }
00290 SUMA_free(Opt->dvec); Opt->dvec = NULL;
00291 } else if (Opt->xform == SUMA_ISO_XFORM_NONE) {
00292
00293 DSET_load(Opt->in_vol);
00294 Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
00295 if (!Opt->dvec) {
00296 SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
00297 SUMA_RETURN(NOPE);
00298 }
00299 EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
00300 DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) ,
00301 MRI_double , Opt->dvec ) ;
00302
00303 PURGE_DSET(Opt->in_vol);
00304 Opt->ninmask = Opt->nvox;
00305 for (i=0; i<Opt->nvox; ++i) {
00306 Opt->mcdatav[i] = Opt->dvec[i];
00307 }
00308 SUMA_free(Opt->dvec); Opt->dvec = NULL;
00309 } else {
00310 SUMA_SL_Err("Bad Opt->xform.");
00311 SUMA_RETURN(NOPE);
00312 }
00313
00314 if ( Opt->ninmask <= 0 ) {
00315 SUMA_SL_Err("No isovolume found!\nNothing to do." );
00316 SUMA_RETURN(NOPE);
00317 }
00318
00319 if (Opt->debug > 0) {
00320 fprintf( SUMA_STDERR, "%s:\nInput dset %s has nvox = %d, nvals = %d",
00321 FuncName, Opt->in_name, Opt->nvox, DSET_NVALS(Opt->in_vol) );
00322 fprintf( SUMA_STDERR, " (%d voxels in mask)\n", Opt->ninmask );
00323 }
00324
00325 SUMA_RETURN(YUP);
00326 }
00327
00328 #ifdef SUMA_IsoSurface_STANDALONE
00329
00330 static char SUMA_ISO_Obj_Types[][10] = { {"Cushin"}, {"Sphere"}, {"Plane"}, {"Cassini"}, {"Blooby"}, {"Chair"}, {"Cyclide"}, {"2 Torus"}, {"mc case"}, {"Drip"} };
00331
00332 void usage_SUMA_IsoSurface (SUMA_GENERIC_ARGV_PARSE *ps)
00333 {
00334 static char FuncName[]={"usage_SUMA_IsoSurface"};
00335 char * s = NULL, *sio=NULL;
00336 int i;
00337 s = SUMA_help_basics();
00338 sio = SUMA_help_IO_Args(ps);
00339 printf ( "\n"
00340 "Usage: A program to perform isosurface extraction from a volume.\n"
00341 " Based on code by Thomas Lewiner (see below).\n"
00342 "\n"
00343 " IsoSurface < -input VOL | -shape S GR >\n"
00344 " < -isoval V | -isorange V0 V1 | -isocmask MASK_COM >\n"
00345 " [< -o_TYPE PREFIX>]\n"
00346 " [< -debug DBG >]\n"
00347 "\n"
00348 " Mandatory parameters:\n"
00349 " You must use one of the following two options:\n"
00350 " -input VOL: Input volume.\n"
00351 " -shape S GR: Built in shape.\n"
00352 " where S is the shape number, \n"
00353 " between 0 and 9 (see below). \n"
00354 " and GR is the grid size (like 64).\n"
00355 " If you use -debug 1 with this option\n"
00356 " a .1D volume called mc_shape*.1D is\n"
00357 " written to disk. Watch the debug output\n"
00358 " for a command suggesting how to turn\n"
00359 " this 1D file into a BRIK volume for viewing\n"
00360 " in AFNI.\n"
00361 " You must use one of the following iso* options:\n"
00362 " -isoval V: Create isosurface where volume = V\n"
00363 " -isorange V0 V1: Create isosurface where V0 <= volume < V1\n"
00364 " -isocmask MASK_COM: Create isosurface where MASK_COM != 0\n"
00365 " For example: -isocmask '-a VOL+orig -expr (1-bool(a-V))' \n"
00366 " is equivalent to using -isoval V. \n"
00367 " NOTE: -isorange and -isocmask are only allowed with -xform mask\n"
00368 " See -xform below for details.\n"
00369 "\n"
00370 " Optional Parameters:\n"
00371 " -xform XFORM: Transform to apply to volume values\n"
00372 " before searching for sign change\n"
00373 " boundary. XFORM can be one of:\n"
00374 " mask: values that meet the iso* conditions\n"
00375 " are set to 1. All other values are set\n"
00376 " to -1. This is the default XFORM.\n"
00377 " shift: subtract V from the dataset and then \n"
00378 " search for 0 isosurface. This has the\n"
00379 " effect of constructing the V isosurface\n"
00380 " if your dataset has a continuum of values.\n"
00381 " This option can only be used with -isoval V.\n"
00382 " none: apply no transforms. This assumes that\n"
00383 " your volume has a continuum of values \n"
00384 " from negative to positive and that you\n"
00385 " are seeking to 0 isosurface.\n"
00386 " This option can only be used with -isoval 0.\n"
00387 " -o_TYPE PREFIX: prefix of output surface.\n"
00388 " where TYPE specifies the format of the surface\n"
00389 " and PREFIX is, well, the prefix.\n"
00390 " TYPE is one of: fs, 1d (or vec), sf, ply.\n"
00391 " Default is: -o_ply \n"
00392 "\n"
00393 "%s\n"
00394 "\n"
00395 " -debug DBG: debug levels of 0 (default), 1, 2, 3.\n"
00396 " This is no Rick Reynolds debug, which is oft nicer\n"
00397 " than the results, but it will do.\n"
00398 "\n"
00399 " Built In Shapes:\n", sio);
00400 for (i=0; i<10;++i) {
00401 printf( " %d: %s\n", i, SUMA_ISO_Obj_Types[i]);
00402 }
00403 printf ( "\n"
00404 " NOTE:\n"
00405 " The code for the heart of this program is a translation of:\n"
00406 " Thomas Lewiner's C++ implementation of the algorithm in:\n"
00407 " Efficient Implementation of Marching Cubes´ Cases with Topological Guarantees\n"
00408 " by Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira and Geovan Tavares \n"
00409 " in Journal of Graphics Tools. \n"
00410 " http://www-sop.inria.fr/prisme/personnel/Thomas.Lewiner/JGT.pdf\n"
00411 "\n"
00412 "%s"
00413 "\n", s);
00414 SUMA_free(s); s = NULL; SUMA_free(sio); sio = NULL;
00415 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
00416 printf(" Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \n");
00417 exit (0);
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_IsoSurface_ParseInput (char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
00430 {
00431 static char FuncName[]={"SUMA_IsoSurface_ParseInput"};
00432 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
00433 int kar, i, ind;
00434 char *outname;
00435 SUMA_Boolean brk = NOPE;
00436 SUMA_Boolean LocalHead = NOPE;
00437
00438 SUMA_ENTRY;
00439
00440 Opt = (SUMA_GENERIC_PROG_OPTIONS_STRUCT *)SUMA_malloc(sizeof(SUMA_GENERIC_PROG_OPTIONS_STRUCT));
00441
00442 kar = 1;
00443 Opt->spec_file = NULL;
00444 Opt->out_prefix = NULL;
00445 Opt->sv_name = NULL;
00446 Opt->N_surf = -1;
00447 Opt->in_name = NULL;
00448 Opt->cmask = NULL;
00449 Opt->MaskMode = SUMA_ISO_UNDEFINED;
00450 for (i=0; i<SUMA_GENERIC_PROG_MAX_SURF; ++i) { Opt->surf_names[i] = NULL; }
00451 outname = NULL;
00452 Opt->in_vol = NULL;
00453 Opt->nvox = -1;
00454 Opt->ninmask = -1;
00455 Opt->mcdatav = NULL;
00456 Opt->debug = 0;
00457 Opt->v0 = 0.0;
00458 Opt->v1 = -1.0;
00459 Opt->dvec = NULL;
00460 Opt->SurfFileType = SUMA_PLY;
00461 Opt->SurfFileFormat = SUMA_ASCII;
00462 Opt->xform = SUMA_ISO_XFORM_MASK;
00463 Opt->obj_type = -1;
00464 Opt->obj_type_res = -1;
00465 Opt->XYZ = NULL;
00466 Opt->in_1D = NULL;
00467 Opt->N_XYZ = 0;
00468 brk = NOPE;
00469 while (kar < argc) {
00470
00471 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
00472 usage_SUMA_IsoSurface(ps);
00473 exit (0);
00474 }
00475
00476 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
00477
00478 if (!brk && (strcmp(argv[kar], "-xform") == 0)) {
00479 kar ++;
00480 if (kar >= argc) {
00481 fprintf (SUMA_STDERR, "need argument after -xform \n");
00482 exit (1);
00483 }
00484 if (!strcmp(argv[kar], "mask")) {
00485 Opt->xform = SUMA_ISO_XFORM_MASK;
00486 } else if (!strcmp(argv[kar], "none")) {
00487 Opt->xform = SUMA_ISO_XFORM_NONE;
00488 } else if (!strcmp(argv[kar], "shift")) {
00489 Opt->xform = SUMA_ISO_XFORM_SHIFT;
00490 }else {
00491 fprintf (SUMA_STDERR, "%s is a bad parameter for -xform option. \n", argv[kar]);
00492 exit (1);
00493 }
00494 brk = YUP;
00495 }
00496
00497
00498 if (!brk && (strcmp(argv[kar], "-debug") == 0)) {
00499 kar ++;
00500 if (kar >= argc) {
00501 fprintf (SUMA_STDERR, "need argument after -debug \n");
00502 exit (1);
00503 }
00504 Opt->debug = atoi(argv[kar]);
00505 if (Opt->debug > 2) { LocalHead = YUP; }
00506 brk = YUP;
00507 }
00508
00509 if (!brk && (strcmp(argv[kar], "-isocmask") == 0)) {
00510 if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
00511 fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
00512 }
00513 kar ++;
00514 if (kar >= argc) {
00515 fprintf (SUMA_STDERR, "need argument after -isocmask \n");
00516 exit (1);
00517 }
00518 Opt->cmask = argv[kar];
00519 Opt->MaskMode = SUMA_ISO_CMASK;
00520 brk = YUP;
00521 }
00522
00523 if (!brk && (strcmp(argv[kar], "-isoval") == 0)) {
00524 if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
00525 fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
00526 }
00527 kar ++;
00528 if (kar >= argc) {
00529 fprintf (SUMA_STDERR, "need argument after -isoval \n");
00530 exit (1);
00531 }
00532 Opt->v0 = atof(argv[kar]);
00533 Opt->MaskMode = SUMA_ISO_VAL;
00534 brk = YUP;
00535 }
00536
00537 if (!brk && (strcmp(argv[kar], "-isorange") == 0)) {
00538 if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
00539 fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
00540 }
00541 kar ++;
00542 if (kar+1 >= argc) {
00543 fprintf (SUMA_STDERR, "need 2 arguments after -isorange \n");
00544 exit (1);
00545 }
00546 Opt->v0 = atof(argv[kar]);kar ++;
00547 Opt->v1 = atof(argv[kar]);
00548 Opt->MaskMode = SUMA_ISO_RANGE;
00549 if (Opt->v1 < Opt->v0) {
00550 fprintf (SUMA_STDERR, "range values wrong. Must have %f <= %f \n", Opt->v0, Opt->v1);
00551 exit (1);
00552 }
00553 brk = YUP;
00554 }
00555 if (!brk && (strcmp(argv[kar], "-input") == 0)) {
00556 kar ++;
00557 if (kar >= argc) {
00558 fprintf (SUMA_STDERR, "need argument after -input \n");
00559 exit (1);
00560 }
00561 Opt->in_name = SUMA_copy_string(argv[kar]);
00562 brk = YUP;
00563 }
00564
00565 if (!brk && (strcmp(argv[kar], "-shape") == 0)) {
00566 kar ++;
00567 if (kar+1 >= argc) {
00568 fprintf (SUMA_STDERR, "need 2 arguments after -shape \n");
00569 exit (1);
00570 }
00571 Opt->obj_type = atoi(argv[kar]); kar ++;
00572 Opt->obj_type_res = atoi(argv[kar]);
00573 if (Opt->obj_type < 0 || Opt->obj_type > 9) {
00574 fprintf (SUMA_STDERR, "Error %s:\nShape number (S) must be between 0 and 9. I have %d\n", FuncName, Opt->obj_type);
00575 exit (1);
00576 }
00577 if (Opt->obj_type_res < 0) {
00578 fprintf (SUMA_STDERR, "Error %s:\nShape grid resolution (GR) must be > 0 . I have %d\n", FuncName, Opt->obj_type_res);
00579 exit (1);
00580 }
00581 brk = YUP;
00582 }
00583
00584 if (!brk && !ps->arg_checked[kar]) {
00585 fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
00586 exit (1);
00587 } else {
00588 brk = NOPE;
00589 kar ++;
00590 }
00591 }
00592
00593
00594 if (ps->o_N_surfnames) {
00595 Opt->out_prefix = SUMA_copy_string(ps->o_surfnames[0]);
00596 Opt->SurfFileType = ps->o_FT[0];
00597 Opt->SurfFileFormat = ps->o_FF[0];
00598 }
00599
00600 if (Opt->in_name && Opt->obj_type >= 0) {
00601 fprintf (SUMA_STDERR,"Error %s:\nOptions -input and -shape are mutually exclusive.\n", FuncName);
00602 exit(1);
00603 }
00604 if (!Opt->in_name && Opt->obj_type < 0) {
00605 fprintf (SUMA_STDERR,"Error %s:\nEither -input or -shape options must be used.\n", FuncName);
00606 exit(1);
00607 }
00608 if (!Opt->out_prefix) Opt->out_prefix = SUMA_copy_string("isosurface_out");
00609 if (Opt->xform == SUMA_ISO_XFORM_NONE) {
00610 if (Opt->v0 != 0) {
00611 fprintf (SUMA_STDERR,"Error %s: Bad %f isovalue\nWith -xform none you can only extract the 0 isosurface.\n(i.e. -isoval 0)\n", FuncName, Opt->v0);
00612 exit(1);
00613 }
00614 if (Opt->MaskMode != SUMA_ISO_VAL) {
00615 fprintf (SUMA_STDERR,"Error %s: \nWith -xform none you can only use -isoval 0\n", FuncName);
00616 exit(1);
00617 }
00618 }
00619 if (Opt->xform == SUMA_ISO_XFORM_SHIFT) {
00620 if (Opt->MaskMode != SUMA_ISO_VAL) {
00621 fprintf (SUMA_STDERR,"Error %s: \nWith -xform shift you can only use -isoval val\n", FuncName);
00622 exit(1);
00623 }
00624 }
00625
00626 SUMA_RETURN(Opt);
00627 }
00628
00629 int main (int argc,char *argv[])
00630 {
00631 static char FuncName[]={"IsoSurface"};
00632 int i;
00633 void *SO_name=NULL;
00634 SUMA_SurfaceObject *SO = NULL;
00635 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
00636 char stmp[200];
00637 SUMA_Boolean exists = NOPE;
00638 SUMA_Boolean LocalHead = NOPE;
00639 SUMA_GENERIC_ARGV_PARSE *ps=NULL;
00640
00641 SUMA_mainENTRY;
00642 SUMA_STANDALONE_INIT;
00643
00644
00645
00646 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
00647 ps = SUMA_Parse_IO_Args(argc, argv, "-o;");
00648
00649 if (argc < 2) {
00650 usage_SUMA_IsoSurface(ps);
00651 exit (1);
00652 }
00653
00654 Opt = SUMA_IsoSurface_ParseInput (argv, argc, ps);
00655
00656 set_suma_debug(Opt->debug);
00657
00658 SO_name = SUMA_Prefix2SurfaceName(Opt->out_prefix, NULL, NULL, Opt->SurfFileType, &exists);
00659 if (exists) {
00660 fprintf(SUMA_STDERR,"Error %s:\nOutput file(s) %s* on disk.\nWill not overwrite.\n", FuncName, Opt->out_prefix);
00661 exit(1);
00662 }
00663
00664 if (Opt->obj_type < 0) {
00665 if (Opt->debug) {
00666 SUMA_S_Note("Creating mask...");
00667 }
00668 if (!SUMA_Get_isosurface_datasets (Opt)) {
00669 SUMA_SL_Err("Failed to get data.");
00670 exit(1);
00671 }
00672
00673 if (Opt->debug > 1) {
00674 if (Opt->debug == 2) {
00675 FILE *fout=fopen("inmaskvec.1D","w");
00676 SUMA_S_Note("Writing masked values...\n");
00677 if (!fout) {
00678 SUMA_SL_Err("Failed to write maskvec");
00679 exit(1);
00680 }
00681 fprintf(fout, "#Col. 0 Voxel Index\n"
00682 "#Col. 1 Is a mask (all values here should be 1)\n" );
00683 for (i=0; i<Opt->nvox; ++i) {
00684 if (Opt->mcdatav[i]) {
00685 fprintf(fout,"%d %.2f\n", i, Opt->mcdatav[i]);
00686 }
00687 }
00688 fclose(fout); fout = NULL;
00689 } else {
00690 FILE *fout=fopen("maskvec.1D","w");
00691 SUMA_S_Note("Writing all mask values...\n");
00692 if (!fout) {
00693 SUMA_S_Err("Failed to write maskvec");
00694 exit(1);
00695 }
00696 fprintf(fout, "#Col. 0 Voxel Index\n"
00697 "#Col. 1 Is in mask ?\n" );
00698 for (i=0; i<Opt->nvox; ++i) {
00699 fprintf(fout,"%d %.2f\n", i, Opt->mcdatav[i]);
00700 }
00701 fclose(fout); fout = NULL;
00702 }
00703 }
00704 } else {
00705 if (Opt->debug) {
00706 SUMA_S_Note("Using built in types...");
00707 }
00708 }
00709
00710 if (!(SO = SUMA_MarchingCubesSurface(Opt))) {
00711 SUMA_S_Err("Failed to create surface.\n");
00712 exit(1);
00713 }
00714
00715
00716 if (!SUMA_Save_Surface_Object (SO_name, SO, Opt->SurfFileType, Opt->SurfFileFormat, NULL)) {
00717 fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
00718 exit (1);
00719 }
00720
00721 if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
00722 if (Opt->dvec) SUMA_free(Opt->dvec); Opt->dvec = NULL;
00723 if (Opt->mcdatav) {SUMA_free(Opt->mcdatav); Opt->mcdatav = NULL;}
00724 if (Opt->in_vol) { DSET_delete( Opt->in_vol); Opt->in_vol = NULL;}
00725 if (Opt->out_prefix) SUMA_free(Opt->out_prefix); Opt->out_prefix = NULL;
00726 if (Opt) SUMA_free(Opt);
00727 if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) {
00728 SUMA_SL_Err("DO Cleanup Failed!");
00729 }
00730 if (SO_name) SUMA_free(SO_name); SO_name = NULL;
00731 if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
00732 exit(0);
00733 }
00734 #endif
