   PPDisplaced.z = PPinInput.z + Displacement.z;
 
   float3 IndexInInput;
-  IndexInInput.x =  tex1Dfetch(tex_PPInputToIndexInputMatrix, 0) * PPDisplaced.x
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 1) * PPDisplaced.y
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 2) * PPDisplaced.z
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 3);
-  IndexInInput.y =  tex1Dfetch(tex_PPInputToIndexInputMatrix, 4) * PPDisplaced.x
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 5) * PPDisplaced.y
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 6) * PPDisplaced.z
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 7);
-  IndexInInput.z =  tex1Dfetch(tex_PPInputToIndexInputMatrix, 8) * PPDisplaced.x
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 9) * PPDisplaced.y
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 10)* PPDisplaced.z
-                  + tex1Dfetch(tex_PPInputToIndexInputMatrix, 11);
+  IndexInInput.x =  c_PPInputToIndexInputMatrix[0] * PPDisplaced.x
+                  + c_PPInputToIndexInputMatrix[1] * PPDisplaced.y
+                  + c_PPInputToIndexInputMatrix[2] * PPDisplaced.z
+                  + c_PPInputToIndexInputMatrix[3];
+  IndexInInput.y =  c_PPInputToIndexInputMatrix[4] * PPDisplaced.x
+                  + c_PPInputToIndexInputMatrix[5] * PPDisplaced.y
+                  + c_PPInputToIndexInputMatrix[6] * PPDisplaced.z
+                  + c_PPInputToIndexInputMatrix[7];
+  IndexInInput.z =  c_PPInputToIndexInputMatrix[8] * PPDisplaced.x
+                  + c_PPInputToIndexInputMatrix[9] * PPDisplaced.y
+                  + c_PPInputToIndexInputMatrix[10]* PPDisplaced.z
+                  + c_PPInputToIndexInputMatrix[11];
 
   float3 ip;
   float  voxel_data;
 
   // matrix multiply
-  ip.x =  tex1Dfetch(tex_matrix, 0) * IndexInInput.x
-        + tex1Dfetch(tex_matrix, 1) * IndexInInput.y
-        + tex1Dfetch(tex_matrix, 2) * IndexInInput.z
-        + tex1Dfetch(tex_matrix, 3);
-  ip.y =  tex1Dfetch(tex_matrix, 4) * IndexInInput.x
-        + tex1Dfetch(tex_matrix, 5) * IndexInInput.y
-        + tex1Dfetch(tex_matrix, 6) * IndexInInput.z
-        + tex1Dfetch(tex_matrix, 7);
-  ip.z =  tex1Dfetch(tex_matrix, 8) * IndexInInput.x
-        + tex1Dfetch(tex_matrix, 9) * IndexInInput.y
-        + tex1Dfetch(tex_matrix, 10) * IndexInInput.z
-        + tex1Dfetch(tex_matrix, 11);
+  ip.x =  c_matrix[0] * IndexInInput.x
+        + c_matrix[1] * IndexInInput.y
+        + c_matrix[2] * IndexInInput.z
+        + c_matrix[3];
+  ip.y =  c_matrix[4] * IndexInInput.x
+        + c_matrix[5] * IndexInInput.y
+        + c_matrix[6] * IndexInInput.z
+        + c_matrix[7];
+  ip.z =  c_matrix[8] * IndexInInput.x
+        + c_matrix[9] * IndexInInput.y
+        + c_matrix[10] * IndexInInput.z
+        + c_matrix[11];
 
   // Change coordinate systems
   ip.z = 1 / ip.z;

  // 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, int3 vol_size)
//{
//  unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
//  unsigned int j = 0;
//  unsigned int k = blockIdx.y * blockDim.y + threadIdx.y;

//  if (i >= vol_size.x || k >= vol_size.z)
//    {
//    return;
//    }

//  // Index row major into the volume
//  long int vol_idx = i + k*vol_size.y*vol_size.x;

//  float3 ip;

//  // matrix multiply
//  ip.x = tex1Dfetch(tex_matrix, 0)*i + tex1Dfetch(tex_matrix, 2)*k + tex1Dfetch(tex_matrix, 3);
//  ip.y = tex1Dfetch(tex_matrix, 4)*i + tex1Dfetch(tex_matrix, 6)*k + tex1Dfetch(tex_matrix, 7);
//  ip.z = tex1Dfetch(tex_matrix, 8)*i + tex1Dfetch(tex_matrix, 10)*k + tex1Dfetch(tex_matrix, 11);

//  // Change coordinate systems
//  ip.z = 1 / ip.z;
//  ip.x = ip.x * ip.z;
//  ip.y = ip.y * ip.z;
//  float dx = tex1Dfetch(tex_matrix, 1)*ip.z;
//  float dy = tex1Dfetch(tex_matrix, 5)*ip.z;

//  // Place it into the volume segment
//  for(; j<vol_size.y; j++)
//    {
//    dev_vol_out[vol_idx] = dev_vol_in[vol_idx] + tex2D(tex_img, ip.x, ip.y);
//    vol_idx+=vol_size.x;
//    ip.x+=dx;
//    ip.y+=dy;
//    }
//}

//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
// K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_


// T E X T U R E S ////////////////////////////////////////////////////////
texture<float, 2, cudaReadModeElementType> tex_img;
texture<float, 3, cudaReadModeElementType> tex_xdvf;
texture<float, 3, cudaReadModeElementType> tex_ydvf;
texture<float, 3, cudaReadModeElementType> tex_zdvf;
///////////////////////////////////////////////////////////////////////////

// CONSTANTS //////////////////////////////////////////////////////////////
__constant__ float c_matrix[12];
__constant__ float c_IndexInputToIndexDVFMatrix[12];
__constant__ float c_PPInputToIndexInputMatrix[12];
__constant__ float c_IndexInputToPPInputMatrix[12];
////////////////////////////////////////////////////////////////////////////

//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
// K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_( S T A R T )_
//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_

__global__
void kernel_warp_back_project_3Dgrid(float *dev_vol_in,
                                     float * dev_vol_out,
                                     int3 vol_size)
{
  unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
  unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;
