diff --git a/src/core/geometry/bounds.rs b/src/core/geometry/bounds.rs new file mode 100644 index 0000000..28f3b1a --- /dev/null +++ b/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 { + 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 centroid(&self) -> Point { + 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) -> 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 corner(&self, corner_index: usize) -> Point { + 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) -> 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 + Sqrt, +{ + 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 { + #[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) + } +} diff --git a/src/core/geometry/cone.rs b/src/core/geometry/cone.rs new file mode 100644 index 0000000..3fb84f3 --- /dev/null +++ b/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()) + } +} diff --git a/src/core/geometry/mod.rs b/src/core/geometry/mod.rs new file mode 100644 index 0000000..71a8e84 --- /dev/null +++ b/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(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); + 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 } +} diff --git a/src/core/geometry/primitives.rs b/src/core/geometry/primitives.rs new file mode 100644 index 0000000..e178b0b --- /dev/null +++ b/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(pub [T; N]); +// N-dimensional location +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +pub struct Point(pub [T; N]); +// N-dimensional surface normal +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +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)) + } + + #[inline] + pub fn average(&self) -> f32 { + let sum: f32 = self.0.iter().sum(); + sum / (N as f32) + } + } + + impl $Struct + 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 $Struct + where + T: Copy, + { + #[inline] + pub fn fill(value: T) -> Self { + Self([value; N]) + } + + #[inline] + pub fn cast(&self) -> $Struct + where + U: 'static + Copy, + T: 'static + Copy + AsPrimitive, + { + $Struct(self.0.map(|c| c.as_())) + } + } + + 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) => { + impl VectorLike for $Struct + where + T: Copy + + Zero + + Add + + Mul + + Sub + + Div + + 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 $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_abs!(Point); +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 + 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 = 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 Normal3 +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)) + } + + 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 Normal3 +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)) + } + + 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 Hash for Vector { + fn hash(&self, state: &mut H) { + for item in self.0.iter() { + item.to_bits().hash(state); + } + } +} + +impl Hash for Point { + fn hash(&self, state: &mut H) { + for item in self.0.iter() { + item.to_bits().hash(state); + } + } +} + +impl Hash for Normal { + fn hash(&self, state: &mut H) { + for item in self.0.iter() { + item.to_bits().hash(state); + } + } +} + +// INTERVAL STUFF +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 From> for Vector { + fn from(pi: Vector) -> Self { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = pi[i].midpoint(); + } + Vector(arr) + } +} + +impl From> for Vector { + fn from(v: Vector) -> Self { + let mut arr = [Interval::default(); N]; + for i in 0..N { + arr[i] = Interval::new(v[i]); + } + Self(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 + Sqrt, +{ + pub fn face_forward(self, v: Vector3) -> Self { + if Vector3::::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 for OctahedralVector { + fn from(v: Vector3f) -> Self { + Self::new(v) + } +} + +impl From 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()) + } +} diff --git a/src/core/geometry/ray.rs b/src/core/geometry/ray.rs new file mode 100644 index 0000000..d6780c8 --- /dev/null +++ b/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>, + 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 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, +} diff --git a/src/core/geometry/traits.rs b/src/core/geometry/traits.rs new file mode 100644 index 0000000..fc48814 --- /dev/null +++ b/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: + Sized + Copy + Index + IndexMut +{ + 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 + + Sub + + Div + + Mul +{ + type Scalar: Copy + Zero + Add + Mul + 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: Sized + Copy { + fn lerp(t: Factor, a: Self, b: Self) -> Self; +} + +impl Lerp for T +where + T: Copy + Sub + Add, + Diff: Mul, + F: Copy, +{ + #[inline(always)] + fn lerp(t: F, a: Self, b: Self) -> Self { + a + (b - a) * t + } +} diff --git a/src/core/scattering.rs b/src/core/scattering.rs new file mode 100644 index 0000000..db49405 --- /dev/null +++ b/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 { + let cos_corr = clamp_t(cos_theta_i, 0., 1.); + let sin2_theta_i = 1. - square(cos_corr); + let sin2_theta_t: Complex = sin2_theta_i / square(eta); + let cos2_theta_t: Complex = (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 +} diff --git a/src/gpu/mod.rs b/src/gpu/mod.rs new file mode 100644 index 0000000..e69de29 diff --git a/src/spectra/color.rs b/src/spectra/color.rs new file mode 100644 index 0000000..5ee7c2a --- /dev/null +++ b/src/spectra/color.rs @@ -0,0 +1,1001 @@ +use std::any::TypeId; +use std::fmt; +use std::ops::{ + Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign, +}; + +use crate::core::geometry::Point2f; +use crate::core::pbrt::Float; +use crate::spectra::Spectrum; +use crate::utils::math::{SquareMatrix, evaluate_polynomial, lerp}; +use once_cell::sync::Lazy; + +use enum_dispatch::enum_dispatch; + +pub trait Triplet { + fn from_triplet(c1: Float, c2: Float, c3: Float) -> Self; +} + +#[derive(Debug, Clone)] +pub struct XYZ { + pub x: Float, + pub y: Float, + pub z: Float, +} + +impl Triplet for XYZ { + fn from_triplet(c1: Float, c2: Float, c3: Float) -> Self { + XYZ::new(c1, c2, c3) + } +} + +impl<'a> IntoIterator for &'a XYZ { + type Item = &'a Float; + type IntoIter = std::array::IntoIter<&'a Float, 3>; + + fn into_iter(self) -> Self::IntoIter { + [&self.x, &self.y, &self.z].into_iter() + } +} + +impl XYZ { + pub fn new(x: Float, y: Float, z: Float) -> Self { + Self { x, y, z } + } + + pub fn x(&self) -> Float { + self.x + } + + pub fn y(&self) -> Float { + self.y + } + + pub fn z(&self) -> Float { + self.z + } + + pub fn average(&self) -> Float { + (self.x + self.y + self.z) / 3.0 + } + + pub fn xy(&self) -> Point2f { + let sum = self.x + self.y + self.z; + Point2f::new(self.x / sum, self.y / sum) + } + + pub fn from_xyy(xy: Point2f, y: Option) -> Self { + let scale: Float; + if let Some(y_val) = y { + scale = y_val; + } else { + scale = 1.; + } + if xy.y() == 0.0 { + return Self { + x: 0.0, + y: 0.0, + z: 0.0, + }; + } + Self::new( + xy.x() * scale / xy.y(), + scale, + (1.0 - xy.x() - xy.y()) * scale / xy.y(), + ) + } +} + +impl Index for XYZ { + type Output = Float; + fn index(&self, index: usize) -> &Self::Output { + debug_assert!(index < 3); + match index { + 0 => &self.x, + 1 => &self.y, + _ => &self.z, + } + } +} + +impl IndexMut for XYZ { + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + debug_assert!(index < 3); + match index { + 0 => &mut self.x, + 1 => &mut self.y, + _ => &mut self.z, + } + } +} + +impl Neg for XYZ { + type Output = Self; + fn neg(self) -> Self { + Self { + x: -self.x, + y: -self.y, + z: -self.z, + } + } +} + +impl Add for XYZ { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Self { + x: self.x + rhs.x, + y: self.y + rhs.y, + z: self.z + rhs.z, + } + } +} + +impl AddAssign for XYZ { + fn add_assign(&mut self, rhs: Self) { + self.x += rhs.x; + self.y += rhs.y; + self.z += rhs.z; + } +} + +impl Sub for XYZ { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + Self { + x: self.x - rhs.x, + y: self.y - rhs.y, + z: self.z - rhs.z, + } + } +} + +impl SubAssign for XYZ { + fn sub_assign(&mut self, rhs: Self) { + self.x -= rhs.x; + self.y -= rhs.y; + self.z -= rhs.z; + } +} + +impl Sub for Float { + type Output = XYZ; + fn sub(self, rhs: XYZ) -> XYZ { + XYZ::new(self - rhs.x, self - rhs.y, self - rhs.z) + } +} + +impl Mul for XYZ { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Self { + x: self.x * rhs.x, + y: self.y * rhs.y, + z: self.z * rhs.z, + } + } +} + +impl MulAssign for XYZ { + fn mul_assign(&mut self, rhs: Self) { + self.x *= rhs.x; + self.y *= rhs.y; + self.z *= rhs.z; + } +} + +impl Mul for XYZ { + type Output = Self; + fn mul(self, rhs: Float) -> Self { + Self { + x: self.x * rhs, + y: self.y * rhs, + z: self.z * rhs, + } + } +} + +impl MulAssign for XYZ { + fn mul_assign(&mut self, rhs: Float) { + self.x *= rhs; + self.y *= rhs; + self.z *= rhs; + } +} + +impl Mul for Float { + type Output = XYZ; + fn mul(self, rhs: XYZ) -> XYZ { + rhs * self + } +} + +impl Div for XYZ { + type Output = Self; + fn div(self, rhs: Self) -> Self { + Self { + x: self.x / rhs.x, + y: self.y / rhs.y, + z: self.z / rhs.z, + } + } +} + +impl DivAssign for XYZ { + fn div_assign(&mut self, rhs: Self) { + self.x /= rhs.x; + self.y /= rhs.y; + self.z /= rhs.z; + } +} + +impl Div for XYZ { + type Output = Self; + fn div(self, rhs: Float) -> Self { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + self * inv + } +} + +impl DivAssign for XYZ { + fn div_assign(&mut self, rhs: Float) { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + *self *= inv; + } +} + +impl fmt::Display for XYZ { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "[ {}, {}, {} ]", self.x, self.y, self.z) + } +} + +#[derive(Debug, Default, Copy, Clone)] +pub struct RGB { + pub r: Float, + pub g: Float, + pub b: Float, +} + +impl Triplet for RGB { + fn from_triplet(c1: Float, c2: Float, c3: Float) -> Self { + RGB::new(c1, c2, c3) + } +} + +impl<'a> IntoIterator for &'a RGB { + type Item = &'a Float; + type IntoIter = std::array::IntoIter<&'a Float, 3>; + + fn into_iter(self) -> Self::IntoIter { + [&self.r, &self.g, &self.b].into_iter() + } +} + +impl RGB { + pub fn new(r: Float, g: Float, b: Float) -> Self { + Self { r, g, b } + } + + pub fn average(&self) -> Float { + (self.r + self.g + self.b) / 3.0 + } + + pub fn max(&self) -> Float { + self.r.max(self.g).max(self.b) + } + + pub fn clamp_zero(rgb: Self) -> Self { + RGB::new(rgb.r.max(0.), rgb.b.max(0.), rgb.b.max(0.)) + } +} + +impl Index for RGB { + type Output = Float; + fn index(&self, index: usize) -> &Self::Output { + debug_assert!(index < 3); + match index { + 0 => &self.r, + 1 => &self.g, + _ => &self.b, + } + } +} + +impl IndexMut for RGB { + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + debug_assert!(index < 3); + match index { + 0 => &mut self.r, + 1 => &mut self.g, + _ => &mut self.b, + } + } +} + +impl Neg for RGB { + type Output = Self; + fn neg(self) -> Self { + Self { + r: -self.r, + g: -self.g, + b: -self.b, + } + } +} + +impl Add for RGB { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Self { + r: self.r + rhs.r, + g: self.g + rhs.g, + b: self.b + rhs.b, + } + } +} + +impl AddAssign for RGB { + fn add_assign(&mut self, rhs: Self) { + self.r += rhs.r; + self.g += rhs.g; + self.b += rhs.b; + } +} + +impl Sub for RGB { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + Self { + r: self.r - rhs.r, + g: self.g - rhs.g, + b: self.b - rhs.b, + } + } +} + +impl SubAssign for RGB { + fn sub_assign(&mut self, rhs: Self) { + self.r -= rhs.r; + self.g -= rhs.g; + self.b -= rhs.b; + } +} + +impl Sub for Float { + type Output = RGB; + fn sub(self, rhs: RGB) -> RGB { + RGB::new(self - rhs.r, self - rhs.g, self - rhs.b) + } +} + +impl Mul for RGB { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Self { + r: self.r * rhs.r, + g: self.g * rhs.g, + b: self.b * rhs.b, + } + } +} + +impl MulAssign for RGB { + fn mul_assign(&mut self, rhs: Self) { + self.r *= rhs.r; + self.g *= rhs.g; + self.b *= rhs.b; + } +} + +// Scalar multiplication (ryz * float) +impl Mul for RGB { + type Output = Self; + fn mul(self, rhs: Float) -> Self { + Self { + r: self.r * rhs, + g: self.g * rhs, + b: self.b * rhs, + } + } +} + +impl MulAssign for RGB { + fn mul_assign(&mut self, rhs: Float) { + self.r *= rhs; + self.g *= rhs; + self.b *= rhs; + } +} + +impl Mul for Float { + type Output = RGB; + fn mul(self, rhs: RGB) -> RGB { + rhs * self + } +} + +impl Div for RGB { + type Output = Self; + fn div(self, rhs: Self) -> Self { + Self { + r: self.r / rhs.r, + g: self.g / rhs.g, + b: self.b / rhs.b, + } + } +} + +impl DivAssign for RGB { + fn div_assign(&mut self, rhs: Self) { + self.r /= rhs.r; + self.g /= rhs.g; + self.b /= rhs.b; + } +} + +impl Div for RGB { + type Output = Self; + fn div(self, rhs: Float) -> Self { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + self * inv + } +} + +impl DivAssign for RGB { + fn div_assign(&mut self, rhs: Float) { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + *self *= inv; + } +} + +impl fmt::Display for RGB { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "[ {}, {}, {} ]", self.r, self.g, self.b) + } +} + +impl Mul for SquareMatrix { + type Output = RGB; + + fn mul(self, v: XYZ) -> RGB { + let r = self[0][0] * v.x + self[0][1] * v.y + self[0][2] * v.z; + let g = self[1][0] * v.x + self[1][1] * v.y + self[1][2] * v.z; + let b = self[2][0] * v.x + self[2][1] * v.y + self[2][2] * v.z; + RGB::new(r, g, b) + } +} + +impl Mul for SquareMatrix { + type Output = XYZ; + fn mul(self, v: RGB) -> XYZ { + let x = self[0][0] * v.r + self[0][1] * v.g + self[0][2] * v.b; + let y = self[1][0] * v.r + self[1][1] * v.g + self[1][2] * v.b; + let z = self[2][0] * v.r + self[2][1] * v.g + self[2][2] * v.b; + XYZ::new(x, y, z) + } +} + +pub trait MatrixMulColor { + fn mul_rgb(&self, v: RGB) -> RGB; + fn mul_xyz(&self, v: XYZ) -> XYZ; +} + +impl MatrixMulColor for SquareMatrix { + fn mul_rgb(&self, v: RGB) -> RGB { + let m = self; + RGB::new( + m[0][0] * v.r + m[0][1] * v.g + m[0][2] * v.b, + m[1][0] * v.r + m[1][1] * v.g + m[1][2] * v.b, + m[2][0] * v.r + m[2][1] * v.g + m[2][2] * v.b, + ) + } + + fn mul_xyz(&self, v: XYZ) -> XYZ { + let m = self; + XYZ::new( + m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z, + m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z, + m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z, + ) + } +} + +pub const RES: usize = 64; +pub type CoefficientArray = [[[[[Float; 3]; RES]; RES]; RES]; 3]; + +#[derive(Debug)] +pub struct RGBToSpectrumTable { + z_nodes: &'static [f32], + coeffs: &'static CoefficientArray, +} + +impl RGBToSpectrumTable { + pub fn srgb() -> Self { + // use crate::core::constants::{RGB_TO_SPECTRUM_Z_NODES, RGB_TO_SPECTRUM_COEFFS}; + // Self::new(&RGB_TO_SPECTRUM_Z_NODES, &RGB_TO_SPECTRUM_COEFFS) + todo!("Link the static constant arrays for sRGB coefficients here") + } +} + +#[derive(Debug, Default, Copy, Clone)] +pub struct RGBSigmoidPolynomial { + c0: Float, + c1: Float, + c2: Float, +} + +impl RGBSigmoidPolynomial { + pub fn new(c0: Float, c1: Float, c2: Float) -> Self { + Self { c0, c1, c2 } + } + + pub fn evaluate(&self, lambda: Float) -> Float { + let eval = match evaluate_polynomial(lambda, &[self.c0, self.c1, self.c2]) { + Some(value) => value, + None => { + panic!("evaluate_polynomial returned None with non-empty coefficients") + } + }; + + Self::s(eval) + } + + pub fn max_value(&self) -> Float { + let lambda = -self.c1 / (2.0 * self.c0); + let result = self.evaluate(360.0).max(self.evaluate(830.0)); + if (360.0..830.0).contains(&lambda) { + return result.max(self.evaluate(lambda)); + } + result + } + + fn s(x: Float) -> Float { + if x.is_infinite() { + if x > 0.0 { return 1.0 } else { return 0.0 } + } + 0.5 + x / (2.0 * (1.0 + (x * x)).sqrt()) + } +} + +impl RGBToSpectrumTable { + pub fn new(z_nodes: &'static [f32], coeffs: &'static CoefficientArray) -> Self { + Self { z_nodes, coeffs } + } + + pub fn to_polynomial(&self, rgb: RGB) -> RGBSigmoidPolynomial { + if rgb[0] == rgb[1] && rgb[1] == rgb[2] { + return RGBSigmoidPolynomial::new( + 0.0, + 0.0, + (rgb[0] - 0.5) / (rgb[0] * (1.0 - rgb[0])).sqrt(), + ); + } + let maxc; + if rgb[0] > rgb[1] { + if rgb[0] > rgb[2] { + maxc = 0; + } else { + maxc = 2; + } + } else if rgb[1] > rgb[2] { + maxc = 1; + } else { + maxc = 2; + } + + let z = rgb[maxc]; + let x = rgb[(maxc + 1) % 3] * (RES - 1) as Float / z; + let y = rgb[(maxc + 2) % 3] * (RES - 1) as Float / z; + + let xi = x.min(RES as Float - 2.0); + let yi = y.min(RES as Float - 2.0); + let zi = crate::core::pbrt::find_interval(RES, |i: usize| self.z_nodes[i] < z); + let dx = (x - xi) as usize; + let dy = (y - yi) as usize; + let dz = (z - self.z_nodes[zi]) / (self.z_nodes[zi + 1] - self.z_nodes[zi]); + let mut c = [0.0; 3]; + #[allow(clippy::needless_range_loop)] + for i in 0..3 { + let co = |dx: usize, dy: usize, dz: usize| { + self.coeffs[maxc][zi as usize + dz][yi as usize + dy][xi as usize + dx][i] + }; + c[i] = lerp( + dz, + lerp( + dy as Float, + lerp(dx as Float, co(0, 0, 0) as Float, co(1, 0, 0)) as Float, + lerp(dx as Float, co(0, 1, 0) as Float, co(1, 1, 0) as Float), + ), + lerp( + dy as Float, + lerp(dx as Float, co(0, 0, 1) as Float, co(1, 0, 1)) as Float, + lerp(dx as Float, co(0, 1, 1) as Float, co(1, 1, 1) as Float), + ), + ); + } + RGBSigmoidPolynomial { + c0: c[0], + c1: c[1], + c2: c[2], + } + } +} + +const LMS_FROM_XYZ: SquareMatrix = SquareMatrix::new([ + [0.8951, 0.2664, -0.1614], + [-0.7502, 1.7135, 0.0367], + [0.0389, -0.0685, 1.0296], +]); + +const XYZ_FROM_LMS: SquareMatrix = SquareMatrix::new([ + [0.986993, -0.147054, 0.159963], + [0.432305, 0.51836, 0.0492912], + [-0.00852866, 0.0400428, 0.968487], +]); + +pub fn white_balance(src_white: Point2f, target_white: Point2f) -> SquareMatrix { + // Find LMS coefficients for source and target white + let src_xyz = XYZ::from_xyy(src_white, None); + let dst_xyz = XYZ::from_xyy(target_white, None); + let src_lms = LMS_FROM_XYZ * src_xyz; + let dst_lms = LMS_FROM_XYZ * dst_xyz; + + // Return white balancing matrix for source and target white + let lms_correct = SquareMatrix::::diag(&[ + dst_lms[0] / src_lms[0], + dst_lms[1] / src_lms[1], + dst_lms[2] / src_lms[2], + ]); + XYZ_FROM_LMS * lms_correct * LMS_FROM_XYZ +} + +#[enum_dispatch] +pub trait ColorEncodingTrait: 'static + Send + Sync + fmt::Debug + fmt::Display { + fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]); + fn to_linear_slice(&self, vin: &[u8], vout: &mut [Float]); + fn to_float_linear(&self, v: Float) -> Float; + fn type_id(&self) -> TypeId { + TypeId::of::() + } +} + +#[enum_dispatch(ColorEncodingTrait)] +#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)] +pub enum ColorEncoding { + Linear(LinearEncoding), + SRGB(SRGBEncoding), +} + +impl fmt::Display for ColorEncoding { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Encoding") + } +} + +#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)] +pub struct LinearEncoding; +impl ColorEncodingTrait for LinearEncoding { + fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]) { + for (i, &v) in vin.iter().enumerate() { + vout[i] = (v.clamp(0.0, 1.0) * 255.0 + 0.5) as u8; + } + } + fn to_linear_slice(&self, vin: &[u8], vout: &mut [Float]) { + for (i, &v) in vin.iter().enumerate() { + vout[i] = v as Float / 255.0; + } + } + fn to_float_linear(&self, v: Float) -> Float { + v + } +} + +impl fmt::Display for LinearEncoding { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Linear") + } +} + +#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)] +pub struct SRGBEncoding; +impl ColorEncodingTrait for SRGBEncoding { + fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]) { + for (i, &v_linear) in vin.iter().enumerate() { + let v = v_linear.clamp(0.0, 1.0); + let v_encoded = if v <= 0.0031308 { + v * 12.92 + } else { + 1.055 * v.powf(1.0 / 2.4) - 0.055 + }; + vout[i] = (v_encoded * 255.0 + 0.5) as u8; + } + } + + fn to_linear_slice(&self, vin: &[u8], vout: &mut [Float]) { + for (i, &v) in vin.iter().enumerate() { + vout[i] = SRGB_TO_LINEAR_LUT[v as usize]; + } + } + + fn to_float_linear(&self, v: Float) -> Float { + let v = v.clamp(0.0, 1.0); + if v <= 0.04045 { + v / 12.92 + } else { + ((v + 0.055) / 1.055).powf(2.4) + } + } +} + +impl fmt::Display for SRGBEncoding { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "sRGB") + } +} + +pub const LINEAR: ColorEncoding = ColorEncoding::Linear(LinearEncoding); +pub const SRGB: ColorEncoding = ColorEncoding::SRGB(SRGBEncoding); + +const SRGB_TO_LINEAR_LUT: [Float; 256] = [ + 0.0000000000, + 0.0003035270, + 0.0006070540, + 0.0009105810, + 0.0012141080, + 0.0015176350, + 0.0018211619, + 0.0021246888, + 0.0024282159, + 0.0027317430, + 0.0030352699, + 0.0033465356, + 0.0036765069, + 0.0040247170, + 0.0043914421, + 0.0047769533, + 0.0051815170, + 0.0056053917, + 0.0060488326, + 0.0065120910, + 0.0069954102, + 0.0074990317, + 0.0080231922, + 0.0085681248, + 0.0091340570, + 0.0097212177, + 0.0103298230, + 0.0109600937, + 0.0116122449, + 0.0122864870, + 0.0129830306, + 0.0137020806, + 0.0144438436, + 0.0152085144, + 0.0159962922, + 0.0168073755, + 0.0176419523, + 0.0185002182, + 0.0193823613, + 0.0202885624, + 0.0212190095, + 0.0221738834, + 0.0231533647, + 0.0241576303, + 0.0251868572, + 0.0262412224, + 0.0273208916, + 0.0284260381, + 0.0295568332, + 0.0307134409, + 0.0318960287, + 0.0331047624, + 0.0343398079, + 0.0356013142, + 0.0368894450, + 0.0382043645, + 0.0395462364, + 0.0409151986, + 0.0423114114, + 0.0437350273, + 0.0451862030, + 0.0466650836, + 0.0481718220, + 0.0497065634, + 0.0512694679, + 0.0528606549, + 0.0544802807, + 0.0561284944, + 0.0578054339, + 0.0595112406, + 0.0612460710, + 0.0630100295, + 0.0648032799, + 0.0666259527, + 0.0684781820, + 0.0703601092, + 0.0722718611, + 0.0742135793, + 0.0761853904, + 0.0781874284, + 0.0802198276, + 0.0822827145, + 0.0843762159, + 0.0865004659, + 0.0886556059, + 0.0908417329, + 0.0930589810, + 0.0953074843, + 0.0975873619, + 0.0998987406, + 0.1022417471, + 0.1046164930, + 0.1070231125, + 0.1094617173, + 0.1119324341, + 0.1144353822, + 0.1169706732, + 0.1195384338, + 0.1221387982, + 0.1247718409, + 0.1274376959, + 0.1301364899, + 0.1328683347, + 0.1356333494, + 0.1384316236, + 0.1412633061, + 0.1441284865, + 0.1470272839, + 0.1499598026, + 0.1529261619, + 0.1559264660, + 0.1589608639, + 0.1620294005, + 0.1651322246, + 0.1682693958, + 0.1714410931, + 0.1746473908, + 0.1778884083, + 0.1811642349, + 0.1844749898, + 0.1878207624, + 0.1912016720, + 0.1946178079, + 0.1980693042, + 0.2015562356, + 0.2050787061, + 0.2086368501, + 0.2122307271, + 0.2158605307, + 0.2195262313, + 0.2232279778, + 0.2269658893, + 0.2307400703, + 0.2345506549, + 0.2383976579, + 0.2422811985, + 0.2462013960, + 0.2501583695, + 0.2541521788, + 0.2581829131, + 0.2622507215, + 0.2663556635, + 0.2704978585, + 0.2746773660, + 0.2788943350, + 0.2831487954, + 0.2874408960, + 0.2917706966, + 0.2961383164, + 0.3005438447, + 0.3049873710, + 0.3094689548, + 0.3139887452, + 0.3185468316, + 0.3231432438, + 0.3277781308, + 0.3324515820, + 0.3371636569, + 0.3419144452, + 0.3467040956, + 0.3515326977, + 0.3564002514, + 0.3613068759, + 0.3662526906, + 0.3712377846, + 0.3762622178, + 0.3813261092, + 0.3864295185, + 0.3915725648, + 0.3967553079, + 0.4019778669, + 0.4072403014, + 0.4125427008, + 0.4178851545, + 0.4232677519, + 0.4286905527, + 0.4341537058, + 0.4396572411, + 0.4452012479, + 0.4507858455, + 0.4564110637, + 0.4620770514, + 0.4677838385, + 0.4735315442, + 0.4793202281, + 0.4851499796, + 0.4910208881, + 0.4969330430, + 0.5028865933, + 0.5088814497, + 0.5149177909, + 0.5209956765, + 0.5271152258, + 0.5332764983, + 0.5394796133, + 0.5457245708, + 0.5520114899, + 0.5583404899, + 0.5647116303, + 0.5711249113, + 0.5775805116, + 0.5840784907, + 0.5906189084, + 0.5972018838, + 0.6038274169, + 0.6104956269, + 0.6172066331, + 0.6239604354, + 0.6307572126, + 0.6375969648, + 0.6444797516, + 0.6514056921, + 0.6583748460, + 0.6653873324, + 0.6724432111, + 0.6795425415, + 0.6866854429, + 0.6938719153, + 0.7011020184, + 0.7083759308, + 0.7156936526, + 0.7230552435, + 0.7304608822, + 0.7379105687, + 0.7454043627, + 0.7529423237, + 0.7605246305, + 0.7681512833, + 0.7758223414, + 0.7835379243, + 0.7912980318, + 0.7991028428, + 0.8069523573, + 0.8148466945, + 0.8227858543, + 0.8307699561, + 0.8387991190, + 0.8468732834, + 0.8549926877, + 0.8631572723, + 0.8713672161, + 0.8796223402, + 0.8879231811, + 0.8962693810, + 0.9046613574, + 0.9130986929, + 0.9215820432, + 0.9301108718, + 0.9386858940, + 0.9473065734, + 0.9559735060, + 0.9646862745, + 0.9734454751, + 0.9822505713, + 0.9911022186, + 1.0000000000, +]; diff --git a/src/spectra/colorspace.rs b/src/spectra/colorspace.rs new file mode 100644 index 0000000..06970ff --- /dev/null +++ b/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, + pub xyz_from_rgb: SquareMatrix, + pub rgb_from_xyz: SquareMatrix, +} + +impl RGBColorSpace { + pub fn new( + r: Point2f, + g: Point2f, + b: Point2f, + illuminant: Spectrum, + rgb_to_spectrum_table: RGBToSpectrumTable, + ) -> Result> { + 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 { + if self == other { + return SquareMatrix::default(); + } + + self.rgb_from_xyz * other.xyz_from_rgb + } + + pub fn srgb() -> &'static Self { + static SRGB_SPACE: Lazy = 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) + } +}