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 let pairs: Vec<[u64; 2]> = chain.segments.iter().map(|segment| [segment.first.face.id, segment.second.face.id]).collect();
37 let _sides = crate::blend::setback::chain_sides(&pairs)?;
38 if chain.closed {
39 blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
40 } else {
41 blend_open_smooth_chain(solid, &chain, radius, chamfer, name)
42 }
43}
44
45pub fn blend_smooth_chain_if_closed(
52 solid: &BrepSolid,
53 seed_edge_id: u64,
54 radius: f64,
55 chamfer: bool,
56 name: Option<&str>,
57) -> Result<BrepSolid, KernelRefusal> {
58 let _operation = crate::blend::BlendOperation::enter();
61 if !(radius > 0.0) || !radius.is_finite() {
62 return Err(KernelRefusal::input(KernelStage::Collect, "radius", "blend: radius must be positive"));
63 }
64 let chain = collect_smooth_chain(solid, seed_edge_id)?;
65 if chain.segments.len() < 2 {
66 return Err(KernelRefusal::internal(KernelStage::Refine, "chain_single_segment", "blend: chain collapsed to a single segment"));
67 }
68 let pairs: Vec<[u64; 2]> = chain.segments.iter().map(|segment| [segment.first.face.id, segment.second.face.id]).collect();
70 let _sides = crate::blend::setback::chain_sides(&pairs)?;
71 if !chain.closed {
72 return Err(
73 KernelRefusal::unsupported(KernelStage::Classify, "open_tangent_chain", "blend: the edge is one arc of an OPEN tangent chain; the unselected \
74 continuation is capped, not blended"),
75 );
76 }
77 blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
78}
79
80fn blend_closed_smooth_chain(
97 solid: &BrepSolid,
98 segments: &[ChainSegment<'_>],
99 radius: f64,
100 chamfer: bool,
101 name: Option<&str>,
102) -> Result<BrepSolid, KernelRefusal> {
103 let mode = std::env::var("BREP_CHAIN_LADDER_MEMO");
104 let (memo_on, verify) = match mode.as_deref() {
105 Ok("0") => (false, false),
106 Ok("verify") => (true, true),
107 _ => (true, false),
108 };
109 if !memo_on || crate::blend::setback::sides().is_some() {
111 return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
112 }
113 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
114 let signed = if chain_ladder_signed_key_on() {
123 let mut signs = Vec::with_capacity(segments.len());
124 for segment in segments {
125 match crate::blend::stations::signed_radii(
126 segment.edge,
127 segment.first.face,
128 segment.first.coedge,
129 segment.second.face,
130 segment.second.coedge,
131 radius,
132 ) {
133 Ok(pair) => signs.push(pair),
134 Err(_) => return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name),
136 }
137 }
138 Some(signs)
139 } else {
140 None
141 };
142 let key = ladder_memo_key(segments, radius, chamfer, bar, signed.as_deref());
143 let cached = LADDER_MEMO.with(|memo| memo.borrow().iter().find(|(stored, _)| *stored == key).map(|(_, refusal)| refusal.clone()));
144 crate::blend::carve::carve_trace(format_args!(
145 "blend chain ladder memo: {} for segments {:?}",
146 if cached.is_some() { "hit" } else { "miss" },
147 segments
148 .iter()
149 .map(|segment| (segment.edge.id, segment.forward, segment.first.face.id, segment.second.face.id))
150 .collect::<Vec<_>>()
151 ));
152 if let (Some(refusal), false) = (&cached, verify) {
153 crate::blend::record_count("count.chain_ladder_memo.hit", 1);
154 return Err(refusal.clone());
155 }
156 let result = blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
157 if let Some(refusal) = cached {
158 let same = matches!(&result, Err(again) if *again == refusal && format!("{again:?}") == format!("{refusal:?}"));
159 assert!(same, "chain ladder memo mismatch: cached {refusal:?}, recomputed {:?}", result.as_ref().err());
160 crate::blend::record_count("count.chain_ladder_memo.verified", 1);
161 return result;
162 }
163 if let Err(refusal) = &result {
164 if matches!(&refusal.class, crate::RefusalClass::NonConvergence { what } if what == "chain_ladder_top") {
165 LADDER_MEMO.with(|memo| {
166 let mut memo = memo.borrow_mut();
167 if memo.len() >= LADDER_MEMO_CAPACITY {
168 memo.remove(0);
169 }
170 memo.push((key, refusal.clone()));
171 });
172 }
173 }
174 result
175}
176
177fn chain_ladder_signed_key_on() -> bool {
181 switch!("BREP_CHAIN_LADDER_SIGNED_KEY", true)
182}
183
184fn chain_ladder_rotation_key_on() -> bool {
188 switch!("BREP_CHAIN_LADDER_ROTATION_KEY", true)
189}
190
191const LADDER_MEMO_CAPACITY: usize = 32;
193
194thread_local! {
195 static LADDER_MEMO: std::cell::RefCell<Vec<(Vec<u64>, KernelRefusal)>> = const { std::cell::RefCell::new(Vec::new()) };
197}
198
199fn ladder_memo_key(segments: &[ChainSegment<'_>], radius: f64, chamfer: bool, bar: f64, signed: Option<&[(f64, f64)]>) -> Vec<u64> {
201 fn curve(key: &mut Vec<u64>, curve: &NurbsCurve) {
202 key.extend([curve.degree as u64, curve.knots.len() as u64, curve.control_points.len() as u64]);
203 key.extend(curve.knots.iter().map(|knot| knot.to_bits()));
204 for point in &curve.control_points {
205 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
206 }
207 }
208 fn face(key: &mut Vec<u64>, face: &FaceRecord, loops: bool) {
209 let surface = &face.surface;
210 key.extend([face.id, face.same_sense as u64, surface.degree_u as u64, surface.degree_v as u64]);
211 key.extend([surface.knots_u.len() as u64, surface.knots_v.len() as u64, surface.control_points.len() as u64]);
212 key.extend(surface.knots_u.iter().chain(&surface.knots_v).map(|knot| knot.to_bits()));
213 for row in &surface.control_points {
214 key.push(row.len() as u64);
215 for point in row {
216 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
217 }
218 }
219 if !loops {
220 return;
221 }
222 key.push(face.loops.len() as u64);
223 for loop_record in &face.loops {
224 key.extend([loop_record.id, loop_record.coedges.len() as u64]);
225 for coedge in &loop_record.coedges {
226 key.extend([coedge.id, coedge.edge_id, coedge.forward as u64]);
227 curve(key, &coedge.pcurve);
228 }
229 }
230 }
231 let mut key = vec![
232 radius.to_bits(),
233 chamfer as u64,
234 bar.to_bits(),
235 crate::blend::station_refinement_on() as u64,
236 crate::blend::chain_arc_stations_on() as u64,
237 wall_verdict_on() as u64,
238 local_continue_on() as u64,
239 segments.len() as u64,
240 signed.is_some() as u64,
241 ];
242 let start = if chain_ladder_rotation_key_on() {
247 key.push(1);
248 segments.iter().enumerate().min_by_key(|(_, segment)| segment.edge.id).map_or(0, |(index, _)| index)
249 } else {
250 0
251 };
252 for index in (0..segments.len()).map(|offset| (start + offset) % segments.len()) {
253 let segment = &segments[index];
254 if let Some((rho1, rho2)) = signed.and_then(|signs| signs.get(index)) {
255 key.extend([rho1.to_bits(), rho2.to_bits()]);
256 }
257 let edge = segment.edge;
258 key.extend([edge.id, edge.start_vertex_id, edge.end_vertex_id, edge.t0.to_bits(), edge.t1.to_bits(), segment.forward as u64]);
259 curve(&mut key, &edge.curve);
260 for mate in [&segment.first, &segment.second] {
261 key.extend([mate.coedge.id, mate.coedge.edge_id, mate.coedge.forward as u64, mate.loop_index as u64, mate.rho.to_bits()]);
262 curve(&mut key, &mate.coedge.pcurve);
263 face(&mut key, mate.face, signed.is_none());
264 }
265 }
266 key
267}
268fn blend_closed_smooth_chain_unmemoized(
281 solid: &BrepSolid,
282 segments: &[ChainSegment<'_>],
283 radius: f64,
284 chamfer: bool,
285 name: Option<&str>,
286) -> Result<BrepSolid, KernelRefusal> {
287 let bar = crate::blend::setback::chain_bar(solid, crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit);
288 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
291 let mut trend: Vec<(usize, f64)> = Vec::new();
293 let first = loop {
294 let per_segment = budget.rung;
295 let rung = build_closed_smooth_chain(
296 solid,
297 segments,
298 radius,
299 chamfer,
300 name,
301 &budget,
302 FoldPolicy::Refuse,
303 bar,
304 );
305 match rung {
306 Ok(Rung::TooCoarse { deviation }) => {
307 crate::blend::carve::carve_trace(format_args!(
308 "blend chain: {} leaves the rails \
309 {deviation:.6e} from the rolling ball, against {bar:.3e}",
310 budget.describe()
311 ));
312 trend.push((per_segment, deviation));
313 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
314 let (best_rung, best, listing) = ladder_trend(&trend);
315 return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
316 "blend: the chain's ladder reached its top, {per_segment} stations per segment, \
317 and no rung passed — its best, {best_rung} stations per segment, left this \
318 collar's fitted rails {best:.6e} from the rolling ball's own contacts against \
319 a bar of {bar:.3e} (rung by rung: {listing}), so there is no wall accurate \
320 enough to build"
321 )));
322 }
323 budget = budget.doubled();
324 }
325 Ok(Rung::StationFailed(failure)) => {
330 crate::blend::carve::carve_trace(format_args!(
331 "blend chain: {per_segment} stations per segment: segment {} position {} \
332 (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
333 failure.segment, failure.position, failure.t, failure.from_position,
334 failure.from_t, failure.seed_uv, failure.error
335 ));
336 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
337 break Err(failure.error);
338 }
339 budget = budget.doubled();
340 }
341 other => break other,
342 }
343 };
344 if let Err(error) = &first {
345 crate::blend::carve::carve_trace(format_args!(
346 "blend chain: the march at {} refused: {error}", budget.describe()
347 ));
348 }
349 match first {
350 Ok(Rung::Built(built)) => Ok(built),
351 Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
352 Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
353 }
354 Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
360 Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
363 carve_ladder(solid, segments, radius, name, bar)
364 .map_err(|carve_error| {
372 if crate::blend::is_wall_fold(&carve_error) {
373 carve_error
374 } else {
375 crate::blend::fold::wall_fold_refusal(
376 Vec::new(),
377 format!(
378 "{} {} fits this edge: the wall folds, and the fold band could not \
379 be carved — {carve_error}",
380 crate::blend::WALL_FOLDS,
381 radius.abs()
382 ),
383 )
384 }
385 })
386 }
387 Err(error) => Err(error),
388 }
389}
390
391fn carve_ladder(
401 solid: &BrepSolid,
402 segments: &[ChainSegment<'_>],
403 radius: f64,
404 name: Option<&str>,
405 bar: f64,
406) -> Result<BrepSolid, KernelRefusal> {
407 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
408 loop {
409 let per_segment = budget.rung;
410 match build_closed_smooth_chain(
411 solid,
412 segments,
413 radius,
414 false,
415 name,
416 &budget,
417 FoldPolicy::Carve,
418 bar,
419 )? {
420 Rung::Built(built) => {
421 crate::blend::carve::carve_trace(format_args!(
422 "blend carve: {} is the first rung inside \
423 {bar:.3e}",
424 budget.describe()
425 ));
426 return Ok(built);
427 }
428 Rung::StationFailed(failure) => return Err(failure.error),
431 Rung::TooCoarse { deviation } => {
432 crate::blend::carve::carve_trace(format_args!(
433 "blend carve: {} leaves the rails \
434 {deviation:.6e} from the rolling ball, against {bar:.3e}",
435 budget.describe()
436 ));
437 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
438 return Err(crate::blend::fold::wall_fold_refusal(
439 Vec::new(),
440 format!(
441 "{} {} fits this edge: the wall folds, and at {per_segment} stations \
442 per segment — the top of the carve's ladder — its fitted rails are \
443 still {deviation:.6e} from the rolling ball's own contacts against a \
444 bar of {bar:.3e}, so there is no wall accurate enough to carve",
445 crate::blend::WALL_FOLDS,
446 radius.abs()
447 ),
448 ));
449 }
450 budget = budget.doubled();
451 }
452 }
453 }
454}
455
456enum Rung {
459 Built(BrepSolid),
460 TooCoarse { deviation: f64 },
461 StationFailed(ChainStationFailure),
464}
465
466fn rail_deviation(
490 samples: &[ChainSample],
491 cr: &NurbsCurve,
492 cs: &NurbsCurve,
493 fit_low: f64,
494 fit_range: f64,
495 budget: &ChainBudget,
496 bar: f64,
497) -> Result<(f64, Vec<usize>), KernelRefusal> {
498 let mut order: Vec<usize> = (0..samples.len())
502 .filter(|&index| (0..=budget.count(samples[index].segment) as isize).contains(&samples[index].position))
503 .collect();
504 order.sort_by(|&a, &b| {
505 let (a, b) = (&samples[a], &samples[b]);
506 a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
507 });
508 let next_of: std::collections::HashMap<usize, f64> = order
509 .windows(2)
510 .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
511 .map(|pair| (pair[0], samples[pair[1]].parameter))
512 .collect();
513 let mut failing = Vec::new();
514 let closest = closest_on_rail;
515 let mut worst: f64 = 0.0;
516 let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
519 for &index in &order {
520 let sample = &samples[index];
521 let Some([p1, p2]) = sample.midpoint else {
522 continue;
523 };
524 let Some(&next) = next_of.get(&index) else {
525 continue;
526 };
527 let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
528 let seed = (global - fit_low) / fit_range;
529 let reach = (next - sample.parameter).abs() / fit_range;
530 let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
531 worst = worst.max(miss1).max(miss2);
532 if miss1.max(miss2) > bar {
533 failing.push(index);
534 }
535 if per_segment_worst.len() <= sample.segment {
536 per_segment_worst.resize(
537 sample.segment + 1,
538 (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
539 );
540 }
541 if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
542 per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
543 }
544 }
545 for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
546 crate::blend::carve::carve_trace(format_args!(
547 " rail miss segment {segment}: worst {miss:.6e} after station {position} \
548 (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
549 p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
550 ));
551 }
552 Ok((worst, failing))
553}
554
555const PIECE_TRACK_CEILING: usize = 4096;
559
560fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
563 let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
564 let (u, _) = crate::projection::curve_foot_newton(curve, target, seed, [low, high],
565 crate::projection::CurveFootPolicy {
566 domain: crate::projection::FootDomain::Clamp,
567 iterations: 40,
568 slope_guard: -1.0, step_limit: reach,
570 parameter_tolerance: 1e-15,
571 mode: crate::projection::FootMode::Backtracked,
572 }).or_refuse(KernelStage::Refine, "curve_foot_newton")?;
573 Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
574}
575
576pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
579
580fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
584 (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
585}
586
587#[allow(clippy::too_many_arguments)]
593pub(super) fn dense_rail_misses(
594 contacts: &crate::blend::chain::DenseContacts,
595 samples: &[ChainSample],
596 cr: &NurbsCurve,
597 cs: &NurbsCurve,
598 low: f64,
599 range: f64,
600 wrap: bool,
601) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
602 let mut readable = Vec::with_capacity(contacts.solved.len());
603 let mut unreadable = Vec::new();
604 for (left, next, fraction, p1, p2) in &contacts.solved {
605 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
606 let global = from + fraction * (to - from);
610 let global = if wrap { global.rem_euclid(1.0) } else { global };
611 let seed = (global - low) / range;
612 let reach = (to - from).abs() / range;
613 let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
616 let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
617 match rail_pair_miss(cr_miss, cs_miss) {
618 Some(miss) => readable.push((*left, miss)),
619 None => unreadable.push(*left),
620 }
621 }
622 Ok((readable, unreadable))
623}
624
625pub(super) fn combine_rail_verdict(
631 halfway: (f64, Vec<usize>),
632 dense: &[(usize, f64)],
633 unsolved: &[usize],
634 bar: f64,
635) -> (f64, Vec<usize>) {
636 let (mut worst, mut failing) = halfway;
637 for &(left, miss) in dense {
638 if !miss.is_finite() {
639 failing.push(left);
640 continue;
641 }
642 if miss > bar {
643 failing.push(left);
644 }
645 worst = worst.max(miss);
646 }
647 failing.extend_from_slice(unsolved);
648 failing.sort_unstable();
649 failing.dedup();
650 (worst, failing)
651}
652
653pub(super) fn model_rounds_limit() -> usize {
660 REFINE_ROUNDS
661}
662
663pub(super) fn model_request_on() -> bool {
664 wall_verdict_on()
665}
666
667
668pub(super) fn wall_verdict_on() -> bool {
671 true
672}
673
674
675
676
677
678pub(super) const WALL_SECTION_PROBES: usize = 8;
681
682pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
703 let j = locate_foot(p, nodes, j, radius)?;
704 sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
705}
706
707const FOOT_TRAVEL: usize = 4;
710
711fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
721 locate_foot_read(p, nodes, j, radius).map(|(node, _)| node)
722}
723
724fn locate_foot_read(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(usize, f64)> {
727 locate_foot_read_prepared(p, nodes, j, radius, None)
728}
729
730fn locate_foot_read_prepared(p: Vec3, nodes: &[Vec3], j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Option<(usize, f64)> {
731 let n = nodes.len();
732 if n < 4 || j == 0 || j + 1 >= n {
733 return None;
734 }
735 let mut tried: Vec<usize> = Vec::new();
736 let mut at = j;
737 loop {
738 let base = (at - 1).min(n - 4);
739 match stencil_solve_prepared(p, nodes, base, at, radius, prepared) {
740 Solve::Read(distance) => return Some((at, distance)),
741 Solve::Unread => return None,
742 Solve::Window(s) => {
743 tried.push(at);
744 let next = (base as f64 + s).round();
745 let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
746 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
747 wall_decline(|| format!(
748 "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
749 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
750 ));
751 return None;
752 }
753 at = next;
754 }
755 }
756 }
757}
758
759fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
762 match stencil_solve(p, nodes, base, j, radius) {
763 Solve::Read(distance) => Some(distance),
764 Solve::Window(_) | Solve::Unread => None,
765 }
766}
767
768enum Solve {
771 Read(f64),
772 Window(f64),
773 Unread,
774}
775
776thread_local! {
777 static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
782}
783
784fn wall_decline(reason: impl FnOnce() -> String) {
785 WALL_DECLINES.with(|log| {
786 if let Some(lines) = log.borrow_mut().as_mut() {
787 lines.push(reason());
788 }
789 });
790}
791
792pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
800 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
801 let read = declared_wall_probe(p, fine, jf, radius);
802 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
803 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
804}
805
806fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
807 stencil_solve_prepared(p, nodes, base, j, radius, None)
808}
809
810fn stencil_solve_prepared(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Solve {
811 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
812 let n = nodes.len();
813 if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || !radius.is_finite() || !finite(p) || !prepared.map_or_else(|| wall_nodes_finite(nodes), |cache| cache.finite) {
814 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
815 return Solve::Unread;
816 }
817 let [p0, d1, d2, d3] = prepared.map_or_else(|| wall_coefficients(nodes, base), |cache| cache.coefficients(nodes, base));
818 let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
819 let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
820 let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
821 let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
822 (c, c1, c2)
823 };
824 let start = j as f64 - base as f64;
825 let initial = prepared.and_then(|cache| cache.start(base, j, &curve));
826 let point_length = p.length();
827 let mut s = start;
828 let mut converged = false;
829 let mut converged_point = None;
830 for iteration in 0..50 {
831 let state = if iteration == 0 { initial.unwrap_or_else(|| WallStart::new(curve(s))) }
832 else { WallStart::new(curve(s)) };
833 let (c, c1, c2) = state.curve;
834 let w = c.sub(p);
835 let (g, slope) = (w.dot(c1), state.tangent_squared + w.dot(c2));
836 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
837 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
838 return Solve::Unread;
839 }
840 let stationary = g.abs() <= 16.0 * f64::EPSILON * (point_length + state.point_length) * state.tangent_length;
853 if stationary {
854 converged = true;
855 converged_point = Some(c);
856 break;
857 }
858 let step = g / slope;
859 let before = s;
860 s -= step;
861 if !((s - start).abs() <= 1.0) {
862 wall_decline(|| {
863 let (c, c1, _) = curve(s);
866 let g_after = c.sub(p).dot(c1);
867 let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
868 format!(
869 "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
870 g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
871 (s - start).abs() - 1.0
872 )
873 });
874 return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
875 }
876 if step.abs() <= 1e-14 * (1.0 + s.abs()) {
877 converged = true;
878 break;
879 }
880 }
881 if !converged {
882 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
883 return Solve::Unread;
884 }
885 let distance = converged_point.unwrap_or_else(|| curve(s).0).sub(p).length() - radius;
886 if distance.is_finite() {
887 Solve::Read(distance)
888 } else {
889 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
890 Solve::Unread
891 }
892}
893
894pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
903 let n = nodes.len();
904 if n < 5 || j == 0 || j + 1 >= n {
905 return None;
906 }
907 let j = locate_foot(p, nodes, j, radius)?;
908 let first = (j - 1).min(n - 4);
909 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
910 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
911 let a = sweep_residual_on(p, nodes, first, j, radius)?;
912 let b = sweep_residual_on(p, nodes, second, j, radius)?;
913 Some((a.abs().max(b.abs()), (a - b).abs()))
914}
915
916fn wall_solve_reuse_on() -> bool {
929 switch!(live "BREP_WALL_SOLVE_REUSE", true)
930}
931
932pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
933 declared_wall_probe_with(p, fine, jf, radius, wall_solve_reuse_on())
934}
935
936fn declared_wall_probe_with(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool) -> Option<(f64, f64)> {
937 declared_wall_probe_prepared(p, fine, jf, radius, reuse, None)
938}
939
940fn declared_wall_probe_prepared(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool, prepared: Option<&PreparedWall<'_>>) -> Option<(f64, f64)> {
941 if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
942 wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
943 return None;
944 }
945 let (reading, spread, at_fine) = if reuse {
946 let read = (|| {
947 let n = fine.len();
948 let (j, a) = locate_foot_read_prepared(p, fine, jf, radius, prepared.map(|c| &c.fine_stencils))?;
949 let first = (j - 1).min(n - 4);
950 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
951 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
952 let b = match stencil_solve_prepared(p, fine, second, j, radius, prepared.map(|c| &c.fine_stencils)) {
953 Solve::Read(distance) => distance,
954 _ => return None,
955 };
956 Some((a.abs().max(b.abs()), (a - b).abs(), a))
957 })();
958 let Some(read) = read else {
959 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
960 return None;
961 };
962 read
963 } else {
964 let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
965 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
966 return None;
967 };
968 let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
969 wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
970 return None;
971 };
972 (reading, spread, at_fine)
973 };
974 let coarse_owned;
975 let coarse = match prepared {
976 Some(cache) => cache.coarse.as_slice(),
977 None => { coarse_owned = wall_coarse(fine); &coarse_owned }
978 };
979 let at_coarse = if reuse {
980 locate_foot_read_prepared(p, coarse, jf / 2, radius, prepared.map(|c| &c.coarse_stencils)).map(|(_, distance)| distance)
981 } else {
982 sweep_residual(p, coarse, jf / 2, radius)
983 };
984 let Some(at_coarse) = at_coarse else {
985 wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
986 return None;
987 };
988 let uncertainty = spread.max((at_fine - at_coarse).abs());
989 uncertainty.is_finite().then_some((reading, uncertainty))
990}
991
992pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
1008 let n = fine.len();
1009 if n < 7 || !(fraction.is_finite()) {
1010 return None;
1011 }
1012 let mut cumulative = Vec::with_capacity(n);
1013 cumulative.push(0.0_f64);
1014 for pair in fine.windows(2) {
1015 let step = pair[1].sub(pair[0]).length();
1016 cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
1017 }
1018 let total = *cumulative.last()?;
1019 if !(total.is_finite() && total > 0.0) {
1020 return None;
1021 }
1022 let target = fraction.clamp(0.0, 1.0) * total;
1023 let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
1024 let even = 2 * ((nearest as f64 / 2.0).round() as usize);
1025 Some(even.clamp(2, (n - 3) & !1))
1026}
1027
1028pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
1035 let first = match declared_wall_probe_why(p, fine, jf, radius) {
1036 Ok(read) => return Ok(read),
1037 Err(why) => why,
1038 };
1039 match centre_path_node(fine, fraction) {
1040 Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
1041 .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
1042 _ => Err(first),
1043 }
1044}
1045
1046pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
1048 match declared_wall_probe(p, fine, jf, radius) {
1049 Some(read) => Some(read),
1050 None => match centre_path_node(fine, fraction) {
1051 Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
1052 _ => None,
1053 },
1054 }
1055}
1056
1057pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
1069 if jf % 2 != 0
1070 || fine.len() < 9
1071 || fine.len() % 2 == 0
1072 || ts.len() != fine.len()
1073 || !ts.iter().all(|t| t.is_finite())
1074 || !ts.windows(2).all(|pair| pair[1] > pair[0])
1075 {
1076 return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
1077 }
1078 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
1079 let read = (|| -> Option<(f64, f64)> {
1080 let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
1081 wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
1082 None
1083 })?;
1084 let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
1085 wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
1086 None
1087 })?;
1088 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
1089 let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
1090 let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
1091 wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
1092 None
1093 })?;
1094 let uncertainty = spread.max((at_fine - at_coarse).abs());
1095 uncertainty.is_finite().then_some((reading, uncertainty))
1096 })();
1097 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
1098 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
1099}
1100
1101fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
1102 let n = nodes.len();
1103 if n < 5 || j == 0 || j + 1 >= n {
1104 return None;
1105 }
1106 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1107 let first = (j - 1).min(n - 4);
1108 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
1109 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
1110 let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
1111 let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
1112 Some((a.abs().max(b.abs()), (a - b).abs()))
1113}
1114
1115fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
1116 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1117 match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
1118 Solve::Read(distance) => Some(distance),
1119 Solve::Window(_) | Solve::Unread => None,
1120 }
1121}
1122
1123fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
1127 let n = nodes.len();
1128 if n < 4 || j == 0 || j + 1 >= n {
1129 return None;
1130 }
1131 let mut tried: Vec<usize> = Vec::new();
1132 let mut at = j;
1133 loop {
1134 let base = (at - 1).min(n - 4);
1135 match stencil_solve_at(p, nodes, ts, base, at, radius) {
1136 Solve::Read(_) => return Some(at),
1137 Solve::Unread => return None,
1138 Solve::Window(t) => {
1139 tried.push(at);
1140 let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
1141 let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
1142 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
1143 wall_decline(|| format!(
1144 "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
1145 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
1146 ));
1147 return None;
1148 }
1149 at = next;
1150 }
1151 }
1152 }
1153}
1154
1155fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
1160 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
1161 let n = nodes.len();
1162 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)) {
1163 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
1164 return Solve::Unread;
1165 }
1166 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
1167 let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
1168 let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
1169 let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
1170 let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
1171 let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
1172 let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
1173 let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
1174 let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
1175 let (a, b, c) = (t - t0, t - t1, t - t2);
1176 let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
1177 let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
1178 let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
1179 (point, first, second)
1180 };
1181 let (low, high) = (ts[j - 1], ts[j + 1]);
1182 let start = ts[j];
1183 let mut t = start;
1184 let mut converged = false;
1185 for _ in 0..50 {
1186 let (c, c1, c2) = curve(t);
1187 let w = c.sub(p);
1188 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
1189 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
1190 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
1191 return Solve::Unread;
1192 }
1193 if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
1194 converged = true;
1195 break;
1196 }
1197 let step = g / slope;
1198 let before = t;
1199 t -= step;
1200 if !(t >= low && t <= high) {
1201 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
1202 return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
1203 }
1204 let node = 0.5 * (high - low);
1209 if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
1210 converged = true;
1211 break;
1212 }
1213 }
1214 if !converged {
1215 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
1216 return Solve::Unread;
1217 }
1218 let distance = curve(t).0.sub(p).length() - radius;
1219 if distance.is_finite() {
1220 Solve::Read(distance)
1221 } else {
1222 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
1223 Solve::Unread
1224 }
1225}
1226
1227pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
1232 readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
1233}
1234
1235#[allow(clippy::too_many_arguments)]
1245pub(super) fn wall_interior_misses(
1246 contacts: &crate::blend::chain::DenseContacts,
1247 samples: &[ChainSample],
1248 surface: &crate::NurbsSurface,
1249 low: f64,
1250 range: f64,
1251 radius: f64,
1252 subdivisions: usize,
1253 probes: usize,
1254 wrap: bool,
1255) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1256 use std::collections::BTreeMap;
1257 let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1260 for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1261 let entry = nodes.entry(*left).or_insert_with(|| {
1262 let mut row = vec![None; subdivisions + 1];
1263 row[0] = Some(samples[*left].station.center);
1264 row[subdivisions] = Some(samples[*next].station.center);
1265 (*next, row)
1266 });
1267 let step = (fraction * subdivisions as f64).round() as usize;
1268 if (1..subdivisions).contains(&step) {
1269 entry.1[step] = contacts.centers.get(index).copied().flatten();
1270 }
1271 }
1272 let mut readable = Vec::new();
1273 let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1274 for (&left, (next, row)) in &nodes {
1275 let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1276 unreadable.push(left);
1277 continue;
1278 };
1279 let (from, to) = (samples[left].parameter, samples[*next].parameter);
1280 let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1281 let mut read = true;
1282 let prepared = PreparedWall::for_interval(¢res);
1283 'probes: for step in (2..subdivisions).step_by(2) {
1284 let global = from + (step as f64 / subdivisions as f64) * (to - from);
1285 let global = if wrap { global.rem_euclid(1.0) } else { global };
1286 let u = (global - low) / range;
1287 for k in 1..probes {
1288 let v = k as f64 / probes as f64;
1289 let reading = surface.evaluate(u, v).ok().and_then(|p| {
1290 match &prepared {
1291 Some(cache) => cache.read(p, step, step as f64 / subdivisions as f64, radius),
1292 None => declared_wall_probe_located(p, ¢res, step, step as f64 / subdivisions as f64, radius),
1293 }
1294 });
1295 match reading {
1296 Some((distance, gap)) => {
1297 declared = declared.max(distance + gap);
1298 worst = worst.max(distance);
1299 spread = spread.max(gap);
1300 }
1301 None => {
1302 read = false;
1303 break 'probes;
1304 }
1305 }
1306 }
1307 }
1308 if read {
1309 readable.push((left, declared, worst, spread));
1310 } else {
1311 unreadable.push(left);
1312 }
1313 }
1314 unreadable.sort_unstable();
1315 unreadable.dedup();
1316 (readable, unreadable)
1317}
1318
1319
1320fn build_closed_smooth_chain(
1321 solid: &BrepSolid,
1322 segments: &[ChainSegment<'_>],
1323 radius: f64,
1324 chamfer: bool,
1325 name: Option<&str>,
1326 budget: &ChainBudget,
1327 fold: FoldPolicy,
1328 bar: f64,
1329) -> Result<Rung, KernelRefusal> {
1330 let samples = match march_chain_typed(segments, radius, false, budget, fold) {
1331 Ok(samples) => samples,
1332 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1333 Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1336 Err(ChainMarchError::Other(error)) => return Err(error),
1337 };
1338 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1339 build_closed_chain_from_samples(solid, segments, radius, chamfer, name, budget, fold, bar, samples, 0)
1340}
1341
1342
1343thread_local! {
1344 static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1348}
1349
1350fn local_continue_on() -> bool {
1356 switch!(live "BREP_CHAIN_LOCAL_CONTINUE", false)
1357}
1358
1359const LOCAL_CONTINUE_MAX_ROUNDS: usize = 4 * REFINE_ROUNDS;
1361
1362const LOCAL_CONTINUE_STALL_ROUNDS: usize = 3;
1364
1365thread_local! {
1366 static CHAIN_LOCAL_BEST: std::cell::Cell<(f64, usize)> = const { std::cell::Cell::new((f64::INFINITY, 0)) };
1368}
1369
1370fn local_rounds_spent(round: usize, deviation: f64) -> bool {
1376 if round == 0 {
1377 CHAIN_LOCAL_BEST.with(|best| best.set((f64::INFINITY, 0)));
1378 return false;
1379 }
1380 let (_, best_round) = CHAIN_LOCAL_BEST.with(|best| {
1381 let (reading, at) = best.get();
1382 if deviation < reading {
1383 best.set((deviation, round));
1384 (deviation, round)
1385 } else {
1386 (reading, at)
1387 }
1388 });
1389 round >= LOCAL_CONTINUE_MAX_ROUNDS || round - best_round >= LOCAL_CONTINUE_STALL_ROUNDS
1390}
1391#[allow(clippy::too_many_arguments)]
1402fn build_closed_chain_from_samples(
1403 solid: &BrepSolid,
1404 segments: &[ChainSegment<'_>],
1405 radius: f64,
1406 chamfer: bool,
1407 name: Option<&str>,
1408 budget: &ChainBudget,
1409 fold: FoldPolicy,
1410 bar: f64,
1411 mut samples: Vec<ChainSample>,
1412 round: usize,
1413) -> Result<Rung, KernelRefusal> {
1414 let per_segment = budget.rung;
1417 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1418
1419 let mut global: Vec<(f64, &ChainSample)> = samples
1421 .iter()
1422 .filter(|sample| (0..budget.count(sample.segment) as isize).contains(&sample.position))
1423 .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1424 .collect();
1425 global.sort_by(|a, b| a.0.total_cmp(&b.0));
1426 let count = global.len();
1427 let degree = FIT_DEGREE;
1428 let mut extended_params = Vec::new();
1429 let mut samples_cr = Vec::new();
1430 let mut samples_mid = Vec::new();
1431 let mut samples_cs = Vec::new();
1432 let mut push = |station: &Station, parameter: f64| {
1433 extended_params.push(parameter);
1434 samples_cr.push(Vec4::from_point(station.p1, 1.0));
1435 samples_cs.push(Vec4::from_point(station.p2, 1.0));
1436 samples_mid.push(Vec4 {
1437 x: station.apex.x * station.weight,
1438 y: station.apex.y * station.weight,
1439 z: station.apex.z * station.weight,
1440 w: station.weight,
1441 });
1442 };
1443 for offset in (1..=degree).rev() {
1444 let (parameter, sample) = &global[count - offset];
1445 push(&sample.station, parameter - 1.0);
1446 }
1447 for (parameter, sample) in &global {
1448 push(&sample.station, *parameter);
1449 }
1450 push(&global[0].1.station, global[0].0 + 1.0);
1453 for offset in 1..=degree {
1454 let (parameter, sample) = &global[offset];
1455 push(&sample.station, parameter + 1.0);
1456 }
1457 let low = extended_params[0];
1458 let high = *extended_params.last().unwrap();
1459 let range = high - low;
1460 let normalized: Vec<f64> = extended_params
1461 .iter()
1462 .map(|parameter| (parameter - low) / range)
1463 .collect();
1464 let seam_low = (0.0 - low) / range;
1465 let seam_high = (1.0 - low) / range;
1466 let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1467 let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1468 let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1469 let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1470 Ok(middle)
1471 };
1472 let cr = fit_row(&samples_cr)?;
1473 let cs = fit_row(&samples_cs)?;
1474 let mid = if chamfer {
1475 None
1476 } else {
1477 Some(fit_row(&samples_mid)?)
1478 };
1479 let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1480 let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1481
1482 let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, budget, bar)?;
1495 let stencil = if round > 0 {
1507 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, LOCAL_RAIL_STENCIL)?)
1508 } else {
1509 None
1510 };
1511 let wall_stencil = if !chamfer && wall_verdict_on() {
1516 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1517 } else {
1518 None
1519 };
1520 let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1521 let contacts = stencil.as_ref().expect("the local path solves the stencil");
1522 let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1523 let reads = dense.len() + unreadable.len() + contacts.unsolved;
1524 let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1525 let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1526 crate::blend::carve::carve_trace(format_args!(
1527 "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1528 {} failing interval(s), {} unsolved contact(s)",
1529 verdict.0, verdict.1.len(), contacts.unsolved
1530 ));
1531 (verdict.0, verdict.1, reads)
1532 } else {
1533 (halfway, halfway_failing, 0)
1534 };
1535 let mut _wall_worst = 0.0_f64;
1541 let mut _wall_spread = 0.0_f64;
1542 let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1543 if let Some(contacts) = wall_stencil.as_ref() {
1544 let (wall, unreadable) =
1545 wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1546 wall_readings = wall.clone();
1547 for &(left, declared, _, uncertainty) in &wall {
1549 if declared > bar {
1550 failing.push(left);
1551 }
1552 _wall_worst = _wall_worst.max(declared);
1553 _wall_spread = _wall_spread.max(uncertainty);
1554 }
1555 failing.extend_from_slice(&unreadable);
1556 failing.sort_unstable();
1557 failing.dedup();
1558 deviation = deviation.max(_wall_worst);
1559 crate::blend::carve::carve_trace(format_args!(
1560 "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1561 interval(s), {} unreadable, bar {bar:.3e}",
1562 wall.len(),
1563 unreadable.len()
1564 ));
1565 }
1566 if deviation > bar || !failing.is_empty() {
1567 let in_segment = |segment: usize| {
1568 samples
1569 .iter()
1570 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1571 .count()
1572 };
1573 let within_ceiling = (0..segments.len()).all(|segment| {
1574 in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1575 <= CHAIN_CARVE_MAX_PER_SEGMENT
1576 });
1577 crate::blend::carve::carve_trace(format_args!(
1578 "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1579 worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1580 samples.len(),
1581 failing.len(),
1582 if within_ceiling { "ok" } else { "reached" }
1583 ));
1584 let rounds_spent = if local_continue_on() {
1586 local_rounds_spent(round, deviation)
1587 } else {
1588 round >= REFINE_ROUNDS
1589 };
1590 if fold != FoldPolicy::Refuse || rounds_spent || !within_ceiling {
1592 return Ok(Rung::TooCoarse { deviation });
1593 }
1594 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1596 match insert_chain_stations(segments, radius, budget, &mut samples, &failing) {
1597 Ok(()) => {}
1598 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1599 Err(ChainMarchError::Incoherent(failure)) => {
1603 crate::blend::carve::carve_trace(format_args!(
1604 "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1605 at t {:.9} does not reach its retained right neighbour: {}",
1606 failure.segment, failure.t, failure.error
1607 ));
1608 return Ok(Rung::TooCoarse { deviation });
1609 }
1610 Err(ChainMarchError::Other(error)) => return Err(error),
1611 }
1612 assign_chain_parameters(&mut samples, segments.len(), budget, false);
1613 return build_closed_chain_from_samples(
1614 solid, segments, radius, chamfer, name, budget, fold, bar, samples, round + 1,
1615 );
1616 }
1617 crate::blend::carve::carve_trace(format_args!(
1618 "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1619 from the rolling ball, inside {bar:.3e}"
1620 ));
1621 let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1635 let stations_now = (0..segments.len())
1636 .map(|segment| {
1637 samples
1638 .iter()
1639 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1640 .count()
1641 })
1642 .max()
1643 .unwrap_or(0);
1644 let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1645 let mut request_detail: Option<String> = None;
1649 let mut request_component = crate::MeasuredComponent::Wall;
1650 let mut request_unread = 0usize;
1651 if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1652 let wall_missing = wall_model_missing(&wall_readings, model);
1653 let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1661 Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1662 None => (Vec::new(), Vec::new()),
1663 };
1664 let rail_over: std::collections::BTreeSet<usize> =
1666 rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1667 let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1668 let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1669 let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1670 missing.sort_unstable();
1671 missing.dedup();
1672 if !missing.is_empty() {
1673 let wall_short = wall_readings
1675 .iter()
1676 .filter(|reading| wall_missing.contains(&reading.0))
1677 .map(|&(_, declared, _, _)| declared)
1678 .fold(0.0_f64, f64::max);
1679 let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1680 let short = wall_short.max(rail_short);
1681 request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1682 request_unread = rail_unread_intervals.len();
1683 request_detail = Some(format!(
1684 "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1685 wall {wall_short:.3e} over {} interval(s)",
1686 rail_over.len(),
1687 rail_unread_intervals.len(),
1688 wall_missing.len()
1689 ));
1690 let within_ceiling = (0..segments.len()).all(|segment| {
1691 let present = samples
1692 .iter()
1693 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1694 .count();
1695 present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1696 <= CHAIN_CARVE_MAX_PER_SEGMENT
1697 });
1698 crate::blend::carve::carve_trace(format_args!(
1699 "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1700 (worst reading + spread {short:.3e}), ceiling {}",
1701 missing.len(),
1702 if within_ceiling { "ok" } else { "reached" }
1703 ));
1704 unmet = if round >= model_rounds_limit() {
1705 Some((short, crate::BudgetReason::RoundsSpent))
1706 } else if !within_ceiling {
1707 Some((short, crate::BudgetReason::StationCeiling))
1708 } else {
1709 let mut refined = samples.clone();
1710 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1712 match insert_chain_stations(segments, radius, budget, &mut refined, &missing) {
1713 Ok(()) => {
1714 assign_chain_parameters(&mut refined, segments.len(), budget, false);
1715 match build_closed_chain_from_samples(
1716 solid, segments, radius, chamfer, name, budget, fold, bar, refined, round + 1,
1717 ) {
1718 Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1719 _ => Some((short, crate::BudgetReason::Rejected)),
1720 }
1721 }
1722 Err(_) => {
1723 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1725 Some((short, crate::BudgetReason::Incoherent))
1726 }
1727 }
1728 };
1729 }
1730 }
1731 if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1736 let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, subdivisions.max(2))?;
1737 let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1738 let mut unreadable = 0usize;
1739 for (left, next, fraction, p1, p2) in &contacts.solved {
1740 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1741 let global = (from + fraction * (to - from)).rem_euclid(1.0);
1742 let seed = (global - low) / range;
1743 let reach = (to - from).abs() / range;
1744 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)) {
1745 Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1746 Some(_) => {}
1747 None => unreadable += 1,
1748 }
1749 }
1750 eprintln!(
1751 "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1752 worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1753 contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1754 );
1755 }
1756 let crease = if fold == FoldPolicy::Carve {
1757 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1767 let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1768 return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1769 "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"
1770 )));
1771 };
1772 Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1773 } else {
1774 None
1775 };
1776
1777 let single_face = |side: usize| -> bool {
1781 let first_id = segments[0].mate(side).face.id;
1782 segments
1783 .iter()
1784 .all(|segment| segment.mate(side).face.id == first_id)
1785 };
1786 let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1787 for side in 0..2 {
1788 if !single_face(side) {
1789 continue;
1790 }
1791 let mut row = Vec::new();
1792 let select = |station: &Station| -> [f64; 2] {
1793 if side == 0 {
1794 station.uv1
1795 } else {
1796 station.uv2
1797 }
1798 };
1799 let mut push_uv = |station: &Station| {
1800 let uv = select(station);
1801 row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1802 };
1803 for offset in (1..=degree).rev() {
1804 push_uv(&global[count - offset].1.station);
1805 }
1806 for (_, sample) in &global {
1807 push_uv(&sample.station);
1808 }
1809 push_uv(&global[0].1.station);
1810 for offset in 1..=degree {
1811 push_uv(&global[offset].1.station);
1812 }
1813 whole_pcurves[side] = Some(fit_row(&row)?);
1814 }
1815
1816 let mut pcurves1 = Vec::with_capacity(segments.len());
1819 let mut pcurves2 = Vec::with_capacity(segments.len());
1820 for segment in 0..segments.len() {
1821 let mut window: Vec<&ChainSample> = samples
1822 .iter()
1823 .filter(|sample| sample.segment == segment)
1824 .collect();
1825 window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1826 let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1827 let low = params[0];
1828 let high = *params.last().unwrap();
1829 let normalized: Vec<f64> = params
1830 .iter()
1831 .map(|parameter| (parameter - low) / (high - low))
1832 .collect();
1833 let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1834 let points: Vec<Vec4> = window
1835 .iter()
1836 .map(|sample| {
1837 let uv = select(&sample.station);
1838 Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1839 })
1840 .collect();
1841 let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1842 Ok(curve)
1846 };
1847 pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1848 pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1849 }
1850
1851 let wall_surface = unmet.as_ref().map(|_| surface.clone());
1852 let built = chain_surgery(
1853 solid,
1854 segments,
1855 ChainRows {
1856 surface,
1857 cr,
1858 cs,
1859 u_domain,
1860 fit_low: low,
1861 fit_range: range,
1862 pcurves1,
1863 pcurves2,
1864 whole_pcurves,
1865 crease,
1866 },
1867 name,
1868 )?;
1869 if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1870 crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1871 name: name.map(str::to_string),
1872 surface,
1873 request: model,
1874 residual,
1875 detail: request_detail,
1876 budget: crate::ApproximationBudget {
1877 reason,
1878 rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1881 rounds_limit: model_rounds_limit(),
1882 stations: stations_now,
1883 station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1884 mechanism: crate::BudgetMechanism::LocalRounds,
1885 measured_component: request_component,
1886 unread: request_unread,
1887 },
1888 });
1889 }
1890 Ok(Rung::Built(built))
1891}
1892
1893pub(super) struct ChainRows {
1894 pub(super) surface: NurbsSurface,
1895 pub(super) cr: NurbsCurve,
1896 pub(super) cs: NurbsCurve,
1897 pub(super) u_domain: [f64; 2],
1898 pub(super) fit_low: f64,
1902 pub(super) fit_range: f64,
1903 pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1906 pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1907 pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1909 pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1912}
1913
1914impl ChainRows {
1915 pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1916 self.fit_low + fit_parameter * self.fit_range
1917 }
1918}
1919
1920pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1924 let (low, high, curve) = fitted;
1925 let to_local = |value: f64| {
1928 let mut best = value;
1929 for candidate in [value, value + 1.0, value - 1.0] {
1930 if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1931 best = candidate;
1932 break;
1933 }
1934 }
1935 ((best - low) / (high - low)).clamp(0.0, 1.0)
1936 };
1937 let mut local_a = to_local(a);
1938 let mut local_b = to_local(b);
1939 if local_a > local_b {
1940 std::mem::swap(&mut local_a, &mut local_b);
1941 }
1942 let epsilon = 1e-9;
1943 let (_, tail) = if local_a > epsilon {
1944 curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1945 } else {
1946 (curve.clone(), curve.clone())
1947 };
1948 let tail = if local_a > epsilon {
1949 tail
1950 } else {
1951 curve.clone()
1952 };
1953 let portion = if local_b < 1.0 - epsilon {
1954 tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1955 } else {
1956 tail
1957 };
1958 Ok(portion)
1959}
1960
1961pub(super) fn cross_edge_at(
1967 solid: &BrepSolid,
1968 face_before: &FaceRecord,
1969 loop_before: usize,
1970 face_after: &FaceRecord,
1971 loop_after: usize,
1972 vertex: u64,
1973 before_edge: u64,
1974 after_edge: u64,
1975) -> Result<Option<u64>, KernelRefusal> {
1976 let mut found: Option<u64> = None;
1977 let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1978 for coedge in &face.loops[loop_index].coedges {
1979 if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1980 continue;
1981 }
1982 let edge = solid
1983 .edges
1984 .iter()
1985 .find(|edge| edge.id == coedge.edge_id)
1986 .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1987 if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1988 if found.is_some() && found != Some(edge.id) {
1989 return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1990 }
1991 found = Some(edge.id);
1992 }
1993 }
1994 Ok(())
1995 };
1996 consider(face_before, loop_before)?;
1997 if !(face_after.id == face_before.id && loop_after == loop_before) {
1998 consider(face_after, loop_after)?;
1999 }
2000 Ok(found)
2001}
2002
2003fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
2006 let period = high - low;
2007 low + ((value - low) / period).round() * period
2008}
2009
2010pub(in crate::blend) fn project_piece_pcurve(
2028 piece: &NurbsCurve,
2029 surface: &NurbsSurface,
2030 start_is_crossing: bool,
2031 end_is_crossing: bool,
2032) -> Result<NurbsCurve, KernelRefusal> {
2033 project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
2034}
2035
2036
2037pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
2042 let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
2043 let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2044 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2045 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2046 let into = move |value: f64, closed: bool, low: f64, high: f64| {
2047 if closed && high > low {
2048 low + (value - low).rem_euclid(high - low)
2049 } else {
2050 value.clamp(low, high)
2051 }
2052 };
2053 Ok(move |fraction: f64| {
2054 let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
2055 let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
2056 seed.iter().all(|value| value.is_finite()).then_some(seed)
2057 })
2058}
2059
2060pub(in crate::blend) fn project_piece_pcurve_from(
2070 piece: &NurbsCurve,
2071 surface: &NurbsSurface,
2072 start_is_crossing: bool,
2073 end_is_crossing: bool,
2074 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2075) -> Result<NurbsCurve, KernelRefusal> {
2076 project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
2077}
2078
2079pub(crate) const CHAIN_PCURVE_FLOOR_WHAT: &str = "chain_pcurve_floor";
2085
2086pub(in crate::blend) fn project_piece_pcurve_within(
2087 piece: &NurbsCurve,
2088 surface: &NurbsSurface,
2089 start_is_crossing: bool,
2090 end_is_crossing: bool,
2091 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2092 refine_near: &[f64],
2093) -> Result<NurbsCurve, KernelRefusal> {
2094 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false)
2095}
2096
2097pub(in crate::blend) fn project_piece_pcurve_within_knot_retry(
2108 piece: &NurbsCurve,
2109 surface: &NurbsSurface,
2110 start_is_crossing: bool,
2111 end_is_crossing: bool,
2112 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2113 refine_near: &[f64],
2114) -> Result<NurbsCurve, KernelRefusal> {
2115 match project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false) {
2116 Err(uniform) if matches!(&uniform.class, crate::RefusalClass::NonConvergence { what } if what == CHAIN_PCURVE_FLOOR_WHAT) => {
2117 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, true).map_err(|_| uniform)
2118 }
2119 other => other,
2120 }
2121}
2122
2123#[allow(clippy::too_many_arguments)]
2127fn project_piece_pcurve_track(
2128 piece: &NurbsCurve,
2129 surface: &NurbsSurface,
2130 start_is_crossing: bool,
2131 end_is_crossing: bool,
2132 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2133 refine_near: &[f64],
2134 knot_aware: bool,
2135) -> Result<NurbsCurve, KernelRefusal> {
2136 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2137 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2138 let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2139 let period = u1 - u0;
2140 let period_v = v1 - v0;
2141 let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
2142 let ceiling = if knot_aware { 2 * PIECE_TRACK_CEILING } else { PIECE_TRACK_CEILING };
2152 let uniform = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
2153 let mut fractions: Vec<f64> = (0..=uniform).map(|section| section as f64 / uniform as f64).collect();
2154 if knot_aware {
2155 let span = d1 - d0;
2158 let mut knots: Vec<f64> = piece.knots.iter().copied().filter(|knot| *knot > d0 && *knot < d1).map(|knot| (knot - d0) / span).collect();
2159 knots.dedup();
2160 fractions.extend(knots);
2161 fractions.sort_by(|a, b| a.total_cmp(b));
2162 fractions.dedup_by(|b, a| (*b - *a).abs() <= 1e-6 / uniform as f64);
2163 }
2164 let dense = fractions.len() - 1;
2165 crate::blend::record_count("count.track.piece_samples", dense as u64 + 1);
2166 let mut proj_u = Vec::with_capacity(dense + 1);
2167 let mut proj_v = Vec::with_capacity(dense + 1);
2168 let mut point3 = Vec::with_capacity(dense + 1);
2169 let mut params = Vec::with_capacity(dense + 1);
2170 let mut previous_u: Option<f64> = None;
2171 let mut previous_foot: Option<[f64; 2]> = None;
2176 let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
2177 let project_timer = crate::blend::BlendTimer::start("timing.track.project_ms");
2178 for section in 0..=dense {
2179 let fraction = fractions[section];
2180 let t = if knot_aware { d0 + (d1 - d0) * fraction } else { d0 + (d1 - d0) * section as f64 / dense as f64 };
2184 let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
2185 let guide = guide.filter(|_| {
2186 true
2187 });
2188 let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(fraction)).into_iter().chain(previous_foot).collect();
2193 let mut projection: Option<crate::projection::SurfaceProjection> = None;
2194 for seed in &seeds {
2195 let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
2196 if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
2197 projection = Some(found);
2198 }
2199 }
2200 let projection = match projection {
2201 Some(projection) => projection,
2202 None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
2203 };
2204 previous_foot = Some([projection.u, projection.v]);
2205 feet.push([projection.u, projection.v]);
2206 let mut u = projection.u;
2207 if let Some(previous) = previous_u {
2208 while u - previous > period * 0.5 {
2209 u -= period;
2210 }
2211 while previous - u > period * 0.5 {
2212 u += period;
2213 }
2214 }
2215 previous_u = Some(u);
2216 proj_u.push(u);
2217 proj_v.push(projection.v);
2218 point3.push(point);
2219 params.push(fraction);
2220 }
2221 project_timer.stop();
2222 if closed_v {
2229 for index in 2..=dense {
2230 while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
2231 proj_v[index] -= period_v;
2232 }
2233 while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
2234 proj_v[index] += period_v;
2235 }
2236 }
2237 while proj_v[0] - proj_v[1] > period_v * 0.5 {
2238 proj_v[0] -= period_v;
2239 }
2240 while proj_v[1] - proj_v[0] > period_v * 0.5 {
2241 proj_v[0] += period_v;
2242 }
2243 }
2244 let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2247 let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
2248 let pinned = |index: usize,
2257 neighbour: usize,
2258 proj_u: &[f64],
2259 proj_v: &[f64]|
2260 -> Result<(f64, f64), KernelRefusal> {
2261 let point = point3[index];
2262 let u_error = surface
2263 .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
2264 .sub(point)
2265 .length();
2266 if closed_v {
2267 let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
2268 let v_error = surface
2269 .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
2270 .sub(point)
2271 .length();
2272 if v_error < u_error {
2273 return Ok((proj_u[index], seam_v));
2274 }
2275 }
2276 Ok((seam_u, proj_v[index]))
2277 };
2278 if start_is_crossing {
2279 let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
2280 proj_u[0] = u;
2281 proj_v[0] = v;
2282 }
2283 if end_is_crossing {
2284 let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
2285 proj_u[dense] = u;
2286 proj_v[dense] = v;
2287 }
2288 let mut shift = 0.0;
2289 while mean_u + shift > u1 + 1e-9 {
2290 shift -= period;
2291 }
2292 while mean_u + shift < u0 - 1e-9 {
2293 shift += period;
2294 }
2295 for u in proj_u.iter_mut() {
2296 *u += shift;
2297 }
2298 if closed_v {
2299 let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2300 let mut shift_v = 0.0;
2301 while mean_v + shift_v > v1 + 1e-9 {
2302 shift_v -= period_v;
2303 }
2304 while mean_v + shift_v < v0 - 1e-9 {
2305 shift_v += period_v;
2306 }
2307 for v in proj_v.iter_mut() {
2308 *v += shift_v;
2309 }
2310 }
2311 let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
2326 let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
2333 let mut track = crate::blend::track_fit::Track {
2334 fractions: params,
2335 uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
2336 feet,
2337 breaks: Vec::new(),
2338 };
2339 let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
2340 let polish_timer = crate::blend::BlendTimer::start("timing.track.polish_ms");
2341 crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
2342 polish_timer.stop();
2343 let crossings_timer = crate::blend::BlendTimer::start("timing.track.crossings_ms");
2344 crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
2345 crossings_timer.stop();
2346 crate::blend::record_count("count.track.breaks", track.breaks.len() as u64);
2347 for &fraction in refine_near {
2348 for _ in 0..2 {
2349 let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
2350 if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], ceiling)?.is_none() {
2351 break;
2352 }
2353 }
2354 }
2355 let fit_timer = crate::blend::BlendTimer::start("timing.track.fit_local_ms");
2362 let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
2363 surface,
2364 &mut track,
2365 &at,
2366 tolerance,
2367 REFINE_ROUNDS,
2368 ceiling,
2369 )?;
2370 fit_timer.stop();
2371 crate::blend::record_count("count.track.fit_sections", fit.sections as u64);
2372 crate::blend::record_count("count.track.local_inserted", inserted as u64);
2373 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2374 eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2375 eprintln!(
2376 "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2377 fit.miss,
2378 fit.samples,
2379 track.fractions.len(),
2380 track.breaks.len(),
2381 fit.on_floor
2382 );
2383 if let Some(worst) = fit.worst {
2388 let parameter = d0 + (d1 - d0) * worst.fraction;
2389 let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2390 .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2391 .map(|pair| [pair[0], pair[1]]).collect();
2392 let mean = (d1 - d0) / spans.len().max(1) as f64;
2393 let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2394 let width = span.map(|index| spans[index][1] - spans[index][0]);
2395 eprintln!(
2396 "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2397 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2398 worst.fraction, parameter, span, spans.len(), width, mean,
2399 worst.bracket[0], worst.bracket[1], worst.interval
2400 );
2401 }
2402 }
2403 if !fit.on_floor {
2410 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2411 for (label, t) in [("start", d0), ("end", d1)] {
2414 if let Ok(point) = piece.evaluate(t) {
2415 let standoff = crate::project_point_to_surface(surface, point).map(|p| p.distance).unwrap_or(f64::NAN);
2416 eprintln!(
2417 "TRACK_FIT closed-chain piece {label} ({:.9}, {:.9}, {:.9}) stands {standoff:.3e} off the carrier",
2418 point.x, point.y, point.z
2419 );
2420 }
2421 }
2422 }
2423 return Err(KernelRefusal::non_convergence(KernelStage::Refine, CHAIN_PCURVE_FLOOR_WHAT, format!(
2424 "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2425 against a floor of {tolerance:.1e}",
2426 crate::blend::PCURVE_OFF_FLOOR,
2427 fit.miss,
2428 fit.samples
2429 )));
2430 }
2431 Ok(fit.curve)
2432}
2433
2434
2435
2436pub(in crate::blend) struct PreparedWall<'a> {
2439 fine: &'a [Vec3],
2440 coarse: Vec<Vec3>,
2441 fine_stencils: StencilPreparation,
2442 coarse_stencils: StencilPreparation,
2443}
2444
2445#[derive(Clone, Copy)]
2446struct WallStart {
2447 curve: (Vec3, Vec3, Vec3),
2448 point_length: f64,
2449 tangent_length: f64,
2450 tangent_squared: f64,
2451}
2452impl WallStart {
2453 fn new(curve: (Vec3, Vec3, Vec3)) -> Self {
2454 Self { point_length: curve.0.length(), tangent_length: curve.1.length(),
2455 tangent_squared: curve.1.dot(curve.1), curve }
2456 }
2457}
2458
2459struct StencilPreparation {
2460 finite: bool,
2461 coefficients: std::cell::RefCell<Vec<Option<[Vec3; 4]>>>,
2462 starts: std::cell::RefCell<Vec<[Option<WallStart>; 2]>>,
2465}
2466impl StencilPreparation {
2467 fn new(nodes: &[Vec3]) -> Self {
2468 Self { finite: wall_nodes_finite(nodes), coefficients: std::cell::RefCell::new(vec![None; nodes.len().saturating_sub(3)]),
2469 starts: std::cell::RefCell::new(vec![[None; 2]; nodes.len().saturating_sub(3)]) }
2470 }
2471 fn start(&self, base: usize, j: usize, curve: &impl Fn(f64) -> (Vec3, Vec3, Vec3)) -> Option<WallStart> {
2472 let offset = j.checked_sub(base)?;
2473 if !(1..=2).contains(&offset) { return None; }
2474 let mut entries = self.starts.borrow_mut();
2475 Some(*entries.get_mut(base)?[offset-1].get_or_insert_with(|| WallStart::new(curve(offset as f64))))
2476 }
2477 fn coefficients(&self, nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2478 let mut entries = self.coefficients.borrow_mut();
2479 *entries[base].get_or_insert_with(|| wall_coefficients(nodes, base))
2480 }
2481}
2482
2483fn wall_nodes_finite(nodes: &[Vec3]) -> bool {
2484 nodes.iter().all(|c| c.x.is_finite() && c.y.is_finite() && c.z.is_finite())
2485}
2486
2487fn wall_coefficients(nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2488 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
2490 let d1 = p1.sub(p0);
2491 let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
2492 let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
2493 [p0, d1, d2, d3]
2494}
2495
2496fn wall_coarse(fine: &[Vec3]) -> Vec<Vec3> {
2497 fine.iter().step_by(2).copied().collect()
2498}
2499
2500fn wall_preparation_on() -> bool {
2501 switch!(live "BREP_WALL_PREPARE", true)
2502}
2503
2504impl<'a> PreparedWall<'a> {
2505 pub(in crate::blend) fn for_interval(fine: &'a [Vec3]) -> Option<Self> {
2506 (wall_preparation_on() && wall_solve_reuse_on()).then(|| Self::new(fine))
2507 }
2508 fn new(fine: &'a [Vec3]) -> Self {
2509 let coarse = wall_coarse(fine);
2510 let fine_stencils = StencilPreparation::new(fine);
2511 let coarse_stencils = StencilPreparation::new(&coarse);
2512 Self { fine, coarse, fine_stencils, coarse_stencils }
2513 }
2514 pub(in crate::blend) fn read(&self, p: Vec3, j: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
2515 let read = |at| declared_wall_probe_prepared(p, self.fine, at, radius, true, Some(self));
2516 match read(j) {
2517 Some(read) => Some(read),
2518 None => match centre_path_node(self.fine, fraction) {
2519 Some(located) if located != j => read(located),
2520 _ => None,
2521 },
2522 }
2523 }
2524}