00001 #include "mrilib.h"
00002 
00003 /*==============================================================================*/
00004 /*========== The following functions are lifted from afni_fimfunc.c ============*/
00005 /*------------------------------------------------------------------------------
00006    Rank-order a float array, with ties getting the average rank.
00007    The output overwrites the input.
00008 --------------------------------------------------------------------------------*/
00009 
00010 static void rank_order_float( int n , float * a )
00011 {
00012    register int ii , ns , n1 , ib ;
00013    static int    nb = 0 ;
00014    static int *   b = NULL ;  /* workspaces */
00015    static float * c = NULL ;
00016    float cs ;
00017 
00018    /*- handle special cases -*/
00019 
00020    if( a == NULL ){
00021       if( b != NULL ){ free(b); free(c); b=NULL ; c=NULL; nb=0; }  /* free workspaces */
00022       return ;
00023    }
00024 
00025    if( n < 1 ) return ;                     /* meaningless input */
00026    if( n == 1 ){ a[0] = 0.0 ; return ; }    /* only one point!? */
00027 
00028    /*- make workspaces, if needed -*/
00029 
00030    if( nb < n ){
00031       if( b != NULL ){ free(b); free(c); }
00032       b  = (int   *) malloc(sizeof(int  )*n) ;
00033       c  = (float *) malloc(sizeof(float)*n) ;
00034       nb = n ;
00035    }
00036 
00037    for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ;
00038 
00039    /*- sort input, carrying b along -*/
00040 
00041    qsort_floatint( n , a , b ) ;  /* see cs_sort_fi.c */
00042 
00043    /* compute ranks into c[] */
00044 
00045    n1 = n-1 ;
00046    for( ii=0 ; ii < n1 ; ii++ ){
00047       if( a[ii] == a[ii+1] ){                  /* handle ties */
00048          cs = 2*ii+1 ; ns = 2 ; ib=ii ; ii++ ;
00049          while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; }
00050          for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ;
00051       }
00052    }
00053 
00054    for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ;
00055 
00056    return ;
00057 }
00058 
00059 /*---------------------------------------------------------------------------
00060    Rank orders a[], subtracts the mean rank, and returns the sum-of-squares
00061 -----------------------------------------------------------------------------*/
00062 
00063 static float spearman_rank_prepare( int n , float * a )
00064 {
00065    register int ii ;
00066    register float rb , rs ;
00067 
00068    rank_order_float( n , a ) ;
00069 
00070    rb = 0.5*(n-1) ; rs=0.0 ;
00071    for( ii=0 ; ii < n ; ii++ ){
00072       a[ii] -= rb ;
00073       rs    += a[ii]*a[ii] ;
00074    }
00075 
00076    return rs ;
00077 }
00078 
00079 /*------------------------------------------------------------------------------*/
00080 
00081 static float quadrant_corr_prepare( int n , float * a )
00082 {
00083    register int ii ;
00084    register float rb , rs ;
00085 
00086    rank_order_float( n , a ) ;
00087 
00088    rb = 0.5*(n-1) ; rs=0.0 ;
00089    for( ii=0 ; ii < n ; ii++ ){
00090       a[ii] = (a[ii] > rb) ? 1.0
00091                            : (a[ii] < rb) ? -1.0 : 0.0 ;
00092       rs    += a[ii]*a[ii] ;
00093    }
00094 
00095    return rs ;
00096 }
00097 
00098 /*-----------------------------------------------------------------------------
00099     To correlate x[] with r[], do
00100       rv = spearman_rank_prepare(n,r) ;
00101     then
00102       corr = spearman_rank_corr(n,x,rv,r) ;
00103     Note that these 2 routines are destructive (r and x are replaced by ranks)
00104 -------------------------------------------------------------------------------*/
00105 
00106 static float spearman_rank_corr( int n , float * x , float rv , float * r )
00107 {
00108    register int ii ;
00109    register float ss ; float xv ;
00110 
00111    xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00112 
00113    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00114 
00115    return ( ss/sqrt(rv*xv) ) ;
00116 }
00117 
00118 /*------------------------------------------------------------------------------*/
00119 
00120 static float quadrant_corr( int n , float * x , float rv , float * r )
00121 {
00122    register int ii ;
00123    register float ss ; float xv ;
00124 
00125    xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ;
00126 
00127    for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
00128 
00129    return ( ss/sqrt(rv*xv) ) ;
00130 }
00131 
00132 /*==============================================================================*/
00133 
00134 #define SPEARMAN 1
00135 #define QUADRANT 2
00136 
00137 int main( int argc , char *argv[] )
00138 {
00139    THD_3dim_dataset *dset ;
00140    int nopt=1 , method=SPEARMAN , do_range=0 , do_autoclip=0 ;
00141    int nvox , nvals , ii,iv,jj ;
00142    float *medar, *var , rv , *corr , cmed,cmad,cbot,ctop , clip_val=0.0 ;
00143    MRI_IMAGE *tsim , *medim ;
00144    int nkeep , *keep ;
00145    float *mkeep , *tkeep ;
00146 
00147    /*----*/
00148 
00149    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00150       printf("Usage: 3dTqual [options] dataset\n"
00151              "Computes a `quality index' for each sub-brick in a 3D+time dataset.\n"
00152              "The output is a 1D time series with the index for each sub-brick.\n"
00153              "The results are written to stdout.\n"
00154              "\n"
00155              "Note that small values of the index are 'good', indicating that\n"
00156              "the sub-brick is not very different from the norm.  The purpose\n"
00157              "of this program is to provide a crude way of screening FMRI\n"
00158              "time series for sporadic abnormal images, such as might be\n"
00159              "caused by large subject head motion or scanner glitches.\n"
00160              "\n"
00161              "Do not take the results of this program too literally.  It\n"
00162              "is intended as a GUIDE to help you find data problems, and no\n"
00163              "more.  It is not an assurance that the dataset is good, and\n"
00164              "it may indicate problems where nothing is wrong.\n"
00165              "\n"
00166              "Sub-bricks with index values much higher than others should be\n"
00167              "examined for problems.  How you determine what 'much higher' means\n"
00168              "is mostly up to you.  I suggest graphical inspection of the indexes\n"
00169              "(cf. EXAMPLE, infra).  As a guide, the program will print (stderr)\n"
00170              "the median quality index and the range median-3.5*MAD .. median+3.5*MAD\n"
00171              "(MAD=Median Absolute Deviation).  Values well outside this range might\n"
00172              "be considered suspect; if the quality index were normally distributed,\n"
00173              "then values outside this range would occur only about 1%% of the time.\n"
00174              "\n"
00175              "OPTIONS:\n"
00176              "  -spearman = Quality index is 1 minus the Spearman (rank)\n"
00177              "               correlation coefficient of each sub-brick\n"
00178              "               with the median sub-brick.\n"
00179              "               [This is the default method.]\n"
00180              "  -quadrant = Similar to -spearman, but using 1 minus the\n"
00181              "               quadrant correlation coefficient as the\n"
00182              "               quality index.\n"
00183              "\n"
00184              "  -autoclip = Clip off low-intensity regions in the median sub-brick,\n"
00185              "  -automask =  so that the correlation is only computed between\n"
00186              "               high-intensity (presumably brain) voxels.  The\n"
00187              "               intensity level is determined the same way that\n"
00188              "               3dClipLevel works.  This prevents the vast number\n"
00189              "               of nearly 0 voxels outside the brain from biasing\n"
00190              "               the correlation coefficient calculations.\n"
00191 #if 1
00192              "\n"
00193              "  -clip val = Clip off values below 'val' in the median sub-brick.\n"
00194 #endif
00195              "\n"
00196              "  -range    = Print the median-3.5*MAD and median+3.5*MAD values\n"
00197              "               out with EACH quality index, so that they\n"
00198              "               can be plotted (cf. Example, infra).\n"
00199              "     Notes: * These values are printed to stderr in any case.\n"
00200              "            * This is only useful for plotting with 1dplot.\n"
00201              "            * The lower value median-3.5*MAD is never allowed\n"
00202              "                to go below 0.\n"
00203              "\n"
00204              "EXAMPLE:\n"
00205              "   3dTqual -range -automask fred+orig | 1dplot -one -stdin\n"
00206              "will calculate the time series of quality indexes and plot them\n"
00207              "to an X11 window, along with the median+/-3.5*MAD bands.\n"
00208              "\n"
00209              "NOTE: cf. program 3dToutcount for a somewhat different quality check.\n"
00210              "\n"
00211              "-- RWCox - Aug 2001\n"
00212             ) ;
00213       exit(0) ;
00214    }
00215 
00216    mainENTRY("3dTqual main"); machdep(); AFNI_logger("3dTqual",argc,argv);
00217    PRINT_VERSION("3dTqual");
00218 
00219    /*-- option processing --*/
00220 
00221    while( nopt < argc && argv[nopt][0] == '-' ){
00222 
00223       if( strcmp(argv[nopt],"-autoclip") == 0 ||
00224           strcmp(argv[nopt],"-automask") == 0   ){
00225 
00226          do_autoclip = 1 ; nopt++ ; continue ;
00227       }
00228 
00229 #if 1
00230       if( strcmp(argv[nopt],"-clip") == 0 ){
00231          do_autoclip = 0 ;
00232          clip_val = strtod(argv[++nopt],NULL) ;
00233          if( clip_val <= 0.0 ){
00234             fprintf(stderr,"** value after -clip is illegal!\n"); exit(1);
00235          }
00236          nopt++ ;continue ;
00237       }
00238 #endif
00239 
00240       if( strcmp(argv[nopt],"-range") == 0 ){
00241          do_range = 1 ; nopt++ ; continue ;
00242       }
00243 
00244       if( strcmp(argv[nopt],"-spearman") == 0 ){
00245          method = SPEARMAN ; nopt++ ; continue ;
00246       }
00247 
00248       if( strcmp(argv[nopt],"-quadrant") == 0 ){
00249          method = QUADRANT ; nopt++ ; continue ;
00250       }
00251 
00252       fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]); exit(1);
00253    }
00254 
00255    /*-- open dataset --*/
00256 
00257    if( nopt >= argc ){
00258       fprintf(stderr,"*** No dataset on command line!?\n"); exit(1);
00259    }
00260 
00261    dset = THD_open_dataset( argv[nopt] ) ;
00262    if( dset == NULL ){
00263       fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
00264    }
00265    if( DSET_NUM_TIMES(dset) < 2 ){
00266       fprintf(stderr,"*** Input dataset is not 3D+time\n"); exit(1);
00267    }
00268    DSET_load(dset) ;
00269    if( !DSET_LOADED(dset) ){
00270       fprintf(stderr,"*** Can't read dataset bricks\n"); exit(1);
00271    }
00272 
00273    /*-- compute median brick, and order it for correlation --*/
00274 
00275    nvox  = DSET_NVOX(dset) ;
00276    nvals = DSET_NVALS(dset) ;
00277 
00278    medim = THD_median_brick( dset ) ; if( medim == NULL ) exit(1) ;
00279    medar = MRI_FLOAT_PTR(medim) ;
00280 
00281    if( do_autoclip ){
00282       clip_val = THD_cliplevel( medim , 0.5 ) ;
00283       fprintf(stderr,"++ Autoclip value = %g\n",clip_val) ;
00284    }
00285 
00286    if( clip_val > 0.0 ){
00287       nkeep = 0 ;
00288       for( ii=0 ; ii < nvox ; ii++ )
00289          if( medar[ii] >= clip_val ) nkeep++ ;
00290       if( nkeep < 9 ){
00291          fprintf(stderr,"*** %d voxels survive the clip: can't continue\n",nkeep) ;
00292          exit(1) ;
00293       } else if( nkeep == nvox ){ /* no clipping */
00294          fprintf(stderr,"++ All %d voxels survive the clip\n",nvox) ;
00295          mkeep = medar ;
00296       } else {                    /* nkeep < nvox */
00297          keep  = (int *)   malloc( sizeof(int)  *nkeep ) ;
00298          mkeep = (float *) malloc( sizeof(float)*nkeep ) ;
00299          tkeep = (float *) malloc( sizeof(float)*nkeep ) ;
00300          for( jj=ii=0 ; ii < nvox ; ii++ )
00301             if( medar[ii] >= clip_val ){
00302                keep[jj] = ii ; mkeep[jj] = medar[ii] ; jj++ ;
00303             }
00304          mri_free(medim) ;
00305          fprintf(stderr,"++ %d out of %d voxels survive the clip\n",nkeep,nvox) ;
00306       }
00307    } else {
00308       nkeep = nvox ; mkeep = medar ; /* no clipping */
00309    }
00310 
00311    switch( method ){
00312      default:
00313      case SPEARMAN:
00314        rv = spearman_rank_prepare( nkeep , mkeep ) ; break ;
00315 
00316      case QUADRANT:
00317        rv = quadrant_corr_prepare( nkeep , mkeep ) ; break ;
00318    }
00319 
00320    /*-- loop over input bricks, and correlate with median --*/
00321 
00322    corr = (float *) malloc( sizeof(float)*nvals ) ;
00323 
00324    for( iv=0 ; iv < nvals ; iv++ ){
00325 
00326       /*- get sub-brick -*/
00327 
00328       tsim = THD_extract_float_brick( iv , dset ) ;
00329       if( tsim == NULL ) exit(1) ;
00330       var = MRI_FLOAT_PTR(tsim) ;
00331 
00332       if( nkeep < nvox ){
00333          for( jj=0 ; jj < nkeep ; jj++ ) tkeep[jj] = var[keep[jj]] ;
00334       } else {
00335          tkeep = var ;
00336       }
00337 
00338       /*- compute correlation -*/
00339 
00340       switch( method ){
00341         default:
00342         case SPEARMAN:
00343            corr[iv] = 1.0-spearman_rank_corr(nkeep,tkeep,rv,mkeep) ; break ;
00344 
00345         case QUADRANT:
00346            corr[iv] = 1.0-quadrant_corr(nkeep,tkeep,rv,mkeep) ; break ;
00347       }
00348 
00349       mri_free(tsim) ; DSET_unload_one(dset,iv) ;
00350 
00351       if( !do_range ) printf( "%g\n" , corr[iv] ) ;
00352 
00353    } /* end of loop over sub-bricks */
00354 
00355    /*-- now compute median and MAD of corr[] --*/
00356 
00357    qmedmad_float( nvals,corr , &cmed,&cmad ) ;
00358 
00359    cbot = cmed - 3.5*cmad ;
00360    ctop = cmed + 3.5*cmad ; if( cbot < 0.0 ) cbot = 0.0 ;
00361 
00362    /*-- now print results (if not already out) --*/
00363 
00364    if( do_range ){
00365       for( iv=0 ; iv < nvals ; iv++ )
00366          printf( "%g  %g  %g\n", corr[iv] , cbot,ctop ) ;
00367    }
00368 
00369    fprintf(stderr,"++ Median=%g  Median-3.5*MAD=%g  Median+3.5*MAD=%g\n",
00370            cmed,cbot,ctop) ;
00371    exit(0) ;
00372 }