use axiolid_contracts::{GeomError, GeomResult};
use axiolid_core::{Point3, Scalar, Vec3};
use axiolid_mesh::TriMesh;
use crate::loft::{self, Frame, Station};
use crate::profile::Rings;
use crate::revolve::rotate;
use crate::sweep;
mod band;
pub use band::MitreBand;
#[derive(Debug, Clone, Copy, PartialEq)]
#[non_exhaustive]
pub enum PathPiece {
Segment {
start: Point3,
end: Point3,
},
Arc {
centre: Point3,
axis: Vec3,
start: Point3,
angle: Scalar,
},
}
impl PathPiece {
#[must_use]
pub fn start_point(&self) -> Point3 {
match *self {
Self::Segment { start, .. } | Self::Arc { start, .. } => start,
}
}
#[must_use]
pub fn end_point(&self) -> Point3 {
self.point(1.0)
}
#[must_use]
pub fn point(&self, s: Scalar) -> Point3 {
match *self {
Self::Segment { start, end } => start + (end - start) * s,
Self::Arc {
centre,
axis,
start,
angle,
} => rotate(start, centre, axis, angle * s),
}
}
fn tangent_at(&self, point: Point3) -> Vec3 {
match *self {
Self::Segment { start, end } => (end - start).normalize_or_zero(),
Self::Arc {
centre,
axis,
angle,
..
} => (axis.cross(point - centre) * angle.signum()).normalize_or_zero(),
}
}
#[must_use]
pub fn start_tangent(&self) -> Vec3 {
self.tangent_at(self.start_point())
}
#[must_use]
pub fn end_tangent(&self) -> Vec3 {
self.tangent_at(self.end_point())
}
#[must_use]
pub fn length(&self) -> Scalar {
match *self {
Self::Segment { start, end } => (end - start).length(),
Self::Arc {
centre,
start,
angle,
..
} => (start - centre).length() * angle.abs(),
}
}
#[must_use]
pub fn reversed(&self) -> Self {
match *self {
Self::Segment { start, end } => Self::Segment {
start: end,
end: start,
},
Self::Arc {
centre,
axis,
angle,
..
} => Self::Arc {
centre,
axis,
start: self.end_point(),
angle: -angle,
},
}
}
#[must_use]
pub fn between(&self, from: Scalar, to: Scalar) -> Self {
match *self {
Self::Segment { .. } => Self::Segment {
start: self.point(from),
end: self.point(to),
},
Self::Arc {
centre,
axis,
angle,
..
} => Self::Arc {
centre,
axis,
start: self.point(from),
angle: angle * (to - from),
},
}
}
fn check(&self, radius: Scalar) -> GeomResult<()> {
let finite = |p: Point3| p.is_finite();
match *self {
Self::Segment { start, end } => {
if !(finite(start) && finite(end)) {
return Err(GeomError::InvalidInput(
"swept disk directrix segment must be finite".to_owned(),
));
}
if start == end {
return Err(GeomError::Degenerate(
"swept disk directrix segment has no length".to_owned(),
));
}
}
Self::Arc {
centre,
axis,
start,
angle,
} => {
if !(finite(centre) && finite(axis) && finite(start) && angle.is_finite()) {
return Err(GeomError::InvalidInput(
"swept disk directrix arc must be finite".to_owned(),
));
}
let bend = (start - centre).length();
if (axis.length() - 1.0).abs() > 1e-9
|| bend == 0.0
|| (start - centre).dot(axis).abs() > 1e-9 * bend
{
return Err(GeomError::InvalidInput(
"swept disk directrix arc needs a unit axis normal to a non-zero \
radius"
.to_owned(),
));
}
if angle == 0.0 || angle.abs() > core::f64::consts::TAU * (1.0 + 1e-12) {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix arc must turn by a non-zero angle of at most \
a full turn, got {angle}"
)));
}
if bend < radius {
return Err(GeomError::InvalidInput(format!(
"swept disk radius {radius} is not below the bend radius {bend}: \
the tube would fold through the bend's axis"
)));
}
if bend == radius {
return Err(horn_torus("bend", radius));
}
}
}
Ok(())
}
fn steps(&self, budget: Scalar) -> GeomResult<usize> {
match *self {
Self::Segment { .. } => Ok(1),
Self::Arc {
centre,
start,
angle,
..
} => crate::revolve::steps((start - centre).length(), angle, budget).map_err(|_| {
GeomError::BudgetExceeded {
resource: "swept disk bend steps",
}
}),
}
}
}
#[must_use]
pub fn joint_tolerance(chord: Scalar, radius: Scalar) -> Scalar {
2.0 * (chord / (8.0 * radius)).min(1.0).asin()
}
fn turn(a: Vec3, b: Vec3) -> Scalar {
a.cross(b).length().atan2(a.dot(b))
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[non_exhaustive]
pub enum MitreSide {
Segment {
corner: Point3,
direction: Vec3,
},
Arc(MitreBand),
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[non_exhaustive]
pub struct ArcMitre {
pub joint: usize,
pub normal: Vec3,
pub sides: [MitreSide; 2],
}
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub struct PipePlan {
pub pieces: Vec<PathPiece>,
pub closed: bool,
pub arc_mitres: Vec<ArcMitre>,
pub twist: Scalar,
}
pub fn swept_disk_along_pieces(
pieces: &[PathPiece],
radius: Scalar,
inner_radius: Option<Scalar>,
fillet_radius: Option<Scalar>,
chord: Scalar,
) -> GeomResult<TriMesh> {
Pipe::new(pieces, false, radius, inner_radius, fillet_radius, chord)?.mesh()
}
pub fn swept_disk_along_loop(
pieces: &[PathPiece],
radius: Scalar,
inner_radius: Option<Scalar>,
fillet_radius: Option<Scalar>,
chord: Scalar,
) -> GeomResult<TriMesh> {
Pipe::new(pieces, true, radius, inner_radius, fillet_radius, chord)?.mesh()
}
pub fn pipe_plan(
pieces: &[PathPiece],
closed: bool,
radius: Scalar,
inner_radius: Option<Scalar>,
fillet_radius: Option<Scalar>,
chord: Scalar,
) -> GeomResult<PipePlan> {
Ok(Pipe::new(pieces, closed, radius, inner_radius, fillet_radius, chord)?.plan())
}
const MAX_REFINEMENTS: usize = 16;
const MAX_TWIST_STEPS: usize = 4096;
const ROUNDING_TWIST: Scalar = 1e-3;
const STRIP_OFFSET: Scalar = 1.0 / 16.0;
struct Pipe {
path: Vec<PathPiece>,
closed: bool,
radius: Scalar,
inner: Option<Scalar>,
chord: Scalar,
mitres: Vec<Option<Vec3>>,
sharpest: Scalar,
bands: Vec<[Option<MitreBand>; 2]>,
twist: Scalar,
rate: Scalar,
}
struct Built {
stations: Vec<Station>,
end: Frame,
turned: Scalar,
}
#[derive(Clone, Copy)]
enum Close<'a> {
Open,
To(&'a Station),
Own,
}
struct Sampled {
points: Vec<Point3>,
tangents: Vec<Vec3>,
twists: Vec<Scalar>,
}
impl Pipe {
fn new(
pieces: &[PathPiece],
closed: bool,
radius: Scalar,
inner: Option<Scalar>,
fillet_radius: Option<Scalar>,
chord: Scalar,
) -> GeomResult<Self> {
sweep::check_radii(radius, inner)?;
if !(chord.is_finite() && chord > 0.0) {
return Err(GeomError::InvalidInput(format!(
"swept disk chord budget must be positive and finite, got {chord}"
)));
}
if let Some(f) = fillet_radius {
if !(f.is_finite() && f > 0.0) {
return Err(GeomError::InvalidInput(format!(
"swept disk fillet radius must be positive and finite, got {f}"
)));
}
}
if pieces.is_empty() {
return Err(GeomError::InvalidInput(
"a swept disk directrix needs at least one piece".to_owned(),
));
}
let tolerance = joint_tolerance(chord, radius);
let path = match fillet_radius {
Some(f) => fillet(pieces, f, radius, tolerance, closed)?,
None => pieces.to_vec(),
};
for piece in &path {
piece.check(radius)?;
}
let m = path.len();
let joints = if closed { m } else { m - 1 };
let mut mitres: Vec<Option<Vec3>> = Vec::with_capacity(joints);
let mut sharpest: Scalar = 0.0;
for k in 0..joints {
let (a, b) = (path[k], path[(k + 1) % m]);
let gap = (b.start_point() - a.end_point()).length();
if gap > 0.25 * chord {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix pieces leave a {gap} gap at joint {k}, more than a \
quarter of the {chord} chord budget"
)));
}
let (u, v) = (a.end_tangent(), b.start_tangent());
let theta = turn(u, v);
if theta <= tolerance {
mitres.push(None);
continue;
}
let straight = matches!(
(a, b),
(PathPiece::Segment { .. }, PathPiece::Segment { .. })
);
if fillet_radius.is_some() && !straight {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix turns a corner of {theta} rad beside an arc at joint \
{k}, above the {tolerance} rad a tangent-continuous joint may turn at this \
chord budget: only corners between two straight segments are filleted"
)));
}
let bisector = u + v;
if theta >= core::f64::consts::PI - 1e-9 || bisector.length() <= 1e-9 {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix reverses at joint {k}: a mitre at half angle would \
hold both legs"
)));
}
sharpest = sharpest.max(theta);
mitres.push(Some(bisector.normalize()));
}
let mut bands: Vec<[Option<MitreBand>; 2]> = vec![[None, None]; m];
for (k, normal) in mitres.iter().enumerate() {
let Some(n) = *normal else { continue };
for (i, end) in [(k, 1), ((k + 1) % m, 0)] {
let PathPiece::Arc {
centre,
axis,
angle,
..
} = path[i]
else {
continue;
};
let (toward, corner, normal) = if end == 1 {
(axis * angle.signum(), path[i].end_point(), n)
} else {
(-axis * angle.signum(), path[i].start_point(), -n)
};
let (back, ahead) = band::reach(centre, toward, corner, normal, radius)
.ok_or_else(|| {
GeomError::InvalidInput(format!(
"swept disk directrix mitres a corner beside an arc at joint {k}, but \
the mitre plane does not clear the bend's tube within a quarter turn \
of the corner: the bend is too tight for a corner this sharp"
))
})?;
bands[i][end] = Some(MitreBand {
centre,
axis: toward,
corner,
normal,
back,
ahead,
});
}
}
for (i, piece) in path.iter().enumerate() {
let PathPiece::Arc { angle, .. } = *piece else {
continue;
};
let need: Scalar = bands[i].iter().flatten().map(|b| b.back).sum();
if need > angle.abs() {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix arc {i} turns {} rad, but the mitre bands beside it \
need {need} rad of it: the mitres reach past the arc",
angle.abs()
)));
}
}
let mut pipe = Self {
path,
closed,
radius,
inner,
chord,
mitres,
sharpest,
bands,
twist: 0.0,
rate: 0.0,
};
pipe.check_segment_reach()?;
if closed {
pipe.twist = pipe.closing_twist();
if pipe.twist.abs() * radius > ROUNDING_TWIST * chord {
let clear: Scalar = (0..m).map(|i| pipe.clear(i).1).sum();
if clear <= 0.0 {
return Err(GeomError::InvalidInput(format!(
"swept disk along a closed directrix: a frame carried round the loop \
comes back turned by {} rad about the tangent, and no part of the \
loop is clear of its mitres to take that twist back",
pipe.twist
)));
}
pipe.rate = pipe.twist / clear;
}
}
Ok(pipe)
}
fn before(&self, i: usize) -> Option<usize> {
match i.checked_sub(1) {
Some(k) => Some(k),
None => self.closed.then(|| self.path.len() - 1),
}
}
fn after(&self, i: usize) -> Option<usize> {
(self.closed || i + 1 < self.path.len()).then_some(i)
}
fn mitre_ends(&self, i: usize) -> [Option<Vec3>; 2] {
[
self.before(i).and_then(|k| self.mitres[k]),
self.after(i).and_then(|k| self.mitres[k]),
]
}
fn check_segment_reach(&self) -> GeomResult<()> {
let radius = self.radius;
for (i, piece) in self.path.iter().enumerate() {
let PathPiece::Segment { start, end } = *piece else {
continue;
};
let ends = self.mitre_ends(i);
if ends == [None, None] {
continue;
}
let u = (end - start).normalize();
let lean = |n: Option<Vec3>| n.map_or(u, |n| n / u.dot(n));
let g = lean(ends[1]) - lean(ends[0]);
let reach = radius * (g - u * u.dot(g)).length();
let length = (end - start).length();
if reach >= length {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix segment {i} is {length} long, but the mitres at its \
ends reach {reach} along it at the disk radius {radius}: the tube would cut \
through itself"
)));
}
}
Ok(())
}
fn hold(&self, u: Vec3, n: Vec3) -> Scalar {
let c = u.dot(n).abs();
self.radius * (1.0 - c * c).max(0.0).sqrt() / c
}
fn clear(&self, i: usize) -> (Scalar, Scalar) {
match self.path[i] {
PathPiece::Segment { start, end } => {
let u = (end - start).normalize();
let [z0, z1] = self
.mitre_ends(i)
.map(|n| n.map_or(0.0, |n| self.hold(u, n)));
let clear = (end - start).length() - z0 - z1;
if clear > 0.0 {
(z0, clear)
} else {
(0.0, 0.0)
}
}
PathPiece::Arc {
centre,
start,
angle,
..
} => {
let banded: Scalar = self.bands[i].iter().flatten().map(|b| b.back).sum();
(
0.0,
(start - centre).length() * (angle.abs() - banded).max(0.0),
)
}
}
}
fn closing_twist(&self) -> Scalar {
let m = self.path.len();
let t0 = self.path[0].start_tangent();
let x0 = seed_reference(t0).normalize();
let mut x = x0;
for (i, piece) in self.path.iter().enumerate() {
if let PathPiece::Arc { axis, angle, .. } = *piece {
x = rotate(x, Point3::ZERO, axis, angle);
}
x = least_rotation(
x,
piece.end_tangent(),
self.path[(i + 1) % m].start_tangent(),
);
}
let x = (x - t0 * t0.dot(x)).normalize();
x.cross(x0).dot(t0).atan2(x.dot(x0))
}
fn plan(&self) -> PipePlan {
let m = self.path.len();
let side = |i: usize, end: usize| -> Option<MitreSide> {
let piece = &self.path[i];
match piece {
PathPiece::Segment { .. } => Some(MitreSide::Segment {
corner: if end == 1 {
piece.end_point()
} else {
piece.start_point()
},
direction: piece.start_tangent(),
}),
PathPiece::Arc { .. } => self.bands[i][end].map(MitreSide::Arc),
}
};
let arc_mitres = self
.mitres
.iter()
.enumerate()
.filter_map(|(k, normal)| {
let normal = (*normal)?;
let sides = [side(k, 1)?, side((k + 1) % m, 0)?];
matches!(sides, [MitreSide::Arc(_), _] | [_, MitreSide::Arc(_)]).then_some(
ArcMitre {
joint: k,
normal,
sides,
},
)
})
.collect();
PipePlan {
pieces: self.path.clone(),
closed: self.closed,
arc_mitres,
twist: self.twist,
}
}
fn mesh(&self) -> GeomResult<TriMesh> {
let rings = sweep::disk_rings(
self.radius,
self.inner,
0.5 * self.chord * (0.5 * self.sharpest).cos(),
)?;
let m = self.path.len();
let mut stations: Vec<Station> = Vec::new();
let mut carried: Option<Frame> = None;
let mut turned: Scalar = 0.0;
for i in 0..m {
let start = self.start_frame(i, carried.as_ref())?;
let closing = self.closed && i + 1 == m;
let close = match (closing, stations.first()) {
(false, _) => Close::Open,
(true, Some(first)) => Close::To(first),
(true, None) => Close::Own,
};
let built = self.piece_stations(i, &rings, &start, stations.last(), close, turned)?;
let skip = usize::from(!stations.is_empty());
let keep = built.stations.len() - usize::from(closing);
stations.extend(built.stations.into_iter().take(keep).skip(skip));
carried = Some(built.end);
turned = built.turned;
}
let mut mesh = loft::loft(&rings, &stations, self.closed)?;
if matches!(crate::extrude::outward_orientation(&mesh), Some(false)) {
for triangle in mesh.indices.as_chunks_mut::<3>().0 {
triangle.swap(1, 2);
}
}
Ok(mesh)
}
fn start_frame(&self, i: usize, carried: Option<&Frame>) -> GeomResult<Frame> {
let piece = &self.path[i];
let tangent = piece.start_tangent();
let seed = carried.map_or_else(
|| seed_reference(tangent),
|f| least_rotation(f.x, f.x.cross(f.y), tangent),
);
Frame::from_reference(piece.start_point(), tangent, seed)
}
fn piece_stations(
&self,
i: usize,
rings: &Rings,
start: &Frame,
joint: Option<&Station>,
close: Close<'_>,
turned: Scalar,
) -> GeomResult<Built> {
match self.path[i] {
PathPiece::Segment { .. } => {
self.segment_stations(i, rings, start, joint, close, turned)
}
PathPiece::Arc { .. } => self.arc_stations(i, rings, start, joint, close, turned),
}
}
fn segment_stations(
&self,
i: usize,
rings: &Rings,
start: &Frame,
joint: Option<&Station>,
close: Close<'_>,
turned: Scalar,
) -> GeomResult<Built> {
let piece = self.path[i];
let length = piece.length();
let tangent = piece.start_tangent();
let (held, clear) = self.clear(i);
let tail = if clear > 0.0 {
length - held - clear
} else {
0.0
};
let twist = self.rate * clear;
let sample = |budget: Scalar| -> GeomResult<Sampled> {
let mut at = vec![0.0];
if twist != 0.0 {
let n = twist_steps(twist, budget, self.radius)?;
if held > 0.0 {
at.push(held);
}
at.extend((1..n).map(|k| held + clear * k as Scalar / n as Scalar));
if tail > 0.0 {
at.push(length - tail);
}
}
at.push(length);
Ok(Sampled {
points: at.iter().map(|s| piece.point(s / length)).collect(),
tangents: vec![tangent; at.len()],
twists: at
.iter()
.map(|s| turned + self.rate * (s - held).clamp(0.0, clear))
.collect(),
})
};
let (frames, placed) = bounded_piece(
sample,
rings,
start.x,
self.radius,
0.5 * self.chord,
joint,
self.mitre_ends(i),
close,
)?;
Ok(Built {
stations: placed,
end: last_frame(frames)?,
turned: turned + twist,
})
}
fn arc_stations(
&self,
i: usize,
rings: &Rings,
start: &Frame,
joint: Option<&Station>,
close: Close<'_>,
turned: Scalar,
) -> GeomResult<Built> {
let piece = self.path[i];
let PathPiece::Arc {
centre,
start: from,
angle,
axis,
} = piece
else {
return Err(not_an_arc());
};
let total = angle.abs();
let along = axis * angle.signum();
let reach_radius = (from - centre).length() + self.radius;
let share = 0.5 * self.chord;
let band_tol = 0.5 * share;
let offset = STRIP_OFFSET * self.chord;
let [first, last] = self.bands[i];
let (head, tail) = (first.map_or(0.0, |b| b.back), last.map_or(0.0, |b| b.back));
let place = |frame: &Frame| loft::place(rings, |p| loft::at(&twisted(frame, turned), p));
let carried = |frame: &Frame, by: Scalar| Frame {
origin: rotate(frame.origin, centre, along, by),
x: rotate(frame.x, Point3::ZERO, along, by),
y: rotate(frame.y, Point3::ZERO, along, by),
};
let mut out: Vec<Station> = Vec::new();
let mut frame = start.clone();
if let Some(band) = first {
let k = self.before(i).ok_or_else(not_an_arc)?;
let shared = self.shared_ring(k, &twisted(start, turned), rings)?;
frame = carried(start, head);
let far = place(&frame);
let strip = matches!(self.path[k], PathPiece::Segment { .. });
let mut toward = band.stations_toward(
&far,
&shared,
strip.then_some(offset),
reach_radius,
band_tol,
)?;
toward.reverse();
out.extend(toward);
out.push(far);
}
let middle = total - head - tail;
let mut twist = 0.0;
if middle > 0.0 {
let part = piece.between(head / total, 1.0 - tail / total);
twist = self.rate * part.length();
let sample = |budget: Scalar| -> GeomResult<Sampled> {
let n = part
.steps(budget)?
.max(twist_steps(twist, budget, self.radius)?);
let points: Vec<Point3> = (0..=n)
.map(|k| part.point(k as Scalar / n as Scalar))
.collect();
let tangents = points.iter().map(|p| part.tangent_at(*p)).collect();
let twists = (0..=n)
.map(|k| turned + twist * k as Scalar / n as Scalar)
.collect();
Ok(Sampled {
points,
tangents,
twists,
})
};
let (joint, close) = (
if first.is_some() { out.last() } else { joint },
if last.is_some() { Close::Open } else { close },
);
let (frames, placed) = bounded_piece(
sample,
rings,
frame.x,
self.radius,
share,
joint,
[None, None],
close,
)?;
let skip = usize::from(first.is_some());
out.extend(placed.into_iter().skip(skip));
frame = last_frame(frames)?;
} else if out.is_empty() {
out.push(place(start));
}
let turned_end = turned + twist;
if let Some(band) = last {
let k = self.after(i).ok_or_else(not_an_arc)?;
let end = carried(&frame, tail);
let target = match close {
Close::To(first) => first.clone(),
Close::Own => out[0].clone(),
Close::Open => {
let next = self.start_frame((k + 1) % self.path.len(), Some(&end))?;
self.shared_ring(k, &twisted(&next, turned_end), rings)?
}
};
let far = out.last().cloned().ok_or_else(not_an_arc)?;
out.extend(band.stations_toward(
&far,
&target,
Some(offset),
reach_radius,
band_tol,
)?);
frame = end;
}
Ok(Built {
stations: out,
end: frame,
turned: turned_end,
})
}
fn shared_ring(&self, k: usize, frame: &Frame, rings: &Rings) -> GeomResult<Station> {
let next = (k + 1) % self.path.len();
let normal = self.mitres[k].ok_or_else(not_an_arc)?;
let square = loft::place(rings, |p| loft::at(frame, p));
match (self.path[k], self.path[next]) {
(PathPiece::Arc { .. }, PathPiece::Arc { .. }) => {
self.bands[next][0].ok_or_else(not_an_arc)?.cut(&square)
}
_ => Ok(mitred(&square, frame, normal)),
}
}
}
fn not_an_arc() -> GeomError {
GeomError::Degenerate("swept disk: a mitre band without its arc or corner".to_owned())
}
fn last_frame(frames: Vec<Frame>) -> GeomResult<Frame> {
frames
.into_iter()
.last()
.ok_or_else(|| GeomError::Degenerate("swept disk: a piece without a station".to_owned()))
}
fn twisted(frame: &Frame, angle: Scalar) -> Frame {
if angle == 0.0 {
return frame.clone();
}
let (s, c) = angle.sin_cos();
Frame {
origin: frame.origin,
x: frame.x * c + frame.y * s,
y: frame.y * c - frame.x * s,
}
}
fn twist_steps(twist: Scalar, budget: Scalar, radius: Scalar) -> GeomResult<usize> {
if twist == 0.0 {
return Ok(1);
}
let per = 2.0 * (1.0 - (budget / radius).min(1.0)).acos();
let n = (twist.abs() / per).ceil();
if n.is_nan() || n > MAX_TWIST_STEPS as Scalar {
return Err(GeomError::BudgetExceeded {
resource: "swept disk twist steps",
});
}
Ok((n as usize).max(1))
}
#[allow(clippy::too_many_arguments)]
fn bounded_piece(
sample: impl Fn(Scalar) -> GeomResult<Sampled>,
rings: &Rings,
seed: Vec3,
reach: Scalar,
share: Scalar,
joint: Option<&Station>,
ends: [Option<Vec3>; 2],
close: Close<'_>,
) -> GeomResult<(Vec<Frame>, Vec<Station>)> {
let mut b = share;
for _ in 0..MAX_REFINEMENTS {
let Sampled {
points,
tangents,
twists,
} = sample(b)?;
let frames = sweep::frames_carried(&points, &tangents, seed)?;
let placing: Vec<Frame> = frames
.iter()
.zip(&twists)
.map(|(f, a)| twisted(f, *a))
.collect();
let mut placed: Vec<Station> = placing
.iter()
.map(|f| loft::place(rings, |p| loft::at(f, p)))
.collect();
let last = placed.len() - 1;
for (k, normal) in [(0, ends[0]), (last, ends[1])] {
if let Some(n) = normal {
placed[k] = mitred(&placed[k], &placing[k], n);
}
}
let shift = joint.map_or(0.0, |j| displacement(j, &placed[0]));
let closing = match close {
Close::Open => 0.0,
Close::To(station) => displacement(station, &placed[last]),
Close::Own => displacement(&placed[0], &placed[last]),
};
let worst = placing
.windows(2)
.zip(placed.windows(2))
.enumerate()
.map(|(k, (f, s))| {
let moved = if k == 0 { shift } else { 0.0 };
let closed = if k + 1 == last { closing } else { 0.0 };
sweep::span_bound(&f[0], &f[1], b, reach)
+ loft::span_deviation(&s[0], &s[1])
+ moved
+ closed
})
.fold(0.0, Scalar::max);
if worst <= share {
return Ok((frames, placed));
}
b *= 0.5;
}
Err(GeomError::BudgetExceeded {
resource: "swept disk directrix refinement",
})
}
fn mitred(station: &Station, frame: &Frame, normal: Vec3) -> Station {
let tangent = frame.x.cross(frame.y);
let along = tangent.dot(normal);
Station {
loops: station
.loops
.iter()
.map(|ring| {
ring.iter()
.map(|q| *q - tangent * ((*q - frame.origin).dot(normal) / along))
.collect()
})
.collect(),
}
}
fn horn_torus(what: &str, radius: Scalar) -> GeomError {
GeomError::InvalidInput(format!(
"swept disk radius {radius} equals the {what} radius: the format rule permits a \
fillet radius equal to the disk radius, but the bend is then a horn torus whose \
inner wall degenerates to one point on the bend's axis, a boundary that is not a \
two-manifold and that no closed mesh bounds without meeting itself there"
))
}
fn displacement(a: &Station, b: &Station) -> Scalar {
a.loops
.iter()
.zip(&b.loops)
.flat_map(|(p, q)| p.iter().zip(q).map(|(u, v)| (*u - *v).length()))
.fold(0.0, Scalar::max)
}
fn least_rotation(v: Vec3, from: Vec3, to: Vec3) -> Vec3 {
let axis = from.cross(to);
let sin = axis.length();
if sin == 0.0 {
return v;
}
rotate(v, Point3::ZERO, axis / sin, turn(from, to))
}
fn seed_reference(t: Vec3) -> Vec3 {
let candidate = if t.x.abs() < 0.9 { Vec3::X } else { Vec3::Y };
candidate - t * t.dot(candidate)
}
fn fillet(
pieces: &[PathPiece],
fillet: Scalar,
radius: Scalar,
tolerance: Scalar,
closed: bool,
) -> GeomResult<Vec<PathPiece>> {
struct Leg {
piece: PathPiece,
dir: Vec3,
trim_start: Scalar,
trim_end: Scalar,
}
let mut legs: Vec<Leg> = Vec::with_capacity(pieces.len());
for piece in pieces {
piece.check(radius)?;
legs.push(Leg {
piece: *piece,
dir: piece.start_tangent(),
trim_start: 0.0,
trim_end: 0.0,
});
}
let m = legs.len();
let mut arcs: Vec<Option<PathPiece>> = vec![None; m];
let first = if closed { 0 } else { 1 };
for k in first..m {
let prev = (k + m - 1) % m;
if prev == k {
continue;
}
let (PathPiece::Segment { .. }, PathPiece::Segment { start: corner, .. }) =
(legs[prev].piece, legs[k].piece)
else {
continue;
};
let (u, v) = (legs[prev].dir, legs[k].dir);
let theta = turn(u, v);
if theta <= tolerance {
continue;
}
if radius > fillet {
return Err(GeomError::InvalidInput(format!(
"swept disk radius {radius} is above the fillet radius {fillet}: the \
tube would fold through the bend's axis"
)));
}
if radius == fillet {
return Err(horn_torus("fillet", radius));
}
let normal = u.cross(v);
if normal.length() <= 1e-12 {
return Err(GeomError::InvalidInput(format!(
"swept disk directrix reverses at corner {k}; no fillet rounds it"
)));
}
let cut = fillet * (0.5 * theta).tan();
legs[prev].trim_end = cut;
legs[k].trim_start = cut;
let tangent_point = corner - u * cut;
let inward = (v - u * u.dot(v)).normalize();
arcs[k] = Some(PathPiece::Arc {
centre: tangent_point + inward * fillet,
axis: normal.normalize(),
start: tangent_point,
angle: theta,
});
}
let mut out = Vec::with_capacity(2 * m);
for (k, leg) in legs.iter().enumerate() {
if k > 0 {
if let Some(arc) = arcs[k] {
out.push(arc);
}
}
let length = leg.piece.length();
let needed = leg.trim_start + leg.trim_end;
if needed == 0.0 {
out.push(leg.piece);
continue;
}
let slack = 1e-9 * length;
if needed > length + slack {
return Err(GeomError::InvalidInput(format!(
"swept disk fillet radius {fillet} does not fit: segment {k} is {length} \
long but its fillets need {needed} of it"
)));
}
if needed < length - slack {
out.push(
leg.piece
.between(leg.trim_start / length, 1.0 - leg.trim_end / length),
);
}
}
if closed {
if let Some(arc) = arcs[0] {
out.push(arc);
}
}
Ok(out)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn a_fillet_is_tangent_to_both_segments() {
let corner = Point3::new(1.0, 0.0, 0.0);
let pieces = [
PathPiece::Segment {
start: Point3::ZERO,
end: corner,
},
PathPiece::Segment {
start: corner,
end: Point3::new(1.0, 1.0, 1.0),
},
];
let path = fillet(&pieces, 0.2, 0.05, 1e-6, false).unwrap();
assert_eq!(path.len(), 3);
for pair in path.windows(2) {
assert!((pair[0].end_point() - pair[1].start_point()).length() < 1e-12);
assert!(turn(pair[0].end_tangent(), pair[1].start_tangent()) < 1e-12);
}
}
#[test]
fn a_mitre_ring_lies_on_both_cylinders_and_on_the_bisector_plane() {
let corner = Point3::new(1.0, 0.2, -0.3);
let (u, v) = (
Vec3::new(1.0, 0.0, 0.0),
Vec3::new(0.2, 0.7, -0.4).normalize(),
);
let pieces = [
PathPiece::Segment {
start: corner - u,
end: corner,
},
PathPiece::Segment {
start: corner,
end: corner + v,
},
];
let (r, inner) = (0.05, 0.03);
let mesh = swept_disk_along_pieces(&pieces, r, Some(inner), None, 1e-3).unwrap();
let per = mesh.positions.len() / 3;
let n = (u + v).normalize();
let off = |p: Point3, axis: Vec3| {
let d = p - corner;
(d - axis * d.dot(axis)).length()
};
for p in &mesh.positions[per..2 * per] {
assert!((*p - corner).dot(n).abs() < 1e-12);
let (a, b) = (off(*p, u), off(*p, v));
assert!((a - b).abs() < 1e-12, "{a} {b}");
assert!((a - r).abs() < 1e-12 || (a - inner).abs() < 1e-12, "{a}");
}
}
#[test]
fn an_arc_walked_backwards_retraces_itself() {
let arc = PathPiece::Arc {
centre: Point3::new(0.0, 1.0, 0.0),
axis: Vec3::Z,
start: Point3::ZERO,
angle: 1.0,
};
let back = arc.reversed();
assert!((back.end_point() - arc.start_point()).length() < 1e-12);
assert!((back.start_tangent() + arc.end_tangent()).length() < 1e-12);
let part = arc.between(0.25, 0.75);
assert!((part.length() - 0.5).abs() < 1e-12);
assert!((part.start_point() - arc.point(0.25)).length() < 1e-12);
}
}