1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------------------*/
4
mytanh(float x)5 static INLINE float mytanh( float x )
6 {
7 register float ex , exi ;
8 if( x > 7.0f ) return 1.0f ;
9 else if( x < -7.0f ) return -1.0f ;
10 ex = exp(x) ; exi = 1.0f/ex ;
11 return (ex-exi)/(ex+exi) ;
12 }
13
14 /*----------------------------------------------------------------------------*/
15 /* Despike the hard way, in place (as in 3dDespike). */
16 /*----------------------------------------------------------------------------*/
17
THD_vectim_despike_L1(MRI_vectim * mrv,int localedit)18 void THD_vectim_despike_L1( MRI_vectim *mrv , int localedit )
19 {
20 int iv , nxyz , ii,jj,kk ;
21 float cut1=2.5,cut2=4.0 , sq2p,sfac , fq , tm,fac ;
22 int corder, nref , nuse ;
23 float **ref ;
24 float c21,ic21 ;
25
26 /*-- initialization --*/
27
28 sq2p = sqrt(0.5*PI) ;
29 sfac = sq2p / 1.4826f ;
30 c21 = cut2 - cut1 ;
31 ic21 = 1.0f / c21 ;
32
33 nuse = mrv->nvals ;
34 nxyz = mrv->nvec ;
35 corder = (int)rintf(nuse/30.0f) ;
36 if( corder < 2 ) corder = 2 ; else if( corder > 50 ) corder = 50 ;
37
38 /* ref functions */
39
40 nref = 2*corder+3 ; /* always odd */
41 ref = (float **)malloc( sizeof(float *) * nref ) ;
42 for( jj=0 ; jj < nref ; jj++ )
43 ref[jj] = (float *)malloc( sizeof(float) * nuse ) ;
44
45 tm = 0.5f * (nuse-1.0f) ; fac = 2.0f / nuse ;
46 for( iv=0 ; iv < nuse ; iv++ ){ /* polynomials */
47 ref[0][iv] = 1.0f ;
48 ref[1][iv] = (iv-tm)*fac ;
49 ref[2][iv] = ref[1][iv] * ref[1][iv] - 0.3333333f ;
50 }
51
52 for( jj=2,kk=1 ; kk <= corder ; kk++,jj+=2 ){ /* sines and cosines */
53 fq = (2.0*PI*kk)/nuse ;
54 for( iv=0 ; iv < nuse ; iv++ ){
55 ref[jj ][iv] = sinf(fq*iv) ;
56 ref[jj+1][iv] = cosf(fq*iv) ;
57 }
58 }
59
60 /*--- loop over voxels and do some work ---*/
61
62 AFNI_OMP_START ;
63 #pragma omp parallel if( nxyz > 99 )
64 { int ii , iv , jj , didit ;
65 float *far , *dar , *var , *fitar , *ssp , *fit , *zar ;
66 float fsig , fq , cls , snew , val ;
67
68 #pragma omp critical (DESPIKE_malloc)
69 { far = (float *) malloc( sizeof(float) * nuse ) ;
70 dar = (float *) malloc( sizeof(float) * nuse ) ;
71 var = (float *) malloc( sizeof(float) * nuse ) ;
72 fitar = (float *) malloc( sizeof(float) * nuse ) ;
73 ssp = (float *) malloc( sizeof(float) * nuse ) ;
74 fit = (float *) malloc( sizeof(float) * nref ) ;
75 }
76
77 #pragma omp for
78 for( ii=0 ; ii < nxyz ; ii++ ){ /* ii = voxel index */
79
80 /*** extract ii-th time series into far[] and dar[] ***/
81
82 zar = VECTIM_PTR(mrv,ii) ;
83 for( iv=0 ; iv < nuse ; iv++ ) far[iv] = dar[iv] = zar[iv] ;
84
85 /*** solve for L1 fit ***/
86
87 cls = cl1_solve( nuse , nref , far , ref , fit,0 ) ; /* the slow part */
88
89 if( cls < 0.0f ) continue ; /* fit failed! */
90
91 for( iv=0 ; iv < nuse ; iv++ ){ /* detrend */
92 val = fit[0] ;
93 for( jj=1 ; jj < nref ; jj+=2 ) /* unrolled (nref is odd) */
94 val += fit[jj] * ref[jj][iv] + fit[jj+1] * ref[jj+1][iv] ;
95
96 fitar[iv] = val ; /* save curve fit value */
97 var[iv] = dar[iv]-val ; /* remove fitted value = resid */
98 far[iv] = fabsf(var[iv]) ; /* abs value of resid */
99 }
100
101 /*** estimate standard deviation of detrended data ***/
102
103 fsig = sq2p * qmed_float(nuse,far) ; /* also mangles far array */
104
105 /*** process time series for spikes, editing data in dar[] ***/
106
107 if( fsig > 0.0f ){ /* data wasn't fit perfectly */
108
109 /* find spikiness for each point in time */
110
111 fq = 1.0f / fsig ;
112 for( iv=0 ; iv < nuse ; iv++ ){
113 ssp[iv] = fq * var[iv] ; /* spikiness s = how many sigma out */
114 }
115 didit = 0 ;
116
117 /* process values of |s| > cut1, editing dar[] */
118
119 for( iv=0 ; iv < nuse ; iv++ ){ /* loop over time points */
120 if( !localedit ){ /** classic 'smash' edit **/
121 if( ssp[iv] > cut1 ){
122 snew = cut1 + c21*mytanh((ssp[iv]-cut1)*ic21) ; /* edit s down */
123 dar[iv] = fitar[iv] + snew*fsig ; didit++ ;
124 } else if( ssp[iv] < -cut1 ){
125 snew = -cut1 + c21*mytanh((ssp[iv]+cut1)*ic21) ; /* edit s up */
126 dar[iv] = fitar[iv] + snew*fsig ; didit++ ;
127 }
128 } else { /** local edit **/
129 if( ssp[iv] >= cut2 || ssp[iv] <= -cut2 ){
130 int iu , id ;
131 for( iu=iv+1 ; iu < nuse ; iu++ ) /* find non-spike above */
132 if( ssp[iu] < cut2 && ssp[iu] > -cut2 ) break ;
133 for( id=iv-1 ; id >= 0 ; id-- ) /* find non-spike below */
134 if( ssp[id] < cut2 && ssp[id] > -cut2 ) break ;
135 switch( (id>=0) + 2*(iu<nuse) ){ /* compute replacement val */
136 case 3: val = 0.5*(dar[iu]+dar[id]); break; /* iu and id OK */
137 case 2: val = dar[iu] ; break; /* only iu OK */
138 case 1: val = dar[id] ; break; /* only id OK */
139 default: val = fitar[iv] ; break; /* shouldn't be */
140 }
141 dar[iv] = val ; didit++ ;
142 }
143 }
144 } /* end of loop over time points */
145
146 /* put edited result back into source */
147
148 if( didit){
149 for( iv=0 ; iv < nuse ; iv++ ) zar[iv] = dar[iv] ;
150 }
151
152 } /* end of processing time series when fsig is positive */
153
154 } /* end of loop over voxels #ii */
155
156 #pragma omp critical (DESPIKE_malloc)
157 { free(fit); free(ssp); free(fitar); free(var); free(dar); free(far); }
158
159 } /* end OpenMP */
160 AFNI_OMP_END ;
161
162 for( jj=0 ; jj < nref ; jj++ ) free(ref[jj]) ;
163 free(ref) ;
164 return ;
165 }
166