From 21a2c0c6740835063c80a596aaab7401889e1793 Mon Sep 17 00:00:00 2001 From: pingupingou Date: Mon, 10 Nov 2025 21:27:12 +0000 Subject: [PATCH] Fixed Traits in geometry, refactored geometry into module, still cant get my head around the exr image stuff, C++ translation is not that linear --- src/geometry/bounds.rs | 272 ++++++++++++++++ src/geometry/cone.rs | 117 +++++++ src/geometry/mod.rs | 100 ++++++ src/geometry/primitives.rs | 632 +++++++++++++++++++++++++++++++++++++ src/geometry/ray.rs | 118 +++++++ src/geometry/traits.rs | 58 ++++ 6 files changed, 1297 insertions(+) create mode 100644 src/geometry/bounds.rs create mode 100644 src/geometry/cone.rs create mode 100644 src/geometry/mod.rs create mode 100644 src/geometry/primitives.rs create mode 100644 src/geometry/ray.rs create mode 100644 src/geometry/traits.rs diff --git a/src/geometry/bounds.rs b/src/geometry/bounds.rs new file mode 100644 index 0000000..d9805a0 --- /dev/null +++ b/src/geometry/bounds.rs @@ -0,0 +1,272 @@ +use super::{Float, NumFloat}; +use super::{Point, Point3, Point2f, Point3f, Vector, Vector2, Vector2f, Vector3, Vector3f}; +use crate::geometry::traits::VectorLike; +use crate::core::pbrt::lerp; +use crate::geometry::{max, min}; +use crate::utils::interval::Interval; +use num_traits::{Num, Bounded}; +use std::mem; +use std::ops::{Add, Div, DivAssign, Mul, Sub}; + +// AABB BOUNDING BOXES + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct Bounds { + pub p_min: Point, + pub p_max: Point, +} + +impl<'a, T, const N: usize> IntoIterator for &'a Bounds { + type Item = &'a Point; + type IntoIter = std::array::IntoIter<&'a Point, 2>; + + fn into_iter(self) -> Self::IntoIter { + [&self.p_min, &self.p_max].into_iter() + } +} + +impl Bounds +where + T: Num + PartialOrd + Copy, +{ + pub fn from_point(p: Point) -> Self { + Self { p_min: p, p_max: p } + } + + pub fn from_points(p1: Point, p2: Point) -> 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) -> 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 { + self.p_max - self.p_min + } + + 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) -> Point { + 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: Sub>, + { + 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) -> Vector + where + Point: Sub>, + Vector: DivAssign, + { + 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 contains(&self, p: Point) -> bool { + (0..N).all(|i| p[i] >= self.p_min[i] && p[i] <= self.p_max[i]) + } + + pub fn contains_exclusive(&self, p: Point) -> 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 Default for Bounds +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 = Bounds; +pub type Bounds2f = Bounds2; +pub type Bounds2i = Bounds2; +pub type Bounds2fi = Bounds2; +pub type Bounds3 = Bounds; +pub type Bounds3i = Bounds3; +pub type Bounds3f = Bounds3; +pub type Bounds3fi = Bounds3; + +impl Bounds3 +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 Bounds3 +where + T: NumFloat + PartialOrd + Copy + Default +{ + pub fn bounding_sphere(&self) -> (Point3, 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, d: Vector3, 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 Bounds2 +where + T: Num + Copy + Default, +{ + pub fn area(&self) -> T { + let d: Vector2 = self.p_max - self.p_min; + d.x() * d.y() + } +} + +impl Bounds3f { + 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) + } +} diff --git a/src/geometry/cone.rs b/src/geometry/cone.rs new file mode 100644 index 0000000..c0b3c60 --- /dev/null +++ b/src/geometry/cone.rs @@ -0,0 +1,117 @@ +use super::{Bounds3f, Float, PI, Point3f, Vector3f, VectorLike}; +use crate::core::pbrt::square; +use crate::utils::math::{degrees, safe_acos, safe_asin, safe_sqrt}; +use crate::utils::transform::Transform; + +#[derive(Debug, Clone)] +pub struct DirectionCone { + w: Vector3f, + 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()) +} diff --git a/src/geometry/mod.rs b/src/geometry/mod.rs new file mode 100644 index 0000000..d32445d --- /dev/null +++ b/src/geometry/mod.rs @@ -0,0 +1,100 @@ +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::primitives::{ + Frame, Normal, Normal3f, Point, Point2f, Point2i, Point2fi, Point3, Point3i, Point3f, Point3fi, Vector, Vector2, Vector2f, + Vector2fi, Vector2i, Vector3, Vector3f, Vector3fi, Vector3i +}; +pub use self::cone::DirectionCone; +pub use self::ray::{Ray, RayDifferential}; +pub use self::traits::{FloatVectorLike, VectorLike}; + +use crate::core::pbrt::{Float, PI, clamp_t, square}; +use num_traits::Float as NumFloat; + +#[inline] +pub fn min(a: T, b: T) -> T { + if a < b { a } else { b } +} + +#[inline] +pub fn max(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); + let result = if sin_theta == 0. { + 1. + } else { + clamp_t(w.x() / sin_theta, -1., 1.) + }; + result +} + +#[inline] +pub fn sin_phi(w: Vector3f) -> Float { + let sin_theta = sin_theta(w); + let result = if sin_theta == 0. { + 0. + } else { + clamp_t(w.y() / sin_theta, -1., 1.) + }; + result +} + +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_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 } +} + + diff --git a/src/geometry/primitives.rs b/src/geometry/primitives.rs new file mode 100644 index 0000000..357d71e --- /dev/null +++ b/src/geometry/primitives.rs @@ -0,0 +1,632 @@ +use super::traits::Tuple; +use super::{Float, NumFloat, PI}; +use crate::utils::interval::Interval; +use crate::utils::math::safe_asin; +use num_traits::{Num, Signed, Zero, FloatConst}; +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, Hash)] +pub struct Vector(pub [T; N]); +// N-dimensional location +#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)] +pub struct Point(pub [T; N]); +// N-dimensional surface normal +#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)] +pub struct Normal(pub [T; N]); + +#[macro_export] +macro_rules! impl_tuple_core { + ($Struct:ident) => { + impl Tuple for $Struct { + #[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 Default for $Struct { + fn default() -> Self { + Self([T::default(); N]) + } + } + + impl $Struct + where + T: Zero + Copy, + { + #[inline] + pub fn zero() -> Self { + Self([T::zero(); N]) + } + } + + impl $Struct { + #[inline] + pub fn floor(&self) -> $Struct { + $Struct(self.0.map(|v| v.floor() as i32)) + } + } + + impl $Struct + where + T: Copy, + { + #[inline] + pub fn fill(value: T) -> Self { + Self([value; N]) + } + } + + impl Index for $Struct { + type Output = T; + #[inline] + fn index(&self, index: usize) -> &Self::Output { + &self.0[index] + } + } + impl IndexMut for $Struct { + #[inline] + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + &mut self.0[index] + } + } + + impl Neg for $Struct + where + T: Neg + 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 Mul for $Struct + where + T: Mul + 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 Mul<$Struct> for Float { + type Output = $Struct; + fn mul(self, rhs: $Struct) -> Self::Output { + rhs * self + } + } + + impl MulAssign for $Struct + where + T: MulAssign + Copy, + { + fn mul_assign(&mut self, rhs: T) { + for i in 0..N { + self.0[i] *= rhs; + } + } + } + + impl Div for $Struct + where + T: Div + 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 DivAssign for $Struct + 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 $Op<$Rhs> for $Lhs + where + T: $Op + Copy, + { + type Output = $Output; + fn $op(self, rhs: $Rhs) -> 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 $OpAssign<$Rhs> for $Lhs + where + T: $OpAssign + Copy, + { + fn $op_assign(&mut self, rhs: $Rhs) { + 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) => { + // This impl block is constrained to only apply when the scalar type `T` is a float. + impl $Struct + where + T: NumFloat, + { + pub fn dot(self, rhs: Self) -> T { + let mut sum = T::zero(); + for i in 0..N { + sum = sum + self[i] * rhs[i]; + } + sum + } + + pub fn norm_squared(&self) -> T { + let mut sum = T::zero(); + for i in 0..N { + sum = sum + self[i] * self[i]; + } + sum + } + + pub fn norm(&self) -> T { + self.norm_squared().sqrt() + } + + pub fn normalize(self) -> Self { + let n = self.norm(); + if n.is_zero() { + self + } else { + self / n + } + } + + pub fn angle_between(self, rhs: Self) -> T { + let dot_product = self.normalize().dot(rhs.normalize()); + let clamped_dot = dot_product.min(T::one()).max(-T::one()); + clamped_dot.acos() + } + + pub fn project_on(self, rhs: Self) -> Self { + let rhs_norm_sq = rhs.norm_squared(); + if rhs_norm_sq.is_zero() { + // This now calls the inherent `zero()` method from `impl_tuple_core!` + Self::zero() + } else { + // Note: This requires Mul to be implemented for the struct, + // which your `impl_scalar_ops!` macro already does. + rhs * (self.dot(rhs) / rhs_norm_sq) + } + } + } + }; +} + +macro_rules! impl_abs { + ($Struct:ident) => { + impl $Struct + 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 $Struct { + pub fn x(&self) -> T { + self.0[0] + } + pub fn y(&self) -> T { + self.0[1] + } + } + impl $Struct { + 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_accessors!(Vector); +impl_accessors!(Point); +impl_accessors!(Normal); + + +impl From> for Normal { + fn from(v: Vector) -> Self { + Self(v.0) + } +} +impl From> for Vector { + fn from(n: Normal) -> Self { + Self(n.0) + } +} + +impl From> for Point { + fn from(v: Vector) -> Self { + Self(v.0) + } +} + +impl From> for Vector { + fn from(n: Point) -> Self { + Self(n.0) + } +} + +impl Point +where + T: NumFloat, +{ + pub fn distance(self, other: Self) -> T { + (self - other).norm() + } + + pub fn distance_squared(self, other: Self) -> T { + (self - other).norm_squared() + } +} + +// Utility aliases and functions +pub type Point2 = Point; +pub type Point2f = Point2; +pub type Point2i = Point2; +pub type Point2fi = Point2; +pub type Point3 = Point; +pub type Point3f = Point3; +pub type Point3i = Point3; +pub type Point3fi = Point3; +pub type Vector2 = Vector; +pub type Vector2f = Vector2; +pub type Vector2i = Vector2; +pub type Vector2fi = Vector2; +pub type Vector3 = Vector; +pub type Vector3f = Vector3; +pub type Vector3i = Vector3; +pub type Vector3fi = Vector3; +pub type Normal3 = Normal; +pub type Normal3f = Normal3; +pub type Normal3i = Normal3; + +impl Vector2 { + pub fn new(x: T, y: T) -> Self { + Self([x, y]) + } +} +impl Point2 { + pub fn new(x: T, y: T) -> Self { + Self([x, y]) + } +} +impl Vector3 { + pub fn new(x: T, y: T, z: T) -> Self { + Self([x, y, z]) + } +} +impl Point3 { + pub fn new(x: T, y: T, z: T) -> Self { + Self([x, y, z]) + } +} +impl Normal3 { + pub fn new(x: T, y: T, z: T) -> Self { + Self([x, y, z]) + } +} + +// Vector operations +impl Vector3 +where + T: Num + Copy + Neg, +{ + 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 Vector3 +where + T: Num + NumFloat + Copy + Neg, +{ + 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)) + } +} + +impl Point { + pub fn new_from_point(p: Point) -> 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, e: Vector) -> 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 { + 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 { + 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 Vector { + pub fn new_from_vector(v: Vector) -> 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, e: Vector) -> 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 { + 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 { + 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 From> for Point { + fn from(pi: Point) -> Self { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = pi[i].midpoint(); + } + Point(arr) + } +} + +impl Mul> for Interval { + type Output = Vector; + fn mul(self, rhs: Vector) -> Self::Output { + rhs * self + } +} + +impl Div> for Interval { + type Output = Vector; + fn div(self, rhs: Vector) -> Self::Output { + let mut result = rhs.0; + for i in 0..N { + result[i] = self / rhs[i]; + } + Vector(result) + } +} + +impl From> for Vector { + fn from(v: Vector) -> Self { + Self(v.0.map(|c| c as f32)) + } +} + +impl From> for Point { + fn from(p: Point) -> Self { + Point(p.0.map(|c| c as Float)) + } +} + +impl Normal3 +where + T: Num + PartialOrd + Copy + Neg, +{ + pub fn face_forward(self, v: Vector3) -> Self { + if Vector3::::from(self).dot(v.into()) < T::zero() { -self } else { self } + } +} + +#[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_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()) + } +} diff --git a/src/geometry/ray.rs b/src/geometry/ray.rs new file mode 100644 index 0000000..a065e79 --- /dev/null +++ b/src/geometry/ray.rs @@ -0,0 +1,118 @@ +use super::{Normal3f, Point3f, Point3fi, Vector3f, VectorLike}; +use crate::core::medium::Medium; +use crate::core::pbrt::{Float, 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>, + pub time: Float, + // We do this instead of creating a trait for Rayable or some gnarly thing like that + pub differential: Option, +} + +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, medium: Option>) -> 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 + } + + pub fn offset_origin(p: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f { + let d: Float = Vector3f::from(n.abs()).dot(p.error().into()); + 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, +} diff --git a/src/geometry/traits.rs b/src/geometry/traits.rs new file mode 100644 index 0000000..391d375 --- /dev/null +++ b/src/geometry/traits.rs @@ -0,0 +1,58 @@ +use std::ops::{Add, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub}; + +use num_traits::{Float as NumFloat, Num, Zero, One, FloatConst}; + +pub trait Tuple: + Sized + Copy + Index + IndexMut +{ + fn data(&self) -> &[T; N]; + + fn data_mut(&mut self) -> &mut [T; N]; + + fn from_array(arr: [T; N]) -> Self; +} + +pub trait VectorLike: + Sized + + Copy + + Add + + Sub + + Div +{ + type Scalar: NumFloat; + + fn dot(self, rhs: Self) -> Self::Scalar; + fn norm_squared(&self) -> Self::Scalar; + + fn norm(&self) -> Self::Scalar { + self.norm_squared().sqrt() + } + + fn normalize(self) -> Self { + let n = self.norm(); + if n.is_zero() { + self + } else { + self / n + } + } + + fn angle_between(self, rhs: Self) -> Self::Scalar { + let dot_product = self.normalize().dot(rhs.normalize()); + let clamped_dot = dot_product.min(Self::Scalar::one()).max(-Self::Scalar::one()); + clamped_dot.acos() + } + + fn project_on(self, rhs: Self) -> Self { + let rhs_norm_sq = rhs.norm_squared(); + if rhs_norm_sq.is_zero() { + Self::zero() + } else { + let scale = self.dot(rhs) / rhs_norm_sq; + rhs * scale + } + } + + fn zero() -> Self; +} +