Refactored, added camera types, filters, films, starting on reflection and scattering

Cargo.toml

@@ -4,5 +4,8 @@ version = "0.1.0"
 edition = "2024"
 
 [dependencies]
+bitflags = "2.10.0"
+bumpalo = "3.19.0"
 num-traits = "0.2.19"
 once_cell = "1.21.3"
+rand = "0.9.2"

src/camera/mod.rs

@@ -0,0 +1,211 @@
+mod perspective;
+mod orthographic;
+mod spherical;
+mod realistic;
+
+pub use perspective::PerspectiveCamera;
+pub use orthographic::OrthographicCamera;
+pub use spherical::SphericalCamera;
+pub use realistic::RealisticCamera;
+
+use crate::core::film::Film;
+use crate::core::medium::Medium;
+use crate::core::pbrt::{Float, lerp, RenderingCoordinateSystem};
+use crate::core::sampler::CameraSample;
+use crate::utils::geometry::{Ray, RayDifferential, Vector3f, Point3f, Normal3f, Dot, Normed};
+use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+use crate::utils::transform::{AnimatedTransform, Transform};
+
+use std::sync::Arc;
+
+#[derive(Debug, Clone)]
+pub struct CameraRay {
+    pub ray: Ray,
+    pub weight: SampledSpectrum,
+}
+
+pub struct CameraTransform {
+    render_from_camera: AnimatedTransform,
+    world_from_render: Transform<Float>,
+}
+
+impl CameraTransform {
+    pub fn from_world(world_from_camera: AnimatedTransform, rendering_space: RenderingCoordinateSystem) -> Self {
+        let world_from_render = match rendering_space {
+            RenderingCoordinateSystem::Camera => {
+                let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.;
+                world_from_camera.interpolate(t_mid)
+            }
+            RenderingCoordinateSystem::CameraWorld => {
+                let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.;
+                let p_camera = world_from_camera.apply_point(Point3f::default(), t_mid);
+                Transform::translate(p_camera - Point3f::new(0., 0., 0.))
+            }
+            RenderingCoordinateSystem::World => Transform::identity(),
+        };
+        let render_from_world = world_from_render.inverse();
+        let rfc = [render_from_world * world_from_camera.start_transform, render_from_world * world_from_camera.end_transform];
+        let render_from_camera = AnimatedTransform::new(&rfc[0], world_from_camera.start_time, &rfc[1], world_from_camera.end_time).expect("Render wrong");
+        Self { render_from_camera, world_from_render }
+
+    }
+
+    pub fn render_from_camera(&self, p: Point3f, time: Float) -> Point3f {
+        self.render_from_camera.apply_point(p, time)
+    }
+
+    pub fn render_from_camera_vector(&self, v: Vector3f, time: Float) -> Vector3f {
+        self.render_from_camera.apply_vector(v, time)
+    }
+
+    pub fn render_from_camera_ray(&self, r: &Ray, t_max: &mut Option<Float>) -> Ray {
+        self.render_from_camera.apply_ray(r, t_max)
+    }
+
+    pub fn camera_from_render(&self, p: Point3f, time: Float) -> Point3f {
+        self.render_from_camera.apply_inverse_point(p, time)
+    }
+
+    pub fn camera_from_render_vector(&self, v: Vector3f, time: Float) -> Vector3f {
+        self.render_from_camera.apply_inverse_vector(v, time)
+    }
+
+    pub fn camera_from_render_normal(&self, n: Normal3f, time: Float) -> Normal3f {
+        self.render_from_camera.apply_inverse_normal(n, time)
+    }
+}
+
+pub struct CameraBase {
+    pub camera_transform: CameraTransform,
+    pub shutter_open: Float,
+    pub shutter_close: Float,
+    pub film: Film,
+    pub medium: Option<Arc<Medium>>,
+    pub min_pos_differential_x: Vector3f,
+    pub min_pos_differential_y: Vector3f,
+    pub min_dir_differential_x: Vector3f,
+    pub min_dir_differential_y: Vector3f,
+}
+
+pub trait CameraTrait {
+    fn base(&self) -> &CameraBase;
+    fn generate_ray(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay>;
+    fn generate_ray_differential(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay> {
+        let mut central_cam_ray = match self.generate_ray(sample, lambda) {
+            Some(cr) => cr,
+            None => return None,
+        };
+
+        let mut rd = RayDifferential::default();
+        let mut rx_found = false;
+        let mut ry_found = false;
+
+        for eps in [0.05, -0.05] {
+            let mut s_shift = sample;
+            s_shift.p_film[0] += eps;
+
+            if let Some(rx_cam_ray) = self.generate_ray(s_shift, lambda) {
+                rd.rx_origin = central_cam_ray.ray.o + (rx_cam_ray.ray.o - central_cam_ray.ray.o) / eps;
+                rd.rx_direction = central_cam_ray.ray.d + (rx_cam_ray.ray.d - central_cam_ray.ray.d) / eps;
+                rx_found = true;
+                break;
+            }
+        }
+
+        for eps in [0.05, -0.05] {
+            let mut s_shift = sample;
+            s_shift.p_film[1] += eps;
+
+            if let Some(ry_cam_ray) = self.generate_ray(s_shift, lambda) {
+                rd.ry_origin = central_cam_ray.ray.o + (ry_cam_ray.ray.o - central_cam_ray.ray.o) / eps;
+                rd.ry_direction = central_cam_ray.ray.d + (ry_cam_ray.ray.d - central_cam_ray.ray.d) / eps;
+                ry_found = true;
+                break
+            }
+        }
+
+        if rx_found && ry_found { 
+            central_cam_ray.ray.differential = Some(rd);
+        }
+
+        Some(central_cam_ray)
+    }
+
+    fn sample_time(&self, u: Float) -> Float { 
+        lerp(u, self.base().shutter_open, self.base().shutter_close)
+    }
+
+    fn approximate_dp_dxy(&self, p: Point3f, n: Normal3f, time: Float, samples_per_pixel: i32, dpdx: &mut Vector3f, dpdy: &mut Vector3f) {
+        let p_camera = self.base().camera_transform.camera_from_render(p, time);
+        let p_camera_vec = p_camera - Point3f::new(0., 0., 0.);
+        let down_z_from_camera = Transform::rotate_from_to(p_camera_vec.normalize(), Vector3f::new(1., 0., 0.));
+        let p_down_z = down_z_from_camera.apply_to_point(p_camera);
+        let n_down_z = down_z_from_camera.apply_to_normal(self.base().camera_transform.camera_from_render_normal(n, time));
+        let d = n_down_z.z() * p_down_z.z();
+        let x_ray = Ray::new(Point3f::new(0., 0., 0.) + self.base().min_pos_differential_x, 
+            Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_x, 
+            None, None);
+        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, None);
+        let tx = -(n_down_z.dot(y_ray.o)) / n_down_z.dot(x_ray.d);
+        let ty = -(n_down_z.dot(x_ray.o) - d) / n_down_z.dot(y_ray.d);
+        let px = x_ray.evaluate(tx);
+        let py = y_ray.evaluate(ty);
+        let spp_scale = 0.125_f32.max((samples_per_pixel as Float).sqrt());
+        *dpdx = spp_scale * self.base().camera_transform.render_from_camera_vector(down_z_from_camera.apply_inverse_vector(px - p_down_z), time);
+        *dpdy = spp_scale * self.base().camera_transform.render_from_camera_vector(down_z_from_camera.apply_inverse_vector(py - p_down_z), time);
+    }
+
+    fn render_from_camera(&self, r: &Ray, t_max: &mut Option<Float>) -> Ray {
+        self.base().camera_transform.render_from_camera_ray(r, t_max)
+    }
+
+}
+
+pub enum Camera {
+    Perspective(PerspectiveCamera),
+    Orthographic(OrthographicCamera),
+    Spherical(SphericalCamera),
+    Realistic(RealisticCamera),
+}
+impl CameraTrait for Camera {
+    fn base(&self) -> &CameraBase {
+        match self {
+            Camera::Perspective(c) => c.base(),
+            Camera::Orthographic(c) => c.base(),
+            Camera::Spherical(c) => c.base(),
+            Camera::Realistic(c) => c.base(),
+        }
+    }
+
+    fn generate_ray(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay> {
+        match self {
+            Camera::Perspective(c) => c.generate_ray(sample, lambda),
+            Camera::Orthographic(c) => c.generate_ray(sample, lambda),
+            Camera::Spherical(c) => c.generate_ray(sample, lambda),
+            Camera::Realistic(c) => c.generate_ray(sample, lambda),
+        }
+    }
+
+    fn generate_ray_differential(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay> {
+        match self {
+            Camera::Perspective(c) => c.generate_ray_differential(sample, lambda),
+            Camera::Orthographic(c) => c.generate_ray_differential(sample, lambda),
+            Camera::Spherical(c) => c.generate_ray_differential(sample, lambda),
+            Camera::Realistic(c) => c.generate_ray_differential(sample, lambda),
+        }
+    }
+}
+
+pub struct LensElementInterface {
+    pub curvature_radius: Float,
+    pub thickness: Float,
+    pub eta: Float,
+    pub aperture_radius: Float,
+}
+
+struct ExitPupilSample {
+    pub p_pupil: Point3f,
+    pub pdf: Float,
+}

src/camera/orthographic.rs

@@ -0,0 +1,96 @@
+use super::{CameraBase, CameraRay, CameraTrait};
+use crate::core::sampler::CameraSample;
+use crate::core::pbrt::{Float, sample_uniform_disk_concentric};
+use crate::core::film::FilmTrait;
+use crate::utils::transform::Transform;
+use crate::utils::geometry::{Bounds2f, Point3f, Vector3f, Ray, RayDifferential, Normed};
+use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+
+pub struct OrthographicCamera {
+    pub base: CameraBase,
+    pub screen_from_camera: Transform<Float>,
+    pub camera_from_raster: Transform<Float>,
+    pub raster_from_screen: Transform<Float>,
+    pub screen_from_raster: Transform<Float>,
+    pub lens_radius: Float, 
+    pub focal_distance: Float,
+    pub dx_camera: Vector3f,
+    pub dy_camera: Vector3f,
+}
+
+impl OrthographicCamera {
+    pub fn new(base: CameraBase, screen_window: Bounds2f, lens_radius: Float, focal_distance: Float) -> Self {
+        let ndc_from_screen: Transform<Float> =  Transform::scale(
+            1./(screen_window.p_max.x() - screen_window.p_min.x()), 
+            1. / (screen_window.p_max.y() - screen_window.p_min.y()), 1.) * 
+            Transform::translate(Vector3f::new(-screen_window.p_min.x(), -screen_window.p_max.y(), 0.));
+        let raster_from_ndc = Transform::scale(base.film.full_resolution().x() as Float, -base.film.full_resolution().y() as Float, 1.);
+        let raster_from_screen = raster_from_ndc * ndc_from_screen;
+        let screen_from_raster = raster_from_screen.inverse();
+        let screen_from_camera = Transform::orthographic(0., 1.);
+        let camera_from_raster = screen_from_camera.inverse() * screen_from_raster;
+        let dx_camera = camera_from_raster.apply_to_vector(Vector3f::new(1., 0., 0.));
+        let dy_camera = camera_from_raster.apply_to_vector(Vector3f::new(0., 1., 0.));
+        let mut base_ortho = base;
+        base_ortho.min_dir_differential_x = Vector3f::new(0., 0., 0.);
+        base_ortho.min_dir_differential_y = Vector3f::new(0., 0., 0.);
+        base_ortho.min_pos_differential_x = dx_camera;
+        base_ortho.min_pos_differential_y = dy_camera;
+        Self { base: base_ortho, screen_from_camera, camera_from_raster, raster_from_screen, screen_from_raster, lens_radius, focal_distance, dx_camera, dy_camera }
+    }
+}
+
+
+impl CameraTrait for OrthographicCamera {
+    fn base(&self) -> &CameraBase {
+        &self.base
+    }
+    fn generate_ray(&self, sample: CameraSample, _lambda: &SampledWavelengths) -> Option<CameraRay> {
+        // Compute raster and camera sample positions
+        let p_film = Point3f::new(sample.p_film.x(), sample.p_film.y(), 0.);
+        let p_camera = self.camera_from_raster.apply_to_point(p_film);
+
+        let mut ray = Ray::new(p_camera, Vector3f::new(0., 0., 1.), Some(self.sample_time(sample.time)), self.base().medium.clone());
+        // Modify ray for depth of field
+        if self.lens_radius > 0. {
+            // Sample point on lens
+            let p_lens = self.lens_radius * sample_uniform_disk_concentric(sample.p_lens);
+
+            // Compute point on plane of focus
+            let ft = self.focal_distance / ray.d.z();
+            let p_focus = ray.evaluate(ft);
+
+            // Update ray for effect of lens
+            ray.o = Point3f::new(p_lens.x(), p_lens.y(), 0.);
+            ray.d = (p_focus - ray.o).normalize();
+        }
+
+        let camera_ray = self.render_from_camera(&ray, &mut None);
+
+        Some(CameraRay{ray: camera_ray, weight: SampledSpectrum::default()})
+    }
+
+    fn generate_ray_differential(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay> {
+        let mut central_cam_ray = match self.generate_ray(sample, lambda) {
+            Some(cr) => cr,
+            None => return None,
+        };
+
+        let mut rd = RayDifferential::default();
+        if self.lens_radius > 0.0 {
+            return self.generate_ray_differential(sample, lambda);
+        } else {
+            let time = self.sample_time(sample.time);
+            let world_dx = self.base.camera_transform.render_from_camera_vector(self.dx_camera, time);
+            let world_dy = self.base.camera_transform.render_from_camera_vector(self.dy_camera, time);
+
+            rd.rx_origin = central_cam_ray.ray.o + world_dx;
+            rd.ry_origin = central_cam_ray.ray.o + world_dy;
+
+            rd.rx_direction = central_cam_ray.ray.d;
+            rd.ry_direction = central_cam_ray.ray.d;
+        }
+        central_cam_ray.ray.differential = Some(rd);
+        Some(central_cam_ray)
+    }
+}

src/camera/perspective.rs

@@ -0,0 +1,82 @@
+use super::{CameraRay, CameraBase, CameraTrait};
+use crate::camera;
+use crate::core::film::FilmTrait;
+use crate::core::filter::FilterTrait;
+use crate::core::pbrt::{Float, sample_uniform_disk_concentric};
+use crate::core::sampler::CameraSample;
+use crate::utils::geometry::{Point2f, Point3f, Bounds2f, Vector3f, Normed, Ray, RayDifferential};
+use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+use crate::utils::transform::Transform;
+
+pub struct PerspectiveCamera {
+    pub base: CameraBase,
+    pub screen_from_camera: Transform<Float>,
+    pub camera_from_raster: Transform<Float>,
+    pub raster_from_screen: Transform<Float>,
+    pub screen_from_raster: Transform<Float>,
+    pub lens_radius: Float, 
+    pub focal_distance: Float,
+    pub dx_camera: Vector3f,
+    pub dy_camera: Vector3f,
+    pub cos_total_width: Float,
+}
+
+impl PerspectiveCamera {
+    pub fn new(base: CameraBase, screen_from_camera: &Transform<Float>, screen_window: Bounds2f, lens_radius: Float, focal_distance: Float) -> Self {
+        let ndc_from_screen: Transform<Float> =  Transform::scale(
+            1./(screen_window.p_max.x() - screen_window.p_min.x()), 
+            1. / (screen_window.p_max.y() - screen_window.p_min.y()), 1.) * 
+            Transform::translate(Vector3f::new(-screen_window.p_min.x(), -screen_window.p_max.y(), 0.));
+        let raster_from_ndc = Transform::scale(base.film.full_resolution().x() as Float, -base.film.full_resolution().y() as Float, 1.);
+        let raster_from_screen = raster_from_ndc * ndc_from_screen;
+        let screen_from_raster = raster_from_screen.inverse();
+        let camera_from_raster = screen_from_camera.inverse() * screen_from_raster;
+        let dx_camera = camera_from_raster.apply_to_point(Point3f::new(1., 0., 0.)) - camera_from_raster.apply_to_point(Point3f::new(0., 0., 0.));
+        let dy_camera = camera_from_raster.apply_to_point(Point3f::new(0., 1., 0.)) - camera_from_raster.apply_to_point(Point3f::new(0., 0., 0.));
+        let radius = base.film.get_filter().radius();
+        let p_corner = Point3f::new(-radius.x(), -radius.y(), 0.);
+        let w_corner_camera = (camera_from_raster.apply_to_point(p_corner) - Point3f::new(0., 0., 0.)).normalize();
+        let cos_total_width = w_corner_camera.z();
+        Self { base, 
+            screen_from_camera: *screen_from_camera, 
+            camera_from_raster, 
+            raster_from_screen, 
+            screen_from_raster, 
+            lens_radius, 
+            focal_distance,
+            dx_camera,
+            dy_camera,
+            cos_total_width,
+        }
+    }
+}
+
+impl CameraTrait for PerspectiveCamera {
+    fn base(&self) -> &CameraBase {
+        &self.base
+    }
+    fn generate_ray(&self, sample: CameraSample, _lambda: &SampledWavelengths) -> Option<CameraRay> {
+        // Compute raster and camera sample positions
+        let p_film = Point3f::new(sample.p_film.x(), sample.p_film.y(), 0.);
+        let p_camera = self.camera_from_raster.apply_to_point(p_film);
+        let p_vector = p_camera - Point3f::new(0., 0., 0.);
+
+        let mut r = Ray::new(Point3f::new(0., 0., 0.), p_vector.normalize(), Some(self.sample_time(sample.time)), self.base().medium.clone());
+        // Modify ray for depth of field
+        if self.lens_radius > 0. {
+            // Sample point on lens
+            let p_lens = self.lens_radius * sample_uniform_disk_concentric(sample.p_lens);
+
+            // Compute point on plane of focus
+            let ft = self.focal_distance / r.d.z();
+            let p_focus = r.evaluate(ft);
+
+            // Update ray for effect of lens
+            r.o = Point3f::new(p_lens.x(), p_lens.y(), 0.);
+            r.d = (p_focus - r.o).normalize();
+        }
+
+        let ray = self.render_from_camera(&r, &mut None);
+        Some(CameraRay{ray, weight: SampledSpectrum::default()})
+    }
+}

src/camera/realistic.rs

@@ -0,0 +1,261 @@
+use super::{CameraBase, LensElementInterface, CameraRay, CameraTrait, ExitPupilSample};
+use crate::core::film::FilmTrait;
+use crate::core::pbrt::{Float, square, lerp};
+use crate::core::sampler::CameraSample;
+use crate::utils::math::quadratic;
+use crate::utils::geometry::{Bounds2f, Point2f, Point3f, Point2i, Normal3f, Vector2i, Vector3f, Ray, Normed, Dot};
+use crate::utils::scattering::refract;
+use crate::utils::spectrum::{SampledWavelengths, SampledSpectrum};
+
+pub struct RealisticCamera {
+    base: CameraBase,
+    focus_distance: Float,
+    set_aperture_diameter: Float,
+    // aperture_image: Image,
+    element_interface: Vec<LensElementInterface>,
+    physical_extent: Bounds2f,
+    exit_pupil_bounds: Vec<Bounds2f>,
+}
+
+impl RealisticCamera {
+    pub fn new(&self, 
+        base: CameraBase, 
+        lens_params: Vec<Float>, 
+        focus_distance: Float, 
+        set_aperture_diameter: Float,
+        ) -> Self {
+
+        let aspect = base.film.full_resolution().x() as Float / base.film.full_resolution().y() as Float;
+        let diagonal = base.film.diagonal();
+        let x = (square(diagonal) / (1.0 + square(diagonal))).sqrt();
+        let y = x * aspect;
+        let physical_extent = Bounds2f::from_points(Point2f::new(-x / 2., -y / 2.), Point2f::new(x / 2., y / 2.));
+        let mut element_interface: Vec<LensElementInterface> = Vec::new();
+
+        for i in (0..lens_params.len()).step_by(4) {
+            let curvature_radius = lens_params[i] / 1000.0;
+            let thickness = lens_params[i + 1] / 1000.0;
+            let eta = lens_params[i + 2];
+            let mut aperture_diameter = lens_params[i + 3] / 1000.0;
+
+            if curvature_radius == 0.0 {
+                aperture_diameter /= 1000.0;
+                if set_aperture_diameter > aperture_diameter {
+                    println!("Aperture is larger than possible")
+                } else {
+                    aperture_diameter = set_aperture_diameter;
+                }
+            }
+            let el_int = LensElementInterface { curvature_radius, thickness, eta, aperture_radius: aperture_diameter / 2.0 };
+            element_interface.push(el_int);
+        }
+
+        let mut exit_pupil_bounds: Vec<Bounds2f> = Vec::new();
+ 
+        let n_samples = 64;
+        for i in 0..64 {
+            let r0 = i  as Float / 64. * base.film.diagonal() / 2.;
+            let r1 = (i + 1) as Float / n_samples as Float * base.film.diagonal() / 2.;
+            exit_pupil_bounds[i] = self.bound_exit_pupil(r0, r1);
+        }
+
+        Self { base, focus_distance, element_interface, physical_extent, set_aperture_diameter, exit_pupil_bounds }
+    }
+
+    pub fn lens_rear_z(&self) -> Float  { self.element_interface.last().unwrap().thickness }
+
+    pub fn lens_front_z(&self) -> Float {
+        let mut z_sum = 0.;
+        for element in &self.element_interface {
+            z_sum += element.thickness;
+        }
+        z_sum
+    }
+
+    pub fn rear_element_radius(&self) -> Float { self.element_interface.last().unwrap().aperture_radius }
+
+    pub fn bound_exit_pupil(&self, film_x_0: Float, film_x_1: Float) -> Bounds2f {
+        let mut pupil_bounds = Bounds2f::default();
+        let n_samples = 1024 * 1024;
+        let rear_radius = self.rear_element_radius();
+        let proj_rear_bounds = Bounds2f::from_points(Point2f::new(-1.5 * rear_radius, -1.5 * rear_radius),
+                                Point2f::new(1.5 * rear_radius, 1.5 * rear_radius));
+
+        let radical_inverse = |x: i32, _y: i64| x as Float;
+        let lens_rear_z = || 1. ;
+        let trace_lenses_from_film = |_ray: Ray, _place: Option<Ray>| true;
+        for i in 0..n_samples {
+            // Find location of sample points on $x$ segment and rear lens element
+            //
+            let p_film = Point3f::new(lerp((i as Float + 0.5) / n_samples as Float, film_x_0, film_x_1), 0., 0.);
+            let u: [Float; 2] = [radical_inverse(0, i), radical_inverse(1, i)];
+            let p_rear = Point3f::new(lerp(u[0], proj_rear_bounds.p_min.x(), proj_rear_bounds.p_max.x()),
+                          lerp(u[1], proj_rear_bounds.p_min.y(), proj_rear_bounds.p_max.y()),
+                          lens_rear_z());
+
+            // 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, None), None) {
+                pupil_bounds = pupil_bounds.union_point(Point2f::new(p_rear.x(), p_rear.y()));
+            }
+        }
+
+        // Return degenerate bounds if no rays made it through the lens system
+        if pupil_bounds.is_degenerate() {
+            print!("Unable to find exit pupil in x = {},{} on film.", film_x_0, film_x_1);
+            return pupil_bounds;
+        }
+
+        pupil_bounds.expand(2. * proj_rear_bounds.diagonal().norm() / (n_samples as Float).sqrt())
+    }    
+
+    pub fn intersect_spherical_element(radius: Float, z_center: Float, ray: &Ray) -> Option<(Float, Normal3f)> {
+        let o = ray.o - Vector3f::new(0.0, 0.0, z_center);
+
+        let a = ray.d.norm_squared();
+        let b = 2.0 * ray.d.dot(o);
+        let c = Vector3f::from(o).norm_squared() - radius * radius;
+
+        let (t0, t1) = quadratic(a, b, c)?;
+
+        let use_closer_t = (ray.d.z() > 0.0) ^ (radius < 0.0);
+        let t = if use_closer_t { t0.min(t1) } else { t0.max(t1) };
+
+        // Intersection is behind the ray's origin.
+        if t < 0.0 {
+            return None;
+        }
+
+        let p_hit_relative = o + ray.d * t;
+        // Ensures the normal points towards the incident ray.
+        let n = Normal3f::from(Vector3f::from(p_hit_relative)).normalize().face_forward(-ray.d);
+
+        Some((t, n))
+    }
+
+    pub fn trace_lenses_from_film(&self, r_camera: &Ray) -> Option<(Float, Ray)> {
+        let mut element_z = 0.; 
+        let weight = 1.;
+        // Transform _r_camera_ from camera to lens system space
+        let mut r_lens = Ray::new(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), None);
+
+        for i in (0..self.element_interface.len() - 1).rev() {
+            let element: &LensElementInterface = &self.element_interface[i];
+            // Update ray from film accounting for interaction with _element_
+            element_z -= element.thickness;
+
+            let is_stop = element.curvature_radius == 0.;
+            let t: Float;
+            let mut n = Normal3f::default();
+            if is_stop {
+                // Compute _t_ at plane of aperture stop
+                t = (element_z - r_lens.o.z()) / r_lens.d.z();
+            } else {
+                // Intersect ray with element to compute _t_ and _n_
+                let radius = element.curvature_radius;
+                let z_center = element_z + element.curvature_radius;
+                if let Some((intersect_t, intersect_n)) = RealisticCamera::intersect_spherical_element(radius, z_center, &r_lens) {
+                    t = intersect_t;
+                    n = intersect_n;
+                } else {
+                    return None; // Ray missed the element
+                }
+            }
+
+            if t < 0. {
+                return None;
+            }
+
+            // Test intersection point against element aperture
+            let p_hit = r_lens.evaluate(t);
+            if square(p_hit.x()) + square(p_hit.y()) > square(element.aperture_radius) {
+                return None;
+            }
+
+            r_lens.o = p_hit;
+
+            // Update ray path for element interface interaction
+            if !is_stop {
+                let eta_i = element.eta;
+                let eta_t = if i > 0 && self.element_interface[i - 1].eta != 0. { self.element_interface[i - 1].eta } else { 1. };
+                let wi = -r_lens.d.normalize();
+                let eta_ratio = eta_t / eta_i;
+
+                // Handle refraction idiomatically
+                if let Some((wt, _final_eta)) = refract(wi, n, eta_ratio) {
+                    r_lens.d = wt;
+                } else {
+                    // Total internal reflection occurred
+                    return None;
+                }
+            }
+        }
+
+        // Transform lens system space ray back to camera space
+        let r_out = Ray::new(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), None);
+
+        Some((weight, r_out))
+    }
+
+    pub fn sample_exit_pupil(&self, p_film: Point2f, u_lens: Point2f) -> Option<ExitPupilSample> {
+        // Find exit pupil bound for sample distance from film center
+        let r_film = (square(p_film.x()) + square(p_film.y())).sqrt();
+        let mut r_index = (r_film / (self.base.film.diagonal() / 2.)) as usize * self.exit_pupil_bounds.len();
+        r_index = (self.exit_pupil_bounds.len() - 1).min(r_index);
+
+        let pupil_bounds = self.exit_pupil_bounds[r_index];
+        if pupil_bounds.is_degenerate() {
+            return None;
+        }
+
+        // Generate sample point inside exit pupil bound
+        let p_lens = pupil_bounds.lerp(u_lens);
+        let pdf = 1. / pupil_bounds.area();
+
+        // Return sample point rotated by angle of _pFilm_ with $+x$ axis
+        let sin_theta = if r_film != 0. { p_film.y() / r_film } else { 0. };
+        let cos_theta = if r_film != 0. { p_film.x() / r_film } else { 1. };
+        let p_pupil = Point3f::new(cos_theta * p_lens.x() - sin_theta * p_lens.y(),
+                       sin_theta * p_lens.x() + cos_theta * p_lens.y(), self.lens_rear_z());
+
+        Some(ExitPupilSample{p_pupil, pdf})
+    }
+}
+
+impl CameraTrait for RealisticCamera {
+    fn base(&self) -> &CameraBase {
+        &self.base
+    }
+
+    fn generate_ray(&self, sample: CameraSample, _lambda: &SampledWavelengths) -> Option<CameraRay> {
+        // Find point on film, _pFilm_, corresponding to _sample.pFilm_
+        let s = Point2f::new(sample.p_film.x() / self.base.film.full_resolution().x() as Float, 
+            sample.p_film.y() / self.base.film.full_resolution().y() as Float);
+        let p_film2 = self.physical_extent.lerp(s);
+        let p_film = Point3f::new(-p_film2.x(), p_film2.y(), 0.);
+
+        // Trace ray from _pFilm_ through lens system
+        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, None);
+        let (weight, mut ray) = self.trace_lenses_from_film(&r_film)?;
+        if weight == 0. {
+            return None
+        }
+
+        // Finish initialization of _RealisticCamera_ ray
+        ray.time = self.sample_time(sample.time);
+        ray.medium = self.base.medium.clone();
+        ray = self.render_from_camera(&ray, &mut None);
+        ray.d = ray.d.normalize();
+
+        // Compute weighting for _RealisticCamera_ ray
+        let cos_theta = r_film.d.normalize().z();
+        let final_weight = weight * cos_theta.powf(4.) / (eps.pdf as Float * square(self.lens_rear_z()));
+
+        Some(CameraRay{ray, weight: SampledSpectrum::new(final_weight)})
+    }
+}

src/camera/spherical.rs

@@ -0,0 +1,51 @@
+use super::{CameraRay, CameraTrait, CameraBase};
+use crate::utils::geometry::{Bounds2f, Point2f, Point3f, Vector3f, Ray, spherical_direction};
+use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+use crate::utils::math::{wrap_equal_area_square, equal_area_square_to_sphere};
+use crate::core::film::FilmTrait;
+use crate::core::pbrt::{Float, PI};
+use crate::core::sampler::CameraSample;
+
+#[derive(PartialEq)]
+pub struct EquiRectangularMapping;
+
+#[derive(PartialEq)]
+pub enum Mapping {
+    EquiRectangular(EquiRectangularMapping),
+}
+
+pub struct SphericalCamera {
+    pub base: CameraBase,
+    pub screen: Bounds2f,
+    pub lens_radius: Float, 
+    pub focal_distance: Float,
+    pub mapping: Mapping,
+}
+
+impl CameraTrait for SphericalCamera {
+    fn base(&self) -> &CameraBase {
+        &self.base
+    }
+    
+    fn generate_ray(&self, sample: CameraSample, _lamdba: &SampledWavelengths) -> Option<CameraRay> {
+        // Compute spherical camera ray direction
+        let mut uv = Point2f::new(sample.p_film.x() / self.base().film.full_resolution().x() as Float,
+                   sample.p_film.y() / self.base().film.full_resolution().y() as Float);
+        let dir: Vector3f;
+        if self.mapping == Mapping::EquiRectangular(EquiRectangularMapping) {
+            // Compute ray direction using equirectangular mapping
+            let theta = PI * uv[1];
+            let phi = 2. * PI * uv[0];
+            dir = spherical_direction(theta.sin(), theta.cos(), phi);
+
+        } else {
+            // Compute ray direction using equal area mapping
+            uv = wrap_equal_area_square(&mut uv);
+            dir = equal_area_square_to_sphere(uv);
+        }
+        std::mem::swap(&mut dir.y(), &mut dir.z());
+
+        let ray = Ray::new(Point3f::new(0., 0., 0.), dir, Some(self.sample_time(sample.time)), self.base().medium.clone());
+        Some(CameraRay{ray: self.render_from_camera(&ray, &mut None), weight: SampledSpectrum::default() })
+    }
+}

src/core/bxdf.rs

@@ -0,0 +1,213 @@
+use bitflags::bitflags;
+use std::fmt;
+use std::ops::Not;
+
+use crate::utils::geometry::{Point2f, Vector3f, abs_cos_theta};
+use crate::utils::spectrum::SampledSpectrum;
+use crate::core::pbrt::{Float, PI};
+use crate::utils::sampling::{uniform_hemisphere_pdf, sample_uniform_hemisphere};
+
+bitflags! {
+    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
+    pub struct BxDFReflTransFlags: u8 {
+        const REFLECTION   = 1 << 0;
+        const TRANSMISSION = 1 << 1;
+        const ALL = Self::REFLECTION.bits() | Self::TRANSMISSION.bits();
+    }
+}
+
+bitflags! {
+    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
+    pub struct BxDFFlags: u8 {
+        const REFLECTION   = 1 << 0;
+        const TRANSMISSION = 1 << 1;
+        const DIFFUSE      = 1 << 2;
+        const GLOSSY       = 1 << 3;
+        const SPECULAR     = 1 << 4;
+
+        // Composite BxDFFlags definitions
+        const DIFFUSE_REFLECTION     = Self::DIFFUSE.bits() | Self::REFLECTION.bits();
+        const DIFFUSE_TRANSMISSION   = Self::DIFFUSE.bits() | Self::TRANSMISSION.bits();
+        const GLOSSY_REFLECTION      = Self::GLOSSY.bits() | Self::REFLECTION.bits();
+        const GLOSSY_TRANSMISSION    = Self::GLOSSY.bits() | Self::TRANSMISSION.bits();
+        const SPECULAR_REFLECTION    = Self::SPECULAR.bits() | Self::REFLECTION.bits();
+        const SPECULAR_TRANSMISSION  = Self::SPECULAR.bits() | Self::TRANSMISSION.bits();
+        
+        const SCATTERING = Self::DIFFUSE.bits() | Self::GLOSSY.bits() | Self::SPECULAR.bits();
+        const ALL = Self::REFLECTION.bits() | Self::TRANSMISSION.bits() | Self::SCATTERING.bits();
+    }
+}
+
+impl BxDFFlags {
+    #[inline]
+    pub fn is_reflective(&self) -> bool { self.contains(Self::REFLECTION) }
+    #[inline]
+    pub fn is_transmissive(&self) -> bool { self.contains(Self::TRANSMISSION) }
+    #[inline]
+    pub fn is_diffuse(&self) -> bool { self.contains(Self::DIFFUSE) }
+    #[inline]
+    pub fn is_glossy(&self) -> bool { self.contains(Self::GLOSSY) }
+    #[inline]
+    pub fn is_specular(&self) -> bool { self.contains(Self::SPECULAR) }
+    #[inline]
+    pub fn is_non_specular(&self) -> bool { self.intersects(Self::DIFFUSE | Self::GLOSSY) }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum TransportMode {
+    Radiance,
+    Importance,
+}
+
+impl Not for TransportMode {
+    type Output = Self;
+    fn not(self) -> Self::Output {
+        match self {
+            TransportMode::Radiance => TransportMode::Importance,
+            TransportMode::Importance => TransportMode::Radiance,
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct BSDFSample {
+    pub f: SampledSpectrum,
+    pub wi: Vector3f,
+    pub pdf: Float,
+    pub flags: BxDFFlags,
+    pub eta: Float,
+    pub pdf_is_proportional: bool,
+}
+
+impl BSDFSample {
+    pub fn new(
+        f: SampledSpectrum, wi: Vector3f, pdf: Float, flags: BxDFFlags,
+        eta: Float, pdf_is_proportional: bool
+    ) -> Self {
+        Self { f, wi, pdf, flags, eta, pdf_is_proportional }
+    }
+
+    pub fn default() -> Self {
+        Self {
+            f: SampledSpectrum::default(),
+            wi: Vector3f::default(),
+            pdf: 0.0,
+            flags: BxDFFlags::empty(),
+            eta: 1.0,
+            pdf_is_proportional: false,
+        }
+    }
+    
+    #[inline] pub fn is_reflective(&self) -> bool { self.flags.is_reflective() }
+    #[inline] pub fn is_transmissive(&self) -> bool { self.flags.is_transmissive() }
+    #[inline] pub fn is_diffuse(&self) -> bool { self.flags.is_diffuse() }
+    #[inline] pub fn is_glossy(&self) -> bool { self.flags.is_glossy() }
+    #[inline] pub fn is_specular(&self) -> bool { self.flags.is_specular() }
+}
+
+pub struct DiffuseBxDF;
+pub struct DiffuseTransmissionBxDF;
+pub struct DielectricBxDF;
+pub struct ThinDielectricBxDF;
+pub struct CoatedDiffuseBxDF;
+pub struct CoatedConductorBxDF;
+pub struct HairBxDF;
+pub struct MeasuredBxDF;
+pub struct ConductorBxDF;
+
+pub trait BxDFTrait {
+    fn flags(&self) -> BxDFFlags;
+
+    fn f(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode) -> SampledSpectrum;
+    
+    fn sample_f(
+        &self, wo: Vector3f, uc: Float, u: Point2f, mode: Option<TransportMode>,
+        sample_flags: Option<BxDFReflTransFlags>
+    ) -> Option<BSDFSample>;
+
+    fn pdf(
+        &self, wo: Vector3f, wi: Vector3f, mode: TransportMode,
+        sample_flags: BxDFReflTransFlags
+    ) -> Float;
+    
+    fn rho_wo(&self, wo: Vector3f, uc: &[Float], u2: &[Point2f]) -> SampledSpectrum {
+        let mut r = SampledSpectrum::new(0.);
+        for i in 0..uc.len() {
+            if let Some(bs) = self.sample_f(wo, uc[i], u2[i], Some(TransportMode::Radiance), Some(BxDFReflTransFlags::ALL)) {
+                r += bs.f * abs_cos_theta(bs.wi) / bs.pdf;
+            }
+        }
+        r / uc.len() as Float
+    }
+
+    fn rho_u(&self, u1: &[Point2f], uc: &[Float], u2: &[Point2f]) -> SampledSpectrum {
+        let mut r = SampledSpectrum::new(0.);
+        for i in 0..uc.len() {
+            let wo = sample_uniform_hemisphere(u1[i]);
+            if wo.z() == 0. {
+                continue
+            }
+            let pdfo = uniform_hemisphere_pdf();
+            if let Some(bs) = self.sample_f(wo, uc[i], u2[i], Some(TransportMode::Radiance), Some(BxDFReflTransFlags::ALL)) {
+                r += bs.f * abs_cos_theta(bs.wi) * abs_cos_theta(wo) / (pdfo * bs.pdf);
+            }
+        }
+        r / (PI * uc.len() as Float)
+    }
+    
+    fn regularize(&mut self) {
+        todo!();
+    }
+}
+
+pub enum BxDF {
+    Diffuse(DiffuseBxDF),
+    DiffuseTransmission(DiffuseTransmissionBxDF),
+    Dielectric(DielectricBxDF),
+    // ThinDielectric(ThinDielectricBxDF),
+    // Hair(HairBxDF),
+    // Measured(MeasuredBxDF),
+    // Conductor(ConductorBxDF),
+}
+
+impl BxDFTrait for BxDF {
+    fn flags(&self) -> BxDFFlags {
+        match self {
+            BxDF::Diffuse(b) => b.flags(),
+            BxDF::DiffuseTransmission(b) => b.flags(),
+            BxDF::Dielectric(b) => b.flags(),
+        }
+    }
+    
+    fn f(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode) -> SampledSpectrum {
+        match self {
+            BxDF::Diffuse(b) => b.f(wo, wi, mode),
+            BxDF::DiffuseTransmission(b) => b.f(wo, wi, mode),
+            BxDF::Dielectric(b) => b.f(wo, wi, mode),
+        }
+    }
+    
+    fn sample_f(&self, wo: Vector3f, uc: Float, u: Point2f, mode: TransportMode, sample_flags: BxDFReflTransFlags) -> Option<BSDFSample> {
+        match self {
+            BxDF::Diffuse(b) => b.sample_f(wo, uc, u, mode, sample_flags),
+            BxDF::DiffuseTransmission(b) => b.sample_f(wo, uc, u, mode, sample_flags),
+            BxDF::Dielectric(b) => b.sample_f(wo, uc, u, mode, sample_flags),
+        }
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode, sample_flags: BxDFReflTransFlags) -> Float {
+        match self {
+            BxDF::Diffuse(b) => b.pdf(wo, wi, mode, sample_flags),
+            BxDF::DiffuseTransmission(b) => b.pdf(wo, wi, mode, sample_flags),
+            BxDF::Dielectric(b) => b.pdf(wo, wi, mode, sample_flags),
+        }
+    }
+    
+    fn regularize(&mut self) {
+        match self {
+            BxDF::Diffuse(b) => b.regularize(),
+            BxDF::DiffuseTransmission(b) => b.regularize(),
+            BxDF::Dielectric(b) => b.regularize(),
+        }
+    }
+}

src/core/camera.rs

@@ -1,76 +0,0 @@
-use crate::core::film::Film;
-use crate::core::geometry::{Point2f, Point3f, Ray, RayDifferential};
-use crate::core::medium::Medium;
-use crate::core::pbrt::{Float, RenderingCoordinateSystem};
-use crate::core::spectrum::SampledSpectrum;
-use crate::core::transform::{Transform, AnimatedTransform};
-
-pub enum Camera {
-    Perspective(PerspectiveCamera),
-    Ortographic(OrtographicCamera),
-    Spherical(SphericalCamera),
-    Realistic(RealisticCamera),
-}
-
-impl Camera {
-    pub fn create(&self, name: str, medium: Medium, camera_transform: CameraTransform, film: Film) -> Float {
-        match self {
-            Camera::Perspective(s) => s.create(name, medium, camera_transform, film),
-            Camera::Ortographic(s) => s.create(name, medium, camera_transform, film),
-            Camera::Spherical(s) => s.create(name, medium, camera_transform, film),
-            Camera::Realistic(s) => s.create(name, medium, camera_transform, film),
-        }
-    }
-}
-
-pub struct CameraSample {
-    p_film: Point2f,
-    p_lens: Point2f,
-    time: Float,
-    filter_weight: Float,
-}
-
-pub struct CameraRay {
-    ray: Ray,
-    weight: SampledSpectrum,
-}
-
-pub struct CameraDifferential {
-    ray: RayDifferential,
-    weight: SampledSpectrum,
-}
-
-pub struct CameraTransform {
-    render_from_camera: AnimatedTransform,
-    world_from_render: Transform<Float>,
-}
-
-impl CameraTransform {
-    pub fn from_world(world_from_camera: AnimatedTransform, rendering_space: RenderingCoordinateSystem) -> Self {
-        let world_from_render = match rendering_space {
-            RenderingCoordinateSystem::Camera => {
-                let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.;
-                world_from_camera.interpolate(t_mid);
-            }
-            RenderingCoordinateSystem::CameraWorld => {
-                let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.;
-                let p_camera = world_from_camera.apply_point(Point3f::default(), t_mid);
-                Transform::translate(p_camera.to_vec());
-            }
-            RenderingCoordinateSystem::World => Transform::identity(),
-        };
-    }
-
-    pub fn render_from_camera(&self, p: Point3f, time: Float) -> Point3f {
-        self.render_from_camera.apply_point(p, time)
-    }
-
-    pub fn camera_from_render(&self, p: Point3f, time: Float) -> Point3f {
-        self.render_from_camera.apply_inverse_point(p, time)
-    }
-}
-
-pub struct PerspectiveCamera;
-pub struct OrtographicCamera;
-pub struct SphericalCamera;
-pub struct RealisticCamera;

src/core/film.rs

@@ -1,9 +1,110 @@
+use crate::core::pbrt::Float;
+use crate::utils::color::RGB;
+use crate::core::filter::Filter;
+use crate::utils::geometry::{Point2i, Point2f, Bounds2f, Bounds2fi};
+use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+
+
+pub struct RGBFilm {
+    pub base: FilmBase,
+}
+pub struct GBufferBFilm {
+    pub base: FilmBase,
+}
+pub struct SpectralFilm {
+    pub base: FilmBase,
+}
+
+pub struct VisibleSurface;
+pub struct PixelSensor;
+pub struct ImageMetadata;
+
+pub struct FilmBase {
+    pub full_resolution: Point2i,
+    pub pixel_bounds: Bounds2fi,
+    pub filter: Filter,
+    pub diagonal: Float,
+    pub sensor: Option<PixelSensor>,
+    pub filename: String,
+}
+
+pub trait FilmTrait {
+    // Must be implemented
+    fn base(&self) -> &FilmBase;
+    fn base_mut(&mut self) -> &mut FilmBase;
+    fn add_sample(&mut self, 
+        _p_filme: Point2i, 
+        _l: SampledSpectrum, 
+        _lambda: &SampledWavelengths, 
+        _visible_surface: Option<&VisibleSurface>, 
+        _weight: Float) {todo!()}
+    fn add_splat(&mut self, _p: Point2f, _v: SampledSpectrum, _lambda: &SampledWavelengths) { todo!() }
+    fn write_image(&self, _splat_scale: Float) { todo!() }
+    fn get_image(&self, _metadata: &ImageMetadata, _splat_scale: Float) { todo!() }
+    fn get_pixel_rgb(&self, _p: Point2i) -> RGB { todo!() }
+    fn reset_pixel(&mut self, _p: Point2i) { todo!() }
+    fn to_output_rgb(&self, _v: SampledSpectrum, _lambda: &SampledWavelengths) -> RGB { todo!() }
+    fn sample_wavelengths(&self) -> &SampledWavelengths { todo!() }
+    fn uses_visible_surface(&self) -> bool { todo!() }
+
+    // Sensible defaults
+    fn full_resolution(&self) -> Point2i { self.base().full_resolution }
+    fn sample_bounds(&self) -> Bounds2f { Bounds2f::default() }
+    fn diagonal(&self) -> Float { self.base().diagonal }
+    fn get_filter(&self) -> &Filter { &self.base().filter }
+    fn get_pixel_sensor(&self) -> Option<&PixelSensor> { self.base().sensor.as_ref() }
+} 
+
 pub enum Film {
     RGB(RGBFilm),
     GBuffer(GBufferBFilm),
     Spectral(SpectralFilm),
 }
 
-pub struct RGBFilm;
-pub struct GBufferBFilm;
-pub struct SpectralFilm;
+impl FilmTrait for RGBFilm {
+    fn base(&self) -> &FilmBase {
+        &self.base
+    }
+
+    fn base_mut(&mut self) -> &mut FilmBase {
+        &mut self.base
+    }
+}
+
+impl FilmTrait for GBufferBFilm {
+    fn base(&self) -> &FilmBase {
+        &self.base
+    }
+
+    fn base_mut(&mut self) -> &mut FilmBase {
+        &mut self.base
+    }
+}
+
+impl FilmTrait for SpectralFilm {
+    fn base(&self) -> &FilmBase {
+        &self.base
+    }
+
+    fn base_mut(&mut self) -> &mut FilmBase {
+        &mut self.base
+    }
+}
+
+impl FilmTrait for Film {
+    fn base(&self) -> &FilmBase {
+        match self {
+            Film::RGB(film) => film.base(),
+            Film::GBuffer(film) => film.base(),
+            Film::Spectral(film) => film.base(),
+        }
+    }
+    
+    fn base_mut(&mut self) -> &mut FilmBase {
+        match self {
+            Film::RGB(film) => film.base_mut(),
+            Film::GBuffer(film) => film.base_mut(),
+            Film::Spectral(film) => film.base_mut(),
+        }
+    }
+}

src/core/filter.rs

@@ -0,0 +1,344 @@
+use crate::utils::geometry::{Vector2f, Point2f, Point2i, Bounds2i, Bounds2f};
+use crate::core::pbrt::{Float, lerp};
+use crate::utils::math::{gaussian, gaussian_integral, sample_tent, windowed_sinc};
+use crate::core::sampler::PiecewiseConstant2D;
+use crate::utils::containers::Array2D;
+
+use std::hash::Hash;
+use rand::Rng;
+
+pub struct FilterSample {
+    p: Point2f,
+    weight: Float,
+}
+pub struct FilterSampler {
+    domain: Bounds2f,
+    distrib: PiecewiseConstant2D,
+    f: Array2D<Float>,
+}
+
+impl FilterSampler {
+    /// A redesigned constructor that takes a closure to avoid initialization issues.
+    pub fn new<F>(radius: Vector2f, resolution: Point2i, evaluate_fn: F) -> Self
+    where
+        F: Fn(Point2f) -> Float,
+    {
+        let domain = Bounds2f::from_points(Point2f::new(-radius.x(), -radius.y()), Point2f::new(radius.x(), radius.y()));
+        let array_bounds = Bounds2i::from_points(Point2i::new(0, 0), resolution);
+        let mut f = Array2D::new(array_bounds);
+        for j in 0..resolution.y() {
+            for i in 0..resolution.x() {
+                let p = domain.lerp(Point2f::new(
+                    (i as Float + 0.5) / resolution.x() as Float,
+                    (j as Float + 0.5) / resolution.y() as Float,
+                ));
+                f[Point2i::new(i, j)] = evaluate_fn(p);
+            }
+        }
+        let distrib = PiecewiseConstant2D::new(&f, domain);
+        Self { domain, f, distrib }
+    }
+
+    /// Samples the filter's distribution.
+    pub fn sample(&self, u: Point2f) -> FilterSample {
+        let (p, pdf, pi) = self.distrib.sample(u);
+        if pdf == 0.0 {
+            return FilterSample { p, weight: 0.0 };
+        }
+        let weight = self.f[pi] / pdf;
+        FilterSample { p, weight }
+    }
+}
+
+pub trait FilterTrait {
+    fn radius(&self) -> Vector2f;
+    fn evaluate(&self, p: Point2f) -> Float;
+    fn integral(&self) -> Float;
+    fn sample(&self, u: Point2f) -> FilterSample;
+}
+
+pub enum Filter {
+    Box(BoxFilter),
+    Gaussian(GaussianFilter),
+    Mitchell(MitchellFilter),
+    LanczosSinc(LanczosSincFilter),
+    Triangle(TriangleFilter),
+}
+
+impl FilterTrait for Filter {
+    fn radius(&self) -> Vector2f {
+        match self {
+            Filter::Box(c) => c.radius(),
+            Filter::Gaussian(c) => c.radius(),
+            Filter::Mitchell(c) => c.radius(),
+            Filter::LanczosSinc(c) => c.radius(),
+            Filter::Triangle(c) => c.radius(),
+        }
+    }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        match self {
+            Filter::Box(c) => c.evaluate(p),
+            Filter::Gaussian(c) => c.evaluate(p),
+            Filter::Mitchell(c) => c.evaluate(p),
+            Filter::LanczosSinc(c) => c.evaluate(p),
+            Filter::Triangle(c) => c.evaluate(p),
+        }
+    }
+
+    fn integral(&self) -> Float {
+        match self {
+            Filter::Box(c) => c.integral(),
+            Filter::Gaussian(c) => c.integral(),
+            Filter::Mitchell(c) => c.integral(),
+            Filter::LanczosSinc(c) => c.integral(),
+            Filter::Triangle(c) => c.integral(),
+        }
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        match self {
+            Filter::Box(c) => c.sample(u),
+            Filter::Gaussian(c) => c.sample(u),
+            Filter::Mitchell(c) => c.sample(u),
+            Filter::LanczosSinc(c) => c.sample(u),
+            Filter::Triangle(c) => c.sample(u),
+        }
+    }
+}
+
+pub struct BoxFilter {
+    pub radius: Vector2f,
+}
+
+impl BoxFilter {
+    pub fn new(radius: Vector2f) -> Self {
+        Self { radius }
+    }
+}
+
+impl FilterTrait for BoxFilter {
+    fn radius(&self) -> Vector2f {
+        self.radius
+    }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        if p.x().abs() <= self.radius.x() && p.y().abs() <= self.radius.y() {
+            1.
+        } else {
+            0.
+        }
+    }
+
+    fn integral(&self) -> Float {
+        (2.0 * self.radius.x()) * (2.0 * self.radius.y())
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        let p = Point2f::new(lerp(u[0], -self.radius.x(), self.radius.x()), lerp(u[1], -self.radius.y(), self.radius.y()));
+        FilterSample { p,  weight: 1.0 }
+    }
+}
+
+pub struct GaussianFilter {
+    pub radius: Vector2f,
+    pub sigma: Float,
+    pub exp_x: Float,
+    pub exp_y: Float,
+    pub sampler: FilterSampler,
+}
+
+impl GaussianFilter {
+    pub fn new(radius: Vector2f, sigma: Float) -> Self {
+        let exp_x = gaussian(radius.x(), 0., sigma);
+        let exp_y = gaussian(radius.y(), 0., sigma);
+
+        let sampler = FilterSampler::new(
+            radius,
+            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
+            |p: Point2f| {
+                (gaussian(p.x(), 0., sigma) - exp_x).max(0.) *
+                    (gaussian(p.y(), 0., sigma) - exp_y).max(0.)
+            }
+        );
+        Self {
+            radius,
+            sigma,
+            exp_x: gaussian(radius.x(), 0., sigma),
+            exp_y: gaussian(radius.y(), 0., sigma),
+            sampler
+        }
+    }
+}
+
+impl FilterTrait for GaussianFilter {
+    fn radius(&self) -> Vector2f {
+        self.radius
+    }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        (gaussian(p.x(), 0.0, self.sigma) - self.exp_x).max(0.0) *
+        (gaussian(p.y(), 0.0, self.sigma) - self.exp_y).max(0.0)
+    }
+
+    fn integral(&self) -> Float {
+        (gaussian_integral(-self.radius.x(), self.radius.x(), 0.0, self.sigma) - 2.0 * self.radius.x() * self.exp_x) *
+        (gaussian_integral(-self.radius.y(), self.radius.y(), 0.0, self.sigma) - 2.0 * self.radius.y() * self.exp_y)
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        self.sampler.sample(u)
+    }
+}
+
+pub struct MitchellFilter {
+    pub radius: Vector2f,
+    pub b: Float,
+    pub c: Float,
+    pub sampler: FilterSampler,
+}
+
+impl MitchellFilter {
+    pub fn new(radius: Vector2f, b: Float, c: Float) -> Self {
+        let sampler = FilterSampler::new(
+            radius,
+            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
+            move |p: Point2f| {
+                let mitchell_1d = |x: Float| {
+                    let x = x.abs();
+                    if x <= 1.0 {
+                        ((12.0 - 9.0 * b - 6.0 * c) * x.powi(3) +
+                         (-18.0 + 12.0 * b + 6.0 * c) * x.powi(2) +
+                         (6.0 - 2.0 * b)) * (1.0 / 6.0)
+                    } else if x <= 2.0 {
+                        ((-b - 6.0 * c) * x.powi(3) +
+                         (6.0 * b + 30.0 * c) * x.powi(2) +
+                         (-12.0 * b - 48.0 * c) * x +
+                         (8.0 * b + 24.0 * c)) * (1.0 / 6.0)
+                    } else {
+                        0.0
+                    }
+                };
+                mitchell_1d(2.0 * p.x() / radius.x()) * mitchell_1d(2.0 * p.y() / radius.y())
+            },
+        );
+        Self { radius, b, c, sampler }
+    }
+    
+    fn mitchell_1d(&self, x: Float) -> Float {
+        let x = x.abs();
+        if x <= 1.0 {
+            ((12.0 - 9.0 * self.b - 6.0 * self.c) * x.powi(3) +
+                (-18.0 + 12.0 * self.b + 6.0 * self.c) * x.powi(2) +
+                (6.0 - 2.0 * self.b)) * (1.0 / 6.0)
+        } else if x <= 2.0 {
+            ((-self.b - 6.0 * self.c) * x.powi(3) +
+                (6.0 * self.b + 30.0 * self.c) * x.powi(2) +
+                (-12.0 * self.b - 48.0 * self.c) * x +
+                (8.0 * self.b + 24.0 * self.c)) * (1.0 / 6.0)
+        } else {
+            0.0
+        }
+    }
+}
+
+impl FilterTrait for MitchellFilter {
+    fn radius(&self) -> Vector2f { self.radius }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        self.mitchell_1d(2.0 * p.x() / self.radius.x()) * self.mitchell_1d(2.0 * p.y() / self.radius.y())
+    }
+    
+    fn integral(&self) -> Float {
+        self.radius.x() * self.radius.y() / 4.0
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        self.sampler.sample(u)
+    }
+}
+
+pub struct LanczosSincFilter {
+    pub radius: Vector2f,
+    pub tau: Float,
+    pub sampler: FilterSampler,
+}
+
+impl LanczosSincFilter {
+    pub fn new(radius: Vector2f, tau: Float) -> Self {
+        let sampler = FilterSampler::new(
+            radius,
+            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
+            move |p: Point2f| {
+                windowed_sinc(p.x(), radius.x(), tau) * windowed_sinc(p.y(), radius.y(), tau)
+            },
+        );
+        Self { radius, tau, sampler }
+    }
+}
+
+impl FilterTrait for LanczosSincFilter {
+    fn radius(&self) -> Vector2f { self.radius }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        windowed_sinc(p.x(), self.radius.x(), self.tau) * windowed_sinc(p.y(), self.radius.y(), self.tau)
+    }
+
+    fn integral(&self) -> Float {
+        let sqrt_samples = 64;
+        let n_samples = sqrt_samples * sqrt_samples;
+        let area = (2.0 * self.radius.x()) * (2.0 * self.radius.y());
+        let mut sum = 0.0;
+        let mut rng = rand::rng();
+
+        for y in 0..sqrt_samples {
+            for x in 0..sqrt_samples {
+                let u = Point2f::new(
+                    (x as Float + rng.random::<Float>()) / sqrt_samples as Float,
+                    (y as Float + rng.random::<Float>()) / sqrt_samples as Float,
+                );
+                let p = Point2f::new(
+                    lerp(u.x(), -self.radius.x(), self.radius.x()),
+                    lerp(u.y(), -self.radius.y(), self.radius.y()),
+                );
+                sum += self.evaluate(p);
+            }
+        }
+        sum / n_samples as Float * area
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        self.sampler.sample(u)
+    }
+}
+
+pub struct TriangleFilter {
+    pub radius: Vector2f,
+}
+
+impl TriangleFilter {
+    pub fn new(radius: Vector2f) -> Self {
+        Self { radius }
+    }
+}
+
+impl FilterTrait for TriangleFilter {
+    fn radius(&self) -> Vector2f { self.radius }
+
+    fn evaluate(&self, p: Point2f) -> Float {
+        (self.radius.x() - p.x().abs()).max(0.0) *
+        (self.radius.y() - p.y().abs()).max(0.0)
+    }
+
+    fn integral(&self) -> Float {
+        self.radius.x().powi(2) * self.radius.y().powi(2)
+    }
+
+    fn sample(&self, u: Point2f) -> FilterSample {
+        let p = Point2f::new(
+            sample_tent(u[0], self.radius.x()),
+            sample_tent(u[1], self.radius.y()),
+        );
+        FilterSample { p, weight: 1.0 }
+    }
+}

src/core/interaction.rs

@@ -1,7 +1,7 @@
-use crate::core::geometry::{Vector3f, Normal3f, Point2f, Point3f, Point3fi, Normed, Ray, Dot};
 use crate::core::pbrt::Float;
 use crate::core::medium::{Medium, MediumInterface};
 use crate::core::light::Light;
+use crate::utils::geometry::{Vector3f, Normal3f, Point2f, Point3f, Point3fi, Normed, Ray, Dot};
 
 use std::any::Any;
 use std::sync::Arc;

src/core/medium.rs

@@ -1,11 +1,17 @@
 use std::sync::Arc;
 
+#[derive(Debug, Clone)]
 pub struct HomogeneousMedium;
+#[derive(Debug, Clone)]
 pub struct GridMedium;
+#[derive(Debug, Clone)]
 pub struct RGBGridMedium;
+#[derive(Debug, Clone)]
 pub struct CloudMedium;
+#[derive(Debug, Clone)]
 pub struct NanoVDBMedium;
 
+#[derive(Debug, Clone)]
 pub enum Medium {
     Homogeneous(HomogeneousMedium),
     Grid(GridMedium),

src/core/mod.rs

@@ -1,18 +1,14 @@
-pub mod camera;
-pub mod color;
-pub mod colorspace;
+pub mod bxdf;
+pub mod cie;
 pub mod film;
-pub mod geometry;
+pub mod filter;
 pub mod integrator;
 pub mod light;
 pub mod interaction;
-pub mod interval;
 pub mod material;
 pub mod medium;
 pub mod pbrt;
 pub mod primitive;
-pub mod quaternion;
+pub mod sampler;
 pub mod shape;
 pub mod texture;
-pub mod transform;
-pub mod spectrum;

src/core/pbrt.rs

@@ -1,10 +1,14 @@
 use num_traits::Num;
 use std::ops::{Add, Mul};
 
+use crate::utils::geometry::{Point2f, Vector2f, Vector3f};
+
 pub type Float = f32;
 
 pub const MACHINE_EPSILON: Float = std::f32::EPSILON * 0.5;
 pub const SHADOW_EPSILON: Float = 0.0001;
+pub const ONE_MINUS_EPSILON: Float = 0.99999994;
+pub const PI: Float = std::f32::consts::PI;
 pub const INV_PI: Float = 0.318_309_886_183_790_671_54;
 pub const INV_2_PI: Float = 0.159_154_943_091_895_335_77;
 pub const INV_4_PI: Float = 0.079_577_471_545_947_667_88;
@@ -41,6 +45,34 @@ where
     u * (a + b)
 }
 
+pub fn sample_uniform_disk_polar(u: Point2f) -> Point2f {
+    let r = u[0].sqrt();
+    let theta = 2. * PI * u[1];
+    Point2f::new(r * theta.cos(), r * theta.sin())
+}
+
+pub fn sample_uniform_disk_concentric(u: Point2f) -> Vector2f {
+    // Map _u_ to $[-1,1]^2$ and handle degeneracy at the origin
+    let u_vector = u - Point2f::new(0., 0.);
+    let u_offset = 2. * u_vector - Vector2f::new(1., 1.);
+    if u_offset.x() == 0. && u_offset.y() == 0. {
+        return Vector2f::new(0., 0.)
+    }
+
+    // Apply concentric mapping to point
+    let theta: Float;
+    let r: Float;
+    if u_offset.x().abs() > u_offset.y().abs() {
+        r = u_offset.x();
+        theta = PI_OVER_4 * (u_offset.y() / u_offset.x());
+    } else {
+        r = u_offset.y();
+        theta = PI_OVER_2 - PI_OVER_4 * (u_offset.x() / u_offset.y());
+    }
+    let s_vector = Point2f::new(theta.cos(), theta.sin()) - Point2f::new(0., 0.);
+    return r * s_vector;
+}
+
 pub fn clamp_t<T>(val: T, low: T, high: T) -> T
 where
     T: PartialOrd,
@@ -144,29 +176,16 @@ where
 
 
 
-#[inline]
-pub fn safe_asin(x: Float) -> Float {
-    if x >= -1.0001 && x <= 1.0001 {
-        clamp_t(x, -1., 1.).asin()
-    } else {
-        panic!("Not valid value for asin")
-    }
-}
-
-
-#[inline]
-pub fn sinx_over_x(x: Float) -> Float {
-    if 1. - x * x == 1. {
-        return 1.;
-    }
-    return x.sin() / x;
-}
-
 #[inline]
 pub fn gamma(n: i32) -> Float {
     return (n as Float * MACHINE_EPSILON) / (1. - n as Float * MACHINE_EPSILON);
 }
 
+#[inline]
+pub fn square<T>(n: T) -> T
+where 
+    T: Mul<Output = T> + Copy { n * n }
+
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum RenderingCoordinateSystem {

src/core/sampler.rs

@@ -0,0 +1,128 @@
+use crate::utils::geometry::{Point2f, Point2i, Vector2f, Bounds2f};
+use crate::core::pbrt::{find_interval, Float, PI, PI_OVER_2, PI_OVER_4, lerp};
+use crate::utils::containers::Array2D;
+
+#[derive(Debug, Clone, Copy)]
+pub struct CameraSample {
+    pub p_film: Point2f,
+    pub p_lens: Point2f,
+    pub time: Float,
+    pub filter_weight: Float,
+}
+
+impl Default for CameraSample {
+    fn default() -> Self {
+        Self { p_film: Point2f::default(), p_lens: Point2f::default(), time: 0.0, filter_weight: 1.0 }
+    }
+}
+
+pub struct PiecewiseConstant1D {
+    pub func: Vec<Float>,
+    pub cdf: Vec<Float>,
+    pub min: Float,
+    pub max: Float,
+    pub func_integral: Float,
+}
+impl PiecewiseConstant1D {
+    pub fn new(f: &[Float]) -> Self {
+        Self::new_with_bounds(f, 0., 1.)
+    }
+    pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
+        assert!(max > min);
+        let n = f.len();
+        let mut func = Vec::with_capacity(n);
+        for &val in f {
+            func.push(val.abs());
+        }
+        let mut cdf = vec![0.; n + 1];
+        for i in 1..=n {
+            debug_assert!(func[i - 1] >= 0.);
+            cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n as Float;
+        }
+        let func_integral = cdf[n];
+        if func_integral == 0. {
+            for i in 1..=n {
+                cdf[i] = i as Float / n as Float;
+            }
+        } else {
+            for i in 1..=n {
+                cdf[i] /= func_integral;
+            }
+        }
+        Self { func, cdf, func_integral, min, max }
+    }
+    pub fn integral(&self) -> Float {
+        self.func_integral
+    }
+    pub fn len(&self) -> usize {
+        self.func.len()
+    }
+    pub fn sample(&self, u: Float) -> (Float, Float, usize) {
+        let o = find_interval(self.cdf.len(), |idx| self.cdf[idx] <= u);
+        let mut du = u - self.cdf[o];
+        if self.cdf[o + 1] - self.cdf[o] > 0. {
+            du /= self.cdf[o + 1] - self.cdf[o];
+        }
+        debug_assert!(!du.is_nan());
+        let value = lerp((o as Float + du) / self.len() as Float, self.min, self.max);
+        let pdf_val = if self.func_integral > 0. { self.func[o] / self.func_integral } else { 0. };
+        (value, pdf_val, o)
+    }
+}
+pub struct PiecewiseConstant2D {
+    pub p_conditional_v: Vec<PiecewiseConstant1D>,
+    pub p_marginal: PiecewiseConstant1D,
+    pub domain: Bounds2f,
+}
+
+impl PiecewiseConstant2D {
+    pub fn new(data: &Array2D<f32>, domain: Bounds2f) -> Self {
+        let resolution = data.size();
+        let nu = resolution.x as usize;
+        let nv = resolution.y as usize;
+        let mut p_conditional_v = Vec::with_capacity(nv);
+        for v in 0..nv {
+            let start = v * nu;
+            let end = start + nu;
+            p_conditional_v.push(PiecewiseConstant1D::new_with_bounds(
+                &data.as_slice()[start..end],
+                domain.p_min.x(),
+                domain.p_max.x(),
+            ));
+        }
+        // Compute marginal sampling distribution from the integrals of the conditional distributions
+        let marginal_func: Vec<f32> = p_conditional_v.iter().map(|p| p.integral()).collect();
+        let p_marginal = PiecewiseConstant1D::new_with_bounds(
+            &marginal_func,
+            domain.p_min.y(),
+            domain.p_max.y(),
+        );
+        Self { p_conditional_v, p_marginal, domain }
+    }
+   
+    pub fn integral(&self) -> f32 {
+        self.p_marginal.integral()
+    }
+    pub fn sample(&self, u: Point2f) -> (Point2f, f32, Point2i) {
+        let (d1, pdf1, off_y) = self.p_marginal.sample(u.y());
+        let (d0, pdf0, off_x) = self.p_conditional_v[off_y].sample(u.x());
+        let pdf = pdf0 * pdf1;
+        let offset = Point2i::new(off_x as i32, off_y as i32);
+        (Point2f::new(d0, d1), pdf, offset)
+    }
+    /// Calculates the PDF at a given point.
+    pub fn pdf(&self, p: Point2f) -> f32 {
+        let p_offset = self.domain.offset(&p);
+        let nu = self.p_conditional_v[0].len();
+        let nv = self.p_marginal.len();
+       
+        let iu = (p_offset.x() * nu as f32).clamp(0.0, nu as f32 - 1.0) as usize;
+        let iv = (p_offset.y() * nv as f32).clamp(0.0, nv as f32 - 1.0) as usize;
+        let integral = self.p_marginal.integral();
+        if integral == 0.0 {
+            0.0
+        } else {
+            self.p_conditional_v[iv].func[iu] / integral
+        }
+    }
+}

src/core/transform.rs

@@ -1,847 +0,0 @@
-use num_traits::{Num, One, Zero, Signed, Float as NumFloat};
-use std::ops::{Add, Div, Index, IndexMut, Mul};
-use std::fmt;
-use std::fmt::{Display, Debug};
-
-use crate::core::geometry::{Vector, Point, Vector3f, Point3f, Point3fi};
-use crate::core::color::{RGB, XYZ};
-use crate::core::pbrt::{gamma, Float};
-use crate::core::quaternion::Quaternion;
-
-#[derive(Copy, Clone)]
-pub struct SquareMatrix<T, const N: usize> {
-    pub m: [[T; N]; N],
-}
-
-impl<T, const N: usize> SquareMatrix<T, N> 
-{
-    pub fn new(data: [[T; N]; N]) -> Self {
-        Self { m: data }
-    }
-
-    fn identity() -> Self where T: Copy + Zero + One {
-        let mut m = [[T::zero(); N]; N];
-        for i in 0..N {
-            m[i][i] = T::one();
-        }
-        Self { m }
-    }
-
-
-    pub fn zero() -> Self where T: Copy + Zero {
-        SquareMatrix { m: [[T::zero(); N]; N] }
-    }
-
-    pub fn is_identity(&self) -> bool where T: Copy + Zero + One + PartialEq {
-        for i in 0..N {
-            for j in 0..N {
-                if i == j {
-                    if self.m[i][j] != T::one() {
-                        return false;
-                    }
-                } else if self.m[i][j] != T::zero() {
-                    return false;
-                }
-            }
-        }
-        true
-    }
-
-    fn trace(&self) -> T 
-        where 
-            T: Zero + Copy
-    {
-        let mut sum = T::zero();
-        for i in 0..N {
-            for j in 0..N {
-                if i == j {
-                    sum = sum + self.m[i][j];
-                }
-            }
-        }
-        sum
-    }
-
-    pub fn transpose(&self) -> Self where T: Copy + Zero {
-        let mut result = SquareMatrix::zero();
-        for i in 0..N {
-            for j in 0..N {
-                result.m[j][i] = self.m[i][j];
-            }
-        }
-        result
-    } 
-}
-
-
-impl<T, const N: usize> SquareMatrix<T, N> 
-    where 
-        T: NumFloat
-{
-    pub fn inverse(&self) -> Option<Self> {
-        if N == 0 {
-            return Some(Self::identity());
-        }
-
-        let mut mat = self.m;
-        let mut inv = Self::identity();
-
-        for i in 0..N {
-            let mut pivot_row = i;
-            for j in (i + 1)..N {
-                if mat[j][i].abs() > mat[pivot_row][i].abs() {
-                    pivot_row = j;
-                }
-            }
-
-            if pivot_row != i {
-                mat.swap(i, pivot_row);
-                inv.m.swap(i, pivot_row);
-            }
-
-            let pivot = mat[i][i];
-            if pivot.is_zero() {
-                return None;
-            }
-
-            for j in i..N {
-                mat[i][j] = mat[i][j] / pivot;
-            }
-
-            for j in 0..N {
-                inv.m[i][j] = inv.m[i][j] / pivot;
-            }
-
-            for j in 0..N {
-                if i != j {
-                    let factor = mat[j][i];
-                    for k in i..N {
-                        mat[j][k] = mat[j][k] - factor * mat[i][k];
-                    }
-                    for k in 0..N {
-                        inv.m[j][k] = inv.m[j][k] - factor * inv.m[i][k];
-                    }
-                }
-            }
-        }
-        Some(inv)
-    }
-
-    pub fn determinant(&self) -> T {
-        let m = &self.m;
-
-        match N {
-            0 => T::one(),
-            1 => m[0][0],
-            2 => m[0][0] * m[1][1] - m[0][1] * m[1][0],
-            3 => {
-                let a = m[0][0]; let b = m[0][1]; let c = m[0][2];
-                let d = m[1][0]; let e = m[1][1]; let f = m[1][2];
-                let g = m[2][0]; let h = m[2][1]; let i = m[2][2];
-
-                a * (e * i - f * h) -
-                b * (d * i - f * g) +
-                c * (d * h - e * g)
-            }
-            4 => {
-                let det3 = |m11: T, m12: T, m13: T,
-                            m21: T, m22: T, m23: T,
-                            m31: T, m32: T, m33: T| -> T {
-                    m11 * (m22 * m33 - m23 * m32) -
-                    m12 * (m21 * m33 - m23 * m31) +
-                    m13 * (m21 * m32 - m22 * m31)
-                };
-
-                let c0 = det3(m[1][1], m[1][2], m[1][3], m[2][1], m[2][2], m[2][3], m[3][1], m[3][2], m[3][3]);
-                let c1 = det3(m[1][0], m[1][2], m[1][3], m[2][0], m[2][2], m[2][3], m[3][0], m[3][2], m[3][3]);
-                let c2 = det3(m[1][0], m[1][1], m[1][3], m[2][0], m[2][1], m[2][3], m[3][0], m[3][1], m[3][3]);
-                let c3 = det3(m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2], m[3][0], m[3][1], m[3][2]);
-
-                m[0][0] * c0 - m[0][1] * c1 + m[0][2] * c2 - m[0][3] * c3
-            }
-            _ => { // Fallback to LU decomposition for N > 4
-                let mut lum = self.m;
-                let mut parity = 0;
-
-                for i in 0..N {
-                    let mut max_row = i;
-                    for k in (i + 1)..N {
-                        if lum[k][i].abs() > lum[max_row][i].abs() {
-                            max_row = k;
-                        }
-                    }
-
-                    if max_row != i {
-                        lum.swap(i, max_row);
-                        parity += 1;
-                    }
-                    
-                     // Singular matrix
-                    if lum[i][i] == T::zero() {
-                        return T::zero();
-                    }
-
-                    // Gaussian elimination
-                    for j in (i + 1)..N {
-                        let factor = lum[j][i] / lum[i][i];
-                        for k in i..N {
-                            lum[j][k] = lum[j][k] - factor * lum[i][k];
-                        }
-                    }
-                }
-
-                let mut det = T::one();
-                for i in 0..N {
-                    det = det * lum[i][i];
-                }
-
-                if parity % 2 != 0 {
-                    det = -det;
-                }
-
-                det
-            }
-        }
-    }
-}
-
-impl<T, const N: usize> Default for SquareMatrix<T, N>
-where
-    T: Copy + Zero + One,
-{
-    fn default() -> Self {
-        Self::identity()
-    }
-}
-
-impl<T, const N: usize> PartialEq for SquareMatrix<T, N>
-where
-    T: PartialEq,
-{
-    fn eq(&self, other: &Self) -> bool {
-        self.m == other.m
-    }
-}
-
-impl<T: Eq, const N: usize> Eq for SquareMatrix<T, N> {}
-
-
-impl<T, const N: usize> Index<usize> for SquareMatrix<T, N> {
-    type Output = [T; N];
-
-    fn index(&self, index: usize) -> &Self::Output {
-        &self.m[index]
-    }
-}
-
-impl<T, const N: usize> PartialOrd for SquareMatrix<T, N>
-where
-    T: PartialOrd,
-{
-    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
-        self.m.iter().flatten().partial_cmp(other.m.iter().flatten())
-    }
-}
-
-impl<T, const N: usize> IndexMut<usize> for SquareMatrix<T, N> {
-    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
-        &mut self.m[index]
-    }
-}
-
-impl<T: Debug, const N: usize> Debug for SquareMatrix<T, N> {
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        f.debug_struct("SquareMatrix")
-         .field("m", &self.m)
-         .finish()
-    }
-}
-
-impl<T, const N: usize> Display for SquareMatrix<T, N>
-where
-    T: Display + Copy,
-{
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        let mut col_widths = [0; N];
-        for row in self.m.iter() {
-            for (j, element) in row.iter().enumerate() {
-                let width = format!("{}", element).len();
-                if width > col_widths[j] {
-                    col_widths[j] = width;
-                }
-            }
-        }
-        
-        for i in 0..N {
-            write!(f, "[")?;
-            for j in 0..N {
-                write!(f, "{: >width$} ", self.m[i][j], width = col_widths[j])?;
-            }
-            writeln!(f, "]")?;
-        }
-        Ok(())
-    }
-}
-
-impl<T, const N: usize> Add for SquareMatrix<T, N>
-where 
-    T: Copy + Add<Output=T> + Zero,
-{
-    type Output = Self;
-    fn add(mut self, rhs: Self) -> Self::Output {
-        for i in 0..N {
-            for j in 0..N {
-                self.m[i][j] = self.m[i][j] + rhs.m[i][j];
-            }
-        }
-        self
-    }
-}
-
-impl<T, const N: usize> Mul for SquareMatrix<T, N>
-where
-    T: Copy + Num,
-{
-    type Output = Self;
-    fn mul(self, rhs: Self) -> Self::Output {
-        let mut result = SquareMatrix::zero();
-        for i in 0..N {
-            for j in 0..N {
-                let mut sum = T::zero();
-                for k in 0..N {
-                    sum = sum + self.m[i][k] * rhs.m[k][j];
-                }
-                result.m[i][j] = sum;
-            }
-        }
-        result
-    }
-}
-
-impl<T, const N: usize> Mul<T> for SquareMatrix<T, N>
-where
-    T: Copy + Mul<Output = T>,
-{
-    type Output = Self;
-    fn mul(mut self, rhs: T) -> Self::Output {
-        for row in self.m.iter_mut() {
-            for element in row.iter_mut() {
-                *element = *element * rhs;
-            }
-        }
-        self
-    }
-}
-
-impl<const N: usize> Mul<SquareMatrix<Float, N>> for Float
-{
-    type Output = SquareMatrix<Float, N>;
-    fn mul(self, rhs: SquareMatrix<Float, N>) -> Self::Output {
-        rhs * self
-    }
-}
-
-impl<T, const N: usize> Div<T> for SquareMatrix<T, N>
-where
-    T: Copy + Div<Output = T>,
-{
-    type Output = Self;
-    fn div(mut self, rhs: T) -> Self::Output {
-        for row in self.m.iter_mut() {
-            for element in row.iter_mut() {
-                *element = *element / rhs;
-            }
-        }
-        self
-    }
-}
-
-impl<T, const N: usize> Mul<Vector<T, N>> for SquareMatrix<T, N>
-where 
-    T: Num + Copy,
-{
-    type Output = Vector<T, N>;
-    fn mul(self, rhs: Vector<T, N>) -> Self::Output {
-        let mut result_arr = [T::zero(); N];
-        for i in 0..N {
-            for j in 0..N {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs.0[j];
-            }
-        }
-        Vector(result_arr)
-    }
-}
-
-impl<T, const N: usize> Mul<Point<T, N>> for SquareMatrix<T, N>
-where
-    T: Num + Copy,
-{
-    type Output = Point<T, N>;
-
-    fn mul(self, rhs: Point<T, N>) -> Self::Output {
-        let mut result_arr = [T::zero(); N];
-        for i in 0..N {
-            for j in 0..N {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs.0[j];
-            }
-        }
-        Point(result_arr)
-    }
-}
-
-impl SquareMatrix<Float, 3> {
-    pub fn transform_to_xyz(&self, rgb: RGB) -> XYZ {
-        let r = rgb[0];
-        let g = rgb[1];
-        let b = rgb[2];
-
-        let x = self.m[0][0] * r + self.m[0][1] * g + self.m[0][2] * b;
-        let y = self.m[1][0] * r + self.m[1][1] * g + self.m[1][2] * b;
-        let z = self.m[2][0] * r + self.m[2][1] * g + self.m[2][2] * b;
-        
-        XYZ::new(x, y, z)
-    }
-
-    // This name is also perfectly clear
-    pub fn transform_to_rgb(&self, xyz: XYZ) -> RGB {
-        let x = xyz[0];
-        let y = xyz[1];
-        let z = xyz[2];
-
-        let r = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
-        let g = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
-        let b = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
-        
-        RGB::new(r, g, b)
-    }
-}
-
-impl Mul<XYZ> for SquareMatrix<Float, 3>
-{
-    type Output = XYZ;
-    fn mul(self, rhs: XYZ) -> Self::Output {
-        let mut result_arr = [0.0; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j];
-            }
-        }
-        XYZ::new(result_arr[0], result_arr[1], result_arr[2])
-    }
-}
-
-impl Mul<RGB> for SquareMatrix<Float, 3>
-{
-    type Output = RGB;
-    fn mul(self, rhs: RGB) -> Self::Output {
-        let mut result_arr = [0.0; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j];
-            }
-        }
-        RGB::new(result_arr[0], result_arr[1], result_arr[2])
-    }
-}
-
-#[derive(Copy, Clone)]
-pub struct Transform<T: num_traits::Float> {
-    m: SquareMatrix<T, 4>,
-    m_inv: SquareMatrix<T, 4>,
-}
-
-impl<T: num_traits::Float + Signed> Transform<T> {
-    pub fn new(m: &SquareMatrix<T, 4>, m_inv: &SquareMatrix<T, 4>) -> Self {
-        Self { m: *m, m_inv: *m_inv }
-    }
-
-    pub fn from_matrix(m: &SquareMatrix<T, 4>) -> Self {
-        let inv = match m.inverse() {
-            Some(inv) => inv ,
-            None => SquareMatrix::zero(),
-        };
-
-        Self { m: *m, m_inv: inv }
-    }
-
-    pub fn identity() -> Self {
-        let m: SquareMatrix<T, 4> = SquareMatrix::identity();
-        let m_inv = m.inverse().expect("Singular matrix");
-        Self { m, m_inv }
-    }
-
-    pub fn is_identity(&self) -> bool {
-        self.m.is_identity()
-    }
-
-    pub fn inv(&self) -> Self {
-        Self { m: self.m_inv, m_inv: self.m}
-    }
-
-    pub fn apply_inverse(&self, p: Point<T, 3>) -> Point<T, 3> {
-        *self * p
-    }
-
-    pub fn apply_inverse_vector(&self, v: Vector<T, 3>) -> Vector<T, 3> {
-        *self * v
-    }
-
-}
-
-impl Transform<Float> {
-    pub fn apply(&self, pi: &Point3fi) -> Point3fi {
-        let p = pi.midpoint();
-        let x = p.x();
-        let y = p.y();
-        let z = p.z();
-        let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z + self.m[0][3];
-        let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z + self.m[1][3];
-        let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z + self.m[2][3];
-        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z + self.m[3][3];
-        let mut p_error = Vector3f::default();
-
-        if pi.is_exact() {
-            p_error[0] = gamma(3) * (self.m[0][0] * x).abs() + (self.m[0][1] * y).abs() +
-                (self.m[0][2] + z).abs() + self.m[0][3].abs();
-            p_error[1] = gamma(3) * (self.m[1][0] * x).abs() + (self.m[1][1] * y).abs() +
-                (self.m[1][2] + z).abs() + self.m[1][3].abs();
-            p_error[2] = gamma(3) * (self.m[2][0] * x).abs() + (self.m[2][1] * y).abs() +
-                (self.m[2][2] + z).abs() + self.m[2][3].abs();
-        }
-
-        if wp == 1. {
-            return Point3fi::new_with_error(Point([xp, yp, zp]), p_error);
-        } else {
-            return Point3fi::new_with_error(Point([xp/wp, yp/wp, zp/wp]), p_error)
-        }
-    }
-
-    pub fn to_quaternion(&self) -> Quaternion {
-        let trace = self.m.trace();
-        let mut quat = Quaternion::default();
-        if trace > 0. {
-            let mut s = (trace + 1.).sqrt();
-            quat.w = 0.5 * s;
-            s = 0.5 / s;
-            quat.v[0] = (self.m[2][1] - self.m[1][2]) * s;
-            quat.v[1] = (self.m[2][1] - self.m[1][2]) * s;
-            quat.v[2] = (self.m[2][1] - self.m[1][2]) * s;
-        } else {
-            let nxt = [1, 2, 0]; 
-            let mut q = [0.; 3];
-            let mut i = 0;
-            if self.m[1][1] > self.m[0][0] {
-                i = 1;
-            }
-            if self.m[2][2] > self.m[i][i] {
-                i = 2;
-            }
-            let j = nxt[i];
-            let k = nxt[j];
-            let mut s = (self.m[i][i] - (self.m[j][j] + self.m[k][k]) + 1.).sqrt();
-            q[i] = 0.5 * s;
-            if s != 0. {
-                s = 0.5 / s;
-            }
-            quat.w = (self.m[k][j] - self.m[j][k]) * s;
-            q[j] = (self.m[j][i] + self.m[i][j]) * s;
-            q[k] = (self.m[k][i] + self.m[i][k]) * s;
-            quat.v[0] = q[0];
-            quat.v[1] = q[1];
-            quat.v[2] = q[2];
-        }
-        quat
-    }
-
-    pub fn translate(delta: Vector3f) -> Self {
-        Transform {
-            m: SquareMatrix {
-                m: [
-                    [1.0, 0.0, 0.0, delta.x()],
-                    [0.0, 1.0, 0.0, delta.y()],
-                    [0.0, 0.0, 1.0, delta.z()],
-                    [0.0, 0.0, 0.0, 1.0],
-                ],
-            },
-            m_inv: SquareMatrix {
-                m: [
-                    [1.0, 0.0, 0.0, -delta.x()],
-                    [0.0, 1.0, 0.0, -delta.y()],
-                    [0.0, 0.0, 1.0, -delta.z()],
-                    [0.0, 0.0, 0.0, 1.0],
-                ],
-            },
-        }
-    }
-
-
-    pub fn scale(x: Float, y: Float, z: Float) -> Self {
-        let m = SquareMatrix::new([
-            [x, 0., 0., 0.],
-            [0., y, 0., 0.],
-            [0., 0., z, 0.],
-            [0., 0., 0., 1.],
-        ]);
-
-        let m_inv = SquareMatrix::new([
-            [1./x, 0., 0., 0.],
-            [0., 1./y, 0., 0.],
-            [0., 0., 1./z, 0.],
-            [0., 0., 0., 1.],
-        ]);
-        Self { m, m_inv }
-    }
-
-    pub fn rotate_x(theta: Float) -> Self {
-        let sin_theta = theta.to_radians().sin();
-        let cos_theta = theta.to_radians().cos();
-        let m = SquareMatrix::new([
-            [1., 0., 0., 0.],
-            [0., cos_theta, -sin_theta, 0.],
-            [0., sin_theta, cos_theta, 0.],
-            [0., 0., 0., 1.],
-        ]);
-        Self { m, m_inv: m.transpose() }
-    }
-
-    pub fn rotate_y(theta: Float) -> Self {
-        let sin_theta = theta.to_radians().sin();
-        let cos_theta = theta.to_radians().cos();
-        let m = SquareMatrix::new([
-            [cos_theta, 0., sin_theta, 0.],
-            [0., 1., 0., 0.],
-            [-sin_theta, 0., cos_theta, 0.],
-            [0., 0., 0., 1.],
-        ]);
-        Self { m, m_inv: m.transpose() }
-    }
-
-    pub fn rotate_z(theta: Float) -> Self {
-        let sin_theta = theta.to_radians().sin();
-        let cos_theta = theta.to_radians().cos();
-        let m = SquareMatrix::new([
-            [cos_theta, -sin_theta, 0., 0.],
-            [sin_theta, cos_theta, 0., 0.],
-            [0., 0., 1., 0.],
-            [0., 0., 0., 1.],
-        ]);
-        Self { m, m_inv: m.transpose() }
-    }
-}
-
-impl<T: num_traits::Float> PartialEq for Transform<T> {
-    fn eq(&self, other: &Self) -> bool {
-        self.m == other.m && self.m_inv == other.m_inv
-    }
-}
-
-impl<T> Mul for Transform<T> 
-where T: NumFloat
-{
-    type Output = Self;
-    fn mul(self, rhs: Self) -> Self::Output {
-        Transform { m: self.m * rhs.m, m_inv: rhs.m_inv * self.m_inv }
-    }
-}
-
-impl Mul<Transform<Float>> for Float
-{
-    type Output = Transform<Float>;
-    fn mul(self, rhs: Transform<Float>) -> Self::Output {
-        Transform { m: rhs.m * self , m_inv: rhs.m_inv * self }
-    }
-}
-
-impl<T> Mul<Point<T, 3>> for Transform<T>
-where 
-    T: NumFloat
-{
-    type Output = Point<T, 3>;
-    fn mul(self, p: Point<T, 3>) -> Self::Output {
-        let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3];
-        let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3];
-        let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3];
-        let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3];
-
-        if wp == T::one() {
-            Point([xp, yp, zp])
-        } else {
-            Point([xp/wp, yp/wp, zp/wp])
-        }
-    }
-}
-
-impl<T> Mul<Vector<T, 3>> for Transform<T>
-where 
-    T: NumFloat
-{
-    type Output = Vector<T, 3>;
-    fn mul(self, p: Vector<T, 3>) -> Self::Output {
-        let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3];
-        let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3];
-        let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3];
-        let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3];
-
-        if wp == T::one() {
-            Vector([xp, yp, zp])
-        } else {
-            Vector([xp/wp, yp/wp, zp/wp])
-        }
-    }
-}
-
-impl From<Quaternion> for Transform<Float> {
-    fn from(q: Quaternion) -> Self {
-        let xx = q.v.x() * q.v.x();
-        let yy = q.v.y() * q.v.y();
-        let zz = q.v.z() * q.v.z();
-        let xy = q.v.x() * q.v.y();
-        let xz = q.v.x() * q.v.z();
-        let yz = q.v.y() * q.v.z();
-        let wx = q.v.x() * q.w;
-        let wy = q.v.y() * q.w;
-        let wz = q.v.z() * q.w;
-
-        let mut m = [[0.0; 4]; 4];
-        
-        m[0][0] = 1. - 2. * (yy + zz);
-        m[0][1] = 2. * (xy - wz);
-        m[0][2] = 2. * (xz + wy);
-
-        m[1][0] = 2. * (xy + wz);
-        m[1][1] = 1. - 2. * (xx + zz);
-        m[1][2] = 2. * (yz - wx);
-
-        m[2][0] = 2. * (xz - wy);
-        m[2][1] = 2. * (yz + wx);
-        m[2][2] = 1. - 2. * (xx + yy);
-
-        m[3][3] = 1.;
-
-        let m_sq = SquareMatrix::new(m);
-        // For a pure rotation, the inverse is the transpose.
-        let m_inv_sq = m_sq.transpose();
-
-        Transform::new(&m_sq, &m_inv_sq)
-    }
-}
-
-pub struct DerivativeTerm {
-    kc: Float,
-    kx: Float,
-    ky: Float,
-    kz: Float,
-}
-
-impl DerivativeTerm {
-    pub fn new(kc: Float, kx: Float, ky: Float, kz: Float) -> Self {
-        Self { kc, kx, ky, kz }
-    }
-}
-
-pub struct AnimatedTransform {
-    pub start_transform: Transform<Float>,
-    pub end_transform: Transform<Float>,
-    pub start_time: Float,
-    pub end_time: Float,
-    actually_animated: bool,
-    t: [Vector3f; 2],
-    r: [Quaternion; 2],
-    s: [SquareMatrix<Float, 4>; 2],
-    has_rotation: bool,
-    c1: [DerivativeTerm; 3],
-    c2: [DerivativeTerm; 3],
-    c3: [DerivativeTerm; 3],
-    c4: [DerivativeTerm; 3],
-    c5: [DerivativeTerm; 3],
-}
-
-impl AnimatedTransform {
-    pub fn interpolate(&self, time: Float) -> Transform<Float> {
-        if !self.actually_animated || time <= self.start_time { return self.start_transform }
-        if time >= self.end_time { return self.end_transform }
-
-        let dt = (time - self.start_time) / (self.end_time - self.start_time);
-        let trans = (1.0 - dt) * self.t[0] + dt * self.t[1];
-
-        let rotate = Quaternion::slerp(dt, self.r[0], self.r[1]);
-
-        let scale = (1.0 - dt) * self.s[0] + dt * self.s[1];
-
-        // Return interpolated matrix as product of interpolated components
-        return Transform::translate(trans) * Transform::from(rotate) * Transform::from_matrix(&scale);
-    }
-
-    pub fn apply_inverse_point(&self, p: Point3f, time: Float) -> Point3f {
-        if !self.actually_animated {
-            return self.start_transform.apply_inverse(p);
-        }
-        return self.interpolate(time).apply_inverse(p)
-    }
-
-    pub fn apply_inverse_vector(&self, v: Vector3f, time: Float) -> Vector3f {
-        if !self.actually_animated {
-            return self.start_transform.apply_inverse_vector(v);
-        }
-        return self.interpolate(time).apply_inverse_vector(v)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    fn assert_matrix_approx_eq<const N: usize>(a: &SquareMatrix<f64, N>, b: &SquareMatrix<f64, N>) {
-        const EPSILON: f64 = 1e-9;
-        for i in 0..N {
-            for j in 0..N {
-                assert!(
-                    (a[i][j] - b[i][j]).abs() < EPSILON,
-                    "Matrices differ at ({},{}): {} vs {}\nLeft:\n{}\nRight:\n{}",
-                    i, j, a[i][j], b[i][j], a, b
-                );
-            }
-        }
-    }
-
-    #[test]
-    fn test_inverse_2x2() {
-        let m = SquareMatrix { m: [[4.0, 7.0], [2.0, 6.0]] };
-        let identity = SquareMatrix::identity();
-        let inv = m.inverse().expect("Matrix should be invertible");
-        let product = m * inv;
-        assert_matrix_approx_eq(&product, &identity);
-    }
-    
-    #[test]
-    fn test_inverse_3x3() {
-        let m = SquareMatrix { m: [[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]] };
-        let identity = SquareMatrix::identity();
-        let inv = m.inverse().expect("Matrix should be invertible");
-        let product = m * inv;
-        let product_inv = inv * m;
-        assert_matrix_approx_eq(&product, &identity);
-        assert_matrix_approx_eq(&product_inv, &identity);
-    }
-
-    #[test]
-    fn test_singular_inverse() {
-        let m = SquareMatrix { m: [[1.0, 2.0], [2.0, 4.0]] }; // Determinant is 0
-        assert!(m.inverse().is_none());
-    }
-
-    #[test]
-    fn test_multiplication_2x2() {
-        let a = SquareMatrix { m: [[1.0, 2.0], [3.0, 4.0]] };
-        let b = SquareMatrix { m: [[2.0, 0.0], [1.0, 2.0]] };
-        let expected = SquareMatrix { m: [[4.0, 4.0], [10.0, 8.0]] };
-        let result = a * b;
-        assert_matrix_approx_eq(&result, &expected);
-    }
-
-    #[test]
-    fn test_determinant_3x3() {
-        let m = SquareMatrix { m: [[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]] };
-        assert_eq!(m.determinant(), 1.0);
-    }
-}

src/lib.rs

@@ -1,4 +1,7 @@
 #![allow(unused_imports, dead_code)]
+#![feature(float_erf)]
 
+mod camera;
 mod core;
 mod lights;
+mod utils;

src/lights/diffuse.rs

@@ -1,6 +1,6 @@
 use crate::core::medium::MediumInterface;
 use crate::core::pbrt::Float;
-use crate::core::spectrum::Spectrum;
+use crate::utils::spectrum::Spectrum;
 
 pub struct DiffuseAreaLight {
     pub l_emit: Spectrum,

src/utils/color.rs

@@ -2,8 +2,8 @@ use std::ops::{Index, IndexMut, Add, AddAssign, Sub, SubAssign, Mul, MulAssign,
 use std::fmt;
 
 use crate::core::pbrt::{Float, lerp};
-use crate::core::geometry::Point2f;
-use crate::core::spectrum::Spectrum;
+use super::geometry::Point2f;
+use super::spectrum::Spectrum;
 
 #[derive(Debug, Clone)]
 pub struct XYZ {

src/utils/colorspace.rs

@@ -1,8 +1,8 @@
 use crate::core::pbrt::Float;
-use crate::core::geometry::Point2f;
-use crate::core::transform::SquareMatrix;
-use crate::core::color::{RGBSigmoidPolynomial, RGBToSpectrumTable, RGB, XYZ};
-use crate::core::spectrum::{DenselySampledSpectrum, Spectrum, SampledSpectrum};
+use super::geometry::Point2f;
+use super::transform::SquareMatrix;
+use super::color::{RGBSigmoidPolynomial, RGBToSpectrumTable, RGB, XYZ};
+use super::spectrum::{DenselySampledSpectrum, Spectrum, SampledSpectrum};
 
 use std::cmp::{Eq, PartialEq};
 use std::sync::Arc;

src/utils/containers.rs

@@ -0,0 +1,75 @@
+use std::collections::hash_map::RandomState;
+use std::hash::{BuildHasher, Hash, Hasher};
+use std::ops::{Index, IndexMut};
+use std::sync::RwLock;
+
+use crate::utils::geometry::{Bounds2i, Bounds3i, Point2i, Point3i, Point3f, Vector3i, Vector3f};
+
+pub struct Array2D<T> {
+    values: Vec<T>,
+    extent: Bounds2i,
+}
+
+impl<T> Array2D<T> {
+    pub fn new(extent: Bounds2i) -> Self
+    where
+        T: Default,
+    {
+        let size = extent.area() as usize;
+        let mut values = Vec::with_capacity(size);
+        values.resize_with(size, T::default);
+        Self { values, extent }
+    }
+
+    pub fn with_default(extent: Bounds2i, default_val: T) -> Self
+    where
+        T: Clone,
+    {
+        let size = extent.area() as usize;
+        let values = vec![default_val; size];
+        Self { values, extent }
+    }
+
+    pub fn x_size(&self) -> i32 {
+        self.extent.p_max.x() - self.extent.p_min.x()
+    }
+
+    pub fn y_size(&self) -> i32 {
+        self.extent.p_max.y() - self.extent.p_min.y()
+    }
+
+    pub fn size(&self) -> usize {
+        self.values.len()
+    }
+    
+    pub fn as_slice(&self) -> &[T] {
+        &self.values
+    }
+
+    pub fn as_mut_slice(&mut self) -> &mut [T] {
+        &mut self.values
+    }
+
+    fn get_index(&self, p: Point2i) -> usize {
+        debug_assert!(p.x() >= self.extent.p_min.x() && p.x() < self.extent.p_max.x());
+        debug_assert!(p.y() >= self.extent.p_min.y() && p.y() < self.extent.p_max.y());
+        let width = self.x_size();
+        let pp = Point2i::new(p.x() - self.extent.p_min.x(), p.y() - self.extent.p_min.y());
+        (pp.y() * width + pp.x()) as usize
+    }
+}
+
+impl<T> Index<Point2i> for Array2D<T> {
+    type Output = T;
+    fn index(&self, p: Point2i) -> &Self::Output {
+        let idx = self.get_index(p);
+        &self.values[idx]
+    }
+}
+
+impl<T> IndexMut<Point2i> for Array2D<T> {
+    fn index_mut(&mut self, p: Point2i) -> &mut Self::Output {
+        let idx = self.get_index(p);
+        &mut self.values[idx]
+    }
+}

src/utils/geometry.rs

@@ -3,10 +3,9 @@ use std::ops::{Sub, SubAssign, Add, AddAssign, Div, DivAssign, Mul, MulAssign, N
 use std::sync::Arc;
 use std::f32::consts::PI;
 
-use crate::core::pbrt::{Float, next_float_up, next_float_down};
-use crate::core::pbrt;
-use crate::core::interval::Interval;
 use crate::core::medium::Medium;
+use crate::core::pbrt::{Float, next_float_up, next_float_down, lerp, clamp_t};
+use super::interval::Interval;
 
 pub trait Tuple<T, const N: usize>:
     Sized + Copy + Index<usize, Output = T> + IndexMut<usize>
@@ -45,6 +44,13 @@ macro_rules! impl_tuple_core {
             #[inline] pub fn zero() -> Self { Self([T::zero(); N]) }
         }
 
+        impl<T, const N: usize> $Struct<T, N> where T: Copy {
+            #[inline]
+            pub fn fill(value: T) -> Self {
+                Self([value; N])
+            }
+        }
+
         impl<T, const N: usize> Index<usize> for $Struct<T, N> {
             type Output = T;
             #[inline] fn index(&self, index: usize) -> &Self::Output { &self.0[index] }
@@ -191,6 +197,8 @@ impl_dot_for!(Normal, Vector);
 impl_dot_for!(Vector, Normal);
 impl_dot_for!(Vector, Vector);
 impl_dot_for!(Normal, Normal);
+impl_dot_for!(Normal, Point);
+impl_dot_for!(Vector, Point);
 
 impl<T: Copy, const N: usize> From<Vector<T, N>> for Normal<T, N> {
     fn from(v: Vector<T, N>) -> Self { Self(v.0) }
@@ -199,6 +207,13 @@ impl<T: Copy, const N: usize> From<Normal<T, N>> for Vector<T, N> {
     fn from(n: Normal<T, N>) -> Self { Self(n.0) }
 }
 
+impl<T: Copy, const N: usize> From<Vector<T, N>> for Point<T, N> {
+    fn from(v: Vector<T, N>) -> Self { Self(v.0) }
+}
+
+impl<T: Copy, const N: usize> From<Point<T, N>> for Vector<T, N> {
+    fn from(n: Point<T, N>) -> Self { Self(n.0) }
+}
 macro_rules! impl_accessors {
     ($Struct:ident) => {
         impl<T: Copy> $Struct<T, 2> {
@@ -421,8 +436,16 @@ where T: Num + PartialOrd + Copy + Neg<Output = T>
     }
 }
 
+pub fn abs_cos_theta(w: Vector3f) -> Float {
+    w.z().abs()
+}
+
 // SPHERICAL GEOMETRY
-pub fn spherical_direction<T: NumFloat>(sin_theta: Float, cos_theta: Float, phi: Float) -> Vector3f {
+pub fn same_hemisphere(w: Vector3f, wp: Vector3f) -> bool {
+    w.z() * wp.z() > 0.
+}
+
+pub fn spherical_direction(sin_theta: Float, cos_theta: Float, phi: Float) -> Vector3f {
     Vector3f::new(sin_theta * phi.cos(), sin_theta * phi.sin(), cos_theta)
 }
 
@@ -431,7 +454,7 @@ pub fn spherical_triangle_area<T: NumFloat>(a: Vector3f, b: Vector3f, c: Vector3
 }
 
 pub fn spherical_theta<T: NumFloat>(v: Vector3f) -> Float {
-    pbrt::clamp_t(v.z(), -1.0, 1.0).acos()
+    clamp_t(v.z(), -1.0, 1.0).acos()
 }
 
 pub fn spherical_phi<T: NumFloat>(v: Vector3f) -> Float {
@@ -504,11 +527,19 @@ where
         d.0.iter().fold(T::one(), |acc, &val| acc * val)
     }
 
+    pub fn expand(&self, delta: T) -> Self {
+        let mut p_min = self.p_min;
+        let mut p_max = self.p_max;
+        p_min = p_min - Vector::fill(delta);
+        p_max = p_max + Vector::fill(delta);
+        Self { p_min, p_max }
+    }
+
     pub fn lerp(&self, t: Point<T, N>) -> Point<T, N>
     {
         let mut results_arr = [T::zero(); N];
         for i in 0..N {
-            results_arr[i] = pbrt::lerp(t[i], self.p_min[i], self.p_max[i])
+            results_arr[i] = lerp(t[i], self.p_min[i], self.p_max[i])
         }
 
         Point(results_arr)
@@ -576,8 +607,13 @@ where
 }
 
 pub type Bounds2<T> = Bounds<T, 2>;
+pub type Bounds2f = Bounds2<Float>;
+pub type Bounds2i = Bounds2<i32>;
+pub type Bounds2fi = Bounds2<Interval>;
 pub type Bounds3<T> = Bounds<T, 3>;
+pub type Bounds3i = Bounds3<i32>;
 pub type Bounds3f = Bounds3<Float>;
+pub type Bounds3fi = Bounds3<Interval>;
 
 impl<T> Bounds3<T> 
 where 
@@ -600,6 +636,16 @@ where
     }
 }
 
+impl<T> Bounds2<T>
+where 
+    T: Num + Copy + Default
+{
+    pub fn area(&self) -> T {
+        let d: Vector2<T> = self.p_max - self.p_min;
+        d.x() * d.y()
+    }
+}
+
 #[derive(Copy, Clone, Default, PartialEq)]
 pub struct Frame {
     pub x: Vector3f,
@@ -632,7 +678,7 @@ impl Frame {
     }
 }
 
-#[derive(Clone, Default)]
+#[derive(Clone, Debug)]
 pub struct Ray {
     pub o: Point3f,
     pub d: Vector3f,
@@ -642,8 +688,24 @@ pub struct Ray {
     pub differential: Option<RayDifferential>,
 }
 
+impl Default for Ray {
+    fn default() -> Self {
+        Self {
+            o: Point3f::new(0.0, 0.0, 0.0),
+            d: Vector3f::new(0.0, 0.0, 0.0),
+            medium: None,
+            time: 0.0,
+            differential: None,
+        }
+    }
+}
+
 impl Ray {
-    pub fn operator(&self, t: Float) -> Point3f {
+    pub fn new(o: Point3f, d: Vector3f, time: Option<Float>, medium: Option<Arc<Medium>>) -> Self {
+        Self { o, d, time: time.unwrap_or_else(|| Self::default().time), medium, ..Self::default() }
+    }
+
+    pub fn evaluate(&self, t: Float) -> Point3f {
         self.o + self.d * t
     }
 

src/utils/math.rs

@@ -0,0 +1,229 @@
+use super::geometry::{Point2f, Vector3f};
+use crate::core::pbrt::{PI, PI_OVER_4, square, lerp, Float, ONE_MINUS_EPSILON, clamp_t};
+use std::mem;
+use std::ops::{Mul, Neg, Add};
+
+#[inline]
+fn fma<T>(a: T, b: T, c: T) -> T 
+where
+    T: Mul<Output = T> + Add<Output = T> + Copy
+{
+    a * b + c
+}
+
+#[inline]
+fn difference_of_products<T>(a: T, b: T, c: T, d: T) -> T 
+where 
+    T: Mul<Output = T> + Add<Output = T> + Neg<Output = T> + Copy
+{
+    let cd = c * d;
+    let difference_of_products = fma(a, b, -cd);
+    let error = fma(-c, d, cd);
+    return difference_of_products + error;
+}
+
+#[inline]
+pub fn safe_asin(x: Float) -> Float {
+    if x >= -1.0001 && x <= 1.0001 {
+        clamp_t(x, -1., 1.).asin()
+    } else {
+        panic!("Not valid value for asin")
+    }
+}
+
+#[inline]
+pub fn safe_acos(x: Float) -> Float {
+    if x >= -1.0001 && x <= 1.0001 {
+        clamp_t(x, -1., 1.).asin()
+    } else {
+        panic!("Not valid value for asin")
+    }
+}
+
+#[inline]
+pub fn sinx_over_x(x: Float) -> Float {
+    if 1. - x * x == 1. {
+        return 1.;
+    }
+    return x.sin() / x;
+}
+
+pub fn windowed_sinc(x: Float, radius: Float, tau: Float) -> Float {
+    if x.abs() > radius {
+        return 0.;
+    }
+    return sinc(x) * sinc(x / tau);
+}
+
+pub fn sinc(x: Float) -> Float {
+    return sinx_over_x(PI * x);
+}
+
+pub fn quadratic(a: Float, b: Float, c: Float) -> Option<(Float, Float)> {
+    if a == 0. {
+        if b == 0. {
+            return None;
+        }
+        let t0 = -c / b;
+        let t1 = - c / b;
+        return Some((t0, t1));
+    }
+
+    let discrim = difference_of_products(b, b, 4. * a, c);
+    if discrim < 0. {
+        return None;
+    }
+    let root_discrim = discrim.sqrt();
+
+    let q = -0.5 * (b + root_discrim.copysign(b));
+    let mut t0 = q / a;
+    let mut t1 = c / q;
+    if t0 > t1 {
+        mem::swap(&mut t0, &mut t1);
+    }
+
+    Some((t0, t1))
+}
+
+pub fn wrap_equal_area_square(uv: &mut Point2f) -> Point2f {
+    if uv[0] < 0. {
+        uv[0] = -uv[0];     // mirror across u = 0
+        uv[1] = 1. - uv[1];  // mirror across v = 0.5
+    } else if uv[0] > 1. 
+    {
+        uv[0] = 2. - uv[0];  // mirror across u = 1
+        uv[1] = 1. - uv[1];  // mirror across v = 0.5
+    }
+    if uv[1] < 0. {
+        uv[0] = 1. - uv[0];  // mirror across u = 0.5
+        uv[1] = -uv[1];     // mirror across v = 0;
+    } else if uv[1] > 1. 
+    {
+        uv[0] = 1. - uv[0];  // mirror across u = 0.5
+        uv[1] = 2. - uv[1];  // mirror across v = 1
+    }
+    *uv
+}
+
+pub fn equal_area_square_to_sphere(p: Point2f) -> Vector3f {
+    assert!(p.x() >= 0. && p.x() <= 1. && p.y() >= 0. && p.y() <= 1.);
+
+    // Transform _p_ to $[-1,1]^2$ and compute absolute values
+    let u = 2. * p.x() - 1.;
+    let v = 2. * p.y() - 1.;
+    let up = u.abs(); 
+    let vp = v.abs();
+
+    // Compute radius _r_ as signed distance from diagonal
+    let signed_distance = 1. - (up + vp);
+    let d = signed_distance.abs();
+    let r = 1. - d;
+
+    // Compute angle $\phi$ for square to sphere mapping
+    let mut phi = if r == 0. { 1. } else  {(vp - up) / r + 1.};
+    phi /= PI_OVER_4;
+
+    // Find $z$ coordinate for spherical direction
+    let z = (1. - square(r)).copysign(signed_distance);
+
+    // Compute $\cos\phi$ and $\sin\phi$ for original quadrant and return vector
+    let cos_phi = phi.cos().copysign(u);
+    let sin_phi = phi.sin().copysign(v);
+    return Vector3f::new(cos_phi * r * (2. - square(r)).sqrt(), 
+        sin_phi * r * (2. - square(r)).sqrt(), z);
+}
+
+pub fn gaussian(x: Float, y: Float, sigma: Float) -> Float
+{
+    (-(x * x + y * y) / (2. * sigma * sigma)).exp()
+}
+
+pub fn gaussian_integral(x0: Float, x1: Float, mu: Float, sigma: Float) -> Float {
+    assert!(sigma > 0.);
+    let sigma_root2 = sigma * 1.414213562373095;
+    0.5 * (((mu - x0) / sigma_root2).erf() - ((mu - x1) / sigma_root2).erf())
+}
+
+pub fn sample_linear(u: Float, a: Float, b: Float) -> Float {
+    assert!(a >= 0. && b >= 0.);
+    if u == 0. && a == 0. {
+        return 0.
+    }
+    let x = u * (a + b) / (a + (lerp(u, square(a), square(b))));
+    x.min(ONE_MINUS_EPSILON)
+}
+
+fn next_float_down(a: f32) -> f32 {
+    if a.is_infinite() && a < 0.0 {
+        return a;
+    }
+    if a == 0.0 {
+        return -0.0;
+    }
+    let bits = a.to_bits();
+    if a > 0.0 {
+        f32::from_bits(bits - 1)
+    } else {
+        f32::from_bits(bits + 1)
+    }
+}
+
+pub fn sample_discrete(weights: &[f32], u: f32, pmf: Option<&mut f32>, u_remapped: Option<&mut f32>,
+) -> Option<usize> {
+    // Handle empty weights for discrete sampling.
+    if weights.is_empty() {
+        if let Some(p) = pmf {
+            *p = 0.0;
+        }
+        return None;
+    }
+
+    // Compute the sum of all weights.
+    let sum_weights: f32 = weights.iter().sum();
+
+    // If the total weight is zero, sampling is not possible.
+    if sum_weights == 0.0 {
+        if let Some(p) = pmf {
+            *p = 0.0;
+        }
+        return None;
+    }
+
+    // Compute rescaled u' sample, ensuring it's strictly less than sum_weights.
+    let mut up = u * sum_weights;
+    if up >= sum_weights {
+        up = next_float_down(up);
+    }
+
+    // Find the offset in weights corresponding to the rescaled sample u'.
+    let mut offset = 0;
+    let mut sum = 0.0;
+    while sum + weights[offset] <= up {
+        sum += weights[offset];
+        offset += 1;
+        // This assertion should hold true due to the guard `up = next_float_down(up)`.
+        debug_assert!(offset < weights.len());
+    }
+
+    // Compute PMF and remapped u value, if requested.
+    if let Some(p) = pmf {
+        *p = weights[offset] / sum_weights;
+    }
+    if let Some(ur) = u_remapped {
+        let weight = weights[offset];
+        // The logic guarantees weight > 0 here if sum_weights > 0.
+        *ur = ((up - sum) / weight).min(ONE_MINUS_EPSILON);
+    }
+
+    Some(offset)
+}
+
+pub fn sample_tent(u: Float, r: Float) -> Float {
+    let mut u_remapped = 0.0;
+    let offset = sample_discrete(&[0.5, 0.5], u, None, Some(&mut u_remapped)).expect("Discrete sampling shouldn't fail");
+    if offset == 0 {
+        return -r + r * sample_linear(u, 0., 1.);
+    } else {
+        return r * sample_linear(u, 1., 0.);
+    }
+}

src/utils/mod.rs

@@ -0,0 +1,11 @@
+pub mod color;
+pub mod colorspace;
+pub mod containers;
+pub mod geometry;
+pub mod interval;
+pub mod math;
+pub mod quaternion;
+pub mod spectrum;
+pub mod transform;
+pub mod scattering;
+pub mod sampling;

src/utils/quaternion.rs

@@ -1,10 +1,11 @@
-use crate::core::geometry::{Vector3f, Dot, Normed};
-use crate::core::pbrt::{safe_asin, sinx_over_x, Float};
 use std::ops::{Index, IndexMut};
 use std::ops::{Add, AddAssign, Sub, SubAssign, Mul, MulAssign, Div, DivAssign, Neg};
 use std::f32::consts::PI;
 
-#[derive(Copy, Clone, PartialEq)]
+use crate::core::pbrt::{safe_asin, sinx_over_x, Float};
+use super::geometry::{Vector3f, Dot, Normed};
+
+#[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Quaternion {
     pub v: Vector3f,
     pub w: Float,
@@ -106,4 +107,12 @@ impl Quaternion {
         return q1 * (1. - t) * sinx_over_x((1. - t) * theta) / sin_theta_over_theta +
                q2 * t * sinx_over_x(t * theta) / sin_theta_over_theta;
     }
+
+    pub fn length(&self) -> Float {
+        self.v.norm()
+    }
+
+    pub fn normalize(&self) -> Self {
+        Quaternion { v: self.v.normalize(), w: self.w}
+    }
 }

src/utils/sampling.rs

@@ -0,0 +1,13 @@
+use super::geometry::{Point2f, Vector3f};
+use crate::core::pbrt::{square, PI, INV_2_PI};
+
+pub fn sample_uniform_hemisphere(u: Point2f) -> Vector3f {
+    let z = u[0];
+    let r = SafeSqrt(1 - square(z));
+    let phi = 2 * PI * u[1];
+    Vector3f::new(r * phi.cos(), r * phi.sin(), z);
+}
+
+pub fn uniform_hemisphere_pdf() -> Float {
+    INV_2_PI
+}

src/utils/scattering.rs

@@ -0,0 +1,31 @@
+use crate::utils::geometry::{Normal3f, Vector3f, Dot};
+use crate::core::pbrt::{Float, square};
+
+pub fn refract(wi: Vector3f, n: Normal3f, eta_ratio: Float) -> Option<(Vector3f, Float)> {
+    // Make local mutable copies of variables that may be changed
+    let mut n_interface = n;
+    let mut eta = eta_ratio;
+
+    let mut cos_theta_i = n_interface.dot(wi);
+
+    // Potentially flip interface orientation for Snell's law
+    if cos_theta_i < 0.0 {
+        eta = 1.0 / eta;
+        cos_theta_i = -cos_theta_i;
+        n_interface = -n_interface;
+    }
+
+    // Compute sin^2(theta_t) using Snell's law
+    let sin2_theta_i = (1.0 - square(cos_theta_i)).max(0.0);
+    let sin2_theta_t = sin2_theta_i / square(eta);
+
+    // Handle total internal reflection
+    if sin2_theta_t >= 1.0 {
+        return None;
+    }
+
+    let cos_theta_t = (1.0 - sin2_theta_t).sqrt();
+
+    let wt = -wi / eta + (cos_theta_i / eta - cos_theta_t) * Vector3f::from(n_interface);
+    Some((wt, eta))
+}

src/utils/spectrum.rs

@@ -3,9 +3,10 @@ use std::ops::{Add, AddAssign, Sub, SubAssign, Mul, MulAssign, Div, DivAssign, N
 use once_cell::sync::Lazy;
 use std::sync::Arc;
 
-use crate::core::colorspace::RGBColorspace;
 use crate::core::pbrt::Float;
-use crate::core::color::{RGB, XYZ, RGBSigmoidPolynomial};
+use crate::core::cie;
+use super::color::{RGB, XYZ, RGBSigmoidPolynomial};
+use super::colorspace::RGBColorspace;
 
 pub const CIE_Y_INTEGRAL: Float = 106.856895;
 pub const N_SPECTRUM_SAMPLES: usize = 1200;
@@ -167,6 +168,13 @@ impl Mul<Float> for SampledSpectrum {
     }
 }
 
+impl Mul<SampledSpectrum> for Float {
+    type Output = SampledSpectrum;
+    fn mul(self, rhs: SampledSpectrum) -> SampledSpectrum {
+        rhs * self
+    }
+}
+
 impl MulAssign<Float> for SampledSpectrum {
     fn mul_assign(&mut self, rhs: Float) {
         for i in 0..N_SPECTRUM_SAMPLES {
@@ -184,6 +192,14 @@ impl DivAssign for SampledSpectrum {
     }
 }
 
+impl Div for SampledSpectrum {
+    type Output = Self;
+    fn div(self, rhs: Self) -> Self::Output {
+        let mut ret = self;
+        ret /= rhs;
+        ret
+    }
+}
 impl Div<Float> for SampledSpectrum {
     type Output = Self;
 
@@ -404,13 +420,8 @@ impl UnboundedRGBSpectrum {
     pub fn new(cs: RGBColorspace, rgb: RGB) -> Self {
         let m = rgb.r.max(rgb.g).max(rgb.b);
         let scale = 2.0 * m;
-        let lambda;
-        if scale == 1.0 {
-            lambda = rgb / scale; 
-        } else {
-            lambda = RGB::new(0.0, 0.0, 0.0);
-        }
-        Self { scale, rsp: cs.to_rgb_coeffs(lambda) }
+        let scaled_rgb = if scale != 0.0 { rgb / scale } else { RGB::new(0.0, 0.0, 0.0) };
+        Self { scale, rsp: cs.to_rgb_coeffs(scaled_rgb) }
     }
     pub fn sample_at(&self, lambda: Float) -> Float {
         self.scale * self.rsp.evaluate(lambda)
@@ -577,15 +588,31 @@ pub struct RGBUnboundedSpectrum(pub RGBSpectrum);
 pub mod spectra {
     use super::*;
     pub static X: Lazy<Spectrum> = Lazy::new(|| {
-        let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_X_DATA);
+        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_X);
+        let dss = DenselySampledSpectrum::from_spectrum(
+            &Spectrum::PiecewiseLinear(pls),
+            LAMBDA_MIN,
+            LAMBDA_MAX,
+        );
         Spectrum::DenselySampled(dss)
     });
     pub static Y: Lazy<Spectrum> = Lazy::new(|| {
-        let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Y_DATA);
+        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_Y);
+        let dss = DenselySampledSpectrum::from_spectrum(
+            &Spectrum::PiecewiseLinear(pls),
+            LAMBDA_MIN,
+            LAMBDA_MAX,
+        );
         Spectrum::DenselySampled(dss)
     });
     pub static Z: Lazy<Spectrum> = Lazy::new(|| {
-        let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Z_DATA);
+        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_Z);
+        let dss = DenselySampledSpectrum::from_spectrum(
+            &Spectrum::PiecewiseLinear(pls),
+            LAMBDA_MIN,
+            LAMBDA_MAX,
+        );
         Spectrum::DenselySampled(dss)
     });
+
  }