Major refactoring

src/core/geometry/bounds.rs

@@ -0,0 +1,344 @@
+use super::{Float, NumFloat};
+use super::{Point, Point2f, Point3, Point3f, Vector, Vector2, Vector2f, Vector3, Vector3f};
+use crate::core::geometry::traits::{Sqrt, VectorLike};
+use crate::core::geometry::{max, min};
+use crate::utils::interval::Interval;
+use crate::utils::math::lerp;
+use num_traits::{Bounded, Num};
+use std::mem;
+use std::ops::{Add, Div, DivAssign, Mul, Sub};
+
+// AABB BOUNDING BOXES
+
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct Bounds<T, const N: usize> {
+    pub p_min: Point<T, N>,
+    pub p_max: Point<T, N>,
+}
+
+impl<'a, T, const N: usize> IntoIterator for &'a Bounds<T, N> {
+    type Item = &'a Point<T, N>;
+    type IntoIter = std::array::IntoIter<&'a Point<T, N>, 2>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        [&self.p_min, &self.p_max].into_iter()
+    }
+}
+
+impl<T, const N: usize> Bounds<T, N>
+where
+    T: Num + PartialOrd + Copy,
+{
+    pub fn from_point(p: Point<T, N>) -> Self {
+        Self { p_min: p, p_max: p }
+    }
+
+    pub fn from_points(p1: Point<T, N>, p2: Point<T, N>) -> Self {
+        let mut p_min_arr = [T::zero(); N];
+        let mut p_max_arr = [T::zero(); N];
+
+        for i in 0..N {
+            if p1[i] < p2[i] {
+                p_min_arr[i] = p1[i];
+                p_max_arr[i] = p2[i];
+            } else {
+                p_min_arr[i] = p2[i];
+                p_max_arr[i] = p1[i];
+            }
+        }
+
+        Self {
+            p_min: Point(p_min_arr),
+            p_max: Point(p_max_arr),
+        }
+    }
+
+    pub fn union_point(self, p: Point<T, N>) -> Self {
+        let mut p_min = self.p_min;
+        let mut p_max = self.p_max;
+        for i in 0..N {
+            p_min[i] = min(p_min[i], p[i]);
+            p_max[i] = max(p_max[i], p[i]);
+        }
+        Self { p_min, p_max }
+    }
+
+    pub fn union(self, b2: Self) -> Self {
+        let mut p_min = self.p_min;
+        let mut p_max = self.p_max;
+        for i in 0..N {
+            p_min[i] = min(p_min[i], b2.p_min[i]);
+            p_max[i] = max(p_max[i], b2.p_max[i]);
+        }
+        Self { p_min, p_max }
+    }
+
+    pub fn diagonal(&self) -> Vector<T, N> {
+        self.p_max - self.p_min
+    }
+
+    pub fn centroid(&self) -> Point<T, N> {
+        let two = T::one() + T::one();
+        self.p_min + (self.diagonal() / two)
+    }
+
+    pub fn volume(&self) -> T {
+        let d = self.diagonal();
+        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] = lerp(t[i], self.p_min[i], self.p_max[i])
+        }
+
+        Point(results_arr)
+    }
+
+    pub fn max_dimension(&self) -> usize
+    where
+        Point<T, N>: Sub<Output = Vector<T, N>>,
+    {
+        let d = self.diagonal();
+        let mut max_dim = 0;
+        let mut max_span = d[0];
+
+        for i in 1..N {
+            if d[i] > max_span {
+                max_span = d[i];
+                max_dim = i;
+            }
+        }
+        max_dim
+    }
+
+    pub fn offset(&self, p: &Point<T, N>) -> Vector<T, N>
+    where
+        Point<T, N>: Sub<Output = Vector<T, N>>,
+        Vector<T, N>: DivAssign<T>,
+    {
+        let mut o = *p - self.p_min;
+        let d = self.diagonal();
+        for i in 0..N {
+            if d[i] > T::zero() {
+                o[i] = o[i] / d[i];
+            }
+        }
+        o
+    }
+
+    pub fn corner(&self, corner_index: usize) -> Point<T, N> {
+        Point(std::array::from_fn(|i| {
+            if (corner_index >> i) & 1 == 1 {
+                self.p_max[i]
+            } else {
+                self.p_min[i]
+            }
+        }))
+    }
+
+    pub fn overlaps(&self, rhs: &Self) -> bool {
+        for i in 0..N {
+            if self.p_max[i] < rhs.p_min[i] || self.p_min[i] > rhs.p_max[i] {
+                return false;
+            }
+        }
+        true
+    }
+
+    pub fn contains(&self, p: Point<T, N>) -> bool {
+        (0..N).all(|i| p[i] >= self.p_min[i] && p[i] <= self.p_max[i])
+    }
+
+    pub fn contains_exclusive(&self, p: Point<T, N>) -> bool {
+        (0..N).all(|i| p[i] >= self.p_min[i] && p[i] < self.p_max[i])
+    }
+
+    pub fn is_empty(&self) -> bool {
+        (0..N).any(|i| self.p_min[i] >= self.p_max[i])
+    }
+
+    pub fn is_degenerate(&self) -> bool {
+        (0..N).any(|i| self.p_min[i] > self.p_max[i])
+    }
+}
+
+impl<T, const N: usize> Default for Bounds<T, N>
+where
+    T: Bounded + Copy,
+{
+    fn default() -> Self {
+        Self {
+            p_min: Point([T::max_value(); N]),
+            p_max: Point([T::min_value(); N]),
+        }
+    }
+}
+
+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
+    T: Num + PartialOrd + Copy + Default,
+{
+    pub fn surface_area(&self) -> T {
+        let d = self.diagonal();
+        let two = T::one() + T::one();
+        two * (d.x() * d.y() + d.x() * d.z() + d.y() * d.z())
+    }
+}
+
+impl<T> Bounds3<T>
+where
+    T: NumFloat + PartialOrd + Copy + Default + Sqrt,
+{
+    pub fn bounding_sphere(&self) -> (Point3<T>, T) {
+        let two = T::one() + T::one();
+        let center = self.p_min + self.diagonal() / two;
+        let radius = if self.contains(center) {
+            center.distance(self.p_max)
+        } else {
+            T::zero()
+        };
+        (center, radius)
+    }
+
+    pub fn insersect(&self, o: Point3<T>, d: Vector3<T>, t_max: T) -> Option<(T, T)> {
+        let mut t0 = T::zero();
+        let mut t1 = t_max;
+
+        for i in 0..3 {
+            let inv_ray_dir = T::one() / d[i];
+            let mut t_near = (self.p_min[i] - o[i]) * inv_ray_dir;
+            let mut t_far = (self.p_max[i] - o[i]) * inv_ray_dir;
+            if t_near > t_far {
+                mem::swap(&mut t_near, &mut t_far);
+            }
+            t0 = if t_near > t0 { t_near } else { t0 };
+            t1 = if t_far < t1 { t_far } else { t1 };
+            if t0 > t1 {
+                return None;
+            }
+        }
+        Some((t0, t1))
+    }
+}
+
+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()
+    }
+}
+
+impl Bounds3f {
+    #[inline(always)]
+    pub fn intersect_p(
+        &self,
+        o: Point3f,
+        ray_t_max: Float,
+        inv_dir: Vector3f,
+        dir_is_neg: &[usize; 3],
+    ) -> Option<(Float, Float)> {
+        let bounds = [&self.p_min, &self.p_max];
+
+        // Check X
+        let mut t_min = (bounds[dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+        let mut t_max = (bounds[1 - dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+
+        // Check Y
+        let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+        let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+
+        if t_min > ty_max || ty_min > t_max {
+            return None;
+        }
+        if ty_min > t_min {
+            t_min = ty_min;
+        }
+        if ty_max < t_max {
+            t_max = ty_max;
+        }
+
+        // Check Z
+        let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+
+        if t_min > tz_max || tz_min > t_max {
+            return None;
+        }
+        if tz_min > t_min {
+            t_min = tz_min;
+        }
+        if tz_max < t_max {
+            t_max = tz_max;
+        }
+
+        if (t_min < ray_t_max) && (t_max > 0.0) {
+            Some((t_min, t_max))
+        } else {
+            None
+        }
+    }
+
+    pub fn intersect_with_inverse(&self, o: Point3f, d: Vector3f, ray_t_max: Float) -> bool {
+        let inv_dir = Vector3::new(1.0 / d.x(), 1.0 / d.y(), 1.0 / d.z());
+        let dir_is_neg: [usize; 3] = [
+            (d.x() < 0.0) as usize,
+            (d.y() < 0.0) as usize,
+            (d.z() < 0.0) as usize,
+        ];
+
+        let bounds = [&self.p_min, &self.p_max];
+
+        // Check for ray intersection against x and y slabs
+        let mut t_min = (bounds[dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+        let mut t_max = (bounds[1 - dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+        let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+        let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+
+        if t_min > ty_max || ty_min > t_max {
+            return false;
+        }
+        if ty_min > t_min {
+            t_min = ty_min;
+        }
+        if ty_max < t_max {
+            t_max = ty_max;
+        }
+
+        let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+
+        if t_min > tz_max || tz_min > t_max {
+            return false;
+        }
+        if tz_min > t_min {
+            t_min = tz_min;
+        }
+        if tz_max < t_max {
+            t_max = tz_max;
+        }
+
+        (t_min < ray_t_max) && (t_max > 0.0)
+    }
+}

src/core/geometry/cone.rs

@@ -0,0 +1,116 @@
+use super::{Bounds3f, Float, PI, Point3f, Vector3f, VectorLike};
+use crate::utils::math::{degrees, safe_acos, safe_asin, safe_sqrt, square};
+use crate::utils::transform::Transform;
+
+#[derive(Debug, Clone)]
+pub struct DirectionCone {
+    pub w: Vector3f,
+    pub cos_theta: Float,
+}
+
+impl Default for DirectionCone {
+    fn default() -> Self {
+        Self {
+            w: Vector3f::zero(),
+            cos_theta: Float::INFINITY,
+        }
+    }
+}
+
+impl DirectionCone {
+    pub fn new(w: Vector3f, cos_theta: Float) -> Self {
+        Self {
+            w: w.normalize(),
+            cos_theta,
+        }
+    }
+
+    pub fn new_from_vector(w: Vector3f) -> Self {
+        Self::new(w, 1.0)
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.cos_theta == Float::INFINITY
+    }
+
+    pub fn entire_sphere() -> Self {
+        Self::new(Vector3f::new(0., 0., 1.), -1.)
+    }
+
+    pub fn closest_vector_income(&self, wt: Vector3f) -> Vector3f {
+        let wp = wt.normalize();
+        let w = self.w;
+        if wp.dot(w) > self.cos_theta {
+            return wp;
+        }
+
+        let sin_theta = -safe_sqrt(1. - self.cos_theta * self.cos_theta);
+        let a = wp.cross(w);
+        self.cos_theta * w
+            + sin_theta / a.norm()
+                * Vector3f::new(
+                    w.x()
+                        * (wp.y() * w.y() + wp.z() * w.z()
+                            - wp.x() * (square(w.y() + square(w.z())))),
+                    w.y()
+                        * (wp.x() * w.x() + wp.z() * w.z()
+                            - wp.y() * (square(w.x() + square(w.z())))),
+                    w.z()
+                        * (wp.x() * w.x() + wp.y() * w.y()
+                            - wp.z() * (square(w.x() + square(w.y())))),
+                )
+    }
+
+    pub fn inside(d: &DirectionCone, w: Vector3f) -> bool {
+        !d.is_empty() && d.w.dot(w.normalize()) > d.cos_theta
+    }
+
+    pub fn bound_subtended_directions(b: &Bounds3f, p: Point3f) -> DirectionCone {
+        let (p_center, radius) = b.bounding_sphere();
+        if p.distance_squared(p_center) < square(radius) {
+            return DirectionCone::entire_sphere();
+        }
+
+        let w = (p_center - p).normalize();
+        let sin2_theta_max = square(radius) / p_center.distance_squared(p);
+        let cos_theta_max = safe_sqrt(1. - sin2_theta_max);
+        DirectionCone::new(w, cos_theta_max)
+    }
+
+    pub fn union(a: &DirectionCone, b: &DirectionCone) -> DirectionCone {
+        if a.is_empty() {
+            return b.clone();
+        }
+        if b.is_empty() {
+            return a.clone();
+        }
+
+        // Handle the cases where one cone is inside the other
+        let theta_a = safe_acos(a.cos_theta);
+        let theta_b = safe_acos(b.cos_theta);
+        let theta_d = a.w.angle_between(b.w);
+
+        if (theta_d + theta_b).min(PI) <= theta_b {
+            return a.clone();
+        }
+        if (theta_d + theta_a).min(PI) <= theta_a {
+            return b.clone();
+        }
+
+        // Compute the spread angle of the merged cone, $\theta_o$
+        let theta_o = (theta_a + theta_d + theta_b) / 2.;
+        if theta_o >= PI {
+            return DirectionCone::entire_sphere();
+        }
+
+        // Find the merged cone's axis and return cone union
+        let theta_r = theta_o - theta_a;
+        let wr = a.w.cross(b.w);
+        if wr.norm_squared() >= 0. {
+            return DirectionCone::entire_sphere();
+        }
+
+        let w = Transform::rotate_around_axis(degrees(theta_r), wr).apply_to_vector(a.w);
+        DirectionCone::new(w, theta_o.cos())
+    }
+}

src/core/geometry/mod.rs

@@ -0,0 +1,122 @@
+pub mod bounds;
+pub mod cone;
+pub mod primitives;
+pub mod ray;
+pub mod traits;
+
+pub use self::bounds::{Bounds, Bounds2f, Bounds2fi, Bounds2i, Bounds3f, Bounds3fi, Bounds3i};
+pub use self::cone::DirectionCone;
+pub use self::primitives::{
+    Frame, Normal, Normal3f, Point, Point2f, Point2fi, Point2i, Point3, Point3f, Point3fi, Point3i,
+    Vector, Vector2, Vector2f, Vector2fi, Vector2i, Vector3, Vector3f, Vector3fi, Vector3i,
+};
+pub use self::ray::{Ray, RayDifferential};
+pub use self::traits::{Lerp, Sqrt, Tuple, VectorLike};
+
+use crate::core::pbrt::{Float, PI, clamp_t};
+use crate::utils::math::square;
+use num_traits::Float as NumFloat;
+
+#[inline]
+pub fn min<T: PartialOrd>(a: T, b: T) -> T {
+    if a < b { a } else { b }
+}
+
+#[inline]
+pub fn max<T: PartialOrd>(a: T, b: T) -> T {
+    if a > b { a } else { b }
+}
+
+#[inline]
+pub fn cos_theta(w: Vector3f) -> Float {
+    w.z()
+}
+#[inline]
+pub fn abs_cos_theta(w: Vector3f) -> Float {
+    w.z().abs()
+}
+#[inline]
+pub fn cos2_theta(w: Vector3f) -> Float {
+    square(w.z())
+}
+#[inline]
+pub fn sin2_theta(w: Vector3f) -> Float {
+    0_f32.max(1. - cos2_theta(w))
+}
+#[inline]
+pub fn sin_theta(w: Vector3f) -> Float {
+    sin2_theta(w).sqrt()
+}
+#[inline]
+pub fn tan_theta(w: Vector3f) -> Float {
+    sin_theta(w) / cos_theta(w)
+}
+#[inline]
+pub fn tan2_theta(w: Vector3f) -> Float {
+    sin2_theta(w) / cos2_theta(w)
+}
+#[inline]
+pub fn cos_phi(w: Vector3f) -> Float {
+    let sin_theta = sin_theta(w);
+    if sin_theta == 0. {
+        1.
+    } else {
+        clamp_t(w.x() / sin_theta, -1., 1.)
+    }
+}
+
+#[inline]
+pub fn sin_phi(w: Vector3f) -> Float {
+    let sin_theta = sin_theta(w);
+    if sin_theta == 0. {
+        0.
+    } else {
+        clamp_t(w.y() / sin_theta, -1., 1.)
+    }
+}
+
+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)
+}
+
+pub fn spherical_triangle_area(a: Vector3f, b: Vector3f, c: Vector3f) -> Float {
+    (2.0 * (a.dot(b.cross(c))).atan2(1.0 + a.dot(b) + a.dot(c) + b.dot(c))).abs()
+}
+
+pub fn spherical_quad_area(a: Vector3f, b: Vector3f, c: Vector3f, d: Vector3f) -> Float {
+    let mut axb = a.cross(b);
+    let mut bxc = b.cross(c);
+    let mut cxd = c.cross(d);
+    let mut dxa = d.cross(a);
+    if axb.norm_squared() == 0.
+        || bxc.norm_squared() == 0.
+        || cxd.norm_squared() == 0.
+        || dxa.norm_squared() == 0.
+    {
+        return 0.;
+    }
+    axb = axb.normalize();
+    bxc = bxc.normalize();
+    cxd = cxd.normalize();
+    dxa = dxa.normalize();
+
+    let alpha = dxa.angle_between(-axb);
+    let beta = axb.angle_between(-bxc);
+    let gamma = bxc.angle_between(-cxd);
+    let delta = cxd.angle_between(-dxa);
+
+    (alpha + beta + gamma + delta - 2. * PI).abs()
+}
+
+pub fn spherical_theta(v: Vector3f) -> Float {
+    clamp_t(v.z(), -1.0, 1.0).acos()
+}
+
+pub fn spherical_phi(v: Vector3f) -> Float {
+    let p = v.y().atan2(v.x());
+    if p < 0.0 { p + 2.0 * PI } else { p }
+}

src/core/geometry/primitives.rs

@@ -0,0 +1,880 @@
+use super::traits::{Sqrt, Tuple, VectorLike};
+use super::{Float, NumFloat, PI, clamp_t};
+use crate::utils::interval::Interval;
+use crate::utils::math::{difference_of_products, quadratic, safe_asin};
+use num_traits::{AsPrimitive, FloatConst, Num, Signed, Zero};
+use std::hash::{Hash, Hasher};
+use std::iter::Sum;
+use std::ops::{
+    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
+};
+
+// N-dimensional displacement
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Vector<T, const N: usize>(pub [T; N]);
+// N-dimensional location
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Point<T, const N: usize>(pub [T; N]);
+// N-dimensional surface normal
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Normal<T, const N: usize>(pub [T; N]);
+
+#[macro_export]
+macro_rules! impl_tuple_core {
+    ($Struct:ident) => {
+        impl<T: Copy, const N: usize> Tuple<T, N> for $Struct<T, N> {
+            #[inline]
+            fn data(&self) -> &[T; N] {
+                &self.0
+            }
+
+            #[inline]
+            fn data_mut(&mut self) -> &mut [T; N] {
+                &mut self.0
+            }
+
+            #[inline]
+            fn from_array(arr: [T; N]) -> Self {
+                Self(arr)
+            }
+        }
+
+        impl<T: Default + Copy, const N: usize> Default for $Struct<T, N> {
+            fn default() -> Self {
+                Self([T::default(); N])
+            }
+        }
+
+        impl<T, const N: usize> $Struct<T, N>
+        where
+            T: Zero + Copy,
+        {
+            #[inline]
+            pub fn zero() -> Self {
+                Self([T::zero(); N])
+            }
+        }
+
+        impl<const N: usize> $Struct<f32, N> {
+            #[inline]
+            pub fn floor(&self) -> $Struct<i32, N> {
+                $Struct(self.0.map(|v| v.floor() as i32))
+            }
+
+            #[inline]
+            pub fn average(&self) -> f32 {
+                let sum: f32 = self.0.iter().sum();
+                sum / (N as f32)
+            }
+        }
+
+        impl<T, const N: usize> $Struct<T, N>
+        where
+            T: Copy + PartialOrd,
+        {
+            #[inline]
+            pub fn min(&self, other: Self) -> Self {
+                let mut out = self.0;
+                for i in 0..N {
+                    if other.0[i] < out[i] {
+                        out[i] = other.0[i];
+                    }
+                }
+                Self(out)
+            }
+
+            #[inline]
+            pub fn max(&self, other: Self) -> Self {
+                let mut out = self.0;
+                for i in 0..N {
+                    if other.0[i] > out[i] {
+                        out[i] = other.0[i]
+                    }
+                }
+                Self(out)
+            }
+
+            #[inline]
+            pub fn max_component_value(&self) -> T {
+                let mut m = self.0[0];
+                for i in 1..N {
+                    if self.0[i] > m {
+                        m = self.0[i];
+                    }
+                }
+                m
+            }
+        }
+
+        impl<T, const N: usize> $Struct<T, N>
+        where
+            T: Copy,
+        {
+            #[inline]
+            pub fn fill(value: T) -> Self {
+                Self([value; N])
+            }
+
+            #[inline]
+            pub fn cast<U>(&self) -> $Struct<U, N>
+            where
+                U: 'static + Copy,
+                T: 'static + Copy + AsPrimitive<U>,
+            {
+                $Struct(self.0.map(|c| c.as_()))
+            }
+        }
+
+        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]
+            }
+        }
+        impl<T, const N: usize> IndexMut<usize> for $Struct<T, N> {
+            #[inline]
+            fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+                &mut self.0[index]
+            }
+        }
+
+        impl<T, const N: usize> Neg for $Struct<T, N>
+        where
+            T: Neg<Output = T> + Copy,
+        {
+            type Output = Self;
+            fn neg(self) -> Self::Output {
+                Self(self.0.map(|c| -c))
+            }
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! impl_scalar_ops {
+    ($Struct:ident) => {
+        impl<T, const N: usize> Mul<T> for $Struct<T, N>
+        where
+            T: Mul<Output = T> + Copy,
+        {
+            type Output = Self;
+            fn mul(self, rhs: T) -> Self::Output {
+                let mut result = self.0;
+                for i in 0..N {
+                    result[i] = result[i] * rhs;
+                }
+                Self(result)
+            }
+        }
+
+        impl<const N: usize> Mul<$Struct<Float, N>> for Float {
+            type Output = $Struct<Float, N>;
+            fn mul(self, rhs: $Struct<Float, N>) -> Self::Output {
+                rhs * self
+            }
+        }
+
+        impl<T, const N: usize> MulAssign<T> for $Struct<T, N>
+        where
+            T: MulAssign + Copy,
+        {
+            fn mul_assign(&mut self, rhs: T) {
+                for i in 0..N {
+                    self.0[i] *= rhs;
+                }
+            }
+        }
+
+        impl<T, const N: usize> Div<T> for $Struct<T, N>
+        where
+            T: Div<Output = T> + Copy,
+        {
+            type Output = Self;
+            fn div(self, rhs: T) -> Self::Output {
+                let mut result = self.0;
+                for i in 0..N {
+                    result[i] = result[i] / rhs;
+                }
+                Self(result)
+            }
+        }
+        impl<T, const N: usize> DivAssign<T> for $Struct<T, N>
+        where
+            T: DivAssign + Copy,
+        {
+            fn div_assign(&mut self, rhs: T) {
+                for i in 0..N {
+                    self.0[i] /= rhs;
+                }
+            }
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! impl_op {
+    ($Op:ident, $op:ident, $Lhs:ident, $Rhs:ident, $Output:ident) => {
+        impl<T, const N: usize> $Op<$Rhs<T, N>> for $Lhs<T, N>
+        where
+            T: $Op<Output = T> + Copy,
+        {
+            type Output = $Output<T, N>;
+            fn $op(self, rhs: $Rhs<T, N>) -> Self::Output {
+                let mut result = self.0;
+                for i in 0..N {
+                    result[i] = $Op::$op(self.0[i], rhs.0[i]);
+                }
+                $Output(result)
+            }
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! impl_op_assign {
+    ($OpAssign:ident, $op_assign:ident, $Lhs:ident, $Rhs:ident) => {
+        impl<T, const N: usize> $OpAssign<$Rhs<T, N>> for $Lhs<T, N>
+        where
+            T: $OpAssign + Copy,
+        {
+            fn $op_assign(&mut self, rhs: $Rhs<T, N>) {
+                for i in 0..N {
+                    $OpAssign::$op_assign(&mut self.0[i], rhs.0[i]);
+                }
+            }
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! impl_float_vector_ops {
+    ($Struct:ident) => {
+        impl<T, const N: usize> VectorLike for $Struct<T, N>
+        where
+            T: Copy
+                + Zero
+                + Add<Output = T>
+                + Mul<Output = T>
+                + Sub<Output = T>
+                + Div<Output = T>
+                + Sqrt,
+        {
+            type Scalar = T;
+            fn dot(self, rhs: Self) -> T {
+                let mut sum = T::zero();
+                for i in 0..N {
+                    sum = sum + self[i] * rhs[i];
+                }
+                sum
+            }
+        }
+    };
+}
+
+macro_rules! impl_abs {
+    ($Struct:ident) => {
+        impl<T, const N: usize> $Struct<T, N>
+        where
+            T: Signed + Copy,
+        {
+            pub fn abs(self) -> Self {
+                let mut result = self.0;
+                for i in 0..N {
+                    result[i] = result[i].abs();
+                }
+                Self(result)
+            }
+        }
+    };
+}
+
+macro_rules! impl_accessors {
+    ($Struct:ident) => {
+        impl<T: Copy> $Struct<T, 2> {
+            pub fn x(&self) -> T {
+                self.0[0]
+            }
+            pub fn y(&self) -> T {
+                self.0[1]
+            }
+        }
+        impl<T: Copy> $Struct<T, 3> {
+            pub fn x(&self) -> T {
+                self.0[0]
+            }
+            pub fn y(&self) -> T {
+                self.0[1]
+            }
+            pub fn z(&self) -> T {
+                self.0[2]
+            }
+        }
+    };
+}
+
+impl_tuple_core!(Vector);
+impl_tuple_core!(Point);
+impl_tuple_core!(Normal);
+
+impl_scalar_ops!(Vector);
+impl_scalar_ops!(Normal);
+
+// Addition
+impl_op!(Add, add, Vector, Vector, Vector);
+impl_op!(Add, add, Point, Vector, Point);
+impl_op!(Add, add, Vector, Point, Point);
+impl_op!(Add, add, Normal, Normal, Normal);
+
+// Subtraction
+impl_op!(Sub, sub, Vector, Vector, Vector);
+impl_op!(Sub, sub, Point, Vector, Point);
+impl_op!(Sub, sub, Point, Point, Vector);
+impl_op!(Sub, sub, Normal, Normal, Normal);
+
+// AddAssign
+impl_op_assign!(AddAssign, add_assign, Vector, Vector);
+impl_op_assign!(AddAssign, add_assign, Point, Vector);
+impl_op_assign!(AddAssign, add_assign, Normal, Normal);
+
+// SubAssign
+impl_op_assign!(SubAssign, sub_assign, Vector, Vector);
+impl_op_assign!(SubAssign, sub_assign, Point, Vector);
+impl_op_assign!(SubAssign, sub_assign, Normal, Normal);
+
+impl_float_vector_ops!(Vector);
+impl_float_vector_ops!(Normal);
+impl_abs!(Vector);
+impl_abs!(Normal);
+impl_abs!(Point);
+impl_accessors!(Vector);
+impl_accessors!(Point);
+impl_accessors!(Normal);
+
+impl<T: Copy, const N: usize> From<Vector<T, N>> for Normal<T, N> {
+    fn from(v: Vector<T, N>) -> Self {
+        Self(v.0)
+    }
+}
+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)
+    }
+}
+
+impl<T, const N: usize> Point<T, N>
+where
+    T: NumFloat + Sqrt,
+{
+    pub fn distance(self, other: Self) -> T {
+        (self - other).norm()
+    }
+
+    pub fn distance_squared(self, other: Self) -> T {
+        (self - other).norm_squared()
+    }
+}
+
+impl Point2f {
+    pub fn invert_bilinear(p: Point2f, vert: &[Point2f]) -> Point2f {
+        let a = vert[0];
+        let b = vert[1];
+        let c = vert[3];
+        let d = vert[2];
+        let e = b - a;
+        let f = d - a;
+        let g = (a - b) + (c - d);
+        let h = p - a;
+
+        let cross2d = |a: Vector2f, b: Vector2f| difference_of_products(a.x(), b.y(), a.y(), b.x());
+
+        let k2 = cross2d(g, f);
+        let k1 = cross2d(e, f) + cross2d(h, g);
+        let k0 = cross2d(h, e);
+
+        // if edges are parallel, this is a linear equation
+        if k2.abs() < 0.001 {
+            if (e.x() * k1 - g.x() * k0).abs() < 1e-5 {
+                return Point2f::new(
+                    (h.y() * k1 + f.y() * k0) / (e.y() * k1 - g.y() * k0),
+                    -k0 / k1,
+                );
+            } else {
+                return Point2f::new(
+                    (h.x() * k1 + f.x() * k0) / (e.x() * k1 - g.x() * k0),
+                    -k0 / k1,
+                );
+            }
+        }
+
+        if let Some((v0, v1)) = quadratic(k2, k1, k0) {
+            let u = (h.x() - f.x() * v0) / (e.x() + g.x() * v0);
+            if !(0.0..=1.).contains(&u) || !(0.0..=1.0).contains(&v0) {
+                return Point2f::new((h.x() - f.x() * v1) / (e.x() + g.x() * v1), v1);
+            }
+            Point2f::new(u, v0)
+        } else {
+            Point2f::zero()
+        }
+    }
+}
+
+// Utility aliases and functions
+pub type Point2<T> = Point<T, 2>;
+pub type Point2f = Point2<Float>;
+pub type Point2i = Point2<i32>;
+pub type Point2fi = Point2<Interval>;
+pub type Point3<T> = Point<T, 3>;
+pub type Point3f = Point3<Float>;
+pub type Point3i = Point3<i32>;
+pub type Point3fi = Point3<Interval>;
+pub type Vector2<T> = Vector<T, 2>;
+pub type Vector2f = Vector2<Float>;
+pub type Vector2i = Vector2<i32>;
+pub type Vector2fi = Vector2<Interval>;
+pub type Vector3<T> = Vector<T, 3>;
+pub type Vector3f = Vector3<Float>;
+pub type Vector3i = Vector3<i32>;
+pub type Vector3fi = Vector3<Interval>;
+pub type Normal3<T> = Normal<T, 3>;
+pub type Normal3f = Normal3<Float>;
+pub type Normal3i = Normal3<i32>;
+
+impl<T: Copy> Vector2<T> {
+    pub fn new(x: T, y: T) -> Self {
+        Self([x, y])
+    }
+}
+impl<T: Copy> Point2<T> {
+    pub fn new(x: T, y: T) -> Self {
+        Self([x, y])
+    }
+}
+impl<T: Copy> Vector3<T> {
+    pub fn new(x: T, y: T, z: T) -> Self {
+        Self([x, y, z])
+    }
+}
+impl<T: Copy> Point3<T> {
+    pub fn new(x: T, y: T, z: T) -> Self {
+        Self([x, y, z])
+    }
+}
+impl<T: Copy> Normal3<T> {
+    pub fn new(x: T, y: T, z: T) -> Self {
+        Self([x, y, z])
+    }
+}
+
+// Vector operations
+impl<T> Vector3<T>
+where
+    T: Num + Copy + Neg<Output = T>,
+{
+    pub fn cross(self, rhs: Self) -> Self {
+        Self([
+            self[1] * rhs[2] - self[2] * rhs[1],
+            self[2] * rhs[0] - self[0] * rhs[2],
+            self[0] * rhs[1] - self[1] * rhs[0],
+        ])
+    }
+}
+
+impl<T> Normal3<T>
+where
+    T: Num + Copy + Neg<Output = T>,
+{
+    pub fn cross(self, rhs: Self) -> Self {
+        Self([
+            self[1] * rhs[2] - self[2] * rhs[1],
+            self[2] * rhs[0] - self[0] * rhs[2],
+            self[0] * rhs[1] - self[1] * rhs[0],
+        ])
+    }
+}
+
+impl<T> Vector3<T>
+where
+    T: Num + NumFloat + Copy + Neg<Output = T>,
+{
+    pub fn coordinate_system(&self) -> (Self, Self)
+    where
+        T: NumFloat,
+    {
+        let v2 = if self[0].abs() > self[1].abs() {
+            Self::new(-self[2], T::zero(), self[0]) / (self[0] * self[0] + self[2] * self[2]).sqrt()
+        } else {
+            Self::new(T::zero(), self[2], -self[1]) / (self[1] * self[1] + self[2] * self[2]).sqrt()
+        };
+        (v2, self.cross(v2))
+    }
+
+    pub fn coordinate_system_from_cpp(&self) -> (Self, Self) {
+        let sign = self.z().copysign(T::one());
+        let a = -T::one() / (sign + self.z());
+        let b = self.x() * self.y() * a;
+        let v2 = Self::new(
+            T::one() + sign * self.x().powi(2) * a,
+            sign * b,
+            -sign * self.x(),
+        );
+        let v3 = Self::new(b, sign + self.y().powi(2) * a, -self.y());
+        (v2, v3)
+    }
+}
+
+impl<T> Normal3<T>
+where
+    T: Num + NumFloat + Copy + Neg<Output = T>,
+{
+    pub fn coordinate_system(&self) -> (Self, Self)
+    where
+        T: NumFloat,
+    {
+        let v2 = if self[0].abs() > self[1].abs() {
+            Self::new(-self[2], T::zero(), self[0]) / (self[0] * self[0] + self[2] * self[2]).sqrt()
+        } else {
+            Self::new(T::zero(), self[2], -self[1]) / (self[1] * self[1] + self[2] * self[2]).sqrt()
+        };
+        (v2, self.cross(v2))
+    }
+
+    pub fn coordinate_system_from_cpp(&self) -> (Self, Self) {
+        let sign = self.z().copysign(T::one());
+        let a = -T::one() / (sign + self.z());
+        let b = self.x() * self.y() * a;
+        let v2 = Self::new(
+            T::one() + sign * self.x().powi(2) * a,
+            sign * b,
+            -sign * self.x(),
+        );
+        let v3 = Self::new(b, sign + self.y().powi(2) * a, -self.y());
+        (v2, v3)
+    }
+}
+
+impl<const N: usize> Hash for Vector<Float, N> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        for item in self.0.iter() {
+            item.to_bits().hash(state);
+        }
+    }
+}
+
+impl<const N: usize> Hash for Point<Float, N> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        for item in self.0.iter() {
+            item.to_bits().hash(state);
+        }
+    }
+}
+
+impl<const N: usize> Hash for Normal<Float, N> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        for item in self.0.iter() {
+            item.to_bits().hash(state);
+        }
+    }
+}
+
+// INTERVAL STUFF
+impl<const N: usize> Point<Interval, N> {
+    pub fn new_from_point(p: Point<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new(p[i]);
+        }
+        Self(arr)
+    }
+
+    pub fn new_with_error(p: Point<Float, N>, e: Vector<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new_from_value_and_error(p[i], e[i]);
+        }
+        Self(arr)
+    }
+
+    pub fn error(&self) -> Vector<Float, N> {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = self[i].width() / 2.0;
+        }
+        Vector(arr)
+    }
+
+    pub fn midpoint(&self) -> Point<Float, N> {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = self[i].midpoint();
+        }
+        Point(arr)
+    }
+
+    pub fn is_exact(&self) -> bool {
+        self.0.iter().all(|interval| interval.width() == 0.0)
+    }
+}
+
+impl<const N: usize> Vector<Interval, N> {
+    pub fn new_from_vector(v: Vector<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new(v[i]);
+        }
+        Self(arr)
+    }
+
+    pub fn new_with_error(v: Vector<Float, N>, e: Vector<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new_from_value_and_error(v[i], e[i]);
+        }
+        Self(arr)
+    }
+
+    pub fn error(&self) -> Vector<Float, N> {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = self[i].width() / 2.0;
+        }
+        Vector(arr)
+    }
+
+    pub fn midpoint(&self) -> Vector<Float, N> {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = self[i].midpoint();
+        }
+        Vector(arr)
+    }
+
+    pub fn is_exact(&self) -> bool {
+        self.0.iter().all(|interval| interval.width() == 0.0)
+    }
+}
+
+impl<const N: usize> From<Point<Interval, N>> for Point<Float, N> {
+    fn from(pi: Point<Interval, N>) -> Self {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = pi[i].midpoint();
+        }
+        Point(arr)
+    }
+}
+
+impl<const N: usize> From<Vector<Interval, N>> for Vector<Float, N> {
+    fn from(pi: Vector<Interval, N>) -> Self {
+        let mut arr = [0.0; N];
+        for i in 0..N {
+            arr[i] = pi[i].midpoint();
+        }
+        Vector(arr)
+    }
+}
+
+impl<const N: usize> From<Vector<Float, N>> for Vector<Interval, N> {
+    fn from(v: Vector<Float, N>) -> Self {
+        let mut arr = [Interval::default(); N];
+        for i in 0..N {
+            arr[i] = Interval::new(v[i]);
+        }
+        Self(arr)
+    }
+}
+
+impl<const N: usize> Mul<Vector<Interval, N>> for Interval {
+    type Output = Vector<Interval, N>;
+    fn mul(self, rhs: Vector<Interval, N>) -> Self::Output {
+        rhs * self
+    }
+}
+
+impl<const N: usize> Div<Vector<Interval, N>> for Interval {
+    type Output = Vector<Interval, N>;
+    fn div(self, rhs: Vector<Interval, N>) -> Self::Output {
+        let mut result = rhs.0;
+        for i in 0..N {
+            result[i] = self / rhs[i];
+        }
+        Vector(result)
+    }
+}
+
+impl<const N: usize> From<Vector<i32, N>> for Vector<f32, N> {
+    fn from(v: Vector<i32, N>) -> Self {
+        Self(v.0.map(|c| c as f32))
+    }
+}
+
+impl<const N: usize> From<Point<i32, N>> for Point<Float, N> {
+    fn from(p: Point<i32, N>) -> Self {
+        Point(p.0.map(|c| c as Float))
+    }
+}
+
+impl<T> Normal3<T>
+where
+    T: Num + PartialOrd + Copy + Neg<Output = T> + Sqrt,
+{
+    pub fn face_forward(self, v: Vector3<T>) -> Self {
+        if Vector3::<T>::from(self).dot(v) < T::zero() {
+            -self
+        } else {
+            self
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
+#[repr(C)]
+pub struct OctahedralVector {
+    x: u16,
+    y: u16,
+}
+
+impl OctahedralVector {
+    pub fn new(mut v: Vector3f) -> Self {
+        v /= v.x().abs() + v.y().abs() + v.z().abs();
+
+        let (x_enc, y_enc) = if v.z() >= 0.0 {
+            (Self::encode(v.x()), Self::encode(v.y()))
+        } else {
+            (
+                Self::encode((1.0 - v.y().abs()) * Self::sign(v.x())),
+                Self::encode((1.0 - v.x().abs()) * Self::sign(v.y())),
+            )
+        };
+
+        Self { x: x_enc, y: y_enc }
+    }
+
+    pub fn to_vector(self) -> Vector3f {
+        let mut v = Vector3f::default();
+
+        // Map [0, 65535] back to [-1, 1]
+        v[0] = -1.0 + 2.0 * (self.x as Float / 65535.0);
+        v[1] = -1.0 + 2.0 * (self.y as Float / 65535.0);
+        v[2] = 1.0 - (v.x().abs() + v.y().abs());
+
+        if v.z() < 0.0 {
+            let xo = v.x();
+            v[0] = (1.0 - v.y().abs()) * Self::sign(xo);
+            v[1] = (1.0 - xo.abs()) * Self::sign(v.y());
+        }
+
+        v.normalize()
+    }
+
+    #[inline]
+    pub fn sign(v: Float) -> Float {
+        1.0.copysign(v)
+    }
+
+    #[inline]
+    pub fn encode(f: Float) -> u16 {
+        (clamp_t((f + 1.0) / 2.0, 0.0, 1.0) * 65535.0).round() as u16
+    }
+}
+
+impl From<Vector3f> for OctahedralVector {
+    fn from(v: Vector3f) -> Self {
+        Self::new(v)
+    }
+}
+
+impl From<OctahedralVector> for Vector3f {
+    fn from(ov: OctahedralVector) -> Self {
+        ov.to_vector()
+    }
+}
+
+#[derive(Copy, Clone, Debug, Default, PartialEq)]
+pub struct Frame {
+    pub x: Vector3f,
+    pub y: Vector3f,
+    pub z: Vector3f,
+}
+
+impl Frame {
+    pub fn new(x: Vector3f, z: Vector3f) -> Self {
+        Self {
+            x,
+            y: z.cross(x),
+            z,
+        }
+    }
+
+    pub fn from_x(x: Vector3f) -> Self {
+        let (y, z) = x.normalize().coordinate_system();
+        Self {
+            x: x.normalize(),
+            y,
+            z,
+        }
+    }
+
+    pub fn from_xz(x: Vector3f, z: Vector3f) -> Self {
+        Self {
+            x,
+            y: z.cross(x),
+            z,
+        }
+    }
+
+    pub fn from_xy(x: Vector3f, y: Vector3f) -> Self {
+        Self {
+            x,
+            y,
+            z: x.cross(y),
+        }
+    }
+
+    pub fn from_y(y: Vector3f) -> Self {
+        let (z, x) = y.normalize().coordinate_system();
+        Self {
+            x,
+            y: y.normalize(),
+            z,
+        }
+    }
+
+    pub fn from_z(z: Vector3f) -> Self {
+        let (x, y) = z.normalize().coordinate_system();
+        Self {
+            x,
+            y,
+            z: z.normalize(),
+        }
+    }
+
+    pub fn to_local(&self, v: Vector3f) -> Vector3f {
+        Vector3f::new(v.dot(self.x), v.dot(self.y), v.dot(self.z))
+    }
+
+    pub fn to_local_normal(&self, n: Normal3f) -> Normal3f {
+        let n: Vector3f = n.into();
+        Normal3f::new(n.dot(self.x), n.dot(self.y), n.dot(self.z))
+    }
+
+    pub fn from_local(&self, v: Vector3f) -> Vector3f {
+        self.x * v.x() + self.y * v.y() + self.z * v.z()
+    }
+
+    pub fn from_local_normal(&self, v: Normal3f) -> Normal3f {
+        Normal3f::from(self.x * v.x() + self.y * v.y() + self.z * v.z())
+    }
+}

src/core/geometry/ray.rs

@@ -0,0 +1,119 @@
+use super::{Normal3f, Point3f, Point3fi, Vector3f, VectorLike};
+use crate::core::medium::Medium;
+use crate::core::pbrt::Float;
+use crate::utils::math::{next_float_down, next_float_up};
+use std::sync::Arc;
+
+#[derive(Clone, Debug)]
+pub struct Ray {
+    pub o: Point3f,
+    pub d: Vector3f,
+    pub medium: Option<Arc<Medium>>,
+    pub time: Float,
+    // We do this instead of creating a trait for Rayable or some gnarly thing like that
+    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 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 at(&self, t: Float) -> Point3f {
+        self.o + self.d * t
+    }
+
+    pub fn offset_origin(p: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f {
+        let d: Float = Vector3f::from(n.abs()).dot(p.error());
+        let normal: Vector3f = Vector3f::from(*n);
+
+        let mut offset = p.midpoint();
+        if w.dot(normal) < 0.0 {
+            offset -= normal * d;
+        } else {
+            offset += normal * d;
+        }
+
+        for i in 0..3 {
+            if n[i] > 0.0 {
+                offset[i] = next_float_up(offset[i]);
+            } else if n[i] < 0.0 {
+                offset[i] = next_float_down(offset[i]);
+            }
+        }
+        offset
+    }
+
+    pub fn spawn(pi: &Point3fi, n: &Normal3f, time: Float, d: Vector3f) -> Ray {
+        let origin = Self::offset_origin(pi, n, &d);
+        Ray {
+            o: origin,
+            d,
+            time,
+            medium: None,
+            differential: None,
+        }
+    }
+
+    pub fn spawn_to_point(p_from: &Point3fi, n: &Normal3f, time: Float, p_to: Point3f) -> Ray {
+        let d = p_to - p_from.midpoint();
+        Self::spawn(p_from, n, time, d)
+    }
+
+    pub fn spawn_to_interaction(
+        p_from: &Point3fi,
+        n_from: &Normal3f,
+        time: Float,
+        p_to: &Point3fi,
+        n_to: &Normal3f,
+    ) -> Ray {
+        let dir_for_offset = p_to.midpoint() - p_from.midpoint();
+        let pf = Self::offset_origin(p_from, n_from, &dir_for_offset);
+        let pt = Self::offset_origin(p_to, n_to, &(pf - p_to.midpoint()));
+
+        let d = pt - pf;
+
+        Ray {
+            o: pf,
+            d,
+            time,
+            medium: None,
+            differential: None,
+        }
+    }
+
+    pub fn scale_differentials(&mut self, s: Float) {
+        if let Some(differential) = &mut self.differential {
+            differential.rx_origin = self.o + (differential.rx_origin - self.o) * s;
+            differential.ry_origin = self.o + (differential.ry_origin - self.o) * s;
+            differential.rx_direction = self.d + (differential.rx_direction - self.d) * s;
+            differential.ry_direction = self.d + (differential.ry_direction - self.d) * s;
+        }
+    }
+}
+
+#[derive(Debug, Default, Copy, Clone)]
+pub struct RayDifferential {
+    pub rx_origin: Point3f,
+    pub ry_origin: Point3f,
+    pub rx_direction: Vector3f,
+    pub ry_direction: Vector3f,
+}

src/core/geometry/traits.rs

@@ -0,0 +1,183 @@
+use crate::core::pbrt::Float;
+use crate::utils::interval::Interval;
+use crate::utils::math::{next_float_down, next_float_up};
+use num_integer::Roots;
+use num_traits::{Float as NumFloat, FloatConst, Num, One, Signed, Zero};
+use std::ops::{Add, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub};
+
+pub trait Tuple<T, const N: usize>:
+    Sized + Copy + Index<usize, Output = T> + IndexMut<usize>
+{
+    fn data(&self) -> &[T; N];
+
+    fn data_mut(&mut self) -> &mut [T; N];
+
+    fn from_array(arr: [T; N]) -> Self;
+
+    #[inline]
+    fn permute(&self, p: [usize; N]) -> Self
+    where
+        T: Copy,
+    {
+        let new_data = p.map(|index| self[index]);
+        Self::from_array(new_data)
+    }
+
+    fn max_component_value(&self) -> T
+    where
+        T: PartialOrd + Copy,
+    {
+        self.data()
+            .iter()
+            .copied()
+            .reduce(|a, b| if a > b { a } else { b })
+            .expect("Cannot get max component of a zero-length tuple")
+    }
+
+    fn min_component_value(&self) -> T
+    where
+        T: PartialOrd + Copy,
+    {
+        self.data()
+            .iter()
+            .copied()
+            .reduce(|a, b| if a < b { a } else { b })
+            .expect("Cannot get min component of a zero-length tuple")
+    }
+
+    fn max_component_index(&self) -> usize
+    where
+        T: PartialOrd,
+    {
+        self.data()
+            .iter()
+            .enumerate()
+            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
+            .map(|(index, _)| index)
+            .unwrap_or(0)
+    }
+
+    fn min_component_index(&self) -> usize
+    where
+        T: PartialOrd,
+    {
+        self.data()
+            .iter()
+            .enumerate()
+            .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
+            .map(|(index, _)| index)
+            .unwrap_or(0)
+    }
+}
+
+pub trait VectorLike:
+    Sized
+    + Copy
+    + Add<Output = Self>
+    + Sub<Output = Self>
+    + Div<Self::Scalar, Output = Self>
+    + Mul<Self::Scalar, Output = Self>
+{
+    type Scalar: Copy + Zero + Add<Output = Self::Scalar> + Mul<Output = Self::Scalar> + Sqrt;
+
+    fn dot(self, rhs: Self) -> Self::Scalar;
+    fn norm_squared(self) -> Self::Scalar {
+        self.dot(self)
+    }
+
+    fn abs_dot(self, rhs: Self) -> Self::Scalar
+    where
+        Self::Scalar: Signed,
+    {
+        self.dot(rhs).abs()
+    }
+
+    fn gram_schmidt(self, rhs: Self) -> Self {
+        self - rhs * self.dot(rhs)
+    }
+
+    fn norm(&self) -> Self::Scalar {
+        self.norm_squared().sqrt()
+    }
+
+    fn normalize(self) -> Self
+    where
+        Self::Scalar: NumFloat,
+    {
+        let n = self.norm();
+        if n.is_zero() { self } else { self / n }
+    }
+
+    fn angle_between(self, rhs: Self) -> Self::Scalar
+    where
+        Self::Scalar: NumFloat,
+    {
+        let dot_product = self.normalize().dot(rhs.normalize());
+        let clamped_dot = dot_product
+            .min(Self::Scalar::one())
+            .max(-Self::Scalar::one());
+        clamped_dot.acos()
+    }
+}
+
+pub trait Sqrt {
+    fn sqrt(self) -> Self;
+}
+
+impl Sqrt for Float {
+    fn sqrt(self) -> Self {
+        self.sqrt()
+    }
+}
+
+impl Sqrt for f64 {
+    fn sqrt(self) -> Self {
+        self.sqrt()
+    }
+}
+
+impl Sqrt for i32 {
+    fn sqrt(self) -> Self {
+        self.isqrt()
+    }
+}
+
+impl Sqrt for u32 {
+    fn sqrt(self) -> Self {
+        self.isqrt()
+    }
+}
+
+impl Sqrt for Interval {
+    fn sqrt(self) -> Self {
+        let low = if self.low < 0.0 {
+            0.0
+        } else {
+            next_float_down(self.low.sqrt())
+        };
+
+        let high = if self.high < 0.0 {
+            0.0
+        } else {
+            next_float_up(self.high.sqrt())
+        };
+
+        Self { low, high }
+    }
+}
+
+pub trait Lerp<Factor = Float>: Sized + Copy {
+    fn lerp(t: Factor, a: Self, b: Self) -> Self;
+}
+
+impl<T, F, Diff> Lerp<F> for T
+where
+    T: Copy + Sub<Output = Diff> + Add<Diff, Output = T>,
+    Diff: Mul<F, Output = Diff>,
+    F: Copy,
+{
+    #[inline(always)]
+    fn lerp(t: F, a: Self, b: Self) -> Self {
+        a + (b - a) * t
+    }
+}

src/core/scattering.rs

@@ -0,0 +1,173 @@
+use crate::core::geometry::{
+    Normal3f, Point2f, Vector2f, Vector3f, VectorLike, abs_cos_theta, cos_phi, cos2_theta, sin_phi,
+    tan2_theta,
+};
+use crate::core::pbrt::{Float, PI, clamp_t};
+use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
+use crate::utils::math::safe_sqrt;
+use crate::utils::math::{lerp, square};
+use crate::utils::sampling::sample_uniform_disk_polar;
+
+use num::complex::Complex;
+
+#[derive(Debug, Default, Clone, Copy)]
+pub struct TrowbridgeReitzDistribution {
+    alpha_x: Float,
+    alpha_y: Float,
+}
+
+impl TrowbridgeReitzDistribution {
+    pub fn new(alpha_x: Float, alpha_y: Float) -> Self {
+        Self { alpha_x, alpha_y }
+    }
+
+    pub fn d(&self, wm: Vector3f) -> Float {
+        let tan2_theta = tan2_theta(wm);
+        if tan2_theta.is_infinite() {
+            return 0.;
+        }
+        let cos4_theta = square(cos2_theta(wm));
+        let e =
+            tan2_theta * (square(cos_phi(wm) / self.alpha_x) + square(sin_phi(wm) / self.alpha_y));
+        1.0 / (PI * self.alpha_x * self.alpha_y * cos4_theta * square(1. + e))
+    }
+
+    pub fn effectively_smooth(&self) -> bool {
+        self.alpha_x.max(self.alpha_y) < 1e-3
+    }
+
+    pub fn lambda(&self, w: Vector3f) -> Float {
+        let tan2_theta = tan2_theta(w);
+        if tan2_theta.is_infinite() {
+            return 0.;
+        }
+        let alpha2 = square(cos_phi(w) * self.alpha_x) + square(sin_phi(w) * self.alpha_y);
+        ((1. + alpha2 * tan2_theta).sqrt() - 1.) / 2.
+    }
+
+    pub fn g(&self, wo: Vector3f, wi: Vector3f) -> Float {
+        1. / (1. + self.lambda(wo) + self.lambda(wi))
+    }
+
+    pub fn g1(&self, w: Vector3f) -> Float {
+        1. / (1. / self.lambda(w))
+    }
+
+    pub fn d_from_w(&self, w: Vector3f, wm: Vector3f) -> Float {
+        self.g1(w) / abs_cos_theta(w) * self.d(wm) * w.dot(wm).abs()
+    }
+
+    pub fn pdf(&self, w: Vector3f, wm: Vector3f) -> Float {
+        self.d_from_w(w, wm)
+    }
+
+    pub fn sample_wm(&self, w: Vector3f, u: Point2f) -> Vector3f {
+        let mut wh = Vector3f::new(self.alpha_x * w.x(), self.alpha_y * w.y(), w.z()).normalize();
+        if wh.z() < 0. {
+            wh = -wh;
+        }
+        let t1 = if wh.z() < 0.99999 {
+            Vector3f::new(0., 0., 1.).cross(wh).normalize()
+        } else {
+            Vector3f::new(1., 0., 0.)
+        };
+        let t2 = wh.cross(t1);
+        let mut p = sample_uniform_disk_polar(u);
+        let h = (1. - square(p.x())).sqrt();
+        p[1] = lerp((1. + wh.z()) / 2., h, p.y());
+        let pz = 0_f32.max(1. - Vector2f::from(p).norm_squared());
+        let nh = p.x() * t1 + p.y() * t2 + pz * wh;
+        Vector3f::new(
+            self.alpha_x * nh.x(),
+            self.alpha_y * nh.y(),
+            nh.z().max(1e-6),
+        )
+        .normalize()
+    }
+
+    pub fn roughness_to_alpha(roughness: Float) -> Float {
+        roughness.sqrt()
+    }
+
+    pub fn regularize(&mut self) {
+        if self.alpha_x < 0.3 {
+            self.alpha_x = clamp_t(2. * self.alpha_x, 0.1, 0.3);
+        }
+        if self.alpha_y < 0.3 {
+            self.alpha_y = clamp_t(2. * self.alpha_y, 0.1, 0.3);
+        }
+    }
+}
+
+pub fn refract(wi: Vector3f, n: Normal3f, eta_ratio: Float) -> Option<(Vector3f, Float)> {
+    let mut n_interface = n;
+    let mut eta = eta_ratio;
+
+    let mut cos_theta_i = Vector3f::from(n_interface).dot(wi);
+
+    if cos_theta_i < 0.0 {
+        eta = 1.0 / eta;
+        cos_theta_i = -cos_theta_i;
+        n_interface = -n_interface;
+    }
+
+    let sin2_theta_i = (1.0 - square(cos_theta_i)).max(0.0_f32);
+    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))
+}
+
+pub fn reflect(wo: Vector3f, n: Normal3f) -> Vector3f {
+    -wo + Vector3f::from(2. * wo.dot(n.into()) * n)
+}
+
+pub fn fr_dielectric(cos_theta_i: Float, eta: Float) -> Float {
+    let mut cos_safe = clamp_t(cos_theta_i, -1., 1.);
+    let mut eta_corr = eta;
+    if cos_theta_i < 0. {
+        eta_corr = 1. / eta_corr;
+        cos_safe = -cos_safe;
+    }
+    let sin2_theta_i = 1. - square(cos_safe);
+    let sin2_theta_t = sin2_theta_i / square(eta_corr);
+    if sin2_theta_t >= 1. {
+        return 1.;
+    }
+    let cos_theta_t = safe_sqrt(1. - sin2_theta_t);
+    let r_parl = (eta_corr * cos_safe - cos_theta_t) / (eta_corr * cos_safe + cos_theta_t);
+    let r_perp = (cos_safe - eta_corr * cos_theta_t) / (cos_safe + eta_corr * cos_theta_t);
+
+    (square(r_parl) + square(r_perp)) / 2.
+}
+
+pub fn fr_complex(cos_theta_i: Float, eta: Complex<Float>) -> Float {
+    let cos_corr = clamp_t(cos_theta_i, 0., 1.);
+    let sin2_theta_i = 1. - square(cos_corr);
+    let sin2_theta_t: Complex<Float> = sin2_theta_i / square(eta);
+    let cos2_theta_t: Complex<Float> = (1. - sin2_theta_t).sqrt();
+
+    let r_parl = (eta * cos_corr - cos2_theta_t) / (eta * cos_corr + cos2_theta_t);
+    let r_perp = (cos_corr - eta * cos2_theta_t) / (cos_corr + eta * cos2_theta_t);
+
+    (r_parl.norm() + r_perp.norm()) / 2.
+}
+
+pub fn fr_complex_from_spectrum(
+    cos_theta_i: Float,
+    eta: SampledSpectrum,
+    k: SampledSpectrum,
+) -> SampledSpectrum {
+    let mut result = SampledSpectrum::default();
+    for i in 0..N_SPECTRUM_SAMPLES {
+        result[i] = fr_complex(cos_theta_i, Complex::new(eta[i], k[i]));
+    }
+    result
+}

src/spectra/colorspace.rs

@@ -0,0 +1,107 @@
+use super::color::{RGB, RGBSigmoidPolynomial, RGBToSpectrumTable, XYZ};
+use crate::core::geometry::Point2f;
+use crate::core::pbrt::Float;
+use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, Spectrum};
+use crate::utils::math::SquareMatrix;
+
+use once_cell::sync::Lazy;
+
+use std::cmp::{Eq, PartialEq};
+use std::error::Error;
+use std::sync::Arc;
+
+#[derive(Debug, Clone)]
+pub struct RGBColorSpace {
+    pub r: Point2f,
+    pub g: Point2f,
+    pub b: Point2f,
+    pub w: Point2f,
+    pub illuminant: Spectrum,
+    pub rgb_to_spectrum_table: Arc<RGBToSpectrumTable>,
+    pub xyz_from_rgb: SquareMatrix<Float, 3>,
+    pub rgb_from_xyz: SquareMatrix<Float, 3>,
+}
+
+impl RGBColorSpace {
+    pub fn new(
+        r: Point2f,
+        g: Point2f,
+        b: Point2f,
+        illuminant: Spectrum,
+        rgb_to_spectrum_table: RGBToSpectrumTable,
+    ) -> Result<Self, Box<dyn Error>> {
+        let w_xyz: XYZ = illuminant.to_xyz();
+        let w = w_xyz.xy();
+        let r_xyz = XYZ::from_xyy(r, Some(1.0));
+        let g_xyz = XYZ::from_xyy(g, Some(1.0));
+        let b_xyz = XYZ::from_xyy(b, Some(1.0));
+        let rgb_values = [
+            [r_xyz.x(), g_xyz.x(), b_xyz.x()],
+            [r_xyz.y(), g_xyz.y(), b_xyz.y()],
+            [r_xyz.z(), g_xyz.z(), g_xyz.z()],
+        ];
+        let rgb = SquareMatrix::new(rgb_values);
+        let c: RGB = rgb.inverse()? * w_xyz;
+        let xyz_from_rgb = rgb * SquareMatrix::diag(&[c.r, c.g, c.b]);
+        let rgb_from_xyz = xyz_from_rgb
+            .inverse()
+            .expect("XYZ from RGB matrix is singular");
+
+        Ok(Self {
+            r,
+            g,
+            b,
+            w,
+            illuminant,
+            rgb_to_spectrum_table: Arc::new(rgb_to_spectrum_table),
+            xyz_from_rgb,
+            rgb_from_xyz,
+        })
+    }
+
+    pub fn to_xyz(&self, rgb: RGB) -> XYZ {
+        self.xyz_from_rgb * rgb
+    }
+
+    pub fn to_rgb(&self, xyz: XYZ) -> RGB {
+        self.rgb_from_xyz * xyz
+    }
+
+    pub fn to_rgb_coeffs(&self, rgb: RGB) -> RGBSigmoidPolynomial {
+        self.rgb_to_spectrum_table.to_polynomial(rgb)
+    }
+
+    pub fn convert_colorspace(&self, other: &RGBColorSpace) -> SquareMatrix<Float, 3> {
+        if self == other {
+            return SquareMatrix::default();
+        }
+
+        self.rgb_from_xyz * other.xyz_from_rgb
+    }
+
+    pub fn srgb() -> &'static Self {
+        static SRGB_SPACE: Lazy<RGBColorSpace> = Lazy::new(|| {
+            let r = Point2f::new(0.64, 0.33);
+            let g = Point2f::new(0.30, 0.60);
+            let b = Point2f::new(0.15, 0.06);
+
+            let illuminant = Spectrum::std_illuminant_d65();
+            let table = RGBToSpectrumTable::srgb();
+
+            RGBColorSpace::new(r, g, b, illuminant, table)
+                .expect("Failed to initialize standard sRGB color space")
+        });
+
+        &SRGB_SPACE
+    }
+}
+
+impl PartialEq for RGBColorSpace {
+    fn eq(&self, other: &Self) -> bool {
+        self.r == other.r
+            && self.g == other.g
+            && self.b == other.b
+            && self.w == other.w
+            && Arc::ptr_eq(&self.rgb_to_spectrum_table, &other.rgb_to_spectrum_table)
+    }
+}