00001 /*********************** 3dMax.c **********************************************/
00002 /* Author: Daniel Glen, 26 Apr 2005 */
00003 #include "mrilib.h"
00004 #include "thd_shear3d.h"
00005 
00006 static int datum                   = MRI_float ;
00007 static void Print_Header_MinMax(int Minflag, int Maxflag, THD_3dim_dataset * dset);
00008 static void Max_func(int Minflag, int Maxflag, int Meanflag, int Countflag, int Posflag,\
00009     int Negflag, int Zeroflag, THD_3dim_dataset * dset, byte *mmm, int mmvox);
00010 static void Max_tsfunc( double tzero , double tdelta ,
00011                          int npts , float ts[] , double ts_mean ,
00012                          double ts_slope , void * ud , int nbriks, float * val ) ;
00013 static float minvalue=1E10, maxvalue=-1E10;
00014 /*! compute the overall minimum and maximum voxel values for a dataset */
00015 int main( int argc , char * argv[] )
00016 {
00017    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00018    int nopt, nbriks;
00019    int slow_flag, quick_flag, min_flag, max_flag, mean_flag, automask,count_flag;
00020    int positive_flag, negative_flag, zero_flag;
00021    byte * mmm=NULL ;
00022    int    mmvox=0 ;
00023    int nxyz, i;
00024    MRI_IMAGE *anat_im = NULL;
00025 
00026 
00027    /*----- Read command line -----*/
00028    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00029       printf("Usage: 3dMax [options] dataset\n"
00030              "Compute maximum and/or minimum voxel values of an input dataset\n"
00031              "\n"
00032              "The output is a number to the console.  The input dataset\n"
00033              "may use a sub-brick selection list, as in program 3dcalc.\n"
00034              "Options :\n"
00035              "  -quick = get the information from the header only (default)\n"
00036              "  -slow = read the whole dataset to find the min and max values\n"
00037              "  -min = print the minimum value in dataset\n"
00038              "  -max = print the minimum value in dataset (default)\n"
00039              "  -mean = print the mean value in dataset (implies slow)\n"
00040              "  -count = print the number of voxels included (implies slow)\n"
00041              "  -positive = include only positive voxel values (implies slow)\n"
00042              "  -negative = include only negative voxel values (implies slow)\n"
00043              "  -zero = include only zero voxel values (implies slow)\n"
00044              "  -non-positive = include only voxel values 0 or negative (implies slow)\n"
00045              "  -non-negative = include only voxel values 0 or greater (implies slow)\n"
00046              "  -non-zero = include only voxel values not equal to 0 (implies slow)\n"
00047              "  -mask dset = use dset as mask to include/exclude voxels\n"
00048              "  -automask = automatically compute mask for dataset\n"
00049              "    Can not be combined with -mask\n"
00050              "  -help = print this help screen\n"
00051            ) ;
00052       printf("\n" MASTER_SHORTHELP_STRING ) ;
00053       exit(0) ;
00054    }
00055 
00056    mainENTRY("3dMax main"); machdep(); AFNI_logger("3dMax",argc,argv);
00057 
00058    nopt = 1 ;
00059 
00060    min_flag  = 0;
00061    max_flag = -1;
00062    mean_flag = 0;
00063    slow_flag = 0;
00064    quick_flag = -1;
00065    automask = 0;
00066    count_flag = 0;
00067    positive_flag = -1;
00068    negative_flag = -1;
00069    zero_flag = -1;
00070 
00071    datum = MRI_float;
00072    while( nopt < argc && argv[nopt][0] == '-' ){
00073       if( strcmp(argv[nopt],"-quick") == 0 ){
00074         quick_flag = 1;
00075         nopt++; continue;
00076       }
00077 
00078       if( strcmp(argv[nopt],"-slow") == 0 ){
00079         slow_flag = 1;
00080         nopt++; continue;
00081       }
00082 
00083       if( strcmp(argv[nopt],"-min") == 0 ){
00084         min_flag = 1;
00085         nopt++; continue;
00086       }
00087 
00088       if( strcmp(argv[nopt],"-max") == 0 ){
00089         max_flag = 1;
00090         nopt++; continue;
00091       }
00092 
00093       if( strcmp(argv[nopt],"-mean") == 0 ){
00094         mean_flag = 1;
00095         nopt++; continue;
00096       }
00097 
00098       if( strcmp(argv[nopt],"-count") == 0 ){
00099         count_flag = 1;
00100         nopt++; continue;
00101       }
00102 
00103       if( strcmp(argv[nopt],"-positive") == 0 ){
00104         if(positive_flag!=-1) {
00105           fprintf(stderr, "***Can not use multiple +/-/0 options");
00106           exit(1) ;
00107         }
00108         positive_flag = 1;
00109         negative_flag = 0;
00110         zero_flag = 0;
00111         nopt++; continue;
00112       }
00113 
00114       if( strcmp(argv[nopt],"-negative") == 0 ){
00115         if(positive_flag!=-1) {
00116           fprintf(stderr, "***Can not use multiple +/-/0 options");
00117           exit(1) ;
00118         }
00119         positive_flag = 0;
00120         negative_flag = 1;
00121         zero_flag = 0;
00122         nopt++; continue;
00123       }
00124 
00125       if( strcmp(argv[nopt],"-zero") == 0 ){
00126         if(positive_flag!=-1) {
00127           fprintf(stderr, "***Can not use multiple +/-/0 options");
00128           exit(1) ;
00129         }
00130         positive_flag = 0;
00131         negative_flag = 0;
00132         zero_flag = 1;
00133         nopt++; continue;
00134       }
00135 
00136       if( strcmp(argv[nopt],"-non-positive") == 0 ){
00137         if(positive_flag!=-1) {
00138           fprintf(stderr, "***Can not use multiple +/-/0 options");
00139           exit(1) ;
00140         }
00141         positive_flag = 0;
00142         negative_flag = 1;
00143         zero_flag = 1;
00144         nopt++; continue;
00145       }
00146       if( strcmp(argv[nopt],"-non-negative") == 0 ){
00147         if(positive_flag!=-1) {
00148           fprintf(stderr, "***Can not use multiple +/-/0 options");
00149           exit(1) ;
00150         }
00151         positive_flag = 1;
00152         negative_flag = 0;
00153         zero_flag = 1;
00154         nopt++; continue;
00155       }
00156 
00157       if( strcmp(argv[nopt],"-non-zero") == 0 ){
00158         if(positive_flag!=-1) {
00159           fprintf(stderr, "***Can not use multiple +/-/0 options");
00160           exit(1) ;
00161         }
00162         positive_flag = 1;
00163         negative_flag = 1;
00164         zero_flag = 0;
00165         nopt++; continue;
00166       }
00167 
00168 
00169       if( strcmp(argv[nopt],"-autoclip") == 0 ||
00170           strcmp(argv[nopt],"-automask") == 0   ){
00171 
00172          if( mmm != NULL ){
00173            fprintf(stderr,"** ERROR: can't use -autoclip/mask with -mask!\n");
00174            exit(1) ;
00175          }
00176          automask = 1 ; nopt++ ; continue ;
00177       }
00178 
00179       if( strcmp(argv[nopt],"-mask") == 0 ){
00180          THD_3dim_dataset * mask_dset ;
00181          if( automask ){
00182            fprintf(stderr,"** ERROR: can't use -mask with -automask!\n");
00183            exit(1) ;
00184          }
00185          mask_dset = THD_open_dataset(argv[++nopt]) ;
00186          if( mask_dset == NULL ){
00187             fprintf(stderr,"** ERROR: can't open -mask dataset!\n"); exit(1);
00188          }
00189          if( mmm != NULL ){
00190             fprintf(stderr,"** ERROR: can't have 2 -mask options!\n"); exit(1);
00191          }
00192          mmm = THD_makemask( mask_dset , 0 , 1.0,-1.0 ) ;
00193          mmvox = DSET_NVOX( mask_dset ) ;
00194 
00195          DSET_delete(mask_dset) ; nopt++ ; continue ;
00196       }
00197 
00198       fprintf(stderr, "*** Error - unknown option %s\n", argv[nopt]);
00199       exit(1);
00200    }
00201 
00202    if(((mmm!=NULL) && (quick_flag))||(automask &&quick_flag)) {
00203       if(quick_flag==1)
00204          fprintf(stderr, "+++ Warning - can't have quick option with mask\n");
00205       quick_flag = 0;
00206       slow_flag = 1;
00207    }
00208 
00209    if(max_flag==-1) {                   /* if max_flag is not set by user,*/
00210      if(min_flag || mean_flag ||count_flag)   /* check if other user options set */
00211          max_flag = 0;
00212       else
00213         max_flag = 1;                  /* otherwise check only for max */
00214      }
00215 
00216    if((mean_flag==1)||(count_flag==1)||(positive_flag!=-1))  /* mean flag or count_flag implies slow */
00217      slow_flag = 1;
00218 
00219    /* check slow and quick options */
00220    if((slow_flag)&&(quick_flag!=1))  /* if user asked for slow give it to him */
00221       quick_flag = 0;
00222    else
00223       quick_flag = 1;
00224 
00225    if((max_flag==0)&&(min_flag==0))   /* if the user only asked for mean */
00226      quick_flag = 0;                  /*  no need to do quick way */
00227 
00228    if((quick_flag) && ((positive_flag==1)||(negative_flag==1)||(zero_flag==1)))
00229      fprintf(stderr, "+++ Warning - ignoring +/-/0 flags for quick computations\n");
00230 
00231    if(positive_flag==-1) {   /* if no +/-/0 options set, allow all voxels */
00232      positive_flag = 1;
00233      negative_flag = 1;
00234      zero_flag = 1;
00235    }
00236 
00237    /*----- read input dataset -----*/
00238 
00239    if( nopt >= argc ){
00240       fprintf(stderr,"*** No input dataset!?\n"); exit(1);
00241    }
00242 
00243    old_dset = THD_open_dataset( argv[nopt] ) ;
00244    if( !ISVALID_DSET(old_dset) ){
00245       fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
00246    }
00247 
00248    nxyz = DSET_NVOX(old_dset) ;
00249    if( mmm != NULL && mmvox != nxyz ){
00250       fprintf(stderr,"** Mask and input datasets not the same size!\n") ;
00251       exit(1) ;
00252    }
00253 
00254    if(automask && mmm == NULL ){
00255       mmm = THD_automask( old_dset ) ;
00256       for(i=0;i<nxyz;i++) {
00257         if(mmm[i]!=0) ++mmvox;
00258       }
00259    }
00260 
00261    if(quick_flag)
00262       Print_Header_MinMax(min_flag, max_flag, old_dset);
00263  
00264    if(slow_flag!=1)
00265       exit(0);
00266 
00267    Max_func(min_flag, max_flag, mean_flag,count_flag,positive_flag, negative_flag, zero_flag,\
00268       old_dset, mmm, mmvox);
00269 
00270    if(mmm!=NULL)
00271      free(mmm);
00272    exit(0);
00273 
00274 /* unused code time series method for extracting data */
00275 #if 0
00276    EDIT_dset_items( old_dset ,
00277                     ADN_ntt    , DSET_NVALS(old_dset) ,
00278                     ADN_ttorg  , 0.0 ,
00279                     ADN_ttdel  , 1.0 ,
00280                     ADN_tunits , UNITS_SEC_TYPE ,
00281                     NULL ) ;
00282    nbriks = 1;
00283 
00284    /*------------- ready to compute new min, max -----------*/
00285    new_dset = MAKER_4D_to_typed_fbuc(
00286                  old_dset ,             /* input dataset */
00287                  "temp" ,               /* output prefix */
00288                  datum ,                /* output datum  */
00289                  0 ,                    /* ignore count  */
00290                  0 ,              /* can't detrend in maker function  KRH 12/02*/
00291                  nbriks ,               /* number of briks */
00292                  Max_tsfunc ,         /* timeseries processor */
00293                  NULL                   /* data for tsfunc */
00294               ) ;
00295    if(min_flag)
00296      printf("%-13.6g ", minvalue); 
00297    if(max_flag)
00298      printf("%-13.6g", maxvalue); 
00299    printf("\n");
00300    exit(0) ;
00301 #endif
00302 }
00303 
00304 /*! Print the minimum and maximum values from the header */
00305 static void
00306 Print_Header_MinMax(Minflag, Maxflag, dset)
00307 int Minflag, Maxflag;
00308 THD_3dim_dataset * dset;
00309 {
00310   int ival, nval_per;
00311   float tf; 
00312   double scaledmin, scaledmax, internalmin, internalmax, overallmin, overallmax;
00313 
00314   overallmin = 1E10;
00315   overallmax = -1E10;
00316 
00317   ENTRY("Print_Header_MinMax");
00318    /* print out stuff for each sub-brick */
00319    nval_per = dset->dblk->nvals ;
00320    for( ival=0 ; ival < nval_per ; ival++ ){
00321       tf = DSET_BRICK_FACTOR(dset,ival) ;
00322       if( ISVALID_STATISTIC(dset->stats) ){
00323          if( tf != 0.0 ){
00324            internalmin = dset->stats->bstat[ival].min/tf;
00325            internalmax = dset->stats->bstat[ival].max/tf;
00326           }
00327          scaledmin = dset->stats->bstat[ival].min;
00328          scaledmax = dset->stats->bstat[ival].max;
00329          if( tf != 0.0 ){
00330             if(internalmin < overallmin)
00331                overallmin = scaledmin;
00332             if(internalmax > overallmax)
00333               overallmax = scaledmax;
00334          }
00335          else {
00336             if(scaledmin < overallmin)
00337                overallmin = scaledmin;
00338             if(scaledmax > overallmax)
00339               overallmax = scaledmax;
00340          }
00341       } 
00342       else {
00343          printf( "No valid statistics in header\n") ;
00344          EXRETURN;
00345       }
00346    }
00347 
00348    if(Minflag)
00349       printf("%-13.6g ", overallmin);
00350    if(Maxflag)
00351       printf("%-13.6g", overallmax);
00352    if( tf != 0.0)
00353       printf(" [*%g]\n",tf) ;
00354    else
00355       printf("\n");
00356 
00357    EXRETURN;
00358 }
00359 
00360 
00361 /*! search whole dataset for minimum and maximum */
00362 /* load all at one time */
00363 static void Max_func(Minflag, Maxflag, Meanflag, Countflag, Posflag, Negflag, Zeroflag, dset, mmm, mmvox)
00364 int Minflag, Maxflag;
00365 THD_3dim_dataset * dset;
00366 byte *mmm;  /* mask pointer */
00367 int mmvox;
00368 {
00369    double overallmin, overallmax, overallmean;
00370    double voxval, fac, sum;
00371    int nvox, npts;
00372    int i,k;
00373 
00374    MRI_IMAGE *data_im = NULL;
00375 
00376    ENTRY("Max_func");
00377 
00378    overallmin = 1E10;
00379    overallmax = -1E10;
00380    sum = 0.0;
00381    npts = 0;
00382    DSET_mallocize (dset);
00383    DSET_load (dset);                    /* load dataset */
00384    npts = 0;                            /* keep track of number of points */
00385    for(i=0;i<dset->dblk->nvals; i++) {  /* for each sub-brik in dataset */
00386       data_im = DSET_BRICK (dset, i);   /* set pointer to the 0th sub-brik of the dataset */
00387       fac = DSET_BRICK_FACTOR(dset, i); /* get scale factor for each sub-brik*/
00388       if(fac==0.0) fac=1.0;
00389       if( mmm != NULL)                  /* masked out */
00390         nvox = mmvox;
00391       else
00392         nvox = data_im->nvox;           /* number of voxels in the sub-brik */
00393 
00394       for(k=0;k<nvox;k++) {
00395              if( mmm != NULL && mmm[k] == 0 ) continue ;  /* masked out */
00396 
00397               switch( data_im->kind ){
00398                case MRI_short:{
00399                   short *ar = mri_data_pointer(data_im) ;
00400                   voxval = ar[k];
00401                }
00402                break ;
00403 
00404                case MRI_byte:{
00405                   byte *ar = mri_data_pointer(data_im) ;
00406                   voxval = ar[k];
00407                }
00408                break ;
00409 
00410                case MRI_float:{
00411                   float *ar = mri_data_pointer(data_im) ;
00412                   voxval = ar[k];
00413                }
00414                break ;
00415 
00416               case MRI_double:{
00417                   double *ar = mri_data_pointer(data_im) ;
00418                   voxval = ar[k];
00419                }
00420                break ;
00421 
00422               case MRI_int:{
00423                   int *ar = mri_data_pointer(data_im) ;
00424                   voxval = ar[k];
00425                }
00426                break ;
00427 
00428               case MRI_complex:{
00429                   complex *ar = mri_data_pointer(data_im) ;
00430                   voxval = CABS(ar[k]);
00431                }
00432                break ;
00433 
00434               case MRI_rgb:{
00435                   byte *ar = mri_data_pointer(data_im) ;
00436                   voxval = 0.299*ar[3*k]+0.587*ar[3*k+1]+0.114*ar[3*k+2];
00437                }
00438                break ;
00439 
00440               default:                          /* unknown type */
00441                  voxval = 0.0;                   /* ignore this voxel */
00442                  k = nvox;                       /* skip to next sub-brik */
00443                  fprintf(stderr,"Unknown type, %s, in sub-brik %d\n", MRI_TYPE_name[data_im->kind], i);
00444                break;
00445             }
00446 
00447              if( mmm == NULL || ((mmm!=NULL) && mmm[k] != 0 )){   /* masked in voxel? */
00448               voxval = voxval * fac;             /* apply scale factor */
00449               if(((voxval<0)&&Negflag)||((voxval==0)&&Zeroflag)||((voxval>0)&&Posflag)) {
00450                  sum += voxval;
00451                  ++npts;            
00452                  if(voxval<overallmin)
00453                     overallmin = voxval;
00454                  if(voxval>overallmax)
00455                     overallmax = voxval;
00456               }
00457              }
00458       }
00459    }
00460    if(Minflag)
00461       printf("%-13.6g ", overallmin);
00462    if(Maxflag)
00463       printf("%-13.6g ", overallmax);
00464    if(Meanflag)
00465      {
00466        overallmean = sum/npts;
00467        printf("%-13.6g ", overallmean);
00468      }
00469    if(Countflag)
00470      printf("%-13d", npts);
00471 
00472    printf("\n");
00473 
00474     mri_free (data_im);
00475     /*    DSET_unload_one (dset, 0);*/
00476     EXRETURN;
00477 }
00478 
00479 /* unused code time series method for extracting data */
00480 #if 0
00481 
00482 /*! search whole dataset for minimum and maximum */
00483 static void Max_tsfunc( double tzero, double tdelta ,
00484                           int npts, float ts[],
00485                           double ts_mean, double ts_slope,
00486                           void * ud, int nbriks, float * val          )
00487 {
00488    static int nvox, ncall;
00489    int i;
00490 
00491   ENTRY("Max_tsfunc"); 
00492   /* ts is input vector data */
00493 
00494    /** is this a "notification"? **/
00495    if( val == NULL ){
00496 
00497       if( npts > 0 ){  /* the "start notification" */
00498 
00499          nvox  = npts ;                       /* keep track of   */
00500          ncall = 0 ;                          /* number of calls */
00501 
00502       } else {  /* the "end notification" */
00503 
00504          /* nothing to do here */
00505       }
00506       return ;
00507    }
00508    ncall++;
00509    for(i=0;i<npts;i++) {
00510      if(ts[i]>maxvalue)
00511        maxvalue = ts[i];
00512      if(ts[i]<minvalue)
00513        minvalue = ts[i];
00514    }
00515 
00516   EXRETURN;
00517 }
00518 #endif