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

  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

  PPDisplaced.z = PPinInput.z + Displacement.z;

  float3 IndexInInput;
  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 =  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;
