Add multiple importance sampling tutorial/demo

res/org/lwjgl/demo/opengl/raytracing/tutorial3/random.glsl

@@ -136,16 +136,15 @@ vec3 isotropic(float rp, float c) {
  *          value
  */
 vec4 randomHemispherePoint(vec3 n, spatialrand rand) {
-  float c = rand.y;
-  return vec4(around(isotropic(rand.x, c), n), ONE_OVER_2PI);
+  return vec4(around(isotropic(rand.x, rand.y), n), ONE_OVER_2PI);
 }
 
 /**
- * Generate a cosine-weighted random vector on the hemisphere around
- * the given normal vector 'n'.
+ * Generate a cosine-weighted random vector on the hemisphere around the
+ * given normal vector 'n'.
  *
- * The probability density of any vector is directly proportional to
- * the cosine of the angle between that vector and the given normal 'n'.
+ * The probability density of any vector is directly proportional to the
+ * cosine of the angle between that vector and the given normal 'n'.
  *
  * http://www.rorydriscoll.com/2009/01/07/better-sampling/
  *
@@ -163,10 +162,10 @@ vec4 randomCosineWeightedHemispherePoint(vec3 n, spatialrand rand) {
 /**
  * Generate Phong-weighted random vector around the given reflection
  * vector 'r'.
- * Since the Phong BRDF has higher values when the outgoing vector
- * is close to the perfect reflection vector of the incoming vector
- * across the normal, we generate directions primarily around that
- * reflection vector.
+ * Since the Phong BRDF has higher values when the outgoing vector is
+ * close to the perfect reflection vector of the incoming vector across
+ * the normal, we generate directions primarily around that reflection
+ * vector.
  *
  * http://blog.tobias-franke.eu/2014/03/30/notes_on_importance_sampling.html
  *
@@ -176,6 +175,6 @@ vec4 randomCosineWeightedHemispherePoint(vec3 n, spatialrand rand) {
  * @returns the Phong-weighted random vector
  */
 vec4 randomPhongWeightedHemispherePoint(vec3 r, float a, spatialrand rand) {
-  float ai = 1.0 / (a + 1.0), c = pow(rand.y, ai);
-  return vec4(around(isotropic(rand.x, c), r), (a + 1.0) * pow(rand.y, a * ai) * ONE_OVER_2PI);
+  float ai = 1.0 / (a + 1.0), pr = (a + 1.0) * pow(rand.y, a * ai) * ONE_OVER_2PI;
+  return vec4(around(isotropic(rand.x, pow(rand.y, ai)), r), pr);
 }

res/org/lwjgl/demo/opengl/raytracing/tutorial3/raytracing.glsl

@@ -118,6 +118,11 @@ ivec2 px;
  * box 'b' and return the values of the parameter 't' at which the ray
  * enters and exists the box, called (tNear, tFar). If there is no
  * intersection then tNear > tFar or tFar < 0.
+ *
+ * @param origin the ray's origin position
+ * @param dir the ray's direction vector
+ * @param b the box to test
+ * @returns (tNear, tFar)
  */
 vec2 intersectBox(vec3 origin, vec3 dir, const box b) {
   vec3 tMin = (b.min - origin) / dir;
@@ -422,7 +427,7 @@ vec3 trace(vec3 origin, vec3 dir) {
       att *= brdf(b, s.xyz, -dir, normal);
     }
     /*
-     * Use the newly generate direction vector as the next one to
+     * Use the newly generated direction vector as the next one to
      * follow.
      */
     dir = s.xyz;

res/org/lwjgl/demo/opengl/raytracing/tutorial4/quad.fs.glsl

@@ -0,0 +1,16 @@
+/*
+ * Copyright LWJGL. All rights reserved.
+ * License terms: http://lwjgl.org/license.php
+ */
+/* The texture we are going to sample */
+uniform sampler2D tex;
+
+/* This comes interpolated from the vertex shader */
+in vec2 texcoord;
+
+out vec4 color;
+
+void main(void) {
+  /* Well, simply sample the texture */
+  color = texture2D(tex, texcoord);
+}

res/org/lwjgl/demo/opengl/raytracing/tutorial4/quad.vs.glsl

@@ -0,0 +1,18 @@
+/*
+ * Copyright LWJGL. All rights reserved.
+ * License terms: http://lwjgl.org/license.php
+ */
+/* The position of the vertex as two-dimensional vector */
+in vec2 vertex;
+
+/* Write interpolated texture coordinate to fragment shader */
+out vec2 texcoord;
+
+void main(void) {
+  gl_Position = vec4(vertex, 0.0, 1.0);
+  /*
+   * Compute texture coordinate by simply
+   * interval-mapping from [-1..+1] to [0..1]
+   */
+  texcoord = vertex * 0.5 + vec2(0.5, 0.5);
+}

res/org/lwjgl/demo/opengl/raytracing/tutorial4/random.glsl

@@ -0,0 +1,227 @@
+/*
+ * Copyright LWJGL. All rights reserved.
+ * License terms: http://lwjgl.org/license.php
+ */
+#version 330 core
+
+#define PI 3.14159265359
+#define TWO_PI 6.28318530718
+#define ONE_OVER_PI (1.0 / PI)
+#define ONE_OVER_2PI (1.0 / TWO_PI)
+
+/*
+ * Define the type of a random variable/vector which is used to compute
+ * spatial values (positions, directions, angles, etc.).
+ * Actually this should have been a typedef, but we can't do this in
+ * GLSL.
+ */
+#define spatialrand vec2
+
+/**
+ * Compute an arbitrary unit vector orthogonal to the unit vector 'v'
+ * which is either of the three base coordinate axes (or the negation of
+ * either of them).
+ *
+ * http://lolengine.net/blog/2013/09/21/picking-orthogonal-vector-combing-coconuts
+ *
+ * @param v the unit base vector to compute an orthogonal unit vector
+ *          from
+ * @returns the unit vector orthogonal to 'v'
+ */
+vec3 orthoBase(vec3 v) {
+  return abs(v.x) > abs(v.z) ? vec3(-v.y, v.x, 0.0) : vec3(0.0, -v.z, v.y);
+}
+
+/**
+ * Compute an arbitrary unit vector orthogonal to any vector 'v'.
+ *
+ * @param v the vector to compute an orthogonal unit vector from
+ * @returns the unit vector orthogonal to 'v'
+ */
+vec3 ortho(vec3 v) {
+  return normalize(cross(v, orthoBase(v)));
+}
+
+/**
+ * http://amindforeverprogramming.blogspot.de/2013/07/random-floats-in-glsl-330.html?showComment=1507064059398#c5427444543794991219
+ */
+uint hash3(uint x, uint y, uint z) {
+  x += x >> 11;
+  x ^= x << 7;
+  x += y;
+  x ^= x << 3;
+  x += z ^ (x >> 14);
+  x ^= x << 6;
+  x += x >> 15;
+  x ^= x << 5;
+  x += x >> 12;
+  x ^= x << 9;
+  return x;
+}
+
+/**
+ * Generate a floating-point pseudo-random number in [0, 1).
+ *
+ * With Monte Carlo sampling we need random numbers to generate random
+ * vectors for shooting rays.
+ * This function takes a vector of three arbitrary floating-point
+ * numbers and based on those numbers returns a pseudo-random number in
+ * the range [0..1).
+ * Since in GLSL we have no other state/entropy than what we input as
+ * uniforms or compute from varying variables (e.g. the invocation id of
+ * the current work item), any "pseudo-random" number will necessarily
+ * only be some kind of hash over the input.
+ * This function however has very very good properties exhibiting no
+ * patterns in the distribution of the random numbers, no matter the
+ * magnitude of the input values.
+ * 
+ * http://amindforeverprogramming.blogspot.de/2013/07/random-floats-in-glsl-330.html
+ *
+ * @param f some input vector of numbers to generate a single
+ *          pseudo-random number from
+ * @returns a single pseudo-random number in [0, 1)
+ */
+float random(vec3 f) {
+  uint mantissaMask = 0x007FFFFFu;
+  uint one = 0x3F800000u;
+  uvec3 u = floatBitsToUint(f);
+  uint h = hash3(u.x, u.y, u.z);
+  return uintBitsToFloat((h & mantissaMask) | one) - 1.0;
+}
+
+/**
+ * Transform the given vector 'v' from its local frame into the
+ * orthonormal basis with +Z = z.
+ *
+ * @param v the vector to transform
+ * @param z the direction to transform +Z into
+ * @returns the vector in the new basis
+ */
+vec3 around(vec3 v, vec3 z) {
+  vec3 t = ortho(z), b = cross(z, t);
+  return t * v.x + b * v.y + z * v.z;
+}
+
+/**
+ * Compute the cartesian coordinates from the values typically computed
+ * when generating sample directions/angles for isotropic BRDFs.
+ *
+ * @param rp this is a pseudo-random number in [0, 1) representing the
+ *           angle `phi` to rotate around the principal vector. For
+ *           isotropic BRDFs where `phi` has the same probability of
+ *           being chosen, this will always be a random number in [0, 1)
+ *           that can directly be used from rand.x given to all sample
+ *           direction generation functions
+ * @param c this is the cosine of theta, the angle in [0, PI/2) giving
+ *          the angle between the principal vector and the one to
+ *          generate. Being the cosine of an angle in [0, PI/2) the
+ *          value 'c' is itself in the range [0, 1)
+ * @returns the cartesian direction vector
+ */
+vec3 isotropic(float rp, float c) {
+  float p = TWO_PI * rp, s = sqrt(1.0 - c*c);
+  return vec3(cos(p) * s, sin(p) * s, c);
+}
+
+/**
+ * Generate a uniformly distributed random vector on the hemisphere
+ * around the given normal vector 'n'.
+ *
+ * http://mathworld.wolfram.com/SpherePointPicking.html
+ *
+ * @param n the normal vector determining the direction of the
+ *          hemisphere
+ * @param rand a vector of two floating-point pseudo-random numbers
+ * @returns the random hemisphere vector plus its probability density
+ *          value
+ */
+vec4 randomHemispherePoint(vec3 n, spatialrand rand) {
+  return vec4(around(isotropic(rand.x, rand.y), n), ONE_OVER_2PI);
+}
+/**
+ * Compute the probability density value of randomHemispherePoint()
+ * generating the vector 'v'.
+ *
+ * @param n the normal vector determining the direction of the
+ *          hemisphere
+ * @param v the v vector as returned by randomHemispherePoint()
+ * @returns pdf(v) for the uniform hemisphere distribution
+ */
+float hemisphereProbability(vec3 n, vec3 v) {
+  return dot(v, n) <= 0.0 ? 0.0 : ONE_OVER_2PI;
+}
+
+/**
+ * Generate a cosine-weighted random vector on the hemisphere around the
+ * given normal vector 'n'.
+ *
+ * The probability density of any vector is directly proportional to the
+ * cosine of the angle between that vector and the given normal 'n'.
+ *
+ * http://www.rorydriscoll.com/2009/01/07/better-sampling/
+ *
+ * @param n the normal vector determining the direction of the
+ *          hemisphere
+ * @param rand a vector of two floating-point pseudo-random numbers
+ * @returns the cosine-weighted random hemisphere vector plus its
+ *          probability density value
+ */
+vec4 randomCosineWeightedHemispherePoint(vec3 n, spatialrand rand) {
+  float c = sqrt(rand.y);
+  return vec4(around(isotropic(rand.x, c), n), c * ONE_OVER_PI);
+}
+
+/**
+ * Generate Phong-weighted random vector around the given reflection
+ * vector 'r'.
+ * Since the Phong BRDF has higher values when the outgoing vector is
+ * close to the perfect reflection vector of the incoming vector across
+ * the normal, we generate directions primarily around that reflection
+ * vector.
+ *
+ * http://blog.tobias-franke.eu/2014/03/30/notes_on_importance_sampling.html
+ *
+ * @param r the direction of perfect reflection
+ * @param a the power to raise the cosine term to in the Phong model
+ * @param rand a vector of two pseudo-random numbers
+ * @returns the Phong-weighted random vector
+ */
+vec4 randomPhongWeightedHemispherePoint(vec3 r, float a, spatialrand rand) {
+  float ai = 1.0 / (a + 1.0), pr = (a + 1.0) * pow(rand.y, a * ai) * ONE_OVER_2PI;
+  return vec4(around(isotropic(rand.x, pow(rand.y, ai)), r), pr);
+}
+
+/**
+ * Generate a random vector uniformly distributed over the disk with 'd'
+ * pointing at its center and with the radius 'r' at a distance of 'd'.
+ *
+ * https://en.wikipedia.org/wiki/Solid_angle#Cone,_spherical_cap,_hemisphere
+ * https://en.wikipedia.org/wiki/Angular_diameter#Formula
+ *
+ * @param n the normalize direction vector towards the disk's center
+ * @param d the distance to the disk
+ * @param r the radius of the disk
+ * @param rand a vector of two pseudo-random numbers
+ * @returns the random disk sample vector
+ */
+vec4 randomDiskPoint(vec3 n, float d, float r, spatialrand rand) {
+  float D = r/d, c = inversesqrt(1 + D * D), pr = ONE_OVER_2PI / (1.0 - c);
+  return vec4(around(isotropic(rand.x, 1.0 - rand.y * (1.0 - c)), n), pr);
+}
+/**
+ * Compute the probability density value of randomDiskPoint()
+ * generating the z vector 'z'.
+ *
+ * https://en.wikipedia.org/wiki/Solid_angle#Cone,_spherical_cap,_hemisphere
+ * https://en.wikipedia.org/wiki/Angular_diameter#Formula
+ *
+ * @param the normalize direction vector towards the disk's center
+ * @param d the distance to the disk
+ * @param r the radius of the disk
+ * @param v the vector as returned by randomDiskPoint()
+ * @returns pdf(v) for the disk distribution
+ */
+float diskProbability(vec3 n, float d, float r, vec3 v) {
+  float D = r/d, c = inversesqrt(1 + D * D), pr = ONE_OVER_2PI / (1.0 - c);
+  return dot(n, v) <= c ? 0.0 : pr;
+}