Fix for ray origin offsets when intersection point is at 0.0.

For some reason the ulp incrementing is unreliable when starting
at zero.  It creates subnormal numbers, and that seems to be an
issue somewhere in the pipeline, ultimately leading to weird
render artifacts.  Not entirely sure why.

This fixes it by avoiding subnormal numbers in the final offset
ray origin.  Left a note suggesting investigating in more detail
at some point.

src/fp_utils.rs

@@ -18,6 +18,8 @@ pub fn increment_ulp(v: f32) -> f32 {
     if (v.is_infinite() && v > 0.0) || v.is_nan() {
         return v;
     }
+
+    // Handle zero
     let v = if v == -0.0 { 0.0 } else { v };
 
     // Increase ulp by 1
@@ -34,13 +36,15 @@ pub fn decrement_ulp(v: f32) -> f32 {
     if (v.is_infinite() && v < 0.0) || v.is_nan() {
         return v;
     }
-    let v = if v == -0.0 { 0.0 } else { v };
+
+    // Handle zero
+    let v = if v == 0.0 { -0.0 } else { v };
 
     // Decrease ulp by 1
-    if v >= 0.0 {
+    if v <= -0.0 {
-        bits_to_f32(f32_to_bits(v) - 1)
-    } else {
         bits_to_f32(f32_to_bits(v) + 1)
+    } else {
+        bits_to_f32(f32_to_bits(v) - 1)
     }
 }
 
@@ -58,22 +62,34 @@ pub fn robust_ray_origin(pos: Point, pos_err: Vector, nor: Normal, ray_dir: Vect
     let p = pos + offset;
 
     // Calculate ulp offsets
-    let x = if offset.x() >= 0.0 {
+    //
-        increment_ulp(p.x())
+    // 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)
     } else {
-        decrement_ulp(p.x())
+        decrement_ulp(p.x() - MIN_POSITIVE)
     };
 
-    let y = if offset.y() >= 0.0 {
+    let y = if nor.y() >= 0.0 {
-        increment_ulp(p.y())
+        increment_ulp(p.y() + MIN_POSITIVE)
     } else {
-        decrement_ulp(p.y())
+        decrement_ulp(p.y() - MIN_POSITIVE)
     };
 
-    let z = if offset.z() >= 0.0 {
+    let z = if nor.z() >= 0.0 {
-        increment_ulp(p.z())
+        increment_ulp(p.z() + MIN_POSITIVE)
     } else {
-        decrement_ulp(p.z())
+        decrement_ulp(p.z() - MIN_POSITIVE)
     };
 
     Point::new(x, y, z)
@@ -91,3 +107,54 @@ fn bits_to_f32(bits: u32) -> f32 {
     use std::mem::transmute_copy;
     unsafe { transmute_copy::<u32, f32>(&bits) }
 }
+
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn inc_ulp() {
+        assert!(increment_ulp(1.0) > 1.0);
+        assert!(increment_ulp(-1.0) > -1.0);
+    }
+
+    #[test]
+    fn dec_ulp() {
+        assert!(decrement_ulp(1.0) < 1.0);
+        assert!(decrement_ulp(-1.0) < -1.0);
+    }
+
+    #[test]
+    fn inc_ulp_zero() {
+        assert!(increment_ulp(0.0) > 0.0);
+        assert!(increment_ulp(0.0) > -0.0);
+        assert!(increment_ulp(-0.0) > 0.0);
+        assert!(increment_ulp(-0.0) > -0.0);
+    }
+
+    #[test]
+    fn dec_ulp_zero() {
+        assert!(decrement_ulp(0.0) < 0.0);
+        assert!(decrement_ulp(0.0) < -0.0);
+        assert!(decrement_ulp(-0.0) < 0.0);
+        assert!(decrement_ulp(-0.0) < -0.0);
+    }
+
+    #[test]
+    fn inc_dec_ulp() {
+        assert_eq!(decrement_ulp(increment_ulp(1.0)), 1.0);
+        assert_eq!(decrement_ulp(increment_ulp(-1.0)), -1.0);
+        assert_eq!(decrement_ulp(increment_ulp(1.2)), 1.2);
+        assert_eq!(decrement_ulp(increment_ulp(-1.2)), -1.2);
+    }
+
+    #[test]
+    fn dec_inc_ulp() {
+        assert_eq!(increment_ulp(decrement_ulp(1.0)), 1.0);
+        assert_eq!(increment_ulp(decrement_ulp(-1.0)), -1.0);
+        assert_eq!(increment_ulp(decrement_ulp(1.2)), 1.2);
+        assert_eq!(increment_ulp(decrement_ulp(-1.2)), -1.2);
+    }
+
+}