1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 /************************************************************************/
8 /******* Hack to remove large spikes from time series data (oog). *******/
9 /******* 30 Aug 2002 - RWCox *******/
10 /************************************************************************/
11
12 #define TFAC 0.1 /* scale factor for -ssave dataset */
13 #define ITFAC 10.0 /* inverse of above */
14
15 /*----------------------------------------------------------------------*/
16
mytanh(float x)17 static INLINE float mytanh( float x )
18 {
19 register float ex , exi ;
20 if( x > 7.0f ) return 1.0f ; /* 03 Sep: check for stupid inputs */
21 else if( x < -7.0f ) return -1.0f ;
22 ex = exp(x) ; exi = 1.0f/ex ;
23 return (ex-exi)/(ex+exi) ;
24 }
25
26 /*============================================================================*/
27 /***************** The following stuff was added 29 Nov 2013 ******************/
28 /*----------------------------------------------------------------------------*/
29
30 #undef SWAP
31 #define SWAP(x,y) (temp=x,x=y,y=temp)
32
33 #undef SORT2
34 #define SORT2(a,b) if(a>b) SWAP(a,b)
35
36 /*--- fast median of 9 values ---*/
37
median9f(float * p)38 static INLINE float median9f(float *p)
39 {
40 register float temp ;
41 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
42 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
43 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
44 SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
45 SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
46 SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
47 SORT2(p[4],p[2]) ; return(p[4]) ;
48 }
49
50 /*--- get the local median and MAD of values vec[j-4 .. j+4] ---*/
51
52 #undef mead9
53 #define mead9(j) \
54 { float qqq[9] ; int jj = (j)-4 ; \
55 if( jj < 0 ) jj = 0; else if( jj+9 > num ) jj = num-9; \
56 qqq[0] = vec[jj+0]; qqq[1] = vec[jj+1]; qqq[2] = vec[jj+2]; \
57 qqq[3] = vec[jj+3]; qqq[4] = vec[jj+4]; qqq[5] = vec[jj+5]; \
58 qqq[6] = vec[jj+6]; qqq[7] = vec[jj+7]; qqq[8] = vec[jj+8]; \
59 med = median9f(qqq); qqq[0] = fabsf(qqq[0]-med); \
60 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
61 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
62 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
63 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
64 mad = median9f(qqq); }
65
66 /*-------------------------------------------------------------------------*/
67 /*! Remove spikes from a time series, in a very simplistic way.
68 Return value is the number of spikes that were squashed [RWCox].
69 *//*-----------------------------------------------------------------------*/
70
DES_despike9(int num,float * vec,float * wks)71 int DES_despike9( int num , float *vec , float *wks )
72 {
73 int ii , nsp ; float *zma,*zme , med,mad,val ;
74
75 if( num < 9 || vec == NULL ) return 0 ;
76
77 zme = wks ; zma = zme + num ;
78
79 for( ii=0 ; ii < num ; ii++ ){
80 mead9(ii) ; zme[ii] = med ; zma[ii] = mad ;
81 }
82 mad = qmed_float(num,zma) ;
83 if( mad <= 0.0f ){ if( wks == NULL ) free(zme); return 0; } /* should not happen */
84 mad *= 6.789f ; /* threshold value */
85
86 for( nsp=ii=0 ; ii < num ; ii++ )
87 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
88
89 return nsp ;
90 }
91 #undef mead9
92
93 /*----------------------------------------------------------------------------*/
94 /* Similar code to the above, but for spans of length 25 instead of 9 */
95
median25f(float * p)96 static INLINE float median25f(float *p)
97 {
98 register float temp ;
99 SORT2(p[0], p[1]) ; SORT2(p[3], p[4]) ; SORT2(p[2], p[4]) ;
100 SORT2(p[2], p[3]) ; SORT2(p[6], p[7]) ; SORT2(p[5], p[7]) ;
101 SORT2(p[5], p[6]) ; SORT2(p[9], p[10]) ; SORT2(p[8], p[10]) ;
102 SORT2(p[8], p[9]) ; SORT2(p[12], p[13]) ; SORT2(p[11], p[13]) ;
103 SORT2(p[11], p[12]) ; SORT2(p[15], p[16]) ; SORT2(p[14], p[16]) ;
104 SORT2(p[14], p[15]) ; SORT2(p[18], p[19]) ; SORT2(p[17], p[19]) ;
105 SORT2(p[17], p[18]) ; SORT2(p[21], p[22]) ; SORT2(p[20], p[22]) ;
106 SORT2(p[20], p[21]) ; SORT2(p[23], p[24]) ; SORT2(p[2], p[5]) ;
107 SORT2(p[3], p[6]) ; SORT2(p[0], p[6]) ; SORT2(p[0], p[3]) ;
108 SORT2(p[4], p[7]) ; SORT2(p[1], p[7]) ; SORT2(p[1], p[4]) ;
109 SORT2(p[11], p[14]) ; SORT2(p[8], p[14]) ; SORT2(p[8], p[11]) ;
110 SORT2(p[12], p[15]) ; SORT2(p[9], p[15]) ; SORT2(p[9], p[12]) ;
111 SORT2(p[13], p[16]) ; SORT2(p[10], p[16]) ; SORT2(p[10], p[13]) ;
112 SORT2(p[20], p[23]) ; SORT2(p[17], p[23]) ; SORT2(p[17], p[20]) ;
113 SORT2(p[21], p[24]) ; SORT2(p[18], p[24]) ; SORT2(p[18], p[21]) ;
114 SORT2(p[19], p[22]) ; SORT2(p[8], p[17]) ; SORT2(p[9], p[18]) ;
115 SORT2(p[0], p[18]) ; SORT2(p[0], p[9]) ; SORT2(p[10], p[19]) ;
116 SORT2(p[1], p[19]) ; SORT2(p[1], p[10]) ; SORT2(p[11], p[20]) ;
117 SORT2(p[2], p[20]) ; SORT2(p[2], p[11]) ; SORT2(p[12], p[21]) ;
118 SORT2(p[3], p[21]) ; SORT2(p[3], p[12]) ; SORT2(p[13], p[22]) ;
119 SORT2(p[4], p[22]) ; SORT2(p[4], p[13]) ; SORT2(p[14], p[23]) ;
120 SORT2(p[5], p[23]) ; SORT2(p[5], p[14]) ; SORT2(p[15], p[24]) ;
121 SORT2(p[6], p[24]) ; SORT2(p[6], p[15]) ; SORT2(p[7], p[16]) ;
122 SORT2(p[7], p[19]) ; SORT2(p[13], p[21]) ; SORT2(p[15], p[23]) ;
123 SORT2(p[7], p[13]) ; SORT2(p[7], p[15]) ; SORT2(p[1], p[9]) ;
124 SORT2(p[3], p[11]) ; SORT2(p[5], p[17]) ; SORT2(p[11], p[17]) ;
125 SORT2(p[9], p[17]) ; SORT2(p[4], p[10]) ; SORT2(p[6], p[12]) ;
126 SORT2(p[7], p[14]) ; SORT2(p[4], p[6]) ; SORT2(p[4], p[7]) ;
127 SORT2(p[12], p[14]) ; SORT2(p[10], p[14]) ; SORT2(p[6], p[7]) ;
128 SORT2(p[10], p[12]) ; SORT2(p[6], p[10]) ; SORT2(p[6], p[17]) ;
129 SORT2(p[12], p[17]) ; SORT2(p[7], p[17]) ; SORT2(p[7], p[10]) ;
130 SORT2(p[12], p[18]) ; SORT2(p[7], p[12]) ; SORT2(p[10], p[18]) ;
131 SORT2(p[12], p[20]) ; SORT2(p[10], p[20]) ; SORT2(p[10], p[12]) ;
132 return (p[12]);
133 }
134
135 /*--- get the local median and MAD of values vec[j-12 .. j+12] ---*/
136
137 #undef mead25
138 #define mead25(j) \
139 { float qqq[25] ; int jj=(j)-12 ; register int pp; \
140 if( jj < 0 ) jj = 0; else if( jj+25 > num ) jj = num-25; \
141 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = vec[jj+pp] ; \
142 med = median25f(qqq) ; \
143 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = fabsf(qqq[pp]-med) ; \
144 mad = median25f(qqq); }
145
DES_despike25(int num,float * vec,float * wks)146 int DES_despike25( int num , float *vec , float *wks )
147 {
148 int ii , nsp ; float *zma,*zme , med,mad,val ;
149
150 if( vec == NULL ) return 0 ;
151 if( num < 25 ) return DES_despike9(num,vec,wks) ;
152
153 zme = wks ; zma = zme + num ;
154
155 for( ii=0 ; ii < num ; ii++ ){
156 mead25(ii) ; zme[ii] = med ; zma[ii] = mad ;
157 }
158 mad = qmed_float(num,zma) ;
159 if( mad <= 0.0f ){ if( wks == NULL ) free(zme); return 0; } /* should not happen */
160 mad *= 6.789f ; /* threshold value */
161
162 for( nsp=ii=0 ; ii < num ; ii++ )
163 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
164
165 return nsp ;
166 }
167 #undef mead25
168 #undef SORT2
169 #undef SWAP
170
171 /*----------------------------------------------------------------------------*/
172
DES_workspace_size(int ntim,int nref)173 int DES_workspace_size( int ntim, int nref )
174 {
175 return 2*ntim ;
176 }
177
178 /*----------------------------------------------------------------------------*/
179
180 static int use_des25 = 0 ;
181
DES_solve(MRI_IMAGE * psinv,float * z,float * coef,float * wks)182 float DES_solve( MRI_IMAGE *psinv , float *z , float *coef , float *wks )
183 {
184 float *psar , *iar , zi ; int ii,jj , ntim,nref ;
185
186 ntim = psinv->ny ;
187 nref = psinv->nx ;
188 psar = MRI_FLOAT_PTR(psinv) ;
189
190 /* step 1: despike the data the simplistic way */
191
192 if( use_des25 )
193 (void) DES_despike25( ntim , z , wks ) ;
194 else
195 (void) DES_despike9( ntim , z , wks ) ;
196
197 /* least squares solve the equations with the modified data */
198
199 for( jj=0 ; jj < nref ; jj++ ) coef[jj] = 0.0f ;
200
201 for( ii=0 ; ii < ntim ; ii++ ){
202 iar = psar + ii*nref ; zi = z[ii] ;
203 for( jj=0 ; jj < nref ; jj++ ) coef[jj] += iar[jj]*zi ;
204 }
205
206 return 0.0f ;
207 }
208
209 /*----------------------------------------------------------------------------*/
210
DES_get_psinv(int ntim,int nref,float ** ref)211 MRI_IMAGE * DES_get_psinv( int ntim , int nref , float **ref )
212 {
213 MRI_IMAGE *refim , *psinv ; float *refar , *jar ; int ii , jj ;
214
215 refim = mri_new(ntim,nref,MRI_float) ;
216 refar = MRI_FLOAT_PTR(refim) ;
217 for( jj=0 ; jj < nref ; jj++ ){
218 jar = refar + jj*ntim ;
219 for( ii=0 ; ii < ntim ; ii++ ) jar[ii] = ref[jj][ii] ;
220 }
221 mri_matrix_psinv_svd(1) ;
222 psinv = mri_matrix_psinv(refim,NULL,0.0f) ;
223 mri_free(refim) ;
224 return psinv ;
225 }
226
227 /*----------------------------------------------------------------------*/
228
main(int argc,char * argv[])229 int main( int argc , char *argv[] )
230 {
231 THD_3dim_dataset *dset , *oset=NULL , *tset=NULL ;
232 int nvals , iv , nxyz , ii,jj,kk , iarg , kz,kzold ;
233 float cut1=2.5,cut2=4.0 , sq2p,sfac , fq ;
234 MRI_IMAGE *flim ;
235 char *prefix="despike" , *tprefix=NULL ;
236
237 int corder=-1 , nref , ignore=0 , polort=2 , nuse , nomask=0 ;
238 int nspike, nbig, nproc ;
239 float **ref ;
240 float c21,ic21 , pspike,pbig ;
241 short *sar , *qar ;
242 byte *tar , *mask=NULL ;
243 float *zar , *yar ;
244 int in_datum , out_datum ;
245 int localedit=0 ; /* 04 Apr 2007 */
246 int verb=1 ;
247
248 int do_NEW = 0 ; /* 29 Nov 2013 */
249 MRI_IMAGE *NEW_psinv=NULL ;
250 int dilate = 4 ; /* 04 Dec 2013 */
251 int ctim = 0 ;
252
253 /*----- Read command line -----*/
254
255 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
256
257 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
258 printf("\n"
259 "Usage: 3dDespike [options] dataset\n"
260 "\n"
261 "Removes 'spikes' from the 3D+time input dataset and writes\n"
262 "a new dataset with the spike values replaced by something\n"
263 "more pleasing to the eye.\n"
264 "\n"
265 "------------------\n"
266 "Outline of Method:\n"
267 "------------------\n"
268 " * L1 fit a smooth-ish curve to each voxel time series\n"
269 " [see -corder option for description of the curve]\n"
270 " [see -NEW option for a different & faster fitting method]\n"
271 " * Compute the MAD of the difference between the curve and\n"
272 " the data time series (the residuals).\n"
273 " * Estimate the standard deviation 'sigma' of the residuals\n"
274 " from the MAD.\n"
275 " * For each voxel value, define s = (value-curve)/sigma.\n"
276 " * Values with s > c1 are replaced with a value that yields\n"
277 " a modified s' = c1+(c2-c1)*tanh((s-c1)/(c2-c1)).\n"
278 " * c1 is the threshold value of s for a 'spike' [default c1=2.5].\n"
279 " * c2 is the upper range of the allowed deviation from the curve:\n"
280 " s=[c1..infinity) is mapped to s'=[c1..c2) [default c2=4].\n"
281 "\n"
282 "An alternative method for replacing the spike value is provided\n"
283 "by the '-localedit' option, and that method is preferred by\n"
284 "many users.\n"
285 "\n"
286 "The input dataset can be stored in short or float formats.\n"
287 "The output dataset will always be stored in floats. [Feb 2017]\n"
288 "\n"
289 "--------\n"
290 "Options:\n"
291 "--------\n"
292 " -ignore I = Ignore the first I points in the time series:\n"
293 " these values will just be copied to the\n"
294 " output dataset [default I=0].\n"
295 " -corder L = Set the curve fit order to L:\n"
296 " the curve that is fit to voxel data v(t) is\n"
297 "\n"
298 " k=L [ (2*PI*k*t) (2*PI*k*t) ]\n"
299 " f(t) = a+b*t+c*t*t + SUM [ d * sin(--------) + e * cos(--------) ]\n"
300 " k=1 [ k ( T ) k ( T ) ]\n"
301 "\n"
302 " where T = duration of time series;\n"
303 " the a,b,c,d,e parameters are chosen to minimize\n"
304 " the sum over t of |v(t)-f(t)| (L1 regression);\n"
305 " this type of fitting is is insensitive to large\n"
306 " spikes in the data. The default value of L is\n"
307 " NT/30, where NT = number of time points.\n"
308 "\n"
309 " -cut c1 c2 = Alter default values for the spike cut values\n"
310 " [default c1=2.5, c2=4.0].\n"
311 "\n"
312 " -prefix pp = Save de-spiked dataset with prefix 'pp'\n"
313 " [default pp='despike']\n"
314 "\n"
315 " -ssave ttt = Save 'spikiness' measure s for each voxel into a\n"
316 " 3D+time dataset with prefix 'ttt' [default=no save]\n"
317 "\n"
318 " -nomask = Process all voxels\n"
319 " [default=use a mask of high-intensity voxels, ]\n"
320 " [as created via '3dAutomask -dilate 4 dataset'].\n"
321 "\n"
322 " -dilate nd = Dilate 'nd' times (as in 3dAutomask). The default\n"
323 " value of 'nd' is 4.\n"
324 "\n"
325 " -q[uiet] = Don't print '++' informational messages.\n"
326 "\n"
327 " -localedit = Change the editing process to the following:\n"
328 " If a voxel |s| value is >= c2, then replace\n"
329 " the voxel value with the average of the two\n"
330 " nearest non-spike (|s| < c2) values; the first\n"
331 " one previous and the first one after.\n"
332 " Note that the c1 cut value is not used here.\n"
333 "\n"
334 " -NEW = Use the 'new' method for computing the fit, which\n"
335 " should be faster than the L1 method for long time\n"
336 " series (200+ time points); however, the results\n"
337 " are similar but NOT identical. [29 Nov 2013]\n"
338 " * You can also make the program use the 'new'\n"
339 " method by setting the environment variable\n"
340 " AFNI_3dDespike_NEW\n"
341 " to the value YES; as in\n"
342 " setenv AFNI_3dDespike_NEW YES (csh)\n"
343 " export AFNI_3dDespike_NEW=YES (bash)\n"
344 " * If this variable is set to YES, you can turn off\n"
345 " the '-NEW' processing by using the '-OLD' option.\n"
346 " -->>* For time series more than 500 points long, the\n"
347 " '-OLD' algorithm is tremendously slow. You should\n"
348 " use the '-NEW' algorith in such cases.\n"
349 " ** At some indeterminate point in the future, the '-NEW'\n"
350 " method will become the default!\n"
351 " -->>* As of 29 Sep 2016, '-NEW' is the default if there\n"
352 " is more than 500 points in the time series dataset.\n"
353 "\n"
354 " -NEW25 = A slightly more aggressive despiking approach than\n"
355 " the '-NEW' method.\n"
356 "\n"
357 "--------\n"
358 "Caveats:\n"
359 "--------\n"
360 "* Despiking may interfere with image registration, since head\n"
361 " movement may produce 'spikes' at the edge of the brain, and\n"
362 " this information would be used in the registration process.\n"
363 " This possibility has not been explored or calibrated.\n"
364 "\n"
365 "* [LATER] Actually, it seems like the registration problem\n"
366 " does NOT happen, and in fact, despiking seems to help!\n"
367 "\n"
368 "* Check your data visually before and after despiking and\n"
369 " registration!\n"
370 ) ;
371
372 PRINT_AFNI_OMP_USAGE("3dDespike",NULL) ;
373 PRINT_COMPILE_DATE ; exit(0) ;
374 }
375
376 /** AFNI package setup and logging **/
377
378 mainENTRY("3dDespike main"); machdep(); AFNI_logger("3dDespike",argc,argv);
379 PRINT_VERSION("3dDespike") ; AUTHOR("RW Cox") ;
380
381 /** parse options **/
382
383 if( AFNI_yesenv("AFNI_3dDespike_NEW") ) do_NEW = 1 ; /* 29 Nov 2013 */
384
385 iarg = 1 ;
386 while( iarg < argc && argv[iarg][0] == '-' ){
387
388 if( strncmp(argv[iarg],"-q",2) == 0 ){ /* 04 Apr 2007 */
389 verb = 0 ; iarg++ ; continue ;
390 }
391 if( strncmp(argv[iarg],"-v",2) == 0 ){
392 verb++ ; iarg++ ; continue ;
393 }
394
395 if( strcmp(argv[iarg],"-NEW") == 0 ){ /* 29 Nov 2013 */
396 do_NEW = 1 ; iarg++ ; continue ;
397 }
398 if( strcmp(argv[iarg],"-NEW25") == 0 ){ /* 29 Sep 2016 */
399 do_NEW = 1 ; use_des25 = 1 ; cut1 = 2.5f ; cut2 = 3.2f ; iarg++ ; continue ;
400 }
401 if( strcmp(argv[iarg],"-OLD") == 0 ){
402 do_NEW = 0 ; iarg++ ; continue ;
403 }
404
405 /** -localedit **/
406
407 if( strcmp(argv[iarg],"-localedit") == 0 ){ /* 04 Apr 2007 */
408 localedit = 1 ; iarg++ ; continue ;
409 }
410
411 /** don't use masking **/
412
413 if( strcmp(argv[iarg],"-nomask") == 0 ){
414 nomask = 1 ; iarg++ ; continue ;
415 }
416
417 /** dilation count [04 Dec 2013] **/
418
419 if( strcmp(argv[iarg],"-dilate") == 0 ){
420 dilate = (int)strtod(argv[++iarg],NULL) ;
421 if( dilate <= 0 ) dilate = 1 ;
422 else if( dilate > 99 ) dilate = 99 ;
423 iarg++ ; continue ;
424 }
425
426 /** output dataset prefix **/
427
428 if( strcmp(argv[iarg],"-prefix") == 0 ){
429 prefix = argv[++iarg] ;
430 if( !THD_filename_ok(prefix) ) ERROR_exit("-prefix is not good");
431 iarg++ ; continue ;
432 }
433
434 /** ratio dataset prefix **/
435
436 if( strcmp(argv[iarg],"-ssave") == 0 ){
437 tprefix = argv[++iarg] ;
438 if( !THD_filename_ok(tprefix) ) ERROR_exit("-ssave prefix is not good");
439 iarg++ ; continue ;
440 }
441
442 /** trigonometric polynomial order **/
443
444 if( strcmp(argv[iarg],"-corder") == 0 ){
445 corder = strtol( argv[++iarg] , NULL , 10 ) ;
446 if( corder < 0 ) ERROR_exit("Illegal value of -corder");
447 iarg++ ; continue ;
448 }
449
450 /** how much to ignore at start **/
451
452 if( strcmp(argv[iarg],"-ignore") == 0 ){
453 ignore = strtol( argv[++iarg] , NULL , 10 ) ;
454 if( ignore < 0 ) ERROR_exit("Illegal value of -ignore");
455 iarg++ ; continue ;
456 }
457
458 /** thresholds for s ratio **/
459
460 if( strcmp(argv[iarg],"-cut") == 0 ){
461 cut1 = strtod( argv[++iarg] , NULL ) ;
462 cut2 = strtod( argv[++iarg] , NULL ) ;
463 if( cut1 < 1.0 || cut2 < cut1+0.5 )
464 ERROR_exit("Illegal values after -cut");
465 iarg++ ; continue ;
466 }
467
468 ERROR_exit("Unknown option: %s",argv[iarg]) ;
469 }
470
471 c21 = cut2-cut1 ; ic21 = 1.0/c21 ;
472
473 /*----- read input dataset -----*/
474
475 if( iarg >= argc ) ERROR_exit("No input dataset!!??");
476
477 dset = THD_open_dataset( argv[iarg] ) ;
478 CHECK_OPEN_ERROR(dset,argv[iarg]) ;
479 in_datum = DSET_BRICK_TYPE(dset,0) ;
480 if( (in_datum != MRI_short && in_datum != MRI_float) || !DSET_datum_constant(dset) )
481 ERROR_exit("Can't process non-short, non-float dataset!") ;
482
483 out_datum = MRI_float ;
484 if( verb && (in_datum == MRI_short) ){
485 INFO_message("Input dataset is in short format, but output will be in float format") ;
486 }
487
488 nvals = DSET_NVALS(dset) ; nuse = nvals - ignore ;
489 if( nuse < 15 )
490 ERROR_exit("Can't use dataset with < 15 time points per voxel!") ;
491
492 if( nuse > 500 && !do_NEW ){
493 if( verb )
494 INFO_message("Switching to '-NEW' method since number of time points = %d > 500",nuse) ;
495 do_NEW = 1 ;
496 }
497 if( use_des25 && nuse <= 99 ){
498 if( verb ){
499 INFO_message("'-NEW25' method was ordered, but need more than 99 time points for that") ;
500 INFO_message(" switching to the '-NEW' method instead") ;
501 }
502 use_des25 = 0 ;
503 }
504
505 if( verb ) INFO_message("ignoring first %d time points, using last %d",ignore,nuse);
506 if( corder > 0 && 4*corder+2 > nuse ){
507 ERROR_exit("-corder %d is too big for NT=%d",corder,nvals) ;
508 } else if( corder < 0 ){
509 corder = rint(nuse/30.0) ; if( corder > 50 && !do_NEW ) corder = 50 ;
510 if( verb ) INFO_message("using %d time points => -corder %d",nuse,corder) ;
511 } else {
512 if( verb ) INFO_message("-corder %d set from command line",corder) ;
513 }
514 nxyz = DSET_NVOX(dset) ;
515 if( verb ) INFO_message("Loading dataset %s",argv[iarg]) ;
516 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
517
518 /*-- create automask --*/
519
520 if( !nomask ){
521 mask = THD_automask( dset ) ;
522 ii = THD_countmask( DSET_NVOX(dset) , mask ) ;
523 if( verb && ii > 0 )
524 INFO_message("%d voxels in the automask [out of %d in dataset]",ii,DSET_NVOX(dset)) ;
525 else if( ii == 0 )
526 ERROR_exit("Nothing to process -- automask is empty :(") ;
527
528 for( ii=0 ; ii < dilate ; ii++ )
529 THD_mask_dilate( DSET_NX(dset), DSET_NY(dset), DSET_NZ(dset), mask, 3, 2 ) ;
530
531 ii = THD_countmask( DSET_NVOX(dset) , mask ) ;
532 if( verb )
533 INFO_message("%d voxels in the dilated automask [out of %d in dataset]",ii,DSET_NVOX(dset)) ;
534 if( ii == 0 )
535 ERROR_exit("Nothing to process -- no voxels in automask?!") ;
536 } else {
537 if( verb ) INFO_message("processing all %d voxels in dataset",DSET_NVOX(dset)) ;
538 }
539
540 /*-- create empty despiked dataset --*/
541
542 oset = EDIT_empty_copy( dset ) ;
543 EDIT_dset_items( oset ,
544 ADN_prefix , prefix ,
545 ADN_brick_fac , NULL ,
546 ADN_datum_all , out_datum ,
547 ADN_none ) ;
548
549 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(oset)) )
550 ERROR_exit("output dataset already exists: %s",DSET_HEADNAME(oset));
551
552 tross_Copy_History( oset , dset ) ;
553 tross_Make_History( "3dDespike" , argc , argv , oset ) ;
554
555 /* copy the ignored bricks, if any */
556
557 for( iv=0 ; iv < ignore ; iv++ ){
558 flim = THD_extract_float_brick(iv,dset) ;
559 EDIT_substitute_brick( oset , iv , MRI_float , MRI_FLOAT_PTR(flim) ) ;
560 mri_clear_and_free(flim) ;
561 }
562
563 /* create rest of the bricks (will be filled with zeros) */
564
565 for( iv=ignore ; iv < nvals ; iv++ )
566 EDIT_substitute_brick( oset , iv , out_datum , NULL ) ;
567
568 /*-- setup to save a threshold statistic dataset, if desired --*/
569
570 if( tprefix != NULL ){
571 float *fac ;
572 tset = EDIT_empty_copy( dset ) ;
573 fac = (float *) malloc( sizeof(float) * nvals ) ;
574 for( ii=0 ; ii < nvals ; ii++ ) fac[ii] = TFAC ;
575 EDIT_dset_items( tset ,
576 ADN_prefix , tprefix ,
577 ADN_brick_fac , fac ,
578 ADN_datum_all , MRI_byte ,
579 ADN_func_type , FUNC_FIM_TYPE ,
580 ADN_none ) ;
581 free(fac) ;
582
583 tross_Copy_History( tset , dset ) ;
584 tross_Make_History( "3dDespike" , argc , argv , tset ) ;
585
586 #if 0
587 if( THD_is_file(DSET_HEADNAME(tset)) )
588 ERROR_exit("-ssave dataset already exists");
589 #endif
590
591 tross_Copy_History( tset , dset ) ;
592 tross_Make_History( "3dDespike" , argc , argv , tset ) ;
593
594 for( iv=0 ; iv < nvals ; iv++ )
595 EDIT_substitute_brick( tset , iv , MRI_byte , NULL ) ;
596 }
597
598 /*-- setup to find spikes --*/
599
600 sq2p = sqrt(0.5*PI) ;
601 sfac = sq2p / 1.4826f ;
602
603 /* make ref functions */
604
605 nref = 2*corder+3 ;
606 ref = (float **) malloc( sizeof(float *) * nref ) ;
607 for( jj=0 ; jj < nref ; jj++ )
608 ref[jj] = (float *) malloc( sizeof(float) * nuse ) ;
609
610 /* r(t) = 1 */
611
612 for( iv=0 ; iv < nuse ; iv++ ) ref[0][iv] = 1.0 ;
613 jj = 1 ;
614
615 /* r(t) = t - tmid */
616
617 { float tm = 0.5 * (nuse-1.0) ; float fac = 2.0 / nuse ;
618 for( iv=0 ; iv < nuse ; iv++ ) ref[1][iv] = (iv-tm)*fac ;
619 jj = 2 ;
620
621 /* r(t) = (t-tmid)**jj [NB: polort==2] */
622
623 for( ; jj <= polort ; jj++ )
624 for( iv=0 ; iv < nuse ; iv++ )
625 ref[jj][iv] = pow( (iv-tm)*fac , (double)jj ) ;
626 }
627
628 for( kk=1 ; kk <= corder ; kk++ ){ /* sines and cosines */
629 fq = (2.0*PI*kk)/nuse ;
630
631 /* r(t) = sin(2*PI*k*t/N) */
632
633 for( iv=0 ; iv < nuse ; iv++ )
634 ref[jj][iv] = sin(fq*iv) ;
635 jj++ ;
636
637 /* r(t) = cos(2*PI*k*t/N) */
638
639 for( iv=0 ; iv < nuse ; iv++ )
640 ref[jj][iv] = cos(fq*iv) ;
641 jj++ ;
642 }
643
644 /****** setup for the NEW solution method [29 Nov 2013] ******/
645
646 if( do_NEW ){
647 NEW_psinv = DES_get_psinv(nuse,nref,ref) ;
648 if( verb )
649 INFO_message("Procesing time series with %s model fit algorithm",
650 (use_des25) ? "NEW25" : "NEW" ) ;
651 } else {
652 if( verb )
653 INFO_message("Procesing time series with OLD model fit algorithm") ;
654 }
655
656 /*--- loop over voxels and do work ---*/
657
658 #define Laplace_t2p(val) ( 1.0 - nifti_stat2cdf( (val), 15, 0.0, 1.4427 , 0.0 ) )
659
660 if( verb ){
661 if( !localedit ){
662 INFO_message("smash edit thresholds: %.1f .. %.1f MADs",cut1*sq2p,cut2*sq2p) ;
663 ININFO_message(" [ %.3f%% .. %.3f%% of normal distribution]",
664 200.0*qg(cut1*sfac) , 200.0*qg(cut2*sfac) ) ;
665 ININFO_message(" [ %.3f%% .. %.3f%% of Laplace distribution]" ,
666 100.0*Laplace_t2p(cut1) , 100.0*Laplace_t2p(cut2) ) ;
667 } else {
668 INFO_message("local edit threshold: %.1f MADS",cut2*sq2p) ;
669 ININFO_message(" [ %.3f%% of normal distribution]",
670 200.0*qg(cut2*sfac) ) ;
671 ININFO_message(" [ %.3f%% of Laplace distribution]",
672 100.0*Laplace_t2p(cut1) ) ;
673 }
674 }
675
676 kzold = -1 ;
677 nspike = 0 ; nbig = 0 ; nproc = 0 ; ctim = NI_clock_time() ;
678
679 /* OpenMP-ized across voxels */
680
681 AFNI_OMP_START ;
682 #pragma omp parallel if( nxyz > 666 )
683 { int ii , iv , iu , id , jj ;
684 float *far , *dar , *var , *fitar , *ssp , *fit , *zar ;
685 short *sar , *qar ; byte *tar ;
686 float fsig , fq , cls , snew , val ;
687 float *NEW_wks=NULL ;
688
689 #pragma omp critical (DESPIKE_malloc)
690 { far = (float *) malloc( sizeof(float) * nvals ) ;
691 dar = (float *) malloc( sizeof(float) * nvals ) ;
692 var = (float *) malloc( sizeof(float) * nvals ) ;
693 fitar = (float *) malloc( sizeof(float) * nvals ) ;
694 ssp = (float *) malloc( sizeof(float) * nvals ) ;
695 fit = (float *) malloc( sizeof(float) * nref ) ;
696 if( do_NEW )
697 NEW_wks = (float *)malloc(sizeof(float)*DES_workspace_size(nuse,nref)) ;
698 }
699
700 #ifdef USE_OMP
701 if( verb )
702 INFO_message("start OpenMP thread #%d",omp_get_thread_num()) ;
703 #endif
704
705 #pragma omp for
706 for( ii=0 ; ii < nxyz ; ii++ ){ /* ii = voxel index */
707
708 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip this voxel */
709
710 #ifndef USE_OMP
711 kz = DSET_index_to_kz(dset,ii) ; /* starting a new slice */
712 if( kz != kzold ){
713 if( verb ){
714 fprintf(stderr, "++ start slice %2d",kz ) ;
715 if( nproc > 0 ){
716 pspike = (100.0*nspike)/nproc ;
717 pbig = (100.0*nbig )/nproc ;
718 fprintf(stderr,
719 "; so far %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]",
720 nproc,nspike,pspike,nbig,pbig ) ;
721 }
722 fprintf(stderr,"\n") ;
723 }
724 kzold = kz ;
725 }
726 #else
727 if( verb && ii % 2345 == 1234 ) fprintf(stderr,".") ;
728 #endif
729
730 /*** extract ii-th time series into far[] ***/
731
732 #if 1 /* 20 Feb 2017 */
733 if( ignore > 0 ){
734 (void)THD_extract_array(ii,dset,0,dar) ;
735 for( iv=0 ; iv < nuse ; iv++ )
736 far[iv] = dar[iv+ignore] ;
737 } else {
738 (void)THD_extract_array(ii,dset,0,far) ;
739 }
740 #else
741 switch( datum ){
742 case MRI_short:
743 for( iv=0 ; iv < nuse ; iv++ ){
744 qar = DSET_ARRAY(dset,iv+ignore) ; /* skip ignored data */
745 far[iv] = (float)qar[ii] ;
746 }
747 break ;
748 case MRI_float:
749 for( iv=0 ; iv < nuse ; iv++ ){
750 zar = DSET_ARRAY(dset,iv+ignore) ;
751 far[iv] = zar[ii] ;
752 }
753 break ;
754 }
755 #endif
756 AAmemcpy(dar,far,sizeof(float)*nuse) ; /* copy time series into dar[] */
757
758 /*** solve for L1 fit ***/
759
760 if( do_NEW )
761 cls = DES_solve( NEW_psinv , far , fit , NEW_wks ) ; /* 29 Nov 2013 */
762 else
763 cls = cl1_solve( nuse , nref , far , ref , fit,0 ) ; /* the slow part */
764
765 if( cls < 0.0f ){ /* fit failed! */
766 #if 0
767 fprintf(stderr,"curve fit fails at voxel %d %d %d\n",
768 DSET_index_to_ix(dset,ii) ,
769 DSET_index_to_jy(dset,ii) ,
770 DSET_index_to_kz(dset,ii) ) ;
771 #endif
772 continue ; /* skip this voxel */
773 }
774
775 for( iv=0 ; iv < nuse ; iv++ ){ /* detrend */
776 val = fit[0]
777 + fit[1]*ref[1][iv] /* quadratic part of curve fit */
778 + fit[2]*ref[2][iv] ;
779 for( jj=3 ; jj < nref ; jj++ ) /* rest of curve fit */
780 val += fit[jj] * ref[jj][iv] ;
781
782 fitar[iv] = val ; /* save curve fit value */
783 var[iv] = dar[iv]-val ; /* remove fitted value = resid */
784 far[iv] = fabsf(var[iv]) ; /* abs value of resid */
785 }
786
787 /*** compute estimate standard deviation of detrended data ***/
788
789 fsig = sq2p * qmed_float(nuse,far) ; /* also mangles far array */
790
791 /*** process time series for spikes, editing data in dar[] ***/
792
793 if( fsig > 0.0f ){ /* data wasn't fit perfectly */
794
795 /* find spikiness for each point in time */
796
797 fq = 1.0f / fsig ;
798 for( iv=0 ; iv < nuse ; iv++ ){
799 ssp[iv] = fq * var[iv] ; /* spikiness s = how many sigma out */
800 }
801
802 /* save spikiness in -ssave datset */
803
804 if( tset != NULL ){
805 for( iv=0 ; iv < nuse ; iv++ ){
806 tar = DSET_ARRAY(tset,iv+ignore) ;
807 snew = ITFAC*fabsf(ssp[iv]) ; /* scale for byte storage */
808 tar[ii] = BYTEIZE(snew) ; /* cf. mrilib.h */
809 }
810 }
811
812 /* process values of |s| > cut1, editing dar[] */
813
814 for( iv=0 ; iv < nuse ; iv++ ){ /* loop over time points */
815 if( !localedit ){ /** classic 'smash' edit **/
816 if( ssp[iv] > cut1 ){
817 snew = cut1 + c21*mytanh((ssp[iv]-cut1)*ic21) ; /* edit s down */
818 dar[iv] = fitar[iv] + snew*fsig ;
819 #pragma omp critical (DESPIKE_counter)
820 { nspike++ ; if( ssp[iv] > cut2 ) nbig++ ; }
821 } else if( ssp[iv] < -cut1 ){
822 snew = -cut1 + c21*mytanh((ssp[iv]+cut1)*ic21) ; /* edit s up */
823 dar[iv] = fitar[iv] + snew*fsig ;
824 #pragma omp critical (DESPIKE_counter)
825 { nspike++ ; if( ssp[iv] < -cut2 ) nbig++ ; }
826 }
827 } else { /** local edit: 04 Apr 2007 **/
828 if( ssp[iv] >= cut2 || ssp[iv] <= -cut2 ){
829 for( iu=iv+1 ; iu < nuse ; iu++ ) /* find non-spike above */
830 if( ssp[iu] < cut2 && ssp[iu] > -cut2 ) break ;
831 for( id=iv-1 ; id >= 0 ; id-- ) /* find non-spike below */
832 if( ssp[id] < cut2 && ssp[id] > -cut2 ) break ;
833 switch( (id>=0) + 2*(iu<nuse) ){ /* compute replacement val */
834 case 3: val = 0.5*(dar[iu]+dar[id]); break; /* iu and id OK */
835 case 2: val = dar[iu] ; break; /* only iu OK */
836 case 1: val = dar[id] ; break; /* only id OK */
837 default: val = fitar[iv] ; break; /* shouldn't be */
838 }
839 dar[iv] = val ;
840 #pragma omp critical (DESPIKE_counter)
841 { nspike++ ; nbig++ ; }
842 }
843 }
844 } /* end of loop over time points */
845 #pragma omp atomic
846 nproc += nuse ; /* number data points processed */
847
848 } /* end of processing time series when fsig is positive */
849
850 /* put dar[] time series (possibly edited above) into output bricks */
851
852 for( iv=0 ; iv < nuse ; iv++ ){
853 zar = DSET_ARRAY(oset,iv+ignore) ; /* output brick */
854 zar[ii] = dar[iv] ; /* original or mutated data */
855 }
856
857 } /* end of loop over voxels #ii */
858
859 #pragma omp critical (DESPIKE_malloc)
860 { free(fit); free(ssp); free(fitar); free(var); free(dar); free(far);
861 if( do_NEW ) free(NEW_wks) ; }
862
863 } /* end OpenMP */
864 AFNI_OMP_END ;
865
866 #ifdef USE_OMP
867 if( verb ) fprintf(stderr,"\n") ;
868 #endif
869 ctim = NI_clock_time() - ctim ;
870 INFO_message( "Elapsed despike time = %s" , nice_time_string(ctim) ) ;
871 if( ctim > 34567 && !do_NEW )
872 ININFO_message("That was SLOW -- try the '-NEW' option for a speedup") ;
873
874 #ifdef USE_OMP
875 if( verb ) fprintf(stderr,"\n") ;
876 #endif
877
878 /*--- finish up ---*/
879
880 if( do_NEW ) mri_free(NEW_psinv) ;
881
882 DSET_delete(dset) ; /* delete input dataset */
883
884 if( verb ){
885 if( nproc > 0 ){
886 pspike = (100.0*nspike)/nproc ;
887 pbig = (100.0*nbig )/nproc ;
888 INFO_message("FINAL: %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]",
889 nproc,nspike,pspike,nbig,pbig ) ;
890 } else {
891 INFO_message("FINAL: no good voxels found to process!!??") ;
892 }
893 }
894
895 /* write results */
896
897 DSET_write(oset) ;
898 if( verb ) WROTE_DSET(oset) ;
899 DSET_delete(oset) ;
900
901 if( tset != NULL ){
902 DSET_write(tset) ;
903 if( verb ) WROTE_DSET(tset) ;
904 DSET_delete(tset) ;
905 }
906
907 exit( THD_get_write_error_count() ) ;
908 }
909