Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
thd_outlier_count.c

00001 #include "mrilib.h"
00002 
00003 /*----------------------------------------------------------------------------------
00004   Inputs: dset (3D+time) and qthr in (0..0.1).
00005   Outputs: *count = array of integer counts of outliers (1 count per sub-brick)
00006            *ctop  = median + 3.5 * MAD of count[]
00007   05 Nov 2001: modified to use THD_extract_array() instead of THD_extract_series()
00008 ------------------------------------------------------------------------------------*/
00009 
00010 void THD_outlier_count( THD_3dim_dataset *dset, float qthr, int **count, int *ctop )
00011 {
00012    int nvals , iv , nxyz , ii , oot , *ccc ;
00013    float alph,fmed,fmad , fbot,ftop ;
00014    MRI_IMAGE * flim ;
00015    float * far , * var , clip_val ;
00016 
00017 ENTRY("THD_outlier_count") ;
00018 
00019    /*-- sanity checks --*/
00020 
00021    if( !ISVALID_DSET(dset) ) EXRETURN ;
00022    if( qthr <= 0.0 || qthr >= 0.1 ) qthr = 0.001 ;
00023 
00024    nvals = DSET_NUM_TIMES(dset) ;
00025    nxyz  = DSET_NVOX(dset) ;
00026    if( nvals < 5 ){ *count = NULL ; *ctop = 0 ; EXRETURN ; }
00027    DSET_load(dset) ;
00028    if( !DSET_LOADED(dset) ){ *count = NULL ; *ctop = 0 ; EXRETURN ; }
00029 
00030    /*-- find clip level [will ignore voxels below this value] --*/
00031 
00032    flim = THD_median_brick( dset ) ;
00033    clip_val = THD_cliplevel( flim , 0.5 ) ;
00034    mri_free(flim) ;
00035 
00036    /*-- setup to count outliers --*/
00037 
00038    alph   = qginv(qthr/nvals) * sqrt(0.5*PI) ;
00039    *count = ccc = (int *) calloc( sizeof(int) , nvals ) ;
00040    var    = (float *) malloc( sizeof(float) * nvals ) ;
00041 
00042    /*--- loop over voxels and count ---*/
00043 
00044    far = (float *) calloc(sizeof(float),nvals+1) ;  /* 05 Nov 2001 */
00045 
00046    for( ii=0 ; ii < nxyz ; ii++ ){
00047 
00048       /*- get time series from voxel #ii -*/
00049 
00050       THD_extract_array( ii , dset , 0 , far ) ;          /* 05 Nov 2001 */
00051       memcpy(var,far,sizeof(float)*nvals ) ;              /* copy it */
00052 
00053       fmed = qmed_float( nvals , far ) ;                  /* median */
00054       if( clip_val > 0.0 && fmed < clip_val ) continue ;  /* below clip? */
00055       for( iv=0 ; iv < nvals ; iv++ )
00056          far[iv] = fabs(far[iv]-fmed) ;
00057       fmad = qmed_float( nvals , far ) ;                  /* MAD */
00058       fbot = fmed - alph*fmad ; ftop = fmed + alph*fmad ; /* inlier range */
00059 
00060       if( fmad > 0.0 ){                                   /* count outliers */
00061          for( iv=0 ; iv < nvals ; iv++ )
00062             if( var[iv] < fbot || var[iv] > ftop ) ccc[iv]++ ;
00063       }
00064 
00065    }
00066 
00067    free(far) ;  /* 05 Nov 2001 */
00068 
00069    for( iv=0 ; iv < nvals ; iv++ ) var[iv] = ccc[iv] ;    /* float-ize counts */
00070    qmedmad_float( nvals,var , &fmed,&fmad ) ; free(var) ; /* median and MAD */
00071    *ctop = (int)(fmed+3.5*fmad+0.499) ;                   /* too much? */
00072 
00073    EXRETURN ;
00074 }
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

thd_outlier_count.c
