Skip to main content

brep_kernel/blending/blend/chain/
closed.rs

1use crate::{KernelRefusal, KernelStage, OrRefuse};
2use super::*;
3
4/// Piece of a support rim assigned to one mate face.
5pub(super) struct RimPiece {
6    pub(super) face_id: u64,
7    pub(super) loop_index: usize,
8    /// Global-parameter window of this piece.
9    pub(super) window: [f64; 2],
10    pub(super) edge_id: u64,
11}
12
13/// Blend a chain of conjugated edges with one rolling-ball blend face
14/// (Golovanov §6.9.5: all conjugated edges processed together).  A CLOSED
15/// chain (stadium rim, T-pipe saddle) welds into a periodic blend; an OPEN
16/// chain (line->arc->line capped by end faces) rolls a clamped blend that
17/// terminates in a transverse edge on each end face.
18pub 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    // One blend operation: the reports of an earlier one, never consumed, are
26    // dropped as it starts (nested calls keep this operation's).
27    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
42/// `blend_smooth_chain` restricted to CLOSED chains: an edge that is one arc
43/// of a seam-split rim (a cyl×cyl saddle) is blended as the whole rim, but an
44/// OPEN tangent chain is not walked — filleting one edge of it must not
45/// fillet the unselected edges it runs into.  Tangent propagation is a
46/// selection policy for the feature layer; the kernel caps an unselected
47/// continuation instead (`network.rs`).
48pub 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    // One blend operation: the reports of an earlier one, never consumed, are
56    // dropped as it starts (nested calls keep this operation's).
57    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
74/// Blend a CLOSED chain of conjugated edges with one rolling-ball blend
75/// face (Golovanov §6.9.5: all conjugated edges processed together).
76///
77/// A wall whose ball-centre curve turns tighter than the ball FOLDS, and the
78/// envelope a plain march sweeps is then not the blend (`blend/fold.rs`).  It
79/// is not refused here any more where the fold is a lens the wall closes on
80/// itself: the chain is re-marched up a measured ladder of budgets and the lens
81/// is carved out along the crease (`blend/carve.rs`).  The re-march is the
82/// whole of the extra cost and only a folding wall pays it — the first march is
83/// the ordinary one, and its refusal is what asks for the second.
84fn blend_closed_smooth_chain(
85    solid: &BrepSolid,
86    segments: &[ChainSegment<'_>],
87    radius: f64,
88    chamfer: bool,
89    name: Option<&str>,
90) -> Result<BrepSolid, KernelRefusal> {
91    let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
92    // The plain collar climbs the same measured ladder the carve does, from the
93    // same first rung: a budget picked once is a resolution request in disguise.
94    let mut per_segment = CHAIN_PER_SEGMENT;
95    let first = loop {
96        let rung = build_closed_smooth_chain(
97            solid,
98            segments,
99            radius,
100            chamfer,
101            name,
102            per_segment,
103            FoldPolicy::Refuse,
104            bar,
105        );
106        match rung {
107            Ok(Rung::TooCoarse { deviation }) => {
108                crate::blend::carve::carve_trace(format_args!(
109                    "blend chain: {per_segment} stations per segment leaves the rails \
110                     {deviation:.6e} from the rolling ball, against {bar:.3e}"
111                ));
112                if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
113                    return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
114                        "blend: at {per_segment} stations per segment — the top of the chain's \
115                         ladder — this collar's fitted rails are still {deviation:.6e} from the \
116                         rolling ball's own contacts against a bar of {bar:.3e}, so there is no \
117                         wall accurate enough to build"
118                    )));
119                }
120                per_segment *= 2;
121            }
122            // The continuation Newton of one station did not converge from its
123            // neighbour at this spacing. Only that typed station failure climbs:
124            // the next rung seeds every station from a neighbour half as far
125            // away. At the top the march's own refusal is returned verbatim.
126            Ok(Rung::StationFailed(failure)) => {
127                crate::blend::carve::carve_trace(format_args!(
128                    "blend chain: {per_segment} stations per segment: segment {} position {} \
129                     (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
130                    failure.segment, failure.position, failure.t, failure.from_position,
131                    failure.from_t, failure.seed_uv, failure.error
132                ));
133                if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
134                    break Err(failure.error);
135                }
136                per_segment *= 2;
137            }
138            other => break other,
139        }
140    };
141    if let Err(error) = &first {
142        crate::blend::carve::carve_trace(format_args!(
143            "blend chain: the march at {per_segment} stations per segment refused: {error}"
144        ));
145    }
146    match first {
147        Ok(Rung::Built(built)) => Ok(built),
148        Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
149            Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
150        }
151        // A fold that reaches a RAIL is not a lens: every section of that wall
152        // is singular somewhere between its contacts, so there is nothing to cut
153        // away and keep, and the ladder below would only re-march it to its top
154        // before the carve refused the same band. Its refusal names the face it
155        // cannot cross, and it is terminal as it stands (`blend/fold.rs`).
156        Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
157        // A CHAMFER section is a straight line and has no ball centre to read
158        // the fold locus against, so its fold stays the refusal it was.
159        Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
160            carve_ladder(solid, segments, radius, name, bar)
161                // A CARVE THAT FAILS IS STILL A FOLD, and the fold is terminal.
162                // If the carve's own refusal did not read as one, the ladder
163                // below would answer a proven-folding centre curve with the
164                // cutter's unrelated complaint — which is the exact failure
165                // `fold.rs` was made terminal to stop. So every exit from this
166                // branch carries the fold's own prefix, with the carve's reason
167                // inside it.
168                .map_err(|carve_error| {
169                    if crate::blend::is_wall_fold(&carve_error) {
170                        carve_error
171                    } else {
172                        crate::blend::fold::wall_fold_refusal(
173                            Vec::new(),
174                            format!(
175                                "{} {} fits this edge: the wall folds, and the fold band could not \
176                                 be carved — {carve_error}",
177                                crate::blend::WALL_FOLDS,
178                                radius.abs()
179                            ),
180                        )
181                    }
182                })
183        }
184        Err(error) => Err(error),
185    }
186}
187
188/// Re-march a folding chain up a ladder of station budgets and carve at the
189/// first rung whose wall IS the rolling ball's to `bar`.
190///
191/// The budget is not chosen: a larger or a tighter fillet needs a different
192/// rung, and a number picked once from one fixture would be a resolution
193/// request in disguise. Each rung is measured by the rolling ball itself — its
194/// contacts re-solved halfway between stations, where an interpolant is worst —
195/// against the fitted rails, and the ladder stops at the first rung inside
196/// `intersection_fit`. A wall that still misses at the top is refused by name.
197fn carve_ladder(
198    solid: &BrepSolid,
199    segments: &[ChainSegment<'_>],
200    radius: f64,
201    name: Option<&str>,
202    bar: f64,
203) -> Result<BrepSolid, KernelRefusal> {
204    let mut per_segment = CHAIN_PER_SEGMENT;
205    loop {
206        match build_closed_smooth_chain(
207            solid,
208            segments,
209            radius,
210            false,
211            name,
212            per_segment,
213            FoldPolicy::Carve,
214            bar,
215        )? {
216            Rung::Built(built) => {
217                crate::blend::carve::carve_trace(format_args!(
218                    "blend carve: {per_segment} stations per segment is the first rung inside \
219                     {bar:.3e}"
220                ));
221                return Ok(built);
222            }
223            // The carve keeps its own contract: a station failure refuses as
224            // it did before the plain ladder learned to climb on it.
225            Rung::StationFailed(failure) => return Err(failure.error),
226            Rung::TooCoarse { deviation } => {
227                crate::blend::carve::carve_trace(format_args!(
228                    "blend carve: {per_segment} stations per segment leaves the rails \
229                     {deviation:.6e} from the rolling ball, against {bar:.3e}"
230                ));
231                if per_segment * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
232                    return Err(crate::blend::fold::wall_fold_refusal(
233                        Vec::new(),
234                        format!(
235                            "{} {} fits this edge: the wall folds, and at {per_segment} stations \
236                             per segment — the top of the carve's ladder — its fitted rails are \
237                             still {deviation:.6e} from the rolling ball's own contacts against a \
238                             bar of {bar:.3e}, so there is no wall accurate enough to carve",
239                            crate::blend::WALL_FOLDS,
240                            radius.abs()
241                        ),
242                    ));
243                }
244                per_segment *= 2;
245            }
246        }
247    }
248}
249
250/// One rung of a chain build: the finished solid, or — for a carve — the
251/// measured reason this budget is not enough.
252enum Rung {
253    Built(BrepSolid),
254    TooCoarse { deviation: f64 },
255    /// A station's continuation Newton found no root at this budget
256    /// (`march::ChainStationFailure`); a finer rung is a nearer seed.
257    StationFailed(ChainStationFailure),
258}
259
260/// The worst distance from the rolling ball's halfway contacts to the fitted
261/// rails `cr` and `cs`, each found by a closest-point Newton from the
262/// parameter the chord map assigns the midpoint.
263///
264/// THE NEWTON STAYS WHERE IT WAS SEEDED. The foot of a halfway contact is
265/// inside the station interval it was solved in, so every step is bounded by
266/// that interval's own length in the fit's parameter, clamped to the rail's
267/// domain, and taken only when it brings the rail closer. A plain Newton stops
268/// wherever the distance is stationary along the curve, which a rail that
269/// turns sharply at a chain junction offers well away from the foot. Measured
270/// on a 20 cube rounded at r = 3 with a bottom edge at r = 4 (the wall folds at
271/// the arcs' rails, so this is the ladder's reading with that refusal set
272/// aside): unguarded, every rung read the distance from the halfway contact
273/// beside a junction to the junction itself — 2.927709e-1 at 24 stations per
274/// segment against a half interval of 0.291667, halving per rung to
275/// 9.128591e-3 at 768 — and on two of the eight mirror-image edges a foot
276/// across the chain, 1.299089e1. Guarded, all eight read 6.961233e-3 at 24
277/// and 6.920006e-6 to 6.920008e-6 at 768, falling fourfold per rung as a
278/// cubic's does. A far foot can read short as well as long, which would pass
279/// a rung that misses.
280///
281/// Also returns every interval over `bar`, named by its left sample's index:
282/// what local refinement inserts a station into.
283fn rail_deviation(
284    samples: &[ChainSample],
285    cr: &NurbsCurve,
286    cs: &NurbsCurve,
287    fit_low: f64,
288    fit_range: f64,
289    per_segment: usize,
290    bar: f64,
291) -> Result<(f64, Vec<usize>), KernelRefusal> {
292    // Each interval runs from a station to the NEXT one along its segment — the
293    // next grid position, a station local refinement inserted, or the
294    // segment's own end station — in the order the stations take.
295    let mut order: Vec<usize> = (0..samples.len())
296        .filter(|&index| (0..=per_segment as isize).contains(&samples[index].position))
297        .collect();
298    order.sort_by(|&a, &b| {
299        let (a, b) = (&samples[a], &samples[b]);
300        a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
301    });
302    let next_of: std::collections::HashMap<usize, f64> = order
303        .windows(2)
304        .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
305        .map(|pair| (pair[0], samples[pair[1]].parameter))
306        .collect();
307    let mut failing = Vec::new();
308    let closest = closest_on_rail;
309    let mut worst: f64 = 0.0;
310    // Each segment's worst halfway contact, for `BREP_BLEND_CARVE_TRACE=1`: where
311    // along the chain a rung misses is what tells a junction from a bend.
312    let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
313    for &index in &order {
314        let sample = &samples[index];
315        let Some([p1, p2]) = sample.midpoint else {
316            continue;
317        };
318        let Some(&next) = next_of.get(&index) else {
319            continue;
320        };
321        let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
322        let seed = (global - fit_low) / fit_range;
323        let reach = (next - sample.parameter).abs() / fit_range;
324        let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
325        worst = worst.max(miss1).max(miss2);
326        if miss1.max(miss2) > bar {
327            failing.push(index);
328        }
329        if per_segment_worst.len() <= sample.segment {
330            per_segment_worst.resize(
331                sample.segment + 1,
332                (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
333            );
334        }
335        if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
336            per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
337        }
338    }
339    for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
340        crate::blend::carve::carve_trace(format_args!(
341            "  rail miss segment {segment}: worst {miss:.6e} after station {position} \
342             (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
343            p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
344        ));
345    }
346    Ok((worst, failing))
347}
348
349/// Most dense INTERVALS (4097 samples) a closed chain's support-piece track may
350/// hold: the existing ceiling its uniform track is built under, and the one
351/// local refinement of that track stays within — not a sample count.
352const PIECE_TRACK_CEILING: usize = 4096;
353
354/// Distance from `target` to a fitted rail, by a Newton that stays where it is
355/// seeded (`rail_deviation`'s reading, shared with the dense diagnostic).
356fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
357    let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
358    let distance = |u: f64| -> Result<f64, KernelRefusal> {
359        Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
360    };
361    let mut u = seed.clamp(low, high);
362    let mut best = distance(u)?;
363    for _ in 0..40 {
364        let derivatives = curve.derivatives_extended(u, 2).or_refuse(KernelStage::Refine, "derivatives_extended")?;
365        let offset = derivatives[0].sub(target);
366        let slope = offset.dot(derivatives[1]);
367        let curvature = derivatives[1].dot(derivatives[1]) + offset.dot(derivatives[2]);
368        // Off a minimum's basin the Newton step points the wrong way; walk
369        // downhill by the bound instead and let the test below shorten it.
370        let mut step = if curvature > 0.0 {
371            (slope / curvature).clamp(-reach, reach)
372        } else {
373            reach.copysign(slope)
374        };
375        let mut trial = (u - step).clamp(low, high);
376        let mut trial_distance = distance(trial)?;
377        for _ in 0..30 {
378            if trial_distance < best {
379                break;
380            }
381            step *= 0.5;
382            trial = (u - step).clamp(low, high);
383            trial_distance = distance(trial)?;
384        }
385        if !(trial_distance < best) {
386            break;
387        }
388        let moved = (trial - u).abs();
389        u = trial;
390        best = trial_distance;
391        if moved < 1e-15 {
392            break;
393        }
394    }
395    Ok(best)
396}
397
398/// Points of every interval the LOCAL-INSERTION path judges a chain rung's
399/// rails on: j/16, j = 1..15 (the halfway point among them).
400pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
401
402/// One dense contact pair's miss: each rail's raw reading is checked BEFORE
403/// they are joined (`f64::max` alone drops a NaN on either side); `None` when
404/// either is unreadable.
405fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
406    (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
407}
408
409/// Each solved dense contact's distance to the fitted rails, read with the same
410/// closest-point Newton and seeds `rail_deviation` uses, keyed by the
411/// interval's left sample: the readable misses, and separately every interval
412/// with an unreadable reading on either rail (a non-finite distance, or a foot
413/// that could not be evaluated).
414#[allow(clippy::too_many_arguments)]
415pub(super) fn dense_rail_misses(
416    contacts: &crate::blend::chain::DenseContacts,
417    samples: &[ChainSample],
418    cr: &NurbsCurve,
419    cs: &NurbsCurve,
420    low: f64,
421    range: f64,
422    wrap: bool,
423) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
424    let mut readable = Vec::with_capacity(contacts.solved.len());
425    let mut unreadable = Vec::new();
426    for (left, next, fraction, p1, p2) in &contacts.solved {
427        let (from, to) = (samples[*left].parameter, samples[*next].parameter);
428        // A closed chain's chord parameter wraps with period 1; an OPEN
429        // chain's runs straight through its ends (its last station's may
430        // overshoot 1, and wrapping it would read the far end).
431        let global = from + fraction * (to - from);
432        let global = if wrap { global.rem_euclid(1.0) } else { global };
433        let seed = (global - low) / range;
434        let reach = (to - from).abs() / range;
435        // A closest foot that cannot be evaluated is an unreadable reading of
436        // this interval, not a refusal of the rung.
437        let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
438        let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
439        match rail_pair_miss(cr_miss, cs_miss) {
440            Some(miss) => readable.push((*left, miss)),
441            None => unreadable.push(*left),
442        }
443    }
444    Ok((readable, unreadable))
445}
446
447/// The local-insertion path's rail verdict: the halfway reading joined with the
448/// dense stencil's, FAIL-CLOSED — an interval with a dense contact over `bar`,
449/// an unreadable distance, or a contact that was not read (unsolved, or
450/// unreadable on either rail), fails; the worst is the largest readable
451/// distance. Intervals are named by their left sample.
452pub(super) fn combine_rail_verdict(
453    halfway: (f64, Vec<usize>),
454    dense: &[(usize, f64)],
455    unsolved: &[usize],
456    bar: f64,
457) -> (f64, Vec<usize>) {
458    let (mut worst, mut failing) = halfway;
459    for &(left, miss) in dense {
460        if !miss.is_finite() {
461            failing.push(left);
462            continue;
463        }
464        if miss > bar {
465            failing.push(left);
466        }
467        worst = worst.max(miss);
468    }
469    failing.extend_from_slice(unsolved);
470    failing.sort_unstable();
471    failing.dedup();
472    (worst, failing)
473}
474
475/// Whether the construction request runs on an accepted plain fillet rung:
476/// always, outside the controls that turn it off on their own thread
477/// (`MODEL_REQUEST_OFF`) or that turn the wall verdict off.
478/// The construction request's round limit: the existing local-refinement
479/// budget, outside the control that sets it to zero on its own thread to
480/// force the typed unmet report.
481pub(super) fn model_rounds_limit() -> usize {
482    REFINE_ROUNDS
483}
484
485pub(super) fn model_request_on() -> bool {
486    wall_verdict_on()
487}
488
489
490/// Whether the fillet wall verdict runs: always, outside the mutation control
491/// that turns it off on its own thread (`WALL_VERDICT_OFF`).
492pub(super) fn wall_verdict_on() -> bool {
493    true
494}
495
496
497
498
499
500/// Section-angle probes across a fillet wall, as fractions of its `v` span:
501/// k/8 for k = 1..7, between (never on) the rails the rail verdict reads.
502pub(super) const WALL_SECTION_PROBES: usize = 8;
503
504/// How far `p` stands from the rolling ball's sweep, read at p's OWN
505/// characteristic: the stationary point of |p − C(s)|, C the cubic through
506/// four consecutive exact ball centres of `nodes` around node `j` (`nodes`
507/// uniform in fraction across one interval, its ends the interval's own
508/// stations), searched within one node of `j`, minus `radius`.
509///
510/// The cubic is an APPROXIMATION of the centre curve, not the exact sweep:
511/// its error grows with the node spacing and the spine's higher derivatives,
512/// and [`wall_probe_reading`] estimates it by a second, overlapping stencil.
513///
514/// A STATIONARY point, not the nearest centre: past a fold's focal point an
515/// envelope point is a MAXIMUM of its distance to the spine, and a global
516/// nearest centre would read a carved wall's lens as inside another ball. A
517/// centre at the probe's own chain parameter is not used either: the wall's
518/// `u` and the stencil's edge parameter name the same section only to first
519/// order, and a spine offset ℓ adds √(ρ² + ℓ²) − ρ to the reading.
520///
521/// `None` — unreadable, which fails its interval — for a non-finite input,
522/// fewer than four nodes, a probe at an interval end, or a solve that leaves
523/// its window or does not converge.
524pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
525    let j = locate_foot(p, nodes, j, radius)?;
526    sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
527}
528
529/// How many nodes the stencil may follow a probe's foot away from the probe's
530/// own node.
531const FOOT_TRAVEL: usize = 4;
532
533/// The node whose stencil holds `p`'s stationary point, starting from node
534/// `j`: the probe's wall `u` and the centres' edge parameter name the same
535/// section only to first order, so the foot can sit more than the one-node
536/// window from `j` (the 15- and 20-degree rows' feet stood up to 1.5 nodes,
537/// one 3.3, from it on f4ed17782's build, every decline there a window exit).
538/// A solve that leaves its window moves the stencil to the node nearest
539/// where it left, at most [`FOOT_TRAVEL`] nodes from `j` and never back to a
540/// node already tried; `None` when the foot is not found that way (or a
541/// solve is unreadable), which fails its interval as before.
542fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
543    let n = nodes.len();
544    if n < 4 || j == 0 || j + 1 >= n {
545        return None;
546    }
547    let mut tried: Vec<usize> = Vec::new();
548    let mut at = j;
549    loop {
550        let base = (at - 1).min(n - 4);
551        match stencil_solve(p, nodes, base, at, radius) {
552            Solve::Read(_) => return Some(at),
553            Solve::Unread => return None,
554            Solve::Window(s) => {
555                tried.push(at);
556                let next = (base as f64 + s).round();
557                let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
558                if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
559                    wall_decline(|| format!(
560                        "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
561                        if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
562                    ));
563                    return None;
564                }
565                at = next;
566            }
567        }
568    }
569}
570
571/// [`sweep_residual`] on the cubic through `nodes[base..base + 4]`, the
572/// stationary point searched from node `j` within one node of it.
573fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
574    match stencil_solve(p, nodes, base, j, radius) {
575        Solve::Read(distance) => Some(distance),
576        Solve::Window(_) | Solve::Unread => None,
577    }
578}
579
580/// One stencil's solve: the distance read, or where (in node units from
581/// `base`) the search left its one-node window, or unreadable.
582enum Solve {
583    Read(f64),
584    Window(f64),
585    Unread,
586}
587
588thread_local! {
589    /// When Some, every decline inside the declared wall reader appends WHY
590    /// (innermost first): set and taken by `declared_wall_probe_why` only, so
591    /// every other caller pays one thread-local read per decline and the
592    /// reader's arithmetic and outcomes are the same either way.
593    static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
594}
595
596fn wall_decline(reason: impl FnOnce() -> String) {
597    WALL_DECLINES.with(|log| {
598        if let Some(lines) = log.borrow_mut().as_mut() {
599            lines.push(reason());
600        }
601    });
602}
603
604/// [`declared_wall_probe`] with, when it declines, the chain of declines
605/// that made it (innermost first): which tier (fine / coarse), which
606/// stencil, and inside the solver the exit (unreadable input, a vanishing
607/// slope, a window exit with the overshot `s`, its offset from the window and
608/// whether the solver's OWN stationarity predicate already held at the
609/// overshot point, a stall, a non-finite distance) or the foot search's
610/// tried nodes. Exactly the same arithmetic as `declared_wall_probe`.
611pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
612    let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
613    let read = declared_wall_probe(p, fine, jf, radius);
614    let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
615    read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
616}
617
618fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
619    let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
620    let n = nodes.len();
621    if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
622        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
623        return Solve::Unread;
624    }
625    // Newton's divided differences on 0, 1, 2, 3, s in node units from `base`.
626    let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
627    let d1 = p1.sub(p0);
628    let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
629    let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
630    let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
631        let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
632        let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
633        let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
634        (c, c1, c2)
635    };
636    let start = j as f64 - base as f64;
637    let mut s = start;
638    let mut converged = false;
639    for _ in 0..50 {
640        let (c, c1, c2) = curve(s);
641        let w = c.sub(p);
642        let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
643        if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
644            wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
645            return Solve::Unread;
646        }
647        // Converged at THIS s once the stationarity condition holds to its own
648        // rounding: g = (C − p)·C′ carries an absolute error of a few
649        // ε·(|p| + |C|)·|C′|, so at a node spacing of 1e-4 the Newton step
650        // cannot fall below ~1e-11 node units and the relative step test
651        // alone never fires — a line-for-line model of this reader read 87 %
652        // to 97 % of the NX and skew walls' fine-stencil probes unreadable
653        // that way (its Newton run past the stop read |g| / (ε·(|p| + |C|)·|C′|)
654        // at most 0.91 there; 16 keeps a margin). Tested before the step, so
655        // the point read is the point whose stationarity was checked; its
656        // distance is off the true stationary one by about
657        // ½·|g|·δs/|C − p| ≤ 8ε·(|p| + |C|)·|C′|/ρ for δs within the
658        // one-node window, far below any bar.
659        let stationary = g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
660        if stationary {
661            converged = true;
662            break;
663        }
664        let step = g / slope;
665        let before = s;
666        s -= step;
667        if !((s - start).abs() <= 1.0) {
668            wall_decline(|| {
669                // The solver's own stationarity predicate, read at the
670                // overshot point (diagnostic only: nothing decides on it).
671                let (c, c1, _) = curve(s);
672                let g_after = c.sub(p).dot(c1);
673                let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
674                format!(
675                    "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
676                     g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
677                    (s - start).abs() - 1.0
678                )
679            });
680            return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
681        }
682        if step.abs() <= 1e-14 * (1.0 + s.abs()) {
683            converged = true;
684            break;
685        }
686    }
687    if !converged {
688        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
689        return Solve::Unread;
690    }
691    let distance = curve(s).0.sub(p).length() - radius;
692    if distance.is_finite() {
693        Solve::Read(distance)
694    } else {
695        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
696        Solve::Unread
697    }
698}
699
700/// One wall probe read on TWO overlapping four-centre stencils around node
701/// `j` (the one [`sweep_residual`] uses and its neighbour): `(worst, spread)`
702/// — the larger magnitude, which keeps ONE stencil's under-read from passing,
703/// and how far apart the two readings are. Both are PROXY readings: when both
704/// stencils under-read alike (a symmetric node set about a foot near the
705/// focal distance) the larger still under-reads and the spread reads zero,
706/// so neither number is an error bound. `None` when either stencil is
707/// unreadable.
708pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
709    let n = nodes.len();
710    if n < 5 || j == 0 || j + 1 >= n {
711        return None;
712    }
713    let j = locate_foot(p, nodes, j, radius)?;
714    let first = (j - 1).min(n - 4);
715    let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
716    let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
717    let a = sweep_residual_on(p, nodes, first, j, radius)?;
718    let b = sweep_residual_on(p, nodes, second, j, radius)?;
719    Some((a.abs().max(b.abs()), (a - b).abs()))
720}
721
722/// The DECLARED wall reading at one probe, the one rule both wall tiers
723/// decide on: `(reading, uncertainty)`, declared = reading + uncertainty.
724/// `fine` are the interval's centres at twice the probe density, `jf` the
725/// probe's (even) node among them; the coarse set is every other node. The
726/// reading is the larger fine two-stencil magnitude; the uncertainty is the
727/// larger of the fine stencils' spread and how far the fine reading moved
728/// from the coarse one — the proxy's error estimated from its own change with
729/// node density, which two stencils agreeing at one density cannot show (two
730/// symmetric stencils that both under-read read a spread of zero). An
731/// ESTIMATE, not a bound: a proxy whose error does not shrink with density
732/// would still hide. `None` when any reading is unreadable.
733pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
734    if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
735        wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
736        return None;
737    }
738    let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
739        wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
740        return None;
741    };
742    let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
743        wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
744        return None;
745    };
746    let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
747    let Some(at_coarse) = sweep_residual(p, &coarse, jf / 2, radius) else {
748        wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
749        return None;
750    };
751    let uncertainty = spread.max((at_fine - at_coarse).abs());
752    uncertainty.is_finite().then_some((reading, uncertainty))
753}
754
755/// The even fine node whose position along the CENTRE PATH (cumulative
756/// chord of `fine`) is nearest `fraction` of it, inside 2..=n-3; `None` when
757/// the path has no finite positive length.
758///
759/// A wall probe sits at a fraction of its interval in the WALL's parameter,
760/// which is chord-like; the exact centres sit at uniform fractions of the
761/// edge's NATIVE parameter. On a regular edge the two agree to first order and
762/// the probe's foot lies within a node or two of its own node. At a
763/// STATIONARY section (the edge's derivative vanishes there) the centre moves
764/// like t^2 in the native parameter, so a probe at chord fraction f has its foot
765/// near native fraction sqrt(f): the probe at node 2 of 32 (f = 1/16) has its
766/// foot near node 8, past `FOOT_TRAVEL`, and the reader declines though the
767/// centre locus is perfectly regular there (zero_tangent_sites' stationary
768/// stadium and tombstone rims refused that way with every reading under the
769/// bar).
770pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
771    let n = fine.len();
772    if n < 7 || !(fraction.is_finite()) {
773        return None;
774    }
775    let mut cumulative = Vec::with_capacity(n);
776    cumulative.push(0.0_f64);
777    for pair in fine.windows(2) {
778        let step = pair[1].sub(pair[0]).length();
779        cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
780    }
781    let total = *cumulative.last()?;
782    if !(total.is_finite() && total > 0.0) {
783        return None;
784    }
785    let target = fraction.clamp(0.0, 1.0) * total;
786    let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
787    let even = 2 * ((nearest as f64 / 2.0).round() as usize);
788    Some(even.clamp(2, (n - 3) & !1))
789}
790
791/// [`declared_wall_probe_why`] from the probe's own node `jf`; when that
792/// declines, once more from the node the probe's `fraction` reaches along the
793/// centre path ([`centre_path_node`]), when that is another node. The same
794/// reader, bars and fail-closed rule decide; a probe it still cannot read is
795/// unread (both declines named). A probe the first read reads is read exactly
796/// as before.
797pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
798    let first = match declared_wall_probe_why(p, fine, jf, radius) {
799        Ok(read) => return Ok(read),
800        Err(why) => why,
801    };
802    match centre_path_node(fine, fraction) {
803        Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
804            .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
805        _ => Err(first),
806    }
807}
808
809/// [`declared_wall_probe_located_why`] as an `Option`.
810pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
811    match declared_wall_probe(p, fine, jf, radius) {
812        Some(read) => Some(read),
813        None => match centre_path_node(fine, fraction) {
814            Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
815            _ => None,
816        },
817    }
818}
819
820/// [`declared_wall_probe_why`] on centres at NON-UNIFORM native abscissae
821/// `ts` (one per node, strictly increasing): the same two-tier rule, the
822/// same two overlapping four-centre stencils, the same foot search
823/// (`FOOT_TRAVEL`), the same 16·ε stationarity predicate and the same
824/// window (one node either side of the start node), with each cubic the
825/// Newton interpolant on the nodes' OWN abscissae instead of the implicit
826/// 0, 1, 2, 3. For a junction of two intervals of unequal native width (local
827/// insertions make them so): the uniform reader would place the two halves'
828/// nodes at one spacing and distort the centre locus. The coarse tier is
829/// every other node at every other abscissa. A caller whose abscissae are
830/// uniform uses the uniform reader, bit for bit.
831pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
832    if jf % 2 != 0
833        || fine.len() < 9
834        || fine.len() % 2 == 0
835        || ts.len() != fine.len()
836        || !ts.iter().all(|t| t.is_finite())
837        || !ts.windows(2).all(|pair| pair[1] > pair[0])
838    {
839        return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
840    }
841    let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
842    let read = (|| -> Option<(f64, f64)> {
843        let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
844            wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
845            None
846        })?;
847        let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
848            wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
849            None
850        })?;
851        let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
852        let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
853        let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
854            wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
855            None
856        })?;
857        let uncertainty = spread.max((at_fine - at_coarse).abs());
858        uncertainty.is_finite().then_some((reading, uncertainty))
859    })();
860    let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
861    read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
862}
863
864fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
865    let n = nodes.len();
866    if n < 5 || j == 0 || j + 1 >= n {
867        return None;
868    }
869    let j = locate_foot_at(p, nodes, ts, j, radius)?;
870    let first = (j - 1).min(n - 4);
871    let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
872    let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
873    let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
874    let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
875    Some((a.abs().max(b.abs()), (a - b).abs()))
876}
877
878fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
879    let j = locate_foot_at(p, nodes, ts, j, radius)?;
880    match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
881        Solve::Read(distance) => Some(distance),
882        Solve::Window(_) | Solve::Unread => None,
883    }
884}
885
886/// [`locate_foot`] on native abscissae: a window exit moves the stencil to
887/// the node whose abscissa is nearest the exit, within `FOOT_TRAVEL` nodes
888/// and never back to a node already tried.
889fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
890    let n = nodes.len();
891    if n < 4 || j == 0 || j + 1 >= n {
892        return None;
893    }
894    let mut tried: Vec<usize> = Vec::new();
895    let mut at = j;
896    loop {
897        let base = (at - 1).min(n - 4);
898        match stencil_solve_at(p, nodes, ts, base, at, radius) {
899            Solve::Read(_) => return Some(at),
900            Solve::Unread => return None,
901            Solve::Window(t) => {
902                tried.push(at);
903                let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
904                let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
905                if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
906                    wall_decline(|| format!(
907                        "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
908                        if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
909                    ));
910                    return None;
911                }
912                at = next;
913            }
914        }
915    }
916}
917
918/// [`stencil_solve`] on the four nodes `base..base + 4` at their native
919/// abscissae: the Newton (divided-difference) cubic in `t`, searched from
920/// node `j`'s abscissa within one node either side; `Window` carries the
921/// abscissa it left at.
922fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
923    let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
924    let n = nodes.len();
925    if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || ts.len() != n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
926        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
927        return Solve::Unread;
928    }
929    let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
930    let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
931    let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
932    let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
933    let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
934    let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
935    let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
936    let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
937    let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
938        let (a, b, c) = (t - t0, t - t1, t - t2);
939        let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
940        let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
941        let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
942        (point, first, second)
943    };
944    let (low, high) = (ts[j - 1], ts[j + 1]);
945    let start = ts[j];
946    let mut t = start;
947    let mut converged = false;
948    for _ in 0..50 {
949        let (c, c1, c2) = curve(t);
950        let w = c.sub(p);
951        let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
952        if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
953            wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
954            return Solve::Unread;
955        }
956        if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
957            converged = true;
958            break;
959        }
960        let step = g / slope;
961        let before = t;
962        t -= step;
963        if !(t >= low && t <= high) {
964            wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
965            return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
966        }
967        // The uniform reader's stop, 1e-14 · (1 + |s|) in node units from
968        // the stencil's base, read in LOCAL node units here (one node = half
969        // the window), so it does not depend on where the native t's origin
970        // is or on its scale.
971        let node = 0.5 * (high - low);
972        if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
973            converged = true;
974            break;
975        }
976    }
977    if !converged {
978        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
979        return Solve::Unread;
980    }
981    let distance = curve(t).0.sub(p).length() - radius;
982    if distance.is_finite() {
983        Solve::Read(distance)
984    } else {
985        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
986        Solve::Unread
987    }
988}
989
990/// The intervals the construction request refines (and, when it cannot,
991/// reports): every interval whose DECLARED wall reading — as
992/// `wall_interior_misses` returns it, `(left, declared, reading,
993/// uncertainty)` — exceeds `model`. The one decision the request makes.
994pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
995    readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
996}
997
998/// Every fillet-wall interval's worst distance from the rolling ball's sweep,
999/// read on the BUILT wall at the even fractions of the centre stencil
1000/// j/`subdivisions` (the centres at twice the probe density) and section
1001/// angles k/`probes` ([`declared_wall_probe`]), keyed by the interval's left
1002/// sample: each readable interval's largest DECLARED reading, largest reading
1003/// and largest uncertainty, and separately
1004/// every interval with an unreadable reading or a contact the stencil did not
1005/// solve. The wall shares the rails' `u` (one interpolation parameter for its
1006/// three rows).
1007#[allow(clippy::too_many_arguments)]
1008pub(super) fn wall_interior_misses(
1009    contacts: &crate::blend::chain::DenseContacts,
1010    samples: &[ChainSample],
1011    surface: &crate::NurbsSurface,
1012    low: f64,
1013    range: f64,
1014    radius: f64,
1015    subdivisions: usize,
1016    probes: usize,
1017    wrap: bool,
1018) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1019    use std::collections::BTreeMap;
1020    // Each interval's centres at fractions 0, 1/k, …, 1: its stations at the
1021    // ends, the stencil's solved balls between.
1022    let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1023    for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1024        let entry = nodes.entry(*left).or_insert_with(|| {
1025            let mut row = vec![None; subdivisions + 1];
1026            row[0] = Some(samples[*left].station.center);
1027            row[subdivisions] = Some(samples[*next].station.center);
1028            (*next, row)
1029        });
1030        let step = (fraction * subdivisions as f64).round() as usize;
1031        if (1..subdivisions).contains(&step) {
1032            entry.1[step] = contacts.centers.get(index).copied().flatten();
1033        }
1034    }
1035    let mut readable = Vec::new();
1036    let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1037    for (&left, (next, row)) in &nodes {
1038        let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1039            unreadable.push(left);
1040            continue;
1041        };
1042        let (from, to) = (samples[left].parameter, samples[*next].parameter);
1043        let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1044        let mut read = true;
1045        'probes: for step in (2..subdivisions).step_by(2) {
1046            let global = from + (step as f64 / subdivisions as f64) * (to - from);
1047            let global = if wrap { global.rem_euclid(1.0) } else { global };
1048            let u = (global - low) / range;
1049            for k in 1..probes {
1050                let v = k as f64 / probes as f64;
1051                let reading = surface.evaluate(u, v).ok().and_then(|p| {
1052                    declared_wall_probe_located(p, &centres, step, step as f64 / subdivisions as f64, radius)
1053                });
1054                match reading {
1055                    Some((distance, gap)) => {
1056                        declared = declared.max(distance + gap);
1057                        worst = worst.max(distance);
1058                        spread = spread.max(gap);
1059                    }
1060                    None => {
1061                        read = false;
1062                        break 'probes;
1063                    }
1064                }
1065            }
1066        }
1067        if read {
1068            readable.push((left, declared, worst, spread));
1069        } else {
1070            unreadable.push(left);
1071        }
1072    }
1073    unreadable.sort_unstable();
1074    unreadable.dedup();
1075    (readable, unreadable)
1076}
1077
1078
1079fn build_closed_smooth_chain(
1080    solid: &BrepSolid,
1081    segments: &[ChainSegment<'_>],
1082    radius: f64,
1083    chamfer: bool,
1084    name: Option<&str>,
1085    per_segment: usize,
1086    fold: FoldPolicy,
1087    bar: f64,
1088) -> Result<Rung, KernelRefusal> {
1089    let samples = match march_chain_typed(segments, radius, false, per_segment, fold) {
1090        Ok(samples) => samples,
1091        Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1092        // Only local refinement reports incoherence; the march never does, so
1093        // this is its refusal exactly as `march_chain` would return it.
1094        Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1095        Err(ChainMarchError::Other(error)) => return Err(error),
1096    };
1097    CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1098    build_closed_chain_from_samples(solid, segments, radius, chamfer, name, per_segment, fold, bar, samples, 0)
1099}
1100
1101
1102thread_local! {
1103    /// The deepest local round any rung of the current ladder step attempted
1104    /// (reset as each rung's march is built from), so a rejected refinement
1105    /// reports the rounds it actually spent.
1106    static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1107}
1108
1109/// Fit, measure and build one rung from its stations. A plain (non-carve)
1110/// rung whose fitted rails miss the halfway contacts by more than `bar`
1111/// first REFINES LOCALLY: one continuation station into each failing interval
1112/// (`march::insert_chain_stations`), up to [`REFINE_ROUNDS`] rounds and never
1113/// past [`CHAIN_CARVE_MAX_PER_SEGMENT`] stations in a segment, then measures
1114/// again. Only a rung whose local rounds are spent is too coarse, and the
1115/// ladder climbs as it always has. A rung that passes at once is built from
1116/// exactly the stations it was built from before.
1117#[allow(clippy::too_many_arguments)]
1118fn build_closed_chain_from_samples(
1119    solid: &BrepSolid,
1120    segments: &[ChainSegment<'_>],
1121    radius: f64,
1122    chamfer: bool,
1123    name: Option<&str>,
1124    per_segment: usize,
1125    fold: FoldPolicy,
1126    bar: f64,
1127    mut samples: Vec<ChainSample>,
1128    round: usize,
1129) -> Result<Rung, KernelRefusal> {
1130    CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1131
1132    // Global rows from the in-segment samples (wrapped-overlap closed fit).
1133    let mut global: Vec<(f64, &ChainSample)> = samples
1134        .iter()
1135        .filter(|sample| (0..per_segment as isize).contains(&sample.position))
1136        .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1137        .collect();
1138    global.sort_by(|a, b| a.0.total_cmp(&b.0));
1139    let count = global.len();
1140    let degree = FIT_DEGREE;
1141    let mut extended_params = Vec::new();
1142    let mut samples_cr = Vec::new();
1143    let mut samples_mid = Vec::new();
1144    let mut samples_cs = Vec::new();
1145    let mut push = |station: &Station, parameter: f64| {
1146        extended_params.push(parameter);
1147        samples_cr.push(Vec4::from_point(station.p1, 1.0));
1148        samples_cs.push(Vec4::from_point(station.p2, 1.0));
1149        samples_mid.push(Vec4 {
1150            x: station.apex.x * station.weight,
1151            y: station.apex.y * station.weight,
1152            z: station.apex.z * station.weight,
1153            w: station.weight,
1154        });
1155    };
1156    for offset in (1..=degree).rev() {
1157        let (parameter, sample) = &global[count - offset];
1158        push(&sample.station, parameter - 1.0);
1159    }
1160    for (parameter, sample) in &global {
1161        push(&sample.station, *parameter);
1162    }
1163    // Close the loop: repeat the first station at parameter 1, then the
1164    // wrap continuation.
1165    push(&global[0].1.station, global[0].0 + 1.0);
1166    for offset in 1..=degree {
1167        let (parameter, sample) = &global[offset];
1168        push(&sample.station, parameter + 1.0);
1169    }
1170    let low = extended_params[0];
1171    let high = *extended_params.last().unwrap();
1172    let range = high - low;
1173    let normalized: Vec<f64> = extended_params
1174        .iter()
1175        .map(|parameter| (parameter - low) / range)
1176        .collect();
1177    let seam_low = (0.0 - low) / range;
1178    let seam_high = (1.0 - low) / range;
1179    let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1180        let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1181        let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1182        let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1183        Ok(middle)
1184    };
1185    let cr = fit_row(&samples_cr)?;
1186    let cs = fit_row(&samples_cs)?;
1187    let mid = if chamfer {
1188        None
1189    } else {
1190        Some(fit_row(&samples_mid)?)
1191    };
1192    let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1193    let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1194
1195    // The CARVE.  A wall that folds carries a lens of parameter inside the
1196    // swept ball volume; the blend is that wall with the lens cut away along
1197    // the crease where the two sheets meet, and the crease becomes an edge of
1198    // the blend with two coedges on this one face.  Traced on the surface that
1199    // was just fitted — not on the ideal envelope — so the two kept sheets
1200    // meet EXACTLY on the curve the trim is built from (`blend/carve.rs`).
1201    // EVERY rung is measured before anything is built on it, carve or not: the
1202    // rails are interpolants through the stations, and between stations they
1203    // can leave the rolling ball's own contacts — and so the carriers they are
1204    // supposed to ride — by far more than the kernel's intersection-fit
1205    // contract. A carve pays for surgery only on an accurate rung; a plain
1206    // collar pays nothing more than the measurement unless it misses.
1207    let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, per_segment, bar)?;
1208    // A rung that local insertion has refined is judged on the dense stencil of
1209    // exact contacts, not the halfway point alone (the NX collar's rung 24 read
1210    // 1.943844e-6 halfway and 2.009538e-6 at 7/16 against 2e-6). A rung that
1211    // passes without insertion is read exactly as it always was.
1212    // A FILLET wall is also judged between its rails, on every rung (carve
1213    // included): the rails can pass while the wall's interior — its weighted
1214    // apex row interpolated between stations — stands off the rolling ball
1215    // (the 20-degree seam-split row's accepted wall read 2.4335e-6 on a 16 x 9
1216    // probe grid and 2.4430e-6 at its worst interval's midpoint, mid-section,
1217    // against 2e-6). A chamfer is not a rolling ball and is not judged
1218    // against one. One stencil of exact contacts serves rails and wall.
1219    let stencil = if round > 0 {
1220        Some(crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, LOCAL_RAIL_STENCIL)?)
1221    } else {
1222        None
1223    };
1224    // The wall's own centres, at twice the rail stencil's density: its probes
1225    // sit at the rail stencil's fractions, and each is read at both densities
1226    // (`declared_wall_probe`). The rails keep their own stencil, read exactly
1227    // as before.
1228    let wall_stencil = if !chamfer && wall_verdict_on() {
1229        Some(crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1230    } else {
1231        None
1232    };
1233    let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1234        let contacts = stencil.as_ref().expect("the local path solves the stencil");
1235        let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1236        let reads = dense.len() + unreadable.len() + contacts.unsolved;
1237        let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1238        let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1239        crate::blend::carve::carve_trace(format_args!(
1240            "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1241             {} failing interval(s), {} unsolved contact(s)",
1242            verdict.0, verdict.1.len(), contacts.unsolved
1243        ));
1244        (verdict.0, verdict.1, reads)
1245    } else {
1246        (halfway, halfway_failing, 0)
1247    };
1248    // The wall between the rails, joined FAIL-CLOSED: an interval whose wall
1249    // reads over `bar`, or cannot be read, fails exactly as a rail interval
1250    // does, and the ladder (local insertion, climb, the existing refusal)
1251    // decides as it always has. Nothing else about the rung changes: a wall
1252    // inside `bar` adds reads, never stations.
1253    let mut _wall_worst = 0.0_f64;
1254    let mut _wall_spread = 0.0_f64;
1255    let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1256    if let Some(contacts) = wall_stencil.as_ref() {
1257        let (wall, unreadable) =
1258            wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1259        wall_readings = wall.clone();
1260        // Acceptance on the DECLARED reading against the unchanged bar.
1261        for &(left, declared, _, uncertainty) in &wall {
1262            if declared > bar {
1263                failing.push(left);
1264            }
1265            _wall_worst = _wall_worst.max(declared);
1266            _wall_spread = _wall_spread.max(uncertainty);
1267        }
1268        failing.extend_from_slice(&unreadable);
1269        failing.sort_unstable();
1270        failing.dedup();
1271        deviation = deviation.max(_wall_worst);
1272        crate::blend::carve::carve_trace(format_args!(
1273            "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1274             interval(s), {} unreadable, bar {bar:.3e}",
1275            wall.len(),
1276            unreadable.len()
1277        ));
1278    }
1279    if deviation > bar || !failing.is_empty() {
1280        let in_segment = |segment: usize| {
1281            samples
1282                .iter()
1283                .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1284                .count()
1285        };
1286        let within_ceiling = (0..segments.len()).all(|segment| {
1287            in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1288                <= CHAIN_CARVE_MAX_PER_SEGMENT
1289        });
1290        crate::blend::carve::carve_trace(format_args!(
1291            "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1292             worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1293            samples.len(),
1294            failing.len(),
1295            if within_ceiling { "ok" } else { "reached" }
1296        ));
1297        if fold != FoldPolicy::Refuse || round >= REFINE_ROUNDS || !within_ceiling {
1298            return Ok(Rung::TooCoarse { deviation });
1299        }
1300        // The insertion is the next round's attempt, whatever it returns.
1301        CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1302        match insert_chain_stations(segments, radius, per_segment, &mut samples, &failing) {
1303            Ok(()) => {}
1304            Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1305            // An inserted station that does not reach its retained right
1306            // neighbour's branch: this rung cannot be refined coherently, so it
1307            // is as coarse as its measured rails say, and the ladder decides.
1308            Err(ChainMarchError::Incoherent(failure)) => {
1309                crate::blend::carve::carve_trace(format_args!(
1310                    "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1311                     at t {:.9} does not reach its retained right neighbour: {}",
1312                    failure.segment, failure.t, failure.error
1313                ));
1314                return Ok(Rung::TooCoarse { deviation });
1315            }
1316            Err(ChainMarchError::Other(error)) => return Err(error),
1317        }
1318        assign_chain_parameters(&mut samples, segments.len(), per_segment, false);
1319        return build_closed_chain_from_samples(
1320            solid, segments, radius, chamfer, name, per_segment, fold, bar, samples, round + 1,
1321        );
1322    }
1323    crate::blend::carve::carve_trace(format_args!(
1324        "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1325         from the rolling ball, inside {bar:.3e}"
1326    ));
1327    // The CONSTRUCTION request, plain fillet rungs only: an ACCEPTED wall is
1328    // refined toward `KernelTolerances::model` — the identity tolerance, the
1329    // same 1e-7 the pcurve construction floor refines to under its looser
1330    // acceptance — locally, in exactly the intervals whose DECLARED reading
1331    // (`declared_wall_probe`: the proxy reading plus its uncertainty from
1332    // stencil spread AND change with node density) is over it: one rule for
1333    // the acceptance verdict, this request and its report. An estimate with
1334    // its uncertainty, NOT a proven bound.
1335    // Acceptance is unchanged: the refined rung must pass it again, and every
1336    // retained station, branch guard and ceiling holds. A wall already inside
1337    // `model` everywhere inserts nothing and builds exactly as it did. When
1338    // refinement cannot run or does not build, THIS accepted rung stands and
1339    // says so, typed (`WallModelReport`); never silently.
1340    let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1341    let stations_now = (0..segments.len())
1342        .map(|segment| {
1343            samples
1344                .iter()
1345                .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1346                .count()
1347        })
1348        .max()
1349        .unwrap_or(0);
1350    let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1351    // What the request found short, named in its report: the detail, which
1352    // component the residual is the reading of, and how many intervals it
1353    // could not read (then the true worst is unknown).
1354    let mut request_detail: Option<String> = None;
1355    let mut request_component = crate::MeasuredComponent::Wall;
1356    let mut request_unread = 0usize;
1357    if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1358        let wall_missing = wall_model_missing(&wall_readings, model);
1359        // The RAILS against the same request: the shipped support curves ARE
1360        // these rails (`chain_surgery` splits `cr`/`cs` into the support
1361        // edges), so a rail standing off the exact construction ships that
1362        // miss on its carrier. Read against the exact contacts the wall
1363        // reading already solved (j/32 of every interval); an interval whose
1364        // rail reads over `model`, or cannot be read, is short exactly as a
1365        // wall interval is. Acceptance (`bar`) is unchanged.
1366        let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1367            Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1368            None => (Vec::new(), Vec::new()),
1369        };
1370        // Per contact, keyed by the interval's left sample: distinct intervals.
1371        let rail_over: std::collections::BTreeSet<usize> =
1372            rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1373        let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1374        let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1375        let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1376        missing.sort_unstable();
1377        missing.dedup();
1378        if !missing.is_empty() {
1379            // The same declared reading, over the intervals short.
1380            let wall_short = wall_readings
1381                .iter()
1382                .filter(|reading| wall_missing.contains(&reading.0))
1383                .map(|&(_, declared, _, _)| declared)
1384                .fold(0.0_f64, f64::max);
1385            let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1386            let short = wall_short.max(rail_short);
1387            request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1388            request_unread = rail_unread_intervals.len();
1389            request_detail = Some(format!(
1390                "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1391                 wall {wall_short:.3e} over {} interval(s)",
1392                rail_over.len(),
1393                rail_unread_intervals.len(),
1394                wall_missing.len()
1395            ));
1396            let within_ceiling = (0..segments.len()).all(|segment| {
1397                let present = samples
1398                    .iter()
1399                    .filter(|sample| sample.segment == segment && (0..per_segment as isize).contains(&sample.position))
1400                    .count();
1401                present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1402                    <= CHAIN_CARVE_MAX_PER_SEGMENT
1403            });
1404            crate::blend::carve::carve_trace(format_args!(
1405                "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1406                 (worst reading + spread {short:.3e}), ceiling {}",
1407                missing.len(),
1408                if within_ceiling { "ok" } else { "reached" }
1409            ));
1410            unmet = if round >= model_rounds_limit() {
1411                Some((short, crate::BudgetReason::RoundsSpent))
1412            } else if !within_ceiling {
1413                Some((short, crate::BudgetReason::StationCeiling))
1414            } else {
1415                let mut refined = samples.clone();
1416                // The insertion is the next round's attempt, whatever it returns.
1417                CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1418                match insert_chain_stations(segments, radius, per_segment, &mut refined, &missing) {
1419                    Ok(()) => {
1420                        assign_chain_parameters(&mut refined, segments.len(), per_segment, false);
1421                        match build_closed_chain_from_samples(
1422                            solid, segments, radius, chamfer, name, per_segment, fold, bar, refined, round + 1,
1423                        ) {
1424                            Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1425                            _ => Some((short, crate::BudgetReason::Rejected)),
1426                        }
1427                    }
1428                    Err(_) => {
1429                        // The insertion itself was the attempt.
1430                        CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1431                        Some((short, crate::BudgetReason::Incoherent))
1432                    }
1433                }
1434            };
1435        }
1436    }
1437    // `BREP_DEBUG_CHAIN_DENSE_RAILS=k` (diagnostic only): the accepted rung's
1438    // rails against exact contacts at k-1 interior fractions of every interval,
1439    // where `rail_deviation` reads only the halfway one. Nothing reads it to
1440    // decide.
1441    if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1442        let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, per_segment, &samples, subdivisions.max(2))?;
1443        let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1444        let mut unreadable = 0usize;
1445        for (left, next, fraction, p1, p2) in &contacts.solved {
1446            let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1447            let global = (from + fraction * (to - from)).rem_euclid(1.0);
1448            let seed = (global - low) / range;
1449            let reach = (to - from).abs() / range;
1450            match rail_pair_miss(closest_on_rail(&cr, *p1, seed, reach).unwrap_or(f64::NAN), closest_on_rail(&cs, *p2, seed, reach).unwrap_or(f64::NAN)) {
1451                Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1452                Some(_) => {}
1453                None => unreadable += 1,
1454            }
1455        }
1456        eprintln!(
1457            "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1458             worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1459            contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1460        );
1461    }
1462    let crease = if fold == FoldPolicy::Carve {
1463        // THE BAR IS THE KERNEL'S OWN CONTRACT FOR THE QUANTITY MEASURED. What
1464        // the fit is held to is "does this coedge's curve-on-surface, pushed
1465        // back to 3-D, still trace the same locus as the edge's 3-D curve?",
1466        // which is `pcurve_consistency` — not `intersection_fit`, which is the
1467        // accuracy of a fitted intersection CURVE and is what the re-marched
1468        // wall itself is held to (`carve_ladder`). Measured on the
1469        // 2026-09-02 collar the crease's pcurves reach 5.4e-6 against that
1470        // 4e-3, so the margin is three decades and the bar is not what decides
1471        // the case either way.
1472        let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1473        let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1474            return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1475                "blend carve: the wall re-marched at {per_segment} stations per segment carries                  no fold to carve, but the march at {CHAIN_PER_SEGMENT} refused one"
1476            )));
1477        };
1478        Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1479    } else {
1480        None
1481    };
1482
1483    // Single-face sides (e.g. one cap around the whole rim) get ONE uv
1484    // fit in the rows' own normalized space — identical parameterization
1485    // to cr/cs, so the closed support piece and its pcurve agree exactly.
1486    let single_face = |side: usize| -> bool {
1487        let first_id = segments[0].mate(side).face.id;
1488        segments
1489            .iter()
1490            .all(|segment| segment.mate(side).face.id == first_id)
1491    };
1492    let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1493    for side in 0..2 {
1494        if !single_face(side) {
1495            continue;
1496        }
1497        let mut row = Vec::new();
1498        let select = |station: &Station| -> [f64; 2] {
1499            if side == 0 {
1500                station.uv1
1501            } else {
1502                station.uv2
1503            }
1504        };
1505        let mut push_uv = |station: &Station| {
1506            let uv = select(station);
1507            row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1508        };
1509        for offset in (1..=degree).rev() {
1510            push_uv(&global[count - offset].1.station);
1511        }
1512        for (_, sample) in &global {
1513            push_uv(&sample.station);
1514        }
1515        push_uv(&global[0].1.station);
1516        for offset in 1..=degree {
1517            push_uv(&global[offset].1.station);
1518        }
1519        whole_pcurves[side] = Some(fit_row(&row)?);
1520    }
1521
1522    // Per-face pcurve fits over each segment's full sample window
1523    // (in-segment + overshoot), in global parameters.
1524    let mut pcurves1 = Vec::with_capacity(segments.len());
1525    let mut pcurves2 = Vec::with_capacity(segments.len());
1526    for segment in 0..segments.len() {
1527        let mut window: Vec<&ChainSample> = samples
1528            .iter()
1529            .filter(|sample| sample.segment == segment)
1530            .collect();
1531        window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1532        let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1533        let low = params[0];
1534        let high = *params.last().unwrap();
1535        let normalized: Vec<f64> = params
1536            .iter()
1537            .map(|parameter| (parameter - low) / (high - low))
1538            .collect();
1539        let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1540            let points: Vec<Vec4> = window
1541                .iter()
1542                .map(|sample| {
1543                    let uv = select(&sample.station);
1544                    Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1545                })
1546                .collect();
1547            let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1548            // Re-express in GLOBAL parameters: the curve's [0,1] domain
1549            // corresponds to [low, high]; keep as-is and remember the
1550            // affine map through the window bounds.
1551            Ok(curve)
1552        };
1553        pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1554        pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1555    }
1556
1557    let wall_surface = unmet.as_ref().map(|_| surface.clone());
1558    let built = chain_surgery(
1559        solid,
1560        segments,
1561        ChainRows {
1562            surface,
1563            cr,
1564            cs,
1565            u_domain,
1566            fit_low: low,
1567            fit_range: range,
1568            pcurves1,
1569            pcurves2,
1570            whole_pcurves,
1571            crease,
1572        },
1573        name,
1574    )?;
1575    if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1576        crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1577            name: name.map(str::to_string),
1578            surface,
1579            request: model,
1580            residual,
1581            detail: request_detail,
1582            budget: crate::ApproximationBudget {
1583                reason,
1584                // Rounds ATTEMPTED (the deepest any refinement of this rung
1585                // reached); the stations and residual are this shipped rung's.
1586                rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1587                rounds_limit: model_rounds_limit(),
1588                stations: stations_now,
1589                station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1590                mechanism: crate::BudgetMechanism::LocalRounds,
1591                measured_component: request_component,
1592                unread: request_unread,
1593            },
1594        });
1595    }
1596    Ok(Rung::Built(built))
1597}
1598
1599pub(super) struct ChainRows {
1600    pub(super) surface: NurbsSurface,
1601    pub(super) cr: NurbsCurve,
1602    pub(super) cs: NurbsCurve,
1603    pub(super) u_domain: [f64; 2],
1604    /// Affine map from the rows' fit space to GLOBAL chain parameters:
1605    /// global = fit_low + p * fit_range (the global fit normalised its
1606    /// wrap-extended parameters onto [0, 1] before interpolating).
1607    pub(super) fit_low: f64,
1608    pub(super) fit_range: f64,
1609    /// Per-segment (window_low, window_high, uv curve over [0,1]) in
1610    /// GLOBAL parameters.
1611    pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1612    pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1613    /// Whole-chain uv fits (rows' fit space) for single-face sides.
1614    pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1615    /// The crease a CARVED wall is cut along: the inner loop of the blend
1616    /// face, one edge carrying two coedges of that same face.
1617    pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1618}
1619
1620impl ChainRows {
1621    pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1622        self.fit_low + fit_parameter * self.fit_range
1623    }
1624}
1625
1626/// Extract the pcurve portion for a global window [a, b] from a
1627/// segment's fitted uv curve, re-parameterised so its domain maps
1628/// affinely onto the portion (matching the split support edge).
1629pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1630    let (low, high, curve) = fitted;
1631    // The chain parameterisation wraps with period 1; pick the period
1632    // copy of the window that overlaps this segment's fit range.
1633    let to_local = |value: f64| {
1634        let mut best = value;
1635        for candidate in [value, value + 1.0, value - 1.0] {
1636            if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1637                best = candidate;
1638                break;
1639            }
1640        }
1641        ((best - low) / (high - low)).clamp(0.0, 1.0)
1642    };
1643    let mut local_a = to_local(a);
1644    let mut local_b = to_local(b);
1645    if local_a > local_b {
1646        std::mem::swap(&mut local_a, &mut local_b);
1647    }
1648    let epsilon = 1e-9;
1649    let (_, tail) = if local_a > epsilon {
1650        curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1651    } else {
1652        (curve.clone(), curve.clone())
1653    };
1654    let tail = if local_a > epsilon {
1655        tail
1656    } else {
1657        curve.clone()
1658    };
1659    let portion = if local_b < 1.0 - epsilon {
1660        tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1661    } else {
1662        tail
1663    };
1664    Ok(portion)
1665}
1666
1667/// The single boundary edge of a side's mate face(s) incident to a chain
1668/// vertex, EXCLUDING the two chain edges meeting there — i.e. the SEAM
1669/// (same face on both segments) or the SPOKE (face changes) that the
1670/// support curve must cross.  `None` when the chain simply flows through
1671/// (same face, no seam — e.g. a stadium cap rim vertex).
1672pub(super) fn cross_edge_at(
1673    solid: &BrepSolid,
1674    face_before: &FaceRecord,
1675    loop_before: usize,
1676    face_after: &FaceRecord,
1677    loop_after: usize,
1678    vertex: u64,
1679    before_edge: u64,
1680    after_edge: u64,
1681) -> Result<Option<u64>, KernelRefusal> {
1682    let mut found: Option<u64> = None;
1683    let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1684        for coedge in &face.loops[loop_index].coedges {
1685            if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1686                continue;
1687            }
1688            let edge = solid
1689                .edges
1690                .iter()
1691                .find(|edge| edge.id == coedge.edge_id)
1692                .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1693            if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1694                if found.is_some() && found != Some(edge.id) {
1695                    return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1696                }
1697                found = Some(edge.id);
1698            }
1699        }
1700        Ok(())
1701    };
1702    consider(face_before, loop_before)?;
1703    if !(face_after.id == face_before.id && loop_after == loop_before) {
1704        consider(face_after, loop_after)?;
1705    }
1706    Ok(found)
1707}
1708
1709/// The nearest periodic image of a seam boundary to `value` — the boundary a
1710/// crossing endpoint must be pinned to, read off its ADJACENT INTERIOR sample.
1711fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1712    let period = high - low;
1713    low + ((value - low) / period).round() * period
1714}
1715
1716/// The mate-face pcurve of a support PIECE, built by projecting the piece's
1717/// 3D curve onto the (analytic) mate surface and fitting.  Because a chain
1718/// support rides EXACTLY on its mate carrier, the projection is exact; the
1719/// periodic track is unwrapped so it never jumps a period across the carrier
1720/// seam, and crossing endpoints are pinned ONTO the seam so the piece and the
1721/// trimmed seam edge close the loop in parameter space.
1722///
1723/// The carrier is closed in EITHER direction, and the crossed seam is not
1724/// always the u meridian.  A cylinder/cone/sphere is closed in u only, so its
1725/// seam is the u0 ≡ u1 meridian — but a TORUS is biperiodic, and a tube welded
1726/// through a plane crosses the torus' v0 ≡ v1 PARALLEL (its outer equator)
1727/// instead.  Pinning u there drags the endpoint a half-carrier away from its
1728/// own 3D vertex (the collar loop then closes on a uv corner that is nowhere
1729/// near the vertex, and the mate face's trim is garbage), so the pinned AXIS is
1730/// chosen GEOMETRICALLY: whichever closed direction's seam actually passes
1731/// through the endpoint.  A carrier closed in u only can only ever take the u
1732/// branch, exactly as before.
1733pub(in crate::blend) fn project_piece_pcurve(
1734    piece: &NurbsCurve,
1735    surface: &NurbsSurface,
1736    start_is_crossing: bool,
1737    end_is_crossing: bool,
1738) -> Result<NurbsCurve, KernelRefusal> {
1739    project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1740}
1741
1742
1743/// A trim pcurve `old` as a GUIDE for [`project_piece_pcurve_from`]: at track
1744/// fraction `f` its uv at the same fraction of its own domain, closed
1745/// directions wrapped and open ones clamped into `surface`'s domain; `None`
1746/// where that uv is not finite.
1747pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
1748    let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
1749    let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1750    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1751    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1752    let into = move |value: f64, closed: bool, low: f64, high: f64| {
1753        if closed && high > low {
1754            low + (value - low).rem_euclid(high - low)
1755        } else {
1756            value.clamp(low, high)
1757        }
1758    };
1759    Ok(move |fraction: f64| {
1760        let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
1761        let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
1762        seed.iter().all(|value| value.is_finite()).then_some(seed)
1763    })
1764}
1765
1766/// [`project_piece_pcurve`], every sample's foot sought from `guide` at its
1767/// fraction (a uv on the branch the piece is known to run on — a support
1768/// trim's own march pcurve) instead of from its neighbour's foot. A foot
1769/// continued from its neighbour stays wherever that neighbour's did: at a
1770/// STATIONARY edge of the carrier (an extrusion of a profile whose derivative
1771/// vanishes at its start, S_u = 0 along u = 0) a foot seeded on that edge has
1772/// no gradient to leave it by, and every later sample stayed there with it
1773/// (the stationary-start prism's refit: one foot, [0, 0.9], for the whole
1774/// rail, 6.0 off it at the far end). Without a guide, as before.
1775pub(in crate::blend) fn project_piece_pcurve_from(
1776    piece: &NurbsCurve,
1777    surface: &NurbsSurface,
1778    start_is_crossing: bool,
1779    end_is_crossing: bool,
1780    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1781) -> Result<NurbsCurve, KernelRefusal> {
1782    project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
1783}
1784
1785/// [`project_piece_pcurve_from`], its track first refined (two bisections)
1786/// in the dense interval holding each fraction of `refine_near` — where a
1787/// caller that read the SHIPPED curve again found it over the floor between
1788/// the fitter's own stencil points — and then fitted at the same floor.
1789pub(in crate::blend) fn project_piece_pcurve_within(
1790    piece: &NurbsCurve,
1791    surface: &NurbsSurface,
1792    start_is_crossing: bool,
1793    end_is_crossing: bool,
1794    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
1795    refine_near: &[f64],
1796) -> Result<NurbsCurve, KernelRefusal> {
1797    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1798    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1799    let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1800    let period = u1 - u0;
1801    let period_v = v1 - v0;
1802    let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
1803    // Dense reference projection; the pcurve is fit from an adaptively
1804    // thinned subset kept as coarse as tolerance allows.
1805    //
1806    // The track is as dense as the PIECE is, not a fixed 96. A carved wall's
1807    // re-marched rail carries hundreds of control points (574 on the
1808    // 2026-09-02 two-torus collar), and a pcurve fitted and checked at 97
1809    // samples of it was accepted 1.19e-3 off the rail it claims.
1810    let dense = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
1811    let mut proj_u = Vec::with_capacity(dense + 1);
1812    let mut proj_v = Vec::with_capacity(dense + 1);
1813    let mut point3 = Vec::with_capacity(dense + 1);
1814    let mut params = Vec::with_capacity(dense + 1);
1815    let mut previous_u: Option<f64> = None;
1816    // The foot each sample continues from: the track is ONE branch of the
1817    // carrier, so after the first sample each inversion starts from its
1818    // neighbour's foot rather than asking for the nearest point of the whole
1819    // surface, which on a carrier that comes back near itself is another sheet.
1820    let mut previous_foot: Option<[f64; 2]> = None;
1821    let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
1822    for section in 0..=dense {
1823        let t = d0 + (d1 - d0) * section as f64 / dense as f64;
1824        let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
1825        let guide = guide.filter(|_| {
1826            true
1827        });
1828        // The guide's seed and the neighbour's foot are both tried and the
1829        // NEARER foot kept: the seeded projector is local and returns a
1830        // stalled Newton's iterate without saying so, so neither seed is
1831        // trusted alone (a foot on its carrier is nearer than any stall).
1832        let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(section as f64 / dense as f64)).into_iter().chain(previous_foot).collect();
1833        let mut projection: Option<crate::projection::SurfaceProjection> = None;
1834        for seed in &seeds {
1835            let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
1836            if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
1837                projection = Some(found);
1838            }
1839        }
1840        let projection = match projection {
1841            Some(projection) => projection,
1842            None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
1843        };
1844        previous_foot = Some([projection.u, projection.v]);
1845        feet.push([projection.u, projection.v]);
1846        let mut u = projection.u;
1847        if let Some(previous) = previous_u {
1848            while u - previous > period * 0.5 {
1849                u -= period;
1850            }
1851            while previous - u > period * 0.5 {
1852                u += period;
1853            }
1854        }
1855        previous_u = Some(u);
1856        proj_u.push(u);
1857        proj_v.push(projection.v);
1858        point3.push(point);
1859        params.push(section as f64 / dense as f64);
1860    }
1861    // A biperiodic carrier needs the SAME unwrap in v, or a track that walks up
1862    // to the v seam comes back as 0 at the far end and the fit swings across
1863    // the whole carrier.  Anchor the chain on the first INTERIOR sample and
1864    // pull index 0 onto it afterwards: an endpoint sitting exactly ON the seam
1865    // inverts to either boundary at the solver's whim, and letting that
1866    // coin flip anchor the chain would shift the whole piece a period.
1867    if closed_v {
1868        for index in 2..=dense {
1869            while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
1870                proj_v[index] -= period_v;
1871            }
1872            while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
1873                proj_v[index] += period_v;
1874            }
1875        }
1876        while proj_v[0] - proj_v[1] > period_v * 0.5 {
1877            proj_v[0] -= period_v;
1878        }
1879        while proj_v[1] - proj_v[0] > period_v * 0.5 {
1880            proj_v[0] += period_v;
1881        }
1882    }
1883    // Interior mean decides which meridian a crossing endpoint sits on, and
1884    // whether the whole track needs a full-period shift into the domain.
1885    let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
1886    let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
1887    // Pin a crossing endpoint onto the seam it actually crosses.  Each
1888    // candidate keeps the OTHER coordinate as projected, so the losing axis
1889    // moves the endpoint bodily off its own 3D vertex — the residual names the
1890    // winner with no surface-type special-casing.  u keeps the whole-track
1891    // mean; a v run can legitimately cover a full period, so its boundary is
1892    // read off the ADJACENT INTERIOR sample instead.  Both pins land in the
1893    // track's own (pre-shift) period, so the shift below carries them along
1894    // with the rest of the track.
1895    let pinned = |index: usize,
1896                  neighbour: usize,
1897                  proj_u: &[f64],
1898                  proj_v: &[f64]|
1899     -> Result<(f64, f64), KernelRefusal> {
1900        let point = point3[index];
1901        let u_error = surface
1902            .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
1903            .sub(point)
1904            .length();
1905        if closed_v {
1906            let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
1907            let v_error = surface
1908                .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
1909                .sub(point)
1910                .length();
1911            if v_error < u_error {
1912                return Ok((proj_u[index], seam_v));
1913            }
1914        }
1915        Ok((seam_u, proj_v[index]))
1916    };
1917    if start_is_crossing {
1918        let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
1919        proj_u[0] = u;
1920        proj_v[0] = v;
1921    }
1922    if end_is_crossing {
1923        let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
1924        proj_u[dense] = u;
1925        proj_v[dense] = v;
1926    }
1927    let mut shift = 0.0;
1928    while mean_u + shift > u1 + 1e-9 {
1929        shift -= period;
1930    }
1931    while mean_u + shift < u0 - 1e-9 {
1932        shift += period;
1933    }
1934    for u in proj_u.iter_mut() {
1935        *u += shift;
1936    }
1937    if closed_v {
1938        let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
1939        let mut shift_v = 0.0;
1940        while mean_v + shift_v > v1 + 1e-9 {
1941            shift_v -= period_v;
1942        }
1943        while mean_v + shift_v < v0 - 1e-9 {
1944            shift_v += period_v;
1945        }
1946        for v in proj_v.iter_mut() {
1947            *v += shift_v;
1948        }
1949    }
1950    // Coarsest interpolation whose fitted pcurve reproduces the projected TRACK
1951    // to the kernel's pcurve refinement floor — the floor every other pcurve fit
1952    // is held to, and the one the shell vector-area bar is built from.
1953    //
1954    // It used to stop at 0.002 absolute ("half the validator's tolerance"),
1955    // which is `pcurve_consistency / 2` on any model under 20 units and so four
1956    // decades looser than any other trim: a torus face carrying a collar's rail
1957    // was accepted 1.19e-3 off it, and that one coedge was the whole of a
1958    // 1.94e-4 shell vector-area residual, 14.5x its face's bar.
1959    //
1960    // The deviation is measured against the projected track, not the piece's
1961    // 3D points: a pcurve lies on its surface by construction, so a piece that
1962    // is itself off the carrier would hold the ladder at its top however fine
1963    // the fit, and that is a RAIL defect for the march to answer, not this fit.
1964    let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
1965    // The ladder is the one verified fitter (`track_fit`), which checks every
1966    // rung at the MIDPOINTS between dense samples, projected on their own —
1967    // checked at the dense samples, the top rung, which interpolates every one
1968    // of them, reads zero by construction and passes whatever lies between (on
1969    // the 574-control-point collar rail it read 9.6e-15) — and breaks the fit,
1970    // C1, where the track crosses a knot line of the carrier.
1971    let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
1972    let mut track = crate::blend::track_fit::Track {
1973        fractions: params,
1974        uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
1975        feet,
1976        breaks: Vec::new(),
1977    };
1978    let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
1979    crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
1980    crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
1981    for &fraction in refine_near {
1982        for _ in 0..2 {
1983            let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
1984            if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], PIECE_TRACK_CEILING)?.is_none() {
1985                break;
1986            }
1987        }
1988    }
1989    // Where the top rung still misses the floor, the track is refined ONLY in
1990    // the dense intervals whose out-of-sample checks missed (local insertion
1991    // concentrates the rail's stations, and a uniform dense track can leave one
1992    // sample per knot span where the rail turns hardest), within the dense
1993    // ceiling the track was built under; every other sample, the pinned ends
1994    // and the breaks are kept. The floor and the refusal below are unchanged.
1995    let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
1996        surface,
1997        &mut track,
1998        &at,
1999        tolerance,
2000        REFINE_ROUNDS,
2001        PIECE_TRACK_CEILING,
2002    )?;
2003    if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2004        eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2005        eprintln!(
2006            "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2007            fit.miss,
2008            fit.samples,
2009            track.fractions.len(),
2010            track.breaks.len(),
2011            fit.on_floor
2012        );
2013        // Where the miss sits, against what the piece is made of: an endpoint
2014        // pinned onto a seam reads at fraction 0 or 1; an interior miss names
2015        // the piece's knot span there and how that span compares with the
2016        // piece's mean (local insertion clusters spans at the crotch).
2017        if let Some(worst) = fit.worst {
2018            let parameter = d0 + (d1 - d0) * worst.fraction;
2019            let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2020                .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2021                .map(|pair| [pair[0], pair[1]]).collect();
2022            let mean = (d1 - d0) / spans.len().max(1) as f64;
2023            let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2024            let width = span.map(|index| spans[index][1] - spans[index][0]);
2025            eprintln!(
2026                "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2027                 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2028                worst.fraction, parameter, span, spans.len(), width, mean,
2029                worst.bracket[0], worst.bracket[1], worst.interval
2030            );
2031        }
2032    }
2033    // A fit off the floor is refused, not handed on as the best rung. Before the
2034    // fitter broke at knot lines this ladder returned its best rung silently,
2035    // and on the 2026-09-02 two-torus collar that was 1.07e-7 to 8.46e-7 off at
2036    // the top rung on every collar that builds (measured 2026-09-17); broken,
2037    // each of them reaches the floor. What still misses is a rail the march tore
2038    // (2.7 to 3.3 on the oversized radii, a body refused or rejected anyway).
2039    if !fit.on_floor {
2040        return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_pcurve_floor", format!(
2041            "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2042             against a floor of {tolerance:.1e}",
2043            crate::blend::PCURVE_OFF_FLOOR,
2044            fit.miss,
2045            fit.samples
2046        )));
2047    }
2048    Ok(fit.curve)
2049}