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(§ions)
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, ¶meters)
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}