Implemented robust ray origin calculation for bounced rays.

We take a small performance hit for this, but given that it's
making things meaningfully more correct I feel like it's more
than worth it.

src/fp_utils.rs

@@ -0,0 +1,93 @@
+//! Utilities for handling floating point precision issues
+//!
+//! This is based on the work in section 3.9 of "Physically Based Rendering:
+//! From Theory to Implementation" 3rd edition by Pharr et al.
+
+use math::{Point, Vector, Normal, dot};
+
+#[inline(always)]
+pub fn fp_gamma(n: u32) -> f32 {
+    use std::f32::EPSILON;
+    let e = EPSILON * 0.5;
+    (e * n as f32) / (1.0 - (e * n as f32))
+}
+
+
+pub fn increment_ulp(v: f32) -> f32 {
+    // Handle special cases
+    if (v.is_infinite() && v > 0.0) || v.is_nan() {
+        return v;
+    }
+    let v = if v == -0.0 { 0.0 } else { v };
+
+    // Increase ulp by 1
+    if v >= 0.0 {
+        bits_to_f32(f32_to_bits(v) + 1)
+    } else {
+        bits_to_f32(f32_to_bits(v) - 1)
+    }
+}
+
+
+pub fn decrement_ulp(v: f32) -> f32 {
+    // Handle special cases
+    if (v.is_infinite() && v < 0.0) || v.is_nan() {
+        return v;
+    }
+    let v = if v == -0.0 { 0.0 } else { v };
+
+    // Decrease ulp by 1
+    if v >= 0.0 {
+        bits_to_f32(f32_to_bits(v) - 1)
+    } else {
+        bits_to_f32(f32_to_bits(v) + 1)
+    }
+}
+
+pub fn robust_ray_origin(pos: Point, pos_err: Vector, nor: Normal, ray_dir: Vector) -> Point {
+    // Get surface normal pointing in the same
+    // direction as ray_dir.
+    let nor = {
+        let nor = nor.into_vector();
+        if dot(nor, ray_dir) >= 0.0 { nor } else { -nor }
+    };
+
+    // Calculate offset point
+    let d = dot(nor.abs(), pos_err);
+    let offset = nor * d;
+    let p = pos + offset;
+
+    // Calculate ulp offsets
+    let x = if offset.x() >= 0.0 {
+        increment_ulp(p.x())
+    } else {
+        decrement_ulp(p.x())
+    };
+
+    let y = if offset.y() >= 0.0 {
+        increment_ulp(p.y())
+    } else {
+        decrement_ulp(p.y())
+    };
+
+    let z = if offset.z() >= 0.0 {
+        increment_ulp(p.z())
+    } else {
+        decrement_ulp(p.z())
+    };
+
+    Point::new(x, y, z)
+}
+
+
+#[inline(always)]
+fn f32_to_bits(v: f32) -> u32 {
+    use std::mem::transmute_copy;
+    unsafe { transmute_copy::<f32, u32>(&v) }
+}
+
+#[inline(always)]
+fn bits_to_f32(bits: u32) -> f32 {
+    use std::mem::transmute_copy;
+    unsafe { transmute_copy::<u32, f32>(&bits) }
+}

src/main.rs

@@ -28,6 +28,7 @@ mod bbox;
 mod boundable;
 mod camera;
 mod color;
+mod fp_utils;
 mod hash;
 mod hilbert;
 mod image;

src/renderer.rs

@@ -10,10 +10,11 @@ use scoped_threadpool::Pool;
 
 use halton;
 
-use algorithm::partition_pair;
 use accel::ACCEL_TRAV_TIME;
+use algorithm::partition_pair;
 use color::{Color, XYZ, SpectralSample, map_0_1_to_wavelength};
 use float4::Float4;
+use fp_utils::robust_ray_origin;
 use hash::hash_u32;
 use hilbert;
 use image::Image;
@@ -477,14 +478,13 @@ impl LightPath {
 
                             // Calculate the shadow ray for testing if the light is
                             // in shadow or not.
-                            // TODO: use proper ray offsets for avoiding self-shadowing
-                            // rather than this hacky stupid stuff.
-                            *ray = Ray::new(
-                                idata.pos + shadow_vec.normalized() * 0.001,
+                            let offset_pos = robust_ray_origin(
+                                idata.pos,
+                                idata.pos_err,
+                                idata.nor_g.normalized(),
                                 shadow_vec,
-                                self.time,
-                                true,
                             );
+                            *ray = Ray::new(offset_pos, shadow_vec, self.time, true);
 
                             // For distant lights
                             if is_infinite {
@@ -518,12 +518,14 @@ impl LightPath {
                             self.next_attentuation_fac = filter.e / pdf;
 
                             // Calculate the ray for this bounce
-                            self.next_bounce_ray = Some(Ray::new(
-                                idata.pos + dir.normalized() * 0.0001,
+                            let offset_pos = robust_ray_origin(
+                                idata.pos,
+                                idata.pos_err,
+                                idata.nor_g.normalized(),
                                 dir,
-                                self.time,
-                                false,
-                            ));
+                            );
+                            self.next_bounce_ray =
+                                Some(Ray::new(offset_pos, dir, self.time, false));
 
                             true
                         } else {

src/surface/mod.rs

@@ -37,6 +37,7 @@ pub enum SurfaceIntersection {
 pub struct SurfaceIntersectionData {
     pub incoming: Vector, // Direction of the incoming ray
     pub pos: Point, // Position of the intersection
+    pub pos_err: Vector, // Error magnitudes of the intersection position
     pub nor: Normal, // Shading normal
     pub nor_g: Normal, // True geometric normal
     pub local_space: Matrix4x4, // Matrix from global space to local space

src/surface/triangle_mesh.rs

@@ -6,6 +6,7 @@ use accel::BVH;
 use bbox::BBox;
 use boundable::Boundable;
 use color::XYZ;
+use fp_utils::fp_gamma;
 use lerp::{lerp, lerp_slice, lerp_slice_with};
 use math::{Point, Matrix4x4, cross};
 use ray::{Ray, AccelRay};
@@ -92,17 +93,29 @@ impl<'a> Surface for TriangleMesh<'a> {
                         };
                         let mat_inv = mat_space.inverse();
                         let tri = (tri.0 * mat_inv, tri.1 * mat_inv, tri.2 * mat_inv);
-                        if let Some((t, _, _, _)) = triangle::intersect_ray(wr, tri) {
+                        if let Some((t, b0, b1, b2)) = triangle::intersect_ray(wr, tri) {
                             if t < r.max_t {
                                 if r.is_occlusion() {
                                     isects[r.id as usize] = SurfaceIntersection::Occlude;
                                     r.mark_done();
                                 } else {
+                                    // Calculate intersection point and error magnitudes
+                                    let pos = ((tri.0.into_vector() * b0)
+                                        + (tri.1.into_vector() * b1)
+                                        + (tri.2.into_vector() * b2)).into_point();
+
+                                    let pos_err = ((tri.0.into_vector().abs() * b0)
+                                        + (tri.1.into_vector().abs() * b1)
+                                        + (tri.2.into_vector().abs() * b2))
+                                        * fp_gamma(7);
+
+                                    // Fill in intersection data
                                     isects[r.id as usize] = SurfaceIntersection::Hit {
                                         intersection_data: SurfaceIntersectionData {
                                             incoming: wr.dir,
                                             t: t,
-                                            pos: wr.orig + (wr.dir * t),
+                                            pos: pos,
+                                            pos_err: pos_err,
                                             nor: cross(tri.0 - tri.1, tri.0 - tri.2)
                                                 .into_normal(), // TODO
                                             nor_g: cross(tri.0 - tri.1, tri.0 - tri.2)

sub_crates/math3d/src/vector.rs

@@ -6,7 +6,7 @@ use std::ops::{Add, Sub, Mul, Div, Neg};
 use float4::Float4;
 
 use super::{DotProduct, CrossProduct};
-use super::{Matrix4x4, Normal};
+use super::{Matrix4x4, Point, Normal};
 
 
 /// A direction vector in 3d homogeneous space.
@@ -36,6 +36,16 @@ impl Vector {
         *self / self.length()
     }
 
+    #[inline(always)]
+    pub fn abs(&self) -> Vector {
+        Vector::new(self.x().abs(), self.y().abs(), self.z().abs())
+    }
+
+    #[inline(always)]
+    pub fn into_point(self) -> Point {
+        Point::new(self.x(), self.y(), self.z())
+    }
+
     #[inline(always)]
     pub fn into_normal(self) -> Normal {
         Normal::new(self.x(), self.y(), self.z())