 *****************/
 #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;
-}
+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;
  // Unbind the volume and matrix textures
  cudaUnbindTexture (tex_vol);
  CUDA_CHECK_ERROR;
  cudaUnbindTexture (tex_matrix);
  CUDA_CHECK_ERROR;
}

// FUNCTION: CUDA_forward_project() //////////////////////////////////
void
CUDA_forward_project(int blockSize[3],
                     float *host_proj,
                     float *dev_proj,
                     double source_position[3],
                     int projections_size[2],
