   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__ 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) / r.d;
-    float t1 = (boxmin.x - r.o.x) * invR.x;
-    float t2 = (boxmax.x - r.o.x) * invR.x;
+    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 - r.o.y) * invR.y;
-    t2 = (boxmax.y - r.o.y) * invR.y;
+    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 - r.o.z) * invR.z;
-    t2 = (boxmax.z - r.o.z) * invR.z;
+    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);
 

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

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

  const float4 o = mul(c_invViewMatrix, make_float4(u, v, -1.0f, 1.0f));
  const float4 d = mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 1.0f, 1.0f));

  Ray ray;
  ray.o = perspectiveDivide(o);
  ray.d = perspectiveDivide(d);
  ray.d = normalize(ray.d);
//  debug[y * width + x] = ray.o.x;
//  debug[y * width + x + 1] = ray.o.y;
//  debug[y * width + x + 2] = ray.o.z;
  float tnear;
  float tfar;
  const bool hit = intersectBox(ray, boxMin, boxMax, &tnear, &tfar);
  if (!hit)
  {
    d_output[y * width + x] = 0;
