Fixed self-intersection bug for coordinates near zero.

This turned out to be a rather interesting one.  The water-tight
ray/triangle intersection algorithm, while very accurate for
finding if there is an intersection with a line segment, is
not as remarkably accurate for determining if that intersection
is within the interval of the ray.

This is because of the coordinate transformation it does
depending on ray direction: for triangles laying flat on one of
the axis planes near zero, that near-zero coordinate can get
transformed to a much less accurate space for testing.  In fact,
generally speaking, beause of the coordinate transform, you can
only rely on the test being as accurate as the least accurate
axis.

The ray-origin offset code was doing offsets based on the
assumption that the error on the major axes are independent, but
as this triangle intersection algorithm shows, you can't actually
depend on that being the case.  So rather than handling triangle
intersection as a special case, I've changed the intersection
position error to be a single float, representing the maximum
possible error on any axis.  This should be robust for any
geometry type added in the future, and also solves the immediate
issue in a correct way.

src/fp_utils.rs

@@ -48,7 +48,7 @@ pub fn decrement_ulp(v: f32) -> f32 {
     }
 }
 
-pub fn robust_ray_origin(pos: Point, pos_err: Vector, nor: Normal, ray_dir: Vector) -> Point {
+pub fn robust_ray_origin(pos: Point, pos_err: f32, nor: Normal, ray_dir: Vector) -> Point {
     // Get surface normal pointing in the same
     // direction as ray_dir.
     let nor = {
@@ -57,39 +57,27 @@ pub fn robust_ray_origin(pos: Point, pos_err: Vector, nor: Normal, ray_dir: Vect
     };
 
     // Calculate offset point
-    let d = dot(nor.abs(), pos_err);
+    let d = dot(nor.abs(), Vector::new(pos_err, pos_err, pos_err));
     let offset = nor * d;
     let p = pos + offset;
 
     // Calculate ulp offsets
-    //
-    // The additiona/subtraction of MIN_POSITIVE is to keep numbers out of the
-    // subnormal range when the original value is 0.0, because that seems to be
-    // causing issues.  Not sure why.  For now this seems like a reasonable
-    // solution because it only affects extremely small numbers near zero
-    // anyway.  But this seems worth investigating at some point.
-    //
-    // TODO: investigate cause of subnormal numbers being a problem, and fix
-    // if possible.  Test case: a horizontal plane at z = 0.0 and four lights
-    // evenly spaced apart at z = 4.0.  Causes strange render artifacts.
-    use std::f32::MIN_POSITIVE;
-
     let x = if nor.x() >= 0.0 {
-        increment_ulp(p.x() + MIN_POSITIVE)
+        increment_ulp(p.x())
     } else {
-        decrement_ulp(p.x() - MIN_POSITIVE)
+        decrement_ulp(p.x())
     };
 
     let y = if nor.y() >= 0.0 {
-        increment_ulp(p.y() + MIN_POSITIVE)
+        increment_ulp(p.y())
     } else {
-        decrement_ulp(p.y() - MIN_POSITIVE)
+        decrement_ulp(p.y())
     };
 
     let z = if nor.z() >= 0.0 {
-        increment_ulp(p.z() + MIN_POSITIVE)
+        increment_ulp(p.z())
     } else {
-        decrement_ulp(p.z() - MIN_POSITIVE)
+        decrement_ulp(p.z())
     };
 
     Point::new(x, y, z)

src/surface/mod.rs

@@ -37,7 +37,8 @@ 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 pos_err: f32, // Error magnitude of the intersection position.  Imagine
+    // a cube centered around `pos` with dimensions of `2 * pos_err`.
     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

@@ -8,7 +8,7 @@ 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 math::{Point, Vector, Matrix4x4, cross};
 use ray::{Ray, AccelRay};
 use shading::surface_closure::{SurfaceClosureUnion, GTRClosure, LambertClosure};
 
@@ -104,10 +104,10 @@ impl<'a> Surface for TriangleMesh<'a> {
                                         + (tri.1.into_vector() * b1)
                                         + (tri.2.into_vector() * b2)).into_point();
 
-                                    let pos_err = ((tri.0.into_vector().abs() * b0)
+                                    let pos_err = (((tri.0.into_vector().abs() * b0)
-                                        + (tri.1.into_vector().abs() * b1)
+                                            + (tri.1.into_vector().abs() * b1)
-                                        + (tri.2.into_vector().abs() * b2))
+                                            + (tri.2.into_vector().abs() * b2))
-                                        * fp_gamma(7);
+                                            * fp_gamma(7)).co.h_max();
 
                                     // Fill in intersection data
                                     isects[r.id as usize] = SurfaceIntersection::Hit {