   long int vol_idx = i + (j + k*vol_size.y)*(vol_size.x);
 
   float3 IndexInDVF;
-  IndexInDVF.x =  tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 0) * i
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 1) * j
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 2) * k
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 3);
-  IndexInDVF.y =  tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 4) * i
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 5) * j
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 6) * k
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 7);
-  IndexInDVF.z =  tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 8) * i
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 9) * j
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 10) * k
-                + tex1Dfetch(tex_IndexInputToIndexDVFMatrix, 11);
+  IndexInDVF.x =  c_IndexInputToIndexDVFMatrix[0] * i
+                + c_IndexInputToIndexDVFMatrix[1] * j
+                + c_IndexInputToIndexDVFMatrix[2] * k
+                + c_IndexInputToIndexDVFMatrix[3];
+  IndexInDVF.y =  c_IndexInputToIndexDVFMatrix[4] * i
+                + c_IndexInputToIndexDVFMatrix[5] * j
+                + c_IndexInputToIndexDVFMatrix[6] * k
+                + c_IndexInputToIndexDVFMatrix[7];
+  IndexInDVF.z =  c_IndexInputToIndexDVFMatrix[8] * i
+                + c_IndexInputToIndexDVFMatrix[9] * j
+                + c_IndexInputToIndexDVFMatrix[10] * k
+                + c_IndexInputToIndexDVFMatrix[11];
 
   // Get each component of the displacement vector by
   // interpolation in the dvf

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

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

  // Change coordinate systems
  ip.z = 1 / ip.z;


// 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;
