Doxygen Source Code Documentation
3dMean.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| int | main (int argc, char *argv[]) |
Function Documentation
|
||||||||||||
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 10 of file 3dMean.c. References addto_args(), ADN_datum_all, ADN_none, ADN_prefix, AFNI_logger(), argc, THD_3dim_dataset::dblk, THD_datablock::diskptr, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_BRIKNAME, DSET_delete, DSET_load, DSET_LOADED, DSET_NVALS, DSET_NX, DSET_NY, DSET_NZ, EDIT_coerce_scale_type(), EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, THD_diskptr::header_name, ISVALID_DSET, machdep(), mainENTRY, malloc, MAX, MCW_vol_amax(), mri_datum_size(), nz, THD_is_file(), THD_load_statistics(), THD_open_dataset(), THD_write_3dim_dataset(), tross_Copy_History(), and tross_Make_History().
00011 {
00012 THD_3dim_dataset * inset , * outset ;
00013 int nx,ny,nz,nxyz,nval , ii,kk , nopt=1, nsum=0 ;
00014 char * prefix = "mean" ;
00015 int datum=-1 , verb=0 , do_sd=0, do_sum=0 , do_sqr=0, firstds=0 ;
00016 float ** sum , fsum;
00017 float ** sd;
00018
00019 int fscale=0 , gscale=0 , nscale=0 ;
00020
00021 /*-- help? --*/
00022
00023 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00024 printf("Usage: 3dMean [options] dset dset ...\n"
00025 "Takes the voxel-by-voxel mean of all input datasets;\n"
00026 "the main reason is to be faster than 3dcalc.\n"
00027 "\n"
00028 "Options [see 3dcalc -help for more details on these]:\n"
00029 " -verbose = Print out some information along the way.\n"
00030 " -prefix ppp = Sets the prefix of the output dataset.\n"
00031 " -datum ddd = Sets the datum of the output dataset.\n"
00032 " -fscale = Force scaling of the output to the maximum integer range.\n"
00033 " -gscale = Same as '-fscale', but also forces each output sub-brick to\n"
00034 " to get the same scaling factor.\n"
00035 " -nscale = Don't do any scaling on output to byte or short datasets.\n"
00036 "\n"
00037 " -sd *OR* = Calculate the standard deviation (variance/n-1) instead\n"
00038 " -stdev of the mean (cannot be used with -sqr or -sum).\n"
00039 "\n"
00040 " -sqr = Average the squares, instead of the values.\n"
00041 " -sum = Just take the sum (don't divide by number of datasets).\n"
00042 "\n"
00043 "N.B.: All input datasets must have the same number of voxels along\n"
00044 " each axis (x,y,z,t).\n"
00045 " * At least 2 input datasets are required.\n"
00046 " * Dataset sub-brick selectors [] are allowed.\n"
00047 " * The output dataset origin, time steps, etc., are taken from the\n"
00048 " first input dataset.\n"
00049 ) ;
00050 exit(0) ;
00051 }
00052
00053 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00054
00055 mainENTRY("3dMean main"); machdep() ;
00056
00057 { int new_argc ; char ** new_argv ;
00058 addto_args( argc , argv , &new_argc , &new_argv ) ;
00059 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00060 }
00061
00062 AFNI_logger("3dMean",argc,argv) ;
00063
00064 /*-- command line options --*/
00065
00066 while( nopt < argc && argv[nopt][0] == '-' ){
00067
00068 if( strcmp(argv[nopt],"-sd") == 0 || strcmp(argv[nopt],"-stdev") == 0 ){
00069 do_sd = 1 ; nopt++ ; continue ;
00070 }
00071
00072 if( strcmp(argv[nopt],"-sum") == 0 ){
00073 do_sum = 1 ; nopt++ ; continue ;
00074 }
00075
00076 if( strcmp(argv[nopt],"-sqr") == 0 ){
00077 do_sqr = 1 ; nopt++ ; continue ;
00078 }
00079
00080 if( strcmp(argv[nopt],"-prefix") == 0 ){
00081 if( ++nopt >= argc ){
00082 fprintf(stderr,"** ERROR: need an argument after -prefix!\n"); exit(1);
00083 }
00084 prefix = argv[nopt] ;
00085 nopt++ ; continue ;
00086 }
00087
00088 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00089 verb++ ; nopt++ ; continue ;
00090 }
00091
00092 if( strcmp(argv[nopt],"-datum") == 0 ){
00093 if( ++nopt >= argc ){
00094 fprintf(stderr,"** ERROR: need an argument after -datum!\n"); exit(1);
00095 }
00096 if( strcmp(argv[nopt],"short") == 0 ){
00097 datum = MRI_short ;
00098 } else if( strcmp(argv[nopt],"float") == 0 ){
00099 datum = MRI_float ;
00100 } else if( strcmp(argv[nopt],"byte") == 0 ){
00101 datum = MRI_byte ;
00102 } else {
00103 fprintf(stderr,"** ERROR -datum of type '%s' not supported in 3dMean!\n",
00104 argv[nopt] ) ;
00105 exit(1) ;
00106 }
00107 nopt++ ; continue ; /* go to next arg */
00108 }
00109
00110 if( strncmp(argv[nopt],"-nscale",6) == 0 ){
00111 gscale = fscale = 0 ; nscale = 1 ;
00112 nopt++ ; continue ;
00113 }
00114
00115 if( strncmp(argv[nopt],"-fscale",6) == 0 ){
00116 fscale = 1 ; nscale = 0 ;
00117 nopt++ ; continue ;
00118 }
00119
00120 if( strncmp(argv[nopt],"-gscale",6) == 0 ){
00121 gscale = fscale = 1 ; nscale = 0 ;
00122 nopt++ ; continue ;
00123 }
00124
00125 fprintf(stderr,"** ERROR: unknown option %s\n",argv[nopt]) ;
00126 exit(1) ;
00127 }
00128
00129 /*-- MSB: check for legal combination of options --*/
00130
00131 if( (do_sd) && ( do_sum || do_sqr )){
00132 fprintf(stderr,"** ERROR: -sd cannot be used with -sqr or -sum \n") ;
00133 exit(1) ;
00134 }
00135
00136 /*-- rest of command line should be datasets --*/
00137
00138 if( nopt >= argc-1 ){
00139 fprintf(stderr,"** ERROR: need at least 2 input datasets!\n") ;
00140 exit(1) ;
00141 }
00142
00143 /*-- loop over datasets --*/
00144
00145 firstds = nopt;
00146 for( ; nopt < argc ; nopt++,nsum++ ){
00147
00148 /*-- input dataset header --*/
00149
00150 inset = THD_open_dataset( argv[nopt] ) ;
00151 if( !ISVALID_DSET(inset) ){
00152 fprintf(stderr,"** ERROR: can't open dataset %s\n",argv[nopt]) ;
00153 exit(1) ;
00154 }
00155
00156 /*-- 1st time thru: make workspace and empty output dataset --*/
00157
00158 if( nsum == 0 ){
00159
00160 nx = DSET_NX(inset) ;
00161 ny = DSET_NY(inset) ;
00162 nz = DSET_NZ(inset) ; nxyz= nx*ny*nz;
00163 nval = DSET_NVALS(inset) ;
00164
00165 sum = (float **) malloc( sizeof(float *)*nval ) ; /* array of sub-bricks */
00166 for( kk=0 ; kk < nval ; kk++ ){
00167 sum[kk] = (float *) malloc(sizeof(float)*nxyz) ; /* kk-th sub-brick */
00168 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] = 0.0f ;
00169 }
00170
00171 outset = EDIT_empty_copy( inset ) ;
00172
00173 if( datum < 0 ) datum = DSET_BRICK_TYPE(inset,0) ;
00174
00175 tross_Copy_History( inset , outset ) ;
00176 tross_Make_History( "3dMean" , argc,argv , outset ) ;
00177
00178 EDIT_dset_items( outset ,
00179 ADN_prefix , prefix ,
00180 ADN_datum_all , datum ,
00181 ADN_none ) ;
00182
00183 if( THD_is_file(outset->dblk->diskptr->header_name) ){
00184 fprintf(stderr,
00185 "*** Output file %s already exists -- cannot continue!\n",
00186 outset->dblk->diskptr->header_name ) ;
00187 exit(1) ;
00188 }
00189
00190 } else { /*-- later: check if dataset matches 1st one --*/
00191
00192 if( DSET_NX(inset) != nx ||
00193 DSET_NY(inset) != ny ||
00194 DSET_NZ(inset) != nz ||
00195 DSET_NVALS(inset) != nval ){
00196
00197 fprintf(stderr,"** ERROR: dataset %s doesn't match 1st one in sizes\n",
00198 argv[nopt]) ;
00199 fprintf(stderr,"** I'm telling your mother about this!\n") ;
00200 exit(1) ;
00201 }
00202 }
00203
00204 /*-- read data from disk --*/
00205
00206 DSET_load(inset) ;
00207 if( !DSET_LOADED(inset) ){
00208 fprintf(stderr,"** ERROR: can't read data from dataset %s\n",argv[nopt]) ;
00209 exit(1) ;
00210 }
00211
00212 if( verb ) fprintf(stderr," ++ read in dataset %s\n",argv[nopt]) ;
00213
00214 /*-- sum dataset values --*/
00215
00216 for( kk=0 ; kk < nval ; kk++ ){
00217
00218 if( verb )
00219 fprintf(stderr," + sub-brick %d [%s]\n",
00220 kk,MRI_TYPE_name[DSET_BRICK_TYPE(inset,kk)] ) ;
00221
00222 switch( DSET_BRICK_TYPE(inset,kk) ){
00223 default:
00224 fprintf(stderr,"ERROR: illegal input sub-brick datum\n") ;
00225 exit(1) ;
00226
00227 case MRI_float:{
00228 float *pp = (float *) DSET_ARRAY(inset,kk) ;
00229 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00230 if( fac == 0.0 ) fac = 1.0 ;
00231 if( do_sqr )
00232 for( ii=0 ; ii < nxyz ; ii++ )
00233 { val = fac * pp[ii] ; sum[kk][ii] += val*val ; }
00234 else
00235 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] += fac * pp[ii] ;
00236 }
00237 break ;
00238
00239 case MRI_short:{
00240 short *pp = (short *) DSET_ARRAY(inset,kk) ;
00241 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00242 if( fac == 0.0 ) fac = 1.0 ;
00243 if( do_sqr )
00244 for( ii=0 ; ii < nxyz ; ii++ )
00245 { val = fac * pp[ii] ; sum[kk][ii] += val*val ; }
00246 else
00247 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] += fac * pp[ii] ;
00248 }
00249 break ;
00250
00251 case MRI_byte:{
00252 byte *pp = (byte *) DSET_ARRAY(inset,kk) ;
00253 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00254 if( fac == 0.0 ) fac = 1.0 ;
00255 if( do_sqr )
00256 for( ii=0 ; ii < nxyz ; ii++ )
00257 { val = fac * pp[ii] ; sum[kk][ii] += val*val ; }
00258 else
00259 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] += fac * pp[ii] ;
00260 }
00261 break ;
00262 }
00263 }
00264
00265 DSET_delete(inset) ;
00266 }
00267
00268 /* scale to be mean instead of sum */
00269
00270 if( !do_sum ){
00271 fsum = 1.0 / nsum ;
00272 for( kk=0 ; kk < nval ; kk++ )
00273 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] *= fsum ;
00274 }
00275
00276 /* MSB: If calculating SD,
00277 loop through datasets again,
00278 subtract each value from the mean, square, and sum up */
00279
00280 if( do_sd ){
00281 for (nopt=firstds, nsum=0; nopt < argc ; nopt++,nsum++ ){
00282
00283 /*-- input dataset header --*/
00284
00285 inset = THD_open_dataset( argv[nopt] ) ;
00286 if( !ISVALID_DSET(inset) ){
00287 fprintf(stderr,"** ERROR: can't open dataset %s\n",argv[nopt]) ;
00288 exit(1) ;
00289 }
00290
00291 /*-- 1st time thru: make workspace to hold sd sums */
00292
00293 if( nsum == 0 ){
00294 sd = (float **) malloc( sizeof(float *)*nval ) ; /* array of sub-bricks */
00295 for( kk=0 ; kk < nval ; kk++ ){
00296 sd[kk] = (float *) malloc(sizeof(float)*nxyz) ; /* kk-th sub-brick */
00297 for( ii=0 ; ii < nxyz ; ii++ ) sd[kk][ii] = 0.0f ;
00298 }
00299 }
00300
00301 /*-- read data from disk --*/
00302
00303 DSET_load(inset) ;
00304 if( !DSET_LOADED(inset) ){
00305 fprintf(stderr,"** ERROR: can't read data from dataset %s\n",argv[nopt]) ;
00306 exit(1) ;
00307 }
00308
00309 if( verb ) fprintf(stderr," ++ read in dataset %s\n",argv[nopt]) ;
00310
00311 /*-- sum dataset values into sd --*/
00312
00313 for( kk=0 ; kk < nval ; kk++ ){
00314
00315 if( verb )
00316 fprintf(stderr," + sub-brick %d [%s]\n",
00317 kk,MRI_TYPE_name[DSET_BRICK_TYPE(inset,kk)] ) ;
00318
00319 switch( DSET_BRICK_TYPE(inset,kk) ){
00320 default:
00321 fprintf(stderr,"ERROR: illegal input sub-brick datum\n") ;
00322 exit(1) ;
00323
00324 case MRI_float:{
00325 float *pp = (float *) DSET_ARRAY(inset,kk) ;
00326 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00327 if( fac == 0.0 ) fac = 1.0 ;
00328 for( ii=0 ; ii < nxyz ; ii++ )
00329 { val = (fac * pp[ii]) - sum[kk][ii]; sd[kk][ii] += val*val ; }
00330 }
00331 break ;
00332
00333 case MRI_short:{
00334 short *pp = (short *) DSET_ARRAY(inset,kk) ;
00335 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00336 if( fac == 0.0 ) fac = 1.0 ;
00337 for( ii=0 ; ii < nxyz ; ii++ )
00338 { val = (fac * pp[ii]) - sum[kk][ii]; sd[kk][ii] += val*val ; }
00339 }
00340 break ;
00341
00342 case MRI_byte:{
00343 byte *pp = (byte *) DSET_ARRAY(inset,kk) ;
00344 float fac = DSET_BRICK_FACTOR(inset,kk) , val ;
00345 if( fac == 0.0 ) fac = 1.0 ;
00346 for( ii=0 ; ii < nxyz ; ii++ )
00347 { val = (fac * pp[ii]) - sum[kk][ii]; sd[kk][ii] += val*val ; }
00348 }
00349 break ;
00350 }
00351 }
00352
00353 DSET_delete(inset) ;
00354 } /* end dataset loop */
00355
00356 /*-- fill up output dataset with sd instead of sum --*/
00357
00358 fsum = 1.0 / (nsum - 1) ;
00359 for( kk=0 ; kk < nval ; kk++ ){
00360 for( ii=0 ; ii < nxyz ; ii++ ) sum[kk][ii] = sqrt(sd[kk][ii]*fsum) ;
00361 free((void *)sd[kk]) ;
00362 }
00363
00364 } /* end sd loop */
00365
00366 switch( datum ){
00367
00368 default:
00369 fprintf(stderr,
00370 "*** Fatal Error ***\n"
00371 "*** Somehow ended up with datum = %d\n",datum) ;
00372 exit(1) ;
00373
00374 case MRI_float:{
00375 for( kk=0 ; kk < nval ; kk++ ){
00376 EDIT_substitute_brick(outset, kk, MRI_float, sum[kk]);
00377 DSET_BRICK_FACTOR(outset, kk) = 0.0;
00378 }
00379 }
00380 break ;
00381
00382 case MRI_byte:
00383 case MRI_short:{
00384 void ** dfim ;
00385 float gtop , fimfac , gtemp ;
00386
00387 if( verb )
00388 fprintf(stderr," ++ Scaling output to type %s brick(s)\n",
00389 MRI_TYPE_name[datum] ) ;
00390
00391 dfim = (void **) malloc(sizeof(void *)*nval) ;
00392
00393 if( gscale ){ /* allow global scaling */
00394 gtop = 0.0 ;
00395 for( kk=0 ; kk < nval ; kk++ ){
00396 gtemp = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
00397 gtop = MAX( gtop , gtemp ) ;
00398 if( gtemp == 0.0 )
00399 fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
00400 }
00401 }
00402
00403 for (kk = 0 ; kk < nval ; kk ++ ) {
00404
00405 if( ! gscale ){
00406 gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, sum[kk] ) ;
00407 if( gtop == 0.0 )
00408 fprintf(stderr," -- Warning: output sub-brick %d is all zeros!\n",kk) ;
00409 }
00410
00411 if( fscale ){
00412 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[datum] / gtop : 0.0 ;
00413 } else if( !nscale ){
00414 fimfac = (gtop > MRI_TYPE_maxval[datum] || (gtop > 0.0 && gtop <= 1.0) )
00415 ? MRI_TYPE_maxval[datum]/ gtop : 0.0 ;
00416 } else {
00417 fimfac = 0.0 ;
00418 }
00419
00420 if( verb ){
00421 if( fimfac != 0.0 )
00422 fprintf(stderr," ++ Sub-brick %d scale factor = %f\n",kk,fimfac) ;
00423 else
00424 fprintf(stderr," ++ Sub-brick %d: no scale factor\n" ,kk) ;
00425 }
00426
00427 dfim[kk] = (void *) malloc( mri_datum_size(datum) * nxyz ) ;
00428 if( dfim[kk] == NULL ){ fprintf(stderr,"*** malloc fails at output\n");exit(1); }
00429
00430 EDIT_coerce_scale_type( nxyz , fimfac ,
00431 MRI_float, sum[kk] , datum,dfim[kk] ) ;
00432 free( sum[kk] ) ;
00433 EDIT_substitute_brick(outset, kk, datum, dfim[kk] );
00434
00435 DSET_BRICK_FACTOR(outset,kk) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
00436 }
00437 }
00438 break ;
00439 }
00440
00441 if( verb ) fprintf(stderr," ++ Computing output statistics\n") ;
00442 THD_load_statistics( outset ) ;
00443
00444 THD_write_3dim_dataset( NULL,NULL , outset , True ) ;
00445 if( verb ) fprintf(stderr," ++ Wrote output: %s\n",DSET_BRIKNAME(outset)) ;
00446
00447 exit(0) ;
00448 }
|
