       ;
     else if (orient02 > 0.) {
       // 0 and 2 match, so swap 1 with 2
-      //swap(x[1], x[2]);
-      vertex_data temp = x[1];
-      x[1] = x[2];
-      x[2] = temp;
+      swap4(n1, n2);
+      swap3(i1, i2);
     } else {
       // 0 and 3 match, so swap 1 with 3
-      //swap(x[1], x[3]);
-      vertex_data temp = x[1];
-      x[1] = x[3];
-      x[3] = temp;
+      swap4(n1, n3);
+      swap3(i1, i3);
     }
 
-    // Compute the intersection point of segment 01 with segment 23
-    // We need to solve this system of equations in the xy plane:
-    // x[0] + t (x[1] - x[0]) = x[2] + u (x[3] - x[2])
-    // which, by algebra, is:
-    // (x[1] - x[0]) t + (x[2] - x[3]) u = x[2] - x[0]
-    // Which we will write abstractly as:
-    // A * X = B
-    // Note that doing this calculation with x[].n instead of x[].v
-    // will give us correct x, y, z values regardless of w.
-    mat2 A = mat2((x[1].n-x[0].n).xy, (x[2].n-x[3].n).xy);
-    vec2 B = (x[2].n - x[0].n).xy;
-    vec2 X = inverse(A) * B; // NOTE: need to handle ill-conditioned case
-    float t = X[0];
-    float u = X[1];
-
-    // We get two answers; one is the point on the front edge which
-    // appears to intersect the back edge in screen coordinates, and
-    // the other is the point on the back edge which appears to inter-
-    // sect the front edge in screen coordinates.  We don't know which
-    // is which, but we can tell the difference by looking at the z
-    // values.
-    // The two intersection points are:
-    vertex_data y0, y1;
-    vec4 weighted_x0 = (1-t) * x[0].n;
-    vec4 weighted_x1 = t * x[1].n;
-    y0.n = weighted_x0 + weighted_x1;
-    y0.d = weighted_x0.w * x[0].d + weighted_x1.w * x[1].d;
-    y0.d /= y0.n.w;
-
-    vec4 weighted_x2 = (1-u) * x[2].n;
-    vec4 weighted_x3 = u * x[3].n;
-    y1.n = weighted_x2 + weighted_x3;
-    y1.d = weighted_x2.w * x[2].d + weighted_x3.w * x[3].d;
-    y1.d /= y1.n.w;
-
-    // Make sure y0 is front and y1 is back
-    if (y0.n.z > y1.n.z) {
-      //swap(y0, y1);
-      vertex_data temp = y0;
-      y0 = y1;
-      y1 = temp;
+    bad_case(n0, n1, n2, n3, i0, i1, i2, i3);
+  }
+}
+
+vec4 vector_inverse(vec4 a)
+{
+  float w = a.w;
+  a.w = 1.0;
+  a /= w;
+  return a;
+}
+
+void intersect_nearclip(out vec4 nc, out vec3 ic,
+                        vec4 px, vec3 ix, vec4 py, vec3 iy)
+{
+  float near = nearfar[0];
+  float t = (near - py.w) / (px.w - py.w);
+  nc = vector_inverse(t * px + (1-t) * py);
+  ic = t * ix + (1-t) * iy;
+}
+
+void main(void)
+{
+  vec4 p0 = gl_in[0].gl_Position;
+  vec4 p1 = gl_in[1].gl_Position;
+  vec4 p2 = gl_in[2].gl_Position;
+  vec4 p3 = gl_in[3].gl_Position;
+
+  // Number of vertices clipped by the near plane
+  int num_clipped = 0;
+
+  // Count how many vertices are clipped
+  float near = nearfar[0];
+  if (p0.w < near)
+    ++num_clipped;
+  if (p1.w < near)
+    ++num_clipped;
+  if (p2.w < near)
+    ++num_clipped;
+  if (p3.w < near)
+    ++num_clipped;
+
+  if (num_clipped == 4)
+    return;
+
+  vec4 n0 = inverted_position[0];
+  vec4 n1 = inverted_position[1];
+  vec4 n2 = inverted_position[2];
+  vec4 n3 = inverted_position[3];
+
+  vec3 i0 = integrand[0];
+  vec3 i1 = integrand[1];
+  vec3 i2 = integrand[2];
+  vec3 i3 = integrand[3];
+
+  if (num_clipped == 0) {
+    handle_tetrahedron(n0, n1, n2, n3, i0, i1, i2, i3);
+    return;
+  }
+
+  // Handle the near clipping plane. This is the "slow" execution path.
+
+  // Move clipped vertices to the front of the array using a "sorting network"
+  if (p0.w >= near && p2.w < near) {
+    swap4(p0, p2);
+    swap4(n0, n2);
+    swap3(i0, i2);
+  }
+  if (p1.w >= near && p3.w < near) {
+    swap4(p1, p3);
+    swap4(n1, n3);
+    swap3(i1, i3);
+  }
+  if (p0.w >= near && p1.w < near) {
+    swap4(p0, p1);
+    swap4(n0, n1);
+    swap3(i0, i1);
+  }
+  if (p2.w >= near && p3.w < near) {
+    swap4(p2, p3);
+    swap4(n2, n3);
+    swap3(i2, i3);
+  }
+  if (p1.w >= near && p2.w < near) {
+    swap4(p1, p2);
+    swap4(n1, n2);
+    swap3(i1, i2);
+  }
+
+  switch (num_clipped) {
+
+  case 1:
+
+    {
+      // One vertex clipped:
+      // Must determine the intersections of the 0-to-i segments
+      // with the near clipping plane, and construct three new tetrahedra
+      // to render.
+
+      vec4 n03;
+      vec3 i03;
+      intersect_nearclip(n03, i03, p0, i0, p3, i3);
+      handle_tetrahedron(n1, n2, n3, n03,
+                         i1, i2, i3, i03);
+
+      vec4 n02;
+      vec3 i02;
+      intersect_nearclip(n02, i02, p0, i0, p2, i2);
+      handle_tetrahedron(n1, n2, n02, n03,
+                         i1, i2, i02, i03);
+
+      vec4 n01;
+      vec3 i01;
+      intersect_nearclip(n01, i01, p0, i0, p1, i1);
+      handle_tetrahedron(n1, n01, n02, n03,
+                         i1, i01, i02, i03);
     }
 
-    // Draw four triangles.
-    // x[0] - x[2] - y
-    outputSimpleVertex(x[0]);
-    outputSimpleVertex(x[2]);
-    outputVertex(y0, y1);
-
-    // x[2] - y - x[1]
-    outputSimpleVertex(x[1]);
-
-    // A null triangle: y - x[1] - y
-    outputVertex(y0, y1);
-
-    // x[1] - y - x[3]
-    outputSimpleVertex(x[3]);
-
-    // y - x[3] - x[0]
-    outputSimpleVertex(x[0]);
-    EndPrimitive();
-  }
-}
+    return;
+
+  case 2:
+
+    {
+      // Two vertices clipped:
+      // Must determine the intersections of the (0,1)-to-(2,3) segments
+      // with the near clipping plane, and construct three new tetrahedra
+      // to render.
+      vec4 n02;
+      vec3 i02;
+      intersect_nearclip(n02, i02, p0, i0, p2, i2);
+
+      vec4 n12;
+      vec3 i12;
+      intersect_nearclip(n12, i12, p1, i1, p2, i2);
+
+      handle_tetrahedron(n2, n3, n02, n12,
+                         i2, i3, i02, i12);
+
+      vec4 n03;
+      vec3 i03;
+      intersect_nearclip(n03, i03, p0, i0, p3, i3);
+
+      handle_tetrahedron(n3, n02, n12, n03,
+                         i3, i02, i12, i03);
+
+      vec4 n13;
+      vec3 i13;
+      intersect_nearclip(n13, i13, p1, i1, p3, i3);
+
+      handle_tetrahedron(n3, n12, n03, n13,
+                         i3, i12, i03, i13);
+    }
+
+    return;
+
+  case 3:
+
+    // Three vertices clipped:
+    // Must replace 0-2 with their projections toward 3 at
+    // the point they intersect the near clipping plane.
+    vec4 n03;
+    vec3 i03;
+    intersect_nearclip(n03, i03, p0, i0, p3, i3);
+
+    vec4 n13;
+    vec3 i13;
+    intersect_nearclip(n13, i13, p1, i1, p3, i3);
+
+    vec4 n23;
+    vec3 i23;
+    intersect_nearclip(n23, i23, p2, i2, p3, i3);
+
+    handle_tetrahedron(n03, n13, n23, n3,
+                       i03, i13, i23, i3);
+
+    return;
+
+  default:
+
+    return;
+
+  }
+}

  // tetrahedron for our rendering calculations!
  // The induced orientations are:
  // Face 0: 3 -> 2 -> 1
  float orient0 = triangle_orientation(x[3].n, x[2].n, x[1].n);
  // Face 1: 0 -> 2 -> 3
  float orient1 = triangle_orientation(x[0].n, x[2].n, x[3].n);
  // Face 2: 3 -> 1 -> 0
  float orient2 = triangle_orientation(x[3].n, x[1].n, x[0].n);
  // Face 3: 0 -> 1 -> 2
  float orient3 = triangle_orientation(x[0].n, x[1].n, x[2].n);

  float orient01 = orient0 * orient1;
  float orient23 = orient2 * orient3;

  return sign(ab.x * ac.y - ab.y * ac.x);
}

void handle_tetrahedron(vec4 n0, vec4 n1, vec4 n2, vec4 n3,
                        vec3 i0, vec3 i1, vec3 i2, vec3 i3)
{
  // For each face of the tetrahedron, an orientation on that face
  // is induced by the orientation of the input tetrahedron.  Compute
