00001 #include "mrilib.h"
00002 
00003 int main( int argc , char * argv[] )
00004 {
00005    THD_3dim_dataset * dset , * oset=NULL ;
00006    int nvals , iv , nxyz , ii,jj , iarg , saveit=0 , oot , ic,cc ;
00007    int * count ;
00008    float qthr=0.001 , alph,fmed,fmad , fbot,ftop,fsig , sq2p,cls ;
00009    MRI_IMAGE * flim ;
00010    float * far , * var ;
00011    byte * mmm=NULL ;
00012    int    mmvox=0 ;
00013    char * prefix=NULL ;
00014    int do_autoclip=0 , npass=0 , do_range=0 ;   /* 12 Aug 2001 */
00015 
00016    int polort=0 , nref ;                        /* 07 Aug 2002 */
00017    float **ref ;
00018    float  *fit ;
00019 
00020    /*----- Read command line -----*/
00021 
00022    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00023       printf("Usage: 3dToutcount [options] dataset\n"
00024              "Calculates number of 'outliers' a 3D+time dataset, at each\n"
00025              "time point, and writes the results to stdout.\n"
00026              "\n"
00027              "Options:\n"
00028              " -mask mset = Only count voxels in the mask dataset.\n"
00029              " -qthr q    = Use 'q' instead of 0.001 in the calculation\n"
00030              "                of alpha (below): 0 < q < 1.\n"
00031              "\n"
00032              " -autoclip }= Clip off 'small' voxels (as in 3dClipLevel);\n"
00033              " -automask }=   you can't use this with -mask!\n"
00034              "\n"
00035              " -range     = Print out median+3.5*MAD of outlier count with\n"
00036              "                each time point; use with 1dplot as in\n"
00037              "                3dToutcount -range fred+orig | 1dplot -stdin -one\n"
00038              " -save ppp  = Make a new dataset, and save the outlier Q in each\n"
00039              "                voxel, where Q is calculated from voxel value v by\n"
00040              "                Q = -log10(qg(abs((v-median)/(sqrt(PI/2)*MAD))))\n"
00041              "             or Q = 0 if v is 'close' to the median (not an outlier).\n"
00042              "                That is, 10**(-Q) is roughly the p-value of value v\n"
00043              "                under the hypothesis that the v's are iid normal.\n"
00044              "              The prefix of the new dataset (float format) is 'ppp'.\n"
00045              "\n"
00046              " -polort nn = Detrend each voxel time series with polynomials of\n"
00047              "                order 'nn' prior to outlier estimation.  Default\n"
00048              "                value of nn=0, which means just remove the median.\n"
00049              "                Detrending is done with L1 regression, not L2.\n"
00050              "\n"
00051              "OUTLIERS are defined as follows:\n"
00052              " * The trend and MAD of each time series are calculated.\n"
00053              "   - MAD = median absolute deviation\n"
00054              "         = median absolute value of time series minus trend.\n"
00055              " * In each time series, points that are 'far away' from the\n"
00056              "    trend are called outliers, where 'far' is defined by\n"
00057              "      alpha * sqrt(PI/2) * MAD\n"
00058              "      alpha = qginv(0.001/N) (inverse of reversed Gaussian CDF)\n"
00059              "      N     = length of time series\n"
00060              " * Some outliers are to be expected, but if a large fraction of the\n"
00061              "    voxels in a volume are called outliers, you should investigate\n"
00062              "    the dataset more fully.\n"
00063              "\n"
00064              "Since the results are written to stdout, you probably want to redirect\n"
00065              "them to a file or another program, as in this example:\n"
00066              "  3dToutcount -automask v1+orig | 1dplot -stdin\n"
00067              "\n"
00068              "NOTE: also see program 3dTqual for a similar quality check.\n"
00069            ) ;
00070       printf("\n" MASTER_SHORTHELP_STRING ) ;
00071       exit(0) ;
00072    }
00073 
00074    mainENTRY("3dToutcount main"); machdep(); AFNI_logger("3dToutcount",argc,argv);
00075    PRINT_VERSION("3dToutcount");
00076 
00077    iarg = 1 ;
00078    while( iarg < argc && argv[iarg][0] == '-' ){
00079 
00080       if( strcmp(argv[iarg],"-autoclip") == 0 ||
00081           strcmp(argv[iarg],"-automask") == 0   ){
00082 
00083          if( mmm != NULL ){
00084            fprintf(stderr,"** ERROR: can't use -autoclip/mask with -mask!\n");
00085            exit(1) ;
00086          }
00087          do_autoclip = 1 ; iarg++ ; continue ;
00088       }
00089 
00090       if( strcmp(argv[iarg],"-range") == 0 ){  /* 12 Aug 2001 */
00091          do_range = 1 ; iarg++ ; continue ;
00092       }
00093 
00094       if( strcmp(argv[iarg],"-save") == 0 ){
00095          prefix = argv[++iarg] ;
00096          if( !THD_filename_ok(prefix) ){
00097             fprintf(stderr,"** ERROR: -save prefix is not good!\n"); exit(1);
00098          }
00099          saveit = 1 ; iarg++ ; continue ;
00100       }
00101 
00102       if( strcmp(argv[iarg],"-qthr") == 0 ){
00103          qthr = strtod( argv[++iarg] , NULL ) ;
00104          if( qthr <= 0.0 || qthr >= 0.999 ){
00105             fprintf(stderr,"** ERROR: -qthr value is illegal!\n"); exit(1);
00106          }
00107          iarg++ ; continue ;
00108       }
00109 
00110       if( strcmp(argv[iarg],"-mask") == 0 ){
00111          int mcount ;
00112          THD_3dim_dataset * mask_dset ;
00113          if( do_autoclip ){
00114            fprintf(stderr,"** ERROR: can't use -mask with -autoclip/mask!\n");
00115            exit(1) ;
00116          }
00117          mask_dset = THD_open_dataset(argv[++iarg]) ;
00118          if( mask_dset == NULL ){
00119             fprintf(stderr,"** ERROR: can't open -mask dataset!\n"); exit(1);
00120          }
00121          if( mmm != NULL ){
00122             fprintf(stderr,"** ERROR: can't have 2 -mask options!\n"); exit(1);
00123          }
00124          mmm = THD_makemask( mask_dset , 0 , 1.0,-1.0 ) ;
00125          mmvox = DSET_NVOX( mask_dset ) ;
00126          mcount = THD_countmask( mmvox , mmm ) ;
00127          fprintf(stderr,"++ %d voxels in the mask\n",mcount) ;
00128          if( mcount <= 5 ){
00129             fprintf(stderr,"** Mask is too small!\n");exit(1);
00130          }
00131          DSET_delete(mask_dset) ; iarg++ ; continue ;
00132       }
00133 
00134       if( strcmp(argv[iarg],"-polort") == 0 ){
00135         polort = strtol( argv[++iarg] , NULL , 10 ) ;
00136         if( polort < 0 || polort > 3){
00137           fprintf(stderr,"** Illegal value of polort!\n"); exit(1);
00138         }
00139         iarg++ ; continue ;
00140       }
00141 
00142       fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00143    }
00144 
00145    /*----- read input dataset -----*/
00146 
00147    if( iarg >= argc ){
00148       fprintf(stderr,"** No input dataset!?\n"); exit(1);
00149    }
00150 
00151    dset = THD_open_dataset( argv[iarg] ) ;
00152    if( !ISVALID_DSET(dset) ){
00153       fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]); exit(1);
00154    }
00155    nvals = DSET_NUM_TIMES(dset) ;
00156    if( nvals < 5 ){
00157       fprintf(stderr,"** Can't use dataset with < 5 time points per voxel!\n") ;
00158       exit(1) ;
00159    }
00160    nxyz = DSET_NVOX(dset) ;
00161    if( mmm != NULL && mmvox != nxyz ){
00162       fprintf(stderr,"** Mask and input datasets not the same size!\n") ;
00163       exit(1) ;
00164    }
00165    DSET_load(dset) ;
00166    if( !DSET_LOADED(dset) ){
00167       fprintf(stderr,"** Can't read dataset from disk!\n") ;
00168       exit(1) ;
00169    }
00170 
00171    /*-- 12 Aug 2001: compute clip, if desired --*/
00172 
00173    if( do_autoclip && mmm == NULL ){
00174       mmm = THD_automask( dset ) ;
00175    }
00176 
00177    /*-- setup to save a new dataset, if desired --*/
00178 
00179    if( saveit ){
00180       float * bar ;
00181       oset = EDIT_empty_copy( dset ) ;
00182       EDIT_dset_items( oset ,
00183                          ADN_prefix    , prefix ,
00184                          ADN_brick_fac , NULL ,
00185                          ADN_datum_all , MRI_float ,
00186                          ADN_func_type , FUNC_FIM_TYPE ,
00187                        ADN_none ) ;
00188 
00189       if( THD_is_file(DSET_HEADNAME(oset)) ){
00190          fprintf(stderr,"** ERROR: -save dataset already exists!\n"); exit(1);
00191       }
00192 
00193       tross_Copy_History( oset , dset ) ;
00194       tross_Make_History( "3dToutcount" , argc , argv , oset ) ;
00195 
00196       for( iv=0 ; iv < nvals ; iv++ ){
00197          EDIT_substitute_brick( oset , iv , MRI_float , NULL ) ;
00198          bar = DSET_ARRAY(oset,iv) ;
00199          for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = 0.0 ;
00200       }
00201    }
00202 
00203    /*-- setup to count --*/
00204 
00205    sq2p  = sqrt(0.5*PI) ;
00206    alph  = qginv(qthr/nvals) * sq2p ;
00207    count = (int *) calloc( sizeof(int) , nvals ) ;
00208    var   = (float *) malloc( sizeof(float) * nvals ) ;
00209 
00210    /* 07 Aug 2002: make polort refs */
00211 
00212    nref = polort+1 ;
00213    ref  = (float **) malloc( sizeof(float *) * nref ) ;
00214    for( jj=0 ; jj < nref ; jj++ )
00215      ref[jj] = (float *) malloc( sizeof(float) * nvals ) ;
00216 
00217    fit = (float *) malloc( sizeof(float) * nref ) ;
00218 
00219    /* r(t) = 1 */
00220 
00221    for( iv=0 ; iv < nvals ; iv++ ) ref[0][iv] = 1.0 ;
00222 
00223    jj = 1 ;
00224    if( polort > 0 ){
00225 
00226      /* r(t) = t - tmid */
00227 
00228      float tm = 0.5 * (nvals-1.0) ; float fac = 2.0 / nvals ;
00229      for( iv=0 ; iv < nvals ; iv++ ) ref[1][iv] = (iv-tm)*fac ;
00230      jj = 2 ;
00231 
00232      /* r(t) = (t-tmid)**jj */
00233 
00234      for( ; jj <= polort ; jj++ )
00235        for( iv=0 ; iv < nvals ; iv++ )
00236          ref[jj][iv] = pow( (iv-tm)*fac , (double)jj ) ;
00237    }
00238 
00239    /*--- loop over voxels and count ---*/
00240 
00241    for( cc=ii=0 ; ii < nxyz ; ii++ ){
00242       if( mmm != NULL && mmm[ii] == 0 ) continue ;  /* masked out */
00243 
00244       npass++ ;
00245 
00246       flim = THD_extract_series( ii , dset , 0 ) ;  /* get data */
00247       far  = MRI_FLOAT_PTR(flim) ;
00248       memcpy(var,far,sizeof(float)*nvals ) ;        /* copy data */
00249 
00250       if( polort == 0 ){                     /* the old way */
00251 
00252         fmed = qmed_float( nvals , far ) ;
00253         for( iv=0 ; iv < nvals ; iv++ ){
00254           var[iv] = var[iv] - fmed ;         /* remove median = resid */
00255           far[iv] = fabs(var[iv]) ;          /* abs value of resid */
00256         }
00257 
00258       } else {                               /* 07 Aug 2002: detrend */
00259 
00260         float val ;
00261         cls = cl1_solve( nvals , nref , far , ref , fit,0 ); /* get fit */
00262         if( cls < 0.0 ) continue ;                          /* bad! should not happen */
00263         for( iv=0 ; iv < nvals ; iv++ ){                    /* detrend */
00264           val = 0.0 ;
00265           for( jj=0 ; jj < nref ; jj++ )                    /* fitted value */
00266             val += fit[jj] * ref[jj][iv] ;
00267 
00268           var[iv] = var[iv]-val ;            /* remove fitted value = resid */
00269           far[iv] = fabs(var[iv]) ;          /* abs value of resid */
00270         }
00271       }
00272 
00273       /* find median of detrended data */
00274 
00275       fmad = qmed_float( nvals , far ) ;
00276       ftop = alph*fmad ; fbot = -ftop ;
00277 
00278       if( fmad > 0.0 ){
00279         if( saveit ) fsig = 1.0/(sq2p*fmad) ;
00280         for( ic=iv=0 ; iv < nvals ; iv++ ){
00281           oot = (var[iv] < fbot || var[iv] > ftop ) ;
00282           if( oot ){ count[iv]++ ; cc++ ; }
00283           if( saveit ){
00284             if( oot ){ far[iv] = -log10qg(fabs(var[iv]*fsig)); ic++; }
00285             else     { far[iv] = 0.0 ; }
00286           }
00287         }
00288 
00289         if( ic > 0 )
00290           THD_insert_series( ii,oset, nvals,MRI_float,far , 0 ) ;
00291       }
00292 
00293       mri_free(flim) ;
00294    }
00295 
00296    if( saveit && cc > 0 ){
00297      DSET_write( oset ) ;
00298      fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(oset)) ;
00299    }
00300 
00301    if( do_range ){
00302       float *ff = (float *)malloc(sizeof(float)*nvals) , cmed,cmad ;
00303       int ctop ;
00304 
00305       for( iv=0 ; iv < nvals ; iv++ ) ff[iv] = count[iv] ;
00306       qmedmad_float( nvals,ff , &cmed,&cmad ) ; free(ff) ;
00307       ctop = (int)(cmed+3.5*cmad+0.499) ;
00308       for( iv=0 ; iv < nvals ; iv++ ) printf("%6d %d\n",count[iv],ctop) ;
00309    } else {
00310       for( iv=0 ; iv < nvals ; iv++ ) printf("%6d\n",count[iv]) ;
00311    }
00312 
00313 #if 0
00314    DSET_unload(dset) ; free(count) ; free(var) ; if(mmm!=NULL)free(mmm) ;
00315 #endif
00316 
00317    if( npass < nxyz )
00318       fprintf(stderr,"++ %d voxels passed mask/clip\n",npass) ;
00319 
00320    exit(0) ;
00321 }