00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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 ;
00035 static int HI_nomit = 0 ;
00036 static float * HI_omit = NULL ;
00037 static int HI_notit = 0 ;
00038 static int HI_doall = 0 ;
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)
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 ) ;
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
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
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 }
00255
00256 DSET_unload_one(dset,iv_fim) ;
00257 }
00258
00259
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
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
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
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
00397
00398 fprintf(stderr,"*** unrecognized option %s\a\n",argv[nopt]) ;
00399 exit(-1) ;
00400
00401 }
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 }