Doxygen Source Code Documentation
afni_slice.c File Reference
#include "afni_warp.h"
Go to the source code of this file.
Defines | |
#define | LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE) |
#define | LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE) |
#define | LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE) |
#define | B2SL_NAME TWO_TWO(AFNI_br2sl_,DTYPE) |
#define | FMAD2_short(a, d1, b, d2, e) (e)=(a)*(d1)+(b)*(d2) |
#define | FMAD2_float FMAD2_short |
#define | FMAD2_byte FMAD2_short |
#define | FMAD2_int FMAD2_short |
#define | FMAD2_double FMAD2_short |
#define | FMAD2_complex(a, d1, b, d2, e) |
#define | FMAD2_rgbyte(a, d1, bb, d2, e) |
#define | FMAD2 TWO_TWO(FMAD2_,DTYPE) |
#define | FMAD4_short(a, d1, b, d2, c, d3, d, d4, e) (e)=(a)*(d1)+(b)*(d2)+(c)*(d3)+(d)*(d4) |
#define | FMAD4_float FMAD4_short |
#define | FMAD4_byte FMAD4_short |
#define | FMAD4_int FMAD4_short |
#define | FMAD4_double FMAD4_short |
#define | FMAD4_complex(a, d1, b, d2, c, d3, d, d4, e) |
#define | FMAD4_rgbyte(a, d1, bb, d2, c, d3, d, d4, e) |
#define | FMAD4 TWO_TWO(FMAD4_,DTYPE) |
#define | FSCAL_short(a, b) (b)*=(a) |
#define | FSCAL_float FSCAL_short |
#define | FSCAL_byte FSCAL_short |
#define | FSCAL_int FSCAL_short |
#define | FSCAL_double FSCAL_short |
#define | FSCAL_complex(a, b) ( (b).r *= (a) , (b).i *= (a) ) |
#define | FSCAL_rgbyte(a, bb) ( (bb).r*= (a) , (bb).g*= (a) , (bb).b*= (a) ) |
#define | FSCAL TWO_TWO(FSCAL_,DTYPE) |
#define | CLIP_OVERFLOW |
#define | FINAL_short(a, b) |
#define | FINAL_byte(a, b) |
#define | FINAL_int(a, b) (b)=((int)((a)+0.5)) |
#define | FINAL_float(a, b) (b)=(a) |
#define | FINAL_double FINAL_float |
#define | FINAL_complex FINAL_float |
#define | FINAL_rgbyte FINAL_float |
#define | FINAL TWO_TWO(FINAL_,DTYPE) |
#define | FZERO_short(b) (b)=0 |
#define | FZERO_byte FZERO_short |
#define | FZERO_int FZERO_short |
#define | FZERO_float(b) (b)=0.0 |
#define | FZERO_double FZERO_float |
#define | FZERO_complex(b) ( (b).r = 0.0 , (b).i = 0.0 ) |
#define | FZERO_rgbyte(bb) ( (bb).r=(bb).g=(bb).g = 0 ) |
#define | FZERO TWO_TWO(FZERO_,DTYPE) |
#define | ZERO_short 0 |
#define | ZERO_byte 0 |
#define | ZERO_int 0 |
#define | ZERO_float 0.0 |
#define | ZERO_double 0.0 |
#define | ZERO_complex complex_zero |
#define | ZERO_rgbyte rgbyte_zero |
#define | ZERO TWO_TWO(ZERO_,DTYPE) |
#define | INTYPE_short float |
#define | INTYPE_float float |
#define | INTYPE_byte float |
#define | INTYPE_int float |
#define | INTYPE_double double |
#define | INTYPE_complex complex |
#define | INTYPE_rgbyte rgbyte |
#define | INTYPE TWO_TWO(INTYPE_,DTYPE) |
#define | ZZZ 1.e-5 |
#define | NONE_ZERO 0 |
#define | X_ZERO 1 |
#define | Y_ZERO 2 |
#define | XY_ZERO 3 |
#define | Z_ZERO 4 |
#define | XZ_ZERO 5 |
#define | YZ_ZERO 6 |
#define | XYZ_ZERO 7 |
#define | OUTADD 100 |
#define | THREEZ(x, y, z) ((fabs(x)<ZZZ) + 2*(fabs(y)<ZZZ) + 4*(fabs(z)<ZZZ)) |
#define | NN_ALOOP_GEN |
#define | NN_ALOOP_PAR(ijk) (bslice[out_ind++] = bold[ ib[ijk]+ob ]) |
#define | NN_BLOOP_XY_ZERO(ijk) (fzk_old += dfzk_inner , zk_old = GFLOOR(fzk_old) , ib[ijk] = IBASE(0,0,zk_old)) |
#define | NN_BLOOP_XZ_ZERO(ijk) (fyj_old += dfyj_inner , yj_old = GFLOOR(fyj_old) , ib[ijk] = IBASE(0,yj_old,0)) |
#define | NN_BLOOP_YZ_ZERO(ijk) (fxi_old += dfxi_inner , xi_old = GFLOOR(fxi_old) , ib[ijk] = IBASE(xi_old,0,0)) |
#define | TEST_OUT_XXX ( fxi_old < fxi_bot || fxi_old > fxi_top ) |
#define | TEST_OUT_YYY ( fyj_old < fyj_bot || fyj_old > fyj_top ) |
#define | TEST_OUT_ZZZ ( fzk_old < fzk_bot || fzk_old > fzk_top ) |
#define | TEST_OUT_ALL (TEST_OUT_XXX || TEST_OUT_YYY || TEST_OUT_ZZZ) |
#define | NN_CLOOP_GEN |
#define | NN_CLOOP_PAR(ijk) (bslice[out_ind++] = (ib[ijk]<0 || ib[ijk]>=ub) ? ZERO : bold[ib[ijk]+ob]) |
#define | NN_ZLOOP (bslice[out_ind++] = ZERO) |
#define | MAKE_IBIG(top) |
#define | IBASE(i, j, k) ((i)+(j)*jstep+(k)*kstep) |
#define | ROUND(qq) ((int)(qq+0.5)) |
#define | FLOOR(qq) ((int)(qq)) |
#define | GFLOOR(qq) ((int)floor(qq)) |
#define | LP_00(x) (1.0-(x)) |
#define | LP_P1(x) (x) |
#define | BP_hh(x) (8.0*((x)*(x))*((x)*(x))) |
#define | BP_00(x) ( ((x)<0.5) ? (1-BP_hh(x)) : ( BP_hh(1-(x))) ) |
#define | BP_P1(x) ( ((x)<0.5) ? ( BP_hh(x)) : (1-BP_hh(1-(x))) ) |
#define | CP_M1(x) (-(x)*((x)-1)*((x)-2)) |
#define | CP_00(x) (3.0*((x)+1)*((x)-1)*((x)-2)) |
#define | CP_P1(x) (-3.0*(x)*((x)+1)*((x)-2)) |
#define | CP_P2(x) ((x)*((x)+1)*((x)-1)) |
#define | CP_FACTOR 4.62962963e-3 |
#define | FXYZTMP(xx, yy, zz) |
#define | OUD_NAME zk |
#define | IND_NAME yj |
#define | IND_nnn nynew |
#define | OUD TWO_TWO(OUD_NAME , _new) |
#define | OUD_bot TWO_TWO(OUD_NAME , _bot) |
#define | OUD_top TWO_TWO(OUD_NAME , _top) |
#define | IND TWO_TWO(IND_NAME , _new) |
#define | IND_bot TWO_TWO(IND_NAME , _bot) |
#define | IND_top TWO_TWO(IND_NAME , _top) |
#define | CXINT(j, k, ff) |
#define | OUD_NAME xi |
#define | IND_NAME zk |
#define | IND_nnn nznew |
#define | OUD TWO_TWO(OUD_NAME , _new) |
#define | OUD_bot TWO_TWO(OUD_NAME , _bot) |
#define | OUD_top TWO_TWO(OUD_NAME , _top) |
#define | IND TWO_TWO(IND_NAME , _new) |
#define | IND_bot TWO_TWO(IND_NAME , _bot) |
#define | IND_top TWO_TWO(IND_NAME , _top) |
#define | OUD_NAME yj |
#define | IND_NAME xi |
#define | IND_nnn nxnew |
#define | OUD TWO_TWO(OUD_NAME , _new) |
#define | OUD_bot TWO_TWO(OUD_NAME , _bot) |
#define | OUD_top TWO_TWO(OUD_NAME , _top) |
#define | IND TWO_TWO(IND_NAME , _new) |
#define | IND_bot TWO_TWO(IND_NAME , _bot) |
#define | IND_top TWO_TWO(IND_NAME , _top) |
Functions | |
void | LMAP_XNAME (THD_linear_mapping *map, int resam_mode, THD_dataxes *old_daxes, DTYPE *bold, THD_dataxes *new_daxes, int xi_fix, DTYPE *bslice) |
void | LMAP_YNAME (THD_linear_mapping *map, int resam_mode, THD_dataxes *old_daxes, DTYPE *bold, THD_dataxes *new_daxes, int yj_fix, DTYPE *bslice) |
void | LMAP_ZNAME (THD_linear_mapping *map, int resam_mode, THD_dataxes *old_daxes, DTYPE *bold, THD_dataxes *new_daxes, int zk_fix, DTYPE *bslice) |
void | B2SL_NAME (int nxx, int nyy, int nzz, int fixed_axis, int fixed_index, DTYPE *bold, DTYPE *bslice) |
Variables | |
complex | complex_zero = { 0.0,0.0 } |
rgbyte | rgbyte_zero = { 0,0,0 } |
int * | ib = NULL |
int | nib = -1 |
Define Documentation
|
Definition at line 44 of file afni_slice.c. |
|
Definition at line 273 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 272 of file afni_slice.c. |
|
Definition at line 274 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
macros for assigning final result from INTYPE a to DTYPE b * Definition at line 103 of file afni_slice.c. |
|
Definition at line 279 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 282 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 278 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 280 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 281 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Value: FMAD4( xwt_m1 , bold[ibase + in_j ## j ## _k ## k -1 ] , \ xwt_00 , bold[ibase + in_j ## j ## _k ## k ] , \ xwt_p1 , bold[ibase + in_j ## j ## _k ## k +1 ] , \ xwt_p2 , bold[ibase + in_j ## j ## _k ## k +2 ] , ff ) |
|
Definition at line 122 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Value: (b) = ( ((a)< 0.0) ? (0) : \ ((a)> 255.0) ? (255) : ((byte)((a)+0.5)) ) Definition at line 109 of file afni_slice.c. |
|
Definition at line 120 of file afni_slice.c. |
|
Definition at line 119 of file afni_slice.c. |
|
Definition at line 118 of file afni_slice.c. |
|
Definition at line 117 of file afni_slice.c. |
|
Definition at line 121 of file afni_slice.c. |
|
Value: (b) = ( ((a)<-32767.0) ? (-32767) : \ ((a)> 32767.0) ? ( 32767) : ((short)((a)+0.5)) ) Definition at line 106 of file afni_slice.c. |
|
Definition at line 262 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 66 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 56 of file afni_slice.c. |
|
Value: ( (e).r = (a)*(d1).r + (b)*(d2).r, \ (e).i = (a)*(d1).i + (b)*(d2).i ) Definition at line 60 of file afni_slice.c. |
|
Definition at line 58 of file afni_slice.c. |
|
Definition at line 55 of file afni_slice.c. |
|
Definition at line 57 of file afni_slice.c. |
|
Value: ( (e).r = (a)*(d1).r + (bb)*(d2).r, \ (e).g = (a)*(d1).g + (bb)*(d2).g, \ (e).b = (a)*(d1).b + (bb)*(d2).b ) Definition at line 63 of file afni_slice.c. |
|
macros for e = a*d1 + b*d2 (a,b floats; d1,d2 DTYPEs) * Definition at line 54 of file afni_slice.c. |
|
Definition at line 85 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 72 of file afni_slice.c. |
|
Value: ( (e).r = (a)*(d1).r + (b)*(d2).r + (c)*(d3).r + (d)*(d4).r, \ (e).i = (a)*(d1).i + (b)*(d2).i + (c)*(d3).i + (d)*(d4).i ) Definition at line 76 of file afni_slice.c. |
|
Definition at line 74 of file afni_slice.c. |
|
Definition at line 71 of file afni_slice.c. |
|
Definition at line 73 of file afni_slice.c. |
|
Value: ( (e).r = (a)*(d1).r + (bb)*(d2).r + (c)*(d3).r + (d)*(d4).r, \ (e).g = (a)*(d1).g + (bb)*(d2).g + (c)*(d3).g + (d)*(d4).g, \ (e).b = (a)*(d1).b + (bb)*(d2).b + (c)*(d3).b + (d)*(d4).b ) Definition at line 80 of file afni_slice.c. |
|
macros for e = a*d1 + b*d2 + c*d3 + d*d3 (a-d floats; d1-d4 DTYPEs) * Definition at line 70 of file afni_slice.c. |
|
Definition at line 97 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 91 of file afni_slice.c. |
|
Definition at line 94 of file afni_slice.c. |
|
Definition at line 93 of file afni_slice.c. |
|
Definition at line 90 of file afni_slice.c. |
|
Definition at line 92 of file afni_slice.c. |
|
Definition at line 96 of file afni_slice.c. |
|
macros to multiply float a times DTYPE b and store the result in b again * Definition at line 89 of file afni_slice.c. |
|
Value: ( fxi_tmp = mt.mat[0][0] * xx \ + mt.mat[0][1] * yy \ + mt.mat[0][2] * zz - vt.xyz[0] , \ fyj_tmp = mt.mat[1][0] * xx \ + mt.mat[1][1] * yy \ + mt.mat[1][2] * zz - vt.xyz[1] , \ fzk_tmp = mt.mat[2][0] * xx \ + mt.mat[2][1] * yy \ + mt.mat[2][2] * zz - vt.xyz[2] ) Definition at line 284 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 133 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 127 of file afni_slice.c. |
|
Definition at line 131 of file afni_slice.c. |
|
Definition at line 130 of file afni_slice.c. |
|
Definition at line 129 of file afni_slice.c. |
|
Definition at line 128 of file afni_slice.c. |
|
Definition at line 132 of file afni_slice.c. |
|
macros for putting a zero into DTYPE b * Definition at line 126 of file afni_slice.c. |
|
Definition at line 263 of file afni_slice.c. |
|
Definition at line 260 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 159 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 154 of file afni_slice.c. |
|
Definition at line 157 of file afni_slice.c. |
|
Definition at line 156 of file afni_slice.c. |
|
Definition at line 153 of file afni_slice.c. |
|
Definition at line 155 of file afni_slice.c. |
|
Definition at line 158 of file afni_slice.c. |
|
macros for intermediate interpolants data type * Definition at line 152 of file afni_slice.c. |
|
macros for function names defined in this file * Definition at line 41 of file afni_slice.c. Referenced by AFNI_dataset_slice(). |
|
Definition at line 42 of file afni_slice.c. Referenced by AFNI_dataset_slice(). |
|
Definition at line 43 of file afni_slice.c. Referenced by AFNI_dataset_slice(). |
|
Definition at line 267 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 268 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Value: do{ if(nib < (top)){ \ if(ib != NULL) free(ib) ; \ ib = (int *) malloc(sizeof(int)*((top)+9)) ; \ if(ib==NULL){ \ fprintf(stderr,"\nmalloc fails in NN reslice!\n");EXIT(1);} \ nib = (top) ; } } while(0) Definition at line 233 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Value: (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \ xi_old = FLOOR(fxi_old) , yj_old = FLOOR(fyj_old) , zk_old = FLOOR(fzk_old) , \ bslice[out_ind++] = bold[ IBASE(xi_old,yj_old,zk_old) ]) Definition at line 183 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 188 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
BLOOP: assign values to the ib array that will be used to index into the input brick for rapid access; there is one BLOOP for each possible parallel axis * Definition at line 194 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 197 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 200 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Value: (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \ bslice[out_ind++] = (TEST_OUT_ALL) ? ZERO : \ bold[IBASE(FLOOR(fxi_old),FLOOR(fyj_old),FLOOR(fzk_old))] ) Definition at line 214 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 219 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
ZLOOP: just assign zero to each output * Definition at line 224 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 167 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 175 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 261 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 209 of file afni_slice.c. |
|
macros to test if the point we want to resample from is outside the old array * Definition at line 205 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 206 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 207 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 176 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 168 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 170 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 174 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 172 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 169 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 173 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 171 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 148 of file afni_slice.c. Referenced by LMAP_XNAME(), LMAP_YNAME(), and LMAP_ZNAME(). |
|
Definition at line 142 of file afni_slice.c. |
|
Definition at line 146 of file afni_slice.c. |
|
Definition at line 145 of file afni_slice.c. |
|
Definition at line 144 of file afni_slice.c. |
|
Definition at line 143 of file afni_slice.c. |
|
Definition at line 147 of file afni_slice.c. |
|
Definition at line 141 of file afni_slice.c. |
|
test and flags for which steps are zero * Definition at line 166 of file afni_slice.c. |
Function Documentation
|
thz = flag that indicates which of the steps df??_inner are zero. If two of them are zero, then the inner loop is parallel to one of the input brick axes, and so data may be pulled out in a very efficient fashion. In such a case, precompute the indexes for the inner loop: the BLOOP macros load array ib, which holds the inner loop computed indexes for each inner loop position. ub = upper bound value for ib array value to still be inside input brick array. ob = outer loop index into input brick array (computed later) * Definition at line 2203 of file afni_slice.c.
02205 { 02206 int ystep = nxx , zstep = nxx*nyy ; 02207 02208 ENTRY("AFNI_br2sl") ; 02209 02210 switch( fixed_axis ){ 02211 02212 case 1:{ 02213 register int out , base , yy,zz , inn ; 02214 02215 out = 0 ; 02216 base = fixed_index ; 02217 for( zz=0 ; zz < nzz ; zz++ ){ 02218 inn = base ; 02219 base += zstep ; 02220 for( yy=0 ; yy < nyy ; yy++ ){ 02221 bslice[out++] = bold[inn] ; inn += ystep ; 02222 } 02223 } 02224 } 02225 break ; 02226 02227 case 2:{ 02228 register int out , base , xx,zz , inn ; 02229 02230 out = 0 ; 02231 base = fixed_index * ystep ; 02232 for( xx=0 ; xx < nxx ; xx++ ){ 02233 inn = base ; 02234 base += 1 ; 02235 for( zz=0 ; zz < nzz ; zz++ ){ 02236 bslice[out++] = bold[inn] ; inn += zstep ; 02237 } 02238 } 02239 } 02240 break ; 02241 02242 case 3:{ 02243 register int out , base , xx,yy , inn ; 02244 02245 base = fixed_index * zstep ; 02246 #ifdef DONT_USE_MEMCPY 02247 out = 0 ; 02248 for( yy=0 ; yy < nyy ; yy++ ){ 02249 inn = base ; 02250 base += ystep ; 02251 for( xx=0 ; xx < nxx ; xx++ ) 02252 bslice[out++] = bold[inn++] ; 02253 } 02254 #else 02255 (void) memcpy( bslice , bold+base , sizeof(DTYPE) * zstep ) ; 02256 #endif 02257 } 02258 break ; 02259 } 02260 02261 EXRETURN ; 02262 } |
|
Definition at line 296 of file afni_slice.c. References THD_linear_mapping::bot, BP_00, BP_P1, CP_00, CP_FACTOR, CP_M1, CP_P1, CP_P2, ENTRY, FINAL, FLOOR, FMAD2, FMAD4, FSCAL, FXYZTMP, FZERO, IBASE, INTYPE, LP_00, LP_P1, m1, MAKE_IBIG, THD_mat33::mat, MAX, THD_linear_mapping::mbac, MIN, NN_ALOOP_GEN, NN_ALOOP_PAR, NN_BLOOP_XY_ZERO, NN_BLOOP_XZ_ZERO, NN_BLOOP_YZ_ZERO, NN_CLOOP_GEN, NN_CLOOP_PAR, NN_ZLOOP, NONE_ZERO, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, OUTADD, RESAM_BLOCK_TYPE, RESAM_CUBIC_TYPE, RESAM_LINEAR_TYPE, RESAM_NN_TYPE, ROUND, STATUS, THD_linear_mapping::svec, TEST_OUT_XXX, TEST_OUT_YYY, TEST_OUT_ZZZ, THREEZ, THD_linear_mapping::top, X_ZERO, XY_ZERO, THD_fvec3::xyz, XYZ_ZERO, XZ_ZERO, Y_ZERO, YZ_ZERO, Z_ZERO, and ZERO.
00299 { 00300 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */ 00301 THD_fvec3 vt = map->svec ; 00302 00303 int xi_new , yj_new , zk_new ; /* voxel indices in new */ 00304 int xi_old , yj_old , zk_old ; /* voxel indices in old */ 00305 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */ 00306 00307 float fxi_top , fyj_top , fzk_top ; /* floating pt. voxel indices */ 00308 float fxi_bot , fyj_bot , fzk_bot ; 00309 float fxi_base , fyj_base , fzk_base ; 00310 float dfxi_outer , dfyj_outer , dfzk_outer ; 00311 float dfxi_inner , dfyj_inner , dfzk_inner ; 00312 00313 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */ 00314 int out_ind , jstep , kstep ; 00315 int nxold,nyold,nzold , nxnew,nynew,nznew ; 00316 00317 ENTRY("AFNI_lmap_to_xslice") ; 00318 00319 /*--- set up ranges ---*/ 00320 00321 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ; 00322 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ; 00323 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ; 00324 00325 if( xi_fix < xi_bot || xi_fix > xi_top ) EXRETURN ; /* map doesn't apply! */ 00326 00327 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ; 00328 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ; 00329 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ; 00330 00331 jstep = nxold ; 00332 kstep = nxold * nyold ; 00333 00334 /* set up base of indices in old */ 00335 00336 xi_new = xi_fix ; 00337 yj_new = yj_bot-1 ; 00338 zk_new = zk_bot-1 ; 00339 00340 fxi_base = mt.mat[0][0] * xi_new 00341 + mt.mat[0][1] * yj_new 00342 + mt.mat[0][2] * zk_new - vt.xyz[0] ; 00343 00344 fyj_base = mt.mat[1][0] * xi_new 00345 + mt.mat[1][1] * yj_new 00346 + mt.mat[1][2] * zk_new - vt.xyz[1] ; 00347 00348 fzk_base = mt.mat[2][0] * xi_new 00349 + mt.mat[2][1] * yj_new 00350 + mt.mat[2][2] * zk_new - vt.xyz[2] ; 00351 00352 dfxi_outer = mt.mat[0][2] ; /* outer loop is in z = 2 */ 00353 dfyj_outer = mt.mat[1][2] ; 00354 dfzk_outer = mt.mat[2][2] ; 00355 00356 dfxi_inner = mt.mat[0][1] ; /* inner loop is in y = 1 */ 00357 dfyj_inner = mt.mat[1][1] ; 00358 dfzk_inner = mt.mat[2][1] ; 00359 00360 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ; 00361 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ; 00362 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ; 00363 00364 switch( resam_mode ){ 00365 00366 default: 00367 case RESAM_NN_TYPE:{ 00368 float fxi_max , fyj_max , fzk_max ; 00369 float fxi_min , fyj_min , fzk_min ; 00370 float fxi_tmp , fyj_tmp , fzk_tmp ; 00371 int any_outside , all_outside ; 00372 00373 #ifdef AFNI_DEBUG 00374 { char str[256] ; 00375 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g", 00376 dfxi_inner,dfyj_inner,dfzk_inner, 00377 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; } 00378 #endif 00379 /** July 15, 1996: 00380 check if all the points are inside the old grid; 00381 if so, can use a version of the resampling loop 00382 that does not need to check each voxel for being 00383 inside -- hopefully, this will execute more quickly **/ 00384 00385 FXYZTMP(xi_new,yj_bot,zk_bot) ; 00386 fxi_max = fxi_min = fxi_tmp ; 00387 fyj_max = fyj_min = fyj_tmp ; 00388 fzk_max = fzk_min = fzk_tmp ; 00389 00390 FXYZTMP(xi_new,yj_top,zk_bot) ; 00391 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 00392 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 00393 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 00394 00395 FXYZTMP(xi_new,yj_bot,zk_top) ; 00396 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 00397 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 00398 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 00399 00400 FXYZTMP(xi_new,yj_top,zk_top) ; 00401 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 00402 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 00403 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 00404 00405 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) || 00406 (fyj_min < fyj_bot) || (fyj_max > fyj_top) || 00407 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ; 00408 00409 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) || 00410 (fyj_max < fyj_bot) || (fyj_min > fyj_top) || 00411 (fzk_max < fzk_bot) || (fzk_min > fzk_top) 00412 : 0 ; 00413 #ifdef AFNI_DEBUG 00414 { char str[256] ; 00415 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max); 00416 STATUS(str) ; 00417 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max); 00418 STATUS(str) ; 00419 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max); 00420 STATUS(str) ; } 00421 #endif 00422 00423 /** redefine the macros specifying loop variables **/ 00424 00425 #undef OUD_NAME 00426 #undef IND_NAME 00427 #undef OUD 00428 #undef OUD_bot 00429 #undef OUD_top 00430 #undef IND 00431 #undef IND_bot 00432 #undef IND_top 00433 #undef IND_nnn 00434 00435 #define OUD_NAME zk /* name of 2nd dimension of output image */ 00436 #define IND_NAME yj /* name of 1st dimension of output image */ 00437 #define IND_nnn nynew /* inner loop image dimension */ 00438 00439 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */ 00440 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */ 00441 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */ 00442 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */ 00443 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */ 00444 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */ 00445 00446 if( all_outside ){ 00447 STATUS("NN resample has all outside") ; 00448 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */ 00449 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */ 00450 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */ 00451 bslice[out_ind++] = ZERO ; 00452 } 00453 } 00454 } else { /* at least some are inside */ 00455 00456 int thz , tho , ob , ub ; 00457 00458 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ; 00459 00460 fxi_top = nxold-0.0001 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */ 00461 fyj_top = nyold-0.0001 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */ 00462 fzk_top = nzold-0.0001 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */ 00463 /* these variables now. */ 00464 00465 /** thz = flag that indicates which of the steps df??_inner are zero. 00466 If two of them are zero, then the inner loop is parallel to 00467 one of the input brick axes, and so data may be pulled 00468 out in a very efficient fashion. In such a case, precompute 00469 the indexes for the inner loop: 00470 00471 the BLOOP macros load array ib, which holds the inner loop 00472 computed indexes for each inner loop position. 00473 ub = upper bound value for ib array value to still be 00474 inside input brick array. 00475 ob = outer loop index into input brick array (computed later) **/ 00476 00477 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */ 00478 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */ 00479 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */ 00480 else /* treatment if outer */ 00481 thz = NONE_ZERO ; /* axes are special */ 00482 00483 #ifdef AFNI_DEBUG 00484 { char str[256] ; 00485 if( any_outside ) sprintf(str,"NN resample has some outside: thz = %d",thz) ; 00486 else sprintf(str,"NN resample has all inside: thz = %d",thz) ; 00487 STATUS(str) ; 00488 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d", 00489 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ; 00490 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d", 00491 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; } 00492 #endif 00493 00494 switch(thz){ 00495 case XY_ZERO: 00496 MAKE_IBIG(IND_top) ; 00497 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ; 00498 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; } 00499 break ; 00500 00501 case XZ_ZERO: 00502 MAKE_IBIG(IND_top) ; 00503 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ; 00504 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; } 00505 break ; 00506 00507 case YZ_ZERO: 00508 MAKE_IBIG(IND_top) ; 00509 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ; 00510 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; } 00511 break ; 00512 } 00513 00514 thz += OUTADD * any_outside ; 00515 00516 STATUS("beginning NN outer loop") ; 00517 00518 /*** outer loop ***/ 00519 00520 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ 00521 00522 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */ 00523 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */ 00524 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */ 00525 00526 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */ 00527 00528 /*** There are 8 cases for the inner loop: 00529 all inside, inner loop not parallel to any input axis 00530 all inside, inner loop parallel to input brick z-axis 00531 all inside, inner loop parallel to input brick y-axis 00532 all inside, inner loop parallel to input brick x-axis 00533 00534 and then the 4 same cases repeated when not all desired 00535 points are inside the input brick. Each of these is 00536 coded separately for efficiency. This is important for 00537 rapid re-display of results during interactive imaging. ***/ 00538 00539 switch(thz){ 00540 case NONE_ZERO: 00541 case X_ZERO: 00542 case Y_ZERO: 00543 case Z_ZERO: 00544 case XYZ_ZERO: 00545 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; } 00546 break ; 00547 00548 case XY_ZERO: 00549 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 00550 ob = IBASE(xi_old,yj_old,0) ; 00551 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 00552 break ; 00553 00554 case XZ_ZERO: 00555 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 00556 ob = IBASE(xi_old,0,zk_old) ; 00557 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 00558 break ; 00559 00560 case YZ_ZERO: 00561 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 00562 ob = IBASE(0,yj_old,zk_old) ; 00563 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 00564 break ; 00565 00566 default: 00567 case NONE_ZERO+OUTADD: 00568 case X_ZERO+OUTADD: 00569 case Y_ZERO+OUTADD: 00570 case Z_ZERO+OUTADD: 00571 case XYZ_ZERO+OUTADD: 00572 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; } 00573 break ; 00574 00575 case XY_ZERO+OUTADD: 00576 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 00577 if( TEST_OUT_XXX || TEST_OUT_YYY ){ 00578 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 00579 } else { 00580 ob = IBASE(xi_old,yj_old,0) ; 00581 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 00582 } 00583 break ; 00584 00585 case XZ_ZERO+OUTADD: 00586 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 00587 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){ 00588 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 00589 } else { 00590 ob = IBASE(xi_old,0,zk_old) ; 00591 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 00592 } 00593 break ; 00594 00595 case YZ_ZERO+OUTADD: 00596 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 00597 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){ 00598 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 00599 } else { 00600 ob = IBASE(0,yj_old,zk_old) ; 00601 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 00602 } 00603 break ; 00604 } 00605 } 00606 } 00607 } 00608 break ; /* end of NN! */ 00609 00610 case RESAM_BLOCK_TYPE: 00611 case RESAM_LINEAR_TYPE:{ 00612 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ; 00613 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 , 00614 f_k00 , f_kp1 , result ; 00615 float frac_xi , frac_yj , frac_zk ; 00616 int ibase , 00617 in_jp1_k00 = jstep , 00618 in_j00_kp1 = kstep , 00619 in_jp1_kp1 = jstep+kstep ; 00620 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 00621 00622 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){ 00623 00624 fxi_old = (fxi_base += dfxi_outer) ; 00625 fyj_old = (fyj_base += dfyj_outer) ; 00626 fzk_old = (fzk_base += dfzk_outer) ; 00627 00628 out_ind = yj_bot + zk_new * nynew ; 00629 00630 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){ 00631 00632 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 00633 00634 if( fxi_old < fxi_bot || fxi_old > fxi_top || 00635 fyj_old < fyj_bot || fyj_old > fyj_top || 00636 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 00637 00638 /*** 00639 bslice[out_ind++] = 0 ; 00640 ***/ 00641 FZERO(bslice[out_ind]) ; out_ind++ ; 00642 continue ; 00643 } 00644 00645 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 00646 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 00647 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 00648 00649 /* use NN if at edges of old brick */ 00650 00651 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 00652 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 00653 00654 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 00655 ROUND(fyj_old) , 00656 ROUND(fzk_old) ) ] ; 00657 continue ; 00658 } 00659 00660 /* compute weights for LINEAR interpolation in each direction */ 00661 00662 if( resam_mode == RESAM_BLOCK_TYPE ){ 00663 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ; 00664 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ; 00665 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ; 00666 } else { 00667 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 00668 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 00669 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 00670 } 00671 00672 /* interpolate in the x direction for each y & z */ 00673 00674 ibase = IBASE(xi_old,yj_old,zk_old) ; 00675 /*** 00676 f_j00_k00 = xwt_00 * bold[ibase] 00677 + xwt_p1 * bold[ibase+1] ; 00678 00679 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 00680 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 00681 00682 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 00683 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 00684 00685 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 00686 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 00687 ***/ 00688 FMAD2( xwt_00 , bold[ibase] , 00689 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 00690 00691 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 00692 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 00693 00694 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 00695 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 00696 00697 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 00698 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 00699 00700 /* interpolate in the y direction for each z */ 00701 00702 /*** 00703 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 00704 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 00705 ***/ 00706 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 00707 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 00708 00709 /* interpolate in the z direction and 00710 put the result into the output array! 00711 (the +0.5 is to force rounding of the result) */ 00712 00713 /*** 00714 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 00715 ***/ 00716 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 00717 FINAL( result , bslice[out_ind] ) ; 00718 out_ind++ ; 00719 00720 } 00721 } 00722 } 00723 break ; 00724 00725 case RESAM_CUBIC_TYPE:{ 00726 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */ 00727 ywt_m1,ywt_00,ywt_p1,ywt_p2 , 00728 zwt_m1,zwt_00,zwt_p1,zwt_p2 ; 00729 00730 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */ 00731 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 , 00732 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 , 00733 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 , 00734 f_km1 , f_k00 , f_kp1 , f_kp2 , result ; 00735 float frac_xi , frac_yj , frac_zk ; 00736 00737 int ibase , /* base index for interpolant */ 00738 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */ 00739 in_jm1_k00 = -jstep , /* 0 (00) */ 00740 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */ 00741 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */ 00742 in_j00_km1 = - kstep , /* steps in j and k indices */ 00743 in_j00_k00 = 0 , 00744 in_j00_kp1 = kstep , 00745 in_j00_kp2 = 2*kstep , 00746 in_jp1_km1 = jstep- kstep , 00747 in_jp1_k00 = jstep , 00748 in_jp1_kp1 = jstep+ kstep , 00749 in_jp1_kp2 =2*jstep+2*kstep , 00750 in_jp2_km1 =2*jstep- kstep , 00751 in_jp2_k00 =2*jstep , 00752 in_jp2_kp1 =2*jstep+ kstep , 00753 in_jp2_kp2 =2*jstep+2*kstep ; 00754 00755 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 00756 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ; 00757 00758 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){ 00759 00760 fxi_old = (fxi_base += dfxi_outer) ; 00761 fyj_old = (fyj_base += dfyj_outer) ; 00762 fzk_old = (fzk_base += dfzk_outer) ; 00763 00764 out_ind = yj_bot + zk_new * nynew ; 00765 00766 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){ 00767 00768 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 00769 00770 /* check if outside old brick */ 00771 00772 if( fxi_old < fxi_bot || fxi_old > fxi_top || 00773 fyj_old < fyj_bot || fyj_old > fyj_top || 00774 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 00775 /*** 00776 bslice[out_ind++] = 0 ; 00777 ***/ 00778 FZERO(bslice[out_ind]) ; out_ind++ ; 00779 continue ; 00780 } 00781 00782 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 00783 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 00784 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 00785 00786 /* use NN if at very edges of old brick */ 00787 00788 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 00789 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 00790 00791 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 00792 ROUND(fyj_old) , 00793 ROUND(fzk_old) ) ] ; 00794 continue ; 00795 } 00796 00797 ibase = IBASE(xi_old,yj_old,zk_old) ; 00798 00799 /* use LINEAR if close to edges of old brick */ 00800 00801 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 || 00802 xi_old == 0 || yj_old == 0 || zk_old == 0 ){ 00803 00804 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 00805 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 00806 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 00807 00808 /*** 00809 f_j00_k00 = xwt_00 * bold[ibase] 00810 + xwt_p1 * bold[ibase+1] ; 00811 00812 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 00813 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 00814 00815 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 00816 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 00817 00818 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 00819 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 00820 ***/ 00821 FMAD2( xwt_00 , bold[ibase] , 00822 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 00823 00824 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 00825 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 00826 00827 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 00828 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 00829 00830 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 00831 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 00832 00833 /*** 00834 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 00835 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 00836 ***/ 00837 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 00838 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 00839 /*** 00840 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 00841 ***/ 00842 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 00843 FINAL( result , bslice[out_ind] ) ; 00844 out_ind++ ; 00845 continue ; 00846 } 00847 00848 /* compute weights for CUBIC interpolation in each direction */ 00849 00850 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ; 00851 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ; 00852 00853 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ; 00854 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ; 00855 00856 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ; 00857 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ; 00858 00859 /* use the ANSI token-merge operator ## to create an interpolating 00860 macro for the x-direction, at each offset in y (j) and z (k) */ 00861 00862 /*** 00863 #define CXINT(j,k) xwt_m1 * bold[ibase + in_j ## j ## _k ## k -1 ] \ 00864 + xwt_00 * bold[ibase + in_j ## j ## _k ## k ] \ 00865 + xwt_p1 * bold[ibase + in_j ## j ## _k ## k +1 ] \ 00866 + xwt_p2 * bold[ibase + in_j ## j ## _k ## k +2 ] 00867 ***/ 00868 00869 #define CXINT(j,k,ff) \ 00870 FMAD4( xwt_m1 , bold[ibase + in_j ## j ## _k ## k -1 ] , \ 00871 xwt_00 , bold[ibase + in_j ## j ## _k ## k ] , \ 00872 xwt_p1 , bold[ibase + in_j ## j ## _k ## k +1 ] , \ 00873 xwt_p2 , bold[ibase + in_j ## j ## _k ## k +2 ] , ff ) 00874 00875 /* interpolate in the x direction for each y & z */ 00876 00877 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ; 00878 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ; 00879 00880 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ; 00881 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ; 00882 00883 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ; 00884 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ; 00885 00886 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ; 00887 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ; 00888 00889 /* interpolate in the y direction for each z */ 00890 00891 /*** 00892 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1 00893 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ; 00894 00895 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00 00896 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ; 00897 00898 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1 00899 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ; 00900 00901 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2 00902 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ; 00903 ***/ 00904 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 , 00905 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ; 00906 00907 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 , 00908 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ; 00909 00910 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 , 00911 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ; 00912 00913 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 , 00914 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ; 00915 00916 /* interpolate in the z direction and 00917 put the result into the output array! 00918 (the +0.5 is to force rounding of the result) */ 00919 /*** 00920 bslice[out_ind++] = 0.5 + CP_FACTOR 00921 * ( zwt_m1 * f_km1 + zwt_00 * f_k00 00922 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ; 00923 ***/ 00924 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 , 00925 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ; 00926 FSCAL(CP_FACTOR,result) ; 00927 FINAL( result , bslice[out_ind] ) ; 00928 out_ind++ ; 00929 00930 } /* end of inner loop */ 00931 } /* end of outer loop */ 00932 } 00933 break ; 00934 00935 } 00936 00937 EXRETURN ; 00938 } |
|
thz = flag that indicates which of the steps df??_inner are zero. If two of them are zero, then the inner loop is parallel to one of the input brick axes, and so data may be pulled out in a very efficient fashion. In such a case, precompute the indexes for the inner loop: the BLOOP macros load array ib, which holds the inner loop computed indexes for each inner loop position. ub = upper bound value for ib array value to still be inside input brick array. ob = outer loop index into input brick array (computed later) * Definition at line 942 of file afni_slice.c. References THD_linear_mapping::bot, BP_00, BP_P1, CP_00, CP_FACTOR, CP_M1, CP_P1, CP_P2, ENTRY, FINAL, FLOOR, FMAD2, FMAD4, FSCAL, FXYZTMP, FZERO, IBASE, INTYPE, LP_00, LP_P1, m1, MAKE_IBIG, THD_mat33::mat, MAX, THD_linear_mapping::mbac, MIN, NN_ALOOP_GEN, NN_ALOOP_PAR, NN_BLOOP_XY_ZERO, NN_BLOOP_XZ_ZERO, NN_BLOOP_YZ_ZERO, NN_CLOOP_GEN, NN_CLOOP_PAR, NN_ZLOOP, NONE_ZERO, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, OUTADD, RESAM_BLOCK_TYPE, RESAM_CUBIC_TYPE, RESAM_LINEAR_TYPE, RESAM_NN_TYPE, ROUND, STATUS, THD_linear_mapping::svec, TEST_OUT_XXX, TEST_OUT_YYY, TEST_OUT_ZZZ, THREEZ, THD_linear_mapping::top, X_ZERO, XY_ZERO, THD_fvec3::xyz, XYZ_ZERO, XZ_ZERO, Y_ZERO, YZ_ZERO, Z_ZERO, and ZERO.
00945 { 00946 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */ 00947 THD_fvec3 vt = map->svec ; 00948 00949 int xi_new , yj_new , zk_new ; /* voxel indices in new */ 00950 int xi_old , yj_old , zk_old ; /* voxel indices in old */ 00951 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */ 00952 00953 float fxi_top , fyj_top , fzk_top ; 00954 float fxi_bot , fyj_bot , fzk_bot ; 00955 float fxi_base , fyj_base , fzk_base ; 00956 float dfxi_outer , dfyj_outer , dfzk_outer ; 00957 float dfxi_inner , dfyj_inner , dfzk_inner ; 00958 00959 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */ 00960 int out_ind , jstep , kstep ; 00961 int nxold,nyold,nzold , nxnew,nynew,nznew ; 00962 00963 ENTRY("AFNI_lmap_to_yslice") ; 00964 00965 /*--- set up ranges ---*/ 00966 00967 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ; 00968 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ; 00969 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ; 00970 00971 if( yj_fix < yj_bot || yj_fix > yj_top ) EXRETURN ; /* map doesn't apply! */ 00972 00973 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ; 00974 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ; 00975 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ; 00976 00977 jstep = nxold ; 00978 kstep = nxold * nyold ; 00979 00980 /* set up base of indices in old */ 00981 00982 xi_new = xi_bot-1 ; 00983 yj_new = yj_fix ; 00984 zk_new = zk_bot-1 ; 00985 00986 fxi_base = mt.mat[0][0] * xi_new 00987 + mt.mat[0][1] * yj_new 00988 + mt.mat[0][2] * zk_new - vt.xyz[0] ; 00989 00990 fyj_base = mt.mat[1][0] * xi_new 00991 + mt.mat[1][1] * yj_new 00992 + mt.mat[1][2] * zk_new - vt.xyz[1] ; 00993 00994 fzk_base = mt.mat[2][0] * xi_new 00995 + mt.mat[2][1] * yj_new 00996 + mt.mat[2][2] * zk_new - vt.xyz[2] ; 00997 00998 dfxi_outer = mt.mat[0][0] ; /* outer loop is in x = 0 */ 00999 dfyj_outer = mt.mat[1][0] ; 01000 dfzk_outer = mt.mat[2][0] ; 01001 01002 dfxi_inner = mt.mat[0][2] ; /* inner loop is in z = 2 */ 01003 dfyj_inner = mt.mat[1][2] ; 01004 dfzk_inner = mt.mat[2][2] ; 01005 01006 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ; 01007 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ; 01008 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ; 01009 01010 switch( resam_mode ){ 01011 01012 default: 01013 case RESAM_NN_TYPE:{ 01014 float fxi_max , fyj_max , fzk_max ; 01015 float fxi_min , fyj_min , fzk_min ; 01016 float fxi_tmp , fyj_tmp , fzk_tmp ; 01017 int any_outside , all_outside ; 01018 01019 #ifdef AFNI_DEBUG 01020 { char str[256] ; 01021 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g", 01022 dfxi_inner,dfyj_inner,dfzk_inner, 01023 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; } 01024 #endif 01025 /** July 15, 1996: 01026 check if all the points are inside the old grid; 01027 if so, can use a version of the resampling loop 01028 that does not need to check each voxel for being 01029 inside -- hopefully, this will execute more quickly **/ 01030 01031 FXYZTMP(xi_bot,yj_new,zk_bot) ; 01032 fxi_max = fxi_min = fxi_tmp ; 01033 fyj_max = fyj_min = fyj_tmp ; 01034 fzk_max = fzk_min = fzk_tmp ; 01035 01036 FXYZTMP(xi_top,yj_new,zk_bot) ; 01037 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01038 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01039 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01040 01041 FXYZTMP(xi_bot,yj_new,zk_top) ; 01042 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01043 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01044 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01045 01046 FXYZTMP(xi_top,yj_new,zk_top) ; 01047 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01048 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01049 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01050 01051 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) || 01052 (fyj_min < fyj_bot) || (fyj_max > fyj_top) || 01053 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ; 01054 01055 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) || 01056 (fyj_max < fyj_bot) || (fyj_min > fyj_top) || 01057 (fzk_max < fzk_bot) || (fzk_min > fzk_top) 01058 : 0 ; 01059 01060 #ifdef AFNI_DEBUG 01061 { char str[256] ; 01062 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max); 01063 STATUS(str) ; 01064 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max); 01065 STATUS(str) ; 01066 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max); 01067 STATUS(str) ; } 01068 #endif 01069 01070 /** redefine the macros specifying loop variables **/ 01071 01072 #undef OUD_NAME 01073 #undef IND_NAME 01074 #undef OUD 01075 #undef OUD_bot 01076 #undef OUD_top 01077 #undef IND 01078 #undef IND_bot 01079 #undef IND_top 01080 #undef IND_nnn 01081 01082 #define OUD_NAME xi /* name of 2nd dimension of output image */ 01083 #define IND_NAME zk /* name of 1st dimension of output image */ 01084 #define IND_nnn nznew /* inner loop image dimension */ 01085 01086 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */ 01087 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */ 01088 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */ 01089 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */ 01090 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */ 01091 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */ 01092 01093 if( all_outside ){ 01094 STATUS("NN resample has all outside") ; 01095 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */ 01096 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */ 01097 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */ 01098 bslice[out_ind++] = ZERO ; 01099 } 01100 } 01101 } else { /* at least some are inside */ 01102 01103 int thz, tho , ob , ub ; 01104 01105 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ; 01106 01107 fxi_top = nxold - 0.01 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */ 01108 fyj_top = nyold - 0.01 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */ 01109 fzk_top = nzold - 0.01 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */ 01110 /* these variables now. */ 01111 01112 /** thz = flag that indicates which of the steps df??_inner are zero. 01113 If two of them are zero, then the inner loop is parallel to 01114 one of the input brick axes, and so data may be pulled 01115 out in a very efficient fashion. In such a case, precompute 01116 the indexes for the inner loop: 01117 01118 the BLOOP macros load array ib, which holds the inner loop 01119 computed indexes for each inner loop position. 01120 ub = upper bound value for ib array value to still be 01121 inside input brick array. 01122 ob = outer loop index into input brick array (computed later) **/ 01123 01124 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */ 01125 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */ 01126 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */ 01127 else /* treatment if outer */ 01128 thz = NONE_ZERO ; /* axes are special */ 01129 01130 #ifdef AFNI_DEBUG 01131 { char str[256] ; 01132 if( any_outside ) sprintf(str,"NN resample has some outside: thz = %d",thz) ; 01133 else sprintf(str,"NN resample has all inside: thz = %d",thz) ; 01134 STATUS(str) ; 01135 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d", 01136 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ; 01137 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d", 01138 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; } 01139 #endif 01140 01141 switch(thz){ 01142 case XY_ZERO: 01143 MAKE_IBIG(IND_top) ; 01144 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ; 01145 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; } 01146 break ; 01147 01148 case XZ_ZERO: 01149 MAKE_IBIG(IND_top) ; 01150 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ; 01151 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; } 01152 break ; 01153 01154 case YZ_ZERO: 01155 MAKE_IBIG(IND_top) ; 01156 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ; 01157 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; } 01158 break ; 01159 } 01160 01161 thz += OUTADD * any_outside ; 01162 01163 STATUS("beginning NN outer loop") ; 01164 01165 /*** outer loop ***/ 01166 01167 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ 01168 01169 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */ 01170 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */ 01171 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */ 01172 01173 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */ 01174 01175 /*** There are 8 cases for the inner loop: 01176 all inside, inner loop not parallel to any input axis 01177 all inside, inner loop parallel to input brick z-axis 01178 all inside, inner loop parallel to input brick y-axis 01179 all inside, inner loop parallel to input brick x-axis 01180 01181 and then the 4 same cases repeated when not all desired 01182 points are inside the input brick. Each of these is 01183 coded separately for efficiency. This is important for 01184 rapid re-display of results during interactive imaging. ***/ 01185 01186 switch(thz){ 01187 case NONE_ZERO: 01188 case X_ZERO: 01189 case Y_ZERO: 01190 case Z_ZERO: 01191 case XYZ_ZERO: 01192 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; } 01193 break ; 01194 01195 case XY_ZERO: 01196 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 01197 ob = IBASE(xi_old,yj_old,0) ; 01198 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01199 break ; 01200 01201 case XZ_ZERO: 01202 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 01203 ob = IBASE(xi_old,0,zk_old) ; 01204 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01205 break ; 01206 01207 case YZ_ZERO: 01208 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 01209 ob = IBASE(0,yj_old,zk_old) ; 01210 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01211 break ; 01212 01213 default: 01214 case NONE_ZERO+OUTADD: 01215 case X_ZERO+OUTADD: 01216 case Y_ZERO+OUTADD: 01217 case Z_ZERO+OUTADD: 01218 case XYZ_ZERO+OUTADD: 01219 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; } 01220 break ; 01221 01222 case XY_ZERO+OUTADD: 01223 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 01224 if( TEST_OUT_XXX || TEST_OUT_YYY ){ 01225 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01226 } else { 01227 ob = IBASE(xi_old,yj_old,0) ; 01228 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01229 } 01230 break ; 01231 01232 case XZ_ZERO+OUTADD: 01233 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 01234 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){ 01235 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01236 } else { 01237 ob = IBASE(xi_old,0,zk_old) ; 01238 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01239 } 01240 break ; 01241 01242 case YZ_ZERO+OUTADD: 01243 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 01244 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){ 01245 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01246 } else { 01247 ob = IBASE(0,yj_old,zk_old) ; 01248 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01249 } 01250 break ; 01251 } 01252 } 01253 } 01254 } 01255 break ; /* end of NN! */ 01256 01257 case RESAM_BLOCK_TYPE: 01258 case RESAM_LINEAR_TYPE:{ 01259 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ; 01260 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 , 01261 f_k00 , f_kp1 , result ; 01262 float frac_xi , frac_yj , frac_zk ; 01263 int ibase , 01264 in_jp1_k00 = jstep , 01265 in_j00_kp1 = kstep , 01266 in_jp1_kp1 = jstep+kstep ; 01267 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 01268 01269 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){ 01270 01271 fxi_old = (fxi_base += dfxi_outer) ; 01272 fyj_old = (fyj_base += dfyj_outer) ; 01273 fzk_old = (fzk_base += dfzk_outer) ; 01274 01275 out_ind = zk_bot + xi_new * nznew ; 01276 01277 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){ 01278 01279 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 01280 01281 if( fxi_old < fxi_bot || fxi_old > fxi_top || 01282 fyj_old < fyj_bot || fyj_old > fyj_top || 01283 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 01284 01285 /*** 01286 bslice[out_ind++] = 0 ; 01287 ***/ 01288 FZERO(bslice[out_ind]) ; out_ind++ ; 01289 continue ; 01290 } 01291 01292 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 01293 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 01294 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 01295 01296 /* use NN if at edges of old brick */ 01297 01298 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 01299 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 01300 01301 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 01302 ROUND(fyj_old) , 01303 ROUND(fzk_old) ) ] ; 01304 continue ; 01305 } 01306 01307 /* compute weights for LINEAR interpolation in each direction */ 01308 01309 if( resam_mode == RESAM_BLOCK_TYPE ){ 01310 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ; 01311 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ; 01312 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ; 01313 } else { 01314 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 01315 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 01316 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 01317 } 01318 01319 /* interpolate in the x direction for each y & z */ 01320 01321 ibase = IBASE(xi_old,yj_old,zk_old) ; 01322 /*** 01323 f_j00_k00 = xwt_00 * bold[ibase] 01324 + xwt_p1 * bold[ibase+1] ; 01325 01326 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 01327 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 01328 01329 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 01330 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 01331 01332 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 01333 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 01334 ***/ 01335 FMAD2( xwt_00 , bold[ibase] , 01336 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 01337 01338 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 01339 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 01340 01341 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 01342 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 01343 01344 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 01345 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 01346 01347 /* interpolate in the y direction for each z */ 01348 /*** 01349 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 01350 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 01351 ***/ 01352 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 01353 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 01354 01355 /* interpolate in the z direction and 01356 put the result into the output array! 01357 (the +0.5 is to force rounding of the result) */ 01358 /*** 01359 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 01360 ***/ 01361 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 01362 FINAL( result , bslice[out_ind] ) ; 01363 out_ind++ ; 01364 } 01365 } 01366 } 01367 break ; 01368 01369 case RESAM_CUBIC_TYPE:{ 01370 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */ 01371 ywt_m1,ywt_00,ywt_p1,ywt_p2 , 01372 zwt_m1,zwt_00,zwt_p1,zwt_p2 ; 01373 01374 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */ 01375 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 , 01376 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 , 01377 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 , 01378 f_km1 , f_k00 , f_kp1 , f_kp2 , result ; 01379 float frac_xi , frac_yj , frac_zk ; 01380 01381 int ibase , /* base index for interpolant */ 01382 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */ 01383 in_jm1_k00 = -jstep , /* 0 (00) */ 01384 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */ 01385 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */ 01386 in_j00_km1 = - kstep , /* steps in j and k indices */ 01387 in_j00_k00 = 0 , 01388 in_j00_kp1 = kstep , 01389 in_j00_kp2 = 2*kstep , 01390 in_jp1_km1 = jstep- kstep , 01391 in_jp1_k00 = jstep , 01392 in_jp1_kp1 = jstep+ kstep , 01393 in_jp1_kp2 =2*jstep+2*kstep , 01394 in_jp2_km1 =2*jstep- kstep , 01395 in_jp2_k00 =2*jstep , 01396 in_jp2_kp1 =2*jstep+ kstep , 01397 in_jp2_kp2 =2*jstep+2*kstep ; 01398 01399 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 01400 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ; 01401 01402 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){ 01403 01404 fxi_old = (fxi_base += dfxi_outer) ; 01405 fyj_old = (fyj_base += dfyj_outer) ; 01406 fzk_old = (fzk_base += dfzk_outer) ; 01407 01408 out_ind = zk_bot + xi_new * nznew ; 01409 01410 for( zk_new=zk_bot ; zk_new <= zk_top ; zk_new++ ){ 01411 01412 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 01413 01414 /* check if outside old brick */ 01415 01416 if( fxi_old < fxi_bot || fxi_old > fxi_top || 01417 fyj_old < fyj_bot || fyj_old > fyj_top || 01418 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 01419 /*** 01420 bslice[out_ind++] = 0 ; 01421 ***/ 01422 FZERO(bslice[out_ind]) ; out_ind++ ; 01423 continue ; 01424 } 01425 01426 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 01427 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 01428 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 01429 01430 /* use NN if at very edges of old brick */ 01431 01432 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 01433 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 01434 01435 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 01436 ROUND(fyj_old) , 01437 ROUND(fzk_old) ) ] ; 01438 continue ; 01439 } 01440 01441 ibase = IBASE(xi_old,yj_old,zk_old) ; 01442 01443 /* use LINEAR if close to edges of old brick */ 01444 01445 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 || 01446 xi_old == 0 || yj_old == 0 || zk_old == 0 ){ 01447 01448 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 01449 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 01450 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 01451 /*** 01452 f_j00_k00 = xwt_00 * bold[ibase] 01453 + xwt_p1 * bold[ibase+1] ; 01454 01455 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 01456 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 01457 01458 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 01459 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 01460 01461 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 01462 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 01463 ***/ 01464 FMAD2( xwt_00 , bold[ibase] , 01465 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 01466 01467 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 01468 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 01469 01470 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 01471 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 01472 01473 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 01474 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 01475 /*** 01476 01477 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 01478 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 01479 ***/ 01480 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 01481 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 01482 /*** 01483 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 01484 ***/ 01485 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 01486 FINAL( result , bslice[out_ind] ) ; 01487 out_ind++ ; 01488 continue ; 01489 } 01490 01491 /* compute weights for CUBIC interpolation in each direction */ 01492 01493 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ; 01494 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ; 01495 01496 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ; 01497 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ; 01498 01499 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ; 01500 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ; 01501 01502 /* interpolate in the x direction for each y & z */ 01503 01504 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ; 01505 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ; 01506 01507 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ; 01508 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ; 01509 01510 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ; 01511 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ; 01512 01513 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ; 01514 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ; 01515 01516 /* interpolate in the y direction for each z */ 01517 /*** 01518 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1 01519 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ; 01520 01521 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00 01522 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ; 01523 01524 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1 01525 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ; 01526 01527 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2 01528 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ; 01529 ***/ 01530 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 , 01531 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ; 01532 01533 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 , 01534 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ; 01535 01536 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 , 01537 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ; 01538 01539 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 , 01540 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ; 01541 01542 /* interpolate in the z direction and 01543 put the result into the output array! 01544 (the +0.5 is to force rounding of the result) */ 01545 /*** 01546 bslice[out_ind++] = 0.5 + CP_FACTOR 01547 * ( zwt_m1 * f_km1 + zwt_00 * f_k00 01548 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ; 01549 ***/ 01550 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 , 01551 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ; 01552 FSCAL(CP_FACTOR,result) ; 01553 FINAL( result , bslice[out_ind] ) ; 01554 out_ind++ ; 01555 01556 } /* end of inner loop */ 01557 } /* end of outer loop */ 01558 } 01559 break ; 01560 } 01561 01562 EXRETURN ; 01563 } |
|
thz = flag that indicates which of the steps df??_inner are zero. If two of them are zero, then the inner loop is parallel to one of the input brick axes, and so data may be pulled out in a very efficient fashion. In such a case, precompute the indexes for the inner loop: the BLOOP macros load array ib, which holds the inner loop computed indexes for each inner loop position. ub = upper bound value for ib array value to still be inside input brick array. ob = outer loop index into input brick array (computed later) * Definition at line 1567 of file afni_slice.c. References THD_linear_mapping::bot, BP_00, BP_P1, CP_00, CP_FACTOR, CP_M1, CP_P1, CP_P2, ENTRY, FINAL, FLOOR, FMAD2, FMAD4, FSCAL, FXYZTMP, FZERO, IBASE, INTYPE, LP_00, LP_P1, m1, MAKE_IBIG, THD_mat33::mat, MAX, THD_linear_mapping::mbac, MIN, NN_ALOOP_GEN, NN_ALOOP_PAR, NN_BLOOP_XY_ZERO, NN_BLOOP_XZ_ZERO, NN_BLOOP_YZ_ZERO, NN_CLOOP_GEN, NN_CLOOP_PAR, NN_ZLOOP, NONE_ZERO, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, OUTADD, RESAM_BLOCK_TYPE, RESAM_CUBIC_TYPE, RESAM_LINEAR_TYPE, RESAM_NN_TYPE, ROUND, STATUS, THD_linear_mapping::svec, TEST_OUT_XXX, TEST_OUT_YYY, TEST_OUT_ZZZ, THREEZ, THD_linear_mapping::top, X_ZERO, XY_ZERO, THD_fvec3::xyz, XYZ_ZERO, XZ_ZERO, Y_ZERO, YZ_ZERO, Z_ZERO, and ZERO.
01570 { 01571 THD_mat33 mt = map->mbac ; /* map from bslice indices to bold */ 01572 THD_fvec3 vt = map->svec ; 01573 01574 int xi_new , yj_new , zk_new ; /* voxel indices in new */ 01575 int xi_old , yj_old , zk_old ; /* voxel indices in old */ 01576 float fxi_old , fyj_old , fzk_old ; /* voxel indices in old */ 01577 01578 float fxi_top , fyj_top , fzk_top ; 01579 float fxi_bot , fyj_bot , fzk_bot ; 01580 float fxi_base , fyj_base , fzk_base ; 01581 float dfxi_outer , dfyj_outer , dfzk_outer ; 01582 float dfxi_inner , dfyj_inner , dfzk_inner ; 01583 01584 int xi_bot,xi_top , yj_bot,yj_top , zk_bot,zk_top ; /* ranges in new */ 01585 int out_ind , jstep , kstep ; 01586 int nxold,nyold,nzold , nxnew,nynew,nznew ; 01587 01588 ENTRY("AFNI_lmap_to_zslice") ; 01589 01590 /*--- set up ranges ---*/ 01591 01592 xi_bot = map->bot.xyz[0] ; xi_top = map->top.xyz[0] ; 01593 yj_bot = map->bot.xyz[1] ; yj_top = map->top.xyz[1] ; 01594 zk_bot = map->bot.xyz[2] ; zk_top = map->top.xyz[2] ; 01595 01596 if( zk_fix < zk_bot || zk_fix > zk_top ) EXRETURN ; /* map doesn't apply! */ 01597 01598 nxold = old_daxes->nxx ; nxnew = new_daxes->nxx ; 01599 nyold = old_daxes->nyy ; nynew = new_daxes->nyy ; 01600 nzold = old_daxes->nzz ; nznew = new_daxes->nzz ; 01601 01602 jstep = nxold ; 01603 kstep = nxold * nyold ; 01604 01605 /* set up base of indices in old */ 01606 01607 xi_new = xi_bot-1 ; 01608 yj_new = yj_bot-1 ; 01609 zk_new = zk_fix ; 01610 01611 fxi_base = mt.mat[0][0] * xi_new 01612 + mt.mat[0][1] * yj_new 01613 + mt.mat[0][2] * zk_new - vt.xyz[0] ; 01614 01615 fyj_base = mt.mat[1][0] * xi_new 01616 + mt.mat[1][1] * yj_new 01617 + mt.mat[1][2] * zk_new - vt.xyz[1] ; 01618 01619 fzk_base = mt.mat[2][0] * xi_new 01620 + mt.mat[2][1] * yj_new 01621 + mt.mat[2][2] * zk_new - vt.xyz[2] ; 01622 01623 dfxi_outer = mt.mat[0][1] ; /* outer loop is in y = 1 */ 01624 dfyj_outer = mt.mat[1][1] ; 01625 dfzk_outer = mt.mat[2][1] ; 01626 01627 dfxi_inner = mt.mat[0][0] ; /* inner loop is in x = 0 */ 01628 dfyj_inner = mt.mat[1][0] ; 01629 dfzk_inner = mt.mat[2][0] ; 01630 01631 fxi_top = nxold - 0.51 ; fxi_bot = -0.49 ; 01632 fyj_top = nyold - 0.51 ; fyj_bot = -0.49 ; 01633 fzk_top = nzold - 0.51 ; fzk_bot = -0.49 ; 01634 01635 switch( resam_mode ){ 01636 01637 default: 01638 case RESAM_NN_TYPE:{ 01639 float fxi_max , fyj_max , fzk_max ; 01640 float fxi_min , fyj_min , fzk_min ; 01641 float fxi_tmp , fyj_tmp , fzk_tmp ; 01642 int any_outside , all_outside ; 01643 01644 #ifdef AFNI_DEBUG 01645 { char str[256] ; 01646 sprintf(str,"NN inner dxyz=%g %g %g outer dxyz=%g %g %g", 01647 dfxi_inner,dfyj_inner,dfzk_inner, 01648 dfxi_outer,dfyj_outer,dfzk_outer ) ; STATUS(str) ; } 01649 #endif 01650 /** July 15, 1996: 01651 check if all the points are inside the old grid; 01652 if so, can use a version of the resampling loop 01653 that does not need to check each voxel for being 01654 inside -- hopefully, this will execute more quickly **/ 01655 01656 FXYZTMP(xi_bot,yj_bot,zk_new) ; 01657 fxi_max = fxi_min = fxi_tmp ; 01658 fyj_max = fyj_min = fyj_tmp ; 01659 fzk_max = fzk_min = fzk_tmp ; 01660 01661 FXYZTMP(xi_top,yj_bot,zk_new) ; 01662 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01663 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01664 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01665 01666 FXYZTMP(xi_bot,yj_top,zk_new) ; 01667 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01668 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01669 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01670 01671 FXYZTMP(xi_top,yj_top,zk_new) ; 01672 fxi_max = MAX(fxi_max,fxi_tmp) ; fxi_min = MIN(fxi_min,fxi_tmp) ; 01673 fyj_max = MAX(fyj_max,fyj_tmp) ; fyj_min = MIN(fyj_min,fyj_tmp) ; 01674 fzk_max = MAX(fzk_max,fzk_tmp) ; fzk_min = MIN(fzk_min,fzk_tmp) ; 01675 01676 any_outside = (fxi_min < fxi_bot) || (fxi_max > fxi_top) || 01677 (fyj_min < fyj_bot) || (fyj_max > fyj_top) || 01678 (fzk_min < fzk_bot) || (fzk_max > fzk_top) ; 01679 01680 all_outside = (any_outside) ? (fxi_max < fxi_bot) || (fxi_min > fxi_top) || 01681 (fyj_max < fyj_bot) || (fyj_min > fyj_top) || 01682 (fzk_max < fzk_bot) || (fzk_min > fzk_top) 01683 : 0 ; 01684 01685 #ifdef AFNI_DEBUG 01686 { char str[256] ; 01687 sprintf(str,"fxi_bot=%g fxi_top=%g fxi_min=%g fxi_max=%g",fxi_bot,fxi_top,fxi_min,fxi_max); 01688 STATUS(str) ; 01689 sprintf(str,"fyj_bot=%g fyj_top=%g fyj_min=%g fyj_max=%g",fyj_bot,fyj_top,fyj_min,fyj_max); 01690 STATUS(str) ; 01691 sprintf(str,"fzk_bot=%g fzk_top=%g fzk_min=%g fzk_max=%g",fzk_bot,fzk_top,fzk_min,fzk_max); 01692 STATUS(str) ; } 01693 #endif 01694 01695 /** redefine the macros specifying loop variables **/ 01696 01697 #undef OUD_NAME 01698 #undef IND_NAME 01699 #undef OUD 01700 #undef OUD_bot 01701 #undef OUD_top 01702 #undef IND 01703 #undef IND_bot 01704 #undef IND_top 01705 #undef IND_nnn 01706 01707 #define OUD_NAME yj /* name of 2nd dimension of output image */ 01708 #define IND_NAME xi /* name of 1st dimension of output image */ 01709 #define IND_nnn nxnew /* inner loop image dimension */ 01710 01711 #define OUD TWO_TWO(OUD_NAME , _new) /* outer loop index name */ 01712 #define OUD_bot TWO_TWO(OUD_NAME , _bot) /* outer loop index min */ 01713 #define OUD_top TWO_TWO(OUD_NAME , _top) /* outer loop index max */ 01714 #define IND TWO_TWO(IND_NAME , _new) /* inner loop index name */ 01715 #define IND_bot TWO_TWO(IND_NAME , _bot) /* inner loop index min */ 01716 #define IND_top TWO_TWO(IND_NAME , _top) /* inner loop index max */ 01717 01718 if( all_outside ){ 01719 STATUS("NN resample has all outside") ; 01720 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ /* all points are */ 01721 out_ind = IND_bot + OUD * IND_nnn ; /* outside input brick */ 01722 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ /* so just load zeros */ 01723 bslice[out_ind++] = ZERO ; 01724 } 01725 } 01726 } else { /* at least some are inside */ 01727 01728 int thz , tho , ob , ub ; 01729 01730 fxi_base += 0.5 ; fyj_base += 0.5 ; fzk_base += 0.5 ; 01731 01732 fxi_top = nxold - 0.01 ; fxi_bot = 0.0 ; /* we can use FLOOR instead */ 01733 fyj_top = nyold - 0.01 ; fyj_bot = 0.0 ; /* of ROUND to find the NN */ 01734 fzk_top = nzold - 0.01 ; fzk_bot = 0.0 ; /* by adding 0.5 to all */ 01735 /* these variables now. */ 01736 01737 /** thz = flag that indicates which of the steps df??_inner are zero. 01738 If two of them are zero, then the inner loop is parallel to 01739 one of the input brick axes, and so data may be pulled 01740 out in a very efficient fashion. In such a case, precompute 01741 the indexes for the inner loop: 01742 01743 the BLOOP macros load array ib, which holds the inner loop 01744 computed indexes for each inner loop position. 01745 ub = upper bound value for ib array value to still be 01746 inside input brick array. 01747 ob = outer loop index into input brick array (computed later) **/ 01748 01749 tho = THREEZ(dfxi_outer,dfyj_outer,dfzk_outer) ; /* 06 Aug 1996: */ 01750 if( tho == XY_ZERO || tho == XZ_ZERO || tho == YZ_ZERO ) /* only allow thz to */ 01751 thz = THREEZ(dfxi_inner,dfyj_inner,dfzk_inner) ; /* indicate special */ 01752 else /* treatment if outer */ 01753 thz = NONE_ZERO ; /* axes are special */ 01754 01755 #ifdef AFNI_DEBUG 01756 { char str[256] ; 01757 if( any_outside ) sprintf(str,"NN resample has some outside: thz = %d",thz) ; 01758 else sprintf(str,"NN resample has all inside: thz = %d",thz) ; 01759 STATUS(str) ; 01760 sprintf(str,"OUD_bot=%d OUD_top=%d nxold=%d nyold=%d nzold=%d", 01761 OUD_bot,OUD_top,nxold,nyold,nzold ) ; STATUS(str) ; 01762 sprintf(str,"IND_bot=%d IND_top=%d nxold=%d nyold=%d nzold=%d", 01763 IND_bot,IND_top,nxold,nyold,nzold ) ; STATUS(str) ; } 01764 #endif 01765 01766 switch(thz){ 01767 case XY_ZERO: 01768 MAKE_IBIG(IND_top) ; 01769 fzk_old = fzk_base + dfzk_outer ; ub = IBASE(0,0,nzold) ; 01770 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XY_ZERO(IND) ; } 01771 break ; 01772 01773 case XZ_ZERO: 01774 MAKE_IBIG(IND_top) ; 01775 fyj_old = fyj_base + dfyj_outer ; ub = IBASE(0,nyold,0) ; 01776 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_XZ_ZERO(IND) ; } 01777 break ; 01778 01779 case YZ_ZERO: 01780 MAKE_IBIG(IND_top) ; 01781 fxi_old = fxi_base + dfxi_outer ; ub = IBASE(nxold,0,0) ; 01782 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_BLOOP_YZ_ZERO(IND) ; } 01783 break ; 01784 } 01785 01786 thz += OUTADD * any_outside ; 01787 01788 STATUS("beginning NN outer loop") ; 01789 01790 /*** outer loop ***/ 01791 01792 for( OUD=OUD_bot ; OUD <= OUD_top ; OUD++ ){ 01793 01794 fxi_old = (fxi_base += dfxi_outer) ; /* floating indexes in */ 01795 fyj_old = (fyj_base += dfyj_outer) ; /* input brick at start */ 01796 fzk_old = (fzk_base += dfzk_outer) ; /* of next inner loop */ 01797 01798 out_ind = IND_bot + OUD * IND_nnn ; /* index into output brick */ 01799 01800 /*** There are 8 cases for the inner loop: 01801 all inside, inner loop not parallel to any input axis 01802 all inside, inner loop parallel to input brick z-axis 01803 all inside, inner loop parallel to input brick y-axis 01804 all inside, inner loop parallel to input brick x-axis 01805 01806 and then the 4 same cases repeated when not all desired 01807 points are inside the input brick. Each of these is 01808 coded separately for efficiency. This is important for 01809 rapid re-display of results during interactive imaging. ***/ 01810 01811 switch(thz){ 01812 case NONE_ZERO: 01813 case X_ZERO: 01814 case Y_ZERO: 01815 case Z_ZERO: 01816 case XYZ_ZERO: 01817 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_GEN ; } 01818 break ; 01819 01820 case XY_ZERO: 01821 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 01822 ob = IBASE(xi_old,yj_old,0) ; 01823 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01824 break ; 01825 01826 case XZ_ZERO: 01827 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 01828 ob = IBASE(xi_old,0,zk_old) ; 01829 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01830 break ; 01831 01832 case YZ_ZERO: 01833 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 01834 ob = IBASE(0,yj_old,zk_old) ; 01835 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ALOOP_PAR(IND) ; } 01836 break ; 01837 01838 default: 01839 case NONE_ZERO+OUTADD: 01840 case X_ZERO+OUTADD: 01841 case Y_ZERO+OUTADD: 01842 case Z_ZERO+OUTADD: 01843 case XYZ_ZERO+OUTADD: 01844 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_GEN ; } 01845 break ; 01846 01847 case XY_ZERO+OUTADD: 01848 xi_old = FLOOR( fxi_old ) ; yj_old = FLOOR( fyj_old ) ; 01849 if( TEST_OUT_XXX || TEST_OUT_YYY ){ 01850 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01851 } else { 01852 ob = IBASE(xi_old,yj_old,0) ; 01853 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01854 } 01855 break ; 01856 01857 case XZ_ZERO+OUTADD: 01858 xi_old = FLOOR( fxi_old ) ; zk_old = FLOOR( fzk_old ) ; 01859 if( TEST_OUT_XXX || TEST_OUT_ZZZ ){ 01860 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01861 } else { 01862 ob = IBASE(xi_old,0,zk_old) ; 01863 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01864 } 01865 break ; 01866 01867 case YZ_ZERO+OUTADD: 01868 yj_old = FLOOR( fyj_old ) ; zk_old = FLOOR( fzk_old ) ; 01869 if( TEST_OUT_YYY || TEST_OUT_ZZZ ){ 01870 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_ZLOOP ; } 01871 } else { 01872 ob = IBASE(0,yj_old,zk_old) ; 01873 for( IND=IND_bot ; IND <= IND_top ; IND++ ){ NN_CLOOP_PAR(IND) ; } 01874 } 01875 break ; 01876 } 01877 } 01878 } 01879 } 01880 break ; /* end of NN! */ 01881 01882 case RESAM_BLOCK_TYPE: 01883 case RESAM_LINEAR_TYPE:{ 01884 float xwt_00,xwt_p1 , ywt_00,ywt_p1 , zwt_00,zwt_p1 ; 01885 INTYPE f_j00_k00 , f_jp1_k00 , f_j00_kp1 , f_jp1_kp1 , 01886 f_k00 , f_kp1 , result ; 01887 float frac_xi , frac_yj , frac_zk ; 01888 int ibase , 01889 in_jp1_k00 = jstep , 01890 in_j00_kp1 = kstep , 01891 in_jp1_kp1 = jstep+kstep ; 01892 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 01893 01894 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){ 01895 01896 fxi_old = (fxi_base += dfxi_outer) ; 01897 fyj_old = (fyj_base += dfyj_outer) ; 01898 fzk_old = (fzk_base += dfzk_outer) ; 01899 01900 out_ind = xi_bot + yj_new * nxnew ; 01901 01902 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){ 01903 01904 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 01905 01906 if( fxi_old < fxi_bot || fxi_old > fxi_top || 01907 fyj_old < fyj_bot || fyj_old > fyj_top || 01908 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 01909 01910 /*** 01911 bslice[out_ind++] = 0 ; 01912 ***/ 01913 FZERO(bslice[out_ind]) ; out_ind++ ; 01914 continue ; 01915 } 01916 01917 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 01918 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 01919 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 01920 01921 /* use NN if at edges of old brick */ 01922 01923 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 01924 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 01925 01926 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 01927 ROUND(fyj_old) , 01928 ROUND(fzk_old) ) ] ; 01929 continue ; 01930 } 01931 01932 /* compute weights for LINEAR interpolation in each direction */ 01933 01934 if( resam_mode == RESAM_BLOCK_TYPE ){ 01935 xwt_00 = BP_00(frac_xi) ; xwt_p1 = BP_P1(frac_xi) ; 01936 ywt_00 = BP_00(frac_yj) ; ywt_p1 = BP_P1(frac_yj) ; 01937 zwt_00 = BP_00(frac_zk) ; zwt_p1 = BP_P1(frac_zk) ; 01938 } else { 01939 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 01940 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 01941 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 01942 } 01943 01944 /* interpolate in the x direction for each y & z */ 01945 01946 ibase = IBASE(xi_old,yj_old,zk_old) ; 01947 /*** 01948 f_j00_k00 = xwt_00 * bold[ibase] 01949 + xwt_p1 * bold[ibase+1] ; 01950 01951 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 01952 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 01953 01954 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 01955 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 01956 01957 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 01958 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 01959 ***/ 01960 FMAD2( xwt_00 , bold[ibase] , 01961 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 01962 01963 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 01964 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 01965 01966 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 01967 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 01968 01969 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 01970 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 01971 01972 /* interpolate in the y direction for each z */ 01973 /*** 01974 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 01975 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 01976 ***/ 01977 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 01978 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 01979 01980 /* interpolate in the z direction and 01981 put the result into the output array! 01982 (the +0.5 is to force rounding of the result) */ 01983 /*** 01984 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 01985 ***/ 01986 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 01987 FINAL( result , bslice[out_ind] ) ; 01988 out_ind++ ; 01989 01990 } 01991 } 01992 } 01993 break ; 01994 01995 case RESAM_CUBIC_TYPE:{ 01996 float xwt_m1,xwt_00,xwt_p1,xwt_p2 , /* interpolation weights */ 01997 ywt_m1,ywt_00,ywt_p1,ywt_p2 , 01998 zwt_m1,zwt_00,zwt_p1,zwt_p2 ; 01999 02000 INTYPE f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1 , /* interpolants */ 02001 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00 , 02002 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1 , 02003 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2 , 02004 f_km1 , f_k00 , f_kp1 , f_kp2 , result ; 02005 float frac_xi , frac_yj , frac_zk ; 02006 02007 int ibase , /* base index for interpolant */ 02008 in_jm1_km1 = -jstep- kstep , /* offsets for -1 (m1) */ 02009 in_jm1_k00 = -jstep , /* 0 (00) */ 02010 in_jm1_kp1 = -jstep+ kstep , /* +1 (p1) */ 02011 in_jm1_kp2 = -jstep+2*kstep , /* +2 (p2) */ 02012 in_j00_km1 = - kstep , /* steps in j and k indices */ 02013 in_j00_k00 = 0 , 02014 in_j00_kp1 = kstep , 02015 in_j00_kp2 = 2*kstep , 02016 in_jp1_km1 = jstep- kstep , 02017 in_jp1_k00 = jstep , 02018 in_jp1_kp1 = jstep+ kstep , 02019 in_jp1_kp2 =2*jstep+2*kstep , 02020 in_jp2_km1 =2*jstep- kstep , 02021 in_jp2_k00 =2*jstep , 02022 in_jp2_kp1 =2*jstep+ kstep , 02023 in_jp2_kp2 =2*jstep+2*kstep ; 02024 02025 int nxold1 = nxold-1 , nyold1 = nyold-1 , nzold1 = nzold-1 ; 02026 int nxold2 = nxold-2 , nyold2 = nyold-2 , nzold2 = nzold-2 ; 02027 02028 for( yj_new=yj_bot ; yj_new <= yj_top ; yj_new++ ){ 02029 02030 fxi_old = (fxi_base += dfxi_outer) ; 02031 fyj_old = (fyj_base += dfyj_outer) ; 02032 fzk_old = (fzk_base += dfzk_outer) ; 02033 02034 out_ind = xi_bot + yj_new * nxnew ; 02035 02036 for( xi_new=xi_bot ; xi_new <= xi_top ; xi_new++ ){ 02037 02038 fxi_old += dfxi_inner ; fyj_old += dfyj_inner ; fzk_old += dfzk_inner ; 02039 02040 /* check if outside old brick */ 02041 02042 if( fxi_old < fxi_bot || fxi_old > fxi_top || 02043 fyj_old < fyj_bot || fyj_old > fyj_top || 02044 fzk_old < fzk_bot || fzk_old > fzk_top ){ /* outside */ 02045 /*** 02046 bslice[out_ind++] = 0 ; 02047 ***/ 02048 FZERO(bslice[out_ind]) ; out_ind++ ; 02049 continue ; 02050 } 02051 02052 xi_old = FLOOR(fxi_old) ; frac_xi = fxi_old - xi_old ; 02053 yj_old = FLOOR(fyj_old) ; frac_yj = fyj_old - yj_old ; 02054 zk_old = FLOOR(fzk_old) ; frac_zk = fzk_old - zk_old ; 02055 02056 /* use NN if at very edges of old brick */ 02057 02058 if( xi_old == nxold1 || yj_old == nyold1 || zk_old == nzold1 || 02059 frac_xi < 0 || frac_yj < 0 || frac_zk < 0 ){ 02060 02061 bslice[out_ind++] = bold[ IBASE( ROUND(fxi_old) , 02062 ROUND(fyj_old) , 02063 ROUND(fzk_old) ) ] ; 02064 continue ; 02065 } 02066 02067 ibase = IBASE(xi_old,yj_old,zk_old) ; 02068 02069 /* use LINEAR if close to edges of old brick */ 02070 02071 if( xi_old == nxold2 || yj_old == nyold2 || zk_old == nzold2 || 02072 xi_old == 0 || yj_old == 0 || zk_old == 0 ){ 02073 02074 xwt_00 = LP_00(frac_xi) ; xwt_p1 = LP_P1(frac_xi) ; 02075 ywt_00 = LP_00(frac_yj) ; ywt_p1 = LP_P1(frac_yj) ; 02076 zwt_00 = LP_00(frac_zk) ; zwt_p1 = LP_P1(frac_zk) ; 02077 /*** 02078 f_j00_k00 = xwt_00 * bold[ibase] 02079 + xwt_p1 * bold[ibase+1] ; 02080 02081 f_jp1_k00 = xwt_00 * bold[ibase+in_jp1_k00] 02082 + xwt_p1 * bold[ibase+in_jp1_k00+1] ; 02083 02084 f_j00_kp1 = xwt_00 * bold[ibase+in_j00_kp1] 02085 + xwt_p1 * bold[ibase+in_j00_kp1+1] ; 02086 02087 f_jp1_kp1 = xwt_00 * bold[ibase+in_jp1_kp1] 02088 + xwt_p1 * bold[ibase+in_jp1_kp1+1] ; 02089 02090 f_k00 = ywt_00 * f_j00_k00 + ywt_p1 * f_jp1_k00 ; 02091 f_kp1 = ywt_00 * f_j00_kp1 + ywt_p1 * f_jp1_kp1 ; 02092 02093 bslice[out_ind++] = zwt_00 * f_k00 + zwt_p1 * f_kp1 + 0.5 ; 02094 ***/ 02095 FMAD2( xwt_00 , bold[ibase] , 02096 xwt_p1 , bold[ibase+1] , f_j00_k00 ) ; 02097 02098 FMAD2( xwt_00 , bold[ibase+in_jp1_k00] , 02099 xwt_p1 , bold[ibase+in_jp1_k00+1] , f_jp1_k00 ) ; 02100 02101 FMAD2( xwt_00 , bold[ibase+in_j00_kp1] , 02102 xwt_p1 , bold[ibase+in_j00_kp1+1] , f_j00_kp1 ) ; 02103 02104 FMAD2( xwt_00 , bold[ibase+in_jp1_kp1] , 02105 xwt_p1 , bold[ibase+in_jp1_kp1+1] , f_jp1_kp1 ) ; 02106 02107 FMAD2( ywt_00 , f_j00_k00 , ywt_p1 , f_jp1_k00 , f_k00 ) ; 02108 FMAD2( ywt_00 , f_j00_kp1 , ywt_p1 , f_jp1_kp1 , f_kp1 ) ; 02109 02110 FMAD2( zwt_00 , f_k00 , zwt_p1 , f_kp1 , result ) ; 02111 FINAL( result , bslice[out_ind] ) ; 02112 out_ind++ ; 02113 continue ; 02114 } 02115 02116 /* compute weights for CUBIC interpolation in each direction */ 02117 02118 xwt_m1 = CP_M1(frac_xi) ; xwt_00 = CP_00(frac_xi) ; 02119 xwt_p1 = CP_P1(frac_xi) ; xwt_p2 = CP_P2(frac_xi) ; 02120 02121 ywt_m1 = CP_M1(frac_yj) ; ywt_00 = CP_00(frac_yj) ; 02122 ywt_p1 = CP_P1(frac_yj) ; ywt_p2 = CP_P2(frac_yj) ; 02123 02124 zwt_m1 = CP_M1(frac_zk) ; zwt_00 = CP_00(frac_zk) ; 02125 zwt_p1 = CP_P1(frac_zk) ; zwt_p2 = CP_P2(frac_zk) ; 02126 02127 /* interpolate in the x direction for each y & z */ 02128 02129 CXINT(m1,m1,f_jm1_km1) ; CXINT(00,m1,f_j00_km1) ; 02130 CXINT(p1,m1,f_jp1_km1) ; CXINT(p2,m1,f_jp2_km1) ; 02131 02132 CXINT(m1,00,f_jm1_k00) ; CXINT(00,00,f_j00_k00) ; 02133 CXINT(p1,00,f_jp1_k00) ; CXINT(p2,00,f_jp2_k00) ; 02134 02135 CXINT(m1,p1,f_jm1_kp1) ; CXINT(00,p1,f_j00_kp1) ; 02136 CXINT(p1,p1,f_jp1_kp1) ; CXINT(p2,p1,f_jp2_kp1) ; 02137 02138 CXINT(m1,p2,f_jm1_kp2) ; CXINT(00,p2,f_j00_kp2) ; 02139 CXINT(p1,p2,f_jp1_kp2) ; CXINT(p2,p2,f_jp2_kp2) ; 02140 02141 /* interpolate in the y direction for each z */ 02142 /*** 02143 f_km1 = ywt_m1 * f_jm1_km1 + ywt_00 * f_j00_km1 02144 + ywt_p1 * f_jp1_km1 + ywt_p2 * f_jp2_km1 ; 02145 02146 f_k00 = ywt_m1 * f_jm1_k00 + ywt_00 * f_j00_k00 02147 + ywt_p1 * f_jp1_k00 + ywt_p2 * f_jp2_k00 ; 02148 02149 f_kp1 = ywt_m1 * f_jm1_kp1 + ywt_00 * f_j00_kp1 02150 + ywt_p1 * f_jp1_kp1 + ywt_p2 * f_jp2_kp1 ; 02151 02152 f_kp2 = ywt_m1 * f_jm1_kp2 + ywt_00 * f_j00_kp2 02153 + ywt_p1 * f_jp1_kp2 + ywt_p2 * f_jp2_kp2 ; 02154 ***/ 02155 FMAD4( ywt_m1 , f_jm1_km1 , ywt_00 , f_j00_km1 , 02156 ywt_p1 , f_jp1_km1 , ywt_p2 , f_jp2_km1 , f_km1 ) ; 02157 02158 FMAD4( ywt_m1 , f_jm1_k00 , ywt_00 , f_j00_k00 , 02159 ywt_p1 , f_jp1_k00 , ywt_p2 , f_jp2_k00 , f_k00 ) ; 02160 02161 FMAD4( ywt_m1 , f_jm1_kp1 , ywt_00 , f_j00_kp1 , 02162 ywt_p1 , f_jp1_kp1 , ywt_p2 , f_jp2_kp1 , f_kp1 ) ; 02163 02164 FMAD4( ywt_m1 , f_jm1_kp2 , ywt_00 , f_j00_kp2 , 02165 ywt_p1 , f_jp1_kp2 , ywt_p2 , f_jp2_kp2 , f_kp2 ) ; 02166 02167 /* interpolate in the z direction and 02168 put the result into the output array! 02169 (the +0.5 is to force rounding of the result) */ 02170 /*** 02171 bslice[out_ind++] = 0.5 + CP_FACTOR 02172 * ( zwt_m1 * f_km1 + zwt_00 * f_k00 02173 +zwt_p1 * f_kp1 + zwt_p2 * f_kp2 ) ; 02174 ***/ 02175 FMAD4( zwt_m1 , f_km1 , zwt_00 , f_k00 , 02176 zwt_p1 , f_kp1 , zwt_p2 , f_kp2 , result ) ; 02177 FSCAL(CP_FACTOR,result) ; 02178 FINAL( result , bslice[out_ind] ) ; 02179 out_ind++ ; 02180 02181 } /* end of inner loop */ 02182 } /* end of outer loop */ 02183 } 02184 break ; 02185 02186 } 02187 02188 EXRETURN ; 02189 } |
Variable Documentation
|
macros for a zero value * Definition at line 137 of file afni_slice.c. |
|
space for precomputed indices * Definition at line 228 of file afni_slice.c. |
|
Definition at line 229 of file afni_slice.c. |
|
Definition at line 139 of file afni_slice.c. |