                              float *dev_proj,
                              float *dev_matrix)
 {
+  // Unbind the volume and matrix textures
+  cudaUnbindTexture (tex_vol);
+  CUDA_CHECK_ERROR;
+  cudaUnbindTexture (tex_matrix);
+  CUDA_CHECK_ERROR;
   // Deallocate memory on the device
   cudaFree (dev_proj);
   CUDA_CHECK_ERROR;

*****************/
#include <cuda.h>
//#include "/home/mvila/NVIDIA_GPU_Computing_SDK/C/common/inc/cutil_math.h"
inline __host__ __device__ float3 operator-(float3 a, float3 b)
{
    return make_float3(a.x - b.x, a.y - b.y, a.z - b.z);
}
inline __host__ __device__ float3 fminf(float3 a, float3 b)
{
	return make_float3(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z));
}
inline __host__ __device__ float3 fmaxf(float3 a, float3 b)
{
	return make_float3(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z));
}
inline __host__ __device__ float dot(float3 a, float3 b)
{ 
    return a.x * b.x + a.y * b.y + a.z * b.z;
}
inline __host__ __device__ float3 operator/(float3 a, float3 b)
{
    return make_float3(a.x / b.x, a.y / b.y, a.z / b.z);
}
inline __host__ __device__ float3 operator/(float3 a, float b)
{
    return make_float3(a.x / b, a.y / b, a.z / b);
}
inline __host__ __device__ float3 operator*(float3 a, float3 b)
{
    return make_float3(a.x * b.x, a.y * b.y, a.z * b.z);
}
inline __host__ __device__ float3 operator*(float b, float3 a)
{
    return make_float3(b * a.x, b * a.y, b * a.z);
}
inline __host__ __device__ float3 operator+(float3 a, float3 b)
{
    return make_float3(a.x + b.x, a.y + b.y, a.z + b.z);
}
inline __host__ __device__ float3 operator+(float3 a, float b)
{
    return make_float3(a.x + b, a.y + b, a.z + b);
}
inline __host__ __device__ void operator+=(float3 &a, float3 b)
{
    a.x += b.x; a.y += b.y; a.z += b.z;
}

// T E X T U R E S ////////////////////////////////////////////////////////
texture<float, cudaTextureType3D, cudaReadModeElementType> tex_vol;

  CUDA_CHECK_ERROR;
  cudaFreeArray ((cudaArray*)dev_vol);
  CUDA_CHECK_ERROR;

}
