diff --git a/src/fp_utils.rs b/src/fp_utils.rs new file mode 100644 index 0000000..3fa662a --- /dev/null +++ b/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::(&v) } +} + +#[inline(always)] +fn bits_to_f32(bits: u32) -> f32 { + use std::mem::transmute_copy; + unsafe { transmute_copy::(&bits) } +} diff --git a/src/main.rs b/src/main.rs index 864d16d..8e4defa 100644 --- a/src/main.rs +++ b/src/main.rs @@ -28,6 +28,7 @@ mod bbox; mod boundable; mod camera; mod color; +mod fp_utils; mod hash; mod hilbert; mod image; diff --git a/src/renderer.rs b/src/renderer.rs index 4e8671a..9259ca8 100644 --- a/src/renderer.rs +++ b/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 { diff --git a/src/surface/mod.rs b/src/surface/mod.rs index 7d80999..35b8b33 100644 --- a/src/surface/mod.rs +++ b/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 diff --git a/src/surface/triangle_mesh.rs b/src/surface/triangle_mesh.rs index 219693f..2618f0a 100644 --- a/src/surface/triangle_mesh.rs +++ b/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) diff --git a/sub_crates/math3d/src/vector.rs b/sub_crates/math3d/src/vector.rs index 9cbc856..2e52980 100644 --- a/sub_crates/math3d/src/vector.rs +++ b/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())