1 #include "mrilib.h"
2
3 /*** This file is intended to be #include-d into another source file ***/
4
5 #include <time.h>
6 #include <sys/types.h>
7 #include <unistd.h>
8
9 #ifdef USE_OMP
10 # include <omp.h>
11 #endif
12
13 /*** Thread-safe FFT functions in csfft_OMP.c are used ***/
14
15 #include "csfft_OMP.c"
16 #define CFFT csfft_cox_OMP
17 #define CFFT_setup csfft_cox_OMP_SETUP()
18
19 /*----------------------------------------------------------------------------*/
20 /* What it says. The dimensions are assumed to be compatible with csfft_cox().
21 The array tar, if not NULL, must be at least as long as MAX(ny,nz).
22 *//*--------------------------------------------------------------------------*/
23
mri_fft3D_inplace(MRI_IMAGE * fim,complex * tar)24 static void mri_fft3D_inplace( MRI_IMAGE *fim , complex *tar )
25 {
26 complex *far , *qar ;
27 int nx,ny,nz,nxy, ii,jj,kk , nn,qq ;
28
29 if( fim == NULL || fim->kind != MRI_complex ) return ;
30 far = MRI_COMPLEX_PTR(fim) ;
31 if( far == NULL ) return ;
32
33 nx = fim->nx ; ny = fim->ny ; nz = fim->nz ; nxy = nx*ny ;
34
35 /* x FFTs */
36
37 CFFT_setup ; /* for the OMP version of csfft,
38 which allows csfft_cox_OMP to be called inside a thread */
39
40 for( qar=far,kk=0 ; kk < nz ; kk++ ){
41 for( jj=0 ; jj < ny ; jj++ , qar+=nx ){
42 CFFT( -1 , nx , qar ) ;
43 }}
44
45 /* y FFTs */
46
47 if( tar == NULL ){ /* create temp space */
48 qar = (complex *)malloc(sizeof(complex)*MAX(ny,nz)) ;
49 } else {
50 qar = tar ; /* user-supplied temp space */
51 }
52
53 for( kk=0 ; kk < nz ; kk++ ){
54 for( ii=0 ; ii < nx ; ii++ ){
55 qq = ii+kk*nxy ;
56 for( jj=0 ; jj < ny ; jj++ ) qar[jj] = far[qq+jj*nx] ;
57 CFFT( -1 , ny , qar ) ;
58 for( jj=0 ; jj < ny ; jj++ ) far[qq+jj*nx] = qar[jj] ;
59 }}
60
61 /* z FFTs */
62
63 for( jj=0 ; jj < ny ; jj++ ){
64 for( ii=0 ; ii < nx ; ii++ ){
65 qq = ii+jj*nx ;
66 for( kk=0 ; kk < nz ; kk++ ) qar[kk] = far[qq+kk*nxy] ;
67 CFFT( -1 , nz , qar ) ;
68 for( kk=0 ; kk < nz ; kk++ ) far[qq+kk*nxy] = qar[kk] ;
69 }}
70
71 if( tar == NULL ) free(qar) ; /* did we create this? */
72 return ;
73 }
74
75 /*----------------------------------------------------------------------------*/
76 /* Given a field size and a kernel function radius,
77 determine the grid size for the FFTs to perform, including
78 some buffering for edge effects of the kernel function.
79 Output grid sizes are even and compatible with csfft_cox().
80 *//*--------------------------------------------------------------------------*/
81
get_convo_field_size(int nx,int ny,int nz,float dx,float dy,float dz,float rr)82 static int_triple get_convo_field_size( int nx, int ny, int nz,
83 float dx, float dy, float dz, float rr )
84 {
85 int_triple ijk = { 0,0,0 } ;
86 int qq , xx,yy,zz ;
87
88 if( rr <= 0.0f ) return ijk ;
89
90 xx = nx + 2*(int)ceilf(rr/dx) ; /* expand for this buffer */
91 yy = ny + 2*(int)ceilf(rr/dy) ;
92 zz = nz + 2*(int)ceilf(rr/dz) ;
93 if( xx < 16 ) xx = 16 ;
94 if( yy < 16 ) yy = 16 ;
95 if( zz < 16 ) zz = 16 ;
96 xx = csfft_nextup_one35(xx) ; /* expand for FFT allowable sizes */
97 yy = csfft_nextup_one35(yy) ; /* (even with at most one) */
98 zz = csfft_nextup_one35(zz) ; /* (factor of 3 and/or 5,) */
99 /* (for efficient FFT-ing) */
100
101 ijk.i = xx ; ijk.j = yy ; ijk.k = zz ; return ijk ;
102 }
103
104 /*----------------------------------------------------------------------------*/
105 /* Create the kernel smoothing function, then FFT it, and clip/shrink it
106 as far as reasonable. The basic method for convolving an input 3D
107 field is then
108 (1) copy it to a larger buffered volume, then FFT it
109 (2) multiply this field by the radial weight function from here
110 (3) FFT back to real space
111 (4) truncate back to desired 3D grid
112 Note that nx,ny,nz must be compatible with csfft_cox(), which
113 will be easiest if get_convo_field_size() is used to set them.
114 *//*--------------------------------------------------------------------------*/
115
116 /* macros for 3D array access */
117
118 #undef WAR
119 #define WAR(i,j,k) war[(i)+(j)*nx+(k)*nxy]
120 #undef QAR
121 #define QAR(i,j,k) qar[(i)+(j)*nxh+(k)*nxyh]
122 #undef FAR
123 #define FAR(i,j,k) far[(i)+(j)*itop+(k)*ijtop]
124
make_radial_convo_weight(int nx,int ny,int nz,float dx,float dy,float dz,float (* kerfun)(float),float rr)125 static MRI_IMAGE * make_radial_convo_weight( int nx , int ny , int nz ,
126 float dx , float dy , float dz ,
127 float (*kerfun)(float) , float rr )
128 {
129 MRI_IMAGE *wim , *fim , *qim ;
130 complex *war ; float *far , *qar ;
131 int nxy=nx*ny , ii,jj,kk , nxh=nx/2,nyh=ny/2,nzh=nz/2 , nxyh=nxh*nyh ;
132 int itop,jtop,ktop,ijtop ;
133 float xx,yy,zz , ftop ;
134
135 if( kerfun == NULL ) return NULL ;
136
137 /* if radius not provided, make it up */
138
139 if( rr <= 0.0f ){
140 float *kar ;
141 itop = MAX(nx,ny) ; itop = MAX(itop,nz) ;
142 xx = MAX(dx,dy) ; xx = MAX(xx,dz) ;
143 kar = (float *)malloc(sizeof(float)*itop) ; ftop = 0.0f ;
144 for( ii=0 ; ii < itop ; ii++ ){
145 kar[ii] = fabsf(kerfun(ii*xx)) ; if( ftop < kar[ii] ) ftop = kar[ii] ;
146 }
147 if( ftop == 0.0f ){ free(kar) ; return NULL ; }
148 ftop *= 0.0001f ;
149 for( ii=itop-1 ; ii > 0 && kar[ii] < ftop ; ii-- ) ; /* nada */
150 free(kar) ;
151 rr = (ii+1)*xx ;
152 }
153
154 /* setup kernel in real space */
155
156 wim = mri_new_vol( nx,ny,nz , MRI_complex ) ;
157 war = MRI_COMPLEX_PTR(wim) ;
158
159 /* fill the array with the kernel in real space */
160
161 for( kk=0 ; kk < nz ; kk++ ){
162 zz = (kk < nzh) ? kk*dz : (nz-kk)*dz ;
163 for( jj=0 ; jj < ny ; jj++ ){
164 yy = (jj < nyh) ? jj*dy : (ny-jj)*dy ;
165 for( ii=0 ; ii < nx ; ii++ ){
166 if( ii==nxh || jj==nyh || kk==nzh ){ /* Nyquist */
167 WAR(ii,jj,kk) = CMPLX(0.0f,0.0f) ;
168 } else {
169 xx = (ii < nxh) ? ii*dx : (nx-ii)*dx ;
170 rr = kerfun( sqrtf(xx*xx+yy*yy+zz*zz) ) ;
171 WAR(ii,jj,kk) = CMPLX(rr,0.0f) ;
172 }
173 }
174 }
175 }
176
177 /* go to FFT space */
178
179 mri_fft3D_inplace( wim , NULL ) ;
180
181 /* create clipped output */
182
183 qim = mri_new_vol( nxh,nyh,nzh , MRI_float ) ;
184 qar = MRI_FLOAT_PTR(qim) ;
185
186 ftop = 0.00001f * fabsf(war[0].r) ; /* largest value to keep */
187
188 /* half size in each direction since is symmetric */
189 /* note we compute sqrt(FFT(wim)) since we are creating the
190 weight for the noise
191 -- which is effectively squared when the ACF is computed/estimated */
192
193 for( kk=0 ; kk < nzh ; kk++ ){
194 for( jj=0 ; jj < nyh ; jj++ ){
195 for( ii=0 ; ii < nxh ; ii++ ){
196 rr = WAR(ii,jj,kk).r ; /* imag part should be very small */
197 if( rr < ftop ) rr = 0.0f ; else rr = sqrtf(rr) ; /* clip here */
198 QAR(ii,jj,kk) = rr ; /* load output */
199 }}}
200
201 mri_free(wim) ; /* done with this */
202
203 /* shrink it in space, if possible */
204
205 MRI_autobbox( qim , NULL,&itop , NULL,&jtop , NULL,&ktop ) ;
206
207 if( itop < nxh-1 ) itop++ ;
208 if( jtop < nyh-1 ) jtop++ ;
209 if( ktop < nzh-1 ) ktop++ ;
210 ijtop = itop*jtop ;
211 fim = mri_new_vol( itop , jtop , ktop , MRI_float ) ;
212 far = MRI_FLOAT_PTR(fim) ;
213 ININFO_message("Kernel image dimensions %d x %d x %d",itop,jtop,ktop) ;
214
215 for( kk=0 ; kk < ktop ; kk++ ){
216 for( jj=0 ; jj < jtop ; jj++ ){
217 for( ii=0 ; ii < itop ; ii++ ){
218 FAR(ii,jj,kk) = QAR(ii,jj,kk) ;
219 }}}
220
221 mri_free(qim) ; return fim ;
222 }
223
224 #undef WAR
225 #undef QAR
226 #undef FAR
227
228 /*----------------------------------------------------------------------------*/
229 /* Convolve with a pre-computed kernel via FFT.
230 nx,ny,nz = dimensions of output
231 wtim = weight image in FFT space
232 tar = FFT workspace - at least MAX(ny,nz) long (or NULL)
233 *//*--------------------------------------------------------------------------*/
234
235 /* more macros for 3D array access */
236
237 #undef CXAR
238 #define CXAR(i,j,k) cxar[(i)+(j)*nx+(k)*nxy]
239
240 #undef WTAR
241 #define WTAR(i,j,k) wtar[(i)+(j)*nxw+(k)*nxyw]
242
243 #undef CXRAN
244 #define CXRAN(i,j,k) \
245 ( zr=ww*zgaussian_sss(xran), zi=ww*zgaussian_sss(xran), CXAR(i,j,k)=CMPLX(zr,zi) )
246
mri_radial_convo(MRI_IMAGE * fim,MRI_IMAGE * wtim,complex * tar)247 static void mri_radial_convo( MRI_IMAGE *fim , MRI_IMAGE *wtim , complex *tar )
248 {
249 MRI_IMAGE *cxim ; complex *cxar ;
250 MRI_IMARR *outar ;
251 float *wtar ;
252 int nxy , nxh,nyh,nzh , nxw,nyw,nzw,nxyw , ii,jj,kk ;
253 float ww , zr, zi ;
254
255 cxim = mri_new_vol( nx,ny,nz , MRI_complex ) ;
256 cxar = MRI_COMPLEX_PTR(cxim) ;
257 wtar = MRI_FLOAT_PTR(wtim) ;
258 nxw = wtim->nx ; nyw = wtim->ny ; nzw = wtim->nz ; nxyw = nxw*nyw ;
259 nxy = nx*ny ;
260
261 /* fill random values */
262
263 for( kk=1 ; kk < nzw ; kk++ ){ /* interior grid */
264 for( jj=1 ; jj < nyw ; jj++ ){
265 for( ii=1 ; ii < nxw ; ii++ ){
266 ww = WTAR(ii,jj,kk) ;
267 if( ww != 0.0f ){ /* fill all 8 reflections */
268 CXRAN(ii, jj, kk) ; CXRAN(nx-ii, jj, kk) ;
269 CXRAN(ii,ny-jj, kk) ; CXRAN(nx-ii,ny-jj, kk) ;
270 CXRAN(ii, jj,nz-kk) ; CXRAN(nx-ii, jj,nz-kk) ;
271 CXRAN(ii,ny-jj,nz-kk) ; CXRAN(nx-ii,ny-jj,nz-kk) ;
272 }
273 }}}
274
275 for( kk=1 ; kk < nzw ; kk++ ){ /* along ii=0 face */
276 for( jj=1 ; jj < nyw ; jj++ ){
277 ww = WTAR(0,jj,kk) ;
278 if( ww != 0.0f ){ /* just 4 reflections */
279 CXRAN(0,jj, kk) ; CXRAN(0,ny-jj, kk) ;
280 CXRAN(0,jj,nz-kk) ; CXRAN(0,ny-jj,nz-kk) ;
281 }
282 }}
283
284 for( jj=1 ; jj < nyw ; jj++ ){ /* along kk=0 face */
285 for( ii=1 ; ii < nxw ; ii++ ){
286 ww = WTAR(ii,jj,0) ;
287 if( ww != 0.0f ){ /* just 4 reflections */
288 CXRAN(ii, jj,0) ; CXRAN(nx-ii, jj,0) ;
289 CXRAN(ii,ny-jj,0) ; CXRAN(nx-ii,ny-jj,0) ;
290 }
291 }}
292
293 for( kk=1 ; kk < nzw ; kk++ ){ /* along jj=0 face */
294 for( ii=1 ; ii < nxw ; ii++ ){
295 ww = WTAR(ii,0,kk) ;
296 if( ww != 0.0f ){ /* just 4 reflections */
297 CXRAN(ii,0, kk) ; CXRAN(nx-ii,0, kk) ;
298 CXRAN(ii,0,nz-kk) ; CXRAN(nx-ii,0,nz-kk) ;
299 }
300 }}
301
302 for( kk=1 ; kk < nzw ; kk++ ){ /* along ii=jj=0 edge */
303 ww = WTAR(0,0,kk) ;
304 if( ww != 0.0f ){ /* just 2 reflections */
305 CXRAN(0,0,kk) ; CXRAN(0,0,nz-kk) ;
306 }
307 }
308
309 for( jj=1 ; jj < nyw ; jj++ ){ /* along ii=kk=0 edge */
310 ww = WTAR(0,jj,0) ;
311 if( ww != 0.0f ){ /* just 2 reflections */
312 CXRAN(0,jj,0) ; CXRAN(0,ny-jj,0) ;
313 }
314 }
315
316 for( ii=1 ; ii < nxw ; ii++ ){ /* along jj=kk=0 edge */
317 ww = WTAR(ii,0,0) ;
318 if( ww != 0.0f ){ /* just 2 reflections */
319 CXRAN(ii,0,0) ; CXRAN(nx-ii,0,0) ;
320 }
321 }
322
323 CXAR(0,0,0) = CMPLX(0.0f,0.0f) ; /* the origin ii=jj=kk=0 */
324 /* set to zero to get zero mean output */
325 /* FFT to real space */
326
327 mri_fft3D_inplace( cxim , tar ) ;
328
329 /* create output and exit */
330
331 outar = mri_complex_to_pair( cxim ) ;
332 mri_free(cxim) ;
333
334 return outar ;
335 }
336
337 #if 0
338 /*----------------------------------------------------------------------------*/
339 /* Main program for testing purposes */
340
341 int main( int argc , char *argv[] )
342 {
343 float rr , vv , parm[3] ;
344 int_triple ijk ;
345 float dxyz=2.0f ; int nxyz=128 ; double ct,wt ;
346 MRI_IMAGE *wim ;
347 THD_3dim_dataset *dset=NULL ;
348 THD_ivec3 nxyz_vec ; THD_fvec3 fxyz_vec , oxyz_vec ;
349 unsigned short xran[3] ; unsigned int gseed ; int ithr=1,nbrik=2 ;
350 int nthr=1 , ith ;
351
352 parm[0] = 0.66f ;
353 parm[1] = 6.62f ;
354 parm[2] = 11.09f ;
355
356 if( argc > 1 ){
357 nbrik = (int)strtod(argv[1],NULL) ;
358 if( nbrik < 2 ) nbrik = 2 ;
359 }
360 if( nbrik%2 == 1 ) nbrik++ ;
361
362 ijk = get_random_field_size( nxyz , nxyz , nxyz ,
363 dxyz , dxyz , dxyz , parm ) ;
364
365 nxyz = ijk.i ;
366 wim = make_radial_weight( nxyz,nxyz,nxyz , dxyz,dxyz,dxyz , parm ) ;
367
368 gseed = ((unsigned int)time(NULL)) + 17*(unsigned int)getpid() ;
369 xran[2] = ( gseed & 0xffff) + (unsigned short)ithr ;
370 xran[1] = ((gseed >> 16) & 0xffff) - (unsigned short)ithr ;
371 xran[0] = 0x330e + (unsigned short)ithr ;
372
373 INFO_message("creating random dataset") ;
374
375 dset = EDIT_empty_copy(NULL) ;
376 nxyz_vec.ijk[0] = nxyz ; nxyz_vec.ijk[1] = nxyz ; nxyz_vec.ijk[2] = nxyz ;
377 fxyz_vec.xyz[0] = dxyz ; fxyz_vec.xyz[1] = dxyz ; fxyz_vec.xyz[2] = dxyz ;
378 oxyz_vec.xyz[0] = 0.0f ; oxyz_vec.xyz[1] = 0.0f ; oxyz_vec.xyz[2] = 0.0f ;
379
380 EDIT_dset_items( dset ,
381 ADN_datum_all , MRI_float ,
382 ADN_nxyz , nxyz_vec ,
383 ADN_xyzdel , fxyz_vec ,
384 ADN_xyzorg , oxyz_vec ,
385 ADN_prefix , "./randomFFT" ,
386 ADN_nvals , nbrik ,
387 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
388 ADN_type , HEAD_ANAT_TYPE ,
389 ADN_view_type , VIEW_ORIGINAL_TYPE ,
390 ADN_func_type , ANAT_EPI_TYPE ,
391 ADN_none ) ;
392
393 #ifdef USE_OMP
394 #pragma omp parallel
395 { if( omp_get_thread_num() == 0 ){
396 nthr = omp_get_num_threads() ;
397 ININFO_message("number of threads = %d",nthr) ;
398 }}
399 #endif
400
401 ct = COX_cpu_time() ; wt = COX_clock_time() ;
402
403 AFNI_OMP_START ;
404 #pragma omp parallel
405 { int ib ; MRI_IMARR *abar; MRI_IMAGE *aim,*bim; complex *tar;
406 tar = (complex *)malloc(sizeof(complex)*(nxyz+nxyz)) ;
407 #pragma omp for
408 for( ib=0 ; ib < nbrik ; ib+=2 ){
409 abar = make_radial_random_field( nxyz,nxyz,nxyz , wim , tar , xran ) ;
410 aim = IMARR_SUBIM(abar,0) ; bim = IMARR_SUBIM(abar,1) ;
411 #pragma omp critical
412 { EDIT_substitute_brick( dset , ib , MRI_float , MRI_FLOAT_PTR(aim) ) ;
413 EDIT_substitute_brick( dset , ib+1 , MRI_float , MRI_FLOAT_PTR(bim) ) ;
414 }
415 FREE_IMARR(abar) ; abar = NULL ;
416 }
417 free(tar) ;
418 }
419 AFNI_OMP_END ;
420
421 ct = COX_cpu_time()-ct ; wt = COX_clock_time()-wt ;
422 ININFO_message("CPU time = %.1f Clock time = %.1f Ratio = %.1f",ct,wt,ct/wt) ;
423
424 ININFO_message("writing dataset") ;
425 DSET_write(dset) ; WROTE_DSET(dset) ; DSET_delete(dset) ; dset = NULL ;
426
427 exit(0) ;
428 }
429 #endif
430