diff --git a/.gitignore b/.gitignore
index b6f4691..445d527 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 *.lock
 *.log
 *.bak
+flip.rs
diff --git a/Cargo.toml b/Cargo.toml
index 53cc07e..e61eccf 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -4,8 +4,18 @@ version = "0.1.0"
 edition = "2024"
 
 [dependencies]
+anyhow = "1.0.100"
 bitflags = "2.10.0"
 bumpalo = "3.19.0"
+exr = "1.73.0"
+half = "2.7.1"
+image = "0.25.8"
+num = "0.4.3"
+num-integer = "0.1.46"
 num-traits = "0.2.19"
 once_cell = "1.21.3"
+qoi = "0.4.1"
 rand = "0.9.2"
+rayon = "1.11.0"
+smallvec = "1.15.1"
+thiserror = "2.0.17"
diff --git a/src/camera/mod.rs b/src/camera/mod.rs
index 4ec8440..a8d5e61 100644
--- a/src/camera/mod.rs
+++ b/src/camera/mod.rs
@@ -219,6 +219,7 @@ pub enum Camera {
     Spherical(SphericalCamera),
     Realistic(RealisticCamera),
 }
+
 impl CameraTrait for Camera {
     fn base(&self) -> &CameraBase {
         match self {
diff --git a/src/core/interaction.rs b/src/core/interaction.rs
index 90691eb..8cdb239 100644
--- a/src/core/interaction.rs
+++ b/src/core/interaction.rs
@@ -4,7 +4,7 @@ use crate::core::medium::{Medium, MediumInterface};
 use crate::core::pbrt::Float;
 use crate::geometry::{Normal3f, Point2f, Point3f, Point3fi, Ray, Vector3f, VectorLike};
 use crate::lights::Light;
-use crate::shapes::ShapeTrait;
+use crate::shapes::Shape;
 
 use std::any::Any;
 use std::sync::Arc;
@@ -106,28 +106,14 @@ pub struct SurfaceInteraction {
     pub dvdx: Float,
     pub dudy: Float,
     pub dvdy: Float,
-    pub shape: Option<Arc<dyn ShapeTrait>>,
+    pub shape: Arc<Shape>,
     pub bsdf: Option<BSDF>,
 }
 
 impl SurfaceInteraction {
-    pub fn set_intersection_properties(
-        &mut self,
-        mtl: Arc<dyn MaterialTrait>,
-        area: Arc<Light>,
-        prim_medium_interface: Option<Arc<MediumInterface>>,
-        ray_medium: Arc<Medium>,
-    ) {
-        self.material = Some(mtl);
-        self.area_light = Some(area);
-        if prim_medium_interface.as_ref().map_or(false, |mi| mi.is_medium_transition()) {
-            self.common.medium_interface = prim_medium_interface;
-        } else {
-            self.common.medium = Some(ray_medium);
-        }
-    }
 }
 
+#[derive(Default, Debug, Clone)]
 pub struct PhaseFunction;
 
 pub struct MediumInteraction {
@@ -234,7 +220,7 @@ impl SurfaceInteraction {
             dudy: 0.0,
             dvdx: 0.0,
             dvdy: 0.0,
-            shape: None,
+            shape: Arc::new(Shape::default()),
             bsdf: None,
         }
     }
@@ -288,8 +274,30 @@ impl SurfaceInteraction {
         si.uv = uv;
         si
     }
+
+    pub fn set_intersection_properties(
+        &mut self,
+        mtl: Arc<dyn MaterialTrait>,
+        area: Arc<Light>,
+        prim_medium_interface: Option<Arc<MediumInterface>>,
+        ray_medium: Arc<Medium>,
+    ) {
+        self.material = Some(mtl);
+        self.area_light = Some(area);
+        if prim_medium_interface.as_ref().map_or(false, |mi| mi.is_medium_transition()) {
+            self.common.medium_interface = prim_medium_interface;
+        } else {
+            self.common.medium = Some(ray_medium);
+        }
+    }
+
 }
 
+// pub enum InteractionEnum {
+//     Surface(SurfaceInteraction),
+//     Medium(MediumInteraction),
+// }
+
 impl Interaction for MediumInteraction {
     fn get_common(&self) -> &InteractionData {
         &self.common
diff --git a/src/core/primitive.rs b/src/core/primitive.rs
index 82e7572..73f836c 100644
--- a/src/core/primitive.rs
+++ b/src/core/primitive.rs
@@ -1,12 +1,13 @@
-use crate::core::interaction::Interaction;
+use crate::core::interaction::{Interaction, SurfaceInteraction};
 use crate::core::pbrt::Float;
 use crate::geometry::{Bounds3f, Ray};
-use crate::shapes::{ShapeIntersection, ShapeTrait};
+use crate::shapes::{ShapeIntersection, Shape};
 use crate::core::material::MaterialTrait;
 use crate::core::medium::MediumInterface;
 use crate::lights::Light;
 use crate::core::texture::{FloatTextureTrait, TextureEvalContext};
 use crate::utils::hash::hash_float;
+use crate::utils::transform::{AnimatedTransform, Transform};
 
 use std::sync::Arc;
 
@@ -18,7 +19,7 @@ pub trait PrimitiveTrait: Send + Sync + std::fmt::Debug {
 
 #[derive(Debug, Clone)]
 pub struct GeometricPrimitive {
-    shape: Arc<dyn ShapeTrait>,
+    shape: Arc<Shape>,
     material: Arc<dyn MaterialTrait>,
     area_light: Arc<Light>,
     medium_interface: Arc<MediumInterface>,
@@ -45,7 +46,9 @@ impl PrimitiveTrait for GeometricPrimitive {
                 return Some(si_next)
             }
         }
-        si.intr_mut().set_intersection_properties(self.material.clone(), self.area_light.clone(), Some(self.medium_interface.clone()), r.medium.clone()?);
+        if let Some(si) = si.intr().downcast_mut::<SurfaceInteraction>() {
+            si.as_ref().set_intersection_properties(self.material.clone(), self.area_light.clone(), Some(self.medium_interface.clone()), r.medium.clone()?);
+        }
         Some(si)
     }
 
@@ -59,13 +62,56 @@ impl PrimitiveTrait for GeometricPrimitive {
 }
 
 #[derive(Debug, Clone)]
-pub struct SimplePrimitive {
-    shape: Arc<dyn ShapeTrait>,
+pub struct SimplePrimitiv {
+    shape: Arc<Shape>,
     material: Arc<dyn MaterialTrait>,
 }
 
-pub struct TransformedPrimitive;
-pub struct AnimatedPrimitive;
+#[derive(Debug, Clone)]
+pub struct TransformedPrimitive {
+    primitive: Arc<dyn PrimitiveTrait>,
+    render_from_primitive: Transform<Float>,
+}
+
+impl PrimitiveTrait for TransformedPrimitive {
+    fn bounds(&self) -> Bounds3f {
+        self.render_from_primitive.apply_to_bounds(self.primitive.bounds())
+    }
+
+    fn intersect(&self, _r: &Ray, _t_max: Option<Float>) -> Option<ShapeIntersection> {
+        todo!()
+    }
+
+    fn intersect_p(&self, _r: &Ray, _t_max: Option<Float>) -> bool {
+        todo!()
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct AnimatedPrimitive {
+    primitive: Arc<dyn PrimitiveTrait>,
+    render_from_primitive: AnimatedTransform
+}
+
+impl PrimitiveTrait for AnimatedPrimitive {
+    fn bounds(&self) -> Bounds3f {
+        self.render_from_primitive.motion_bounds(&self.primitive.bounds())
+    }
+
+    fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+        let interp_render_from_primitive = self.render_from_primitive.interpolate(r.time);
+        let ray = interp_render_from_primitive.apply_inverse_ray(r, t_max);
+        let mut si = self.primitive.intersect(r, t_max) else { return None };
+        let new_render = interp_render_from_primitive.apply_to_interaction(si?.intr());
+        let transform = new_render.as_any().downcast_mut::<SurfaceInteraction>().expect(
+        "Failed to downcast Interaction to SurfaceInteraction. This should not happen."
+        );
+
+        *si?.intr = new_render;
+        si
+    }
+
+}
 pub struct BVHAggregatePrimitive;
 pub struct KdTreeAggregate;
 
@@ -77,6 +123,7 @@ pub enum Primitive {
     KdTree(KdTreeAggregate),
 }
 
+
 impl Primitive {
     // pub fn bounds(&self) -> Bounds3f {
     //     match self {
diff --git a/src/geometry/primitives.rs b/src/geometry/primitives.rs
index 825f103..350c016 100644
--- a/src/geometry/primitives.rs
+++ b/src/geometry/primitives.rs
@@ -656,6 +656,16 @@ impl<const N: usize> From<Vector<Interval, N>> for Vector<Float, N> {
     }
 }
 
+impl<const N: usize> From<Vector<Float, N>> for Vector<Interval, N> {
+    fn from(v: Vector<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new(v[i]);
+        }
+        Self(arr)
+    }
+}
+
 impl<const N: usize> Mul<Vector<Interval, N>> for Interval {
     type Output = Vector<Interval, N>;
     fn mul(self, rhs: Vector<Interval, N>) -> Self::Output {
diff --git a/src/lights/diffuse.rs b/src/lights/diffuse.rs
index 2c2a4bc..ed0b706 100644
--- a/src/lights/diffuse.rs
+++ b/src/lights/diffuse.rs
@@ -1,14 +1,14 @@
 use crate::core::medium::MediumInterface;
 use crate::core::pbrt::Float;
-use crate::shapes::ShapeTrait;
+use crate::shapes::Shape;
 use crate::utils::spectrum::Spectrum;
 use crate::utils::transform::Transform;
 
 use std::sync::Arc;
 
-pub struct DiffuseAreaLight<'a> {
+pub struct DiffuseAreaLight {
     pub l_emit: Spectrum,
-    pub shape: Arc<&'a dyn ShapeTrait>,
+    pub shape: Arc<Shape>,
     pub two_sided: bool,
     pub area: Float,
     pub flags: u8,
diff --git a/src/shapes/bilinear.rs b/src/shapes/bilinear.rs
index badf952..9ad0566 100644
--- a/src/shapes/bilinear.rs
+++ b/src/shapes/bilinear.rs
@@ -1,7 +1,7 @@
 use super::{
     BilinearIntersection, BilinearPatchShape, Bounds3f, DirectionCone, Interaction, Normal3f,
     Point2f, Point3f, Point3fi, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
-    ShapeTrait, SurfaceInteraction, Vector3f,
+    SurfaceInteraction, Vector3f,
 };
 use crate::core::pbrt::{Float, clamp_t, gamma, lerp};
 use crate::geometry::{Tuple, VectorLike, spherical_quad_area};
@@ -500,14 +500,12 @@ impl BilinearPatchShape {
         }
         n
     }
-}
 
-impl ShapeTrait for BilinearPatchShape {
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         self.area
     }
 
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         let data = self.get_data();
         if data.p00 == data.p10
             || data.p10 == data.p11
@@ -554,19 +552,19 @@ impl ShapeTrait for BilinearPatchShape {
         DirectionCone::new(n_avg.into(), clamp_t(cos_theta, -1., 1.))
     }
 
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         let data = self.get_data();
         Bounds3f::from_points(data.p00, data.p01).union(Bounds3f::from_points(data.p10, data.p11))
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let t_max_val = t_max?;
         let data = self.get_data();
         if let Some(bilinear_hit) = self.intersect_bilinear_patch(ray, t_max_val, &data) {
             let intr = self.interaction_from_intersection(&data, bilinear_hit.uv, ray.time, -ray.d);
 
             Some(ShapeIntersection {
-                intr,
+                intr: Box::new(intr),
                 t_hit: bilinear_hit.t,
             })
         } else {
@@ -574,14 +572,14 @@ impl ShapeTrait for BilinearPatchShape {
         }
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         let t_max_val = t_max.unwrap_or(Float::INFINITY);
         let data = self.get_data();
         self.intersect_bilinear_patch(ray, t_max_val, &data)
             .is_some()
     }
 
-    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
         let data = self.get_data();
         // Sample bilinear patch parametric coordinate (u, v)
         let (uv, pdf) = self.sample_parametric_coords(&data, u);
@@ -613,7 +611,7 @@ impl ShapeTrait for BilinearPatchShape {
         })
     }
 
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
         let data = self.get_data();
         let v00 = (data.p00 - ctx.p()).normalize();
         let v10 = (data.p10 - ctx.p()).normalize();
@@ -630,7 +628,7 @@ impl ShapeTrait for BilinearPatchShape {
         }
     }
 
-    fn pdf(&self, intr: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, intr: &dyn Interaction) -> Float {
         let Some(si) = intr.as_any().downcast_ref::<SurfaceInteraction>() else {
             return 0.;
         };
@@ -660,7 +658,7 @@ impl ShapeTrait for BilinearPatchShape {
         if cross == 0. { 0. } else { param_pdf / cross }
     }
 
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let ray = ctx.spawn_ray(wi);
         let Some(isect) = self.intersect(&ray, None) else {
             return 0.;
@@ -677,7 +675,7 @@ impl ShapeTrait for BilinearPatchShape {
             || data.mesh.image_distribution.is_some()
             || spherical_quad_area(v00, v10, v01, v11) <= Self::MIN_SPHERICAL_SAMPLE_AREA;
         if use_area_sampling {
-            let isect_pdf = self.pdf(Arc::new(&isect.intr));
+            let isect_pdf = self.pdf(isect.intr.as_ref());
             let distsq = ctx.p().distance_squared(isect.intr.p());
             let absdot = Vector3f::from(isect.intr.n()).abs_dot(-wi);
             if absdot == 0. {
diff --git a/src/shapes/curves.rs b/src/shapes/curves.rs
index 7dc517b..b3a21a0 100644
--- a/src/shapes/curves.rs
+++ b/src/shapes/curves.rs
@@ -8,18 +8,18 @@ use crate::utils::transform::look_at;
 use super::{
     Bounds3f, CurveCommon, CurveShape, CurveType, DirectionCone, Float, Interaction, Normal3f,
     Point2f, Point3f, Point3fi, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
-    ShapeTrait, SurfaceInteraction, Transform, Vector2f, Vector3f, VectorLike,
+    SurfaceInteraction, Transform, Vector2f, Vector3f, VectorLike,
 };
 use std::sync::Arc;
 
-struct IntersectionContext<'a> {
-    ray: &'a Ray,
-    object_from_ray: &'a Transform<Float>,
-    common: &'a CurveCommon<'a>,
+struct IntersectionContext {
+    ray: Ray,
+    object_from_ray: Arc<Transform<Float>>,
+    common: CurveCommon,
 }
 
-impl<'a> CurveShape<'a> {
-    pub fn new(common: CurveCommon<'a>, u_min: Float, u_max: Float) -> Self {
+impl CurveShape {
+    pub fn new(common: CurveCommon, u_min: Float, u_max: Float) -> Self {
         Self {
             common,
             u_min,
@@ -76,9 +76,9 @@ impl<'a> CurveShape<'a> {
         };
 
         let context = IntersectionContext {
-            ray: &ray,
-            object_from_ray: &ray_from_object.inverse(),
-            common: &self.common,
+            ray: ray,
+            object_from_ray: Arc::new(ray_from_object.inverse()),
+            common: self.common.clone(),
         };
 
         self.recursive_intersect(&context, t_max, &cp, self.u_min, self.u_max, max_depth)
@@ -257,12 +257,10 @@ impl<'a> CurveShape<'a> {
             flip_normal,
         );
 
-        Some(ShapeIntersection { intr, t_hit })
+        Some(ShapeIntersection { intr: Box::new(intr), t_hit })
     }
-}
 
-impl ShapeTrait for CurveShape<'_> {
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         let cs_span = self.common.cp_obj;
         let obj_bounds = bound_cubic_bezier(&cs_span, self.u_min, self.u_max);
         let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]);
@@ -273,11 +271,11 @@ impl ShapeTrait for CurveShape<'_> {
             .apply_to_bounds(obj_bounds_expand)
     }
 
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         DirectionCone::entire_sphere()
     }
 
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         let cp_obj = cubic_bezier_control_points(&self.common.cp_obj, self.u_min, self.u_max);
         let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]);
         let width1 = lerp(self.u_max, self.common.width[0], self.common.width[1]);
@@ -289,28 +287,28 @@ impl ShapeTrait for CurveShape<'_> {
         approx_length * avg_width
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY))
             .is_some()
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY))
     }
 
-    fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
         todo!()
     }
 
-    fn pdf_from_context(&self, _ctx: &ShapeSampleContext, _wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, _ctx: &ShapeSampleContext, _wi: Vector3f) -> Float {
         todo!()
     }
 
-    fn sample(&self, _u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, _u: Point2f) -> Option<ShapeSample> {
         todo!()
     }
 
-    fn sample_from_context(&self, _ctx: &ShapeSampleContext, _u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, _ctx: &ShapeSampleContext, _u: Point2f) -> Option<ShapeSample> {
         todo!()
     }
 }
diff --git a/src/shapes/cylinder.rs b/src/shapes/cylinder.rs
index 4c8fb96..8ff78bb 100644
--- a/src/shapes/cylinder.rs
+++ b/src/shapes/cylinder.rs
@@ -1,7 +1,7 @@
 use super::{
     Bounds3f, CylinderShape, DirectionCone, Float, Interaction, Normal3f, PI, Point2f, Point3f,
     Point3fi, QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
-    ShapeTrait, SurfaceInteraction, Transform, Vector3f, Vector3fi,
+    SurfaceInteraction, Transform, Vector3f, Vector3fi,
 };
 use crate::core::pbrt::{gamma, lerp};
 use crate::geometry::{Sqrt, Tuple, VectorLike};
@@ -10,10 +10,10 @@ use crate::utils::math::{difference_of_products, square};
 use std::mem;
 use std::sync::Arc;
 
-impl<'a> CylinderShape<'a> {
+impl CylinderShape {
     pub fn new(
-        render_from_object: &'a Transform<Float>,
-        object_from_render: &'a Transform<Float>,
+        render_from_object: Arc<Transform<Float>>,
+        object_from_render:  Arc<Transform<Float>>,
         reverse_orientation: bool,
         radius: Float,
         z_min: Float,
@@ -25,7 +25,7 @@ impl<'a> CylinderShape<'a> {
             z_min,
             z_max,
             phi_max,
-            render_from_object,
+            render_from_object: render_from_object.clone(),
             object_from_render,
             reverse_orientation,
             transform_swap_handedness: render_from_object.swaps_handedness(),
@@ -168,14 +168,12 @@ impl<'a> CylinderShape<'a> {
         );
         surf_point
     }
-}
 
-impl ShapeTrait for CylinderShape<'_> {
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         (self.z_max - self.z_min) * self.radius * self.phi_max
     }
 
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         self.render_from_object
             .apply_to_bounds(Bounds3f::from_points(
                 Point3f::new(-self.radius, -self.radius, self.z_min),
@@ -183,11 +181,11 @@ impl ShapeTrait for CylinderShape<'_> {
             ))
     }
 
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         DirectionCone::entire_sphere()
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
@@ -197,7 +195,7 @@ impl ShapeTrait for CylinderShape<'_> {
         }
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         if let Some(t) = t_max {
             self.basic_intersect(ray, t).is_some()
         } else {
@@ -205,11 +203,11 @@ impl ShapeTrait for CylinderShape<'_> {
         }
     }
 
-    fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
         1. / self.area()
     }
 
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let ray = ctx.spawn_ray(wi);
         if let Some(isect) = self.intersect(&ray, None) {
             let n = isect.intr.n();
@@ -224,7 +222,7 @@ impl ShapeTrait for CylinderShape<'_> {
         }
     }
 
-    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
         let z = lerp(u[0], self.z_min, self.z_max);
         let phi = u[1] * self.phi_max;
         let mut p_obj = Point3f::new(self.radius * phi.cos(), self.radius * phi.sin(), z);
@@ -253,7 +251,7 @@ impl ShapeTrait for CylinderShape<'_> {
         })
     }
 
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
         let mut ss = self.sample(u)?;
         let intr = Arc::make_mut(&mut ss.intr);
         intr.get_common_mut().time = ctx.time;
diff --git a/src/shapes/disk.rs b/src/shapes/disk.rs
index fe3e99f..f435419 100644
--- a/src/shapes/disk.rs
+++ b/src/shapes/disk.rs
@@ -1,21 +1,21 @@
 use super::{
     Bounds3f, DirectionCone, DiskShape, Float, Interaction, Normal3f, PI, Point2f, Point3f,
     Point3fi, QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
-    ShapeTrait, SurfaceInteraction, Transform, Vector3f,
+    SurfaceInteraction, Transform, Vector3f,
 };
 use crate::geometry::VectorLike;
 use crate::utils::math::square;
 use crate::utils::sampling::sample_uniform_disk_concentric;
 use std::sync::Arc;
 
-impl<'a> DiskShape<'a> {
+impl DiskShape {
     pub fn new(
         radius: Float,
         inner_radius: Float,
         height: Float,
         phi_max: Float,
-        render_from_object: &'a Transform<Float>,
-        object_from_render: &'a Transform<Float>,
+        render_from_object: Arc<Transform<Float>>,
+        object_from_render: Arc<Transform<Float>>,
         reverse_orientation: bool,
     ) -> Self {
         Self {
@@ -23,7 +23,7 @@ impl<'a> DiskShape<'a> {
             inner_radius,
             height,
             phi_max,
-            render_from_object,
+            render_from_object: render_from_object.clone(),
             object_from_render,
             reverse_orientation,
             transform_swap_handedness: render_from_object.swaps_handedness(),
@@ -101,21 +101,19 @@ impl<'a> DiskShape<'a> {
         );
         surf_point
     }
-}
 
-impl ShapeTrait for DiskShape<'_> {
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         self.phi_max * 0.5 * (square(self.radius) - square(self.inner_radius))
     }
 
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         self.render_from_object
             .apply_to_bounds(Bounds3f::from_points(
                 Point3f::new(-self.radius, -self.radius, self.height),
                 Point3f::new(self.radius, self.radius, self.height),
             ))
     }
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         let mut n = self
             .render_from_object
             .apply_to_normal(Normal3f::new(0., 0., 1.));
@@ -125,7 +123,7 @@ impl ShapeTrait for DiskShape<'_> {
         DirectionCone::new_from_vector(Vector3f::from(n))
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
@@ -135,7 +133,7 @@ impl ShapeTrait for DiskShape<'_> {
         }
     }
 
-    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
         let pd = sample_uniform_disk_concentric(u);
         let p_obj = Point3f::new(pd.x() * self.radius, pd.y() * self.radius, self.height);
         let pi = self
@@ -163,7 +161,7 @@ impl ShapeTrait for DiskShape<'_> {
         })
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         if let Some(t) = t_max {
             self.basic_intersect(ray, t).is_some()
         } else {
@@ -171,7 +169,7 @@ impl ShapeTrait for DiskShape<'_> {
         }
     }
 
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
         let mut ss = self.sample(u)?;
         let intr = Arc::make_mut(&mut ss.intr);
         intr.get_common_mut().time = ctx.time;
@@ -187,11 +185,11 @@ impl ShapeTrait for DiskShape<'_> {
         return Some(ss);
     }
 
-    fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
         1. / self.area()
     }
 
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let ray = ctx.spawn_ray(wi);
         if let Some(isect) = self.intersect(&ray, None) {
             let n = isect.intr.n();
diff --git a/src/shapes/mod.rs b/src/shapes/mod.rs
index b3d60bb..6ab6847 100644
--- a/src/shapes/mod.rs
+++ b/src/shapes/mod.rs
@@ -9,46 +9,46 @@ use crate::core::interaction::{Interaction, MediumInteraction, SurfaceInteractio
 use crate::core::pbrt::{Float, PI};
 use crate::geometry::{
     Bounds3f, DirectionCone, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3f,
-    Vector3fi, VectorLike,
+    Vector3fi, VectorLike
 };
 use crate::utils::math::{next_float_down, next_float_up};
 use crate::utils::transform::Transform;
-use std::sync::Arc;
+use std::sync::{Arc, Mutex};
 
 #[derive(Debug, Clone)]
-pub struct SphereShape<'a> {
+pub struct SphereShape {
     radius: Float,
     z_min: Float,
     z_max: Float,
     theta_z_min: Float,
     theta_z_max: Float,
     phi_max: Float,
-    render_from_object: &'a Transform<Float>,
-    object_from_render: &'a Transform<Float>,
+    render_from_object: Arc<Transform<Float>>,
+    object_from_render: Arc<Transform<Float>>,
     reverse_orientation: bool,
     transform_swap_handedness: bool,
 }
 
 #[derive(Debug, Clone)]
-pub struct CylinderShape<'a> {
+pub struct CylinderShape {
     radius: Float,
     z_min: Float,
     z_max: Float,
     phi_max: Float,
-    render_from_object: &'a Transform<Float>,
-    object_from_render: &'a Transform<Float>,
+    render_from_object: Arc<Transform<Float>>,
+    object_from_render: Arc<Transform<Float>>,
     reverse_orientation: bool,
     transform_swap_handedness: bool,
 }
 
 #[derive(Debug, Clone)]
-pub struct DiskShape<'a> {
+pub struct DiskShape {
     radius: Float,
     inner_radius: Float,
     height: Float,
     phi_max: Float,
-    render_from_object: &'a Transform<Float>,
-    object_from_render: &'a Transform<Float>,
+    render_from_object: Arc<Transform<Float>>,
+    object_from_render: Arc<Transform<Float>>,
     reverse_orientation: bool,
     transform_swap_handedness: bool,
 }
@@ -75,28 +75,28 @@ pub enum CurveType {
 }
 
 #[derive(Debug, Clone)]
-pub struct CurveCommon<'a> {
+pub struct CurveCommon {
     curve_type: CurveType,
     cp_obj: [Point3f; 4],
     width: [Float; 2],
     n: [Normal3f; 2],
     normal_angle: Float,
     inv_sin_normal_angle: Float,
-    render_from_object: &'a Transform<Float>,
-    object_from_render: &'a Transform<Float>,
+    render_from_object: Arc<Transform<Float>>,
+    object_from_render: Arc<Transform<Float>>,
     reverse_orientation: bool,
     transform_swap_handedness: bool,
 }
 
-impl<'a> CurveCommon<'a> {
+impl CurveCommon {
     pub fn new(
         c: &[Point3f],
         w0: Float,
         w1: Float,
         curve_type: CurveType,
         norm: &[Vector3f],
-        render_from_object: &'a Transform<Float>,
-        object_from_render: &'a Transform<Float>,
+        render_from_object: Arc<Transform<Float>>,
+        object_from_render: Arc<Transform<Float>>,
         reverse_orientation: bool,
     ) -> Self {
         let transform_swap_handedness = render_from_object.swaps_handedness();
@@ -133,8 +133,8 @@ impl<'a> CurveCommon<'a> {
 }
 
 #[derive(Debug, Clone)]
-pub struct CurveShape<'a> {
-    common: CurveCommon<'a>,
+pub struct CurveShape {
+    common: CurveCommon,
     u_min: Float,
     u_max: Float,
 }
@@ -142,21 +142,23 @@ pub struct CurveShape<'a> {
 // Define Intersection objects. This only varies for
 #[derive(Debug, Clone)]
 pub struct ShapeIntersection {
-    intr: SurfaceInteraction,
-    t_hit: Float,
+    pub intr: Box<SurfaceInteraction>,
+    pub t_hit: Float,
 }
 
 impl ShapeIntersection {
     pub fn new(intr: SurfaceInteraction, t_hit: Float) -> Self {
-        Self { intr, t_hit }
+        Self { intr: Box::new(intr), t_hit }
     }
 
-    pub fn intr(&self) -> &SurfaceInteraction {
-        &self.intr
+    pub fn as_ref(&self) -> &dyn Interaction { &*self.intr }
+
+    pub fn intr(&self) -> &dyn Interaction {
+        &*self.intr
     }
 
     pub fn intr_mut(&mut self) -> &mut SurfaceInteraction {
-        &mut self.intr
+        &mut *self.intr
     }
 
     pub fn t_hit(&self) -> Float {
@@ -266,24 +268,124 @@ impl ShapeSampleContext {
     }
 }
 
-#[derive(Debug, Clone)]
-pub enum Shape<'a> {
-    Sphere(SphereShape<'a>),
-    Cylinder(CylinderShape<'a>),
-    Disk(DiskShape<'a>),
+#[derive(Default, Debug, Clone)]
+pub enum Shape {
+    #[default]
+    None,
+    Sphere(SphereShape),
+    Cylinder(CylinderShape),
+    Disk(DiskShape),
     Triangle(TriangleShape),
     BilinearPatch(BilinearPatchShape),
-    Curve(CurveShape<'a>),
+    Curve(CurveShape),
 }
 
-pub trait ShapeTrait: Send + Sync + std::fmt::Debug {
-    fn bounds(&self) -> Bounds3f;
-    fn normal_bounds(&self) -> DirectionCone;
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection>;
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool;
-    fn area(&self) -> Float;
-    fn sample(&self, u: Point2f) -> Option<ShapeSample>;
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample>;
-    fn pdf(&self, interaction: Arc<&dyn Interaction>) -> Float;
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float;
+impl<'a> Shape {
+    pub fn bounds(&self) -> Bounds3f {
+        match self {
+            Shape::None => Bounds3f::default(),
+            Shape::Sphere(s) => s.bounds(),
+            Shape::Cylinder(c) => c.bounds(),
+            Shape::Disk(d) => d.bounds(),
+            Shape::Triangle(t) => t.bounds(),
+            Shape::BilinearPatch(b) => b.bounds(),
+            Shape::Curve(c) => c.bounds(),
+        }
+    }
+
+    pub fn normal_bounds(&self) -> DirectionCone {
+        match self {
+            Shape::None => DirectionCone::entire_sphere(),
+            Shape::Sphere(s) => s.normal_bounds(),
+            Shape::Cylinder(c) => c.normal_bounds(),
+            Shape::Disk(d) => d.normal_bounds(),
+            Shape::Triangle(t) => t.normal_bounds(),
+            Shape::BilinearPatch(b) => b.normal_bounds(),
+            Shape::Curve(c) => c.normal_bounds(),
+        }
+    }
+
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+        match self {
+            Shape::None => None,
+            Shape::Sphere(s) => s.intersect(ray, t_max),
+            Shape::Cylinder(c) => c.intersect(ray, t_max),
+            Shape::Disk(d) => d.intersect(ray, t_max),
+            Shape::Triangle(t) => t.intersect(ray, t_max),
+            Shape::BilinearPatch(b) => b.intersect(ray, t_max),
+            Shape::Curve(c) => c.intersect(ray, t_max),
+        }
+    }
+
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+        match self {
+            Shape::None => false,
+            Shape::Sphere(s) => s.intersect_p(ray, t_max),
+            Shape::Cylinder(c) => c.intersect_p(ray, t_max),
+            Shape::Disk(d) => d.intersect_p(ray, t_max),
+            Shape::Triangle(t) => t.intersect_p(ray, t_max),
+            Shape::BilinearPatch(b) => b.intersect_p(ray, t_max),
+            Shape::Curve(c) => c.intersect_p(ray, t_max),
+        }
+    }
+
+    pub fn area(&self) -> Float {
+        match self {
+            Shape::None => 0.,
+            Shape::Sphere(s) => s.area(),
+            Shape::Cylinder(c) => c.area(),
+            Shape::Disk(d) => d.area(),
+            Shape::Triangle(t) => t.area(),
+            Shape::BilinearPatch(b) => b.area(),
+            Shape::Curve(c) => c.area(),
+        }
+    }
+
+    pub fn pdf(&self, interaction: Arc<&dyn Interaction>) -> Float {
+        match self {
+            Shape::None => 0.,
+            Shape::Sphere(s) => s.pdf(interaction),
+            Shape::Cylinder(c) => c.pdf(interaction),
+            Shape::Disk(d) => d.pdf(interaction),
+            Shape::Triangle(t) => t.pdf(interaction),
+            Shape::BilinearPatch(b) => b.pdf(interaction),
+            Shape::Curve(c) => c.pdf(interaction),
+        }
+    }
+
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+        match self {
+            Shape::None => 0.,
+            Shape::Sphere(s) => s.pdf_from_context(ctx, wi),
+            Shape::Cylinder(c) => c.pdf_from_context(ctx, wi),
+            Shape::Disk(d) => d.pdf_from_context(ctx, wi),
+            Shape::Triangle(t) => t.pdf_from_context(ctx, wi),
+            Shape::BilinearPatch(b) => b.pdf_from_context(ctx, wi),
+            Shape::Curve(c) => c.pdf_from_context(ctx, wi),
+        }
+    }
+
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+        match self {
+            Shape::None => None,
+            Shape::Sphere(s) => s.sample(u),
+            Shape::Cylinder(c) => c.sample(u),
+            Shape::Disk(d) => d.sample(u),
+            Shape::Triangle(t) => t.sample(u),
+            Shape::BilinearPatch(b) => b.sample(u),
+            Shape::Curve(c) => c.sample(u),
+        }
+    }
+
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+        match self {
+            Shape::None => None,
+            Shape::Sphere(s) => s.sample_from_context(ctx, u),
+            Shape::Cylinder(c) => c.sample_from_context(ctx, u),
+            Shape::Disk(d) => d.sample_from_context(ctx, u),
+            Shape::Triangle(t) => t.sample_from_context(ctx, u),
+            Shape::BilinearPatch(b) => b.sample_from_context(ctx, u),
+            Shape::Curve(c) => c.sample_from_context(ctx, u),
+        }
+    }
 }
diff --git a/src/shapes/sphere.rs b/src/shapes/sphere.rs
index f08b9ea..558380a 100644
--- a/src/shapes/sphere.rs
+++ b/src/shapes/sphere.rs
@@ -1,6 +1,6 @@
 use super::{
     Bounds3f, DirectionCone, Float, Interaction, Normal3f, PI, Point2f, Point3f, Point3fi,
-    QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait,
+    QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
     SphereShape, SurfaceInteraction, Transform, Vector3f, Vector3fi,
 };
 use crate::core::pbrt::{clamp_t, gamma};
@@ -12,10 +12,10 @@ use crate::utils::sampling::sample_uniform_sphere;
 use std::mem;
 use std::sync::Arc;
 
-impl<'a> SphereShape<'a> {
+impl SphereShape {
     pub fn new(
-        render_from_object: &'a Transform<Float>,
-        object_from_render: &'a Transform<Float>,
+        render_from_object: Arc<Transform<Float>>,
+        object_from_render: Arc<Transform<Float>>,
         reverse_orientation: bool,
         radius: Float,
         z_min: Float,
@@ -26,8 +26,8 @@ impl<'a> SphereShape<'a> {
         let theta_z_max = clamp_t(z_max.min(z_max) / radius, -1., 1.).acos();
         let phi_max = radians(clamp_t(phi_max, 0., 360.0));
         Self {
-            render_from_object,
-            object_from_render,
+            render_from_object: render_from_object.clone(),
+            object_from_render: object_from_render.clone(),
             radius,
             z_min: clamp_t(z_min.min(z_max), -radius, radius),
             z_max: clamp_t(z_min.max(z_max), -radius, radius),
@@ -193,10 +193,8 @@ impl<'a> SphereShape<'a> {
         // self.render_from_object.apply_to_point(surf_point)
         surf_point
     }
-}
 
-impl ShapeTrait for SphereShape<'_> {
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         self.render_from_object
             .apply_to_bounds(Bounds3f::from_points(
                 Point3f::new(-self.radius, -self.radius, self.z_min),
@@ -204,19 +202,19 @@ impl ShapeTrait for SphereShape<'_> {
             ))
     }
 
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         DirectionCone::entire_sphere()
     }
 
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         self.phi_max * self.radius * (self.z_max - self.z_min)
     }
 
-    fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
         1. / self.area()
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
@@ -226,7 +224,7 @@ impl ShapeTrait for SphereShape<'_> {
         }
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         if let Some(t) = t_max {
             self.basic_intersect(ray, t).is_some()
         } else {
@@ -234,7 +232,7 @@ impl ShapeTrait for SphereShape<'_> {
         }
     }
 
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let p_center = self
             .object_from_render
             .apply_to_point(Point3f::new(0., 0., 0.));
@@ -262,7 +260,7 @@ impl ShapeTrait for SphereShape<'_> {
         1. / (2. * PI * one_minus_cos_theta_max)
     }
 
-    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
         let p_obj = Point3f::new(0., 0., 0.) + self.radius * sample_uniform_sphere(u);
         let mut p_obj_vec = Vector3f::from(p_obj);
         p_obj_vec *= self.radius / p_obj.distance(Point3f::zero());
@@ -294,7 +292,7 @@ impl ShapeTrait for SphereShape<'_> {
         })
     }
 
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
         let p_center = self.render_from_object.apply_to_point(Point3f::zero());
         let p_origin = ctx.offset_ray_origin_from_point(p_center);
         if p_origin.distance_squared(p_center) <= square(self.radius) {
diff --git a/src/shapes/triangle.rs b/src/shapes/triangle.rs
index f7d4e60..3fe3d35 100644
--- a/src/shapes/triangle.rs
+++ b/src/shapes/triangle.rs
@@ -1,6 +1,6 @@
 use super::{
     Bounds3f, DirectionCone, Float, Interaction, Normal3f, Point2f, Point3f, Point3fi, Ray,
-    ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait, SurfaceInteraction,
+    ShapeIntersection, ShapeSample, ShapeSampleContext, SurfaceInteraction,
     TriangleIntersection, TriangleShape, Vector2f, Vector3f,
 };
 use crate::core::pbrt::gamma;
@@ -323,15 +323,13 @@ impl TriangleShape {
         isect.dndu = dndu;
         isect.dndv = dndv;
     }
-}
 
-impl ShapeTrait for TriangleShape {
-    fn bounds(&self) -> Bounds3f {
+    pub fn bounds(&self) -> Bounds3f {
         let [p0, p1, p2] = self.get_data().vertices;
         Bounds3f::from_points(p0, p1).union_point(p2)
     }
 
-    fn normal_bounds(&self) -> DirectionCone {
+    pub fn normal_bounds(&self) -> DirectionCone {
         let data = self.get_data();
         let mut n = data.normal;
         if let Some([n0, n1, n2]) = data.normals {
@@ -342,15 +340,15 @@ impl ShapeTrait for TriangleShape {
         DirectionCone::new_from_vector(Vector3f::from(n))
     }
 
-    fn area(&self) -> Float {
+    pub fn area(&self) -> Float {
         self.get_data().area
     }
 
-    fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
         1. / self.area()
     }
 
-    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
+    pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let solid_angle = self.solid_angle(ctx.p());
         if solid_angle < Self::MIN_SPHERICAL_SAMPLE_AREA
             || solid_angle > Self::MAX_SPHERICAL_SAMPLE_AREA
@@ -387,7 +385,7 @@ impl ShapeTrait for TriangleShape {
         pdf
     }
 
-    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
         let data = self.get_data();
         let [p0, p1, p2] = data.vertices;
         let solid_angle = self.solid_angle(ctx.p());
@@ -472,7 +470,7 @@ impl ShapeTrait for TriangleShape {
         })
     }
 
-    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
+    pub fn sample(&self, u: Point2f) -> Option<ShapeSample> {
         let data = self.get_data();
         let [p0, p1, p2] = data.vertices;
         let [uv0, uv1, uv2] = data.uvs;
@@ -503,15 +501,15 @@ impl ShapeTrait for TriangleShape {
         })
     }
 
-    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY))
             .map(|ti| {
                 let intr = self.interaction_from_intersection(ti, ray.time, -ray.d);
-                ShapeIntersection { intr, t_hit: ti.t }
+                ShapeIntersection { intr: Box::new(intr), t_hit: ti.t }
             })
     }
 
-    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
         self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY))
             .is_some()
     }
diff --git a/src/utils/interval.rs b/src/utils/interval.rs
index 7b9dab0..329cc15 100644
--- a/src/utils/interval.rs
+++ b/src/utils/interval.rs
@@ -44,6 +44,21 @@ impl Interval {
     pub fn is_empty(&self) -> bool {
         self.low > self.high
     }
+
+    pub fn abs(&self) -> Self {
+    if self.low >= 0.0 {
+        return *self;
+    }
+
+    if self.high < 0.0 {
+        return -(*self);
+    }
+
+    Self {
+        low: 0.0,
+        high: next_float_up((-self.low).max(self.high)),
+    }
+}
 }
 
 impl Default for Interval {
@@ -65,6 +80,20 @@ impl Add for Interval {
     }
 }
 
+impl Add<Interval> for Float {
+    type Output = Interval;
+    fn add(self, rhs: Interval) -> Self::Output {
+        Interval::new(self) + rhs
+    }
+}
+
+impl Sub<Interval> for Float {
+    type Output = Interval;
+    fn sub(self, rhs: Interval) -> Self::Output {
+        Interval::new(self) - rhs
+    }
+}
+
 impl Sub for Interval {
     type Output = Self;
     fn sub(self, rhs: Self) -> Self::Output {
diff --git a/src/utils/transform.rs b/src/utils/transform.rs
index 318e760..952b3b2 100644
--- a/src/utils/transform.rs
+++ b/src/utils/transform.rs
@@ -2,10 +2,12 @@ use num_traits::Float as NumFloat;
 use std::error::Error;
 use std::fmt::{self, Display};
 use std::ops::{Add, Div, Index, IndexMut, Mul};
+use std::sync::Arc;
 
 use super::color::{RGB, XYZ};
 use super::math::{SquareMatrix, safe_acos};
 use super::quaternion::Quaternion;
+use crate::core::interaction::{SurfaceInteraction, MediumInteraction, Interaction, InteractionData};
 use crate::core::pbrt::{Float, gamma};
 use crate::geometry::{
     Bounds3f, Normal, Normal3f, Point, Point3f, Point3fi, Ray, Vector, Vector3f, Vector3fi,
@@ -51,8 +53,13 @@ impl<T: NumFloat> Transform<T> {
     }
 
     pub fn apply_inverse_vector(&self, v: Vector<T, 3>) -> Vector<T, 3> {
-        self.clone() * v
-    }
+        let x = v.x();
+        let y = v.y();
+        let z = v.z();
+        Vector::<T, 3>::new(self.m_inv[0][0] * x + self.m_inv[0][1] * y + self.m_inv[0][2] * z,
+                          self.m_inv[1][0] * x + self.m_inv[1][1] * y + self.m_inv[1][2] * z,
+                          self.m_inv[2][0] * x + self.m_inv[2][1] * y + self.m_inv[2][2] * z)
+    }    
 
     pub fn apply_inverse_normal(&self, n: Normal<T, 3>) -> Normal<T, 3> {
         self.clone() * n
@@ -208,6 +215,120 @@ impl Transform<Float> {
         }
     }
 
+    pub fn apply_to_interaction(&self, inter: &dyn Interaction) -> Box<dyn Interaction> {
+        if let Some(si) = inter.as_any().downcast_ref::<SurfaceInteraction>() {
+            let mut ret = si.clone();
+            ret.common.pi = self.apply_to_interval(&si.common.pi);
+            let n = self.apply_to_normal(si.common.n);
+            ret.common.wo = self.apply_to_vector(si.common.wo).normalize();
+
+            ret.dpdu = self.apply_to_vector(si.dpdu);
+            ret.dpdv = self.apply_to_vector(si.dpdv);
+            ret.dndu = self.apply_to_normal(si.dndu);
+            ret.dndv = self.apply_to_normal(si.dndv);
+            ret.dpdx = self.apply_to_vector(si.dpdx);
+            ret.dpdy = self.apply_to_vector(si.dpdy);
+
+            let shading_n = self.apply_to_normal(si.shading.n);
+            ret.shading.n = shading_n.normalize();
+            ret.shading.dpdu = self.apply_to_vector(si.shading.dpdu);
+            ret.shading.dpdv = self.apply_to_vector(si.shading.dpdv);
+            ret.shading.dndu = self.apply_to_normal(si.shading.dndu);
+            ret.shading.dndv = self.apply_to_normal(si.shading.dndv);
+
+            ret.common.n = n.normalize().face_forward(shading_n.into());
+
+            Box::new(ret)
+        } else if let Some(mi) = inter.as_any().downcast_ref::<MediumInteraction>() {
+        let ret = MediumInteraction {
+            common: InteractionData {
+                pi: self.apply_to_interval(&mi.common.pi),
+                n: self.apply_to_normal(mi.common.n).normalize(),
+                wo: self.apply_to_vector(mi.common.wo).normalize(),
+                time: mi.common.time,
+                medium_interface: mi.common.medium_interface.clone(),
+                medium: mi.common.medium.clone(),
+            },
+            medium: mi.medium.clone(),
+            phase: mi.phase.clone(),
+        };
+        Box::new(ret)
+        } else {
+            panic!("Unhandled Interaction type in Transform::apply_to_interaction");
+        }
+    }
+
+    pub fn apply_inverse_interval(&self, p: &Point3fi) -> Point3fi {
+        let x = Float::from(p.x());
+        let y = Float::from(p.y());
+        let z = Float::from(p.z());
+        let m_inv = &self.m_inv;
+
+        // Compute transformed coordinates from point
+        let xp = (m_inv[0][0] * x + m_inv[0][1] * y) + (m_inv[0][2] * z + m_inv[0][3]);
+        let yp = (m_inv[1][0] * x + m_inv[1][1] * y) + (m_inv[1][2] * z + m_inv[1][3]);
+        let zp = (m_inv[2][0] * x + m_inv[2][1] * y) + (m_inv[2][2] * z + m_inv[2][3]);
+        let wp = (m_inv[3][0] * x + m_inv[3][1] * y) + (m_inv[3][2] * z + m_inv[3][3]);
+
+        // Compute absolute error for transformed point
+        let p_out_error: Vector3f;
+        let g3 = gamma(3);
+
+        if p.is_exact() {
+            p_out_error = Vector3f::new(
+                g3 * (m_inv[0][0] * x).abs() + (m_inv[0][1] * y).abs() + (m_inv[0][2] * z).abs(),
+                g3 * (m_inv[1][0] * x).abs() + (m_inv[1][1] * y).abs() + (m_inv[1][2] * z).abs(),
+                g3 * (m_inv[2][0] * x).abs() + (m_inv[2][1] * y).abs() + (m_inv[2][2] * z).abs(),
+            );
+        } else {
+            let p_in_error = p.error();
+            let g3_plus_1 = g3 + 1.0;
+            
+            p_out_error = Vector3f::new(
+                g3_plus_1 * (m_inv[0][0].abs() * p_in_error.x() +
+                             m_inv[0][1].abs() * p_in_error.y() +
+                             m_inv[0][2].abs() * p_in_error.z()) +
+                g3 * ((m_inv[0][0] * x).abs() + (m_inv[0][1] * y).abs() +
+                      (m_inv[0][2] * z).abs() + m_inv[0][3].abs()),
+                
+                g3_plus_1 * (m_inv[1][0].abs() * p_in_error.x() +
+                             m_inv[1][1].abs() * p_in_error.y() +
+                             m_inv[1][2].abs() * p_in_error.z()) +
+                g3 * ((m_inv[1][0] * x).abs() + (m_inv[1][1] * y).abs() +
+                      (m_inv[1][2] * z).abs() + m_inv[1][3].abs()),
+                
+                g3_plus_1 * (m_inv[2][0].abs() * p_in_error.x() +
+                             m_inv[2][1].abs() * p_in_error.y() +
+                             m_inv[2][2].abs() * p_in_error.z()) +
+                g3 * ((m_inv[2][0] * x).abs() + (m_inv[2][1] * y).abs() +
+                      (m_inv[2][2] * z).abs() + m_inv[2][3].abs()),
+            );
+        }
+
+        if wp == 1.0 {
+            Point3fi::new_with_error(Point3f::new(xp, yp, zp), p_out_error)
+        } else {
+            Point3fi::new_with_error(Point3f::new(xp/wp, yp/wp, zp/wp), p_out_error)
+        }
+    }
+
+    pub fn apply_inverse_ray(&self, r: &Ray, t_max: Option<Float>) -> (Ray, Float) {
+        let mut o = self.apply_inverse_interval(&Point3fi::new_from_point(r.o));
+        let d = self.apply_inverse_vector(r.d);
+        // Offset ray origin to edge of error bounds and compute _tMax_
+        let mut t = 0.;
+        let length_squared = d.norm_squared();
+        if length_squared > 0. {
+            let o_error = Vector3f::new(o.x().width() / 2., o.y().width() / 2., o.z().width() / 2.);
+            let dt = d.abs().dot(o_error) / length_squared;
+            o = o + Vector3fi::from(d * dt);
+            if let Some(t_max) = t_max {
+                t = t_max - dt;
+            }
+        }
+        (Ray::new(Point3f::from(o), d, Some(r.time), r.medium.clone()), t)
+    }
+
     pub fn to_quaternion(&self) -> Quaternion {
         let trace = self.m.trace();
         let mut quat = Quaternion::default();
@@ -580,6 +701,7 @@ impl DerivativeTerm {
     }
 }
 
+#[derive(Debug, Clone)]
 pub struct AnimatedTransform {
     pub start_transform: Transform<Float>,
     pub end_transform: Transform<Float>,
@@ -1800,6 +1922,14 @@ impl AnimatedTransform {
         t.apply_to_ray(r, t_max)
     }
 
+    pub fn apply_interaction(&self, si: &dyn Interaction) -> Box<dyn Interaction> {
+        if !self.actually_animated {
+            return self.start_transform.apply_to_interaction(si)
+        }
+        let t = self.interpolate(si.time());
+        t.apply_to_interaction(si)
+    }
+
     pub fn interpolate(&self, time: Float) -> Transform<Float> {
         if !self.actually_animated || time <= self.start_time {
             return self.start_transform;
@@ -1840,6 +1970,13 @@ impl AnimatedTransform {
         }
         return self.interpolate(time).apply_inverse_normal(n);
     }
+
+    pub fn motion_bounds(&self, b: &Bounds3f) -> Bounds3f {
+        if !self.actually_animated {
+            return self.start_transform.apply_to_bounds(*b)
+        }
+        return self.start_transform.apply_to_bounds(*b).union(self.end_transform.apply_to_bounds(*b))
+    }
 }
 
 pub fn look_at(