1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 /*---------------------------------------------------------------------------*/
8
GA_gcd(int m,int n)9 int GA_gcd( int m , int n ) /* Euclid's Greatest Common Denominator */
10 {
11 while( m > 0 ){
12 if( n > m ){ int t=m; m=n; n=t; } /* swap */
13 m -= n ;
14 }
15 return n ;
16 }
17
18 /*---------------------------------------------------------------------------*/
19
GA_find_relprime_fixed(int n)20 int GA_find_relprime_fixed( int n ) /* find number relatively prime to n */
21 {
22 int dj , n5=n/5 ;
23 if( n5 < 2 ) return 1 ;
24 for( dj=n5 ; GA_gcd(n,dj) > 1 ; dj++ ) ; /*nada*/
25 return dj ;
26 }
27
28 /*---------------------------------------------------------------------------*/
29 /*! Smooth an image with a given method to a given radius.
30 Assumes the dx,dy,dz parameters in the image struct are correct! */
31
GA_smooth(MRI_IMAGE * im,int meth,float rad)32 MRI_IMAGE * GA_smooth( MRI_IMAGE *im , int meth , float rad )
33 {
34 MRI_IMAGE *om=NULL ;
35
36 ENTRY("GA_smooth") ;
37
38 if( im == NULL || rad <= 0.0f ) RETURN(NULL) ;
39
40 #undef CALLME
41 #define CALLME(inp,out) (out) = GA_smooth( (inp), meth,rad )
42 if( ISVECTIM(im) ){ VECTORME(im,om) ; RETURN(om) ; }
43
44 if( im->kind != MRI_float ) RETURN(NULL) ;
45
46 RAND_ROUND ;
47 switch( meth ){
48 default:
49 case GA_SMOOTH_GAUSSIAN:
50 om = mri_to_float(im) ;
51 rad = FWHM_TO_SIGMA(rad) ; /* convert rad from FWHM to st.dev. */
52 FIR_blur_volume_3d( om->nx , om->ny , om->nz ,
53 om->dx , om->dy , om->dz ,
54 MRI_FLOAT_PTR(om) , rad,rad,rad ) ;
55 break ;
56
57 case GA_SMOOTH_MEDIAN:{
58 float d;
59 d = MIN(im->dx,im->dy) ; d = MIN(d,im->dz) ;
60 if( rad <= 1.01f*d ) rad = 1.01f*d ;
61 if( rad > 0.0f ){ mri_medianfilter_usedxyz(1); }
62 else { mri_medianfilter_usedxyz(0); rad=1.01f; }
63 om = mri_medianfilter( im , rad , NULL , 0 ) ;
64 }
65 break ;
66 }
67
68 RETURN(om) ;
69 }
70
71 /*---------------------------------------------------------------------------*/
72
73 /* for interpolation access to the (i,j,k) element of an array */
74
75 #undef FAR
76 #define FAR(i,j,k) far[(i)+(j)*nx+(k)*nxy]
77
78 #undef FARJK
79 #define FARJK(j,k) (far+(j)*nx+(k)*nxy)
80
81 /* clip value mm to range 0..nn */
82
83 #undef CLIP
84 #define CLIP(mm,nn) if(mm < 0)mm=0; else if(mm > nn)mm=nn
85
86 static float outval = 0.0f ; /* value for 'outside' voxels */
GA_set_outval(float v)87 void GA_set_outval( float v ){ outval = v; } /* 28 Feb 2007 */
GA_get_outval(void)88 float GA_get_outval(void){ return outval; } /* 10 Dec 2008 */
89
90 #undef ISTINY
91 #define ISTINY(a) (fabsf(a) < 0.0001f)
92
93 /*---------------------------------------------------------------*/
94 /*! Interpolate an image at npp (index) points, using NN method. */
95
GA_interp_NN(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)96 void GA_interp_NN( MRI_IMAGE *fim ,
97 int npp, float *ip, float *jp, float *kp, float *vv )
98 {
99 ENTRY("GA_interp_NN") ;
100
101 RAND_ROUND ;
102
103 AFNI_OMP_START ;
104 #pragma omp parallel if(npp > 6666)
105 {
106 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , ii,jj,kk , pp ;
107 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
108 float *far = MRI_FLOAT_PTR(fim) ;
109 #pragma omp for
110 for( pp=0 ; pp < npp ; pp++ ){
111 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
112 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
113 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
114
115 ii = (int)(xx+0.5f) ; jj = (int)(yy+0.5f) ; kk = (int)(zz+0.5f) ;
116 vv[pp] = FAR(ii,jj,kk) ;
117 }
118 } /* end OpenMP */
119 AFNI_OMP_END ;
120
121 EXRETURN ;
122 }
123
124 /*---------------------------------------------------------------------------*/
125 /*! Interpolate an image at npp (index) points, using linear method. */
126
GA_interp_linear(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)127 void GA_interp_linear( MRI_IMAGE *fim ,
128 int npp, float *ip, float *jp, float *kp, float *vv )
129 {
130 ENTRY("GA_interp_linear") ;
131
132 #if 0
133 if( fim == NULL || ip == NULL || jp == NULL || kp == NULL || vv == NULL )
134 ERROR_message("NULL pointer on entry to GA_interp_linear :-(") ;
135 #endif
136
137 RAND_ROUND ; /* 26 Feb 2020 */
138
139 AFNI_OMP_START ;
140 #pragma omp parallel if(npp > 4444)
141 {
142 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
143 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
144 float fx,fy,fz ;
145 float *far = MRI_FLOAT_PTR(fim) ;
146 int nx1=nx-1,ny1=ny-1,nz1=nz-1 ;
147 float ix,jy,kz ;
148 int ix_00,ix_p1 ; /* interpolation indices */
149 int jy_00,jy_p1 ;
150 int kz_00,kz_p1 ;
151 float wt_00,wt_p1 ; /* interpolation weights */
152 float f_j00_k00, f_jp1_k00, f_j00_kp1, f_jp1_kp1, f_k00, f_kp1 ;
153
154 #pragma omp for
155 for( pp=0 ; pp < npp ; pp++ ){
156 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
157 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
158 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
159
160 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
161 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
162 kz = floorf(zz) ; fz = zz - kz ;
163
164 ix_00 = ix ; ix_p1 = ix_00+1 ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ;
165 jy_00 = jy ; jy_p1 = jy_00+1 ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ;
166 kz_00 = kz ; kz_p1 = kz_00+1 ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ;
167
168 wt_00 = 1.0f-fx ; wt_p1 = fx ; /* weights for ix_00 and ix_p1 points */
169
170 #undef XINT
171 #define XINT(j,k) wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k)
172
173 /* interpolate to location ix+fx at each jy,kz level */
174
175 f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
176 f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
177
178 /* interpolate to jy+fy at each kz level */
179
180 wt_00 = 1.0f-fy ; wt_p1 = fy ;
181 f_k00 = wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 ;
182 f_kp1 = wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 ;
183
184 /* interpolate to kz+fz to get output */
185
186 vv[pp] = (1.0f-fz) * f_k00 + fz * f_kp1 ;
187 }
188
189 } /* end OpenMP */
190 AFNI_OMP_END ;
191
192 EXRETURN ;
193 }
194
195 /*---------------------------------------------------------------------------*/
196
197 /* define Lagrange cubic interpolation polynomials */
198
199 #undef P_M1
200 #undef P_00
201 #undef P_P1
202 #undef P_P2
203 #undef P_FACTOR
204 #define P_M1(x) (-(x)*((x)-1)*((x)-2))
205 #define P_00(x) (3*((x)+1)*((x)-1)*((x)-2))
206 #define P_P1(x) (-3*(x)*((x)+1)*((x)-2))
207 #define P_P2(x) ((x)*((x)+1)*((x)-1))
208 #define P_FACTOR 4.62962963e-3 /* 1/216 = final scaling factor */
209
210 /*------------------------------------------------------------------*/
211 /*! Interpolate an image at npp (index) points, using cubic method. */
212
GA_interp_cubic(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)213 void GA_interp_cubic( MRI_IMAGE *fim ,
214 int npp, float *ip, float *jp, float *kp, float *vv )
215 {
216 ENTRY("GA_interp_cubic") ;
217
218 RAND_ROUND ; /* 26 Feb 2020 */
219
220 AFNI_OMP_START ;
221 #pragma omp parallel if(npp > 2222)
222 {
223 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
224 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
225 float fx,fy,fz ;
226 float *far = MRI_FLOAT_PTR(fim) ;
227 int nx1=nx-1,ny1=ny-1,nz1=nz-1, ix,jy,kz ;
228
229 int ix_m1,ix_00,ix_p1,ix_p2 ; /* interpolation indices */
230 int jy_m1,jy_00,jy_p1,jy_p2 ;
231 int kz_m1,kz_00,kz_p1,kz_p2 ;
232 float wt_m1,wt_00,wt_p1,wt_p2 ; /* interpolation weights */
233 float f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, /* interpolants */
234 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00,
235 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1,
236 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2,
237 f_km1 , f_k00 , f_kp1 , f_kp2 ;
238
239 #pragma omp for
240 for( pp=0 ; pp < npp ; pp++ ){
241 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
242 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
243 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
244
245 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
246 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
247 kz = floorf(zz) ; fz = zz - kz ;
248
249 if( ISTINY(fx) && ISTINY(fy) && ISTINY(fz) ){
250 CLIP(ix,nx1); CLIP(jy,ny1); CLIP(kz,nz1);
251 vv[pp] = FAR(ix,jy,kz) ; continue ;
252 }
253
254 ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ;
255 CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
256
257 jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ;
258 CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
259
260 kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ;
261 CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
262
263 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ; /* interpolation weights */
264 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
265
266 #undef XINT
267 #define XINT(j,k) wt_m1*FAR(ix_m1,j,k)+wt_00*FAR(ix_00,j,k) \
268 +wt_p1*FAR(ix_p1,j,k)+wt_p2*FAR(ix_p2,j,k)
269
270 /* interpolate to location ix+fx at each jy,kz level */
271
272 f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ;
273 f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ;
274 f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ;
275 f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ;
276 f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ;
277 f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ;
278 f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ;
279 f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ;
280
281 /* interpolate to jy+fy at each kz level */
282
283 wt_m1 = P_M1(fy) ; wt_00 = P_00(fy) ;
284 wt_p1 = P_P1(fy) ; wt_p2 = P_P2(fy) ;
285
286 f_km1 = wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
287 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 ;
288 f_k00 = wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
289 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 ;
290 f_kp1 = wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
291 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 ;
292 f_kp2 = wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
293 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 ;
294
295 /* interpolate to kz+fz to get output */
296
297 wt_m1 = P_M1(fz) ; wt_00 = P_00(fz) ;
298 wt_p1 = P_P1(fz) ; wt_p2 = P_P2(fz) ;
299
300 vv[pp] = P_FACTOR * ( wt_m1 * f_km1 + wt_00 * f_k00
301 + wt_p1 * f_kp1 + wt_p2 * f_kp2 ) ;
302 }
303
304 } /* end OpenMP */
305 AFNI_OMP_END ;
306
307 EXRETURN ;
308 }
309
310 /*---------------------------------------------------------------------------*/
311 /* define variance preserving interpolation functions */
312
313 #undef BP
314 #define BP 0.7611164839f /* = (11+sqrt(33))/22 */
315 #undef PSI
316 #define PSI (xb>-0.001f && xb<0.001f) \
317 ? (1.2222222f+0.52008386f*xb)*xb \
318 : 0.16666667f*xb \
319 * sqrtf( ( 28.0f*xxx*(xxx-1.0f)+10.0f+(6.0f-12.0f*xxx)*phi ) \
320 / (xb*xb) )
321 #undef VPWT
322 #define VPWT(x) \
323 { float xxx=(x) , xb , px , psi ; \
324 float phi = sqrtf(1.0f-8.0f*xxx*(xxx-1.0f)) ; \
325 wt_m1 = wt_p2 = 0.25f * (1.0f-phi) ; px = 0.25f + 0.08333333f*phi ; \
326 xb = xxx-BP ; psi = PSI ; wt_00 = px - psi ; \
327 xxx = 1.0f-xxx ; xb = xxx-BP ; psi = PSI ; wt_p1 = px - psi ; }
328
329 /*------------------------------------------------------------------*/
330 /*! Interpolate an image at npp (index) points, using VP method. */
331
GA_interp_varp1(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)332 void GA_interp_varp1( MRI_IMAGE *fim ,
333 int npp, float *ip, float *jp, float *kp, float *vv )
334 {
335 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
336 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
337 float fx,fy,fz ;
338 float *far = MRI_FLOAT_PTR(fim) ;
339 int nx1=nx-1,ny1=ny-1,nz1=nz-1, ix,jy,kz ;
340
341 int ix_m1,ix_00,ix_p1,ix_p2 ; /* interpolation indices */
342 int jy_m1,jy_00,jy_p1,jy_p2 ;
343 int kz_m1,kz_00,kz_p1,kz_p2 ;
344 float wt_m1,wt_00,wt_p1,wt_p2 ; /* interpolation weights */
345 float f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, /* interpolants */
346 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00,
347 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1,
348 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2,
349 f_km1 , f_k00 , f_kp1 , f_kp2 ;
350
351 ENTRY("GA_interp_varp1") ;
352
353 for( pp=0 ; pp < npp ; pp++ ){
354 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
355 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
356 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
357
358 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
359 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
360 kz = floorf(zz) ; fz = zz - kz ;
361
362 ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ;
363 CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
364
365 jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ;
366 CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
367
368 kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ;
369 CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
370
371 VPWT(fx) ; /* interpolation weights in x direction */
372
373 #undef XINT
374 #define XINT(j,k) wt_m1*FAR(ix_m1,j,k)+wt_00*FAR(ix_00,j,k) \
375 +wt_p1*FAR(ix_p1,j,k)+wt_p2*FAR(ix_p2,j,k)
376
377 /* interpolate to location ix+fx at each jy,kz level */
378
379 f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ;
380 f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ;
381 f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ;
382 f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ;
383 f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ;
384 f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ;
385 f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ;
386 f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ;
387
388 /* interpolate to jy+fy at each kz level */
389
390 VPWT(fy) ; /* interpolation weights in y direction */
391
392 f_km1 = wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
393 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 ;
394 f_k00 = wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
395 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 ;
396 f_kp1 = wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
397 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 ;
398 f_kp2 = wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
399 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 ;
400
401 /* interpolate to kz+fz to get output */
402
403 VPWT(fz) ; /* interpolation weights in z direction */
404
405 vv[pp] = wt_m1 * f_km1 + wt_00 * f_k00
406 + wt_p1 * f_kp1 + wt_p2 * f_kp2 ;
407 }
408 EXRETURN ;
409 }
410
411 /*============================================================================*/
412 /* Interpolation with weighted (tapered) sinc in 3D.
413 Funcion GA_interp_wsinc5s() uses a spherical mask (really slow).
414 Funcion GA_interp_wsinc5p() uses a cubical mask (pretty slow).
415 Function setup_wsinc5() allows control via environment variables (once):
416 AFNI_WSINC5_TAPERCUT = between 0 and 0.8, sets start point for tapering
417 (0 = most tapering, 0.8 = least tapering)
418 AFNI_WSINC5_RADIUS = integer from 3 to 21 (bigger = slower)
419 AFNI_WSINC5_TAPERFUN = 'H' for Hamming (2 terms); otherwise, uses 3 terms
420 AFNI_WSINC5_SPHERICAL = 'Y' for spherical mask ; otherwise uses cubical
421 AFNI_WSINC5_SILENT = 'Y' to avoid the message
422 *//*--------------------------------------------------------------------------*/
423
424 #undef PIF
425 #define PIF 3.1415927f /* PI in float */
426
427 static float WCUT = 0.0f ; /* cutoff point for taper function */
428 static float WCUTI = 1.00f ; /* = 1.0/(1.0-WCUT) */
429 static int IRAD = 5 ; /* radius of sinc window */
430 static int IRAD1 = 4 ; /* IRAD - 1 */
431 static float WRAD = 5.001f ; /* float form of IRAD (+ epsilon) */
432 static int WFUN = 0 ; /* window chooser: 0 = M3(x) ; 1 = HW(x) */
433 static int WSHAP = 0 ; /* 0 = cubical ; 1 = spherical */
434
435 #define IRAD_MAX 21 /* largest IRAD allowed */
436
437 /*------------------------------------------------------------------*/
438
setup_wsinc5(void)439 static void setup_wsinc5(void)
440 {
441 char *eee ; float val ;
442
443 eee = getenv("AFNI_WSINC5_TAPERCUT") ;
444 WCUT = 0.0f ;
445 if( eee != NULL ){
446 val = (float)strtod(eee,NULL) ;
447 if( val >= 0.0f && val <= 0.8f ) WCUT = val ;
448 }
449 WCUTI = 1.0f / (1.0f - WCUT) ;
450
451 eee = getenv("AFNI_WSINC5_RADIUS") ;
452 IRAD = 5 ;
453 if( eee != NULL ){
454 val = (float)strtod(eee,NULL) ;
455 if( val >= 3.0f && val <= (IRAD_MAX+0.9f) ) IRAD = (int)val ;
456 }
457 WRAD = 0.001f + (float)IRAD ; IRAD1 = IRAD - 1 ;
458
459 eee = getenv("AFNI_WSINC5_TAPERFUN") ;
460 WFUN = (eee != NULL && toupper(*eee) == 'H' ) ;
461
462 eee = getenv("AFNI_WSINC5_SPHERICAL") ;
463 WSHAP = (eee != NULL && toupper(*eee) == 'Y' ) ;
464
465 eee = getenv("AFNI_WSINC5_SILENT") ;
466 if( eee == NULL || toupper(*eee) != 'Y' ){
467 INFO_message("wsinc5 interpolation setup:") ;
468 ININFO_message(" taper function = %s",(WFUN)?"Hamming":"Min sidelobe 3 term");
469 ININFO_message(" taper cut point = %.3f",WCUT) ;
470 ININFO_message(" window radius = %d voxels",IRAD) ;
471 ININFO_message(" window shape = %s",(WSHAP)?"Spherical":"Cubical") ;
472 ININFO_message(" The above can be altered via the AFNI_WSINC5_* environment variables.") ;
473 ININFO_message(" (To avoid this message, 'setenv AFNI_WSINC5_SILENT YES'.)") ;
474 }
475
476 return ;
477 }
478
479 /*------------------------------------------------------------------*/
480
481 /* sinc function = sin(PI*x)/(PI*x) [N.B.: x will always be >= 0] */
482
483 #undef sinc
484 #define sinc(x) ( ((x)>0.01f) ? sinf(PIF*(x))/(PIF*(x)) \
485 : 1.0f - 1.6449341f*(x)*(x) )
486
487 /*------------------------------------------------------------------*/
488
489 /* HW(x) = Hamming Window = minimum sidelobe 2 term window */
490
491 #undef HW
492 #define HW(x) (0.53836f+0.46164f*cosf(PIF*(x)))
493
494 /*------------------------------------------------------------------*/
495
496 /* M3(x) = minimum sidelobe 3 term window (has no catchy name, alas) */
497
498 #undef M3
499 #define M3(x) (0.4243801f+0.4973406f*cosf(PIF*(x))+0.0782793f*cosf(PIF*(x)*2.0f))
500
501 /*------------------------------------------------------------------*/
502
503 /* Weight (taper) function, declining from ww(WCUT)=1 to ww(1)=0 */
504 /* Note that the input to ww will always be between WCUT and 1. */
505
506 #undef ww
507 #define ww(x) ( (WFUN) ? HW( ((x)-WCUT)*WCUTI ) : M3( ((x)-WCUT)*WCUTI ) )
508
509 /*---------------------------------------------------------------------------*/
510
511 #define UNROLL /* unroll some loops */
512
513 /*---------------------------------------------------------------------------*/
514 /*! Spherical windowed sinc interpolation (really slow). */
515
GA_interp_wsinc5s(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)516 void GA_interp_wsinc5s( MRI_IMAGE *fim ,
517 int npp, float *ip, float *jp, float *kp, float *vv )
518 {
519 static MCW_cluster *smask=NULL ; static int nmask=0 ;
520 static short *di=NULL , *dj=NULL , *dk=NULL ;
521
522 ENTRY("GA_interp_wsinc5s") ;
523
524 /*----- first time in: build spherical mask -----*/
525
526 if( smask == NULL ){
527 char *eee ;
528 smask = MCW_spheremask( 1.0f,1.0f,1.0f , WRAD ) ;
529 nmask = smask->num_pt ;
530 di = smask->i ;
531 dj = smask->j ;
532 dk = smask->k ;
533 eee = getenv("AFNI_WSINC5_SILENT") ;
534 if( eee == NULL || toupper(*eee) != 'Y' )
535 ININFO_message(" wsinc5 SPHERE(%d) mask has %d points",IRAD,nmask) ;
536 }
537
538 /*----- loop over points -----*/
539 AFNI_OMP_START ;
540 #pragma omp parallel if(npp > 444)
541 {
542 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
543 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
544 float fx,fy,fz ;
545 float *far = MRI_FLOAT_PTR(fim) ;
546 int nx1=nx-1,ny1=ny-1,nz1=nz-1, ix,jy,kz ;
547 float xw,yw,zw,rr , sum,wsum,wt ;
548 int iq,jq,kq , qq , ddi,ddj,ddk ;
549 float xsin[1+2*IRAD_MAX] , ysin[1+2*IRAD_MAX] , zsin[1+2*IRAD_MAX] ;
550
551 #pragma omp for
552 for( pp=0 ; pp < npp ; pp++ ){
553 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
554 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
555 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
556
557 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
558 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
559 kz = floorf(zz) ; fz = zz - kz ;
560
561 if( ISTINY(fx) && ISTINY(fy) && ISTINY(fz) ){
562 CLIP(ix,nx1); CLIP(jy,ny1); CLIP(kz,nz1);
563 vv[pp] = FAR(ix,jy,kz) ; continue ;
564 }
565
566 /*- compute sinc at all points plus/minus 5 indexes from current locale -*/
567
568 for( qq=-IRAD ; qq <= IRAD ; qq++ ){
569 xw = fabsf(fx-qq) ; xsin[qq+IRAD] = sinc(xw) ;
570 yw = fabsf(fy-qq) ; ysin[qq+IRAD] = sinc(yw) ;
571 zw = fabsf(fz-qq) ; zsin[qq+IRAD] = sinc(zw) ;
572 }
573
574 for( wsum=sum=0.0f,qq=0 ; qq < nmask ; qq++ ){
575 ddi = di[qq] ; ddj = dj[qq] ; ddk = dk[qq] ;
576 iq = ix + ddi ; CLIP(iq,nx1) ; xw = fx - (float)ddi ;
577 jq = jy + ddj ; CLIP(jq,ny1) ; yw = fy - (float)ddj ;
578 kq = kz + ddk ; CLIP(kq,nz1) ; zw = fz - (float)ddk ;
579 rr = sqrtf(xw*xw+yw*yw+zw*zw) / WRAD ; if( rr >= 1.0f ) continue ;
580 wt = xsin[ddi+IRAD] * ysin[ddj+IRAD] * zsin[ddk+IRAD] ;
581 if( rr > WCUT ) wt *= ww(rr) ;
582 wsum += wt ; sum += FAR(iq,jq,kq) * wt ;
583 }
584
585 vv[pp] = sum / wsum ;
586 }
587 } /* end OpenMP */
588 AFNI_OMP_END ;
589
590 EXRETURN ;
591 }
592
593 /*---------------------------------------------------------------------------*/
594 /*! Cubical (tensor product) windowed sinc interpolation (pretty slow). */
595
GA_interp_wsinc5p(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)596 void GA_interp_wsinc5p( MRI_IMAGE *fim ,
597 int npp, float *ip, float *jp, float *kp, float *vv )
598 {
599 static int first=1 ;
600
601 ENTRY("GA_interp_wsinc5p") ;
602
603 if( first ){
604 char *eee = getenv("AFNI_WSINC5_SILENT") ;
605 if( eee == NULL || toupper(*eee) != 'Y' )
606 ININFO_message(" wsinc5 CUBE(%d) mask has %d points",IRAD,8*IRAD*IRAD*IRAD) ;
607 first = 0 ;
608 }
609
610 AFNI_OMP_START ;
611 #pragma omp parallel if(npp > 444)
612 {
613 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
614 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
615 float fx,fy,fz ;
616 float *far = MRI_FLOAT_PTR(fim) , *farjk ;
617 int nx1=nx-1,ny1=ny-1,nz1=nz-1, ix,jy,kz ;
618
619 float xw,yw,zw,rr , sum,wsum,wfac,wt ;
620 int iq,jq,kq,iqp , qq,jj,kk , ddi,ddj,ddk ;
621 float xsin[2*IRAD_MAX] , ysin[2*IRAD_MAX] , zsin[2*IRAD_MAX] ;
622 float wtt[2*IRAD_MAX] , fjk[2*IRAD_MAX][2*IRAD_MAX] , fk[2*IRAD_MAX] ;
623 int iqq[2*IRAD_MAX] ;
624
625 /*----- loop over points -----*/
626
627
628 #pragma omp for
629 for( pp=0 ; pp < npp ; pp++ ){
630 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
631 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
632 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
633
634 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
635 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
636 kz = floorf(zz) ; fz = zz - kz ;
637
638 if( ISTINY(fx) && ISTINY(fy) && ISTINY(fz) ){
639 CLIP(ix,nx1); CLIP(jy,ny1); CLIP(kz,nz1);
640 vv[pp] = FAR(ix,jy,kz); continue;
641 }
642
643 /*- x interpolations -*/
644
645 for( wsum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){ /* setup weights */
646 xw = fabsf(fx-qq) ; wt = sinc(xw) ;
647 xw /= WRAD ; if( xw > WCUT ) wt *= ww(xw) ;
648 wtt[qq+IRAD1] = wt ; wsum += wt ;
649 iq = ix+qq ; CLIP(iq,nx1) ; iqq[qq+IRAD1] = iq ;
650 }
651 wfac = wsum ;
652
653 for( jj=-IRAD1 ; jj <= IRAD ; jj++ ){
654 jq = jy+jj ; CLIP(jq,ny1) ;
655 for( kk=-IRAD1 ; kk <= IRAD ; kk++ ){
656 kq = kz+kk ; CLIP(kq,nz1) ;
657 #ifndef UNROLL
658 for( sum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){
659 iq = iqq[qq+IRAD1] ; sum += FAR(iq,jq,kq) * wtt[qq+IRAD1] ;
660 }
661 #else
662 # define FW(i) farjk[iqq[i]]*wtt[i]
663 farjk = FARJK(jq,kq) ;
664 switch( IRAD ){
665
666 default: /* not actually needed */
667 for( sum=0.0f,qq=-IRAD1 ; qq < IRAD ; qq+=2 ){ /* unrolled by 2 */
668 iq = iqq[qq+IRAD1] ; iqp = iqq[qq+IRAD] ;
669 sum += farjk[iq] * wtt[qq+IRAD1]
670 +farjk[iqp] * wtt[qq+IRAD ] ;
671 }
672 break ;
673
674 /* all possible cases from 3..21 are manually unrolled */
675 /* adding up FW(0)..FW(2*IRAD-1) */
676
677 case 3:
678 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5) ;
679 break ;
680
681 case 4:
682 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7) ;
683 break ;
684
685 case 5:
686 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9) ;
687 break ;
688
689 case 6:
690 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
691 +FW(10)+FW(11) ;
692 break ;
693
694 case 7:
695 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
696 +FW(10)+FW(11)+FW(12)+FW(13) ;
697 break ;
698
699 case 8:
700 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
701 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15) ;
702 break ;
703
704 case 9:
705 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
706 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17) ;
707 break ;
708
709 case 10:
710 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
711 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
712 +FW(18)+FW(19) ;
713 break ;
714
715 case 11:
716 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
717 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
718 +FW(18)+FW(19)+FW(20)+FW(21) ;
719 break ;
720
721 case 12:
722 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
723 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
724 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23) ;
725 break ;
726
727 case 13:
728 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
729 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
730 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25) ;
731 break ;
732
733 case 14:
734 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
735 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
736 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
737 +FW(26)+FW(27) ;
738 break ;
739
740 case 15:
741 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
742 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
743 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
744 +FW(26)+FW(27)+FW(28)+FW(29) ;
745 break ;
746
747 case 16:
748 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
749 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
750 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
751 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31) ;
752 break ;
753
754 case 17:
755 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
756 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
757 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
758 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31)+FW(32)+FW(33) ;
759 break ;
760
761 case 18:
762 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
763 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
764 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
765 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31)+FW(32)+FW(33)
766 +FW(34)+FW(35) ;
767 break ;
768
769 case 19:
770 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
771 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
772 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
773 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31)+FW(32)+FW(33)
774 +FW(34)+FW(35)+FW(36)+FW(37) ;
775 break ;
776
777 case 20:
778 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
779 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
780 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
781 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31)+FW(32)+FW(33)
782 +FW(34)+FW(35)+FW(36)+FW(37)+FW(38)+FW(39)+FW(40) ;
783 break ;
784
785 case 21:
786 sum = FW(0)+FW(1)+FW(2)+FW(3)+FW(4)+FW(5)+FW(6)+FW(7)+FW(8)+FW(9)
787 +FW(10)+FW(11)+FW(12)+FW(13)+FW(14)+FW(15)+FW(16)+FW(17)
788 +FW(18)+FW(19)+FW(20)+FW(21)+FW(22)+FW(23)+FW(24)+FW(25)
789 +FW(26)+FW(27)+FW(28)+FW(29)+FW(30)+FW(31)+FW(32)+FW(33)
790 +FW(34)+FW(35)+FW(36)+FW(37)+FW(38)+FW(39)+FW(40)+FW(41) ;
791 break ;
792
793 } /* end of switch on IRAD */
794
795 # undef FW
796 #endif /* UNROLL */
797
798 fjk[jj+IRAD1][kk+IRAD1] = sum ;
799 }
800 }
801
802 /*- y interpolations -*/
803
804 for( wsum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){
805 yw = fabsf(fy - qq) ; wt = sinc(yw) ;
806 yw /= WRAD ; if( yw > WCUT ) wt *= ww(yw) ;
807 wtt[qq+IRAD1] = wt ; wsum += wt ;
808 }
809 wfac *= wsum ;
810
811 for( kk=-IRAD1 ; kk <= IRAD ; kk++ ){
812 #ifndef UNROLL
813 for( sum=0.0f,jj=-IRAD1 ; jj <= IRAD ; jj++ ){
814 sum += wtt[jj+IRAD1]*fjk[jj+IRAD1][kk+IRAD1] ;
815 }
816 #else
817 for( sum=0.0f,jj=-IRAD1 ; jj < IRAD ; jj+=2 ){ /* unrolled by 2 */
818 sum += wtt[jj+IRAD1]*fjk[jj+IRAD1][kk+IRAD1]
819 +wtt[jj+IRAD ]*fjk[jj+IRAD ][kk+IRAD1] ;
820 }
821 #endif
822 fk[kk+IRAD1] = sum ;
823 }
824
825 /*- z interpolation -*/
826
827 for( wsum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){
828 zw = fabsf(fz - qq) ; wt = sinc(zw) ;
829 zw /= WRAD ; if( zw > WCUT ) wt *= ww(zw) ;
830 wtt[qq+IRAD1] = wt ; wsum += wt ;
831 }
832 wfac *= wsum ;
833
834 #ifndef UNROLL
835 for( sum=0.0f,kk=-IRAD1 ; kk <= IRAD ; kk++ ){
836 sum += wtt[kk+IRAD1] * fk[kk+IRAD1] ;
837 }
838 #else
839 for( sum=0.0f,kk=-IRAD1 ; kk < IRAD ; kk+=2 ){ /* unrolled by 2 */
840 sum += wtt[kk+IRAD1] * fk[kk+IRAD1]
841 +wtt[kk+IRAD ] * fk[kk+IRAD ] ;
842 }
843 #endif
844
845 vv[pp] = sum / wfac ;
846 }
847
848 } /* end OpenMP */
849 AFNI_OMP_END ;
850
851 EXRETURN ;
852 }
853
854 /*---------------------------------------------------------------------------*/
855 /* Master function for windowed sinc interpolation. */
856
GA_interp_wsinc5(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)857 void GA_interp_wsinc5( MRI_IMAGE *fim ,
858 int npp, float *ip, float *jp, float *kp, float *vv )
859 {
860 static int first = 1 ;
861 ENTRY("GA_interp_wsinc5") ;
862
863 if( first ){ setup_wsinc5() ; first = 0 ; }
864
865 if( WSHAP ) GA_interp_wsinc5s( fim,npp,ip,jp,kp,vv ) ; /* spherical */
866 else GA_interp_wsinc5p( fim,npp,ip,jp,kp,vv ) ; /* not spherical */
867
868 EXRETURN ;
869 }
870
871 /*---------------------------------------------------------------------------*/
872 /* Square (tensor product) windowed sinc interpolation in 2D. */
873
874 #undef FAR2D
875 #define FAR2D(i,j) far[(i)+(j)*nx]
876
877 #undef FARJ2D
878 #define FARJ2D(j,k) (far+(j)*nx)
879
GA_interp_wsinc5_2D(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * vv)880 void GA_interp_wsinc5_2D( MRI_IMAGE *fim ,
881 int npp, float *ip, float *jp, float *vv )
882 {
883 ENTRY("GA_interp_wsinc5_2D") ;
884
885 AFNI_OMP_START ;
886 #pragma omp parallel if(npp > 444)
887 {
888 int nx=fim->nx , ny=fim->ny , nxy=nx*ny , pp ;
889 float nxh=nx-0.501f , nyh=ny-0.501f , xx,yy ;
890 float fx,fy ;
891 float *far = MRI_FLOAT_PTR(fim) , *farj ;
892 int nx1=nx-1,ny1=ny-1, ix,jy ;
893
894 float xw,yw,rr , sum,wsum,wfac,wt ;
895 int iq,jq , qq,jj ;
896 float xsin[2*IRAD_MAX] , ysin[2*IRAD_MAX] ;
897 float wtt [2*IRAD_MAX] , fj [2*IRAD_MAX] ;
898 int iqq[2*IRAD_MAX] ;
899
900 /*----- loop over points -----*/
901
902 #pragma omp for
903 for( pp=0 ; pp < npp ; pp++ ){
904 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
905 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
906
907 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
908 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
909
910 if( ISTINY(fx) && ISTINY(fy) ){
911 CLIP(ix,nx1); CLIP(jy,ny1); vv[pp] = FAR2D(ix,jy); continue;
912 }
913
914 /*- x interpolations -*/
915
916 for( wsum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){ /* setup weights */
917 xw = fabsf(fx-qq) ; wt = sinc(xw) ;
918 xw /= WRAD ; if( xw > WCUT ) wt *= ww(xw) ;
919 wtt[qq+IRAD1] = wt ; wsum += wt ;
920 iq = ix+qq ; CLIP(iq,nx1) ; iqq[qq+IRAD1] = iq ;
921 }
922 wfac = wsum ;
923
924 for( jj=-IRAD1 ; jj <= IRAD ; jj++ ){
925 jq = jy+jj ; CLIP(jq,ny1) ;
926 for( sum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){
927 iq = iqq[qq+IRAD1] ; sum += FAR2D(iq,jq) * wtt[qq+IRAD1] ;
928 }
929 fj[jj+IRAD1] = sum ;
930 }
931
932 /*- y interpolations -*/
933
934 for( wsum=0.0f,qq=-IRAD1 ; qq <= IRAD ; qq++ ){
935 yw = fabsf(fy - qq) ; wt = sinc(yw) ;
936 yw /= WRAD ; if( yw > WCUT ) wt *= ww(yw) ;
937 wtt[qq+IRAD1] = wt ; wsum += wt ;
938 }
939 wfac *= wsum ;
940
941 for( sum=0.0f,jj=-IRAD1 ; jj <= IRAD ; jj++ ){
942 sum += wtt[jj+IRAD1]*fj[jj+IRAD1] ;
943 }
944 vv[pp] = sum / wfac ;
945 }
946
947 } /* end OpenMP */
948 AFNI_OMP_END ;
949
950 EXRETURN ;
951 }
952
953 /*===========================================================================*/
954 /* define quintic interpolation polynomials (Lagrange) */
955
956 #undef Q_M2
957 #undef Q_M1
958 #undef Q_00
959 #undef Q_P1
960 #undef Q_P2
961 #undef Q_P3
962 #define Q_M2(x) (x*(x*x-1.0)*(2.0-x)*(x-3.0)*0.008333333)
963 #define Q_M1(x) (x*(x*x-4.0)*(x-1.0)*(x-3.0)*0.041666667)
964 #define Q_00(x) ((x*x-4.0)*(x*x-1.0)*(3.0-x)*0.083333333)
965 #define Q_P1(x) (x*(x*x-4.0)*(x+1.0)*(x-3.0)*0.083333333)
966 #define Q_P2(x) (x*(x*x-1.0)*(x+2.0)*(3.0-x)*0.041666667)
967 #define Q_P3(x) (x*(x*x-1.0)*(x*x-4.0)*0.008333333)
968
969 /*--------------------------------------------------------------------*/
970 /*! Interpolate an image at npp (index) points, using quintic method. */
971
GA_interp_quintic(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,float * vv)972 void GA_interp_quintic( MRI_IMAGE *fim ,
973 int npp, float *ip, float *jp, float *kp, float *vv )
974 {
975 ENTRY("GA_interp_quintic") ;
976
977 RAND_ROUND ; /* 26 Feb 2020 */
978
979 AFNI_OMP_START ;
980 #pragma omp parallel if(npp > 1111)
981 {
982 int nx=fim->nx , ny=fim->ny , nz=fim->nz , nxy=nx*ny , pp ;
983 float nxh=nx-0.501f , nyh=ny-0.501f , nzh=nz-0.501f , xx,yy,zz ;
984 float fx,fy,fz ;
985 float *far = MRI_FLOAT_PTR(fim) ;
986 int nx1=nx-1,ny1=ny-1,nz1=nz-1, ix,jy,kz ;
987 int ix_m2,ix_m1,ix_00,ix_p1,ix_p2,ix_p3 ; /* interpolation indices */
988 int jy_m2,jy_m1,jy_00,jy_p1,jy_p2,jy_p3 ; /* (input image) */
989 int kz_m2,kz_m1,kz_00,kz_p1,kz_p2,kz_p3 ;
990
991 float wt_m2,wt_m1,wt_00,wt_p1,wt_p2,wt_p3 ; /* interpolation weights */
992
993 float f_jm2_km2, f_jm1_km2, f_j00_km2, f_jp1_km2, f_jp2_km2, f_jp3_km2,
994 f_jm2_km1, f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, f_jp3_km1,
995 f_jm2_k00, f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00, f_jp3_k00,
996 f_jm2_kp1, f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1, f_jp3_kp1,
997 f_jm2_kp2, f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2, f_jp3_kp2,
998 f_jm2_kp3, f_jm1_kp3, f_j00_kp3, f_jp1_kp3, f_jp2_kp3, f_jp3_kp3,
999 f_km2 , f_km1 , f_k00 , f_kp1 , f_kp2 , f_kp3 ;
1000
1001 #pragma omp for
1002 for( pp=0 ; pp < npp ; pp++ ){
1003 xx = ip[pp] ; if( xx < -0.499f || xx > nxh ){ vv[pp]=outval; continue; }
1004 yy = jp[pp] ; if( yy < -0.499f || yy > nyh ){ vv[pp]=outval; continue; }
1005 zz = kp[pp] ; if( zz < -0.499f || zz > nzh ){ vv[pp]=outval; continue; }
1006
1007 ix = floorf(xx) ; fx = xx - ix ; /* integer and */
1008 jy = floorf(yy) ; fy = yy - jy ; /* fractional coords */
1009 kz = floorf(zz) ; fz = zz - kz ;
1010
1011 if( ISTINY(fx) && ISTINY(fy) && ISTINY(fz) ){
1012 CLIP(ix,nx1); CLIP(jy,ny1); CLIP(kz,nz1);
1013 vv[pp] = FAR(ix,jy,kz) ; continue ;
1014 }
1015
1016 /* compute indexes from which to interpolate (-2,-1,0,+1,+2,+3),
1017 but clipped to lie inside input image volume */
1018
1019 ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ;
1020 CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
1021 ix_m2 = ix-2 ; ix_p3 = ix+3 ;
1022 CLIP(ix_m2,nx1) ; CLIP(ix_p3,nx1) ;
1023
1024 jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ;
1025 CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
1026 jy_m2 = jy-2 ; jy_p3 = jy+3 ;
1027 CLIP(jy_m2,ny1) ; CLIP(jy_p3,ny1) ;
1028
1029 kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ;
1030 CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
1031 kz_m2 = kz-2 ; kz_p3 = kz+3 ;
1032 CLIP(kz_m2,nz1) ; CLIP(kz_p3,nz1) ;
1033
1034 wt_m1 = Q_M1(fx) ; wt_00 = Q_00(fx) ; /* interpolation weights */
1035 wt_p1 = Q_P1(fx) ; wt_p2 = Q_P2(fx) ; /* in x-direction */
1036 wt_m2 = Q_M2(fx) ; wt_p3 = Q_P3(fx) ;
1037
1038 #undef XINT
1039 #define XINT(j,k) wt_m2*FAR(ix_m2,j,k)+wt_m1*FAR(ix_m1,j,k) \
1040 +wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k) \
1041 +wt_p2*FAR(ix_p2,j,k)+wt_p3*FAR(ix_p3,j,k)
1042
1043 /* interpolate to location ix+fx at each jy,kz level */
1044
1045 f_jm2_km2 = XINT(jy_m2,kz_m2) ; f_jm1_km2 = XINT(jy_m1,kz_m2) ;
1046 f_j00_km2 = XINT(jy_00,kz_m2) ; f_jp1_km2 = XINT(jy_p1,kz_m2) ;
1047 f_jp2_km2 = XINT(jy_p2,kz_m2) ; f_jp3_km2 = XINT(jy_p3,kz_m2) ;
1048
1049 f_jm2_km1 = XINT(jy_m2,kz_m1) ; f_jm1_km1 = XINT(jy_m1,kz_m1) ;
1050 f_j00_km1 = XINT(jy_00,kz_m1) ; f_jp1_km1 = XINT(jy_p1,kz_m1) ;
1051 f_jp2_km1 = XINT(jy_p2,kz_m1) ; f_jp3_km1 = XINT(jy_p3,kz_m1) ;
1052
1053 f_jm2_k00 = XINT(jy_m2,kz_00) ; f_jm1_k00 = XINT(jy_m1,kz_00) ;
1054 f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
1055 f_jp2_k00 = XINT(jy_p2,kz_00) ; f_jp3_k00 = XINT(jy_p3,kz_00) ;
1056
1057 f_jm2_kp1 = XINT(jy_m2,kz_p1) ; f_jm1_kp1 = XINT(jy_m1,kz_p1) ;
1058 f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
1059 f_jp2_kp1 = XINT(jy_p2,kz_p1) ; f_jp3_kp1 = XINT(jy_p3,kz_p1) ;
1060
1061 f_jm2_kp2 = XINT(jy_m2,kz_p2) ; f_jm1_kp2 = XINT(jy_m1,kz_p2) ;
1062 f_j00_kp2 = XINT(jy_00,kz_p2) ; f_jp1_kp2 = XINT(jy_p1,kz_p2) ;
1063 f_jp2_kp2 = XINT(jy_p2,kz_p2) ; f_jp3_kp2 = XINT(jy_p3,kz_p2) ;
1064
1065 f_jm2_kp3 = XINT(jy_m2,kz_p3) ; f_jm1_kp3 = XINT(jy_m1,kz_p3) ;
1066 f_j00_kp3 = XINT(jy_00,kz_p3) ; f_jp1_kp3 = XINT(jy_p1,kz_p3) ;
1067 f_jp2_kp3 = XINT(jy_p2,kz_p3) ; f_jp3_kp3 = XINT(jy_p3,kz_p3) ;
1068
1069 /* interpolate to jy+fy at each kz level */
1070
1071 wt_m1 = Q_M1(fy) ; wt_00 = Q_00(fy) ; wt_p1 = Q_P1(fy) ;
1072 wt_p2 = Q_P2(fy) ; wt_m2 = Q_M2(fy) ; wt_p3 = Q_P3(fy) ;
1073
1074 f_km2 = wt_m2 * f_jm2_km2 + wt_m1 * f_jm1_km2 + wt_00 * f_j00_km2
1075 + wt_p1 * f_jp1_km2 + wt_p2 * f_jp2_km2 + wt_p3 * f_jp3_km2 ;
1076
1077 f_km1 = wt_m2 * f_jm2_km1 + wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
1078 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 + wt_p3 * f_jp3_km1 ;
1079
1080 f_k00 = wt_m2 * f_jm2_k00 + wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
1081 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 + wt_p3 * f_jp3_k00 ;
1082
1083 f_kp1 = wt_m2 * f_jm2_kp1 + wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
1084 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 + wt_p3 * f_jp3_kp1 ;
1085
1086 f_kp2 = wt_m2 * f_jm2_kp2 + wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
1087 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 + wt_p3 * f_jp3_kp2 ;
1088
1089 f_kp3 = wt_m2 * f_jm2_kp3 + wt_m1 * f_jm1_kp3 + wt_00 * f_j00_kp3
1090 + wt_p1 * f_jp1_kp3 + wt_p2 * f_jp2_kp3 + wt_p3 * f_jp3_kp3 ;
1091
1092 /* interpolate to kz+fz to get output */
1093
1094 wt_m1 = Q_M1(fz) ; wt_00 = Q_00(fz) ; wt_p1 = Q_P1(fz) ;
1095 wt_p2 = Q_P2(fz) ; wt_m2 = Q_M2(fz) ; wt_p3 = Q_P3(fz) ;
1096
1097 vv[pp] = wt_m2 * f_km2 + wt_m1 * f_km1 + wt_00 * f_k00
1098 + wt_p1 * f_kp1 + wt_p2 * f_kp2 + wt_p3 * f_kp3 ;
1099 }
1100
1101 } /* end OpenMP */
1102 AFNI_OMP_END ;
1103
1104 EXRETURN ;
1105 }
1106
1107 /*===========================================================================*/
1108 /*--- Stuff for storing sub-BLOKs of data points for localized cost funcs ---*/
1109 /*--- such as lpc and lpa ---*/
1110 /*---------------------------------------------------------------------------*/
1111
1112 /* FAS Test: is abs(a) <= s ?? */
1113
1114 #undef FAS
1115 #define FAS(a,s) ( (a) <= (s) && (a) >= -(s) )
1116
1117 /**** Macros to test if a point is inside a given BLOK type ****/
1118
1119 /** define inside of a ball; is point (a,b,c) inside? **/
1120 /** volume of ball = 4*PI/3 * rad**3 = 4.1888 * rad**3 **/
1121 /** In this test below, siz = rad^2 **/
1122
1123 #define GA_BLOK_inside_ball(a,b,c,siz) \
1124 ( ((a)*(a)+(b)*(b)+(c)*(c)) <= (siz) )
1125
1126 /** define inside of a cube **/
1127 /** volume of cube = 8 * siz**3, with siz = rad **/
1128 /** lattice vectors = [2*siz,0,0] [0,2*siz,0] [0,0,2*siz] **/
1129
1130 #define GA_BLOK_inside_cube(a,b,c,siz) \
1131 ( FAS((a),(siz)) && FAS((b),(siz)) && FAS((c),(siz)) )
1132
1133 /** define inside of a rhombic dodecahedron (RHDD) **/
1134 /** volume of RHDD = 2 * siz**3, with siz = rad **/
1135 /** lattice vectors = [siz,siz,0] [0,siz,siz] [siz,0,siz] **/
1136
1137 #define GA_BLOK_inside_rhdd(a,b,c,siz) \
1138 ( FAS((a)+(b),(siz)) && FAS((a)-(b),(siz)) && \
1139 FAS((a)+(c),(siz)) && FAS((a)-(c),(siz)) && \
1140 FAS((b)+(c),(siz)) && FAS((b)-(c),(siz)) )
1141
1142 /** define inside of a truncated octahedron (TOHD) **/
1143 /** volume of TOHD = 4 * siz**, with siz = rad3 **/
1144 /** lattice vectors = [-siz,siz,siz] [siz,-siz,siz] [siz,siz,-siz] **/
1145
1146 #define GA_BLOK_inside_tohd(a,b,c,siz) \
1147 ( FAS((a),(siz)) && FAS((b),(siz)) && FAS((c),(siz)) && \
1148 FAS((a)+(b)+(c),1.5f*(siz)) && FAS((a)-(b)+(c),1.5f*(siz)) && \
1149 FAS((a)+(b)-(c),1.5f*(siz)) && FAS((a)-(b)-(c),1.5f*(siz)) )
1150
1151 /** define inside of an arbitrary blok type **/
1152
1153 #define GA_BLOK_inside(bt,a,b,c,s) \
1154 ( ((bt)==GA_BLOK_BALL) ? GA_BLOK_inside_ball((a),(b),(c),(s)) \
1155 : ((bt)==GA_BLOK_CUBE) ? GA_BLOK_inside_cube((a),(b),(c),(s)) \
1156 : ((bt)==GA_BLOK_RHDD) ? GA_BLOK_inside_rhdd((a),(b),(c),(s)) \
1157 : ((bt)==GA_BLOK_TOHD) ? GA_BLOK_inside_tohd((a),(b),(c),(s)) \
1158 : 0 )
1159
1160 /** add 1 value to a dynamically allocated integer array;
1161 note that we expand it by 50% if we have overflowed current space **/
1162
1163 #define GA_BLOK_ADDTO_intar(nar,nal,ar,val) \
1164 do{ if( (nar) == (nal) ){ \
1165 (nal) = 1.5*(nal)+16; (ar) = (int *)realloc((ar),sizeof(int)*(nal)); \
1166 } \
1167 (ar)[(nar)++] = (val); \
1168 } while(0)
1169
1170 /** truncate dynamically allocated integer array down to final size **/
1171
1172 #define GA_BLOK_CLIP_intar(nar,nal,ar) \
1173 do{ if( (nar) < (nal) && (nar) > 0 ){ \
1174 (nal) = (nar); (ar) = (int *)realloc((ar),sizeof(int)*(nal)); \
1175 }} while(0)
1176
1177 /*----------------------------------------------------------------------------*/
1178 /*! Fill a struct with list of points contained in sub-bloks of the base.
1179
1180 - nx,ny,nz = 3D grid dimensions
1181 - dx,dy,dz = 3D grid spacings
1182 - npt = number of points stored in im,jm,km
1183 - im,jm,km = 3D indexes of points to blok-ize
1184 (can be NULL, in which case all nx*ny*nz points are used)
1185 - bloktype = one of GA_BLOK_BALL, GA_BLOK_CUBE, GA_BLOK_RHDD, GA_BLOK_TOHD
1186 - blokrad = radius parameter for the bloks to be built
1187 - minel = minimum number of points to put in a blok
1188 (if 0, function will pick a value)
1189 - shfac = shrinkage factor -- normally, bloks don't overlap much, but
1190 you can specify shfac to change the lattice size:
1191 < 1 makes them overlap more, and > 1 makes them spaced apart
1192 - verb = whether to print out some verbosity stuff FYI
1193
1194 I admit this code is tricksy. I'll try to comment it better now [Jan 2021].
1195 *//*--------------------------------------------------------------------------*/
1196
create_GA_BLOK_set(int nx,int ny,int nz,float dx,float dy,float dz,int npt,float * im,float * jm,float * km,int bloktype,float blokrad,int minel,float shfac,int verb)1197 GA_BLOK_set * create_GA_BLOK_set( int nx , int ny , int nz ,
1198 float dx , float dy , float dz ,
1199 int npt , float *im, float *jm, float *km,
1200 int bloktype , float blokrad , int minel ,
1201 float shfac , int verb )
1202 {
1203 GA_BLOK_set *gbs ;
1204 float dxp,dyp,dzp , dxq,dyq,dzq , dxr,dyr,dzr , xt,yt,zt ;
1205 float xx,yy,zz , uu,vv,ww , siz ;
1206 THD_mat33 latmat , invlatmat ; THD_fvec3 pqr , xyz ;
1207 int pb,pt , qb,qt , rb,rt , pp,qq,rr , nblok,nball , ii , nxy ;
1208 int aa,bb,cc , dd,ss , np,nq,nr,npq , *nelm,*nalm,**elm , ntot,nsav,ndup ;
1209
1210 ENTRY("create_GA_BLOK_set") ;
1211
1212 /* check for badness */
1213
1214 if( nx < 3 || ny < 3 || nz < 1 ) RETURN(NULL) ; /* nonsense */
1215 if( dx <= 0.0f ) dx = 1.0f ; /* fixes */
1216 if( dy <= 0.0f ) dy = 1.0f ;
1217 if( dz <= 0.0f ) dz = 1.0f ;
1218
1219 /* correct for crazy input shrinkage factor */
1220
1221 if( shfac < 0.2f || shfac > 5.0f ) shfac = 1.0f ;
1222
1223 /* Normal operation specifies only a subset of full 3D space
1224 to be used for matching and so to be included in the BLOKs.
1225 If code below is used, then all of 3D space will be used instead. */
1226
1227 if( npt <= 0 || im == NULL || jm == NULL || km == NULL ){
1228 im = jm = km = NULL ; npt = 0 ;
1229 }
1230
1231 /*------------ Mark type of blok being stored ------------*/
1232
1233 /* Create lattice vectors (dxp,dyp,dzp) etc, to generate translated bloks:
1234 The (p,q,r)-th blok - for integral p,q,r - is centered at (x,y,z) offset
1235 (dxp,dyp,dzp)*p + (dxq,dyq,dzq)*q + (dxr,dyr,dzr)*r
1236 Also set the 'siz' parameter for the blok, to test for inclusion;
1237 what 'siz' is with respect to blokrad depends on the blok shape in use */
1238
1239 switch( bloktype ){
1240
1241 /* balls go on a hexagonal close packed lattice,
1242 but with lattice spacing reduced to avoid gaps
1243 (of course, then the balls overlap -- c'est la geometrie) */
1244
1245 case GA_BLOK_BALL:{
1246 float s3=1.73205f , /* sqrt(3) */
1247 s6=2.44949f , /* sqrt(6) */
1248 a =blokrad*0.866025f ; /* shrink spacing to avoid gaps */
1249 siz = blokrad*blokrad ;
1250 /* hexagonal close packing basis vectors for sphere of radius a */
1251 a *= shfac ;
1252 dxp = 2.0f * a ; dyp = 0.0f ; dzp = 0.0f ;
1253 dxq = a ; dyq = a * s3; dzq = 0.0f ;
1254 dxr = a ; dyr = a / s3; dzr = a * 0.666667f*s6 ;
1255 }
1256 break ;
1257
1258 /* cubes go on a simple cubical lattice, spaced so faces touch */
1259
1260 case GA_BLOK_CUBE:{
1261 float a = blokrad ;
1262 siz = a ; a *= shfac ;
1263 dxp = 2*a ; dyp = 0.0f; dzp = 0.0f ;
1264 dxq = 0.0f; dyq = 2*a ; dzq = 0.0f ;
1265 dxr = 0.0f; dyr = 0.0f; dzr = 2*a ;
1266 }
1267 break ;
1268
1269 /* rhombic dodecahedra go on a FCC lattice,
1270 spaced so that faces touch (i.e., no volumetric overlap) */
1271
1272 case GA_BLOK_RHDD:{
1273 float a = blokrad ;
1274 siz = a ; a *= shfac ;
1275 dxp = a ; dyp = a ; dzp = 0.0f ;
1276 dxq = 0.0f; dyq = a ; dzq = a ;
1277 dxr = a ; dyr = 0.0f; dzr = a ;
1278 }
1279 break ;
1280
1281 /* truncated octahedra go on a BCC lattice,
1282 spaced so that faces touch (i.e., no volumetric overlap) */
1283
1284 case GA_BLOK_TOHD:{
1285 float a = blokrad ;
1286 siz = a ; a *= shfac ;
1287 dxp = -a ; dyp = a ; dzp = a ;
1288 dxq = a ; dyq = -a ; dzq = a ;
1289 dxr = a ; dyr = a ; dzr = -a ;
1290 }
1291 break ;
1292
1293 default: RETURN(NULL) ; /** should not happen == stupid caller! **/
1294 }
1295
1296 /*** find range of (p,q,r) indexes needed to cover volume,
1297 by checking out all 7 corners besides (0,0,0) (where p=q=r=0) ***/
1298
1299 /* latmat = lattice matrix
1300 = 3x3 matrix that takes (p,q,r) to blok center (x,y,z)
1301 = set from d[xyz][pqr] values set above for each blok type.
1302 So invlatmat is a 3x3 matrix that takes (x,y,z) to (p,q,r).
1303 We apply invlatmat to each corner of the 3D volume, and
1304 find the smallest and largest (p,q,r) values that happen.
1305 Those min/max values give the range of (p,q,r)
1306 that has to be scanned for entry into the blokset. */
1307
1308 LOAD_MAT( latmat, dxp , dxq , dxr ,
1309 dyp , dyq , dyr ,
1310 dzp , dzq , dzr ) ; invlatmat = MAT_INV(latmat) ;
1311
1312 /*----- corner 0 -----*/
1313 xt = (nx-1)*dx ; yt = (ny-1)*dy ; zt = (nz-1)*dz ; /* xyz top values */
1314 pb = pt = qb = qt = rb = rt = 0 ; /* initialize (p,q,r) bot, top values */
1315
1316 /*----- corner 1 -----*/
1317 LOAD_FVEC3(xyz , xt,0.0f,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
1318 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1319 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1320 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1321
1322 /*----- corner 2 -----*/
1323 LOAD_FVEC3(xyz , xt,yt,0.0f ) ; pqr = MATVEC( invlatmat , xyz ) ;
1324 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1325 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1326 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1327
1328 /*----- corner 3 -----*/
1329 LOAD_FVEC3(xyz , xt,0.0f,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
1330 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1331 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1332 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1333
1334 /*----- corner 4 -----*/
1335 LOAD_FVEC3(xyz , xt,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
1336 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1337 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1338 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1339
1340 /*----- corner 5 -----*/
1341 LOAD_FVEC3(xyz , 0.0f,yt,0.0f ); pqr = MATVEC( invlatmat , xyz ) ;
1342 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1343 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1344 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1345
1346 /*----- corner 6 -----*/
1347 LOAD_FVEC3(xyz , 0.0f,0.0f,zt ); pqr = MATVEC( invlatmat , xyz ) ;
1348 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1349 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1350 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1351
1352 /*----- corner 7 -----*/
1353 LOAD_FVEC3(xyz , 0.0f,yt,zt ) ; pqr = MATVEC( invlatmat , xyz ) ;
1354 pp = (int)floorf( pqr.xyz[0] ) ; pb = MIN(pb,pp) ; pp++ ; pt = MAX(pt,pp) ;
1355 qq = (int)floorf( pqr.xyz[1] ) ; qb = MIN(qb,qq) ; qq++ ; qt = MAX(qt,qq) ;
1356 rr = (int)floorf( pqr.xyz[2] ) ; rb = MIN(rb,rr) ; rr++ ; rt = MAX(rt,rr) ;
1357
1358 /*** Lattice index range is (p,q,r) = (pb..pt,qb..qt,rb..rt) inclusive ***/
1359
1360 np = pt-pb+1 ; /* number of p values to consider */
1361 nq = qt-qb+1 ; npq = np*nq ;
1362 nr = rt-rb+1 ;
1363 nblok = npq*nr ; /* total number of bloks to consider */
1364 /* not all of them will actually get used */
1365
1366 /*----- Now have list of bloks, so put points into each blok list -----*/
1367
1368 /* create empty lists of points (calloc == all contents are zero/NULL) */
1369
1370 nelm = (int *) calloc(sizeof(int) ,nblok) ; /* # pts in each blok */
1371 nalm = (int *) calloc(sizeof(int) ,nblok) ; /* # malloc-ed in each blok */
1372 elm = (int **)calloc(sizeof(int *),nblok) ; /* list of pts in each blok */
1373
1374 nxy = nx*ny ; if( npt == 0 ) npt = nxy*nz ; /* npt = all of 3D space?? */
1375
1376 /*--- loop over dataset control points to put them into bloks ---*/
1377
1378 for( ndup=ntot=ii=0 ; ii < npt ; ii++ ){
1379
1380 if( im != NULL ){ /* 3D index of point #ii comes from arrays */
1381
1382 pp = (int)im[ii]; qq = (int)jm[ii]; rr = (int)km[ii]; /* xyz indexes */
1383
1384 } else { /* 3D index is computed from ii, as we use all points */
1385
1386 pp = ii%nx ; rr = ii/nxy ; qq = (ii-rr*nxy)/nx ;
1387
1388 }
1389
1390 ss = ii ; /* index in 1D array of points to be matched */
1391 /* this is NOT the index in the dataset, unless im==NULL */
1392 /* it is the index into the point list given by (im,jm,km) */
1393 /* so the index in the 3D dataset is pp+qq*nx+rr*nxy */
1394
1395 xx = pp*dx ; yy = qq*dy ; zz = rr*dz ; /* xyz spatial coordinates inside */
1396 LOAD_FVEC3( xyz , xx,yy,zz ) ; /* grid of index (pp,qq,rr) */
1397
1398 pqr = MATVEC( invlatmat , xyz ) ; /* float lattice coordinates */
1399 pp = (int)floorf(pqr.xyz[0]+.499f) ; /* round to integer lattice coords */
1400 qq = (int)floorf(pqr.xyz[1]+.499f) ; /* [note pp,qq,rr have changed from] */
1401 rr = (int)floorf(pqr.xyz[2]+.499f) ; /* [grid indexes to lattice indexes] */
1402
1403 nsav = 0 ; /* nsav = num times point #ii at (xx,yy,zz) is saved to a blok */
1404
1405 /* searching bloks surrounding (pp,qq,rr) for inclusion of (xx,yy,zz) */
1406 /* a point MIGHT be in more than one blok (not common if shfac >= 1) */
1407
1408 for( cc=rr-1 ; cc <= rr+1 ; cc++ ){ /* cc is near rr */
1409 if( cc < rb || cc > rt ) continue ; /* out of range */
1410 for( bb=qq-1 ; bb <= qq+1 ; bb++ ){ /* bb is near qq */
1411 if( bb < qb || bb > qt ) continue ;
1412 for( aa=pp-1 ; aa <= pp+1 ; aa++ ){ /* aa is near pp */
1413 if( aa < pb || aa > pt ) continue ;
1414
1415 LOAD_FVEC3( pqr , aa,bb,cc ) ; /* compute center of this */
1416 xyz = MATVEC( latmat , pqr ) ; /* blok into xyz vector */
1417 uu = xx - xyz.xyz[0] ; /* xyz coords of point #ii (xx,yy,zz) */
1418 vv = yy - xyz.xyz[1] ; /* relative to blok center (xyz vector) */
1419 ww = zz - xyz.xyz[2] ;
1420 if( GA_BLOK_inside( bloktype , uu,vv,ww , siz ) ){ /* add to blok */
1421 dd = (aa-pb) + (bb-qb)*np + (cc-rb)*npq ; /* blok index */
1422 GA_BLOK_ADDTO_intar( nelm[dd], nalm[dd], elm[dd], ss ) ;
1423 ntot++ ; /* total number of points saved into bloks */
1424 nsav++ ; /* how many times THIS point was saved into a blok */
1425 }
1426 }
1427 }
1428 }
1429 if( nsav > 1 ) ndup++ ; /* count of points that had duplicate saves */
1430 }
1431
1432 /**----- compute minimum number of points saved per blok? -----**/
1433
1434 if( minel < 9 ){
1435 for( minel=dd=0 ; dd < nblok ; dd++ ) minel = MAX(minel,nelm[dd]) ;
1436 minel = (int)(0.456*minel)+1 ;
1437 }
1438
1439 /**----- now cast out bloks that have too few points,
1440 and truncate those arrays that pass the threshold;
1441 nsav will now be the number of saved bloks -----**/
1442
1443 for( nsav=dd=0 ; dd < nblok ; dd++ ){
1444 if( nelm[dd] < minel ){ /* too empty == destroy blok */
1445 if( elm[dd] != NULL ){ free(elm[dd]); elm[dd] = NULL; }
1446 nelm[dd] = 0 ;
1447 } else { /* clip point list back to used array size */
1448 GA_BLOK_CLIP_intar( nelm[dd] , nalm[dd] , elm[dd] ) ; nsav++ ;
1449 }
1450 }
1451 free(nalm) ; /* counts of allocated space per blok -- no longer needed */
1452
1453 if( nsav == 0 ){ /* didn't find any blok arrays to keep!? */
1454 ERROR_message("create_GA_BLOK_set can't get bloks with at least %d elements",minel);
1455 free(nelm) ; free(elm) ; RETURN(NULL) ;
1456 }
1457
1458 /*----- create output struct from all of the above -----*/
1459
1460 gbs = (GA_BLOK_set *)malloc(sizeof(GA_BLOK_set)) ;
1461 gbs->num = nsav ;
1462 gbs->nelm = (int *) calloc(sizeof(int) ,nsav) ;
1463 gbs->elm = (int **)calloc(sizeof(int *),nsav) ;
1464 for( ntot=nsav=dd=0 ; dd < nblok ; dd++ ){
1465 if( nelm[dd] > 0 && elm[dd] != NULL ){
1466 gbs->nelm[nsav] = nelm[dd] ; ntot += nelm[dd] ;
1467 gbs->elm [nsav] = elm[dd] ; nsav++ ;
1468 }
1469 }
1470 free(nelm) ; free(elm) ;
1471
1472 /*--- set some other parameters in the blok set ---*/
1473
1474 gbs->ppow = AFNI_numenv("AFNI_LPC_POWER") ; /* 28 Jan 2021 */
1475 if( gbs->ppow <= 0.0f || gbs->ppow > 2.0f ) gbs->ppow = 1.0f ;
1476
1477 gbs->nx = nx ; gbs->ny = ny ; gbs->nz = nz ; /* Biden Day 3 */
1478 gbs->dx = dx ; gbs->dy = dy ; gbs->dz = dz ;
1479
1480 if( verb )
1481 ININFO_message("%d total points stored in %d '%s(%g)' bloks (%d duplicates)",
1482 ntot , gbs->num , GA_BLOK_STRING(bloktype) , blokrad , ndup ) ;
1483
1484 RETURN(gbs) ;
1485 }
1486
1487 /*---------------------------------------------------------------------------*/
1488 /*! Return the vector of individual blok correlations, for further analysis.
1489 Each value in the returned vector corresponds to one blok in gbs. Note
1490 that the number of points in the 'vm' arrays isn't needed as an input,
1491 since that is implicitly encoded in gbs. This is for Paul Taylor.
1492 *//*-------------------------------------------------------------------------*/
1493
GA_pearson_vector(GA_BLOK_set * gbs,float * avm,float * bvm,float * wvm)1494 floatvec * GA_pearson_vector( GA_BLOK_set *gbs ,
1495 float *avm, float *bvm, float *wvm )
1496 {
1497 int nblok , nelm , *elm , dd , ii,jj , nm ;
1498 float xv,yv,xy,xm,ym,vv,ww,ws , wt ;
1499 floatvec *pv=NULL ; float *pvar ;
1500
1501 if( gbs == NULL || avm == NULL || bvm == NULL ) return NULL ;
1502
1503 nblok = gbs->num ; if( nblok < 1 ) return NULL ;
1504
1505 MAKE_floatvec( pv , nblok ) ; pvar = pv->ar ;
1506
1507 /* loop over bloks */
1508
1509 for( dd=0 ; dd < nblok ; dd++ ){
1510 pvar[dd] = 0.0f ; /* default */
1511 nelm = gbs->nelm[dd] ; if( nelm < 9 ) continue ; /* skip it */
1512 elm = gbs->elm[dd] ; /* array of indexes in avm (etc.) to use */
1513
1514 RAND_ROUND ; /* 26 Feb 2020 */
1515
1516 if( wvm == NULL ){ /*** unweighted correlation ***/
1517 xv=yv=xy=xm=ym=0.0f ;
1518 for( ii=0 ; ii < nelm ; ii++ ){ /* compute means */
1519 jj = elm[ii] ;
1520 xm += avm[jj] ; ym += bvm[jj] ;
1521 }
1522 xm /= nelm ; ym /= nelm ;
1523 for( ii=0 ; ii < nelm ; ii++ ){ /* compute (co)variances */
1524 jj = elm[ii] ;
1525 vv = avm[jj]-xm ; ww = bvm[jj]-ym ;
1526 xv += vv*vv ; yv += ww*ww ; xy += vv*ww ;
1527 }
1528
1529 } else { /*** weighted correlation ***/
1530 xv=yv=xy=xm=ym=ws=0.0f ;
1531 for( ii=0 ; ii < nelm ; ii++ ){ /* compute weighted means */
1532 jj = elm[ii] ;
1533 wt = wvm[jj] ; ws += wt ;
1534 xm += avm[jj]*wt ; ym += bvm[jj]*wt ;
1535 }
1536 xm /= ws ; ym /= ws ;
1537 for( ii=0 ; ii < nelm ; ii++ ){ /* compute weighted (co)variances */
1538 jj = elm[ii] ;
1539 wt = wvm[jj] ; vv = avm[jj]-xm ; ww = bvm[jj]-ym ;
1540 xv += wt*vv*vv ; yv += wt*ww*ww ; xy += wt*vv*ww ;
1541 }
1542 }
1543
1544 if( xv > 0.0f && yv > 0.0f ) pvar[dd] = xy/sqrtf(xv*yv) ; /* correlation */
1545 }
1546
1547 return pv ;
1548 }
1549
1550 /*---------------------------------------------------------------------------*/
1551 /*! Return a volume from the pearson vector showing with the
1552 voxel populated with 'their' blok's pearson correlation. [Biden day 3]
1553 This is for Paul Taylor.
1554 *//*-------------------------------------------------------------------------*/
1555
GA_pearson_image(GA_setup * stup,floatvec * pv)1556 MRI_IMAGE * GA_pearson_image( GA_setup *stup , floatvec *pv )
1557 {
1558 MRI_IMAGE *gim ; float *gar , *pvar ;
1559 GA_BLOK_set *gbs ;
1560 float *ima=NULL , *jma=NULL , *kma=NULL ;
1561 int nx,ny,nz,nxy,nxyz,dd,nelm,ii,jj , pp,qq,rr,npt ;
1562 int *elm ;
1563
1564 if( stup == NULL || pv == NULL ) return NULL ; /* dumdum user */
1565
1566 gbs = stup->blokset ; if( gbs == NULL ) return NULL ;
1567
1568 nx = gbs->nx ; ny = gbs->ny ; nz = gbs->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
1569 if( nx < 2 || ny < 2 || nz < 2 || nxyz < 99 ) return NULL ;
1570
1571 /* point indexes available? */
1572
1573 if( stup->im != NULL ) ima = stup->im->ar ;
1574 if( stup->jm != NULL ) jma = stup->jm->ar ;
1575 if( stup->km != NULL ) kma = stup->km->ar ;
1576 npt = stup->npt_match ;
1577 if( ima == NULL || jma == NULL || kma == NULL ){ /* no point index arrays */
1578 ima = jma = kma = NULL ; npt = nxyz ; /* so use all of 3D space */
1579 }
1580
1581 /* create output 3D volume */
1582
1583 gim = mri_new_vol(nx,ny,nz,MRI_float) ; /* zero filled */
1584 if( gim == NULL ) return NULL ; /* should never happen */
1585 gim->dx = gbs->dx ; gim->dy = gbs->dy ; gim->dz = gbs->dz ;
1586 gar = MRI_FLOAT_PTR(gim) ;
1587 pvar = pv->ar ;
1588
1589 /* loop over bloks */
1590
1591 /** INFO_message("GA_pearson_image: npt=%d ima=%s",npt,(ima==NULL)?"NULL":"non-NULL") ; **/
1592 for( dd=0 ; dd < gbs->num ; dd++ ){
1593 nelm = gbs->nelm[dd] ;
1594 elm = gbs->elm[dd] ; /* array of indexes for this blok */
1595 /** ININFO_message(" blok[%d] has %d points with value %g",dd,nelm,pvar[dd]) ; **/
1596 for( ii=0 ; ii < nelm ; ii++ ){ /* push values to image */
1597 /* jj = index into 1D matching arrays (ima,jma,kma) */
1598 jj = elm[ii] ;
1599 /* but 1D matching arrays are not usually all of 3D space! */
1600 if( ima != NULL ){
1601 pp = (int)ima[jj]; qq = (int)jma[jj]; rr = (int)kma[jj]; /* xyz indexes */
1602 jj = pp + qq*nx + rr*nxy ; /* now jj = index into 3D space */
1603 /** ININFO_message(" [%d] is at (%d,%d,%d) = %d",ii,pp,qq,rr,jj) ; **/
1604 }
1605 if( jj >=0 && jj < nxyz ) gar[jj] = pvar[dd] ; /* store into 3D volume */
1606 }
1607 }
1608
1609 return gim ;
1610 }
1611
1612 /*======================== End of BLOK-iness functionality ==================*/
1613
1614 /*---------------------------------------------------------------------------*/
1615
GA_indexwarp(MRI_IMAGE * inim,int interp_code,MRI_IMAGE * wpim)1616 MRI_IMAGE * GA_indexwarp( MRI_IMAGE *inim, int interp_code, MRI_IMAGE *wpim )
1617 {
1618 float_triple delta = {0.0f,0.0f,0.0f} ;
1619
1620 return GA_indexwarp_plus( inim , interp_code , wpim , delta , NULL ) ;
1621 }
1622
1623 /*---------------------------------------------------------------------------*/
1624 /* Output image will be on the same grid as the input, of course.
1625 Points not in the mask will be set to zero.
1626
1627 Input image format ==> output format:
1628 ------------------ -------------
1629 fvect ==> fvect
1630 rgb ==> rgb
1631 complex ==> complex
1632 other ==> float
1633 *//*-------------------------------------------------------------------------*/
1634
GA_indexwarp_plus(MRI_IMAGE * inim,int interp_code,MRI_IMAGE * wpim,float_triple delta,byte * mask)1635 MRI_IMAGE * GA_indexwarp_plus( MRI_IMAGE *inim, int interp_code,
1636 MRI_IMAGE *wpim, float_triple delta, byte *mask )
1637 {
1638 MRI_IMAGE *outim=NULL , *fim , *iim,*jjm,*kkm ;
1639 float *outar , *far , *iar,*jar,*kar ;
1640 float ffmin=0.0f , ffmax=0.0f , delx,dely,delz ;
1641 int ii,jj , nvox,nx,ny,nz,nxy , do_clip ;
1642 byte *mmm=NULL ; int nmask=0 ; float *miar,*mjar,*mkar,*moutar ;
1643
1644 ENTRY("GA_indexwarp_plus") ;
1645
1646 if( inim == NULL || wpim == NULL || wpim->kind != MRI_fvect ) RETURN(NULL);
1647 if( mri_data_pointer(inim) == NULL ||
1648 mri_data_pointer(wpim) == NULL || wpim->vdim != 3 ) RETURN(NULL);
1649 if( inim->nx != wpim->nx ||
1650 inim->ny != wpim->ny || inim->nz != wpim->nz ) RETURN(NULL);
1651
1652 /*- if input is itself a vector, use recursion to process each sub-image -*/
1653 /*- (for usage of VECTORME macro, via CALLME, see mrilib.h) -*/
1654
1655 #undef CALLME
1656 #define CALLME(inee,outee) \
1657 outee = GA_indexwarp_plus( (inee),interp_code,wpim,delta,mask )
1658
1659 if( ISVECTIM(inim) ){ VECTORME(inim,outim) ; RETURN(outim) ; }
1660
1661 /*------------------ here, input image is scalar-valued ------------------*/
1662 /* (convert to float type, if needed) */
1663
1664 fim = (inim->kind == MRI_float ) ? inim : mri_to_float(inim) ;
1665 far = MRI_FLOAT_PTR(fim) ;
1666
1667 outim = mri_new_conforming( fim , MRI_float ) ;
1668 outar = MRI_FLOAT_PTR(outim) ; /* is zero filled */
1669
1670 nx = fim->nx; ny = fim->ny; nz = fim->nz; nxy = nx*ny; nvox = nx*ny*nz;
1671
1672 /* component images of 3D warp */
1673
1674 iim = mri_fvect_subimage( wpim , 0 ) ;
1675 jjm = mri_fvect_subimage( wpim , 1 ) ;
1676 kkm = mri_fvect_subimage( wpim , 2 ) ;
1677 #if 0
1678 if( kkm == NULL ){ /* no k-direction info */
1679 kkm = mri_new_conforming( fim , MRI_float ) ; /* so make some up */
1680 kar = MRI_FLOAT_PTR(kkm) ;
1681 for( ii=0 ; ii < nvox ; ii++ ) kar[ii] = (ii/nxy) ;
1682 }
1683 #endif
1684
1685 /* check mask */
1686
1687 if( mask != NULL ){
1688 mmm = mask ; nmask = THD_countmask( nvox , mmm ) ;
1689 if( nmask < 1 ) nmask = nvox ;
1690 } else {
1691 nmask = nvox ; /* no mask ==> do them all */
1692 }
1693
1694 /* indexes at which to calculate output volume */
1695
1696 iar = MRI_FLOAT_PTR(iim); jar = MRI_FLOAT_PTR(jjm); kar = MRI_FLOAT_PTR(kkm);
1697
1698 /* make subset that fits the mask */
1699
1700 if( nmask == nvox ){
1701 miar = iar ; mjar = jar ; mkar = kar ; moutar = outar ;
1702 } else {
1703 miar = (float *)malloc(sizeof(float)*nmask) ;
1704 mjar = (float *)malloc(sizeof(float)*nmask) ;
1705 mkar = (float *)malloc(sizeof(float)*nmask) ;
1706 moutar = (float *)malloc(sizeof(float)*nmask) ;
1707 for( ii=jj=0 ; ii < nvox ; ii++ ){
1708 if( mmm[ii] ){
1709 miar[jj] = iar[ii] ; mjar[jj] = jar[ii] ; mkar[jj] = kar[ii] ; jj++ ;
1710 }
1711 }
1712 }
1713
1714 /* shift by delta */
1715
1716 delx = delta.a ; dely = delta.b ; delz = delta.c ;
1717 if( delx != 0.0f || dely != 0.0f || delz != 0.0f ){
1718 for( jj=0 ; jj < nmask ; jj++ ){
1719 miar[jj] += delx ; mjar[jj] += dely ; mkar[jj] += delz ;
1720 }
1721 }
1722
1723 /* compute bounds on input to apply to output */
1724
1725 do_clip = ( interp_code != MRI_NN && interp_code != MRI_LINEAR ) ;
1726 if( do_clip ){
1727 ffmin = ffmax = far[0] ;
1728 for( ii=1 ; ii < nvox ; ii++ ){
1729 if( far[ii] < ffmin ) ffmin = far[ii] ;
1730 else if( far[ii] > ffmax ) ffmax = far[ii] ;
1731 }
1732 }
1733
1734 /*-- the actual interpolation work is outsourced --*/
1735
1736 switch( interp_code ){
1737
1738 case MRI_NN:
1739 GA_interp_NN ( fim , nmask,miar,mjar,mkar,moutar ) ;
1740 break ;
1741
1742 case MRI_LINEAR:
1743 GA_interp_linear ( fim , nmask,miar,mjar,mkar,moutar ) ;
1744 break ;
1745
1746 case MRI_CUBIC:
1747 GA_interp_cubic ( fim , nmask,miar,mjar,mkar,moutar ) ;
1748 break ;
1749
1750 case MRI_VARP1:
1751 GA_interp_varp1 ( fim , nmask,miar,mjar,mkar,moutar ) ;
1752 break ;
1753
1754 case MRI_WSINC5:
1755 GA_interp_wsinc5 ( fim , nmask,miar,mjar,mkar,moutar ) ;
1756 break ;
1757
1758 default:
1759 case MRI_QUINTIC:
1760 GA_interp_quintic( fim , nmask,miar,mjar,mkar,moutar ) ;
1761 break ;
1762 }
1763
1764 /* apply the bounds */
1765
1766 if( do_clip ){
1767 for( jj=0 ; jj < nmask ; jj++ ){
1768 if( moutar[jj] < ffmin ) moutar[jj] = ffmin ;
1769 else if( moutar[jj] > ffmax ) moutar[jj] = ffmax ;
1770 }
1771 }
1772
1773 /* copy subset values in moutar back to outar image */
1774
1775 if( nmask < nvox ){
1776 for( ii=jj=0 ; ii < nvox ; ii++ ){
1777 if( mmm[ii] ) outar[ii] = moutar[jj++] ;
1778 }
1779 free(moutar) ; free(mkar) ; free(mjar) ; free(miar) ;
1780 }
1781
1782 /*--- done! ---*/
1783
1784 mri_free(kkm) ; mri_free(jjm) ; mri_free(iim) ;
1785 if( fim != inim ) mri_free(fim) ;
1786
1787 RETURN(outim) ;
1788 }
1789
1790 /*---------------------------------------------------------------------------*/
1791 /*! Apply a matrix to a set of warp vectors (in place). */
1792
GA_affine_edit_warp(mat44 aff,MRI_IMAGE * wpim)1793 void GA_affine_edit_warp( mat44 aff , MRI_IMAGE *wpim )
1794 {
1795 int ii , nvox ;
1796 float *war , aa,bb,cc ;
1797
1798 ENTRY("GA_affine_edit_warp") ;
1799
1800 if( !ISVALID_MAT44(aff) || wpim == NULL ) EXRETURN ;
1801
1802 if( wpim->kind != MRI_fvect ||
1803 mri_data_pointer(wpim) == NULL || wpim->vdim != 3 ) EXRETURN ;
1804
1805 nvox = wpim->nvox ;
1806 war = (float *)mri_data_pointer(wpim) ;
1807
1808 for( ii=0 ; ii < nvox ; ii++ ){
1809 aa = war[3*ii ] ;
1810 bb = war[3*ii+1] ;
1811 cc = war[3*ii+2] ;
1812 MAT44_VEC( aff , aa,bb,cc , war[3*ii],war[3*ii+1],war[3*ii+2] ) ;
1813 }
1814
1815 EXRETURN ;
1816 }
1817
1818 /*---------------------------------------------------------------------------*/
1819 /* B(A(x)) */
1820
1821 #if 0
1822 MRI_IMAGE * GA_compose_indexwarp( MRI_IMAGE *awpim , MRI_IMAGE *bwpim )
1823 {
1824
1825 ENTRY("GA_compose_indexwarp") ;
1826
1827 if( awpim == NULL || awpim->kind != MRI_fvect ) RETURN(NULL);
1828 if( bwpim == NULL || bwpim->kind != MRI_fvect ) RETURN(NULL);
1829 if( mri_data_pointer(awpim) == NULL || awpim->vdim != 3 ) RETURN(NULL) ;
1830 if( mri_data_pointer(bwpim) == NULL || bwpim->vdim != 3 ) RETURN(NULL) ;
1831 if( awpim->nx != bwpim->nx ||
1832 awpim->ny != bwpim->ny || awpim->nz != bwpim->nz ) RETURN(NULL);
1833
1834 cwpim = GA_indexwarp( bwpim
1835 #endif
1836