1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 void osfilt3_func( int num , float *vec ) ;
10 void median3_func( int num , float *vec ) ;
11 void linear3_func( int num , float *vec ) ;
12
13 static int Nadapt=9 , Nadhalf=4 ;
14 void adaptN_func ( int num , float *vec ) ; /* 03 Oct 2014 */
15
16 /*-- 01 & 03 Mar 2001: linear filtering functions --*/
17
18 void linear_filter_extend( int , float * , int , float * ) ;
19 void linear_filter_zero ( int , float * , int , float * ) ;
20 void linear_filter_trend ( int , float * , int , float * ) ;
21 float * hamming_window ( int ) ;
22 float * blackman_window( int ) ;
23 float * custom_filter( int ) ;
24 char custom_file[256] = "\0"; /*KRH 01/03, needed for custom filter */
25 int custom_ntaps = 0;
26
27 static float af=0.15 , bf=0.70 , cf=0.15 ;
28
29 /*--------------------------------------------------------------------------*/
30
main(int argc,char * argv[])31 int main( int argc , char *argv[] )
32 {
33 void (*smth)(int,float *) = linear3_func ; /* default filter */
34
35 char prefix[256] = "smooth" ;
36 int new_datum = ILLEGAL_TYPE , old_datum ;
37 int nopt ;
38 THD_3dim_dataset *old_dset , *new_dset ;
39 int ii,kk , nxyz , ntime , use_fac , ityp , nbytes ;
40 void *new_brick ;
41
42 byte **bptr = NULL ; /* one of these will be the array of */
43 short **sptr = NULL ; /* pointers to input dataset sub-bricks */
44 float **fptr = NULL ; /* (depending on input datum type) */
45
46 byte **new_bptr = NULL ; /* one of these will be the array of */
47 short **new_sptr = NULL ; /* pointers to output dataset sub-bricks */
48 float **new_fptr = NULL ; /* (depending on output datum type) */
49
50 float *fxar = NULL ; /* array loaded from input dataset */
51 float *fac = NULL ; /* array of brick scaling factors */
52 float *faci = NULL ;
53
54 #define BLACKMAN 1
55 #define HAMMING 2
56 #define CUSTOM 3
57
58 #define EXTEND 77
59 #define ZERO 78
60 #define TREND 79
61
62 int ntap=0 ; /* 01 Mar 2001 */
63 float *ftap=NULL ;
64 int nwin=0,nfil=EXTEND ; /* 03 Mar 2001 */
65
66 void (*lfil)(int,float *,int,float *) = linear_filter_extend ;
67 float * (*lwin)(int) = NULL ;
68
69 /* start of code */
70
71 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
72 printf("Usage: 3dTsmooth [options] dataset\n"
73 "Smooths each voxel time series in a 3D+time dataset and produces\n"
74 "as output a new 3D+time dataset (e.g., lowpass filter in time).\n"
75 "\n"
76 "*** Also see program 3dBandpass ***\n"
77 "\n"
78 "General Options:\n"
79 " -prefix ppp = Sets the prefix of the output dataset to be 'ppp'.\n"
80 " [default = 'smooth']\n"
81 " -datum type = Coerce output dataset to be stored as the given type.\n"
82 " [default = input data type]\n"
83 "\n"
84 "Three Point Filtering Options [07 July 1999]\n"
85 "--------------------------------------------\n"
86 "The following options define the smoothing filter to be used.\n"
87 "All these filters use 3 input points to compute one output point:\n"
88 " Let a = input value before the current point\n"
89 " b = input value at the current point\n"
90 " c = input value after the current point\n"
91 " [at the left end, a=b; at the right end, c=b]\n"
92 "\n"
93 " -lin = 3 point linear filter: 0.15*a + 0.70*b + 0.15*c\n"
94 " [This is the default smoother]\n"
95 " -med = 3 point median filter: median(a,b,c)\n"
96 " -osf = 3 point order statistics filter:\n"
97 " 0.15*min(a,b,c) + 0.70*median(a,b,c) + 0.15*max(a,b,c)\n"
98 "\n"
99 " -3lin m = 3 point linear filter: 0.5*(1-m)*a + m*b + 0.5*(1-m)*c\n"
100 " Here, 'm' is a number strictly between 0 and 1.\n"
101 "\n"
102 "General Linear Filtering Options [03 Mar 2001]\n"
103 "----------------------------------------------\n"
104 " -hamming N = Use N point Hamming or Blackman windows.\n"
105 " -blackman N (N must be odd and bigger than 1.)\n"
106 " -custom coeff_filename.1D (odd # of coefficients must be in a \n"
107 " single column in ASCII file)\n"
108 " (-custom added Jan 2003)\n"
109 " WARNING: If you use long filters, you do NOT want to include the\n"
110 " large early images in the program. Do something like\n"
111 " 3dTsmooth -hamming 13 'fred+orig[4..$]'\n"
112 " to eliminate the first 4 images (say).\n"
113 " The following options determing how the general filters treat\n"
114 " time points before the beginning and after the end:\n"
115 " -EXTEND = BEFORE: use the first value; AFTER: use the last value\n"
116 " -ZERO = BEFORE and AFTER: use zero\n"
117 " -TREND = compute a linear trend, and extrapolate BEFORE and AFTER\n"
118 " The default is -EXTEND. These options do NOT affect the operation\n"
119 " of the 3 point filters described above, which always use -EXTEND.\n"
120 "\n"
121 "Adaptive Mean Filtering option [03 Oct 2014]\n"
122 "--------------------------------------------\n"
123 " -adaptive N = use adaptive mean filtering of width N\n"
124 " (where N must be odd and bigger than 3).\n"
125 " * This filter is similar to the 'AdptMean9'\n"
126 " 1D filter in the AFNI GUI, except that the\n"
127 " end points are treated differently.\n"
128 ) ;
129 printf("\n" MASTER_SHORTHELP_STRING ) ;
130 PRINT_COMPILE_DATE ; exit(0) ;
131 }
132
133 mainENTRY("3dTsmooth main"); machdep(); AFNI_logger("3dTsmooth",argc,argv);
134 PRINT_VERSION("3dTsmooth") ;
135
136 /* parse options */
137
138 nopt = 1 ;
139
140 while( nopt < argc && argv[nopt][0] == '-' ){
141
142 if( strcmp(argv[nopt],"-adaptive") == 0 ){ /* 03 Oct 2014 */
143 if( ++nopt >= argc )
144 ERROR_exit("need value after '%s'",argv[nopt]-1) ;
145 Nadapt = (int)strtod(argv[nopt],NULL) ;
146 if( Nadapt < 5 || Nadapt%2 != 1 )
147 ERROR_exit("illegal value after '%s'",argv[nopt-1]) ;
148 Nadhalf = Nadapt/2 ;
149 INFO_message("Adaptive mean filtering (width=%d) will be done",Nadapt) ;
150 smth = adaptN_func ;
151 nopt++ ; continue ;
152 }
153
154 if( strcmp(argv[nopt],"-EXTEND") == 0 ){ /* 03 Mar 2001 */
155 nfil = EXTEND ; lfil = linear_filter_extend ;
156 nopt++ ; continue ;
157 }
158
159 if( strcmp(argv[nopt],"-ZERO") == 0 ){ /* 03 Mar 2001 */
160 nfil = ZERO ; lfil = linear_filter_zero ;
161 nopt++ ; continue ;
162 }
163
164 if( strcmp(argv[nopt],"-TREND") == 0 ){ /* 03 Mar 2001 */
165 nfil = TREND ; lfil = linear_filter_trend ;
166 nopt++ ; continue ;
167 }
168
169 if( strcmp(argv[nopt],"-hamming") == 0 ){
170 if( ++nopt >= argc ) ERROR_exit("Illegal -hamming: no value!") ;
171 ntap = (int) strtod(argv[nopt],NULL) ;
172 if( ntap < 3 || ntap%2 != 1 ) ERROR_exit("Illegal -hamming: bad value!") ;
173 nwin = HAMMING ; lwin = hamming_window ;
174 INFO_message("Hamming window filter (width=%d) will be done",ntap) ;
175 nopt++ ; continue ;
176 }
177
178 if( strcmp(argv[nopt],"-blackman") == 0 ){
179 if( ++nopt >= argc ) ERROR_exit("Illegal -blackman: no value!") ;
180 ntap = (int) strtod(argv[nopt],NULL) ;
181 if( ntap < 3 || ntap%2 != 1 ) ERROR_exit("Illegal -blackman: bad value!") ;
182 nwin = BLACKMAN ; lwin = blackman_window ;
183 INFO_message("Blackman window filter (width=%d) will be done",ntap) ;
184 nopt++ ; continue ;
185 }
186
187 if( strcmp(argv[nopt],"-custom") == 0 ){
188 if( ++nopt >= argc ) ERROR_exit("Illegal -blackman: no filename!") ;
189 strcpy(custom_file, argv[nopt]) ;
190 nwin = CUSTOM ; lwin = custom_filter ;
191 ntap = 1;
192 INFO_message("Custom window filter will be done") ;
193 nopt++ ; continue ;
194 }
195
196 if( strcmp(argv[nopt],"-prefix") == 0 ){
197 if( ++nopt >= argc ) ERROR_exit("Illegal -prefix: no name!") ;
198 strcpy(prefix,argv[nopt]) ;
199 if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal -prefix: bad name!") ;
200 nopt++ ; continue ;
201 }
202
203 if( strcmp(argv[nopt],"-datum") == 0 ){
204 if( ++nopt >= argc ) ERROR_exit("Illegal -datum: no type name!") ;
205 if( strcmp(argv[nopt],"short") == 0 ){
206 new_datum = MRI_short ;
207 } else if( strcmp(argv[nopt],"float") == 0 ){
208 new_datum = MRI_float ;
209 } else if( strcmp(argv[nopt],"byte") == 0 ){
210 new_datum = MRI_byte ;
211 } else {
212 ERROR_exit("Illegal -datum: bad type name!") ;
213 }
214 nopt++ ; continue ;
215 }
216
217 if( strcmp(argv[nopt],"-lin") == 0 ){
218 bf = 0.70 ; af = cf = 0.15 ;
219 smth = linear3_func ;
220 INFO_message("3 point linear filtering will be done") ;
221 nopt++ ; continue ;
222 }
223
224 if( strcmp(argv[nopt],"-med") == 0 ){
225 smth = median3_func ;
226 INFO_message("3 point median filtering will be done") ;
227 nopt++ ; continue ;
228 }
229
230 if( strcmp(argv[nopt],"-osf") == 0 ){
231 smth = osfilt3_func ;
232 INFO_message("3 point order-statistic filtering will be done") ;
233 nopt++ ; continue ;
234 }
235
236 if( strcmp(argv[nopt],"-3lin") == 0 ){
237 if( ++nopt >= argc ) ERROR_exit("Illegal -3lin: no value!") ;
238 bf = strtod( argv[nopt] , NULL ) ;
239 if( bf <= 0.0 || bf >= 1.0 ) ERROR_exit("Illegal -3lin: bad value") ;
240 af = cf = 0.5*(1.0-bf) ;
241 smth = linear3_func ;
242 INFO_message("3 point linear filtering will be done") ;
243 nopt++ ; continue ;
244 }
245
246 ERROR_exit("Unknown option: %s",argv[nopt]) ;
247
248 } /* end of loop over options */
249
250 if( nopt >= argc ) ERROR_exit("No input dataset filename!?") ;
251
252 /*---------- open dataset ----------*/
253
254 old_dset = THD_open_dataset( argv[nopt] ) ;
255 if( old_dset == NULL ) ERROR_exit("Can't open dataset %s",argv[nopt]) ;
256
257 ntime = DSET_NVALS(old_dset) ;
258 nxyz = DSET_NVOX(old_dset) ;
259
260 if( ntime < 4 )
261 ERROR_exit(" Can't smooth dataset with less than 4 time points!") ;
262
263 DSET_load(old_dset) ; CHECK_LOAD_ERROR(old_dset) ;
264
265 old_datum = DSET_BRICK_TYPE(old_dset,0) ;
266 if( new_datum < 0 ) new_datum = old_datum ;
267
268 switch( old_datum ){ /* pointer type depends on input datum type */
269
270 /** create array of pointers into old dataset sub-bricks **/
271
272 /*--------- input is bytes ----------*/
273 /* voxel #i at time #k is bptr[k][i] */
274 /* for i=0..nxyz-1 and k=0..ntime-1. */
275
276 case MRI_byte:
277 bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
278 for( kk=0 ; kk < ntime ; kk++ )
279 bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk) ;
280 break ;
281
282 /*--------- input is shorts ---------*/
283 /* voxel #i at time #k is sptr[k][i] */
284 /* for i=0..nxyz-1 and k=0..ntime-1. */
285
286 case MRI_short:
287 sptr = (short **) malloc( sizeof(short *) * ntime ) ;
288 for( kk=0 ; kk < ntime ; kk++ )
289 sptr[kk] = (short *) DSET_ARRAY(old_dset,kk) ;
290 break ;
291
292 /*--------- input is floats ---------*/
293 /* voxel #i at time #k is fptr[k][i] */
294 /* for i=0..nxyz-1 and k=0..ntime-1. */
295
296 case MRI_float:
297 fptr = (float **) malloc( sizeof(float *) * ntime ) ;
298 for( kk=0 ; kk < ntime ; kk++ )
299 fptr[kk] = (float *) DSET_ARRAY(old_dset,kk) ;
300 break ;
301
302 } /* end of switch on input type */
303
304 /*---- allocate space for 1 voxel timeseries ----*/
305
306 fxar = (float *) malloc( sizeof(float) * ntime ) ; /* voxel timeseries */
307
308 /*--- get scaling factors for sub-bricks ---*/
309
310 fac = (float *) malloc( sizeof(float) * ntime ) ; /* factors */
311
312 use_fac = 0 ;
313 for( kk=0 ; kk < ntime ; kk++ ){
314 fac[kk] = DSET_BRICK_FACTOR(old_dset,kk) ;
315 if( fac[kk] != 0.0 ) use_fac++ ;
316 else fac[kk] = 1.0 ;
317 }
318 if( !use_fac ){
319 free(fac) ; fac = NULL ;
320 } else {
321 faci = (float *) malloc( sizeof(float) * ntime ) ;
322 for( kk=0 ; kk < ntime ; kk++ ) faci[kk] = 1.0 / fac[kk] ;
323 }
324
325 /*---------------------- make a new dataset ----------------------*/
326
327 new_dset = EDIT_empty_copy( old_dset ) ;
328
329 tross_Copy_History( old_dset , new_dset ) ;
330 tross_Make_History( "3dTsmooth" , argc,argv , new_dset ) ;
331
332 /*-- edit some of its internal parameters --*/
333
334 EDIT_dset_items( new_dset ,
335 ADN_prefix , prefix ,
336 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
337 ADN_datum_all , new_datum ,
338 ADN_none ) ;
339
340 /*-- make brick(s) for this dataset --*/
341
342 switch( new_datum ){
343 case MRI_byte:
344 new_bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
345 break ;
346
347 case MRI_short:
348 new_sptr = (short **) malloc( sizeof(short *) * ntime ) ;
349 break ;
350
351 case MRI_float:
352 new_fptr = (float **) malloc( sizeof(float *) * ntime ) ;
353 break ;
354 }
355
356 for( kk=0 ; kk < ntime ; kk++ ){
357 ityp = DSET_BRICK_TYPE(new_dset,kk) ; /* type of data */
358 nbytes = DSET_BRICK_BYTES(new_dset,kk) ; /* how much data */
359 new_brick = malloc( nbytes ) ; /* make room */
360
361 if( new_brick == NULL )
362 ERROR_exit("Can't allocate memory for output dataset!") ;
363
364 EDIT_substitute_brick( new_dset , kk , ityp , new_brick ) ;
365
366 switch( new_datum ){
367 case MRI_byte: new_bptr[kk] = (byte * ) new_brick ; break ;
368 case MRI_short: new_sptr[kk] = (short *) new_brick ; break ;
369 case MRI_float: new_fptr[kk] = (float *) new_brick ; break ;
370 }
371 }
372
373 if( lwin != NULL && ntap > 0 ){ /* 03 Mar 2001 */
374 ftap = lwin(ntap) ;
375 if( lfil == NULL ) lfil = linear_filter_extend ;
376 if( nwin == CUSTOM ) ntap = custom_ntaps ;
377 }
378
379 /*----------------------------------------------------*/
380 /*----- Setup has ended. Now do some real work. -----*/
381
382 /***** loop over voxels *****/
383
384 for( ii=0 ; ii < nxyz ; ii++ ){ /* 1 time series at a time */
385
386 /*** load data from input dataset, depending on type ***/
387
388 switch( old_datum ){
389 case MRI_byte:
390 for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = bptr[kk][ii] ;
391 break ;
392
393 case MRI_short:
394 for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = sptr[kk][ii] ;
395 break ;
396
397 case MRI_float:
398 for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = fptr[kk][ii] ;
399 break ;
400 } /* end of switch over input type */
401
402 if( use_fac )
403 for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] *= fac[kk] ;
404
405 /* do smoothing */
406
407 if( ftap != NULL )
408 lfil( ntap,ftap , ntime,fxar ) ; /* 01 Mar 2001 */
409 else
410 smth( ntime , fxar ) ; /* 3 point smoother */
411
412 /*** put data into output dataset ***/
413
414 switch( new_datum ){
415
416 case MRI_byte:
417 if( use_fac )
418 for( kk=0 ; kk < ntime ; kk++ ) new_bptr[kk][ii] = (byte)(fxar[kk] * faci[kk]) ;
419 else
420 for( kk=0 ; kk < ntime ; kk++ ) new_bptr[kk][ii] = (byte) fxar[kk] ;
421 break ;
422
423 case MRI_short:
424 if( use_fac )
425 for( kk=0 ; kk < ntime ; kk++ ) new_sptr[kk][ii] = (short)(fxar[kk] * faci[kk]) ;
426 else
427 for( kk=0 ; kk < ntime ; kk++ ) new_sptr[kk][ii] = (short) fxar[kk] ;
428 break ;
429
430 case MRI_float:
431 if( use_fac )
432 for( kk=0 ; kk < ntime ; kk++ ) new_fptr[kk][ii] = (float)(fxar[kk] * faci[kk]) ;
433 else
434 for( kk=0 ; kk < ntime ; kk++ ) new_fptr[kk][ii] = (float) fxar[kk] ;
435 break ;
436 }
437 } /* end of loop over voxels */
438
439 DSET_unload(old_dset) ; free(ftap) ;
440
441 if (DSET_write(new_dset) != False) {
442 fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
443 exit(0) ;
444 } else {
445 ERROR_exit("3dTsmooth: Failed to write output dataset :-(" ) ;
446 }
447
448 }
449
450 /*--------------- Order Statistics Filter ----------------*/
451
osfilt3_func(int num,float * vec)452 void osfilt3_func( int num , float *vec )
453 {
454 int ii ;
455 float aa,bb,cc ;
456
457 bb = vec[0] ; cc = vec[1] ; vec[0] = OSFILT(bb,bb,cc) ;
458 for( ii=1 ; ii < num-1 ; ii++ ){
459 aa = bb ; bb = cc ; cc = vec[ii+1] ;
460 vec[ii] = OSFILT(aa,bb,cc) ; /* see mrilib.h */
461 }
462 vec[num-1] = OSFILT(bb,cc,cc) ;
463
464 return ;
465 }
466
467 /*--------------- Median of 3 Filter ----------------*/
468
median3_func(int num,float * vec)469 void median3_func( int num , float *vec )
470 {
471 int ii ;
472 float aa,bb,cc ;
473
474 bb = vec[0] ; cc = vec[1] ;
475 for( ii=1 ; ii < num-1 ; ii++ ){
476 aa = bb ; bb = cc ; cc = vec[ii+1] ;
477 vec[ii] = MEDIAN(aa,bb,cc) ; /* see mrilib.h */
478 }
479
480 return ;
481 }
482
483 /*--------------- Linear Filter ----------------*/
484
485 #define LSUM(a,b,c) af*(a)+bf*(b)+cf*(c)
486
linear3_func(int num,float * vec)487 void linear3_func( int num , float *vec )
488 {
489 int ii ;
490 float aa,bb,cc ;
491
492 bb = vec[0] ; cc = vec[1] ; vec[0] = LSUM(bb,bb,cc) ;
493 for( ii=1 ; ii < num-1 ; ii++ ){
494 aa = bb ; bb = cc ; cc = vec[ii+1] ;
495 vec[ii] = LSUM(aa,bb,cc) ;
496 }
497 vec[num-1] = LSUM(bb,cc,cc) ;
498
499 return ;
500 }
501
502 /*-------------------------------------------------------------------------*/
503
linear_filter_extend(int ntap,float * wt,int npt,float * x)504 void linear_filter_extend( int ntap , float *wt , int npt , float *x )
505 {
506 int ii , nt2=(ntap-1)/2 , jj ;
507 float sum ;
508 static int nfar=0 ;
509 static float *far=NULL ;
510
511 if( npt > nfar ){
512 if(far != NULL) free(far) ;
513 far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
514 }
515
516 #undef XX
517 #define XX(i) ( ((i)<0) ? far[0] : ((i)>npt-1) ? far[npt-1] : far[i] )
518
519 memcpy( far , x , sizeof(float)*npt ) ;
520
521 for( ii=0 ; ii < npt ; ii++ ){
522 for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
523 x[ii] = sum ;
524 }
525
526 return ;
527 }
528
529 /*-------------------------------------------------------------------------*/
530
linear_filter_zero(int ntap,float * wt,int npt,float * x)531 void linear_filter_zero( int ntap , float *wt , int npt , float *x )
532 {
533 int ii , nt2=(ntap-1)/2 , jj ;
534 float sum ;
535 static int nfar=0 ;
536 static float *far=NULL ;
537
538 if( npt > nfar ){
539 if(far != NULL) free(far) ;
540 far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
541 }
542
543 #undef XX
544 #define XX(i) ( ((i)<0 || (i)>npt-1) ? 0 : far[i] )
545
546 memcpy( far , x , sizeof(float)*npt ) ;
547
548 for( ii=0 ; ii < npt ; ii++ ){
549 for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
550 x[ii] = sum ;
551 }
552
553 return ;
554 }
555
556 /*-------------------------------------------------------------------------*/
557
linear_filter_trend(int ntap,float * wt,int npt,float * x)558 void linear_filter_trend( int ntap , float *wt , int npt , float *x )
559 {
560 int ii , nt2=(ntap-1)/2 , jj ;
561 float sum ;
562 static int nfar=0 ;
563 static float *far=NULL ;
564 float a=0.0,b=0.0 ;
565
566 if( npt > nfar ){
567 if(far != NULL) free(far) ;
568 far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
569 }
570
571 #undef XX
572 #define XX(i) ( ((i)<0 || (i)>npt-1) ? (a+b*(i)) : far[i] )
573
574 memcpy( far , x , sizeof(float)*npt ) ;
575
576 get_linear_trend( npt,far , &a,&b ) ; /* cf. thd_detrend.c */
577
578 for( ii=0 ; ii < npt ; ii++ ){
579 for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
580 x[ii] = sum ;
581 }
582
583 return ;
584 }
585
586 /*-------------------------------------------------------------------------*/
587
custom_filter(int dummy)588 float * custom_filter( int dummy )
589 {
590 MRI_IMAGE * filter_data=NULL;
591 float * filter_coefficients = NULL;
592
593 filter_data = mri_read_1D(custom_file);
594 if MRI_IS_1D(filter_data) {
595 custom_ntaps = filter_data->nx;
596 if (custom_ntaps % 2) {
597 filter_coefficients = MRI_FLOAT_PTR(filter_data);
598 } else {
599 /*ERROR, EVEN NUMBER OF FILTER TAPS */
600 }
601 }
602
603 return filter_coefficients ;
604 }
605
606 /*-------------------------------------------------------------------------*/
607
hamming_window(int ntap)608 float * hamming_window( int ntap )
609 {
610 float * wt , tau , t , sum ;
611 int ii , nt2=(ntap-1)/2 ;
612
613 if( ntap < 3 ) return NULL ;
614
615 wt = (float *) calloc(sizeof(float),ntap) ;
616 tau = nt2 + 1.0 ;
617
618 for( sum=0.0,ii=0 ; ii <= 2*nt2 ; ii++ ){
619 t = PI*(ii-nt2)/tau ;
620 wt[ii] = 0.54 + 0.46*cos(t) ;
621 sum += wt[ii] ;
622 }
623 sum = 1.0 / sum ;
624 for( ii=0 ; ii < ntap ; ii++ ) wt[ii] *= sum ;
625
626 return wt ;
627 }
628
629 /*-------------------------------------------------------------------------*/
630
blackman_window(int ntap)631 float * blackman_window( int ntap )
632 {
633 float * wt , tau , t , sum ;
634 int ii , nt2=(ntap-1)/2 ;
635
636 if( ntap < 3 ) return NULL ;
637
638 wt = (float *) calloc(sizeof(float),ntap) ;
639 tau = nt2 + 1.0 ;
640
641 for( sum=0.0,ii=0 ; ii <= 2*nt2 ; ii++ ){
642 t = PI*(ii-nt2)/tau ;
643 wt[ii] = 0.42 + 0.5*cos(t) + 0.08*cos(2*t) ;
644 sum += wt[ii] ;
645 }
646 sum = 1.0 / sum ;
647 for( ii=0 ; ii < ntap ; ii++ ) wt[ii] *= sum ;
648
649 return wt ;
650 }
651
652 /*--------------------------------------------------------------------------*/
653
adaptive_weighted_mean(int num,float * x)654 static float adaptive_weighted_mean( int num , float *x )
655 {
656 float med,mad, wt,wsum, xsum ; int ii ;
657
658 if( num <= 0 || x == NULL ) return (0.0f) ;
659 else if( num == 1 ) return (x[0]) ;
660 else if( num == 2 ) return (0.5f*(x[0]+x[1])) ;
661
662 qmedmad_float( num , x , &med , &mad ) ;
663 if( mad <= 0.0f ) return (med) ;
664
665 wsum = xsum = 0.0f ; mad = 0.4567f / mad ;
666 for( ii=0 ; ii < num ; ii++ ){
667 wt = mad*fabsf(x[ii]-med); wt = 1.0f / (1.0f+wt*wt*wt); wsum += wt;
668 xsum += wt * x[ii] ;
669 }
670 return (xsum/wsum) ;
671 }
672
673 /*--------------------------------------------------------------------------*/
674
adaptN_func(int num,float * vec)675 void adaptN_func( int num , float *vec ) /* 03 Oct 2014 */
676 {
677 static float *x=NULL ;
678 float *nv ; int ii,jj,kk,nn , n1=num-1 ;
679
680 nv = (float *)malloc(sizeof(float)*num) ;
681 if( x == NULL ) x = (float *)malloc(sizeof(float)*Nadapt) ;
682
683 for( ii=0 ; ii < num ; ii++ ){
684
685 for( nn=0,jj=-Nadhalf ; jj <= Nadhalf ; jj++ ){
686 kk = ii+jj ; if( kk < 0 ) continue ; if( kk > n1 ) break ;
687 x[nn++] = vec[kk] ;
688 }
689
690 nv[ii] = adaptive_weighted_mean( nn , x ) ;
691 }
692
693 memcpy(vec,nv,sizeof(float)*num) ; free(nv) ; return ;
694 }
695