1 #include "mrilib.h"
2
3 /*------------------- 08 Oct 2010: functions for despiking ----------------*/
4
5 #undef SWAP
6 #define SWAP(x,y) (temp=x,x=y,y=temp)
7
8 #undef SORT2
9 #define SORT2(a,b) if(a>b) SWAP(a,b)
10
11 /*--- fast median of 9 values ---*/
12
median9f(float * p)13 static INLINE float median9f(float *p)
14 {
15 register float temp ;
16 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
17 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
18 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
19 SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
20 SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
21 SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
22 SORT2(p[4],p[2]) ; return(p[4]) ;
23 }
24 #undef SORT2
25 #undef SWAP
26
27 /*--- get the local median and MAD of values vec[j-4 .. j+4] ---*/
28
29 #undef mead9
30 #define mead9(j) \
31 { float qqq[9] ; int jj = (j)-4 ; \
32 if( jj < 0 ) jj = 0; else if( jj+8 >= num ) jj = num-9; \
33 qqq[0] = vec[jj+0]; qqq[1] = vec[jj+1]; qqq[2] = vec[jj+2]; \
34 qqq[3] = vec[jj+3]; qqq[4] = vec[jj+4]; qqq[5] = vec[jj+5]; \
35 qqq[6] = vec[jj+6]; qqq[7] = vec[jj+7]; qqq[8] = vec[jj+8]; \
36 med = median9f(qqq); qqq[0] = fabsf(qqq[0]-med); \
37 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
38 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
39 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
40 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
41 mad = median9f(qqq); }
42
43
44 /*-------------------------------------------------------------------------*/
45 /*! Return value is the number of spikes that were squashed
46 (bigger spikes count more).
47 *//*-----------------------------------------------------------------------*/
48
MYY_despike9(int num,float * vec)49 float MYY_despike9( int num , float *vec )
50 {
51 int ii ; float *zma,*zme , med,mad,val , spk=0.0f;
52
53 if( num < 9 || vec == NULL ) return 0 ;
54 zme = (float *)malloc(sizeof(float)*num) ;
55 zma = (float *)malloc(sizeof(float)*num) ;
56
57 for( ii=0 ; ii < num ; ii++ ){
58 mead9(ii) ; zme[ii] = med ; zma[ii] = mad ;
59 }
60 mad = qmed_float(num,zma) ; free(zma) ;
61 if( mad <= 0.0f ){ free(zme); return 0; } /* should not happen */
62 mad *= 7.777f ; /* threshold value */
63
64 for( ii=0 ; ii < num ; ii++ ){
65 val = fabsf(vec[ii]-zme[ii]) / mad ;
66 if( val > 1.0f ){ vec[ii] = zme[ii] ; spk += val ; }
67 }
68
69 free(zme) ; return spk ;
70 }
71 #undef mead9
72
73 /*----------------------------------------------------------------------*/
74 /* Count spikes.
75 Then count steps as spikes in the first difference
76 -- after isolated spikes were squashed.
77 *//*--------------------------------------------------------------------*/
78
THD_countspikes(int num,float * vec)79 float_pair THD_countspikes( int num , float *vec )
80 {
81 float_pair spair = {0.0f,0.0f} ;
82 int n1=num-1 , ii ;
83 float *del , spk,stp ;
84
85 if( num < 11 || vec == NULL ) return spair ;
86
87 spk = MYY_despike9( num , vec ) ;
88
89 del = (float *)malloc(sizeof(float)*n1) ;
90 for( ii=0 ; ii < n1 ; ii++ ) del[ii] = vec[ii+1]-vec[ii] ;
91 stp = MYY_despike9( n1 , del ) ;
92 free(del) ;
93
94 spair.a = spk ; spair.b = stp ; return spair ;
95 }
96
97 /*----------------------------------------------------------------------------*/
98
main(int argc,char * argv[])99 int main( int argc , char * argv[] )
100 {
101 int do_automask=1 ;
102 THD_3dim_dataset *inset=NULL ;
103 byte *mask=NULL ;
104 int mask_nx=0,mask_ny=0,mask_nz=0,nmask , verb=0 , nx,ny,nz,nvox , kk,ntime ;
105 float_pair spair ; float spk=0.0f , stp=0.0f ;
106 MRI_vectim *mrv ;
107 int nopt=1 ;
108
109 /*-- help? --*/
110
111 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
112 printf(
113 "Usage: 3dCountSpikes dataset\n"
114 "\n"
115 "Does this program do something useful?\n"
116 "If it does, that circumstance is completely accidental.\n"
117 ) ;
118 exit(0) ;
119 }
120
121 /*-- startup --*/
122
123 mainENTRY("3dCountSpikes"); machdep();
124
125 /*-- input --*/
126
127 inset = THD_open_dataset(argv[nopt]) ;
128 CHECK_OPEN_ERROR(inset,argv[nopt]) ; nopt++ ;
129 ntime = DSET_NVALS(inset) ;
130 if( ntime < 11 ) ERROR_exit("Input dataset is too short!") ;
131
132 nx = DSET_NX(inset); ny = DSET_NY(inset); nz = DSET_NZ(inset); nvox = nx*ny*nz;
133
134 if( mask != NULL ){
135 if( mask_nx != nx || mask_ny != ny || mask_nz != nz )
136 ERROR_exit("-mask dataset grid doesn't match input dataset") ;
137
138 } else if( do_automask && nvox > nx ){
139 mask = THD_automask( inset ) ;
140 if( mask == NULL )
141 ERROR_message("Can't create -automask from input dataset?") ;
142 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
143 if( verb ) INFO_message("Number of voxels in automask = %d",nmask);
144 if( nmask < 1 ) ERROR_exit("Automask is too small to process") ;
145
146 } else {
147 mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;
148 memset(mask,1,sizeof(byte)*nvox) ;
149 if( verb ) INFO_message("No mask ==> processing all %d voxels",nvox);
150 }
151
152 /* convert input dataset to a vectim, which is more fun */
153
154 mrv = THD_dset_to_vectim( inset , mask , 0 ) ;
155 if( mrv == NULL ) ERROR_exit("Can't load time series data!?") ;
156 DSET_unload(inset) ;
157
158 for( kk=0 ; kk < mrv->nvec ; kk++ ){
159 spair = THD_countspikes( mrv->nvals , VECTIM_PTR(mrv,kk) ) ;
160 spk += spair.a ; stp += spair.b ;
161 }
162
163 printf("#spikiness = %.0f #steppiness = %.0f\n",spk,stp) ;
164 exit(0) ;
165 }
166