00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #undef MAIN
00008 
00009 /*********************************************************************
00010    These are the template routines to fill in a slice of data
00011    in a new image array from an old image brick array, warped on demand.
00012 
00013    To create actual routines for this purpose, you must compile
00014    the file with the preprocessor symbol DTYPE set to one of
00015    the following types:
00016 
00017       byte short int float double complex rgbyte
00018 
00019       cc -c -DDTYPE=short afni_slice.c
00020       mv -f afni_slice.o afni_slice_short.o
00021 
00022    In the example above, the resulting routines will be named
00023 
00024       AFNI_lmap_to_xslice_short
00025       AFNI_lmap_to_yslice_short
00026       AFNI_lmap_to_zslice_short
00027       AFNI_br2sl_short
00028 
00029    and will take as input pointers bold and bslice, pointers of
00030    type "short *".
00031 *********************************************************************/
00032 
00033 #ifndef DTYPE
00034 #error "Cannot compile, since DTYPE is undefined."
00035 #endif
00036 
00037 #include "afni_warp.h"
00038 
00039 /** macros for function names defined in this file **/
00040 
00041 #define LMAP_XNAME TWO_TWO(AFNI_lmap_to_xslice_,DTYPE)
00042 #define LMAP_YNAME TWO_TWO(AFNI_lmap_to_yslice_,DTYPE)
00043 #define LMAP_ZNAME TWO_TWO(AFNI_lmap_to_zslice_,DTYPE)
00044 #define B2SL_NAME  TWO_TWO(AFNI_br2sl_,DTYPE)
00045 
00046 /*******************************************************************
00047      To add a new DTYPE, you must add definitions in mrilib.h,
00048      and define the following macros for the new type:
00049        FMAD2, FMAD4, FSCAL, FINAL, FZERO, and INTYPE
00050 ********************************************************************/
00051 
00052 /** macros for e = a*d1 + b*d2 (a,b floats; d1,d2 DTYPEs) **/
00053 
00054 #define FMAD2_short(a,d1,b,d2,e)   (e)=(a)*(d1)+(b)*(d2)
00055 #define FMAD2_float                FMAD2_short
00056 #define FMAD2_byte                 FMAD2_short
00057 #define FMAD2_int                  FMAD2_short
00058 #define FMAD2_double               FMAD2_short
00059 
00060 #define FMAD2_complex(a,d1,b,d2,e) ( (e).r = (a)*(d1).r + (b)*(d2).r, \
00061                                      (e).i = (a)*(d1).i + (b)*(d2).i   )
00062 
00063 #define FMAD2_rgbyte(a,d1,bb,d2,e) ( (e).r = (a)*(d1).r + (bb)*(d2).r, \
00064                                      (e).g = (a)*(d1).g + (bb)*(d2).g, \
00065                                      (e).b = (a)*(d1).b + (bb)*(d2).b   )
00066 #define FMAD2 TWO_TWO(FMAD2_,DTYPE)
00067 
00068 /** macros for e = a*d1 + b*d2 + c*d3 + d*d3 (a-d floats; d1-d4 DTYPEs) **/
00069 
00070 #define FMAD4_short(a,d1,b,d2,c,d3,d,d4,e)   (e)=(a)*(d1)+(b)*(d2)+(c)*(d3)+(d)*(d4)
00071 #define FMAD4_float                          FMAD4_short
00072 #define FMAD4_byte                           FMAD4_short
00073 #define FMAD4_int                            FMAD4_short
00074 #define FMAD4_double                         FMAD4_short
00075 
00076 #define FMAD4_complex(a,d1,b,d2,c,d3,d,d4,e)                              \
00077              ( (e).r = (a)*(d1).r + (b)*(d2).r + (c)*(d3).r + (d)*(d4).r, \
00078                (e).i = (a)*(d1).i + (b)*(d2).i + (c)*(d3).i + (d)*(d4).i   )
00079 
00080 #define FMAD4_rgbyte(a,d1,bb,d2,c,d3,d,d4,e)                               \
00081              ( (e).r = (a)*(d1).r + (bb)*(d2).r + (c)*(d3).r + (d)*(d4).r, \
00082                (e).g = (a)*(d1).g + (bb)*(d2).g + (c)*(d3).g + (d)*(d4).g, \
00083                (e).b = (a)*(d1).b + (bb)*(d2).b + (c)*(d3).b + (d)*(d4).b   )
00084 
00085 #define FMAD4 TWO_TWO(FMAD4_,DTYPE)
00086 
00087 /** macros to multiply float a times DTYPE b and store the result in b again **/
00088 
00089 #define FSCAL_short(a,b)           (b)*=(a)
00090 #define FSCAL_float                FSCAL_short
00091 #define FSCAL_byte                 FSCAL_short
00092 #define FSCAL_int                  FSCAL_short
00093 #define FSCAL_double               FSCAL_short
00094 #define FSCAL_complex(a,b)         ( (b).r *= (a) , (b).i *= (a) )
00095 
00096 #define FSCAL_rgbyte(a,bb)         ( (bb).r*= (a) , (bb).g*= (a) , (bb).b*= (a) )
00097 #define FSCAL TWO_TWO(FSCAL_,DTYPE)
00098 
00099 /** macros for assigning final result from INTYPE a to DTYPE b **/
00100 
00101    /* 18 Nov 1998: modify FINAL_short and FINAL_byte to prevent overflow */
00102 
00103 #define CLIP_OVERFLOW
00104 #ifdef  CLIP_OVERFLOW
00105 
00106 #  define FINAL_short(a,b) (b) = ( ((a)<-32767.0) ? (-32767) : \
00107                                    ((a)> 32767.0) ? ( 32767) : ((short)((a)+0.5)) )
00108 
00109 #  define FINAL_byte(a,b)  (b) = ( ((a)<   0.0) ? (0)   : \
00110                                    ((a)> 255.0) ? (255) : ((byte)((a)+0.5)) )
00111 
00112 #else
00113 # define FINAL_short(a,b)           (b)=((short)((a)+0.5))
00114 # define FINAL_byte(a,b)            (b)=((byte)((a)+0.5))
00115 #endif
00116 
00117 #define FINAL_int(a,b)             (b)=((int)((a)+0.5))
00118 #define FINAL_float(a,b)           (b)=(a)
00119 #define FINAL_double               FINAL_float
00120 #define FINAL_complex              FINAL_float
00121 #define FINAL_rgbyte               FINAL_float
00122 #define FINAL TWO_TWO(FINAL_,DTYPE)
00123 
00124 /** macros for putting a zero into DTYPE b **/
00125 
00126 #define FZERO_short(b)             (b)=0
00127 #define FZERO_byte                 FZERO_short
00128 #define FZERO_int                  FZERO_short
00129 #define FZERO_float(b)             (b)=0.0
00130 #define FZERO_double               FZERO_float
00131 #define FZERO_complex(b)           ( (b).r = 0.0 , (b).i = 0.0 )
00132 #define FZERO_rgbyte(bb)           ( (bb).r=(bb).g=(bb).g = 0 )
00133 #define FZERO TWO_TWO(FZERO_,DTYPE)
00134 
00135 /** macros for a zero value **/
00136 
00137 static complex complex_zero = { 0.0,0.0 } ;
00138 
00139 static rgbyte  rgbyte_zero  = { 0,0,0 } ;
00140 
00141 #define ZERO_short    0
00142 #define ZERO_byte     0
00143 #define ZERO_int      0
00144 #define ZERO_float    0.0
00145 #define ZERO_double   0.0
00146 #define ZERO_complex  complex_zero
00147 #define ZERO_rgbyte   rgbyte_zero
00148 #define ZERO          TWO_TWO(ZERO_,DTYPE)
00149 
00150 /** macros for intermediate interpolants data type **/
00151 
00152 #define INTYPE_short    float
00153 #define INTYPE_float    float
00154 #define INTYPE_byte     float
00155 #define INTYPE_int      float
00156 #define INTYPE_double   double
00157 #define INTYPE_complex  complex
00158 #define INTYPE_rgbyte   rgbyte
00159 #define INTYPE TWO_TWO(INTYPE_,DTYPE)
00160 
00161 /**-------------------------------------------------------------------------**/
00162 /**----- macros for quickly doing NN interpolation along parallel axes -----**/
00163 
00164 /** test and flags for which steps are zero **/
00165 
00166 #define ZZZ           1.e-5  /* effectively zero */
00167 #define NONE_ZERO     0
00168 #define X_ZERO        1
00169 #define Y_ZERO        2
00170 #define XY_ZERO       3
00171 #define Z_ZERO        4
00172 #define XZ_ZERO       5
00173 #define YZ_ZERO       6
00174 #define XYZ_ZERO      7
00175 #define OUTADD      100
00176 #define THREEZ(x,y,z) ((fabs(x)<ZZZ) + 2*(fabs(y)<ZZZ) + 4*(fabs(z)<ZZZ))
00177 
00178 /** ALOOP: inner loop statements for the actual NN assignments,
00179              and there is no need to check for being outside the input brick
00180     _GEN:  for the general (non-parallel) case
00181     _PAR:  for the case when the inner loop is parallel to an input brick axis **/
00182 
00183 #define NN_ALOOP_GEN \
00184  (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \
00185   xi_old = FLOOR(fxi_old) , yj_old = FLOOR(fyj_old) , zk_old = FLOOR(fzk_old) , \
00186   bslice[out_ind++] = bold[ IBASE(xi_old,yj_old,zk_old) ])
00187 
00188 #define NN_ALOOP_PAR(ijk) (bslice[out_ind++] = bold[ ib[ijk]+ob ])
00189 
00190 /** BLOOP: assign values to the ib array that will be used to
00191            index into the input brick for rapid access;
00192            there is one BLOOP for each possible parallel axis **/
00193 
00194 #define NN_BLOOP_XY_ZERO(ijk) \
00195  (fzk_old += dfzk_inner , zk_old = GFLOOR(fzk_old) , ib[ijk] = IBASE(0,0,zk_old))
00196 
00197 #define NN_BLOOP_XZ_ZERO(ijk) \
00198  (fyj_old += dfyj_inner , yj_old = GFLOOR(fyj_old) , ib[ijk] = IBASE(0,yj_old,0))
00199 
00200 #define NN_BLOOP_YZ_ZERO(ijk) \
00201  (fxi_old += dfxi_inner , xi_old = GFLOOR(fxi_old) , ib[ijk] = IBASE(xi_old,0,0))
00202 
00203 /** macros to test if the point we want to resample from is outside the old array **/
00204 
00205 #define TEST_OUT_XXX ( fxi_old < fxi_bot || fxi_old > fxi_top )
00206 #define TEST_OUT_YYY ( fyj_old < fyj_bot || fyj_old > fyj_top )
00207 #define TEST_OUT_ZZZ ( fzk_old < fzk_bot || fzk_old > fzk_top )
00208 
00209 #define TEST_OUT_ALL (TEST_OUT_XXX || TEST_OUT_YYY || TEST_OUT_ZZZ)
00210 
00211 /** CLOOP: like the ALOOP, but for when we must test
00212            if the desired point may be outside the input brick array **/
00213 
00214 #define NN_CLOOP_GEN                                                      \
00215  (fxi_old += dfxi_inner , fyj_old += dfyj_inner , fzk_old += dfzk_inner , \
00216   bslice[out_ind++] = (TEST_OUT_ALL) ? ZERO :                             \
00217                       bold[IBASE(FLOOR(fxi_old),FLOOR(fyj_old),FLOOR(fzk_old))] )
00218 
00219 #define NN_CLOOP_PAR(ijk) \
00220  (bslice[out_ind++] = (ib[ijk]<0 || ib[ijk]>=ub) ? ZERO : bold[ib[ijk]+ob])
00221 
00222 /** ZLOOP: just assign zero to each output **/
00223 
00224 #define NN_ZLOOP (bslice[out_ind++] = ZERO)
00225 
00226 /** space for precomputed indices **/
00227 
00228 static int * ib = NULL ;
00229 static int  nib = -1 ;
00230 
00231 /** macro to make the ib array as big as we need it **/
00232 
00233 #define MAKE_IBIG(top) do{ if(nib < (top)){                                 \
00234                               if(ib != NULL) free(ib) ;                     \
00235                               ib  = (int *) malloc(sizeof(int)*((top)+9)) ; \
00236                               if(ib==NULL){                                 \
00237                                  fprintf(stderr,"\nmalloc fails in NN reslice!\n");EXIT(1);} \
00238                               nib = (top) ; } } while(0)
00239 
00240 /*---------------------------------------------------------------------
00241    routine to apply a linear mapping to a dataset and fill in
00242    a single fixed-x slice from a 3D dataset:
00243 
00244      map        = linear mapping from old to new
00245                    (bot & top field contain index limits in new)
00246      resam_mode = type of interpolation to do in old
00247 
00248      old_daxes  = axis information of old dataset
00249      bold       = pointer to old dataset's 3D data brick array
00250 
00251      new_daxes  = axis information of new dataset
00252      xi_fix     = fixed x index in new dataset
00253      bslice     = pointer to nynew * nznew array to hold slice
00254 
00255   Cognate routines for yslice (output is nznew * nxnew: Y-ZX) and
00256                        zslize (output is nxnew * nynew: Z-XY)
00257   follow this one directly.
00258 -----------------------------------------------------------------------*/
00259 
00260 #define IBASE(i,j,k) ((i)+(j)*jstep+(k)*kstep)
00261 #define ROUND(qq)    ((int)(qq+0.5))
00262 #define FLOOR(qq)    ((int)(qq))          /* cheap and fast */
00263 #define GFLOOR(qq)   ((int)floor(qq))     /* good and slow */
00264 
00265 /* define linear interpolation polynomials */
00266 
00267 #define LP_00(x) (1.0-(x))
00268 #define LP_P1(x) (x)
00269 
00270 /* define blocky interpolation functions */
00271 
00272 #define BP_hh(x) (8.0*((x)*(x))*((x)*(x)))
00273 #define BP_00(x) ( ((x)<0.5) ? (1-BP_hh(x)) : (  BP_hh(1-(x))) )
00274 #define BP_P1(x) ( ((x)<0.5) ? (  BP_hh(x)) : (1-BP_hh(1-(x))) )
00275 
00276 /* define cubic interpolation polynomials */
00277 
00278 #define CP_M1(x)  (-(x)*((x)-1)*((x)-2))
00279 #define CP_00(x)  (3.0*((x)+1)*((x)-1)*((x)-2))
00280 #define CP_P1(x)  (-3.0*(x)*((x)+1)*((x)-2))
00281 #define CP_P2(x)  ((x)*((x)+1)*((x)-1))
00282 #define CP_FACTOR  4.62962963e-3   /* 1/216 = final scaling factor */
00283 
00284 #define FXYZTMP(xx,yy,zz)                            \
00285        ( fxi_tmp =   mt.mat[0][0] * xx               \
00286                    + mt.mat[0][1] * yy               \
00287                    + mt.mat[0][2] * zz - vt.xyz[0] , \
00288          fyj_tmp =   mt.mat[1][0] * xx               \
00289                    + mt.mat[1][1] * yy               \
00290                    + mt.mat[1][2] * zz - vt.xyz[1] , \
00291          fzk_tmp =   mt.mat[2][0] * xx               \
00292                    + mt.mat[2][1] * yy               \
00293                    + mt.mat[2][2] * zz - vt.xyz[2] )
00294 
00295 
00296 void LMAP_XNAME( THD_linear_mapping * map , int resam_mode ,
00297                  THD_dataxes * old_daxes , DTYPE * bold ,
00298                  THD_dataxes * new_daxes , int xi_fix , DTYPE * bslice )
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 }
00939 
00940 /*--------------------------------------------------------------------------*/
00941 
00942 void LMAP_YNAME( THD_linear_mapping * map , int resam_mode ,
00943                  THD_dataxes * old_daxes , DTYPE * bold ,
00944                  THD_dataxes * new_daxes , int yj_fix , DTYPE * bslice )
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 }
01564 
01565 /*----------------------------------------------------------------------------*/
01566 
01567 void LMAP_ZNAME( THD_linear_mapping * map , int resam_mode ,
01568                  THD_dataxes * old_daxes , DTYPE * bold ,
01569                  THD_dataxes * new_daxes , int zk_fix , DTYPE * bslice )
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 }
02190 
02191 /*---------------------------------------------------------------------
02192    Routine to copy data from a brick array into a slice:
02193      nxx,nyy,nzz = brick dimensions
02194      fixed_axis  = 1, 2, or 3 (x, y, or z)
02195      fixed_index = subscript of chosen slice in the fixed_axis direction
02196      bold        = pointer to brick array
02197      bslice      = pointer to area to get slice from brick array
02198                      fixed_axis = 1 --> nyy * nzz
02199                                   2 --> nzz * nxx
02200                                   3 --> nxx * nyy
02201 -----------------------------------------------------------------------*/
02202 
02203 void B2SL_NAME( int nxx, int nyy, int nzz ,
02204                 int fixed_axis , int fixed_index , DTYPE * bold , DTYPE * bslice )
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 }