 
   // Place it into the volume
   dev_vol_out[vol_idx] = dev_vol_in[vol_idx] + voxel_data;
+}
+
+
+__global__
+void kernel_optim(float *dev_vol_in, float *dev_vol_out)
+{
+  // Allocate a few shared buffers
+  extern __shared__ float shared[];
+  float* accumulators = &shared[0];
+
+  // Compute the starting position of the shared memory available to the current thread
+  unsigned int startIndex = (threadIdx.x + (threadIdx.y * ( blockDim.x + threadIdx.z * blockDim.y) ) ) * 8;
+
+  unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
+  unsigned int j = (blockIdx.y * blockDim.y + threadIdx.y)*8;
+  unsigned int k = blockIdx.z * blockDim.z + threadIdx.z;
+
+  if (i >= c_vol_size.x || j >= c_vol_size.y || k >= c_vol_size.z)
+    {
+    return;
+    }
+  unsigned int range_j = max(c_vol_size.y - j, 8);
+
+  // Index row major into the volume
+  long int vol_idx = i + (j + k*c_vol_size.y)*(c_vol_size.x);
+
+  float3 ip;
+
+  // Initialize the shared memory with zeros
+  for(unsigned int l=0; l<range_j; l++)
+    accumulators[startIndex + l]=0;
+
+  // Process all the projections in the stack
+  for (unsigned int proj = 0; proj<c_proj_size.z; proj++)
+    {
+    // matrix multiply
+    ip = matrix_multiply(make_float3(i,j,k), &(c_matrices[12*proj]));
+
+    // Change coordinate systems
+    ip.z = 1 / ip.z;
+    ip.x = ip.x * ip.z;
+    ip.y = ip.y * ip.z;
+    float dx = c_matrices[12*proj +  1]*ip.z;
+    float dy = c_matrices[12*proj +  5]*ip.z;
+
+    // Walk voxels and accumulate their back projection in shared memory
+    for(unsigned int l=0; l<range_j; l++)
+      {
+      accumulators[startIndex + l] += tex2DLayered(tex_proj, ip.x, ip.y, proj);
+      ip.x+=dx;
+      ip.y+=dy;
+      }
+    }
+
+  // Write the accumulated results into global memory
+  for(unsigned int l=0; l<range_j; l++)
+    {
+    dev_vol_out[vol_idx] = dev_vol_in[vol_idx] + accumulators[startIndex + l];
+    vol_idx+=c_vol_size.x;
+    }
 }
 
 //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_

  float3 ip;
  float  voxel_data = 0;

  for (unsigned int proj = 0; proj<c_projSize.z; proj++)
    {
    // matrix multiply
    ip = matrix_multiply(make_float3(i,j,k), &(c_matrices[12*proj]));


// Constant memory
__constant__ float c_matrices[1024 * 12]; //Can process stacks of at most 1024 projections
__constant__ int3 c_proj_size;
__constant__ int3 c_vol_size;

//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
