From 4e9bfd6e79031b20b714028959f874666c2cdf39 Mon Sep 17 00:00:00 2001
From: Nathan Vegdahl <cessen@cessen.com>
Date: Wed, 21 Jun 2017 01:36:08 -0700
Subject: [PATCH] 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 | 93 ++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 80 insertions(+), 13 deletions(-)

diff --git a/src/fp_utils.rs b/src/fp_utils.rs
index 3fa662a..4cf469a 100644
--- a/src/fp_utils.rs
+++ b/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 {
-        bits_to_f32(f32_to_bits(v) - 1)
-    } else {
+    if v <= -0.0 {
         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 {
-        increment_ulp(p.y())
+    let y = if nor.y() >= 0.0 {
+        increment_ulp(p.y() + MIN_POSITIVE)
     } else {
-        decrement_ulp(p.y())
+        decrement_ulp(p.y() - MIN_POSITIVE)
     };
 
-    let z = if offset.z() >= 0.0 {
-        increment_ulp(p.z())
+    let z = if nor.z() >= 0.0 {
+        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);
+    }
+
+}