00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 static float min_float( int n , float *x )
00012 {
00013 float m=0.0 ; int i ;
00014 if( n < 1 || x == NULL ) return m ;
00015 m = x[0] ;
00016 for( i=1 ; i < n ; i++ ) if( m > x[i] ) m = x[i] ;
00017 return m ;
00018 }
00019
00020
00021
00022 static float max_float( int n , float *x )
00023 {
00024 float m=0.0 ; int i ;
00025 if( n < 1 || x == NULL ) return m ;
00026 m = x[0] ;
00027 for( i=1 ; i < n ; i++ ) if( m < x[i] ) m = x[i] ;
00028 return m ;
00029 }
00030
00031
00032
00033
00034
00035 int main( int argc , char * argv[] )
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 ;
00043 int div = -1 , div_bot,div_top , drange=0;
00044 float data_bot=666.0 , data_top=-666.0 ;
00045 int indump = 0 ;
00046 int pslice=-1 , qslice=-1 , nxy,nz ;
00047 int quiet=0 ;
00048
00049 int medianit = 0 ;
00050 int maxit = 0 ;
00051 int minit = 0 ;
00052 float *exar ;
00053 char *sname = "Average" ;
00054 int self_mask = 0 ;
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
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
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
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 ){
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 ){
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 ){
00264 maxit = 1 ; medianit = 0 ; minit = 0 ;
00265 narg++ ; continue ;
00266 }
00267
00268 if( strncmp(argv[narg],"-min",4) == 0 ){
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" ; }
00278 if( minit ){ sigmait = 0; sname = "Min" ; }
00279
00280
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
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 ;
00299
00300 if( miv > 0 ){
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
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 ){
00426 int kz , ibot ;
00427 mcount = 0 ;
00428 for( kz=0 ; kz < nz ; kz++ ){
00429 ibot = kz*nxy ;
00430
00431 if( kz >= pslice && kz <= qslice ){
00432 for( ii=0 ; ii < nxy ; ii++ )
00433 if( mmm[ii+ibot] ) mcount++ ;
00434 } else {
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
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 ) ;
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 }
