diff --git a/src/light/rectangle_light.rs b/src/light/rectangle_light.rs
index db01072..e399d68 100644
--- a/src/light/rectangle_light.rs
+++ b/src/light/rectangle_light.rs
@@ -284,8 +284,9 @@ impl<'a> Surface for RectangleLight<'a> {
             let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv;
 
             // Test against two triangles that make up the light
+            let ray_pre = triangle::RayTriPrecompute::new(dir);
             for tri in &[(p1, p2, p3), (p3, p4, p1)] {
-                if let Some((t, b0, b1, b2)) = triangle::intersect_ray(orig, dir, max_t, *tri) {
+                if let Some((t, b0, b1, b2)) = triangle::intersect_ray(orig, ray_pre, max_t, *tri) {
                     if t < max_t {
                         if rays.is_occlusion(ray_idx) {
                             isects[ray_idx] = SurfaceIntersection::Occlude;
diff --git a/src/surface/triangle.rs b/src/surface/triangle.rs
index 5f0a9f6..4aed3a3 100644
--- a/src/surface/triangle.rs
+++ b/src/surface/triangle.rs
@@ -5,6 +5,45 @@ use crate::{
     math::{Point, Vector},
 };
 
+#[derive(Debug, Copy, Clone)]
+pub struct RayTriPrecompute {
+    i: (usize, usize, usize),
+    s: (f32, f32, f32),
+}
+
+impl RayTriPrecompute {
+    pub fn new(ray_dir: Vector) -> RayTriPrecompute {
+        // Calculate the permuted dimension indices for the new ray space.
+        let (xi, yi, zi) = {
+            let xabs = ray_dir.x().abs();
+            let yabs = ray_dir.y().abs();
+            let zabs = ray_dir.z().abs();
+
+            if xabs > yabs && xabs > zabs {
+                (1, 2, 0)
+            } else if yabs > zabs {
+                (2, 0, 1)
+            } else {
+                (0, 1, 2)
+            }
+        };
+
+        let dir_x = ray_dir.get_n(xi);
+        let dir_y = ray_dir.get_n(yi);
+        let dir_z = ray_dir.get_n(zi);
+
+        // Calculate shear constants.
+        let sx = dir_x / dir_z;
+        let sy = dir_y / dir_z;
+        let sz = 1.0 / dir_z;
+
+        RayTriPrecompute {
+            i: (xi, yi, zi),
+            s: (sx, sy, sz),
+        }
+    }
+}
+
 /// Intersects `ray` with `tri`, returning `Some((t, b0, b1, b2))`, or `None`
 /// if no intersection.
 ///
@@ -18,45 +57,21 @@ use crate::{
 /// Intersection" by Woop et al.
 pub fn intersect_ray(
     ray_orig: Point,
-    ray_dir: Vector,
+    ray_pre: RayTriPrecompute,
     ray_max_t: f32,
     tri: (Point, Point, Point),
 ) -> Option<(f32, f32, f32, f32)> {
-    // Calculate the permuted dimension indices for the new ray space.
-    let (xi, yi, zi) = {
-        let xabs = ray_dir.x().abs();
-        let yabs = ray_dir.y().abs();
-        let zabs = ray_dir.z().abs();
-
-        if xabs > yabs && xabs > zabs {
-            (1, 2, 0)
-        } else if yabs > zabs {
-            (2, 0, 1)
-        } else {
-            (0, 1, 2)
-        }
-    };
-
-    let dir_x = ray_dir.get_n(xi);
-    let dir_y = ray_dir.get_n(yi);
-    let dir_z = ray_dir.get_n(zi);
-
-    // Calculate shear constants.
-    let sx = dir_x / dir_z;
-    let sy = dir_y / dir_z;
-    let sz = 1.0 / dir_z;
-
     // Calculate vertices in ray space.
     let p0 = tri.0 - ray_orig;
     let p1 = tri.1 - ray_orig;
     let p2 = tri.2 - ray_orig;
 
-    let p0x = p0.get_n(xi) - (sx * p0.get_n(zi));
-    let p0y = p0.get_n(yi) - (sy * p0.get_n(zi));
-    let p1x = p1.get_n(xi) - (sx * p1.get_n(zi));
-    let p1y = p1.get_n(yi) - (sy * p1.get_n(zi));
-    let p2x = p2.get_n(xi) - (sx * p2.get_n(zi));
-    let p2y = p2.get_n(yi) - (sy * p2.get_n(zi));
+    let p0x = p0.get_n(ray_pre.i.0) - (ray_pre.s.0 * p0.get_n(ray_pre.i.2));
+    let p0y = p0.get_n(ray_pre.i.1) - (ray_pre.s.1 * p0.get_n(ray_pre.i.2));
+    let p1x = p1.get_n(ray_pre.i.0) - (ray_pre.s.0 * p1.get_n(ray_pre.i.2));
+    let p1y = p1.get_n(ray_pre.i.1) - (ray_pre.s.1 * p1.get_n(ray_pre.i.2));
+    let p2x = p2.get_n(ray_pre.i.0) - (ray_pre.s.0 * p2.get_n(ray_pre.i.2));
+    let p2y = p2.get_n(ray_pre.i.1) - (ray_pre.s.1 * p2.get_n(ray_pre.i.2));
 
     // Calculate scaled barycentric coordinates.
     let mut e0 = (p1x * p2y) - (p1y * p2x);
@@ -82,9 +97,9 @@ pub fn intersect_ray(
     }
 
     // Calculate t of hitpoint.
-    let p0z = sz * p0.get_n(zi);
-    let p1z = sz * p1.get_n(zi);
-    let p2z = sz * p2.get_n(zi);
+    let p0z = ray_pre.s.2 * p0.get_n(ray_pre.i.2);
+    let p1z = ray_pre.s.2 * p1.get_n(ray_pre.i.2);
+    let p2z = ray_pre.s.2 * p2.get_n(ray_pre.i.2);
     let t_scaled = (e0 * p0z) + (e1 * p1z) + (e2 * p2z);
 
     // Check if the hitpoint t is within ray min/max t.
diff --git a/src/surface/triangle_mesh.rs b/src/surface/triangle_mesh.rs
index 967f90c..43388a8 100644
--- a/src/surface/triangle_mesh.rs
+++ b/src/surface/triangle_mesh.rs
@@ -183,6 +183,7 @@ impl<'a> Surface for TriangleMesh<'a> {
                     let mut hit_tri = unsafe { std::mem::uninitialized() };
                     let mut hit_tri_indices = unsafe { std::mem::uninitialized() };
                     let mut hit_tri_data = unsafe { std::mem::uninitialized() };
+                    let ray_pre = triangle::RayTriPrecompute::new(rays.dir(ray_idx));
                     for tri_idx in idx_range.clone() {
                         let tri_indices = self.indices[tri_idx];
 
@@ -229,7 +230,7 @@ impl<'a> Surface for TriangleMesh<'a> {
                         // Test ray against triangle
                         if let Some((t, b0, b1, b2)) = triangle::intersect_ray(
                             rays.orig(ray_idx),
-                            rays.dir(ray_idx),
+                            ray_pre,
                             rays.max_t(ray_idx),
                             tri,
                         ) {