   float4 dst = make_float4(0.0f);
   float t = tnear;
   float3 pos = eyeRay.o + eyeRay.d * tnear;
-  float3 step = eyeRay.d * tstep;
+  const float3 step = eyeRay.d * tstep;
 
   for (int i=0; i<maxSteps; ++i)
   {
-    const float sample = tex3D(volTexture,
-                               (pos.x + volSize.x) / (volSize.x * 2),
-                               (pos.y + volSize.y) / (volSize.y * 2),
-                               (pos.z + volSize.z) / (volSize.z * 2));
+    const float3 texCoord = make_float3((pos.x + volSizeHalf.x) / (volSizeHalf.x * 2),
+                                        (pos.y + volSizeHalf.y) / (volSizeHalf.y * 2),
+                                        (pos.z + volSizeHalf.z) / (volSizeHalf.z * 2));
+    const float sample = volume(texCoord);
 
     // Post-classification transfer-function lookup.
     float4 src = tex1D(tfTexture, sample);
+
+    // Local illumination.
+    if (lighting && (src.w > 0.1))
+    {
+      const float3 Ka = make_float3(0.0f, 0.0f, 0.0f);
+      const float3 Kd = make_float3(0.8f, 0.8f, 0.8f);
+      const float3 Ks = make_float3(0.8f, 0.8f, 0.8f);
+      const float shininess = 1000.0f;
+      src = phong(src, texCoord, L, H, Ka, Kd, Ks, shininess);
+    }
 
     // "under" operator for back-to-front blending
     //dst = lerp(dst, src, src.w);

  // Inverse modelview-projection matrix.
  vvMatrix mvp, pr;
  getModelviewMatrix(&mvp);
  getProjectionMatrix(&pr);
  mvp.multiplyPost(&pr);
  mvp.invert();


__device__ float4 mul(const matrix4x4& M, const float4& v)
{
    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)
