Doxygen Source Code Documentation
3dmaskave.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | INRANGE(i) ( !drange || ( mfac*bar[i] >= data_bot && mfac*bar[i] <= data_top ) ) |
#define | GOOD(i) ( mmm[i] && INRANGE(i) ) |
Functions | |
float | min_float (int n, float *x) |
float | max_float (int n, float *x) |
int | main (int argc, char *argv[]) |
Define Documentation
|
|
|
|
Function Documentation
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 35 of file 3dmaskave.c. References AFNI_logger(), argc, BYTEIZE, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_delete, DSET_index_to_ix, DSET_index_to_jy, DSET_index_to_kz, DSET_load, DSET_NVALS, DSET_NVOX, DSET_NX, DSET_NY, DSET_NZ, DSET_unload, free, machdep(), mainENTRY, malloc, MASTER_SHORTHELP_STRING, max_float(), mc, min_float(), mmm, nz, PRINT_VERSION, qmed_float(), SHORTIZE, SQR, strtod(), and THD_open_dataset().
00036 { 00037 int narg , nvox , ii , mcount , iv , mc ; 00038 THD_3dim_dataset * mask_dset=NULL , * input_dset=NULL ; 00039 float mask_bot=666.0 , mask_top=-666.0 ; 00040 byte * mmm = NULL ; 00041 int dumpit = 0 , sigmait = 0 ; 00042 int miv = 0 ; /* 06 Aug 1998 */ 00043 int div = -1 , div_bot,div_top , drange=0; /* 16 Sep 1998 */ 00044 float data_bot=666.0 , data_top=-666.0 ; 00045 int indump = 0 ; /* 19 Aug 1999 */ 00046 int pslice=-1 , qslice=-1 , nxy,nz ; /* 15 Sep 1999 */ 00047 int quiet=0 ; /* 23 Nov 1999 */ 00048 00049 int medianit = 0 ; /* 06 Jul 2003 */ 00050 int maxit = 0 ; /* 24 Feb 2005 */ 00051 int minit = 0 ; /* 25 Feb 2005 */ 00052 float *exar ; /* 06 Jul 2003 */ 00053 char *sname = "Average" ; /* 06 Jul 2003 */ 00054 int self_mask = 0 ; /* 06 Dec 2004 */ 00055 00056 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ 00057 printf("Usage: 3dmaskave [options] dataset\n" 00058 "Computes average of all voxels in the input dataset\n" 00059 "which satisfy the criterion in the options list.\n" 00060 "If no options are given, then all voxels are included.\n" 00061 "Options:\n" 00062 " -mask mset Means to use the dataset 'mset' as a mask:\n" 00063 " Only voxels with nonzero values in 'mset'\n" 00064 " will be averaged from 'dataset'. Note\n" 00065 " that the mask dataset and the input dataset\n" 00066 " must have the same number of voxels.\n" 00067 " SPECIAL CASE: If 'mset' is the string 'SELF',\n" 00068 " then the input dataset will be\n" 00069 " used to mask itself. That is,\n" 00070 " only nonzero voxels from the\n" 00071 " #miv sub-brick will be used.\n" 00072 " -mindex miv Means to use sub-brick #'miv' from the mask\n" 00073 " dataset. If not given, miv=0.\n" 00074 " -mrange a b Means to further restrict the voxels from\n" 00075 " 'mset' so that only those mask values\n" 00076 " between 'a' and 'b' (inclusive) will\n" 00077 " be used. If this option is not given,\n" 00078 " all nonzero values from 'mset' are used.\n" 00079 " Note that if a voxel is zero in 'mset', then\n" 00080 " it won't be included, even if a < 0 < b.\n" 00081 "\n" 00082 " -dindex div Means to use sub-brick #'div' from the dataset.\n" 00083 " If not given, all sub-bricks will be processed.\n" 00084 " -drange a b Means to only include voxels from the dataset whose\n" 00085 " values fall in the range 'a' to 'b' (inclusive).\n" 00086 " Otherwise, all voxel values are included.\n" 00087 "\n" 00088 " -slices p q Means to only included voxels from the dataset\n" 00089 " whose slice numbers are in the range 'p' to 'q'\n" 00090 " (inclusive). Slice numbers range from 0 to\n" 00091 " NZ-1, where NZ can be determined from the output\n" 00092 " of program 3dinfo. The default is to include\n" 00093 " data from all slices.\n" 00094 " [There is no provision for geometrical voxel]\n" 00095 " [selection except in the slice (z) direction]\n" 00096 "\n" 00097 " -sigma Means to compute the standard deviation as well\n" 00098 " as the mean.\n" 00099 " -median Means to compute the median instead of the mean.\n" 00100 " -max Means to compute the max instead of the mean.\n" 00101 " -min Means to compute the min instead of the mean.\n" 00102 " (-sigma is ignored with -median, -max, or -min)\n" 00103 " -dump Means to print out all the voxel values that\n" 00104 " go into the average.\n" 00105 " -udump Means to print out all the voxel values that\n" 00106 " go into the average, UNSCALED by any internal\n" 00107 " factors.\n" 00108 " N.B.: the scale factors for a sub-brick\n" 00109 " can be found using program 3dinfo.\n" 00110 " -indump Means to print out the voxel indexes (i,j,k) for\n" 00111 " each dumped voxel. Has no effect if -dump\n" 00112 " or -udump is not also used.\n" 00113 " N.B.: if nx,ny,nz are the number of voxels in\n" 00114 " each direction, then the array offset\n" 00115 " in the brick corresponding to (i,j,k)\n" 00116 " is i+j*nx+k*nx*ny.\n" 00117 " -q or\n" 00118 " -quiet Means to print only the minimal results.\n" 00119 " This is useful if you want to create a *.1D file.\n" 00120 "\n" 00121 "The output is printed to stdout (the terminal), and can be\n" 00122 "saved to a file using the usual redirection operation '>'.\n" 00123 ) ; 00124 00125 printf("\n" MASTER_SHORTHELP_STRING ) ; 00126 00127 exit(0) ; 00128 } 00129 00130 mainENTRY("3dmaskave main"); machdep(); AFNI_logger("3dmaskave",argc,argv); 00131 PRINT_VERSION("3dmaskave") ; 00132 00133 /* scan argument list */ 00134 00135 narg = 1 ; 00136 while( narg < argc && argv[narg][0] == '-' ){ 00137 00138 if( strcmp(argv[narg],"-q") == 0 || strcmp(argv[narg],"-quiet") == 0 ){ 00139 quiet++ ; 00140 narg++ ; continue ; 00141 } 00142 00143 if( strncmp(argv[narg],"-mask",5) == 0 ){ 00144 if( mask_dset != NULL || self_mask ){ 00145 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ; 00146 } 00147 if( narg+1 >= argc ){ 00148 fprintf(stderr,"*** -mask option requires a following argument!\n") ; 00149 exit(1) ; 00150 } 00151 narg++ ; 00152 if( strcmp(argv[narg],"SELF") == 0 ){ 00153 self_mask = 1 ; 00154 } else { 00155 mask_dset = THD_open_dataset( argv[narg] ) ; 00156 if( mask_dset == NULL ){ 00157 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ; 00158 } 00159 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){ 00160 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n") ; 00161 exit(1) ; 00162 } else if( DSET_BRICK_TYPE(mask_dset,0) == MRI_rgb ){ 00163 fprintf(stderr,"*** Cannot deal with rgb-valued mask dataset!\n") ; 00164 exit(1) ; 00165 } 00166 } 00167 narg++ ; continue ; 00168 } 00169 00170 /* 06 Aug 1998 */ 00171 00172 if( strncmp(argv[narg],"-mindex",5) == 0 ){ 00173 if( narg+1 >= argc ){ 00174 fprintf(stderr,"*** -mindex option needs 1 following argument!\n") ; exit(1) ; 00175 } 00176 miv = (int) strtod( argv[++narg] , NULL ) ; 00177 if( miv < 0 ){ 00178 fprintf(stderr,"*** -mindex value is negative!\n") ; exit(1) ; 00179 } 00180 narg++ ; continue ; 00181 } 00182 00183 if( strncmp(argv[narg],"-mrange",5) == 0 ){ 00184 if( narg+2 >= argc ){ 00185 fprintf(stderr,"*** -mrange option requires 2 following arguments!\n") ; exit(1) ; 00186 } 00187 mask_bot = strtod( argv[++narg] , NULL ) ; 00188 mask_top = strtod( argv[++narg] , NULL ) ; 00189 if( mask_top < mask_top ){ 00190 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ; 00191 } 00192 narg++ ; continue ; 00193 } 00194 00195 /* 16 Sep 1998 */ 00196 00197 if( strncmp(argv[narg],"-dindex",5) == 0 ){ 00198 if( narg+1 >= argc ){ 00199 fprintf(stderr,"*** -dindex option needs 1 following argument!\n") ; exit(1) ; 00200 } 00201 div = (int) strtod( argv[++narg] , NULL ) ; 00202 if( div < 0 ){ 00203 fprintf(stderr,"*** -dindex value is negative!\n") ; exit(1) ; 00204 } 00205 narg++ ; continue ; 00206 } 00207 00208 if( strncmp(argv[narg],"-drange",5) == 0 ){ 00209 if( narg+2 >= argc ){ 00210 fprintf(stderr, 00211 "*** -drange option requires 2 following arguments!\n") ; 00212 exit(1) ; 00213 } 00214 data_bot = strtod( argv[++narg] , NULL ) ; 00215 data_top = strtod( argv[++narg] , NULL ) ; 00216 if( data_top < data_top ){ 00217 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ; 00218 } 00219 drange = 1 ; 00220 narg++ ; continue ; 00221 } 00222 00223 if( strncmp(argv[narg],"-slices",5) == 0 ){ /* 15 Sep 1999 */ 00224 if( narg+2 >= argc ){ 00225 fprintf(stderr, 00226 "*** -slices option requires 2 following arguments!\n") ; 00227 exit(1) ; 00228 } 00229 pslice = (int) strtod( argv[++narg] , NULL ) ; 00230 qslice = (int) strtod( argv[++narg] , NULL ) ; 00231 if( pslice < 0 || qslice < 0 || qslice < pslice ){ 00232 fprintf(stderr, "*** Illegal values after -slices!\n") ; 00233 exit(1) ; 00234 } 00235 narg++ ; continue ; 00236 } 00237 00238 if( strncmp(argv[narg],"-dump",5) == 0 ){ 00239 dumpit = 1 ; 00240 narg++ ; continue ; 00241 } 00242 00243 if( strncmp(argv[narg],"-udump",5) == 0 ){ 00244 dumpit = 2 ; 00245 narg++ ; continue ; 00246 } 00247 00248 if( strncmp(argv[narg],"-indump",5) == 0 ){ /* 19 Aug 1999 */ 00249 indump = 1 ; 00250 narg++ ; continue ; 00251 } 00252 00253 if( strncmp(argv[narg],"-sigma",5) == 0 ){ 00254 sigmait = 2 ; 00255 narg++ ; continue ; 00256 } 00257 00258 if( strncmp(argv[narg],"-median",5) == 0 ){ 00259 medianit = 1 ; maxit = 0 ; minit = 0 ; 00260 narg++ ; continue ; 00261 } 00262 00263 if( strncmp(argv[narg],"-max",4) == 0 ){ /* 24 Feb 2005 */ 00264 maxit = 1 ; medianit = 0 ; minit = 0 ; 00265 narg++ ; continue ; 00266 } 00267 00268 if( strncmp(argv[narg],"-min",4) == 0 ){ /* 25 Feb 2005 */ 00269 maxit = 0 ; medianit = 0 ; minit = 1 ; 00270 narg++ ; continue ; 00271 } 00272 00273 fprintf(stderr,"*** Unknown option: %s\n",argv[narg]) ; exit(1) ; 00274 } 00275 00276 if( medianit ){ sigmait = 0; sname = "Median"; } 00277 if( maxit ){ sigmait = 0; sname = "Max" ; } /* 24 Feb 2005 */ 00278 if( minit ){ sigmait = 0; sname = "Min" ; } /* 25 Feb 2005 */ 00279 00280 /* should have one more argument */ 00281 00282 if( narg >= argc ){ 00283 fprintf(stderr,"*** No input dataset!?\n") ; exit(1) ; 00284 } 00285 00286 #if 0 00287 if( dumpit && mask_dset == NULL ){ 00288 fprintf(stderr,"*** Can't use dump option without -mask!\n") ; exit(1) ; 00289 } 00290 #endif 00291 00292 /* read input dataset */ 00293 00294 input_dset = THD_open_dataset( argv[narg] ) ; 00295 if( input_dset == NULL ){ 00296 fprintf(stderr,"*** Cannot open input dataset!\n") ; exit(1) ; 00297 } 00298 if( self_mask ) mask_dset = input_dset ; /* 06 Dec 2004 */ 00299 00300 if( miv > 0 ){ /* 06 Aug 1998 */ 00301 if( mask_dset == NULL ){ 00302 fprintf(stderr,"*** -mindex option used without -mask!\n") ; exit(1) ; 00303 } 00304 if( miv >= DSET_NVALS(mask_dset) ){ 00305 fprintf(stderr,"*** -mindex value is too large!\n") ; exit(1) ; 00306 } 00307 } 00308 00309 if( DSET_BRICK_TYPE(input_dset,0) == MRI_complex ){ 00310 fprintf(stderr,"*** Cannot deal with complex-valued input dataset!\n") ; exit(1) ; 00311 } 00312 00313 if( div >= DSET_NVALS(input_dset) ){ 00314 fprintf(stderr,"*** Not enough sub-bricks in dataset for -dindex %d!\n",div) ; 00315 exit(1) ; 00316 } 00317 00318 if( pslice >= 0 ){ 00319 nxy = DSET_NX(input_dset) * DSET_NY(input_dset) ; 00320 nz = DSET_NZ(input_dset) ; 00321 if( qslice >= nz ){ 00322 fprintf(stderr, 00323 "*** There are only %d slices in the input dataset!\n",nz) ; 00324 exit(1) ; 00325 } 00326 00327 if( pslice == 0 && qslice == nz-1 ) 00328 fprintf(stderr,"+++ -slice option says to use all slices!?\n") ; 00329 } 00330 00331 nvox = DSET_NVOX(input_dset) ; 00332 00333 /* make a byte mask from mask dataset */ 00334 00335 mmm = (byte *) malloc( sizeof(byte) * nvox ) ; 00336 if( mmm == NULL ){ 00337 fprintf(stderr,"*** Cannot malloc workspace!\n") ; exit(1) ; 00338 } 00339 00340 if( mask_dset != NULL ){ 00341 if( DSET_NVOX(mask_dset) != nvox ){ 00342 fprintf(stderr,"*** Input and mask datasets are not same dimensions!\n") ; 00343 exit(1) ; 00344 } 00345 DSET_load(mask_dset) ; 00346 if( DSET_ARRAY(mask_dset,miv) == NULL ){ 00347 fprintf(stderr,"*** Cannot read in mask dataset BRIK!\n") ; exit(1) ; 00348 } 00349 00350 switch( DSET_BRICK_TYPE(mask_dset,miv) ){ 00351 default: 00352 fprintf(stderr,"*** Cannot deal with mask dataset datum!\n") ; exit(1) ; 00353 00354 case MRI_short:{ 00355 short mbot , mtop ; 00356 short * mar = (short *) DSET_ARRAY(mask_dset,miv) ; 00357 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ; 00358 if( mfac == 0.0 ) mfac = 1.0 ; 00359 if( mask_bot <= mask_top ){ 00360 mbot = SHORTIZE(mask_bot/mfac) ; 00361 mtop = SHORTIZE(mask_top/mfac) ; 00362 } else { 00363 mbot = (short) -MRI_TYPE_maxval[MRI_short] ; 00364 mtop = (short) MRI_TYPE_maxval[MRI_short] ; 00365 } 00366 for( mcount=0,ii=0 ; ii < nvox ; ii++ ) 00367 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; } 00368 else { mmm[ii] = 0 ; } 00369 } 00370 break ; 00371 00372 case MRI_byte:{ 00373 byte mbot , mtop ; 00374 byte * mar = (byte *) DSET_ARRAY(mask_dset,miv) ; 00375 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ; 00376 if( mfac == 0.0 ) mfac = 1.0 ; 00377 if( mask_bot <= mask_top ){ 00378 mbot = BYTEIZE(mask_bot/mfac) ; 00379 mtop = BYTEIZE(mask_top/mfac) ; 00380 if( mtop == 0 ){ 00381 fprintf(stderr,"*** Illegal mask range for mask dataset of bytes.\n") ; exit(1) ; 00382 } 00383 } else { 00384 mbot = 0 ; 00385 mtop = (byte) MRI_TYPE_maxval[MRI_short] ; 00386 } 00387 for( mcount=0,ii=0 ; ii < nvox ; ii++ ) 00388 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; } 00389 else { mmm[ii] = 0 ; } 00390 } 00391 break ; 00392 00393 case MRI_float:{ 00394 float mbot , mtop ; 00395 float * mar = (float *) DSET_ARRAY(mask_dset,miv) ; 00396 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ; 00397 if( mfac == 0.0 ) mfac = 1.0 ; 00398 if( mask_bot <= mask_top ){ 00399 mbot = (float) (mask_bot/mfac) ; 00400 mtop = (float) (mask_top/mfac) ; 00401 } else { 00402 mbot = -WAY_BIG ; 00403 mtop = WAY_BIG ; 00404 } 00405 for( mcount=0,ii=0 ; ii < nvox ; ii++ ) 00406 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; } 00407 else { mmm[ii] = 0 ; } 00408 } 00409 break ; 00410 } 00411 if( !self_mask ) DSET_unload(mask_dset) ; 00412 00413 if( mcount == 0 ){ 00414 fprintf(stderr,"*** No voxels survive the masking operation\n"); exit(1); 00415 } 00416 00417 fprintf(stderr,"+++ %d voxels survive the mask\n",mcount) ; 00418 00419 } else { 00420 mcount = nvox ; 00421 memset( mmm , 1 , mcount ) ; 00422 fprintf(stderr,"+++ %d voxels in the entire dataset (no mask)\n",mcount) ; 00423 } 00424 00425 if( pslice >= 0 ){ /* 15 Sep 1999 */ 00426 int kz , ibot ; 00427 mcount = 0 ; 00428 for( kz=0 ; kz < nz ; kz++ ){ /* loop over all slices */ 00429 ibot = kz*nxy ; /* base index for this slice */ 00430 00431 if( kz >= pslice && kz <= qslice ){ /* keepers => recount */ 00432 for( ii=0 ; ii < nxy ; ii++ ) 00433 if( mmm[ii+ibot] ) mcount++ ; 00434 } else { /* throw them back */ 00435 for( ii=0 ; ii < nxy ; ii++ ) 00436 mmm[ii+ibot] = 0 ; 00437 } 00438 } 00439 00440 if( mcount == 0 ){ 00441 fprintf(stderr,"*** No voxels survive the slicing operation\n"); exit(1); 00442 } 00443 00444 fprintf(stderr,"+++ %d voxels survive the slicing\n",mcount) ; 00445 } 00446 00447 if( mcount < 2 && sigmait ){ 00448 fprintf(stderr,"+++ [too few voxels; cannot compute sigma]\n") ; 00449 sigmait = 0 ; 00450 } 00451 00452 DSET_load(input_dset) ; 00453 if( DSET_ARRAY(input_dset,0) == NULL ){ 00454 fprintf(stderr,"*** Cannot read in input dataset BRIK!\n") ; exit(1) ; 00455 } 00456 00457 /* loop over input sub-bricks */ 00458 00459 if( div < 0 ){ 00460 div_bot = 0 ; div_top = DSET_NVALS(input_dset) ; 00461 } else { 00462 div_bot = div ; div_top = div+1 ; 00463 } 00464 00465 exar = (float *) malloc( sizeof(float) * mcount ) ; /* 06 Jul 2003 */ 00466 00467 for( iv=div_bot ; iv < div_top ; iv++ ){ 00468 00469 switch( DSET_BRICK_TYPE(input_dset,iv) ){ 00470 00471 default: 00472 printf("*** Illegal sub-brick datum at %d\n",iv) ; 00473 break ; 00474 00475 #define INRANGE(i) ( !drange || ( mfac*bar[i] >= data_bot && mfac*bar[i] <= data_top ) ) 00476 #define GOOD(i) ( mmm[i] && INRANGE(i) ) 00477 00478 case MRI_short:{ 00479 short * bar = (short *) DSET_ARRAY(input_dset,iv) ; 00480 double sum = 0.0 , sigma = 0.0 ; 00481 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ; 00482 if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ; 00483 00484 if( dumpit ){ 00485 int noscal = (dumpit==2) || (mfac==1.0) ; 00486 printf("+++ Dump for sub-brick %d:\n",iv) ; 00487 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){ 00488 if( noscal ) printf(" %d",bar[ii]) ; 00489 else printf(" %g",bar[ii]*mfac) ; 00490 00491 if( indump ) 00492 printf(" (%d,%d,%d)", 00493 DSET_index_to_ix(input_dset,ii), 00494 DSET_index_to_jy(input_dset,ii), 00495 DSET_index_to_kz(input_dset,ii) ) ; 00496 printf("\n") ; 00497 } 00498 } 00499 00500 for( ii=mc=0 ; ii < nvox ; ii++ ) 00501 if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; } 00502 if( mc > 0 ) sum = sum / mc ; 00503 00504 if( sigmait && mc > 1 ){ 00505 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ; 00506 sigma = mfac * sqrt( sigma/(mc-1) ) ; 00507 } 00508 if( medianit ) sum = qmed_float( mc , exar ) ; 00509 else if( maxit ) sum = max_float ( mc , exar ) ; 00510 else if( minit ) sum = min_float ( mc , exar ) ; 00511 sum = mfac * sum ; 00512 00513 if( dumpit ) printf("+++ %s = %g",sname,sum) ; 00514 else printf("%g",sum) ; 00515 00516 if( sigmait ){ 00517 if( dumpit ) printf(" Sigma = %g",sigma) ; 00518 else printf(" %g",sigma) ; 00519 } 00520 if( !quiet ) printf(" [%d voxels]\n",mc) ; 00521 else printf("\n") ; 00522 } 00523 break ; 00524 00525 case MRI_byte:{ 00526 byte * bar = (byte *) DSET_ARRAY(input_dset,iv) ; 00527 double sum = 0.0 , sigma = 0.0 ; 00528 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ; 00529 if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ; 00530 00531 if( dumpit ){ 00532 int noscal = (dumpit==2) || (mfac==1.0) ; 00533 printf("+++ Dump for sub-brick %d:\n",iv) ; 00534 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){ 00535 if( noscal ) printf(" %d",bar[ii]) ; 00536 else printf(" %g",bar[ii]*mfac) ; 00537 00538 if( indump ) 00539 printf(" (%d,%d,%d)", 00540 DSET_index_to_ix(input_dset,ii), 00541 DSET_index_to_jy(input_dset,ii), 00542 DSET_index_to_kz(input_dset,ii) ) ; 00543 printf("\n") ; 00544 } 00545 } 00546 00547 for( ii=mc=0 ; ii < nvox ; ii++ ) 00548 if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; } 00549 if( mc > 0 ) sum = sum / mc ; 00550 00551 if( sigmait && mc > 1 ){ 00552 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ; 00553 sigma = mfac * sqrt( sigma/(mc-1) ) ; 00554 } 00555 if( medianit ) sum = qmed_float( mc , exar ) ; 00556 else if( maxit ) sum = max_float ( mc , exar ) ; 00557 else if( minit ) sum = min_float ( mc , exar ) ; 00558 sum = mfac * sum ; 00559 00560 if( dumpit ) printf("+++ %s = %g",sname,sum) ; 00561 else printf("%g",sum) ; 00562 00563 if( sigmait ){ 00564 if( dumpit ) printf(" Sigma = %g",sigma) ; 00565 else printf(" %g",sigma) ; 00566 } 00567 if( !quiet ) printf(" [%d voxels]\n",mc) ; 00568 else printf("\n") ; 00569 } 00570 break ; 00571 00572 case MRI_float:{ 00573 float * bar = (float *) DSET_ARRAY(input_dset,iv) ; 00574 double sum = 0.0 , sigma = 0.0 ; 00575 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ; 00576 if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ; 00577 00578 if( dumpit ){ 00579 int noscal = (dumpit==2) || (mfac==1.0) ; 00580 printf("+++ Dump for sub-brick %d:\n",iv) ; 00581 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){ 00582 if( noscal ) printf(" %g",bar[ii]) ; 00583 else printf(" %g",bar[ii]*mfac) ; 00584 00585 if( indump ) 00586 printf(" (%d,%d,%d)", 00587 DSET_index_to_ix(input_dset,ii), 00588 DSET_index_to_jy(input_dset,ii), 00589 DSET_index_to_kz(input_dset,ii) ) ; 00590 printf("\n") ; 00591 } 00592 } 00593 00594 for( ii=mc=0 ; ii < nvox ; ii++ ) 00595 if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; } 00596 if( mc > 0 ) sum = sum / mc ; 00597 00598 if( sigmait && mc > 1 ){ 00599 for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ; 00600 sigma = mfac * sqrt( sigma/(mc-1) ) ; 00601 } 00602 if( medianit ) sum = qmed_float( mc , exar ) ; 00603 else if( maxit ) sum = max_float ( mc , exar ) ; 00604 else if( minit ) sum = min_float ( mc , exar ) ; 00605 sum = mfac * sum ; 00606 00607 if( dumpit ) printf("+++ %s = %g",sname,sum) ; 00608 else printf("%g",sum) ; 00609 00610 if( sigmait ){ 00611 if( dumpit ) printf(" Sigma = %g",sigma) ; 00612 else printf(" %g",sigma) ; 00613 } 00614 if( !quiet ) printf(" [%d voxels]\n",mc) ; 00615 else printf("\n") ; 00616 } 00617 break ; 00618 } 00619 } 00620 free(exar) ; DSET_delete(input_dset) ; 00621 exit(0) ; 00622 } |
|
Definition at line 22 of file 3dmaskave.c. References i. Referenced by main().
|
|
Definition at line 11 of file 3dmaskave.c. References i. Referenced by main().
|