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 let n = nodes.len();
672 if n < 4 || j == 0 || j + 1 >= n {
673 return None;
674 }
675 let mut tried: Vec<usize> = Vec::new();
676 let mut at = j;
677 loop {
678 let base = (at - 1).min(n - 4);
679 match stencil_solve(p, nodes, base, at, radius) {
680 Solve::Read(_) => return Some(at),
681 Solve::Unread => return None,
682 Solve::Window(s) => {
683 tried.push(at);
684 let next = (base as f64 + s).round();
685 let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
686 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
687 wall_decline(|| format!(
688 "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
689 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
690 ));
691 return None;
692 }
693 at = next;
694 }
695 }
696 }
697}
698
699fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
702 match stencil_solve(p, nodes, base, j, radius) {
703 Solve::Read(distance) => Some(distance),
704 Solve::Window(_) | Solve::Unread => None,
705 }
706}
707
708enum Solve {
711 Read(f64),
712 Window(f64),
713 Unread,
714}
715
716thread_local! {
717 static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
722}
723
724fn wall_decline(reason: impl FnOnce() -> String) {
725 WALL_DECLINES.with(|log| {
726 if let Some(lines) = log.borrow_mut().as_mut() {
727 lines.push(reason());
728 }
729 });
730}
731
732pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
740 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
741 let read = declared_wall_probe(p, fine, jf, radius);
742 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
743 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
744}
745
746fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
747 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
748 let n = nodes.len();
749 if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
750 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
751 return Solve::Unread;
752 }
753 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
755 let d1 = p1.sub(p0);
756 let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
757 let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
758 let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
759 let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
760 let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
761 let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
762 (c, c1, c2)
763 };
764 let start = j as f64 - base as f64;
765 let mut s = start;
766 let mut converged = false;
767 for _ in 0..50 {
768 let (c, c1, c2) = curve(s);
769 let w = c.sub(p);
770 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
771 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
772 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
773 return Solve::Unread;
774 }
775 let stationary = g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
788 if stationary {
789 converged = true;
790 break;
791 }
792 let step = g / slope;
793 let before = s;
794 s -= step;
795 if !((s - start).abs() <= 1.0) {
796 wall_decline(|| {
797 let (c, c1, _) = curve(s);
800 let g_after = c.sub(p).dot(c1);
801 let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
802 format!(
803 "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
804 g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
805 (s - start).abs() - 1.0
806 )
807 });
808 return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
809 }
810 if step.abs() <= 1e-14 * (1.0 + s.abs()) {
811 converged = true;
812 break;
813 }
814 }
815 if !converged {
816 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
817 return Solve::Unread;
818 }
819 let distance = curve(s).0.sub(p).length() - radius;
820 if distance.is_finite() {
821 Solve::Read(distance)
822 } else {
823 wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
824 Solve::Unread
825 }
826}
827
828pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
837 let n = nodes.len();
838 if n < 5 || j == 0 || j + 1 >= n {
839 return None;
840 }
841 let j = locate_foot(p, nodes, j, radius)?;
842 let first = (j - 1).min(n - 4);
843 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
844 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
845 let a = sweep_residual_on(p, nodes, first, j, radius)?;
846 let b = sweep_residual_on(p, nodes, second, j, radius)?;
847 Some((a.abs().max(b.abs()), (a - b).abs()))
848}
849
850pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
862 if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
863 wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
864 return None;
865 }
866 let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
867 wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
868 return None;
869 };
870 let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
871 wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
872 return None;
873 };
874 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
875 let Some(at_coarse) = sweep_residual(p, &coarse, jf / 2, radius) else {
876 wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
877 return None;
878 };
879 let uncertainty = spread.max((at_fine - at_coarse).abs());
880 uncertainty.is_finite().then_some((reading, uncertainty))
881}
882
883pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
899 let n = fine.len();
900 if n < 7 || !(fraction.is_finite()) {
901 return None;
902 }
903 let mut cumulative = Vec::with_capacity(n);
904 cumulative.push(0.0_f64);
905 for pair in fine.windows(2) {
906 let step = pair[1].sub(pair[0]).length();
907 cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
908 }
909 let total = *cumulative.last()?;
910 if !(total.is_finite() && total > 0.0) {
911 return None;
912 }
913 let target = fraction.clamp(0.0, 1.0) * total;
914 let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
915 let even = 2 * ((nearest as f64 / 2.0).round() as usize);
916 Some(even.clamp(2, (n - 3) & !1))
917}
918
919pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
926 let first = match declared_wall_probe_why(p, fine, jf, radius) {
927 Ok(read) => return Ok(read),
928 Err(why) => why,
929 };
930 match centre_path_node(fine, fraction) {
931 Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
932 .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
933 _ => Err(first),
934 }
935}
936
937pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
939 match declared_wall_probe(p, fine, jf, radius) {
940 Some(read) => Some(read),
941 None => match centre_path_node(fine, fraction) {
942 Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
943 _ => None,
944 },
945 }
946}
947
948pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
960 if jf % 2 != 0
961 || fine.len() < 9
962 || fine.len() % 2 == 0
963 || ts.len() != fine.len()
964 || !ts.iter().all(|t| t.is_finite())
965 || !ts.windows(2).all(|pair| pair[1] > pair[0])
966 {
967 return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
968 }
969 let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
970 let read = (|| -> Option<(f64, f64)> {
971 let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
972 wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
973 None
974 })?;
975 let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
976 wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
977 None
978 })?;
979 let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
980 let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
981 let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
982 wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
983 None
984 })?;
985 let uncertainty = spread.max((at_fine - at_coarse).abs());
986 uncertainty.is_finite().then_some((reading, uncertainty))
987 })();
988 let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
989 read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
990}
991
992fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
993 let n = nodes.len();
994 if n < 5 || j == 0 || j + 1 >= n {
995 return None;
996 }
997 let j = locate_foot_at(p, nodes, ts, j, radius)?;
998 let first = (j - 1).min(n - 4);
999 let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
1000 let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
1001 let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
1002 let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
1003 Some((a.abs().max(b.abs()), (a - b).abs()))
1004}
1005
1006fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
1007 let j = locate_foot_at(p, nodes, ts, j, radius)?;
1008 match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
1009 Solve::Read(distance) => Some(distance),
1010 Solve::Window(_) | Solve::Unread => None,
1011 }
1012}
1013
1014fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
1018 let n = nodes.len();
1019 if n < 4 || j == 0 || j + 1 >= n {
1020 return None;
1021 }
1022 let mut tried: Vec<usize> = Vec::new();
1023 let mut at = j;
1024 loop {
1025 let base = (at - 1).min(n - 4);
1026 match stencil_solve_at(p, nodes, ts, base, at, radius) {
1027 Solve::Read(_) => return Some(at),
1028 Solve::Unread => return None,
1029 Solve::Window(t) => {
1030 tried.push(at);
1031 let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
1032 let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
1033 if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
1034 wall_decline(|| format!(
1035 "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
1036 if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
1037 ));
1038 return None;
1039 }
1040 at = next;
1041 }
1042 }
1043 }
1044}
1045
1046fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
1051 let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
1052 let n = nodes.len();
1053 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)) {
1054 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
1055 return Solve::Unread;
1056 }
1057 let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
1058 let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
1059 let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
1060 let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
1061 let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
1062 let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
1063 let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
1064 let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
1065 let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
1066 let (a, b, c) = (t - t0, t - t1, t - t2);
1067 let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
1068 let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
1069 let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
1070 (point, first, second)
1071 };
1072 let (low, high) = (ts[j - 1], ts[j + 1]);
1073 let start = ts[j];
1074 let mut t = start;
1075 let mut converged = false;
1076 for _ in 0..50 {
1077 let (c, c1, c2) = curve(t);
1078 let w = c.sub(p);
1079 let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
1080 if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
1081 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
1082 return Solve::Unread;
1083 }
1084 if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
1085 converged = true;
1086 break;
1087 }
1088 let step = g / slope;
1089 let before = t;
1090 t -= step;
1091 if !(t >= low && t <= high) {
1092 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
1093 return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
1094 }
1095 let node = 0.5 * (high - low);
1100 if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
1101 converged = true;
1102 break;
1103 }
1104 }
1105 if !converged {
1106 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
1107 return Solve::Unread;
1108 }
1109 let distance = curve(t).0.sub(p).length() - radius;
1110 if distance.is_finite() {
1111 Solve::Read(distance)
1112 } else {
1113 wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
1114 Solve::Unread
1115 }
1116}
1117
1118pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
1123 readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
1124}
1125
1126#[allow(clippy::too_many_arguments)]
1136pub(super) fn wall_interior_misses(
1137 contacts: &crate::blend::chain::DenseContacts,
1138 samples: &[ChainSample],
1139 surface: &crate::NurbsSurface,
1140 low: f64,
1141 range: f64,
1142 radius: f64,
1143 subdivisions: usize,
1144 probes: usize,
1145 wrap: bool,
1146) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1147 use std::collections::BTreeMap;
1148 let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1151 for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1152 let entry = nodes.entry(*left).or_insert_with(|| {
1153 let mut row = vec![None; subdivisions + 1];
1154 row[0] = Some(samples[*left].station.center);
1155 row[subdivisions] = Some(samples[*next].station.center);
1156 (*next, row)
1157 });
1158 let step = (fraction * subdivisions as f64).round() as usize;
1159 if (1..subdivisions).contains(&step) {
1160 entry.1[step] = contacts.centers.get(index).copied().flatten();
1161 }
1162 }
1163 let mut readable = Vec::new();
1164 let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1165 for (&left, (next, row)) in &nodes {
1166 let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1167 unreadable.push(left);
1168 continue;
1169 };
1170 let (from, to) = (samples[left].parameter, samples[*next].parameter);
1171 let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1172 let mut read = true;
1173 'probes: for step in (2..subdivisions).step_by(2) {
1174 let global = from + (step as f64 / subdivisions as f64) * (to - from);
1175 let global = if wrap { global.rem_euclid(1.0) } else { global };
1176 let u = (global - low) / range;
1177 for k in 1..probes {
1178 let v = k as f64 / probes as f64;
1179 let reading = surface.evaluate(u, v).ok().and_then(|p| {
1180 declared_wall_probe_located(p, ¢res, step, step as f64 / subdivisions as f64, radius)
1181 });
1182 match reading {
1183 Some((distance, gap)) => {
1184 declared = declared.max(distance + gap);
1185 worst = worst.max(distance);
1186 spread = spread.max(gap);
1187 }
1188 None => {
1189 read = false;
1190 break 'probes;
1191 }
1192 }
1193 }
1194 }
1195 if read {
1196 readable.push((left, declared, worst, spread));
1197 } else {
1198 unreadable.push(left);
1199 }
1200 }
1201 unreadable.sort_unstable();
1202 unreadable.dedup();
1203 (readable, unreadable)
1204}
1205
1206
1207fn build_closed_smooth_chain(
1208 solid: &BrepSolid,
1209 segments: &[ChainSegment<'_>],
1210 radius: f64,
1211 chamfer: bool,
1212 name: Option<&str>,
1213 budget: &ChainBudget,
1214 fold: FoldPolicy,
1215 bar: f64,
1216) -> Result<Rung, KernelRefusal> {
1217 let samples = match march_chain_typed(segments, radius, false, budget, fold) {
1218 Ok(samples) => samples,
1219 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1220 Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1223 Err(ChainMarchError::Other(error)) => return Err(error),
1224 };
1225 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1226 build_closed_chain_from_samples(solid, segments, radius, chamfer, name, budget, fold, bar, samples, 0)
1227}
1228
1229
1230thread_local! {
1231 static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1235}
1236
1237fn local_continue_on() -> bool {
1243 std::env::var("BREP_CHAIN_LOCAL_CONTINUE").as_deref() == Ok("1")
1244}
1245
1246const LOCAL_CONTINUE_MAX_ROUNDS: usize = 4 * REFINE_ROUNDS;
1248
1249const LOCAL_CONTINUE_STALL_ROUNDS: usize = 3;
1251
1252thread_local! {
1253 static CHAIN_LOCAL_BEST: std::cell::Cell<(f64, usize)> = const { std::cell::Cell::new((f64::INFINITY, 0)) };
1255}
1256
1257fn local_rounds_spent(round: usize, deviation: f64) -> bool {
1263 if round == 0 {
1264 CHAIN_LOCAL_BEST.with(|best| best.set((f64::INFINITY, 0)));
1265 return false;
1266 }
1267 let (_, best_round) = CHAIN_LOCAL_BEST.with(|best| {
1268 let (reading, at) = best.get();
1269 if deviation < reading {
1270 best.set((deviation, round));
1271 (deviation, round)
1272 } else {
1273 (reading, at)
1274 }
1275 });
1276 round >= LOCAL_CONTINUE_MAX_ROUNDS || round - best_round >= LOCAL_CONTINUE_STALL_ROUNDS
1277}
1278#[allow(clippy::too_many_arguments)]
1289fn build_closed_chain_from_samples(
1290 solid: &BrepSolid,
1291 segments: &[ChainSegment<'_>],
1292 radius: f64,
1293 chamfer: bool,
1294 name: Option<&str>,
1295 budget: &ChainBudget,
1296 fold: FoldPolicy,
1297 bar: f64,
1298 mut samples: Vec<ChainSample>,
1299 round: usize,
1300) -> Result<Rung, KernelRefusal> {
1301 let per_segment = budget.rung;
1304 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1305
1306 let mut global: Vec<(f64, &ChainSample)> = samples
1308 .iter()
1309 .filter(|sample| (0..budget.count(sample.segment) as isize).contains(&sample.position))
1310 .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1311 .collect();
1312 global.sort_by(|a, b| a.0.total_cmp(&b.0));
1313 let count = global.len();
1314 let degree = FIT_DEGREE;
1315 let mut extended_params = Vec::new();
1316 let mut samples_cr = Vec::new();
1317 let mut samples_mid = Vec::new();
1318 let mut samples_cs = Vec::new();
1319 let mut push = |station: &Station, parameter: f64| {
1320 extended_params.push(parameter);
1321 samples_cr.push(Vec4::from_point(station.p1, 1.0));
1322 samples_cs.push(Vec4::from_point(station.p2, 1.0));
1323 samples_mid.push(Vec4 {
1324 x: station.apex.x * station.weight,
1325 y: station.apex.y * station.weight,
1326 z: station.apex.z * station.weight,
1327 w: station.weight,
1328 });
1329 };
1330 for offset in (1..=degree).rev() {
1331 let (parameter, sample) = &global[count - offset];
1332 push(&sample.station, parameter - 1.0);
1333 }
1334 for (parameter, sample) in &global {
1335 push(&sample.station, *parameter);
1336 }
1337 push(&global[0].1.station, global[0].0 + 1.0);
1340 for offset in 1..=degree {
1341 let (parameter, sample) = &global[offset];
1342 push(&sample.station, parameter + 1.0);
1343 }
1344 let low = extended_params[0];
1345 let high = *extended_params.last().unwrap();
1346 let range = high - low;
1347 let normalized: Vec<f64> = extended_params
1348 .iter()
1349 .map(|parameter| (parameter - low) / range)
1350 .collect();
1351 let seam_low = (0.0 - low) / range;
1352 let seam_high = (1.0 - low) / range;
1353 let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1354 let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1355 let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1356 let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1357 Ok(middle)
1358 };
1359 let cr = fit_row(&samples_cr)?;
1360 let cs = fit_row(&samples_cs)?;
1361 let mid = if chamfer {
1362 None
1363 } else {
1364 Some(fit_row(&samples_mid)?)
1365 };
1366 let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1367 let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1368
1369 let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, budget, bar)?;
1382 let stencil = if round > 0 {
1394 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, LOCAL_RAIL_STENCIL)?)
1395 } else {
1396 None
1397 };
1398 let wall_stencil = if !chamfer && wall_verdict_on() {
1403 Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1404 } else {
1405 None
1406 };
1407 let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1408 let contacts = stencil.as_ref().expect("the local path solves the stencil");
1409 let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1410 let reads = dense.len() + unreadable.len() + contacts.unsolved;
1411 let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1412 let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1413 crate::blend::carve::carve_trace(format_args!(
1414 "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1415 {} failing interval(s), {} unsolved contact(s)",
1416 verdict.0, verdict.1.len(), contacts.unsolved
1417 ));
1418 (verdict.0, verdict.1, reads)
1419 } else {
1420 (halfway, halfway_failing, 0)
1421 };
1422 let mut _wall_worst = 0.0_f64;
1428 let mut _wall_spread = 0.0_f64;
1429 let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1430 if let Some(contacts) = wall_stencil.as_ref() {
1431 let (wall, unreadable) =
1432 wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1433 wall_readings = wall.clone();
1434 for &(left, declared, _, uncertainty) in &wall {
1436 if declared > bar {
1437 failing.push(left);
1438 }
1439 _wall_worst = _wall_worst.max(declared);
1440 _wall_spread = _wall_spread.max(uncertainty);
1441 }
1442 failing.extend_from_slice(&unreadable);
1443 failing.sort_unstable();
1444 failing.dedup();
1445 deviation = deviation.max(_wall_worst);
1446 crate::blend::carve::carve_trace(format_args!(
1447 "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1448 interval(s), {} unreadable, bar {bar:.3e}",
1449 wall.len(),
1450 unreadable.len()
1451 ));
1452 }
1453 if deviation > bar || !failing.is_empty() {
1454 let in_segment = |segment: usize| {
1455 samples
1456 .iter()
1457 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1458 .count()
1459 };
1460 let within_ceiling = (0..segments.len()).all(|segment| {
1461 in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1462 <= CHAIN_CARVE_MAX_PER_SEGMENT
1463 });
1464 crate::blend::carve::carve_trace(format_args!(
1465 "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1466 worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1467 samples.len(),
1468 failing.len(),
1469 if within_ceiling { "ok" } else { "reached" }
1470 ));
1471 let rounds_spent = if local_continue_on() {
1473 local_rounds_spent(round, deviation)
1474 } else {
1475 round >= REFINE_ROUNDS
1476 };
1477 if fold != FoldPolicy::Refuse || rounds_spent || !within_ceiling {
1479 return Ok(Rung::TooCoarse { deviation });
1480 }
1481 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1483 match insert_chain_stations(segments, radius, budget, &mut samples, &failing) {
1484 Ok(()) => {}
1485 Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1486 Err(ChainMarchError::Incoherent(failure)) => {
1490 crate::blend::carve::carve_trace(format_args!(
1491 "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1492 at t {:.9} does not reach its retained right neighbour: {}",
1493 failure.segment, failure.t, failure.error
1494 ));
1495 return Ok(Rung::TooCoarse { deviation });
1496 }
1497 Err(ChainMarchError::Other(error)) => return Err(error),
1498 }
1499 assign_chain_parameters(&mut samples, segments.len(), budget, false);
1500 return build_closed_chain_from_samples(
1501 solid, segments, radius, chamfer, name, budget, fold, bar, samples, round + 1,
1502 );
1503 }
1504 crate::blend::carve::carve_trace(format_args!(
1505 "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1506 from the rolling ball, inside {bar:.3e}"
1507 ));
1508 let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1522 let stations_now = (0..segments.len())
1523 .map(|segment| {
1524 samples
1525 .iter()
1526 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1527 .count()
1528 })
1529 .max()
1530 .unwrap_or(0);
1531 let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1532 let mut request_detail: Option<String> = None;
1536 let mut request_component = crate::MeasuredComponent::Wall;
1537 let mut request_unread = 0usize;
1538 if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1539 let wall_missing = wall_model_missing(&wall_readings, model);
1540 let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1548 Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1549 None => (Vec::new(), Vec::new()),
1550 };
1551 let rail_over: std::collections::BTreeSet<usize> =
1553 rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1554 let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1555 let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1556 let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1557 missing.sort_unstable();
1558 missing.dedup();
1559 if !missing.is_empty() {
1560 let wall_short = wall_readings
1562 .iter()
1563 .filter(|reading| wall_missing.contains(&reading.0))
1564 .map(|&(_, declared, _, _)| declared)
1565 .fold(0.0_f64, f64::max);
1566 let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1567 let short = wall_short.max(rail_short);
1568 request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1569 request_unread = rail_unread_intervals.len();
1570 request_detail = Some(format!(
1571 "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1572 wall {wall_short:.3e} over {} interval(s)",
1573 rail_over.len(),
1574 rail_unread_intervals.len(),
1575 wall_missing.len()
1576 ));
1577 let within_ceiling = (0..segments.len()).all(|segment| {
1578 let present = samples
1579 .iter()
1580 .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1581 .count();
1582 present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1583 <= CHAIN_CARVE_MAX_PER_SEGMENT
1584 });
1585 crate::blend::carve::carve_trace(format_args!(
1586 "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1587 (worst reading + spread {short:.3e}), ceiling {}",
1588 missing.len(),
1589 if within_ceiling { "ok" } else { "reached" }
1590 ));
1591 unmet = if round >= model_rounds_limit() {
1592 Some((short, crate::BudgetReason::RoundsSpent))
1593 } else if !within_ceiling {
1594 Some((short, crate::BudgetReason::StationCeiling))
1595 } else {
1596 let mut refined = samples.clone();
1597 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1599 match insert_chain_stations(segments, radius, budget, &mut refined, &missing) {
1600 Ok(()) => {
1601 assign_chain_parameters(&mut refined, segments.len(), budget, false);
1602 match build_closed_chain_from_samples(
1603 solid, segments, radius, chamfer, name, budget, fold, bar, refined, round + 1,
1604 ) {
1605 Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1606 _ => Some((short, crate::BudgetReason::Rejected)),
1607 }
1608 }
1609 Err(_) => {
1610 CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1612 Some((short, crate::BudgetReason::Incoherent))
1613 }
1614 }
1615 };
1616 }
1617 }
1618 if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1623 let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, subdivisions.max(2))?;
1624 let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1625 let mut unreadable = 0usize;
1626 for (left, next, fraction, p1, p2) in &contacts.solved {
1627 let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1628 let global = (from + fraction * (to - from)).rem_euclid(1.0);
1629 let seed = (global - low) / range;
1630 let reach = (to - from).abs() / range;
1631 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)) {
1632 Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1633 Some(_) => {}
1634 None => unreadable += 1,
1635 }
1636 }
1637 eprintln!(
1638 "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1639 worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1640 contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1641 );
1642 }
1643 let crease = if fold == FoldPolicy::Carve {
1644 let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1654 let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1655 return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1656 "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"
1657 )));
1658 };
1659 Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1660 } else {
1661 None
1662 };
1663
1664 let single_face = |side: usize| -> bool {
1668 let first_id = segments[0].mate(side).face.id;
1669 segments
1670 .iter()
1671 .all(|segment| segment.mate(side).face.id == first_id)
1672 };
1673 let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1674 for side in 0..2 {
1675 if !single_face(side) {
1676 continue;
1677 }
1678 let mut row = Vec::new();
1679 let select = |station: &Station| -> [f64; 2] {
1680 if side == 0 {
1681 station.uv1
1682 } else {
1683 station.uv2
1684 }
1685 };
1686 let mut push_uv = |station: &Station| {
1687 let uv = select(station);
1688 row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1689 };
1690 for offset in (1..=degree).rev() {
1691 push_uv(&global[count - offset].1.station);
1692 }
1693 for (_, sample) in &global {
1694 push_uv(&sample.station);
1695 }
1696 push_uv(&global[0].1.station);
1697 for offset in 1..=degree {
1698 push_uv(&global[offset].1.station);
1699 }
1700 whole_pcurves[side] = Some(fit_row(&row)?);
1701 }
1702
1703 let mut pcurves1 = Vec::with_capacity(segments.len());
1706 let mut pcurves2 = Vec::with_capacity(segments.len());
1707 for segment in 0..segments.len() {
1708 let mut window: Vec<&ChainSample> = samples
1709 .iter()
1710 .filter(|sample| sample.segment == segment)
1711 .collect();
1712 window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1713 let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1714 let low = params[0];
1715 let high = *params.last().unwrap();
1716 let normalized: Vec<f64> = params
1717 .iter()
1718 .map(|parameter| (parameter - low) / (high - low))
1719 .collect();
1720 let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1721 let points: Vec<Vec4> = window
1722 .iter()
1723 .map(|sample| {
1724 let uv = select(&sample.station);
1725 Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1726 })
1727 .collect();
1728 let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1729 Ok(curve)
1733 };
1734 pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1735 pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1736 }
1737
1738 let wall_surface = unmet.as_ref().map(|_| surface.clone());
1739 let built = chain_surgery(
1740 solid,
1741 segments,
1742 ChainRows {
1743 surface,
1744 cr,
1745 cs,
1746 u_domain,
1747 fit_low: low,
1748 fit_range: range,
1749 pcurves1,
1750 pcurves2,
1751 whole_pcurves,
1752 crease,
1753 },
1754 name,
1755 )?;
1756 if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1757 crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1758 name: name.map(str::to_string),
1759 surface,
1760 request: model,
1761 residual,
1762 detail: request_detail,
1763 budget: crate::ApproximationBudget {
1764 reason,
1765 rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1768 rounds_limit: model_rounds_limit(),
1769 stations: stations_now,
1770 station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1771 mechanism: crate::BudgetMechanism::LocalRounds,
1772 measured_component: request_component,
1773 unread: request_unread,
1774 },
1775 });
1776 }
1777 Ok(Rung::Built(built))
1778}
1779
1780pub(super) struct ChainRows {
1781 pub(super) surface: NurbsSurface,
1782 pub(super) cr: NurbsCurve,
1783 pub(super) cs: NurbsCurve,
1784 pub(super) u_domain: [f64; 2],
1785 pub(super) fit_low: f64,
1789 pub(super) fit_range: f64,
1790 pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1793 pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1794 pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1796 pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1799}
1800
1801impl ChainRows {
1802 pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1803 self.fit_low + fit_parameter * self.fit_range
1804 }
1805}
1806
1807pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1811 let (low, high, curve) = fitted;
1812 let to_local = |value: f64| {
1815 let mut best = value;
1816 for candidate in [value, value + 1.0, value - 1.0] {
1817 if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1818 best = candidate;
1819 break;
1820 }
1821 }
1822 ((best - low) / (high - low)).clamp(0.0, 1.0)
1823 };
1824 let mut local_a = to_local(a);
1825 let mut local_b = to_local(b);
1826 if local_a > local_b {
1827 std::mem::swap(&mut local_a, &mut local_b);
1828 }
1829 let epsilon = 1e-9;
1830 let (_, tail) = if local_a > epsilon {
1831 curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1832 } else {
1833 (curve.clone(), curve.clone())
1834 };
1835 let tail = if local_a > epsilon {
1836 tail
1837 } else {
1838 curve.clone()
1839 };
1840 let portion = if local_b < 1.0 - epsilon {
1841 tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1842 } else {
1843 tail
1844 };
1845 Ok(portion)
1846}
1847
1848pub(super) fn cross_edge_at(
1854 solid: &BrepSolid,
1855 face_before: &FaceRecord,
1856 loop_before: usize,
1857 face_after: &FaceRecord,
1858 loop_after: usize,
1859 vertex: u64,
1860 before_edge: u64,
1861 after_edge: u64,
1862) -> Result<Option<u64>, KernelRefusal> {
1863 let mut found: Option<u64> = None;
1864 let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1865 for coedge in &face.loops[loop_index].coedges {
1866 if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1867 continue;
1868 }
1869 let edge = solid
1870 .edges
1871 .iter()
1872 .find(|edge| edge.id == coedge.edge_id)
1873 .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1874 if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1875 if found.is_some() && found != Some(edge.id) {
1876 return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1877 }
1878 found = Some(edge.id);
1879 }
1880 }
1881 Ok(())
1882 };
1883 consider(face_before, loop_before)?;
1884 if !(face_after.id == face_before.id && loop_after == loop_before) {
1885 consider(face_after, loop_after)?;
1886 }
1887 Ok(found)
1888}
1889
1890fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1893 let period = high - low;
1894 low + ((value - low) / period).round() * period
1895}
1896
1897pub(in crate::blend) fn project_piece_pcurve(
1915 piece: &NurbsCurve,
1916 surface: &NurbsSurface,
1917 start_is_crossing: bool,
1918 end_is_crossing: bool,
1919) -> Result<NurbsCurve, KernelRefusal> {
1920 project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1921}
1922
1923
1924pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
1929 let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
1930 let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1931 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1932 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1933 let into = move |value: f64, closed: bool, low: f64, high: f64| {
1934 if closed && high > low {
1935 low + (value - low).rem_euclid(high - low)
1936 } else {
1937 value.clamp(low, high)
1938 }
1939 };
1940 Ok(move |fraction: f64| {
1941 let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
1942 let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
1943 seed.iter().all(|value| value.is_finite()).then_some(seed)
1944 })
1945}
1946
1947pub(in crate::blend) fn project_piece_pcurve_from(
1957 piece: &NurbsCurve,
1958 surface: &NurbsSurface,
1959 start_is_crossing: bool,
1960 end_is_crossing: bool,
1961 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1962) -> Result<NurbsCurve, KernelRefusal> {
1963 project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
1964}
1965
1966pub(crate) const CHAIN_PCURVE_FLOOR_WHAT: &str = "chain_pcurve_floor";
1972
1973pub(in crate::blend) fn project_piece_pcurve_within(
1974 piece: &NurbsCurve,
1975 surface: &NurbsSurface,
1976 start_is_crossing: bool,
1977 end_is_crossing: bool,
1978 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1979 refine_near: &[f64],
1980) -> Result<NurbsCurve, KernelRefusal> {
1981 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false)
1982}
1983
1984pub(in crate::blend) fn project_piece_pcurve_within_knot_retry(
1995 piece: &NurbsCurve,
1996 surface: &NurbsSurface,
1997 start_is_crossing: bool,
1998 end_is_crossing: bool,
1999 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2000 refine_near: &[f64],
2001) -> Result<NurbsCurve, KernelRefusal> {
2002 match project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false) {
2003 Err(uniform) if matches!(&uniform.class, crate::RefusalClass::NonConvergence { what } if what == CHAIN_PCURVE_FLOOR_WHAT) => {
2004 project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, true).map_err(|_| uniform)
2005 }
2006 other => other,
2007 }
2008}
2009
2010#[allow(clippy::too_many_arguments)]
2014fn project_piece_pcurve_track(
2015 piece: &NurbsCurve,
2016 surface: &NurbsSurface,
2017 start_is_crossing: bool,
2018 end_is_crossing: bool,
2019 guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2020 refine_near: &[f64],
2021 knot_aware: bool,
2022) -> Result<NurbsCurve, KernelRefusal> {
2023 let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2024 let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2025 let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2026 let period = u1 - u0;
2027 let period_v = v1 - v0;
2028 let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
2029 let ceiling = if knot_aware { 2 * PIECE_TRACK_CEILING } else { PIECE_TRACK_CEILING };
2039 let uniform = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
2040 let mut fractions: Vec<f64> = (0..=uniform).map(|section| section as f64 / uniform as f64).collect();
2041 if knot_aware {
2042 let span = d1 - d0;
2045 let mut knots: Vec<f64> = piece.knots.iter().copied().filter(|knot| *knot > d0 && *knot < d1).map(|knot| (knot - d0) / span).collect();
2046 knots.dedup();
2047 fractions.extend(knots);
2048 fractions.sort_by(|a, b| a.total_cmp(b));
2049 fractions.dedup_by(|b, a| (*b - *a).abs() <= 1e-6 / uniform as f64);
2050 }
2051 let dense = fractions.len() - 1;
2052 crate::blend::record_count("count.track.piece_samples", dense as u64 + 1);
2053 let mut proj_u = Vec::with_capacity(dense + 1);
2054 let mut proj_v = Vec::with_capacity(dense + 1);
2055 let mut point3 = Vec::with_capacity(dense + 1);
2056 let mut params = Vec::with_capacity(dense + 1);
2057 let mut previous_u: Option<f64> = None;
2058 let mut previous_foot: Option<[f64; 2]> = None;
2063 let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
2064 let project_timer = crate::blend::BlendTimer::start("timing.track.project_ms");
2065 for section in 0..=dense {
2066 let fraction = fractions[section];
2067 let t = if knot_aware { d0 + (d1 - d0) * fraction } else { d0 + (d1 - d0) * section as f64 / dense as f64 };
2071 let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
2072 let guide = guide.filter(|_| {
2073 true
2074 });
2075 let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(fraction)).into_iter().chain(previous_foot).collect();
2080 let mut projection: Option<crate::projection::SurfaceProjection> = None;
2081 for seed in &seeds {
2082 let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
2083 if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
2084 projection = Some(found);
2085 }
2086 }
2087 let projection = match projection {
2088 Some(projection) => projection,
2089 None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
2090 };
2091 previous_foot = Some([projection.u, projection.v]);
2092 feet.push([projection.u, projection.v]);
2093 let mut u = projection.u;
2094 if let Some(previous) = previous_u {
2095 while u - previous > period * 0.5 {
2096 u -= period;
2097 }
2098 while previous - u > period * 0.5 {
2099 u += period;
2100 }
2101 }
2102 previous_u = Some(u);
2103 proj_u.push(u);
2104 proj_v.push(projection.v);
2105 point3.push(point);
2106 params.push(fraction);
2107 }
2108 project_timer.stop();
2109 if closed_v {
2116 for index in 2..=dense {
2117 while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
2118 proj_v[index] -= period_v;
2119 }
2120 while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
2121 proj_v[index] += period_v;
2122 }
2123 }
2124 while proj_v[0] - proj_v[1] > period_v * 0.5 {
2125 proj_v[0] -= period_v;
2126 }
2127 while proj_v[1] - proj_v[0] > period_v * 0.5 {
2128 proj_v[0] += period_v;
2129 }
2130 }
2131 let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2134 let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
2135 let pinned = |index: usize,
2144 neighbour: usize,
2145 proj_u: &[f64],
2146 proj_v: &[f64]|
2147 -> Result<(f64, f64), KernelRefusal> {
2148 let point = point3[index];
2149 let u_error = surface
2150 .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
2151 .sub(point)
2152 .length();
2153 if closed_v {
2154 let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
2155 let v_error = surface
2156 .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
2157 .sub(point)
2158 .length();
2159 if v_error < u_error {
2160 return Ok((proj_u[index], seam_v));
2161 }
2162 }
2163 Ok((seam_u, proj_v[index]))
2164 };
2165 if start_is_crossing {
2166 let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
2167 proj_u[0] = u;
2168 proj_v[0] = v;
2169 }
2170 if end_is_crossing {
2171 let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
2172 proj_u[dense] = u;
2173 proj_v[dense] = v;
2174 }
2175 let mut shift = 0.0;
2176 while mean_u + shift > u1 + 1e-9 {
2177 shift -= period;
2178 }
2179 while mean_u + shift < u0 - 1e-9 {
2180 shift += period;
2181 }
2182 for u in proj_u.iter_mut() {
2183 *u += shift;
2184 }
2185 if closed_v {
2186 let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2187 let mut shift_v = 0.0;
2188 while mean_v + shift_v > v1 + 1e-9 {
2189 shift_v -= period_v;
2190 }
2191 while mean_v + shift_v < v0 - 1e-9 {
2192 shift_v += period_v;
2193 }
2194 for v in proj_v.iter_mut() {
2195 *v += shift_v;
2196 }
2197 }
2198 let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
2213 let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
2220 let mut track = crate::blend::track_fit::Track {
2221 fractions: params,
2222 uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
2223 feet,
2224 breaks: Vec::new(),
2225 };
2226 let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
2227 let polish_timer = crate::blend::BlendTimer::start("timing.track.polish_ms");
2228 crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
2229 polish_timer.stop();
2230 let crossings_timer = crate::blend::BlendTimer::start("timing.track.crossings_ms");
2231 crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
2232 crossings_timer.stop();
2233 crate::blend::record_count("count.track.breaks", track.breaks.len() as u64);
2234 for &fraction in refine_near {
2235 for _ in 0..2 {
2236 let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
2237 if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], ceiling)?.is_none() {
2238 break;
2239 }
2240 }
2241 }
2242 let fit_timer = crate::blend::BlendTimer::start("timing.track.fit_local_ms");
2249 let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
2250 surface,
2251 &mut track,
2252 &at,
2253 tolerance,
2254 REFINE_ROUNDS,
2255 ceiling,
2256 )?;
2257 fit_timer.stop();
2258 crate::blend::record_count("count.track.fit_sections", fit.sections as u64);
2259 crate::blend::record_count("count.track.local_inserted", inserted as u64);
2260 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2261 eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2262 eprintln!(
2263 "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2264 fit.miss,
2265 fit.samples,
2266 track.fractions.len(),
2267 track.breaks.len(),
2268 fit.on_floor
2269 );
2270 if let Some(worst) = fit.worst {
2275 let parameter = d0 + (d1 - d0) * worst.fraction;
2276 let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2277 .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2278 .map(|pair| [pair[0], pair[1]]).collect();
2279 let mean = (d1 - d0) / spans.len().max(1) as f64;
2280 let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2281 let width = span.map(|index| spans[index][1] - spans[index][0]);
2282 eprintln!(
2283 "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2284 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2285 worst.fraction, parameter, span, spans.len(), width, mean,
2286 worst.bracket[0], worst.bracket[1], worst.interval
2287 );
2288 }
2289 }
2290 if !fit.on_floor {
2297 if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2298 for (label, t) in [("start", d0), ("end", d1)] {
2301 if let Ok(point) = piece.evaluate(t) {
2302 let standoff = crate::project_point_to_surface(surface, point).map(|p| p.distance).unwrap_or(f64::NAN);
2303 eprintln!(
2304 "TRACK_FIT closed-chain piece {label} ({:.9}, {:.9}, {:.9}) stands {standoff:.3e} off the carrier",
2305 point.x, point.y, point.z
2306 );
2307 }
2308 }
2309 }
2310 return Err(KernelRefusal::non_convergence(KernelStage::Refine, CHAIN_PCURVE_FLOOR_WHAT, format!(
2311 "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2312 against a floor of {tolerance:.1e}",
2313 crate::blend::PCURVE_OFF_FLOOR,
2314 fit.miss,
2315 fit.samples
2316 )));
2317 }
2318 Ok(fit.curve)
2319}