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