Fixed sampling of very small rectangle lights.

The sampling method used before is numerically unstable for very
small lights.  That sampling method is still used for large/close
lights, since it works very well for that.  But for small/distant
lights a simpler and numerically stable method is used.

src/light/rectangle_light.rs

@@ -4,9 +4,10 @@ use bbox::BBox;
 use boundable::Boundable;
 use color::{XYZ, SpectralSample, Color};
 use lerp::lerp_slice;
-use math::{Vector, Normal, Point, Matrix4x4, cross};
+use math::{Vector, Normal, Point, Matrix4x4, cross, dot};
 use ray::{Ray, AccelRay};
-use sampling::{spherical_triangle_solid_angle, uniform_sample_spherical_triangle};
+use sampling::{spherical_triangle_solid_angle, uniform_sample_spherical_triangle,
+               triangle_surface_area, uniform_sample_triangle};
 use shading::surface_closure::{SurfaceClosureUnion, EmitClosure};
 use shading::SurfaceShader;
 use surface::{Surface, SurfaceIntersection, SurfaceIntersectionData, triangle};
@@ -14,6 +15,9 @@ use surface::{Surface, SurfaceIntersection, SurfaceIntersectionData, triangle};
 use super::SurfaceLight;
 
 
+const SIMPLE_SAMPLING_THRESHOLD: f32 = 0.01;
+
+
 #[derive(Copy, Clone, Debug)]
 pub struct RectangleLight<'a> {
     dimensions: &'a [(f32, f32)],
@@ -50,13 +54,12 @@ impl<'a> RectangleLight<'a> {
         space: &Matrix4x4,
         arr: Point,
         sample_dir: Vector,
-        sample_u: f32,
+        hit_point: Point,
-        sample_v: f32,
         wavelength: f32,
         time: f32,
     ) -> f32 {
         // We're not using these, silence warnings
-        let _ = (sample_dir, sample_u, sample_v, wavelength);
+        let _ = wavelength;
 
         let dim = lerp_slice(self.dimensions, time);
 
@@ -78,8 +81,19 @@ impl<'a> RectangleLight<'a> {
         let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3);
         let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3);
 
+        // World-space surface normal
+        let normal = Normal::new(0.0, 0.0, 1.0) * space_inv;
+
+        // PDF
+        if (area_1 + area_2) < SIMPLE_SAMPLING_THRESHOLD {
+            let area = triangle_surface_area(p2, p1, p3) + triangle_surface_area(p4, p1, p3);
+            (hit_point - arr).length2() /
+                dot(sample_dir.normalized(), normal.into_vector().normalized()).abs() /
+                area
+        } else {
             1.0 / (area_1 + area_2)
         }
+    }
 
     // fn outgoing(
     //     &self,
@@ -117,7 +131,8 @@ impl<'a> SurfaceLight for RectangleLight<'a> {
         let dim = lerp_slice(self.dimensions, time);
         let col = lerp_slice(self.colors, time);
 
-        let surface_area_inv: f64 = 1.0 / (dim.0 as f64 * dim.1 as f64);
+        let surface_area: f64 = dim.0 as f64 * dim.1 as f64;
+        let surface_area_inv: f64 = 1.0 / surface_area;
 
         // Get the four corners of the rectangle, transformed into world space
         let space_inv = space.inverse();
@@ -126,19 +141,69 @@ impl<'a> SurfaceLight for RectangleLight<'a> {
         let p3 = Point::new(dim.0 * -0.5, dim.1 * -0.5, 0.0) * space_inv;
         let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv;
 
-        // Get the four corners of the rectangle, projected on to the unit
+        // Get the four corners of the rectangle relative to arr.
-        // sphere centered around arr.
+        let lp1 = p1 - arr;
-        let sp1 = (p1 - arr).normalized();
+        let lp2 = p2 - arr;
-        let sp2 = (p2 - arr).normalized();
+        let lp3 = p3 - arr;
-        let sp3 = (p3 - arr).normalized();
+        let lp4 = p4 - arr;
-        let sp4 = (p4 - arr).normalized();
+
+        // Four corners projected on to the unit sphere.
+        let sp1 = lp1.normalized();
+        let sp2 = lp2.normalized();
+        let sp3 = lp3.normalized();
+        let sp4 = lp4.normalized();
 
         // Get the solid angles of the rectangle split into two triangles
         let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3);
         let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3);
 
+        // Calculate world-space surface normal
+        let normal = Normal::new(0.0, 0.0, 1.0) * space_inv;
+
+        if (area_1 + area_2) < SIMPLE_SAMPLING_THRESHOLD {
+            // Simple sampling for more distant lights
+            let surface_area_1 = triangle_surface_area(p2, p1, p3);
+            let surface_area_2 = triangle_surface_area(p4, p1, p3);
+            let sample_point = {
+                // Select which triangle to sample
+                let threshhold = surface_area_1 / (surface_area_1 + surface_area_2);
+                if u < threshhold {
+                    uniform_sample_triangle(
+                        p2.into_vector(),
+                        p1.into_vector(),
+                        p3.into_vector(),
+                        v,
+                        u / threshhold,
+                    )
+                } else {
+                    uniform_sample_triangle(
+                        p4.into_vector(),
+                        p1.into_vector(),
+                        p3.into_vector(),
+                        v,
+                        (u - threshhold) / (1.0 - threshhold),
+                    )
+                }
+            }.into_point();
+            let shadow_vec = sample_point - arr;
+            let spectral_sample =
+                (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength);
+            let pdf = (sample_point - arr).length2() /
+                dot(shadow_vec.normalized(), normal.into_vector().normalized()).abs() /
+                (surface_area_1 + surface_area_2);
+            let point_err = 0.0001; // TODO: this is a hack, do properly.
+            (spectral_sample, (sample_point, normal, point_err), pdf)
+        } else {
+            // Sophisticated sampling for close lights.
+
             // Normalize the solid angles for selection purposes
-        let prob_1 = area_1 / (area_1 + area_2);
+            let prob_1 = if area_1.is_infinite() {
+                1.0
+            } else if area_2.is_infinite() {
+                0.0
+            } else {
+                area_1 / (area_1 + area_2)
+            };
             let prob_2 = 1.0 - prob_1;
 
             // Select one of the triangles and sample it
@@ -161,12 +226,12 @@ impl<'a> SurfaceLight for RectangleLight<'a> {
                 sample_point_local.set_z(0.0);
             }
             let sample_point = sample_point_local * space_inv;
-        let normal = Normal::new(0.0, 0.0, 1.0) * space_inv;
             let point_err = 0.0001; // TODO: this is a hack, do properly.
 
             // Calculate pdf and light energy
             let pdf = 1.0 / (area_1 + area_2); // PDF of the ray direction being sampled
-        let spectral_sample = (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength);
+            let spectral_sample =
+                (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength);
 
             (
                 spectral_sample,
@@ -174,6 +239,7 @@ impl<'a> SurfaceLight for RectangleLight<'a> {
                 pdf as f32,
             )
         }
+    }
 
     fn is_delta(&self) -> bool {
         false
@@ -239,8 +305,7 @@ impl<'a> Surface for RectangleLight<'a> {
                                     &xform,
                                     wr.orig,
                                     wr.dir,
-                                    0.0,
+                                    pos,
-                                    0.0,
                                     wr.wavelength,
                                     r.time,
                                 ),

src/sampling/mod.rs

@@ -2,4 +2,5 @@ mod monte_carlo;
 
 pub use self::monte_carlo::{square_to_circle, cosine_sample_hemisphere, uniform_sample_hemisphere,
                             uniform_sample_sphere, uniform_sample_cone, uniform_sample_cone_pdf,
+                            uniform_sample_triangle, triangle_surface_area,
                             spherical_triangle_solid_angle, uniform_sample_spherical_triangle};

src/sampling/monte_carlo.rs

@@ -4,7 +4,7 @@ use std::f32::consts::FRAC_PI_4 as QPI_32;
 use std::f32::consts::PI as PI_32;
 use std::f64::consts::PI as PI_64;
 
-use math::{Vector, dot};
+use math::{Vector, Point, cross, dot};
 
 
 /// Maps the unit square to the unit circle.
@@ -80,6 +80,22 @@ pub fn uniform_sample_cone_pdf(cos_theta_max: f64) -> f64 {
     1.0 / (2.0 * PI_64 * (1.0 - cos_theta_max))
 }
 
+/// Generates a uniform sample on a triangle given two uniform random
+/// variables i and j in [0, 1].
+pub fn uniform_sample_triangle(va: Vector, vb: Vector, vc: Vector, i: f32, j: f32) -> Vector {
+    let isqrt = i.sqrt();
+    let a = 1.0 - isqrt;
+    let b = isqrt * (1.0 - j);
+    let c = j * isqrt;
+
+    (va * a) + (vb * b) + (vc * c)
+}
+
+/// Calculates the surface area of a triangle.
+pub fn triangle_surface_area(p0: Point, p1: Point, p2: Point) -> f32 {
+    0.5 * cross(p1 - p0, p2 - p0).length()
+}
+
 /// Calculates the projected solid angle of a spherical triangle.
 ///
 /// A, B, and C are the points of the triangle on a unit sphere.