1 #include "mrilib.h"
2
3 /*** For inclusion into 3dBallMatch.c ***/
4
5 #define MINRAD 7 /* min radius for ball mask (voxels) */
6 #define USE_PM /* use +1/-1 instead of +1/0 for masks */
7 #define USE_FMASK /* automask the correlation map */
8 #define USE_SQRTF /* weight voxel CM by sqrtf(correlation) */
9
10 /*---------------------------------------------------------------------------*/
11 /* Modified 12 Nov 2020 to use the fftnf_OMP routine (fftn_omp.c),
12 which is FFTs that have been made thread-safe for OpenMP usage.
13 Not that any individual FFT is parallelized, but you can call
14 the FFT function from multiple threads. Which is what
15 function mri_fft_3Dconvolve_OMP() does.
16 Unlike the original csfft.c function or the original fftn.c function,
17 which are not thread safe.
18 *//*-------------------------------------------------------------------------*/
19 #define USE_FFTN /* undef this to avoid using fftnf_OMP and OpenMP */
20 /*---------------------------------------------------------------------------*/
21 #ifdef USE_FFTN
22
23 /* internal prototypes */
24 static MRI_IMAGE * mri_fft_3Dconvolve_OMP( MRI_IMAGE *aim, MRI_IMAGE *bim ) ;
25 extern int fftn_nextup_one35(int) ;
26 extern int fftn_nextup_even (int) ;
27 # define NEXTUP(x) fftn_nextup_even(x)
28
29 # define CONVO(a,b) mri_fft_3Dconvolve_OMP( (a) , (b) )
30
31 #else /*** not using fftn -- use a library function ***/
32
33 # define CONVO(a,b) mri_fft_3Dconvolve( (a) , (b) )
34 # define NEXTUP(x) csfft_nextup_one35(x)
35
36 #endif
37 /*---------------------------------------------------------------------------*/
38
39 /*---------------------------------------------------------------------------*/
40 /* Make a ball mask, radius brad mm, grid spacings given by dx,dy,dz */
41 /*---------------------------------------------------------------------------*/
42
make_ball_mask_float(float brad,float dx,float dy,float dz)43 MRI_IMAGE * make_ball_mask_float( float brad, float dx, float dy, float dz )
44 {
45 MRI_IMAGE *bim ; float *bmask ;
46 float xq,yq,zq,rq ;
47 int ii,jj,kk, irad,jrad,krad, idim,jdim,kdim, ic,jc,kc, ijk ;
48
49 if( brad <= 0.0f ) return NULL ;
50 if( dx == 0.0f ) dx = 1.0f ; else dx = fabsf(dx) ;
51 if( dy == 0.0f ) dy = 1.0f ; else dy = fabsf(dy) ;
52 if( dz == 0.0f ) dz = 1.0f ; else dz = fabsf(dz) ;
53
54 /* radius of ball in voxels */
55
56 irad = (int)rintf(1.01f*brad/dx) ;
57 jrad = (int)rintf(1.01f*brad/dy) ;
58 krad = (int)rintf(1.01f*brad/dz) ;
59
60 if( irad < MINRAD || jrad < MINRAD || krad < MINRAD ) return NULL ;
61
62 /* size of the box that holds the ball */
63
64 idim = 2*irad+1 ; ic = irad ;
65 jdim = 2*jrad+1 ; jc = jrad ;
66 kdim = 2*krad+1 ; kc = krad ;
67 rq = brad*brad ;
68
69 /* make the box image for the ball mask */
70
71 bim = mri_new_vol( idim,jdim,kdim , MRI_float ) ;
72 bmask = MRI_FLOAT_PTR(bim) ;
73 bim->dx = dx ; bim->dy = dy ; bim->dz = dz ;
74
75 for( ijk=kk=0 ; kk < kdim ; kk++ ){
76 zq = (kk-kc)*dz ; zq = zq*zq ;
77 for( jj=0 ; jj < jdim ; jj++ ){
78 yq = (jj-jc)*dy ; yq = zq + yq*yq ;
79 for( ii=0 ; ii < idim ; ii++,ijk++ ){
80 xq = (ii-ic)*dx ; xq = yq + xq*xq ;
81 #ifndef USE_PM
82 bmask[ijk] = ( xq <= rq ) ? 1.0f : 0.0f ; /* in=1 out=0 */
83 #else
84 bmask[ijk] = ( xq <= rq ) ? 1.0f : -1.0f ; /* in=1 out=-1 */
85 #endif
86 }}}
87
88 return bim ;
89 }
90
91 /*---------------------------------------------------------------------------*/
92 /* Find the center of mass of the overlapation */
93 /*---------------------------------------------------------------------------*/
94
MRI_ball_mask_overlapation(MRI_IMAGE * aim,float brad)95 int_triple MRI_ball_mask_overlapation( MRI_IMAGE *aim , float brad )
96 {
97 MRI_IMAGE *bim ;
98 MRI_IMAGE *fim ; float *far ; byte *fmask=NULL ;
99 int irad,jrad,krad , ii,jj,kk,ijk , nx,ny,nz,nxy ;
100 int_triple ijkout = { -666 , -666 , -666 } ; /* error return */
101 float xc,yc,zc,fc,ff,ftop ;
102
103 if( aim == NULL || brad <= 0.0f ) return ijkout ; /* bad inputs */
104
105 /* make a binary mask with a ball of radius brad mm */
106
107 /* ININFO_message("making ball mask") ; */
108 bim = make_ball_mask_float( brad , aim->dx , aim->dy , aim->dz ) ;
109 if( bim == NULL ) return ijkout ;
110
111 /* dimensions */
112
113 irad = (bim->nx-1)/2 ; /* dimensions of mask in voxels */
114 jrad = (bim->ny-1)/2 ;
115 krad = (bim->nz-1)/2 ;
116
117 /* convolve ball mask with the input anatomical mask */
118
119 /* INFO_message("about to CONVO") ; */
120 fim = CONVO( aim , bim ) ; /* see top of file for definition of CONVO */
121
122 /* INFO_message("about to mri_free(bim)") ; */
123 mri_free(bim) ;
124 if( fim == NULL ) return ijkout ;
125
126 /* find automask of the convolved image? */
127 #ifdef USE_FMASK
128 THD_automask_set_cheapo(1) ;
129 fmask = mri_automask_image( fim ) ;
130 if( fmask == NULL ){ mri_free(fim) ; return ijkout ; } /* bad */
131 #endif
132
133 far = MRI_FLOAT_PTR(fim) ;
134 nx = fim->nx; ny = fim->ny; nz = fim->nz; nxy = nx*ny;
135 fc = xc = yc = zc = ftop = 0.0f ;
136
137 /* find the max convolution value */
138
139 for( kk=krad ; kk < nz ; kk++ ){
140 for( jj=jrad ; jj < ny ; jj++ ){
141 for( ii=irad ; ii < nx ; ii++ ){
142 ijk = ii + jj*nx + kk*nxy ;
143 #ifdef USE_FMASK
144 if( !fmask[ijk] ) far[ijk] = 0.0f ; /* use the automask */
145 #endif
146 if( far[ijk] > ftop ) ftop = far[ijk] ; /* find the max */
147 }}}
148 #ifdef USE_FMASK
149 /* INFO_message("about to free(fmask)") ; */
150 free(fmask) ;
151 #endif
152 if( ftop == 0.0f ){
153 mri_free(fim) ; return ijkout ; /* bad */
154 }
155 ftop = 0.09f*ftop ; /* threshold for using voxel */
156
157 /* compute CM */
158
159 for( kk=krad ; kk < nz ; kk++ ){
160 for( jj=jrad ; jj < ny ; jj++ ){
161 for( ii=irad ; ii < nx ; ii++ ){
162 ijk = ii + jj*nx + kk*nxy ;
163 ff = far[ijk] ;
164 if( ff >= ftop ){
165 #ifdef USE_SQRTF
166 ff = sqrtf(ff); /* downweight */
167 #endif
168 xc += ff*ii; yc += ff*jj; zc += ff*kk; fc += ff;
169 }
170 }}}
171 /* INFO_message("about to mri_free(fim)") ; */
172 mri_free(fim) ;
173 if( fc == 0.0f ) return ijkout ; /* no positive values? */
174
175 /* subtract off ball radius from CM to get location in original dataset */
176
177 ijkout.i = (int)rintf(xc/fc) - irad ;
178 ijkout.j = (int)rintf(yc/fc) - jrad ;
179 ijkout.k = (int)rintf(zc/fc) - krad ;
180
181 return ijkout ;
182 }
183
184 /*---------------------------------------------------------------------------*/
185 /* Find the 'best' overlap with a ball of radius brad mm. */
186 /*---------------------------------------------------------------------------*/
187
THD_ball_mask_overlapation(THD_3dim_dataset * dset,float brad)188 int_triple THD_ball_mask_overlapation( THD_3dim_dataset *dset , float brad )
189 {
190 MRI_IMAGE *aim ; byte *aar ; int ii ;
191 int_triple ijkout = { -666 , -666 , -666 } ;
192 float dxyz ; int npeel ;
193
194 if( !ISVALID_DSET(dset) || brad <= 0.0f ) return ijkout ;
195
196 /* autmask the dataset */
197
198 dxyz = cbrtf(fabsf( DSET_DX(dset)*DSET_DY(dset)*DSET_DZ(dset) )) ;
199 if( dxyz <= 0.0f ) dxyz = 1.0f ;
200 npeel = (int)rintf( 0.222f*brad/dxyz ) ;
201 if( npeel <= 0 ) npeel = 1 ;
202 else if( npeel > 31 ) npeel = 31 ;
203 THD_automask_set_peelcounts(npeel,17) ;
204 THD_automask_set_clipfrac(0.135f) ;
205 THD_automask_extclip(1) ;
206 THD_automask_set_onlypos(1) ;
207 /* INFO_message("automasking dataset") ; */
208 THD_automask_verbose(0) ;
209 aar = THD_automask(dset) ; /* byte array with in=1 out=0 */
210 THD_automask_set_onlypos(0) ;
211
212 /* stuff automask into an MRI_IMAGE */
213
214 #ifndef USE_PM
215 /* use the 0-1 mask */
216 aim = mri_empty_conforming( DSET_BRICK(dset,0) , MRI_byte ) ;
217 mri_set_data_pointer( aim , aar ) ;
218 #else
219 /* convert to a -1/+1 mask */
220 { float *nar ; int nvox ;
221 aim = mri_new_conforming( DSET_BRICK(dset,0), MRI_float ) ;
222 nar = MRI_FLOAT_PTR(aim) ; nvox = aim->nvox ;
223 for( ii=0 ; ii < nvox ; ii++ ){
224 nar[ii] = ( aar[ii] ) ? +1.0f : -1.0f ; /* in=+1 out=-1 */
225 }
226 free(aar) ;
227 }
228 #endif
229 aim->dx = DSET_DX(dset); aim->dy = DSET_DY(dset); aim->dz = DSET_DZ(dset);
230
231 /* compute the 'best' mask overlap location and return it */
232
233 ijkout = MRI_ball_mask_overlapation( aim , brad ) ;
234 mri_free(aim) ;
235 return ijkout ;
236 }
237
238 /*===========================================================================*/
239 /*===========================================================================*/
240 /***** Below, new code for using a spheroid mask instead of a ball mask. ****/
241 /*===========================================================================*/
242 /*===========================================================================*/
243
244 /*---------------------------------------------------------------------------*/
245 /* Make a spheroid (ellipsoid of rotation) mask. Inputs are:
246 * ab = quintuple of floats
247 arad = radius along axis of rotation
248 ax,ay,az = 3-vector defining axis of rotation
249 brad = radius along 2 axis perpendicular to rotation axis
250 brad == 0 resets to brad = arad
251 * de = triple of floats
252 dx,dy,dz = grid spacings to use along each axis
253 (any 0 value is replaced by 1)
254 * spheroid is prolate if arad > brad ; oblate if arad < brad
255 *//*-------------------------------------------------------------------------*/
256
make_spheroid_mask_float(float_quint ab,float_triple de)257 MRI_IMAGE * make_spheroid_mask_float( float_quint ab , float_triple de )
258 {
259 MRI_IMAGE *bim ; float *bmask ;
260 float yq,zq , xqq,aqq,bqq,cqq,crad , xc,yc,zc , arq,brq , yac,zac ;
261 int ii,jj,kk, irad,jrad,krad, idim,jdim,kdim, ic,jc,kc, ijk ;
262 float arad, ax, ay, az ;
263 float brad, dx, dy, dz , dmin ;
264
265 /* unpack input structs into scalars */
266
267 arad = ab.a ; ax = ab.b ; ay = ab.c ; az = ab.d ; brad = ab.e ;
268 dx = de.a ; dy = de.b ; dz = de.c ;
269
270 if( dx == 0.0f ) dx = 1.0f ; else dx = fabsf(dx) ;
271 if( dy == 0.0f ) dy = 1.0f ; else dy = fabsf(dy) ;
272 if( dz == 0.0f ) dz = 1.0f ; else dz = fabsf(dz) ;
273 dmin = MIN(dx,dy) ; if( dz < dmin ) dmin = dz ;
274
275 if( arad <= 0.0f ) return NULL ; /* why do I have to deal with this crap? */
276
277 /* are we in the spherical case? */
278
279 if( brad == 0.0f || fabsf(arad-brad) <= dmin )
280 return make_ball_mask_float( arad , dx,dy,dz ) ;
281
282 /* OK, not a sphere/ball */
283
284 arq = 1.0f / (arad*arad) ; /* for convenience below, */
285 brq = 1.0f / (brad*brad) ; /* when computing in-or-out-ness of a location */
286
287 /* normalize (ax,ay,az) to unit length */
288
289 aqq = sqrtf(ax*ax + ay*ay + az*az) ; /* length of vector */
290 if( aqq == 0.0f ){
291 ax = 1.0f ; ay = az = 0.0f ; /* hopeless user */
292 } else if( aqq != 1.0f ){
293 ax = ax/aqq ; ay = ay/aqq ; az = az/aqq ; /* normalize */
294 }
295
296 /* 'radius' of mask in voxels [does not depend on orientation] */
297
298 crad = MAX(arad,brad) ; /* largest possible distance from center */
299 irad = (int)rintf(1.01f*crad/dx) ;
300 jrad = (int)rintf(1.01f*crad/dy) ;
301 krad = (int)rintf(1.01f*crad/dz) ;
302
303 if( irad < MINRAD || jrad < MINRAD || krad < MINRAD ) return NULL ;
304
305 /* size of the box that holds the ball */
306
307 idim = 2*irad+1 ; ic = irad ; /* dimensions are odd integers */
308 jdim = 2*jrad+1 ; jc = jrad ;
309 kdim = 2*krad+1 ; kc = krad ;
310
311 /* make the box image for the ball mask */
312
313 bim = mri_new_vol( idim,jdim,kdim , MRI_float ) ;
314 bmask = MRI_FLOAT_PTR(bim) ;
315 bim->dx = dx ; bim->dy = dy ; bim->dz = dz ;
316
317 /* loop over (xc,yc,zc) points in the grid;
318 compute aqq = length squared of (xc,yc,zc) projected along (ax,ay,az)
319 compute bqq = length squared of (xc,yc,zc) projected perp to (ax,ay,az)
320 = length squared of (xc,yc,zc) minus aqq (cf. Pythagoras)
321 then, a point is inside the spheroid if aqq/arad^2 + bqq/brad^2 <= 1 */
322
323 for( ijk=kk=0 ; kk < kdim ; kk++ ){
324 zc = (kk-kc)*dz ; zq = zc*zc ; zac = zc*az ;
325 for( jj=0 ; jj < jdim ; jj++ ){
326 yc = (jj-jc)*dy ; yq = yc*yc + zq ; yac = yc*ay + zac ;
327 for( ii=0 ; ii < idim ; ii++,ijk++ ){
328 xc = (ii-ic)*dx; xqq = xc*xc + yq; /* xqq = L^2 of (xc,yc,zc) */
329
330 aqq = xc*ax + yac ; aqq = aqq*aqq ; /* aqq = L^2 of (xc,yc,zc)
331 along (ax,ay,az) */
332 bqq = xqq - aqq ; /* bqq = L^2 of (xc,yc,zc)
333 perp to (ax,ay,az) */
334 cqq = aqq*arq + bqq*brq ; /* cqq = aqq/arad^2 + bqq/brad^2 */
335 #ifndef USE_PM
336 bmask[ijk] = ( cqq <= 1.0f ) ? 1.0f : 0.0f ; /* in=1 out=0 */
337 #else
338 bmask[ijk] = ( cqq <= 1.0f ) ? 1.0f : -1.0f ; /* in=1 out=-1 */
339 #endif
340 }}} /* end of 3D loop */
341
342 /* at this point, could use MRI_autobbox and mri_cut_3D to trim
343 image down to exclude any planes of all zeros around the edges;
344 HOWEVER, I want all boxes for a given spheroid size (arad,brad)
345 to be the same, for ease of combining results from them */
346
347 return bim ;
348 }
349
350 /*---------------------------------------------------------------------------*/
351 /* Find the center of mass of the overlapation */
352 /*---------------------------------------------------------------------------*/
353
354 typedef struct { int i,j,k ; float a ; } int3float1 ;
355
MRI_spheroid_overlapation1(MRI_IMAGE * aim,float_quint abin)356 int3float1 MRI_spheroid_overlapation1( MRI_IMAGE *aim , float_quint abin )
357 {
358 MRI_IMAGE *bim ;
359 MRI_IMAGE *fim ; float *far ; byte *fmask=NULL ;
360 int irad,jrad,krad , ii,jj,kk,ijk , nx,ny,nz,nvox ;
361 int3float1 ijkout = { -666 , -666 , -666 , 0.0f } ; /* error return */
362 float xc,yc,zc,fc,ff,ftop ;
363 float_quint ab ; float_triple de ;
364
365 if( aim == NULL ) return ijkout ; /* why me? */
366
367 de.a = aim->dx ; de.b = aim->dy ; de.c = aim->dz ;
368
369 /* make a binary mask for the desired spheroid */
370
371 bim = make_spheroid_mask_float( abin , de ) ;
372 if( bim == NULL ) return ijkout ;
373
374 /* dimensions */
375
376 irad = (bim->nx-1)/2 ; /* 'radii' of ball in voxels */
377 jrad = (bim->ny-1)/2 ;
378 krad = (bim->nz-1)/2 ;
379
380 /* convolve ball mask with the input anatomical mask */
381
382 /* INFO_message("about to CONVO") ; */
383 fim = CONVO( aim , bim ) ; /* see top of file for definition of CONVO */
384
385 /* INFO_message("about to mri_free(bim)") ; */
386 mri_free(bim) ;
387 if( fim == NULL ) return ijkout ;
388
389 /* find automask of the convolved image? */
390 #ifdef USE_FMASK
391 fmask = mri_automask_image( fim ) ;
392 if( fmask == NULL ){ mri_free(fim) ; return ijkout ; } /* bad */
393 #endif
394
395 far = MRI_FLOAT_PTR(fim) ;
396 nx = fim->nx ; ny = fim->ny ; nz = fim->nz ; nvox = nx*ny*nz ;
397 fc = xc = yc = zc = ftop = 0.0f ;
398
399 /* find the max convolution value */
400
401 for( ijk=0 ; ijk < nvox ; ijk++ ){
402 #ifdef USE_FMASK
403 if( !fmask[ijk] ) far[ijk] = 0.0f ; /* use the automask */
404 #endif
405 if( far[ijk] > ftop ) ftop = far[ijk] ; /* find the max */
406 }
407 #ifdef USE_FMASK
408 /* INFO_message("about to free(fmask)") ; */
409 free(fmask) ;
410 #endif
411 if( ftop == 0.0f ){
412 mri_free(fim) ; return ijkout ; /* bad */
413 }
414 ftop = 0.09f*ftop ; /* threshold for using voxel */
415
416 /* compute CM */
417
418 for( ijk=kk=0 ; kk < nz ; kk++ ){
419 for( jj=0 ; jj < ny ; jj++ ){
420 for( ii=0 ; ii < nx ; ii++,ijk++ ){
421 ff = far[ijk] ;
422 if( ff >= ftop ){
423 #ifdef USE_SQRTF
424 ff = sqrtf(ff); /* downweight */
425 #endif
426 xc += ff*ii; yc += ff*jj; zc += ff*kk; fc += ff;
427 }
428 }}}
429 /* INFO_message("about to mri_free(fim)") ; */
430 mri_free(fim) ;
431 if( fc == 0.0f ) return ijkout ; /* no positive values? */
432
433 /* subtract off ball radius from CM to get location in original dataset */
434
435 ijkout.i = (int)rintf(xc/fc) - irad ;
436 ijkout.j = (int)rintf(yc/fc) - jrad ;
437 ijkout.k = (int)rintf(zc/fc) - krad ;
438 ijkout.a = fc ;
439
440 /*
441 ININFO_message("xc=%g yc=%g zc=%g irad=%d jrad=%d krad=%d",
442 xc/fc , yc/fc , zc/fc , irad,jrad,krad ) ;
443 */
444
445 return ijkout ;
446 }
447
448 /*---------------------------------------------------------------------------*/
449 /* Find the 'best' overlap with various spheroids. */
450 /*---------------------------------------------------------------------------*/
451
THD_spheroid_overlapation(THD_3dim_dataset * dset,float arad,float_triple avec,float brad,float_triple bvec,float dth,int nthbot,int nthtop)452 int_triple THD_spheroid_overlapation( THD_3dim_dataset *dset ,
453 float arad , float_triple avec ,
454 float brad , float_triple bvec ,
455 float dth , int nthbot, int nthtop )
456 {
457 MRI_IMAGE *aim ; byte *aar ; int ii ;
458 int_triple ijkout = { -666 , -666 , -666 } ;
459 float dx,dy,dz,dxyz, dmin, crad ; int npeel ;
460 float_quint ab ; float_triple de ;
461 float ax,ay,az,aqq , bx,by,bz,bqq ;
462 mat33 mrot ; int irot ; float cx,cy,cz ;
463 int3float1 ijka ;
464 float isum,jsum,ksum,fsum,fff ;
465
466 if( !ISVALID_DSET(dset) || arad <= 0.0f ) return ijkout ;
467
468 dx = DSET_DX(dset); if( dx == 0.0f ) dx = 1.0f ; else dx = fabsf(dx) ;
469 dy = DSET_DY(dset); if( dy == 0.0f ) dy = 1.0f ; else dy = fabsf(dy) ;
470 dz = DSET_DZ(dset); if( dz == 0.0f ) dz = 1.0f ; else dz = fabsf(dz) ;
471 dmin = MIN(dx,dy) ; if( dz < dmin ) dmin = dz ;
472 de.a = dx ; de.b = dy ; de.c = dz ;
473
474 /* spherical case? */
475
476 if( brad == 0.0f || fabsf(arad-brad) <= 3.0f*dmin )
477 return THD_ball_mask_overlapation( dset , arad ) ;
478
479 /* autmask the dataset */
480
481 crad = sqrtf(arad*brad) ;
482 dxyz = cbrtf(dx*dy*dz) ;
483 npeel = (int)rintf( 0.266f*crad/dxyz ) ;
484 if( npeel <= 0 ) npeel = 1 ;
485 THD_automask_set_peelcounts(npeel,17) ;
486 THD_automask_set_clipfrac(0.135f) ;
487 THD_automask_extclip(1) ;
488 aar = THD_automask(dset) ; /* byte array with in=1 out=0 */
489
490 /* stuff automask into an MRI_IMAGE */
491
492 #ifndef USE_PM
493 /* use the 0-1 mask */
494 aim = mri_empty_conforming( DSET_BRICK(dset,0) , MRI_byte ) ;
495 mri_set_data_pointer( aim , aar ) ;
496 #else
497 /* convert to a -1/+1 mask */
498 { float *nar ; int nvox ;
499 aim = mri_new_conforming( DSET_BRICK(dset,0), MRI_float ) ;
500 nar = MRI_FLOAT_PTR(aim) ; nvox = aim->nvox ;
501 for( ii=0 ; ii < nvox ; ii++ ){
502 nar[ii] = ( aar[ii] ) ? +1.0f : -1.0f ; /* in=+1 out=-1 */
503 }
504 free(aar) ;
505 }
506 #endif
507 aim->dx = dx ; aim->dy = dy ; aim->dz = dz ;
508
509 /* normalize (ax,ay,az) to unit length */
510
511 ax = avec.a ; ay = avec.b ; az = avec.c ;
512 aqq = sqrtf(ax*ax + ay*ay + az*az) ; /* length of vector */
513 if( aqq == 0.0f ){
514 ay = 1.0f ; ax = az = 0.0f ; /* hopeless user */
515 } else if( aqq != 1.0f ){
516 ax = ax/aqq ; ay = ay/aqq ; az = az/aqq ; /* normalize */
517 }
518
519 /* normalize (bx,by,bz) to unit length */
520
521 bx = bvec.a ; by = bvec.b ; bz = bvec.c ;
522 bqq = sqrtf(bx*bx + by*by + bz*bz) ; /* length of vector */
523 if( bqq == 0.0f ){
524 bx = 1.0f ; by = bz = 0.0f ; /* hopeless user */
525 } else if( bqq != 1.0f ){
526 bx = bx/bqq ; by = by/bqq ; bz = bz/bqq ; /* normalize */
527 }
528
529 /* make sure bvec is perpendicular to avec */
530
531 aqq = ax*bx + ay*by + az*bz ; /* dot product of the 2 vectors */
532
533 if( fabsf(aqq) > 0.99f ){ /* a and b vectors are parallel :( */
534 bx = drand48() ; by = drand48() ; bz = drand48() ; /* so make up a */
535 bqq = sqrtf(bx*bx + by*by + bz*bz) ; /* new b vector */
536 bx = bx/bqq ; by = by/bqq ; bz = bz/bqq ;
537 aqq = ax*bx + ay*by + az*bz ;
538 }
539
540 if( fabsf(aqq) > 0.0001f ){ /* a and b are not perpendicular */
541 bx = bx - aqq*ax ; by = by - aqq*ay ; bz = bz - aqq*az ;
542 bqq = sqrtf(bx*bx + by*by + bz*bz) ;
543 bx = bx/bqq ; by = by/bqq ; bz = bz/bqq ; /* normalize */
544 }
545
546 /* loop over rotations of (ax,ay,az) axis about the (bx,by,bz) axis */
547
548 isum = jsum = ksum = fsum = 0.0f ;
549
550 for( irot=nthbot ; irot <= nthtop ; irot++ ){
551
552 /* rotate (ax,ay,az) to (cx,cy,cz) */
553 mrot = THD_mat33_generic_rotation( irot*dth , bx,by,bz ) ;
554 MAT33_VEC( mrot , ax,ay,az , cx,cy,cz ) ;
555
556 /* compute 'best' location for this orientation */
557
558 ab.a = arad; ab.b = cx; ab.c = cy; ab.d = cz; ab.e = brad;
559 ijka = MRI_spheroid_overlapation1( aim , ab ) ;
560
561 /*
562 INFO_message("spheroid overlapation: theta=%g axis=%g %g %g bxis=%g %g %g",
563 irot*dth , cx,cy,cz , bx,by,bz) ;
564 ININFO_message(" ijka = %d %d %d %g",ijka.i,ijka.j,ijka.k,ijka.a) ;
565 ININFO_message(" rotation matrix:\n"
566 " %13.6f %13.6f %13.6f\n"
567 " %13.6f %13.6f %13.6f\n"
568 " %13.6f %13.6f %13.6f" ,
569 mrot.m[0][0] , mrot.m[0][1] , mrot.m[0][2] ,
570 mrot.m[1][0] , mrot.m[1][1] , mrot.m[1][2] ,
571 mrot.m[2][0] , mrot.m[2][1] , mrot.m[2][2] ) ;
572 */
573
574 /* make running sum */
575 fff = ijka.a ;
576 if( fff > 0.0f ){
577 isum += fff * ijka.i ;
578 jsum += fff * ijka.j ;
579 ksum += fff * ijka.k ;
580 fsum += fff ;
581 }
582 }
583
584 /* compute the 'best' mask overlap location and return it */
585
586 mri_free(aim) ;
587
588 if( fsum > 0.0f ){
589 ijkout.i = (int)rintf(isum/fsum) ;
590 ijkout.j = (int)rintf(jsum/fsum) ;
591 ijkout.k = (int)rintf(ksum/fsum) ;
592 }
593
594 return ijkout ;
595 }
596
597 /*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*/
598 /*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*/
599
600 #ifdef USE_FFTN
601 /*----------------------------------------------------------------------------*/
602 /*----------------------------------------------------------------------------*/
603 /*** Some specialized OpenMP-ized 3D FFT functions [Oct 2020] ***/
604
605 /*----------------------------------------------------------------------------*/
606 /*----------------------------------------------------------------------------*/
607
608 #ifdef USE_OMP
609 # include <omp.h>
610 #endif
611 #include "Aomp.h"
612
613 #include "fftn_OMP.c"
614
615 /*----------------------------------------------------------------------------*/
616 /* Replaces mri_fft_3D() in mri_fft_complex.c -- note that alt arg is ignored
617 FFT lengths are in Lxx, Lyy, Lzz; however,
618 Lxx = 0 ==> no FFT in that direction (etc.).
619 *//*--------------------------------------------------------------------------*/
620
mri_fft_3D_OMP(int Sign,MRI_IMAGE * inim,int Lxx,int Lyy,int Lzz,int alt)621 static MRI_IMAGE * mri_fft_3D_OMP( int Sign, MRI_IMAGE *inim,
622 int Lxx,int Lyy,int Lzz, int alt )
623 {
624 MRI_IMAGE *outim ;
625 int ii,jj,kk , nx,ny,nxy,nz , nbig , fx,fy,fz,fxy ;
626 complex *car , *far ;
627 AO_DEFINE_ARRAY(complex,cbig) ;
628
629 if( inim->kind != MRI_complex ) return NULL ;
630
631 /* input data and its dimensions */
632
633 car = MRI_COMPLEX_PTR(inim) ;
634 nx = inim->nx ; ny = inim->ny ; nz = inim->nz ; nxy = nx*ny ;
635
636 /* output dimensions and data */
637
638 fx = (Lxx == 0) ? nx : (Lxx > nx) ? NEXTUP(Lxx) : NEXTUP(nx);
639 fy = (Lyy == 0) ? ny : (Lyy > ny) ? NEXTUP(Lyy) : NEXTUP(ny);
640 fz = (Lzz == 0) ? nz : (Lzz > nz) ? NEXTUP(Lzz) : NEXTUP(nz);
641 fxy = fx*fy ;
642
643 outim = mri_new_vol( fx,fy,fz , MRI_complex ) ; /* zero filled */
644 far = MRI_COMPLEX_PTR(outim) ;
645
646 /* buffer space */
647
648 nbig = MAX(fx,fy) ; nbig = MAX(nbig,fz) ; nbig = 4*nbig + 512 ;
649 #pragma omp parallel
650 { AO_RESIZE_ARRAY(complex,cbig,nbig) ;
651 /* fprintf(stderr,"resize cbig#%d to %d\n",AOth,nbig) ; */
652 }
653
654 /* copy input data into output image */
655
656 /* INFO_message("about to copy data from car to far") ; */
657 for( kk=0 ; kk < nz ; kk++ )
658 for( jj=0 ; jj < ny ; jj++ )
659 memcpy( far + jj*fx + kk*fxy, car + jj*nx + kk*nxy, sizeof(complex)*nx );
660
661 /* x-direction FFTs */
662
663 AFNI_OMP_START ;
664 if( Lxx > 1 ){
665 /* INFO_message("start x FFTs") ; */
666 #pragma omp parallel
667 { int kk,koff , jj,joff , ii ; complex *cbig ;
668 cbig = AO_VALUE(cbig) ;
669 #pragma omp for
670 for( kk=0 ; kk < fz ; kk++ ){
671 koff = kk*fxy ;
672 for( jj=0 ; jj < fy ; jj++ ){
673 /* ININFO_message("jj=%d kk=%d",jj,kk) ; */
674 joff = koff + jj*fx ;
675 for( ii=0 ; ii < fx ; ii++ ) cbig[ii] = far[ii+joff] ;
676 csfft_OMP( Sign , fx , cbig ) ;
677 for( ii=0 ; ii < fx ; ii++ ) far[ii+joff] = cbig[ii] ;
678 }
679 }
680 }
681 /* INFO_message("%d x FFTs(%d) done",fy*fz,fx) ; */
682 }
683 AFNI_OMP_END ;
684
685 /* y-direction FFTs */
686
687 AFNI_OMP_START ;
688 if( Lyy > 1 ){
689 /* INFO_message("start y FFTs") ; */
690 #pragma omp parallel
691 { int kk,koff , jj,joff , ii ; complex *cbig ;
692 cbig = AO_VALUE(cbig) ;
693 #pragma omp for
694 for( kk=0 ; kk < fz ; kk++ ){
695 koff = kk*fxy ;
696 for( ii=0 ; ii < fx ; ii++ ){
697 /* ININFO_message("ii=%d kk=%d",ii,kk) ; */
698 joff = koff + ii ;
699 for( jj=0 ; jj < fy ; jj++ ) cbig[jj] = far[jj*fx+joff] ; /* copy data */
700 csfft_OMP( Sign , fy , cbig ) ; /* FFT in buffer */
701 for( jj=0 ; jj < fy ; jj++ ) far[jj*fx+joff] = cbig[jj] ; /* copy back */
702 }
703 }
704 }
705 /* INFO_message("%d y FFTs(%d) done",fx*fz,fy) ; */
706 }
707 AFNI_OMP_END ;
708
709 /* z-direction FFTs */
710
711 AFNI_OMP_START ;
712 if( Lzz > 1 ){
713 /* INFO_message("start z FFTs: jj=0..%d ii=0..%d fz=%d",fy-1,fx-1,fz) ; */
714 #pragma omp parallel
715 { int kk,koff , jj,joff , ii ; complex *cbig ;
716 cbig = AO_VALUE(cbig) ;
717 #pragma omp for
718 for( jj=0 ; jj < fy ; jj++ ){
719 joff = jj*fx ;
720 for( ii=0 ; ii < fx ; ii++ ){
721 /* ININFO_message("ii=%d jj=%d",ii,jj) ; */
722 koff = joff + ii ;
723 for( kk=0 ; kk < fz ; kk++ ) cbig[kk] = far[kk*fxy+koff] ;
724 csfft_OMP( Sign , fz , cbig ) ;
725 for( kk=0 ; kk < fz ; kk++ ) far[kk*fxy+koff] = cbig[kk] ;
726 }
727 }
728 }
729 /* INFO_message("%d z FFTs(%d) done",fx*fy,fz) ; */
730 }
731 AFNI_OMP_END ;
732
733 /* INFO_message("all FFTs done") ; */
734 MRI_COPY_AUX(outim,inim) ; return outim ;
735 }
736
737 /*----------------------------------------------------------------------------*/
738 /* Replaces mri_fft_3Dconvolve() in mri_fft_complex.c [Oct 2020]
739 Convolve (via FFT) image aim with bim.
740 Note output image will be at least as big as than the sum of the two sizes.
741 *//*--------------------------------------------------------------------------*/
742
mri_fft_3Dconvolve_OMP(MRI_IMAGE * aim,MRI_IMAGE * bim)743 static MRI_IMAGE * mri_fft_3Dconvolve_OMP( MRI_IMAGE *aim , MRI_IMAGE *bim )
744 {
745 MRI_IMAGE *outim=NULL ;
746 MRI_IMAGE *paim , *pbim , *faim , *fbim ;
747 int nxa,nya,nza , nxb,nyb,nzb , Lxx,Lyy,Lzz , Lxyz,ii ;
748 complex ac , bc , qc ;
749 complex *acar , *bcar ;
750 float linv ;
751
752 if( aim == NULL || bim == NULL ) return NULL ;
753
754 /* input dimensions */
755
756 nxa = aim->nx ; nya = aim->ny ; nza = aim->nz ;
757 nxb = bim->nx ; nyb = bim->ny ; nzb = bim->nz ;
758
759 /* FFT and output dimensions (sum, bumped up for FFT effiency) */
760
761 Lxx = (nxa > 1 && nxb > 1) ? NEXTUP(nxa+nxb) : 0 ;
762 Lyy = (nya > 1 && nyb > 1) ? NEXTUP(nya+nyb) : 0 ;
763 Lzz = (nza > 1 && nzb > 1) ? NEXTUP(nza+nzb) : 0 ;
764
765 /* at this time, we don't allow for convolving a 3D image with a 1D
766 or 2D image, for example, which is possible but more complicated */
767
768 if( Lxx == 0 || Lyy == 0 || Lzz == 0 ) return NULL ;
769
770 /* 1) convert A image to complex
771 2) zero pad it to fit the FFT size
772 3) FFT that
773 Then repeat these steps for the B image */
774
775 faim = mri_to_complex( aim ) ; /* 1) */
776 paim = mri_zeropad_3D( 0,Lxx-nxa , 0,Lyy-nya , 0,Lzz-nza , faim ) ; /* 2) */
777 /* INFO_message("about to mri_free(faim)") ; */
778 mri_free(faim) ;
779 faim = mri_fft_3D_OMP( -1 , paim , Lxx,Lyy,Lzz , 0 ) ; /* 3) */
780 /* INFO_message("about to mri_free(paim)") ; */
781 mri_free(paim) ;
782 acar = MRI_COMPLEX_PTR(faim) ;
783
784 fbim = mri_to_complex( bim ) ; /* 1) */
785 pbim = mri_zeropad_3D( 0,Lxx-nxb , 0,Lyy-nyb , 0,Lzz-nzb , fbim ) ; /* 2) */
786 /* INFO_message("about to mri_free(fbim)") ; */
787 mri_free(fbim) ;
788 fbim = mri_fft_3D_OMP( -1 , pbim , Lxx,Lyy,Lzz , 0 ) ; /* 3) */
789 /* INFO_message("about to mri_free(pbim)") ; */
790 mri_free(pbim) ;
791 bcar = MRI_COMPLEX_PTR(fbim) ;
792
793 /* multiply+scale FFTs, store back in faim/acar */
794
795 Lxyz = Lxx * Lyy * Lzz ;
796 linv = 10.f / (float)Lxyz ; /* scaling for inverse FFT */
797 for( ii=0 ; ii < Lxyz ; ii++ ){
798 ac = acar[ii] ;
799 bc = bcar[ii] ;
800 qc.r = (ac.r * bc.r - ac.i * bc.i) * linv ; /* complex */
801 qc.i = (ac.r * bc.i + ac.i * bc.r) * linv ; /* multiply */
802 acar[ii] = qc ;
803 }
804 /* INFO_message("about to mri_free(fbim) again") ; */
805 mri_free(fbim) ;
806
807 /* inverse FFT back to 'real' space */
808
809 fbim = mri_fft_3D_OMP( +1 , faim , Lxx,Lyy,Lzz , 0 ) ;
810 /* INFO_message("about to mri_free(faim) again") ; */
811 mri_free(faim) ;
812
813 /* convert to float-valued image (from complex FFT) and return */
814
815 outim = mri_complex_to_real( fbim ) ;
816 /* INFO_message("about to mri_free(faim) yet again") ; */
817 mri_free(fbim) ;
818 return outim ;
819 }
820
821 #endif /* USE_FFTN */
822