Skip to main content

axiolid_construct/
pipe.rs

1//! Swept disks along chains of straight segments and circular arcs (#232,
2//! #245, #248).
3//!
4//! A pipe run, a rebar or a cable tray edge is a disk swept along legs
5//! joined by bends or corners: authoring formats write it as a disk swept
6//! along a composite of lines and arcs, or along a polyline whose corners
7//! are mitred or rounded to a fillet radius.
8//! [`swept_disk_along_pieces`] meshes all of them so that every point of
9//! the exact tube lies within the chord budget of the triangles.
10//!
11//! # The bound
12//!
13//! Each piece is bounded on its own, with the split [`crate::sweep`] uses
14//! for one arc (#231): half the budget `c` chords the disk, so a ring of
15//! `n` points inscribed in the circle of radius `r` lies within
16//! `r (1 - cos(pi/n)) <= c/2` of it, and every point of the exact tube is
17//! within `c/2` of the tube the ring polygon traces when it is carried
18//! rigidly along the exact centreline by the piece's own frames (each
19//! exact section circle and its polygon lie in the same plane). The
20//! other half bounds that polygon tube against the triangles:
21//!
22//! - Along a segment the frames do not turn, so the polygon tube is a
23//!   prism whose walls are the mesh's own planar rectangles: zero.
24//! - Along an arc of radius `R` the frames are rotation-minimising, which
25//!   along a circle means they turn with it: consecutive stations are
26//!   rotations of one ring about the arc's axis, and the ring's plane
27//!   contains that axis. The polygon tube between two stations is then the
28//!   revolution [`crate::revolve`] bounds, and its triangles are within
29//!   `(R + r)(1 - cos(phi/2))` of it for a step `phi`: the sagitta of the
30//!   outer side of the bend, not of the centreline. Its wall quads are
31//!   planar trapezoids, so their twist adds nothing; the span bound
32//!   [`crate::sweep`] computes, directrix sagitta plus the section's reach
33//!   plus [`crate::loft`]'s quad bound, is exactly that sum here. The step
34//!   halves until every span fits `c/2`; more than 4096 steps on a bend
35//!   is refused with `BudgetExceeded`.
36//!
37//! # Joints
38//!
39//! Consecutive pieces share one station, so the walls meet watertight and
40//! wound one way: the incoming piece's end section. The outgoing piece's
41//! own start frame is that frame turned by the least rotation that takes
42//! the incoming tangent to the outgoing one, moved to its start point, and
43//! its stations follow from it. Where the pieces are tangent continuous
44//! and meet, the two frames coincide. Where they turn by `theta` or leave
45//! a gap `g`, every vertex of the shared station lies at most
46//! `delta <= g + 2 r sin(theta/2)` from the station the outgoing piece's
47//! bound was proved for; moving a triangle's vertices by at most `delta`
48//! moves each of its points by at most `delta` (same barycentric weights),
49//! so the first span of the outgoing piece adds the measured `delta` to its
50//! bound and is refined until the sum fits.
51//!
52//! A joint turning by at most
53//!
54//! `theta_max = 2 asin(c / (8 r))` (pi once `c >= 8 r`)
55//!
56//! is tangent continuous at the requested accuracy: its shared section
57//! leans off the outgoing piece's own by at most a quarter of the budget,
58//! which that piece's first span pays as above. A gap above `c/4` is
59//! refused.
60//!
61//! # Mitres
62//!
63//! A sharper joint between two straight segments is a corner, and a disk
64//! swept round a corner is mitred at half angle (#245): both tubes are cut
65//! by the plane through the corner whose normal `n = (u + v) / |u + v|`
66//! bisects the incoming and outgoing tangents `u` and `v`. Reflection in
67//! that plane maps `u` to `-v`, so it swaps the two axis lines and with
68//! them the two cylinders of radius `r` (and the two bores); it fixes the
69//! plane, so both cylinders cut it in the same ellipse. The outgoing frame
70//! is the incoming one turned by the least rotation about `u x v`, which
71//! agrees with that reflection on the plane square to `u`. So profile point
72//! `k` of the incoming end section and of the outgoing start section, each
73//! projected along its own axis onto the mitre plane, land on one point:
74//! the two pieces share that mitre ring, and every vertex of it lies on
75//! both exact cylinders. Each wall quad of a segment joins two generator
76//! lines of the prism over the ring polygon, so it is a planar trapezoid,
77//! exactly a face of that prism cut by the end planes: the walls add
78//! nothing, and the span bound measures them as zero.
79//!
80//! The section pays instead. A point `w + s u` of the exact cylinder (`w`
81//! square to `u`, `|w| = r`; the mitre plane cuts its generator at
82//! `s_cut(w) = -w.n / u.n`, at most `r tan(theta/2)`) lies within
83//! `(1 - l) r <= e = r (1 - cos(pi/m))` of the prism point `w' + s u` at
84//! its azimuth, `w' = l w` on the ring polygon of `m` points. When that
85//! point lies past the plane, the plane crosses the prism's generator at
86//! `l s_cut(w)`, which moves the point by at most
87//! `(1 - l) s_cut(w) <= (1 - l) r tan(theta/2)` along `u`, square to the
88//! first move: the distance is at most `e / cos(theta/2)`, attained where
89//! the ellipse's major axis ends. The rings are therefore chorded to
90//! `(c/2) cos(theta_m/2)` for the sharpest mitre `theta_m` (only finer
91//! elsewhere), and a mitred pipe keeps every point of the exact tube
92//! within `c` of its triangles, as a smooth one does. A hollow disk's bore
93//! is mitred by the same plane and bounded the same way.
94//!
95//! Along a segment of length `L` the tube runs, on the generator at `w`,
96//! from its start plane to its end plane; both cuts are linear in `w`, so
97//! the shortest generator is `L - r |g_perp|` with
98//! `g = n1 / u.n1 - n0 / u.n0` (`n0 = u` at a square end), and it must be
99//! positive (one mitre: `L > r tan(theta/2)`). Refused by name: a segment
100//! whose mitres meet inside the tube, which would cut through itself, and a
101//! reversal (`theta = pi`, where the mitre plane holds both legs).
102//!
103//! # Mitres beside an arc
104//!
105//! A corner beside an arc is mitred at half angle too (#248), by the same
106//! bisector plane, but the plane cuts the bend's torus in a different curve
107//! than the tube on its other side: no ring lies on both exact surfaces,
108//! and the exact solid (each tube continued to the plane and cut by it) has
109//! a planar face between the two section curves. The end of the bend is
110//! meshed as a [`MitreBand`] (see its module): from a square section
111//! `back` short of the corner, each profile vertex walks along its own
112//! circle of the torus to the plane, so the band's vertices lie on the
113//! exact torus, and where the shared ring is not the band's own cut, the
114//! band stops just short of the plane and a strip, nearly in the plane,
115//! joins it to that ring: the mesh of the planar face. The shared ring is
116//! the outgoing piece's: a segment's mitre ring as above (so a segment
117//! beside an arc keeps the proof of the previous section, its ring on its
118//! own cylinder), or, between two arcs, the outgoing band's cut. Nothing
119//! proves the band and the strip within the budget by construction; the
120//! compiler certifies them against the mesh and reports that bound beside
121//! the proven one ([`pipe_plan`] gives it the bands). Refused by name: a
122//! corner whose plane does not clear the bend's tube within a quarter turn
123//! of it, and bands that need more of an arc than it turns.
124//!
125//! # Fillets
126//!
127//! A fillet radius `f` rounds every corner between two straight segments
128//! (the swept disk's `fillet_radius`): a corner turning by
129//! `theta` is replaced by the arc of radius `f` tangent to both segments,
130//! which cuts `f tan(theta/2)` off each. Segments shortened to nothing are
131//! dropped (two fillets then meet tangentially). Refused by name: a fillet
132//! that needs more of a segment than it has, a corner beside an arc (only
133//! corners between straight segments are rounded), a disk radius above the
134//! fillet radius, or any bend radius, whose tube folds through the bend's
135//! axis on its inside, and a reversal.
136//!
137//! A disk radius equal to the fillet (or bend) radius stays refused by name
138//! (#248), although the format rule permits it (fillet radius at least the
139//! disk radius): the bend is then a horn torus whose inner wall pinches to
140//! one point on the bend's axis, where every section touches. That boundary
141//! is not a two-manifold; a closed mesh within the budget would have to
142//! either meet itself there (a non-manifold vertex, degenerate triangles
143//! along the pinch) or stand off the pinch and misstate it.
144//!
145//! # Closed directrices
146//!
147//! [`swept_disk_along_loop`] sweeps a loop: the last piece ends where the
148//! first starts, and that joint is mitred, filleted or tangent like any
149//! other (#248). The walls close on the first station, so there are no
150//! caps. A frame carried round a loop that leaves its plane does not come
151//! back to itself: least rotations at the joints and rotation-minimising
152//! frames along the bends return it turned about the first tangent by the
153//! loop's holonomy `tau`, in `(-pi, pi]` ([`PipePlan::twist`], the
154//! measure reported). The pipe turns its frames back by `tau`, distributed
155//! uniformly over the loop's clear length: every segment's length outside
156//! the mitre zones at its ends (`r tan(theta/2)` from each mitred corner,
157//! where the plane reaches along it), and every bend's outside its mitre
158//! bands. In the zones and bands the frame is held, so a mitre still shares
159//! one ring on both cylinders and its walls stay planar. The twist is part
160//! of the proven bound: [`crate::sweep`]'s span bound measures the turn
161//! between two frames, twist included, and a twisted segment gets stations
162//! until each span's twist fits (more than 4096 is refused with
163//! `BudgetExceeded`). The station closing the loop is the first one; the
164//! last span pays its measured distance from the station it was proved for.
165//! Refused by name: a loop whose frame comes back turned when no part of it
166//! is clear to take the twist.
167
168use axiolid_contracts::{GeomError, GeomResult};
169use axiolid_core::{Point3, Scalar, Vec3};
170use axiolid_mesh::TriMesh;
171
172use crate::loft::{self, Frame, Station};
173use crate::profile::Rings;
174use crate::revolve::rotate;
175use crate::sweep;
176
177mod band;
178pub use band::MitreBand;
179
180/// One smooth piece of a swept disk's directrix.
181#[derive(Debug, Clone, Copy, PartialEq)]
182#[non_exhaustive]
183pub enum PathPiece {
184    /// A straight segment from `start` to `end`.
185    Segment {
186        /// First point, in the direction of travel.
187        start: Point3,
188        /// Last point.
189        end: Point3,
190    },
191    /// A circular arc: `start` turned by `angle` radians, right-handed,
192    /// about the line through `centre` along the unit `axis`. `start -
193    /// centre` must be perpendicular to `axis`; a negative `angle` turns
194    /// the other way.
195    Arc {
196        /// Centre of the arc's circle.
197        centre: Point3,
198        /// Unit normal of the arc's plane.
199        axis: Vec3,
200        /// First point, in the direction of travel.
201        start: Point3,
202        /// Signed turn, at most a full turn either way.
203        angle: Scalar,
204    },
205}
206
207impl PathPiece {
208    /// The first point, in the direction of travel.
209    #[must_use]
210    pub fn start_point(&self) -> Point3 {
211        match *self {
212            Self::Segment { start, .. } | Self::Arc { start, .. } => start,
213        }
214    }
215
216    /// The last point, in the direction of travel.
217    #[must_use]
218    pub fn end_point(&self) -> Point3 {
219        self.point(1.0)
220    }
221
222    /// The point a fraction `s` of the way along, by length.
223    #[must_use]
224    pub fn point(&self, s: Scalar) -> Point3 {
225        match *self {
226            Self::Segment { start, end } => start + (end - start) * s,
227            Self::Arc {
228                centre,
229                axis,
230                start,
231                angle,
232            } => rotate(start, centre, axis, angle * s),
233        }
234    }
235
236    /// The unit tangent, in the direction of travel, at `point` on the
237    /// piece; zero for a degenerate piece.
238    fn tangent_at(&self, point: Point3) -> Vec3 {
239        match *self {
240            Self::Segment { start, end } => (end - start).normalize_or_zero(),
241            Self::Arc {
242                centre,
243                axis,
244                angle,
245                ..
246            } => (axis.cross(point - centre) * angle.signum()).normalize_or_zero(),
247        }
248    }
249
250    /// The unit tangent at the start, in the direction of travel.
251    #[must_use]
252    pub fn start_tangent(&self) -> Vec3 {
253        self.tangent_at(self.start_point())
254    }
255
256    /// The unit tangent at the end, in the direction of travel.
257    #[must_use]
258    pub fn end_tangent(&self) -> Vec3 {
259        self.tangent_at(self.end_point())
260    }
261
262    /// The length along the piece.
263    #[must_use]
264    pub fn length(&self) -> Scalar {
265        match *self {
266            Self::Segment { start, end } => (end - start).length(),
267            Self::Arc {
268                centre,
269                start,
270                angle,
271                ..
272            } => (start - centre).length() * angle.abs(),
273        }
274    }
275
276    /// The same piece walked the other way.
277    #[must_use]
278    pub fn reversed(&self) -> Self {
279        match *self {
280            Self::Segment { start, end } => Self::Segment {
281                start: end,
282                end: start,
283            },
284            Self::Arc {
285                centre,
286                axis,
287                angle,
288                ..
289            } => Self::Arc {
290                centre,
291                axis,
292                start: self.end_point(),
293                angle: -angle,
294            },
295        }
296    }
297
298    /// The part between the fractions `from` and `to` of its length, in
299    /// that order (`to < from` walks it backwards).
300    #[must_use]
301    pub fn between(&self, from: Scalar, to: Scalar) -> Self {
302        match *self {
303            Self::Segment { .. } => Self::Segment {
304                start: self.point(from),
305                end: self.point(to),
306            },
307            Self::Arc {
308                centre,
309                axis,
310                angle,
311                ..
312            } => Self::Arc {
313                centre,
314                axis,
315                start: self.point(from),
316                angle: angle * (to - from),
317            },
318        }
319    }
320
321    /// Refuse a piece with no direction or a malformed arc.
322    fn check(&self, radius: Scalar) -> GeomResult<()> {
323        let finite = |p: Point3| p.is_finite();
324        match *self {
325            Self::Segment { start, end } => {
326                if !(finite(start) && finite(end)) {
327                    return Err(GeomError::InvalidInput(
328                        "swept disk directrix segment must be finite".to_owned(),
329                    ));
330                }
331                if start == end {
332                    return Err(GeomError::Degenerate(
333                        "swept disk directrix segment has no length".to_owned(),
334                    ));
335                }
336            }
337            Self::Arc {
338                centre,
339                axis,
340                start,
341                angle,
342            } => {
343                if !(finite(centre) && finite(axis) && finite(start) && angle.is_finite()) {
344                    return Err(GeomError::InvalidInput(
345                        "swept disk directrix arc must be finite".to_owned(),
346                    ));
347                }
348                let bend = (start - centre).length();
349                if (axis.length() - 1.0).abs() > 1e-9
350                    || bend == 0.0
351                    || (start - centre).dot(axis).abs() > 1e-9 * bend
352                {
353                    return Err(GeomError::InvalidInput(
354                        "swept disk directrix arc needs a unit axis normal to a non-zero \
355                         radius"
356                            .to_owned(),
357                    ));
358                }
359                if angle == 0.0 || angle.abs() > core::f64::consts::TAU * (1.0 + 1e-12) {
360                    return Err(GeomError::InvalidInput(format!(
361                        "swept disk directrix arc must turn by a non-zero angle of at most \
362                         a full turn, got {angle}"
363                    )));
364                }
365                if bend < radius {
366                    return Err(GeomError::InvalidInput(format!(
367                        "swept disk radius {radius} is not below the bend radius {bend}: \
368                         the tube would fold through the bend's axis"
369                    )));
370                }
371                if bend == radius {
372                    return Err(horn_torus("bend", radius));
373                }
374            }
375        }
376        Ok(())
377    }
378
379    /// Equal steps to chord budget `budget`: one for a segment, for an arc
380    /// as many as keep each step's sagitta within the budget.
381    fn steps(&self, budget: Scalar) -> GeomResult<usize> {
382        match *self {
383            Self::Segment { .. } => Ok(1),
384            Self::Arc {
385                centre,
386                start,
387                angle,
388                ..
389            } => crate::revolve::steps((start - centre).length(), angle, budget).map_err(|_| {
390                GeomError::BudgetExceeded {
391                    resource: "swept disk bend steps",
392                }
393            }),
394        }
395    }
396}
397
398/// The largest turn a joint may make and still be tangent continuous at
399/// chord budget `chord` for a disk of radius `radius`: `2 asin(c / (8 r))`
400/// (see the module notes).
401#[must_use]
402pub fn joint_tolerance(chord: Scalar, radius: Scalar) -> Scalar {
403    2.0 * (chord / (8.0 * radius)).min(1.0).asin()
404}
405
406/// The angle between two unit vectors, accurate when it is small.
407fn turn(a: Vec3, b: Vec3) -> Scalar {
408    a.cross(b).length().atan2(a.dot(b))
409}
410
411/// One side of a corner mitred beside an arc (#248).
412#[derive(Debug, Clone, Copy, PartialEq)]
413#[non_exhaustive]
414pub enum MitreSide {
415    /// A straight segment: near the corner its tube is the cylinder of the
416    /// disk radius about the line through `corner` along `direction`.
417    Segment {
418        /// The segment's end at the corner.
419        corner: Point3,
420        /// Its unit direction of travel.
421        direction: Vec3,
422    },
423    /// An arc: near the corner its tube is the band.
424    Arc(MitreBand),
425}
426
427/// A corner mitred beside an arc (#248): the mitre plane and the tubes on
428/// its two sides.
429#[derive(Debug, Clone, Copy, PartialEq)]
430#[non_exhaustive]
431pub struct ArcMitre {
432    /// The joint, between piece `joint` of [`PipePlan::pieces`] and the
433    /// next one (the first, for the joint that closes a loop).
434    pub joint: usize,
435    /// Unit normal of the mitre plane: the bisector of the incoming and the
436    /// outgoing tangent.
437    pub normal: Vec3,
438    /// The incoming side, then the outgoing one.
439    pub sides: [MitreSide; 2],
440}
441
442/// How [`swept_disk_along_pieces`] and [`swept_disk_along_loop`] sweep a
443/// directrix (#248): what a caller needs to bound the parts of the mesh no
444/// construction proof covers.
445#[derive(Debug, Clone, PartialEq)]
446#[non_exhaustive]
447pub struct PipePlan {
448    /// The pieces swept, fillets included.
449    pub pieces: Vec<PathPiece>,
450    /// Whether the last piece joins the first.
451    pub closed: bool,
452    /// Every corner mitred beside an arc, in joint order.
453    pub arc_mitres: Vec<ArcMitre>,
454    /// For a loop, the angle in `(-pi, pi]` about the first piece's start
455    /// tangent by which a frame carried once round the loop comes back
456    /// turned (its holonomy); the pipe turns its frames back by it,
457    /// uniformly over the loop's clear length (see the module notes). Zero
458    /// for an open directrix.
459    pub twist: Scalar,
460}
461
462/// Sweep a disk, optionally hollow, along a chain of segments and arcs so
463/// that every point of the exact tube lies within `chord` of the mesh.
464///
465/// `fillet_radius` rounds each corner between two segments; without it,
466/// each corner is mitred at half angle (#245, beside an arc #248). See the
467/// module notes for the bound, the joint rule, the mitres, and what is
468/// refused; a corner beside an arc is not proven but certified by the
469/// caller ([`pipe_plan`]). The pieces run in the direction of travel, each
470/// starting where the last ended.
471pub fn swept_disk_along_pieces(
472    pieces: &[PathPiece],
473    radius: Scalar,
474    inner_radius: Option<Scalar>,
475    fillet_radius: Option<Scalar>,
476    chord: Scalar,
477) -> GeomResult<TriMesh> {
478    Pipe::new(pieces, false, radius, inner_radius, fillet_radius, chord)?.mesh()
479}
480
481/// [`swept_disk_along_pieces`] round a loop (#248): the last piece ends
482/// where the first starts, that joint is mitred, filleted or tangent like
483/// any other, and the walls close on the first station without caps. A
484/// frame that comes back turned is twisted back uniformly over the loop's
485/// clear length (see the module notes and [`PipePlan::twist`]).
486pub fn swept_disk_along_loop(
487    pieces: &[PathPiece],
488    radius: Scalar,
489    inner_radius: Option<Scalar>,
490    fillet_radius: Option<Scalar>,
491    chord: Scalar,
492) -> GeomResult<TriMesh> {
493    Pipe::new(pieces, true, radius, inner_radius, fillet_radius, chord)?.mesh()
494}
495
496/// How [`swept_disk_along_pieces`] (`closed` false) or
497/// [`swept_disk_along_loop`] (`closed` true) sweeps these pieces at this
498/// budget, refused exactly as they refuse it.
499pub fn pipe_plan(
500    pieces: &[PathPiece],
501    closed: bool,
502    radius: Scalar,
503    inner_radius: Option<Scalar>,
504    fillet_radius: Option<Scalar>,
505    chord: Scalar,
506) -> GeomResult<PipePlan> {
507    Ok(Pipe::new(pieces, closed, radius, inner_radius, fillet_radius, chord)?.plan())
508}
509
510/// Rounds of halving a piece's directrix budget before it is refused.
511const MAX_REFINEMENTS: usize = 16;
512
513/// Most stations the twist of one piece may need.
514const MAX_TWIST_STEPS: usize = 4096;
515
516/// A closing twist whose turn at the disk's rim is below this fraction of
517/// the chord budget is rounding: the last span pays it as a measured
518/// shift instead of the loop being twisted.
519const ROUNDING_TWIST: Scalar = 1e-3;
520
521/// How far short of the mitre plane a band stops before its strip, as a
522/// fraction of the chord budget.
523const STRIP_OFFSET: Scalar = 1.0 / 16.0;
524
525/// A directrix checked and planned for sweeping.
526struct Pipe {
527    path: Vec<PathPiece>,
528    closed: bool,
529    radius: Scalar,
530    inner: Option<Scalar>,
531    chord: Scalar,
532    /// The mitre plane's unit normal at each joint, `None` where the joint
533    /// is tangent continuous at this budget. Joint `k` is between pieces
534    /// `k` and `k + 1` (the first, closing a loop).
535    mitres: Vec<Option<Vec3>>,
536    /// The sharpest mitre's turn.
537    sharpest: Scalar,
538    /// Each piece's mitre bands at its start and its end.
539    bands: Vec<[Option<MitreBand>; 2]>,
540    /// The closing twist, and the rate it is turned back at per unit of
541    /// clear length (zero where it is rounding).
542    twist: Scalar,
543    rate: Scalar,
544}
545
546/// A piece's stations, the untwisted frame at its end, and the twist
547/// reached there.
548struct Built {
549    stations: Vec<Station>,
550    end: Frame,
551    turned: Scalar,
552}
553
554/// Where a piece's last station must land.
555#[derive(Clone, Copy)]
556enum Close<'a> {
557    /// Anywhere: an open end, or a joint the next piece pays for.
558    Open,
559    /// On this station, the first of the loop.
560    To(&'a Station),
561    /// On the piece's own first station: a loop of one piece.
562    Own,
563}
564
565/// A piece sampled for a budget: points, the exact unit tangent at each,
566/// and the twist about it.
567struct Sampled {
568    points: Vec<Point3>,
569    tangents: Vec<Vec3>,
570    twists: Vec<Scalar>,
571}
572
573impl Pipe {
574    fn new(
575        pieces: &[PathPiece],
576        closed: bool,
577        radius: Scalar,
578        inner: Option<Scalar>,
579        fillet_radius: Option<Scalar>,
580        chord: Scalar,
581    ) -> GeomResult<Self> {
582        sweep::check_radii(radius, inner)?;
583        if !(chord.is_finite() && chord > 0.0) {
584            return Err(GeomError::InvalidInput(format!(
585                "swept disk chord budget must be positive and finite, got {chord}"
586            )));
587        }
588        if let Some(f) = fillet_radius {
589            if !(f.is_finite() && f > 0.0) {
590                return Err(GeomError::InvalidInput(format!(
591                    "swept disk fillet radius must be positive and finite, got {f}"
592                )));
593            }
594        }
595        if pieces.is_empty() {
596            return Err(GeomError::InvalidInput(
597                "a swept disk directrix needs at least one piece".to_owned(),
598            ));
599        }
600        let tolerance = joint_tolerance(chord, radius);
601        let path = match fillet_radius {
602            Some(f) => fillet(pieces, f, radius, tolerance, closed)?,
603            None => pieces.to_vec(),
604        };
605        for piece in &path {
606            piece.check(radius)?;
607        }
608        let m = path.len();
609        let joints = if closed { m } else { m - 1 };
610        let mut mitres: Vec<Option<Vec3>> = Vec::with_capacity(joints);
611        let mut sharpest: Scalar = 0.0;
612        for k in 0..joints {
613            let (a, b) = (path[k], path[(k + 1) % m]);
614            let gap = (b.start_point() - a.end_point()).length();
615            if gap > 0.25 * chord {
616                return Err(GeomError::InvalidInput(format!(
617                    "swept disk directrix pieces leave a {gap} gap at joint {k}, more than a \
618                     quarter of the {chord} chord budget"
619                )));
620            }
621            let (u, v) = (a.end_tangent(), b.start_tangent());
622            let theta = turn(u, v);
623            if theta <= tolerance {
624                mitres.push(None);
625                continue;
626            }
627            let straight = matches!(
628                (a, b),
629                (PathPiece::Segment { .. }, PathPiece::Segment { .. })
630            );
631            if fillet_radius.is_some() && !straight {
632                return Err(GeomError::InvalidInput(format!(
633                    "swept disk directrix turns a corner of {theta} rad beside an arc at joint \
634                     {k}, above the {tolerance} rad a tangent-continuous joint may turn at this \
635                     chord budget: only corners between two straight segments are filleted"
636                )));
637            }
638            let bisector = u + v;
639            if theta >= core::f64::consts::PI - 1e-9 || bisector.length() <= 1e-9 {
640                return Err(GeomError::InvalidInput(format!(
641                    "swept disk directrix reverses at joint {k}: a mitre at half angle would \
642                     hold both legs"
643                )));
644            }
645            sharpest = sharpest.max(theta);
646            mitres.push(Some(bisector.normalize()));
647        }
648        let mut bands: Vec<[Option<MitreBand>; 2]> = vec![[None, None]; m];
649        for (k, normal) in mitres.iter().enumerate() {
650            let Some(n) = *normal else { continue };
651            for (i, end) in [(k, 1), ((k + 1) % m, 0)] {
652                let PathPiece::Arc {
653                    centre,
654                    axis,
655                    angle,
656                    ..
657                } = path[i]
658                else {
659                    continue;
660                };
661                // Toward the corner: along the travel at the arc's end,
662                // against it at its start.
663                let (toward, corner, normal) = if end == 1 {
664                    (axis * angle.signum(), path[i].end_point(), n)
665                } else {
666                    (-axis * angle.signum(), path[i].start_point(), -n)
667                };
668                let (back, ahead) = band::reach(centre, toward, corner, normal, radius)
669                    .ok_or_else(|| {
670                        GeomError::InvalidInput(format!(
671                            "swept disk directrix mitres a corner beside an arc at joint {k}, but \
672                         the mitre plane does not clear the bend's tube within a quarter turn \
673                         of the corner: the bend is too tight for a corner this sharp"
674                        ))
675                    })?;
676                bands[i][end] = Some(MitreBand {
677                    centre,
678                    axis: toward,
679                    corner,
680                    normal,
681                    back,
682                    ahead,
683                });
684            }
685        }
686        for (i, piece) in path.iter().enumerate() {
687            let PathPiece::Arc { angle, .. } = *piece else {
688                continue;
689            };
690            let need: Scalar = bands[i].iter().flatten().map(|b| b.back).sum();
691            if need > angle.abs() {
692                return Err(GeomError::InvalidInput(format!(
693                    "swept disk directrix arc {i} turns {} rad, but the mitre bands beside it \
694                     need {need} rad of it: the mitres reach past the arc",
695                    angle.abs()
696                )));
697            }
698        }
699        let mut pipe = Self {
700            path,
701            closed,
702            radius,
703            inner,
704            chord,
705            mitres,
706            sharpest,
707            bands,
708            twist: 0.0,
709            rate: 0.0,
710        };
711        pipe.check_segment_reach()?;
712        if closed {
713            pipe.twist = pipe.closing_twist();
714            if pipe.twist.abs() * radius > ROUNDING_TWIST * chord {
715                let clear: Scalar = (0..m).map(|i| pipe.clear(i).1).sum();
716                if clear <= 0.0 {
717                    return Err(GeomError::InvalidInput(format!(
718                        "swept disk along a closed directrix: a frame carried round the loop \
719                         comes back turned by {} rad about the tangent, and no part of the \
720                         loop is clear of its mitres to take that twist back",
721                        pipe.twist
722                    )));
723                }
724                pipe.rate = pipe.twist / clear;
725            }
726        }
727        Ok(pipe)
728    }
729
730    /// The joint before piece `i`, if any.
731    fn before(&self, i: usize) -> Option<usize> {
732        match i.checked_sub(1) {
733            Some(k) => Some(k),
734            None => self.closed.then(|| self.path.len() - 1),
735        }
736    }
737
738    /// The joint after piece `i`, if any.
739    fn after(&self, i: usize) -> Option<usize> {
740        (self.closed || i + 1 < self.path.len()).then_some(i)
741    }
742
743    /// The mitre normals at the start and the end of piece `i`.
744    fn mitre_ends(&self, i: usize) -> [Option<Vec3>; 2] {
745        [
746            self.before(i).and_then(|k| self.mitres[k]),
747            self.after(i).and_then(|k| self.mitres[k]),
748        ]
749    }
750
751    /// A segment's two end planes must not meet inside the tube.
752    fn check_segment_reach(&self) -> GeomResult<()> {
753        let radius = self.radius;
754        for (i, piece) in self.path.iter().enumerate() {
755            let PathPiece::Segment { start, end } = *piece else {
756                continue;
757            };
758            let ends = self.mitre_ends(i);
759            if ends == [None, None] {
760                continue;
761            }
762            let u = (end - start).normalize();
763            let lean = |n: Option<Vec3>| n.map_or(u, |n| n / u.dot(n));
764            let g = lean(ends[1]) - lean(ends[0]);
765            let reach = radius * (g - u * u.dot(g)).length();
766            let length = (end - start).length();
767            if reach >= length {
768                return Err(GeomError::InvalidInput(format!(
769                    "swept disk directrix segment {i} is {length} long, but the mitres at its \
770                     ends reach {reach} along it at the disk radius {radius}: the tube would cut \
771                     through itself"
772                )));
773            }
774        }
775        Ok(())
776    }
777
778    /// How far a mitre plane with unit normal `n` through a segment's end
779    /// reaches along its direction `u` at the disk radius:
780    /// `r tan(theta/2)`.
781    fn hold(&self, u: Vec3, n: Vec3) -> Scalar {
782        let c = u.dot(n).abs();
783        self.radius * (1.0 - c * c).max(0.0).sqrt() / c
784    }
785
786    /// Piece `i`'s clear part, where a loop's twist is turned back: a
787    /// segment's length outside the mitre zones at its ends, as the length
788    /// held before it and its own length (both zero when the zones meet);
789    /// a bend's length outside its bands, with nothing held before it.
790    fn clear(&self, i: usize) -> (Scalar, Scalar) {
791        match self.path[i] {
792            PathPiece::Segment { start, end } => {
793                let u = (end - start).normalize();
794                let [z0, z1] = self
795                    .mitre_ends(i)
796                    .map(|n| n.map_or(0.0, |n| self.hold(u, n)));
797                let clear = (end - start).length() - z0 - z1;
798                if clear > 0.0 {
799                    (z0, clear)
800                } else {
801                    (0.0, 0.0)
802                }
803            }
804            PathPiece::Arc {
805                centre,
806                start,
807                angle,
808                ..
809            } => {
810                let banded: Scalar = self.bands[i].iter().flatten().map(|b| b.back).sum();
811                (
812                    0.0,
813                    (start - centre).length() * (angle.abs() - banded).max(0.0),
814                )
815            }
816        }
817    }
818
819    /// The holonomy of the loop: a frame carried round it (turning with
820    /// each bend, by the least rotation at each joint) comes back to the
821    /// first piece's start turned by this angle about its tangent, in
822    /// `(-pi, pi]`; turning the returned frame on by it gives the first.
823    fn closing_twist(&self) -> Scalar {
824        let m = self.path.len();
825        let t0 = self.path[0].start_tangent();
826        let x0 = seed_reference(t0).normalize();
827        let mut x = x0;
828        for (i, piece) in self.path.iter().enumerate() {
829            if let PathPiece::Arc { axis, angle, .. } = *piece {
830                x = rotate(x, Point3::ZERO, axis, angle);
831            }
832            x = least_rotation(
833                x,
834                piece.end_tangent(),
835                self.path[(i + 1) % m].start_tangent(),
836            );
837        }
838        let x = (x - t0 * t0.dot(x)).normalize();
839        x.cross(x0).dot(t0).atan2(x.dot(x0))
840    }
841
842    fn plan(&self) -> PipePlan {
843        let m = self.path.len();
844        let side = |i: usize, end: usize| -> Option<MitreSide> {
845            let piece = &self.path[i];
846            match piece {
847                PathPiece::Segment { .. } => Some(MitreSide::Segment {
848                    corner: if end == 1 {
849                        piece.end_point()
850                    } else {
851                        piece.start_point()
852                    },
853                    direction: piece.start_tangent(),
854                }),
855                PathPiece::Arc { .. } => self.bands[i][end].map(MitreSide::Arc),
856            }
857        };
858        let arc_mitres = self
859            .mitres
860            .iter()
861            .enumerate()
862            .filter_map(|(k, normal)| {
863                let normal = (*normal)?;
864                let sides = [side(k, 1)?, side((k + 1) % m, 0)?];
865                matches!(sides, [MitreSide::Arc(_), _] | [_, MitreSide::Arc(_)]).then_some(
866                    ArcMitre {
867                        joint: k,
868                        normal,
869                        sides,
870                    },
871                )
872            })
873            .collect();
874        PipePlan {
875            pieces: self.path.clone(),
876            closed: self.closed,
877            arc_mitres,
878            twist: self.twist,
879        }
880    }
881
882    fn mesh(&self) -> GeomResult<TriMesh> {
883        let rings = sweep::disk_rings(
884            self.radius,
885            self.inner,
886            0.5 * self.chord * (0.5 * self.sharpest).cos(),
887        )?;
888        let m = self.path.len();
889        let mut stations: Vec<Station> = Vec::new();
890        let mut carried: Option<Frame> = None;
891        let mut turned: Scalar = 0.0;
892        for i in 0..m {
893            let start = self.start_frame(i, carried.as_ref())?;
894            let closing = self.closed && i + 1 == m;
895            let close = match (closing, stations.first()) {
896                (false, _) => Close::Open,
897                (true, Some(first)) => Close::To(first),
898                (true, None) => Close::Own,
899            };
900            let built = self.piece_stations(i, &rings, &start, stations.last(), close, turned)?;
901            // Consecutive pieces share one station; a loop's last one is
902            // its first.
903            let skip = usize::from(!stations.is_empty());
904            let keep = built.stations.len() - usize::from(closing);
905            stations.extend(built.stations.into_iter().take(keep).skip(skip));
906            carried = Some(built.end);
907            turned = built.turned;
908        }
909        let mut mesh = loft::loft(&rings, &stations, self.closed)?;
910        // As `revolve::revolve` (#221): settle the winding from the built
911        // mesh rather than from the frames' handedness.
912        if matches!(crate::extrude::outward_orientation(&mesh), Some(false)) {
913            for triangle in mesh.indices.as_chunks_mut::<3>().0 {
914                triangle.swap(1, 2);
915            }
916        }
917        Ok(mesh)
918    }
919
920    /// Piece `i`'s own start frame: the incoming end frame turned by the
921    /// least rotation onto its tangent.
922    fn start_frame(&self, i: usize, carried: Option<&Frame>) -> GeomResult<Frame> {
923        let piece = &self.path[i];
924        let tangent = piece.start_tangent();
925        let seed = carried.map_or_else(
926            || seed_reference(tangent),
927            |f| least_rotation(f.x, f.x.cross(f.y), tangent),
928        );
929        Frame::from_reference(piece.start_point(), tangent, seed)
930    }
931
932    fn piece_stations(
933        &self,
934        i: usize,
935        rings: &Rings,
936        start: &Frame,
937        joint: Option<&Station>,
938        close: Close<'_>,
939        turned: Scalar,
940    ) -> GeomResult<Built> {
941        match self.path[i] {
942            PathPiece::Segment { .. } => {
943                self.segment_stations(i, rings, start, joint, close, turned)
944            }
945            PathPiece::Arc { .. } => self.arc_stations(i, rings, start, joint, close, turned),
946        }
947    }
948
949    /// A segment's stations: its two ends (mitred where a corner is), and
950    /// in a twisted loop the ends of its mitre zones and enough stations
951    /// between them for the twist.
952    fn segment_stations(
953        &self,
954        i: usize,
955        rings: &Rings,
956        start: &Frame,
957        joint: Option<&Station>,
958        close: Close<'_>,
959        turned: Scalar,
960    ) -> GeomResult<Built> {
961        let piece = self.path[i];
962        let length = piece.length();
963        let tangent = piece.start_tangent();
964        let (held, clear) = self.clear(i);
965        let tail = if clear > 0.0 {
966            length - held - clear
967        } else {
968            0.0
969        };
970        let twist = self.rate * clear;
971        let sample = |budget: Scalar| -> GeomResult<Sampled> {
972            let mut at = vec![0.0];
973            if twist != 0.0 {
974                let n = twist_steps(twist, budget, self.radius)?;
975                if held > 0.0 {
976                    at.push(held);
977                }
978                at.extend((1..n).map(|k| held + clear * k as Scalar / n as Scalar));
979                if tail > 0.0 {
980                    at.push(length - tail);
981                }
982            }
983            at.push(length);
984            Ok(Sampled {
985                points: at.iter().map(|s| piece.point(s / length)).collect(),
986                tangents: vec![tangent; at.len()],
987                twists: at
988                    .iter()
989                    .map(|s| turned + self.rate * (s - held).clamp(0.0, clear))
990                    .collect(),
991            })
992        };
993        let (frames, placed) = bounded_piece(
994            sample,
995            rings,
996            start.x,
997            self.radius,
998            0.5 * self.chord,
999            joint,
1000            self.mitre_ends(i),
1001            close,
1002        )?;
1003        Ok(Built {
1004            stations: placed,
1005            end: last_frame(frames)?,
1006            turned: turned + twist,
1007        })
1008    }
1009
1010    /// An arc's stations: a mitre band at each end beside a mitred corner,
1011    /// and the part between proven like any bend.
1012    fn arc_stations(
1013        &self,
1014        i: usize,
1015        rings: &Rings,
1016        start: &Frame,
1017        joint: Option<&Station>,
1018        close: Close<'_>,
1019        turned: Scalar,
1020    ) -> GeomResult<Built> {
1021        let piece = self.path[i];
1022        let PathPiece::Arc {
1023            centre,
1024            start: from,
1025            angle,
1026            axis,
1027        } = piece
1028        else {
1029            return Err(not_an_arc());
1030        };
1031        let total = angle.abs();
1032        let along = axis * angle.signum();
1033        let reach_radius = (from - centre).length() + self.radius;
1034        let share = 0.5 * self.chord;
1035        let band_tol = 0.5 * share;
1036        let offset = STRIP_OFFSET * self.chord;
1037        let [first, last] = self.bands[i];
1038        let (head, tail) = (first.map_or(0.0, |b| b.back), last.map_or(0.0, |b| b.back));
1039        let place = |frame: &Frame| loft::place(rings, |p| loft::at(&twisted(frame, turned), p));
1040        let carried = |frame: &Frame, by: Scalar| Frame {
1041            origin: rotate(frame.origin, centre, along, by),
1042            x: rotate(frame.x, Point3::ZERO, along, by),
1043            y: rotate(frame.y, Point3::ZERO, along, by),
1044        };
1045        let mut out: Vec<Station> = Vec::new();
1046        let mut frame = start.clone();
1047        if let Some(band) = first {
1048            let k = self.before(i).ok_or_else(not_an_arc)?;
1049            let shared = self.shared_ring(k, &twisted(start, turned), rings)?;
1050            frame = carried(start, head);
1051            let far = place(&frame);
1052            // The shared ring is this band's own cut unless a segment
1053            // comes in.
1054            let strip = matches!(self.path[k], PathPiece::Segment { .. });
1055            let mut toward = band.stations_toward(
1056                &far,
1057                &shared,
1058                strip.then_some(offset),
1059                reach_radius,
1060                band_tol,
1061            )?;
1062            toward.reverse();
1063            out.extend(toward);
1064            out.push(far);
1065        }
1066        let middle = total - head - tail;
1067        let mut twist = 0.0;
1068        if middle > 0.0 {
1069            let part = piece.between(head / total, 1.0 - tail / total);
1070            twist = self.rate * part.length();
1071            let sample = |budget: Scalar| -> GeomResult<Sampled> {
1072                let n = part
1073                    .steps(budget)?
1074                    .max(twist_steps(twist, budget, self.radius)?);
1075                let points: Vec<Point3> = (0..=n)
1076                    .map(|k| part.point(k as Scalar / n as Scalar))
1077                    .collect();
1078                let tangents = points.iter().map(|p| part.tangent_at(*p)).collect();
1079                let twists = (0..=n)
1080                    .map(|k| turned + twist * k as Scalar / n as Scalar)
1081                    .collect();
1082                Ok(Sampled {
1083                    points,
1084                    tangents,
1085                    twists,
1086                })
1087            };
1088            let (joint, close) = (
1089                if first.is_some() { out.last() } else { joint },
1090                if last.is_some() { Close::Open } else { close },
1091            );
1092            let (frames, placed) = bounded_piece(
1093                sample,
1094                rings,
1095                frame.x,
1096                self.radius,
1097                share,
1098                joint,
1099                [None, None],
1100                close,
1101            )?;
1102            let skip = usize::from(first.is_some());
1103            out.extend(placed.into_iter().skip(skip));
1104            frame = last_frame(frames)?;
1105        } else if out.is_empty() {
1106            out.push(place(start));
1107        }
1108        let turned_end = turned + twist;
1109        if let Some(band) = last {
1110            let k = self.after(i).ok_or_else(not_an_arc)?;
1111            let end = carried(&frame, tail);
1112            let target = match close {
1113                Close::To(first) => first.clone(),
1114                Close::Own => out[0].clone(),
1115                Close::Open => {
1116                    let next = self.start_frame((k + 1) % self.path.len(), Some(&end))?;
1117                    self.shared_ring(k, &twisted(&next, turned_end), rings)?
1118                }
1119            };
1120            let far = out.last().cloned().ok_or_else(not_an_arc)?;
1121            out.extend(band.stations_toward(
1122                &far,
1123                &target,
1124                Some(offset),
1125                reach_radius,
1126                band_tol,
1127            )?);
1128            frame = end;
1129        }
1130        Ok(Built {
1131            stations: out,
1132            end: frame,
1133            turned: turned_end,
1134        })
1135    }
1136
1137    /// The ring shared across mitred joint `k`, from the outgoing piece's
1138    /// start frame `frame` (twist included): a segment's mitre ring, or
1139    /// after a segment an arc's section projected the same way (the
1140    /// incoming segment's own ring, #245); between two arcs the outgoing
1141    /// band's cut.
1142    fn shared_ring(&self, k: usize, frame: &Frame, rings: &Rings) -> GeomResult<Station> {
1143        let next = (k + 1) % self.path.len();
1144        let normal = self.mitres[k].ok_or_else(not_an_arc)?;
1145        let square = loft::place(rings, |p| loft::at(frame, p));
1146        match (self.path[k], self.path[next]) {
1147            (PathPiece::Arc { .. }, PathPiece::Arc { .. }) => {
1148                self.bands[next][0].ok_or_else(not_an_arc)?.cut(&square)
1149            }
1150            _ => Ok(mitred(&square, frame, normal)),
1151        }
1152    }
1153}
1154
1155/// The refusal for a plan that does not hold together; never reached.
1156fn not_an_arc() -> GeomError {
1157    GeomError::Degenerate("swept disk: a mitre band without its arc or corner".to_owned())
1158}
1159
1160fn last_frame(frames: Vec<Frame>) -> GeomResult<Frame> {
1161    frames
1162        .into_iter()
1163        .last()
1164        .ok_or_else(|| GeomError::Degenerate("swept disk: a piece without a station".to_owned()))
1165}
1166
1167/// `frame` turned by `angle` about its own tangent.
1168fn twisted(frame: &Frame, angle: Scalar) -> Frame {
1169    if angle == 0.0 {
1170        return frame.clone();
1171    }
1172    let (s, c) = angle.sin_cos();
1173    Frame {
1174        origin: frame.origin,
1175        x: frame.x * c + frame.y * s,
1176        y: frame.y * c - frame.x * s,
1177    }
1178}
1179
1180/// Spans a twist of `twist` needs for the turn between two frames to move
1181/// the disk's rim by at most `budget`: `r (1 - cos(phi/2)) <= budget`.
1182fn twist_steps(twist: Scalar, budget: Scalar, radius: Scalar) -> GeomResult<usize> {
1183    if twist == 0.0 {
1184        return Ok(1);
1185    }
1186    let per = 2.0 * (1.0 - (budget / radius).min(1.0)).acos();
1187    let n = (twist.abs() / per).ceil();
1188    if n.is_nan() || n > MAX_TWIST_STEPS as Scalar {
1189        return Err(GeomError::BudgetExceeded {
1190            resource: "swept disk twist steps",
1191        });
1192    }
1193    Ok((n as usize).max(1))
1194}
1195
1196/// A piece's frames (untwisted) and stations, its directrix sampled until
1197/// each span's bound fits `share`; the first span also carries the
1198/// measured distance from the shared `joint` station to the piece's own
1199/// first station, and the last the distance from its last station to the
1200/// one it `close`s on. An end with a mitre normal in `ends` (start, end)
1201/// has its station projected along the piece's tangent onto that mitre
1202/// plane.
1203#[allow(clippy::too_many_arguments)]
1204fn bounded_piece(
1205    sample: impl Fn(Scalar) -> GeomResult<Sampled>,
1206    rings: &Rings,
1207    seed: Vec3,
1208    reach: Scalar,
1209    share: Scalar,
1210    joint: Option<&Station>,
1211    ends: [Option<Vec3>; 2],
1212    close: Close<'_>,
1213) -> GeomResult<(Vec<Frame>, Vec<Station>)> {
1214    let mut b = share;
1215    for _ in 0..MAX_REFINEMENTS {
1216        let Sampled {
1217            points,
1218            tangents,
1219            twists,
1220        } = sample(b)?;
1221        let frames = sweep::frames_carried(&points, &tangents, seed)?;
1222        let placing: Vec<Frame> = frames
1223            .iter()
1224            .zip(&twists)
1225            .map(|(f, a)| twisted(f, *a))
1226            .collect();
1227        let mut placed: Vec<Station> = placing
1228            .iter()
1229            .map(|f| loft::place(rings, |p| loft::at(f, p)))
1230            .collect();
1231        let last = placed.len() - 1;
1232        for (k, normal) in [(0, ends[0]), (last, ends[1])] {
1233            if let Some(n) = normal {
1234                placed[k] = mitred(&placed[k], &placing[k], n);
1235            }
1236        }
1237        let shift = joint.map_or(0.0, |j| displacement(j, &placed[0]));
1238        let closing = match close {
1239            Close::Open => 0.0,
1240            Close::To(station) => displacement(station, &placed[last]),
1241            Close::Own => displacement(&placed[0], &placed[last]),
1242        };
1243        let worst = placing
1244            .windows(2)
1245            .zip(placed.windows(2))
1246            .enumerate()
1247            .map(|(k, (f, s))| {
1248                let moved = if k == 0 { shift } else { 0.0 };
1249                let closed = if k + 1 == last { closing } else { 0.0 };
1250                sweep::span_bound(&f[0], &f[1], b, reach)
1251                    + loft::span_deviation(&s[0], &s[1])
1252                    + moved
1253                    + closed
1254            })
1255            .fold(0.0, Scalar::max);
1256        if worst <= share {
1257            return Ok((frames, placed));
1258        }
1259        b *= 0.5;
1260    }
1261    Err(GeomError::BudgetExceeded {
1262        resource: "swept disk directrix refinement",
1263    })
1264}
1265
1266/// `station`, square to `frame`'s tangent at its origin, projected along
1267/// that tangent onto the mitre plane through the origin with unit normal
1268/// `normal` (see the module notes).
1269fn mitred(station: &Station, frame: &Frame, normal: Vec3) -> Station {
1270    let tangent = frame.x.cross(frame.y);
1271    let along = tangent.dot(normal);
1272    Station {
1273        loops: station
1274            .loops
1275            .iter()
1276            .map(|ring| {
1277                ring.iter()
1278                    .map(|q| *q - tangent * ((*q - frame.origin).dot(normal) / along))
1279                    .collect()
1280            })
1281            .collect(),
1282    }
1283}
1284
1285/// The refusal of a disk radius equal to a bend's or fillet's radius
1286/// (see the module notes).
1287fn horn_torus(what: &str, radius: Scalar) -> GeomError {
1288    GeomError::InvalidInput(format!(
1289        "swept disk radius {radius} equals the {what} radius: the format rule permits a \
1290         fillet radius equal to the disk radius, but the bend is then a horn torus whose \
1291         inner wall degenerates to one point on the bend's axis, a boundary that is not a \
1292         two-manifold and that no closed mesh bounds without meeting itself there"
1293    ))
1294}
1295
1296/// The largest distance between corresponding vertices of two stations.
1297fn displacement(a: &Station, b: &Station) -> Scalar {
1298    a.loops
1299        .iter()
1300        .zip(&b.loops)
1301        .flat_map(|(p, q)| p.iter().zip(q).map(|(u, v)| (*u - *v).length()))
1302        .fold(0.0, Scalar::max)
1303}
1304
1305/// `v` turned by the least rotation that takes unit `from` to unit `to`.
1306fn least_rotation(v: Vec3, from: Vec3, to: Vec3) -> Vec3 {
1307    let axis = from.cross(to);
1308    let sin = axis.length();
1309    if sin == 0.0 {
1310        return v;
1311    }
1312    rotate(v, Point3::ZERO, axis / sin, turn(from, to))
1313}
1314
1315/// Any direction perpendicular to unit `t`.
1316fn seed_reference(t: Vec3) -> Vec3 {
1317    let candidate = if t.x.abs() < 0.9 { Vec3::X } else { Vec3::Y };
1318    candidate - t * t.dot(candidate)
1319}
1320
1321/// Round every corner between two segments that turns by more than
1322/// `tolerance` with an arc of radius `fillet` (see the module notes); in a
1323/// loop, the corner where it closes too.
1324fn fillet(
1325    pieces: &[PathPiece],
1326    fillet: Scalar,
1327    radius: Scalar,
1328    tolerance: Scalar,
1329    closed: bool,
1330) -> GeomResult<Vec<PathPiece>> {
1331    // Each segment keeps its direction even when the fillets either side
1332    // use it up, so the corner after it is still measured against it.
1333    struct Leg {
1334        piece: PathPiece,
1335        dir: Vec3,
1336        trim_start: Scalar,
1337        trim_end: Scalar,
1338    }
1339    let mut legs: Vec<Leg> = Vec::with_capacity(pieces.len());
1340    for piece in pieces {
1341        piece.check(radius)?;
1342        legs.push(Leg {
1343            piece: *piece,
1344            dir: piece.start_tangent(),
1345            trim_start: 0.0,
1346            trim_end: 0.0,
1347        });
1348    }
1349    let m = legs.len();
1350    // The fillet at the corner where leg `k` starts.
1351    let mut arcs: Vec<Option<PathPiece>> = vec![None; m];
1352    let first = if closed { 0 } else { 1 };
1353    for k in first..m {
1354        let prev = (k + m - 1) % m;
1355        if prev == k {
1356            continue;
1357        }
1358        let (PathPiece::Segment { .. }, PathPiece::Segment { start: corner, .. }) =
1359            (legs[prev].piece, legs[k].piece)
1360        else {
1361            continue;
1362        };
1363        let (u, v) = (legs[prev].dir, legs[k].dir);
1364        let theta = turn(u, v);
1365        if theta <= tolerance {
1366            continue;
1367        }
1368        if radius > fillet {
1369            return Err(GeomError::InvalidInput(format!(
1370                "swept disk radius {radius} is above the fillet radius {fillet}: the \
1371                 tube would fold through the bend's axis"
1372            )));
1373        }
1374        if radius == fillet {
1375            return Err(horn_torus("fillet", radius));
1376        }
1377        let normal = u.cross(v);
1378        if normal.length() <= 1e-12 {
1379            return Err(GeomError::InvalidInput(format!(
1380                "swept disk directrix reverses at corner {k}; no fillet rounds it"
1381            )));
1382        }
1383        let cut = fillet * (0.5 * theta).tan();
1384        legs[prev].trim_end = cut;
1385        legs[k].trim_start = cut;
1386        let tangent_point = corner - u * cut;
1387        let inward = (v - u * u.dot(v)).normalize();
1388        arcs[k] = Some(PathPiece::Arc {
1389            centre: tangent_point + inward * fillet,
1390            axis: normal.normalize(),
1391            start: tangent_point,
1392            angle: theta,
1393        });
1394    }
1395    let mut out = Vec::with_capacity(2 * m);
1396    for (k, leg) in legs.iter().enumerate() {
1397        if k > 0 {
1398            if let Some(arc) = arcs[k] {
1399                out.push(arc);
1400            }
1401        }
1402        let length = leg.piece.length();
1403        let needed = leg.trim_start + leg.trim_end;
1404        if needed == 0.0 {
1405            out.push(leg.piece);
1406            continue;
1407        }
1408        // A fillet using up exactly its segment meets the next one
1409        // tangentially; rounding in `tan` may overshoot by an ulp or two.
1410        let slack = 1e-9 * length;
1411        if needed > length + slack {
1412            return Err(GeomError::InvalidInput(format!(
1413                "swept disk fillet radius {fillet} does not fit: segment {k} is {length} \
1414                 long but its fillets need {needed} of it"
1415            )));
1416        }
1417        if needed < length - slack {
1418            out.push(
1419                leg.piece
1420                    .between(leg.trim_start / length, 1.0 - leg.trim_end / length),
1421            );
1422        }
1423    }
1424    if closed {
1425        if let Some(arc) = arcs[0] {
1426            out.push(arc);
1427        }
1428    }
1429    Ok(out)
1430}
1431
1432#[cfg(test)]
1433mod tests {
1434    use super::*;
1435
1436    #[test]
1437    fn a_fillet_is_tangent_to_both_segments() {
1438        let corner = Point3::new(1.0, 0.0, 0.0);
1439        let pieces = [
1440            PathPiece::Segment {
1441                start: Point3::ZERO,
1442                end: corner,
1443            },
1444            PathPiece::Segment {
1445                start: corner,
1446                end: Point3::new(1.0, 1.0, 1.0),
1447            },
1448        ];
1449        let path = fillet(&pieces, 0.2, 0.05, 1e-6, false).unwrap();
1450        assert_eq!(path.len(), 3);
1451        for pair in path.windows(2) {
1452            assert!((pair[0].end_point() - pair[1].start_point()).length() < 1e-12);
1453            assert!(turn(pair[0].end_tangent(), pair[1].start_tangent()) < 1e-12);
1454        }
1455    }
1456
1457    #[test]
1458    fn a_mitre_ring_lies_on_both_cylinders_and_on_the_bisector_plane() {
1459        let corner = Point3::new(1.0, 0.2, -0.3);
1460        let (u, v) = (
1461            Vec3::new(1.0, 0.0, 0.0),
1462            Vec3::new(0.2, 0.7, -0.4).normalize(),
1463        );
1464        let pieces = [
1465            PathPiece::Segment {
1466                start: corner - u,
1467                end: corner,
1468            },
1469            PathPiece::Segment {
1470                start: corner,
1471                end: corner + v,
1472            },
1473        ];
1474        let (r, inner) = (0.05, 0.03);
1475        let mesh = swept_disk_along_pieces(&pieces, r, Some(inner), None, 1e-3).unwrap();
1476        // One station per vertex: start, mitre, end.
1477        let per = mesh.positions.len() / 3;
1478        let n = (u + v).normalize();
1479        let off = |p: Point3, axis: Vec3| {
1480            let d = p - corner;
1481            (d - axis * d.dot(axis)).length()
1482        };
1483        for p in &mesh.positions[per..2 * per] {
1484            assert!((*p - corner).dot(n).abs() < 1e-12);
1485            let (a, b) = (off(*p, u), off(*p, v));
1486            assert!((a - b).abs() < 1e-12, "{a} {b}");
1487            assert!((a - r).abs() < 1e-12 || (a - inner).abs() < 1e-12, "{a}");
1488        }
1489    }
1490
1491    #[test]
1492    fn an_arc_walked_backwards_retraces_itself() {
1493        let arc = PathPiece::Arc {
1494            centre: Point3::new(0.0, 1.0, 0.0),
1495            axis: Vec3::Z,
1496            start: Point3::ZERO,
1497            angle: 1.0,
1498        };
1499        let back = arc.reversed();
1500        assert!((back.end_point() - arc.start_point()).length() < 1e-12);
1501        assert!((back.start_tangent() + arc.end_tangent()).length() < 1e-12);
1502        let part = arc.between(0.25, 0.75);
1503        assert!((part.length() - 0.5).abs() < 1e-12);
1504        assert!((part.start_point() - arc.point(0.25)).length() < 1e-12);
1505    }
1506}