 
 __device__ float4 mul(const matrix4x4& M, const float4& v)
 {
-    float4 r;
-    r.x = dot(v, M.m[0]);
-    r.y = dot(v, M.m[1]);
-    r.z = dot(v, M.m[2]);
-    r.w = dot(v, M.m[3]);
-    return r;
+    float4 result;
+    result.x = M.m[0][0] * v.x + M.m[0][1] * v.y + M.m[0][2] * v.z + M.m[0][3] * v.w;
+    result.y = M.m[1][0] * v.x + M.m[1][1] * v.y + M.m[1][2] * v.z + M.m[1][3] * v.w;
+    result.z = M.m[2][0] * v.x + M.m[2][1] * v.y + M.m[2][2] * v.z + M.m[2][3] * v.w;
+    result.w = M.m[3][0] * v.x + M.m[3][1] * v.y + M.m[3][2] * v.z + M.m[3][3] * v.w;
+
+//    result.x = M.m[0][0] * v.x + M.m[1][0] * v.y + M.m[2][0] * v.z + M.m[3][0] * v.w;
+//    result.y = M.m[0][1] * v.x + M.m[1][1] * v.y + M.m[2][1] * v.z + M.m[3][1] * v.w;
+//    result.z = M.m[0][2] * v.x + M.m[1][2] * v.y + M.m[2][2] * v.z + M.m[3][2] * v.w;
+//    result.w = M.m[0][3] * v.x + M.m[1][3] * v.y + M.m[2][3] * v.z + M.m[3][3] * v.w;
+
+    return result;
 }
 
 __device__ float3 perspectiveDivide(const float4& v)

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

  render<true><<<gridSize, blockSize>>>(d_output, width, height, 2.0f / (float)numSlices, volSize);
  cudaGLUnmapBufferObject(pbo);

  glBindBuffer(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);


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

__device__ float4 phong(const float4& classification, const float3& pos,
                        const float3& L, const float3& H,
                        const float3& Ka, const float3& Kd, const float3& Ks,
                        const float shininess)
{
  float3 N;
  float3 sample1;
  float3 sample2;
  const float DELTA = 0.01f;
  sample1.x = volume(pos - make_float3(DELTA, 0.0f, 0.0f));
  sample2.x = volume(pos + make_float3(DELTA, 0.0f, 0.0f));
  sample1.y = volume(pos - make_float3(0.0f, DELTA, 0.0f));
  sample2.y = volume(pos + make_float3(0.0f, DELTA, 0.0f));
  sample1.z = volume(pos - make_float3(0.0f, 0.0f, DELTA));
  sample2.z = volume(pos + make_float3(0.0f, 0.0f, DELTA));

  N = normalize(sample2 - sample1);
  const float diffuse = fabsf(dot(L, N));
  const float specular = powf(dot(H, N), shininess);

  const float3 c = make_float3(classification);
  float3 tmp = Ka * c + Kd * diffuse * c;
  if (specular > 0.0f)
  {
    tmp += Ks * specular * c;
  }
  return make_float4(tmp.x, tmp.y, tmp.z, classification.w);
}

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

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