Add the orbital-explorer shaders that Paul sent me.

1 files changed, 636 insertions, 0 deletions

diff --git a/shaders/orbital_explorer.shader_test b/shaders/orbital_explorer.shader_test
new file mode 100644
index 0000000..70ea084
--- /dev/null
+++ b/shaders/orbital_explorer.shader_test
@@ -0,0 +1,636 @@
+# Shader_runner test illustrating slow compilation of geometry shaders
+# in mesa/i965.
+#
+# I've made the following changes from the original shaders:
+# - Change shaderDepth from a texture to a uniform float
+#
+# The shaders in this test come from
+# https://github.com/bjthinks/orbital-explorer, which contains this
+# copyright notice:
+#
+# This file is part of the Electron Orbital Explorer. The Electron
+# Orbital Explorer is distributed under the Simplified BSD License
+# (also called the "BSD 2-Clause License"), in hopes that these
+# rendering techniques might be used by other programmers in
+# applications such as scientific visualization, video gaming, and so
+# on. If you find value in this software and use its technologies for
+# another purpose, I would love to hear back from you at bjthinks (at)
+# gmail (dot) com. If you improve this software and agree to release
+# your modifications under the below license, I encourage you to fork
+# the development tree on github and push your modifications. The
+# Electron Orbital Explorer's development URL is:
+# https://github.com/bjthinks/orbital-explorer
+# (This paragraph is not part of the software license and may be
+# removed.)
+#
+# Copyright (c) 2013, Brian W. Johnson
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# + Redistributions of source code must retain the above copyright
+#   notice, this list of conditions and the following disclaimer.
+#
+# + Redistributions in binary form must reproduce the above copyright
+#   notice, this list of conditions and the following disclaimer in
+#   the documentation and/or other materials provided with the
+#   distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+
+[require]
+GLSL >= 1.50
+
+[vertex shader]
+#version 150
+
+uniform mat4 modelViewProjMatrix;
+
+in vec4 position;
+in vec3 rim;
+out vec4 inverted_position;
+out vec3 integrand;
+
+// Note sure what this transform is called, but since doing it twice
+// gives back the original vector, we'll call it a "vector inverse"
+// for now.
+// (x, y, z, w) --> (x/w, y/w, z/w, 1/w)
+vec4 vector_inverse(vec4 a)
+{
+  float w = a.w;
+  a.w = 1.0;
+  a /= w;
+  return a;
+}
+
+// Calculate the position and inverted position of the vertex, setting
+// the output varyings "gl_Position" and "inverted_position".
+void calculate_position()
+{
+  vec4 pos = modelViewProjMatrix * position;
+  gl_Position = pos;
+  inverted_position = vector_inverse(pos);
+}
+
+// Apply brightness adjustment and color rotation to the input color,
+// setting the output varying "integrand" to transformed color
+// coordinates that are sensible to integrate.
+void calculate_color()
+{
+  integrand = rim;
+}
+
+void main(void)
+{
+  calculate_position();
+  calculate_color();
+}
+
+
+[geometry shader]
+#version 150
+
+layout(lines_adjacency) in;
+layout(triangle_strip,max_vertices=21) out;
+
+uniform vec2 nearfar;
+
+in vec4 inverted_position[4];
+in vec3 integrand[4];
+
+noperspective out float one_over_w_front;
+noperspective out float one_over_w_back;
+noperspective out vec3 integrand_over_w_front;
+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
+// in the xy plane (ignoring z and w)
+float triangle_orientation(vec4 a, vec4 b, vec4 c)
+{
+  vec2 ab = b.xy - a.xy;
+  vec2 ac = c.xy - a.xy;
+  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
+  // whether those orientations are equal to or opposite from the
+  // orientation in screen coordinates.
+  // The orientation will be the same for front faces and reversed
+  // for back faces (or vice versa, depending on whether the
+  // tetrahedron is left or right handed).
+  // Surprisingly, we DON'T need to know the handedness of the input
+  // tetrahedron for our rendering calculations!
+  // The induced orientations are:
+  // Face 0: 3 -> 2 -> 1
+  float orient0 = triangle_orientation(n3, n2, n1);
+  // Face 1: 0 -> 2 -> 3
+  float orient1 = triangle_orientation(n0, n2, n3);
+  // Face 2: 3 -> 1 -> 0
+  float orient2 = triangle_orientation(n3, n1, n0);
+  // Face 3: 0 -> 1 -> 2
+  float orient3 = triangle_orientation(n0, n1, n2);
+
+  float orient01 = orient0 * orient1;
+  float orient23 = orient2 * orient3;
+  float orient = orient01 * orient23;
+
+  // For now, discard tetrahedra where a face is seen exactly edge-on.
+  // The sign() function should return 0.0 in that case, but exercise
+  // an abundance of caution with floating point math.
+  if (orient > -0.5 && orient < 0.5)
+    return;
+
+  // We call a tetrahedron good if its projection onto the canvas
+  // has a triangular convex hull with one vertex in the interior.
+  // It is conversely bad if the convex hull is a quadrilateral.
+  // A tetrahedron is good iff an odd number of faces have flipped
+  // orientation, i.e. there are either three front and one rear
+  // faces or vice versa.
+  bool good = orient < 0.;
+
+  if (good) {
+    // The vertex in the center is the one whose orientN value is
+    // of a different sign.  Figure out the majority sign and which
+    // vertex is in the middle.
+    if (orient01 > 0.) {
+      // 0 and 1 have the same sign; either 2 or 3 is opposite
+      float majority = orient0;
+      if (orient2 * majority < 0.) {
+        swap4(n0, n2);
+        swap3(i0, i2);
+      } else {
+        swap4(n0, n3);
+        swap3(i0, i3);
+      }
+    } else {
+      // 2 and 3 have the same sign; either 0 or 1 is opposite
+      float majority = orient2;
+      if (orient0 * majority < 0.)
+        ;
+      else {
+        swap4(n0, n1);
+        swap3(i0, i1);
+      }
+    }
+
+    good_case(n0, n1, n2, n3,
+              i0, i1, i2, i3);
+
+  } else {
+    float orient02 = orient0 * orient2;
+
+    // We assume here that two of the orientations are positive, and
+    // two are negative.  The other possibility, that all four have
+    // the same sign, would indicate an error; it's unclear if that
+    // can happen at all, but it's probably worth checking for.
+    if (orient01 > 0. && orient23 > 0. && orient02 > 0.)
+         return;
+
+    // Reorder so that vertices 0 and 1 have the same sign
+    if (orient01 > 0.)
+      // 0 and 1 already match, do nothng
+      ;
+    else if (orient02 > 0.) {
+      // 0 and 2 match, so swap 1 with 2
+      swap4(n1, n2);
+      swap3(i1, i2);
+    } else {
+      // 0 and 3 match, so swap 1 with 3
+      swap4(n1, n3);
+      swap3(i1, i3);
+    }
+
+    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);
+    }
+
+    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;
+
+  }
+}
+
+[fragment shader]
+#version 150
+
+noperspective in float one_over_w_front;
+noperspective in float one_over_w_back;
+noperspective in vec3 integrand_over_w_front;
+noperspective in vec3 integrand_over_w_back;
+noperspective in vec2 texPosition;
+out vec3 integratedValue;
+uniform float solidDepth; // was: uniform sampler2D solidDepth;
+uniform vec2 nearfar;
+uniform bool depth_obscuration;
+
+// Extract the w value (which is the pre-projection z value) of the
+// solid object from the depth buffer that was used in the solid
+// object rendering pass.
+float compute_solid_w()
+{
+  // The depth buffer from the solid object rendering pass stores
+  // something called z_w.
+  float z_w = solidDepth; // was: texture(solidDepth, texPosition).x;
+
+  // We also need the "near" and "far" values from the frustum matrix,
+  // which are passed in as uniforms.
+  float n = nearfar[0];
+  float f = nearfar[1];
+
+  // Also, z_d = z_c / w_c
+  // From the frustum matrix,
+  // z_c = -(f + n) / (f - n) * z - 2 * f * n / (f - n) * w
+  // w_c = -z
+  // Substituting,
+  // z_d = z_c / w_c = (f + n) / (f - n) + 2 * f * n / (f - n) * w / z
+  // But pre-frustum w = 1, so
+  // z_d = (f + n) / (f - n) + 2 * f * n / (f - n) / z
+  // z_d * (f - n) = f + n + 2 * f * n / z
+  // z_d * (f - n) - f - n = 2 * f * n / z
+  // z = 2 * f * n / (z_d * (f - n) - f - n)
+  // w = -z = -2 * f * n / (z_d * (f - n) - f - n)
+  // Substitute z_d = 2 * z_w - 1
+  // w = -2 * f * n / ((2 * z_w - 1) * (f - n) - f - n)
+  // Simplify RHS
+  // w = -2 * f * n / (2 * z_w * (f - n) - 2 * f)
+  // w = f * n / (f - z_w * (f - n))
+  return f * n / (f - z_w * (f - n));
+}
+
+void main(void)
+{
+  float w_front = 1.0 / one_over_w_front;
+  float w_back = 1.0 / one_over_w_back;
+
+  vec3 integrand_front = w_front * integrand_over_w_front;
+  vec3 integrand_back = w_back * integrand_over_w_back;
+
+  float solid_w = compute_solid_w();
+
+  if (solid_w <= w_front)
+    discard;
+
+  if (solid_w < w_back) {
+    float t = (solid_w - w_front) / (w_back - w_front);
+    integrand_back = (1-t) * integrand_front + t * integrand_back;
+    w_back = solid_w;
+  }
+
+  float depth = w_back - w_front;
+  vec3 integrand_middle = (integrand_front + integrand_back) / 2.0;
+  if (!depth_obscuration) {
+    integratedValue = depth * integrand_middle;
+  } else {
+    // Assume the z-component of the integrand represents magnitude,
+    // and calculate the falloff of all components based on its
+    // integral along the line of sight.
+    // TODO
+  }
+}
+
+[vertex data]
+position/float/4  rim/float/3
+-1 -1 -2 1        1 1 1
+-1  1 -3 1        1 1 1
+ 1 -1 -3 1        1 1 1
+ 1  1 -2 1        1 1 1
+
+[test]
+# Frustum with n=1 f=10 l=-1 r=1 b=-1 t=1
+uniform mat4 modelViewProjMatrix 1 0 0 0  0 1 0 0  0 0 -1.222 -1  0 0 -2.222 0
+uniform vec2 nearfar 1 10
+uniform float solidDepth 1.0
+uniform int depth_obscuration 0
+clear color 0 0 0 0
+clear
+draw arrays GL_LINES_ADJACENCY 0 4