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 /**************************************************************************/
10 /***** prototypes for FIR blurring functions at the end of this file ******/
11
12 void fir_blurx( int m, float *wt,int nx, int ny, int nz, float *f ) ;
13 void fir_blury( int m, float *wt,int nx, int ny, int nz, float *f ) ;
14 void fir_blurz( int m, float *wt,int nx, int ny, int nz, float *f ) ;
15
16 static void fir_gaussian_load( int m, float dx, float *ff ) ; /* 03 Apr 2007 */
17
18 #undef FIR_MAX
19 #define FIR_MAX 35 /** max length of FIR filter to use instead of FFTs **/
20
21 #undef DO_INFO
22 #define DO_INFO 0 /** for debugging and timing **/
23
24 /*! If set to 0, EDIT_blur_volume() will not use the fir_blur? functions. */
25
26 static int allow_fir = 1 ;
27
28 /*! Controls whether EDIT_blur_volume() will use the fir_blur? functions. */
29
EDIT_blur_allow_fir(int i)30 void EDIT_blur_allow_fir( int i ){ allow_fir = i; }
31
32 /*************************************************************************/
33 /*! Routine to blur a 3D volume with a Gaussian, using FFTs.
34 If the blurring parameter (sigma) is small enough, will use
35 real-space FIR filter instead. This function actually just
36 calls EDIT_blur_volume_3d() to do the real work.
37 **************************************************************************/
38
EDIT_blur_volume(int nx,int ny,int nz,float dx,float dy,float dz,int ftype,void * vfim,float sigma)39 void EDIT_blur_volume( int nx, int ny, int nz,
40 float dx, float dy, float dz,
41 int ftype , void * vfim , float sigma )
42 {
43 EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, ftype, vfim,
44 sigma, sigma, sigma);
45 }
46
47 /**************************************************************************/
48 /*! The following slightly modified version of EDIT_blur_volume allows
49 independent specification of Gaussian filter widths along the three
50 perpendicular axes.
51 - BDW - 21 Feb 1997
52 - RWC - 04 Feb 2005: use fir_blur? function if sigma? is small
53 - also see EDIT_blur_allow_fir() and FIR_blur_volume_3d()
54 --------------------------------------------------------------------------*/
55
EDIT_blur_volume_3d(int nx,int ny,int nz,float dx,float dy,float dz,int ftype,void * vfim,float sigmax,float sigmay,float sigmaz)56 void EDIT_blur_volume_3d( int nx, int ny, int nz,
57 float dx, float dy, float dz,
58 int ftype , void *vfim ,
59 float sigmax, float sigmay, float sigmaz )
60 {
61 int jj,kk , nxy , base,nby2 ;
62 float dk , aa , k , fac ;
63 register int ii , nup ;
64
65 complex *cx = NULL ;
66 float *gg = NULL ; int ncg=0 ;
67
68 byte *bfim = NULL ; /* pointers to data array vfim */
69 short *sfim = NULL ;
70 float *ffim = NULL ;
71 complex *cfim = NULL ;
72
73 float fbot=0.0,ftop=0.0 ; /* 10 Jan 2003: for clipping results */
74 int nxyz ; /* number of voxels */
75
76 int fir_m , fir_num=0 ; /* 03 Oct 2005: for fir_blur? filtering */
77 float fir_wt[FIR_MAX+1] ;
78 int all_fir=0 ; /* 06 Oct 2005 */
79
80 int as_float = 0 ; /* process as float (to apply FIR) 15 Jun 2009 [rickr] */
81 /* - convert short/byte to float, blur, convert back */
82 int dtype = ftype ; /* processing datatype (ftype, if no conversion) */
83
84 double sfac = AFNI_numenv("AFNI_BLUR_FIRFAC") ;
85 /* <2 -> <=0 to allow for truncated blurs [8 May 2019 rickr] */
86 /* */
87 /* This can be abused for a "fast ANATICOR", for example. */
88 /* Since sigma = 0.4246609 * fwhm, consider using: */
89 /* sfac = 1/(2*.0.4246609) = 1.17741 */
90 /* That number of sigmas should match the half width at half max, */
91 /* which should terminate the blur just after a half height. */
92 /* */
93 /* Or use 2*FWHM and sfac = 1.17741/2 = 0.588705 to make it more flat,*/
94 /* with a min contribution of ~0.84, rather than 0.5, yet limiting */
95 /* the output to the same HWHM radius (e.g. FWHM=80mm with sfac=0.589 */
96 /* results in a fairly flat blur out to a radius of ~20 mm). */
97 if( sfac <= 0.0 ) sfac = 2.5 ;
98
99 /***---------- initialize ----------***/
100
101 ENTRY("EDIT_blur_volume_3d") ;
102
103 if( vfim == NULL ) EXRETURN ; /* no data? */
104
105 if( sigmax <= 0.0 && sigmay <= 0.0 && sigmaz <= 0.0 ) EXRETURN ;
106
107 if( dx <= 0.0 ) dx = 1.0 ; /* 03 Oct 2005: regularize grid sizes */
108 if( dy <= 0.0 ) dy = dx ;
109 if( dz <= 0.0 ) dz = dx ;
110
111 switch( ftype ){ /* cast pointer to correct type */
112 default: EXRETURN ;
113 case MRI_short: sfim = (short *) vfim ; break ;
114 case MRI_float: ffim = (float *) vfim ; break ;
115 case MRI_byte: bfim = (byte *) vfim ; break ;
116 case MRI_complex: cfim = (complex *) vfim ; break ;
117 }
118 nxy = nx * ny ; nxyz = nxy * nz ;
119
120 /*** 10 Jan 2003: find bot and top of data input */
121
122 if( allow_fir ){
123 ii = (int) ceil( sfac * sigmax / dx ) ;
124 jj = (int) ceil( sfac * sigmay / dy ) ;
125 kk = (int) ceil( sfac * sigmaz / dz ) ;
126 if( ii <= FIR_MAX && jj <= FIR_MAX && kk <= FIR_MAX ) all_fir = 1 ;
127
128 /* try to process most data using Finite Impule Response, not with Fourier
129 interpolation (leads to ringing problems both inside and outside a mask)
130 - set AFNI_BLUR_INTS_AS_OLD=YES to use old method 15 Jun 2009 [rickr] */
131 if( (ftype == MRI_short || ftype == MRI_byte) &&
132 ! AFNI_yesenv("AFNI_BLUR_INTS_AS_OLD") ) {
133 ffim = (float *)malloc(nxyz * sizeof(float));
134 as_float = 1; /* flag, redundant but clear */
135 dtype = MRI_float; /* process data as float */
136 }
137 }
138 if( ftype != MRI_float ) all_fir = 0 ; /* 17 Nov 2005: oopsie */
139
140 if( !all_fir ){
141 switch( ftype ){
142 default:
143 fbot = ftop = 0.0 ; /* for complex */
144 break ;
145
146 case MRI_short:
147 fbot = ftop = sfim[0] ;
148 for( ii=1 ; ii < nxyz ; ii++ ) {
149 if( sfim[ii] < fbot ) fbot = sfim[ii] ;
150 else if( sfim[ii] > ftop ) ftop = sfim[ii] ;
151 if( as_float ) ffim[ii] = (float)sfim[ii] ; /* copy data as int */
152 }
153 break ;
154
155 case MRI_float:
156 fbot = ftop = ffim[0] ;
157 for( ii=1 ; ii < nxyz ; ii++ )
158 if( ffim[ii] < fbot ) fbot = ffim[ii] ;
159 else if( ffim[ii] > ftop ) ftop = ffim[ii] ;
160 break ;
161
162 case MRI_byte:
163 fbot = ftop = bfim[0] ;
164 for( ii=1 ; ii < nxyz ; ii++ ) {
165 if( bfim[ii] < fbot ) fbot = bfim[ii] ;
166 else if( bfim[ii] > ftop ) ftop = bfim[ii] ;
167 if( as_float ) ffim[ii] = (float)bfim[ii] ; /* copy data as int */
168 }
169 break ;
170 }
171 }
172
173 /*** do x-direction ***/
174
175 /** 03 Oct 2005: perhaps do the x-blur in real-space? **/
176
177 if( nx < 2 || sigmax <= 0.0 ){
178 STATUS("skipping x blur") ; fir_num++ ; goto DO_Y_BLUR ;
179 }
180
181 fir_m = (int) ceil( sfac * sigmax / dx ) ;
182 if( allow_fir && dtype == MRI_float && fir_m <= FIR_MAX ){
183 STATUS("start x FIR") ;
184 if( fir_m < 1 ) fir_m = 1 ;
185 fir_gaussian_load( fir_m , dx/sigmax , fir_wt ) ;
186 #if 0
187 fac = fir_wt[0] = 1.0f ;
188 for( ii=1 ; ii <= fir_m ; ii++ ){
189 fir_wt[ii] = exp(-0.5*(ii*dx)*(ii*dx)/(sigmax*sigmax)) ;
190 fac += 2.0f * fir_wt[ii] ;
191 }
192 fac = 1.0f / fac ;
193 for( ii=0 ; ii <= fir_m ; ii++ ) fir_wt[ii] *= fac ;
194 #endif
195 fir_blurx( fir_m , fir_wt , nx,ny,nz , ffim ) ;
196 fir_num++ ; goto DO_Y_BLUR ;
197 }
198
199 STATUS("start x FFTs") ;
200
201 aa = sigmax * sigmax * 0.5 ;
202 nup = nx + (int)(3.0 * sigmax / dx) ; /* min FFT length */
203 nup = csfft_nextup_one35(nup) ; nby2 = nup / 2 ;
204
205 if(DO_INFO) INFO_message("x FFT(%d)",nup) ;
206
207 if( nup > ncg ){
208 cx = (complex *)realloc(cx,sizeof(complex)*nup) ;
209 gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ;
210 }
211
212 dk = (2.0*PI) / (nup * dx) ;
213 fac = 1.0 / nup ;
214 gg[0] = fac ;
215 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
216
217 /** July 20: double up on FFTs **/
218 /** Feb 09: extend to other data types besides shorts;
219 doubling up does not apply to complex data! **/
220
221 switch( dtype ){
222 case MRI_short:{
223 register short *qfim ;
224 for( kk=0 ; kk < nz ; kk++ ){
225 for( jj=0 ; jj < ny ; jj+=2 ){
226 base = jj*nx + kk*nxy ;
227 qfim = sfim + base ;
228 if( jj == ny-1 )
229 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
230 else
231 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
232 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
233 csfft_cox( -1 , nup , cx ) ;
234 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
235 csfft_cox( 1 , nup , cx ) ;
236 if( jj == ny-1 )
237 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
238 else
239 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
240 }
241 }
242 }
243 break ;
244
245 case MRI_float:{
246 register float *qfim ;
247 for( kk=0 ; kk < nz ; kk++ ){
248 for( jj=0 ; jj < ny ; jj+=2 ){
249 base = jj*nx + kk*nxy ;
250 qfim = ffim + base ;
251 if( jj == ny-1 )
252 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
253 else
254 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
255 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
256 csfft_cox( -1 , nup , cx ) ;
257 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
258 csfft_cox( 1 , nup , cx ) ;
259 if( jj == ny-1 )
260 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
261 else
262 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
263 }
264 }
265 }
266 break ;
267
268 case MRI_byte:{
269 register byte *qfim ;
270 for( kk=0 ; kk < nz ; kk++ ){
271 for( jj=0 ; jj < ny ; jj+=2 ){
272 base = jj*nx + kk*nxy ;
273 qfim = bfim + base ;
274 if( jj == ny-1 )
275 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
276 else
277 for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
278 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
279 csfft_cox( -1 , nup , cx ) ;
280 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
281 csfft_cox( 1 , nup , cx ) ;
282 if( jj == ny-1 )
283 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
284 else
285 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
286 }
287 }
288 }
289 break ;
290
291 case MRI_complex:{
292 register complex *qfim ;
293 for( kk=0 ; kk < nz ; kk++ ){
294 for( jj=0 ; jj < ny ; jj++ ){
295 base = jj*nx + kk*nxy ;
296 qfim = cfim + base ;
297 for( ii=0 ; ii<nx ; ii++) { cx[ii] = qfim[ii] ; }
298 for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
299 csfft_cox( -1 , nup , cx ) ;
300 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
301 csfft_cox( 1 , nup , cx ) ;
302 for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii] ; }
303 }
304 }
305 }
306 break ;
307 }
308
309 /*** do y-direction ***/
310 DO_Y_BLUR:
311
312 /** 03 Oct 2005: perhaps do the y-blur in real-space? **/
313
314 if( ny < 2 || sigmay <= 0.0 ){
315 STATUS("skip y blur") ; fir_num++ ; goto DO_Z_BLUR ;
316 }
317
318 fir_m = (int) ceil( sfac * sigmay / dy ) ;
319 if( allow_fir && dtype == MRI_float && fir_m <= FIR_MAX ){
320 STATUS("start y FIR") ;
321 if( fir_m < 1 ) fir_m = 1 ;
322 fir_gaussian_load( fir_m , dy/sigmay , fir_wt ) ;
323 #if 0
324 fac = fir_wt[0] = 1.0f ;
325 for( ii=1 ; ii <= fir_m ; ii++ ){
326 fir_wt[ii] = exp(-0.5*(ii*dy)*(ii*dy)/(sigmay*sigmay)) ;
327 fac += 2.0f * fir_wt[ii] ;
328 }
329 fac = 1.0f / fac ;
330 for( ii=0 ; ii <= fir_m ; ii++ ) fir_wt[ii] *= fac ;
331 #endif
332 fir_blury( fir_m , fir_wt , nx,ny,nz , ffim ) ;
333 fir_num++ ; goto DO_Z_BLUR ;
334 }
335
336 STATUS("start y FFTs") ;
337
338 aa = sigmay * sigmay * 0.5 ;
339 nup = ny + (int)(3.0 * sigmay / dy) ; /* min FFT length */
340 nup = csfft_nextup_one35(nup) ; nby2 = nup / 2 ;
341
342 if(DO_INFO) INFO_message("y FFT(%d)",nup) ;
343
344 if( nup > ncg ){
345 cx = (complex *)realloc(cx,sizeof(complex)*nup) ;
346 gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ;
347 }
348
349 dk = (2.0*PI) / (nup * dy) ;
350 fac = 1.0 / nup ;
351 gg[0] = fac ;
352 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
353
354 switch( dtype ){
355 case MRI_short:{
356 register short *qfim ;
357 for( kk=0 ; kk < nz ; kk++ ){
358 for( jj=0 ; jj < nx ; jj+=2 ){
359 base = jj + kk*nxy ;
360 qfim = sfim + base ;
361 if( jj == nx-1 )
362 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
363 else
364 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
365 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
366 csfft_cox( -1 , nup , cx ) ;
367 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
368 csfft_cox( 1 , nup , cx ) ;
369 if( jj == nx-1 )
370 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
371 else
372 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
373 }
374 }
375 }
376 break ;
377
378 case MRI_byte:{
379 register byte *qfim ;
380 for( kk=0 ; kk < nz ; kk++ ){
381 for( jj=0 ; jj < nx ; jj+=2 ){
382 base = jj + kk*nxy ;
383 qfim = bfim + base ;
384 if( jj == nx-1 )
385 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
386 else
387 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
388 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
389 csfft_cox( -1 , nup , cx ) ;
390 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
391 csfft_cox( 1 , nup , cx ) ;
392 if( jj == nx-1 )
393 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
394 else
395 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
396 }
397 }
398 }
399 break ;
400
401 case MRI_float:{
402 register float *qfim ;
403 for( kk=0 ; kk < nz ; kk++ ){
404 for( jj=0 ; jj < nx ; jj+=2 ){
405 base = jj + kk*nxy ;
406 qfim = ffim + base ;
407 if( jj == nx-1 )
408 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
409 else
410 for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
411 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
412 csfft_cox( -1 , nup , cx ) ;
413 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
414 csfft_cox( 1 , nup , cx ) ;
415 if( jj == nx-1 )
416 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
417 else
418 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
419 }
420 }
421 }
422 break ;
423
424 case MRI_complex:{
425 register complex *qfim ;
426 for( kk=0 ; kk < nz ; kk++ ){
427 for( jj=0 ; jj < nx ; jj++ ){
428 base = jj + kk*nxy ;
429 qfim = cfim + base ;
430 for( ii=0 ; ii<ny ; ii++){ cx[ii] = qfim[ii*nx] ; }
431 for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
432 csfft_cox( -1 , nup , cx ) ;
433 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
434 csfft_cox( 1 , nup , cx ) ;
435 for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii] ; }
436 }
437 }
438 }
439 break ;
440 }
441
442 /*** do z-direction ***/
443 DO_Z_BLUR:
444
445 /** 03 Oct 2005: perhaps do the z-blur in real-space? **/
446
447 /* sigmay to sizmaz, noted by Patryk on MB 21 July, 2011 [rickr] */
448 if( nz < 2 || sigmaz <= 0.0 ){
449 STATUS("skip z blur") ; fir_num++ ; goto ALL_DONE_NOW ;
450 }
451
452 fir_m = (int) ceil( sfac * sigmaz / dz ) ;
453 if( allow_fir && dtype == MRI_float && fir_m <= FIR_MAX ){
454 STATUS("start z FIR") ;
455 if( fir_m < 1 ) fir_m = 1 ;
456 fir_gaussian_load( fir_m , dz/sigmaz , fir_wt ) ;
457 #if 0
458 fac = fir_wt[0] = 1.0f ;
459 for( ii=1 ; ii <= fir_m ; ii++ ){
460 fir_wt[ii] = exp(-0.5*(ii*dz)*(ii*dz)/(sigmaz*sigmaz)) ;
461 fac += 2.0f * fir_wt[ii] ;
462 }
463 fac = 1.0f / fac ;
464 for( ii=0 ; ii <= fir_m ; ii++ ) fir_wt[ii] *= fac ;
465 #endif
466 fir_blurz( fir_m , fir_wt , nx,ny,nz , ffim ) ;
467 fir_num++ ; goto ALL_DONE_NOW ;
468 }
469
470 STATUS("start z FFTs") ;
471
472 aa = sigmaz * sigmaz * 0.5 ;
473 nup = nz + (int)(3.0 * sigmaz / dz) ; /* min FFT length */
474 nup = csfft_nextup_one35(nup) ; nby2 = nup / 2 ;
475
476 if(DO_INFO) INFO_message("z FFT(%d)",nup) ;
477
478 if( nup > ncg ){
479 cx = (complex *)realloc(cx,sizeof(complex)*nup) ;
480 gg = (float *)realloc(gg,sizeof(float )*nup) ; ncg = nup ;
481 }
482
483 dk = (2.0*PI) / (nup * dz) ;
484 fac = 1.0 / nup ;
485 gg[0] = fac ;
486 for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
487
488 switch( dtype ){
489 case MRI_short:{
490 register short *qfim ;
491 for( kk=0 ; kk < ny ; kk++ ){
492 for( jj=0 ; jj < nx ; jj+=2 ){
493 base = jj + kk*nx ;
494 qfim = sfim + base ;
495 if( jj == nx-1 )
496 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
497 else
498 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
499 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
500 csfft_cox( -1 , nup , cx ) ;
501 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
502 csfft_cox( 1 , nup , cx ) ;
503 if( jj == nx-1 )
504 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
505 else
506 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
507 }
508 }
509 }
510 break ;
511
512 case MRI_float:{
513 register float *qfim ;
514 for( kk=0 ; kk < ny ; kk++ ){
515 for( jj=0 ; jj < nx ; jj+=2 ){
516 base = jj + kk*nx ;
517 qfim = ffim + base ;
518 if( jj == nx-1 )
519 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
520 else
521 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
522 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
523 csfft_cox( -1 , nup , cx ) ;
524 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
525 csfft_cox( 1 , nup , cx ) ;
526 if( jj == nx-1 )
527 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
528 else
529 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
530 }
531 }
532 }
533 break ;
534
535 case MRI_byte:{
536 register byte *qfim ;
537 for( kk=0 ; kk < ny ; kk++ ){
538 for( jj=0 ; jj < nx ; jj+=2 ){
539 base = jj + kk*nx ;
540 qfim = bfim + base ;
541 if( jj == nx-1 )
542 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
543 else
544 for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
545 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
546 csfft_cox( -1 , nup , cx ) ;
547 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
548 csfft_cox( 1 , nup , cx ) ;
549 if( jj == nx-1 )
550 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
551 else
552 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
553 }
554 }
555 }
556 break ;
557
558 case MRI_complex:{
559 register complex *qfim ;
560 for( kk=0 ; kk < ny ; kk++ ){
561 for( jj=0 ; jj < nx ; jj++ ){
562 base = jj + kk*nx ;
563 qfim = cfim + base ;
564 for( ii=0 ; ii<nz ; ii++){ cx[ii] = qfim[ii*nxy] ; }
565 for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
566 csfft_cox( -1 , nup , cx ) ;
567 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
568 csfft_cox( 1 , nup , cx ) ;
569 for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii] ; }
570 }
571 }
572 }
573 break ;
574 }
575
576 /*** 10 Jan 2003: clip data to bot and top found above ***/
577 /*** to minimize Gibbs ringing artifacts ***/
578
579 ALL_DONE_NOW:
580
581 if( !all_fir && (fir_num < 3 || as_float) ){
582 STATUS("clipping results") ;
583 switch( ftype ){
584
585 case MRI_short:
586 for( ii=0 ; ii < nxyz ; ii++ ) {
587 if( as_float ) { /* return floats to rounded shorts */
588 if( ffim[ii] >= 0.0f ) sfim[ii] = (short)(ffim[ii]+0.5) ;
589 else sfim[ii] = (short)(ffim[ii]-0.5) ;
590 }
591 if( sfim[ii] < fbot ) sfim[ii] = fbot ;
592 else if( sfim[ii] > ftop ) sfim[ii] = ftop ;
593 }
594 break ;
595
596 case MRI_float:
597 for( ii=0 ; ii < nxyz ; ii++ )
598 if( ffim[ii] < fbot ) ffim[ii] = fbot ;
599 else if( ffim[ii] > ftop ) ffim[ii] = ftop ;
600 break ;
601
602 case MRI_byte:
603 for( ii=0 ; ii < nxyz ; ii++ ) {
604 if( as_float ) bfim[ii] = (byte)(ffim[ii]+0.5) ; /* return to byte */
605 if( bfim[ii] < fbot ) bfim[ii] = fbot ;
606 else if( bfim[ii] > ftop ) bfim[ii] = ftop ;
607 }
608 break ;
609 }
610 }
611
612 /*** done! ***/
613
614 if( cx != NULL ) free(cx) ;
615 if( gg != NULL ) free(gg) ;
616 if( as_float && ffim ) free(ffim) ;
617 EXRETURN ;
618 }
619
620 #if 0
621 /*-------------------------------------------------------------------*/
622 /*! Partition of unity function:
623
624 n=+infinity
625 sum partu(x-n) = 1 identically
626 n=-infinity
627
628 and support of partu(x) is interval [-1,1]
629 *//*-----------------------------------------------------------------*/
630
631 static INLINE float partu( float x )
632 {
633 register float ax ;
634 ax = fabsf(x) ;
635 if( ax >= 1.0f ) return 0.0f ;
636 return (x*x*(20.0f*ax+10.0f)+4.0f*ax+1.0f) ;
637 }
638 #endif
639
640 /*-------------------------------------------------------------------*/
641 /*! Function to blur a 3D volume in-place with a symmetric FIR filter
642 along the x-direction.
643 - m = stencil size (+/- m points in each direction; m >= 1)
644 - wt = array of weights
645 - nx,ny,nz = dimensions of 3D array
646 - f = 3D array
647
648 f_out[i,j,k] = wt[0] * f_in[i,j,k]
649 + wt[1] *(f_in[i+1,j,k]+f_in[i-1,j,k))
650 + wt[2] *(f_in[i+2,j,k]+f_in[i-2,j,k))
651 + ...
652 + wt[m] *(f_in[i+m,j,k]+f_in[i-m,j,k))
653
654 Similar routines for blurring along the y- and z-directions
655 are fir_blury() and fir_blurz().
656
657 -- RWCox - 03 Oct 2005 - trying for some speedup for Daniel Glen
658 ---------------------------------------------------------------------*/
659
fir_blurx(int m,float * wt,int nx,int ny,int nz,float * f)660 void fir_blurx( int m, float *wt,int nx, int ny, int nz, float *f )
661 {
662 int ii,jj,kk,qq , nxy=nx*ny , off ;
663 float *r , wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ;
664
665 ENTRY("fir_blurx") ;
666 if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);}
667
668 if( m < 1 || wt == NULL || nx < (m+1) || f == NULL ) EXRETURN ;
669 if( ny <= 0 || nz <= 0 ) EXRETURN ;
670 switch(m){ /** assign weights to variables not arrays **/
671 case 7:
672 wt7 = wt[7]; /* let cases fall through to next case to assign weights */
673 case 6:
674 wt6 = wt[6];
675 case 5:
676 wt5 = wt[5];
677 case 4:
678 wt4 = wt[4];
679 case 3:
680 wt3 = wt[3];
681 case 2:
682 wt2 = wt[2];
683 case 1:
684 wt1 = wt[1];
685 case 0:
686 wt0 = wt[0];
687 default:
688 break ;
689 }
690
691 /* 1 row workspace, with m-long buffers at each end
692 (so that the i-th element of the row is in r[i+m]) */
693
694 r = (float *)calloc(sizeof(float),(nx+2*m)) ;
695
696 switch( m ){
697
698 default:
699 for( kk=0 ; kk < nz ; kk++ ){
700 for( jj=0 ; jj < ny ; jj++ ){
701 off = jj*nx + kk*nxy ; ff = f+off ; /* ff = ptr to this row */
702 memcpy( r+m , ff , sizeof(float)*nx ) ; /* copy row into workspace */
703 r[m-1] = r[m+1] ; r[nx+m] = r[nx+m-2] ; /* mirror at ends */
704
705 for( ii=0 ; ii < nx ; ii++ ){ /* filter at ii-th location */
706 sum = wt[0]*r[ii+m] ;
707 for( qq=1 ; qq <= m ; qq++ )
708 sum += wt[qq] * ( r[ii+m-qq] + r[ii+m+qq] ) ;
709 ff[ii] = sum ; /* save result back in input array */
710 }
711 }}
712 break ;
713
714 /** for the cases below, the innermost loop
715 (qq, above) is completely unrolled for speedup **/
716
717 #undef M
718 #define M 7
719 case 7:
720 for( kk=0 ; kk < nz ; kk++ ){
721 for( jj=0 ; jj < ny ; jj++ ){
722 off = jj*nx + kk*nxy ; ff = f+off ;
723 memcpy( r+m , ff , sizeof(float)*nx ) ;
724 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
725
726 for( ii=0 ; ii < nx ; ii++ )
727 ff[ii] = wt7*(r[ii ]+r[ii+14])
728 +wt6*(r[ii+1]+r[ii+13])
729 +wt5*(r[ii+2]+r[ii+12])
730 +wt4*(r[ii+3]+r[ii+11])
731 +wt3*(r[ii+4]+r[ii+10])
732 +wt2*(r[ii+5]+r[ii+ 9])
733 +wt1*(r[ii+6]+r[ii+ 8])+wt0*r[ii+7] ;
734 }}
735 break ;
736
737 #undef M
738 #define M 6
739 case 6:
740 for( kk=0 ; kk < nz ; kk++ ){
741 for( jj=0 ; jj < ny ; jj++ ){
742 off = jj*nx + kk*nxy ; ff = f+off ;
743 memcpy( r+m , ff , sizeof(float)*nx ) ;
744 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
745
746 for( ii=0 ; ii < nx ; ii++ )
747 ff[ii] = wt6*(r[ii ]+r[ii+12])
748 +wt5*(r[ii+1]+r[ii+11])
749 +wt4*(r[ii+2]+r[ii+10])
750 +wt3*(r[ii+3]+r[ii+ 9])
751 +wt2*(r[ii+4]+r[ii+ 8])
752 +wt1*(r[ii+5]+r[ii+ 7])+wt0*r[ii+6] ;
753 }}
754 break ;
755
756 #undef M
757 #define M 5
758 case 5:
759 for( kk=0 ; kk < nz ; kk++ ){
760 for( jj=0 ; jj < ny ; jj++ ){
761
762 off = jj*nx + kk*nxy ; ff = f+off ;
763 memcpy( r+m , ff , sizeof(float)*nx ) ;
764 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
765
766 for( ii=0 ; ii < nx ; ii++ )
767 ff[ii] = wt5*(r[ii ]+r[ii+10])
768 +wt4*(r[ii+1]+r[ii+ 9])
769 +wt3*(r[ii+2]+r[ii+ 8])
770 +wt2*(r[ii+3]+r[ii+ 7])
771 +wt1*(r[ii+4]+r[ii+ 6])+wt0*r[ii+5] ;
772 }}
773 break ;
774
775 #undef M
776 #define M 4
777 case 4:
778 for( kk=0 ; kk < nz ; kk++ ){
779 for( jj=0 ; jj < ny ; jj++ ){
780
781 off = jj*nx + kk*nxy ; ff = f+off ;
782 memcpy( r+m , ff , sizeof(float)*nx ) ;
783 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
784
785 for( ii=0 ; ii < nx ; ii++ )
786 ff[ii] = wt4*(r[ii ]+r[ii+ 8])
787 +wt3*(r[ii+1]+r[ii+ 7])
788 +wt2*(r[ii+2]+r[ii+ 6])
789 +wt1*(r[ii+3]+r[ii+ 5])+wt0*r[ii+4] ;
790 }}
791 break ;
792
793 #undef M
794 #define M 3
795 case 3:
796 for( kk=0 ; kk < nz ; kk++ ){
797 for( jj=0 ; jj < ny ; jj++ ){
798
799 off = jj*nx + kk*nxy ; ff = f+off ;
800 memcpy( r+m , ff , sizeof(float)*nx ) ;
801 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
802
803 for( ii=0 ; ii < nx ; ii++ )
804 ff[ii] = wt3*(r[ii ]+r[ii+ 6])
805 +wt2*(r[ii+1]+r[ii+ 5])
806 +wt1*(r[ii+2]+r[ii+ 4])+wt0*r[ii+3] ;
807 }}
808 break ;
809
810 #undef M
811 #define M 2
812 case 2:
813 for( kk=0 ; kk < nz ; kk++ ){
814 for( jj=0 ; jj < ny ; jj++ ){
815
816 off = jj*nx + kk*nxy ; ff = f+off ;
817 memcpy( r+m , ff , sizeof(float)*nx ) ;
818 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
819
820 for( ii=0 ; ii < nx ; ii++ )
821 ff[ii] = wt2*(r[ii ]+r[ii+ 4])
822 +wt1*(r[ii+1]+r[ii+ 3])+wt0*r[ii+2] ;
823 }}
824 break ;
825
826 #undef M
827 #define M 1
828 case 1:
829 for( kk=0 ; kk < nz ; kk++ ){
830 for( jj=0 ; jj < ny ; jj++ ){
831
832 off = jj*nx + kk*nxy ; ff = f+off ;
833 memcpy( r+m , ff , sizeof(float)*nx ) ;
834 r[M-1] = r[M+1] ; r[nx+M] = r[nx+M-2] ;
835
836 for( ii=0 ; ii < nx ; ii++ )
837 ff[ii] = wt1*(r[ii]+r[ii+2])+wt0*r[ii+1] ;
838 }}
839 break ;
840
841 } /* end of switch on m */
842
843 free((void *)r) ; EXRETURN ;
844 }
845
846 /*-------------------------------------------------------------------*/
847 #undef D
848 #define D nx /* stride along y-axis */
849
850 /*! Similar to fir_blurx(), but along the y-axis.
851 For further comments, see fir_blurx() source code. */
852
fir_blury(int m,float * wt,int nx,int ny,int nz,float * f)853 void fir_blury( int m, float *wt,int nx, int ny, int nz, float *f )
854 {
855 int ii,jj,kk,qq , nxy=nx*ny , off ;
856 float *r, wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ;
857 float *rr, *ss;
858 int ny2m = ny+2*m;
859
860 ENTRY("fir_blury") ;
861 if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);}
862
863 if( m < 1 || wt == NULL || ny < (m+1) || f == NULL ) EXRETURN ;
864 if( nx <= 0 || nz <= 0 ) EXRETURN ;
865 switch(m){ /**assign weights to variables not arrays **/
866 case 7:
867 wt7 = wt[7]; /* let cases fall through to next case to assign weights */
868 case 6:
869 wt6 = wt[6];
870 case 5:
871 wt5 = wt[5];
872 case 4:
873 wt4 = wt[4];
874 case 3:
875 wt3 = wt[3];
876 case 2:
877 wt2 = wt[2];
878 case 1:
879 wt1 = wt[1];
880 case 0:
881 wt0 = wt[0];
882 default:
883 break ;
884 }
885
886 if( nx < 512 ) goto SMALLIMAGE;
887
888 /* In this function, for each value of kk (z index), we extract a
889 2D (y,x) slice, with m-long buffers on each side in the y-direction.
890 The purpose of this is to get multiple lines of y-direction data into
891 the CPU cache, to speed up processing (a lot). For the x-axis, this
892 was unneeded, since the x-rows are contiguous in memory. For data at 256x256
893 this 2D extract/process/insert trick was nugatory. However, for 512x512 data
894 this trick becomes important.
895 The same method is used for the z-axis in fir_blurz(). */
896
897 /* macro to access the input data 2D slice: (i,j) = (x,y) indexes */
898
899 #undef RR
900 #define RR(i,j) rr[(j)+m+(i)*ny2m] /*** 0 <= i <= nx-1 ; -m <= m <= ny-1+m ***/
901
902 /* macro to access the output data 2D slice */
903
904 #undef SS
905 #define SS(i,k) ss[(k)+(i)*ny]
906
907 rr = (float *)calloc(sizeof(float),ny2m*nx) ; /* ny2m = ny+2*m */
908 ss = (float *)malloc(sizeof(float)*ny *nx) ;
909
910 for( kk=0 ; kk < nz ; kk++ ){ /* loop in z-direction (an xy slice at a time) */
911 off = kk*nxy ; ff = f+off ; /* ff = ptr to start of this 2D slice */
912
913 /* load data into 2D (y+2m,x) slice from 3D (x,y,z) array;
914 inner loop is over ii so as to access in the most contiguous way */
915
916 for( jj=0 ; jj < ny ; jj++ ){
917 for( ii=0 ; ii < nx ; ii++ ) RR(ii,jj) = ff[ii+D*jj] ; /* D = nx here */
918 }
919 for( ii=0 ; ii < nx ; ii++ ){
920 RR(ii,-1) = RR(ii,1) ; RR(ii,ny) = RR(ii,ny-2) ; /* edge reflection - */
921 /* only 1 point reflected*/
922 }
923
924 /* filter data in RR along y-direction, put into 2D SS array */
925
926 switch(m){ /** for small m, unroll the inner loop for speed **/
927
928 default:
929 for( ii=0 ; ii < nx ; ii++ ){
930 for( jj=0 ; jj < ny ; jj++ ){
931 sum = wt[0]*RR(ii,jj) ;
932 for( qq=1 ; qq <= m ; qq++ )
933 sum += wt[qq] * ( RR(ii,jj+qq) + RR(ii,jj-qq) ) ;
934 SS(ii,jj) = sum ;
935 }}
936 break ;
937
938 case 7:
939 for( ii=0 ; ii < nx ; ii++ ){
940 for( jj=0 ; jj < ny ; jj++ ){
941 SS(ii,jj) = wt7 * ( RR(ii,jj+7) + RR(ii,jj-7) )
942 + wt6 * ( RR(ii,jj+6) + RR(ii,jj-6) )
943 + wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) )
944 + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) )
945 + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) )
946 + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
947 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
948 + wt0 * RR(ii,jj) ;
949 }}
950 break ;
951
952 case 6:
953 for( ii=0 ; ii < nx ; ii++ ){
954 for( jj=0 ; jj < ny ; jj++ ){
955 SS(ii,jj) = wt6 * ( RR(ii,jj+6) + RR(ii,jj-6) )
956 + wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) )
957 + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) )
958 + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) )
959 + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
960 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
961 + wt0 * RR(ii,jj) ;
962 }}
963 break ;
964
965 case 5:
966 for( ii=0 ; ii < nx ; ii++ ){
967 for( jj=0 ; jj < ny ; jj++ ){
968 SS(ii,jj) = wt5 * ( RR(ii,jj+5) + RR(ii,jj-5) )
969 + wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) )
970 + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) )
971 + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
972 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
973 + wt0 * RR(ii,jj) ;
974 }}
975 break ;
976
977 case 4:
978 for( ii=0 ; ii < nx ; ii++ ){
979 for( jj=0 ; jj < ny ; jj++ ){
980 SS(ii,jj) = wt4 * ( RR(ii,jj+4) + RR(ii,jj-4) )
981 + wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) )
982 + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
983 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
984 + wt0 * RR(ii,jj) ;
985 }}
986 break ;
987
988 case 3:
989 for( ii=0 ; ii < nx ; ii++ ){
990 for( jj=0 ; jj < ny ; jj++ ){
991 SS(ii,jj) = wt3 * ( RR(ii,jj+3) + RR(ii,jj-3) )
992 + wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
993 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
994 + wt0 * RR(ii,jj) ;
995 }}
996 break ;
997
998 case 2:
999 for( ii=0 ; ii < nx ; ii++ ){
1000 for( jj=0 ; jj < ny ; jj++ ){
1001 SS(ii,jj) = wt2 * ( RR(ii,jj+2) + RR(ii,jj-2) )
1002 + wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
1003 + wt0 * RR(ii,jj) ;
1004 }}
1005 break ;
1006
1007 case 1:
1008 for( ii=0 ; ii < nx ; ii++ ){
1009 for( jj=0 ; jj < ny ; jj++ ){
1010 SS(ii,jj) = wt1 * ( RR(ii,jj+1) + RR(ii,jj-1) )
1011 + wt0 * RR(ii,jj) ;
1012 }}
1013 break ;
1014
1015 } /* end of special cases of m */
1016
1017 /* put SS array back into input 3D array;
1018 again, inner loop over ii for most contiguous access to f[] array */
1019
1020 for( jj=0 ; jj < ny ; jj++ ){
1021 for( ii=0 ; ii < nx ; ii++ ) ff[ii+D*jj] = SS(ii,jj) ;
1022 }
1023
1024 } /* end of loop over y-direction (zz) */
1025
1026 /*** finito, cara mia mine, oh, oh, oh, each time we part, my heart wants to die...***/
1027
1028 free((void *)ss) ; free((void *)rr) ; EXRETURN ;
1029
1030 /* for small images (nx<512), use slice as is and don't use reslicing trick*/
1031
1032 SMALLIMAGE:
1033
1034 r = (float *)calloc(sizeof(float),(ny+2*m)) ;
1035
1036 switch( m ){
1037
1038 default:
1039 for( kk=0 ; kk < nz ; kk++ ){
1040 for( ii=0 ; ii < nx ; ii++ ){
1041 off = ii + kk*nxy ; ff = f+off ;
1042 for( jj=0 ; jj < ny ; jj++ ) r[jj+m] = ff[D*jj] ;
1043 r[m-1] = r[m+1] ; r[ny+m] = r[ny+m-2] ;
1044
1045 for( jj=0 ; jj < ny ; jj++ ){
1046 sum = wt[0]*r[jj+m] ;
1047 for( qq=1 ; qq <= m ; qq++ )
1048 sum += wt[qq] * ( r[jj+m-qq] + r[jj+m+qq] ) ;
1049 ff[D*jj] = sum ;
1050 }
1051 }}
1052 break ;
1053
1054 #undef M
1055 #define M 7
1056 case 7:
1057 for( kk=0 ; kk < nz ; kk++ ){
1058 for( ii=0 ; ii < nx ; ii++ ){
1059 off = ii + kk*nxy ; ff = f+off ;
1060 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1061 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1062
1063 for( jj=0 ; jj < ny ; jj++ )
1064 ff[D*jj] = wt7*(r[jj ]+r[jj+14])
1065 +wt6*(r[jj+1]+r[jj+13])
1066 +wt5*(r[jj+2]+r[jj+12])
1067 +wt4*(r[jj+3]+r[jj+11])
1068 +wt3*(r[jj+4]+r[jj+10])
1069 +wt2*(r[jj+5]+r[jj+ 9])
1070 +wt1*(r[jj+6]+r[jj+ 8])+wt0*r[jj+7] ;
1071 }}
1072 break ;
1073
1074 #undef M
1075 #define M 6
1076 case 6:
1077 for( kk=0 ; kk < nz ; kk++ ){
1078 for( ii=0 ; ii < nx ; ii++ ){
1079 off = ii + kk*nxy ; ff = f+off ;
1080 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1081 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1082
1083 for( jj=0 ; jj < ny ; jj++ )
1084 ff[D*jj] = wt6*(r[jj ]+r[jj+12])
1085 +wt5*(r[jj+1]+r[jj+11])
1086 +wt4*(r[jj+2]+r[jj+10])
1087 +wt3*(r[jj+3]+r[jj+ 9])
1088 +wt2*(r[jj+4]+r[jj+ 8])
1089 +wt1*(r[jj+5]+r[jj+ 7])+wt0*r[jj+6] ;
1090 }}
1091 break ;
1092
1093 #undef M
1094 #define M 5
1095 case 5:
1096 for( kk=0 ; kk < nz ; kk++ ){
1097 for( ii=0 ; ii < nx ; ii++ ){
1098 off = ii + kk*nxy ; ff = f+off ;
1099 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1100 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1101
1102 for( jj=0 ; jj < ny ; jj++ )
1103 ff[D*jj] = wt5*(r[jj ]+r[jj+10])
1104 +wt4*(r[jj+1]+r[jj+ 9])
1105 +wt3*(r[jj+2]+r[jj+ 8])
1106 +wt2*(r[jj+3]+r[jj+ 7])
1107 +wt1*(r[jj+4]+r[jj+ 6])+wt0*r[jj+5] ;
1108 }}
1109 break ;
1110
1111 #undef M
1112 #define M 4
1113 case 4:
1114 for( kk=0 ; kk < nz ; kk++ ){
1115 for( ii=0 ; ii < nx ; ii++ ){
1116 off = ii + kk*nxy ; ff = f+off ;
1117 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1118 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1119
1120 for( jj=0 ; jj < ny ; jj++ )
1121 ff[D*jj] = wt4*(r[jj ]+r[jj+ 8])
1122 +wt3*(r[jj+1]+r[jj+ 7])
1123 +wt2*(r[jj+2]+r[jj+ 6])
1124 +wt1*(r[jj+3]+r[jj+ 5])+wt0*r[jj+4] ;
1125 }}
1126 break ;
1127
1128 #undef M
1129 #define M 3
1130 case 3:
1131 for( kk=0 ; kk < nz ; kk++ ){
1132 for( ii=0 ; ii < nx ; ii++ ){
1133 off = ii + kk*nxy ; ff = f+off ;
1134 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1135 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1136
1137 for( jj=0 ; jj < ny ; jj++ )
1138 ff[D*jj] = wt3*(r[jj ]+r[jj+ 6])
1139 +wt2*(r[jj+1]+r[jj+ 5])
1140 +wt1*(r[jj+2]+r[jj+ 4])+wt0*r[jj+3] ;
1141 }}
1142 break ;
1143
1144 #undef M
1145 #define M 2
1146 case 2:
1147 for( kk=0 ; kk < nz ; kk++ ){
1148 for( ii=0 ; ii < nx ; ii++ ){
1149
1150 off = ii + kk*nxy ; ff = f+off ;
1151 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1152 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1153
1154 for( jj=0 ; jj < ny ; jj++ )
1155 ff[D*jj] = wt2*(r[jj ]+r[jj+ 4])
1156 +wt1*(r[jj+1]+r[jj+ 3])+wt0*r[jj+2] ;
1157 }}
1158 break ;
1159
1160 #undef M
1161 #define M 1
1162 case 1:
1163 for( kk=0 ; kk < nz ; kk++ ){
1164 for( ii=0 ; ii < nx ; ii++ ){
1165 off = ii + kk*nxy ; ff = f+off ;
1166 for( jj=0 ; jj < ny ; jj++ ) r[jj+M] = ff[D*jj] ;
1167 r[M-1] = r[M+1] ; r[ny+M] = r[ny+M-2] ;
1168
1169 for( jj=0 ; jj < ny ; jj++ )
1170 ff[D*jj] = wt1*(r[jj]+r[jj+2])+wt0*r[jj+1] ;
1171 }}
1172 break ;
1173
1174 } /* end of switch on m */
1175
1176 free((void *)r) ; EXRETURN ;
1177 }
1178
1179 /*-------------------------------------------------------------------*/
1180 #undef D
1181 #define D nxy /* stride along z-axis */
1182
1183 /*! Similar to fir_blurx(), but along z-axis. */
1184
fir_blurz(int m,float * wt,int nx,int ny,int nz,float * f)1185 void fir_blurz( int m, float *wt,int nx, int ny, int nz, float *f )
1186 {
1187 int ii,jj,kk,qq , nxy=nx*ny , off ;
1188 float *rr,*ss , wt0=0.0,wt1=0.0,wt2=0.0,wt3=0.0,wt4=0.0,wt5=0.0,wt6=0.0,wt7=0.0 , sum , *ff ;
1189 int nz2m = nz+2*m ;
1190
1191 ENTRY("fir_blurz") ;
1192 if(PRINT_TRACING){char str[256];sprintf(str,"m=%d",m);STATUS(str);}
1193
1194 if( m < 1 || wt == NULL || nz < (m+1) || f == NULL ) EXRETURN ;
1195 if( nxy <= 0 ) EXRETURN ;
1196
1197 /* In this function, for each value of jj (y index), we extract a
1198 2D (z,x) slice, with m-long buffers on each side in the z-direction.
1199 The purpose of this is to get multiple lines of z-direction data into
1200 the CPU cache, to speed up processing (a lot). For the x-axis, this
1201 was unneeded, since the x-rows are contiguous in memory. For the
1202 y-axis, this trick might help, but only if a single (x,y) plane
1203 doesn't fit into cache. For nx=ny=256, 1 plane is 256 KB, so I
1204 decided that this 2D extract/process/insert trick was nugatory. */
1205
1206 switch(m){ /**assign weights to variables not arrays **/
1207 case 7:
1208 wt7 = wt[7]; /* let cases fall through to next case to assign weights */
1209 case 6:
1210 wt6 = wt[6];
1211 case 5:
1212 wt5 = wt[5];
1213 case 4:
1214 wt4 = wt[4];
1215 case 3:
1216 wt3 = wt[3];
1217 case 2:
1218 wt2 = wt[2];
1219 case 1:
1220 wt1 = wt[1];
1221 case 0:
1222 wt0 = wt[0];
1223 default:
1224 break ;
1225 }
1226
1227 /* macro to access the input data 2D slice: (i,k) = (x,z) indexes */
1228
1229 #undef RR
1230 #define RR(i,k) rr[(k)+m+(i)*nz2m] /*** 0 <= i <= nx-1 ; -m <= k <= nz-1+m ***/
1231
1232 /* macro to access the output data 2D slice */
1233
1234 #undef SS
1235 #define SS(i,k) ss[(k)+(i)*nz]
1236
1237 rr = (float *)calloc(sizeof(float),nz2m*nx) ; /* nz2m = nz+2*m */
1238 ss = (float *)malloc(sizeof(float)*nz *nx) ;
1239
1240 for( jj=0 ; jj < ny ; jj++ ){ /* loop in y-direction */
1241 off = jj*nx ; ff = f+off ; /* ff = ptr to start of this 2D slice */
1242
1243 /* load data into 2D (z,x) slice from 3D (x,y,z) array;
1244 inner loop is over ii so as to access in the most contiguous way */
1245
1246 for( kk=0 ; kk < nz ; kk++ ){
1247 for( ii=0 ; ii < nx ; ii++ ) RR(ii,kk) = ff[ii+D*kk] ;
1248 }
1249 for( ii=0 ; ii < nx ; ii++ ){
1250 RR(ii,-1) = RR(ii,1) ; RR(ii,nz) = RR(ii,nz-2) ; /* edge reflection */
1251 }
1252
1253 /* filter data in RR along z-direction, put into 2D SS array */
1254
1255 switch(m){ /** for small m, unroll the inner loop for speed **/
1256
1257 default:
1258 for( ii=0 ; ii < nx ; ii++ ){
1259 for( kk=0 ; kk < nz ; kk++ ){
1260 sum = wt[0]*RR(ii,kk) ;
1261 for( qq=1 ; qq <= m ; qq++ )
1262 sum += wt[qq] * ( RR(ii,kk+qq) + RR(ii,kk-qq) ) ;
1263 SS(ii,kk) = sum ;
1264 }}
1265 break ;
1266
1267 case 7:
1268 for( ii=0 ; ii < nx ; ii++ ){
1269 for( kk=0 ; kk < nz ; kk++ ){
1270 SS(ii,kk) = wt7 * ( RR(ii,kk+7) + RR(ii,kk-7) )
1271 + wt6 * ( RR(ii,kk+6) + RR(ii,kk-6) )
1272 + wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) )
1273 + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) )
1274 + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) )
1275 + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1276 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1277 + wt0 * RR(ii,kk) ;
1278 }}
1279 break ;
1280
1281 case 6:
1282 for( ii=0 ; ii < nx ; ii++ ){
1283 for( kk=0 ; kk < nz ; kk++ ){
1284 SS(ii,kk) = wt6 * ( RR(ii,kk+6) + RR(ii,kk-6) )
1285 + wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) )
1286 + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) )
1287 + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) )
1288 + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1289 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1290 + wt0 * RR(ii,kk) ;
1291 }}
1292 break ;
1293
1294 case 5:
1295 for( ii=0 ; ii < nx ; ii++ ){
1296 for( kk=0 ; kk < nz ; kk++ ){
1297 SS(ii,kk) = wt5 * ( RR(ii,kk+5) + RR(ii,kk-5) )
1298 + wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) )
1299 + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) )
1300 + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1301 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1302 + wt0 * RR(ii,kk) ;
1303 }}
1304 break ;
1305
1306 case 4:
1307 for( ii=0 ; ii < nx ; ii++ ){
1308 for( kk=0 ; kk < nz ; kk++ ){
1309 SS(ii,kk) = wt4 * ( RR(ii,kk+4) + RR(ii,kk-4) )
1310 + wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) )
1311 + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1312 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1313 + wt0 * RR(ii,kk) ;
1314 }}
1315 break ;
1316
1317 case 3:
1318 for( ii=0 ; ii < nx ; ii++ ){
1319 for( kk=0 ; kk < nz ; kk++ ){
1320 SS(ii,kk) = wt3 * ( RR(ii,kk+3) + RR(ii,kk-3) )
1321 + wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1322 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1323 + wt0 * RR(ii,kk) ;
1324 }}
1325 break ;
1326
1327 case 2:
1328 for( ii=0 ; ii < nx ; ii++ ){
1329 for( kk=0 ; kk < nz ; kk++ ){
1330 SS(ii,kk) = wt2 * ( RR(ii,kk+2) + RR(ii,kk-2) )
1331 + wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1332 + wt0 * RR(ii,kk) ;
1333 }}
1334 break ;
1335
1336 case 1:
1337 for( ii=0 ; ii < nx ; ii++ ){
1338 for( kk=0 ; kk < nz ; kk++ ){
1339 SS(ii,kk) = wt1 * ( RR(ii,kk+1) + RR(ii,kk-1) )
1340 + wt0 * RR(ii,kk) ;
1341 }}
1342 break ;
1343
1344 } /* end of special cases of m */
1345
1346 /* put SS array back into input 3D array;
1347 again, inner loop over ii for most contiguous access to f[] array */
1348
1349 for( kk=0 ; kk < nz ; kk++ ){
1350 for( ii=0 ; ii < nx ; ii++ ) ff[ii+D*kk] = SS(ii,kk) ;
1351 }
1352
1353 } /* end of loop over y-direction (jj) */
1354
1355 /*** finito, cara mia ***/
1356
1357 free((void *)ss) ; free((void *)rr) ; EXRETURN ;
1358 }
1359
1360 /*-------------------------------------------------------------------*/
1361 /*! Like EDIT_blur_volume(), using FIR only (no FFTs)
1362 and only for float arrays.
1363 ---------------------------------------------------------------------*/
1364
FIR_blur_volume(int nx,int ny,int nz,float dx,float dy,float dz,float * ffim,float sigma)1365 void FIR_blur_volume( int nx, int ny, int nz,
1366 float dx, float dy, float dz,
1367 float *ffim , float sigma )
1368 {
1369 if( ffim != NULL && sigma > 0.0 )
1370 FIR_blur_volume_3d(nx,ny,nz, dx,dy,dz, ffim, sigma,sigma,sigma) ;
1371 }
1372
1373 /*-------------------------------------------------------------------*/
1374 /* Gaussian blur in real space, of a float volume; like
1375 EDIT_blur_volume_3d(), but using FIR only.
1376 - nx,ny,nz = dimensions of array
1377 - dx,dy,dz = grid step sizes
1378 - ffim = array
1379 - sigmax = stdev for blur along x; if 0, no blurring in x; etc.
1380 ---------------------------------------------------------------------*/
1381
FIR_blur_volume_3d(int nx,int ny,int nz,float dx,float dy,float dz,float * ffim,float sigmax,float sigmay,float sigmaz)1382 void FIR_blur_volume_3d( int nx, int ny, int nz,
1383 float dx, float dy, float dz,
1384 float *ffim ,
1385 float sigmax, float sigmay, float sigmaz )
1386 {
1387 int fir_m , ii ;
1388 float *fir_wt=NULL , fac ;
1389
1390 double sfac = AFNI_numenv("AFNI_BLUR_FIRFAC") ;
1391 if( sfac < 2.0 ) sfac = 2.5 ;
1392
1393 ENTRY("FIR_blur_volume_3d") ;
1394
1395 if( ffim == NULL ) EXRETURN ;
1396 if( sigmax <= 0.0 && sigmay <= 0.0 && sigmaz <= 0.0 ) EXRETURN ;
1397
1398 if( dx <= 0.0 ) dx = 1.0 ;
1399 if( dy <= 0.0 ) dy = dx ;
1400 if( dz <= 0.0 ) dz = dx ;
1401
1402 /*-- blur along x --*/
1403
1404 if( sigmax > 0.0 && nx > 1 ){
1405 fir_m = (int) ceil( sfac * sigmax / dx ) ; /* about the 5% level */
1406 if( fir_m < 1 ) fir_m = 1 ;
1407 if( fir_m > nx/2 ) fir_m = nx/2 ;
1408 fir_wt = (float *)malloc(sizeof(float)*(fir_m+1)) ;
1409 fir_gaussian_load( fir_m , dx/sigmax , fir_wt ) ;
1410 fir_blurx( fir_m , fir_wt , nx,ny,nz , ffim ) ;
1411 }
1412
1413 /*-- blur along y --*/
1414
1415 if( sigmay > 0.0 && ny > 1 ){
1416 fir_m = (int) ceil( sfac * sigmay / dy ) ;
1417 if( fir_m < 1 ) fir_m = 1 ;
1418 if( fir_m > ny/2 ) fir_m = ny/2 ;
1419 fir_wt = (float *)realloc(fir_wt,sizeof(float)*(fir_m+1)) ;
1420 fir_gaussian_load( fir_m , dy/sigmay , fir_wt ) ;
1421 fir_blury( fir_m , fir_wt , nx,ny,nz , ffim ) ;
1422 }
1423
1424 /*-- blur along z --*/
1425
1426 if( sigmaz > 0.0 && nz > 1 ){
1427 fir_m = (int) ceil( sfac * sigmaz / dz ) ;
1428 if( fir_m < 1 ) fir_m = 1 ;
1429 if( fir_m > nz/2 ) fir_m = nz/2 ;
1430 fir_wt = (float *)realloc(fir_wt,sizeof(float)*(fir_m+1)) ;
1431 fir_gaussian_load( fir_m , dz/sigmaz , fir_wt ) ;
1432 fir_blurz( fir_m , fir_wt , nx,ny,nz , ffim ) ;
1433 }
1434
1435 if( fir_wt != NULL ) free(fir_wt) ;
1436 EXRETURN ;
1437 }
1438
1439 /*-----------------------------------------------------------------------------*/
1440 /*! Blurring, applied to a 2D RGB image. */
1441
mri_rgb_blur2D(float sig,MRI_IMAGE * im)1442 MRI_IMAGE * mri_rgb_blur2D( float sig , MRI_IMAGE *im )
1443 {
1444 MRI_IMARR *imar ;
1445 MRI_IMAGE *rim , *gim , *bim , *newim ;
1446
1447 ENTRY("mri_rgb_blur2D") ;
1448
1449 if( im == NULL || im->kind != MRI_rgb || sig <= 0.0f ) RETURN(NULL) ;
1450
1451 imar = mri_rgb_to_3float(im) ;
1452 rim = IMARR_SUBIM(imar,0) ;
1453 gim = IMARR_SUBIM(imar,1) ;
1454 bim = IMARR_SUBIM(imar,2) ;
1455
1456 FIR_blur_volume_3d( rim->nx,rim->ny,1 , 1.0f,1.0f,1.0f ,
1457 MRI_FLOAT_PTR(rim) , sig,sig,0.0f ) ;
1458 FIR_blur_volume_3d( gim->nx,gim->ny,1 , 1.0f,1.0f,1.0f ,
1459 MRI_FLOAT_PTR(gim) , sig,sig,0.0f ) ;
1460 FIR_blur_volume_3d( bim->nx,bim->ny,1 , 1.0f,1.0f,1.0f ,
1461 MRI_FLOAT_PTR(bim) , sig,sig,0.0f ) ;
1462
1463 newim = mri_3to_rgb(rim,gim,bim) ; MRI_COPY_AUX(newim,im) ;
1464 DESTROY_IMARR(imar) ; RETURN(newim) ;
1465 }
1466
1467 /*-----------------------------------------------------------------------------*/
1468 /*! Blurring, applied to a 2D byte image. */
1469
mri_byte_blur2D(float sig,MRI_IMAGE * im)1470 MRI_IMAGE * mri_byte_blur2D( float sig , MRI_IMAGE *im )
1471 {
1472 MRI_IMAGE *rim, *newim ;
1473
1474 ENTRY("mri_byte_blur2D") ;
1475 if( im == NULL || im->kind != MRI_byte || sig <= 0.0f ) RETURN(NULL) ;
1476
1477 rim = mri_to_mri( MRI_float, im ) ;
1478 FIR_blur_volume_3d( rim->nx,rim->ny,1 , 1.0f,1.0f,1.0f ,
1479 MRI_FLOAT_PTR(rim) , sig,sig,0.0f ) ;
1480
1481 newim = mri_to_mri( MRI_byte , rim );
1482 MRI_COPY_AUX(newim,im) ;
1483 mri_free(rim);
1484
1485 RETURN(newim) ;
1486 }
1487
1488 /*-----------------------------------------------------------------------------*/
1489 /*! Blurring, applied to a 2D float image. */
1490
mri_float_blur2D(float sig,MRI_IMAGE * im)1491 MRI_IMAGE * mri_float_blur2D( float sig , MRI_IMAGE *im )
1492 {
1493 MRI_IMAGE *newim ;
1494
1495 ENTRY("mri_float_blur2D") ;
1496
1497 if( im == NULL || im->kind != MRI_float || sig <= 0.0f ) RETURN(NULL) ;
1498 newim = mri_to_float(im) ;
1499 FIR_blur_volume_3d( newim->nx,newim->ny,1 , 1.0f,1.0f,1.0f ,
1500 MRI_FLOAT_PTR(newim) , sig,sig,0.0f ) ;
1501 RETURN(newim) ;
1502 }
1503
1504 /*-----------------------------------------------------------------------------*/
1505 /*! Blurring, applied to a 3D float image (sig is in voxels, not mm). */
1506
mri_float_blur3D(float sig,MRI_IMAGE * im)1507 MRI_IMAGE * mri_float_blur3D( float sig , MRI_IMAGE *im )
1508 {
1509 MRI_IMAGE *newim ;
1510
1511 ENTRY("mri_float_blur3D") ;
1512
1513 if( im == NULL || im->kind != MRI_float || sig <= 0.0f ) RETURN(NULL) ;
1514 newim = mri_to_float(im) ;
1515 FIR_blur_volume_3d( newim->nx,newim->ny,newim->nz , 1.0f,1.0f,1.0f ,
1516 MRI_FLOAT_PTR(newim) , sig,sig,sig ) ;
1517 RETURN(newim) ;
1518 }
1519
1520 /*-----------------------------------------------------------------------------*/
1521
1522 #undef NG
1523 #define NG 11
fir_gaussian_load(int m,float dx,float * ff)1524 static void fir_gaussian_load( int m, float dx, float *ff )
1525 {
1526 int ii ; float fac ; double xx ;
1527
1528 if( m < 0 || dx <= 0.0f || ff == NULL ) return ;
1529
1530 if(DO_INFO) INFO_message("fir_gaussian_load(%d,%g)",m,dx) ;
1531
1532 if( AFNI_yesenv("AFNI_BLUR_FIROLD") ){ /** the olde waye **/
1533 fac = ff[0] = 1.0f ;
1534 for( ii=1 ; ii <= m ; ii++ ){
1535 xx = ii*dx ;
1536 ff[ii] = exp(-0.5*xx*xx ) ; /* center of each cell */
1537 fac += 2.0f * ff[ii] ;
1538 }
1539 } else { /** the newe bettere waye **/
1540 double sum , dxj , ee ; int jj ;
1541 fac = 0.0f ; dxj = dx/(2.0*NG) ; /* subdivision of each cell */
1542 for( ii=0 ; ii <= m ; ii++ ){
1543 sum = 0.0 ;
1544 for( jj=-NG ; jj <= NG ; jj++ ){ /* sum over subdivisions */
1545 xx = ii*dx + jj*dxj ;
1546 ee = exp( -0.5*xx*xx ) ;
1547 if( jj==-NG || jj==NG ) sum += 0.5*ee ; else sum += ee ;
1548 }
1549 ff[ii] = sum ;
1550 if( ii==0 ) fac += sum ; else fac += 2.0*sum ;
1551 }
1552 }
1553
1554 fac = 1.0f / fac ;
1555 for( ii=0 ; ii <= m ; ii++ ) ff[ii] *= fac ;
1556 return ;
1557 }
1558
1559 /**************************************************************************/
1560 /**************************************************************************/
1561 /**************************************************************************/
1562 #if 0
1563 /*------------------------------------------------------------------------*/
1564 /* Below is a test program that can be used to evaluate the relative
1565 speed of FFT and FIR blurring. If compiled into 'tblur', then a
1566 command line would be
1567 'tblur 240 13.0'
1568 meaning blur a 240x240x240 volume with sigma=13.0. The output is
1569 a line indicating the input parameters and the CPU time ratio
1570 (FFT time)/(FIR time).
1571
1572 Compilation should be something like so:
1573
1574 make tblur.o
1575 cc -o tblur tblur.o -L. -L/usr/X11R6/lib -lmri -lXt -lm
1576 --------------------------------------------------------------------------*/
1577 #include "mrilib.h"
1578
1579 int main( int argc , char *argv[] )
1580 {
1581 int nx , nxxx , ii , jj ;
1582 float sigma , *gar , *far ;
1583 double c1,c2,c3 , cg,cf ;
1584
1585 EDIT_blur_allow_fir(0) ; /* turn off auto-FIR in EDIT_blur_volume() */
1586
1587 /* get command line params */
1588
1589 if( argc < 3 ) exit(0) ;
1590 nx = (int)strtod( argv[1] , NULL) ; if( nx < 32 ) nx = 128 ;
1591 sigma = (float)strtod( argv[2] , NULL ) ; if( sigma <= 0.0 ) exit(0) ;
1592
1593 /* allocate volumes */
1594
1595 nxxx = nx*nx*nx ;
1596 gar = (float *)malloc( sizeof(float)*nxxx ) ;
1597 far = (float *)malloc( sizeof(float)*nxxx ) ;
1598
1599 /* load volume for FFT and then blur 5 times */
1600
1601 for( ii=0 ; ii < nxxx ; ii++ ) gar[ii] = cos(0.37382*ii) ;
1602 c1 = COX_cpu_time() ;
1603 for( jj=0 ; jj < 5 ; jj++ )
1604 EDIT_blur_volume( nx,nx,nx , 1.0,1.0,1.0 , MRI_float,gar , sigma ) ;
1605 c2 = COX_cpu_time() ; cg = c2-c1 ;
1606 /** printf("nx=%d sigma=%.2f Gaussian cpu=%.3f\n",nx,sigma,cg) ; **/
1607
1608 /* load volume for FIR and then blur 5 times */
1609
1610 for( ii=0 ; ii < nxxx ; ii++ ) far[ii] = cos(0.37382*ii) ;
1611 c1 = COX_cpu_time() ;
1612 for( jj=0 ; jj < 5 ; jj++ )
1613 FIR_blur_volume( nx,nx,nx , 1.0,1.0,1.0 , far , sigma ) ;
1614 c2 = COX_cpu_time() ; cf = c2-c1 ;
1615 /** printf("nx=%d sigma=%.2f FIR_blur cpu=%.3f\n",nx,sigma,cf) ; **/
1616
1617 /* output CPU time ratio */
1618
1619 printf("ratio: %d %.2f %.2f\n",nx,sigma,cg/cf) ;
1620
1621 /* compute mean difference between the 2 approaches */
1622
1623 /**
1624 c1 = 0.0 ;
1625 for( ii=0 ; ii < nxxx ; ii++ ) c1 += fabs(far[ii]-gar[ii]) ;
1626 c1 = c1 / nxxx ;
1627 printf("mean abs diff = %g\n",c1) ;
1628 **/
1629
1630 exit(0) ;
1631 }
1632 #endif
1633