Doxygen Source Code Documentation
thd_outlier_count.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| void | THD_outlier_count (THD_3dim_dataset *dset, float qthr, int **count, int *ctop) |
Function Documentation
|
||||||||||||||||||||
|
Definition at line 10 of file thd_outlier_count.c. References calloc, DSET_load, DSET_LOADED, DSET_NUM_TIMES, DSET_NVOX, ENTRY, far, free, ISVALID_DSET, malloc, mri_free(), qginv(), qmed_float(), qmedmad_float(), THD_cliplevel(), THD_extract_array(), THD_median_brick(), and var. Referenced by T3D_check_outliers().
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 }
|
