   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

  cudaBindTextureToArray(tex_zdvf, (cudaArray*) DVFcomponentArrays[2], channelDesc);
  CUDA_CHECK_ERROR;

  ///////////////////////////////////////
  // Copy matrices, bind them to textures

  float *dev_IndexInputToPPInput;
  cudaMalloc( (void**)&dev_IndexInputToPPInput, 12*sizeof(float) );
  cudaMemcpy (dev_IndexInputToPPInput, IndexInputToPPInputMatrix, 12*sizeof(float), cudaMemcpyHostToDevice);
  cudaBindTexture (0, tex_IndexInputToPPInputMatrix, dev_IndexInputToPPInput, 12*sizeof(float) );

  float *dev_IndexInputToIndexDVF;
  cudaMalloc( (void**)&dev_IndexInputToIndexDVF, 12*sizeof(float) );
  cudaMemcpy (dev_IndexInputToIndexDVF, IndexInputToIndexDVFMatrix, 12*sizeof(float), cudaMemcpyHostToDevice);
  cudaBindTexture (0, tex_IndexInputToIndexDVFMatrix, dev_IndexInputToIndexDVF, 12*sizeof(float) );

  float *dev_PPInputToIndexInput;
  cudaMalloc( (void**)&dev_PPInputToIndexInput, 12*sizeof(float) );
  cudaMemcpy (dev_PPInputToIndexInput, PPInputToIndexInputMatrix, 12*sizeof(float), cudaMemcpyHostToDevice);
  cudaBindTexture (0, tex_PPInputToIndexInputMatrix, dev_PPInputToIndexInput, 12*sizeof(float) );

  // The optimized version runs when only one of the axis of the detector is
  // parallel to the y axis of the volume

//  if(fabs(matrix[5])<1e-10 && fabs(matrix[9])<1e-10)
//    {
//    // Thread Block Dimensions
//    const int tBlock_x = 32;
//    const int tBlock_y = 16;

//    // Each segment gets 1 thread
//    unsigned int  blocksInX = (vol_size[0]-1)/tBlock_x + 1;
//    unsigned int  blocksInY = (vol_size[2]-1)/tBlock_y + 1;
//    dim3 dimGrid  = dim3(blocksInX, blocksInY);
//    dim3 dimBlock = dim3(tBlock_x, tBlock_y, 1);

//    // Note: cbi->img AND cbi->matrix are passed via texture memory
//    //-------------------------------------
//    kernel_optim <<< dimGrid, dimBlock >>> ( dev_vol_in,
//                                             dev_vol_out,
//                                             make_int3(vol_size[0], vol_size[1], vol_size[2]) );
//    }
//  else
//    {
//    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;

//    if(CUDA_VERSION<4000 || GetCudaComputeCapability(device).first<=1)
//      {
//      dim3 dimGrid  = dim3(blocksInX, blocksInY*blocksInZ);
//      dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);


//      // Note: cbi->img AND cbi->matrix are passed via texture memory
//      //-------------------------------------
//      kernel <<< dimGrid, dimBlock >>> ( dev_vol_in,
//                                         dev_vol_out,
//                                         make_int3(vol_size[0], vol_size[1], vol_size[2]),
//                                         blocksInY );
//      }
//    else
//      {
      dim3 dimGrid  = dim3(blocksInX, blocksInY, blocksInZ);
      dim3 dimBlock = dim3(tBlock_x, tBlock_y, tBlock_z);


      // Note: cbi->img AND cbi->matrix are passed via texture memory
      //-------------------------------------
      kernel_warp_back_project_3Dgrid <<< dimGrid, dimBlock >>> ( dev_vol_in,
                                                dev_vol_out,
                                                make_int3(vol_size[0], vol_size[1], vol_size[2]));
//      }

//    }
  CUDA_CHECK_ERROR;

  // Unbind the image and projection matrix textures
  cudaUnbindTexture (tex_xdvf);
  cudaUnbindTexture (tex_ydvf);
  cudaUnbindTexture (tex_zdvf);
  CUDA_CHECK_ERROR;
  cudaUnbindTexture (tex_img);
  CUDA_CHECK_ERROR;
  cudaUnbindTexture (tex_matrix);
  CUDA_CHECK_ERROR;
  cudaUnbindTexture (tex_IndexInputToPPInputMatrix);
  cudaUnbindTexture (tex_IndexInputToIndexDVFMatrix);
  cudaUnbindTexture (tex_PPInputToIndexInputMatrix);
  CUDA_CHECK_ERROR;

  // Cleanup
  cudaFreeArray ((cudaArray*) DVFcomponentArrays[0]);
  cudaFreeArray ((cudaArray*) DVFcomponentArrays[1]);
  cudaFreeArray ((cudaArray*) DVFcomponentArrays[2]);
  delete[] DVFcomponentArrays;
  CUDA_CHECK_ERROR;
  cudaFreeArray ((cudaArray*)array_img);
  CUDA_CHECK_ERROR;
  cudaFree (dev_matrix);
  CUDA_CHECK_ERROR;
  cudaFree (dev_IndexInputToPPInput);
  cudaFree (dev_IndexInputToIndexDVF);
  cudaFree (dev_PPInputToIndexInput);
  CUDA_CHECK_ERROR;

  // Destroy CUBLAS context

  cudaMemcpyToArray( array_img, 0, 0, dev_img, img_size[0] * img_size[1] * sizeof(float), cudaMemcpyDeviceToDevice);
  cudaBindTextureToArray( tex_img, (cudaArray*)array_img, channelDesc);

  // Extent stuff, will be used for each component extraction
  cudaExtent dvfExtent = make_cudaExtent(dvf_size[0], dvf_size[1], dvf_size[2]);

