00001 /*------------------------------------------------------------------------
00002   This program takes the mean of a bunch of datasets, voxel-by-voxel.
00003   Much of this code was taken from 3dcalc.c [I wrote it, I can steal it].
00004   30 Jan 2001 - RWCox
00005   20 Dec 2004 - MSB adds -sd option (cleaned up by RWC a little)
00006 --------------------------------------------------------------------------*/
00007 
00008 #include "mrilib.h"
00009 
00010 int main( int argc , char * argv[] )
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 }