diff --git a/src/lib.rs b/src/lib.rs index 15535c0..71b2fe9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,6 +4,7 @@ pub mod mathx; pub mod mesh_loader; pub mod mesh_tracer; pub mod riemann; +pub mod shape; pub mod tube; pub mod types; pub mod utils; diff --git a/src/shape/cylinder.rs b/src/shape/cylinder.rs new file mode 100644 index 0000000..baf21a9 --- /dev/null +++ b/src/shape/cylinder.rs @@ -0,0 +1,154 @@ +use glam::{Vec3, Vec3Swizzles as _}; + +use crate::types::Ray; + +/// Цилиндр с центром в начале координат и осью вдоль оси Y. +pub struct YCylinder { + pub half_length: f32, + pub radius: f32, +} + +impl YCylinder { + /// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии YCylinder). + fn flip_ray(ray: Ray) -> Ray { + Ray { + pos: ray.pos * ray.dir.signum(), + dir: ray.dir.abs(), + } + } + + pub fn is_inside(&self, pt: Vec3) -> bool { + let r = f32::hypot(pt.x, pt.z); + pt.y.abs() < self.half_length && r < self.radius + } + + pub fn trace_into(&self, ray: Ray) -> Option { + let ray = Self::flip_ray(ray); + + // 1. ray.pos.y + t * ray.dir.y = −half_length + let t_cap_in = (-self.half_length - ray.pos.y) / ray.dir.y; + let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y; + + // 2. (ray.pos.x + t * ray.dir.x)² + (ray.pos.z + t * ray.dir.z)² = radius² + let pos = ray.pos.xz(); + let dir = ray.dir.xz(); + if dir.length_squared() < 1e-3 { + if pos.length_squared() >= self.radius.powi(2) { + return None; + } + return Some(t_cap_in).filter(|&t| t > 0.); + } + let (t_side_in, t_side_out) = solve_quadratic( + dir.length_squared(), + pos.dot(dir), + pos.length_squared() - self.radius.powi(2), + )?; + let t = f32::max(t_cap_in, t_side_in); + if t < 0. { + return None; + } + if t >= t_cap_out || t >= t_side_out { + return None; + } + Some(t) + } + + pub fn trace_out_of(&self, ray: Ray) -> Option { + let ray = Self::flip_ray(ray); + let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y; + let pos = ray.pos.xz(); + let dir = ray.dir.xz(); + if dir.length_squared() < 1e-3 { + return Some(t_cap_out); + } + let (_t_side_in, t_side_out) = solve_quadratic( + dir.length_squared(), + pos.dot(dir), + pos.length_squared() - self.radius.powi(2), + ) + .expect("the ray starts inside and is not along the axis so *has* to cross the side"); + Some(t_side_out) + } +} + +fn solve_quadratic(a: f32, half_b: f32, c: f32) -> Option<(f32, f32)> { + let base = -half_b / a; + let d = base * base - c / a; + if d < 0. { + None + } else { + let δ = d.sqrt(); + Some((base - δ, base + δ)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::types::ray; + use approx::assert_abs_diff_eq; + use glam::vec3; + + #[test] + fn test_cylinder() { + assert_eq!( + YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))), + ray(vec3(2., 3., 2.), vec3(4., 5., 4.)), + ); + assert_eq!( + YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))), + ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)), + ); + assert_eq!( + YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))), + ray(vec3(2., -3., 2.), vec3(4., 5., 4.)), + ); + assert_eq!( + YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))), + ray(vec3(2., 3., 2.), vec3(4., 0., 4.)), + ); + + let r = YCylinder { + half_length: 3., + radius: 2., + }; + + assert_eq!(r.trace_into(ray(vec3(3., 4., 3.), vec3(0., -1., 0.))), None); + assert_eq!(r.trace_into(ray(vec3(1., 4., 1.), vec3(0., -1., 0.))), Some(1.)); + assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None); + assert_abs_diff_eq!( + r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))).unwrap(), + 1.5857864 + ); + assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None); + assert_abs_diff_eq!( + r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))).unwrap(), + 0.7928932 + ); + assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None); + + assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None); + assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None); + assert_eq!( + r.trace_into(ray(vec3(1.4142135, 3., 1.4142135), vec3(-1., -1., -1.))), + Some(0.) + ); + assert_eq!( + r.trace_into(ray(vec3(1.4142135, -3., 1.4142135), vec3(-1., 1., -1.))), + Some(0.) + ); + + assert_abs_diff_eq!( + r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))).unwrap(), + 1.4142135 + ); + assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.)); + assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.)); + assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.)); + assert_abs_diff_eq!( + r.trace_out_of(ray(vec3(1.4142135, 3., 1.4142135), vec3(1., 1., 1.))) + .unwrap(), + 0. + ); + } +} diff --git a/src/shape/mod.rs b/src/shape/mod.rs new file mode 100644 index 0000000..2803b51 --- /dev/null +++ b/src/shape/mod.rs @@ -0,0 +1,5 @@ +pub mod cylinder; +pub mod rect; + +pub use cylinder::YCylinder; +pub use rect::Rect; diff --git a/src/shape/rect.rs b/src/shape/rect.rs new file mode 100644 index 0000000..f71284c --- /dev/null +++ b/src/shape/rect.rs @@ -0,0 +1,90 @@ +use glam::Vec3; + +use crate::types::Ray; + +pub struct Rect { + pub size: Vec3, +} + +impl Rect { + /// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect). + fn flip_ray(ray: Ray) -> Ray { + Ray { + pos: ray.pos * ray.dir.signum(), + dir: ray.dir.abs(), + } + } + + pub fn is_inside(&self, pt: Vec3) -> bool { + pt.abs().cmplt(self.size).all() + } + + pub fn trace_into(&self, ray: Ray) -> Option { + let ray = Self::flip_ray(ray); + // ray.pos.x + t * ray.dir.x = −size.x + let ts = (-self.size - ray.pos) / ray.dir; + let t = ts.max_element(); + let pt = ray.pos + t * ray.dir; + if t < 0.0 { + return None; + } + if pt.cmpgt(self.size).any() { + return None; + } + Some(t) + } + + pub fn trace_out_of(&self, ray: Ray) -> Option { + let ray = Self::flip_ray(ray); + // ray.pos.x + t * ray.dir.x = +size.x + let ts = (self.size - ray.pos) / ray.dir; + let t = ts.min_element(); + Some(t) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::types::ray; + use glam::vec3; + + #[test] + fn test_rect() { + assert_eq!( + Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))), + ray(vec3(2., 3., 2.), vec3(4., 5., 4.)), + ); + assert_eq!( + Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))), + ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)), + ); + assert_eq!( + Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))), + ray(vec3(2., -3., 2.), vec3(4., 5., 4.)), + ); + assert_eq!( + Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))), + ray(vec3(2., 3., 2.), vec3(4., 0., 4.)), + ); + + let r = Rect { size: vec3(2., 3., 2.) }; + + assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None); + assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))), Some(1.)); + assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None); + assert_eq!(r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))), Some(0.5)); + assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None); + + assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None); + assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None); + assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(-1., -1., -1.))), Some(0.)); + assert_eq!(r.trace_into(ray(vec3(2., -3., 2.), vec3(-1., 1., -1.))), Some(0.)); + + assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))), Some(2.)); + assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.)); + assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.)); + assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.)); + assert_eq!(r.trace_out_of(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), Some(0.)); + } +} diff --git a/src/tube/coords.rs b/src/tube/coords.rs index 657420b..6516b86 100644 --- a/src/tube/coords.rs +++ b/src/tube/coords.rs @@ -1,9 +1,10 @@ use glam::{vec3, Mat3, Vec3}; use crate::riemann::Metric; +use crate::shape::YCylinder; use crate::types::{Location, Ray}; -use super::{Tube, YCylinder}; +use super::Tube; pub trait FlatCoordinateSystem { fn flat_to_global(&self, v: T) -> T; diff --git a/src/tube/mod.rs b/src/tube/mod.rs index 7bd5a70..e47d300 100644 --- a/src/tube/mod.rs +++ b/src/tube/mod.rs @@ -1,4 +1,4 @@ -use glam::{bool, f32, Mat3, Vec3, Vec3Swizzles}; +use glam::{f32, Mat3, Vec3}; use crate::ifaces::{DebugTraceable, RayPath, Traceable}; use coords::{FlatCoordinateSystem, InnerCS, OuterCS}; @@ -226,234 +226,3 @@ impl DebugTraceable for Space { panic!("tracing didn't terminate"); } } - -struct Rect { - pub size: Vec3, -} - -impl Rect { - /// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect). - fn flip_ray(ray: Ray) -> Ray { - Ray { - pos: ray.pos * ray.dir.signum(), - dir: ray.dir.abs(), - } - } - - fn is_inside(&self, pt: Vec3) -> bool { - pt.abs().cmplt(self.size).all() - } - - fn trace_into(&self, ray: Ray) -> Option { - let ray = Self::flip_ray(ray); - // ray.pos.x + t * ray.dir.x = −size.x - let ts = (-self.size - ray.pos) / ray.dir; - let t = ts.max_element(); - let pt = ray.pos + t * ray.dir; - if t < 0.0 { - return None; - } - if pt.cmpgt(self.size).any() { - return None; - } - Some(t) - } - - fn trace_out_of(&self, ray: Ray) -> Option { - let ray = Self::flip_ray(ray); - // ray.pos.x + t * ray.dir.x = +size.x - let ts = (self.size - ray.pos) / ray.dir; - let t = ts.min_element(); - Some(t) - } -} - -fn solve_quadratic(a: f32, half_b: f32, c: f32) -> Option<(f32, f32)> { - let base = -half_b / a; - let d = base * base - c / a; - if d < 0. { - None - } else { - let δ = d.sqrt(); - Some((base - δ, base + δ)) - } -} - -/// Цилиндр с центром в начале координат и осью вдоль оси Y. -struct YCylinder { - pub half_length: f32, - pub radius: f32, -} - -impl YCylinder { - /// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии YCylinder). - fn flip_ray(ray: Ray) -> Ray { - Ray { - pos: ray.pos * ray.dir.signum(), - dir: ray.dir.abs(), - } - } - - fn is_inside(&self, pt: Vec3) -> bool { - let r = f32::hypot(pt.x, pt.z); - pt.y.abs() < self.half_length && r < self.radius - } - - fn trace_into(&self, ray: Ray) -> Option { - let ray = Self::flip_ray(ray); - - // 1. ray.pos.y + t * ray.dir.y = −half_length - let t_cap_in = (-self.half_length - ray.pos.y) / ray.dir.y; - let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y; - - // 2. (ray.pos.x + t * ray.dir.x)² + (ray.pos.z + t * ray.dir.z)² = radius² - let pos = ray.pos.xz(); - let dir = ray.dir.xz(); - if dir.length_squared() < 1e-3 { - if pos.length_squared() >= self.radius.powi(2) { - return None; - } - return Some(t_cap_in).filter(|&t| t > 0.); - } - let (t_side_in, t_side_out) = solve_quadratic( - dir.length_squared(), - pos.dot(dir), - pos.length_squared() - self.radius.powi(2), - )?; - let t = f32::max(t_cap_in, t_side_in); - if t < 0. { - return None; - } - if t >= t_cap_out || t >= t_side_out { - return None; - } - Some(t) - } - - fn trace_out_of(&self, ray: Ray) -> Option { - let ray = Self::flip_ray(ray); - let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y; - let pos = ray.pos.xz(); - let dir = ray.dir.xz(); - if dir.length_squared() < 1e-3 { - return Some(t_cap_out); - } - let (_t_side_in, t_side_out) = solve_quadratic( - dir.length_squared(), - pos.dot(dir), - pos.length_squared() - self.radius.powi(2), - ) - .expect("the ray starts inside and is not along the axis so *has* to cross the side"); - Some(t_side_out) - } -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::types::ray; - use approx::assert_abs_diff_eq; - use glam::vec3; - - #[test] - fn test_rect() { - assert_eq!( - Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))), - ray(vec3(2., 3., 2.), vec3(4., 5., 4.)), - ); - assert_eq!( - Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))), - ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)), - ); - assert_eq!( - Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))), - ray(vec3(2., -3., 2.), vec3(4., 5., 4.)), - ); - assert_eq!( - Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))), - ray(vec3(2., 3., 2.), vec3(4., 0., 4.)), - ); - - let r = Rect { size: vec3(2., 3., 2.) }; - - assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None); - assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))), Some(1.)); - assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None); - assert_eq!(r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))), Some(0.5)); - assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None); - - assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None); - assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None); - assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(-1., -1., -1.))), Some(0.)); - assert_eq!(r.trace_into(ray(vec3(2., -3., 2.), vec3(-1., 1., -1.))), Some(0.)); - - assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))), Some(2.)); - assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.)); - assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.)); - assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.)); - assert_eq!(r.trace_out_of(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), Some(0.)); - } - - #[test] - fn test_cylinder() { - assert_eq!( - YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))), - ray(vec3(2., 3., 2.), vec3(4., 5., 4.)), - ); - assert_eq!( - YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))), - ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)), - ); - assert_eq!( - YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))), - ray(vec3(2., -3., 2.), vec3(4., 5., 4.)), - ); - assert_eq!( - YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))), - ray(vec3(2., 3., 2.), vec3(4., 0., 4.)), - ); - - let r = YCylinder { - half_length: 3., - radius: 2., - }; - - assert_eq!(r.trace_into(ray(vec3(3., 4., 3.), vec3(0., -1., 0.))), None); - assert_eq!(r.trace_into(ray(vec3(1., 4., 1.), vec3(0., -1., 0.))), Some(1.)); - assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None); - assert_abs_diff_eq!( - r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))).unwrap(), - 1.5857864 - ); - assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None); - assert_abs_diff_eq!( - r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))).unwrap(), - 0.7928932 - ); - assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None); - - assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None); - assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None); - assert_eq!( - r.trace_into(ray(vec3(1.4142135, 3., 1.4142135), vec3(-1., -1., -1.))), - Some(0.) - ); - assert_eq!( - r.trace_into(ray(vec3(1.4142135, -3., 1.4142135), vec3(-1., 1., -1.))), - Some(0.) - ); - - assert_abs_diff_eq!( - r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))).unwrap(), - 1.4142135 - ); - assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.)); - assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.)); - assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.)); - assert_abs_diff_eq!( - r.trace_out_of(ray(vec3(1.4142135, 3., 1.4142135), vec3(1., 1., 1.))) - .unwrap(), - 0. - ); - } -}