     tnear = 0.0f;
   }
 
+  // Calc hits with clip sphere.
+  float tsnear;
+  float tsfar;
+  if (clipSphere)
+  {
+    intersectSphere(ray, sphereCenter, sphereRadius, &tsnear, &tsfar);
+  }
+
+  // Calc hits with clip plane.
+  float tpnear;
+  float nddot;
+  if (clipPlane)
+  {
+    intersectPlane(ray, planeNormal, planeDist, &nddot, &tpnear);
+  }
+
   float4 dst = make_float4(0.0f);
   float t = tnear;
-  float3 pos = eyeRay.o + eyeRay.d * tnear;
-  const float3 step = eyeRay.d * tstep;
+  float3 pos = ray.o + ray.d * tnear;
+  const float3 step = ray.d * tstep;
+
+  float maxIntensity = 0.0f;
+  float minIntensity = FLT_MAX;
 
   for (int i=0; i<maxSteps; ++i)
   {
-    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));
+    // Test for clipping.
+    if ((
+        (clipSphere && (t >= tsnear) && (t <= tsfar)) // Sphere.
+     || (clipPlane && (((t <= tpnear) && (nddot >= 0.0f))
+         || ((t >= tpnear) && (nddot < 0.0f)))) // Plane.
+    ))
+    {
+      t += tstep;
+      if (t > tfar)
+      {
+        break;
+      }
+      pos += step;
+      continue;
+    }
+
+    const float3 texCoord = calcTexCoord(pos, volSizeHalf);
     const float sample = volume(texCoord);
 
     // Post-classification transfer-function lookup.
-    float4 src = tex1D(tfTexture, sample);
+    float4 src;
+
+    if (mipMode == 0)
+    {
+      src = tex1D(tfTexture, sample);
+    }
+    else if (mipMode == 1 && (sample > maxIntensity))
+    {
+      dst = tex1D(tfTexture, sample);
+      maxIntensity = sample;
+    }
+    else if (mipMode == 2 && (sample  < minIntensity))
+    {
+      dst = tex1D(tfTexture, sample);
+      minIntensity = sample;
+    }
 
     // Local illumination.
     if (lighting && (src.w > 0.1))

    src.y *= src.w;
    src.z *= src.w;

    if (frontToBack)
    {
      dst = dst + src * (1.0f - dst.w);
    }

  return tex3D(volTexture, pos.x, pos.y, pos.z);
}

__device__ float3 calcTexCoord(const float3& pos, const float3& volSizeHalf)
{
  return make_float3((pos.x + volSizeHalf.x) / (volSizeHalf.x * 2.0f),
                     (pos.y + volSizeHalf.y) / (volSizeHalf.y * 2.0f),
                     (pos.z + volSizeHalf.z) / (volSizeHalf.z * 2.0f));
}

__device__ void solveQuadraticEquation(const float A, const float B, const float C,
                                       float* tnear, float* tfar)
{
  const float discrim = B * B - 4.0f * A * C;
  if (discrim < 0.0f)
  {
    *tnear = -1.0f;
    *tfar = -1.0f;
  }
  const float rootDiscrim = __fsqrt_rn(discrim);
  float q;
  if (B < 0)
  {
    q = -0.5f * (B - rootDiscrim);
  }
  else
  {
    q = -0.5f * (B + rootDiscrim);
  }
  *tnear = q / A;
  *tfar = C / q;
  if (*tnear > *tfar)
  {
    float tmp = *tnear;
    *tnear = *tfar;
    *tfar = tmp;
  }
}

__device__ bool intersectBox(const Ray& ray, const float3& boxmin, const float3& boxmax,
                             float* tnear, float* tfar)
{
    // compute intersection of ray with all six bbox planes
    float3 invR = make_float3(1.0f, 1.0f, 1.0f) / ray.d;
    float t1 = (boxmin.x - ray.o.x) * invR.x;
    float t2 = (boxmax.x - ray.o.x) * invR.x;
    float tmin = fminf(t1, t2);
    float tmax = fmaxf(t1, t2);

    t1 = (boxmin.y - ray.o.y) * invR.y;
    t2 = (boxmax.y - ray.o.y) * invR.y;
    tmin = fmaxf(fminf(t1, t2), tmin);
    tmax = fminf(fmaxf(t1, t2), tmax);

    t1 = (boxmin.z - ray.o.z) * invR.z;
    t2 = (boxmax.z - ray.o.z) * invR.z;
    tmin = fmaxf(fminf(t1, t2), tmin);
    tmax = fminf(fmaxf(t1, t2), tmax);

