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 /*** NOT 7D SAFE ***/
10
11 /**************************************************************************/
12 /*
13 mode = -1 for forward transform
14 = +1 for inverse (including scaling)
15
16 taper = fraction of data to taper (0 to 1)
17 */
18
mri_fft_complex(int mode,float taper,MRI_IMAGE * im)19 void mri_fft_complex( int mode , float taper , MRI_IMAGE *im )
20 {
21 float *rbuf , *ibuf , *xtap , *ytap ;
22 complex *cxim ;
23 int ii , jj , npix , jbase , nx,ny ;
24
25 if( im->kind != MRI_complex ){
26 fprintf( stderr , "mri_fft_complex only works on complex images!\n" ) ;
27 MRI_FATAL_ERROR ;
28 }
29
30 if( ! MRI_IS_2D(im) ){
31 fprintf(stderr,"mri_fft_complex only works on 2D images!\n") ;
32 MRI_FATAL_ERROR ;
33 }
34
35 /*** set up buffers ***/
36
37 npix = im->nx * im->ny ; /* number of pixels */
38 rbuf = (float *)malloc( sizeof(float) * npix ) ; /* real and imag buffs */
39 ibuf = (float *)malloc( sizeof(float) * npix ) ;
40 cxim = mri_data_pointer( im ) ; /* easy acces to im */
41
42 for( ii=0 ; ii < npix ; ii++ ){
43 rbuf[ii] = cxim[ii].r ;
44 ibuf[ii] = cxim[ii].i ;
45 }
46
47 /*** taper buffers, if desired ***/
48
49 if( taper > 0.0 && taper <= 1.0 ){
50 nx = im->nx ;
51 ny = im->ny ;
52 xtap = mri_setup_taper( nx , taper ) ;
53
54 /***
55 printf( "taper" ) ;
56 for( ii=0 ; ii < nx ; ii++ ){
57 if( (ii%5) == 0 ) printf("\n") ;
58 printf( "%12.4e " , xtap[ii] ) ;
59 }
60 printf("\n") ;
61 ***/
62
63 if( nx == ny ) ytap = xtap ;
64 else ytap = mri_setup_taper( ny , taper ) ;
65
66 for( jj=0 ; jj < ny ; jj++ ){
67 jbase = jj * nx ;
68 for( ii=0 ; ii < nx ; ii++ ){
69 rbuf[ii] *= xtap[ii] * ytap[jj] ;
70 ibuf[ii] *= xtap[ii] * ytap[jj] ;
71 }
72 }
73 free( xtap ) ;
74 if( ytap != xtap ) free(ytap) ;
75 }
76
77 /*** FFT buffers and copy them back to original image ***/
78
79 cfft2d_cox( mode , im->nx , im->ny , rbuf,ibuf ) ;
80
81 for( ii=0 ; ii < npix ; ii++ ){
82 cxim[ii].r = rbuf[ii] ;
83 cxim[ii].i = ibuf[ii] ;
84 }
85
86 return ;
87 }
88
89 /***********************************************************************/
90
mri_setup_taper(int nx,float taper)91 float *mri_setup_taper( int nx , float taper )
92 {
93 register int ii ;
94 int ntap ;
95 float *tap ;
96 float phi ;
97
98 tap = (float *)malloc( sizeof(float) * nx ) ; /* make array */
99
100 for( ii=0 ; ii < nx ; ii++ ) tap[ii] = 1.0 ; /* default 1's */
101
102 ntap = (int) (nx * 0.5 * taper ) ; /* # pts on each end */
103
104 if( ntap == 0 ){ /* special case of no points: taper is tiny */
105 tap[0] = tap[nx-1] = 0.5 ;
106 return tap ;
107 }
108
109 phi = PI / ntap ;
110 for( ii=0 ; ii < ntap ; ii++ ){
111 tap[ii] = 0.54 - 0.46 * cos( ii*phi ) ;
112 tap[nx-1-ii] = tap[ii] ;
113 }
114
115 return tap ;
116 }
117
118 /*----------------------------------------------------------------------------*/
119 /* macro to alternate signs in workspace array */
120
121 #undef ALTERN
122 #define ALTERN(nn) \
123 do{ register int qq; \
124 for( qq=1; qq<(nn); qq+=2 ){ cbig[qq].r=-cbig[qq].r; cbig[qq].i=-cbig[qq].i; } \
125 } while(0)
126
127 /*******
128 Note that thd_ballcorr.c has OpenMP compatible versions of these
129 3D FFT functions, which could make your life better (or at least faster).
130 ********/
131 /*----------------------------------------------------------------------------*/
132 /* FFT lengths are in Lxx, Lyy, Lzz; however,
133 Lxx = 0 ==> no FFT in that direction (etc.).
134 *//*--------------------------------------------------------------------------*/
135
mri_fft_3D(int Sign,MRI_IMAGE * inim,int Lxx,int Lyy,int Lzz,int alt)136 MRI_IMAGE * mri_fft_3D( int Sign, MRI_IMAGE *inim,
137 int Lxx,int Lyy,int Lzz, int alt )
138 {
139 MRI_IMAGE *outim ;
140 int ii,jj,kk , nx,ny,nxy,nz , nbig , fx,fy,fz,fxy , joff,koff ;
141 complex *cbig , *car , *far ;
142
143 if( inim->kind != MRI_complex ) return NULL ;
144
145 /* input data and its dimensions */
146
147 car = MRI_COMPLEX_PTR(inim) ;
148 nx = inim->nx ; ny = inim->ny ; nz = inim->nz ; nxy = nx*ny ;
149
150 /* output dimensions and data */
151
152 fx = (Lxx == 0) ? nx : (Lxx > nx) ? csfft_nextup_one35(Lxx) : csfft_nextup_one35(nx);
153 fy = (Lyy == 0) ? ny : (Lyy > ny) ? csfft_nextup_one35(Lyy) : csfft_nextup_one35(ny);
154 fz = (Lzz == 0) ? nz : (Lzz > nz) ? csfft_nextup_one35(Lzz) : csfft_nextup_one35(nz);
155 fxy = fx*fy ;
156
157 outim = mri_new_vol( fx,fy,fz , MRI_complex ) ; /* zero filled */
158 far = MRI_COMPLEX_PTR(outim) ;
159
160 /* buffer space */
161
162 nbig = MAX(fx,fy) ; nbig = MAX(nbig,fz) ; nbig = 4*nbig + 512 ;
163 cbig = (complex *)malloc(sizeof(complex)*nbig) ;
164
165 /* copy input data into output image */
166
167 for( kk=0 ; kk < nz ; kk++ )
168 for( jj=0 ; jj < ny ; jj++ )
169 memcpy( far + jj*fx + kk*fxy, car + jj*nx + kk*nxy, sizeof(complex)*nx );
170
171 /* x-direction FFTs */
172
173 if( Lxx > 1 ){
174 for( kk=0 ; kk < fz ; kk++ ){
175 koff = kk*fxy ;
176 for( jj=0 ; jj < fy ; jj++ ){
177 joff = koff + jj*fx ;
178 for( ii=0 ; ii < fx ; ii++ ) cbig[ii] = far[ii+joff] ;
179 if( alt > 0 ) ALTERN(fx) ;
180 csfft_cox( Sign , fx , cbig ) ;
181 if( alt < 0 ) ALTERN(fx) ;
182 for( ii=0 ; ii < fx ; ii++ ) far[ii+joff] = cbig[ii] ;
183 }
184 }
185 }
186
187 /* y-direction FFTs */
188
189 if( Lyy > 1 ){
190 for( kk=0 ; kk < fz ; kk++ ){
191 koff = kk*fxy ;
192 for( ii=0 ; ii < fx ; ii++ ){
193 joff = koff + ii ;
194 for( jj=0 ; jj < fy ; jj++ ) cbig[jj] = far[jj*fx+joff] ; /* copy data */
195 if( alt > 0 ) ALTERN(fy) ;
196 csfft_cox( Sign , fy , cbig ) ; /* FFT in buffer */
197 if( alt < 0 ) ALTERN(fy) ;
198 for( jj=0 ; jj < fy ; jj++ ) far[jj*fx+joff] = cbig[jj] ; /* copy back */
199 }
200 }
201 }
202
203 /* z-direction FFTs */
204
205 if( Lzz > 1 ){
206 for( jj=0 ; jj < fy ; jj++ ){
207 joff = jj*fx ;
208 for( ii=0 ; ii < fx ; ii++ ){
209 koff = joff + ii ;
210 for( kk=0 ; kk < fz ; kk++ ) cbig[kk] = far[kk*fxy+koff] ;
211 if( alt > 0 ) ALTERN(fz) ;
212 csfft_cox( Sign , fz , cbig ) ;
213 if( alt < 0 ) ALTERN(fz) ;
214 for( kk=0 ; kk < fz ; kk++ ) far[kk*fxy+koff] = cbig[kk] ;
215 }
216 }
217 }
218
219 free(cbig) ; MRI_COPY_AUX(outim,inim) ; return outim ;
220 }
221
222 /*----------------------------------------------------------------------------*/
223 /* Convolve (via FFT) image aim with bim. [Sep 2020]
224 Note output image will be at least as big as than the sum of the two sizes.
225 *//*--------------------------------------------------------------------------*/
226
mri_fft_3Dconvolve(MRI_IMAGE * aim,MRI_IMAGE * bim)227 MRI_IMAGE * mri_fft_3Dconvolve( MRI_IMAGE *aim , MRI_IMAGE *bim )
228 {
229 MRI_IMAGE *outim=NULL ;
230 MRI_IMAGE *paim , *pbim , *faim , *fbim ;
231 int nxa,nya,nza , nxb,nyb,nzb , Lxx,Lyy,Lzz , Lxyz,ii ;
232 complex ac , bc , qc ;
233 complex *acar , *bcar ;
234 float linv ;
235
236 if( aim == NULL || bim == NULL ) return NULL ;
237
238 /* input dimensions */
239
240 nxa = aim->nx ; nya = aim->ny ; nza = aim->nz ;
241 nxb = bim->nx ; nyb = bim->ny ; nzb = bim->nz ;
242
243 /* FFT and output dimensions (sum, bumped up for FFT effiency) */
244
245 Lxx = (nxa > 1 && nxb > 1) ? csfft_nextup_one35(nxa+nxb) : 0 ;
246 Lyy = (nya > 1 && nyb > 1) ? csfft_nextup_one35(nya+nyb) : 0 ;
247 Lzz = (nza > 1 && nzb > 1) ? csfft_nextup_one35(nza+nzb) : 0 ;
248
249 /* at this time, we don't allow for convolving a 3D image with a 1D
250 or 2D image, for example, which is possible but more complicated */
251
252 if( Lxx == 0 || Lyy == 0 || Lzz == 0 ) return NULL ;
253
254 /* 1) convert A image to complex
255 2) zero pad it to fit the FFT size
256 3) FFT that
257 Then repeat these steps for the B image */
258
259 faim = mri_to_complex( aim ) ; /* 1) */
260 paim = mri_zeropad_3D( 0,Lxx-nxa , 0,Lyy-nya , 0,Lzz-nza , faim ) ; /* 2) */
261 mri_free(faim) ;
262 faim = mri_fft_3D( -1 , paim , Lxx,Lyy,Lzz , 0 ) ; /* 3) */
263 mri_free(paim) ;
264 acar = MRI_COMPLEX_PTR(faim) ;
265
266 fbim = mri_to_complex( bim ) ; /* 1) */
267 pbim = mri_zeropad_3D( 0,Lxx-nxb , 0,Lyy-nyb , 0,Lzz-nzb , fbim ) ; /* 2) */
268 mri_free(fbim) ;
269 fbim = mri_fft_3D( -1 , pbim , Lxx,Lyy,Lzz , 0 ) ; /* 3) */
270 mri_free(pbim) ;
271 bcar = MRI_COMPLEX_PTR(fbim) ;
272
273 /* multiply+scale FFTs, store back in faim/acar */
274
275 Lxyz = Lxx * Lyy * Lzz ;
276 linv = 10.f / (float)Lxyz ; /* scaling for inverse FFT */
277 for( ii=0 ; ii < Lxyz ; ii++ ){
278 ac = acar[ii] ;
279 bc = bcar[ii] ;
280 qc.r = (ac.r * bc.r - ac.i * bc.i) * linv ; /* complex */
281 qc.i = (ac.r * bc.i + ac.i * bc.r) * linv ; /* multiply */
282 acar[ii] = qc ;
283 }
284 mri_free(fbim) ;
285
286 /* inverse FFT back to 'real' space */
287
288 fbim = mri_fft_3D( +1 , faim , Lxx,Lyy,Lzz , 0 ) ;
289 mri_free(faim) ;
290
291 /* convert to float-valued image (from complex FFT) and return */
292
293 outim = mri_complex_to_real( fbim ) ;
294 mri_free(fbim) ;
295 return outim ;
296 }
297