1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #undef SL_DEBUG /* 17 Oct 2006: some extra debugging stuff */
8
9 #undef MAIN
10
11 /*********************************************************************
12 These are the template routines to fill in a slice of data
13 in a new image array from an old image brick array, warped on demand.
14
15 To create actual routines for this purpose, you must compile
16 the file with the preprocessor symbol DTYPE set to one of
17 the following types:
18
19 byte short int float double complex rgbyte rgba
20
21 cc -c -DDTYPE=short afni_slice.c
22 mv -f afni_slice.o afni_slice_short.o
23
24 In the example above, the resulting routines will be named
25
26 AFNI_lmap_to_xslice_short
27 AFNI_lmap_to_yslice_short
28 AFNI_lmap_to_zslice_short
29 AFNI_br2sl_short
30
31 and will take as input pointers bold and bslice, pointers of
32 type "short *".
33 *********************************************************************/
34
35 #ifndef DTYPE
36 #error "Cannot compile, since DTYPE is undefined."
37 #endif
38
39 #include "afni_warp.h"
40
41 /** macros for function names defined in this file **/
42
43 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
44 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
45 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
46 #define B2SL_NAME TWO_TWO(AFNI_br2sl_,DTYPE)
47
48 #define XSTRING(qqq) STRING(qqq)
49 #define STRING(nnn) #nnn
50
51 /*******************************************************************
52 To add a new DTYPE, you must add definitions in mrilib.h,
53 and define the following macros for the new type:
54 FMAD2, FMAD4, FSCAL, FINAL, FZERO, and INTYPE
55 ********************************************************************/
56
57 /** macros for e = a*d1 + b*d2 (a,b floats; d1,d2 DTYPEs) **/
58
59 #define FMAD2_short(a,d1,b,d2,e) (e)=(a)*(d1)+(b)*(d2)
60 #define FMAD2_float FMAD2_short
61 #define FMAD2_byte FMAD2_short
62 #define FMAD2_int FMAD2_short
63 #define FMAD2_double FMAD2_short
64
65 #define FMAD2_complex(a,d1,b,d2,e) ( (e).r = (a)*(d1).r + (b)*(d2).r, \
66 (e).i = (a)*(d1).i + (b)*(d2).i )
67
68 #define FMAD2_rgbyte(a,d1,bb,d2,e) ( (e).r = (a)*(d1).r + (bb)*(d2).r, \
69 (e).g = (a)*(d1).g + (bb)*(d2).g, \
70 (e).b = (a)*(d1).b + (bb)*(d2).b )
71
72 #define FMAD2_rgba(p,d1,bb,d2,e) ( (e).r = (p)*(d1).r + (bb)*(d2).r, \
73 (e).g = (p)*(d1).g + (bb)*(d2).g, \
74 (e).b = (p)*(d1).b + (bb)*(d2).b, \
75 (e).a = (p)*(d1).a + (bb)*(d2).a )
76
77 #define FMAD2 TWO_TWO(FMAD2_,DTYPE)
78
79 /** macros for e = a*d1 + b*d2 + c*d3 + d*d3 (a-d floats; d1-d4 DTYPEs) **/
80
81 #define FMAD4_short(a,d1,b,d2,c,d3,d,d4,e) (e)=(a)*(d1)+(b)*(d2)+(c)*(d3)+(d)*(d4)
82 #define FMAD4_float FMAD4_short
83 #define FMAD4_byte FMAD4_short
84 #define FMAD4_int FMAD4_short
85 #define FMAD4_double FMAD4_short
86
87 #define FMAD4_complex(a,d1,b,d2,c,d3,d,d4,e) \
88 ( (e).r = (a)*(d1).r + (b)*(d2).r + (c)*(d3).r + (d)*(d4).r, \
89 (e).i = (a)*(d1).i + (b)*(d2).i + (c)*(d3).i + (d)*(d4).i )
90
91 #define FMAD4_rgbyte(a,d1,bb,d2,c,d3,d,d4,e) \
92 ( (e).r = (a)*(d1).r + (bb)*(d2).r + (c)*(d3).r + (d)*(d4).r, \
93 (e).g = (a)*(d1).g + (bb)*(d2).g + (c)*(d3).g + (d)*(d4).g, \
94 (e).b = (a)*(d1).b + (bb)*(d2).b + (c)*(d3).b + (d)*(d4).b )
95
96 #define FMAD4_rgba(p,d1,bb,d2,c,d3,d,d4,e) \
97 ( (e).r = (p)*(d1).r + (bb)*(d2).r + (c)*(d3).r + (d)*(d4).r, \
98 (e).g = (p)*(d1).g + (bb)*(d2).g + (c)*(d3).g + (d)*(d4).g, \
99 (e).b = (p)*(d1).b + (bb)*(d2).b + (c)*(d3).b + (d)*(d4).b, \
100 (e).a = (p)*(d1).a + (bb)*(d2).a + (c)*(d3).a + (d)*(d4).a )
101
102 #define FMAD4 TWO_TWO(FMAD4_,DTYPE)
103
104 /** macros to multiply float a times DTYPE b and store the result in b again **/
105
106 #define FSCAL_short(a,b) (b)*=(a)
107 #define FSCAL_float FSCAL_short
108 #define FSCAL_byte FSCAL_short
109 #define FSCAL_int FSCAL_short
110 #define FSCAL_double FSCAL_short
111 #define FSCAL_complex(a,b) ( (b).r *= (a) , (b).i *= (a) )
112
113 #define FSCAL_rgbyte(a,bb) ( (bb).r*= (a) , (bb).g*= (a) , (bb).b*= (a) )
114
115 #define FSCAL_rgba(p,bb) ( (bb).r*=(p), (bb).g*=(p), (bb).b*=(p), (bb).a*=(p) )
116
117 #define FSCAL TWO_TWO(FSCAL_,DTYPE)
118
119 /** macros for assigning final result from INTYPE a to DTYPE b **/
120
121 /* 18 Nov 1998: modify FINAL_short and FINAL_byte to prevent overflow */
122
123 #define CLIP_OVERFLOW
124 #ifdef CLIP_OVERFLOW
125
126 # define FINAL_short(a,b) (b) = ( ((a)<-32767.0) ? (-32767) : \
127 ((a)> 32767.0) ? ( 32767) : ((short)((a)+0.5)) )
128
129 # define FINAL_byte(a,b) (b) = ( ((a)< 0.0) ? (0) : \
130 ((a)> 255.0) ? (255) : ((byte)((a)+0.5)) )
131
132 #else
133 # define FINAL_short(a,b) (b)=((short)((a)+0.5))
134 # define FINAL_byte(a,b) (b)=((byte)((a)+0.5))
135 #endif
136
137 #define FINAL_int(a,b) (b)=((int)((a)+0.5))
138 #define FINAL_float(a,b) (b)=(a)
139 #define FINAL_double FINAL_float
140 #define FINAL_complex FINAL_float
141 #define FINAL_rgbyte FINAL_float
142 #define FINAL_rgba FINAL_float
143 #define FINAL TWO_TWO(FINAL_,DTYPE)
144
145 /** macros for putting a zero into DTYPE b **/
146
147 #define FZERO_short(b) (b)=0
148 #define FZERO_byte FZERO_short
149 #define FZERO_int FZERO_short
150 #define FZERO_float(b) (b)=0.0
151 #define FZERO_double FZERO_float
152 #define FZERO_complex(b) ( (b).r = 0.0 , (b).i = 0.0 )
153 #define FZERO_rgbyte(bb) ( (bb).r=(bb).g=(bb).b = 0 )
154 #define FZERO_rgba(bb) ( (bb).r=(bb).g=(bb).b=(bb).a = 0 )
155 #define FZERO TWO_TWO(FZERO_,DTYPE)
156
157 /** macros for a zero value **/
158
159 static complex complex_zero = { 0.0f,0.0f } ;
160 static rgbyte rgbyte_zero = { 0,0,0 } ;
161 static rgba rgba_zero = { 0,0,0,0 } ;
162
163 #define ZERO_short 0
164 #define ZERO_byte 0
165 #define ZERO_int 0
166 #define ZERO_float 0.0f
167 #define ZERO_double 0.0
168 #define ZERO_complex complex_zero
169 #define ZERO_rgbyte rgbyte_zero
170 #define ZERO_rgba rgba_zero
171 #define ZERO TWO_TWO(ZERO_,DTYPE)
172
173 /** macros for intermediate interpolants data type **/
174
175 #define INTYPE_short float
176 #define INTYPE_float float
177 #define INTYPE_byte float
178 #define INTYPE_int float
179 #define INTYPE_double double
180 #define INTYPE_complex complex
181 #define INTYPE_rgbyte rgbyte
182 #define INTYPE_rgba rgba
183 #define INTYPE TWO_TWO(INTYPE_,DTYPE)
184
185 /**-------------------------------------------------------------------------**/
186 /**----- macros for quickly doing NN interpolation along parallel axes -----**/
187
188 /** test and flags for which steps are zero **/
189
190 #define ZZZ 1.e-5 /* effectively zero */
191 #define NONE_ZERO 0
192 #define X_ZERO 1
193 #define Y_ZERO 2
194 #define XY_ZERO 3
195 #define Z_ZERO 4
196 #define XZ_ZERO 5
197 #define YZ_ZERO 6
198 #define XYZ_ZERO 7
199 #define OUTADD 100
200 #define THREEZ(x,y,z) ((fabs(x)<ZZZ) + 2*(fabs(y)<ZZZ) + 4*(fabs(z)<ZZZ))
201
202 /** ALOOP: inner loop statements for the actual NN assignments,
203 and there is no need to check for being outside the input brick
204 _GEN: for the general (non-parallel) case
205 _PAR: for the case when the inner loop is parallel to an input brick axis **/
206
207 #ifndef SL_DEBUG
208
209 #define NN_ALOOP_GEN \
210 (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \
211 xi_old = FLOOR(fxi_old) , yj_old = FLOOR(fyj_old) , zk_old = FLOOR(fzk_old) , \
212 bslice[out_ind++] = bold[ IBASE(xi_old,yj_old,zk_old) ])
213
214 #else
215
216 #define NN_ALOOP_GEN \
217 do{ \
218 fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \
219 xi_old = FLOOR(fxi_old) , yj_old = FLOOR(fyj_old) , zk_old = FLOOR(fzk_old) ; \
220 if(PRINT_TRACING){ char str[256]; \
221 sprintf(str,"out_ind=%d xi_old=%d yj_old=%d zk_old=%d",out_ind,xi_old,yj_old,zk_old); \
222 STATUS(str) ; } \
223 bslice[out_ind++] = bold[ IBASE(xi_old,yj_old,zk_old) ] ; \
224 } while(0)
225
226 #endif /* SL_DEBUG */
227
228 #define NN_ALOOP_PAR(ijk) (bslice[out_ind++] = bold[ ib[ijk]+ob ])
229
230 /** BLOOP: assign values to the ib array that will be used to
231 index into the input brick for rapid access;
232 there is one BLOOP for each possible parallel axis **/
233
234 #define NN_BLOOP_XY_ZERO(ijk) \
235 (fzk_old += dfzk_inner , zk_old = GFLOOR(fzk_old) , ib[ijk] = IBASE(0,0,zk_old))
236
237 #define NN_BLOOP_XZ_ZERO(ijk) \
238 (fyj_old += dfyj_inner , yj_old = GFLOOR(fyj_old) , ib[ijk] = IBASE(0,yj_old,0))
239
240 #define NN_BLOOP_YZ_ZERO(ijk) \
241 (fxi_old += dfxi_inner , xi_old = GFLOOR(fxi_old) , ib[ijk] = IBASE(xi_old,0,0))
242
243 /** macros to test if the point we want to resample from is outside the old array **/
244
245 #define TEST_OUT_XXX ( fxi_old < fxi_bot || fxi_old > fxi_top )
246 #define TEST_OUT_YYY ( fyj_old < fyj_bot || fyj_old > fyj_top )
247 #define TEST_OUT_ZZZ ( fzk_old < fzk_bot || fzk_old > fzk_top )
248
249 #define TEST_OUT_ALL (TEST_OUT_XXX || TEST_OUT_YYY || TEST_OUT_ZZZ)
250
251 /** CLOOP: like the ALOOP, but for when we must test
252 if the desired point may be outside the input brick array **/
253
254 #define NN_CLOOP_GEN \
255 (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \
256 bslice[out_ind++] = (TEST_OUT_ALL) ? ZERO : \
257 bold[IBASE(FLOOR(fxi_old),FLOOR(fyj_old),FLOOR(fzk_old))] )
258
259 #define NN_CLOOP_PAR(ijk) \
260 (bslice[out_ind++] = (ib[ijk]<0 || ib[ijk]>=ub) ? ZERO : bold[ib[ijk]+ob])
261
262 /** ZLOOP: just assign zero to each output **/
263
264 #define NN_ZLOOP (bslice[out_ind++] = ZERO)
265
266 /** space for precomputed indices **/
267
268 static int * ib = NULL ;
269 static int nib = -1 ;
270
271 /** macro to make the ib array as big as we need it **/
272
273 #define MAKE_IBIG(top) do{ if(nib < (top)){ \
274 if(ib != NULL) free(ib) ; \
275 ib = (int *) malloc(sizeof(int)*((top)+9)) ; \
276 if(ib==NULL){ \
277 fprintf(stderr,"\nmalloc fails in NN reslice!\n");EXIT(1);} \
278 nib = (top) ; } } while(0)
279
280 /*---------------------------------------------------------------------
281 routine to apply a linear mapping to a dataset and fill in
282 a single fixed-x slice from a 3D dataset:
283
284 map = linear mapping from old to new
285 (bot & top field contain index limits in new)
286 resam_mode = type of interpolation to do in old
287
288 old_daxes = axis information of old dataset
289 bold = pointer to old dataset's 3D data brick array
290
291 new_daxes = axis information of new dataset
292 xi_fix = fixed x index in new dataset
293 bslice = pointer to nynew * nznew array to hold slice
294
295 Cognate routines for yslice (output is nznew * nxnew: Y-ZX) and
296 zslize (output is nxnew * nynew: Z-XY)
297 follow this one directly.
298 -----------------------------------------------------------------------*/
299
300 #define IBASE(i,j,k) ((i)+(j)*jstep+(k)*kstep)
301 #define ROUND(qq) ((int)(qq+0.5))
302 #define FLOOR(qq) ((int)(qq)) /* cheap and fast */
303 #define GFLOOR(qq) ((int)floor(qq)) /* good and slow */
304
305 /* define linear interpolation polynomials */
306
307 #define LP_00(x) (1.0-(x))
308 #define LP_P1(x) (x)
309
310 /* define blocky interpolation functions */
311
312 #define BP_hh(x) (8.0*((x)*(x))*((x)*(x)))
313 #define BP_00(x) ( ((x)<0.5) ? (1-BP_hh(x)) : ( BP_hh(1-(x))) )
314 #define BP_P1(x) ( ((x)<0.5) ? ( BP_hh(x)) : (1-BP_hh(1-(x))) )
315
316 /* define cubic interpolation polynomials */
317
318 #define CP_M1(x) (-(x)*((x)-1)*((x)-2))
319 #define CP_00(x) (3.0*((x)+1)*((x)-1)*((x)-2))
320 #define CP_P1(x) (-3.0*(x)*((x)+1)*((x)-2))
321 #define CP_P2(x) ((x)*((x)+1)*((x)-1))
322 #define CP_FACTOR 4.62962963e-3 /* 1/216 = final scaling factor */
323
324 #define FXYZTMP(xx,yy,zz) \
325 ( fxi_tmp = mt.mat[0][0] * xx \
326 + mt.mat[0][1] * yy \
327 + mt.mat[0][2] * zz - vt.xyz[0] , \
328 fyj_tmp = mt.mat[1][0] * xx \
329 + mt.mat[1][1] * yy \
330 + mt.mat[1][2] * zz - vt.xyz[1] , \
331 fzk_tmp = mt.mat[2][0] * xx \
332 + mt.mat[2][1] * yy \
333 + mt.mat[2][2] * zz - vt.xyz[2] )
334
335
LMAP_XNAME(THD_linear_mapping * map,int resam_mode,THD_dataxes * old_daxes,DTYPE * bold,THD_dataxes * new_daxes,int xi_fix,DTYPE * bslice)336 void LMAP_XNAME( THD_linear_mapping * map , int resam_mode ,
337 THD_dataxes * old_daxes , DTYPE * bold ,
338 THD_dataxes * new_daxes , int xi_fix , DTYPE * bslice )
339 {
340 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */
341 THD_fvec3 vt = map->svec ;
342
343 int xi_new , yj_new , zk_new ; /* voxel indices in new */
344 int xi_old , yj_old , zk_old ; /* voxel indices in old */
345 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */
346
347 float fxi_top , fyj_top , fzk_top ; /* floating pt. voxel indices */
348 float fxi_bot , fyj_bot , fzk_bot ;
349 float fxi_base , fyj_base , fzk_base ;
350 float dfxi_outer , dfyj_outer , dfzk_outer ;
351 float dfxi_inner , dfyj_inner , dfzk_inner ;
352
353 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */
354 int out_ind , jstep , kstep ;
355 int nxold,nyold,nzold , nxnew,nynew,nznew ;
356
357 #if 0
358 ENTRY("AFNI_lmap_to_xslice") ;
359 #else
360 ENTRY( XSTRING(LMAP_XNAME) ) ;
361 #endif
362
363 /*--- set up ranges ---*/
364
365 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ;
366 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ;
367 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ;
368
369 if( xi_fix < xi_bot || xi_fix > xi_top ) EXRETURN ; /* map doesn't apply! */
370
371 if( bold==NULL || bslice==NULL ){ STATUS("NULL pointers?!"); EXRETURN; }
372
373 #ifdef SL_DEBUG
374 if(PRINT_TRACING)
375 { char str[256];
376 sprintf(str,"xi_bot=%d xi_fix=%d xi_top=%d",xi_bot,xi_fix,xi_top);
377 STATUS(str) ;
378 sprintf(str,"yj_bot=%d yj_top=%d zk_bot=%d zk_top=%d",
379 yj_bot,yj_top,zk_bot,zk_top) ; STATUS(str) ;
380 }
381 #endif
382
383 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ;
384 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ;
385 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ;
386
387 #ifdef SL_DEBUG
388 if(PRINT_TRACING)
389 { char str[256];
390 sprintf(str,"OLD: nx=%d ny=%d nz=%d NEW: nx=%d ny=%d nz=%d",
391 nxold,nyold,nzold , nxnew,nynew,nznew ) ; STATUS(str); }
392 #endif
393
394 jstep = nxold ;
395 kstep = nxold * nyold ;
396
397 /* set up base of indices in old */
398
399 xi_new = xi_fix ;
400 yj_new = yj_bot-1 ;
401 zk_new = zk_bot-1 ;
402
403 fxi_base = mt.mat[0][0] * xi_new
404 + mt.mat[0][1] * yj_new
405 + mt.mat[0][2] * zk_new - vt.xyz[0] ;
406
407 fyj_base = mt.mat[1][0] * xi_new
408 + mt.mat[1][1] * yj_new
409 + mt.mat[1][2] * zk_new - vt.xyz[1] ;
410
411 fzk_base = mt.mat[2][0] * xi_new
412 + mt.mat[2][1] * yj_new
413 + mt.mat[2][2] * zk_new - vt.xyz[2] ;
414
415 dfxi_outer = mt.mat[0][2] ; /* outer loop is in z = 2 */
416 dfyj_outer = mt.mat[1][2] ;
417 dfzk_outer = mt.mat[2][2] ;
418
419 dfxi_inner = mt.mat[0][1] ; /* inner loop is in y = 1 */
420 dfyj_inner = mt.mat[1][1] ;
421 dfzk_inner = mt.mat[2][1] ;
422
423 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ;
424 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ;
425 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ;
426
427 #ifdef SL_DEBUG
428 if(PRINT_TRACING)
429 { char str[256] ;
430 sprintf(str,"switch(resam_mode=%d)",resam_mode); STATUS(str); }
431 #endif
432
433 switch( resam_mode ){
434
435 default:
436 case RESAM_NN_TYPE:{
437 float fxi_max , fyj_max , fzk_max ;
438 float fxi_min , fyj_min , fzk_min ;
439 float fxi_tmp , fyj_tmp , fzk_tmp ;
440 int any_outside , all_outside ;
441
442 #ifdef SL_DEBUG
443 if(PRINT_TRACING)
444 { char str[256] ;
445 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g",
446 dfxi_inner,dfyj_inner,dfzk_inner,
447 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; }
448 #endif
449 /** July 15, 1996:
450 check if all the points are inside the old grid;
451 if so, can use a version of the resampling loop
452 that does not need to check each voxel for being
453 inside -- hopefully, this will execute more quickly **/
454
455 FXYZTMP(xi_new,yj_bot,zk_bot) ;
456 fxi_max = fxi_min = fxi_tmp ;
457 fyj_max = fyj_min = fyj_tmp ;
458 fzk_max = fzk_min = fzk_tmp ;
459
460 FXYZTMP(xi_new,yj_top,zk_bot) ;
461 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
462 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
463 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
464
465 FXYZTMP(xi_new,yj_bot,zk_top) ;
466 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
467 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
468 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
469
470 FXYZTMP(xi_new,yj_top,zk_top) ;
471 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
472 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
473 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
474
475 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) ||
476 (fyj_min < fyj_bot) || (fyj_max > fyj_top) ||
477 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ;
478
479 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) ||
480 (fyj_max < fyj_bot) || (fyj_min > fyj_top) ||
481 (fzk_max < fzk_bot) || (fzk_min > fzk_top)
482 : 0 ;
483 #ifdef SL_DEBUG
484 if(PRINT_TRACING)
485 { char str[256] ;
486 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max);
487 STATUS(str) ;
488 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max);
489 STATUS(str) ;
490 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max);
491 STATUS(str) ; }
492 #endif
493
494 /** redefine the macros specifying loop variables **/
495
496 #undef OUD_NAME
497 #undef IND_NAME
498 #undef OUD
499 #undef OUD_bot
500 #undef OUD_top
501 #undef IND
502 #undef IND_bot
503 #undef IND_top
504 #undef IND_nnn
505
506 #define OUD_NAME zk /* name of 2nd dimension of output image */
507 #define IND_NAME yj /* name of 1st dimension of output image */
508 #define IND_nnn nynew /* inner loop image dimension */
509
510 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */
511 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */
512 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */
513 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */
514 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */
515 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */
516
517 if( all_outside ){
518 #ifdef SL_DEBUG
519 STATUS("NN resample has all outside") ;
520 #endif
521 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */
522 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */
523 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */
524 bslice[out_ind++] = ZERO ;
525 }
526 }
527 } else { /* at least some are inside */
528
529 int thz , tho , ob , ub=0 ;
530
531 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ;
532
533 fxi_top = nxold-0.0001 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */
534 fyj_top = nyold-0.0001 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */
535 fzk_top = nzold-0.0001 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */
536 /* these variables now. */
537
538 /** thz = flag that indicates which of the steps df??_inner are zero.
539 If two of them are zero, then the inner loop is parallel to
540 one of the input brick axes, and so data may be pulled
541 out in a very efficient fashion. In such a case, precompute
542 the indexes for the inner loop:
543
544 the BLOOP macros load array ib, which holds the inner loop
545 computed indexes for each inner loop position.
546 ub = upper bound value for ib array value to still be
547 inside input brick array.
548 ob = outer loop index into input brick array (computed later) **/
549
550 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */
551 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */
552 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */
553 else /* treatment if outer */
554 thz = NONE_ZERO ; /* axes are special */
555
556 #ifdef SL_DEBUG
557 if(PRINT_TRACING)
558 { char str[256] ;
559 if( any_outside ) sprintf(str,"NN resample has some outside: thz=%d tho=%d",thz,tho);
560 else sprintf(str,"NN resample has all inside: thz=%d thd=%d",thz,tho);
561 STATUS(str) ;
562 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d",
563 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ;
564 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d",
565 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; }
566 #endif
567
568 switch(thz){
569 case XY_ZERO:
570 MAKE_IBIG(IND_top) ;
571 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ;
572 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; }
573 break ;
574
575 case XZ_ZERO:
576 MAKE_IBIG(IND_top) ;
577 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ;
578 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; }
579 break ;
580
581 case YZ_ZERO:
582 MAKE_IBIG(IND_top) ;
583 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ;
584 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; }
585 break ;
586 }
587
588 thz += OUTADD * any_outside ;
589
590 #ifdef SL_DEBUG
591 STATUS("beginning NN outer loop") ;
592 #endif
593
594 /*** outer loop ***/
595
596 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){
597
598 #ifdef SL_DEBUG
599 if(PRINT_TRACING){
600 char str[256] ; sprintf(str,"index=%d",OUD_bot); STATUS(str); }
601 #endif
602
603 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */
604 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */
605 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */
606
607 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */
608
609 /*** There are 8 cases for the inner loop:
610 all inside, inner loop not parallel to any input axis
611 all inside, inner loop parallel to input brick z-axis
612 all inside, inner loop parallel to input brick y-axis
613 all inside, inner loop parallel to input brick x-axis
614
615 and then the 4 same cases repeated when not all desired
616 points are inside the input brick. Each of these is
617 coded separately for efficiency. This is important for
618 rapid re-display of results during interactive imaging. ***/
619
620 switch(thz){
621 case NONE_ZERO:
622 case X_ZERO:
623 case Y_ZERO:
624 case Z_ZERO:
625 case XYZ_ZERO:
626 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; }
627 break ;
628
629 case XY_ZERO:
630 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
631 ob = IBASE(xi_old,yj_old,0) ;
632 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
633 break ;
634
635 case XZ_ZERO:
636 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
637 ob = IBASE(xi_old,0,zk_old) ;
638 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
639 break ;
640
641 case YZ_ZERO:
642 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
643 ob = IBASE(0,yj_old,zk_old) ;
644 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
645 break ;
646
647 default:
648 case NONE_ZERO+OUTADD:
649 case X_ZERO+OUTADD:
650 case Y_ZERO+OUTADD:
651 case Z_ZERO+OUTADD:
652 case XYZ_ZERO+OUTADD:
653 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; }
654 break ;
655
656 case XY_ZERO+OUTADD:
657 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
658 if( TEST_OUT_XXX || TEST_OUT_YYY ){
659 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
660 } else {
661 ob = IBASE(xi_old,yj_old,0) ;
662 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
663 }
664 break ;
665
666 case XZ_ZERO+OUTADD:
667 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
668 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){
669 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
670 } else {
671 ob = IBASE(xi_old,0,zk_old) ;
672 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
673 }
674 break ;
675
676 case YZ_ZERO+OUTADD:
677 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
678 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){
679 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
680 } else {
681 ob = IBASE(0,yj_old,zk_old) ;
682 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
683 }
684 break ;
685 }
686 }
687 }
688 }
689 break ; /* end of NN! */
690
691 case RESAM_BLOCK_TYPE:
692 case RESAM_LINEAR_TYPE:{
693 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ;
694 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 ,
695 f_k00 , f_kp1 , result ;
696 float frac_xi , frac_yj , frac_zk ;
697 int ibase ,
698 in_jp1_k00 = jstep ,
699 in_j00_kp1 = kstep ,
700 in_jp1_kp1 = jstep+kstep ;
701 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
702
703 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){
704
705 fxi_old = (fxi_base += dfxi_outer) ;
706 fyj_old = (fyj_base += dfyj_outer) ;
707 fzk_old = (fzk_base += dfzk_outer) ;
708
709 out_ind = yj_bot + zk_new * nynew ;
710
711 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){
712
713 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
714
715 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
716 fyj_old < fyj_bot || fyj_old > fyj_top ||
717 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
718
719 /***
720 bslice[out_ind++] = 0 ;
721 ***/
722 FZERO(bslice[out_ind]) ; out_ind++ ;
723 continue ;
724 }
725
726 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
727 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
728 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
729
730 /* use NN if at edges of old brick */
731
732 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
733 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
734
735 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
736 ROUND(fyj_old) ,
737 ROUND(fzk_old) ) ] ;
738 continue ;
739 }
740
741 /* compute weights for LINEAR interpolation in each direction */
742
743 if( resam_mode == RESAM_BLOCK_TYPE ){
744 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ;
745 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ;
746 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ;
747 } else {
748 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
749 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
750 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
751 }
752
753 /* interpolate in the x direction for each y & z */
754
755 ibase = IBASE(xi_old,yj_old,zk_old) ;
756 /***
757 f_j00_k00 = xwt_00 * bold[ibase]
758 + xwt_p1 * bold[ibase+1] ;
759
760 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
761 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
762
763 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
764 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
765
766 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
767 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
768 ***/
769 FMAD2( xwt_00 , bold[ibase] ,
770 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
771
772 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
773 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
774
775 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
776 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
777
778 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
779 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
780
781 /* interpolate in the y direction for each z */
782
783 /***
784 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
785 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
786 ***/
787 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
788 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
789
790 /* interpolate in the z direction and
791 put the result into the output array!
792 (the +0.5 is to force rounding of the result) */
793
794 /***
795 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
796 ***/
797 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
798 FINAL( result , bslice[out_ind] ) ;
799 out_ind++ ;
800
801 }
802 }
803 }
804 break ;
805
806 case RESAM_CUBIC_TYPE:{
807 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */
808 ywt_m1,ywt_00,ywt_p1,ywt_p2 ,
809 zwt_m1,zwt_00,zwt_p1,zwt_p2 ;
810
811 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */
812 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 ,
813 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 ,
814 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 ,
815 f_km1 , f_k00 , f_kp1 , f_kp2 , result ;
816 float frac_xi , frac_yj , frac_zk ;
817
818 int ibase , /* base index for interpolant */
819 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */
820 in_jm1_k00 = -jstep , /* 0 (00) */
821 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */
822 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */
823 in_j00_km1 = - kstep , /* steps in j and k indices */
824 in_j00_k00 = 0 ,
825 in_j00_kp1 = kstep ,
826 in_j00_kp2 = 2*kstep ,
827 in_jp1_km1 = jstep- kstep ,
828 in_jp1_k00 = jstep ,
829 in_jp1_kp1 = jstep+ kstep ,
830 in_jp1_kp2 = jstep+2*kstep , /* no '2*' [PT: 27 Jul 2021] */
831 in_jp2_km1 =2*jstep- kstep ,
832 in_jp2_k00 =2*jstep ,
833 in_jp2_kp1 =2*jstep+ kstep ,
834 in_jp2_kp2 =2*jstep+2*kstep ;
835
836 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
837 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ;
838
839 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){
840
841 fxi_old = (fxi_base += dfxi_outer) ;
842 fyj_old = (fyj_base += dfyj_outer) ;
843 fzk_old = (fzk_base += dfzk_outer) ;
844
845 out_ind = yj_bot + zk_new * nynew ;
846
847 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){
848
849 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
850
851 /* check if outside old brick */
852
853 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
854 fyj_old < fyj_bot || fyj_old > fyj_top ||
855 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
856 /***
857 bslice[out_ind++] = 0 ;
858 ***/
859 FZERO(bslice[out_ind]) ; out_ind++ ;
860 continue ;
861 }
862
863 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
864 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
865 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
866
867 /* use NN if at very edges of old brick */
868
869 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
870 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
871
872 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
873 ROUND(fyj_old) ,
874 ROUND(fzk_old) ) ] ;
875 continue ;
876 }
877
878 ibase = IBASE(xi_old,yj_old,zk_old) ;
879
880 /* use LINEAR if close to edges of old brick */
881
882 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 ||
883 xi_old == 0 || yj_old == 0 || zk_old == 0 ){
884
885 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
886 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
887 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
888
889 /***
890 f_j00_k00 = xwt_00 * bold[ibase]
891 + xwt_p1 * bold[ibase+1] ;
892
893 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
894 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
895
896 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
897 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
898
899 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
900 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
901 ***/
902 FMAD2( xwt_00 , bold[ibase] ,
903 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
904
905 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
906 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
907
908 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
909 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
910
911 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
912 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
913
914 /***
915 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
916 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
917 ***/
918 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
919 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
920 /***
921 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
922 ***/
923 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
924 FINAL( result , bslice[out_ind] ) ;
925 out_ind++ ;
926 continue ;
927 }
928
929 /* compute weights for CUBIC interpolation in each direction */
930
931 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ;
932 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ;
933
934 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ;
935 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ;
936
937 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ;
938 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ;
939
940 /* use the ANSI token-merge operator ## to create an interpolating
941 macro for the x-direction, at each offset in y (j) and z (k) */
942
943 /***
944 #define CXINT(j,k) xwt_m1 * bold[ibase + in_j ## j ## _k ## k -1 ] \
945 + xwt_00 * bold[ibase + in_j ## j ## _k ## k ] \
946 + xwt_p1 * bold[ibase + in_j ## j ## _k ## k +1 ] \
947 + xwt_p2 * bold[ibase + in_j ## j ## _k ## k +2 ]
948 ***/
949
950 #define CXINT(j,k,ff) \
951 FMAD4( xwt_m1 , bold[ibase + in_j ## j ## _k ## k -1 ] , \
952 xwt_00 , bold[ibase + in_j ## j ## _k ## k ] , \
953 xwt_p1 , bold[ibase + in_j ## j ## _k ## k +1 ] , \
954 xwt_p2 , bold[ibase + in_j ## j ## _k ## k +2 ] , ff )
955
956 /* interpolate in the x direction for each y & z */
957
958 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ;
959 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ;
960
961 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ;
962 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ;
963
964 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ;
965 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ;
966
967 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ;
968 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ;
969
970 /* interpolate in the y direction for each z */
971
972 /***
973 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1
974 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ;
975
976 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00
977 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ;
978
979 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1
980 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ;
981
982 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2
983 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ;
984 ***/
985 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 ,
986 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ;
987
988 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 ,
989 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ;
990
991 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 ,
992 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ;
993
994 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 ,
995 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ;
996
997 /* interpolate in the z direction and
998 put the result into the output array!
999 (the +0.5 is to force rounding of the result) */
1000 /***
1001 bslice[out_ind++] = 0.5 + CP_FACTOR
1002 * ( zwt_m1 * f_km1 + zwt_00 * f_k00
1003 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ;
1004 ***/
1005 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 ,
1006 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ;
1007 FSCAL(CP_FACTOR,result) ;
1008 FINAL( result , bslice[out_ind] ) ;
1009 out_ind++ ;
1010
1011 } /* end of inner loop */
1012 } /* end of outer loop */
1013 }
1014 break ;
1015
1016 }
1017
1018 EXRETURN ;
1019 }
1020
1021 /*--------------------------------------------------------------------------*/
1022
LMAP_YNAME(THD_linear_mapping * map,int resam_mode,THD_dataxes * old_daxes,DTYPE * bold,THD_dataxes * new_daxes,int yj_fix,DTYPE * bslice)1023 void LMAP_YNAME( THD_linear_mapping * map , int resam_mode ,
1024 THD_dataxes * old_daxes , DTYPE * bold ,
1025 THD_dataxes * new_daxes , int yj_fix , DTYPE * bslice )
1026 {
1027 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */
1028 THD_fvec3 vt = map->svec ;
1029
1030 int xi_new , yj_new , zk_new ; /* voxel indices in new */
1031 int xi_old , yj_old , zk_old ; /* voxel indices in old */
1032 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */
1033
1034 float fxi_top , fyj_top , fzk_top ;
1035 float fxi_bot , fyj_bot , fzk_bot ;
1036 float fxi_base , fyj_base , fzk_base ;
1037 float dfxi_outer , dfyj_outer , dfzk_outer ;
1038 float dfxi_inner , dfyj_inner , dfzk_inner ;
1039
1040 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */
1041 int out_ind , jstep , kstep ;
1042 int nxold,nyold,nzold , nxnew,nynew,nznew ;
1043
1044 #if 0
1045 ENTRY("AFNI_lmap_to_yslice") ;
1046 #else
1047 ENTRY( XSTRING(LMAP_YNAME) ) ;
1048 #endif
1049
1050 /*--- set up ranges ---*/
1051
1052 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ;
1053 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ;
1054 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ;
1055
1056 if( yj_fix < yj_bot || yj_fix > yj_top ) EXRETURN ; /* map doesn't apply! */
1057
1058 if( bold==NULL || bslice==NULL ){ STATUS("NULL pointers?!"); EXRETURN; }
1059
1060 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ;
1061 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ;
1062 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ;
1063
1064 jstep = nxold ;
1065 kstep = nxold * nyold ;
1066
1067 /* set up base of indices in old */
1068
1069 xi_new = xi_bot-1 ;
1070 yj_new = yj_fix ;
1071 zk_new = zk_bot-1 ;
1072
1073 fxi_base = mt.mat[0][0] * xi_new
1074 + mt.mat[0][1] * yj_new
1075 + mt.mat[0][2] * zk_new - vt.xyz[0] ;
1076
1077 fyj_base = mt.mat[1][0] * xi_new
1078 + mt.mat[1][1] * yj_new
1079 + mt.mat[1][2] * zk_new - vt.xyz[1] ;
1080
1081 fzk_base = mt.mat[2][0] * xi_new
1082 + mt.mat[2][1] * yj_new
1083 + mt.mat[2][2] * zk_new - vt.xyz[2] ;
1084
1085 dfxi_outer = mt.mat[0][0] ; /* outer loop is in x = 0 */
1086 dfyj_outer = mt.mat[1][0] ;
1087 dfzk_outer = mt.mat[2][0] ;
1088
1089 dfxi_inner = mt.mat[0][2] ; /* inner loop is in z = 2 */
1090 dfyj_inner = mt.mat[1][2] ;
1091 dfzk_inner = mt.mat[2][2] ;
1092
1093 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ;
1094 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ;
1095 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ;
1096
1097 switch( resam_mode ){
1098
1099 default:
1100 case RESAM_NN_TYPE:{
1101 float fxi_max , fyj_max , fzk_max ;
1102 float fxi_min , fyj_min , fzk_min ;
1103 float fxi_tmp , fyj_tmp , fzk_tmp ;
1104 int any_outside , all_outside ;
1105
1106 #ifdef SL_DEBUG
1107 if(PRINT_TRACING)
1108 { char str[256] ;
1109 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g",
1110 dfxi_inner,dfyj_inner,dfzk_inner,
1111 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; }
1112 #endif
1113 /** July 15, 1996:
1114 check if all the points are inside the old grid;
1115 if so, can use a version of the resampling loop
1116 that does not need to check each voxel for being
1117 inside -- hopefully, this will execute more quickly **/
1118
1119 FXYZTMP(xi_bot,yj_new,zk_bot) ;
1120 fxi_max = fxi_min = fxi_tmp ;
1121 fyj_max = fyj_min = fyj_tmp ;
1122 fzk_max = fzk_min = fzk_tmp ;
1123
1124 FXYZTMP(xi_top,yj_new,zk_bot) ;
1125 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1126 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1127 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1128
1129 FXYZTMP(xi_bot,yj_new,zk_top) ;
1130 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1131 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1132 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1133
1134 FXYZTMP(xi_top,yj_new,zk_top) ;
1135 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1136 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1137 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1138
1139 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) ||
1140 (fyj_min < fyj_bot) || (fyj_max > fyj_top) ||
1141 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ;
1142
1143 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) ||
1144 (fyj_max < fyj_bot) || (fyj_min > fyj_top) ||
1145 (fzk_max < fzk_bot) || (fzk_min > fzk_top)
1146 : 0 ;
1147
1148 #ifdef SL_DEBUG
1149 if(PRINT_TRACING)
1150 { char str[256] ;
1151 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max);
1152 STATUS(str) ;
1153 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max);
1154 STATUS(str) ;
1155 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max);
1156 STATUS(str) ; }
1157 #endif
1158
1159 /** redefine the macros specifying loop variables **/
1160
1161 #undef OUD_NAME
1162 #undef IND_NAME
1163 #undef OUD
1164 #undef OUD_bot
1165 #undef OUD_top
1166 #undef IND
1167 #undef IND_bot
1168 #undef IND_top
1169 #undef IND_nnn
1170
1171 #define OUD_NAME xi /* name of 2nd dimension of output image */
1172 #define IND_NAME zk /* name of 1st dimension of output image */
1173 #define IND_nnn nznew /* inner loop image dimension */
1174
1175 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */
1176 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */
1177 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */
1178 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */
1179 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */
1180 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */
1181
1182 if( all_outside ){
1183 #ifdef SL_DEBUG
1184 STATUS("NN resample has all outside") ;
1185 #endif
1186 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */
1187 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */
1188 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */
1189 bslice[out_ind++] = ZERO ;
1190 }
1191 }
1192 } else { /* at least some are inside */
1193
1194 int thz, tho , ob , ub=0 ;
1195
1196 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ;
1197
1198 fxi_top = nxold - 0.01 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */
1199 fyj_top = nyold - 0.01 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */
1200 fzk_top = nzold - 0.01 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */
1201 /* these variables now. */
1202
1203 /** thz = flag that indicates which of the steps df??_inner are zero.
1204 If two of them are zero, then the inner loop is parallel to
1205 one of the input brick axes, and so data may be pulled
1206 out in a very efficient fashion. In such a case, precompute
1207 the indexes for the inner loop:
1208
1209 the BLOOP macros load array ib, which holds the inner loop
1210 computed indexes for each inner loop position.
1211 ub = upper bound value for ib array value to still be
1212 inside input brick array.
1213 ob = outer loop index into input brick array (computed later) **/
1214
1215 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */
1216 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */
1217 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */
1218 else /* treatment if outer */
1219 thz = NONE_ZERO ; /* axes are special */
1220
1221 #ifdef SL_DEBUG
1222 if(PRINT_TRACING)
1223 { char str[256] ;
1224 if( any_outside ) sprintf(str,"NN resample has some outside: thz = %d",thz) ;
1225 else sprintf(str,"NN resample has all inside: thz = %d",thz) ;
1226 STATUS(str) ;
1227 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d",
1228 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ;
1229 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d",
1230 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; }
1231 #endif
1232
1233 switch(thz){
1234 case XY_ZERO:
1235 MAKE_IBIG(IND_top) ;
1236 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ;
1237 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; }
1238 break ;
1239
1240 case XZ_ZERO:
1241 MAKE_IBIG(IND_top) ;
1242 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ;
1243 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; }
1244 break ;
1245
1246 case YZ_ZERO:
1247 MAKE_IBIG(IND_top) ;
1248 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ;
1249 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; }
1250 break ;
1251 }
1252
1253 thz += OUTADD * any_outside ;
1254
1255 #ifdef SL_DEBUG
1256 STATUS("beginning NN outer loop") ;
1257 #endif
1258
1259 /*** outer loop ***/
1260
1261 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){
1262
1263 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */
1264 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */
1265 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */
1266
1267 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */
1268
1269 /*** There are 8 cases for the inner loop:
1270 all inside, inner loop not parallel to any input axis
1271 all inside, inner loop parallel to input brick z-axis
1272 all inside, inner loop parallel to input brick y-axis
1273 all inside, inner loop parallel to input brick x-axis
1274
1275 and then the 4 same cases repeated when not all desired
1276 points are inside the input brick. Each of these is
1277 coded separately for efficiency. This is important for
1278 rapid re-display of results during interactive imaging. ***/
1279
1280 switch(thz){
1281 case NONE_ZERO:
1282 case X_ZERO:
1283 case Y_ZERO:
1284 case Z_ZERO:
1285 case XYZ_ZERO:
1286 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; }
1287 break ;
1288
1289 case XY_ZERO:
1290 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
1291 ob = IBASE(xi_old,yj_old,0) ;
1292 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1293 break ;
1294
1295 case XZ_ZERO:
1296 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
1297 ob = IBASE(xi_old,0,zk_old) ;
1298 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1299 break ;
1300
1301 case YZ_ZERO:
1302 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
1303 ob = IBASE(0,yj_old,zk_old) ;
1304 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1305 break ;
1306
1307 default:
1308 case NONE_ZERO+OUTADD:
1309 case X_ZERO+OUTADD:
1310 case Y_ZERO+OUTADD:
1311 case Z_ZERO+OUTADD:
1312 case XYZ_ZERO+OUTADD:
1313 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; }
1314 break ;
1315
1316 case XY_ZERO+OUTADD:
1317 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
1318 if( TEST_OUT_XXX || TEST_OUT_YYY ){
1319 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
1320 } else {
1321 ob = IBASE(xi_old,yj_old,0) ;
1322 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
1323 }
1324 break ;
1325
1326 case XZ_ZERO+OUTADD:
1327 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
1328 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){
1329 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
1330 } else {
1331 ob = IBASE(xi_old,0,zk_old) ;
1332 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
1333 }
1334 break ;
1335
1336 case YZ_ZERO+OUTADD:
1337 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
1338 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){
1339 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
1340 } else {
1341 ob = IBASE(0,yj_old,zk_old) ;
1342 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
1343 }
1344 break ;
1345 }
1346 }
1347 }
1348 }
1349 break ; /* end of NN! */
1350
1351 case RESAM_BLOCK_TYPE:
1352 case RESAM_LINEAR_TYPE:{
1353 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ;
1354 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 ,
1355 f_k00 , f_kp1 , result ;
1356 float frac_xi , frac_yj , frac_zk ;
1357 int ibase ,
1358 in_jp1_k00 = jstep ,
1359 in_j00_kp1 = kstep ,
1360 in_jp1_kp1 = jstep+kstep ;
1361 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
1362
1363 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){
1364
1365 fxi_old = (fxi_base += dfxi_outer) ;
1366 fyj_old = (fyj_base += dfyj_outer) ;
1367 fzk_old = (fzk_base += dfzk_outer) ;
1368
1369 out_ind = zk_bot + xi_new * nznew ;
1370
1371 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){
1372
1373 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
1374
1375 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
1376 fyj_old < fyj_bot || fyj_old > fyj_top ||
1377 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
1378
1379 /***
1380 bslice[out_ind++] = 0 ;
1381 ***/
1382 FZERO(bslice[out_ind]) ; out_ind++ ;
1383 continue ;
1384 }
1385
1386 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
1387 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
1388 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
1389
1390 /* use NN if at edges of old brick */
1391
1392 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
1393 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
1394
1395 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
1396 ROUND(fyj_old) ,
1397 ROUND(fzk_old) ) ] ;
1398 continue ;
1399 }
1400
1401 /* compute weights for LINEAR interpolation in each direction */
1402
1403 if( resam_mode == RESAM_BLOCK_TYPE ){
1404 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ;
1405 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ;
1406 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ;
1407 } else {
1408 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
1409 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
1410 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
1411 }
1412
1413 /* interpolate in the x direction for each y & z */
1414
1415 ibase = IBASE(xi_old,yj_old,zk_old) ;
1416 /***
1417 f_j00_k00 = xwt_00 * bold[ibase]
1418 + xwt_p1 * bold[ibase+1] ;
1419
1420 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
1421 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
1422
1423 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
1424 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
1425
1426 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
1427 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
1428 ***/
1429 FMAD2( xwt_00 , bold[ibase] ,
1430 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
1431
1432 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
1433 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
1434
1435 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
1436 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
1437
1438 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
1439 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
1440
1441 /* interpolate in the y direction for each z */
1442 /***
1443 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
1444 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
1445 ***/
1446 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
1447 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
1448
1449 /* interpolate in the z direction and
1450 put the result into the output array!
1451 (the +0.5 is to force rounding of the result) */
1452 /***
1453 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
1454 ***/
1455 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
1456 FINAL( result , bslice[out_ind] ) ;
1457 out_ind++ ;
1458 }
1459 }
1460 }
1461 break ;
1462
1463 case RESAM_CUBIC_TYPE:{
1464 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */
1465 ywt_m1,ywt_00,ywt_p1,ywt_p2 ,
1466 zwt_m1,zwt_00,zwt_p1,zwt_p2 ;
1467
1468 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */
1469 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 ,
1470 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 ,
1471 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 ,
1472 f_km1 , f_k00 , f_kp1 , f_kp2 , result ;
1473 float frac_xi , frac_yj , frac_zk ;
1474
1475 int ibase , /* base index for interpolant */
1476 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */
1477 in_jm1_k00 = -jstep , /* 0 (00) */
1478 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */
1479 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */
1480 in_j00_km1 = - kstep , /* steps in j and k indices */
1481 in_j00_k00 = 0 ,
1482 in_j00_kp1 = kstep ,
1483 in_j00_kp2 = 2*kstep ,
1484 in_jp1_km1 = jstep- kstep ,
1485 in_jp1_k00 = jstep ,
1486 in_jp1_kp1 = jstep+ kstep ,
1487 in_jp1_kp2 = jstep+2*kstep , /* no '2*' [PT: 27 Jul 2021] */
1488 in_jp2_km1 =2*jstep- kstep ,
1489 in_jp2_k00 =2*jstep ,
1490 in_jp2_kp1 =2*jstep+ kstep ,
1491 in_jp2_kp2 =2*jstep+2*kstep ;
1492
1493 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
1494 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ;
1495
1496 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){
1497
1498 fxi_old = (fxi_base += dfxi_outer) ;
1499 fyj_old = (fyj_base += dfyj_outer) ;
1500 fzk_old = (fzk_base += dfzk_outer) ;
1501
1502 out_ind = zk_bot + xi_new * nznew ;
1503
1504 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){
1505
1506 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
1507
1508 /* check if outside old brick */
1509
1510 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
1511 fyj_old < fyj_bot || fyj_old > fyj_top ||
1512 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
1513 /***
1514 bslice[out_ind++] = 0 ;
1515 ***/
1516 FZERO(bslice[out_ind]) ; out_ind++ ;
1517 continue ;
1518 }
1519
1520 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
1521 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
1522 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
1523
1524 /* use NN if at very edges of old brick */
1525
1526 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
1527 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
1528
1529 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
1530 ROUND(fyj_old) ,
1531 ROUND(fzk_old) ) ] ;
1532 continue ;
1533 }
1534
1535 ibase = IBASE(xi_old,yj_old,zk_old) ;
1536
1537 /* use LINEAR if close to edges of old brick */
1538
1539 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 ||
1540 xi_old == 0 || yj_old == 0 || zk_old == 0 ){
1541
1542 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
1543 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
1544 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
1545 /***
1546 f_j00_k00 = xwt_00 * bold[ibase]
1547 + xwt_p1 * bold[ibase+1] ;
1548
1549 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
1550 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
1551
1552 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
1553 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
1554
1555 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
1556 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
1557 ***/
1558 FMAD2( xwt_00 , bold[ibase] ,
1559 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
1560
1561 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
1562 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
1563
1564 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
1565 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
1566
1567 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
1568 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
1569 /***
1570
1571 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
1572 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
1573 ***/
1574 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
1575 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
1576 /***
1577 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
1578 ***/
1579 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
1580 FINAL( result , bslice[out_ind] ) ;
1581 out_ind++ ;
1582 continue ;
1583 }
1584
1585 /* compute weights for CUBIC interpolation in each direction */
1586
1587 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ;
1588 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ;
1589
1590 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ;
1591 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ;
1592
1593 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ;
1594 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ;
1595
1596 /* interpolate in the x direction for each y & z */
1597
1598 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ;
1599 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ;
1600
1601 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ;
1602 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ;
1603
1604 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ;
1605 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ;
1606
1607 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ;
1608 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ;
1609
1610 /* interpolate in the y direction for each z */
1611 /***
1612 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1
1613 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ;
1614
1615 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00
1616 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ;
1617
1618 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1
1619 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ;
1620
1621 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2
1622 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ;
1623 ***/
1624 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 ,
1625 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ;
1626
1627 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 ,
1628 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ;
1629
1630 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 ,
1631 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ;
1632
1633 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 ,
1634 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ;
1635
1636 /* interpolate in the z direction and
1637 put the result into the output array!
1638 (the +0.5 is to force rounding of the result) */
1639 /***
1640 bslice[out_ind++] = 0.5 + CP_FACTOR
1641 * ( zwt_m1 * f_km1 + zwt_00 * f_k00
1642 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ;
1643 ***/
1644 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 ,
1645 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ;
1646 FSCAL(CP_FACTOR,result) ;
1647 FINAL( result , bslice[out_ind] ) ;
1648 out_ind++ ;
1649
1650 } /* end of inner loop */
1651 } /* end of outer loop */
1652 }
1653 break ;
1654 }
1655
1656 EXRETURN ;
1657 }
1658
1659 /*----------------------------------------------------------------------------*/
1660
1661 #undef DBG_ZZZ
1662
LMAP_ZNAME(THD_linear_mapping * map,int resam_mode,THD_dataxes * old_daxes,DTYPE * bold,THD_dataxes * new_daxes,int zk_fix,DTYPE * bslice)1663 void LMAP_ZNAME( THD_linear_mapping * map , int resam_mode ,
1664 THD_dataxes * old_daxes , DTYPE * bold ,
1665 THD_dataxes * new_daxes , int zk_fix , DTYPE * bslice )
1666 {
1667 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */
1668 THD_fvec3 vt = map->svec ;
1669
1670 int xi_new , yj_new , zk_new ; /* voxel indices in new */
1671 int xi_old , yj_old , zk_old ; /* voxel indices in old */
1672 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */
1673
1674 float fxi_top , fyj_top , fzk_top ;
1675 float fxi_bot , fyj_bot , fzk_bot ;
1676 float fxi_base , fyj_base , fzk_base ;
1677 float dfxi_outer , dfyj_outer , dfzk_outer ;
1678 float dfxi_inner , dfyj_inner , dfzk_inner ;
1679
1680 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */
1681 int out_ind , jstep , kstep ;
1682 int nxold,nyold,nzold , nxnew,nynew,nznew ;
1683
1684 #if 0
1685 ENTRY("AFNI_lmap_to_zslice") ;
1686 #else
1687 ENTRY( XSTRING(LMAP_ZNAME) ) ;
1688 #endif
1689
1690 #ifdef DBG_ZZZ
1691 INFO_message("+++++ Enter %s",XSTRING(LMAP_ZNAME)) ;
1692 ININFO_message(" zk_fix = %d",zk_fix) ;
1693 DUMP_MAT33(" + mt ",mt) ;
1694 DUMP_FVEC3(" + vt ",vt) ;
1695 DUMP_FVEC3(" + bot",map->bot) ;
1696 DUMP_FVEC3(" + top",map->top) ;
1697 #endif
1698
1699 /*--- set up ranges ---*/
1700
1701 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ;
1702 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ;
1703 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ;
1704
1705 if( zk_fix < zk_bot || zk_fix > zk_top ) EXRETURN ; /* map doesn't apply! */
1706
1707 if( bold==NULL || bslice==NULL ){ STATUS("NULL pointers?!"); EXRETURN; }
1708
1709 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ;
1710 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ;
1711 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ;
1712
1713 jstep = nxold ;
1714 kstep = nxold * nyold ;
1715
1716 /* set up base of indices in old */
1717
1718 xi_new = xi_bot-1 ;
1719 yj_new = yj_bot-1 ;
1720 zk_new = zk_fix ;
1721
1722 fxi_base = mt.mat[0][0] * xi_new
1723 + mt.mat[0][1] * yj_new
1724 + mt.mat[0][2] * zk_new - vt.xyz[0] ;
1725
1726 fyj_base = mt.mat[1][0] * xi_new
1727 + mt.mat[1][1] * yj_new
1728 + mt.mat[1][2] * zk_new - vt.xyz[1] ;
1729
1730 fzk_base = mt.mat[2][0] * xi_new
1731 + mt.mat[2][1] * yj_new
1732 + mt.mat[2][2] * zk_new - vt.xyz[2] ;
1733
1734 #ifdef DBG_ZZZ
1735 ININFO_message(" f??_base = %f %f %f",fxi_base,fyj_base,fzk_base) ;
1736 #endif
1737
1738 dfxi_outer = mt.mat[0][1] ; /* outer loop is in y = 1 */
1739 dfyj_outer = mt.mat[1][1] ;
1740 dfzk_outer = mt.mat[2][1] ;
1741
1742 #ifdef DBG_ZZZ
1743 ININFO_message(" df??_outer = %f %f %f",dfxi_outer,dfyj_outer,dfzk_outer) ;
1744 #endif
1745
1746 dfxi_inner = mt.mat[0][0] ; /* inner loop is in x = 0 */
1747 dfyj_inner = mt.mat[1][0] ;
1748 dfzk_inner = mt.mat[2][0] ;
1749
1750 #ifdef DBG_ZZZ
1751 ININFO_message(" df??_inner = %f %f %f",dfxi_inner,dfyj_inner,dfzk_inner) ;
1752 #endif
1753
1754 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ;
1755 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ;
1756 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ;
1757
1758 switch( resam_mode ){
1759
1760 default:
1761 case RESAM_NN_TYPE:{
1762 float fxi_max , fyj_max , fzk_max ;
1763 float fxi_min , fyj_min , fzk_min ;
1764 float fxi_tmp , fyj_tmp , fzk_tmp ;
1765 int any_outside , all_outside ;
1766
1767 #ifdef SL_DEBUG
1768 if(PRINT_TRACING)
1769 { char str[256] ;
1770 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g",
1771 dfxi_inner,dfyj_inner,dfzk_inner,
1772 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; }
1773 #endif
1774 /** July 15, 1996:
1775 check if all the points are inside the old grid;
1776 if so, can use a version of the resampling loop
1777 that does not need to check each voxel for being
1778 inside -- hopefully, this will execute more quickly **/
1779
1780 FXYZTMP(xi_bot,yj_bot,zk_new) ;
1781 fxi_max = fxi_min = fxi_tmp ;
1782 fyj_max = fyj_min = fyj_tmp ;
1783 fzk_max = fzk_min = fzk_tmp ;
1784
1785 FXYZTMP(xi_top,yj_bot,zk_new) ;
1786 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1787 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1788 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1789
1790 FXYZTMP(xi_bot,yj_top,zk_new) ;
1791 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1792 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1793 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1794
1795 FXYZTMP(xi_top,yj_top,zk_new) ;
1796 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ;
1797 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ;
1798 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ;
1799
1800 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) ||
1801 (fyj_min < fyj_bot) || (fyj_max > fyj_top) ||
1802 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ;
1803
1804 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) ||
1805 (fyj_max < fyj_bot) || (fyj_min > fyj_top) ||
1806 (fzk_max < fzk_bot) || (fzk_min > fzk_top)
1807 : 0 ;
1808
1809 #ifdef DBG_ZZZ
1810 ININFO_message(" any_outside=%d all_outside=%d",any_outside,all_outside) ;
1811 #endif
1812
1813 #ifdef SL_DEBUG
1814 if(PRINT_TRACING)
1815 { char str[256] ;
1816 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max);
1817 STATUS(str) ;
1818 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max);
1819 STATUS(str) ;
1820 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max);
1821 STATUS(str) ; }
1822 #endif
1823
1824 /** redefine the macros specifying loop variables **/
1825
1826 #undef OUD_NAME
1827 #undef IND_NAME
1828 #undef OUD
1829 #undef OUD_bot
1830 #undef OUD_top
1831 #undef IND
1832 #undef IND_bot
1833 #undef IND_top
1834 #undef IND_nnn
1835
1836 #define OUD_NAME yj /* name of 2nd dimension of output image */
1837 #define IND_NAME xi /* name of 1st dimension of output image */
1838 #define IND_nnn nxnew /* inner loop image dimension */
1839
1840 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */
1841 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */
1842 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */
1843 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */
1844 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */
1845 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */
1846
1847 if( all_outside ){
1848 #ifdef SL_DEBUG
1849 STATUS("NN resample has all outside") ;
1850 #endif
1851 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */
1852 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */
1853 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */
1854 bslice[out_ind++] = ZERO ;
1855 }
1856 }
1857 } else { /* at least some are inside */
1858
1859 int thz , tho , ob , ub=0 ;
1860
1861 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ;
1862
1863 #ifdef DBG_ZZZ
1864 ININFO_message(" new f??_base = %f %f %f",fxi_base,fyj_base,fzk_base) ;
1865 #endif
1866
1867 fxi_top = nxold - 0.01 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */
1868 fyj_top = nyold - 0.01 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */
1869 fzk_top = nzold - 0.01 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */
1870 /* these variables now. */
1871
1872 /** thz = flag that indicates which of the steps df??_inner are zero.
1873 If two of them are zero, then the inner loop is parallel to
1874 one of the input brick axes, and so data may be pulled
1875 out in a very efficient fashion. In such a case, precompute
1876 the indexes for the inner loop:
1877
1878 the BLOOP macros load array ib, which holds the inner loop
1879 computed indexes for each inner loop position.
1880 ub = upper bound value for ib array value to still be
1881 inside input brick array.
1882 ob = outer loop index into input brick array (computed later) **/
1883
1884 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */
1885 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */
1886 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */
1887 else /* treatment if outer */
1888 thz = NONE_ZERO ; /* axes are special */
1889
1890 #ifdef SL_DEBUG
1891 if(PRINT_TRACING)
1892 { char str[256] ;
1893 if( any_outside ) sprintf(str,"NN resample has some outside: thz = %d",thz) ;
1894 else sprintf(str,"NN resample has all inside: thz = %d",thz) ;
1895 STATUS(str) ;
1896 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d",
1897 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ;
1898 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d",
1899 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; }
1900 #endif
1901
1902 switch(thz){
1903 case XY_ZERO:
1904 MAKE_IBIG(IND_top) ;
1905 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ;
1906 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; }
1907 break ;
1908
1909 case XZ_ZERO:
1910 MAKE_IBIG(IND_top) ;
1911 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ;
1912 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; }
1913 break ;
1914
1915 case YZ_ZERO:
1916 MAKE_IBIG(IND_top) ;
1917 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ;
1918 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; }
1919 break ;
1920 }
1921
1922 thz += OUTADD * any_outside ;
1923
1924 #ifdef SL_DEBUG
1925 STATUS("beginning NN outer loop") ;
1926 #endif
1927
1928 /*** outer loop ***/
1929
1930 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){
1931
1932 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */
1933 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */
1934 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */
1935
1936 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */
1937
1938 /*** There are 8 cases for the inner loop:
1939 all inside, inner loop not parallel to any input axis
1940 all inside, inner loop parallel to input brick z-axis
1941 all inside, inner loop parallel to input brick y-axis
1942 all inside, inner loop parallel to input brick x-axis
1943
1944 and then the 4 same cases repeated when not all desired
1945 points are inside the input brick. Each of these is
1946 coded separately for efficiency. This is important for
1947 rapid re-display of results during interactive imaging. ***/
1948
1949 switch(thz){
1950 case NONE_ZERO:
1951 case X_ZERO:
1952 case Y_ZERO:
1953 case Z_ZERO:
1954 case XYZ_ZERO:
1955 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; }
1956 break ;
1957
1958 case XY_ZERO:
1959 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
1960 ob = IBASE(xi_old,yj_old,0) ;
1961 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1962 break ;
1963
1964 case XZ_ZERO:
1965 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
1966 ob = IBASE(xi_old,0,zk_old) ;
1967 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1968 break ;
1969
1970 case YZ_ZERO:
1971 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
1972 ob = IBASE(0,yj_old,zk_old) ;
1973 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; }
1974 break ;
1975
1976 default:
1977 case NONE_ZERO+OUTADD:
1978 case X_ZERO+OUTADD:
1979 case Y_ZERO+OUTADD:
1980 case Z_ZERO+OUTADD:
1981 case XYZ_ZERO+OUTADD:
1982 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; }
1983 break ;
1984
1985 case XY_ZERO+OUTADD:
1986 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ;
1987 if( TEST_OUT_XXX || TEST_OUT_YYY ){
1988 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
1989 } else {
1990 ob = IBASE(xi_old,yj_old,0) ;
1991 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
1992 }
1993 break ;
1994
1995 case XZ_ZERO+OUTADD:
1996 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ;
1997 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){
1998 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
1999 } else {
2000 ob = IBASE(xi_old,0,zk_old) ;
2001 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
2002 }
2003 break ;
2004
2005 case YZ_ZERO+OUTADD:
2006 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ;
2007 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){
2008 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; }
2009 } else {
2010 ob = IBASE(0,yj_old,zk_old) ;
2011 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; }
2012 }
2013 break ;
2014 }
2015 }
2016 }
2017 }
2018 break ; /* end of NN! */
2019
2020 case RESAM_BLOCK_TYPE:
2021 case RESAM_LINEAR_TYPE:{
2022 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ;
2023 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 ,
2024 f_k00 , f_kp1 , result ;
2025 float frac_xi , frac_yj , frac_zk ;
2026 int ibase ,
2027 in_jp1_k00 = jstep ,
2028 in_j00_kp1 = kstep ,
2029 in_jp1_kp1 = jstep+kstep ;
2030 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
2031
2032 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){
2033
2034 fxi_old = (fxi_base += dfxi_outer) ;
2035 fyj_old = (fyj_base += dfyj_outer) ;
2036 fzk_old = (fzk_base += dfzk_outer) ;
2037
2038 out_ind = xi_bot + yj_new * nxnew ;
2039
2040 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){
2041
2042 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
2043
2044 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
2045 fyj_old < fyj_bot || fyj_old > fyj_top ||
2046 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
2047
2048 /***
2049 bslice[out_ind++] = 0 ;
2050 ***/
2051 FZERO(bslice[out_ind]) ; out_ind++ ;
2052 continue ;
2053 }
2054
2055 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
2056 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
2057 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
2058
2059 /* use NN if at edges of old brick */
2060
2061 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
2062 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
2063
2064 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
2065 ROUND(fyj_old) ,
2066 ROUND(fzk_old) ) ] ;
2067 continue ;
2068 }
2069
2070 /* compute weights for LINEAR interpolation in each direction */
2071
2072 if( resam_mode == RESAM_BLOCK_TYPE ){
2073 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ;
2074 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ;
2075 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ;
2076 } else {
2077 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
2078 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
2079 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
2080 }
2081
2082 /* interpolate in the x direction for each y & z */
2083
2084 ibase = IBASE(xi_old,yj_old,zk_old) ;
2085 /***
2086 f_j00_k00 = xwt_00 * bold[ibase]
2087 + xwt_p1 * bold[ibase+1] ;
2088
2089 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
2090 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
2091
2092 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
2093 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
2094
2095 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
2096 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
2097 ***/
2098 FMAD2( xwt_00 , bold[ibase] ,
2099 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
2100
2101 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
2102 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
2103
2104 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
2105 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
2106
2107 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
2108 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
2109
2110 /* interpolate in the y direction for each z */
2111 /***
2112 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
2113 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
2114 ***/
2115 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
2116 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
2117
2118 /* interpolate in the z direction and
2119 put the result into the output array!
2120 (the +0.5 is to force rounding of the result) */
2121 /***
2122 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
2123 ***/
2124 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
2125 FINAL( result , bslice[out_ind] ) ;
2126 out_ind++ ;
2127
2128 }
2129 }
2130 }
2131 break ;
2132
2133 case RESAM_CUBIC_TYPE:{
2134 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */
2135 ywt_m1,ywt_00,ywt_p1,ywt_p2 ,
2136 zwt_m1,zwt_00,zwt_p1,zwt_p2 ;
2137
2138 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */
2139 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 ,
2140 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 ,
2141 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 ,
2142 f_km1 , f_k00 , f_kp1 , f_kp2 , result ;
2143 float frac_xi , frac_yj , frac_zk ;
2144
2145 int ibase , /* base index for interpolant */
2146 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */
2147 in_jm1_k00 = -jstep , /* 0 (00) */
2148 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */
2149 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */
2150 in_j00_km1 = - kstep , /* steps in j and k indices */
2151 in_j00_k00 = 0 ,
2152 in_j00_kp1 = kstep ,
2153 in_j00_kp2 = 2*kstep ,
2154 in_jp1_km1 = jstep- kstep ,
2155 in_jp1_k00 = jstep ,
2156 in_jp1_kp1 = jstep+ kstep ,
2157 in_jp1_kp2 = jstep+2*kstep , /* no '2*' [27 Jul 2021 rickr] */
2158 in_jp2_km1 =2*jstep- kstep ,
2159 in_jp2_k00 =2*jstep ,
2160 in_jp2_kp1 =2*jstep+ kstep ,
2161 in_jp2_kp2 =2*jstep+2*kstep ;
2162
2163 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ;
2164 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ;
2165
2166 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){
2167
2168 fxi_old = (fxi_base += dfxi_outer) ;
2169 fyj_old = (fyj_base += dfyj_outer) ;
2170 fzk_old = (fzk_base += dfzk_outer) ;
2171
2172 out_ind = xi_bot + yj_new * nxnew ;
2173
2174 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){
2175
2176 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ;
2177
2178 /* check if outside old brick */
2179
2180 if( fxi_old < fxi_bot || fxi_old > fxi_top ||
2181 fyj_old < fyj_bot || fyj_old > fyj_top ||
2182 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */
2183 /***
2184 bslice[out_ind++] = 0 ;
2185 ***/
2186 FZERO(bslice[out_ind]) ; out_ind++ ;
2187 continue ;
2188 }
2189
2190 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ;
2191 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ;
2192 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ;
2193
2194 /* use NN if at very edges of old brick */
2195
2196 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 ||
2197 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){
2198
2199 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) ,
2200 ROUND(fyj_old) ,
2201 ROUND(fzk_old) ) ] ;
2202 continue ;
2203 }
2204
2205 ibase = IBASE(xi_old,yj_old,zk_old) ;
2206
2207 /* use LINEAR if close to edges of old brick */
2208
2209 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 ||
2210 xi_old == 0 || yj_old == 0 || zk_old == 0 ){
2211
2212 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ;
2213 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ;
2214 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ;
2215 /***
2216 f_j00_k00 = xwt_00 * bold[ibase]
2217 + xwt_p1 * bold[ibase+1] ;
2218
2219 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00]
2220 + xwt_p1 * bold[ibase+in_jp1_k00+1] ;
2221
2222 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1]
2223 + xwt_p1 * bold[ibase+in_j00_kp1+1] ;
2224
2225 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1]
2226 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ;
2227
2228 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ;
2229 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ;
2230
2231 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ;
2232 ***/
2233 FMAD2( xwt_00 , bold[ibase] ,
2234 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ;
2235
2236 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] ,
2237 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ;
2238
2239 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] ,
2240 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ;
2241
2242 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] ,
2243 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ;
2244
2245 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ;
2246 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ;
2247
2248 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ;
2249 FINAL( result , bslice[out_ind] ) ;
2250 out_ind++ ;
2251 continue ;
2252 }
2253
2254 /* compute weights for CUBIC interpolation in each direction */
2255
2256 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ;
2257 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ;
2258
2259 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ;
2260 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ;
2261
2262 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ;
2263 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ;
2264
2265 /* interpolate in the x direction for each y & z */
2266
2267 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ;
2268 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ;
2269
2270 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ;
2271 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ;
2272
2273 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ;
2274 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ;
2275
2276 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ;
2277 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ;
2278
2279 /* interpolate in the y direction for each z */
2280 /***
2281 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1
2282 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ;
2283
2284 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00
2285 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ;
2286
2287 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1
2288 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ;
2289
2290 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2
2291 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ;
2292 ***/
2293 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 ,
2294 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ;
2295
2296 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 ,
2297 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ;
2298
2299 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 ,
2300 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ;
2301
2302 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 ,
2303 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ;
2304
2305 /* interpolate in the z direction and
2306 put the result into the output array!
2307 (the +0.5 is to force rounding of the result) */
2308 /***
2309 bslice[out_ind++] = 0.5 + CP_FACTOR
2310 * ( zwt_m1 * f_km1 + zwt_00 * f_k00
2311 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ;
2312 ***/
2313 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 ,
2314 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ;
2315 FSCAL(CP_FACTOR,result) ;
2316 FINAL( result , bslice[out_ind] ) ;
2317 out_ind++ ;
2318
2319 } /* end of inner loop */
2320 } /* end of outer loop */
2321 }
2322 break ;
2323
2324 }
2325
2326 EXRETURN ;
2327 }
2328
2329 /*---------------------------------------------------------------------
2330 Routine to copy data from a brick array into a slice:
2331 nxx,nyy,nzz = brick dimensions
2332 fixed_axis = 1, 2, or 3 (x, y, or z)
2333 fixed_index = subscript of chosen slice in the fixed_axis direction
2334 bold = pointer to brick array
2335 bslice = pointer to area to get slice from brick array
2336 fixed_axis = 1 --> nyy * nzz
2337 2 --> nzz * nxx
2338 3 --> nxx * nyy
2339 -----------------------------------------------------------------------*/
2340
B2SL_NAME(int nxx,int nyy,int nzz,int fixed_axis,int fixed_index,DTYPE * bold,DTYPE * bslice)2341 void B2SL_NAME( int nxx, int nyy, int nzz ,
2342 int fixed_axis , int fixed_index , DTYPE * bold , DTYPE * bslice )
2343 {
2344 int ystep = nxx , zstep = nxx*nyy ;
2345
2346 #if 0
2347 ENTRY("AFNI_br2sl") ;
2348 #else
2349 ENTRY(XSTRING(B2SL_NAME)) ;
2350 #endif
2351
2352 if( bold == NULL || bslice == NULL ) EXRETURN ;
2353
2354 switch( fixed_axis ){
2355
2356 case 1:{
2357 register int out , base , yy,zz , inn ;
2358
2359 out = 0 ;
2360 base = fixed_index ;
2361 for( zz=0 ; zz < nzz ; zz++ ){
2362 inn = base ;
2363 base += zstep ;
2364 for( yy=0 ; yy < nyy ; yy++ ){
2365 bslice[out++] = bold[inn] ; inn += ystep ;
2366 }
2367 }
2368 }
2369 break ;
2370
2371 case 2:{
2372 register int out , base , xx,zz , inn ;
2373
2374 out = 0 ;
2375 base = fixed_index * ystep ;
2376 for( xx=0 ; xx < nxx ; xx++ ){
2377 inn = base ;
2378 base += 1 ;
2379 for( zz=0 ; zz < nzz ; zz++ ){
2380 bslice[out++] = bold[inn] ; inn += zstep ;
2381 }
2382 }
2383 }
2384 break ;
2385
2386 case 3:{
2387 register int out , base , xx,yy , inn ;
2388
2389 base = fixed_index * zstep ;
2390 #if 0
2391 out = 0 ;
2392 for( yy=0 ; yy < nyy ; yy++ ){
2393 inn = base ;
2394 base += ystep ;
2395 for( xx=0 ; xx < nxx ; xx++ )
2396 bslice[out++] = bold[inn++] ;
2397 }
2398 #else
2399 (void) memcpy( bslice , bold+base , sizeof(DTYPE) * zstep ) ;
2400 #endif
2401 }
2402 break ;
2403 }
2404
2405 EXRETURN ;
2406 }
2407