241 lines
6.6 KiB
Rust
241 lines
6.6 KiB
Rust
use glam::{Vec3, Vec3Swizzles as _};
|
||
|
||
use crate::types::Ray;
|
||
|
||
pub struct Cylinder {
|
||
pub center: Vec3,
|
||
pub semiaxis: Vec3,
|
||
pub radius: f32,
|
||
}
|
||
|
||
impl Cylinder {
|
||
/// Split a vector into a component along the axis and one orthogonal to it.
|
||
fn split(&self, dir: Vec3) -> (f32, Vec3) {
|
||
let along = dir.dot(self.semiaxis) / self.semiaxis.length_squared();
|
||
(along, dir - along * self.semiaxis)
|
||
}
|
||
|
||
fn cap_inout(p: f32, d: f32) -> (f32, f32) {
|
||
((-d.signum() - p) / d.abs(), (d.signum() - p) / d.abs())
|
||
}
|
||
|
||
pub fn is_inside(&self, pt: Vec3) -> bool {
|
||
let (along, ortho) = self.split(pt - self.center);
|
||
along.abs() < 1. && ortho.length_squared() < self.radius.powi(2)
|
||
}
|
||
|
||
fn trace_inout(&self, ray: Ray) -> Option<(f32, f32)> {
|
||
let (pos_along, pos_ortho) = self.split(ray.pos - self.center);
|
||
let (dir_along, dir_ortho) = self.split(ray.dir);
|
||
let (t_cap_in, t_cap_out) = Self::cap_inout(pos_along, dir_along);
|
||
if dir_ortho.length_squared() < 1e-3 {
|
||
if pos_ortho.length_squared() >= self.radius.powi(2) {
|
||
return None;
|
||
}
|
||
return Some((t_cap_in, t_cap_out));
|
||
}
|
||
let (t_side_in, t_side_out) = solve_quadratic(
|
||
dir_ortho.length_squared(),
|
||
pos_ortho.dot(dir_ortho),
|
||
pos_ortho.length_squared() - self.radius.powi(2),
|
||
)?;
|
||
let t_in = f32::max(t_cap_in, t_side_in);
|
||
let t_out = f32::min(t_cap_out, t_side_out);
|
||
if t_out <= t_in {
|
||
return None;
|
||
}
|
||
Some((t_in, t_out))
|
||
}
|
||
|
||
pub fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||
let (t, _) = self.trace_inout(ray)?;
|
||
if t < 0. {
|
||
return None;
|
||
}
|
||
Some(t)
|
||
}
|
||
|
||
pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||
let (_, t) = self
|
||
.trace_inout(ray)
|
||
.expect("the ray starts inside so *has* to cross the boundary");
|
||
Some(t)
|
||
}
|
||
}
|
||
|
||
/// Цилиндр с центром в начале координат и осью вдоль оси 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<f32> {
|
||
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-6 * ray.dir.length_squared() {
|
||
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<f32> {
|
||
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;
|
||
use rand::{Rng, SeedableRng};
|
||
|
||
#[test]
|
||
fn test_split() {
|
||
let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
|
||
let cyl = Cylinder {
|
||
center: vec3(1., 2., 3.),
|
||
semiaxis: vec3(4., 5., 6.),
|
||
radius: 7.,
|
||
};
|
||
for _ in 0..100 {
|
||
let dir = vec3(rng.gen(), rng.gen(), rng.gen());
|
||
let (along, ortho) = cyl.split(dir);
|
||
assert_abs_diff_eq!(along * cyl.semiaxis + ortho, dir, epsilon = 1e-5);
|
||
assert_abs_diff_eq!(cyl.semiaxis.dot(ortho), 0., epsilon = 1e-5);
|
||
}
|
||
}
|
||
|
||
#[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_eq!(
|
||
YCylinder {
|
||
half_length: 300.,
|
||
radius: 50.
|
||
}
|
||
.trace_into(ray(vec3(-125., 375., 0.), vec3(3., -11., 0.) / 1024.)),
|
||
Some(25600.)
|
||
);
|
||
|
||
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.
|
||
);
|
||
}
|
||
}
|