diff --git a/shared/src/cameras/orthographic.rs b/shared/src/cameras/orthographic.rs
index c60b8fd..8a9ab3a 100644
--- a/shared/src/cameras/orthographic.rs
+++ b/shared/src/cameras/orthographic.rs
@@ -96,7 +96,7 @@ impl CameraTrait for OrthographicCamera {
             p_camera,
             Vector3f::new(0., 0., 1.),
             Some(self.sample_time(sample.time)),
-            &self.base().medium,
+            self.base().medium,
         );
         if self.lens_radius > 0. {
             let p_lens_vec =
diff --git a/shared/src/cameras/realistic.rs b/shared/src/cameras/realistic.rs
index 6cdf29a..8b738e7 100644
--- a/shared/src/cameras/realistic.rs
+++ b/shared/src/cameras/realistic.rs
@@ -132,7 +132,7 @@ impl RealisticCamera {
             Point3f::new(0., x, self.lens_front_z() + 1.),
             Vector3f::new(0., 0., -1.),
             None,
-            &Ptr::null(),
+            Ptr::null(),
         );
         let Some((_, r_film)) = self.trace_lenses_from_film(&r_scene) else {
             panic!(
@@ -144,7 +144,7 @@ impl RealisticCamera {
             Point3f::new(x, 0., self.lens_rear_z() - 1.),
             Vector3f::new(0., 0., 1.),
             None,
-            &Ptr::null(),
+            Ptr::null(),
         );
         let Some((_, r_scene)) = self.trace_lenses_from_film(&r_film) else {
             panic!(
@@ -209,7 +209,7 @@ impl RealisticCamera {
             // Expand pupil bounds if ray makes it through the lens system
             if !pupil_bounds.contains(Point2f::new(p_rear.x(), p_rear.y()))
                 && trace_lenses_from_film(
-                    Ray::new(p_film, p_rear - p_film, None, &Ptr::null()),
+                    Ray::new(p_film, p_rear - p_film, None, Ptr::null()),
                     None,
                 )
             {
@@ -270,7 +270,7 @@ impl RealisticCamera {
             Point3f::new(r_camera.o.x(), r_camera.o.y(), -r_camera.o.z()),
             Vector3f::new(r_camera.d.x(), r_camera.d.y(), -r_camera.d.z()),
             Some(r_camera.time),
-            &Ptr::null(),
+            Ptr::null(),
         );
 
         for i in (0..self.n_elements - 1).rev() {
@@ -337,7 +337,7 @@ impl RealisticCamera {
             Point3f::new(r_lens.o.x(), r_lens.o.y(), -r_lens.o.z()),
             Vector3f::new(r_lens.d.x(), r_lens.d.y(), -r_lens.d.z()),
             Some(r_lens.time),
-            &Ptr::null(),
+            Ptr::null(),
         );
 
         Some((weight, r_out))
@@ -415,7 +415,7 @@ impl CameraTrait for RealisticCamera {
         let eps = self.sample_exit_pupil(Point2f::new(p_film.x(), p_film.y()), sample.p_lens)?;
 
         let p_pupil = Point3f::new(0., 0., 0.);
-        let r_film = Ray::new(p_film, p_pupil - p_film, None, &Ptr::null());
+        let r_film = Ray::new(p_film, p_pupil - p_film, None, Ptr::null());
         let (weight, mut ray) = self.trace_lenses_from_film(&r_film)?;
         if weight == 0. {
             return None;
diff --git a/shared/src/cameras/spherical.rs b/shared/src/cameras/spherical.rs
index 467f0b8..6caaf11 100644
--- a/shared/src/cameras/spherical.rs
+++ b/shared/src/cameras/spherical.rs
@@ -54,7 +54,7 @@ impl CameraTrait for SphericalCamera {
             Point3f::new(0., 0., 0.),
             dir,
             Some(self.sample_time(sample.time)),
-            &self.base().medium,
+            self.base().medium,
         );
         Some(CameraRay {
             ray: self.render_from_camera(&ray, &mut None),
diff --git a/shared/src/core/camera.rs b/shared/src/core/camera.rs
index b192f7d..bfbe04d 100644
--- a/shared/src/core/camera.rs
+++ b/shared/src/core/camera.rs
@@ -227,13 +227,13 @@ pub trait CameraTrait {
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_x,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_x,
             None,
-            &Ptr::default(),
+            Ptr::default(),
         );
         let y_ray = Ray::new(
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_y,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_y,
             None,
-            &Ptr::default(),
+            Ptr::default(),
         );
         let n_down = Vector3f::from(n_down_z);
         let tx = -(n_down.dot(y_ray.o.into())) / n_down.dot(x_ray.d);
diff --git a/shared/src/core/geometry/ray.rs b/shared/src/core/geometry/ray.rs
index 1daf9c0..1851ea5 100644
--- a/shared/src/core/geometry/ray.rs
+++ b/shared/src/core/geometry/ray.rs
@@ -30,12 +30,12 @@ impl Default for Ray {
 }
 
 impl Ray {
-    pub fn new(o: Point3f, d: Vector3f, time: Option<Float>, medium: &Medium) -> Self {
+    pub fn new(o: Point3f, d: Vector3f, time: Option<Float>, medium: Ptr<Medium>) -> Self {
         Self {
             o,
             d,
             time: time.unwrap_or_else(|| Self::default().time),
-            medium: Ptr::from(medium),
+            medium,
             ..Self::default()
         }
     }
diff --git a/shared/src/core/shape.rs b/shared/src/core/shape.rs
index af160f0..228529f 100644
--- a/shared/src/core/shape.rs
+++ b/shared/src/core/shape.rs
@@ -118,7 +118,7 @@ impl ShapeSampleContext {
     }
 
     pub fn spawn_ray(&self, w: Vector3f) -> Ray {
-        Ray::new(self.offset_ray_origin(w), w, Some(self.time), &Ptr::null())
+        Ray::new(self.offset_ray_origin(w), w, Some(self.time), Ptr::null())
     }
 }
 
diff --git a/shared/src/shapes/sphere.rs b/shared/src/shapes/sphere.rs
index 0d11077..3918196 100644
--- a/shared/src/shapes/sphere.rs
+++ b/shared/src/shapes/sphere.rs
@@ -57,7 +57,7 @@ impl SphereShape {
         phi_max: Float,
     ) -> Self {
         let theta_z_min = clamp(z_min.min(z_max) / radius, -1., 1.).acos();
-        let theta_z_max = clamp(z_max.min(z_max) / radius, -1., 1.).acos();
+        let theta_z_max = clamp(z_min.max(z_max) / radius, -1., 1.).acos();
         let phi_max = radians(clamp(phi_max, 0., 360.0));
         Self {
             render_from_object: render_from_object.clone(),
diff --git a/shared/src/shapes/triangle.rs b/shared/src/shapes/triangle.rs
index 8ba2d97..6e795de 100644
--- a/shared/src/shapes/triangle.rs
+++ b/shared/src/shapes/triangle.rs
@@ -78,7 +78,6 @@ impl TriangleShape {
         Some([mesh.s[v0], mesh.s[v1], mesh.s[v2]])
     }
 
-
     fn get_uvs(&self) -> Option<[Point2f; 3]> {
         let mesh = self.mesh();
         if mesh.s.is_empty() {
@@ -107,13 +106,125 @@ impl TriangleShape {
 
     fn intersect_triangle(
         &self,
-        _ray: &Ray,
-        _t_max: Float,
-        _p0: Point3f,
-        _p1: Point3f,
-        _p2: Point3f,
+        ray: &Ray,
+        t_max: Float,
+        p0: Point3f,
+        p1: Point3f,
+        p2: Point3f,
     ) -> Option<TriangleIntersection> {
-        todo!()
+        if (p2 - p0).cross(p1 - p0).norm_squared() == 0. {
+            return None;
+        }
+
+        // Transform triangle vertices to ray coordinate space
+        // Translate vertices based on ray origin
+        let mut p0t = p0 - Vector3f::from(ray.o);
+        let mut p1t = p1 - Vector3f::from(ray.o);
+        let mut p2t = p2 - Vector3f::from(ray.o);
+
+        // Permute components of triangle vertices and ray direction
+        let kz = ray.d.abs().max_component_index();
+        let mut kx = kz + 1;
+        if kx == 3 {
+            kx = 0;
+        }
+        let mut ky = kx + 1;
+        if ky == 3 {
+            ky = 0;
+        }
+        let d = ray.d.permute([kx, ky, kz]);
+        p0t = p0t.permute([kx, ky, kz]);
+        p1t = p1t.permute([kx, ky, kz]);
+        p2t = p2t.permute([kx, ky, kz]);
+
+        // Apply shear transformation to translated vertex positions
+        let sx = -d.x() / d.z();
+        let sy = -d.y() / d.z();
+        let sz = 1. / d.z();
+        p0t[0] += sx * p0t[2];
+        p0t[1] += sy * p0t[2];
+        p1t[0] += sx * p1t[2];
+        p1t[1] += sy * p1t[2];
+        p2t[0] += sx * p2t[2];
+        p2t[1] += sy * p2t[2];
+
+        // Compute edge function coefficients _e0_, _e1_, and _e2_
+        let e0 = difference_of_products(p1t.x(), p2t.y(), p1t.y(), p2t.x());
+        let e1 = difference_of_products(p2t.x(), p0t.y(), p2t.y(), p0t.x());
+        let e2 = difference_of_products(p0t.x(), p1t.y(), p0t.y(), p1t.x());
+
+        // Fall back to double-precision test at triangle edges
+        // if sizeof(Float) == sizeof(float) && (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f)) {
+        //     double p2txp1ty = (double)p2t.x * (double)p1t.y;
+        //     double p2typ1tx = (double)p2t.y * (double)p1t.x;
+        //     e0 = (float)(p2typ1tx - p2txp1ty);
+        //     double p0txp2ty = (double)p0t.x * (double)p2t.y;
+        //     double p0typ2tx = (double)p0t.y * (double)p2t.x;
+        //     e1 = (float)(p0typ2tx - p0txp2ty);
+        //     double p1txp0ty = (double)p1t.x * (double)p0t.y;
+        //     double p1typ0tx = (double)p1t.y * (double)p0t.x;
+        //     e2 = (float)(p1typ0tx - p1txp0ty);
+        // }
+
+        // Perform triangle edge and determinant tests
+        if (e0 < 0. || e1 < 0. || e2 < 0.) && (e0 > 0. || e1 > 0. || e2 > 0.) {
+            return None;
+        }
+        let det = e0 + e1 + e2;
+        if det == 0. {
+            return None;
+        }
+
+        // Compute scaled hit distance to triangle and test against ray $t$ range
+        p0t[2] *= sz;
+        p1t[2] *= sz;
+        p2t[2] *= sz;
+        let t_scaled = e0 * p0t.z() + e1 * p1t.z() + e2 * p2t.z();
+        if det < 0. && (t_scaled >= 0. || t_scaled < t_max * det) {
+            return None;
+        } else if det > 0. && (t_scaled <= 0. || t_scaled > t_max * det) {
+            return None;
+        }
+
+        // Compute barycentric coordinates and $t$ value for triangle intersection
+        let inv_det = 1. / det;
+        let b0 = e0 * inv_det;
+        let b1 = e1 * inv_det;
+        let b2 = e2 * inv_det;
+        let t = t_scaled * inv_det;
+
+        debug_assert!(t.is_finite());
+
+        // Ensure that computed triangle $t$ is conservatively greater than zero
+        // Compute $\delta_z$ term for triangle $t$ error bounds
+        let maxZt = Vector3f::new(p0t.z(), p1t.z(), p2t.z())
+            .abs()
+            .max_component_value();
+        let deltaZ = gamma(3) * maxZt;
+
+        // Compute $\delta_x$ and $\delta_y$ terms for triangle $t$ error bounds
+        let maxXt = Vector3f::new(p0t.x(), p1t.x(), p2t.x())
+            .abs()
+            .max_component_value();
+        let maxYt = Vector3f::new(p0t.y(), p1t.y(), p2t.y())
+            .abs()
+            .max_component_value();
+        let deltaX = gamma(5) * (maxXt + maxZt);
+        let deltaY = gamma(5) * (maxYt + maxZt);
+
+        // Compute $\delta_e$ term for triangle $t$ error bounds
+        let deltaE = 2. * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
+
+        // Compute $\delta_t$ term for triangle $t$ error bounds and check _t_
+        let maxE = Vector3f::new(e0, e1, e2).abs().max_component_value();
+        let deltaT =
+            3. * (gamma(3) * maxE * maxZt + deltaE * maxZt + deltaZ * maxE) * inv_det.abs();
+        if t <= deltaT {
+            return None;
+        }
+
+        // Return _TriangleIntersection_ for intersection
+        return Some(TriangleIntersection { b0, b1, b2, t });
     }
 
     fn interaction_from_intersection(
diff --git a/shared/src/utils/math.rs b/shared/src/utils/math.rs
index c86045a..400cbc8 100644
--- a/shared/src/utils/math.rs
+++ b/shared/src/utils/math.rs
@@ -113,7 +113,7 @@ pub fn safe_asin<T: NumFloat>(x: T) -> T {
 #[inline]
 pub fn safe_acos(x: Float) -> Float {
     if (-1.001..1.001).contains(&x) {
-        clamp(x, -1., 1.).asin()
+        clamp(x, -1., 1.).acos()
     } else {
         panic!("Not valid value for acos")
     }
@@ -542,9 +542,9 @@ pub fn next_float_up(v: Float) -> Float {
 
     let mut ui = float_to_bits(v);
     if v >= 0.0 {
-        ui += 1;
+        ui = ui.wrapping_add(1);
     } else {
-        ui -= 1;
+        ui = ui.wrapping_sub(1);
     }
     bits_to_float(ui)
 }
@@ -558,9 +558,9 @@ pub fn next_float_down(v: Float) -> Float {
 
     let mut ui = float_to_bits(v);
     if v > 0.0 {
-        ui -= 1;
+        ui = ui.wrapping_sub(1);
     } else {
-        ui += 1;
+        ui = ui.wrapping_add(1);
     }
     bits_to_float(ui)
 }
diff --git a/shared/src/utils/sampling.rs b/shared/src/utils/sampling.rs
index 315d691..aaa6200 100644
--- a/shared/src/utils/sampling.rs
+++ b/shared/src/utils/sampling.rs
@@ -774,6 +774,9 @@ impl PiecewiseConstant1D {
 
     pub fn find_interval(&self, u: Float) -> usize {
         let n = self.func.len();
+        if n == 0 {
+            return 0
+        }
         let mut size = n;
         let mut first = 0usize;
         while size > 0 {
diff --git a/shared/src/utils/transform.rs b/shared/src/utils/transform.rs
index c757d9d..dd86b9f 100644
--- a/shared/src/utils/transform.rs
+++ b/shared/src/utils/transform.rs
@@ -165,7 +165,7 @@ impl TransformGeneric<Float> {
                 *t -= dt;
             }
         }
-        Ray::new(o.into(), r.d, Some(r.time), &*r.medium)
+        Ray::new(o.into(), r.d, Some(r.time), r.medium)
     }
 
     pub fn apply_to_interval(&self, pi: &Point3fi) -> Point3fi {
@@ -366,7 +366,7 @@ impl TransformGeneric<Float> {
                 t = t_max - dt;
             }
         }
-        (Ray::new(Point3f::from(o), d, Some(r.time), &*r.medium), t)
+        (Ray::new(Point3f::from(o), d, Some(r.time), r.medium), t)
     }
 
     pub fn to_quaternion(self) -> Quaternion {