/* * $Id: spm_vol_utils.c 4452 2011-09-02 10:45:26Z guillaume $ * John Ashburner */ #define TINY 5e-2 #ifdef SPM_UNSIGNED_CHAR #define RESAMPLE resample_uchar #define RESAMPLE_D resample_d_uchar #define SLICE slice_uchar #define RESAMPLE_0 resample_uchar_0 #define RESAMPLE_1 resample_uchar_1 #define RESAMPLE_D_1 resample_d_uchar_1 #define RESAMPLE_POLY resample_uchar_poly #define RESAMPLE_D_POLY resample_d_uchar_poly #define SLICE_0 slice_uchar_0 #define SLICE_1 slice_uchar_1 #define SLICE_POLY slice_uchar_poly #define PLANE plane_uchar #define GET(x) (x) #define IMAGE_DTYPE unsigned char #endif #ifdef SPM_SIGNED_CHAR #define RESAMPLE resample_schar #define RESAMPLE_D resample_d_schar #define SLICE slice_schar #define RESAMPLE_0 resample_schar_0 #define RESAMPLE_1 resample_schar_1 #define RESAMPLE_D_1 resample_d_schar_1 #define RESAMPLE_POLY resample_schar_poly #define RESAMPLE_D_POLY resample_d_schar_poly #define SLICE_0 slice_schar_0 #define SLICE_1 slice_schar_1 #define SLICE_POLY slice_schar_poly #define PLANE plane_schar #define GET(x) (x) #define IMAGE_DTYPE signed char #endif #ifdef SPM_SIGNED_SHORT #ifdef SPM_BYTESWAP #define GET(x) getshort(x) #define RESAMPLE resample_short_s #define RESAMPLE_D resample_d_short_s #define SLICE slice_short_s #define RESAMPLE_0 resample_short_s_0 #define RESAMPLE_1 resample_short_s_1 #define RESAMPLE_D_1 resample_d_short_s_1 #define RESAMPLE_POLY resample_short_s_poly #define RESAMPLE_D_POLY resample_d_short_s_poly #define SLICE_0 slice_short_s_0 #define SLICE_1 slice_short_s_1 #define SLICE_POLY slice_short_s_poly #define PLANE plane_short_s #else #define GET(x) (x) #define RESAMPLE resample_short #define RESAMPLE_D resample_d_short #define SLICE slice_short #define RESAMPLE_0 resample_short_0 #define RESAMPLE_1 resample_short_1 #define RESAMPLE_D_1 resample_d_short_1 #define RESAMPLE_POLY resample_short_poly #define RESAMPLE_D_POLY resample_d_short_poly #define SLICE_0 slice_short_0 #define SLICE_1 slice_short_1 #define SLICE_POLY slice_short_poly #define PLANE plane_short #endif #define IMAGE_DTYPE short int #endif #ifdef SPM_UNSIGNED_SHORT #ifdef SPM_BYTESWAP #define GET(x) getushort(x) #define RESAMPLE resample_ushort_s #define RESAMPLE_D resample_d_ushort_s #define SLICE slice_ushort_s #define RESAMPLE_0 resample_ushort_s_0 #define RESAMPLE_1 resample_ushort_s_1 #define RESAMPLE_D_1 resample_d_ushort_s_1 #define RESAMPLE_POLY resample_ushort_s_poly #define RESAMPLE_D_POLY resample_d_ushort_s_poly #define SLICE_0 slice_ushort_s_0 #define SLICE_1 slice_ushort_s_1 #define SLICE_POLY slice_ushort_s_poly #define PLANE plane_ushort_s #else #define GET(x) (x) #define RESAMPLE resample_ushort #define RESAMPLE_D resample_d_ushort #define SLICE slice_ushort #define RESAMPLE_0 resample_ushort_0 #define RESAMPLE_1 resample_ushort_1 #define RESAMPLE_D_1 resample_d_ushort_1 #define RESAMPLE_POLY resample_ushort_poly #define RESAMPLE_D_POLY resample_d_ushort_poly #define SLICE_0 slice_ushort_0 #define SLICE_1 slice_ushort_1 #define SLICE_POLY slice_ushort_poly #define PLANE plane_ushort #endif #define IMAGE_DTYPE unsigned short int #endif #ifdef SPM_SIGNED_INT #ifdef SPM_BYTESWAP #define GET(x) getint(x) #define RESAMPLE resample_int_s #define RESAMPLE_D resample_d_int_s #define SLICE slice_int_s #define RESAMPLE_0 resample_int_s_0 #define RESAMPLE_1 resample_int_s_1 #define RESAMPLE_D_1 resample_d_int_s_1 #define RESAMPLE_POLY resample_int_s_poly #define RESAMPLE_D_POLY resample_d_int_s_poly #define SLICE_0 slice_int_s_0 #define SLICE_1 slice_int_s_1 #define SLICE_POLY slice_int_s_poly #define PLANE plane_int_s #else #define GET(x) (x) #define RESAMPLE resample_int #define RESAMPLE_D resample_d_int #define SLICE slice_int #define RESAMPLE_0 resample_int_0 #define RESAMPLE_1 resample_int_1 #define RESAMPLE_D_1 resample_d_int_1 #define RESAMPLE_POLY resample_int_poly #define RESAMPLE_D_POLY resample_d_int_poly #define SLICE_0 slice_int_0 #define SLICE_1 slice_int_1 #define SLICE_POLY slice_int_poly #define PLANE plane_int #endif #define IMAGE_DTYPE int #endif #ifdef SPM_UNSIGNED_INT #ifdef SPM_BYTESWAP #define GET(x) getuint(x) #define RESAMPLE resample_uint_s #define RESAMPLE_D resample_d_uint_s #define SLICE slice_uint_s #define RESAMPLE_0 resample_uint_s_0 #define RESAMPLE_1 resample_uint_s_1 #define RESAMPLE_D_1 resample_d_uint_s_1 #define RESAMPLE_POLY resample_uint_s_poly #define RESAMPLE_D_POLY resample_d_uint_s_poly #define SLICE_0 slice_uint_s_0 #define SLICE_1 slice_uint_s_1 #define SLICE_POLY slice_uint_s_poly #define PLANE plane_uint_s #else #define GET(x) (x) #define RESAMPLE resample_uint #define RESAMPLE_D resample_d_uint #define SLICE slice_uint #define RESAMPLE_0 resample_uint_0 #define RESAMPLE_1 resample_uint_1 #define RESAMPLE_D_1 resample_d_uint_1 #define RESAMPLE_POLY resample_uint_poly #define RESAMPLE_D_POLY resample_d_uint_poly #define SLICE_0 slice_uint_0 #define SLICE_1 slice_uint_1 #define SLICE_POLY slice_uint_poly #define PLANE plane_uint #endif #define IMAGE_DTYPE unsigned int #endif #ifdef SPM_FLOAT #ifdef SPM_BYTESWAP #define GET(x) getfloat(x) #define RESAMPLE resample_float_s #define RESAMPLE_D resample_d_float_s #define SLICE slice_float_s #define RESAMPLE_0 resample_float_s_0 #define RESAMPLE_1 resample_float_s_1 #define RESAMPLE_D_1 resample_d_float_s_1 #define RESAMPLE_POLY resample_float_s_poly #define RESAMPLE_D_POLY resample_d_float_s_poly #define SLICE_0 slice_float_s_0 #define SLICE_1 slice_float_s_1 #define SLICE_POLY slice_float_s_poly #define PLANE plane_float_s #else #define GET(x) (x) #define RESAMPLE resample_float #define RESAMPLE_D resample_d_float #define SLICE slice_float #define RESAMPLE_0 resample_float_0 #define RESAMPLE_1 resample_float_1 #define RESAMPLE_D_1 resample_d_float_1 #define RESAMPLE_POLY resample_float_poly #define RESAMPLE_D_POLY resample_d_float_poly #define SLICE_0 slice_float_0 #define SLICE_1 slice_float_1 #define SLICE_POLY slice_float_poly #define PLANE plane_float #endif #define IMAGE_DTYPE float #endif #ifdef SPM_DOUBLE #ifdef SPM_BYTESWAP #define GET(x) getdouble(x) #define RESAMPLE resample_double_s #define RESAMPLE_D resample_d_double_s #define SLICE slice_double_s #define RESAMPLE_0 resample_double_s_0 #define RESAMPLE_1 resample_double_s_1 #define RESAMPLE_D_1 resample_d_double_s_1 #define RESAMPLE_POLY resample_double_s_poly #define RESAMPLE_D_POLY resample_d_double_s_poly #define SLICE_0 slice_double_s_0 #define SLICE_1 slice_double_s_1 #define SLICE_POLY slice_double_s_poly #define PLANE plane_double_s #else #define GET(x) (x) #define RESAMPLE resample_double #define RESAMPLE_D resample_d_double #define SLICE slice_double #define RESAMPLE_0 resample_double_0 #define RESAMPLE_1 resample_double_1 #define RESAMPLE_D_1 resample_d_double_1 #define RESAMPLE_POLY resample_double_poly #define RESAMPLE_D_POLY resample_d_double_poly #define SLICE_0 slice_double_0 #define SLICE_1 slice_double_1 #define SLICE_POLY slice_double_poly #define PLANE plane_double #endif #define IMAGE_DTYPE double #endif #include #include #define RINT(A) floor((A)+0.5) #include "spm_make_lookup.h" #include "spm_getdata.h" static void (*make_lookup)() = make_lookup_poly, (*make_lookup_grad)() = make_lookup_poly_grad; /* Zero order hold resampling - nearest neighbour */ void RESAMPLE_0(m,vol,out,x,y,z,xdim,ydim,zdim,background, scale,offset) int m, xdim,ydim,zdim; double out[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { int i; for (i=0; i=0 && xcoord=0 && ycoord=0 && zcoord=-TINY && zi=-TINY && yi=-TINY && xi=xdim-1) ? ((offx=0),xdim-1) : ((offx=1 ),xcoord)); ycoord = (ycoord < 0) ? ((offy=0),0) : ((ycoord>=ydim-1) ? ((offy=0),ydim-1) : ((offy=xdim),ycoord)); zcoord = (zcoord < 0) ? ((offz=0),0) : ((zcoord>=zdim-1) ? ((offz=0),zdim-1) : ((offz=1 ),zcoord)); off1 = xcoord + xdim*ycoord; off2 = off1+offy; k222 = GET(vol[zcoord ][off1]); k122 = GET(vol[zcoord ][off1+offx]); k212 = GET(vol[zcoord ][off2]); k112 = GET(vol[zcoord ][off2+offx]); k221 = GET(vol[zcoord+offz][off1]); k121 = GET(vol[zcoord+offz][off1+offx]); k211 = GET(vol[zcoord+offz][off2]); k111 = GET(vol[zcoord+offz][off2+offx]); /* resampled pixel value (trilinear interpolation) */ out[i] = (((k222*dx2 + k122*dx1)*dy2 + (k212*dx2 + k112*dx1)*dy1)*scale[zcoord ] + offset[zcoord ])*dz2 + (((k221*dx2 + k121*dx1)*dy2 + (k211*dx2 + k111*dx1)*dy1)*scale[zcoord+offz] + offset[zcoord+offz])*dz1; } else out[i] = background; } } /* First order hold resampling - trilinear interpolation */ void RESAMPLE_D_1(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim,background, scale,offset) int m, xdim,ydim,zdim; double out[],gradx[],grady[],gradz[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { int i; for (i=0; i=-TINY && zi=-TINY && yi=-TINY && xi=xdim-1) ? ((offx=0),xdim-1) : ((offx=1 ),xcoord)); ycoord = (ycoord < 0) ? ((offy=0),0) : ((ycoord>=ydim-1) ? ((offy=0),ydim-1) : ((offy=xdim),ycoord)); zcoord = (zcoord < 0) ? ((offz=0),0) : ((zcoord>=zdim-1) ? ((offz=0),zdim-1) : ((offz=1 ),zcoord)); off1 = xcoord + xdim*ycoord; off2 = off1+offy; k222 = GET(vol[zcoord ][off1]); k122 = GET(vol[zcoord ][off1+offx]); k212 = GET(vol[zcoord ][off2]); k112 = GET(vol[zcoord ][off2+offx]); k221 = GET(vol[zcoord+offz][off1]); k121 = GET(vol[zcoord+offz][off1+offx]); k211 = GET(vol[zcoord+offz][off2]); k111 = GET(vol[zcoord+offz][off2+offx]); /* resampled pixel value (trilinear interpolation) and gradients old code: gradx[i] = scale*(((k111 - k211 )*dy1 + (k121 - k221 )*dy2)*dz1 + ((k112 - k212 )*dy1 + (k122 - k222 )*dy2)*dz2); grady[i] = scale*(((k111*dx1 + k211*dx2) - (k121*dx1 + k221*dx2) )*dz1 + ((k112*dx1 + k212*dx2) - (k122*dx1 + k222*dx2) )*dz2); gradz[i] = scale*(((k111*dx1 + k211*dx2)*dy1 + (k121*dx1 + k221*dx2)*dy2) - ((k112*dx1 + k212*dx2)*dy1 + (k122*dx1 + k222*dx2)*dy2) ); out[i] = scale*(((k111*dx1 + k211*dx2)*dy1 + (k121*dx1 + k221*dx2)*dy2)*dz1 + ((k112*dx1 + k212*dx2)*dy1 + (k122*dx1 + k222*dx2)*dy2)*dz2) + offset; */ gradx[i] = (((k111 - k211)*dy1 + (k121 - k221)*dy2)*scale[zcoord+offz])*dz1 + (((k112 - k212)*dy1 + (k122 - k222)*dy2)*scale[zcoord ])*dz2; k111 = (k111*dx1 + k211*dx2)*scale[zcoord+offz]+offset[zcoord+offz]; k121 = (k121*dx1 + k221*dx2)*scale[zcoord+offz]+offset[zcoord+offz]; k112 = (k112*dx1 + k212*dx2)*scale[zcoord ]+offset[zcoord ]; k122 = (k122*dx1 + k222*dx2)*scale[zcoord ]+offset[zcoord ]; grady[i] = (k111 - k121)*dz1 + (k112 - k122)*dz2; k111 = k111*dy1 + k121*dy2; k112 = k112*dy1 + k122*dy2; gradz[i] = k111 - k112; out[i] = k111*dz1 + k112*dz2; } else { out[i] = background; gradx[i] = 0.0; grady[i] = 0.0; gradz[i] = 0.0; } } } /* Sinc resampling */ void RESAMPLE_POLY(m,vol,out,x,y,z,xdim,ydim,zdim, q,background, scale,offset) int m, xdim,ydim,zdim, q; double out[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { int i; int dx1, dy1, dz1; static double tablex[255], tabley[255], tablez[255]; for (i=0; i=1-TINY && z[i]<=zdim+TINY && y[i]>=1-TINY && y[i]<=ydim+TINY && x[i]>=1-TINY && x[i]<=xdim+TINY) { double dat=0.0, *tp1, *tp1end, *tp2end, *tp3end; make_lookup(x[i], q, xdim, &dx1, tablex, &tp3end); make_lookup(y[i], q, ydim, &dy1, tabley, &tp2end); make_lookup(z[i], q, zdim, &dz1, tablez, &tp1end); tp1 = tablez; dy1 *= xdim; while(tp1 <= tp1end) { IMAGE_DTYPE *dp2 = &vol[dz1][dy1]; double dat2 = 0.0, *tp2 = tabley; while (tp2 <= tp2end) { register double dat3 = 0.0, *tp3 = tablex; register IMAGE_DTYPE *dp3 = dp2 + dx1; while(tp3 <= tp3end) dat3 += GET(*(dp3++)) * *(tp3++); dat2 += dat3 * *(tp2++); dp2 += xdim; } dat += (dat2*scale[dz1]+offset[dz1]) * *(tp1++); dz1 ++; } out[i] = dat; } else out[i] = background; } } /* Sinc resampling */ void RESAMPLE_D_POLY(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim, q,background, scale,offset) int m, xdim,ydim,zdim, q; double out[],gradx[],grady[],gradz[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { int i; int dx1, dy1, dz1; static double tablex[255], tabley[255], tablez[255]; static double dtablex[255], dtabley[255], dtablez[255]; for (i=0; i=1-TINY && z[i]<=zdim+TINY && y[i]>=1-TINY && y[i]<=ydim+TINY && x[i]>=1-TINY && x[i]<=xdim+TINY) { double dat=0.0, datx = 0.0, daty = 0.0, datz = 0.0, *tp1, *tp1end, *tp2end, *tp3end, *dp1; make_lookup_grad(x[i], q, xdim, &dx1, tablex, dtablex, &tp3end); make_lookup_grad(y[i], q, ydim, &dy1, tabley, dtabley, &tp2end); make_lookup_grad(z[i], q, zdim, &dz1, tablez, dtablez, &tp1end); tp1 = tablez; dp1 = dtablez; dy1 *= xdim; while(tp1 <= tp1end) { IMAGE_DTYPE *d2 = &vol[dz1][dy1]; double dat2 = 0.0, *tp2 = tabley; double dat2x = 0.0, *dp2 = dtabley, dat2y = 0.0; while (tp2 <= tp2end) { register IMAGE_DTYPE *d3 = d2 + dx1; register double dat3 = 0.0, *tp3 = tablex; register double dat3x = 0.0, *dp3 = dtablex; while(tp3 <= tp3end) { dat3x += GET(*(d3 )) * *(dp3++); dat3 += GET(*(d3++)) * *(tp3++); } dat2x += dat3x * *(tp2 ); dat2 += dat3 * *(tp2++); dat2y += dat3 * *(dp2++); d2 += xdim; } datx += (dat2x*scale[dz1]) * *(tp1 ); daty += (dat2y*scale[dz1]) * *(tp1 ); dat2 = dat2*scale[dz1]+offset[dz1]; dat += dat2 * *(tp1++); datz += dat2 * *(dp1++); dz1 ++; } out[i] = dat; gradx[i] = datx; grady[i] = daty; gradz[i] = datz; } else { out[i] = background; gradx[i] = 0.0; grady[i] = 0.0; gradz[i] = 0.0; } } } /* Zero order hold resampling - nearest neighbour */ int SLICE_0(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, background, scale,offset) double mat[16], background, scale[],offset[]; double image[]; IMAGE_DTYPE *vol[]; int xdim1, ydim1, xdim2, ydim2, zdim2; { double y, x2, y2, z2, s2, dx3=mat[0], dy3=mat[1], dz3=mat[2], ds3=mat[3]; int t = 0; x2 = mat[12] + 0*mat[8]; y2 = mat[13] + 0*mat[9]; z2 = mat[14] + 0*mat[10]; s2 = mat[15] + 0*mat[11]; for(y=1; y<=ydim1; y++) { double x; double x3 = x2 + y*mat[4]; double y3 = y2 + y*mat[5]; double z3 = z2 + y*mat[6]; double s3 = s2 + y*mat[7]; for(x=1; x<=xdim1; x++) { int ix4, iy4, iz4; s3 += ds3; if (s3 == 0.0) return(-1); ix4 = floor(((x3 += dx3)/s3)-0.5); iy4 = floor(((y3 += dy3)/s3)-0.5); iz4 = floor(((z3 += dz3)/s3)-0.5); if (iz4>=0 && iz4=0 && iy4=0 && ix4=-TINY && z4=-TINY && y4=-TINY && x4=xdim2-1) ? ((offx=0),xdim2-1) : ((offx=1 ),ix4)); iy4 = (iy4 < 0) ? ((offy=0),0) : ((iy4>=ydim2-1) ? ((offy=0),ydim2-1) : ((offy=xdim2),iy4)); iz4 = (iz4 < 0) ? ((offz=0),0) : ((iz4>=zdim2-1) ? ((offz=0),zdim2-1) : ((offz=1 ),iz4)); off1 = ix4 + xdim2*iy4; off2 = off1+offy; k222 = GET(vol[iz4 ][off1]); k122 = GET(vol[iz4 ][off1+offx]); k212 = GET(vol[iz4 ][off2]); k112 = GET(vol[iz4 ][off2+offx]); k221 = GET(vol[iz4+offz][off1]); k121 = GET(vol[iz4+offz][off1+offx]); k211 = GET(vol[iz4+offz][off2]); k111 = GET(vol[iz4+offz][off2+offx]); /* resampled pixel value (trilinear interpolation) */ image[t] = (((k222*dx2 + k122*dx1)*dy2 + (k212*dx2 + k112*dx1)*dy1)*scale[iz4 ] + offset[iz4 ])*dz2 + (((k221*dx2 + k121*dx1)*dy2 + (k211*dx2 + k111*dx1)*dy1)*scale[iz4+offz] + offset[iz4+offz])*dz1; } else image[t] = background; t++; } } return(0); } /* Sinc resampling */ int SLICE_POLY(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, q, background, scale,offset) int ydim1,xdim1, xdim2,ydim2,zdim2, q; double image[], mat[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { int dx1, dy1, dz1; static double tablex[255], tabley[255], tablez[255]; double y, x2, y2, z2, s2, dx3=mat[0], dy3=mat[1], dz3=mat[2], ds3=mat[3]; x2 = mat[12] + 0*mat[8]; y2 = mat[13] + 0*mat[9]; z2 = mat[14] + 0*mat[10]; s2 = mat[15] + 0*mat[11]; for(y=1; y<=ydim1; y++) { double x; double x3 = x2 + y*mat[4]; double y3 = y2 + y*mat[5]; double z3 = z2 + y*mat[6]; double s3 = s2 + y*mat[7]; for(x=1; x<=xdim1; x++) { double x4,y4,z4; s3 += ds3; if (s3 == 0.0) return(-1); x4=(x3 += dx3)/s3; y4=(y3 += dy3)/s3; z4=(z3 += dz3)/s3; if (z4>=1-TINY && z4<=zdim2+TINY && y4>=1-TINY && y4<=ydim2+TINY && x4>=1-TINY && x4<=xdim2+TINY) { double dat=0.0, *tp1, *tp1end, *tp2end, *tp3end; make_lookup(x4, q, xdim2, &dx1, tablex, &tp3end); make_lookup(y4, q, ydim2, &dy1, tabley, &tp2end); make_lookup(z4, q, zdim2, &dz1, tablez, &tp1end); tp1 = tablez; dy1 *= xdim2; while(tp1 <= tp1end) { IMAGE_DTYPE *dp2 = &vol[dz1][dy1]; double dat2 = 0.0, *tp2 = tabley; while (tp2 <= tp2end) { register double dat3 = 0.0, *tp3 = tablex; register IMAGE_DTYPE *dp3 = dp2 + dx1; while(tp3 <= tp3end) dat3 += GET(*(dp3++)) * *(tp3++); dat2 += dat3 * *(tp2++); dp2 += xdim2; } dat += (dat2*scale[dz1]+offset[dz1]) * *(tp1++); dz1 ++; } *(image++) = dat; } else *(image++) = background; } } return(0); } /* simple extraction of transverse plane */ void PLANE(p,image,vol,xdim,ydim,scale,offset) int p, xdim,ydim; double image[], scale[],offset[]; IMAGE_DTYPE *vol[]; { int n = xdim*ydim, i; IMAGE_DTYPE *ptr = vol[p-1]; for(i=0; i 1-t && mat[0+0*4] < 1+t && mat[0+1*4] > -t && mat[0+1*4] < t && mat[0+2*4] > -t && mat[0+2*4] < t && mat[0+3*4] > -t && mat[0+3*4] < t && mat[1+0*4] > -t && mat[1+0*4] < t && mat[1+1*4] > 1-t && mat[1+1*4] < 1+t && mat[1+2*4] > -t && mat[1+2*4] < t && mat[1+3*4] > -t && mat[1+3*4] < t && mat[2+0*4] > -t && mat[2+0*4] < t && mat[2+1*4] > -t && mat[2+1*4] < t && mat[2+2*4] > 1-t && mat[2+2*4] < 1+t && fabs(RINT(mat[2+3*4])-mat[2+3*4]) -t && mat[3+0*4] < t && mat[3+1*4] > -t && mat[3+1*4] < t && mat[3+2*4] > -t && mat[3+2*4] < t && mat[3+3*4] > 1-t && mat[3+3*4] < 1+t && xdim1 == xdim2 && ydim1 == ydim2 && mat[2+3*4]>=1.0 && mat[2+3*4]<=zdim2) { int p; p = RINT(mat[2+3*4]); PLANE(p,image,vol,xdim2,ydim2,scale,offset); return(0); } else { if (hold<0) { hold=abs(hold); make_lookup = make_lookup_sinc; return(SLICE_POLY(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, hold+1, background, scale,offset)); } else { make_lookup = make_lookup_poly; } if (hold == 0) return(SLICE_0(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, background, scale,offset)); if (hold == 1) return(SLICE_1(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, background, scale,offset)); else return(SLICE_POLY(mat, image, xdim1, ydim1, vol, xdim2, ydim2, zdim2, hold+1, background, scale,offset)); } } /* Resample image */ void RESAMPLE(m,vol,out,x,y,z,xdim,ydim,zdim, hold, background, scale,offset) int m, xdim,ydim,zdim, hold; double out[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { if (hold<0) { hold=abs(hold); make_lookup = make_lookup_sinc; RESAMPLE_POLY(m,vol,out,x,y,z,xdim,ydim,zdim, hold+1, background, scale,offset); } else { make_lookup = make_lookup_poly; if (hold == 0) RESAMPLE_0(m,vol,out,x,y,z,xdim,ydim,zdim, background, scale,offset); else if (hold == 1) RESAMPLE_1(m,vol,out,x,y,z,xdim,ydim,zdim, background, scale,offset); else RESAMPLE_POLY(m,vol,out,x,y,z,xdim,ydim,zdim, hold+1, background, scale,offset); } } /* Resample image and derivatives */ void RESAMPLE_D(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim, hold, background, scale,offset) int m, xdim,ydim,zdim, hold; double out[],gradx[],grady[],gradz[], x[], y[], z[], background, scale[],offset[]; IMAGE_DTYPE *vol[]; { if (hold<0) { hold=abs(hold); make_lookup_grad = make_lookup_sinc_grad; RESAMPLE_D_POLY(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim, hold+1, background, scale,offset); } else { make_lookup_grad = make_lookup_poly_grad; if (hold == 1) RESAMPLE_D_1(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim, background, scale,offset); else RESAMPLE_D_POLY(m,vol,out,gradx,grady,gradz,x,y,z,xdim,ydim,zdim, hold+1, background, scale,offset); } }