1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------------------------
4 Inputs: dset (3D+time) and qthr in (0..0.1).
5 Outputs: *count = array of integer counts of outliers (1 count per sub-brick)
6 *ctop = median + 3.5 * MAD of count[]
7 05 Nov 2001: modified to use THD_extract_array() instead of THD_extract_series()
8 ------------------------------------------------------------------------------------*/
9
THD_outlier_count(THD_3dim_dataset * dset,float qthr,int ** count,int * ctop)10 void THD_outlier_count( THD_3dim_dataset *dset, float qthr, int **count, int *ctop )
11 {
12 int nvals , iv , nxyz , ii , oot , *ccc ;
13 float alph,fmed,fmad , fbot,ftop ;
14 MRI_IMAGE * flim ;
15 float * far , * var , clip_val ;
16
17 ENTRY("THD_outlier_count") ;
18
19 /*-- sanity checks --*/
20
21 if( !ISVALID_DSET(dset) ) EXRETURN ;
22 if( qthr <= 0.0 || qthr >= 0.1 ) qthr = 0.001 ;
23
24 nvals = DSET_NUM_TIMES(dset) ;
25 nxyz = DSET_NVOX(dset) ;
26 if( nvals < 5 ){ *count = NULL ; *ctop = 0 ; EXRETURN ; }
27 DSET_load(dset) ;
28 if( !DSET_LOADED(dset) ){ *count = NULL ; *ctop = 0 ; EXRETURN ; }
29
30 /*-- find clip level [will ignore voxels below this value] --*/
31
32 flim = THD_median_brick( dset ) ;
33 clip_val = THD_cliplevel( flim , 0.5 ) ;
34 mri_free(flim) ;
35
36 /*-- setup to count outliers --*/
37
38 alph = qginv(qthr/nvals) * sqrt(0.5*PI) ;
39 *count = ccc = (int *) calloc( sizeof(int) , nvals ) ;
40 var = (float *) malloc( sizeof(float) * nvals ) ;
41
42 /*--- loop over voxels and count ---*/
43
44 far = (float *) calloc(sizeof(float),nvals+1) ; /* 05 Nov 2001 */
45
46 for( ii=0 ; ii < nxyz ; ii++ ){
47
48 /*- get time series from voxel #ii -*/
49
50 THD_extract_array( ii , dset , 0 , far ) ; /* 05 Nov 2001 */
51 memcpy(var,far,sizeof(float)*nvals ) ; /* copy it */
52
53 fmed = qmed_float( nvals , far ) ; /* median */
54 if( clip_val > 0.0 && fmed < clip_val ) continue ; /* below clip? */
55 for( iv=0 ; iv < nvals ; iv++ )
56 far[iv] = fabs(far[iv]-fmed) ;
57 fmad = qmed_float( nvals , far ) ; /* MAD */
58 fbot = fmed - alph*fmad ; ftop = fmed + alph*fmad ; /* inlier range */
59
60 if( fmad > 0.0 ){ /* count outliers */
61 for( iv=0 ; iv < nvals ; iv++ )
62 if( var[iv] < fbot || var[iv] > ftop ) ccc[iv]++ ;
63 }
64
65 }
66
67 free(far) ; /* 05 Nov 2001 */
68
69 for( iv=0 ; iv < nvals ; iv++ ) var[iv] = ccc[iv] ; /* float-ize counts */
70 qmedmad_float( nvals,var , &fmed,&fmad ) ; free(var) ; /* median and MAD */
71 *ctop = (int)(fmed+3.5*fmad+0.499) ; /* too much? */
72
73 EXRETURN ;
74 }
75