1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #endif
5
6 /*--------------------------------------------------------------------------*/
7 /* adaptive = downweight things a long ways from median */
8
adaptive_weighted_mean(int num,float * x)9 float adaptive_weighted_mean( int num , float *x )
10 {
11 float med,mad, wt,wsum, xsum ; int ii ;
12
13 ENTRY("adaptive_weighted_mean") ;
14
15 if( num <= 0 || x == NULL ) RETURN (0.0f) ;
16 else if( num == 1 ) RETURN (x[0]) ;
17 else if( num == 2 ) RETURN (0.5f*(x[0]+x[1])) ;
18
19 qmedmad_float( num , x , &med , &mad ) ;
20 if( mad <= 0.0f ) RETURN (med) ;
21
22 wsum = xsum = 0.0f ; mad = 0.56789f / mad ;
23 for( ii=0 ; ii < num ; ii++ ){
24 wt = mad*fabsf(x[ii]-med); wt = 1.0f / (1.0f+wt*wt*wt); wsum += wt;
25 xsum += wt * x[ii] ;
26 }
27 RETURN (xsum/wsum) ;
28 }
29
30 /*--------------------------------------------------------------------------*/
31
32 static int nXX=0 , nHH_adapt=0 ;
33 static float *aXX=NULL ;
34 static float *tvv=NULL , *uvv=NULL ;
35 static int ntv=0 ;
36
setup_adaptive_filter(int hwid,int num)37 void setup_adaptive_filter( int hwid , int num )
38 {
39 ENTRY("setup_adaptive_filter") ;
40 if( hwid > 0 && num > 0 && num < 4*hwid+1 ){
41 ERROR_message(
42 "Illegal setup_adaptive_filter: 4*hwid+1=%d > num=%d",4*hwid+1,num) ;
43 EXRETURN ;
44 }
45 if( hwid > 0 ){
46 nHH_adapt = hwid ; nXX = 2*nHH_adapt+1 ;
47 aXX = (float *)realloc(aXX,sizeof(float)*nXX) ;
48 } else if( hwid <= 0 && nHH_adapt > 0 ){ /* cleanup */
49 nHH_adapt = nXX = 0 ;
50 if( aXX != NULL ){ free(aXX) ; aXX = NULL ; }
51 }
52 if( hwid > 0 && num > 0 ){
53 ntv = num + nXX ;
54 tvv = (float *)realloc(tvv,sizeof(float)*ntv) ;
55 uvv = (float *)realloc(uvv,sizeof(float)*ntv) ;
56 } else if( hwid <= 0 || num <= 0 ){
57 if( tvv != NULL ) free(tvv) ;
58 if( uvv != NULL ) free(uvv) ;
59 tvv = uvv = NULL ; ntv = 0 ;
60 }
61 EXRETURN ;
62 }
63
64 /*--------------------------------------------------------------------------*/
65
66 static int dosub=0 ;
67
FILTER_adaptive(int num,float * vec)68 void FILTER_adaptive( int num , float *vec )
69 {
70 int ii,jj,kk , nn1=num-1 , nn2=num-2 , nph = num+nHH_adapt ;
71
72 ENTRY("FILTER_adaptive") ;
73
74 if( num < 2 || vec == NULL || nHH_adapt <= 0 ) EXRETURN ;
75
76 if( aXX == NULL ){ aXX = (float *)malloc(sizeof(float)*nXX) ; }
77
78 if( ntv < num+nXX ){
79 ntv = num+nXX ;
80 tvv = (float *)realloc(tvv,sizeof(float)*ntv) ;
81 uvv = (float *)realloc(uvv,sizeof(float)*ntv) ;
82 }
83
84 memcpy( tvv+nHH_adapt , vec , sizeof(float)*num ) ;
85 for( ii=0 ; ii < nHH_adapt ; ii++ ){ /* end reflections */
86 tvv[ii] = vec[nHH_adapt-ii] ;
87 tvv[nph+ii] = vec[nn2-ii] ;
88 }
89
90 for( ii=0 ; ii < num ; ii++ ){
91 memcpy( aXX , tvv+ii , sizeof(float)*nXX ) ;
92 uvv[ii] = adaptive_weighted_mean( nXX , aXX ) ;
93 }
94
95 if( dosub ){
96 for( ii=0 ; ii < num ; ii++ ) vec[ii] -= uvv[ii] ;
97 } else {
98 memcpy(vec,uvv,sizeof(float)*num) ;
99 }
100
101 EXRETURN ;
102 }
103
104 /*--------------------------------------------------------------------------*/
105
106 static int polort=-1 ;
107
FILTER_polort(int num,float * vec)108 void FILTER_polort( int num , float *vec )
109 {
110 THD_generic_detrend_LSQ( num , vec , polort , 0,NULL,NULL ) ;
111 return ;
112 }
113
114 /*----------------------------------------------------------------------------*/
115 /* Despiking filter */
116
117 #undef SWAP
118 #define SWAP(x,y) (temp=x,x=y,y=temp)
119
120 #undef SORT2
121 #define SORT2(a,b) if(a>b) SWAP(a,b)
122
123 static float *deswks = NULL ;
124
125 /*--- fast median of 9 values ---*/
126
median9f(float * p)127 static INLINE float median9f(float *p)
128 {
129 register float temp ;
130 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
131 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
132 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
133 SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
134 SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
135 SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
136 SORT2(p[4],p[2]) ; return(p[4]) ;
137 }
138
139 /*--- get the local median and MAD of values vec[j-4 .. j+4] ---*/
140
141 #undef mead9
142 #define mead9(j) \
143 { float qqq[9] ; int jj = (j)-4 ; \
144 if( jj < 0 ) jj = 0; else if( jj+9 > num ) jj = num-9; \
145 qqq[0] = vec[jj+0]; qqq[1] = vec[jj+1]; qqq[2] = vec[jj+2]; \
146 qqq[3] = vec[jj+3]; qqq[4] = vec[jj+4]; qqq[5] = vec[jj+5]; \
147 qqq[6] = vec[jj+6]; qqq[7] = vec[jj+7]; qqq[8] = vec[jj+8]; \
148 med = median9f(qqq); qqq[0] = fabsf(qqq[0]-med); \
149 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
150 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
151 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
152 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
153 mad = median9f(qqq); }
154
155 /*-------------------------------------------------------------------------*/
156 /*! Remove spikes from a time series, in a very simplistic way.
157 Return value is the number of spikes that were squashed [RWCox].
158 *//*-----------------------------------------------------------------------*/
159
DES_despike9(int num,float * vec)160 void DES_despike9( int num , float *vec )
161 {
162 int ii , nsp ; float *zma,*zme , med,mad,val ;
163
164 if( num < 9 || vec == NULL ) return ; /* nada */
165
166 if( deswks == NULL ) deswks = (float *)malloc(sizeof(float)*(4*num)) ;
167
168 zme = deswks ; zma = zme + num ;
169
170 for( ii=0 ; ii < num ; ii++ ){
171 mead9(ii) ; zme[ii] = med ; zma[ii] = mad ;
172 }
173 mad = qmed_float(num,zma) ;
174 if( mad <= 0.0f ) return ;
175 mad *= 6.789f ; /* threshold value */
176
177 for( nsp=ii=0 ; ii < num ; ii++ )
178 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
179
180 return ;
181 }
182 #undef mead9
183
184 /*----------------------------------------------------------------------------*/
185 /* Similar code to the above, but for spans of length 25 instead of 9 */
186
median25f(float * p)187 static INLINE float median25f(float *p)
188 {
189 register float temp ;
190 SORT2(p[0], p[1]) ; SORT2(p[3], p[4]) ; SORT2(p[2], p[4]) ;
191 SORT2(p[2], p[3]) ; SORT2(p[6], p[7]) ; SORT2(p[5], p[7]) ;
192 SORT2(p[5], p[6]) ; SORT2(p[9], p[10]) ; SORT2(p[8], p[10]) ;
193 SORT2(p[8], p[9]) ; SORT2(p[12], p[13]) ; SORT2(p[11], p[13]) ;
194 SORT2(p[11], p[12]) ; SORT2(p[15], p[16]) ; SORT2(p[14], p[16]) ;
195 SORT2(p[14], p[15]) ; SORT2(p[18], p[19]) ; SORT2(p[17], p[19]) ;
196 SORT2(p[17], p[18]) ; SORT2(p[21], p[22]) ; SORT2(p[20], p[22]) ;
197 SORT2(p[20], p[21]) ; SORT2(p[23], p[24]) ; SORT2(p[2], p[5]) ;
198 SORT2(p[3], p[6]) ; SORT2(p[0], p[6]) ; SORT2(p[0], p[3]) ;
199 SORT2(p[4], p[7]) ; SORT2(p[1], p[7]) ; SORT2(p[1], p[4]) ;
200 SORT2(p[11], p[14]) ; SORT2(p[8], p[14]) ; SORT2(p[8], p[11]) ;
201 SORT2(p[12], p[15]) ; SORT2(p[9], p[15]) ; SORT2(p[9], p[12]) ;
202 SORT2(p[13], p[16]) ; SORT2(p[10], p[16]) ; SORT2(p[10], p[13]) ;
203 SORT2(p[20], p[23]) ; SORT2(p[17], p[23]) ; SORT2(p[17], p[20]) ;
204 SORT2(p[21], p[24]) ; SORT2(p[18], p[24]) ; SORT2(p[18], p[21]) ;
205 SORT2(p[19], p[22]) ; SORT2(p[8], p[17]) ; SORT2(p[9], p[18]) ;
206 SORT2(p[0], p[18]) ; SORT2(p[0], p[9]) ; SORT2(p[10], p[19]) ;
207 SORT2(p[1], p[19]) ; SORT2(p[1], p[10]) ; SORT2(p[11], p[20]) ;
208 SORT2(p[2], p[20]) ; SORT2(p[2], p[11]) ; SORT2(p[12], p[21]) ;
209 SORT2(p[3], p[21]) ; SORT2(p[3], p[12]) ; SORT2(p[13], p[22]) ;
210 SORT2(p[4], p[22]) ; SORT2(p[4], p[13]) ; SORT2(p[14], p[23]) ;
211 SORT2(p[5], p[23]) ; SORT2(p[5], p[14]) ; SORT2(p[15], p[24]) ;
212 SORT2(p[6], p[24]) ; SORT2(p[6], p[15]) ; SORT2(p[7], p[16]) ;
213 SORT2(p[7], p[19]) ; SORT2(p[13], p[21]) ; SORT2(p[15], p[23]) ;
214 SORT2(p[7], p[13]) ; SORT2(p[7], p[15]) ; SORT2(p[1], p[9]) ;
215 SORT2(p[3], p[11]) ; SORT2(p[5], p[17]) ; SORT2(p[11], p[17]) ;
216 SORT2(p[9], p[17]) ; SORT2(p[4], p[10]) ; SORT2(p[6], p[12]) ;
217 SORT2(p[7], p[14]) ; SORT2(p[4], p[6]) ; SORT2(p[4], p[7]) ;
218 SORT2(p[12], p[14]) ; SORT2(p[10], p[14]) ; SORT2(p[6], p[7]) ;
219 SORT2(p[10], p[12]) ; SORT2(p[6], p[10]) ; SORT2(p[6], p[17]) ;
220 SORT2(p[12], p[17]) ; SORT2(p[7], p[17]) ; SORT2(p[7], p[10]) ;
221 SORT2(p[12], p[18]) ; SORT2(p[7], p[12]) ; SORT2(p[10], p[18]) ;
222 SORT2(p[12], p[20]) ; SORT2(p[10], p[20]) ; SORT2(p[10], p[12]) ;
223 return (p[12]);
224 }
225
226 /*--- get the local median and MAD of values vec[j-12 .. j+12] ---*/
227
228 #undef mead25
229 #define mead25(j) \
230 { float qqq[25] ; int jj=(j)-12 ; register int pp; \
231 if( jj < 0 ) jj = 0; else if( jj+25 > num ) jj = num-25; \
232 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = vec[jj+pp] ; \
233 med = median25f(qqq) ; \
234 for( pp=0 ; pp < 25 ; pp++ ) qqq[pp] = fabsf(qqq[pp]-med) ; \
235 mad = median25f(qqq); }
236
DES_despike25(int num,float * vec)237 void DES_despike25( int num , float *vec )
238 {
239 int ii , nsp ; float *zma,*zme , med,mad,val ;
240
241 if( deswks == NULL ) deswks = (float *)malloc(sizeof(float)*(4*num)) ;
242
243 if( vec == NULL ) return ;
244 if( num < 25 ) { DES_despike9(num,vec) ; return ; }
245
246 zme = deswks ; zma = zme + num ;
247
248 for( ii=0 ; ii < num ; ii++ ){
249 mead25(ii) ; zme[ii] = med ; zma[ii] = mad ;
250 }
251 mad = qmed_float(num,zma) ;
252 if( mad <= 0.0f ) return ;
253 mad *= 6.234f ; /* threshold value */
254
255 for( nsp=ii=0 ; ii < num ; ii++ )
256 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
257
258 return ;
259 }
260 #undef mead25
261 #undef SORT2
262 #undef SWAP
263
264 /*--------------- Generalized to arbitrary widths ----------------*/
265
266 #undef meadHH
267 #define meadHH(j) \
268 { int jj=(j)-hwid ; register int pp; \
269 if( jj < 0 ) jj = 0; else if( jj+nqqq > num ) jj = num-nqqq; \
270 for( pp=0 ; pp < nqqq ; pp++ ) qqq[pp] = vec[jj+pp] ; \
271 qmedmad_float(nqqq,qqq,&med,&mad) ; \
272 }
273
DES_despikeHH(int num,float * vec,int hwid)274 void DES_despikeHH( int num , float *vec , int hwid )
275 {
276 int ii , nsp ; float *zma,*zme , med,mad,val ;
277 static float *qqq=NULL ; static int nhold=0 , nqqq=0 ;
278
279 ENTRY("DES_despikeHH") ;
280
281 if( vec == NULL || num < 9 || hwid < 4 ) EXRETURN ;
282
283 if( deswks == NULL ) deswks = (float *)malloc(sizeof(float)*(4*num)) ;
284
285 if( hwid == 12 ){
286 DES_despike25(num,vec) ; EXRETURN ;
287 } else if( hwid == 4 ){
288 DES_despike9(num,vec) ; EXRETURN ;
289 }
290
291 if( hwid > nhold ){
292 nhold = hwid ;
293 qqq = (float *)realloc(qqq,sizeof(float)*(2*nhold+1)) ;
294 }
295 nqqq = 2*hwid+1 ;
296
297 zme = deswks ; zma = zme + num ;
298
299 for( ii=0 ; ii < num ; ii++ ){
300 meadHH(ii) ; zme[ii] = med ; zma[ii] = mad ;
301 }
302 mad = qmed_float(num,zma) ;
303 if( mad <= 0.0f ) EXRETURN ;
304 mad *= 6.234f ; /* threshold value */
305
306 for( nsp=ii=0 ; ii < num ; ii++ )
307 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
308
309 EXRETURN ;
310 }
311 #undef meadHH
312 #undef SORT2
313 #undef SWAP
314
315 static int nHH_despike=0 ;
316
FILTER_despikeHH(int num,float * vec)317 void FILTER_despikeHH( int num , float *vec )
318 {
319 DES_despikeHH( num , vec , nHH_despike ) ;
320 return ;
321 }
322
323 /*--------------------------------------------------------------------------*/
324 /* Array of filter functions to apply to time series */
325
326 static int nffunc=0 ;
327 static generic_func **ffunc=NULL ;
328
329 #define ADD_FILTER(ff) \
330 do{ ffunc = (generic_func **)realloc(ffunc,sizeof(generic_func)*(nffunc+1)) ; \
331 ffunc[nffunc++] = (generic_func *)(ff) ; \
332 } while(0)
333
334 /*----------------------------------------------------------------------------*/
335
336 static int vstep = 6666 ;
337
vstep_print(void)338 static void vstep_print(void) /* pacifier */
339 {
340 static char xx[10] = "0123456789" ; static int vn=0 ;
341 fprintf(stderr , "%c" , xx[vn%10] ) ;
342 if( vn%10 == 9) fprintf(stderr,".") ; vn++ ;
343 }
344
345 /*--------------------------------------------------------------------------*/
346
FILTER_tsfunc(double tzero,double tdelta,int npts,float * ts,double ts_mean,double ts_slope,void * ud,int nbriks,float * val)347 void FILTER_tsfunc( double tzero , double tdelta ,
348 int npts , float *ts , double ts_mean ,
349 double ts_slope , void *ud , int nbriks, float *val )
350 {
351 int qq ; static int ncalls=0 ;
352
353 if( ts == NULL || val == NULL ) return ; /* setup/cleanup calls */
354
355 memcpy(val,ts,sizeof(float)*npts) ;
356
357 for( qq=1 ; qq < npts ; qq++ ) if( val[qq] != val[0] ) break ;
358 if( qq == npts ) return ;
359
360 ncalls++ ; if( ncalls%vstep == 77 ) vstep_print() ;
361
362 for( qq=0 ; qq < nffunc ; qq++ )
363 AFNI_CALL_VOID_2ARG( ffunc[qq] , int,npts , float *,val ) ;
364
365 return ;
366 }
367
368 /*--------------------------------------------------------------------------*/
369
main(int argc,char * argv[])370 int main( int argc , char *argv[] )
371 {
372 THD_3dim_dataset *old_dset=NULL , *new_dset=NULL ;
373 char *prefix = "Filtered" ;
374 int hwid=0 ;
375 int nvals , nopt ;
376
377 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
378 printf(
379 "\n"
380 "3dTfilter takes as input a dataset, filters the time series in\n"
381 "each voxel as ordered by the user, and outputs a new dataset.\n"
382 "The data in each voxel is processed separately.\n"
383 "\n"
384 "The user (you?) specifies the filter functions to apply.\n"
385 "They are applied in the order given on the command line:\n"
386 " -filter rank -filter adaptive:7\n"
387 "means to do the following operations\n"
388 " (1) turn the data into ranks\n"
389 " (2) apply the adaptive mean filter to the ranks\n"
390 "\n"
391 "Notes:\n"
392 "------\n"
393 "** This program is a work in progress, and more capabilities\n"
394 " will be added as time allows, as the need arises, and as\n"
395 " the author's whims bubble to the surface of his febrile brain.\n"
396 "\n"
397 "** This program is for people who have Sisu.\n"
398 "\n"
399 "Options:\n"
400 "--------\n"
401 "\n"
402 " -input inputdataset\n"
403 "\n"
404 " -prefix outputdataset\n"
405 "\n"
406 " -filter FunctionName\n"
407 " At least one '-filter' option is required!\n"
408 " The FunctionName values that you can give are:\n"
409 "\n"
410 " rank = smallest value is replaced by 0,\n"
411 " next smallest value by 1, and so forth.\n"
412 " ** This filter is pretty useless.\n"
413 "\n"
414 " adaptive:H = adaptive mean filter with half-width of\n"
415 " 'H' time points (H > 0).\n"
416 " ** At most one 'adaptive' filter can be used!\n"
417 " ** The filter 'footprint' is 2*H+1 points.\n"
418 " ** This filter does local smoothing over the\n"
419 " 'footprint', with values far away from\n"
420 " the local median being weighted less.\n"
421 "\n"
422 " adetrend:H = apply adaptive mean filter with half-width\n"
423 " of 'H' time points to get a local baseline,\n"
424 " then subtract this baseline from the actual\n"
425 " data, to provide an adaptive detrending.\n"
426 " ** At most one 'adaptive' OR 'adetrend' filter\n"
427 " can be used.\n"
428 "\n"
429 " despike = apply the 'NEW25' despiking algorithm, as in\n"
430 " program 3dDespike.\n"
431 "\n"
432 " despike:H = apply the despiking algorithm over a window\n"
433 " of half-with 'H' time points (667 > H > 3).\n"
434 " ** H=12 is the same as 'despike'.\n"
435 " ** At most one 'despike' filter can be used.\n"
436 "\n"
437 " detrend:P = (least squares) detrend with polynomials of up\n"
438 " order 'P' for P=0, 1, 2, ....\n"
439 " ** At most one 'detrend' filter can be used!\n"
440 " ** You can use both '-adetrend' and '-detrend',\n"
441 " but I don't know why you would try this.\n"
442 "\n"
443 "Example:\n"
444 "--------\n"
445 " 3dTfilter -input fred.nii -prefix fred.af.nii -filter adaptive:7\n"
446 "\n"
447 "-------\n"
448 "Author: The Programmer with No Name\n"
449 "-------\n"
450 "\n"
451 ) ;
452 exit(0) ;
453 }
454
455 /* bureaucracy */
456
457 mainENTRY("3dTfilter main"); machdep(); AFNI_logger("3dTfilter",argc,argv);
458 PRINT_VERSION("3dTfilter"); AUTHOR("Thorby Baslim");
459
460
461 /*--- scan command line for options ---*/
462
463 nopt = 1 ;
464 while( nopt < argc && argv[nopt][0] == '-' ){
465
466 /*-- prefix --*/
467
468 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
469 if( ++nopt >= argc ) ERROR_exit("%s needs an argument!",argv[nopt-1]);
470 prefix = strdup(argv[nopt]) ;
471 if( !THD_filename_ok(prefix) )
472 ERROR_exit("%s is not a valid prefix!",prefix);
473 nopt++ ; continue ;
474 }
475
476 if( strcasecmp(argv[nopt],"-input") == 0 ||
477 strcasecmp(argv[nopt],"-inset") == 0 ){
478 if( ++nopt >= argc ) ERROR_exit("%s needs an argument!",argv[nopt-1]);
479 if( old_dset != NULL ) ERROR_exit("you can't have 2 input datasets!") ;
480 old_dset = THD_open_dataset(argv[nopt]) ;
481 CHECK_OPEN_ERROR(old_dset,argv[nopt]) ;
482 nopt++ ; continue ;
483 }
484
485 if( strcasecmp(argv[nopt],"-filter") == 0 ){
486 if( ++nopt >= argc ) ERROR_exit("%s needs an argument!",argv[nopt-1]);
487
488 if( strcasecmp(argv[nopt],"rank") == 0 ){
489 ADD_FILTER(rank_order_float) ; /* external function */
490 INFO_message("Filter #%d = rank",nffunc) ;
491
492 } else if( strncasecmp(argv[nopt],"adaptive:",9) == 0 ){
493 char *cpt=argv[nopt]+9 ;
494 if( hwid > 0 )
495 ERROR_exit("You can't use more than one 'adaptive' filter :(") ;
496 if( !isdigit(*cpt) )
497 ERROR_exit("'%s' is not a valid 'adaptive' filter name",argv[nopt]) ;
498 hwid = (int)strtod(cpt,NULL) ;
499 if( hwid > 29 )
500 WARNING_message("Very long filter '%s' will be very slow",argv[nopt]) ;
501 else if( hwid <= 0 )
502 ERROR_exit("'%s' is not a legal 'adaptive' filter name",argv[nopt]) ;
503 ADD_FILTER(FILTER_adaptive) ; dosub = 0 ;
504 INFO_message("Filter #%d = adaptive:%d",nffunc,hwid) ;
505
506 } else if( strncasecmp(argv[nopt],"adetrend:",9) == 0 ){
507 char *cpt=argv[nopt]+9 ;
508 if( hwid > 0 )
509 ERROR_exit("You can't use more than one 'adaptive' filter :(") ;
510 if( !isdigit(*cpt) )
511 ERROR_exit("'%s' is not a valid 'adetrend' filter name",argv[nopt]) ;
512 hwid = (int)strtod(cpt,NULL) ;
513 if( hwid > 29 )
514 WARNING_message("Very long filter '%s' will be very slow",argv[nopt]) ;
515 else if( hwid <= 0 )
516 ERROR_exit("'%s' is not a legal 'adaptive' filter name",argv[nopt]) ;
517 ADD_FILTER(FILTER_adaptive) ; dosub = 1 ;
518 INFO_message("Filter #%d = adetrend:%d",nffunc,hwid) ;
519
520 } else if( strncasecmp(argv[nopt],"detrend:",8) == 0 ){
521 char *cpt=argv[nopt]+8 ;
522 if( polort > 0 )
523 ERROR_exit("You can't use more than one 'detrend' filter :(") ;
524 if( !isdigit(*cpt) )
525 ERROR_exit("'%s' is not a valid 'detrend' filter name",argv[nopt]) ;
526 polort = (int)strtod(cpt,NULL) ;
527 if( polort < 0 )
528 ERROR_exit("'%s' is not a legal 'detrend' filter name",argv[nopt]) ;
529 ADD_FILTER(FILTER_polort) ;
530 INFO_message("Filter #%d = detrend:%d",nffunc,polort) ;
531
532 } else if( strcasecmp(argv[nopt],"despike") == 0 ){
533 if( nHH_despike > 0 )
534 ERROR_exit("You can't use more than one 'despike' filter :(") ;
535 ADD_FILTER(FILTER_despikeHH) ; nHH_despike = 12 ;
536 INFO_message("Filter #%d = despike:12",nffunc) ;
537
538 } else if( strncasecmp(argv[nopt],"despike:",8) == 0 ){
539 char *cpt=argv[nopt]+8 ;
540 if( nHH_despike > 0 )
541 ERROR_exit("You can't use more than one 'despike' filter :(") ;
542 if( !isdigit(*cpt) )
543 ERROR_exit("'%s' is not a valid 'despike' filter name",argv[nopt]) ;
544 nHH_despike = (int)strtod(cpt,NULL) ;
545 if( nHH_despike < 4 || nHH_despike > 666 )
546 ERROR_exit("'%s' is not a legal 'despike' filter name",argv[nopt]) ;
547 if( nHH_despike > 39 )
548 WARNING_message("Very long filter '%s' will be very slow",argv[nopt]) ;
549 ADD_FILTER(FILTER_despikeHH) ;
550 INFO_message("Filter #%d = despike:%d",nffunc,nHH_despike) ;
551
552 } else {
553 ERROR_exit("Unkown filter type '%s'",argv[nopt]) ;
554 }
555 nopt++ ; continue ;
556 }
557
558 ERROR_exit("Unknown option: '%s'",argv[nopt]) ;
559 }
560
561 if( nffunc == 0 ) ERROR_exit("No -filter options given !? :(") ;
562
563 if( old_dset == NULL ){
564 if( nopt >= argc ) ERROR_exit("no input dataset?") ;
565 old_dset = THD_open_dataset(argv[nopt]) ;
566 CHECK_OPEN_ERROR(old_dset,argv[nopt]) ;
567 }
568
569 nvals = DSET_NVALS(old_dset) ;
570 if( nvals < 9 )
571 ERROR_exit("Input dataset too short to filter: nvals=%d < 9",nvals) ;
572
573 if( 4*hwid+1 > nvals )
574 ERROR_exit(
575 "adaptive filter half-width %d is too big for data length %d",
576 hwid,nvals) ;
577
578 if( 4*nHH_despike > nvals )
579 ERROR_exit(
580 "despike filter half-width %d is too big for data length %d",
581 nHH_despike,nvals) ;
582
583 if( hwid > 0 ) setup_adaptive_filter( hwid , nvals ) ;
584
585 INFO_message("Load input dataset") ;
586
587 DSET_load(old_dset) ; CHECK_LOAD_ERROR(old_dset) ;
588
589 /** do the work **/
590
591 fprintf(stderr,"++ Voxel processing: ") ;
592 vstep = DSET_NVOX(old_dset) / 50 ;
593
594 new_dset = MAKER_4D_to_typed_fbuc(
595 old_dset , /* input dataset */
596 prefix , /* output prefix */
597 MRI_float , /* output datum */
598 0 , /* ignore count */
599 0 , /* don't detrend */
600 nvals , /* number of briks */
601 FILTER_tsfunc , /* timeseries processor */
602 NULL, /* data for tsfunc */
603 NULL, /* mask */
604 0 /* Allow auto scaling of output */
605 ) ;
606
607 fprintf(stderr,"!\n") ;
608
609 DSET_unload(old_dset) ;
610
611 if( new_dset != NULL ){
612 tross_Copy_History( old_dset , new_dset ) ;
613 tross_Make_History( "3dTfilter" , argc,argv , new_dset ) ;
614 if( DSET_NUM_TIMES(old_dset) > 1 )
615 EDIT_dset_items( new_dset ,
616 ADN_ntt , DSET_NVALS(old_dset) ,
617 ADN_ttorg , DSET_TIMEORIGIN(old_dset) ,
618 ADN_ttdel , DSET_TR(old_dset) ,
619 ADN_tunits , UNITS_SEC_TYPE ,
620 NULL ) ;
621 DSET_write( new_dset ) ; WROTE_DSET( new_dset ) ;
622 } else {
623 ERROR_exit("Unable to compute output dataset!\n") ;
624 }
625
626 exit(0) ;
627 }
628