 
 typedef struct
 {
-  float4 m[4];
+  float m[4][4];
 } matrix4x4;
 
-__constant__ matrix4x4 c_invViewMatrix;  // inverse view matrix
+__constant__ matrix4x4 c_invViewMatrix;
 
 struct Ray
 {
   float3 o;
   float3 d;
 };
+
+__device__ float volume(const float x, const float y, const float z)
+{
+  return tex3D(volTexture, x, y, z);
+}
+
+__device__ float volume(const float3& pos)
+{
+  return tex3D(volTexture, pos.x, pos.y, pos.z);
+}
 
 __device__ bool intersectBox(const Ray& r, const float3& boxmin, const float3& boxmax,
                              float *tnear, float *tfar)


__device__ uint rgbaFloatToInt(float3 rgb)
{
    float4 rgba = make_float4(rgb.x, rgb.y, rgb.z, 1.0f);
    return rgbaFloatToInt(rgba);
}

template<bool frontToBack>
__global__ void render(uint *d_output, const uint width, const uint height, const float dist,
                       const float3 volSize)
{
  const int maxSteps = 10000;
  const float tstep = dist;
  const float opacityThreshold = 0.95f;
  const float3 boxMin = make_float3(-volSize.x, -volSize.y, -volSize.z);
  const float3 boxMax = make_float3(volSize.x, volSize.y, volSize.z);

  const uint x = blockIdx.x * blockDim.x + threadIdx.x;
  const uint y = blockIdx.y * blockDim.y + threadIdx.y;

    volSize = make_float3(vd->vox[0] * inv, vd->vox[1] * inv, 1.0f);
  }

  bool isOrtho = pr.isProjOrtho();

  vvVector3 eye;
  getEyePosition(&eye);
  eye.multiply(&invMV);

  vvVector3 origin;

  vvVector3 normal;
  getObjNormal(normal, origin, eye, invMV, isOrtho);

  const float3 N = make_float3(normal[0], normal[1], normal[2]);

  const float3 L(-N);

  // Viewing direction.
  const float3 V(-N);

  // Half way vector.
  const float3 H = normalize(L + V);
  render<true, true><<<gridSize, blockSize>>>(d_output, width, height, 2.0f / (float)numSlices, volSize * 0.5f, L, H);
  cudaGLUnmapBufferObject(pbo);

  glBindBuffer(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
