Skip to main content

brep_kernel/construction/sweep_topology/
sweep.rs

1use super::*;
2
3/// Sweep a CLOSED PLANAR profile loop along a path curve (§3.3/§5.7),
4/// TRANSPLANTING the profile onto the path — see [`SectionPlacement`] for the
5/// other reading, and for why the choice changes the solid.
6///
7/// Reuses the loft builder rather than a bespoke swept surface: the path is
8/// sampled at `STATIONS` uniform stations, a rotation-minimizing frame is
9/// propagated along it (double-reflection RMF, Wang et al. 2008 — the frame
10/// does NOT spin at inflections the way raw Frenet does), a rigidly
11/// transformed copy of the profile is placed at each station (identical
12/// degree/knots/weights, only the control points moved — which guarantees
13/// loft's per-section compatibility), and the stations are lofted through
14/// to produce the tube plus planar end caps.
15///
16/// Guards return a clear `Err` on an open/non-planar profile, a degenerate
17/// path tangent, or a loft failure. A self-intersecting result (path
18/// curvature radius smaller than the profile extent) is OUT OF SCOPE — the
19/// caller is responsible for keeping the tube from folding onto itself.
20pub fn sweep_profile_along_path(
21    profile: &[NurbsCurve],
22    path: &NurbsCurve,
23    name: Option<&str>,
24) -> Result<BrepSolid, String> {
25    // 32 stations is the original fixed sampling — golden parity pins the
26    // emitted geometry to it.  The helix variant raises the count with the
27    // turn count instead, hence the shared `_stations` core.
28    sweep_profile_along_path_stations(
29        profile,
30        path,
31        name,
32        32,
33        0.0,
34        None,
35        SectionPlacement::Transplant,
36    )
37}
38
39/// Twisted path sweep (§5.7): identical to [`sweep_profile_along_path`], but
40/// the profile additionally ROTATES about the path tangent, linearly in ARC
41/// LENGTH, from 0 at the sweep start to `twist_angle` radians (right-handed
42/// about the tangent) at the end.  The arc-length fraction comes from the
43/// sampled station polyline, not the raw path parameter, so a non-uniformly
44/// parameterized path still twists uniformly in space.
45///
46/// STATION LAW: 16 stations per quarter turn of twist, floored at the path
47/// sweep's 32 and capped at 1024 (the loft's dense interpolation solve is
48/// O(stations³) per control column — the same cap the helix uses).  At
49/// 16/quarter-turn the inter-station twist step is Δφ ≈ 5.6°, so the cubic
50/// v-interpolation error on a profile point circling at radius r is
51/// ≈ r·Δφ⁴/384 ≈ 2.4·10⁻¹⁰·r — far below any geometric tolerance.  The cap
52/// holds that density up to 16 full turns (|twist| = 32π = 1024/16 quarter
53/// turns); a larger twist would silently alias under the cap, so it is
54/// REJECTED with an honest error instead.  Everything else the path sweep
55/// documents (profile validity, self-intersection being the caller's
56/// responsibility) applies unchanged; twisting about the profile's own
57/// centroid adds no new radial extent, so no extra collision guard exists
58/// to compute here.
59pub fn sweep_profile_twisted(
60    profile: &[NurbsCurve],
61    path: &NurbsCurve,
62    twist_angle: f64,
63    name: Option<&str>,
64) -> Result<BrepSolid, String> {
65    use std::f64::consts::{FRAC_PI_2, TAU};
66
67    if !twist_angle.is_finite() {
68        return Err("sweep_profile_twisted: twist angle must be finite".into());
69    }
70    // 16 turns is where the 1024-station cap meets 16 stations/quarter-turn;
71    // beyond it the cap would degrade the twist sampling density silently.
72    const MAX_TURNS: f64 = 16.0;
73    if twist_angle.abs() > MAX_TURNS * TAU {
74        return Err(format!(
75            "sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
76             limit the 1024-station cap can resolve at 16 stations per quarter turn; \
77             split the sweep or reduce the twist",
78            twist_angle.abs() / TAU
79        ));
80    }
81    let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
82    let stations = (quarter_turns * 16).clamp(32, 1024);
83    sweep_profile_along_path_stations(
84        profile,
85        path,
86        name,
87        stations,
88        twist_angle,
89        None,
90        SectionPlacement::Transplant,
91    )
92    .map_err(|error| format!("sweep_profile_twisted: {error}"))
93}
94
95/// HOW a section is placed at each station — the two things "sweep a profile
96/// along a path" can mean, and the reason they are not the same shape.
97///
98/// `Transplant` MOVES the profile onto the path: the anchor origin lands on the
99/// station point and the profile plane becomes the station's NORMAL plane. The
100/// profile's own position and its angle to the path are both discarded — every
101/// sweep comes out centred on the path and square to it. That is the classic
102/// swept-surface reading, and it is what `SWP` builds.
103///
104/// `Rigid` CARRIES the profile: the section at station `k` is the profile moved
105/// by the same rigid motion the PATH undergoes between its start and station
106/// `k`,
107///
108/// ```text
109///     section_k(x) = P_k + R_k · (x − P_0)
110/// ```
111///
112/// where `R_k` is the rotation taking the start frame `(T_0, r_0, s_0)` to the
113/// station's `(T_k, r_k, s_k)`. The profile keeps its drawn POSITION and its
114/// drawn ANGLE to the path; what changes along the sweep is only what the path
115/// itself does. Two consequences are the whole point of the mode:
116///   - a STRAIGHT path makes every `R_k` the identity, so the section merely
117///     translates and the result is the oblique prism `extrude_profile_brep`
118///     builds from the same profile and vector — `pathAlign` and `translate`
119///     agree on a straight path rather than merely having the same volume;
120///   - a circular ARC makes `R_k` the rotation about the ARC'S OWN CENTRE AXIS
121///     (for a planar path the RMF's transport rotation IS that rotation), so
122///     the profile is carried round the pivot and the end cap arrives at the
123///     angle the start cap had to the tangent. A full-circle path degenerates
124///     to exactly a revolve of the profile about that axis.
125///
126/// `R_k = F_k · F_0ᵀ` does NOT depend on which perpendicular `r_0` the RMF was
127/// seeded with: re-seeding rolls `F_0` and every `F_k` by the same angle about
128/// the tangent, and the two rolls cancel in the product. That matters because
129/// `Transplant` maps `pu → r_k` directly, so under IT the profile's roll about
130/// the path is whatever `tangents[0].perpendicular()` happened to return — an
131/// arbitrary orientation no caller can predict. `Rigid` has no such freedom.
132#[derive(Debug, Clone, Copy, PartialEq, Eq)]
133pub enum SectionPlacement {
134    /// Move the profile onto the path, square to it (the classic sweep).
135    Transplant,
136    /// Carry the profile by the path's own rigid motion (see above).
137    Rigid,
138}
139
140/// The placement anchor a path sweep transplants its profile with: the plane
141/// frame `(origin = boundary-sample centroid, normal, pu, pv)` the station
142/// loop maps profile points through (`local = p − origin` → `station + ri·(local
143/// ·pu) + si·(local·pv)`). Extracted as data so a HOLE loop can sweep with its
144/// OUTER loop's anchor — sweeping each loop with its OWN centroid would
145/// re-center every loop onto the path and lose the hole's in-plane offset.
146#[derive(Debug, Clone, Copy)]
147pub struct ProfileAnchor {
148    pub origin: Vec3,
149    pub normal: Vec3,
150    pub pu: Vec3,
151    pub pv: Vec3,
152}
153
154/// [`sweep_profile_along_path`] with an explicit placement anchor (see
155/// [`ProfileAnchor`]) — the hole-loop cutter path: the swept loop is validated
156/// as usual but PLACED in its outer loop's frame.
157pub fn sweep_profile_along_path_anchored(
158    profile: &[NurbsCurve],
159    path: &NurbsCurve,
160    name: Option<&str>,
161    anchor: ProfileAnchor,
162) -> Result<BrepSolid, String> {
163    sweep_profile_along_path_stations(
164        profile,
165        path,
166        name,
167        32,
168        0.0,
169        Some(anchor),
170        SectionPlacement::Transplant,
171    )
172}
173
174/// [`sweep_profile_twisted`] with an explicit placement anchor: the hole loop
175/// twists about the SAME path axis as its outer loop (shared anchor), so the
176/// cutter stays registered with the outer wall through the whole twist.
177pub fn sweep_profile_twisted_anchored(
178    profile: &[NurbsCurve],
179    path: &NurbsCurve,
180    twist_angle: f64,
181    name: Option<&str>,
182    anchor: ProfileAnchor,
183) -> Result<BrepSolid, String> {
184    use std::f64::consts::{FRAC_PI_2, TAU};
185
186    if !twist_angle.is_finite() {
187        return Err("sweep_profile_twisted: twist angle must be finite".into());
188    }
189    const MAX_TURNS: f64 = 16.0;
190    if twist_angle.abs() > MAX_TURNS * TAU {
191        return Err(format!(
192            "sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
193             limit the 1024-station cap can resolve at 16 stations per quarter turn; \
194             split the sweep or reduce the twist",
195            twist_angle.abs() / TAU
196        ));
197    }
198    let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
199    let stations = (quarter_turns * 16).clamp(32, 1024);
200    sweep_profile_along_path_stations(
201        profile,
202        path,
203        name,
204        stations,
205        twist_angle,
206        Some(anchor),
207        SectionPlacement::Transplant,
208    )
209    .map_err(|error| format!("sweep_profile_twisted: {error}"))
210}
211
212/// Sweep a CLOSED PLANAR profile along a CHAIN of path curves joined
213/// head-to-tail — the multi-segment entry point, with the twist and anchor
214/// options the single-curve family exposes as separate functions folded into
215/// one signature (the feature layer picks all three per loop, so splitting them
216/// four ways here would only push the same match into the caller).
217///
218/// `segment_names` name the chain's segments for the joint refusal; a short or
219/// empty slice degrades to `<unnamed>` rather than failing.
220///
221/// `corner_advice` is appended to the cornered-joint refusal and is the CALLER's
222/// to write, because the builder has two callers with different answers: path
223/// sweep sends the user to SW (a different feature), while SW's own `pathAlign`
224/// sends them to its `translate` mode. A builder that named one feature would be
225/// telling half its users to switch to the feature they are already in.
226///
227/// A ONE-segment chain delegates to the single-curve builders unchanged, so the
228/// geometry every existing path sweep emits is untouched by this entry point
229/// existing — the multi-segment sampler is reached only by a path that actually
230/// has a joint.
231///
232/// Joints must be tangent-continuous; see [`MAX_JOINT_TANGENT_BREAK`] for why a
233/// corner is refused here instead of rounded off, and where such a path belongs.
234pub fn sweep_profile_along_chain(
235    profile: &[NurbsCurve],
236    chain: &[NurbsCurve],
237    segment_names: &[String],
238    twist_angle: f64,
239    name: Option<&str>,
240    anchor: Option<ProfileAnchor>,
241    placement_mode: SectionPlacement,
242    corner_advice: &str,
243) -> Result<BrepSolid, String> {
244    use std::f64::consts::{FRAC_PI_2, TAU};
245
246    if chain.is_empty() {
247        return Err("sweepSolid: path chain is empty".into());
248    }
249    // Single segment under `Transplant`: the existing builders, byte for byte.
250    // `Rigid` has no single-curve wrapper to delegate to and falls through to the
251    // chain sampler, which samples a 1-element chain over exactly the domain the
252    // single-curve sampler would (no joints, the whole budget on one segment).
253    if let ([single], SectionPlacement::Transplant) = (chain, placement_mode) {
254        return match (anchor, twist_angle == 0.0) {
255            (None, true) => sweep_profile_along_path(profile, single, name),
256            (Some(anchor), true) => {
257                sweep_profile_along_path_anchored(profile, single, name, anchor)
258            }
259            (None, false) => sweep_profile_twisted(profile, single, twist_angle, name),
260            (Some(anchor), false) => {
261                sweep_profile_twisted_anchored(profile, single, twist_angle, name, anchor)
262            }
263        };
264    }
265
266    // Everything else — any multi-segment chain, and a single segment under
267    // `Rigid` — takes the same station law the single-curve twisted builder
268    // uses, so density scales with twist identically on a chained path.
269    if !twist_angle.is_finite() {
270        return Err("sweep_profile_along_chain: twist angle must be finite".into());
271    }
272    const MAX_TURNS: f64 = 16.0;
273    if twist_angle.abs() > MAX_TURNS * TAU {
274        return Err(format!(
275            "sweep_profile_along_chain: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
276             limit the 1024-station cap can resolve at 16 stations per quarter turn; \
277             split the sweep or reduce the twist",
278            twist_angle.abs() / TAU
279        ));
280    }
281    let stations = if twist_angle == 0.0 {
282        32
283    } else {
284        let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
285        (quarter_turns * 16).clamp(32, 1024)
286    };
287    let _ = name;
288    sweep_profile_along_chain_stations(
289        profile,
290        chain,
291        segment_names,
292        stations,
293        twist_angle,
294        anchor,
295        placement_mode,
296        corner_advice,
297    )
298}
299
300/// Validate a closed planar profile loop and derive its placement anchor —
301/// the path sweep's §1 block, extracted bit-identically: 16 samples per curve,
302/// closure at `tolerance`, Newell normal, boundary-sample-mean origin,
303/// planarity at `tolerance * 100`, `pu = np.perpendicular()`, `pv = np × pu`.
304pub fn profile_anchor(profile: &[NurbsCurve]) -> Result<ProfileAnchor, String> {
305    let tolerance = 1e-6;
306    if profile.len() < 2 {
307        return Err("sweepSolid: profile needs at least 2 curves forming a closed loop".into());
308    }
309    let mut samples = Vec::new();
310    for (index, curve) in profile.iter().enumerate() {
311        let [start, end] = curve.domain()?;
312        let next = &profile[(index + 1) % profile.len()];
313        let next_start = next.domain()?[0];
314        if curve
315            .evaluate(end)?
316            .sub(next.evaluate(next_start)?)
317            .length()
318            > tolerance
319        {
320            return Err(format!(
321                "sweepSolid: profile is not closed at curve {index}"
322            ));
323        }
324        for sample in 0..16 {
325            samples.push(curve.evaluate(start + (end - start) * sample as f64 / 16.0)?);
326        }
327    }
328    let mut normal = Vec3::default();
329    let mut centroid = Vec3::default();
330    for index in 0..samples.len() {
331        let point = samples[index];
332        let next = samples[(index + 1) % samples.len()];
333        normal.x += (point.y - next.y) * (point.z + next.z);
334        normal.y += (point.z - next.z) * (point.x + next.x);
335        normal.z += (point.x - next.x) * (point.y + next.y);
336        centroid = centroid.add(point);
337    }
338    let np = normal
339        .normalized()
340        .map_err(|_| "sweepSolid: profile is degenerate (zero enclosed area)".to_string())?;
341    let origin = centroid.scale(1.0 / samples.len() as f64);
342    if samples
343        .iter()
344        .any(|point| point.sub(origin).dot(np).abs() > tolerance * 100.0)
345    {
346        return Err("sweepSolid: profile is not planar".into());
347    }
348    let pu = np.perpendicular()?;
349    let pv = np.cross(pu).normalized()?;
350    Ok(ProfileAnchor {
351        origin,
352        normal: np,
353        pu,
354        pv,
355    })
356}
357
358/// Core of the path sweep with an explicit station count.  Every consumer of
359/// the station count (path sampling, RMF propagation, section placement)
360/// derives from the one `stations` argument so the density scales as a unit.
361/// `twist_angle` (radians; 0 for the untwisted variants) rotates the placed
362/// profile about the path tangent linearly in sampled arc length — the
363/// `twist_angle == 0.0` fast path leaves the RMF axes bit-identical, so the
364/// untwisted callers keep golden parity.
365fn sweep_profile_along_path_stations(
366    profile: &[NurbsCurve],
367    path: &NurbsCurve,
368    name: Option<&str>,
369    stations: usize,
370    twist_angle: f64,
371    anchor: Option<ProfileAnchor>,
372    placement_mode: SectionPlacement,
373) -> Result<BrepSolid, String> {
374    // Loft carries no face names; the app stamps them onto the emitted face
375    // order.  Accept `name` for ABI symmetry with the other builders.
376    let _ = name;
377    let tolerance = 1e-6;
378    if stations < 2 {
379        return Err("sweepSolid: need at least 2 stations".into());
380    }
381
382    // --- 1. Validate the profile: closed + planar; derive (origin, np, pu, pv).
383    //        A caller-supplied anchor OVERRIDES the placement frame (the swept
384    //        loop is still validated against its own plane), so a hole loop
385    //        rides the path in its outer loop's frame instead of re-centering.
386    //        Resolved HERE, before the path is touched, so a bad profile is
387    //        still reported ahead of a bad path exactly as it always was.
388    let placement = resolve_placement(profile, anchor, placement_mode)?;
389
390    // --- 2. Sample the path; require a non-degenerate tangent at every station.
391    let [t0, t1] = path.domain()?;
392    if (t1 - t0).abs() <= tolerance {
393        return Err("sweepSolid: path domain is degenerate".into());
394    }
395    let mut points = Vec::with_capacity(stations);
396    let mut tangents = Vec::with_capacity(stations);
397    for index in 0..stations {
398        let t = t0 + (t1 - t0) * index as f64 / (stations - 1) as f64;
399        let derivatives = path.derivatives(t, 1)?;
400        let tangent = derivatives[1]
401            .normalized()
402            .map_err(|_| format!("sweepSolid: path tangent is degenerate at station {index}"))?;
403        points.push(derivatives[0]);
404        tangents.push(tangent);
405    }
406
407    sweep_sections_through_samples(
408        profile,
409        &points,
410        &tangents,
411        twist_angle,
412        placement,
413        placement_mode,
414    )
415}
416
417/// The placement frame a sweep transplants its profile with: the caller's
418/// `anchor` when it supplied one (a hole loop riding its outer loop's frame),
419/// otherwise the profile's own. Either way the profile is VALIDATED — closed and
420/// planar — because `profile_anchor` is what does that validating, and skipping
421/// it for an anchored loop would let an open hole loop through.
422///
423/// Called by each SAMPLER before it touches the path, so the "bad profile" error
424/// still precedes the "bad path" one, as it did when the two lived in one
425/// function.
426fn resolve_placement(
427    profile: &[NurbsCurve],
428    anchor: Option<ProfileAnchor>,
429    placement_mode: SectionPlacement,
430) -> Result<ProfileAnchor, String> {
431    let computed = profile_anchor(profile)?;
432    Ok(match placement_mode {
433        // RIGID never transplants, so there is no frame to borrow: every loop
434        // keeps its own position already, which is the whole reason the anchor
435        // exists under `Transplant`. Taking a caller's anchor here would silently
436        // swap in another loop's NORMAL for the fold-back guard below.
437        SectionPlacement::Rigid => computed,
438        SectionPlacement::Transplant => anchor.unwrap_or(computed),
439    })
440}
441
442/// The largest tangent break, in radians, a chained path may have at a joint.
443///
444/// This is a REAL geometric boundary, not a tuning knob. The builder skins
445/// consecutive stations by lofting through them, so a corner between two
446/// stations is rendered as a smooth blend across the corner — the tube cuts it,
447/// and near a sharp one the section sweeps through itself. A path sweep can
448/// honestly build a G1 chain (line→fillet→line, spline pieces, an arc train)
449/// and cannot honestly build a cornered polyline. `Sweep` (SW) is the feature
450/// for a cornered path: it builds one prism per segment and unions them, which
451/// is exactly the construction that gives a corner a real mitre.
452///
453/// ~1.15° (0.02 rad). Loose enough to absorb the tangent disagreement of two
454/// curves a sketch chained on endpoint coincidence, tight enough that anything
455/// a user would call a corner is refused rather than silently rounded off.
456///
457/// This is INDEPENDENT of the chainer's join tolerance, and deliberately so.
458/// `common::chain_path_segments` joins on POSITION (`1e-5 * scale`); this gate
459/// measures DIRECTION. Two collinear segments meeting with a positional gap have
460/// a zero tangent break, and two segments meeting exactly at a point can still
461/// break 90°. So a run the chainer happily orders head-to-tail may still be
462/// refused here — which is not an inconsistency but the SW/SWP split itself: the
463/// chainer's job is to decide what order the picks form a run in, and this gate's
464/// job is to decide whether that run is one this builder can skin.
465const MAX_JOINT_TANGENT_BREAK: f64 = 0.02;
466
467/// Sweep a CLOSED PLANAR profile along a CHAIN of path curves, joined
468/// head-to-tail, as one continuous tube.
469///
470/// The chain is sampled as a single trajectory and handed to the shared core,
471/// so ONE rotation-minimizing frame is propagated across the whole path: a joint
472/// is just another pair of adjacent stations to the RMF, which is what makes the
473/// profile arrive at the far end in the orientation the near end implies rather
474/// than snapping at each segment.
475///
476/// STATION BUDGET. `stations` are distributed across the segments in proportion
477/// to estimated ARC LENGTH (chord sums over 16 samples per curve), floored at 2
478/// per segment so a short fillet between two long lines still carries its
479/// curvature. The shared endpoint of two adjacent segments is sampled ONCE, not
480/// twice — a duplicated station would land in the RMF's coincident-step branch
481/// and waste a section on zero advance.
482///
483/// JOINTS must be tangent-continuous to within [`MAX_JOINT_TANGENT_BREAK`];
484/// a break beyond it is refused, naming the two segments and the angle.
485fn sweep_profile_along_chain_stations(
486    profile: &[NurbsCurve],
487    chain: &[NurbsCurve],
488    segment_names: &[String],
489    stations: usize,
490    twist_angle: f64,
491    anchor: Option<ProfileAnchor>,
492    placement_mode: SectionPlacement,
493    corner_advice: &str,
494) -> Result<BrepSolid, String> {
495    let tolerance = 1e-6;
496    if chain.is_empty() {
497        return Err("sweepSolid: path chain is empty".into());
498    }
499    if stations < 2 {
500        return Err("sweepSolid: need at least 2 stations".into());
501    }
502    // Profile first, matching the single-curve sampler's order of complaint.
503    let placement = resolve_placement(profile, anchor, placement_mode)?;
504
505    // --- Per-segment arc-length estimate, and the domain each one is sampled
506    //     over. A degenerate domain is refused here rather than producing a
507    //     silently collapsed station block.
508    let mut domains = Vec::with_capacity(chain.len());
509    let mut lengths = Vec::with_capacity(chain.len());
510    for (index, curve) in chain.iter().enumerate() {
511        let [t0, t1] = curve.domain()?;
512        if (t1 - t0).abs() <= tolerance {
513            return Err(format!(
514                "sweepSolid: path segment {index} has a degenerate domain"
515            ));
516        }
517        let mut length = 0.0;
518        let mut previous = curve.evaluate(t0)?;
519        for sample in 1..=16 {
520            let point = curve.evaluate(t0 + (t1 - t0) * sample as f64 / 16.0)?;
521            length += point.sub(previous).length();
522            previous = point;
523        }
524        domains.push([t0, t1]);
525        lengths.push(length);
526    }
527    let total_length: f64 = lengths.iter().sum();
528    if total_length <= tolerance {
529        return Err("sweepSolid: path chain has zero length".into());
530    }
531
532    // --- Station budget per segment: proportional to arc length, floored at 2
533    //     (a segment needs its two ends), and the remainder handed to the
534    //     longest segment so the totals land exactly on `stations`.
535    let mut budget: Vec<usize> = lengths
536        .iter()
537        .map(|length| {
538            let share = (stations as f64 * length / total_length).round() as usize;
539            share.max(2)
540        })
541        .collect();
542    // Each interior joint shares one station between its two segments, so the
543    // emitted total is sum(budget) - (segments - 1).
544    let joints = chain.len() - 1;
545    let emitted: usize = budget.iter().sum::<usize>() - joints;
546    // The 2-per-segment floor means a chain of n segments emits at least n+1
547    // stations, so a long enough chain cannot be held under the loft's station
548    // cap by any distribution. The loft's interpolation solve is O(stations³)
549    // per control column, so silently accepting such a path would look like a
550    // hang. Refuse with the count instead — the same honesty the twist cap uses.
551    const MAX_STATIONS: usize = 1024;
552    if emitted > MAX_STATIONS {
553        return Err(format!(
554            "sweepSolid: a {}-segment path needs at least {emitted} stations, past the \
555             {MAX_STATIONS}-station cap the loft solve can carry; sweep it in fewer, \
556             longer pieces",
557            chain.len()
558        ));
559    }
560    if emitted < stations {
561        let longest = lengths
562            .iter()
563            .enumerate()
564            .max_by(|a, b| a.1.total_cmp(b.1))
565            .map(|(index, _)| index)
566            .unwrap_or(0);
567        budget[longest] += stations - emitted;
568    }
569
570    // --- Sample the chain as ONE trajectory, checking tangent continuity at
571    //     every joint as it is crossed.
572    let mut points: Vec<Vec3> = Vec::new();
573    let mut tangents: Vec<Vec3> = Vec::new();
574    for (index, curve) in chain.iter().enumerate() {
575        let [t0, t1] = domains[index];
576        let count = budget[index];
577
578        // The JOINT check runs at this segment's START parameter, which is the
579        // point the previous segment ended on — comparing the tangent the path
580        // arrives with against the one it departs with. Sampling the next
581        // station instead would fold the segment's own curvature into the
582        // measured break and let a genuine corner through on a curved segment.
583        if index > 0 {
584            let departing = curve.derivatives(t0, 1)?[1].normalized().map_err(|_| {
585                format!("sweepSolid: path tangent is degenerate at the start of segment {index}")
586            })?;
587            let arriving = *tangents
588                .last()
589                .expect("a previous segment emitted stations");
590            // Both unit, so the dot is the cosine of the break.
591            let break_angle = arriving.dot(departing).clamp(-1.0, 1.0).acos();
592            if break_angle > MAX_JOINT_TANGENT_BREAK {
593                let previous_name = segment_names
594                    .get(index - 1)
595                    .map(String::as_str)
596                    .unwrap_or("<unnamed>");
597                let name = segment_names
598                    .get(index)
599                    .map(String::as_str)
600                    .unwrap_or("<unnamed>");
601                return Err(format!(
602                    "the path must be tangent-continuous: segments '{previous_name}' and \
603                     '{name}' meet at a {:.1}° corner. The profile is skinned between sampled \
604                     stations along the path, so a corner would be ROUNDED OFF rather than \
605                     mitred. {corner_advice}",
606                    break_angle.to_degrees()
607                ));
608            }
609        }
610
611        // Skip the first sample of every segment after the first: that station
612        // is the joint, already emitted as the previous segment's last.
613        let first = usize::from(index > 0);
614        for sample in first..count {
615            let t = t0 + (t1 - t0) * sample as f64 / (count - 1) as f64;
616            let derivatives = curve.derivatives(t, 1)?;
617            let tangent = derivatives[1].normalized().map_err(|_| {
618                format!("sweepSolid: path tangent is degenerate on segment {index}")
619            })?;
620            points.push(derivatives[0]);
621            tangents.push(tangent);
622        }
623    }
624
625    sweep_sections_through_samples(
626        profile,
627        &points,
628        &tangents,
629        twist_angle,
630        placement,
631        placement_mode,
632    )
633}
634
635/// The part of the sweep that does not care HOW the stations were sampled:
636/// validate the profile, distribute the twist over the sampled polyline,
637/// propagate the rotation-minimizing frame, place a profile copy per station
638/// and loft through them.
639///
640/// Split out so a MULTI-SEGMENT path can share it
641/// ([`sweep_profile_along_chain_stations`]). The single-curve sampler above is
642/// deliberately left calling this with its OWN sampling rather than being
643/// re-expressed as a one-element chain: two code paths that agree today are not
644/// the same thing as one code path, and the untwisted single-curve geometry is
645/// pinned by volume tests that must not move.
646///
647/// The RMF is propagated over the station sequence as a whole, so a chain's
648/// frame carries across a joint exactly as it carries across any other pair of
649/// adjacent stations — there is no per-segment restart to reconcile.
650fn sweep_sections_through_samples(
651    profile: &[NurbsCurve],
652    points: &[Vec3],
653    tangents: &[Vec3],
654    twist_angle: f64,
655    placement: ProfileAnchor,
656    placement_mode: SectionPlacement,
657) -> Result<BrepSolid, String> {
658    let tolerance = 1e-6;
659    let stations = points.len();
660    if stations < 2 {
661        return Err("sweepSolid: need at least 2 stations".into());
662    }
663
664    let ProfileAnchor { origin, pu, pv, .. } = placement;
665
666    // --- 2b. Twist distribution: cumulative ARC-LENGTH fractions over the
667    //         sampled station polyline (chord sums), so the twist advances
668    //         uniformly in space even on a non-uniformly parameterized path.
669    //         Only computed when a twist is actually requested — the
670    //         `twist_angle == 0.0` path must stay bit-identical to the
671    //         pre-twist builder.
672    let twist_fractions: Option<Vec<f64>> = if twist_angle != 0.0 {
673        let mut cumulative = vec![0.0; stations];
674        let mut total = 0.0;
675        for index in 1..stations {
676            total += points[index].sub(points[index - 1]).length();
677            cumulative[index] = total;
678        }
679        if total <= tolerance {
680            return Err("sweepSolid: path has zero length; cannot distribute the twist".into());
681        }
682        for length in &mut cumulative {
683            *length /= total;
684        }
685        Some(cumulative)
686    } else {
687        None
688    };
689
690    // --- 3. Rotation-minimizing frames via the double-reflection method.
691    let mut r_axes = Vec::with_capacity(stations);
692    let mut s_axes = Vec::with_capacity(stations);
693    let r0 = tangents[0].perpendicular()?; // any unit vector ⟂ T0
694    s_axes.push(tangents[0].cross(r0).normalized()?);
695    r_axes.push(r0);
696    for index in 0..stations - 1 {
697        let t_next = tangents[index + 1];
698        let v1 = points[index + 1].sub(points[index]);
699        let c1 = v1.dot(v1);
700        let r_candidate = if c1 <= 1e-18 {
701            // Coincident stations: carry the reference axis forward unchanged.
702            r_axes[index]
703        } else {
704            // First reflection across the plane bisecting the step vector.
705            let reflected_r = r_axes[index].sub(v1.scale(2.0 / c1 * v1.dot(r_axes[index])));
706            let reflected_t = tangents[index].sub(v1.scale(2.0 / c1 * v1.dot(tangents[index])));
707            // Second reflection across the plane bisecting the tangents.
708            let v2 = t_next.sub(reflected_t);
709            let c2 = v2.dot(v2);
710            if c2 <= 1e-18 {
711                reflected_r
712            } else {
713                reflected_r.sub(v2.scale(2.0 / c2 * v2.dot(reflected_r)))
714            }
715        };
716        // Re-orthogonalize against the new tangent to shed floating drift.
717        let r_next = r_candidate
718            .sub(t_next.scale(r_candidate.dot(t_next)))
719            .normalized()
720            .map_err(|_| format!("sweepSolid: frame degenerated at station {index}"))?;
721        s_axes.push(t_next.cross(r_next).normalized()?);
722        r_axes.push(r_next);
723    }
724
725    // --- 3b. A CLOSED trajectory cannot be capped: the first and last sections
726    //         land on top of each other, and the loft's own complaint about that
727    //         ("end sections coincide") says nothing about the path. Name it here
728    //         instead, with the two things that do work. The check is on the
729    //         PATH's first and last station POINTS, which is enough: a path that
730    //         comes back to where it started puts the two caps in the same place
731    //         whatever direction it arrives from — coincident when it also
732    //         arrives pointing the same way, interpenetrating when it does not,
733    //         and neither is a solid worth building.
734    //
735    //         The scale is the path's own travel, so this is a proportion, not a
736    //         length: ends within a thousandth of the distance travelled are the
737    //         same place for capping purposes. On a circular arc that band is the
738    //         last ~0.36° before closure, so a 359° sweep still builds.
739    {
740        let travel: f64 = points
741            .windows(2)
742            .map(|pair| pair[1].sub(pair[0]).length())
743            .sum();
744        let closure = points[stations - 1].sub(points[0]).length();
745        if travel > tolerance && closure <= 1e-3 * travel {
746            return Err(format!(
747                "sweepSolid: the path returns to where it started ({closure:.3e} apart after \
748                 travelling {travel:.3e}), so the sweep's two end caps would land on top of each \
749                 other. A closed path swept this way IS a revolution — build it with Revolve \
750                 about the same axis, or sweep the run in two halves and union them"
751            ));
752        }
753    }
754
755    // --- 4. Place a transformed copy of the profile at each station, by the mode
756    //        the caller asked for (see `SectionPlacement`).
757    //
758    //   TRANSPLANT maps the profile's own plane frame onto the station frame,
759    //       world = P_k + r_k·(l·pu) + s_k·(l·pv),   l = x − origin
760    //   which lands the anchor origin ON the path and the profile square to it.
761    //   The profile's distance from the path and its angle to the path are both
762    //   discarded — there is no term along the tangent to carry them.
763    //
764    //   RIGID applies the path's own motion to the profile where it stands,
765    //       world = P_k + a·T_k + b·r_k + c·s_k,     (a,b,c) = (l·T_0, l·r_0, l·s_0)
766    //                                                l = x − P_0
767    //   which is `P_k + R_k·(x − P_0)` for `R_k = F_k·F_0ᵀ`, written in the frame
768    //   basis rather than as a matrix. The tangent term `a·T_k` is exactly what
769    //   TRANSPLANT drops, and carrying it is what preserves both the profile's
770    //   offset from the path and its angle to it. At station 0 the frames
771    //   coincide and the section IS the profile, so the START cap lies in the
772    //   plane it was drawn in.
773    let path_start = points[0];
774    let (t0, r0, s0) = (tangents[0], r_axes[0], s_axes[0]);
775    let rigid = placement_mode == SectionPlacement::Rigid;
776    let mut sections: Vec<Vec<NurbsCurve>> = Vec::with_capacity(stations);
777    // Only the RIGID guard below reads these; TRANSPLANT leaves them empty.
778    let mut station_points: Vec<Vec<Vec3>> = Vec::with_capacity(if rigid { stations } else { 0 });
779    let mut section_normals: Vec<Vec3> = Vec::with_capacity(if rigid { stations } else { 0 });
780    let mut anchor_track: Vec<Vec3> = Vec::with_capacity(if rigid { stations } else { 0 });
781    for station in 0..stations {
782        let station_origin = points[station];
783        // Rotate the RMF axes about the tangent by the station's twist angle
784        // (Rodrigues on vectors ⟂ the tangent: r' = r·cosφ + s·sinφ,
785        // s' = s·cosφ − r·sinφ, since s = t × r and t × s = −r).
786        let (ri, si) = match &twist_fractions {
787            Some(fractions) => {
788                let phi = twist_angle * fractions[station];
789                let (sin_phi, cos_phi) = phi.sin_cos();
790                let r = r_axes[station];
791                let s = s_axes[station];
792                (
793                    r.scale(cos_phi).add(s.scale(sin_phi)),
794                    s.scale(cos_phi).sub(r.scale(sin_phi)),
795                )
796            }
797            None => (r_axes[station], s_axes[station]),
798        };
799        let tk = tangents[station];
800        let place = |x: Vec3| -> Vec3 {
801            if rigid {
802                let local = x.sub(path_start);
803                station_origin
804                    .add(tk.scale(local.dot(t0)))
805                    .add(ri.scale(local.dot(r0)))
806                    .add(si.scale(local.dot(s0)))
807            } else {
808                let local = x.sub(origin);
809                station_origin
810                    .add(ri.scale(local.dot(pu)))
811                    .add(si.scale(local.dot(pv)))
812            }
813        };
814        if rigid {
815            // The section's own normal takes the same rotation (a direction, so
816            // `R_k` without the translation), and the anchor origin's track is
817            // the profile's bulk motion for the near-parallel guard.
818            let n0 = placement.normal;
819            section_normals.push(
820                tk.scale(n0.dot(t0))
821                    .add(ri.scale(n0.dot(r0)))
822                    .add(si.scale(n0.dot(s0))),
823            );
824            anchor_track.push(place(origin));
825            station_points.push(Vec::with_capacity(
826                profile.iter().map(|c| c.control_points.len()).sum(),
827            ));
828        }
829        let mut section = Vec::with_capacity(profile.len());
830        for curve in profile {
831            let control_points = curve
832                .control_points
833                .iter()
834                .map(|point| {
835                    let weight = point.w;
836                    let euclidean = Vec3::new(point.x / weight, point.y / weight, point.z / weight);
837                    let world = place(euclidean);
838                    if rigid {
839                        station_points[station].push(world);
840                    }
841                    Vec4 {
842                        x: world.x * weight,
843                        y: world.y * weight,
844                        z: world.z * weight,
845                        w: weight,
846                    }
847                })
848                .collect();
849            section.push(NurbsCurve::new(
850                curve.degree,
851                curve.knots.clone(),
852                control_points,
853            )?);
854        }
855        sections.push(section);
856    }
857
858    // --- 4b. RIGID only: refuse a sweep whose sections fold through each other.
859    //
860    // TRANSPLANT keeps every section square to the path, so it can only fold
861    // where the path's curvature radius drops below the profile's extent — the
862    // case this builder has always documented as the caller's to avoid. RIGID
863    // can fold for a second, much more reachable reason: a turning path carries
864    // the profile around the TURN'S OWN AXIS, so a profile straddling that axis
865    // has one half advancing while the other retreats. That is not a tolerance
866    // question, it is a sign question, so it is answered exactly.
867    //
868    // The placement map is AFFINE in `x`, so the advance is affine in `x` too and
869    // its extremes over the section lie on the control points. Checking the hull
870    // therefore bounds the advance over the curves AND over the region they
871    // enclose — this is a proof, not a sampling.
872    if rigid {
873        let mut most_positive = 0.0f64;
874        let mut most_negative = 0.0f64;
875        let mut longest_step = 0.0f64;
876        for station in 0..stations - 1 {
877            let normal = section_normals[station];
878            for (before, after) in station_points[station]
879                .iter()
880                .zip(&station_points[station + 1])
881            {
882                let step = after.sub(*before);
883                longest_step = longest_step.max(step.length());
884                let advance = step.dot(normal);
885                most_positive = most_positive.max(advance);
886                most_negative = most_negative.min(advance);
887            }
888        }
889        // Scale-relative, and far below any geometric tolerance: this band only
890        // has to separate a sign from floating noise, never a small motion from
891        // a large one.
892        let noise = 1e-9 * longest_step.max(tolerance);
893        if most_positive > noise && most_negative < -noise {
894            return Err(
895                "sweepSolid: the profile sweeps back through itself — part of it advances along \
896                 the path while part of it retreats. A turning path carries the profile around \
897                 the turn's own axis, so a profile that STRADDLES that axis folds into itself; \
898                 move the profile clear of the axis, or sweep the run in pieces"
899                    .into(),
900            );
901        }
902        if most_positive <= noise && most_negative >= -noise {
903            return Err(
904                "sweepSolid: the path does not advance through the profile — it runs inside the \
905                 profile plane, so the sweep encloses no volume"
906                    .into(),
907            );
908        }
909        // The near-parallel case, on the anchor origin's own track: the same 0.1
910        // threshold `extrude_profile_brep` refuses a sliver at, so a STRAIGHT
911        // path is accepted here exactly when the translational sweep of the same
912        // profile and vector is. A station whose origin barely moves is skipped —
913        // an origin sitting on the turn axis means the profile straddles it, and
914        // the fold-back guard above has already spoken.
915        for station in 0..stations - 1 {
916            let step = anchor_track[station + 1].sub(anchor_track[station]);
917            let length = step.length();
918            if length <= noise {
919                continue;
920            }
921            if step.dot(section_normals[station]).abs() / length < 0.1 {
922                return Err(
923                    "sweepSolid: the path is nearly parallel to the profile plane, which sweeps a \
924                     sliver rather than a solid"
925                        .into(),
926                );
927            }
928        }
929    }
930
931    // --- 5. Loft through the swept sections (side walls + planar end caps).
932    loft_profile_brep(&sections)
933        .map_err(|error| format!("sweepSolid: loft through swept sections failed: {error}"))
934}
935
936/// Helical sweep (§5.7): sweep a CLOSED PLANAR profile loop along a helix of
937/// `helix_radius` about the axis through `axis_origin` with direction
938/// `axis_direction`, rising `pitch` per revolution for `turns` revolutions.
939///
940/// The helix is transcendental, not exactly NURBS-representable, so the path
941/// is FIT: points (R·cosθ, R·sinθ, pitch·θ/2π) in the axis frame are sampled
942/// uniformly in θ and globally interpolated with a cubic (`interpolate_curve`
943/// — the same machinery every other fitted path here uses).  DENSITY: 64
944/// samples per turn, capped at 1025 total nodes.  Cubic interpolation error
945/// on a circle of radius R with node spacing Δθ is ≈ R·Δθ⁴/384: ~2·10⁻⁷·R at
946/// 64/turn, and still ~6·10⁻⁵·R at the cap's worst case (16/turn at the
947/// 64-turn limit) — orders below any profile a caller could sweep without
948/// self-intersecting.  Uniform-in-θ parameters are exact chord-length for a
949/// helix (constant speed), which is what the averaged-knot interpolation
950/// assumes.
951///
952/// The fitted path then drives the EXISTING path-sweep core with 32 stations
953/// per turn (min 32, capped at 1024 — the loft's dense interpolation solve is
954/// O(stations³) per control column, so the cap trades per-turn density, never
955/// correctness, at high turn counts).
956///
957/// GUARDS beyond the path sweep's own: the path sweep documents
958/// self-intersection as caller responsibility, so the helix variant — which
959/// knows its curvature analytically — rejects the two garbage modes itself:
960///   • fold-over: the helix curvature radius (R² + c²)/R (c = pitch/2π) must
961///     exceed the profile's max extent about its centroid, or the tube folds
962///     through itself on the inner side (the torus tube-radius > major-radius
963///     failure, pitch-relaxed);
964///   • coil collision (turns ≥ 1): the normal gap between consecutive coils,
965///     pitch·2πR/√((2πR)² + pitch²), must exceed the profile diameter.
966/// Both use the profile's max sample distance from its centroid — conservative
967/// for asymmetric profiles (extent in a harmless direction still counts), but
968/// a false reject beats silent garbage.
969///
970/// KNOWN LIMITATION: whole-turn shallow helixes (pitch ≲ 0.63·R) are rejected
971/// by the loft's cap-plane guard (`|n·axis| ≥ 0.1` — the end-to-end axis is
972/// purely axial while the cap normal is nearly tangential); the error
973/// propagates honestly rather than being worked around.
974pub fn sweep_profile_helix(
975    profile: &[NurbsCurve],
976    axis_origin: Vec3,
977    axis_direction: Vec3,
978    helix_radius: f64,
979    pitch: f64,
980    turns: f64,
981    name: Option<&str>,
982) -> Result<BrepSolid, String> {
983    use std::f64::consts::TAU;
984
985    // --- 1. Validate the helix parameters with honest errors.
986    let w = axis_direction
987        .normalized()
988        .map_err(|_| "sweep_profile_helix: axis direction is degenerate".to_string())?;
989    if !(helix_radius.is_finite() && helix_radius > 0.0) {
990        return Err("sweep_profile_helix: helix radius must be positive".into());
991    }
992    if !(pitch.is_finite() && pitch > 0.0) {
993        return Err("sweep_profile_helix: pitch must be positive".into());
994    }
995    if !(turns.is_finite() && turns > 0.0) {
996        return Err("sweep_profile_helix: turns must be positive".into());
997    }
998    // 64 turns bounds the loft's O(stations³) interpolation solve; beyond it
999    // the station cap would silently degrade per-turn density anyway.
1000    if turns > 64.0 {
1001        return Err("sweep_profile_helix: turns must be at most 64".into());
1002    }
1003
1004    // --- 2. Profile extent about its centroid, sampled exactly like the path
1005    //        sweep derives its placement origin (boundary-sample mean), so the
1006    //        extent is measured about the point that actually rides the path.
1007    let mut samples = Vec::new();
1008    for curve in profile {
1009        let [start, end] = curve.domain()?;
1010        for sample in 0..16 {
1011            samples.push(curve.evaluate(start + (end - start) * sample as f64 / 16.0)?);
1012        }
1013    }
1014    if !samples.is_empty() {
1015        let mut centroid = Vec3::default();
1016        for point in &samples {
1017            centroid = centroid.add(*point);
1018        }
1019        let origin = centroid.scale(1.0 / samples.len() as f64);
1020        let extent = samples
1021            .iter()
1022            .map(|point| point.sub(origin).length())
1023            .fold(0.0, f64::max);
1024        let c = pitch / TAU; // axial rise per radian
1025                             // Fold-over: profile reaches past the helix's center of curvature.
1026        let curvature_radius = (helix_radius * helix_radius + c * c) / helix_radius;
1027        if extent >= curvature_radius {
1028            return Err(format!(
1029                "sweep_profile_helix: profile extent {extent:.6} reaches the helix \
1030                 curvature radius {curvature_radius:.6}; the tube would fold through \
1031                 itself — increase the helix radius or pitch, or shrink the profile"
1032            ));
1033        }
1034        // Coil collision: only possible once the sweep spans a full revolution.
1035        if turns >= 1.0 {
1036            let circumference = TAU * helix_radius;
1037            let turn_length = (circumference * circumference + pitch * pitch).sqrt();
1038            let coil_gap = pitch * circumference / turn_length;
1039            if coil_gap <= 2.0 * extent {
1040                return Err(format!(
1041                    "sweep_profile_helix: consecutive turns would self-intersect — \
1042                     coil gap {coil_gap:.6} does not clear the profile diameter {:.6}; \
1043                     increase the pitch or shrink the profile",
1044                    2.0 * extent
1045                ));
1046            }
1047        }
1048    }
1049
1050    // --- 3. Fit the helical path (see the density rationale in the fn docs).
1051    let u = w.perpendicular()?;
1052    let v = w.cross(u).normalized()?;
1053    let total_angle = turns * TAU;
1054    let rise = pitch / TAU;
1055    let count = ((turns * 64.0).ceil() as usize + 1).clamp(9, 1025);
1056    let mut points = Vec::with_capacity(count);
1057    let mut parameters = Vec::with_capacity(count);
1058    for index in 0..count {
1059        let s = index as f64 / (count - 1) as f64;
1060        let theta = total_angle * s;
1061        points.push(
1062            axis_origin
1063                .add(u.scale(helix_radius * theta.cos()))
1064                .add(v.scale(helix_radius * theta.sin()))
1065                .add(w.scale(rise * theta)),
1066        );
1067        parameters.push(s);
1068    }
1069    let path = interpolate_curve(&points, 3, &parameters)
1070        .map_err(|error| format!("sweep_profile_helix: helix path fit failed: {error}"))?;
1071
1072    // --- 4. Drive the existing sweep core; its (or the loft's) failures
1073    //        propagate with helix context prepended.
1074    let stations = ((turns * 32.0).ceil() as usize).clamp(32, 1024);
1075    // The helix builder places the profile ON the fitted helix (the coil's
1076    // section is drawn once and carried round), which is `Transplant` by
1077    // construction — the radius comes from the helix argument, not the profile's
1078    // own position.
1079    sweep_profile_along_path_stations(
1080        profile,
1081        &path,
1082        name,
1083        stations,
1084        0.0,
1085        None,
1086        SectionPlacement::Transplant,
1087    )
1088        .map_err(|error| format!("sweep_profile_helix: {error}"))
1089}