        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;

//#define SHMLOAD
#define VOLTEX3D 1 // undef, 1 or 3
//#define PITCHED
#define FLOATDATA
//#define CONSTLOAD
//#define THREADPERVOXEL
//#define CONSTDATA
//#define FLOATIMG
#define EARLYRAYTERM

const int Repetitions = 1;

            switch(principal)
            {
                case 0:
                    v = tex3D(tex_raw, c_vox[2]-slice-1, c_vox[0]-vidx-1, c_vox[1]+iPosY-line-1);
                    break;
                case 1:
                    v = tex3D(tex_raw, line-iPosY, slice, c_vox[1]-vidx-1);
                    break;
                case 2:
                    v = tex3D(tex_raw, vidx, c_vox[0]+iPosY-line-1, slice);
