Doxygen Source Code Documentation
SUMA_IsoSurface.h File Reference
Go to the source code of this file.
Enumeration Type Documentation
|
|
Header for functions in SUMA_Isosurface.c 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 Definition at line 14 of file SUMA_IsoSurface.h.
|
|
|
Definition at line 15 of file SUMA_IsoSurface.h.
|
Function Documentation
|
|
contains code shamelessly stolen from Rick who stole it from Bob. Definition at line 166 of file SUMA_IsoSurface.c. References SUMA_GENERIC_PROG_OPTIONS_STRUCT::cmask, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_load, DSET_NVALS, DSET_NVOX, SUMA_GENERIC_PROG_OPTIONS_STRUCT::dvec, EDIT_coerce_scale_type(), EDT_calcmask(), free, i, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_name, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, ISVALID_DSET, malloc, SUMA_GENERIC_PROG_OPTIONS_STRUCT::MaskMode, SUMA_GENERIC_PROG_OPTIONS_STRUCT::mcdatav, SUMA_GENERIC_PROG_OPTIONS_STRUCT::ninmask, SUMA_GENERIC_PROG_OPTIONS_STRUCT::nvox, PURGE_DSET, SUMA_Boolean, SUMA_ENTRY, SUMA_free, SUMA_ISO_CMASK, SUMA_ISO_RANGE, SUMA_ISO_VAL, SUMA_ISO_XFORM_MASK, SUMA_ISO_XFORM_NONE, SUMA_ISO_XFORM_SHIFT, SUMA_malloc, SUMA_RETURN, SUMA_SL_Crit, SUMA_SL_Err, THD_countmask(), THD_open_dataset(), SUMA_GENERIC_PROG_OPTIONS_STRUCT::v0, SUMA_GENERIC_PROG_OPTIONS_STRUCT::v1, and SUMA_GENERIC_PROG_OPTIONS_STRUCT::xform.
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 /* here's the second order grand theft */
00204 int clen = strlen( Opt->cmask );
00205 char * cmd;
00206 byte *bmask;
00207
00208 /* save original cmask command, as EDT_calcmask() is destructive */
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 ); /* free EDT_calcmask() string */
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 /* load the dset */
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) , /* input */
00245 MRI_double , Opt->dvec ) ; /* output */
00246 /* no need for data in input volume anymore */
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; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
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 /* load the dset */
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) , /* input */
00283 MRI_double , Opt->dvec ) ; /* output */
00284 /* no need for data in input volume anymore */
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; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
00291 } else if (Opt->xform == SUMA_ISO_XFORM_NONE) {
00292 /* load the dset */
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) , /* input */
00301 MRI_double , Opt->dvec ) ; /* output */
00302 /* no need for data in input volume anymore */
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; /* this vector is not even created in SUMA_ISO_CMASK mode ...*/
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 }
|
|
|
A function version of the program mc by Thomas Lewiner see main.c in ./MarchingCubes Definition at line 47 of file SUMA_IsoSurface.c. References THD_3dim_dataset::daxes, SUMA_GENERIC_PROG_OPTIONS_STRUCT::debug, DSET_DX, DSET_DY, DSET_DZ, DSET_NX, DSET_NY, DSET_NZ, DSET_XORG, DSET_YORG, DSET_ZORG, free, i, SUMA_GENERIC_PROG_OPTIONS_STRUCT::in_vol, SUMA_GENERIC_PROG_OPTIONS_STRUCT::mcdatav, SUMA_GENERIC_PROG_OPTIONS_STRUCT::obj_type, SUMA_GENERIC_PROG_OPTIONS_STRUCT::obj_type_res, SUMA_ENTRY, SUMA_FreeNewSOOpt(), SUMA_malloc, SUMA_NewNewSOOpt(), SUMA_NewSO(), SUMA_RETURN, SUMA_SET_MC_DATA, SUMA_SL_Crit, SUMA_THD_3dmm_to_dicomm(), THD_dataxes::xxorient, THD_fvec3::xyz, THD_dataxes::yyorient, and THD_dataxes::zzorient. Referenced by SUMA_GimmeSomeSOs().
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 /* built in */
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; /* change from index coordinates to mm DICOM, next three lines are equivalent of SUMA_THD_3dfind_to_3dmm*/
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 /* change mm to RAI coords */
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 /* built in */
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 }
|
