BREP_kernel 0.11.0

A boundary representation (BREP) geometry kernel for building CAD applications.
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
//! Bounded reconstruction from stored carrier geometry and planar trim images.
//! Exact analytic, isoparametric and rational ruled-plane sections supplement
//! planar images when the inherited image itself misses the paired carrier.
//! Candidates must reach settled vertices, stay on both carriers, and pass
//! movement and nonincident face clearance checks. No station fit is used.

use crate::{BrepSolid, KernelRefusal, KernelStage, KernelTolerances, NurbsCurve, NurbsSurface, OrRefuse, Vec3, Vec4};
use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};

/// An edge this pass replaced, with the curve it replaced.
pub(in crate::boolean) struct ReplacedEdge {
    pub edge_id: u64,
    pub previous: NurbsCurve,
    pub t0: f64,
    pub t1: f64,
    /// Worst sampled distance to either carrier, before and after.
    pub residual_before: f64,
    pub residual_after: f64,
}

/// Put back every edge in `replaced`.
pub(in crate::boolean) fn restore_edges(solid: &mut BrepSolid, replaced: Vec<ReplacedEdge>) {
    for item in replaced {
        if let Some(edge) = solid.edges.iter_mut().find(|edge| edge.id == item.edge_id) {
            edge.curve = item.previous;
            edge.t0 = item.t0;
            edge.t1 = item.t1;
        }
    }
}

const STATIONS: usize = 33;

thread_local! {
    static RESIDUAL_QUERIES: std::cell::Cell<u64> = const { std::cell::Cell::new(0) };
    static RESIDUAL_WORK: std::cell::Cell<(u64, u64)> = const { std::cell::Cell::new((0, 0)) };
}


/// Worst distance from `curve` over `[t0, t1]` to the carriers in `surfaces`.
fn worst_residual(curve: &NurbsCurve, t0: f64, t1: f64, surfaces: &[&NurbsSurface]) -> Result<f64, KernelRefusal> {
    RESIDUAL_QUERIES.with(|queries| queries.set(queries.get() + 1));
    // Union folds repeatedly revisit edges beside a touched face, even when
    // their geometry and both carriers are unchanged. Cache the reading,
    // never the repair decision: policy bars and candidate checks still run.
    let key = crate::mass_properties::carrier_residual_identity(curve, t0, t1, surfaces);
    crate::reading_memo::memoized(crate::reading_memo::Reading::CarrierResidual, key,
        || worst_residual_uncached(curve, t0, t1, surfaces).ok())
        .map(Ok).unwrap_or_else(|| worst_residual_uncached(curve, t0, t1, surfaces))
}

fn worst_residual_uncached(curve: &NurbsCurve, t0: f64, t1: f64, surfaces: &[&NurbsSurface]) -> Result<f64, KernelRefusal> {
    let stations = carrier_stations(curve, t0, t1);
    RESIDUAL_WORK.with(|work| {
        let (calls, projections) = work.get();
        work.set((calls + 1, projections + (stations.len() * surfaces.len()) as u64));
    });
    let mut worst = 0.0f64;
    for t in stations {
        let point = curve.evaluate(t).or_refuse(KernelStage::Sew, "evaluate")?;
        for surface in surfaces {
            let projection = crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Sew, "project_point_to_surface")?;
            worst = worst.max(projection.distance);
        }
    }
    Ok(worst)
}

/// Uniform stations alone alias inherited polygon vertices. Include each
/// represented knot span's quarters in both entry and candidate checks.
pub(super) fn carrier_stations(curve: &NurbsCurve, t0: f64, t1: f64) -> Vec<f64> {
    let mut breaks = vec![t0, t1];
    breaks.extend(curve.knots.iter().copied().filter(|&t| t > t0 && t < t1));
    breaks.sort_by(f64::total_cmp);
    breaks.dedup();
    let mut stations: Vec<_> = (0..STATIONS)
        .map(|i| t0 + (t1 - t0) * i as f64 / (STATIONS - 1) as f64).collect();
    for span in breaks.windows(2) {
        for i in 0..=4 { stations.push(span[0] + (span[1] - span[0]) * i as f64 / 4.0); }
    }
    stations.sort_by(f64::total_cmp);
    stations.dedup();
    stations
}

/// Restrict the nearest-point witness to the represented edge interval.
/// A parent curve can contain the complementary branch of a section; using
/// that parent as the displacement witness would make a branch swap read zero.
pub(super) fn represented_curve(edge: &crate::topology::EdgeRecord) -> Option<NurbsCurve> {
    super::image::restricted_curve(&edge.curve,edge.t0,edge.t1)
}

/// Exact carrier constructions, without station interpolation. A constant
/// trim coordinate supports a represented isocurve; recognized analytic
/// surface pairs supply their existing exact sections. Unsupported pairs stay
/// unchanged. Restriction must reach both settled vertices: never bend these
/// curves by pinning their controls.
fn exact_carrier_candidates(
    surfaces: [&NurbsSurface; 2],
    uses: &[(usize, NurbsCurve, bool)],
    faces: &[&crate::FaceRecord],
    start: Vec3,
    end: Vec3,
    bar: f64,
) -> Vec<NurbsCurve> {
    let mut full = crate::intersect_analytic_pair(surfaces[0], surfaces[1], bar).unwrap_or_default();
    if let Ok(Some(curve))=crate::imprint::planar_carrier_section(surfaces[0],surfaces[1],bar) {full.push(curve);}
    for (plane,ruled) in [(surfaces[0],surfaces[1]),(surfaces[1],surfaces[0])] {
        if let Some(pair)=crate::imprint::plane_section_graph(plane,ruled).filter(|pair|plane_supports_curve(plane,&pair.curve,bar*0.1)) {
            // The algebraic construction also supplies the carrier chart.
            // Check its endpoint images before handing the 3D branch to the
            // existing restriction, displacement and clearance guards.
            if let Ok([a,b])=pair.pcurve.domain() {
                if [a,b].into_iter().all(|t| pair.pcurve.evaluate(t).ok().and_then(|uv|ruled.evaluate(uv.x,uv.y).ok())
                    .zip(pair.curve.evaluate(t).ok()).is_some_and(|(p,q)|p.sub(q).length()<=bar*0.1)) {full.push(pair.curve);}
            }
        }
        if let (Ok(a),Ok(b))=(crate::project_point_to_surface(ruled,start),crate::project_point_to_surface(ruled,end)) {
            if let Some(pair)=crate::imprint::plane_section_graph_range(plane,ruled,[a.u.min(b.u),a.u.max(b.u)]).filter(|pair|plane_supports_curve(plane,&pair.curve,bar*0.1)) {full.push(pair.curve);}
        }
    }
    for (face_index, trim, _) in uses {
        let surface = &faces[*face_index].surface;
        if trim.control_points.is_empty() || trim.control_points.iter().any(|p| !(p.w > 0.0)) { continue; }
        let first = trim.control_points[0];
        let u = first.x / first.w;
        let v = first.y / first.w;
        // Exact constant coordinates only. A UV-near trim is not proof of
        // physical boundary support on a large or reparameterized carrier.
        if trim.control_points.iter().all(|p| p.x / p.w == u) {
            if let Ok(curve) = surface.iso_curve_u(u) { full.push(curve); }
        }
        if trim.control_points.iter().all(|p| p.y / p.w == v) {
            if let Ok(curve) = surface.iso_curve_v(v) { full.push(curve); }
        }
    }
    full.into_iter().filter_map(|curve| {
        let a = crate::project_point_to_curve(&curve, start).ok()?;
        let b = crate::project_point_to_curve(&curve, end).ok()?;
        if a.distance > bar || b.distance > bar { return None; }
        if start.sub(end).length() <= bar {
            let [d0,d1] = curve.domain().ok()?;
            if curve.evaluate(d0).ok()?.sub(start).length() <= bar && curve.evaluate(d1).ok()?.sub(end).length() <= bar {
                return Some(curve);
            }
            return None;
        }
        if (a.u - b.u).abs() < 1e-12 { return None; }
        let [d0, d1] = curve.domain().ok()?;
        let lo = a.u.min(b.u);
        let hi = a.u.max(b.u);
        let mut part = curve;
        if hi < d1 { part = part.split(hi).ok()?.0; }
        if lo > d0 { part = part.split(lo).ok()?.1; }
        if a.u > b.u { part = part.reversed().ok()?; }
        Some(part)
    }).collect()
}

/// The exact rational image through a supported bilinear patch. Approximate
/// affine recognition is only a scope gate: it never removes a warp term.
pub(super) fn affine_image(surface: &NurbsSurface, pcurve: &NurbsCurve) -> Result<Option<NurbsCurve>, KernelRefusal> {
    if !surface.is_affine().or_refuse(KernelStage::Sew,"is_affine")? {return Ok(None);}
    Ok(super::image::bilinear_image(surface,pcurve))
}

/// A linear chart is a proposal, not proof of an approximately affine
/// surface. Compose it through the complete bilinear net and bound its
/// same-parameter difference from the proposed section over every span.
pub(super) fn plane_supports_curve(surface:&NurbsSurface,curve:&NurbsCurve,bar:f64)->bool {
    if !surface.is_affine().unwrap_or(false) {return false;}
    let Ok([a,b])=curve.domain() else {return false;};
    let Ok(fit)=crate::fit_pcurve_on_surface_range(surface,curve,a,b,true,bar) else {return false;};
    // Range fitting uses a normalized trim parameter. Restore the known
    // affine correspondence before a same-parameter whole-image comparison.
    let Ok([x,y])=fit.curve.domain() else {return false;};
    let mut knots=fit.curve.knots.clone();
    for t in &mut knots { *t=if *t==x {a} else if *t==y {b} else {a+(b-a)*(*t-x)/(y-x)}; }
    let Ok(pc)=NurbsCurve::new(fit.curve.degree,knots,fit.curve.control_points) else {return false;};
    let Some(image)=super::image::bilinear_image(surface,&pc) else {return false;};
    super::image::difference_bound(curve,&image).is_some_and(|d|d<=bar)
}

/// Move a curve's two end control points onto `start` and `end`.
fn pin_ends(curve: &NurbsCurve, start: Vec3, end: Vec3) -> Result<NurbsCurve, KernelRefusal> {
    let mut controls = curve.control_points.clone();
    let last = controls.len() - 1;
    controls[0] = Vec4::from_point(start, controls[0].w);
    controls[last] = Vec4::from_point(end, controls[last].w);
    NurbsCurve::new(curve.degree, curve.knots.clone(), controls).or_refuse(KernelStage::Sew, "NurbsCurve::new")
}

/// A one-span rational quadratic — the form every circular, elliptic and
/// other conic arc trim takes — read as the CONIC it lies on, so an arc can
/// be re-cut to new ends without leaving its circle. Moving one end control
/// point of such an arc (`pin_ends`) bends it (3.07e-5 off the offset ball
/// for `BadBoolean`'s arc 40, record §9); re-deriving the Bézier form on the
/// same conic from the new ends does not.
pub(super) struct Conic {
    h: [Vec4; 3],
    k0: f64,
    k1: f64,
}

impl Conic {
    pub(super) fn of(curve: &NurbsCurve) -> Option<Conic> {
        if curve.degree != 2 || curve.control_points.len() != 3 || curve.knots.len() != 6 {
            return None;
        }
        let k = &curve.knots;
        if (k[0] - k[2]).abs() > 0.0 || (k[3] - k[5]).abs() > 0.0 || !(k[3] > k[0]) {
            return None;
        }
        Some(Conic { h: [curve.control_points[0], curve.control_points[1], curve.control_points[2]], k0: k[0], k1: k[3] })
    }

    /// Homogeneous point and its derivative at Bézier parameter `u`, for ANY
    /// real `u` (the conic continues past the arc's ends).
    fn homogeneous(&self, u: f64) -> (Vec4, Vec4) {
        let b = [(1.0 - u) * (1.0 - u), 2.0 * u * (1.0 - u), u * u];
        let db = [-2.0 * (1.0 - u), 2.0 - 4.0 * u, 2.0 * u];
        let mut point = Vec4 { x: 0.0, y: 0.0, z: 0.0, w: 0.0 };
        let mut derivative = Vec4 { x: 0.0, y: 0.0, z: 0.0, w: 0.0 };
        for index in 0..3 {
            let h = self.h[index];
            point.x += b[index] * h.x;
            point.y += b[index] * h.y;
            point.z += b[index] * h.z;
            point.w += b[index] * h.w;
            derivative.x += db[index] * h.x;
            derivative.y += db[index] * h.y;
            derivative.z += db[index] * h.z;
            derivative.w += db[index] * h.w;
        }
        (point, derivative)
    }

    fn point(&self, u: f64) -> Option<Vec3> {
        let (p, _) = self.homogeneous(u);
        (p.w.abs() > 1e-300).then(|| Vec3::new(p.x / p.w, p.y / p.w, p.z / p.w))
    }

    fn tangent(&self, u: f64) -> Option<Vec3> {
        let (p, d) = self.homogeneous(u);
        if p.w.abs() <= 1e-300 {
            return None;
        }
        let w2 = p.w * p.w;
        Some(Vec3::new((d.x * p.w - p.x * d.w) / w2, (d.y * p.w - p.y * d.w) / w2, (d.z * p.w - p.z * d.w) / w2))
    }

    /// The parameter of the conic point nearest `target`, from `seed`.
    pub(super) fn nearest(&self, target: Vec3, seed: f64) -> Option<f64> {
        let mut u = seed;
        for _ in 0..50 {
            let point = self.point(u)?;
            let tangent = self.tangent(u)?;
            let speed2 = tangent.dot(tangent);
            if !(speed2 > 1e-300) {
                return None;
            }
            let step = point.sub(target).dot(tangent) / speed2;
            u -= step;
            if !u.is_finite() || u.abs() > 1e3 {
                return None;
            }
            if step.abs() < 1e-16 {
                break;
            }
        }
        Some(u)
    }

    /// The same conic as a one-span rational quadratic from `ua` to `ub`,
    /// with the old parameter values as its knots.
    pub(super) fn arc(&self, ua: f64, ub: f64) -> Option<NurbsCurve> {
        if !(ub > ua) {
            return None;
        }
        let a = self.point(ua)?;
        let b = self.point(ub)?;
        let ta = self.tangent(ua)?;
        let tb = self.tangent(ub)?;
        // Middle control point: where the end tangents meet (least squares
        // in 3D; the conic is planar so the lines meet to rounding).
        let d = a.sub(b);
        let aa = ta.dot(ta);
        let bb = tb.dot(tb);
        let ab = ta.dot(tb);
        let det = aa * bb - ab * ab;
        if !(det.abs() > 1e-14 * aa * bb) {
            return None;
        }
        let da = d.dot(ta);
        let db = d.dot(tb);
        // a + s ta = b + r tb  ->  [aa -ab; ab -bb] [s r]^T = [-da; -db]
        let s = (-da * bb + ab * db) / det;
        let m = a.add(ta.scale(s));
        // Weight from a third conic point: its barycentric position in the
        // triangle A M B fixes both the Bézier parameter and the weight.
        let q = self.point(0.5 * (ua + ub))?;
        let e1 = m.sub(a);
        let e2 = b.sub(a);
        let n = e1.cross(e2);
        let n2 = n.dot(n);
        if !(n2 > 1e-300) {
            return None;
        }
        let qa = q.sub(a);
        let beta = qa.cross(e2).dot(n) / n2; // coordinate on M
        let gamma = e1.cross(qa).dot(n) / n2; // coordinate on B
        let alpha = 1.0 - beta - gamma;
        if !(alpha > 0.0 && beta > 0.0 && gamma > 0.0) {
            return None;
        }
        let ratio = (gamma / alpha).sqrt();
        let u = ratio / (1.0 + ratio);
        let w = beta * (1.0 - u) / (2.0 * u * alpha);
        if !(w.is_finite() && w > 0.0) {
            return None;
        }
        let ka = self.k0 + ua * (self.k1 - self.k0);
        let kb = self.k0 + ub * (self.k1 - self.k0);
        NurbsCurve::new(2, vec![ka, ka, ka, kb, kb, kb], vec![Vec4::from_point(a, 1.0), Vec4::from_point(m, w), Vec4::from_point(b, 1.0)]).ok()
    }
}

/// The candidate for one planar use: the trim's image re-cut to the edge's
/// vertices — on its own conic for a rational arc, by restriction for any
/// other curve whose image reaches both vertices — or, failing both, the
/// image with its ends pinned (which a later bound then judges).
fn candidate_from_image(image: &NurbsCurve, start: Vec3, end: Vec3, bar: f64) -> Result<(NurbsCurve, &'static str), KernelRefusal> {
    if let Some(conic) = Conic::of(image) {
        let (Some(ua), Some(ub)) = (conic.nearest(start, 0.0), conic.nearest(end, 1.0)) else {
            return Ok((pin_ends(image, start, end)?, "pinned"));
        };
        let on_conic = conic.point(ua).map_or(f64::INFINITY, |p| p.sub(start).length()).max(conic.point(ub).map_or(f64::INFINITY, |p| p.sub(end).length()));
        if on_conic <= bar {
            if let Some(arc) = conic.arc(ua, ub) {
                return Ok((pin_ends(&arc, start, end)?, "conic"));
            }
        }
        return Ok((pin_ends(image, start, end)?, "pinned"));
    }
    let [d0, d1] = image.domain().or_refuse(KernelStage::Sew, "domain")?;
    let epsilon = ((d1 - d0).abs().max(1.0) * 1e-10).max(2e-9);
    let pa = crate::project_point_to_curve(image, start).or_refuse(KernelStage::Sew, "project_point_to_curve")?;
    let pb = crate::project_point_to_curve(image, end).or_refuse(KernelStage::Sew, "project_point_to_curve")?;
    if pa.distance <= bar && pb.distance <= bar && pb.u > pa.u + 2.0 * epsilon {
        let mut restricted = image.clone();
        if pb.u < d1 - epsilon {
            restricted = restricted.split(pb.u).or_refuse(KernelStage::Sew, "split")?.0;
        }
        if pa.u > d0 + epsilon {
            restricted = restricted.split(pa.u).or_refuse(KernelStage::Sew, "split")?.1;
        }
        return Ok((pin_ends(&restricted, start, end)?, "restricted"));
    }
    Ok((pin_ends(image, start, end)?, "pinned"))
}

/// Reconstruct every two-use edge of `solid` that stands off one of its
/// carriers by more than the joint bar, from a planar use's trim. Returns the
/// edges replaced, with their previous curves.
pub(in crate::boolean) fn reconstruct_edges_from_carriers(solid: &mut BrepSolid, tolerance: f64) -> Result<Vec<ReplacedEdge>, KernelRefusal> {
    reconstruct_edges_from_carriers_in(solid, tolerance, None)
}

/// [`reconstruct_edges_from_carriers`] over the edges with a use on a face
/// in `scope` (`None` reads every edge): a boolean passes the faces it
/// touched; an edge between two untouched faces stands where the operand
/// put it.
pub(in crate::boolean) fn reconstruct_edges_from_carriers_in(
    solid: &mut BrepSolid,
    tolerance: f64,
    scope: Option<&HashSet<u64>>,
) -> Result<Vec<ReplacedEdge>, KernelRefusal> {
    let queries_before = RESIDUAL_QUERIES.with(|queries| queries.get());
    let work_before = RESIDUAL_WORK.with(|work| work.get());
    let policy = KernelTolerances::for_solid(solid, 1e-7);
    let bar = policy.model.min(tolerance.max(1e-12));
    let carrier_band = policy.pcurve_consistency;
    let debug = std::env::var("BREP_DEBUG_BOOL").is_ok() || std::env::var("BREP_DEBUG_JOINTS").is_ok();
    if std::env::var("BREP_RECONSTRUCT_EDGES").as_deref() == Ok("0") {
        // Diagnostic A/B hatch, not a policy (the record's before column).
        return Ok(Vec::new());
    }

    let faces: Vec<&crate::FaceRecord> = solid.shells.iter().flat_map(|shell| &shell.faces).collect();
    // edge id -> (face index, pcurve, forward) for every use.
    let mut uses: HashMap<u64, Vec<(usize, NurbsCurve, bool)>> = HashMap::default();
    for (face_index, face) in faces.iter().enumerate() {
        for loop_record in &face.loops {
            for coedge in &loop_record.coedges {
                uses.entry(coedge.edge_id).or_default().push((face_index, coedge.pcurve.clone(), coedge.forward));
            }
        }
    }
    let points: HashMap<u64, Vec3> = solid.vertices.iter().map(|vertex| (vertex.id, vertex.point)).collect();

    let mut candidates = Vec::<(usize, NurbsCurve, f64, f64)>::new();
    let mut declined = Vec::<String>::new();
    let mut off_edges = 0usize;
    for (edge_index, edge) in solid.edges.iter().enumerate() {
        if edge.degenerate || !(edge.t1 > edge.t0) {
            continue;
        }
        let Some(edge_uses) = uses.get(&edge.id) else { continue };
        if edge_uses.len() != 2 || edge_uses[0].0 == edge_uses[1].0 {
            continue;
        }
        if scope.is_some_and(|scope| edge_uses.iter().all(|(face_index, _, _)| !scope.contains(&faces[*face_index].id))) {
            continue;
        }
        let surfaces = [&faces[edge_uses[0].0].surface, &faces[edge_uses[1].0].surface];
        let before = worst_residual(&edge.curve, edge.t0, edge.t1, &surfaces)?;
        if before <= bar {
            continue;
        }
        off_edges += 1;
        if before > carrier_band {
            declined.push(format!("edge {} stands {before:.3e} off a carrier, over the {carrier_band:.3e} contract", edge.id));
            continue;
        }
        let (Some(&start), Some(&end)) = (points.get(&edge.start_vertex_id), points.get(&edge.end_vertex_id)) else { continue };
        let mut best: Option<(NurbsCurve, f64, String)> = None;
        let mut planar_uses = 0usize;
        for (face_index, pcurve, forward) in edge_uses {
            let Some(image) = affine_image(&faces[*face_index].surface, pcurve)? else { continue };
            planar_uses += 1;
            let image = if *forward { image } else { image.reversed().or_refuse(KernelStage::Sew, "reversed")? };
            let [i0, i1] = image.domain().or_refuse(KernelStage::Sew, "domain")?;
            let image_start = image.evaluate(i0).or_refuse(KernelStage::Sew, "evaluate")?;
            let image_end = image.evaluate(i1).or_refuse(KernelStage::Sew, "evaluate")?;
            let end_miss = image_start.sub(start).length().max(image_end.sub(end).length());
            if end_miss > carrier_band {
                declined.push(format!("edge {}: face {}'s trim image ends {end_miss:.3e} from the edge's vertices", edge.id, faces[*face_index].id));
                continue;
            }
            let (pinned, how) = candidate_from_image(&image, start, end, bar)?;
            let [i0, i1] = pinned.domain().or_refuse(KernelStage::Sew, "domain")?;
            let after = worst_residual(&pinned, i0, i1, &surfaces)?;
            // Displacement from the old edge: the distance from each sample
            // of the candidate to the old curve. Geometric, not at equal
            // parameter fractions — a re-cut conic carries the trim's
            // parameterisation, not the old edge's, and equal fractions read
            // that as a 0.29 move on an arc that lies on both carriers to 5e-12.
            let mut moved = 0.0f64;
            for step in 0..STATIONS {
                let fraction = step as f64 / (STATIONS - 1) as f64;
                let new_point = pinned.evaluate(i0 + (i1 - i0) * fraction).or_refuse(KernelStage::Sew, "evaluate")?;
                let nearest = crate::project_point_to_curve(&edge.curve, new_point).or_refuse(KernelStage::Sew, "project_point_to_curve")?;
                moved = moved.max(nearest.distance);
            }
            let note = format!("face {} trim image ({how}): on carriers to {after:.3e}, ends from {end_miss:.3e}, moved {moved:.3e}", faces[*face_index].id);
            if after > bar || moved > carrier_band {
                declined.push(format!("edge {} ({before:.3e} off): {note}", edge.id));
                continue;
            }
            if best.as_ref().map_or(true, |(_, best_after, _)| after < *best_after) {
                best = Some((pinned, after, note));
            }
        }
        if best.is_none() {
            let Some(previous)=represented_curve(edge) else {continue;};
            for mut candidate in exact_carrier_candidates(surfaces, edge_uses, &faces, start, end, bar) {
                if start.sub(end).length() <= bar {
                    if let (Ok(a), Ok(b)) = (candidate.domain().and_then(|d| candidate.derivatives(d[0],1)), edge.curve.derivatives(edge.t0,1)) {
                        if a[1].dot(b[1]) < 0.0 {
                            let Ok(reverse) = candidate.reversed() else { continue };
                            candidate = reverse;
                        }
                    }
                }
                let Ok([c0, c1]) = candidate.domain() else { continue };
                let Ok(after) = worst_residual(&candidate, c0, c1, &surfaces) else { continue };
                if after > bar { continue; }
                // Bidirectional geometric movement avoids accepting a remote
                // analytic branch or the long complementary arc. The entry
                // residual supplies the repair budget, capped by the existing
                // representation contract; it never enlarges acceptance.
                let limit = (8.0 * before).min(carrier_band);
                let mut moved = 0.0f64;
                for (from, a, b, to) in [(&candidate, c0, c1, &previous), (&previous, edge.t0, edge.t1, &candidate)] {
                    for t in carrier_stations(from, a, b) {
                        let Ok(point) = from.evaluate(t) else { moved = f64::INFINITY; break };
                        let Ok(projection) = crate::project_point_to_curve(to, point) else { moved = f64::INFINITY; break };
                        moved = moved.max(projection.distance);
                    }
                }
                if moved <= limit && best.as_ref().map_or(true, |(_, r, _)| after < *r)
                    && super::clearance::motion_clear(edge, &candidate, &faces, &solid.edges) {
                    best = Some((candidate, after, format!("exact carrier section, moved {moved:.3e}, limit {limit:.3e}")));
                }
            }
        }
        if best.as_ref().is_some_and(|(candidate,_,_)|!super::clearance::motion_clear(edge,candidate,&faces,&solid.edges)) {
            declined.push(format!("edge {}: complete old-to-new image motion cannot exclude a nonincident face",edge.id));
            best=None;
        }
        match best {
            Some((curve, after, note)) => {
                if debug {
                    eprintln!("edges: reconstructed edge {} ({} controls, {before:.3e} off) from {note}", edge.id, edge.curve.control_points.len());
                }
                candidates.push((edge_index, curve, before, after));
            }
            None if planar_uses == 0 => declined.push(format!("edge {} ({before:.3e} off): both carriers curved, no planar trim to take", edge.id)),
            None => {}
        }
    }

    if crate::imprint::profile_enabled() {
        let work = RESIDUAL_WORK.with(|work| work.get());
        eprintln!("assemble.reconstruct.work queries={} readings={} projections={}", RESIDUAL_QUERIES.with(|queries| queries.get()) - queries_before, work.0 - work_before.0, work.1 - work_before.1);
    }
    let mut replaced = Vec::new();
    for (edge_index, curve, residual_before, residual_after) in candidates {
        let edge = &mut solid.edges[edge_index];
        let [t0, t1] = curve.domain().or_refuse(KernelStage::Sew, "domain")?;
        let previous = std::mem::replace(&mut edge.curve, curve);
        let (previous_t0, previous_t1) = (edge.t0, edge.t1);
        edge.t0 = t0;
        edge.t1 = t1;
        replaced.push(ReplacedEdge { edge_id: edge.id, previous, t0: previous_t0, t1: previous_t1, residual_before, residual_after });
    }
    if debug {
        eprintln!("edges: {off_edges} edges over the {bar:.1e} bar; {} reconstructed, {} declined", replaced.len(), declined.len());
        for line in &declined {
            eprintln!("edges:   declined {line}");
        }
    }
    Ok(replaced)
}