/*! File to contain routines for creating isosurfaces.
NOTE: MarchingCube code was translated from Thomas Lewiner's C++
implementation of the paper:
Efficient Implementation of Marching Cubes´ Cases with Topological Guarantees
by Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira and Geovan Tavares
in Journal of Graphics Tools.
http://www-sop.inria.fr/prisme/personnel/Thomas.Lewiner/JGT.pdf
*/
#include "SUMA_suma.h"
#include "MarchingCubes/MarchingCubes.h"
#undef STAND_ALONE
#if defined SUMA_IsoSurface_STANDALONE
#define STAND_ALONE
#endif
#ifdef STAND_ALONE
/* these global variables must be declared even if they will not be used by this main */
SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/
SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures
SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
SUMA_DO *SUMAg_DOv = NULL; /*!< Global pointer to Displayable Object structure vector*/
int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */
#else
extern SUMA_CommonFields *SUMAg_CF;
extern SUMA_DO *SUMAg_DOv;
extern SUMA_SurfaceViewer *SUMAg_SVv;
extern int SUMAg_N_SVv;
extern int SUMAg_N_DOv;
#endif
/*! a macro version of MarchingCubes' set_data */
#define SUMA_SET_MC_DATA(mcb, val, i, j, k) { \
mcb->data[ i + j*mcb->size_x + k*mcb->size_x*mcb->size_y] = val; \
}
/*!
A function version of the program mc by Thomas Lewiner
see main.c in ./MarchingCubes
*/
SUMA_SurfaceObject *SUMA_MarchingCubesSurface(SUMA_GENERIC_PROG_OPTIONS_STRUCT * Opt)
{
static char FuncName[]={"SUMA_MarchingCubesSurface"};
SUMA_SurfaceObject *SO=NULL;
int nxx, nyy, nzz, cnt, i, j, k, *FaceSetList=NULL;
float *NodeList=NULL;
SUMA_NEW_SO_OPT *nsoopt = NULL;
THD_fvec3 fv, iv;
MCB *mcp ;
SUMA_ENTRY;
if (Opt->obj_type < 0) {
nxx = DSET_NX(Opt->in_vol);
nyy = DSET_NY(Opt->in_vol);
nzz = DSET_NZ(Opt->in_vol);
if (Opt->debug) {
fprintf(SUMA_STDERR,"%s:\nNxx=%d\tNyy=%d\tNzz=%d\n", FuncName, nxx, nyy, nzz);
}
mcp = MarchingCubes(-1, -1, -1);
set_resolution( mcp, nxx, nyy, nzz ) ;
init_all(mcp) ;
if (Opt->debug) fprintf(SUMA_STDERR,"%s:\nSetting data...\n", FuncName);
cnt = 0;
for( k = 0 ; k < mcp->size_z ; k++ ) {
for( j = 0 ; j < mcp->size_y ; j++ ) {
for( i = 0 ; i < mcp->size_x ; i++ ) {
SUMA_SET_MC_DATA ( mcp, Opt->mcdatav[cnt], i, j, k);
++cnt;
}
}
}
} else {
/* built in */
nxx = nyy = nzz = Opt->obj_type_res;
mcp = MarchingCubes(-1, -1, -1);
set_resolution( mcp, nxx, nyy, nzz) ;
init_all(mcp) ;
compute_data( *mcp , Opt->obj_type) ;
}
if (Opt->debug) fprintf(SUMA_STDERR,"%s:\nrunning MarchingCubes...\n", FuncName);
run(mcp) ;
clean_temps(mcp) ;
if (Opt->debug > 1) {
fprintf(SUMA_STDERR,"%s:\nwriting out NodeList and FaceSetList...\n", FuncName);
write1Dmcb(mcp);
}
if (Opt->debug) {
fprintf(SUMA_STDERR,"%s:\nNow creating SO...\n", FuncName);
}
NodeList = (float *)SUMA_malloc(sizeof(float)*3*mcp->nverts);
FaceSetList = (int *)SUMA_malloc(sizeof(int)*3*mcp->ntrigs);
if (!NodeList || !FaceSetList) {
SUMA_SL_Crit("Failed to allocate!");
SUMA_RETURN(SO);
}
nsoopt = SUMA_NewNewSOOpt();
if (Opt->obj_type < 0) {
nsoopt->LargestBoxSize = -1;
if (Opt->debug) {
fprintf(SUMA_STDERR,"%s:\nCopying vertices, changing to DICOM \nOrig:(%f %f %f) \nD:(%f %f %f)...\n",
FuncName, DSET_XORG(Opt->in_vol), DSET_YORG(Opt->in_vol), DSET_ZORG(Opt->in_vol),
DSET_DX(Opt->in_vol), DSET_DY(Opt->in_vol), DSET_DZ(Opt->in_vol));
}
for ( i = 0; i < mcp->nverts; i++ ) {
j = 3*i; /* change from index coordinates to mm DICOM, next three lines are equivalent of SUMA_THD_3dfind_to_3dmm*/
fv.xyz[0] = DSET_XORG(Opt->in_vol) + mcp->vertices[i].x * DSET_DX(Opt->in_vol);
fv.xyz[1] = DSET_YORG(Opt->in_vol) + mcp->vertices[i].y * DSET_DY(Opt->in_vol);
fv.xyz[2] = DSET_ZORG(Opt->in_vol) + mcp->vertices[i].z * DSET_DZ(Opt->in_vol);
/* change mm to RAI coords */
iv = SUMA_THD_3dmm_to_dicomm( Opt->in_vol->daxes->xxorient, Opt->in_vol->daxes->yyorient, Opt->in_vol->daxes->zzorient, fv );
NodeList[j ] = iv.xyz[0];
NodeList[j+1] = iv.xyz[1];
NodeList[j+2] = iv.xyz[2];
}
for ( i = 0; i < mcp->ntrigs; i++ ) {
j = 3*i;
FaceSetList[j ] = mcp->triangles[i].v3;
FaceSetList[j+1] = mcp->triangles[i].v2;
FaceSetList[j+2] = mcp->triangles[i].v1;
}
} else {
nsoopt->LargestBoxSize = 100;
/* built in */
for ( i = 0; i < mcp->nverts; i++ ) {
j = 3*i;
NodeList[j ] = mcp->vertices[i].x;
NodeList[j+1] = mcp->vertices[i].y;
NodeList[j+2] = mcp->vertices[i].z;
}
for ( i = 0; i < mcp->ntrigs; i++ ) {
j = 3*i;
FaceSetList[j ] = mcp->triangles[i].v3;
FaceSetList[j+1] = mcp->triangles[i].v2;
FaceSetList[j+2] = mcp->triangles[i].v1;
}
}
SO = SUMA_NewSO(&NodeList, mcp->nverts, &FaceSetList, mcp->ntrigs, nsoopt);
if (Opt->obj_type < 0) {
/* not sure if anything needs to be done here ...*/
} else {
if (Opt->obj_type == 0) SO->normdir = 1;
else SO->normdir = -1;
}
if (Opt->debug) {
fprintf(SUMA_STDERR,"%s:\nCleaning mcp...\n", FuncName);
}
clean_all(mcp) ;
free(mcp);
nsoopt=SUMA_FreeNewSOOpt(nsoopt);
SUMA_RETURN(SO);
}
/*!
contains code shamelessly stolen from Rick who stole it from Bob.
*/
SUMA_Boolean SUMA_Get_isosurface_datasets (SUMA_GENERIC_PROG_OPTIONS_STRUCT * Opt)
{
static char FuncName[]={"SUMA_Get_isosurface_datasets"};
int i;
SUMA_ENTRY;
Opt->in_vol = THD_open_dataset( Opt->in_name );
if (!ISVALID_DSET(Opt->in_vol)) {
if (!Opt->in_name) {
SUMA_SL_Err("NULL input volume.");
SUMA_RETURN(NOPE);
} else {
SUMA_SL_Err("invalid volume.");
SUMA_RETURN(NOPE);
}
} else if ( DSET_BRICK_TYPE(Opt->in_vol, 0) == MRI_complex) {
SUMA_SL_Err("Can't do complex data.");
SUMA_RETURN(NOPE);
}
Opt->nvox = DSET_NVOX( Opt->in_vol );
if (DSET_NVALS( Opt->in_vol) != 1) {
SUMA_SL_Err("Input volume can only have one sub-brick in it.\nUse [.] selectors to choose sub-brick needed.");
SUMA_RETURN(NOPE);
}
Opt->mcdatav = (double *)SUMA_malloc(sizeof(double)*Opt->nvox);
if (!Opt->mcdatav) {
SUMA_SL_Crit("Failed to allocate for maskv");
SUMA_RETURN(NOPE);
}
if (Opt->xform == SUMA_ISO_XFORM_MASK) {
switch (Opt->MaskMode) {
case SUMA_ISO_CMASK:
if (Opt->cmask) {
/* here's the second order grand theft */
int clen = strlen( Opt->cmask );
char * cmd;
byte *bmask;
/* save original cmask command, as EDT_calcmask() is destructive */
cmd = (char *)malloc((clen + 1) * sizeof(char));
strcpy( cmd, Opt->cmask);
bmask = EDT_calcmask( cmd, &Opt->ninmask );
free( cmd ); /* free EDT_calcmask() string */
if ( bmask == NULL ) {
SUMA_SL_Err("Failed to compute mask from -cmask option");
SUMA_free(Opt->mcdatav); Opt->mcdatav=NULL;
SUMA_RETURN(NOPE);
}
if ( Opt->ninmask != Opt->nvox ) {
SUMA_SL_Err("Input and cmask datasets do not have "
"the same dimensions\n" );
SUMA_free(Opt->mcdatav); Opt->mcdatav=NULL;
SUMA_RETURN(NOPE);
}
Opt->ninmask = THD_countmask( Opt->ninmask, bmask );
for (i=0; i<Opt->nvox; ++i) if (bmask[i]) Opt->mcdatav[i] = (double)bmask[i]; else Opt->mcdatav[i] = -1;
free(bmask);bmask=NULL;
} else {
SUMA_SL_Err("NULL cmask"); SUMA_RETURN(NOPE);
}
break;
case SUMA_ISO_VAL:
case SUMA_ISO_RANGE:
/* load the dset */
DSET_load(Opt->in_vol);
Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
if (!Opt->dvec) {
SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
SUMA_RETURN(NOPE);
}
EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) , /* input */
MRI_double , Opt->dvec ) ; /* output */
/* no need for data in input volume anymore */
PURGE_DSET(Opt->in_vol);
Opt->ninmask = 0;
if (Opt->MaskMode == SUMA_ISO_VAL) {
for (i=0; i<Opt->nvox; ++i) {
if (Opt->dvec[i] == Opt->v0) {
Opt->mcdatav[i] = 1; ++ Opt->ninmask;
} else Opt->mcdatav[i] = -1;
}
} else if (Opt->MaskMode == SUMA_ISO_RANGE) {
for (i=0; i<Opt->nvox; ++i) {
if (Opt->dvec[i] >= Opt->v0 && Opt->dvec[i] < Opt->v1) {
Opt->mcdatav[i] = 1; ++ Opt->ninmask;
} else Opt->mcdatav[i] = -1;
}
} else {
SUMA_SL_Err("Bad Miracle.");
SUMA_RETURN(NOPE);
}
SUMA_free(Opt->dvec); Opt->dvec = NULL; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
break;
default:
SUMA_SL_Err("Unexpected value of MaskMode");
SUMA_RETURN(NOPE);
break;
}
} else if (Opt->xform == SUMA_ISO_XFORM_SHIFT) {
/* load the dset */
DSET_load(Opt->in_vol);
Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
if (!Opt->dvec) {
SUMA_SL_Crit("Failed to allocate for dvec.\nOh misery.");
SUMA_RETURN(NOPE);
}
EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) , /* input */
MRI_double , Opt->dvec ) ; /* output */
/* no need for data in input volume anymore */
PURGE_DSET(Opt->in_vol);
Opt->ninmask = Opt->nvox;
for (i=0; i<Opt->nvox; ++i) {
Opt->mcdatav[i] = Opt->dvec[i] - Opt->v0;
}
SUMA_free(Opt->dvec); Opt->dvec = NULL; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
} else if (Opt->xform == SUMA_ISO_XFORM_NONE) {
/* load the dset */
DSET_load(Opt->in_vol);
Opt->dvec = (double *)SUMA_malloc(sizeof(double) * Opt->nvox);
if (!Opt->dvec) {
SUMA_SL_Crit("Faile to allocate for dvec.\nOh misery.");
SUMA_RETURN(NOPE);
}
EDIT_coerce_scale_type( Opt->nvox , DSET_BRICK_FACTOR(Opt->in_vol,0) ,
DSET_BRICK_TYPE(Opt->in_vol,0), DSET_ARRAY(Opt->in_vol, 0) , /* input */
MRI_double , Opt->dvec ) ; /* output */
/* no need for data in input volume anymore */
PURGE_DSET(Opt->in_vol);
Opt->ninmask = Opt->nvox;
for (i=0; i<Opt->nvox; ++i) {
Opt->mcdatav[i] = Opt->dvec[i];
}
SUMA_free(Opt->dvec); Opt->dvec = NULL; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
} else {
SUMA_SL_Err("Bad Opt->xform.");
SUMA_RETURN(NOPE);
}
if ( Opt->ninmask <= 0 ) {
SUMA_SL_Err("No isovolume found!\n Nothing to do." );
SUMA_RETURN(NOPE);
}
if (Opt->debug > 0) {
fprintf( SUMA_STDERR, "%s:\nInput dset %s has nvox = %d, nvals = %d",
FuncName, Opt->in_name, Opt->nvox, DSET_NVALS(Opt->in_vol) );
fprintf( SUMA_STDERR, " (%d voxels in mask)\n", Opt->ninmask );
}
SUMA_RETURN(YUP);
}
#ifdef SUMA_IsoSurface_STANDALONE
static char SUMA_ISO_Obj_Types[][10] = { {"Cushin"}, {"Sphere"}, {"Plane"}, {"Cassini"}, {"Blooby"}, {"Chair"}, {"Cyclide"}, {"2 Torus"}, {"mc case"}, {"Drip"} };
void usage_SUMA_IsoSurface (SUMA_GENERIC_ARGV_PARSE *ps)
{
static char FuncName[]={"usage_SUMA_IsoSurface"};
char * s = NULL, *sio=NULL;
int i;
s = SUMA_help_basics();
sio = SUMA_help_IO_Args(ps);
printf ( "\n"
"Usage: A program to perform isosurface extraction from a volume.\n"
" Based on code by Thomas Lewiner (see below).\n"
"\n"
" IsoSurface < -input VOL | -shape S GR >\n"
" < -isoval V | -isorange V0 V1 | -isocmask MASK_COM >\n"
" [< -o_TYPE PREFIX>]\n"
" [< -debug DBG >]\n"
"\n"
" Mandatory parameters:\n"
" You must use one of the following two options:\n"
" -input VOL: Input volume.\n"
" -shape S GR: Built in shape.\n"
" where S is the shape number, \n"
" between 0 and 9 (see below). \n"
" and GR is the grid size (like 64).\n"
" If you use -debug 1 with this option\n"
" a .1D volume called mc_shape*.1D is\n"
" written to disk. Watch the debug output\n"
" for a command suggesting how to turn\n"
" this 1D file into a BRIK volume for viewing\n"
" in AFNI.\n"
" You must use one of the following iso* options:\n"
" -isoval V: Create isosurface where volume = V\n"
" -isorange V0 V1: Create isosurface where V0 <= volume < V1\n"
" -isocmask MASK_COM: Create isosurface where MASK_COM != 0\n"
" For example: -isocmask '-a VOL+orig -expr (1-bool(a-V))' \n"
" is equivalent to using -isoval V. \n"
" NOTE: -isorange and -isocmask are only allowed with -xform mask\n"
" See -xform below for details.\n"
"\n"
" Optional Parameters:\n"
" -xform XFORM: Transform to apply to volume values\n"
" before searching for sign change\n"
" boundary. XFORM can be one of:\n"
" mask: values that meet the iso* conditions\n"
" are set to 1. All other values are set\n"
" to -1. This is the default XFORM.\n"
" shift: subtract V from the dataset and then \n"
" search for 0 isosurface. This has the\n"
" effect of constructing the V isosurface\n"
" if your dataset has a continuum of values.\n"
" This option can only be used with -isoval V.\n"
" none: apply no transforms. This assumes that\n"
" your volume has a continuum of values \n"
" from negative to positive and that you\n"
" are seeking to 0 isosurface.\n"
" This option can only be used with -isoval 0.\n"
" -o_TYPE PREFIX: prefix of output surface.\n"
" where TYPE specifies the format of the surface\n"
" and PREFIX is, well, the prefix.\n"
" TYPE is one of: fs, 1d (or vec), sf, ply.\n"
" Default is: -o_ply \n"
"\n"
"%s\n"
"\n"
" -debug DBG: debug levels of 0 (default), 1, 2, 3.\n"
" This is no Rick Reynolds debug, which is oft nicer\n"
" than the results, but it will do.\n"
"\n"
" Built In Shapes:\n", sio);
for (i=0; i<10;++i) {
printf( " %d: %s\n", i, SUMA_ISO_Obj_Types[i]);
}
printf ( "\n"
" NOTE:\n"
" The code for the heart of this program is a translation of:\n"
" Thomas Lewiner's C++ implementation of the algorithm in:\n"
" Efficient Implementation of Marching Cubes´ Cases with Topological Guarantees\n"
" by Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira and Geovan Tavares \n"
" in Journal of Graphics Tools. \n"
" http://www-sop.inria.fr/prisme/personnel/Thomas.Lewiner/JGT.pdf\n"
"\n"
"%s"
"\n", s);
SUMA_free(s); s = NULL; SUMA_free(sio); sio = NULL;
s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
printf(" Ziad S. Saad SSCC/NIMH/NIH saadz@mail.nih.gov \n");
exit (0);
}
/*!
\brief parse the arguments for SurfSmooth program
\param argv (char *)
\param argc (int)
\return Opt (SUMA_GENERIC_PROG_OPTIONS_STRUCT *) options structure.
To free it, use
SUMA_free(Opt->out_name);
SUMA_free(Opt);
*/
SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_IsoSurface_ParseInput (char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
{
static char FuncName[]={"SUMA_IsoSurface_ParseInput"};
SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
int kar, i, ind;
char *outname;
SUMA_Boolean brk = NOPE;
SUMA_Boolean LocalHead = NOPE;
SUMA_ENTRY;
Opt = (SUMA_GENERIC_PROG_OPTIONS_STRUCT *)SUMA_malloc(sizeof(SUMA_GENERIC_PROG_OPTIONS_STRUCT));
kar = 1;
Opt->spec_file = NULL;
Opt->out_prefix = NULL;
Opt->sv_name = NULL;
Opt->N_surf = -1;
Opt->in_name = NULL;
Opt->cmask = NULL;
Opt->MaskMode = SUMA_ISO_UNDEFINED;
for (i=0; i<SUMA_GENERIC_PROG_MAX_SURF; ++i) { Opt->surf_names[i] = NULL; }
outname = NULL;
Opt->in_vol = NULL;
Opt->nvox = -1;
Opt->ninmask = -1;
Opt->mcdatav = NULL;
Opt->debug = 0;
Opt->v0 = 0.0;
Opt->v1 = -1.0;
Opt->dvec = NULL;
Opt->SurfFileType = SUMA_PLY;
Opt->SurfFileFormat = SUMA_ASCII;
Opt->xform = SUMA_ISO_XFORM_MASK;
Opt->obj_type = -1;
Opt->obj_type_res = -1;
Opt->XYZ = NULL;
Opt->in_1D = NULL;
Opt->N_XYZ = 0;
brk = NOPE;
while (kar < argc) { /* loop accross command ine options */
/*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
usage_SUMA_IsoSurface(ps);
exit (0);
}
SUMA_SKIP_COMMON_OPTIONS(brk, kar);
if (!brk && (strcmp(argv[kar], "-xform") == 0)) {
kar ++;
if (kar >= argc) {
fprintf (SUMA_STDERR, "need argument after -xform \n");
exit (1);
}
if (!strcmp(argv[kar], "mask")) {
Opt->xform = SUMA_ISO_XFORM_MASK;
} else if (!strcmp(argv[kar], "none")) {
Opt->xform = SUMA_ISO_XFORM_NONE;
} else if (!strcmp(argv[kar], "shift")) {
Opt->xform = SUMA_ISO_XFORM_SHIFT;
}else {
fprintf (SUMA_STDERR, "%s is a bad parameter for -xform option. \n", argv[kar]);
exit (1);
}
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-debug") == 0)) {
kar ++;
if (kar >= argc) {
fprintf (SUMA_STDERR, "need argument after -debug \n");
exit (1);
}
Opt->debug = atoi(argv[kar]);
if (Opt->debug > 2) { LocalHead = YUP; }
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-isocmask") == 0)) {
if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
}
kar ++;
if (kar >= argc) {
fprintf (SUMA_STDERR, "need argument after -isocmask \n");
exit (1);
}
Opt->cmask = argv[kar];
Opt->MaskMode = SUMA_ISO_CMASK;
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-isoval") == 0)) {
if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
}
kar ++;
if (kar >= argc) {
fprintf (SUMA_STDERR, "need argument after -isoval \n");
exit (1);
}
Opt->v0 = atof(argv[kar]);
Opt->MaskMode = SUMA_ISO_VAL;
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-isorange") == 0)) {
if (Opt->MaskMode != SUMA_ISO_UNDEFINED) {
fprintf (SUMA_STDERR, "only one masking mode (-iso*) allowed.\n");
}
kar ++;
if (kar+1 >= argc) {
fprintf (SUMA_STDERR, "need 2 arguments after -isorange \n");
exit (1);
}
Opt->v0 = atof(argv[kar]);kar ++;
Opt->v1 = atof(argv[kar]);
Opt->MaskMode = SUMA_ISO_RANGE;
if (Opt->v1 < Opt->v0) {
fprintf (SUMA_STDERR, "range values wrong. Must have %f <= %f \n", Opt->v0, Opt->v1);
exit (1);
}
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-input") == 0)) {
kar ++;
if (kar >= argc) {
fprintf (SUMA_STDERR, "need argument after -input \n");
exit (1);
}
Opt->in_name = SUMA_copy_string(argv[kar]);
brk = YUP;
}
if (!brk && (strcmp(argv[kar], "-shape") == 0)) {
kar ++;
if (kar+1 >= argc) {
fprintf (SUMA_STDERR, "need 2 arguments after -shape \n");
exit (1);
}
Opt->obj_type = atoi(argv[kar]); kar ++;
Opt->obj_type_res = atoi(argv[kar]);
if (Opt->obj_type < 0 || Opt->obj_type > 9) {
fprintf (SUMA_STDERR, "Error %s:\nShape number (S) must be between 0 and 9. I have %d\n", FuncName, Opt->obj_type);
exit (1);
}
if (Opt->obj_type_res < 0) {
fprintf (SUMA_STDERR, "Error %s:\nShape grid resolution (GR) must be > 0 . I have %d\n", FuncName, Opt->obj_type_res);
exit (1);
}
brk = YUP;
}
if (!brk && !ps->arg_checked[kar]) {
fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
exit (1);
} else {
brk = NOPE;
kar ++;
}
}
/* transfer some options to Opt from ps. Clunky because this is retrofitting */
if (ps->o_N_surfnames) {
Opt->out_prefix = SUMA_copy_string(ps->o_surfnames[0]);
Opt->SurfFileType = ps->o_FT[0];
Opt->SurfFileFormat = ps->o_FF[0];
}
if (Opt->in_name && Opt->obj_type >= 0) {
fprintf (SUMA_STDERR,"Error %s:\nOptions -input and -shape are mutually exclusive.\n", FuncName);
exit(1);
}
if (!Opt->in_name && Opt->obj_type < 0) {
fprintf (SUMA_STDERR,"Error %s:\nEither -input or -shape options must be used.\n", FuncName);
exit(1);
}
if (!Opt->out_prefix) Opt->out_prefix = SUMA_copy_string("isosurface_out");
if (Opt->xform == SUMA_ISO_XFORM_NONE) {
if (Opt->v0 != 0) {
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);
exit(1);
}
if (Opt->MaskMode != SUMA_ISO_VAL) {
fprintf (SUMA_STDERR,"Error %s: \nWith -xform none you can only use -isoval 0\n", FuncName);
exit(1);
}
}
if (Opt->xform == SUMA_ISO_XFORM_SHIFT) {
if (Opt->MaskMode != SUMA_ISO_VAL) {
fprintf (SUMA_STDERR,"Error %s: \nWith -xform shift you can only use -isoval val\n", FuncName);
exit(1);
}
}
SUMA_RETURN(Opt);
}
int main (int argc,char *argv[])
{/* Main */
static char FuncName[]={"IsoSurface"};
int i;
void *SO_name=NULL;
SUMA_SurfaceObject *SO = NULL;
SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
char stmp[200];
SUMA_Boolean exists = NOPE;
SUMA_Boolean LocalHead = NOPE;
SUMA_GENERIC_ARGV_PARSE *ps=NULL;
SUMA_STANDALONE_INIT;
SUMA_mainENTRY;
/* Allocate space for DO structure */
SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
ps = SUMA_Parse_IO_Args(argc, argv, "-o;");
if (argc < 2) {
usage_SUMA_IsoSurface(ps);
exit (1);
}
Opt = SUMA_IsoSurface_ParseInput (argv, argc, ps);
set_suma_debug(Opt->debug);
SO_name = SUMA_Prefix2SurfaceName(Opt->out_prefix, NULL, NULL, Opt->SurfFileType, &exists);
if (exists) {
fprintf(SUMA_STDERR,"Error %s:\nOutput file(s) %s* on disk.\nWill not overwrite.\n", FuncName, Opt->out_prefix);
exit(1);
}
if (Opt->obj_type < 0) {
if (Opt->debug) {
SUMA_S_Note("Creating mask...");
}
if (!SUMA_Get_isosurface_datasets (Opt)) {
SUMA_SL_Err("Failed to get data.");
exit(1);
}
if (Opt->debug > 1) {
if (Opt->debug == 2) {
FILE *fout=fopen("inmaskvec.1D","w");
SUMA_S_Note("Writing masked values...\n");
if (!fout) {
SUMA_SL_Err("Failed to write maskvec");
exit(1);
}
fprintf(fout, "#Col. 0 Voxel Index\n"
"#Col. 1 Is a mask (all values here should be 1)\n" );
for (i=0; i<Opt->nvox; ++i) {
if (Opt->mcdatav[i]) {
fprintf(fout,"%d %.2f\n", i, Opt->mcdatav[i]);
}
}
fclose(fout); fout = NULL;
} else {
FILE *fout=fopen("maskvec.1D","w");
SUMA_S_Note("Writing all mask values...\n");
if (!fout) {
SUMA_S_Err("Failed to write maskvec");
exit(1);
}
fprintf(fout, "#Col. 0 Voxel Index\n"
"#Col. 1 Is in mask ?\n" );
for (i=0; i<Opt->nvox; ++i) {
fprintf(fout,"%d %.2f\n", i, Opt->mcdatav[i]);
}
fclose(fout); fout = NULL;
}
}
} else {
if (Opt->debug) {
SUMA_S_Note("Using built in types...");
}
}
/* Now call Marching Cube functions */
if (!(SO = SUMA_MarchingCubesSurface(Opt))) {
SUMA_S_Err("Failed to create surface.\n");
exit(1);
}
/* write the surface to disk */
if (!SUMA_Save_Surface_Object (SO_name, SO, Opt->SurfFileType, Opt->SurfFileFormat, NULL)) {
fprintf (SUMA_STDERR,"Error %s: Failed to write surface object.\n", FuncName);
exit (1);
}
if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
if (Opt->dvec) SUMA_free(Opt->dvec); Opt->dvec = NULL;
if (Opt->mcdatav) {SUMA_free(Opt->mcdatav); Opt->mcdatav = NULL;}
if (Opt->in_vol) { DSET_delete( Opt->in_vol); Opt->in_vol = NULL;}
if (Opt->out_prefix) SUMA_free(Opt->out_prefix); Opt->out_prefix = NULL;
if (Opt) SUMA_free(Opt);
if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) {
SUMA_SL_Err("DO Cleanup Failed!");
}
if (SO_name) SUMA_free(SO_name); SO_name = NULL;
if (!SUMA_Free_CommonFields(SUMAg_CF)) SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
exit(0);
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1