   // 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);
+  float orient0 = triangle_orientation(n3, n2, n1);
   // Face 1: 0 -> 2 -> 3
-  float orient1 = triangle_orientation(x[0].n, x[2].n, x[3].n);
+  float orient1 = triangle_orientation(n0, n2, n3);
   // Face 2: 3 -> 1 -> 0
-  float orient2 = triangle_orientation(x[3].n, x[1].n, x[0].n);
+  float orient2 = triangle_orientation(n3, n1, n0);
   // Face 3: 0 -> 1 -> 2
-  float orient3 = triangle_orientation(x[0].n, x[1].n, x[2].n);
+  float orient3 = triangle_orientation(n0, n1, n2);
 
   float orient01 = orient0 * orient1;
   float orient23 = orient2 * orient3;

      ;
    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;
    } 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;
    }

    // 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;
    }

    // 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();
  }
}

noperspective out vec3 integrand_over_w_back;
noperspective out vec2 texPosition;

void swap4(inout vec4 a, inout vec4 b)
{
  vec4 temp = a;
  a = b;
  b = temp;
}

void swap3(inout vec3 a, inout vec3 b)
{
  vec3 temp = a;
  a = b;
  b = temp;
}

void outputVertex(vec4 nfront, vec4 nback, vec3 ifront, vec3 iback)
{
  one_over_w_front = nfront.w;
  one_over_w_back  = nback.w;
  integrand_over_w_front = ifront * nfront.w;
  integrand_over_w_back = iback * nback.w;
  gl_Position = vec4(nfront.xy, 0.0, 1.0);
  texPosition = (nfront.xy + 1.0) / 2.0;
  EmitVertex();
}

void outputSimpleVertex(vec4 n, vec3 i)
{
  one_over_w_front = one_over_w_back  = n.w;
  integrand_over_w_front = integrand_over_w_back = i * n.w;
  gl_Position = vec4(n.xy, 0.0, 1.0);
  texPosition = (n.xy + 1.0) / 2.0;
  EmitVertex();
}

void good_case(vec4 n0, vec4 n1, vec4 n2, vec4 n3,
               vec3 i0, vec3 i1, vec3 i2, vec3 i3)
{
  // We need the point on the big face that intersects the middle point
  // in screen cooordinates.  To do this, we need to solve the following
  // system of equations in the xy plane:
  // n0 = n1 + t (n2 - n1) + u (n3 - n1)
  // Rearranging, we get:
  // (n2 - n1) t + (n3 - n1) u = n0 - n1
  // Which we write abstractly as:
  // A * X = B, where
  mat2 A = mat2((n2 - n1).xy, (n3 - n1).xy);
  vec2 B = (n0 - n1).xy;
  vec2 X = inverse(A) * B;
  float t = X[0];
  float u = X[1];
  vec4 weighted_x1 = n1 * (1 - t - u);
  vec4 weighted_x2 = n2 * t;
  vec4 weighted_x3 = n3 * u;
  vec4 ny = weighted_x1 + weighted_x2 + weighted_x3;
  vec3 iy = weighted_x1.w * i1 + weighted_x2.w * i2 + weighted_x3.w * i3;
  iy /= ny.w;
  if (n0.z > ny.z) {
    swap4(n0, ny);
    swap3(i0, iy);
  }

  // Draw three triangles.
  // x[1] - x[2] - y
  outputSimpleVertex(n1, i1);
  outputSimpleVertex(n2, i2);
  outputVertex(n0, ny, i0, iy);

  // x[2] - y - x[3]
  outputSimpleVertex(n3, i3);

  // y - x[3] - x[1]
  outputSimpleVertex(n1, i1);
  EndPrimitive();
}

void bad_case(vec4 n0, vec4 n1, vec4 n2, vec4 n3,
              vec3 i0, vec3 i1, vec3 i2, vec3 i3)
{
  // Compute the intersection point of segment 01 with segment 23
  // We need to solve this system of equations in the xy plane:
  // n0 + t (n1 - n0) = n2 + u (n3 - n2)
  // which, by algebra, is:
  // (n1 - n0) t + (n2 - n3) u = n2 - n0
  // 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((n1 - n0).xy, (n2 - n3).xy);
  vec2 B = (n2 - n0).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:
  vec4 weighted_x0 = n0 * (1 - t);
  vec4 weighted_x1 = n1 * t;
  vec4 ny0 = weighted_x0 + weighted_x1;
  vec3 iy0 = weighted_x0.w * i0 + weighted_x1.w * i1;
  iy0 /= ny0.w;

  vec4 weighted_x2 = n2 * (1 - u);
  vec4 weighted_x3 = n3 * u;
  vec4 ny1 = weighted_x2 + weighted_x3;
  vec3 iy1 = weighted_x2.w * i2 + weighted_x3.w * i3;
  iy1 /= ny1.w;

  // Make sure y0 is front and y1 is back
  if (ny0.z > ny1.z) {
    swap4(ny0, ny1);
    swap3(iy0, iy1);
  }

  // Draw four triangles.
  // x[0] - x[2] - y
  outputSimpleVertex(n0, i0);
  outputSimpleVertex(n2, i2);
  outputVertex(ny0, ny1, iy0, iy1);

  // x[2] - y - x[1]
  outputSimpleVertex(n1, i1);

  // A null triangle: y - x[1] - y
  outputVertex(ny0, ny1, iy0, iy1);

  // x[1] - y - x[3]
  outputSimpleVertex(n3, i3);

  // y - x[3] - x[0]
  outputSimpleVertex(n0, i0);
  EndPrimitive();
}

// Return +/- 1, depending on the orientation of triangle a -> b -> c
