00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008   This program generates a histogram for the input AFNI dataset.
00009 
00010   Mod:   2 June 1997
00011          Added option to generate histogram for only those voxels which are
00012          above the operator specified threshold value.
00013 
00014   Mod:   5 Dec 2000, by Vinai Roopchansingh, to include -mask option
00015 
00016   Mod:   16 Feb 2005, RWCox: remove threshold stuff, and add -doall.
00017 */
00018 
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <math.h>
00022 
00023 #include "mrilib.h"
00024 
00025 static EDIT_options HI_edopt ;
00026 
00027 #define NBIN_SPECIAL 65536
00028 #define BIG_NUMBER 9.999e+37
00029 
00030 static int   HI_nopt ;
00031 static int   HI_nbin = 100 ;
00032 static int   HI_log  = 0 ;
00033 
00034 static int     HI_dind  = -1 ;   /* 23 Sep 1998 */
00035 static int     HI_nomit = 0 ;
00036 static float * HI_omit  = NULL ;
00037 static int     HI_notit = 0 ;
00038 static int     HI_doall = 0 ;    /* 16 Feb 2005 */
00039 
00040 static byte  * HI_mask      = NULL ;
00041 static int     HI_mask_nvox = 0 ;
00042 static int     HI_mask_hits = 0 ;
00043 
00044 static double  HI_min = BIG_NUMBER;
00045 static double  HI_max = -BIG_NUMBER;
00046 
00047 #define KEEP(x) ( (HI_nomit==0) ? 1 :  \
00048                   (HI_nomit==1) ? ((x) != HI_omit[0]) : HI_keep(x) )
00049 
00050 void HI_read_opts( int , char ** ) ;
00051 #define HI_syntax(str) \
00052   do{ fprintf(stderr,"\n** ERROR: %s\a\n",str) ; exit(1) ; } while(0)
00053 
00054 /*---------------------------------------------------------------------------------*/
00055 
00056 int HI_keep(float x)  /* check if value x is in the omitted list */
00057 {
00058    register int ii ;
00059    for( ii=0 ; ii < HI_nomit ; ii++ ) if( x == HI_omit[ii] ) return 0 ;
00060    return 1 ;
00061 }
00062 
00063 /*---------------------------------------------------------------------------------*/
00064 
00065 int main( int argc , char * argv[] )
00066 {
00067    int iarg ;
00068    THD_3dim_dataset *dset ;
00069 
00070    int nx,ny,nz , nxyz , ii , kk , nopt , nbin ;
00071    float fbot , ftop ;
00072    int *fbin=NULL , *tbin=NULL ;
00073    float df , dfi ;
00074 
00075    float fimfac;
00076    int iv_fim, fim_type;
00077    byte    *bfim = NULL ;
00078    short   *sfim = NULL ;
00079    float   *ffim = NULL ;
00080    void    *vfim = NULL ;
00081    long cumfbin, cumtbin;
00082    int iv_bot , iv_top ;
00083    float vbot , vtop ;
00084 
00085    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){
00086       printf("Compute histogram of 3D Dataset\n"
00087              "Usage: 3dhistog [editing options] [histogram options] dataset\n"
00088              "\n"
00089              "The editing options are the same as in 3dmerge\n"
00090              " (i.e., the options starting with '-1').\n"
00091              "\n"
00092              "The histogram options are:\n"
00093              "  -nbin #   Means to use '#' bins [default=100]\n"
00094              "            Special Case: for short or byte dataset bricks,\n"
00095              "                          set '#' to zero to have the number\n"
00096              "                          of bins set by the brick range.\n"
00097              "  -dind i   Means to take data from sub-brick #i, rather than #0\n"
00098              "  -omit x   Means to omit the value 'x' from the count;\n"
00099              "              -omit can be used more than once to skip multiple values.\n"
00100              "  -mask m   Means to use dataset 'm' to determine which voxels to use\n"
00101              "  -doall    Means to include all sub-bricks in the calculation;\n"
00102              "              otherwise, only sub-brick #0 (or that from -dind) is used.\n"
00103              "  -notit    Means to leave the title line off the output.\n"
00104              "  -log10    Output log10() of the counts, instead of the count values.\n"
00105              "  -min x    Means specify minimum of histogram.\n"
00106              "  -max x    Means specify maximum of histogram.\n"
00107              "\n"
00108              "The histogram is written to stdout.  Use redirection '>' if you\n"
00109              "want to save it to a file.  The format is a title line, then\n"
00110              "three numbers printed per line:\n"
00111              "  bottom-of-interval  count-in-interval  cumulative-count\n"
00112              "\n"
00113              "-- by RW Cox (V Roopchansingh added the -mask option)\n"
00114          ) ;
00115 
00116       printf("\n" MASTER_SHORTHELP_STRING ) ;
00117 
00118       exit(0) ;
00119    }
00120 
00121    HI_read_opts( argc , argv ) ;
00122    nopt = HI_nopt ;
00123 
00124    if( nopt >=  argc ) HI_syntax("no dset argument?") ;
00125 
00126    fbin = (int *) calloc( sizeof(int) , HI_nbin ) ;
00127    if( fbin == NULL ) HI_syntax("can't allocate histogram array!") ;
00128 
00129    iarg = nopt ;
00130    dset = THD_open_dataset( argv[iarg] ) ;
00131    if( dset == NULL ){
00132      fprintf(stderr,"** ERROR: Can't open dataset %s\n",argv[iarg]) ;
00133      exit(1) ;
00134    }
00135 
00136    if( (HI_mask_nvox > 0) && (HI_mask_nvox != DSET_NVOX(dset)) )
00137      HI_syntax("mask and input dataset bricks don't match in size!") ;
00138 
00139    DSET_mallocize( dset ) ;
00140    DSET_load( dset ) ;
00141    if( !DSET_LOADED(dset) ){
00142      fprintf(stderr,"** ERROR: can't load dataset %s\n",argv[iarg]) ;
00143      exit(1) ;
00144    }
00145    EDIT_one_dataset( dset , &HI_edopt ) ;  /* edit value in memory */
00146 
00147    nx = dset->daxes->nxx ;
00148    ny = dset->daxes->nyy ;
00149    nz = dset->daxes->nzz ; nxyz = nx * ny * nz ;
00150 
00151    if( HI_doall ){
00152      iv_bot = 0 ; iv_top = DSET_NVALS(dset)-1 ;
00153      if( !THD_datum_constant(dset->dblk) ){
00154        fprintf(stderr,
00155          "** ERROR: Dataset %s doesn't have same datum type in all sub-bricks!\n",
00156          argv[iarg]) ;
00157        exit(1) ;
00158      }
00159    } else {
00160      iv_bot = (HI_dind >= 0) ? HI_dind
00161                              : DSET_IS_MASTERED(dset) ? 0
00162                                                       : DSET_PRINCIPAL_VALUE(dset) ;
00163      iv_top = iv_bot ;
00164      if( iv_bot < 0 || iv_bot >= DSET_NVALS(dset) ){
00165        fprintf(stderr,"*** Sub-brick index %d out of range for dataset %s\n",
00166                iv_bot , argv[iarg] ) ;
00167        exit(1) ;
00168      }
00169    }
00170    fim_type = DSET_BRICK_TYPE(dset,iv_bot) ;
00171 
00172    /* find global min and max of data in all used bricks */
00173 
00174    fbot = BIG_NUMBER ; ftop = -fbot ;
00175    for( iv_fim=iv_bot ; iv_fim <= iv_top ; iv_fim++ ){
00176      vbot = mri_min( DSET_BRICK(dset,iv_fim) ) ;
00177      vtop = mri_max( DSET_BRICK(dset,iv_fim) ) ;
00178      fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ; if (fimfac == 0.0)  fimfac = 1.0;
00179      vbot *= fimfac ; vtop *= fimfac ;
00180      if( vbot < fbot ) fbot = vbot;
00181      if( vtop > ftop ) ftop = vtop;
00182    }
00183 
00184    if(HI_min != BIG_NUMBER)
00185      fbot = HI_min;
00186 
00187    if(HI_max != -BIG_NUMBER)
00188      ftop = HI_max;
00189 
00190    if( fbot >= ftop ){
00191      fprintf(stderr,"** ERROR: all values in dataset are = %f!\n",fbot) ;
00192      exit(1) ;
00193    }
00194    switch( fim_type ){
00195       default:
00196         fprintf(stderr,"** ERROR: can't process data of this type!\n") ;
00197       exit(1) ;
00198 
00199       case MRI_byte:
00200       case MRI_short:
00201         nbin = (int)(ftop-fbot+1.0) ;
00202         if( nbin > HI_nbin ) nbin = HI_nbin ;
00203         if( nbin < 2       ) nbin = 2 ;
00204       break ;
00205 
00206       case MRI_float:
00207         nbin = (HI_nbin==NBIN_SPECIAL) ? 100 : HI_nbin ;
00208       break ;
00209    }
00210    df  = (ftop-fbot) / (nbin-1) ;
00211    dfi = 1.0 / df ;
00212 
00213    /* loop over all bricks and accumulate histogram */
00214 
00215    for( iv_fim=iv_bot ; iv_fim <= iv_top ; iv_fim++ ){
00216      fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ;
00217      if (fimfac == 0.0)  fimfac = 1.0;
00218      vfim = DSET_ARRAY(dset,iv_fim) ;
00219 
00220      switch( fim_type ){
00221 
00222        case MRI_short:{
00223          short *fim = (short *)vfim ;
00224          for( ii=0 ; ii < nxyz ; ii++ ){
00225            if( KEEP(fim[ii]*fimfac) && (HI_mask == NULL || HI_mask[ii]) ){
00226              kk = (int)( (fim[ii]-fbot)*dfi ) ;
00227              fbin[kk]++ ;
00228            }
00229          }
00230       }
00231       break ;
00232 
00233        case MRI_byte:{
00234          byte *fim = (byte *)vfim ;
00235          for( ii=0 ; ii < nxyz ; ii++ ){
00236            if( KEEP(fim[ii]*fimfac) && (HI_mask == NULL || HI_mask[ii]) ){
00237              kk = (int)( (fim[ii]-fbot)*dfi ) ;
00238              fbin[kk]++ ;
00239            }
00240          }
00241        }
00242        break ;
00243 
00244        case MRI_float:{
00245          float *fim = (float *)vfim ;
00246          for( ii=0 ; ii < nxyz ; ii++ ){
00247            if( KEEP(fim[ii]*fimfac) && (HI_mask == NULL || HI_mask[ii]) ){
00248              kk = (int)( (fim[ii]-fbot)*dfi ) ;
00249              fbin[kk]++ ;
00250            }
00251          }
00252        }
00253        break ;
00254      }  /* end of switch on data brick type */
00255 
00256      DSET_unload_one(dset,iv_fim) ;
00257    }
00258 
00259    /*** print something ***/
00260 
00261    cumfbin = 0;
00262 
00263    if( HI_log ){
00264      if( ! HI_notit )
00265        printf ("%12s %13s %13s\n", "#Magnitude", "Log_Freq", "Log_Cum_Freq");
00266 
00267      for( kk=0 ; kk < nbin ; kk++ ){
00268        cumfbin += fbin[kk];
00269        printf ("%12.6f %13.6f %13.6f\n",
00270                (fbot+kk*df)*fimfac,
00271                log10((double)fbin[kk]+1.0), log10((double)cumfbin+1.0));
00272      }
00273    } else {
00274      if( ! HI_notit )
00275        printf ("%12s %13s %13s\n",  "#Magnitude", "Freq", "Cum_Freq");
00276 
00277      for( kk=0 ; kk < nbin ; kk++ ){
00278        cumfbin += fbin[kk];
00279        printf ("%12.6f %13d %13ld\n",
00280                (fbot+kk*df)*fimfac, fbin[kk], cumfbin);
00281      }
00282    }
00283    exit(0) ;
00284 }
00285 
00286 
00287 /*--------------------------------------------------------------------
00288    read the arguments, and load the global variables
00289 ----------------------------------------------------------------------*/
00290 
00291 #ifdef HIDEBUG
00292 #  define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
00293 #  define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
00294 #  define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
00295 #else
00296 #  define DUMP1
00297 #  define DUMP2
00298 #  define DUMP3
00299 #endif
00300 
00301 void HI_read_opts( int argc , char * argv[] )
00302 {
00303    int nopt = 1 ;
00304    float val ;
00305    int  ival , kk ;
00306 
00307    INIT_EDOPT( &HI_edopt ) ;
00308 
00309    while( nopt < argc && argv[nopt][0] == '-' ){
00310 
00311       /**** check for editing options ****/
00312 
00313       ival = EDIT_check_argv( argc , argv , nopt , &HI_edopt ) ;
00314       if( ival > 0 ){ nopt += ival ; continue ; }
00315 
00316       if( strncmp(argv[nopt],"-nbin",5) == 0 ){
00317         HI_nbin = strtol( argv[++nopt] , NULL , 10 ) ;
00318         if( HI_nbin < 10 && HI_nbin != 0 ) HI_syntax("illegal value of -nbin!") ;
00319         if( HI_nbin == 0 ) HI_nbin = NBIN_SPECIAL ;
00320         nopt++ ; continue ;
00321       }
00322 
00323       if( strncmp(argv[nopt],"-dind",5) == 0 ){
00324         HI_dind = strtol( argv[++nopt] , NULL , 10 ) ;
00325         if( HI_dind < 0 ) HI_syntax("illegal value of -dind!") ;
00326         nopt++ ; continue ;
00327       }
00328 
00329       if( strncmp(argv[nopt],"-doall",6) == 0 ){
00330         HI_doall = 1 ;
00331         nopt++ ; continue ;
00332       }
00333 
00334       if( strncmp(argv[nopt],"-omit",5) == 0 ){
00335          char *cpt ; float val ;
00336          val = strtod( argv[++nopt] , &cpt ) ;
00337          if( cpt != NULL && *cpt != '\0' ) HI_syntax("illegal value of -omit!") ;
00338          HI_nomit++ ;
00339          if( HI_nomit == 1 )
00340             HI_omit = (float *) malloc( sizeof(float) ) ;
00341          else
00342             HI_omit = (float *) realloc( HI_omit , sizeof(float)*HI_nomit ) ;
00343          HI_omit[HI_nomit-1] = val ;
00344         nopt++ ; continue ;
00345       }
00346 
00347       if( strncmp(argv[nopt],"-notit",5) == 0 ){
00348          HI_notit = 1 ;
00349          nopt++ ; continue ;
00350       }
00351 
00352       if( strncmp(argv[nopt],"-log10",5) == 0 ){
00353          HI_log = 1 ;
00354          nopt++ ; continue ;
00355       }
00356 
00357       /* ----- -mask ----- */
00358 
00359       if( strncmp(argv[nopt],"-mask",5) == 0 )
00360       {
00361           THD_3dim_dataset * mset ; int ii,mc ;
00362           nopt++ ;
00363 
00364           if( nopt >= argc ) HI_syntax("need argument after -mask!") ;
00365 
00366           mset = THD_open_dataset( argv[nopt] ) ;
00367           if( mset == NULL ) HI_syntax("can't open -mask dataset!") ;
00368           HI_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00369           if( HI_mask == NULL ) HI_syntax("can't use -mask dataset!") ;
00370 
00371           HI_mask_nvox = DSET_NVOX(mset) ;
00372           DSET_delete(mset) ;
00373 
00374           for( ii=mc=0 ; ii < HI_mask_nvox ; ii++ ) if( HI_mask[ii] ) mc++ ;
00375 
00376           if( mc == 0 ) HI_syntax("mask is all zeros!") ;
00377 
00378           fprintf(stderr,"++ %d voxels in mask\n",mc) ;
00379           HI_mask_hits = mc ;
00380           nopt++ ; continue ;
00381 
00382       }
00383 
00384       if( strncmp(argv[nopt],"-min",4) == 0 ){
00385          HI_min = strtod( argv[++nopt] , NULL ) ;
00386          nopt++ ; continue ;
00387       }
00388 
00389 
00390       if( strncmp(argv[nopt],"-max",4) == 0 ){
00391          HI_max = strtod( argv[++nopt] , NULL ) ;
00392          nopt++ ; continue ;
00393       }
00394 
00395 
00396       /**** unknown switch ****/
00397 
00398       fprintf(stderr,"*** unrecognized option %s\a\n",argv[nopt]) ;
00399       exit(-1) ;
00400 
00401    }  /* end of loop over options */
00402 
00403 #ifdef HIDEBUG
00404 printf("*** finished with options\n") ;
00405 #endif
00406 
00407    if( HI_doall ){
00408      if( HI_dind >= 0 ){
00409        fprintf(stderr,"** WARNING: -dind ignored with -doall!\n") ;
00410        HI_dind = -1 ;
00411      }
00412    }
00413 
00414    HI_nopt = nopt ;
00415    return ;
00416 }