   int device;
   cudaGetDevice(&device);
 
-  // Thread Block Dimensions
-  const int tBlock_x = 16;
-  const int tBlock_y = 4;
-  const int tBlock_z = 4;
-
-  // Each element in the volume (each voxel) gets 1 thread
-  unsigned int  blocksInX = (vol_size[0]-1)/tBlock_x + 1;
-  unsigned int  blocksInY = (vol_size[1]-1)/tBlock_y + 1;
-  unsigned int  blocksInZ = (vol_size[2]-1)/tBlock_z + 1;
+  // Choose which kernel to use. The optimized version runs only when
+  // one of the axes of the detector is parallel to the y axis of the volume
+  bool runKernelOptim = true;
+  for (unsigned int proj = 0; proj<proj_size[2]; proj++)
+    {
+    if(fabs(matrices[12*proj + 5])>1e-10 || fabs(matrices[12*proj + 9])>1e-10)
+      runKernelOptim = false;
+    }
 
   // Run kernels. Note: Projection data is passed via texture memory,
   // transform matrix is passed via constant memory
-  if(CUDA_VERSION<4000 || GetCudaComputeCapability(device).first<=1)
-    {
-    dim3 dimGrid  = dim3(blocksInX, blocksInY*blocksInZ);
-    dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);
-
-    kernel <<< dimGrid, dimBlock >>> ( dev_vol_in,
-                                       dev_vol_out,
-                                       blocksInY );
-    }
-  else
-    {
+  if (runKernelOptim)
+    {
+    cudaDeviceProp properties;
+    cudaGetDeviceProperties(&properties, device);
+
+    // Thread Block Dimensions
+    size_t sharedMemPerThread = sizeof(float) * 8;
+    const int tBlock_x = 16;
+    const int tBlock_y = 4;
+    int maximumZSizeOfBlocks = properties.sharedMemPerBlock / ( sharedMemPerThread * tBlock_x * tBlock_y );
+    int tBlock_z = min(maximumZSizeOfBlocks, 4);
+
+    // Each segment gets 1 thread
+    unsigned int  blocksInX = (vol_size[0]-1)/tBlock_x + 1;
+    unsigned int  blocksInY = (vol_size[1]-1)/(tBlock_y * 8) + 1;
+    unsigned int  blocksInZ = (vol_size[2]-1)/tBlock_z + 1;
     dim3 dimGrid  = dim3(blocksInX, blocksInY, blocksInZ);
     dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);
     CUDA_CHECK_ERROR;
 
-    kernel_3Dgrid <<< dimGrid, dimBlock >>> ( dev_vol_in,
-                                              dev_vol_out);
+    std::cout << "Using blocks of size " << tBlock_x << "x" << tBlock_y << "x" << tBlock_z << std::endl;
+    std::cout << "Using grid of size " << blocksInX << "x" << blocksInY << "x" << blocksInZ << std::endl;
+
+    // Note: cbi->img AND cbi->matrix are passed via texture memory
+    //-------------------------------------
+    kernel_optim <<< dimGrid, dimBlock, sharedMemPerThread * tBlock_x * tBlock_y * tBlock_z >>> ( dev_vol_in,
+                                                                                                  dev_vol_out);
+    }
+  else
+    {
+    // Thread Block Dimensions
+    const int tBlock_x = 16;
+    const int tBlock_y = 4;
+    const int tBlock_z = 4;
+
+    // Each element in the volume (each voxel) gets 1 thread
+    unsigned int  blocksInX = (vol_size[0]-1)/tBlock_x + 1;
+    unsigned int  blocksInY = (vol_size[1]-1)/tBlock_y + 1;
+    unsigned int  blocksInZ = (vol_size[2]-1)/tBlock_z + 1;
+
+    if(CUDA_VERSION<4000 || GetCudaComputeCapability(device).first<=1)
+      {
+      dim3 dimGrid  = dim3(blocksInX, blocksInY*blocksInZ);
+      dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);
+
+      kernel <<< dimGrid, dimBlock >>> ( dev_vol_in,
+                                         dev_vol_out,
+                                         blocksInY );
+      }
+    else
+      {
+      dim3 dimGrid  = dim3(blocksInX, blocksInY, blocksInZ);
+      dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);
+      CUDA_CHECK_ERROR;
+
+      kernel_3Dgrid <<< dimGrid, dimBlock >>> ( dev_vol_in,
+                                                dev_vol_out);
+      }
     }
 
   // Unbind the image and projection matrix textures


  // Place it into the volume
  dev_vol_out[vol_idx] = dev_vol_in[vol_idx] + voxel_data;
}

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

  float3 ip;
  float  voxel_data = 0;

  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]));
