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(
90 solid: &BrepSolid,
91 segments: &[ChainSegment<'_>],
92 radius: f64,
93 chamfer: bool,
94 name: Option<&str>,
95) -> Result<BrepSolid, KernelRefusal> {
96 let mode = std::env::var("BREP_CHAIN_LADDER_MEMO");
97 let (memo_on, verify) = match mode.as_deref() {
98 Ok("0") => (false, false),
99 Ok("verify") => (true, true),
100 _ => (true, false),
101 };
102 if !memo_on {
103 return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
104 }
105 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
106 let signed = if chain_ladder_signed_key_on() {
114 let mut signs = Vec::with_capacity(segments.len());
115 for segment in segments {
116 match crate::blend::stations::signed_radii(
117 segment.edge,
118 segment.first.face,
119 segment.first.coedge,
120 segment.second.face,
121 segment.second.coedge,
122 radius,
123 ) {
124 Ok(pair) => signs.push(pair),
125 Err(_) => return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name),
127 }
128 }
129 Some(signs)
130 } else {
131 None
132 };
133 let key = ladder_memo_key(segments, radius, chamfer, bar, signed.as_deref());
134 let cached = LADDER_MEMO.with(|memo| memo.borrow().iter().find(|(stored, _)| *stored == key).map(|(_, refusal)| refusal.clone()));
135 if let (Some(refusal), false) = (&cached, verify) {
136 crate::blend::record_count("count.chain_ladder_memo.hit", 1);
137 return Err(refusal.clone());
138 }
139 let result = blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
140 if let Some(refusal) = cached {
141 let same = matches!(&result, Err(again) if *again == refusal && format!("{again:?}") == format!("{refusal:?}"));
142 assert!(same, "chain ladder memo mismatch: cached {refusal:?}, recomputed {:?}", result.as_ref().err());
143 crate::blend::record_count("count.chain_ladder_memo.verified", 1);
144 return result;
145 }
146 if let Err(refusal) = &result {
147 if matches!(&refusal.class, crate::RefusalClass::NonConvergence { what } if what == "chain_ladder_top") {
148 LADDER_MEMO.with(|memo| {
149 let mut memo = memo.borrow_mut();
150 if memo.len() >= LADDER_MEMO_CAPACITY {
151 memo.remove(0);
152 }
153 memo.push((key, refusal.clone()));
154 });
155 }
156 }
157 result
158}
159
160fn chain_ladder_signed_key_on() -> bool {
164 switch!("BREP_CHAIN_LADDER_SIGNED_KEY", false)
165}
166
167const LADDER_MEMO_CAPACITY: usize = 32;
169
170thread_local! {
171 static LADDER_MEMO: std::cell::RefCell<Vec<(Vec<u64>, KernelRefusal)>> = const { std::cell::RefCell::new(Vec::new()) };
173}
174
175fn ladder_memo_key(segments: &[ChainSegment<'_>], radius: f64, chamfer: bool, bar: f64, signed: Option<&[(f64, f64)]>) -> Vec<u64> {
177 fn curve(key: &mut Vec<u64>, curve: &NurbsCurve) {
178 key.extend([curve.degree as u64, curve.knots.len() as u64, curve.control_points.len() as u64]);
179 key.extend(curve.knots.iter().map(|knot| knot.to_bits()));
180 for point in &curve.control_points {
181 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
182 }
183 }
184 fn face(key: &mut Vec<u64>, face: &FaceRecord, loops: bool) {
185 let surface = &face.surface;
186 key.extend([face.id, face.same_sense as u64, surface.degree_u as u64, surface.degree_v as u64]);
187 key.extend([surface.knots_u.len() as u64, surface.knots_v.len() as u64, surface.control_points.len() as u64]);
188 key.extend(surface.knots_u.iter().chain(&surface.knots_v).map(|knot| knot.to_bits()));
189 for row in &surface.control_points {
190 key.push(row.len() as u64);
191 for point in row {
192 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
193 }
194 }
195 if !loops {
196 return;
197 }
198 key.push(face.loops.len() as u64);
199 for loop_record in &face.loops {
200 key.extend([loop_record.id, loop_record.coedges.len() as u64]);
201 for coedge in &loop_record.coedges {
202 key.extend([coedge.id, coedge.edge_id, coedge.forward as u64]);
203 curve(key, &coedge.pcurve);
204 }
205 }
206 }
207 let mut key = vec![
208 radius.to_bits(),
209 chamfer as u64,
210 bar.to_bits(),
211 crate::blend::station_refinement_on() as u64,
212 crate::blend::chain_arc_stations_on() as u64,
213 wall_verdict_on() as u64,
214 local_continue_on() as u64,
215 segments.len() as u64,
216 signed.is_some() as u64,
217 ];
218 for (index, segment) in segments.iter().enumerate() {
219 if let Some((rho1, rho2)) = signed.and_then(|signs| signs.get(index)) {
220 key.extend([rho1.to_bits(), rho2.to_bits()]);
221 }
222 let edge = segment.edge;
223 key.extend([edge.id, edge.start_vertex_id, edge.end_vertex_id, edge.t0.to_bits(), edge.t1.to_bits(), segment.forward as u64]);
224 curve(&mut key, &edge.curve);
225 for mate in [&segment.first, &segment.second] {
226 key.extend([mate.coedge.id, mate.coedge.edge_id, mate.coedge.forward as u64, mate.loop_index as u64, mate.rho.to_bits()]);
227 curve(&mut key, &mate.coedge.pcurve);
228 face(&mut key, mate.face, signed.is_none());
229 }
230 }
231 key
232}
233fn blend_closed_smooth_chain_unmemoized(
246 solid: &BrepSolid,
247 segments: &[ChainSegment<'_>],
248 radius: f64,
249 chamfer: bool,
250 name: Option<&str>,
251) -> Result<BrepSolid, KernelRefusal> {
252 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
253 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
256 let mut trend: Vec<(usize, f64)> = Vec::new();
258 let first = loop {
259 let per_segment = budget.rung;
260 let rung = build_closed_smooth_chain(
261 solid,
262 segments,
263 radius,
264 chamfer,
265 name,
266 &budget,
267 FoldPolicy::Refuse,
268 bar,
269 );
270 match rung {
271 Ok(Rung::TooCoarse { deviation }) => {
272 crate::blend::carve::carve_trace(format_args!(
273 "blend chain: {} leaves the rails \
274 {deviation:.6e} from the rolling ball, against {bar:.3e}",
275 budget.describe()
276 ));
277 trend.push((per_segment, deviation));
278 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
279 let (best_rung, best, listing) = ladder_trend(&trend);
280 return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
281 "blend: the chain's ladder reached its top, {per_segment} stations per segment, \
282 and no rung passed — its best, {best_rung} stations per segment, left this \
283 collar's fitted rails {best:.6e} from the rolling ball's own contacts against \
284 a bar of {bar:.3e} (rung by rung: {listing}), so there is no wall accurate \
285 enough to build"
286 )));
287 }
288 budget = budget.doubled();
289 }
290 Ok(Rung::StationFailed(failure)) => {
295 crate::blend::carve::carve_trace(format_args!(
296 "blend chain: {per_segment} stations per segment: segment {} position {} \
297 (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
298 failure.segment, failure.position, failure.t, failure.from_position,
299 failure.from_t, failure.seed_uv, failure.error
300 ));
301 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
302 break Err(failure.error);
303 }
304 budget = budget.doubled();
305 }
306 other => break other,
307 }
308 };
309 if let Err(error) = &first {
310 crate::blend::carve::carve_trace(format_args!(
311 "blend chain: the march at {} refused: {error}", budget.describe()
312 ));
313 }
314 match first {
315 Ok(Rung::Built(built)) => Ok(built),
316 Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
317 Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
318 }
319 Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
325 Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
328 carve_ladder(solid, segments, radius, name, bar)
329 .map_err(|carve_error| {
337 if crate::blend::is_wall_fold(&carve_error) {
338 carve_error
339 } else {
340 crate::blend::fold::wall_fold_refusal(
341 Vec::new(),
342 format!(
343 "{} {} fits this edge: the wall folds, and the fold band could not \
344 be carved — {carve_error}",
345 crate::blend::WALL_FOLDS,
346 radius.abs()
347 ),
348 )
349 }
350 })
351 }
352 Err(error) => Err(error),
353 }
354}
355
356fn carve_ladder(
366 solid: &BrepSolid,
367 segments: &[ChainSegment<'_>],
368 radius: f64,
369 name: Option<&str>,
370 bar: f64,
371) -> Result<BrepSolid, KernelRefusal> {
372 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
373 loop {
374 let per_segment = budget.rung;
375 match build_closed_smooth_chain(
376 solid,
377 segments,
378 radius,
379 false,
380 name,
381 &budget,
382 FoldPolicy::Carve,
383 bar,
384 )? {
385 Rung::Built(built) => {
386 crate::blend::carve::carve_trace(format_args!(
387 "blend carve: {} is the first rung inside \
388 {bar:.3e}",
389 budget.describe()
390 ));
391 return Ok(built);
392 }
393 Rung::StationFailed(failure) => return Err(failure.error),
396 Rung::TooCoarse { deviation } => {
397 crate::blend::carve::carve_trace(format_args!(
398 "blend carve: {} leaves the rails \
399 {deviation:.6e} from the rolling ball, against {bar:.3e}",
400 budget.describe()
401 ));
402 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
403 return Err(crate::blend::fold::wall_fold_refusal(
404 Vec::new(),
405 format!(
406 "{} {} fits this edge: the wall folds, and at {per_segment} stations \
407 per segment — the top of the carve's ladder — its fitted rails are \
408 still {deviation:.6e} from the rolling ball's own contacts against a \
409 bar of {bar:.3e}, so there is no wall accurate enough to carve",
410 crate::blend::WALL_FOLDS,
411 radius.abs()
412 ),
413 ));
414 }
415 budget = budget.doubled();
416 }
417 }
418 }
419}
420
421enum Rung {
424 Built(BrepSolid),
425 TooCoarse { deviation: f64 },
426 StationFailed(ChainStationFailure),
429}
430
431fn rail_deviation(
455 samples: &[ChainSample],
456 cr: &NurbsCurve,
457 cs: &NurbsCurve,
458 fit_low: f64,
459 fit_range: f64,
460 budget: &ChainBudget,
461 bar: f64,
462) -> Result<(f64, Vec<usize>), KernelRefusal> {
463 let mut order: Vec<usize> = (0..samples.len())
467 .filter(|&index| (0..=budget.count(samples[index].segment) as isize).contains(&samples[index].position))
468 .collect();
469 order.sort_by(|&a, &b| {
470 let (a, b) = (&samples[a], &samples[b]);
471 a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
472 });
473 let next_of: std::collections::HashMap<usize, f64> = order
474 .windows(2)
475 .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
476 .map(|pair| (pair[0], samples[pair[1]].parameter))
477 .collect();
478 let mut failing = Vec::new();
479 let closest = closest_on_rail;
480 let mut worst: f64 = 0.0;
481 let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
484 for &index in &order {
485 let sample = &samples[index];
486 let Some([p1, p2]) = sample.midpoint else {
487 continue;
488 };
489 let Some(&next) = next_of.get(&index) else {
490 continue;
491 };
492 let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
493 let seed = (global - fit_low) / fit_range;
494 let reach = (next - sample.parameter).abs() / fit_range;
495 let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
496 worst = worst.max(miss1).max(miss2);
497 if miss1.max(miss2) > bar {
498 failing.push(index);
499 }
500 if per_segment_worst.len() <= sample.segment {
501 per_segment_worst.resize(
502 sample.segment + 1,
503 (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
504 );
505 }
506 if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
507 per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
508 }
509 }
510 for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
511 crate::blend::carve::carve_trace(format_args!(
512 " rail miss segment {segment}: worst {miss:.6e} after station {position} \
513 (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
514 p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
515 ));
516 }
517 Ok((worst, failing))
518}
519
520const PIECE_TRACK_CEILING: usize = 4096;
524
525fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
528 let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
529 let (u, _) = crate::projection::curve_foot_newton(curve, target, seed, [low, high],
530 crate::projection::CurveFootPolicy {
531 domain: crate::projection::FootDomain::Clamp,
532 iterations: 40,
533 slope_guard: -1.0, step_limit: reach,
535 parameter_tolerance: 1e-15,
536 mode: crate::projection::FootMode::Backtracked,
537 }).or_refuse(KernelStage::Refine, "curve_foot_newton")?;
538 Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
539}
540
541pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
544
545fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
549 (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
550}
551
552#[allow(clippy::too_many_arguments)]
558pub(super) fn dense_rail_misses(
559 contacts: &crate::blend::chain::DenseContacts,
560 samples: &[ChainSample],
561 cr: &NurbsCurve,
562 cs: &NurbsCurve,
563 low: f64,
564 range: f64,
565 wrap: bool,
566) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
567 let mut readable = Vec::with_capacity(contacts.solved.len());
568 let mut unreadable = Vec::new();
569 for (left, next, fraction, p1, p2) in &contacts.solved {
570 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
571 let global = from + fraction * (to - from);
575 let global = if wrap { global.rem_euclid(1.0) } else { global };
576 let seed = (global - low) / range;
577 let reach = (to - from).abs() / range;
578 let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
581 let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
582 match rail_pair_miss(cr_miss, cs_miss) {
583 Some(miss) => readable.push((*left, miss)),
584 None => unreadable.push(*left),
585 }
586 }
587 Ok((readable, unreadable))
588}
589
590pub(super) fn combine_rail_verdict(
596 halfway: (f64, Vec<usize>),
597 dense: &[(usize, f64)],
598 unsolved: &[usize],
599 bar: f64,
600) -> (f64, Vec<usize>) {
601 let (mut worst, mut failing) = halfway;
602 for &(left, miss) in dense {
603 if !miss.is_finite() {
604 failing.push(left);
605 continue;
606 }
607 if miss > bar {
608 failing.push(left);
609 }
610 worst = worst.max(miss);
611 }
612 failing.extend_from_slice(unsolved);
613 failing.sort_unstable();
614 failing.dedup();
615 (worst, failing)
616}
617
618pub(super) fn model_rounds_limit() -> usize {
625 REFINE_ROUNDS
626}
627
628pub(super) fn model_request_on() -> bool {
629 wall_verdict_on()
630}
631
632
633pub(super) fn wall_verdict_on() -> bool {
636 true
637}
638
639
640
641
642
643pub(super) const WALL_SECTION_PROBES: usize = 8;
646
647pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
668 let j = locate_foot(p, nodes, j, radius)?;
669 sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
670}
671
672const FOOT_TRAVEL: usize = 4;
675
676fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
686 locate_foot_read(p, nodes, j, radius).map(|(node, _)| node)
687}
688
689fn locate_foot_read(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(usize, f64)> {
692 locate_foot_read_prepared(p, nodes, j, radius, None)
693}
694
695fn locate_foot_read_prepared(p: Vec3, nodes: &[Vec3], j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Option<(usize, f64)> {
696 let n = nodes.len();
697 if n < 4 || j == 0 || j + 1 >= n {
698 return None;
699 }
700 let mut tried: Vec<usize> = Vec::new();
701 let mut at = j;
702 loop {
703 let base = (at - 1).min(n - 4);
704 match stencil_solve_prepared(p, nodes, base, at, radius, prepared) {
705 Solve::Read(distance) => return Some((at, distance)),
706 Solve::Unread => return None,
707 Solve::Window(s) => {
708 tried.push(at);
709 let next = (base as f64 + s).round();
710 let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
711 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
712 wall_decline(|| format!(
713 "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
714 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
715 ));
716 return None;
717 }
718 at = next;
719 }
720 }
721 }
722}
723
724fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
727 match stencil_solve(p, nodes, base, j, radius) {
728 Solve::Read(distance) => Some(distance),
729 Solve::Window(_) | Solve::Unread => None,
730 }
731}
732
733enum Solve {
736 Read(f64),
737 Window(f64),
738 Unread,
739}
740
741thread_local! {
742 static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
747}
748
749fn wall_decline(reason: impl FnOnce() -> String) {
750 WALL_DECLINES.with(|log| {
751 if let Some(lines) = log.borrow_mut().as_mut() {
752 lines.push(reason());
753 }
754 });
755}
756
757pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
765 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
766 let read = declared_wall_probe(p, fine, jf, radius);
767 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
768 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
769}
770
771fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
772 stencil_solve_prepared(p, nodes, base, j, radius, None)
773}
774
775fn stencil_solve_prepared(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Solve {
776 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
777 let n = nodes.len();
778 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) {
779 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
780 return Solve::Unread;
781 }
782 let [p0, d1, d2, d3] = prepared.map_or_else(|| wall_coefficients(nodes, base), |cache| cache.coefficients(nodes, base));
783 let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
784 let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
785 let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
786 let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
787 (c, c1, c2)
788 };
789 let start = j as f64 - base as f64;
790 let initial = prepared.and_then(|cache| cache.start(base, j, &curve));
791 let point_length = p.length();
792 let mut s = start;
793 let mut converged = false;
794 let mut converged_point = None;
795 for iteration in 0..50 {
796 let state = if iteration == 0 { initial.unwrap_or_else(|| WallStart::new(curve(s))) }
797 else { WallStart::new(curve(s)) };
798 let (c, c1, c2) = state.curve;
799 let w = c.sub(p);
800 let (g, slope) = (w.dot(c1), state.tangent_squared + w.dot(c2));
801 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
802 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
803 return Solve::Unread;
804 }
805 let stationary = g.abs() <= 16.0 * f64::EPSILON * (point_length + state.point_length) * state.tangent_length;
818 if stationary {
819 converged = true;
820 converged_point = Some(c);
821 break;
822 }
823 let step = g / slope;
824 let before = s;
825 s -= step;
826 if !((s - start).abs() <= 1.0) {
827 wall_decline(|| {
828 let (c, c1, _) = curve(s);
831 let g_after = c.sub(p).dot(c1);
832 let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
833 format!(
834 "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
835 g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
836 (s - start).abs() - 1.0
837 )
838 });
839 return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
840 }
841 if step.abs() <= 1e-14 * (1.0 + s.abs()) {
842 converged = true;
843 break;
844 }
845 }
846 if !converged {
847 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
848 return Solve::Unread;
849 }
850 let distance = converged_point.unwrap_or_else(|| curve(s).0).sub(p).length() - radius;
851 if distance.is_finite() {
852 Solve::Read(distance)
853 } else {
854 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
855 Solve::Unread
856 }
857}
858
859pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
868 let n = nodes.len();
869 if n < 5 || j == 0 || j + 1 >= n {
870 return None;
871 }
872 let j = locate_foot(p, nodes, j, radius)?;
873 let first = (j - 1).min(n - 4);
874 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
875 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
876 let a = sweep_residual_on(p, nodes, first, j, radius)?;
877 let b = sweep_residual_on(p, nodes, second, j, radius)?;
878 Some((a.abs().max(b.abs()), (a - b).abs()))
879}
880
881fn wall_solve_reuse_on() -> bool {
894 switch!(live "BREP_WALL_SOLVE_REUSE", true)
895}
896
897pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
898 declared_wall_probe_with(p, fine, jf, radius, wall_solve_reuse_on())
899}
900
901fn declared_wall_probe_with(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool) -> Option<(f64, f64)> {
902 declared_wall_probe_prepared(p, fine, jf, radius, reuse, None)
903}
904
905fn declared_wall_probe_prepared(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool, prepared: Option<&PreparedWall<'_>>) -> Option<(f64, f64)> {
906 if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
907 wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
908 return None;
909 }
910 let (reading, spread, at_fine) = if reuse {
911 let read = (|| {
912 let n = fine.len();
913 let (j, a) = locate_foot_read_prepared(p, fine, jf, radius, prepared.map(|c| &c.fine_stencils))?;
914 let first = (j - 1).min(n - 4);
915 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
916 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
917 let b = match stencil_solve_prepared(p, fine, second, j, radius, prepared.map(|c| &c.fine_stencils)) {
918 Solve::Read(distance) => distance,
919 _ => return None,
920 };
921 Some((a.abs().max(b.abs()), (a - b).abs(), a))
922 })();
923 let Some(read) = read else {
924 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
925 return None;
926 };
927 read
928 } else {
929 let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
930 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
931 return None;
932 };
933 let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
934 wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
935 return None;
936 };
937 (reading, spread, at_fine)
938 };
939 let coarse_owned;
940 let coarse = match prepared {
941 Some(cache) => cache.coarse.as_slice(),
942 None => { coarse_owned = wall_coarse(fine); &coarse_owned }
943 };
944 let at_coarse = if reuse {
945 locate_foot_read_prepared(p, coarse, jf / 2, radius, prepared.map(|c| &c.coarse_stencils)).map(|(_, distance)| distance)
946 } else {
947 sweep_residual(p, coarse, jf / 2, radius)
948 };
949 let Some(at_coarse) = at_coarse else {
950 wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
951 return None;
952 };
953 let uncertainty = spread.max((at_fine - at_coarse).abs());
954 uncertainty.is_finite().then_some((reading, uncertainty))
955}
956
957pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
973 let n = fine.len();
974 if n < 7 || !(fraction.is_finite()) {
975 return None;
976 }
977 let mut cumulative = Vec::with_capacity(n);
978 cumulative.push(0.0_f64);
979 for pair in fine.windows(2) {
980 let step = pair[1].sub(pair[0]).length();
981 cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
982 }
983 let total = *cumulative.last()?;
984 if !(total.is_finite() && total > 0.0) {
985 return None;
986 }
987 let target = fraction.clamp(0.0, 1.0) * total;
988 let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
989 let even = 2 * ((nearest as f64 / 2.0).round() as usize);
990 Some(even.clamp(2, (n - 3) & !1))
991}
992
993pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
1000 let first = match declared_wall_probe_why(p, fine, jf, radius) {
1001 Ok(read) => return Ok(read),
1002 Err(why) => why,
1003 };
1004 match centre_path_node(fine, fraction) {
1005 Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
1006 .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
1007 _ => Err(first),
1008 }
1009}
1010
1011pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
1013 match declared_wall_probe(p, fine, jf, radius) {
1014 Some(read) => Some(read),
1015 None => match centre_path_node(fine, fraction) {
1016 Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
1017 _ => None,
1018 },
1019 }
1020}
1021
1022pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
1034 if jf % 2 != 0
1035 || fine.len() < 9
1036 || fine.len() % 2 == 0
1037 || ts.len() != fine.len()
1038 || !ts.iter().all(|t| t.is_finite())
1039 || !ts.windows(2).all(|pair| pair[1] > pair[0])
1040 {
1041 return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
1042 }
1043 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
1044 let read = (|| -> Option<(f64, f64)> {
1045 let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
1046 wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
1047 None
1048 })?;
1049 let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
1050 wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
1051 None
1052 })?;
1053 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
1054 let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
1055 let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
1056 wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
1057 None
1058 })?;
1059 let uncertainty = spread.max((at_fine - at_coarse).abs());
1060 uncertainty.is_finite().then_some((reading, uncertainty))
1061 })();
1062 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
1063 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
1064}
1065
1066fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
1067 let n = nodes.len();
1068 if n < 5 || j == 0 || j + 1 >= n {
1069 return None;
1070 }
1071 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1072 let first = (j - 1).min(n - 4);
1073 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
1074 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
1075 let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
1076 let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
1077 Some((a.abs().max(b.abs()), (a - b).abs()))
1078}
1079
1080fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
1081 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1082 match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
1083 Solve::Read(distance) => Some(distance),
1084 Solve::Window(_) | Solve::Unread => None,
1085 }
1086}
1087
1088fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
1092 let n = nodes.len();
1093 if n < 4 || j == 0 || j + 1 >= n {
1094 return None;
1095 }
1096 let mut tried: Vec<usize> = Vec::new();
1097 let mut at = j;
1098 loop {
1099 let base = (at - 1).min(n - 4);
1100 match stencil_solve_at(p, nodes, ts, base, at, radius) {
1101 Solve::Read(_) => return Some(at),
1102 Solve::Unread => return None,
1103 Solve::Window(t) => {
1104 tried.push(at);
1105 let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
1106 let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
1107 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
1108 wall_decline(|| format!(
1109 "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
1110 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
1111 ));
1112 return None;
1113 }
1114 at = next;
1115 }
1116 }
1117 }
1118}
1119
1120fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
1125 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
1126 let n = nodes.len();
1127 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)) {
1128 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
1129 return Solve::Unread;
1130 }
1131 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
1132 let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
1133 let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
1134 let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
1135 let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
1136 let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
1137 let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
1138 let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
1139 let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
1140 let (a, b, c) = (t - t0, t - t1, t - t2);
1141 let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
1142 let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
1143 let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
1144 (point, first, second)
1145 };
1146 let (low, high) = (ts[j - 1], ts[j + 1]);
1147 let start = ts[j];
1148 let mut t = start;
1149 let mut converged = false;
1150 for _ in 0..50 {
1151 let (c, c1, c2) = curve(t);
1152 let w = c.sub(p);
1153 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
1154 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
1155 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
1156 return Solve::Unread;
1157 }
1158 if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
1159 converged = true;
1160 break;
1161 }
1162 let step = g / slope;
1163 let before = t;
1164 t -= step;
1165 if !(t >= low && t <= high) {
1166 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
1167 return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
1168 }
1169 let node = 0.5 * (high - low);
1174 if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
1175 converged = true;
1176 break;
1177 }
1178 }
1179 if !converged {
1180 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
1181 return Solve::Unread;
1182 }
1183 let distance = curve(t).0.sub(p).length() - radius;
1184 if distance.is_finite() {
1185 Solve::Read(distance)
1186 } else {
1187 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
1188 Solve::Unread
1189 }
1190}
1191
1192pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
1197 readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
1198}
1199
1200#[allow(clippy::too_many_arguments)]
1210pub(super) fn wall_interior_misses(
1211 contacts: &crate::blend::chain::DenseContacts,
1212 samples: &[ChainSample],
1213 surface: &crate::NurbsSurface,
1214 low: f64,
1215 range: f64,
1216 radius: f64,
1217 subdivisions: usize,
1218 probes: usize,
1219 wrap: bool,
1220) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1221 use std::collections::BTreeMap;
1222 let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1225 for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1226 let entry = nodes.entry(*left).or_insert_with(|| {
1227 let mut row = vec![None; subdivisions + 1];
1228 row[0] = Some(samples[*left].station.center);
1229 row[subdivisions] = Some(samples[*next].station.center);
1230 (*next, row)
1231 });
1232 let step = (fraction * subdivisions as f64).round() as usize;
1233 if (1..subdivisions).contains(&step) {
1234 entry.1[step] = contacts.centers.get(index).copied().flatten();
1235 }
1236 }
1237 let mut readable = Vec::new();
1238 let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1239 for (&left, (next, row)) in &nodes {
1240 let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1241 unreadable.push(left);
1242 continue;
1243 };
1244 let (from, to) = (samples[left].parameter, samples[*next].parameter);
1245 let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1246 let mut read = true;
1247 let prepared = PreparedWall::for_interval(¢res);
1248 'probes: for step in (2..subdivisions).step_by(2) {
1249 let global = from + (step as f64 / subdivisions as f64) * (to - from);
1250 let global = if wrap { global.rem_euclid(1.0) } else { global };
1251 let u = (global - low) / range;
1252 for k in 1..probes {
1253 let v = k as f64 / probes as f64;
1254 let reading = surface.evaluate(u, v).ok().and_then(|p| {
1255 match &prepared {
1256 Some(cache) => cache.read(p, step, step as f64 / subdivisions as f64, radius),
1257 None => declared_wall_probe_located(p, ¢res, step, step as f64 / subdivisions as f64, radius),
1258 }
1259 });
1260 match reading {
1261 Some((distance, gap)) => {
1262 declared = declared.max(distance + gap);
1263 worst = worst.max(distance);
1264 spread = spread.max(gap);
1265 }
1266 None => {
1267 read = false;
1268 break 'probes;
1269 }
1270 }
1271 }
1272 }
1273 if read {
1274 readable.push((left, declared, worst, spread));
1275 } else {
1276 unreadable.push(left);
1277 }
1278 }
1279 unreadable.sort_unstable();
1280 unreadable.dedup();
1281 (readable, unreadable)
1282}
1283
1284
1285fn build_closed_smooth_chain(
1286 solid: &BrepSolid,
1287 segments: &[ChainSegment<'_>],
1288 radius: f64,
1289 chamfer: bool,
1290 name: Option<&str>,
1291 budget: &ChainBudget,
1292 fold: FoldPolicy,
1293 bar: f64,
1294) -> Result<Rung, KernelRefusal> {
1295 let samples = match march_chain_typed(segments, radius, false, budget, fold) {
1296 Ok(samples) => samples,
1297 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1298 Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1301 Err(ChainMarchError::Other(error)) => return Err(error),
1302 };
1303 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1304 build_closed_chain_from_samples(solid, segments, radius, chamfer, name, budget, fold, bar, samples, 0)
1305}
1306
1307
1308thread_local! {
1309 static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1313}
1314
1315fn local_continue_on() -> bool {
1321 switch!(live "BREP_CHAIN_LOCAL_CONTINUE", false)
1322}
1323
1324const LOCAL_CONTINUE_MAX_ROUNDS: usize = 4 * REFINE_ROUNDS;
1326
1327const LOCAL_CONTINUE_STALL_ROUNDS: usize = 3;
1329
1330thread_local! {
1331 static CHAIN_LOCAL_BEST: std::cell::Cell<(f64, usize)> = const { std::cell::Cell::new((f64::INFINITY, 0)) };
1333}
1334
1335fn local_rounds_spent(round: usize, deviation: f64) -> bool {
1341 if round == 0 {
1342 CHAIN_LOCAL_BEST.with(|best| best.set((f64::INFINITY, 0)));
1343 return false;
1344 }
1345 let (_, best_round) = CHAIN_LOCAL_BEST.with(|best| {
1346 let (reading, at) = best.get();
1347 if deviation < reading {
1348 best.set((deviation, round));
1349 (deviation, round)
1350 } else {
1351 (reading, at)
1352 }
1353 });
1354 round >= LOCAL_CONTINUE_MAX_ROUNDS || round - best_round >= LOCAL_CONTINUE_STALL_ROUNDS
1355}
1356#[allow(clippy::too_many_arguments)]
1367fn build_closed_chain_from_samples(
1368 solid: &BrepSolid,
1369 segments: &[ChainSegment<'_>],
1370 radius: f64,
1371 chamfer: bool,
1372 name: Option<&str>,
1373 budget: &ChainBudget,
1374 fold: FoldPolicy,
1375 bar: f64,
1376 mut samples: Vec<ChainSample>,
1377 round: usize,
1378) -> Result<Rung, KernelRefusal> {
1379 let per_segment = budget.rung;
1382 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1383
1384 let mut global: Vec<(f64, &ChainSample)> = samples
1386 .iter()
1387 .filter(|sample| (0..budget.count(sample.segment) as isize).contains(&sample.position))
1388 .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1389 .collect();
1390 global.sort_by(|a, b| a.0.total_cmp(&b.0));
1391 let count = global.len();
1392 let degree = FIT_DEGREE;
1393 let mut extended_params = Vec::new();
1394 let mut samples_cr = Vec::new();
1395 let mut samples_mid = Vec::new();
1396 let mut samples_cs = Vec::new();
1397 let mut push = |station: &Station, parameter: f64| {
1398 extended_params.push(parameter);
1399 samples_cr.push(Vec4::from_point(station.p1, 1.0));
1400 samples_cs.push(Vec4::from_point(station.p2, 1.0));
1401 samples_mid.push(Vec4 {
1402 x: station.apex.x * station.weight,
1403 y: station.apex.y * station.weight,
1404 z: station.apex.z * station.weight,
1405 w: station.weight,
1406 });
1407 };
1408 for offset in (1..=degree).rev() {
1409 let (parameter, sample) = &global[count - offset];
1410 push(&sample.station, parameter - 1.0);
1411 }
1412 for (parameter, sample) in &global {
1413 push(&sample.station, *parameter);
1414 }
1415 push(&global[0].1.station, global[0].0 + 1.0);
1418 for offset in 1..=degree {
1419 let (parameter, sample) = &global[offset];
1420 push(&sample.station, parameter + 1.0);
1421 }
1422 let low = extended_params[0];
1423 let high = *extended_params.last().unwrap();
1424 let range = high - low;
1425 let normalized: Vec<f64> = extended_params
1426 .iter()
1427 .map(|parameter| (parameter - low) / range)
1428 .collect();
1429 let seam_low = (0.0 - low) / range;
1430 let seam_high = (1.0 - low) / range;
1431 let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1432 let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1433 let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1434 let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1435 Ok(middle)
1436 };
1437 let cr = fit_row(&samples_cr)?;
1438 let cs = fit_row(&samples_cs)?;
1439 let mid = if chamfer {
1440 None
1441 } else {
1442 Some(fit_row(&samples_mid)?)
1443 };
1444 let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1445 let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1446
1447 let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, budget, bar)?;
1460 let stencil = if round > 0 {
1472 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, LOCAL_RAIL_STENCIL)?)
1473 } else {
1474 None
1475 };
1476 let wall_stencil = if !chamfer && wall_verdict_on() {
1481 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1482 } else {
1483 None
1484 };
1485 let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1486 let contacts = stencil.as_ref().expect("the local path solves the stencil");
1487 let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1488 let reads = dense.len() + unreadable.len() + contacts.unsolved;
1489 let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1490 let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1491 crate::blend::carve::carve_trace(format_args!(
1492 "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1493 {} failing interval(s), {} unsolved contact(s)",
1494 verdict.0, verdict.1.len(), contacts.unsolved
1495 ));
1496 (verdict.0, verdict.1, reads)
1497 } else {
1498 (halfway, halfway_failing, 0)
1499 };
1500 let mut _wall_worst = 0.0_f64;
1506 let mut _wall_spread = 0.0_f64;
1507 let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1508 if let Some(contacts) = wall_stencil.as_ref() {
1509 let (wall, unreadable) =
1510 wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1511 wall_readings = wall.clone();
1512 for &(left, declared, _, uncertainty) in &wall {
1514 if declared > bar {
1515 failing.push(left);
1516 }
1517 _wall_worst = _wall_worst.max(declared);
1518 _wall_spread = _wall_spread.max(uncertainty);
1519 }
1520 failing.extend_from_slice(&unreadable);
1521 failing.sort_unstable();
1522 failing.dedup();
1523 deviation = deviation.max(_wall_worst);
1524 crate::blend::carve::carve_trace(format_args!(
1525 "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1526 interval(s), {} unreadable, bar {bar:.3e}",
1527 wall.len(),
1528 unreadable.len()
1529 ));
1530 }
1531 if deviation > bar || !failing.is_empty() {
1532 let in_segment = |segment: usize| {
1533 samples
1534 .iter()
1535 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1536 .count()
1537 };
1538 let within_ceiling = (0..segments.len()).all(|segment| {
1539 in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1540 <= CHAIN_CARVE_MAX_PER_SEGMENT
1541 });
1542 crate::blend::carve::carve_trace(format_args!(
1543 "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1544 worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1545 samples.len(),
1546 failing.len(),
1547 if within_ceiling { "ok" } else { "reached" }
1548 ));
1549 let rounds_spent = if local_continue_on() {
1551 local_rounds_spent(round, deviation)
1552 } else {
1553 round >= REFINE_ROUNDS
1554 };
1555 if fold != FoldPolicy::Refuse || rounds_spent || !within_ceiling {
1557 return Ok(Rung::TooCoarse { deviation });
1558 }
1559 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1561 match insert_chain_stations(segments, radius, budget, &mut samples, &failing) {
1562 Ok(()) => {}
1563 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1564 Err(ChainMarchError::Incoherent(failure)) => {
1568 crate::blend::carve::carve_trace(format_args!(
1569 "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1570 at t {:.9} does not reach its retained right neighbour: {}",
1571 failure.segment, failure.t, failure.error
1572 ));
1573 return Ok(Rung::TooCoarse { deviation });
1574 }
1575 Err(ChainMarchError::Other(error)) => return Err(error),
1576 }
1577 assign_chain_parameters(&mut samples, segments.len(), budget, false);
1578 return build_closed_chain_from_samples(
1579 solid, segments, radius, chamfer, name, budget, fold, bar, samples, round + 1,
1580 );
1581 }
1582 crate::blend::carve::carve_trace(format_args!(
1583 "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1584 from the rolling ball, inside {bar:.3e}"
1585 ));
1586 let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1600 let stations_now = (0..segments.len())
1601 .map(|segment| {
1602 samples
1603 .iter()
1604 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1605 .count()
1606 })
1607 .max()
1608 .unwrap_or(0);
1609 let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1610 let mut request_detail: Option<String> = None;
1614 let mut request_component = crate::MeasuredComponent::Wall;
1615 let mut request_unread = 0usize;
1616 if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1617 let wall_missing = wall_model_missing(&wall_readings, model);
1618 let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1626 Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1627 None => (Vec::new(), Vec::new()),
1628 };
1629 let rail_over: std::collections::BTreeSet<usize> =
1631 rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1632 let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1633 let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1634 let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1635 missing.sort_unstable();
1636 missing.dedup();
1637 if !missing.is_empty() {
1638 let wall_short = wall_readings
1640 .iter()
1641 .filter(|reading| wall_missing.contains(&reading.0))
1642 .map(|&(_, declared, _, _)| declared)
1643 .fold(0.0_f64, f64::max);
1644 let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1645 let short = wall_short.max(rail_short);
1646 request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1647 request_unread = rail_unread_intervals.len();
1648 request_detail = Some(format!(
1649 "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1650 wall {wall_short:.3e} over {} interval(s)",
1651 rail_over.len(),
1652 rail_unread_intervals.len(),
1653 wall_missing.len()
1654 ));
1655 let within_ceiling = (0..segments.len()).all(|segment| {
1656 let present = samples
1657 .iter()
1658 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1659 .count();
1660 present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1661 <= CHAIN_CARVE_MAX_PER_SEGMENT
1662 });
1663 crate::blend::carve::carve_trace(format_args!(
1664 "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1665 (worst reading + spread {short:.3e}), ceiling {}",
1666 missing.len(),
1667 if within_ceiling { "ok" } else { "reached" }
1668 ));
1669 unmet = if round >= model_rounds_limit() {
1670 Some((short, crate::BudgetReason::RoundsSpent))
1671 } else if !within_ceiling {
1672 Some((short, crate::BudgetReason::StationCeiling))
1673 } else {
1674 let mut refined = samples.clone();
1675 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1677 match insert_chain_stations(segments, radius, budget, &mut refined, &missing) {
1678 Ok(()) => {
1679 assign_chain_parameters(&mut refined, segments.len(), budget, false);
1680 match build_closed_chain_from_samples(
1681 solid, segments, radius, chamfer, name, budget, fold, bar, refined, round + 1,
1682 ) {
1683 Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1684 _ => Some((short, crate::BudgetReason::Rejected)),
1685 }
1686 }
1687 Err(_) => {
1688 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1690 Some((short, crate::BudgetReason::Incoherent))
1691 }
1692 }
1693 };
1694 }
1695 }
1696 if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1701 let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, subdivisions.max(2))?;
1702 let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1703 let mut unreadable = 0usize;
1704 for (left, next, fraction, p1, p2) in &contacts.solved {
1705 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1706 let global = (from + fraction * (to - from)).rem_euclid(1.0);
1707 let seed = (global - low) / range;
1708 let reach = (to - from).abs() / range;
1709 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)) {
1710 Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1711 Some(_) => {}
1712 None => unreadable += 1,
1713 }
1714 }
1715 eprintln!(
1716 "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1717 worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1718 contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1719 );
1720 }
1721 let crease = if fold == FoldPolicy::Carve {
1722 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1732 let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1733 return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1734 "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"
1735 )));
1736 };
1737 Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1738 } else {
1739 None
1740 };
1741
1742 let single_face = |side: usize| -> bool {
1746 let first_id = segments[0].mate(side).face.id;
1747 segments
1748 .iter()
1749 .all(|segment| segment.mate(side).face.id == first_id)
1750 };
1751 let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1752 for side in 0..2 {
1753 if !single_face(side) {
1754 continue;
1755 }
1756 let mut row = Vec::new();
1757 let select = |station: &Station| -> [f64; 2] {
1758 if side == 0 {
1759 station.uv1
1760 } else {
1761 station.uv2
1762 }
1763 };
1764 let mut push_uv = |station: &Station| {
1765 let uv = select(station);
1766 row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1767 };
1768 for offset in (1..=degree).rev() {
1769 push_uv(&global[count - offset].1.station);
1770 }
1771 for (_, sample) in &global {
1772 push_uv(&sample.station);
1773 }
1774 push_uv(&global[0].1.station);
1775 for offset in 1..=degree {
1776 push_uv(&global[offset].1.station);
1777 }
1778 whole_pcurves[side] = Some(fit_row(&row)?);
1779 }
1780
1781 let mut pcurves1 = Vec::with_capacity(segments.len());
1784 let mut pcurves2 = Vec::with_capacity(segments.len());
1785 for segment in 0..segments.len() {
1786 let mut window: Vec<&ChainSample> = samples
1787 .iter()
1788 .filter(|sample| sample.segment == segment)
1789 .collect();
1790 window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1791 let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1792 let low = params[0];
1793 let high = *params.last().unwrap();
1794 let normalized: Vec<f64> = params
1795 .iter()
1796 .map(|parameter| (parameter - low) / (high - low))
1797 .collect();
1798 let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1799 let points: Vec<Vec4> = window
1800 .iter()
1801 .map(|sample| {
1802 let uv = select(&sample.station);
1803 Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1804 })
1805 .collect();
1806 let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1807 Ok(curve)
1811 };
1812 pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1813 pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1814 }
1815
1816 let wall_surface = unmet.as_ref().map(|_| surface.clone());
1817 let built = chain_surgery(
1818 solid,
1819 segments,
1820 ChainRows {
1821 surface,
1822 cr,
1823 cs,
1824 u_domain,
1825 fit_low: low,
1826 fit_range: range,
1827 pcurves1,
1828 pcurves2,
1829 whole_pcurves,
1830 crease,
1831 },
1832 name,
1833 )?;
1834 if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1835 crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1836 name: name.map(str::to_string),
1837 surface,
1838 request: model,
1839 residual,
1840 detail: request_detail,
1841 budget: crate::ApproximationBudget {
1842 reason,
1843 rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1846 rounds_limit: model_rounds_limit(),
1847 stations: stations_now,
1848 station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1849 mechanism: crate::BudgetMechanism::LocalRounds,
1850 measured_component: request_component,
1851 unread: request_unread,
1852 },
1853 });
1854 }
1855 Ok(Rung::Built(built))
1856}
1857
1858pub(super) struct ChainRows {
1859 pub(super) surface: NurbsSurface,
1860 pub(super) cr: NurbsCurve,
1861 pub(super) cs: NurbsCurve,
1862 pub(super) u_domain: [f64; 2],
1863 pub(super) fit_low: f64,
1867 pub(super) fit_range: f64,
1868 pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1871 pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1872 pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1874 pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1877}
1878
1879impl ChainRows {
1880 pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1881 self.fit_low + fit_parameter * self.fit_range
1882 }
1883}
1884
1885pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1889 let (low, high, curve) = fitted;
1890 let to_local = |value: f64| {
1893 let mut best = value;
1894 for candidate in [value, value + 1.0, value - 1.0] {
1895 if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1896 best = candidate;
1897 break;
1898 }
1899 }
1900 ((best - low) / (high - low)).clamp(0.0, 1.0)
1901 };
1902 let mut local_a = to_local(a);
1903 let mut local_b = to_local(b);
1904 if local_a > local_b {
1905 std::mem::swap(&mut local_a, &mut local_b);
1906 }
1907 let epsilon = 1e-9;
1908 let (_, tail) = if local_a > epsilon {
1909 curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1910 } else {
1911 (curve.clone(), curve.clone())
1912 };
1913 let tail = if local_a > epsilon {
1914 tail
1915 } else {
1916 curve.clone()
1917 };
1918 let portion = if local_b < 1.0 - epsilon {
1919 tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1920 } else {
1921 tail
1922 };
1923 Ok(portion)
1924}
1925
1926pub(super) fn cross_edge_at(
1932 solid: &BrepSolid,
1933 face_before: &FaceRecord,
1934 loop_before: usize,
1935 face_after: &FaceRecord,
1936 loop_after: usize,
1937 vertex: u64,
1938 before_edge: u64,
1939 after_edge: u64,
1940) -> Result<Option<u64>, KernelRefusal> {
1941 let mut found: Option<u64> = None;
1942 let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1943 for coedge in &face.loops[loop_index].coedges {
1944 if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1945 continue;
1946 }
1947 let edge = solid
1948 .edges
1949 .iter()
1950 .find(|edge| edge.id == coedge.edge_id)
1951 .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1952 if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1953 if found.is_some() && found != Some(edge.id) {
1954 return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1955 }
1956 found = Some(edge.id);
1957 }
1958 }
1959 Ok(())
1960 };
1961 consider(face_before, loop_before)?;
1962 if !(face_after.id == face_before.id && loop_after == loop_before) {
1963 consider(face_after, loop_after)?;
1964 }
1965 Ok(found)
1966}
1967
1968fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1971 let period = high - low;
1972 low + ((value - low) / period).round() * period
1973}
1974
1975pub(in crate::blend) fn project_piece_pcurve(
1993 piece: &NurbsCurve,
1994 surface: &NurbsSurface,
1995 start_is_crossing: bool,
1996 end_is_crossing: bool,
1997) -> Result<NurbsCurve, KernelRefusal> {
1998 project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1999}
2000
2001
2002pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
2007 let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
2008 let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2009 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2010 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2011 let into = move |value: f64, closed: bool, low: f64, high: f64| {
2012 if closed && high > low {
2013 low + (value - low).rem_euclid(high - low)
2014 } else {
2015 value.clamp(low, high)
2016 }
2017 };
2018 Ok(move |fraction: f64| {
2019 let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
2020 let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
2021 seed.iter().all(|value| value.is_finite()).then_some(seed)
2022 })
2023}
2024
2025pub(in crate::blend) fn project_piece_pcurve_from(
2035 piece: &NurbsCurve,
2036 surface: &NurbsSurface,
2037 start_is_crossing: bool,
2038 end_is_crossing: bool,
2039 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2040) -> Result<NurbsCurve, KernelRefusal> {
2041 project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
2042}
2043
2044pub(crate) const CHAIN_PCURVE_FLOOR_WHAT: &str = "chain_pcurve_floor";
2050
2051pub(in crate::blend) fn project_piece_pcurve_within(
2052 piece: &NurbsCurve,
2053 surface: &NurbsSurface,
2054 start_is_crossing: bool,
2055 end_is_crossing: bool,
2056 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2057 refine_near: &[f64],
2058) -> Result<NurbsCurve, KernelRefusal> {
2059 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false)
2060}
2061
2062pub(in crate::blend) fn project_piece_pcurve_within_knot_retry(
2073 piece: &NurbsCurve,
2074 surface: &NurbsSurface,
2075 start_is_crossing: bool,
2076 end_is_crossing: bool,
2077 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2078 refine_near: &[f64],
2079) -> Result<NurbsCurve, KernelRefusal> {
2080 match project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false) {
2081 Err(uniform) if matches!(&uniform.class, crate::RefusalClass::NonConvergence { what } if what == CHAIN_PCURVE_FLOOR_WHAT) => {
2082 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, true).map_err(|_| uniform)
2083 }
2084 other => other,
2085 }
2086}
2087
2088#[allow(clippy::too_many_arguments)]
2092fn project_piece_pcurve_track(
2093 piece: &NurbsCurve,
2094 surface: &NurbsSurface,
2095 start_is_crossing: bool,
2096 end_is_crossing: bool,
2097 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2098 refine_near: &[f64],
2099 knot_aware: bool,
2100) -> Result<NurbsCurve, KernelRefusal> {
2101 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2102 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2103 let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2104 let period = u1 - u0;
2105 let period_v = v1 - v0;
2106 let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
2107 let ceiling = if knot_aware { 2 * PIECE_TRACK_CEILING } else { PIECE_TRACK_CEILING };
2117 let uniform = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
2118 let mut fractions: Vec<f64> = (0..=uniform).map(|section| section as f64 / uniform as f64).collect();
2119 if knot_aware {
2120 let span = d1 - d0;
2123 let mut knots: Vec<f64> = piece.knots.iter().copied().filter(|knot| *knot > d0 && *knot < d1).map(|knot| (knot - d0) / span).collect();
2124 knots.dedup();
2125 fractions.extend(knots);
2126 fractions.sort_by(|a, b| a.total_cmp(b));
2127 fractions.dedup_by(|b, a| (*b - *a).abs() <= 1e-6 / uniform as f64);
2128 }
2129 let dense = fractions.len() - 1;
2130 crate::blend::record_count("count.track.piece_samples", dense as u64 + 1);
2131 let mut proj_u = Vec::with_capacity(dense + 1);
2132 let mut proj_v = Vec::with_capacity(dense + 1);
2133 let mut point3 = Vec::with_capacity(dense + 1);
2134 let mut params = Vec::with_capacity(dense + 1);
2135 let mut previous_u: Option<f64> = None;
2136 let mut previous_foot: Option<[f64; 2]> = None;
2141 let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
2142 let project_timer = crate::blend::BlendTimer::start("timing.track.project_ms");
2143 for section in 0..=dense {
2144 let fraction = fractions[section];
2145 let t = if knot_aware { d0 + (d1 - d0) * fraction } else { d0 + (d1 - d0) * section as f64 / dense as f64 };
2149 let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
2150 let guide = guide.filter(|_| {
2151 true
2152 });
2153 let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(fraction)).into_iter().chain(previous_foot).collect();
2158 let mut projection: Option<crate::projection::SurfaceProjection> = None;
2159 for seed in &seeds {
2160 let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
2161 if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
2162 projection = Some(found);
2163 }
2164 }
2165 let projection = match projection {
2166 Some(projection) => projection,
2167 None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
2168 };
2169 previous_foot = Some([projection.u, projection.v]);
2170 feet.push([projection.u, projection.v]);
2171 let mut u = projection.u;
2172 if let Some(previous) = previous_u {
2173 while u - previous > period * 0.5 {
2174 u -= period;
2175 }
2176 while previous - u > period * 0.5 {
2177 u += period;
2178 }
2179 }
2180 previous_u = Some(u);
2181 proj_u.push(u);
2182 proj_v.push(projection.v);
2183 point3.push(point);
2184 params.push(fraction);
2185 }
2186 project_timer.stop();
2187 if closed_v {
2194 for index in 2..=dense {
2195 while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
2196 proj_v[index] -= period_v;
2197 }
2198 while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
2199 proj_v[index] += period_v;
2200 }
2201 }
2202 while proj_v[0] - proj_v[1] > period_v * 0.5 {
2203 proj_v[0] -= period_v;
2204 }
2205 while proj_v[1] - proj_v[0] > period_v * 0.5 {
2206 proj_v[0] += period_v;
2207 }
2208 }
2209 let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2212 let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
2213 let pinned = |index: usize,
2222 neighbour: usize,
2223 proj_u: &[f64],
2224 proj_v: &[f64]|
2225 -> Result<(f64, f64), KernelRefusal> {
2226 let point = point3[index];
2227 let u_error = surface
2228 .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
2229 .sub(point)
2230 .length();
2231 if closed_v {
2232 let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
2233 let v_error = surface
2234 .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
2235 .sub(point)
2236 .length();
2237 if v_error < u_error {
2238 return Ok((proj_u[index], seam_v));
2239 }
2240 }
2241 Ok((seam_u, proj_v[index]))
2242 };
2243 if start_is_crossing {
2244 let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
2245 proj_u[0] = u;
2246 proj_v[0] = v;
2247 }
2248 if end_is_crossing {
2249 let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
2250 proj_u[dense] = u;
2251 proj_v[dense] = v;
2252 }
2253 let mut shift = 0.0;
2254 while mean_u + shift > u1 + 1e-9 {
2255 shift -= period;
2256 }
2257 while mean_u + shift < u0 - 1e-9 {
2258 shift += period;
2259 }
2260 for u in proj_u.iter_mut() {
2261 *u += shift;
2262 }
2263 if closed_v {
2264 let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2265 let mut shift_v = 0.0;
2266 while mean_v + shift_v > v1 + 1e-9 {
2267 shift_v -= period_v;
2268 }
2269 while mean_v + shift_v < v0 - 1e-9 {
2270 shift_v += period_v;
2271 }
2272 for v in proj_v.iter_mut() {
2273 *v += shift_v;
2274 }
2275 }
2276 let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
2291 let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
2298 let mut track = crate::blend::track_fit::Track {
2299 fractions: params,
2300 uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
2301 feet,
2302 breaks: Vec::new(),
2303 };
2304 let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
2305 let polish_timer = crate::blend::BlendTimer::start("timing.track.polish_ms");
2306 crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
2307 polish_timer.stop();
2308 let crossings_timer = crate::blend::BlendTimer::start("timing.track.crossings_ms");
2309 crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
2310 crossings_timer.stop();
2311 crate::blend::record_count("count.track.breaks", track.breaks.len() as u64);
2312 for &fraction in refine_near {
2313 for _ in 0..2 {
2314 let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
2315 if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], ceiling)?.is_none() {
2316 break;
2317 }
2318 }
2319 }
2320 let fit_timer = crate::blend::BlendTimer::start("timing.track.fit_local_ms");
2327 let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
2328 surface,
2329 &mut track,
2330 &at,
2331 tolerance,
2332 REFINE_ROUNDS,
2333 ceiling,
2334 )?;
2335 fit_timer.stop();
2336 crate::blend::record_count("count.track.fit_sections", fit.sections as u64);
2337 crate::blend::record_count("count.track.local_inserted", inserted as u64);
2338 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2339 eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2340 eprintln!(
2341 "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2342 fit.miss,
2343 fit.samples,
2344 track.fractions.len(),
2345 track.breaks.len(),
2346 fit.on_floor
2347 );
2348 if let Some(worst) = fit.worst {
2353 let parameter = d0 + (d1 - d0) * worst.fraction;
2354 let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2355 .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2356 .map(|pair| [pair[0], pair[1]]).collect();
2357 let mean = (d1 - d0) / spans.len().max(1) as f64;
2358 let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2359 let width = span.map(|index| spans[index][1] - spans[index][0]);
2360 eprintln!(
2361 "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2362 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2363 worst.fraction, parameter, span, spans.len(), width, mean,
2364 worst.bracket[0], worst.bracket[1], worst.interval
2365 );
2366 }
2367 }
2368 if !fit.on_floor {
2375 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2376 for (label, t) in [("start", d0), ("end", d1)] {
2379 if let Ok(point) = piece.evaluate(t) {
2380 let standoff = crate::project_point_to_surface(surface, point).map(|p| p.distance).unwrap_or(f64::NAN);
2381 eprintln!(
2382 "TRACK_FIT closed-chain piece {label} ({:.9}, {:.9}, {:.9}) stands {standoff:.3e} off the carrier",
2383 point.x, point.y, point.z
2384 );
2385 }
2386 }
2387 }
2388 return Err(KernelRefusal::non_convergence(KernelStage::Refine, CHAIN_PCURVE_FLOOR_WHAT, format!(
2389 "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2390 against a floor of {tolerance:.1e}",
2391 crate::blend::PCURVE_OFF_FLOOR,
2392 fit.miss,
2393 fit.samples
2394 )));
2395 }
2396 Ok(fit.curve)
2397}
2398
2399
2400
2401pub(in crate::blend) struct PreparedWall<'a> {
2404 fine: &'a [Vec3],
2405 coarse: Vec<Vec3>,
2406 fine_stencils: StencilPreparation,
2407 coarse_stencils: StencilPreparation,
2408}
2409
2410#[derive(Clone, Copy)]
2411struct WallStart {
2412 curve: (Vec3, Vec3, Vec3),
2413 point_length: f64,
2414 tangent_length: f64,
2415 tangent_squared: f64,
2416}
2417impl WallStart {
2418 fn new(curve: (Vec3, Vec3, Vec3)) -> Self {
2419 Self { point_length: curve.0.length(), tangent_length: curve.1.length(),
2420 tangent_squared: curve.1.dot(curve.1), curve }
2421 }
2422}
2423
2424struct StencilPreparation {
2425 finite: bool,
2426 coefficients: std::cell::RefCell<Vec<Option<[Vec3; 4]>>>,
2427 starts: std::cell::RefCell<Vec<[Option<WallStart>; 2]>>,
2430}
2431impl StencilPreparation {
2432 fn new(nodes: &[Vec3]) -> Self {
2433 Self { finite: wall_nodes_finite(nodes), coefficients: std::cell::RefCell::new(vec![None; nodes.len().saturating_sub(3)]),
2434 starts: std::cell::RefCell::new(vec![[None; 2]; nodes.len().saturating_sub(3)]) }
2435 }
2436 fn start(&self, base: usize, j: usize, curve: &impl Fn(f64) -> (Vec3, Vec3, Vec3)) -> Option<WallStart> {
2437 let offset = j.checked_sub(base)?;
2438 if !(1..=2).contains(&offset) { return None; }
2439 let mut entries = self.starts.borrow_mut();
2440 Some(*entries.get_mut(base)?[offset-1].get_or_insert_with(|| WallStart::new(curve(offset as f64))))
2441 }
2442 fn coefficients(&self, nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2443 let mut entries = self.coefficients.borrow_mut();
2444 *entries[base].get_or_insert_with(|| wall_coefficients(nodes, base))
2445 }
2446}
2447
2448fn wall_nodes_finite(nodes: &[Vec3]) -> bool {
2449 nodes.iter().all(|c| c.x.is_finite() && c.y.is_finite() && c.z.is_finite())
2450}
2451
2452fn wall_coefficients(nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2453 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
2455 let d1 = p1.sub(p0);
2456 let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
2457 let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
2458 [p0, d1, d2, d3]
2459}
2460
2461fn wall_coarse(fine: &[Vec3]) -> Vec<Vec3> {
2462 fine.iter().step_by(2).copied().collect()
2463}
2464
2465fn wall_preparation_on() -> bool {
2466 switch!(live "BREP_WALL_PREPARE", true)
2467}
2468
2469impl<'a> PreparedWall<'a> {
2470 pub(in crate::blend) fn for_interval(fine: &'a [Vec3]) -> Option<Self> {
2471 (wall_preparation_on() && wall_solve_reuse_on()).then(|| Self::new(fine))
2472 }
2473 fn new(fine: &'a [Vec3]) -> Self {
2474 let coarse = wall_coarse(fine);
2475 let fine_stencils = StencilPreparation::new(fine);
2476 let coarse_stencils = StencilPreparation::new(&coarse);
2477 Self { fine, coarse, fine_stencils, coarse_stencils }
2478 }
2479 pub(in crate::blend) fn read(&self, p: Vec3, j: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
2480 let read = |at| declared_wall_probe_prepared(p, self.fine, at, radius, true, Some(self));
2481 match read(j) {
2482 Some(read) => Some(read),
2483 None => match centre_path_node(self.fine, fraction) {
2484 Some(located) if located != j => read(located),
2485 _ => None,
2486 },
2487 }
2488 }
2489}