        to = intImg->height;
 
    // compute data for determining upper left image corner of each slice and copy it to device
-   vvVector4 start, inc;
-   findSlicePosition(firstSlice, &start, NULL);
-   findSlicePosition(firstSlice+sliceStep, &inc, NULL);
-   inc.sub(&start);
-   vvVector4 cur = start;
+   vvVector4 start, end;
+   findSlicePosition(firstSlice, &start, &end);
+   vvVector4 sinc, einc;
+   findSlicePosition(firstSlice+sliceStep, &sinc, &einc);
+   sinc.sub(&start);
+   einc.sub(&end);
+   vvVector4 scur = start;
+   vvVector4 ecur = end;
+#if defined(VOLTEX3D) && VOLTEX3D==1
+   const int p = principal;
+#else
+   const int p = 2;
+#endif
    for(int slice=firstSlice; slice != lastSlice; slice += sliceStep)
    {
-       h_start[slice].x = cur.e[0] / cur.e[3];
-       h_start[slice].y = cur.e[1] / cur.e[3];
-       cur.add(&inc);
+       if(sliceInterpol)
+       {
+           const float sx = scur.e[0]/scur.e[3];
+           const float sy = scur.e[1]/scur.e[3];
+           const float ex = ecur.e[0]/ecur.e[3];
+           const float ey = ecur.e[1]/ecur.e[3];
+
+           h_start[slice].x = int(floor(sx));
+           h_start[slice].y = int(floor(sy));
+
+           h_stop[slice].x = int(ceil(ex));
+           h_stop[slice].y = int(ceil(ey));
+
+           switch(p)
+           {
+               case 0:
+                   h_tcStep[slice].x = -1.f/(ex-sx);
+                   h_tcStep[slice].y = -1.f/(ey-sy);
+
+                   h_tcStart[slice].x = 1.f + (h_start[slice].x - sx + 0.5f)*h_tcStep[slice].x;
+                   h_tcStart[slice].y = 1.f + (h_start[slice].y - sy + 0.5f)*h_tcStep[slice].y;
+
+                   h_tc3[slice] = 1.f-(slice+0.5f)*1.f/vd->vox[principal];
+                   break;
+                case 1:
+                   h_tcStep[slice].x = -1.f/(ex-sx);
+                   h_tcStep[slice].y = 1.f/(ey-sy);
+
+                   h_tcStart[slice].x = 1.f + (h_start[slice].x - sx + 0.5f)*h_tcStep[slice].x;
+                   h_tcStart[slice].y = (h_start[slice].y - sy + 0.5f)*h_tcStep[slice].y;
+
+                   h_tc3[slice] = (slice+0.5f)*1.f/vd->vox[principal];
+                   break;
+                case 2:
+                   h_tcStep[slice].x = 1.f/(ex-sx);
+                   h_tcStep[slice].y = -1.f/(ey-sy);
+
+                   h_tcStart[slice].x = (h_start[slice].x - sx + 0.5f)*h_tcStep[slice].x;
+                   h_tcStart[slice].y = 1.f + (h_start[slice].y - sy + 0.5f)*h_tcStep[slice].y;
+
+                   h_tc3[slice] = (slice+0.5f)*1.f/vd->vox[principal];
+                   break;
+           }
+
+           ecur.add(&sinc);
+       }
+       else
+       {
+           h_start[slice].x = int(scur.e[0] / scur.e[3] + 0.5f);
+           h_start[slice].y = int(scur.e[1] / scur.e[3] + 0.5f);
+       }
+       scur.add(&sinc);
    }
 
    bool ok = true;
    vvCuda::checkError(&ok, cudaMemcpyToSymbol(c_start, h_start, sizeof(h_start)), "cudaMemcpy start");
+   if(sliceInterpol)
+   {
+       vvCuda::checkError(&ok, cudaMemcpyToSymbol(c_stop, h_stop, sizeof(h_stop)), "cudaMemcpy stop");
+       vvCuda::checkError(&ok, cudaMemcpyToSymbol(c_tcStart, h_tcStart, sizeof(h_tcStart)), "cudaMemcpy tcStart");
+       vvCuda::checkError(&ok, cudaMemcpyToSymbol(c_tcStep, h_tcStep, sizeof(h_tcStep)), "cudaMemcpy tcStep");
+       vvCuda::checkError(&ok, cudaMemcpyToSymbol(c_tc3, h_tc3, sizeof(h_tc3)), "cudaMemcpy tc3");
+   }
 
 #ifdef VOLTEX3D
-   tex_raw.normalized = false;
-   tex_raw.filterMode = cudaFilterModePoint;
+   tex_raw.normalized = sliceInterpol;
+   tex_raw.filterMode = sliceInterpol ? cudaFilterModeLinear : cudaFilterModePoint;
    tex_raw.addressMode[0] = cudaAddressModeClamp;
    tex_raw.addressMode[1] = cudaAddressModeClamp;
    tex_raw.addressMode[2] = cudaAddressModeClamp;

    return NULL;
}

CompositionFunction selectCompositionWithPrincipalAndSliceStep(int principal, int sliceStep)
{
    switch(principal)
    {
        case 0:
            return selectCompositionWithSliceStep<0>(sliceStep);
        case 1:
            return selectCompositionWithSliceStep<1>(sliceStep);
        case 2:
            return selectCompositionWithSliceStep<2>(sliceStep);
        default:
            assert("principal axis out of range" == NULL);


#include "vvcuda.h"
#include "vvcudapar.h"

const int MAX_SLICES = 1600;

__constant__ int   c_vox[5];
__constant__ int2 c_start[MAX_SLICES];
static int2 h_start[MAX_SLICES];
__constant__ int2 c_stop[MAX_SLICES];
static int2 h_stop[MAX_SLICES];
__constant__ float2 c_tcStart[MAX_SLICES];
static float2 h_tcStart[MAX_SLICES];
__constant__ float2 c_tcStep[MAX_SLICES];
static float2 h_tcStep[MAX_SLICES];
__constant__ float c_tc3[MAX_SLICES];
static float h_tc3[MAX_SLICES];


const int MaxCompositeSlices = MAX_SLICES;

//#define SHMCLASS
//#define NOOP
