Doxygen Source Code Documentation
3dfractionize.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Functions | |
THD_fvec3 | AFNI_backward_warp_vector (THD_warp *warp, THD_fvec3 old_fv) |
int | main (int argc, char *argv[]) |
Function Documentation
|
Definition at line 546 of file 3dfractionize.c. References THD_linear_mapping::bot, MATVEC_SUB, THD_linear_mapping::mbac, THD_warp::rig_bod, THD_linear_mapping::svec, THD_warp::tal_12, THD_linear_mapping::top, THD_warp::type, THD_affine_warp::warp, THD_talairach_12_warp::warp, WARP_AFFINE_TYPE, WARP_TALAIRACH_12_TYPE, and THD_fvec3::xyz.
00547 { 00548 THD_fvec3 new_fv ; 00549 00550 if( warp == NULL ) return old_fv ; 00551 00552 switch( warp->type ){ 00553 00554 default: new_fv = old_fv ; break ; 00555 00556 case WARP_TALAIRACH_12_TYPE:{ 00557 THD_linear_mapping map ; 00558 int iw ; 00559 00560 /* test if input is in bot..top of each defined map */ 00561 00562 for( iw=0 ; iw < 12 ; iw++ ){ 00563 map = warp->tal_12.warp[iw] ; 00564 00565 if( old_fv.xyz[0] >= map.bot.xyz[0] && 00566 old_fv.xyz[1] >= map.bot.xyz[1] && 00567 old_fv.xyz[2] >= map.bot.xyz[2] && 00568 old_fv.xyz[0] <= map.top.xyz[0] && 00569 old_fv.xyz[1] <= map.top.xyz[1] && 00570 old_fv.xyz[2] <= map.top.xyz[2] ) break ; /* leave loop */ 00571 } 00572 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ; 00573 } 00574 break ; 00575 00576 case WARP_AFFINE_TYPE:{ 00577 THD_linear_mapping map = warp->rig_bod.warp ; 00578 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ; 00579 } 00580 break ; 00581 00582 } 00583 return new_fv ; 00584 } |
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 15 of file 3dfractionize.c. References ADN_brick_fac, ADN_datum_all, ADN_func_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, AFNI_backward_warp_vector(), AFNI_logger(), argc, THD_3dim_dataset::daxes, THD_3dim_dataset::dblk, THD_datablock::diskptr, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_BRIKNAME, DSET_delete, DSET_load, DSET_LOADED, DSET_write, EDIT_BRICK_FACTOR, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, FUNC_FIM_TYPE, THD_diskptr::header_name, i, ISFUNC, LOAD_FVEC3, machdep(), mainENTRY, malloc, MIN, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, PRINT_VERSION, qsort_int(), realloc, strtod(), THD_3dmm_to_dicomm(), THD_dicomm_to_3dmm(), THD_filename_ok(), THD_is_file(), THD_open_dataset(), THD_open_one_dataset(), tross_Copy_History(), tross_Make_History(), UNLOAD_FVEC3, VIEW_ORIGINAL_TYPE, THD_3dim_dataset::view_type, THD_3dim_dataset::warp, x2, THD_dataxes::xxdel, THD_dataxes::xxorg, y1, THD_dataxes::yydel, THD_dataxes::yyorg, z1, THD_dataxes::zzdel, and THD_dataxes::zzorg.
00016 { 00017 char *prefix = "fractionize" ; 00018 THD_3dim_dataset *tset=NULL , *iset=NULL , *dset=NULL , *wset=NULL ; 00019 int iarg=1 ; 00020 int nxin,nyin,nzin , nxyin ; 00021 float dxin,dyin,dzin , xorgin,yorgin,zorgin , clip=0.0 ; 00022 int nxout,nyout,nzout , nxyout , nvoxout ; 00023 float dxout,dyout,dzout , xorgout,yorgout,zorgout ; 00024 float f1,f2,f , g1,g2,g , h1,h2,h , sx,sy,sz , tx,ty,tz ; 00025 float x1,x2 , y1,y2 , z1,z2 , xx1,xx2 , yy1,yy2 , zz1,zz2 ; 00026 int i,ip , j,jp , k,kp , iv,jv,kv , vtype , ijk ; 00027 float *voxout ; 00028 byte *bin = NULL ; 00029 short *sin = NULL , sclip=0 ; 00030 float *fin = NULL , fclip=0.00001 ; 00031 void *voxin = NULL ; 00032 THD_fvec3 vv ; 00033 00034 THD_warp *warp=NULL ; /* 14 Oct 1999 */ 00035 00036 int do_vote=0 ; /* 18 Oct 1999 */ 00037 int *vote_val = NULL ; 00038 int nvote_val , ivote , voter , vote_print=0 ; 00039 byte *vote_bout = NULL ; 00040 short *vote_sout = NULL ; 00041 float *vote_best = NULL ; 00042 00043 /*-- help? --*/ 00044 00045 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ 00046 printf("Usage: 3dfractionize [options]\n" 00047 "\n" 00048 "* For each voxel in the output dataset, computes the fraction\n" 00049 " of it that is occupied by nonzero voxels from the input.\n" 00050 "* The fraction is stored as a short in the range 0..10000,\n" 00051 " indicating fractions running from 0..1.\n" 00052 "* The template dataset is used only to define the output grid;\n" 00053 " its brick(s) will not be read into memory. (The same is\n" 00054 " true of the warp dataset, if it is used.)\n" 00055 "* The actual values stored in the input dataset are irrelevant,\n" 00056 " except in that they are zero or nonzero (UNLESS the -preserve\n" 00057 " option is used).\n" 00058 "\n" 00059 "The purpose of this program is to allow the resampling of a mask\n" 00060 "dataset (the input) from a fine grid to a coarse grid (defined by\n" 00061 "the template). When you are using the output, you will probably\n" 00062 "want to threshold the mask so that voxels with a tiny occupancy\n" 00063 "fraction aren't used. This can be done in 3dmaskave, by using\n" 00064 "3calc, or with the '-clip' option below.\n" 00065 "\n" 00066 "Options are [the first 2 are 'mandatory options']:\n" 00067 " -template tset = Use dataset 'tset' as a template for the output.\n" 00068 " The output dataset will be on the same grid as\n" 00069 " this dataset.\n" 00070 "\n" 00071 " -input iset = Use dataset 'iset' for the input.\n" 00072 " Only the sub-brick #0 of the input is used.\n" 00073 " You can use the sub-brick selection technique\n" 00074 " described in '3dcalc -help' to choose the\n" 00075 " desired sub-brick from a multi-brick dataset.\n" 00076 "\n" 00077 " -prefix ppp = Use 'ppp' for the prefix of the output.\n" 00078 " [default prefix = 'fractionize']\n" 00079 "\n" 00080 " -clip fff = Clip off voxels that are less than 'fff' occupied.\n" 00081 " 'fff' can be a number between 0.0 and 1.0, meaning\n" 00082 " the fraction occupied, can be a number between 1.0\n" 00083 " and 100.0, meaning the percent occupied, or can be\n" 00084 " a number between 100.0 and 10000.0, meaning the\n" 00085 " direct output value to use as a clip level.\n" 00086 " ** Some sort of clipping is desirable; otherwise,\n" 00087 " an output voxel that is barely overlapped by a\n" 00088 " single nonzero input voxel will enter the mask.\n" 00089 " [default clip = 0.0]\n" 00090 "\n" 00091 " -warp wset = If this option is used, 'wset' is a dataset that\n" 00092 " provides a transformation (warp) from +orig\n" 00093 " coordinates to the coordinates of 'iset'.\n" 00094 " In this case, the output dataset will be in\n" 00095 " +orig coordinates rather than the coordinates\n" 00096 " of 'iset'. With this option:\n" 00097 " ** 'tset' must be in +orig coordinates\n" 00098 " ** 'iset' must be in +acpc or +tlrc coordinates\n" 00099 " ** 'wset' must be in the same coordinates as 'iset'\n" 00100 "\n" 00101 " -preserve = When this option is used, the program will copy\n" 00102 " or the nonzero values of input voxels to the output\n" 00103 " -vote dataset, rather than create a fractional mask.\n" 00104 " Since each output voxel might be overlapped\n" 00105 " by more than one input voxel, the program 'votes'\n" 00106 " for which input value to preserve. For example,\n" 00107 " if input voxels with value=1 occupy 10%% of an\n" 00108 " output voxel, and inputs with value=2 occupy 20%%\n" 00109 " of the same voxel, then the output value in that\n" 00110 " voxel will be set to 2 (provided that 20%% is >=\n" 00111 " to the clip fraction).\n" 00112 " ** Voting can only be done on short-valued datasets,\n" 00113 " or on byte-valued datasets.\n" 00114 " ** Voting is a relatively time-consuming option,\n" 00115 " since a separate loop is made through the\n" 00116 " input dataset for each distinct value found.\n" 00117 " ** Combining this with the -warp option does NOT\n" 00118 " make a general +trlc to +orig transformer!\n" 00119 " This is because for any value to survive the\n" 00120 " vote, its fraction in the output voxel must be\n" 00121 " >= clip fraction, regardless of other values\n" 00122 " present in the output voxel.\n" 00123 "\n" 00124 "Example usage:\n" 00125 " 3dfractionize -template a+orig -input b+tlrc -warp anat+tlrc -clip 0.2\n" 00126 "\n" 00127 "This program will also work in going from a coarse grid to a fine grid,\n" 00128 "but it isn't clear that this capability has any purpose.\n" 00129 "-- RWCox - February 1999\n" 00130 " - October 1999: added -warp and -preserve options\n" 00131 ) ; 00132 exit(0) ; 00133 } 00134 00135 mainENTRY("3dfractionize main"); machdep(); 00136 AFNI_logger("3dfractionize",argc,argv); PRINT_VERSION("3dfractionize"); 00137 00138 /*-- read command line args --*/ 00139 00140 while( iarg < argc && argv[iarg][0] == '-' ){ 00141 00142 if( strcmp(argv[iarg],"-clip") == 0 ){ 00143 clip = strtod( argv[++iarg] , NULL ) ; 00144 00145 if( clip <= 1.0 ) sclip = (short)(10000.0 * clip + 0.49999) ; 00146 else if( clip <= 100.0 ) sclip = (short)( 100.0 * clip + 0.49999) ; 00147 else if( clip <= 10000.0 ) sclip = (short)(clip) ; 00148 else sclip = -1 ; 00149 00150 if( sclip < 0 || sclip > 10000 ){ 00151 fprintf(stderr,"** Illegal clip value!\n") ; exit(1) ; 00152 } 00153 00154 fclip = 0.0001 * sclip ; 00155 iarg++ ; continue ; 00156 } 00157 00158 if( strcmp(argv[iarg],"-template") == 0 || strcmp(argv[iarg],"-tset") == 0 ){ 00159 if( tset != NULL ){ 00160 fprintf(stderr,"** Can't have more than one -template argument!\n") ; 00161 exit(1) ; 00162 } 00163 tset = THD_open_one_dataset( argv[++iarg] ) ; 00164 if( tset == NULL ){ 00165 fprintf(stderr,"** Can't open template %s\n",argv[iarg]) ; 00166 exit(1) ; 00167 } 00168 iarg++ ; continue ; 00169 } 00170 00171 if( strcmp(argv[iarg],"-input") == 0 || strcmp(argv[iarg],"-iset") == 0 ){ 00172 if( iset != NULL ){ 00173 fprintf(stderr,"** Can't have more than one -input argument!\n") ; 00174 exit(1) ; 00175 } 00176 iset = THD_open_dataset( argv[++iarg] ) ; 00177 if( iset == NULL ){ 00178 fprintf(stderr,"** Can't open input %s\n",argv[iarg]) ; 00179 exit(1) ; 00180 } 00181 iarg++ ; continue ; 00182 } 00183 00184 if( strcmp(argv[iarg],"-warp") == 0 || strcmp(argv[iarg],"-wset") == 0 ){ 00185 if( wset != NULL ){ 00186 fprintf(stderr,"** Can't have more than one -warp argument!\n") ; 00187 exit(1) ; 00188 } 00189 wset = THD_open_dataset( argv[++iarg] ) ; 00190 if( wset == NULL ){ 00191 fprintf(stderr,"** Can't open warp %s\n",argv[iarg]) ; 00192 exit(1) ; 00193 } 00194 iarg++ ; continue ; 00195 } 00196 00197 if( strcmp(argv[iarg],"-prefix") == 0 ){ 00198 prefix = argv[++iarg] ; 00199 iarg++ ; continue ; 00200 } 00201 00202 if( strcmp(argv[iarg],"-preserve") == 0 || strcmp(argv[iarg],"-vote") == 0 ){ 00203 do_vote = 1 ; 00204 iarg++ ; continue ; 00205 } 00206 00207 fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ; 00208 } 00209 00210 /*-- check inputs for sanity --*/ 00211 00212 if( tset == NULL ){ 00213 fprintf(stderr,"** No template dataset?\n") ; exit(1) ; 00214 } 00215 if( iset == NULL ){ 00216 fprintf(stderr,"** No input dataset?\n") ; exit(1) ; 00217 } 00218 if( !THD_filename_ok(prefix) ){ 00219 fprintf(stderr,"** Illegal prefix?\n") ; exit(1) ; 00220 } 00221 00222 /*- 14 Oct 1999: additional checks with the -warp option -*/ 00223 00224 if( wset != NULL ){ 00225 00226 if( tset->view_type != VIEW_ORIGINAL_TYPE ){ 00227 fprintf(stderr,"** Template is not in +orig view - this is illegal!\n"); 00228 exit(1); 00229 } 00230 00231 if( iset->view_type == VIEW_ORIGINAL_TYPE ){ 00232 fprintf(stderr,"** Input is in +orig view - this is illegal!\n"); exit(1); 00233 } 00234 00235 if( wset != NULL && wset->view_type != iset->view_type ){ 00236 fprintf(stderr,"** Warp and Input are not in same view - this is illegal!\n"); 00237 exit(1); 00238 } 00239 00240 warp = wset->warp ; 00241 if( warp == NULL ){ 00242 fprintf(stderr,"** Warp dataset doesn't contain a warp transformation!\n"); 00243 exit(1) ; 00244 } 00245 } 00246 00247 /*- 18 Oct 1999: check input for legality if we are voting --*/ 00248 00249 if( do_vote ){ 00250 vtype = DSET_BRICK_TYPE(iset,0) ; 00251 if( vtype != MRI_short && vtype != MRI_byte ){ 00252 fprintf(stderr,"** -preserve option requires short- or byte-valued input dataset!\n"); 00253 exit(1); 00254 } 00255 } 00256 00257 /*-- start to create output dataset --*/ 00258 00259 dset = EDIT_empty_copy( tset ) ; 00260 00261 tross_Copy_History( iset , dset ) ; /* 14 Oct 1999 */ 00262 tross_Make_History( "3dfractionize" , argc,argv , dset ) ; 00263 00264 EDIT_dset_items( dset , 00265 ADN_prefix , prefix , 00266 ADN_nvals , 1 , 00267 ADN_ntt , 0 , 00268 ADN_brick_fac , NULL , 00269 ADN_datum_all , MRI_short , 00270 ADN_none ) ; 00271 00272 if( ISFUNC(dset) ) 00273 EDIT_dset_items( dset , ADN_func_type,FUNC_FIM_TYPE , ADN_none ) ; 00274 00275 if( THD_is_file(dset->dblk->diskptr->header_name) ){ 00276 fprintf(stderr, 00277 "** Output file %s already exists -- cannot continue!\n", 00278 dset->dblk->diskptr->header_name ) ; 00279 exit(1) ; 00280 } 00281 00282 DSET_delete( tset ) ; /* don't need template any more */ 00283 00284 /* load the input dataset */ 00285 00286 DSET_load( iset ) ; 00287 if( ! DSET_LOADED(iset) ){ 00288 fprintf(stderr,"** Can't read input dataset brick!\n") ; exit(1) ; 00289 } 00290 voxin = DSET_ARRAY(iset,0) ; vtype = DSET_BRICK_TYPE(iset,0) ; 00291 switch( vtype ){ 00292 default: fprintf(stderr,"** Illegal brick type in input dataset!\n"); exit(1); 00293 00294 case MRI_byte: bin = (byte * ) voxin ; break ; 00295 case MRI_short: sin = (short *) voxin ; break ; 00296 case MRI_float: fin = (float *) voxin ; break ; 00297 } 00298 00299 /*-- setup voxel index (iv,jv,kv) to coordinate (x,y,z) calculations --*/ 00300 00301 nxin =iset->daxes->nxx ; nyin =iset->daxes->nyy ; nzin =iset->daxes->nzz ; 00302 dxin =iset->daxes->xxdel; dyin =iset->daxes->yydel; dzin =iset->daxes->zzdel; 00303 xorgin =iset->daxes->xxorg; yorgin =iset->daxes->yyorg; zorgin =iset->daxes->zzorg; 00304 00305 nxout =dset->daxes->nxx ; nyout =dset->daxes->nyy ; nzout =dset->daxes->nzz ; 00306 dxout =dset->daxes->xxdel; dyout =dset->daxes->yydel; dzout =dset->daxes->zzdel; 00307 xorgout=dset->daxes->xxorg; yorgout=dset->daxes->yyorg; zorgout=dset->daxes->zzorg; 00308 00309 /* voxel (i,j,k) center is at (xorg+i*dx,yorg+j*dy,zorg+k*dz) */ 00310 00311 nxyout = nxout * nyout ; nvoxout = nxyout * nzout ; 00312 00313 voxout = (float *) malloc( sizeof(float) * nvoxout ) ; 00314 00315 nxyin = nxin * nyin ; 00316 00317 /*-- if voting, do some setup --*/ 00318 00319 if( do_vote ){ 00320 int nvoxin = nxyin * nzin ; 00321 00322 /* 1: find all distinct nonzero values in dataset */ 00323 00324 vote_val = (int *) malloc( sizeof(int) ) ; 00325 nvote_val = 0 ; 00326 for( iv=0 ; iv < nvoxin ; iv++ ){ 00327 kv = (vtype == MRI_short) ? sin[iv] : bin[iv] ; /* voxel value */ 00328 if( kv == 0 ) continue ; /* skip zeroes */ 00329 for( jv=0 ; jv < nvote_val ; jv++ ) /* find in list */ 00330 if( vote_val[jv] == kv ) break ; 00331 if( jv < nvote_val ) continue ; /* skip if already in list */ 00332 00333 vote_val = (int *) realloc( vote_val , sizeof(int)*(nvote_val+1) ) ; 00334 vote_val[nvote_val++] = kv ; 00335 } 00336 if( nvote_val == 0 ){ 00337 fprintf(stderr,"** Input dataset is all zero!\n") ; exit(1) ; 00338 } 00339 fprintf(stderr,"++ Found %d distinct nonzero values in input.\n",nvote_val) ; 00340 00341 if( nvote_val > 1 ) /* arrange into ascending value */ 00342 qsort_int( nvote_val , vote_val ) ; 00343 00344 /* 2: make a receptacle for the voting results */ 00345 00346 if( vtype == MRI_byte ){ 00347 vote_bout = (byte *) malloc( sizeof(byte) * nvoxout ) ; 00348 memset( vote_bout , 0 , sizeof(byte) * nvoxout ) ; 00349 } else { 00350 vote_sout = (short *) malloc( sizeof(short) * nvoxout ) ; 00351 memset( vote_sout , 0 , sizeof(short) * nvoxout ) ; 00352 } 00353 00354 /* 3: make a place to hold the best fraction yet */ 00355 00356 vote_best = (float *) malloc( sizeof(float) * nvoxout ) ; 00357 for( i=0 ; i < nvoxout ; i++ ) vote_best[i] = 0.0 ; 00358 00359 /* 4: if needed, attach the correct brick factor to the output */ 00360 00361 f = DSET_BRICK_FACTOR(iset,0) ; 00362 if( f != 0.0 && f != 1.0 ) EDIT_BRICK_FACTOR(dset,0,f) ; 00363 00364 /* 5: if input is bytes, make the output be that too */ 00365 00366 if( vtype == MRI_byte ) 00367 EDIT_dset_items( dset , ADN_datum_all , MRI_byte , ADN_none ) ; 00368 } 00369 00370 /*-- loop over values to vote on --*/ 00371 00372 ivote = 0 ; 00373 do{ 00374 00375 if( do_vote ){ 00376 voter = vote_val[ivote] ; /* who gets to vote this time */ 00377 if( ivote%10 == 9 ){ 00378 if( !vote_print ) fprintf(stderr,"++ Voting:") ; 00379 fprintf(stderr,"%d",(ivote/10)%10) ; 00380 vote_print++ ; 00381 } 00382 } 00383 00384 for( i=0 ; i < nvoxout ; i++ ) voxout[i] = 0.0 ; /* fractions */ 00385 00386 /*-- loop over input voxels --*/ 00387 00388 for( kv=0 ; kv < nzin ; kv++ ){ 00389 00390 z1 = zorgin + dzin * (kv-0.5) ; z2 = zorgin + dzin * (kv+0.49999) ; 00391 00392 for( jv=0 ; jv < nyin ; jv++ ){ 00393 00394 y1 = yorgin + dyin * (jv-0.5) ; y2 = yorgin + dyin * (jv+0.49999) ; 00395 00396 for( iv=0 ; iv < nxin ; iv++ ){ 00397 00398 ijk = iv + jv*nxin + kv*nxyin ; /* 1D index of voxel (iv,jv,kv) */ 00399 00400 /* decide if we use this voxel, based on its value */ 00401 00402 if( do_vote ){ 00403 if( vtype == MRI_short && sin[ijk] != voter ) continue ; 00404 else if( vtype == MRI_byte && bin[ijk] != voter ) continue ; 00405 } else { 00406 if( vtype == MRI_short && sin[ijk] == 0 ) continue ; 00407 else if( vtype == MRI_byte && bin[ijk] == 0 ) continue ; 00408 else if( vtype == MRI_float && fin[ijk] == 0 ) continue ; 00409 } 00410 00411 x1 = xorgin + dxin * (iv-0.5) ; x2 = xorgin + dxin * (iv+0.49999) ; 00412 00413 /* input voxel (iv,jv,kv) spans coordinates [x1,x2] X [y1,y2] X [z1,z2] */ 00414 00415 /* transform these corner coordinates to output dataset grid coordinates */ 00416 00417 LOAD_FVEC3(vv , x1,y1,z1) ; 00418 vv = THD_3dmm_to_dicomm( iset , vv ); /* transpose from iset to Dicom */ 00419 vv = AFNI_backward_warp_vector( warp , vv ); /* transform to +orig ???*/ 00420 vv = THD_dicomm_to_3dmm( dset , vv ); /* and back from Dicom to dset */ 00421 UNLOAD_FVEC3(vv , xx1,yy1,zz1) ; 00422 00423 LOAD_FVEC3(vv , x2,y2,z2) ; 00424 vv = THD_3dmm_to_dicomm( iset , vv ) ; 00425 vv = AFNI_backward_warp_vector( warp , vv ) ; 00426 vv = THD_dicomm_to_3dmm( dset , vv ) ; 00427 UNLOAD_FVEC3(vv , xx2,yy2,zz2) ; 00428 00429 /* [xx1,xx2] X [yy1,yy2] X [zz1,zz2] is now in coordinates of output dataset */ 00430 00431 /* compute indices into output dataset voxel (keeping fractions) */ 00432 00433 f1 = (xx1-xorgout)/dxout + 0.49999 ; f2 = (xx2-xorgout)/dxout + 0.49999 ; 00434 if( f1 > f2 ){ tx = f1 ; f1 = f2 ; f2 = tx ; } 00435 if( f1 >= nxout || f2 <= 0.0 ) continue ; 00436 if( f1 < 0.0 ) f1 = 0.0 ; if( f2 >= nxout ) f2 = nxout - 0.001 ; 00437 00438 g1 = (yy1-yorgout)/dyout + 0.49999 ; g2 = (yy2-yorgout)/dyout + 0.49999 ; 00439 if( g1 > g2 ){ ty = g1 ; g1 = g2 ; g2 = ty ; } 00440 if( g1 >= nyout || g2 <= 0.0 ) continue ; 00441 if( g1 < 0.0 ) g1 = 0.0 ; if( g2 >= nyout ) g2 = nyout - 0.001 ; 00442 00443 h1 = (zz1-zorgout)/dzout + 0.49999 ; h2 = (zz2-zorgout)/dzout + 0.49999 ; 00444 if( h1 > h2 ){ tz = h1 ; h1 = h2 ; h2 = tz ; } 00445 if( h1 >= nzout || h2 <= 0.0 ) continue ; 00446 if( h1 < 0.0 ) h1 = 0.0 ; if( h2 >= nzout ) h2 = nzout - 0.001 ; 00447 00448 /* input voxel covers voxels [f1,f2] X [g1,g2] X [h1,h2] in the output */ 00449 00450 /* For example, [6.3,7.2] X [9.3,9.6] X [11.7,13.4], which must be */ 00451 /* distributed into these voxels: */ 00452 /* (6,9,11), (7,9,11), (6,9,12), (7,9,12), (6,9,13), and (7,9,13) */ 00453 00454 for( f=f1 ; f < f2 ; f = ip ){ 00455 i = (int) f ; ip = i+1 ; tx = MIN(ip,f2) ; sx = tx - f ; 00456 for( g=g1 ; g < g2 ; g = jp ){ 00457 j = (int) g ; jp = j+1 ; ty = MIN(jp,g2) ; sy = ty - g ; 00458 for( h=h1 ; h < h2 ; h = kp ){ 00459 k = (int) h ; kp = k+1 ; tz = MIN(kp,h2) ; sz = tz - h ; 00460 voxout[ i + j*nxout + k * nxyout ] += sx * sy * sz ; 00461 00462 } 00463 } 00464 } 00465 00466 }}} /* end of loops over voxels */ 00467 00468 /* at this point, voxout[i] contains the fraction that output voxel [i] 00469 is occupied by input voxels that were allowed to vote (or were nonzero) */ 00470 00471 /* if not voting, we are done; otherwise, we must tally the count so far */ 00472 00473 if( do_vote ){ 00474 for( i=0 ; i < nvoxout ; i++ ){ 00475 if( voxout[i] > vote_best[i] && voxout[i] >= fclip ){ /* wins */ 00476 vote_best[i] = voxout[i] ; 00477 if( vtype == MRI_byte ) vote_bout[i] = (byte) voter ; 00478 else if( vtype == MRI_short ) vote_sout[i] = (short) voter ; 00479 } 00480 } 00481 } 00482 00483 /* if we are voting, loop back for the next ballot */ 00484 00485 } while( do_vote && ++ivote < nvote_val ) ; /* end of voting loop */ 00486 00487 /*- throw out some trash -*/ 00488 00489 DSET_delete(iset) ; 00490 if( wset != NULL ) DSET_delete(wset) ; 00491 00492 if( do_vote ){ 00493 free(voxout) ; free(vote_best) ; free(vote_val) ; 00494 if( vote_print ) fprintf(stderr,"\n") ; 00495 } 00496 00497 /*-- if not voting, compute the output brick --*/ 00498 00499 if( ! do_vote ){ 00500 sin = (short *) malloc( sizeof(short) * nvoxout ) ; 00501 if( sin == NULL ){ 00502 fprintf(stderr,"** Can't malloc output brick!\n") ; exit(1) ; 00503 } 00504 00505 for( i=ijk=0 ; i < nvoxout ; i++ ){ 00506 sin[i] = (short) (10000.0 * voxout[i] + 0.49999) ; 00507 if( sin[i] < sclip ) sin[i] = 0 ; 00508 else if( sin[i] > 10000 ) sin[i] = 10000 ; 00509 00510 if( sin[i] != 0 ) ijk++ ; 00511 } 00512 00513 free(voxout) ; 00514 EDIT_substitute_brick( dset , 0 , MRI_short , sin ) ; 00515 00516 fprintf(stderr,"-- Writing %d nonzero mask voxels to dataset %s\n", 00517 ijk , DSET_BRIKNAME(dset) ) ; 00518 00519 /*-- if voting, output brick will be in vote_bout or vote_sout --*/ 00520 00521 } else { 00522 if( vtype == MRI_byte ){ 00523 EDIT_substitute_brick( dset , 0 , MRI_byte , vote_bout ) ; 00524 for( i=ijk=0 ; i < nvoxout ; i++ ) 00525 if( vote_bout[i] != 0 ) ijk++ ; 00526 } else { 00527 EDIT_substitute_brick( dset , 0 , MRI_short , vote_sout ) ; 00528 for( i=ijk=0 ; i < nvoxout ; i++ ) 00529 if( vote_sout[i] != 0 ) ijk++ ; 00530 } 00531 00532 fprintf(stderr,"-- Writing %d nonzero voted voxels to dataset %s\n", 00533 ijk , DSET_BRIKNAME(dset) ) ; 00534 } 00535 00536 DSET_write(dset) ; 00537 exit(0) ; 00538 } |