1use crate::{KernelRefusal, KernelStage, OrRefuse};
2use super::*;
3
4pub(super) struct RimPiece {
6 pub(super) face_id: u64,
7 pub(super) loop_index: usize,
8 pub(super) window: [f64; 2],
10 pub(super) edge_id: u64,
11}
12
13pub fn blend_smooth_chain(
19 solid: &BrepSolid,
20 seed_edge_id: u64,
21 radius: f64,
22 chamfer: bool,
23 name: Option<&str>,
24) -> Result<BrepSolid, KernelRefusal> {
25 let _operation = crate::blend::BlendOperation::enter();
28 if !(radius > 0.0) || !radius.is_finite() {
29 return Err(KernelRefusal::input(KernelStage::Collect, "radius", "blend: radius must be positive"));
30 }
31 let chain = collect_smooth_chain(solid, seed_edge_id)?;
32 if chain.segments.len() < 2 {
33 return Err(KernelRefusal::internal(KernelStage::Refine, "chain_single_segment", "blend: chain collapsed to a single segment"));
34 }
35 if chain.closed {
36 blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
37 } else {
38 blend_open_smooth_chain(solid, &chain, radius, chamfer, name)
39 }
40}
41
42pub fn blend_smooth_chain_if_closed(
49 solid: &BrepSolid,
50 seed_edge_id: u64,
51 radius: f64,
52 chamfer: bool,
53 name: Option<&str>,
54) -> Result<BrepSolid, KernelRefusal> {
55 let _operation = crate::blend::BlendOperation::enter();
58 if !(radius > 0.0) || !radius.is_finite() {
59 return Err(KernelRefusal::input(KernelStage::Collect, "radius", "blend: radius must be positive"));
60 }
61 let chain = collect_smooth_chain(solid, seed_edge_id)?;
62 if chain.segments.len() < 2 {
63 return Err(KernelRefusal::internal(KernelStage::Refine, "chain_single_segment", "blend: chain collapsed to a single segment"));
64 }
65 if !chain.closed {
66 return Err(
67 KernelRefusal::unsupported(KernelStage::Classify, "open_tangent_chain", "blend: the edge is one arc of an OPEN tangent chain; the unselected \
68 continuation is capped, not blended"),
69 );
70 }
71 blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
72}
73
74fn blend_closed_smooth_chain(
85 solid: &BrepSolid,
86 segments: &[ChainSegment<'_>],
87 radius: f64,
88 chamfer: bool,
89 name: Option<&str>,
90) -> Result<BrepSolid, KernelRefusal> {
91 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
92 let mut per_segment = CHAIN_PER_SEGMENT;
95 let first = loop {
96 let rung = build_closed_smooth_chain(
97 solid,
98 segments,
99 radius,
100 chamfer,
101 name,
102 per_segment,
103 FoldPolicy::Refuse,
104 bar,
105 );
106 match rung {
107 Ok(Rung::TooCoarse { deviation }) => {
108 crate::blend::carve::carve_trace(format_args!(
109 "blend chain: {per_segment} stations per segment leaves the rails \
110 {deviation:.6e} from the rolling ball, against {bar:.3e}"
111 ));
112 if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
113 return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
114 "blend: at {per_segment} stations per segment — the top of the chain's \
115 ladder — this collar's fitted rails are still {deviation:.6e} from the \
116 rolling ball's own contacts against a bar of {bar:.3e}, so there is no \
117 wall accurate enough to build"
118 )));
119 }
120 per_segment *= 2;
121 }
122 Ok(Rung::StationFailed(failure)) => {
127 crate::blend::carve::carve_trace(format_args!(
128 "blend chain: {per_segment} stations per segment: segment {} position {} \
129 (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
130 failure.segment, failure.position, failure.t, failure.from_position,
131 failure.from_t, failure.seed_uv, failure.error
132 ));
133 if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
134 break Err(failure.error);
135 }
136 per_segment *= 2;
137 }
138 other => break other,
139 }
140 };
141 if let Err(error) = &first {
142 crate::blend::carve::carve_trace(format_args!(
143 "blend chain: the march at {per_segment} stations per segment refused: {error}"
144 ));
145 }
146 match first {
147 Ok(Rung::Built(built)) => Ok(built),
148 Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
149 Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
150 }
151 Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
157 Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
160 carve_ladder(solid, segments, radius, name, bar)
161 .map_err(|carve_error| {
169 if crate::blend::is_wall_fold(&carve_error) {
170 carve_error
171 } else {
172 crate::blend::fold::wall_fold_refusal(
173 Vec::new(),
174 format!(
175 "{} {} fits this edge: the wall folds, and the fold band could not \
176 be carved — {carve_error}",
177 crate::blend::WALL_FOLDS,
178 radius.abs()
179 ),
180 )
181 }
182 })
183 }
184 Err(error) => Err(error),
185 }
186}
187
188fn carve_ladder(
198 solid: &BrepSolid,
199 segments: &[ChainSegment<'_>],
200 radius: f64,
201 name: Option<&str>,
202 bar: f64,
203) -> Result<BrepSolid, KernelRefusal> {
204 let mut per_segment = CHAIN_PER_SEGMENT;
205 loop {
206 match build_closed_smooth_chain(
207 solid,
208 segments,
209 radius,
210 false,
211 name,
212 per_segment,
213 FoldPolicy::Carve,
214 bar,
215 )? {
216 Rung::Built(built) => {
217 crate::blend::carve::carve_trace(format_args!(
218 "blend carve: {per_segment} stations per segment is the first rung inside \
219 {bar:.3e}"
220 ));
221 return Ok(built);
222 }
223 Rung::StationFailed(failure) => return Err(failure.error),
226 Rung::TooCoarse { deviation } => {
227 crate::blend::carve::carve_trace(format_args!(
228 "blend carve: {per_segment} stations per segment leaves the rails \
229 {deviation:.6e} from the rolling ball, against {bar:.3e}"
230 ));
231 if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
232 return Err(crate::blend::fold::wall_fold_refusal(
233 Vec::new(),
234 format!(
235 "{} {} fits this edge: the wall folds, and at {per_segment} stations \
236 per segment — the top of the carve's ladder — its fitted rails are \
237 still {deviation:.6e} from the rolling ball's own contacts against a \
238 bar of {bar:.3e}, so there is no wall accurate enough to carve",
239 crate::blend::WALL_FOLDS,
240 radius.abs()
241 ),
242 ));
243 }
244 per_segment *= 2;
245 }
246 }
247 }
248}
249
250enum Rung {
253 Built(BrepSolid),
254 TooCoarse { deviation: f64 },
255 StationFailed(ChainStationFailure),
258}
259
260fn rail_deviation(
284 samples: &[ChainSample],
285 cr: &NurbsCurve,
286 cs: &NurbsCurve,
287 fit_low: f64,
288 fit_range: f64,
289 per_segment: usize,
290 bar: f64,
291) -> Result<(f64, Vec<usize>), KernelRefusal> {
292 let mut order: Vec<usize> = (0..samples.len())
296 .filter(|&index| (0..=per_segment as isize).contains(&samples[index].position))
297 .collect();
298 order.sort_by(|&a, &b| {
299 let (a, b) = (&samples[a], &samples[b]);
300 a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
301 });
302 let next_of: std::collections::HashMap<usize, f64> = order
303 .windows(2)
304 .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
305 .map(|pair| (pair[0], samples[pair[1]].parameter))
306 .collect();
307 let mut failing = Vec::new();
308 let closest = closest_on_rail;
309 let mut worst: f64 = 0.0;
310 let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
313 for &index in &order {
314 let sample = &samples[index];
315 let Some([p1, p2]) = sample.midpoint else {
316 continue;
317 };
318 let Some(&next) = next_of.get(&index) else {
319 continue;
320 };
321 let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
322 let seed = (global - fit_low) / fit_range;
323 let reach = (next - sample.parameter).abs() / fit_range;
324 let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
325 worst = worst.max(miss1).max(miss2);
326 if miss1.max(miss2) > bar {
327 failing.push(index);
328 }
329 if per_segment_worst.len() <= sample.segment {
330 per_segment_worst.resize(
331 sample.segment + 1,
332 (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
333 );
334 }
335 if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
336 per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
337 }
338 }
339 for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
340 crate::blend::carve::carve_trace(format_args!(
341 " rail miss segment {segment}: worst {miss:.6e} after station {position} \
342 (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
343 p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
344 ));
345 }
346 Ok((worst, failing))
347}
348
349const PIECE_TRACK_CEILING: usize = 4096;
353
354fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
357 let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
358 let distance = |u: f64| -> Result<f64, KernelRefusal> {
359 Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
360 };
361 let mut u = seed.clamp(low, high);
362 let mut best = distance(u)?;
363 for _ in 0..40 {
364 let derivatives = curve.derivatives_extended(u, 2).or_refuse(KernelStage::Refine, "derivatives_extended")?;
365 let offset = derivatives[0].sub(target);
366 let slope = offset.dot(derivatives[1]);
367 let curvature = derivatives[1].dot(derivatives[1]) + offset.dot(derivatives[2]);
368 let mut step = if curvature > 0.0 {
371 (slope / curvature).clamp(-reach, reach)
372 } else {
373 reach.copysign(slope)
374 };
375 let mut trial = (u - step).clamp(low, high);
376 let mut trial_distance = distance(trial)?;
377 for _ in 0..30 {
378 if trial_distance < best {
379 break;
380 }
381 step *= 0.5;
382 trial = (u - step).clamp(low, high);
383 trial_distance = distance(trial)?;
384 }
385 if !(trial_distance < best) {
386 break;
387 }
388 let moved = (trial - u).abs();
389 u = trial;
390 best = trial_distance;
391 if moved < 1e-15 {
392 break;
393 }
394 }
395 Ok(best)
396}
397
398pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
401
402fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
406 (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
407}
408
409#[allow(clippy::too_many_arguments)]
415pub(super) fn dense_rail_misses(
416 contacts: &crate::blend::chain::DenseContacts,
417 samples: &[ChainSample],
418 cr: &NurbsCurve,
419 cs: &NurbsCurve,
420 low: f64,
421 range: f64,
422 wrap: bool,
423) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
424 let mut readable = Vec::with_capacity(contacts.solved.len());
425 let mut unreadable = Vec::new();
426 for (left, next, fraction, p1, p2) in &contacts.solved {
427 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
428 let global = from + fraction * (to - from);
432 let global = if wrap { global.rem_euclid(1.0) } else { global };
433 let seed = (global - low) / range;
434 let reach = (to - from).abs() / range;
435 let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
438 let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
439 match rail_pair_miss(cr_miss, cs_miss) {
440 Some(miss) => readable.push((*left, miss)),
441 None => unreadable.push(*left),
442 }
443 }
444 Ok((readable, unreadable))
445}
446
447pub(super) fn combine_rail_verdict(
453 halfway: (f64, Vec<usize>),
454 dense: &[(usize, f64)],
455 unsolved: &[usize],
456 bar: f64,
457) -> (f64, Vec<usize>) {
458 let (mut worst, mut failing) = halfway;
459 for &(left, miss) in dense {
460 if !miss.is_finite() {
461 failing.push(left);
462 continue;
463 }
464 if miss > bar {
465 failing.push(left);
466 }
467 worst = worst.max(miss);
468 }
469 failing.extend_from_slice(unsolved);
470 failing.sort_unstable();
471 failing.dedup();
472 (worst, failing)
473}
474
475pub(super) fn model_rounds_limit() -> usize {
482 REFINE_ROUNDS
483}
484
485pub(super) fn model_request_on() -> bool {
486 wall_verdict_on()
487}
488
489
490pub(super) fn wall_verdict_on() -> bool {
493 true
494}
495
496
497
498
499
500pub(super) const WALL_SECTION_PROBES: usize = 8;
503
504pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
525 let j = locate_foot(p, nodes, j, radius)?;
526 sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
527}
528
529const FOOT_TRAVEL: usize = 4;
532
533fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
543 let n = nodes.len();
544 if n < 4 || j == 0 || j + 1 >= n {
545 return None;
546 }
547 let mut tried: Vec<usize> = Vec::new();
548 let mut at = j;
549 loop {
550 let base = (at - 1).min(n - 4);
551 match stencil_solve(p, nodes, base, at, radius) {
552 Solve::Read(_) => return Some(at),
553 Solve::Unread => return None,
554 Solve::Window(s) => {
555 tried.push(at);
556 let next = (base as f64 + s).round();
557 let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
558 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
559 wall_decline(|| format!(
560 "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
561 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
562 ));
563 return None;
564 }
565 at = next;
566 }
567 }
568 }
569}
570
571fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
574 match stencil_solve(p, nodes, base, j, radius) {
575 Solve::Read(distance) => Some(distance),
576 Solve::Window(_) | Solve::Unread => None,
577 }
578}
579
580enum Solve {
583 Read(f64),
584 Window(f64),
585 Unread,
586}
587
588thread_local! {
589 static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
594}
595
596fn wall_decline(reason: impl FnOnce() -> String) {
597 WALL_DECLINES.with(|log| {
598 if let Some(lines) = log.borrow_mut().as_mut() {
599 lines.push(reason());
600 }
601 });
602}
603
604pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
612 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
613 let read = declared_wall_probe(p, fine, jf, radius);
614 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
615 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
616}
617
618fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
619 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
620 let n = nodes.len();
621 if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
622 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
623 return Solve::Unread;
624 }
625 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
627 let d1 = p1.sub(p0);
628 let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
629 let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
630 let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
631 let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
632 let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
633 let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
634 (c, c1, c2)
635 };
636 let start = j as f64 - base as f64;
637 let mut s = start;
638 let mut converged = false;
639 for _ in 0..50 {
640 let (c, c1, c2) = curve(s);
641 let w = c.sub(p);
642 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
643 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
644 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
645 return Solve::Unread;
646 }
647 let stationary = g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
660 if stationary {
661 converged = true;
662 break;
663 }
664 let step = g / slope;
665 let before = s;
666 s -= step;
667 if !((s - start).abs() <= 1.0) {
668 wall_decline(|| {
669 let (c, c1, _) = curve(s);
672 let g_after = c.sub(p).dot(c1);
673 let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
674 format!(
675 "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
676 g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
677 (s - start).abs() - 1.0
678 )
679 });
680 return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
681 }
682 if step.abs() <= 1e-14 * (1.0 + s.abs()) {
683 converged = true;
684 break;
685 }
686 }
687 if !converged {
688 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
689 return Solve::Unread;
690 }
691 let distance = curve(s).0.sub(p).length() - radius;
692 if distance.is_finite() {
693 Solve::Read(distance)
694 } else {
695 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
696 Solve::Unread
697 }
698}
699
700pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
709 let n = nodes.len();
710 if n < 5 || j == 0 || j + 1 >= n {
711 return None;
712 }
713 let j = locate_foot(p, nodes, j, radius)?;
714 let first = (j - 1).min(n - 4);
715 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
716 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
717 let a = sweep_residual_on(p, nodes, first, j, radius)?;
718 let b = sweep_residual_on(p, nodes, second, j, radius)?;
719 Some((a.abs().max(b.abs()), (a - b).abs()))
720}
721
722pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
734 if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
735 wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
736 return None;
737 }
738 let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
739 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
740 return None;
741 };
742 let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
743 wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
744 return None;
745 };
746 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
747 let Some(at_coarse) = sweep_residual(p, &coarse, jf / 2, radius) else {
748 wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
749 return None;
750 };
751 let uncertainty = spread.max((at_fine - at_coarse).abs());
752 uncertainty.is_finite().then_some((reading, uncertainty))
753}
754
755pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
771 let n = fine.len();
772 if n < 7 || !(fraction.is_finite()) {
773 return None;
774 }
775 let mut cumulative = Vec::with_capacity(n);
776 cumulative.push(0.0_f64);
777 for pair in fine.windows(2) {
778 let step = pair[1].sub(pair[0]).length();
779 cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
780 }
781 let total = *cumulative.last()?;
782 if !(total.is_finite() && total > 0.0) {
783 return None;
784 }
785 let target = fraction.clamp(0.0, 1.0) * total;
786 let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
787 let even = 2 * ((nearest as f64 / 2.0).round() as usize);
788 Some(even.clamp(2, (n - 3) & !1))
789}
790
791pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
798 let first = match declared_wall_probe_why(p, fine, jf, radius) {
799 Ok(read) => return Ok(read),
800 Err(why) => why,
801 };
802 match centre_path_node(fine, fraction) {
803 Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
804 .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
805 _ => Err(first),
806 }
807}
808
809pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
811 match declared_wall_probe(p, fine, jf, radius) {
812 Some(read) => Some(read),
813 None => match centre_path_node(fine, fraction) {
814 Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
815 _ => None,
816 },
817 }
818}
819
820pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
832 if jf % 2 != 0
833 || fine.len() < 9
834 || fine.len() % 2 == 0
835 || ts.len() != fine.len()
836 || !ts.iter().all(|t| t.is_finite())
837 || !ts.windows(2).all(|pair| pair[1] > pair[0])
838 {
839 return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
840 }
841 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
842 let read = (|| -> Option<(f64, f64)> {
843 let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
844 wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
845 None
846 })?;
847 let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
848 wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
849 None
850 })?;
851 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
852 let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
853 let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
854 wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
855 None
856 })?;
857 let uncertainty = spread.max((at_fine - at_coarse).abs());
858 uncertainty.is_finite().then_some((reading, uncertainty))
859 })();
860 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
861 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
862}
863
864fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
865 let n = nodes.len();
866 if n < 5 || j == 0 || j + 1 >= n {
867 return None;
868 }
869 let j = locate_foot_at(p, nodes, ts, j, radius)?;
870 let first = (j - 1).min(n - 4);
871 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
872 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
873 let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
874 let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
875 Some((a.abs().max(b.abs()), (a - b).abs()))
876}
877
878fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
879 let j = locate_foot_at(p, nodes, ts, j, radius)?;
880 match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
881 Solve::Read(distance) => Some(distance),
882 Solve::Window(_) | Solve::Unread => None,
883 }
884}
885
886fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
890 let n = nodes.len();
891 if n < 4 || j == 0 || j + 1 >= n {
892 return None;
893 }
894 let mut tried: Vec<usize> = Vec::new();
895 let mut at = j;
896 loop {
897 let base = (at - 1).min(n - 4);
898 match stencil_solve_at(p, nodes, ts, base, at, radius) {
899 Solve::Read(_) => return Some(at),
900 Solve::Unread => return None,
901 Solve::Window(t) => {
902 tried.push(at);
903 let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
904 let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
905 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
906 wall_decline(|| format!(
907 "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
908 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
909 ));
910 return None;
911 }
912 at = next;
913 }
914 }
915 }
916}
917
918fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
923 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
924 let n = nodes.len();
925 if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || ts.len() != n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
926 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
927 return Solve::Unread;
928 }
929 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
930 let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
931 let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
932 let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
933 let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
934 let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
935 let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
936 let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
937 let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
938 let (a, b, c) = (t - t0, t - t1, t - t2);
939 let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
940 let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
941 let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
942 (point, first, second)
943 };
944 let (low, high) = (ts[j - 1], ts[j + 1]);
945 let start = ts[j];
946 let mut t = start;
947 let mut converged = false;
948 for _ in 0..50 {
949 let (c, c1, c2) = curve(t);
950 let w = c.sub(p);
951 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
952 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
953 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
954 return Solve::Unread;
955 }
956 if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
957 converged = true;
958 break;
959 }
960 let step = g / slope;
961 let before = t;
962 t -= step;
963 if !(t >= low && t <= high) {
964 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
965 return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
966 }
967 let node = 0.5 * (high - low);
972 if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
973 converged = true;
974 break;
975 }
976 }
977 if !converged {
978 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
979 return Solve::Unread;
980 }
981 let distance = curve(t).0.sub(p).length() - radius;
982 if distance.is_finite() {
983 Solve::Read(distance)
984 } else {
985 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
986 Solve::Unread
987 }
988}
989
990pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
995 readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
996}
997
998#[allow(clippy::too_many_arguments)]
1008pub(super) fn wall_interior_misses(
1009 contacts: &crate::blend::chain::DenseContacts,
1010 samples: &[ChainSample],
1011 surface: &crate::NurbsSurface,
1012 low: f64,
1013 range: f64,
1014 radius: f64,
1015 subdivisions: usize,
1016 probes: usize,
1017 wrap: bool,
1018) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1019 use std::collections::BTreeMap;
1020 let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1023 for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1024 let entry = nodes.entry(*left).or_insert_with(|| {
1025 let mut row = vec![None; subdivisions + 1];
1026 row[0] = Some(samples[*left].station.center);
1027 row[subdivisions] = Some(samples[*next].station.center);
1028 (*next, row)
1029 });
1030 let step = (fraction * subdivisions as f64).round() as usize;
1031 if (1..subdivisions).contains(&step) {
1032 entry.1[step] = contacts.centers.get(index).copied().flatten();
1033 }
1034 }
1035 let mut readable = Vec::new();
1036 let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1037 for (&left, (next, row)) in &nodes {
1038 let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1039 unreadable.push(left);
1040 continue;
1041 };
1042 let (from, to) = (samples[left].parameter, samples[*next].parameter);
1043 let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1044 let mut read = true;
1045 'probes: for step in (2..subdivisions).step_by(2) {
1046 let global = from + (step as f64 / subdivisions as f64) * (to - from);
1047 let global = if wrap { global.rem_euclid(1.0) } else { global };
1048 let u = (global - low) / range;
1049 for k in 1..probes {
1050 let v = k as f64 / probes as f64;
1051 let reading = surface.evaluate(u, v).ok().and_then(|p| {
1052 declared_wall_probe_located(p, ¢res, step, step as f64 / subdivisions as f64, radius)
1053 });
1054 match reading {
1055 Some((distance, gap)) => {
1056 declared = declared.max(distance + gap);
1057 worst = worst.max(distance);
1058 spread = spread.max(gap);
1059 }
1060 None => {
1061 read = false;
1062 break 'probes;
1063 }
1064 }
1065 }
1066 }
1067 if read {
1068 readable.push((left, declared, worst, spread));
1069 } else {
1070 unreadable.push(left);
1071 }
1072 }
1073 unreadable.sort_unstable();
1074 unreadable.dedup();
1075 (readable, unreadable)
1076}
1077
1078
1079fn build_closed_smooth_chain(
1080 solid: &BrepSolid,
1081 segments: &[ChainSegment<'_>],
1082 radius: f64,
1083 chamfer: bool,
1084 name: Option<&str>,
1085 per_segment: usize,
1086 fold: FoldPolicy,
1087 bar: f64,
1088) -> Result<Rung, KernelRefusal> {
1089 let samples = match march_chain_typed(segments, radius, false, per_segment, fold) {
1090 Ok(samples) => samples,
1091 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1092 Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1095 Err(ChainMarchError::Other(error)) => return Err(error),
1096 };
1097 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1098 build_closed_chain_from_samples(solid, segments, radius, chamfer, name, per_segment, fold, bar, samples, 0)
1099}
1100
1101
1102thread_local! {
1103 static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1107}
1108
1109#[allow(clippy::too_many_arguments)]
1118fn build_closed_chain_from_samples(
1119 solid: &BrepSolid,
1120 segments: &[ChainSegment<'_>],
1121 radius: f64,
1122 chamfer: bool,
1123 name: Option<&str>,
1124 per_segment: usize,
1125 fold: FoldPolicy,
1126 bar: f64,
1127 mut samples: Vec<ChainSample>,
1128 round: usize,
1129) -> Result<Rung, KernelRefusal> {
1130 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1131
1132 let mut global: Vec<(f64, &ChainSample)> = samples
1134 .iter()
1135 .filter(|sample| (0..per_segment as isize).contains(&sample.position))
1136 .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1137 .collect();
1138 global.sort_by(|a, b| a.0.total_cmp(&b.0));
1139 let count = global.len();
1140 let degree = FIT_DEGREE;
1141 let mut extended_params = Vec::new();
1142 let mut samples_cr = Vec::new();
1143 let mut samples_mid = Vec::new();
1144 let mut samples_cs = Vec::new();
1145 let mut push = |station: &Station, parameter: f64| {
1146 extended_params.push(parameter);
1147 samples_cr.push(Vec4::from_point(station.p1, 1.0));
1148 samples_cs.push(Vec4::from_point(station.p2, 1.0));
1149 samples_mid.push(Vec4 {
1150 x: station.apex.x * station.weight,
1151 y: station.apex.y * station.weight,
1152 z: station.apex.z * station.weight,
1153 w: station.weight,
1154 });
1155 };
1156 for offset in (1..=degree).rev() {
1157 let (parameter, sample) = &global[count - offset];
1158 push(&sample.station, parameter - 1.0);
1159 }
1160 for (parameter, sample) in &global {
1161 push(&sample.station, *parameter);
1162 }
1163 push(&global[0].1.station, global[0].0 + 1.0);
1166 for offset in 1..=degree {
1167 let (parameter, sample) = &global[offset];
1168 push(&sample.station, parameter + 1.0);
1169 }
1170 let low = extended_params[0];
1171 let high = *extended_params.last().unwrap();
1172 let range = high - low;
1173 let normalized: Vec<f64> = extended_params
1174 .iter()
1175 .map(|parameter| (parameter - low) / range)
1176 .collect();
1177 let seam_low = (0.0 - low) / range;
1178 let seam_high = (1.0 - low) / range;
1179 let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1180 let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1181 let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1182 let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1183 Ok(middle)
1184 };
1185 let cr = fit_row(&samples_cr)?;
1186 let cs = fit_row(&samples_cs)?;
1187 let mid = if chamfer {
1188 None
1189 } else {
1190 Some(fit_row(&samples_mid)?)
1191 };
1192 let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1193 let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1194
1195 let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, per_segment, bar)?;
1208 let stencil = if round > 0 {
1220 Some(crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, LOCAL_RAIL_STENCIL)?)
1221 } else {
1222 None
1223 };
1224 let wall_stencil = if !chamfer && wall_verdict_on() {
1229 Some(crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1230 } else {
1231 None
1232 };
1233 let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1234 let contacts = stencil.as_ref().expect("the local path solves the stencil");
1235 let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1236 let reads = dense.len() + unreadable.len() + contacts.unsolved;
1237 let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1238 let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1239 crate::blend::carve::carve_trace(format_args!(
1240 "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1241 {} failing interval(s), {} unsolved contact(s)",
1242 verdict.0, verdict.1.len(), contacts.unsolved
1243 ));
1244 (verdict.0, verdict.1, reads)
1245 } else {
1246 (halfway, halfway_failing, 0)
1247 };
1248 let mut _wall_worst = 0.0_f64;
1254 let mut _wall_spread = 0.0_f64;
1255 let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1256 if let Some(contacts) = wall_stencil.as_ref() {
1257 let (wall, unreadable) =
1258 wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1259 wall_readings = wall.clone();
1260 for &(left, declared, _, uncertainty) in &wall {
1262 if declared > bar {
1263 failing.push(left);
1264 }
1265 _wall_worst = _wall_worst.max(declared);
1266 _wall_spread = _wall_spread.max(uncertainty);
1267 }
1268 failing.extend_from_slice(&unreadable);
1269 failing.sort_unstable();
1270 failing.dedup();
1271 deviation = deviation.max(_wall_worst);
1272 crate::blend::carve::carve_trace(format_args!(
1273 "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1274 interval(s), {} unreadable, bar {bar:.3e}",
1275 wall.len(),
1276 unreadable.len()
1277 ));
1278 }
1279 if deviation > bar || !failing.is_empty() {
1280 let in_segment = |segment: usize| {
1281 samples
1282 .iter()
1283 .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1284 .count()
1285 };
1286 let within_ceiling = (0..segments.len()).all(|segment| {
1287 in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1288 <= CHAIN_CARVE_MAX_PER_SEGMENT
1289 });
1290 crate::blend::carve::carve_trace(format_args!(
1291 "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1292 worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1293 samples.len(),
1294 failing.len(),
1295 if within_ceiling { "ok" } else { "reached" }
1296 ));
1297 if fold != FoldPolicy::Refuse || round >= REFINE_ROUNDS || !within_ceiling {
1298 return Ok(Rung::TooCoarse { deviation });
1299 }
1300 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1302 match insert_chain_stations(segments, radius, per_segment, &mut samples, &failing) {
1303 Ok(()) => {}
1304 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1305 Err(ChainMarchError::Incoherent(failure)) => {
1309 crate::blend::carve::carve_trace(format_args!(
1310 "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1311 at t {:.9} does not reach its retained right neighbour: {}",
1312 failure.segment, failure.t, failure.error
1313 ));
1314 return Ok(Rung::TooCoarse { deviation });
1315 }
1316 Err(ChainMarchError::Other(error)) => return Err(error),
1317 }
1318 assign_chain_parameters(&mut samples, segments.len(), per_segment, false);
1319 return build_closed_chain_from_samples(
1320 solid, segments, radius, chamfer, name, per_segment, fold, bar, samples, round + 1,
1321 );
1322 }
1323 crate::blend::carve::carve_trace(format_args!(
1324 "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1325 from the rolling ball, inside {bar:.3e}"
1326 ));
1327 let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1341 let stations_now = (0..segments.len())
1342 .map(|segment| {
1343 samples
1344 .iter()
1345 .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1346 .count()
1347 })
1348 .max()
1349 .unwrap_or(0);
1350 let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1351 let mut request_detail: Option<String> = None;
1355 let mut request_component = crate::MeasuredComponent::Wall;
1356 let mut request_unread = 0usize;
1357 if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1358 let wall_missing = wall_model_missing(&wall_readings, model);
1359 let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1367 Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1368 None => (Vec::new(), Vec::new()),
1369 };
1370 let rail_over: std::collections::BTreeSet<usize> =
1372 rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1373 let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1374 let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1375 let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1376 missing.sort_unstable();
1377 missing.dedup();
1378 if !missing.is_empty() {
1379 let wall_short = wall_readings
1381 .iter()
1382 .filter(|reading| wall_missing.contains(&reading.0))
1383 .map(|&(_, declared, _, _)| declared)
1384 .fold(0.0_f64, f64::max);
1385 let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1386 let short = wall_short.max(rail_short);
1387 request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1388 request_unread = rail_unread_intervals.len();
1389 request_detail = Some(format!(
1390 "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1391 wall {wall_short:.3e} over {} interval(s)",
1392 rail_over.len(),
1393 rail_unread_intervals.len(),
1394 wall_missing.len()
1395 ));
1396 let within_ceiling = (0..segments.len()).all(|segment| {
1397 let present = samples
1398 .iter()
1399 .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1400 .count();
1401 present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1402 <= CHAIN_CARVE_MAX_PER_SEGMENT
1403 });
1404 crate::blend::carve::carve_trace(format_args!(
1405 "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1406 (worst reading + spread {short:.3e}), ceiling {}",
1407 missing.len(),
1408 if within_ceiling { "ok" } else { "reached" }
1409 ));
1410 unmet = if round >= model_rounds_limit() {
1411 Some((short, crate::BudgetReason::RoundsSpent))
1412 } else if !within_ceiling {
1413 Some((short, crate::BudgetReason::StationCeiling))
1414 } else {
1415 let mut refined = samples.clone();
1416 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1418 match insert_chain_stations(segments, radius, per_segment, &mut refined, &missing) {
1419 Ok(()) => {
1420 assign_chain_parameters(&mut refined, segments.len(), per_segment, false);
1421 match build_closed_chain_from_samples(
1422 solid, segments, radius, chamfer, name, per_segment, fold, bar, refined, round + 1,
1423 ) {
1424 Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1425 _ => Some((short, crate::BudgetReason::Rejected)),
1426 }
1427 }
1428 Err(_) => {
1429 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1431 Some((short, crate::BudgetReason::Incoherent))
1432 }
1433 }
1434 };
1435 }
1436 }
1437 if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1442 let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, subdivisions.max(2))?;
1443 let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1444 let mut unreadable = 0usize;
1445 for (left, next, fraction, p1, p2) in &contacts.solved {
1446 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1447 let global = (from + fraction * (to - from)).rem_euclid(1.0);
1448 let seed = (global - low) / range;
1449 let reach = (to - from).abs() / range;
1450 match rail_pair_miss(closest_on_rail(&cr, *p1, seed, reach).unwrap_or(f64::NAN), closest_on_rail(&cs, *p2, seed, reach).unwrap_or(f64::NAN)) {
1451 Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1452 Some(_) => {}
1453 None => unreadable += 1,
1454 }
1455 }
1456 eprintln!(
1457 "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1458 worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1459 contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1460 );
1461 }
1462 let crease = if fold == FoldPolicy::Carve {
1463 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1473 let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1474 return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1475 "blend carve: the wall re-marched at {per_segment} stations per segment carries no fold to carve, but the march at {CHAIN_PER_SEGMENT} refused one"
1476 )));
1477 };
1478 Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1479 } else {
1480 None
1481 };
1482
1483 let single_face = |side: usize| -> bool {
1487 let first_id = segments[0].mate(side).face.id;
1488 segments
1489 .iter()
1490 .all(|segment| segment.mate(side).face.id == first_id)
1491 };
1492 let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1493 for side in 0..2 {
1494 if !single_face(side) {
1495 continue;
1496 }
1497 let mut row = Vec::new();
1498 let select = |station: &Station| -> [f64; 2] {
1499 if side == 0 {
1500 station.uv1
1501 } else {
1502 station.uv2
1503 }
1504 };
1505 let mut push_uv = |station: &Station| {
1506 let uv = select(station);
1507 row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1508 };
1509 for offset in (1..=degree).rev() {
1510 push_uv(&global[count - offset].1.station);
1511 }
1512 for (_, sample) in &global {
1513 push_uv(&sample.station);
1514 }
1515 push_uv(&global[0].1.station);
1516 for offset in 1..=degree {
1517 push_uv(&global[offset].1.station);
1518 }
1519 whole_pcurves[side] = Some(fit_row(&row)?);
1520 }
1521
1522 let mut pcurves1 = Vec::with_capacity(segments.len());
1525 let mut pcurves2 = Vec::with_capacity(segments.len());
1526 for segment in 0..segments.len() {
1527 let mut window: Vec<&ChainSample> = samples
1528 .iter()
1529 .filter(|sample| sample.segment == segment)
1530 .collect();
1531 window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1532 let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1533 let low = params[0];
1534 let high = *params.last().unwrap();
1535 let normalized: Vec<f64> = params
1536 .iter()
1537 .map(|parameter| (parameter - low) / (high - low))
1538 .collect();
1539 let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1540 let points: Vec<Vec4> = window
1541 .iter()
1542 .map(|sample| {
1543 let uv = select(&sample.station);
1544 Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1545 })
1546 .collect();
1547 let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1548 Ok(curve)
1552 };
1553 pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1554 pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1555 }
1556
1557 let wall_surface = unmet.as_ref().map(|_| surface.clone());
1558 let built = chain_surgery(
1559 solid,
1560 segments,
1561 ChainRows {
1562 surface,
1563 cr,
1564 cs,
1565 u_domain,
1566 fit_low: low,
1567 fit_range: range,
1568 pcurves1,
1569 pcurves2,
1570 whole_pcurves,
1571 crease,
1572 },
1573 name,
1574 )?;
1575 if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1576 crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1577 name: name.map(str::to_string),
1578 surface,
1579 request: model,
1580 residual,
1581 detail: request_detail,
1582 budget: crate::ApproximationBudget {
1583 reason,
1584 rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1587 rounds_limit: model_rounds_limit(),
1588 stations: stations_now,
1589 station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1590 mechanism: crate::BudgetMechanism::LocalRounds,
1591 measured_component: request_component,
1592 unread: request_unread,
1593 },
1594 });
1595 }
1596 Ok(Rung::Built(built))
1597}
1598
1599pub(super) struct ChainRows {
1600 pub(super) surface: NurbsSurface,
1601 pub(super) cr: NurbsCurve,
1602 pub(super) cs: NurbsCurve,
1603 pub(super) u_domain: [f64; 2],
1604 pub(super) fit_low: f64,
1608 pub(super) fit_range: f64,
1609 pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1612 pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1613 pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1615 pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1618}
1619
1620impl ChainRows {
1621 pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1622 self.fit_low + fit_parameter * self.fit_range
1623 }
1624}
1625
1626pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1630 let (low, high, curve) = fitted;
1631 let to_local = |value: f64| {
1634 let mut best = value;
1635 for candidate in [value, value + 1.0, value - 1.0] {
1636 if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1637 best = candidate;
1638 break;
1639 }
1640 }
1641 ((best - low) / (high - low)).clamp(0.0, 1.0)
1642 };
1643 let mut local_a = to_local(a);
1644 let mut local_b = to_local(b);
1645 if local_a > local_b {
1646 std::mem::swap(&mut local_a, &mut local_b);
1647 }
1648 let epsilon = 1e-9;
1649 let (_, tail) = if local_a > epsilon {
1650 curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1651 } else {
1652 (curve.clone(), curve.clone())
1653 };
1654 let tail = if local_a > epsilon {
1655 tail
1656 } else {
1657 curve.clone()
1658 };
1659 let portion = if local_b < 1.0 - epsilon {
1660 tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1661 } else {
1662 tail
1663 };
1664 Ok(portion)
1665}
1666
1667pub(super) fn cross_edge_at(
1673 solid: &BrepSolid,
1674 face_before: &FaceRecord,
1675 loop_before: usize,
1676 face_after: &FaceRecord,
1677 loop_after: usize,
1678 vertex: u64,
1679 before_edge: u64,
1680 after_edge: u64,
1681) -> Result<Option<u64>, KernelRefusal> {
1682 let mut found: Option<u64> = None;
1683 let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1684 for coedge in &face.loops[loop_index].coedges {
1685 if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1686 continue;
1687 }
1688 let edge = solid
1689 .edges
1690 .iter()
1691 .find(|edge| edge.id == coedge.edge_id)
1692 .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1693 if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1694 if found.is_some() && found != Some(edge.id) {
1695 return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1696 }
1697 found = Some(edge.id);
1698 }
1699 }
1700 Ok(())
1701 };
1702 consider(face_before, loop_before)?;
1703 if !(face_after.id == face_before.id && loop_after == loop_before) {
1704 consider(face_after, loop_after)?;
1705 }
1706 Ok(found)
1707}
1708
1709fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1712 let period = high - low;
1713 low + ((value - low) / period).round() * period
1714}
1715
1716pub(in crate::blend) fn project_piece_pcurve(
1734 piece: &NurbsCurve,
1735 surface: &NurbsSurface,
1736 start_is_crossing: bool,
1737 end_is_crossing: bool,
1738) -> Result<NurbsCurve, KernelRefusal> {
1739 project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1740}
1741
1742
1743pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
1748 let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
1749 let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1750 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1751 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1752 let into = move |value: f64, closed: bool, low: f64, high: f64| {
1753 if closed && high > low {
1754 low + (value - low).rem_euclid(high - low)
1755 } else {
1756 value.clamp(low, high)
1757 }
1758 };
1759 Ok(move |fraction: f64| {
1760 let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
1761 let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
1762 seed.iter().all(|value| value.is_finite()).then_some(seed)
1763 })
1764}
1765
1766pub(in crate::blend) fn project_piece_pcurve_from(
1776 piece: &NurbsCurve,
1777 surface: &NurbsSurface,
1778 start_is_crossing: bool,
1779 end_is_crossing: bool,
1780 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1781) -> Result<NurbsCurve, KernelRefusal> {
1782 project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
1783}
1784
1785pub(in crate::blend) fn project_piece_pcurve_within(
1790 piece: &NurbsCurve,
1791 surface: &NurbsSurface,
1792 start_is_crossing: bool,
1793 end_is_crossing: bool,
1794 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1795 refine_near: &[f64],
1796) -> Result<NurbsCurve, KernelRefusal> {
1797 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1798 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1799 let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1800 let period = u1 - u0;
1801 let period_v = v1 - v0;
1802 let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
1803 let dense = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
1811 let mut proj_u = Vec::with_capacity(dense + 1);
1812 let mut proj_v = Vec::with_capacity(dense + 1);
1813 let mut point3 = Vec::with_capacity(dense + 1);
1814 let mut params = Vec::with_capacity(dense + 1);
1815 let mut previous_u: Option<f64> = None;
1816 let mut previous_foot: Option<[f64; 2]> = None;
1821 let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
1822 for section in 0..=dense {
1823 let t = d0 + (d1 - d0) * section as f64 / dense as f64;
1824 let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
1825 let guide = guide.filter(|_| {
1826 true
1827 });
1828 let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(section as f64 / dense as f64)).into_iter().chain(previous_foot).collect();
1833 let mut projection: Option<crate::projection::SurfaceProjection> = None;
1834 for seed in &seeds {
1835 let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
1836 if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
1837 projection = Some(found);
1838 }
1839 }
1840 let projection = match projection {
1841 Some(projection) => projection,
1842 None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
1843 };
1844 previous_foot = Some([projection.u, projection.v]);
1845 feet.push([projection.u, projection.v]);
1846 let mut u = projection.u;
1847 if let Some(previous) = previous_u {
1848 while u - previous > period * 0.5 {
1849 u -= period;
1850 }
1851 while previous - u > period * 0.5 {
1852 u += period;
1853 }
1854 }
1855 previous_u = Some(u);
1856 proj_u.push(u);
1857 proj_v.push(projection.v);
1858 point3.push(point);
1859 params.push(section as f64 / dense as f64);
1860 }
1861 if closed_v {
1868 for index in 2..=dense {
1869 while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
1870 proj_v[index] -= period_v;
1871 }
1872 while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
1873 proj_v[index] += period_v;
1874 }
1875 }
1876 while proj_v[0] - proj_v[1] > period_v * 0.5 {
1877 proj_v[0] -= period_v;
1878 }
1879 while proj_v[1] - proj_v[0] > period_v * 0.5 {
1880 proj_v[0] += period_v;
1881 }
1882 }
1883 let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
1886 let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
1887 let pinned = |index: usize,
1896 neighbour: usize,
1897 proj_u: &[f64],
1898 proj_v: &[f64]|
1899 -> Result<(f64, f64), KernelRefusal> {
1900 let point = point3[index];
1901 let u_error = surface
1902 .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
1903 .sub(point)
1904 .length();
1905 if closed_v {
1906 let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
1907 let v_error = surface
1908 .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
1909 .sub(point)
1910 .length();
1911 if v_error < u_error {
1912 return Ok((proj_u[index], seam_v));
1913 }
1914 }
1915 Ok((seam_u, proj_v[index]))
1916 };
1917 if start_is_crossing {
1918 let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
1919 proj_u[0] = u;
1920 proj_v[0] = v;
1921 }
1922 if end_is_crossing {
1923 let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
1924 proj_u[dense] = u;
1925 proj_v[dense] = v;
1926 }
1927 let mut shift = 0.0;
1928 while mean_u + shift > u1 + 1e-9 {
1929 shift -= period;
1930 }
1931 while mean_u + shift < u0 - 1e-9 {
1932 shift += period;
1933 }
1934 for u in proj_u.iter_mut() {
1935 *u += shift;
1936 }
1937 if closed_v {
1938 let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
1939 let mut shift_v = 0.0;
1940 while mean_v + shift_v > v1 + 1e-9 {
1941 shift_v -= period_v;
1942 }
1943 while mean_v + shift_v < v0 - 1e-9 {
1944 shift_v += period_v;
1945 }
1946 for v in proj_v.iter_mut() {
1947 *v += shift_v;
1948 }
1949 }
1950 let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
1965 let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
1972 let mut track = crate::blend::track_fit::Track {
1973 fractions: params,
1974 uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
1975 feet,
1976 breaks: Vec::new(),
1977 };
1978 let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
1979 crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
1980 crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
1981 for &fraction in refine_near {
1982 for _ in 0..2 {
1983 let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
1984 if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], PIECE_TRACK_CEILING)?.is_none() {
1985 break;
1986 }
1987 }
1988 }
1989 let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
1996 surface,
1997 &mut track,
1998 &at,
1999 tolerance,
2000 REFINE_ROUNDS,
2001 PIECE_TRACK_CEILING,
2002 )?;
2003 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2004 eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2005 eprintln!(
2006 "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2007 fit.miss,
2008 fit.samples,
2009 track.fractions.len(),
2010 track.breaks.len(),
2011 fit.on_floor
2012 );
2013 if let Some(worst) = fit.worst {
2018 let parameter = d0 + (d1 - d0) * worst.fraction;
2019 let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2020 .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2021 .map(|pair| [pair[0], pair[1]]).collect();
2022 let mean = (d1 - d0) / spans.len().max(1) as f64;
2023 let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2024 let width = span.map(|index| spans[index][1] - spans[index][0]);
2025 eprintln!(
2026 "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2027 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2028 worst.fraction, parameter, span, spans.len(), width, mean,
2029 worst.bracket[0], worst.bracket[1], worst.interval
2030 );
2031 }
2032 }
2033 if !fit.on_floor {
2040 return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_pcurve_floor", format!(
2041 "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2042 against a floor of {tolerance:.1e}",
2043 crate::blend::PCURVE_OFF_FLOOR,
2044 fit.miss,
2045 fit.samples
2046 )));
2047 }
2048 Ok(fit.curve)
2049}