#![expect(
clippy::many_single_char_names,
reason = "ray and frame algebra (o, p, q, a, c, d, r, s, t, z) mirrors the derivation in the plan; longer names would hide the formulas"
)]
#![expect(
clippy::pub_underscore_fields,
clippy::used_underscore_binding,
reason = "`_pad` is part of the specified 80-byte GPU layout and is read only to check it is zero"
)]
use crate::{geometry::plane::GpuFacetPlane, optics::raytracer::camera::Ray};
use glam::Vec3;
mod interval;
use interval::{ball_interval, frustum_interval, wedge_interval};
pub const MAX_TOOL_PRIMITIVES: usize = 128;
#[repr(u32)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum ToolKind {
Ball = 0,
Cylinder = 1,
Frustum = 2,
Bicone = 3,
}
impl TryFrom<u32> for ToolKind {
type Error = u32;
fn try_from(value: u32) -> Result<Self, Self::Error> {
match value {
0 => Ok(Self::Ball),
1 => Ok(Self::Cylinder),
2 => Ok(Self::Frustum),
3 => Ok(Self::Bicone),
other => Err(other),
}
}
}
#[repr(u32)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum ToolSweep {
None = 0,
AlongAxis = 1,
AcrossAxis = 2,
}
impl TryFrom<u32> for ToolSweep {
type Error = u32;
fn try_from(value: u32) -> Result<Self, Self::Error> {
match value {
0 => Ok(Self::None),
1 => Ok(Self::AlongAxis),
2 => Ok(Self::AcrossAxis),
other => Err(other),
}
}
}
#[repr(C, align(16))]
#[derive(Debug, Clone, Copy, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ToolPrimitive {
pub kind: u32,
pub sweep_kind: u32,
pub _pad: [u32; 2],
pub origin: [f32; 4],
pub axis: [f32; 4],
pub profile: [f32; 4],
pub sweep_dir: [f32; 4],
}
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub enum ToolPrimitiveError {
NonFinite,
NonPositiveRadius,
AxisNotUnit,
UnknownKind(u32),
UnknownSweep(u32),
PaddingNotZero,
}
impl std::fmt::Display for ToolPrimitiveError {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Self::NonFinite => f.write_str("tool primitive has a non-finite value"),
Self::NonPositiveRadius => {
f.write_str("tool primitive has an invalid radius or length")
}
Self::AxisNotUnit => f.write_str("tool primitive axis is not a unit vector"),
Self::UnknownKind(k) => write!(f, "unknown tool kind {k}"),
Self::UnknownSweep(k) => write!(f, "unknown tool sweep {k}"),
Self::PaddingNotZero => f.write_str("tool primitive padding is not zero"),
}
}
}
impl std::error::Error for ToolPrimitiveError {}
#[derive(Debug, Clone, Copy)]
pub struct StoneGeometry<'a> {
pub planes: &'a [GpuFacetPlane],
pub tools: &'a [ToolPrimitive],
}
impl<'a> StoneGeometry<'a> {
#[must_use]
pub const fn planes_only(planes: &'a [GpuFacetPlane]) -> Self {
Self { planes, tools: &[] }
}
#[must_use]
pub const fn is_convex(&self) -> bool {
self.tools.is_empty()
}
#[must_use]
pub const fn facet_count(&self) -> usize {
self.planes.len() + self.tools.len()
}
}
#[derive(Debug, Clone, Copy)]
enum Piece {
Ball {
c: Vec3,
r: f32,
},
Frustum {
c: Vec3,
a: Vec3,
hl: f32,
r0: f32,
r1: f32,
},
Wedge {
c: Vec3,
a: Vec3,
d: Vec3,
hl: f32,
s: f32,
r0: f32,
r1: f32,
},
}
const MAX_PIECES: usize = 6;
struct Pieces {
items: [Piece; MAX_PIECES],
len: usize,
}
impl Pieces {
const fn new() -> Self {
Self {
items: [Piece::Ball {
c: Vec3::ZERO,
r: 0.0,
}; MAX_PIECES],
len: 0,
}
}
fn push(&mut self, p: Piece) {
debug_assert!(
self.len < MAX_PIECES,
"tool decomposes into too many pieces"
);
if self.len < MAX_PIECES {
self.items[self.len] = p;
self.len += 1;
}
}
fn as_slice(&self) -> &[Piece] {
&self.items[..self.len]
}
}
#[derive(Clone, Copy)]
struct Segment {
z_lo: f32,
z_hi: f32,
r_lo: f32,
r_hi: f32,
}
impl Piece {
fn interval(&self, ray: Ray) -> Option<(f32, f32)> {
match *self {
Self::Ball { c, r } => ball_interval(ray, c, r),
Self::Frustum { c, a, hl, r0, r1 } => frustum_interval(ray, c, a, hl, r0, r1),
Self::Wedge {
c,
a,
d,
hl,
s,
r0,
r1,
} => wedge_interval(ray, c, a, d, hl, s, r0, r1),
}
}
#[expect(
clippy::manual_midpoint,
reason = "pinned reference arithmetic: plain f32 `0.5 * (a + b)` (`f32::midpoint` widens through f64)"
)]
#[allow(
clippy::suboptimal_flops,
reason = "pinned reference arithmetic: unfused products (`mul_add` rounds once and moves the last bit)"
)]
fn sd_normal(&self, p: Vec3) -> (f32, Vec3) {
match *self {
Self::Ball { c, r } => {
let q = p - c;
let len = q.length();
(len - r, if len > 1e-12 { q / len } else { Vec3::Z })
}
Self::Frustum { c, a, hl, r0, r1 } => {
let q = p - c;
let z = q.dot(a);
let perp = q - a * z;
let rho = perp.length();
let rho_hat = if rho > 1e-12 {
perp / rho
} else {
a.any_orthonormal_vector()
};
let (cap_sd, cap_n) = if z >= 0.0 { (z - hl, a) } else { (-z - hl, -a) };
let dr = r1 - r0;
let len = (4.0 * hl).mul_add(hl, dr * dr).sqrt();
let rm = 0.5 * (r0 + r1);
let side_sd = z.mul_add(-dr, (rho - rm) * 2.0 * hl) / len;
let side_n = rho_hat * (2.0 * hl / len) - a * (dr / len);
if side_sd > cap_sd {
(side_sd, side_n)
} else {
(cap_sd, cap_n)
}
}
Self::Wedge {
c,
a,
d,
hl,
s,
r0,
r1,
} => {
let q = p - c;
let w = a.cross(d);
let z = q.dot(a);
let x = q.dot(d);
let y = q.dot(w);
let k = (r1 - r0) / (2.0 * hl);
let rm = 0.5 * (r0 + r1);
let kappa = 1.0 / k.mul_add(k, 1.0).sqrt();
let mut best = if z >= 0.0 { (z - hl, a) } else { (-z - hl, -a) };
let along = if x >= 0.0 { (x - s, d) } else { (-x - s, -d) };
if along.0 > best.0 {
best = along;
}
let slant_pos = ((y - rm - k * z) * kappa, (w - a * k) * kappa);
if slant_pos.0 > best.0 {
best = slant_pos;
}
let slant_neg = ((-y - rm - k * z) * kappa, (-w - a * k) * kappa);
if slant_neg.0 > best.0 {
best = slant_neg;
}
best
}
}
}
}
impl ToolPrimitive {
const fn base(
kind: ToolKind,
centre: Vec3,
axis: Vec3,
rim: f32,
half_length: f32,
r: [f32; 2],
) -> Self {
Self {
kind: kind as u32,
sweep_kind: ToolSweep::None as u32,
_pad: [0; 2],
origin: [centre.x, centre.y, centre.z, rim],
axis: [axis.x, axis.y, axis.z, half_length],
profile: [r[0], r[1], 0.0, 0.0],
sweep_dir: [0.0; 4],
}
}
#[must_use]
pub const fn ball(centre: Vec3, radius: f32) -> Self {
Self::base(
ToolKind::Ball,
centre,
Vec3::Z,
radius,
0.0,
[radius, radius],
)
}
#[must_use]
pub const fn cylinder(centre: Vec3, axis: Vec3, radius: f32, half_length: f32) -> Self {
Self::base(
ToolKind::Cylinder,
centre,
axis,
radius,
half_length,
[radius, radius],
)
}
#[must_use]
#[expect(
clippy::manual_midpoint,
reason = "pinned reference arithmetic: plain f32 `0.5 * (r_neg + r_pos)`; `f32::midpoint` widens through f64"
)]
pub fn frustum(centre: Vec3, axis: Vec3, r_neg: f32, r_pos: f32, half_length: f32) -> Self {
Self::base(
ToolKind::Frustum,
centre,
axis,
0.5 * (r_neg + r_pos),
half_length,
[r_neg, r_pos],
)
}
#[must_use]
pub const fn bicone(centre: Vec3, axis: Vec3, rim_radius: f32, half_length: f32) -> Self {
Self::base(
ToolKind::Bicone,
centre,
axis,
rim_radius,
half_length,
[0.0, 0.0],
)
}
#[must_use]
pub fn with_sweep(mut self, sweep: ToolSweep, half_stroke: f32, dir: Vec3) -> Self {
self.sweep_kind = sweep as u32;
self.profile[2] = half_stroke;
self.sweep_dir = if sweep == ToolSweep::AcrossAxis {
[dir.x, dir.y, dir.z, 0.0]
} else {
[0.0; 4]
};
self
}
#[must_use]
pub fn kind(&self) -> Option<ToolKind> {
ToolKind::try_from(self.kind).ok()
}
#[must_use]
pub fn sweep(&self) -> Option<ToolSweep> {
ToolSweep::try_from(self.sweep_kind).ok()
}
pub fn validate(&self) -> Result<(), ToolPrimitiveError> {
let all = self
.origin
.iter()
.chain(&self.axis)
.chain(&self.profile)
.chain(&self.sweep_dir);
if !all.into_iter().all(|v| v.is_finite()) {
return Err(ToolPrimitiveError::NonFinite);
}
let kind = ToolKind::try_from(self.kind).map_err(ToolPrimitiveError::UnknownKind)?;
let sweep =
ToolSweep::try_from(self.sweep_kind).map_err(ToolPrimitiveError::UnknownSweep)?;
if self._pad != [0; 2] || self.profile[3] != 0.0 || self.sweep_dir[3] != 0.0 {
return Err(ToolPrimitiveError::PaddingNotZero);
}
let radii = [self.origin[3], self.profile[0], self.profile[1]];
let hl = self.axis[3];
let stroke = self.profile[2];
if radii.iter().any(|&r| r < 0.0)
|| radii.iter().all(|&r| r <= 0.0)
|| hl < 0.0
|| (kind != ToolKind::Ball && hl <= 0.0)
|| (kind == ToolKind::Ball && hl != 0.0)
|| stroke < 0.0
{
return Err(ToolPrimitiveError::NonPositiveRadius);
}
let axis = self.axis_vec();
if (axis.length() - 1.0).abs() > 1e-4 {
return Err(ToolPrimitiveError::AxisNotUnit);
}
if sweep == ToolSweep::AcrossAxis {
let dir = self.sweep_dir_vec();
if (dir.length() - 1.0).abs() > 1e-4 || dir.dot(axis).abs() > 1e-4 {
return Err(ToolPrimitiveError::AxisNotUnit);
}
}
Ok(())
}
const fn centre(&self) -> Vec3 {
Vec3::new(self.origin[0], self.origin[1], self.origin[2])
}
const fn axis_vec(&self) -> Vec3 {
Vec3::new(self.axis[0], self.axis[1], self.axis[2])
}
const fn sweep_dir_vec(&self) -> Vec3 {
Vec3::new(self.sweep_dir[0], self.sweep_dir[1], self.sweep_dir[2])
}
fn segments(&self, kind: ToolKind) -> ([Segment; 2], usize, f32, f32) {
let hl = self.axis[3];
let rim = self.origin[3];
let (rn, rp) = (self.profile[0], self.profile[1]);
let seg = |z_lo, z_hi, r_lo, r_hi| Segment {
z_lo,
z_hi,
r_lo,
r_hi,
};
match kind {
ToolKind::Ball | ToolKind::Cylinder => (
[seg(-hl, hl, rim, rim), seg(0.0, 0.0, 0.0, 0.0)],
1,
hl,
rim,
),
ToolKind::Frustum => {
let peak = if rp >= rn { (hl, rp) } else { (-hl, rn) };
(
[seg(-hl, hl, rn, rp), seg(0.0, 0.0, 0.0, 0.0)],
1,
peak.0,
peak.1,
)
}
ToolKind::Bicone => ([seg(-hl, 0.0, rn, rim), seg(0.0, hl, rim, rp)], 2, 0.0, rim),
}
}
#[expect(
clippy::manual_midpoint,
reason = "pinned reference arithmetic: plain f32 `0.5 * (z_lo + z_hi)`; `f32::midpoint` widens through f64"
)]
fn pieces(&self) -> Option<Pieces> {
let kind = self.kind()?;
let sweep = self.sweep()?;
let c = self.centre();
let a = self.axis_vec();
let stroke = self.profile[2];
let mut out = Pieces::new();
let swept = sweep != ToolSweep::None && stroke > 0.0;
if kind == ToolKind::Ball {
let r = self.origin[3];
if swept {
let dir = if sweep == ToolSweep::AlongAxis {
a
} else {
self.sweep_dir_vec()
};
out.push(Piece::Ball {
c: c - dir * stroke,
r,
});
out.push(Piece::Ball {
c: c + dir * stroke,
r,
});
out.push(Piece::Frustum {
c,
a: dir,
hl: stroke,
r0: r,
r1: r,
});
} else {
out.push(Piece::Ball { c, r });
}
return Some(out);
}
let (segs, n, peak_z, peak_r) = self.segments(kind);
let seg_piece = |s: &Segment, shift: Vec3| Piece::Frustum {
c: c + a * (0.5 * (s.z_lo + s.z_hi)) + shift,
a,
hl: 0.5 * (s.z_hi - s.z_lo),
r0: s.r_lo,
r1: s.r_hi,
};
for s in &segs[..n] {
if !swept {
out.push(seg_piece(s, Vec3::ZERO));
} else if sweep == ToolSweep::AlongAxis {
let shift = if s.r_hi >= s.r_lo {
-a * stroke
} else {
a * stroke
};
out.push(seg_piece(s, shift));
} else {
let dir = self.sweep_dir_vec();
out.push(seg_piece(s, -dir * stroke));
out.push(seg_piece(s, dir * stroke));
out.push(Piece::Wedge {
c: c + a * (0.5 * (s.z_lo + s.z_hi)),
a,
d: dir,
hl: 0.5 * (s.z_hi - s.z_lo),
s: stroke,
r0: s.r_lo,
r1: s.r_hi,
});
}
}
if swept && sweep == ToolSweep::AlongAxis {
out.push(Piece::Frustum {
c: c + a * peak_z,
a,
hl: stroke,
r0: peak_r,
r1: peak_r,
});
}
Some(out)
}
#[must_use]
pub fn contains(&self, p: Vec3) -> bool {
self.pieces()
.is_some_and(|ps| ps.as_slice().iter().any(|pc| pc.sd_normal(p).0 <= 0.0))
}
#[cfg(test)]
fn boundary_gap(&self, p: Vec3) -> f32 {
self.pieces().map_or(0.0, |ps| {
ps.as_slice()
.iter()
.map(|pc| pc.sd_normal(p).0.abs())
.fold(f32::INFINITY, f32::min)
})
}
#[must_use]
pub fn outward_normal(&self, p: Vec3) -> Vec3 {
let Some(ps) = self.pieces() else {
return Vec3::Z;
};
let mut best = (f32::INFINITY, Vec3::Z);
for pc in ps.as_slice() {
let (sd, n) = pc.sd_normal(p);
if sd.abs() < best.0 {
best = (sd.abs(), n);
}
}
best.1
}
#[must_use]
pub fn ray_interval(&self, ray: Ray) -> Option<(f32, f32)> {
let ps = self.pieces()?;
let mut hull: Option<(f32, f32)> = None;
for pc in ps.as_slice() {
if let Some((c, d)) = pc.interval(ray) {
hull = Some(hull.map_or((c, d), |(lo, hi)| (lo.min(c), hi.max(d))));
}
}
hull
}
}
#[cfg(test)]
mod tests;