             d.z += w*c.z;
             d.w -= w;
         }
+        sf = sb;
+    }
+};
+
+__device__ bool outsideBounds(int2 p, int2 from1, int2 to1, int2 from2, int2 to2)
+{
+    return (p.x < from1.x && p.x < from2.x)
+        || (p.y < from1.y && p.y < from2.y)
+        || (p.x >= to1.x && p.x >= to2.x)
+        || (p.y >= to1.y && p.y >= to2.y);
+}
+
+__device__ bool fullyInsideIntersection(int2 p1, int2 p2, int2 from1, int2 to1, int2 from2, int2 to2)
+{
+    return p1.x >= from1.x && p1.x >= from2.x
+        && p1.y >= from1.y && p1.y >= from2.y
+        && p2.x < to1.x && p2.x < to2.x
+        && p2.y < to1.y && p2.y < to2.y;
+}
+
+#ifdef VOLTEX3D
+template<typename Scalar, int BPV, typename Pixel, int sliceStep, int principal, bool earlyRayTerm, bool preInt>
+__global__ void compositeSlicesBilinear(
+      uchar4 * __restrict__ img, int width, int height,
+#ifdef PITCHED
+      const cudaPitchedPtr pvoxels,
+#else
+      const Scalar * __restrict__ voxels,
+#endif
+      int firstSlice, int lastSlice,
+      int2 from, int2 to)
+{
+#ifdef PATCHES
+    Ray<Pixel, principal, preInt> ray;
+    int2 p(coord(from));
+
+    // composite slices for this image line
+    for (int sliceb=firstSlice; sliceStep>0 ? sliceb<lastSlice : sliceb>lastSlice; sliceb += sliceStep*SliceStack)
+    {
+        if(earlyRayTerm && isOpaque(ray.d))
+            break;
+
+        int last = lastSlice;
+        if(sliceStep > 0 && last>sliceb+SliceStack)
+            last = sliceb+SliceStack;
+        if(sliceStep < 0 && last<sliceb-SliceStack)
+            last = sliceb-SliceStack;
+
+        if(outsideBounds(p, c_start[sliceb], c_stop[sliceb],
+                    c_start[last-sliceStep], c_stop[last-sliceStep]))
+            continue;
+
+#ifdef UNROLL
+        const int2 p1 = make_int2(blockDim.x*blockIdx.x+from.x,
+                blockDim.y*blockIdx.y+from.y);
+        const int2 p2 = make_int2(blockDim.x*blockIdx.x+from.x+blockDim.x-1,
+                blockDim.y*blockIdx.y+from.y+blockDim.y-1);
+        if(fullyInsideIntersection(p1, p2,
+                    c_start[sliceb], c_stop[sliceb],
+                    c_start[last-sliceStep], c_stop[last-sliceStep]))
+        {
+#pragma unroll 4
+            for (int slice=sliceb; sliceStep>0 ? slice<last : slice>last; slice += sliceStep)
+            {
+
+                // compute upper left image corner
+                const int iPosY = c_start[slice].y;
+                const int iPosX = c_start[slice].x;
+                const int vidx = p.x - iPosX;
+                ray.accumulate(vidx, p.y-iPosY, slice);
+            }
+        }
+        else
+#endif
+        {
+            for (int slice=sliceb; sliceStep>0 ? slice<last : slice>last; slice += sliceStep)
+            {
+                if(earlyRayTerm && isOpaque(ray.d))
+                    break;
+
+                // compute upper left image corner
+                const int iPosY = c_start[slice].y;
+
+                const int iPosX = c_start[slice].x;
+                const int endX = c_stop[slice].x;
+
+                if(p.y < iPosY)
+                    continue;
+                if(p.y >= c_stop[slice].y)
+                    continue;
+
+                // Traverse intermediate image pixels which correspond to the current slice.
+                // 1 is subtracted from each loop counter to remain inside of the volume boundaries:
+                if(p.x<iPosX)
+                    continue;
+                if(p.x>=endX)
+                    continue;
+
+                const int vidx = p.x - iPosX;
+                ray.accumulate(vidx, p.y-iPosY, slice);
+            }
+        }
     }
 
     // copy pixel to intermediate image
     if(p.x >= from.x && p.x < to.x && p.y >= from.y && p.y < to.y)
     {
         uchar4 *dest = img + p.y*width+p.x;
-        blend(dest, d);
+        blend(dest, ray.d);
     }
 #else
     // this block's line from the intermediate image

#include "vvcudasw.h"

const int MAX_SLICES = 1600;

__constant__ int   c_vox[5];
__constant__ int2 c_start[MAX_SLICES];

#ifdef NOSHMEM
    for (int ix=threadIdx.x+from.x; ix<to.x; ix+=blockDim.x)
    {
        Ray<Pixel, principal, preInt> ray;

        // composite slices for this image line
        for (int slice=firstSlice; slice!=lastSlice; slice += sliceStep)
        {
            if(earlyRayTerm && isOpaque(ray.d))
                break;

            // compute upper left image corner
