diff --git a/src/fp_utils.rs b/src/fp_utils.rs
index 4cf469a..c2ee8f9 100644
--- a/src/fp_utils.rs
+++ b/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)
diff --git a/src/surface/mod.rs b/src/surface/mod.rs
index 35b8b33..5fee6f2 100644
--- a/src/surface/mod.rs
+++ b/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
diff --git a/src/surface/triangle_mesh.rs b/src/surface/triangle_mesh.rs
index 168cc40..aae8354 100644
--- a/src/surface/triangle_mesh.rs
+++ b/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)
-                                        + (tri.1.into_vector().abs() * b1)
-                                        + (tri.2.into_vector().abs() * b2))
-                                        * fp_gamma(7);
+                                    let pos_err = (((tri.0.into_vector().abs() * b0)
+                                            + (tri.1.into_vector().abs() * b1)
+                                            + (tri.2.into_vector().abs() * b2))
+                                            * fp_gamma(7)).co.h_max();
 
                                     // Fill in intersection data
                                     isects[r.id as usize] = SurfaceIntersection::Hit {