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(
88 solid: &BrepSolid,
89 segments: &[ChainSegment<'_>],
90 radius: f64,
91 chamfer: bool,
92 name: Option<&str>,
93) -> Result<BrepSolid, KernelRefusal> {
94 let mode = std::env::var("BREP_CHAIN_LADDER_MEMO");
95 let (memo_on, verify) = match mode.as_deref() {
96 Ok("0") => (false, false),
97 Ok("verify") => (true, true),
98 _ => (true, false),
99 };
100 if !memo_on {
101 return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
102 }
103 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
104 let key = ladder_memo_key(segments, radius, chamfer, bar);
105 let cached = LADDER_MEMO.with(|memo| memo.borrow().iter().find(|(stored, _)| *stored == key).map(|(_, refusal)| refusal.clone()));
106 if let (Some(refusal), false) = (&cached, verify) {
107 crate::blend::record_count("count.chain_ladder_memo.hit", 1);
108 return Err(refusal.clone());
109 }
110 let result = blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
111 if let Some(refusal) = cached {
112 let same = matches!(&result, Err(again) if *again == refusal && format!("{again:?}") == format!("{refusal:?}"));
113 assert!(same, "chain ladder memo mismatch: cached {refusal:?}, recomputed {:?}", result.as_ref().err());
114 crate::blend::record_count("count.chain_ladder_memo.verified", 1);
115 return result;
116 }
117 if let Err(refusal) = &result {
118 if matches!(&refusal.class, crate::RefusalClass::NonConvergence { what } if what == "chain_ladder_top") {
119 LADDER_MEMO.with(|memo| {
120 let mut memo = memo.borrow_mut();
121 if memo.len() >= LADDER_MEMO_CAPACITY {
122 memo.remove(0);
123 }
124 memo.push((key, refusal.clone()));
125 });
126 }
127 }
128 result
129}
130
131const LADDER_MEMO_CAPACITY: usize = 32;
133
134thread_local! {
135 static LADDER_MEMO: std::cell::RefCell<Vec<(Vec<u64>, KernelRefusal)>> = const { std::cell::RefCell::new(Vec::new()) };
137}
138
139fn ladder_memo_key(segments: &[ChainSegment<'_>], radius: f64, chamfer: bool, bar: f64) -> Vec<u64> {
141 fn curve(key: &mut Vec<u64>, curve: &NurbsCurve) {
142 key.extend([curve.degree as u64, curve.knots.len() as u64, curve.control_points.len() as u64]);
143 key.extend(curve.knots.iter().map(|knot| knot.to_bits()));
144 for point in &curve.control_points {
145 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
146 }
147 }
148 fn face(key: &mut Vec<u64>, face: &FaceRecord) {
149 let surface = &face.surface;
150 key.extend([face.id, face.same_sense as u64, surface.degree_u as u64, surface.degree_v as u64]);
151 key.extend([surface.knots_u.len() as u64, surface.knots_v.len() as u64, surface.control_points.len() as u64]);
152 key.extend(surface.knots_u.iter().chain(&surface.knots_v).map(|knot| knot.to_bits()));
153 for row in &surface.control_points {
154 key.push(row.len() as u64);
155 for point in row {
156 key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
157 }
158 }
159 key.push(face.loops.len() as u64);
160 for loop_record in &face.loops {
161 key.extend([loop_record.id, loop_record.coedges.len() as u64]);
162 for coedge in &loop_record.coedges {
163 key.extend([coedge.id, coedge.edge_id, coedge.forward as u64]);
164 curve(key, &coedge.pcurve);
165 }
166 }
167 }
168 let mut key = vec![
169 radius.to_bits(),
170 chamfer as u64,
171 bar.to_bits(),
172 crate::blend::station_refinement_on() as u64,
173 crate::blend::chain_arc_stations_on() as u64,
174 wall_verdict_on() as u64,
175 local_continue_on() as u64,
176 segments.len() as u64,
177 ];
178 for segment in segments {
179 let edge = segment.edge;
180 key.extend([edge.id, edge.start_vertex_id, edge.end_vertex_id, edge.t0.to_bits(), edge.t1.to_bits(), segment.forward as u64]);
181 curve(&mut key, &edge.curve);
182 for mate in [&segment.first, &segment.second] {
183 key.extend([mate.coedge.id, mate.coedge.edge_id, mate.coedge.forward as u64, mate.loop_index as u64, mate.rho.to_bits()]);
184 curve(&mut key, &mate.coedge.pcurve);
185 face(&mut key, mate.face);
186 }
187 }
188 key
189}
190fn blend_closed_smooth_chain_unmemoized(
203 solid: &BrepSolid,
204 segments: &[ChainSegment<'_>],
205 radius: f64,
206 chamfer: bool,
207 name: Option<&str>,
208) -> Result<BrepSolid, KernelRefusal> {
209 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
210 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
213 let mut trend: Vec<(usize, f64)> = Vec::new();
215 let first = loop {
216 let per_segment = budget.rung;
217 let rung = build_closed_smooth_chain(
218 solid,
219 segments,
220 radius,
221 chamfer,
222 name,
223 &budget,
224 FoldPolicy::Refuse,
225 bar,
226 );
227 match rung {
228 Ok(Rung::TooCoarse { deviation }) => {
229 crate::blend::carve::carve_trace(format_args!(
230 "blend chain: {} leaves the rails \
231 {deviation:.6e} from the rolling ball, against {bar:.3e}",
232 budget.describe()
233 ));
234 trend.push((per_segment, deviation));
235 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
236 let (best_rung, best, listing) = ladder_trend(&trend);
237 return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
238 "blend: the chain's ladder reached its top, {per_segment} stations per segment, \
239 and no rung passed — its best, {best_rung} stations per segment, left this \
240 collar's fitted rails {best:.6e} from the rolling ball's own contacts against \
241 a bar of {bar:.3e} (rung by rung: {listing}), so there is no wall accurate \
242 enough to build"
243 )));
244 }
245 budget = budget.doubled();
246 }
247 Ok(Rung::StationFailed(failure)) => {
252 crate::blend::carve::carve_trace(format_args!(
253 "blend chain: {per_segment} stations per segment: segment {} position {} \
254 (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
255 failure.segment, failure.position, failure.t, failure.from_position,
256 failure.from_t, failure.seed_uv, failure.error
257 ));
258 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
259 break Err(failure.error);
260 }
261 budget = budget.doubled();
262 }
263 other => break other,
264 }
265 };
266 if let Err(error) = &first {
267 crate::blend::carve::carve_trace(format_args!(
268 "blend chain: the march at {} refused: {error}", budget.describe()
269 ));
270 }
271 match first {
272 Ok(Rung::Built(built)) => Ok(built),
273 Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
274 Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
275 }
276 Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
282 Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
285 carve_ladder(solid, segments, radius, name, bar)
286 .map_err(|carve_error| {
294 if crate::blend::is_wall_fold(&carve_error) {
295 carve_error
296 } else {
297 crate::blend::fold::wall_fold_refusal(
298 Vec::new(),
299 format!(
300 "{} {} fits this edge: the wall folds, and the fold band could not \
301 be carved — {carve_error}",
302 crate::blend::WALL_FOLDS,
303 radius.abs()
304 ),
305 )
306 }
307 })
308 }
309 Err(error) => Err(error),
310 }
311}
312
313fn carve_ladder(
323 solid: &BrepSolid,
324 segments: &[ChainSegment<'_>],
325 radius: f64,
326 name: Option<&str>,
327 bar: f64,
328) -> Result<BrepSolid, KernelRefusal> {
329 let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
330 loop {
331 let per_segment = budget.rung;
332 match build_closed_smooth_chain(
333 solid,
334 segments,
335 radius,
336 false,
337 name,
338 &budget,
339 FoldPolicy::Carve,
340 bar,
341 )? {
342 Rung::Built(built) => {
343 crate::blend::carve::carve_trace(format_args!(
344 "blend carve: {} is the first rung inside \
345 {bar:.3e}",
346 budget.describe()
347 ));
348 return Ok(built);
349 }
350 Rung::StationFailed(failure) => return Err(failure.error),
353 Rung::TooCoarse { deviation } => {
354 crate::blend::carve::carve_trace(format_args!(
355 "blend carve: {} leaves the rails \
356 {deviation:.6e} from the rolling ball, against {bar:.3e}",
357 budget.describe()
358 ));
359 if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
360 return Err(crate::blend::fold::wall_fold_refusal(
361 Vec::new(),
362 format!(
363 "{} {} fits this edge: the wall folds, and at {per_segment} stations \
364 per segment — the top of the carve's ladder — its fitted rails are \
365 still {deviation:.6e} from the rolling ball's own contacts against a \
366 bar of {bar:.3e}, so there is no wall accurate enough to carve",
367 crate::blend::WALL_FOLDS,
368 radius.abs()
369 ),
370 ));
371 }
372 budget = budget.doubled();
373 }
374 }
375 }
376}
377
378enum Rung {
381 Built(BrepSolid),
382 TooCoarse { deviation: f64 },
383 StationFailed(ChainStationFailure),
386}
387
388fn rail_deviation(
412 samples: &[ChainSample],
413 cr: &NurbsCurve,
414 cs: &NurbsCurve,
415 fit_low: f64,
416 fit_range: f64,
417 budget: &ChainBudget,
418 bar: f64,
419) -> Result<(f64, Vec<usize>), KernelRefusal> {
420 let mut order: Vec<usize> = (0..samples.len())
424 .filter(|&index| (0..=budget.count(samples[index].segment) as isize).contains(&samples[index].position))
425 .collect();
426 order.sort_by(|&a, &b| {
427 let (a, b) = (&samples[a], &samples[b]);
428 a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
429 });
430 let next_of: std::collections::HashMap<usize, f64> = order
431 .windows(2)
432 .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
433 .map(|pair| (pair[0], samples[pair[1]].parameter))
434 .collect();
435 let mut failing = Vec::new();
436 let closest = closest_on_rail;
437 let mut worst: f64 = 0.0;
438 let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
441 for &index in &order {
442 let sample = &samples[index];
443 let Some([p1, p2]) = sample.midpoint else {
444 continue;
445 };
446 let Some(&next) = next_of.get(&index) else {
447 continue;
448 };
449 let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
450 let seed = (global - fit_low) / fit_range;
451 let reach = (next - sample.parameter).abs() / fit_range;
452 let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
453 worst = worst.max(miss1).max(miss2);
454 if miss1.max(miss2) > bar {
455 failing.push(index);
456 }
457 if per_segment_worst.len() <= sample.segment {
458 per_segment_worst.resize(
459 sample.segment + 1,
460 (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
461 );
462 }
463 if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
464 per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
465 }
466 }
467 for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
468 crate::blend::carve::carve_trace(format_args!(
469 " rail miss segment {segment}: worst {miss:.6e} after station {position} \
470 (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
471 p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
472 ));
473 }
474 Ok((worst, failing))
475}
476
477const PIECE_TRACK_CEILING: usize = 4096;
481
482fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
485 let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
486 let distance = |u: f64| -> Result<f64, KernelRefusal> {
487 Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
488 };
489 let mut u = seed.clamp(low, high);
490 let mut best = distance(u)?;
491 for _ in 0..40 {
492 let derivatives = curve.derivatives_extended(u, 2).or_refuse(KernelStage::Refine, "derivatives_extended")?;
493 let offset = derivatives[0].sub(target);
494 let slope = offset.dot(derivatives[1]);
495 let curvature = derivatives[1].dot(derivatives[1]) + offset.dot(derivatives[2]);
496 let mut step = if curvature > 0.0 {
499 (slope / curvature).clamp(-reach, reach)
500 } else {
501 reach.copysign(slope)
502 };
503 let mut trial = (u - step).clamp(low, high);
504 let mut trial_distance = distance(trial)?;
505 for _ in 0..30 {
506 if trial_distance < best {
507 break;
508 }
509 step *= 0.5;
510 trial = (u - step).clamp(low, high);
511 trial_distance = distance(trial)?;
512 }
513 if !(trial_distance < best) {
514 break;
515 }
516 let moved = (trial - u).abs();
517 u = trial;
518 best = trial_distance;
519 if moved < 1e-15 {
520 break;
521 }
522 }
523 Ok(best)
524}
525
526pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
529
530fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
534 (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
535}
536
537#[allow(clippy::too_many_arguments)]
543pub(super) fn dense_rail_misses(
544 contacts: &crate::blend::chain::DenseContacts,
545 samples: &[ChainSample],
546 cr: &NurbsCurve,
547 cs: &NurbsCurve,
548 low: f64,
549 range: f64,
550 wrap: bool,
551) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
552 let mut readable = Vec::with_capacity(contacts.solved.len());
553 let mut unreadable = Vec::new();
554 for (left, next, fraction, p1, p2) in &contacts.solved {
555 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
556 let global = from + fraction * (to - from);
560 let global = if wrap { global.rem_euclid(1.0) } else { global };
561 let seed = (global - low) / range;
562 let reach = (to - from).abs() / range;
563 let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
566 let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
567 match rail_pair_miss(cr_miss, cs_miss) {
568 Some(miss) => readable.push((*left, miss)),
569 None => unreadable.push(*left),
570 }
571 }
572 Ok((readable, unreadable))
573}
574
575pub(super) fn combine_rail_verdict(
581 halfway: (f64, Vec<usize>),
582 dense: &[(usize, f64)],
583 unsolved: &[usize],
584 bar: f64,
585) -> (f64, Vec<usize>) {
586 let (mut worst, mut failing) = halfway;
587 for &(left, miss) in dense {
588 if !miss.is_finite() {
589 failing.push(left);
590 continue;
591 }
592 if miss > bar {
593 failing.push(left);
594 }
595 worst = worst.max(miss);
596 }
597 failing.extend_from_slice(unsolved);
598 failing.sort_unstable();
599 failing.dedup();
600 (worst, failing)
601}
602
603pub(super) fn model_rounds_limit() -> usize {
610 REFINE_ROUNDS
611}
612
613pub(super) fn model_request_on() -> bool {
614 wall_verdict_on()
615}
616
617
618pub(super) fn wall_verdict_on() -> bool {
621 true
622}
623
624
625
626
627
628pub(super) const WALL_SECTION_PROBES: usize = 8;
631
632pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
653 let j = locate_foot(p, nodes, j, radius)?;
654 sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
655}
656
657const FOOT_TRAVEL: usize = 4;
660
661fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
671 locate_foot_read(p, nodes, j, radius).map(|(node, _)| node)
672}
673
674fn locate_foot_read(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(usize, f64)> {
677 locate_foot_read_prepared(p, nodes, j, radius, None)
678}
679
680fn locate_foot_read_prepared(p: Vec3, nodes: &[Vec3], j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Option<(usize, f64)> {
681 let n = nodes.len();
682 if n < 4 || j == 0 || j + 1 >= n {
683 return None;
684 }
685 let mut tried: Vec<usize> = Vec::new();
686 let mut at = j;
687 loop {
688 let base = (at - 1).min(n - 4);
689 match stencil_solve_prepared(p, nodes, base, at, radius, prepared) {
690 Solve::Read(distance) => return Some((at, distance)),
691 Solve::Unread => return None,
692 Solve::Window(s) => {
693 tried.push(at);
694 let next = (base as f64 + s).round();
695 let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
696 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
697 wall_decline(|| format!(
698 "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
699 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
700 ));
701 return None;
702 }
703 at = next;
704 }
705 }
706 }
707}
708
709fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
712 match stencil_solve(p, nodes, base, j, radius) {
713 Solve::Read(distance) => Some(distance),
714 Solve::Window(_) | Solve::Unread => None,
715 }
716}
717
718enum Solve {
721 Read(f64),
722 Window(f64),
723 Unread,
724}
725
726thread_local! {
727 static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
732}
733
734fn wall_decline(reason: impl FnOnce() -> String) {
735 WALL_DECLINES.with(|log| {
736 if let Some(lines) = log.borrow_mut().as_mut() {
737 lines.push(reason());
738 }
739 });
740}
741
742pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
750 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
751 let read = declared_wall_probe(p, fine, jf, radius);
752 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
753 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
754}
755
756fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
757 stencil_solve_prepared(p, nodes, base, j, radius, None)
758}
759
760fn stencil_solve_prepared(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Solve {
761 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
762 let n = nodes.len();
763 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) {
764 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
765 return Solve::Unread;
766 }
767 let [p0, d1, d2, d3] = prepared.map_or_else(|| wall_coefficients(nodes, base), |cache| cache.coefficients(nodes, base));
768 let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
769 let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
770 let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
771 let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
772 (c, c1, c2)
773 };
774 let start = j as f64 - base as f64;
775 let mut s = start;
776 let mut converged = false;
777 for _ in 0..50 {
778 let (c, c1, c2) = curve(s);
779 let w = c.sub(p);
780 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
781 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
782 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
783 return Solve::Unread;
784 }
785 let stationary = g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
798 if stationary {
799 converged = true;
800 break;
801 }
802 let step = g / slope;
803 let before = s;
804 s -= step;
805 if !((s - start).abs() <= 1.0) {
806 wall_decline(|| {
807 let (c, c1, _) = curve(s);
810 let g_after = c.sub(p).dot(c1);
811 let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
812 format!(
813 "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
814 g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
815 (s - start).abs() - 1.0
816 )
817 });
818 return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
819 }
820 if step.abs() <= 1e-14 * (1.0 + s.abs()) {
821 converged = true;
822 break;
823 }
824 }
825 if !converged {
826 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
827 return Solve::Unread;
828 }
829 let distance = curve(s).0.sub(p).length() - radius;
830 if distance.is_finite() {
831 Solve::Read(distance)
832 } else {
833 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
834 Solve::Unread
835 }
836}
837
838pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
847 let n = nodes.len();
848 if n < 5 || j == 0 || j + 1 >= n {
849 return None;
850 }
851 let j = locate_foot(p, nodes, j, radius)?;
852 let first = (j - 1).min(n - 4);
853 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
854 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
855 let a = sweep_residual_on(p, nodes, first, j, radius)?;
856 let b = sweep_residual_on(p, nodes, second, j, radius)?;
857 Some((a.abs().max(b.abs()), (a - b).abs()))
858}
859
860fn wall_solve_reuse_on() -> bool {
873 std::env::var("BREP_WALL_SOLVE_REUSE").map_or(true, |value| value != "0" && value != "off")
874}
875
876pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
877 declared_wall_probe_with(p, fine, jf, radius, wall_solve_reuse_on())
878}
879
880fn declared_wall_probe_with(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool) -> Option<(f64, f64)> {
881 declared_wall_probe_prepared(p, fine, jf, radius, reuse, None)
882}
883
884fn declared_wall_probe_prepared(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool, prepared: Option<&PreparedWall<'_>>) -> Option<(f64, f64)> {
885 if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
886 wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
887 return None;
888 }
889 let (reading, spread, at_fine) = if reuse {
890 let read = (|| {
891 let n = fine.len();
892 let (j, a) = locate_foot_read_prepared(p, fine, jf, radius, prepared.map(|c| &c.fine_stencils))?;
893 let first = (j - 1).min(n - 4);
894 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
895 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
896 let b = match stencil_solve_prepared(p, fine, second, j, radius, prepared.map(|c| &c.fine_stencils)) {
897 Solve::Read(distance) => distance,
898 _ => return None,
899 };
900 Some((a.abs().max(b.abs()), (a - b).abs(), a))
901 })();
902 let Some(read) = read else {
903 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
904 return None;
905 };
906 read
907 } else {
908 let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
909 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
910 return None;
911 };
912 let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
913 wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
914 return None;
915 };
916 (reading, spread, at_fine)
917 };
918 let coarse_owned;
919 let coarse = match prepared {
920 Some(cache) => cache.coarse.as_slice(),
921 None => { coarse_owned = wall_coarse(fine); &coarse_owned }
922 };
923 let at_coarse = if reuse {
924 locate_foot_read_prepared(p, coarse, jf / 2, radius, prepared.map(|c| &c.coarse_stencils)).map(|(_, distance)| distance)
925 } else {
926 sweep_residual(p, coarse, jf / 2, radius)
927 };
928 let Some(at_coarse) = at_coarse else {
929 wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
930 return None;
931 };
932 let uncertainty = spread.max((at_fine - at_coarse).abs());
933 uncertainty.is_finite().then_some((reading, uncertainty))
934}
935
936pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
952 let n = fine.len();
953 if n < 7 || !(fraction.is_finite()) {
954 return None;
955 }
956 let mut cumulative = Vec::with_capacity(n);
957 cumulative.push(0.0_f64);
958 for pair in fine.windows(2) {
959 let step = pair[1].sub(pair[0]).length();
960 cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
961 }
962 let total = *cumulative.last()?;
963 if !(total.is_finite() && total > 0.0) {
964 return None;
965 }
966 let target = fraction.clamp(0.0, 1.0) * total;
967 let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
968 let even = 2 * ((nearest as f64 / 2.0).round() as usize);
969 Some(even.clamp(2, (n - 3) & !1))
970}
971
972pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
979 let first = match declared_wall_probe_why(p, fine, jf, radius) {
980 Ok(read) => return Ok(read),
981 Err(why) => why,
982 };
983 match centre_path_node(fine, fraction) {
984 Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
985 .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
986 _ => Err(first),
987 }
988}
989
990pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
992 match declared_wall_probe(p, fine, jf, radius) {
993 Some(read) => Some(read),
994 None => match centre_path_node(fine, fraction) {
995 Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
996 _ => None,
997 },
998 }
999}
1000
1001pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
1013 if jf % 2 != 0
1014 || fine.len() < 9
1015 || fine.len() % 2 == 0
1016 || ts.len() != fine.len()
1017 || !ts.iter().all(|t| t.is_finite())
1018 || !ts.windows(2).all(|pair| pair[1] > pair[0])
1019 {
1020 return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
1021 }
1022 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
1023 let read = (|| -> Option<(f64, f64)> {
1024 let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
1025 wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
1026 None
1027 })?;
1028 let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
1029 wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
1030 None
1031 })?;
1032 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
1033 let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
1034 let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
1035 wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
1036 None
1037 })?;
1038 let uncertainty = spread.max((at_fine - at_coarse).abs());
1039 uncertainty.is_finite().then_some((reading, uncertainty))
1040 })();
1041 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
1042 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
1043}
1044
1045fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
1046 let n = nodes.len();
1047 if n < 5 || j == 0 || j + 1 >= n {
1048 return None;
1049 }
1050 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1051 let first = (j - 1).min(n - 4);
1052 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
1053 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
1054 let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
1055 let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
1056 Some((a.abs().max(b.abs()), (a - b).abs()))
1057}
1058
1059fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
1060 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1061 match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
1062 Solve::Read(distance) => Some(distance),
1063 Solve::Window(_) | Solve::Unread => None,
1064 }
1065}
1066
1067fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
1071 let n = nodes.len();
1072 if n < 4 || j == 0 || j + 1 >= n {
1073 return None;
1074 }
1075 let mut tried: Vec<usize> = Vec::new();
1076 let mut at = j;
1077 loop {
1078 let base = (at - 1).min(n - 4);
1079 match stencil_solve_at(p, nodes, ts, base, at, radius) {
1080 Solve::Read(_) => return Some(at),
1081 Solve::Unread => return None,
1082 Solve::Window(t) => {
1083 tried.push(at);
1084 let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
1085 let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
1086 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
1087 wall_decline(|| format!(
1088 "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
1089 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
1090 ));
1091 return None;
1092 }
1093 at = next;
1094 }
1095 }
1096 }
1097}
1098
1099fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
1104 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
1105 let n = nodes.len();
1106 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)) {
1107 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
1108 return Solve::Unread;
1109 }
1110 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
1111 let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
1112 let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
1113 let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
1114 let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
1115 let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
1116 let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
1117 let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
1118 let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
1119 let (a, b, c) = (t - t0, t - t1, t - t2);
1120 let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
1121 let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
1122 let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
1123 (point, first, second)
1124 };
1125 let (low, high) = (ts[j - 1], ts[j + 1]);
1126 let start = ts[j];
1127 let mut t = start;
1128 let mut converged = false;
1129 for _ in 0..50 {
1130 let (c, c1, c2) = curve(t);
1131 let w = c.sub(p);
1132 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
1133 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
1134 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
1135 return Solve::Unread;
1136 }
1137 if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
1138 converged = true;
1139 break;
1140 }
1141 let step = g / slope;
1142 let before = t;
1143 t -= step;
1144 if !(t >= low && t <= high) {
1145 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
1146 return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
1147 }
1148 let node = 0.5 * (high - low);
1153 if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
1154 converged = true;
1155 break;
1156 }
1157 }
1158 if !converged {
1159 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
1160 return Solve::Unread;
1161 }
1162 let distance = curve(t).0.sub(p).length() - radius;
1163 if distance.is_finite() {
1164 Solve::Read(distance)
1165 } else {
1166 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
1167 Solve::Unread
1168 }
1169}
1170
1171pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
1176 readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
1177}
1178
1179#[allow(clippy::too_many_arguments)]
1189pub(super) fn wall_interior_misses(
1190 contacts: &crate::blend::chain::DenseContacts,
1191 samples: &[ChainSample],
1192 surface: &crate::NurbsSurface,
1193 low: f64,
1194 range: f64,
1195 radius: f64,
1196 subdivisions: usize,
1197 probes: usize,
1198 wrap: bool,
1199) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1200 use std::collections::BTreeMap;
1201 let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1204 for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1205 let entry = nodes.entry(*left).or_insert_with(|| {
1206 let mut row = vec![None; subdivisions + 1];
1207 row[0] = Some(samples[*left].station.center);
1208 row[subdivisions] = Some(samples[*next].station.center);
1209 (*next, row)
1210 });
1211 let step = (fraction * subdivisions as f64).round() as usize;
1212 if (1..subdivisions).contains(&step) {
1213 entry.1[step] = contacts.centers.get(index).copied().flatten();
1214 }
1215 }
1216 let mut readable = Vec::new();
1217 let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1218 for (&left, (next, row)) in &nodes {
1219 let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1220 unreadable.push(left);
1221 continue;
1222 };
1223 let (from, to) = (samples[left].parameter, samples[*next].parameter);
1224 let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1225 let mut read = true;
1226 let prepared = PreparedWall::for_interval(¢res);
1227 'probes: for step in (2..subdivisions).step_by(2) {
1228 let global = from + (step as f64 / subdivisions as f64) * (to - from);
1229 let global = if wrap { global.rem_euclid(1.0) } else { global };
1230 let u = (global - low) / range;
1231 for k in 1..probes {
1232 let v = k as f64 / probes as f64;
1233 let reading = surface.evaluate(u, v).ok().and_then(|p| {
1234 match &prepared {
1235 Some(cache) => cache.read(p, step, step as f64 / subdivisions as f64, radius),
1236 None => declared_wall_probe_located(p, ¢res, step, step as f64 / subdivisions as f64, radius),
1237 }
1238 });
1239 match reading {
1240 Some((distance, gap)) => {
1241 declared = declared.max(distance + gap);
1242 worst = worst.max(distance);
1243 spread = spread.max(gap);
1244 }
1245 None => {
1246 read = false;
1247 break 'probes;
1248 }
1249 }
1250 }
1251 }
1252 if read {
1253 readable.push((left, declared, worst, spread));
1254 } else {
1255 unreadable.push(left);
1256 }
1257 }
1258 unreadable.sort_unstable();
1259 unreadable.dedup();
1260 (readable, unreadable)
1261}
1262
1263
1264fn build_closed_smooth_chain(
1265 solid: &BrepSolid,
1266 segments: &[ChainSegment<'_>],
1267 radius: f64,
1268 chamfer: bool,
1269 name: Option<&str>,
1270 budget: &ChainBudget,
1271 fold: FoldPolicy,
1272 bar: f64,
1273) -> Result<Rung, KernelRefusal> {
1274 let samples = match march_chain_typed(segments, radius, false, budget, fold) {
1275 Ok(samples) => samples,
1276 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1277 Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1280 Err(ChainMarchError::Other(error)) => return Err(error),
1281 };
1282 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1283 build_closed_chain_from_samples(solid, segments, radius, chamfer, name, budget, fold, bar, samples, 0)
1284}
1285
1286
1287thread_local! {
1288 static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1292}
1293
1294fn local_continue_on() -> bool {
1300 std::env::var("BREP_CHAIN_LOCAL_CONTINUE").as_deref() == Ok("1")
1301}
1302
1303const LOCAL_CONTINUE_MAX_ROUNDS: usize = 4 * REFINE_ROUNDS;
1305
1306const LOCAL_CONTINUE_STALL_ROUNDS: usize = 3;
1308
1309thread_local! {
1310 static CHAIN_LOCAL_BEST: std::cell::Cell<(f64, usize)> = const { std::cell::Cell::new((f64::INFINITY, 0)) };
1312}
1313
1314fn local_rounds_spent(round: usize, deviation: f64) -> bool {
1320 if round == 0 {
1321 CHAIN_LOCAL_BEST.with(|best| best.set((f64::INFINITY, 0)));
1322 return false;
1323 }
1324 let (_, best_round) = CHAIN_LOCAL_BEST.with(|best| {
1325 let (reading, at) = best.get();
1326 if deviation < reading {
1327 best.set((deviation, round));
1328 (deviation, round)
1329 } else {
1330 (reading, at)
1331 }
1332 });
1333 round >= LOCAL_CONTINUE_MAX_ROUNDS || round - best_round >= LOCAL_CONTINUE_STALL_ROUNDS
1334}
1335#[allow(clippy::too_many_arguments)]
1346fn build_closed_chain_from_samples(
1347 solid: &BrepSolid,
1348 segments: &[ChainSegment<'_>],
1349 radius: f64,
1350 chamfer: bool,
1351 name: Option<&str>,
1352 budget: &ChainBudget,
1353 fold: FoldPolicy,
1354 bar: f64,
1355 mut samples: Vec<ChainSample>,
1356 round: usize,
1357) -> Result<Rung, KernelRefusal> {
1358 let per_segment = budget.rung;
1361 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1362
1363 let mut global: Vec<(f64, &ChainSample)> = samples
1365 .iter()
1366 .filter(|sample| (0..budget.count(sample.segment) as isize).contains(&sample.position))
1367 .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1368 .collect();
1369 global.sort_by(|a, b| a.0.total_cmp(&b.0));
1370 let count = global.len();
1371 let degree = FIT_DEGREE;
1372 let mut extended_params = Vec::new();
1373 let mut samples_cr = Vec::new();
1374 let mut samples_mid = Vec::new();
1375 let mut samples_cs = Vec::new();
1376 let mut push = |station: &Station, parameter: f64| {
1377 extended_params.push(parameter);
1378 samples_cr.push(Vec4::from_point(station.p1, 1.0));
1379 samples_cs.push(Vec4::from_point(station.p2, 1.0));
1380 samples_mid.push(Vec4 {
1381 x: station.apex.x * station.weight,
1382 y: station.apex.y * station.weight,
1383 z: station.apex.z * station.weight,
1384 w: station.weight,
1385 });
1386 };
1387 for offset in (1..=degree).rev() {
1388 let (parameter, sample) = &global[count - offset];
1389 push(&sample.station, parameter - 1.0);
1390 }
1391 for (parameter, sample) in &global {
1392 push(&sample.station, *parameter);
1393 }
1394 push(&global[0].1.station, global[0].0 + 1.0);
1397 for offset in 1..=degree {
1398 let (parameter, sample) = &global[offset];
1399 push(&sample.station, parameter + 1.0);
1400 }
1401 let low = extended_params[0];
1402 let high = *extended_params.last().unwrap();
1403 let range = high - low;
1404 let normalized: Vec<f64> = extended_params
1405 .iter()
1406 .map(|parameter| (parameter - low) / range)
1407 .collect();
1408 let seam_low = (0.0 - low) / range;
1409 let seam_high = (1.0 - low) / range;
1410 let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1411 let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1412 let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1413 let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1414 Ok(middle)
1415 };
1416 let cr = fit_row(&samples_cr)?;
1417 let cs = fit_row(&samples_cs)?;
1418 let mid = if chamfer {
1419 None
1420 } else {
1421 Some(fit_row(&samples_mid)?)
1422 };
1423 let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1424 let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1425
1426 let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, budget, bar)?;
1439 let stencil = if round > 0 {
1451 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, LOCAL_RAIL_STENCIL)?)
1452 } else {
1453 None
1454 };
1455 let wall_stencil = if !chamfer && wall_verdict_on() {
1460 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1461 } else {
1462 None
1463 };
1464 let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1465 let contacts = stencil.as_ref().expect("the local path solves the stencil");
1466 let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1467 let reads = dense.len() + unreadable.len() + contacts.unsolved;
1468 let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1469 let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1470 crate::blend::carve::carve_trace(format_args!(
1471 "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1472 {} failing interval(s), {} unsolved contact(s)",
1473 verdict.0, verdict.1.len(), contacts.unsolved
1474 ));
1475 (verdict.0, verdict.1, reads)
1476 } else {
1477 (halfway, halfway_failing, 0)
1478 };
1479 let mut _wall_worst = 0.0_f64;
1485 let mut _wall_spread = 0.0_f64;
1486 let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1487 if let Some(contacts) = wall_stencil.as_ref() {
1488 let (wall, unreadable) =
1489 wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1490 wall_readings = wall.clone();
1491 for &(left, declared, _, uncertainty) in &wall {
1493 if declared > bar {
1494 failing.push(left);
1495 }
1496 _wall_worst = _wall_worst.max(declared);
1497 _wall_spread = _wall_spread.max(uncertainty);
1498 }
1499 failing.extend_from_slice(&unreadable);
1500 failing.sort_unstable();
1501 failing.dedup();
1502 deviation = deviation.max(_wall_worst);
1503 crate::blend::carve::carve_trace(format_args!(
1504 "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1505 interval(s), {} unreadable, bar {bar:.3e}",
1506 wall.len(),
1507 unreadable.len()
1508 ));
1509 }
1510 if deviation > bar || !failing.is_empty() {
1511 let in_segment = |segment: usize| {
1512 samples
1513 .iter()
1514 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1515 .count()
1516 };
1517 let within_ceiling = (0..segments.len()).all(|segment| {
1518 in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1519 <= CHAIN_CARVE_MAX_PER_SEGMENT
1520 });
1521 crate::blend::carve::carve_trace(format_args!(
1522 "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1523 worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1524 samples.len(),
1525 failing.len(),
1526 if within_ceiling { "ok" } else { "reached" }
1527 ));
1528 let rounds_spent = if local_continue_on() {
1530 local_rounds_spent(round, deviation)
1531 } else {
1532 round >= REFINE_ROUNDS
1533 };
1534 if fold != FoldPolicy::Refuse || rounds_spent || !within_ceiling {
1536 return Ok(Rung::TooCoarse { deviation });
1537 }
1538 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1540 match insert_chain_stations(segments, radius, budget, &mut samples, &failing) {
1541 Ok(()) => {}
1542 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1543 Err(ChainMarchError::Incoherent(failure)) => {
1547 crate::blend::carve::carve_trace(format_args!(
1548 "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1549 at t {:.9} does not reach its retained right neighbour: {}",
1550 failure.segment, failure.t, failure.error
1551 ));
1552 return Ok(Rung::TooCoarse { deviation });
1553 }
1554 Err(ChainMarchError::Other(error)) => return Err(error),
1555 }
1556 assign_chain_parameters(&mut samples, segments.len(), budget, false);
1557 return build_closed_chain_from_samples(
1558 solid, segments, radius, chamfer, name, budget, fold, bar, samples, round + 1,
1559 );
1560 }
1561 crate::blend::carve::carve_trace(format_args!(
1562 "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1563 from the rolling ball, inside {bar:.3e}"
1564 ));
1565 let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1579 let stations_now = (0..segments.len())
1580 .map(|segment| {
1581 samples
1582 .iter()
1583 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1584 .count()
1585 })
1586 .max()
1587 .unwrap_or(0);
1588 let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1589 let mut request_detail: Option<String> = None;
1593 let mut request_component = crate::MeasuredComponent::Wall;
1594 let mut request_unread = 0usize;
1595 if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1596 let wall_missing = wall_model_missing(&wall_readings, model);
1597 let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1605 Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1606 None => (Vec::new(), Vec::new()),
1607 };
1608 let rail_over: std::collections::BTreeSet<usize> =
1610 rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1611 let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1612 let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1613 let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1614 missing.sort_unstable();
1615 missing.dedup();
1616 if !missing.is_empty() {
1617 let wall_short = wall_readings
1619 .iter()
1620 .filter(|reading| wall_missing.contains(&reading.0))
1621 .map(|&(_, declared, _, _)| declared)
1622 .fold(0.0_f64, f64::max);
1623 let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1624 let short = wall_short.max(rail_short);
1625 request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1626 request_unread = rail_unread_intervals.len();
1627 request_detail = Some(format!(
1628 "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1629 wall {wall_short:.3e} over {} interval(s)",
1630 rail_over.len(),
1631 rail_unread_intervals.len(),
1632 wall_missing.len()
1633 ));
1634 let within_ceiling = (0..segments.len()).all(|segment| {
1635 let present = samples
1636 .iter()
1637 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1638 .count();
1639 present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1640 <= CHAIN_CARVE_MAX_PER_SEGMENT
1641 });
1642 crate::blend::carve::carve_trace(format_args!(
1643 "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1644 (worst reading + spread {short:.3e}), ceiling {}",
1645 missing.len(),
1646 if within_ceiling { "ok" } else { "reached" }
1647 ));
1648 unmet = if round >= model_rounds_limit() {
1649 Some((short, crate::BudgetReason::RoundsSpent))
1650 } else if !within_ceiling {
1651 Some((short, crate::BudgetReason::StationCeiling))
1652 } else {
1653 let mut refined = samples.clone();
1654 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1656 match insert_chain_stations(segments, radius, budget, &mut refined, &missing) {
1657 Ok(()) => {
1658 assign_chain_parameters(&mut refined, segments.len(), budget, false);
1659 match build_closed_chain_from_samples(
1660 solid, segments, radius, chamfer, name, budget, fold, bar, refined, round + 1,
1661 ) {
1662 Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1663 _ => Some((short, crate::BudgetReason::Rejected)),
1664 }
1665 }
1666 Err(_) => {
1667 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1669 Some((short, crate::BudgetReason::Incoherent))
1670 }
1671 }
1672 };
1673 }
1674 }
1675 if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1680 let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, subdivisions.max(2))?;
1681 let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1682 let mut unreadable = 0usize;
1683 for (left, next, fraction, p1, p2) in &contacts.solved {
1684 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1685 let global = (from + fraction * (to - from)).rem_euclid(1.0);
1686 let seed = (global - low) / range;
1687 let reach = (to - from).abs() / range;
1688 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)) {
1689 Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1690 Some(_) => {}
1691 None => unreadable += 1,
1692 }
1693 }
1694 eprintln!(
1695 "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1696 worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1697 contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1698 );
1699 }
1700 let crease = if fold == FoldPolicy::Carve {
1701 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1711 let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1712 return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1713 "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"
1714 )));
1715 };
1716 Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1717 } else {
1718 None
1719 };
1720
1721 let single_face = |side: usize| -> bool {
1725 let first_id = segments[0].mate(side).face.id;
1726 segments
1727 .iter()
1728 .all(|segment| segment.mate(side).face.id == first_id)
1729 };
1730 let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1731 for side in 0..2 {
1732 if !single_face(side) {
1733 continue;
1734 }
1735 let mut row = Vec::new();
1736 let select = |station: &Station| -> [f64; 2] {
1737 if side == 0 {
1738 station.uv1
1739 } else {
1740 station.uv2
1741 }
1742 };
1743 let mut push_uv = |station: &Station| {
1744 let uv = select(station);
1745 row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1746 };
1747 for offset in (1..=degree).rev() {
1748 push_uv(&global[count - offset].1.station);
1749 }
1750 for (_, sample) in &global {
1751 push_uv(&sample.station);
1752 }
1753 push_uv(&global[0].1.station);
1754 for offset in 1..=degree {
1755 push_uv(&global[offset].1.station);
1756 }
1757 whole_pcurves[side] = Some(fit_row(&row)?);
1758 }
1759
1760 let mut pcurves1 = Vec::with_capacity(segments.len());
1763 let mut pcurves2 = Vec::with_capacity(segments.len());
1764 for segment in 0..segments.len() {
1765 let mut window: Vec<&ChainSample> = samples
1766 .iter()
1767 .filter(|sample| sample.segment == segment)
1768 .collect();
1769 window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1770 let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1771 let low = params[0];
1772 let high = *params.last().unwrap();
1773 let normalized: Vec<f64> = params
1774 .iter()
1775 .map(|parameter| (parameter - low) / (high - low))
1776 .collect();
1777 let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1778 let points: Vec<Vec4> = window
1779 .iter()
1780 .map(|sample| {
1781 let uv = select(&sample.station);
1782 Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1783 })
1784 .collect();
1785 let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1786 Ok(curve)
1790 };
1791 pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1792 pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1793 }
1794
1795 let wall_surface = unmet.as_ref().map(|_| surface.clone());
1796 let built = chain_surgery(
1797 solid,
1798 segments,
1799 ChainRows {
1800 surface,
1801 cr,
1802 cs,
1803 u_domain,
1804 fit_low: low,
1805 fit_range: range,
1806 pcurves1,
1807 pcurves2,
1808 whole_pcurves,
1809 crease,
1810 },
1811 name,
1812 )?;
1813 if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1814 crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1815 name: name.map(str::to_string),
1816 surface,
1817 request: model,
1818 residual,
1819 detail: request_detail,
1820 budget: crate::ApproximationBudget {
1821 reason,
1822 rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1825 rounds_limit: model_rounds_limit(),
1826 stations: stations_now,
1827 station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1828 mechanism: crate::BudgetMechanism::LocalRounds,
1829 measured_component: request_component,
1830 unread: request_unread,
1831 },
1832 });
1833 }
1834 Ok(Rung::Built(built))
1835}
1836
1837pub(super) struct ChainRows {
1838 pub(super) surface: NurbsSurface,
1839 pub(super) cr: NurbsCurve,
1840 pub(super) cs: NurbsCurve,
1841 pub(super) u_domain: [f64; 2],
1842 pub(super) fit_low: f64,
1846 pub(super) fit_range: f64,
1847 pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1850 pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1851 pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1853 pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1856}
1857
1858impl ChainRows {
1859 pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1860 self.fit_low + fit_parameter * self.fit_range
1861 }
1862}
1863
1864pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1868 let (low, high, curve) = fitted;
1869 let to_local = |value: f64| {
1872 let mut best = value;
1873 for candidate in [value, value + 1.0, value - 1.0] {
1874 if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1875 best = candidate;
1876 break;
1877 }
1878 }
1879 ((best - low) / (high - low)).clamp(0.0, 1.0)
1880 };
1881 let mut local_a = to_local(a);
1882 let mut local_b = to_local(b);
1883 if local_a > local_b {
1884 std::mem::swap(&mut local_a, &mut local_b);
1885 }
1886 let epsilon = 1e-9;
1887 let (_, tail) = if local_a > epsilon {
1888 curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1889 } else {
1890 (curve.clone(), curve.clone())
1891 };
1892 let tail = if local_a > epsilon {
1893 tail
1894 } else {
1895 curve.clone()
1896 };
1897 let portion = if local_b < 1.0 - epsilon {
1898 tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1899 } else {
1900 tail
1901 };
1902 Ok(portion)
1903}
1904
1905pub(super) fn cross_edge_at(
1911 solid: &BrepSolid,
1912 face_before: &FaceRecord,
1913 loop_before: usize,
1914 face_after: &FaceRecord,
1915 loop_after: usize,
1916 vertex: u64,
1917 before_edge: u64,
1918 after_edge: u64,
1919) -> Result<Option<u64>, KernelRefusal> {
1920 let mut found: Option<u64> = None;
1921 let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1922 for coedge in &face.loops[loop_index].coedges {
1923 if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1924 continue;
1925 }
1926 let edge = solid
1927 .edges
1928 .iter()
1929 .find(|edge| edge.id == coedge.edge_id)
1930 .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1931 if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1932 if found.is_some() && found != Some(edge.id) {
1933 return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1934 }
1935 found = Some(edge.id);
1936 }
1937 }
1938 Ok(())
1939 };
1940 consider(face_before, loop_before)?;
1941 if !(face_after.id == face_before.id && loop_after == loop_before) {
1942 consider(face_after, loop_after)?;
1943 }
1944 Ok(found)
1945}
1946
1947fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1950 let period = high - low;
1951 low + ((value - low) / period).round() * period
1952}
1953
1954pub(in crate::blend) fn project_piece_pcurve(
1972 piece: &NurbsCurve,
1973 surface: &NurbsSurface,
1974 start_is_crossing: bool,
1975 end_is_crossing: bool,
1976) -> Result<NurbsCurve, KernelRefusal> {
1977 project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1978}
1979
1980
1981pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
1986 let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
1987 let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1988 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1989 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1990 let into = move |value: f64, closed: bool, low: f64, high: f64| {
1991 if closed && high > low {
1992 low + (value - low).rem_euclid(high - low)
1993 } else {
1994 value.clamp(low, high)
1995 }
1996 };
1997 Ok(move |fraction: f64| {
1998 let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
1999 let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
2000 seed.iter().all(|value| value.is_finite()).then_some(seed)
2001 })
2002}
2003
2004pub(in crate::blend) fn project_piece_pcurve_from(
2014 piece: &NurbsCurve,
2015 surface: &NurbsSurface,
2016 start_is_crossing: bool,
2017 end_is_crossing: bool,
2018 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2019) -> Result<NurbsCurve, KernelRefusal> {
2020 project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
2021}
2022
2023pub(crate) const CHAIN_PCURVE_FLOOR_WHAT: &str = "chain_pcurve_floor";
2029
2030pub(in crate::blend) fn project_piece_pcurve_within(
2031 piece: &NurbsCurve,
2032 surface: &NurbsSurface,
2033 start_is_crossing: bool,
2034 end_is_crossing: bool,
2035 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2036 refine_near: &[f64],
2037) -> Result<NurbsCurve, KernelRefusal> {
2038 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false)
2039}
2040
2041pub(in crate::blend) fn project_piece_pcurve_within_knot_retry(
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 match project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false) {
2060 Err(uniform) if matches!(&uniform.class, crate::RefusalClass::NonConvergence { what } if what == CHAIN_PCURVE_FLOOR_WHAT) => {
2061 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, true).map_err(|_| uniform)
2062 }
2063 other => other,
2064 }
2065}
2066
2067#[allow(clippy::too_many_arguments)]
2071fn project_piece_pcurve_track(
2072 piece: &NurbsCurve,
2073 surface: &NurbsSurface,
2074 start_is_crossing: bool,
2075 end_is_crossing: bool,
2076 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2077 refine_near: &[f64],
2078 knot_aware: bool,
2079) -> Result<NurbsCurve, KernelRefusal> {
2080 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2081 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2082 let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2083 let period = u1 - u0;
2084 let period_v = v1 - v0;
2085 let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
2086 let ceiling = if knot_aware { 2 * PIECE_TRACK_CEILING } else { PIECE_TRACK_CEILING };
2096 let uniform = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
2097 let mut fractions: Vec<f64> = (0..=uniform).map(|section| section as f64 / uniform as f64).collect();
2098 if knot_aware {
2099 let span = d1 - d0;
2102 let mut knots: Vec<f64> = piece.knots.iter().copied().filter(|knot| *knot > d0 && *knot < d1).map(|knot| (knot - d0) / span).collect();
2103 knots.dedup();
2104 fractions.extend(knots);
2105 fractions.sort_by(|a, b| a.total_cmp(b));
2106 fractions.dedup_by(|b, a| (*b - *a).abs() <= 1e-6 / uniform as f64);
2107 }
2108 let dense = fractions.len() - 1;
2109 crate::blend::record_count("count.track.piece_samples", dense as u64 + 1);
2110 let mut proj_u = Vec::with_capacity(dense + 1);
2111 let mut proj_v = Vec::with_capacity(dense + 1);
2112 let mut point3 = Vec::with_capacity(dense + 1);
2113 let mut params = Vec::with_capacity(dense + 1);
2114 let mut previous_u: Option<f64> = None;
2115 let mut previous_foot: Option<[f64; 2]> = None;
2120 let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
2121 let project_timer = crate::blend::BlendTimer::start("timing.track.project_ms");
2122 for section in 0..=dense {
2123 let fraction = fractions[section];
2124 let t = if knot_aware { d0 + (d1 - d0) * fraction } else { d0 + (d1 - d0) * section as f64 / dense as f64 };
2128 let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
2129 let guide = guide.filter(|_| {
2130 true
2131 });
2132 let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(fraction)).into_iter().chain(previous_foot).collect();
2137 let mut projection: Option<crate::projection::SurfaceProjection> = None;
2138 for seed in &seeds {
2139 let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
2140 if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
2141 projection = Some(found);
2142 }
2143 }
2144 let projection = match projection {
2145 Some(projection) => projection,
2146 None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
2147 };
2148 previous_foot = Some([projection.u, projection.v]);
2149 feet.push([projection.u, projection.v]);
2150 let mut u = projection.u;
2151 if let Some(previous) = previous_u {
2152 while u - previous > period * 0.5 {
2153 u -= period;
2154 }
2155 while previous - u > period * 0.5 {
2156 u += period;
2157 }
2158 }
2159 previous_u = Some(u);
2160 proj_u.push(u);
2161 proj_v.push(projection.v);
2162 point3.push(point);
2163 params.push(fraction);
2164 }
2165 project_timer.stop();
2166 if closed_v {
2173 for index in 2..=dense {
2174 while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
2175 proj_v[index] -= period_v;
2176 }
2177 while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
2178 proj_v[index] += period_v;
2179 }
2180 }
2181 while proj_v[0] - proj_v[1] > period_v * 0.5 {
2182 proj_v[0] -= period_v;
2183 }
2184 while proj_v[1] - proj_v[0] > period_v * 0.5 {
2185 proj_v[0] += period_v;
2186 }
2187 }
2188 let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2191 let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
2192 let pinned = |index: usize,
2201 neighbour: usize,
2202 proj_u: &[f64],
2203 proj_v: &[f64]|
2204 -> Result<(f64, f64), KernelRefusal> {
2205 let point = point3[index];
2206 let u_error = surface
2207 .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
2208 .sub(point)
2209 .length();
2210 if closed_v {
2211 let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
2212 let v_error = surface
2213 .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
2214 .sub(point)
2215 .length();
2216 if v_error < u_error {
2217 return Ok((proj_u[index], seam_v));
2218 }
2219 }
2220 Ok((seam_u, proj_v[index]))
2221 };
2222 if start_is_crossing {
2223 let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
2224 proj_u[0] = u;
2225 proj_v[0] = v;
2226 }
2227 if end_is_crossing {
2228 let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
2229 proj_u[dense] = u;
2230 proj_v[dense] = v;
2231 }
2232 let mut shift = 0.0;
2233 while mean_u + shift > u1 + 1e-9 {
2234 shift -= period;
2235 }
2236 while mean_u + shift < u0 - 1e-9 {
2237 shift += period;
2238 }
2239 for u in proj_u.iter_mut() {
2240 *u += shift;
2241 }
2242 if closed_v {
2243 let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2244 let mut shift_v = 0.0;
2245 while mean_v + shift_v > v1 + 1e-9 {
2246 shift_v -= period_v;
2247 }
2248 while mean_v + shift_v < v0 - 1e-9 {
2249 shift_v += period_v;
2250 }
2251 for v in proj_v.iter_mut() {
2252 *v += shift_v;
2253 }
2254 }
2255 let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
2270 let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
2277 let mut track = crate::blend::track_fit::Track {
2278 fractions: params,
2279 uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
2280 feet,
2281 breaks: Vec::new(),
2282 };
2283 let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
2284 let polish_timer = crate::blend::BlendTimer::start("timing.track.polish_ms");
2285 crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
2286 polish_timer.stop();
2287 let crossings_timer = crate::blend::BlendTimer::start("timing.track.crossings_ms");
2288 crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
2289 crossings_timer.stop();
2290 crate::blend::record_count("count.track.breaks", track.breaks.len() as u64);
2291 for &fraction in refine_near {
2292 for _ in 0..2 {
2293 let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
2294 if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], ceiling)?.is_none() {
2295 break;
2296 }
2297 }
2298 }
2299 let fit_timer = crate::blend::BlendTimer::start("timing.track.fit_local_ms");
2306 let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
2307 surface,
2308 &mut track,
2309 &at,
2310 tolerance,
2311 REFINE_ROUNDS,
2312 ceiling,
2313 )?;
2314 fit_timer.stop();
2315 crate::blend::record_count("count.track.fit_sections", fit.sections as u64);
2316 crate::blend::record_count("count.track.local_inserted", inserted as u64);
2317 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2318 eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2319 eprintln!(
2320 "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2321 fit.miss,
2322 fit.samples,
2323 track.fractions.len(),
2324 track.breaks.len(),
2325 fit.on_floor
2326 );
2327 if let Some(worst) = fit.worst {
2332 let parameter = d0 + (d1 - d0) * worst.fraction;
2333 let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2334 .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2335 .map(|pair| [pair[0], pair[1]]).collect();
2336 let mean = (d1 - d0) / spans.len().max(1) as f64;
2337 let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2338 let width = span.map(|index| spans[index][1] - spans[index][0]);
2339 eprintln!(
2340 "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2341 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2342 worst.fraction, parameter, span, spans.len(), width, mean,
2343 worst.bracket[0], worst.bracket[1], worst.interval
2344 );
2345 }
2346 }
2347 if !fit.on_floor {
2354 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2355 for (label, t) in [("start", d0), ("end", d1)] {
2358 if let Ok(point) = piece.evaluate(t) {
2359 let standoff = crate::project_point_to_surface(surface, point).map(|p| p.distance).unwrap_or(f64::NAN);
2360 eprintln!(
2361 "TRACK_FIT closed-chain piece {label} ({:.9}, {:.9}, {:.9}) stands {standoff:.3e} off the carrier",
2362 point.x, point.y, point.z
2363 );
2364 }
2365 }
2366 }
2367 return Err(KernelRefusal::non_convergence(KernelStage::Refine, CHAIN_PCURVE_FLOOR_WHAT, format!(
2368 "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2369 against a floor of {tolerance:.1e}",
2370 crate::blend::PCURVE_OFF_FLOOR,
2371 fit.miss,
2372 fit.samples
2373 )));
2374 }
2375 Ok(fit.curve)
2376}
2377
2378
2379
2380pub(in crate::blend) struct PreparedWall<'a> {
2383 fine: &'a [Vec3],
2384 coarse: Vec<Vec3>,
2385 fine_stencils: StencilPreparation,
2386 coarse_stencils: StencilPreparation,
2387}
2388
2389struct StencilPreparation {
2390 finite: bool,
2391 coefficients: std::cell::RefCell<Vec<Option<[Vec3; 4]>>>,
2392}
2393impl StencilPreparation {
2394 fn new(nodes: &[Vec3]) -> Self {
2395 Self { finite: wall_nodes_finite(nodes), coefficients: std::cell::RefCell::new(vec![None; nodes.len().saturating_sub(3)]) }
2396 }
2397 fn coefficients(&self, nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2398 let mut entries = self.coefficients.borrow_mut();
2399 *entries[base].get_or_insert_with(|| wall_coefficients(nodes, base))
2400 }
2401}
2402
2403fn wall_nodes_finite(nodes: &[Vec3]) -> bool {
2404 nodes.iter().all(|c| c.x.is_finite() && c.y.is_finite() && c.z.is_finite())
2405}
2406
2407fn wall_coefficients(nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2408 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
2410 let d1 = p1.sub(p0);
2411 let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
2412 let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
2413 [p0, d1, d2, d3]
2414}
2415
2416fn wall_coarse(fine: &[Vec3]) -> Vec<Vec3> {
2417 fine.iter().step_by(2).copied().collect()
2418}
2419
2420fn wall_preparation_on() -> bool {
2421 std::env::var("BREP_WALL_PREPARE").ok().as_deref() != Some("0")
2422}
2423
2424impl<'a> PreparedWall<'a> {
2425 pub(in crate::blend) fn for_interval(fine: &'a [Vec3]) -> Option<Self> {
2426 (wall_preparation_on() && wall_solve_reuse_on()).then(|| Self::new(fine))
2427 }
2428 fn new(fine: &'a [Vec3]) -> Self {
2429 let coarse = wall_coarse(fine);
2430 let fine_stencils = StencilPreparation::new(fine);
2431 let coarse_stencils = StencilPreparation::new(&coarse);
2432 Self { fine, coarse, fine_stencils, coarse_stencils }
2433 }
2434 pub(in crate::blend) fn read(&self, p: Vec3, j: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
2435 let read = |at| declared_wall_probe_prepared(p, self.fine, at, radius, true, Some(self));
2436 match read(j) {
2437 Some(read) => Some(read),
2438 None => match centre_path_node(self.fine, fraction) {
2439 Some(located) if located != j => read(located),
2440 _ => None,
2441 },
2442 }
2443 }
2444}