rscad_core/feature/

edge.rs

//! Analytic edge resolution: turn a pair of leaf faces into a verified edge
//! curve on the **final solid**, or a typed rejection.
//!
//! Supported curves (v1): straight lines (plane ∩ plane, incl. cube edges)
//! and full circles (cylinder/cone rims, plane ∩ cylinder-side with the plane
//! ⊥ axis, plane ∩ sphere, plane ∩ cone-side ⊥ axis). Everything else is
//! rejected with a reason a frontend can show ([`EdgeReject::message`]).
//!
//! Leaf transforms may scale round surfaces arbitrarily: axisymmetric scale
//! about the leaf's own axis keeps the curves circular in the evaluation
//! frame ([`axisymmetric_scale`]); anisotropic scale makes them elliptical,
//! so resolution retries in the leaf's local frame where they are exact
//! circles again, and the result carries that frame
//! ([`CircularEdge::frame`] — the world curve is its affine image). Only
//! pairs mixing two *differently* anisotropic round leaves reject
//! ([`EdgeReject::NonUniformScale`]).
//!
//! Resolution is *verified* against the final solid, not just the two
//! surfaces: sample points along the candidate curve must lie on the actual
//! boundary (membership probes straddling the edge bisector), and the
//! surviving side of each face is found by probing — so an edge fully
//! consumed by a boolean resolves to nothing, and a partially consumed
//! straight edge is trimmed to its longest surviving run. Circles must
//! survive at every sample (no partial-arc features in v1).
//!
//! Convexity is classified geometrically from the probed face directions:
//! with `d̂A` pointing along face A's surviving surface away from the edge,
//! the edge is convex iff `d̂A·n̂B < 0` (the faces bend away from each other).
//! Both faces must agree. Straight concave edges are rejected in v1; concave
//! *circles* (boss-base rings) resolve fine — fillet lowering uses them.

use crate::feature::FaceTag;
use crate::feature::geom::{CsgExpr, CsgLeaves, GeomLeaf, LeafPrim};
use crate::prelude_::*;

/// Membership-probe offset as a fraction of the model's AABB diagonal.
const PROBE_FRAC: f64 = 1.0e-4;
/// Samples along a straight edge for existence verification.
const LINE_SAMPLES: usize = 9;
/// Samples around a circle for existence verification (all must pass).
const CIRCLE_SAMPLES: usize = 8;
/// Bisection iterations when trimming a partially surviving line.
const TRIM_ITERS: usize = 12;
/// Relative tolerance for axis-alignment / similarity checks.
const ALIGN_EPS: f64 = 1.0e-7;

// ── Public result types ──────────────────────────────────────────

/// Whether material bends away from (convex) or wraps around (concave) an
/// edge, on the final solid.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Convexity {
    Convex,
    Concave,
}

/// Why an edge candidate did not resolve. Frontends surface
/// [`message`](Self::message) verbatim.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum EdgeReject {
    /// Proxy/unknown faces, or a face tag that no longer matches its leaf.
    Unsupported,
    /// The surfaces intersect in something other than a line or full circle
    /// (oblique plane∩cylinder ellipse, cylinder∩cylinder, …).
    UnsupportedCurve,
    /// Round faces of two leaves under mismatched anisotropic scales: no
    /// single frame makes both surfaces canonical, so their intersection is
    /// not a line/circle in any frame. (A single anisotropic leaf — or one
    /// paired with planes — resolves in its own frame instead.)
    NonUniformScale,
    /// The faces don't meet in a crease on the final solid (miss, parallel,
    /// flush, tangent, or fully consumed by booleans).
    NotAnEdge,
    /// The circle exists geometrically but part of it is cut away.
    PartialCircle,
    /// Straight concave creases are not supported in v1.
    StraightConcave,
}

impl EdgeReject {
    pub fn message(self) -> &'static str {
        match self {
            Self::Unsupported => "Edge not supported on this surface",
            Self::UnsupportedCurve => "Edge curve not supported (line/circle only)",
            Self::NonUniformScale => "Mismatched anisotropic scales not supported",
            Self::NotAnEdge => "No edge between these faces",
            Self::PartialCircle => "Partial circular edges not supported",
            Self::StraightConcave => "Concave straight edges not supported",
        }
    }
}

/// Resolved edge geometry, in the evaluation frame.
#[derive(Clone, Copy, Debug)]
pub enum EdgeGeom {
    Straight(StraightEdge),
    Circular(CircularEdge),
}

/// A verified straight edge on the final solid, in the evaluation frame.
#[derive(Clone, Copy, Debug)]
pub struct StraightEdge {
    pub start: DVec3,
    pub end: DVec3,
    /// Final-solid outward normals of the two faces along the edge.
    pub na: DVec3,
    pub nb: DVec3,
    /// Unit in-face directions away from the edge, along each surviving face.
    pub da: DVec3,
    pub db: DVec3,
    /// How far each face extends from the edge along `da`/`db` (setback room).
    pub extent_a: f64,
    pub extent_b: f64,
}

impl StraightEdge {
    pub fn dir(&self) -> DVec3 {
        (self.end - self.start).normalize_or_zero()
    }
}

/// A verified full-circle edge (rim), in `frame` coordinates.
///
/// The meridian frame at any rim point has x = outward radial (away from the
/// axis) and y = `axis`. `axis` is oriented as the **planar face's final
/// outward normal**, so the planar face is always `y = 0` with normal `+y`
/// in meridian coordinates.
///
/// `frame` is identity for rims that are circles in the evaluation frame.
/// Rims that are *elliptical* in the evaluation frame (anisotropically
/// scaled leaves) resolve in the owning leaf's local frame, where they are
/// exact circles; `frame` maps that frame to the evaluation frame, so the
/// world curve is the frame image of this circle. All fields (center, axis,
/// radius, meridian data, extents) are frame-local.
#[derive(Clone, Copy, Debug)]
pub struct CircularEdge {
    /// Frame the circle lives in (frame → evaluation coordinates).
    pub frame: DAffine3,
    pub center: DVec3,
    /// Unit circle-plane normal — the planar face's final outward normal.
    pub axis: DVec3,
    pub radius: f64,
    /// Angular origin of the curved leaf's tessellation: its local `+X`
    /// image (unit, ⊥ `axis`). Masks anchor their polygon phase here so
    /// their chords coincide with the leaf's facets in the mesh booleans.
    pub ref_dir: DVec3,
    /// Segment count of the curved leaf's render tessellation (`None` for
    /// spheres, which have no prism phase to match).
    pub leaf_segments: Option<usize>,
    /// Final outward normal of the curved face, meridian coords.
    pub nb2: DVec2,
    /// In-face directions away from the rim, meridian coords. Face A is the
    /// planar face (`da2.y == 0`).
    pub da2: DVec2,
    pub db2: DVec2,
    /// Setback room along `da2`/`db2`.
    pub extent_a: f64,
    pub extent_b: f64,
}

impl CircularEdge {
    /// Meridian frame at angle `theta` from `ref_dir` (leaf tessellation
    /// phase): (outward radial, rim point).
    pub fn radial_at(&self, theta: f64) -> (DVec3, DVec3) {
        let (u, v) = (self.ref_dir, self.ref_dir.cross(self.axis));
        let r_hat = theta.cos() * u + theta.sin() * v;
        (r_hat, self.center + self.radius * r_hat)
    }
}

/// Distance from `p` to the segment `[a, b]`.
fn point_segment_distance(p: DVec3, a: DVec3, b: DVec3) -> f64 {
    let ab = b - a;
    let len_sq = ab.length_squared();
    if len_sq < 1.0e-24 {
        return (p - a).length();
    }
    let t = ((p - a).dot(ab) / len_sq).clamp(0.0, 1.0);
    (p - (a + t * ab)).length()
}

/// Samples used to approximate distances to frame-mapped (elliptical) rims.
const ELLIPSE_DISTANCE_SEGMENTS: usize = 64;

impl EdgeGeom {
    /// Shortest distance from evaluation-frame `p` to the edge curve.
    /// Analytic for lines and true circles; frame-mapped rims (world
    /// ellipses) fall back to a fine polyline approximation.
    pub fn distance_to(&self, p: DVec3) -> f64 {
        match self {
            EdgeGeom::Straight(e) => point_segment_distance(p, e.start, e.end),
            EdgeGeom::Circular(c) if c.frame == DAffine3::IDENTITY => {
                let rel = p - c.center;
                let axial = rel.dot(c.axis);
                let radial = rel - axial * c.axis;
                let rho = radial.length();
                if rho < 1.0e-12 {
                    // On the axis: every rim point is equidistant.
                    return (axial * axial + c.radius * c.radius).sqrt();
                }
                let nearest = c.center + radial / rho * c.radius;
                (p - nearest).length()
            }
            EdgeGeom::Circular(_) => {
                let pts = self.polyline(ELLIPSE_DISTANCE_SEGMENTS);
                (0..pts.len())
                    .map(|i| point_segment_distance(p, pts[i], pts[(i + 1) % pts.len()]))
                    .fold(f64::INFINITY, f64::min)
            }
        }
    }

    /// Evaluation-frame polyline for viewport highlighting (closed for
    /// circles; frame-mapped rims trace the world ellipse).
    pub fn polyline(&self, segments: usize) -> Vec<DVec3> {
        match self {
            EdgeGeom::Straight(e) => vec![e.start, e.end],
            EdgeGeom::Circular(c) => (0..segments.max(3))
                .map(|i| {
                    let theta = core::f64::consts::TAU * i as f64 / segments.max(3) as f64;
                    c.frame.transform_point3(c.radial_at(theta).1)
                })
                .collect(),
        }
    }

    /// Re-express a resolution done in `frame` coordinates in the evaluation
    /// frame. Straight edges map exactly (lines stay lines; the in-face
    /// directions are re-derived orthogonal to the mapped edge). Circles keep
    /// their frame-local data and compose `frame` (the world curve is the
    /// frame image — an ellipse in general).
    fn into_world_from(self, frame: &DAffine3) -> EdgeGeom {
        match self {
            EdgeGeom::Straight(e) => {
                let m = frame.matrix3;
                let n_it = m.inverse().transpose();
                let na = (n_it * e.na).normalize_or_zero();
                let nb = (n_it * e.nb).normalize_or_zero();
                let start = frame.transform_point3(e.start);
                let end = frame.transform_point3(e.end);
                let u = (end - start).normalize_or_zero();
                // The edge line lies in both face planes, so u ⊥ na, nb: the
                // world in-face direction is u × n, signed by the mapped one.
                let redir = |n: DVec3, d_frame: DVec3| {
                    let c = u.cross(n).normalize_or_zero();
                    if c.dot(m * d_frame) >= 0.0 { c } else { -c }
                };
                EdgeGeom::Straight(StraightEdge {
                    start,
                    end,
                    na,
                    nb,
                    da: redir(na, e.da),
                    db: redir(nb, e.db),
                    extent_a: (m * (e.da * e.extent_a)).length(),
                    extent_b: (m * (e.db * e.extent_b)).length(),
                })
            }
            EdgeGeom::Circular(c) => EdgeGeom::Circular(CircularEdge {
                frame: *frame * c.frame,
                ..c
            }),
        }
    }
}

// ── Leaf helpers ─────────────────────────────────────────────────

/// Face tags a leaf primitive exposes (candidate enumeration).
pub fn leaf_face_tags(prim: &LeafPrim) -> &'static [FaceTag] {
    match prim {
        LeafPrim::Cube { .. } => &[
            FaceTag::CubePosX,
            FaceTag::CubeNegX,
            FaceTag::CubePosY,
            FaceTag::CubeNegY,
            FaceTag::CubePosZ,
            FaceTag::CubeNegZ,
        ],
        LeafPrim::Cylinder { .. } => &[FaceTag::CylSide, FaceTag::CylTop, FaceTag::CylBottom],
        LeafPrim::Cone { .. } => &[FaceTag::ConeSide, FaceTag::ConeBottom],
        LeafPrim::Sphere { .. } => &[FaceTag::Sphere],
        LeafPrim::Aabb { .. } => &[],
    }
}

/// Per-leaf orientation: `true` when the leaf sits under an odd number of
/// `Difference` subtrahend positions, so its surfaces bound the final solid
/// with *inverted* normals (hole walls).
pub fn leaf_flipped(csg: &CsgLeaves) -> Vec<bool> {
    fn walk(expr: &CsgExpr, flipped: bool, out: &mut Vec<(usize, bool)>) {
        match expr {
            CsgExpr::Leaf(i) => out.push((*i, flipped)),
            CsgExpr::Union(cs) | CsgExpr::Inter(cs) => {
                cs.iter().for_each(|c| walk(c, flipped, out));
            }
            CsgExpr::Diff(cs) => {
                if let Some((first, rest)) = cs.split_first() {
                    walk(first, flipped, out);
                    rest.iter().for_each(|c| walk(c, !flipped, out));
                }
            }
        }
    }
    let mut pairs = Vec::with_capacity(csg.leaves.len());
    walk(&csg.expr, false, &mut pairs);
    let mut flags = vec![false; csg.leaves.len()];
    for (i, f) in pairs {
        flags[i] = f;
    }
    flags
}

/// Evaluation-frame AABB over every leaf (probe scaling).
pub fn csg_aabb(csg: &CsgLeaves) -> Option<(DVec3, DVec3)> {
    let mut bounds: Option<(DVec3, DVec3)> = None;
    for leaf in &csg.leaves {
        let (lo, hi) = match leaf.prim {
            LeafPrim::Cube { size } => (DVec3::ZERO, size),
            LeafPrim::Cylinder { radius, height, .. } | LeafPrim::Cone { radius, height, .. } => (
                DVec3::new(-radius, 0.0, -radius),
                DVec3::new(radius, height, radius),
            ),
            LeafPrim::Sphere { radius } => (DVec3::splat(-radius), DVec3::splat(radius)),
            LeafPrim::Aabb { min, max } => (min, max),
        };
        for i in 0..8 {
            let corner = DVec3::new(
                if i & 1 == 0 { lo.x } else { hi.x },
                if i & 2 == 0 { lo.y } else { hi.y },
                if i & 4 == 0 { lo.z } else { hi.z },
            );
            let p = leaf.world.transform_point3(corner);
            bounds = Some(match bounds {
                None => (p, p),
                Some((mn, mx)) => (mn.min(p), mx.max(p)),
            });
        }
    }
    bounds
}

/// Radial/axial scale factors `(s_r, s_h)` of an affine that keeps a local
/// Y-axisymmetric surface a surface of revolution: the Gram matrix `MᵀM`
/// must be `diag(s_r², s_h², s_r²)` — the image frame stays orthogonal with
/// equal X/Z stretch, so cross-section circles map to circles. `None` when
/// the linear part shears or squashes the cross-section anisotropically
/// (circles become ellipses). Tolerances are relative per entry, so extreme
/// radial/axial aspect ratios don't loosen the radial check.
fn axisymmetric_scale(world: &DAffine3) -> Option<(f64, f64)> {
    let m = world.matrix3;
    let g = m.transpose() * m;
    let r2 = (g.x_axis.x + g.z_axis.z) / 2.0;
    let h2 = g.y_axis.y;
    if r2 <= 0.0 || h2 <= 0.0 {
        return None;
    }
    // |g_xy|, |g_yz| ≤ √(g_xx·g_yy) by Cauchy–Schwarz: bound them mixed.
    let mix = ALIGN_EPS * (r2 * h2).sqrt();
    let ok = (g.x_axis.x - r2).abs() < ALIGN_EPS * r2
        && g.x_axis.z.abs() < ALIGN_EPS * r2
        && g.x_axis.y.abs() < mix
        && g.y_axis.z.abs() < mix;
    ok.then(|| (r2.sqrt(), h2.sqrt()))
}

/// Uniform scale factor of an affine similarity, or `None` when the linear
/// part shears/scales non-uniformly (spheres stop being spheres).
fn similarity_scale(world: &DAffine3) -> Option<f64> {
    let (s_r, s_h) = axisymmetric_scale(world)?;
    let (r2, h2) = (s_r * s_r, s_h * s_h);
    ((r2 - h2).abs() < ALIGN_EPS * r2.max(h2)).then(|| ((2.0 * r2 + h2) / 3.0).sqrt())
}

// ── Face surfaces ────────────────────────────────────────────────

/// In-plane extent of a planar face, for clipping edge candidates.
#[derive(Clone, Copy, Debug)]
enum PlaneExtent {
    /// Affine image of the unit square: `origin + u·eu + v·ev`, u,v ∈ [0,1].
    Quad { origin: DVec3, eu: DVec3, ev: DVec3 },
    /// Disc about `center` (uniform-scale radius).
    Disc { center: DVec3, radius: f64 },
}

/// Analytic surface of one leaf face in the evaluation frame. Round surfaces
/// only exist under similarity transforms (checked at construction).
#[derive(Clone, Copy, Debug)]
enum Surface {
    Plane {
        point: DVec3,
        /// Leaf-orientation outward unit normal (flip applied by callers).
        normal: DVec3,
        extent: PlaneExtent,
    },
    Cylinder {
        base: DVec3,
        axis: DVec3,
        radius: f64,
        height: f64,
    },
    Cone {
        base: DVec3,
        axis: DVec3,
        radius: f64,
        height: f64,
    },
    Sphere {
        center: DVec3,
        radius: f64,
    },
}

impl Surface {
    /// Leaf-orientation outward normal at a point on the surface.
    fn normal_at(&self, p: DVec3) -> DVec3 {
        match self {
            Surface::Plane { normal, .. } => *normal,
            Surface::Cylinder { base, axis, .. } => {
                let rel = p - *base;
                (rel - rel.dot(*axis) * *axis).normalize_or_zero()
            }
            Surface::Cone {
                base,
                axis,
                radius,
                height,
            } => {
                let rel = p - *base;
                let r_hat = (rel - rel.dot(*axis) * *axis).normalize_or_zero();
                if r_hat == DVec3::ZERO {
                    return *axis;
                }
                (r_hat + radius / height * *axis).normalize()
            }
            Surface::Sphere { center, .. } => (p - *center).normalize_or_zero(),
        }
    }

    /// Project a nearby point back onto the surface (survival probes walk
    /// along curved faces).
    fn project(&self, p: DVec3) -> DVec3 {
        match self {
            Surface::Plane { point, normal, .. } => p - (p - *point).dot(*normal) * *normal,
            Surface::Cylinder {
                base, axis, radius, ..
            } => {
                let rel = p - *base;
                let y = rel.dot(*axis);
                let radial = rel - y * *axis;
                let r_hat = radial.normalize_or_zero();
                *base + y * *axis + *radius * r_hat
            }
            Surface::Cone {
                base,
                axis,
                radius,
                height,
            } => {
                let rel = p - *base;
                let y = rel.dot(*axis).clamp(0.0, *height);
                let radial = rel - rel.dot(*axis) * *axis;
                let r_hat = radial.normalize_or_zero();
                let rho = radius * (1.0 - y / height);
                *base + y * *axis + rho * r_hat
            }
            Surface::Sphere { center, radius } => {
                *center + (p - *center).normalize_or_zero() * *radius
            }
        }
    }
}

/// Surface of `(leaf, face)` in the evaluation frame, or a rejection.
fn face_surface(leaf: &GeomLeaf, face: FaceTag) -> Result<Surface, EdgeReject> {
    let w = &leaf.world;
    let normal = |n_local: DVec3| (leaf.world_inv.matrix3.transpose() * n_local).normalize();
    let quad = |axis: usize, positive: bool, size: DVec3| {
        let mut corner = DVec3::ZERO;
        corner[axis] = if positive { size[axis] } else { 0.0 };
        let (u, v) = ((axis + 1) % 3, (axis + 2) % 3);
        let mut eu = DVec3::ZERO;
        eu[u] = size[u];
        let mut ev = DVec3::ZERO;
        ev[v] = size[v];
        PlaneExtent::Quad {
            origin: w.transform_point3(corner),
            eu: w.matrix3 * eu,
            ev: w.matrix3 * ev,
        }
    };
    let similarity = || similarity_scale(w).ok_or(EdgeReject::NonUniformScale);
    let axisym = || axisymmetric_scale(w).ok_or(EdgeReject::NonUniformScale);
    match (&leaf.prim, face) {
        (LeafPrim::Cube { size }, tag) => {
            let (axis, positive) = match tag {
                FaceTag::CubePosX => (0, true),
                FaceTag::CubeNegX => (0, false),
                FaceTag::CubePosY => (1, true),
                FaceTag::CubeNegY => (1, false),
                FaceTag::CubePosZ => (2, true),
                FaceTag::CubeNegZ => (2, false),
                _ => return Err(EdgeReject::Unsupported),
            };
            let mut n_local = DVec3::ZERO;
            n_local[axis] = if positive { 1.0 } else { -1.0 };
            let mut on_face = DVec3::ZERO;
            on_face[axis] = if positive { size[axis] } else { 0.0 };
            Ok(Surface::Plane {
                point: w.transform_point3(on_face),
                normal: normal(n_local),
                extent: quad(axis, positive, *size),
            })
        }
        (LeafPrim::Cylinder { radius, height, .. }, FaceTag::CylTop | FaceTag::CylBottom) => {
            let (s_r, _) = axisym()?;
            let y = if face == FaceTag::CylTop {
                *height
            } else {
                0.0
            };
            let n_local = if face == FaceTag::CylTop {
                DVec3::Y
            } else {
                -DVec3::Y
            };
            let center = w.transform_point3(DVec3::new(0.0, y, 0.0));
            Ok(Surface::Plane {
                point: center,
                normal: normal(n_local),
                extent: PlaneExtent::Disc {
                    center,
                    radius: radius * s_r,
                },
            })
        }
        (LeafPrim::Cylinder { radius, height, .. }, FaceTag::CylSide) => {
            let (s_r, s_h) = axisym()?;
            Ok(Surface::Cylinder {
                base: w.transform_point3(DVec3::ZERO),
                axis: (w.matrix3 * DVec3::Y).normalize(),
                radius: radius * s_r,
                height: height * s_h,
            })
        }
        (LeafPrim::Cone { radius, .. }, FaceTag::ConeBottom) => {
            let (s_r, _) = axisym()?;
            let center = w.transform_point3(DVec3::ZERO);
            Ok(Surface::Plane {
                point: center,
                normal: normal(-DVec3::Y),
                extent: PlaneExtent::Disc {
                    center,
                    radius: radius * s_r,
                },
            })
        }
        (LeafPrim::Cone { radius, height, .. }, FaceTag::ConeSide) => {
            let (s_r, s_h) = axisym()?;
            Ok(Surface::Cone {
                base: w.transform_point3(DVec3::ZERO),
                axis: (w.matrix3 * DVec3::Y).normalize(),
                radius: radius * s_r,
                height: height * s_h,
            })
        }
        (LeafPrim::Sphere { radius }, FaceTag::Sphere) => {
            let s = similarity()?;
            Ok(Surface::Sphere {
                center: w.transform_point3(DVec3::ZERO),
                radius: radius * s,
            })
        }
        _ => Err(EdgeReject::Unsupported),
    }
}

// ── Probing ──────────────────────────────────────────────────────

struct Prober<'a> {
    csg: &'a CsgLeaves,
    /// Membership-probe offset (model-scaled).
    delta: f64,
}

impl<'a> Prober<'a> {
    fn new(csg: &'a CsgLeaves) -> Self {
        let diag = csg_aabb(csg)
            .map(|(mn, mx)| (mx - mn).length())
            .unwrap_or(1.0)
            .max(1.0e-6);
        Self {
            csg,
            delta: PROBE_FRAC * diag,
        }
    }

    /// `q` lies on the final solid's boundary crease between two faces whose
    /// summed outward normal is `bisector` (unnormalized OK).
    fn edge_exists(&self, q: DVec3, bisector: DVec3) -> bool {
        let b = bisector.normalize_or_zero();
        if b == DVec3::ZERO {
            return false;
        }
        self.csg.contains(q - self.delta * b) && !self.csg.contains(q + self.delta * b)
    }

    /// Whether `surface` survives on the final solid just off `q` along `d`
    /// (walk δ along the face, then straddle it by δ/2).
    fn face_survives(&self, surface: &Surface, n_final_sign: f64, q: DVec3, d: DVec3) -> bool {
        let p = surface.project(q + self.delta * d);
        let n = surface.normal_at(p) * n_final_sign;
        if n == DVec3::ZERO {
            return false;
        }
        let h = 0.5 * self.delta;
        self.csg.contains(p - h * n) && !self.csg.contains(p + h * n)
    }

    /// Pick the in-face direction (out of `±candidate`) along which the face
    /// survives. Exactly one side must survive.
    fn orient(
        &self,
        surface: &Surface,
        n_sign: f64,
        q: DVec3,
        candidate: DVec3,
    ) -> Result<DVec3, EdgeReject> {
        let pos = self.face_survives(surface, n_sign, q, candidate);
        let neg = self.face_survives(surface, n_sign, q, -candidate);
        match (pos, neg) {
            (true, false) => Ok(candidate),
            (false, true) => Ok(-candidate),
            _ => Err(EdgeReject::NotAnEdge),
        }
    }
}

/// Convexity from the probed face directions: convex iff each face's
/// in-face direction points *behind* the other face (`d·n < 0`). The two
/// faces must agree.
fn classify(da_dot_nb: f64, db_dot_na: f64) -> Result<Convexity, EdgeReject> {
    let eps = 1.0e-9;
    if da_dot_nb < -eps && db_dot_na < -eps {
        Ok(Convexity::Convex)
    } else if da_dot_nb > eps && db_dot_na > eps {
        Ok(Convexity::Concave)
    } else {
        Err(EdgeReject::NotAnEdge) // tangent or inconsistent corner
    }
}

// ── Extent clipping ──────────────────────────────────────────────

/// Interval of `t` where `q0 + t·dir` stays inside the extent. `dir` need
/// not be unit. Returns `None` when the line misses entirely.
fn clip_to_extent(extent: &PlaneExtent, q0: DVec3, dir: DVec3) -> Option<(f64, f64)> {
    match extent {
        PlaneExtent::Quad { origin, eu, ev } => {
            // Dual (gram) solve for the affine params of a point in the plane.
            let g00 = eu.dot(*eu);
            let g01 = eu.dot(*ev);
            let g11 = ev.dot(*ev);
            let det = g00 * g11 - g01 * g01;
            if det.abs() < 1.0e-18 {
                return None;
            }
            let params = |p: DVec3| {
                let r = p - *origin;
                let (bu, bv) = (r.dot(*eu), r.dot(*ev));
                ((bu * g11 - bv * g01) / det, (bv * g00 - bu * g01) / det)
            };
            let (u0, v0) = params(q0);
            let (u1, v1) = params(q0 + dir);
            let (du, dv) = (u1 - u0, v1 - v0);
            // Edges lie exactly on the face boundary (u or v = 0/1), so the
            // constant-parameter test needs slack for accumulated transform
            // error, or every rotated cube edge would clip to nothing.
            const PARAM_EPS: f64 = 1.0e-9;
            let mut lo = f64::NEG_INFINITY;
            let mut hi = f64::INFINITY;
            for (start, delta) in [(u0, du), (v0, dv)] {
                if delta.abs() < 1.0e-12 {
                    if !(-PARAM_EPS..=1.0 + PARAM_EPS).contains(&start) {
                        return None;
                    }
                    continue;
                }
                let (t0, t1) = ((0.0 - start) / delta, (1.0 - start) / delta);
                lo = lo.max(t0.min(t1));
                hi = hi.min(t0.max(t1));
            }
            (lo < hi).then_some((lo, hi))
        }
        PlaneExtent::Disc { center, radius } => {
            let rel = q0 - *center;
            let a = dir.dot(dir);
            if a < 1.0e-18 {
                return (rel.length() <= *radius).then_some((f64::NEG_INFINITY, f64::INFINITY));
            }
            let b = 2.0 * rel.dot(dir);
            let c = rel.dot(rel) - radius * radius;
            let disc = b * b - 4.0 * a * c;
            (disc > 0.0).then(|| {
                let sq = disc.sqrt();
                ((-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a))
            })
        }
    }
}

/// How far the face extends from `q` along unit `d` (setback room).
fn extent_along(extent: &PlaneExtent, q: DVec3, d: DVec3) -> f64 {
    clip_to_extent(extent, q, d).map_or(0.0, |(_, hi)| hi.max(0.0))
}

// ── Resolution ───────────────────────────────────────────────────

/// Resolve the edge where faces `(leaf_a, face_a)` and `(leaf_b, face_b)`
/// meet on the final solid. See the module docs for what resolves and why
/// things reject.
///
/// Anisotropically scaled round leaves make the evaluation-frame attempt
/// reject with [`EdgeReject::NonUniformScale`]; the resolution is then
/// retried in each round face's leaf frame (where that leaf's surfaces are
/// canonical), and the result is mapped back — circles keep the frame
/// ([`CircularEdge::frame`]), straight cap creases map exactly.
pub fn resolve_edge(
    csg: &CsgLeaves,
    a: (usize, FaceTag),
    b: (usize, FaceTag),
) -> Result<(EdgeGeom, Convexity), EdgeReject> {
    let world = resolve_edge_framed(csg, a, b);
    let Err(EdgeReject::NonUniformScale) = world else {
        return world;
    };
    let mut err = EdgeReject::NonUniformScale;
    let mut tried: Option<usize> = None;
    for (li, tag) in [a, b] {
        // Only round-face leaves define a useful frame (cube faces are
        // planes under any affine and never reject on scale).
        let round = matches!(
            tag,
            FaceTag::CylSide
                | FaceTag::CylTop
                | FaceTag::CylBottom
                | FaceTag::ConeSide
                | FaceTag::ConeBottom
                | FaceTag::Sphere
        );
        if !round || tried == Some(li) {
            continue;
        }
        tried = Some(li);
        let Some(leaf) = csg.leaves.get(li) else {
            continue;
        };
        let frame = leaf.world;
        let reframed = reframe(csg, &leaf.world_inv, &frame);
        match resolve_edge_framed(&reframed, a, b) {
            Ok((geom, convexity)) => return Ok((geom.into_world_from(&frame), convexity)),
            Err(e) => err = e,
        }
    }
    Err(err)
}

/// The same solid expressed in the frame `F` (world = `F` · frame coords):
/// every leaf transform is left-multiplied by `F⁻¹`.
fn reframe(csg: &CsgLeaves, f_inv: &DAffine3, f: &DAffine3) -> CsgLeaves {
    CsgLeaves {
        leaves: csg
            .leaves
            .iter()
            .map(|l| GeomLeaf {
                prim: l.prim,
                world: *f_inv * l.world,
                // (F⁻¹ · W)⁻¹ = W⁻¹ · F
                world_inv: l.world_inv * *f,
            })
            .collect(),
        expr: csg.expr.clone(),
        paths: csg.paths.clone(),
    }
}

/// One resolution attempt with all geometry read in the given CSG's own
/// coordinates (the evaluation frame, or a leaf frame after [`reframe`]).
fn resolve_edge_framed(
    csg: &CsgLeaves,
    a: (usize, FaceTag),
    b: (usize, FaceTag),
) -> Result<(EdgeGeom, Convexity), EdgeReject> {
    if a == b {
        return Err(EdgeReject::NotAnEdge);
    }
    let flipped = leaf_flipped(csg);
    let leaf_a = csg.leaves.get(a.0).ok_or(EdgeReject::Unsupported)?;
    let leaf_b = csg.leaves.get(b.0).ok_or(EdgeReject::Unsupported)?;
    let surf_a = face_surface(leaf_a, a.1)?;
    let surf_b = face_surface(leaf_b, b.1)?;
    let sign_a = if flipped[a.0] { -1.0 } else { 1.0 };
    let sign_b = if flipped[b.0] { -1.0 } else { 1.0 };
    let prober = Prober::new(csg);

    // Same-leaf cylinder/cone rims come from parameters directly (the generic
    // plane∩surface path would see the cap plane flush with the wall's end).
    if a.0 == b.0
        && let Some(circle) = same_leaf_rim(leaf_a, a.1, b.1)?
    {
        let plane_first = matches!(
            (a.1, b.1),
            (
                FaceTag::CylTop | FaceTag::CylBottom | FaceTag::ConeBottom,
                _
            )
        );
        let (plane, curved, sp, sc) = if plane_first {
            (&surf_a, &surf_b, sign_a, sign_b)
        } else {
            (&surf_b, &surf_a, sign_b, sign_a)
        };
        return resolve_circle(&prober, circle, plane, curved, leaf_a, sp, sc, true);
    }

    match (&surf_a, &surf_b) {
        (Surface::Plane { .. }, Surface::Plane { .. }) => {
            resolve_line(&prober, &surf_a, &surf_b, sign_a, sign_b)
        }
        (Surface::Plane { .. }, _) => {
            let circle = plane_curved_circle(&surf_a, &surf_b)?;
            resolve_circle(
                &prober, circle, &surf_a, &surf_b, leaf_b, sign_a, sign_b, false,
            )
        }
        (_, Surface::Plane { .. }) => {
            let circle = plane_curved_circle(&surf_b, &surf_a)?;
            resolve_circle(
                &prober, circle, &surf_b, &surf_a, leaf_a, sign_b, sign_a, false,
            )
        }
        _ => Err(EdgeReject::UnsupportedCurve),
    }
}

/// Same-leaf rim circle from primitive parameters (`None`: not a rim pair).
fn same_leaf_rim(
    leaf: &GeomLeaf,
    ta: FaceTag,
    tb: FaceTag,
) -> Result<Option<(DVec3, DVec3, f64)>, EdgeReject> {
    use FaceTag as F;
    let pair = |x: FaceTag, y: FaceTag| (ta == x && tb == y) || (ta == y && tb == x);
    let (y_local, radius_local) = match &leaf.prim {
        LeafPrim::Cylinder { radius, height, .. } => {
            if pair(F::CylSide, F::CylTop) {
                (*height, *radius)
            } else if pair(F::CylSide, F::CylBottom) {
                (0.0, *radius)
            } else {
                return if pair(F::CylTop, F::CylBottom) {
                    Err(EdgeReject::NotAnEdge)
                } else {
                    Ok(None)
                };
            }
        }
        LeafPrim::Cone { radius, .. } if pair(F::ConeSide, F::ConeBottom) => (0.0, *radius),
        _ => return Ok(None),
    };
    let (s_r, _) = axisymmetric_scale(&leaf.world).ok_or(EdgeReject::NonUniformScale)?;
    let center = leaf.world.transform_point3(DVec3::new(0.0, y_local, 0.0));
    let axis = (leaf.world.matrix3 * DVec3::Y).normalize();
    Ok(Some((center, axis, radius_local * s_r)))
}

/// Circle where a plane meets a curved surface, or a typed rejection.
/// Returns `(center, plane_normal_of_circle, radius)` (axis unsigned).
fn plane_curved_circle(
    plane: &Surface,
    curved: &Surface,
) -> Result<(DVec3, DVec3, f64), EdgeReject> {
    let Surface::Plane { point, normal, .. } = plane else {
        return Err(EdgeReject::UnsupportedCurve);
    };
    match curved {
        Surface::Cylinder {
            base,
            axis,
            radius,
            height,
        } => {
            if normal.dot(*axis).abs() < 1.0 - ALIGN_EPS {
                return Err(EdgeReject::UnsupportedCurve); // oblique → ellipse
            }
            // Closed range: a plane flush with a cap still meets the wall in
            // that cap's rim circle (boss bases, stacked parts) — the
            // existence probes decide whether the crease is real.
            let y = (*point - *base).dot(*axis);
            let tol = 1.0e-9 * height.max(1.0);
            if !(-tol..=height + tol).contains(&y) {
                return Err(EdgeReject::NotAnEdge);
            }
            Ok((*base + y.clamp(0.0, *height) * *axis, *axis, *radius))
        }
        Surface::Cone {
            base,
            axis,
            radius,
            height,
        } => {
            if normal.dot(*axis).abs() < 1.0 - ALIGN_EPS {
                return Err(EdgeReject::UnsupportedCurve);
            }
            let y = (*point - *base).dot(*axis);
            let tol = 1.0e-9 * height.max(1.0);
            // Open at the apex (a point, not a circle); closed at the base.
            if y < -tol || y >= *height - tol {
                return Err(EdgeReject::NotAnEdge);
            }
            let y = y.max(0.0);
            Ok((*base + y * *axis, *axis, radius * (1.0 - y / height)))
        }
        Surface::Sphere { center, radius } => {
            let d = (*center - *point).dot(*normal);
            if d.abs() >= *radius {
                return Err(EdgeReject::NotAnEdge);
            }
            let r = (radius * radius - d * d).sqrt();
            if r < 1.0e-9 * radius {
                return Err(EdgeReject::NotAnEdge); // tangent point
            }
            Ok((*center - d * *normal, *normal, r))
        }
        _ => Err(EdgeReject::UnsupportedCurve),
    }
}

/// Resolve a straight plane∩plane edge: intersect, clip to both faces,
/// verify against the final solid, trim to the surviving run.
fn resolve_line(
    prober: &Prober,
    surf_a: &Surface,
    surf_b: &Surface,
    sign_a: f64,
    sign_b: f64,
) -> Result<(EdgeGeom, Convexity), EdgeReject> {
    let (
        Surface::Plane {
            point: pa,
            normal: na_leaf,
            extent: ext_a,
        },
        Surface::Plane {
            point: pb,
            normal: nb_leaf,
            extent: ext_b,
        },
    ) = (surf_a, surf_b)
    else {
        return Err(EdgeReject::UnsupportedCurve);
    };
    let na = sign_a * *na_leaf;
    let nb = sign_b * *nb_leaf;
    let u = na_leaf.cross(*nb_leaf);
    let u_len = u.length();
    if u_len < 1.0e-9 {
        return Err(EdgeReject::NotAnEdge); // parallel/coplanar faces
    }
    let u_hat = u / u_len;
    // Point on both planes: q = (d1·(n2×u) − d2·(n1×u)) / |u|².
    let (d1, d2) = (pa.dot(*na_leaf), pb.dot(*nb_leaf));
    let q0 = (d1 * nb_leaf.cross(u) - d2 * na_leaf.cross(u)) / (u_len * u_len);
    debug_assert!((q0.dot(*na_leaf) - d1).abs() < 1.0e-6);
    debug_assert!((q0.dot(*nb_leaf) - d2).abs() < 1.0e-6);

    // Clip to both faces' extents.
    let (a_lo, a_hi) = clip_to_extent(ext_a, q0, u_hat).ok_or(EdgeReject::NotAnEdge)?;
    let (b_lo, b_hi) = clip_to_extent(ext_b, q0, u_hat).ok_or(EdgeReject::NotAnEdge)?;
    let (lo, hi) = (a_lo.max(b_lo), a_hi.min(b_hi));
    if !(hi - lo).is_finite() || hi - lo < prober.delta {
        return Err(EdgeReject::NotAnEdge);
    }

    // Verify samples on the final solid; keep the longest surviving run.
    let bisector = na + nb;
    if bisector.length_squared() < 1.0e-12 {
        return Err(EdgeReject::NotAnEdge); // flush faces, no crease
    }
    let at = |t: f64| q0 + t * u_hat;
    let sample_t: Vec<f64> = (0..LINE_SAMPLES)
        .map(|i| {
            let f = (i as f64 + 0.5) / LINE_SAMPLES as f64;
            lo + f * (hi - lo)
        })
        .collect();
    let alive: Vec<bool> = sample_t
        .iter()
        .map(|&t| prober.edge_exists(at(t), bisector))
        .collect();
    let (run_start, run_len) = longest_run(&alive).ok_or(EdgeReject::NotAnEdge)?;

    // Trim the run's outer bounds: clipped ends when fully alive, else
    // bisect between the last alive and first dead sample.
    let exists = |t: f64| prober.edge_exists(at(t), bisector);
    let t_start = if run_start == 0 {
        lo
    } else {
        bisect(sample_t[run_start], sample_t[run_start - 1], exists)
    };
    let run_end = run_start + run_len - 1;
    let t_end = if run_end + 1 == LINE_SAMPLES {
        hi
    } else {
        bisect(sample_t[run_end], sample_t[run_end + 1], exists)
    };
    if t_end - t_start < prober.delta {
        return Err(EdgeReject::NotAnEdge);
    }

    // Orient the in-face directions at the run midpoint and classify.
    let q_mid = at(0.5 * (t_start + t_end));
    let da = prober.orient(surf_a, sign_a, q_mid, u_hat.cross(na).normalize())?;
    let db = prober.orient(surf_b, sign_b, q_mid, u_hat.cross(nb).normalize())?;
    let convexity = classify(da.dot(nb), db.dot(na))?;
    if convexity == Convexity::Concave {
        return Err(EdgeReject::StraightConcave);
    }

    let edge = StraightEdge {
        start: at(t_start),
        end: at(t_end),
        na,
        nb,
        da,
        db,
        extent_a: extent_along(ext_a, q_mid, da),
        extent_b: extent_along(ext_b, q_mid, db),
    };
    Ok((EdgeGeom::Straight(edge), convexity))
}

/// Resolve a circular rim: verify every sample, orient the meridian
/// directions, classify convexity.
#[allow(clippy::too_many_arguments)]
fn resolve_circle(
    prober: &Prober,
    (center, axis_raw, radius): (DVec3, DVec3, f64),
    plane: &Surface,
    curved: &Surface,
    curved_leaf: &GeomLeaf,
    sign_plane: f64,
    sign_curved: f64,
    same_leaf: bool,
) -> Result<(EdgeGeom, Convexity), EdgeReject> {
    let Surface::Plane {
        normal: np_leaf,
        extent: ext_p,
        ..
    } = plane
    else {
        return Err(EdgeReject::UnsupportedCurve);
    };
    // The circle axis, oriented as the planar face's FINAL outward normal.
    let np_final = sign_plane * *np_leaf;
    let axis = if axis_raw.dot(np_final) >= 0.0 {
        axis_raw
    } else {
        -axis_raw
    };
    // Tessellation phase of the curved leaf: its local +X, ⊥ the axis. Masks
    // anchored here share the leaf's chord positions in mesh booleans.
    let ref_dir = {
        let mx = curved_leaf.world.matrix3 * DVec3::X;
        let perp = mx - mx.dot(axis) * axis;
        let len = perp.length();
        if len > 1.0e-12 {
            perp / len
        } else {
            plane_basis(axis).0
        }
    };
    let leaf_segments = match curved_leaf.prim {
        LeafPrim::Cylinder { segments, .. } | LeafPrim::Cone { segments, .. } => Some(segments),
        _ => None,
    };
    let (bu, bv) = plane_basis(axis);

    // Containment in the planar face (skipped for same-leaf rims, where the
    // rim IS the disc boundary) + existence at every sample.
    let mut nb2_acc = DVec2::ZERO;
    for i in 0..CIRCLE_SAMPLES {
        let theta = core::f64::consts::TAU * i as f64 / CIRCLE_SAMPLES as f64;
        let r_hat = theta.cos() * bu + theta.sin() * bv;
        let q = center + radius * r_hat;
        if !same_leaf {
            let inside = match ext_p {
                PlaneExtent::Quad { .. } => {
                    clip_to_extent(ext_p, q, r_hat).is_some_and(|(lo, hi)| lo <= 0.0 && hi >= 0.0)
                }
                PlaneExtent::Disc {
                    center: c,
                    radius: r,
                } => (q - *c).length() <= *r + 1.0e-9,
            };
            if !inside {
                return Err(EdgeReject::PartialCircle);
            }
        }
        let nb3 = curved.normal_at(q) * sign_curved;
        if !prober.edge_exists(q, np_final + nb3) {
            return Err(EdgeReject::PartialCircle);
        }
        nb2_acc += DVec2::new(nb3.dot(r_hat), nb3.dot(axis));
    }
    let nb2 = (nb2_acc / CIRCLE_SAMPLES as f64).normalize_or_zero();
    if nb2 == DVec2::ZERO {
        return Err(EdgeReject::NotAnEdge);
    }

    // Orientation probes at θ = 0 (meridian frame: x = radial, y = axis).
    let r_hat = bu;
    let q = center + radius * r_hat;
    let da3 = prober.orient(plane, sign_plane, q, r_hat)?;
    let da2 = DVec2::new(da3.dot(r_hat), 0.0);
    // Curved-face meridian tangent candidates (along the surface's own axis,
    // which may oppose the circle axis — e.g. a bottom-cap rim).
    let tangent = match curved {
        Surface::Cylinder { axis: ca, .. } => *ca,
        Surface::Cone {
            axis: ca,
            radius,
            height,
            ..
        } => {
            let k = radius / height;
            (*ca - k * r_hat).normalize()
        }
        Surface::Sphere { center: c, .. } => {
            let n = (q - *c).normalize_or_zero();
            // Meridian tangent ⊥ sphere normal, in the (r̂, axis) plane.
            (axis - axis.dot(n) * n).normalize_or_zero()
        }
        Surface::Plane { .. } => return Err(EdgeReject::UnsupportedCurve),
    };
    if tangent == DVec3::ZERO {
        return Err(EdgeReject::NotAnEdge);
    }
    let db3 = prober.orient(curved, sign_curved, q, tangent)?;
    let db2 = DVec2::new(db3.dot(r_hat), db3.dot(axis)).normalize_or_zero();

    // Meridian normals: planar face is +y by construction.
    let na2 = DVec2::Y;
    let convexity = classify(da2.dot(nb2), db2.dot(na2))?;

    // Setback room: planar face along da3; curved face along its meridian.
    let extent_a = extent_along(ext_p, q, da3);
    let extent_b = match curved {
        Surface::Cylinder {
            base,
            axis: ca,
            height,
            ..
        } => {
            let y = (q - *base).dot(*ca);
            if db3.dot(*ca) >= 0.0 { height - y } else { y }
        }
        Surface::Cone { radius, height, .. } => (radius * radius + height * height).sqrt(),
        Surface::Sphere { radius, .. } => *radius,
        Surface::Plane { .. } => 0.0,
    };

    let edge = CircularEdge {
        frame: DAffine3::IDENTITY,
        center,
        axis,
        radius,
        ref_dir,
        leaf_segments,
        nb2,
        da2,
        db2,
        extent_a,
        extent_b: extent_b.max(0.0),
    };
    Ok((EdgeGeom::Circular(edge), convexity))
}

// ── Small utilities ──────────────────────────────────────────────

/// Longest run of `true` in `flags`: `(start, len)`.
fn longest_run(flags: &[bool]) -> Option<(usize, usize)> {
    let mut best: Option<(usize, usize)> = None;
    let mut cur: Option<(usize, usize)> = None;
    for (i, &f) in flags.iter().enumerate() {
        cur = match (cur, f) {
            (None, true) => Some((i, 1)),
            (Some((s, n)), true) => Some((s, n + 1)),
            (_, false) => None,
        };
        if let Some((_, n)) = cur
            && best.is_none_or(|(_, bn)| n > bn)
        {
            best = cur;
        }
    }
    best
}

/// Bisect between a passing `t_ok` and failing `t_bad` parameter.
fn bisect(mut t_ok: f64, mut t_bad: f64, exists: impl Fn(f64) -> bool) -> f64 {
    for _ in 0..TRIM_ITERS {
        let mid = 0.5 * (t_ok + t_bad);
        if exists(mid) {
            t_ok = mid;
        } else {
            t_bad = mid;
        }
    }
    t_ok
}

/// Orthonormal basis of the plane with unit normal `n`.
pub(crate) fn plane_basis(n: DVec3) -> (DVec3, DVec3) {
    let helper = if n.y.abs() < 0.9 { DVec3::Y } else { DVec3::X };
    let u = n.cross(helper).normalize();
    (u, n.cross(u))
}

#[cfg(test)]
mod tests {
    use super::*;

    fn cube(size: [f64; 3]) -> DynamicObject {
        DynamicObject::Cube { size }
    }

    fn cylinder(radius: f64, height: f64) -> DynamicObject {
        DynamicObject::Cylinder {
            radius,
            height,
            segments: 32,
        }
    }

    fn translate(offset: [f64; 3], child: DynamicObject) -> DynamicObject {
        DynamicObject::Translate {
            offset,
            child: Box::new(child),
        }
    }

    fn leaves(obj: &DynamicObject) -> CsgLeaves {
        CsgLeaves::from_dynamic(obj, DAffine3::IDENTITY, None)
    }

    /// Plate with a through-hole: leaf 0 = plate cube, leaf 1 = hole cylinder.
    fn plate_with_hole() -> DynamicObject {
        DynamicObject::Difference {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate([2.0, -0.5, 2.0], cylinder(0.5, 2.0)),
            ],
        }
    }

    /// Boss on a plate: leaf 0 = plate, leaf 1 = boss cylinder on top.
    fn boss_on_plate() -> DynamicObject {
        DynamicObject::Union {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate([2.0, 1.0, 2.0], cylinder(0.5, 1.0)),
            ],
        }
    }

    #[test]
    fn cube_adjacent_faces_resolve_convex_lines() {
        let csg = leaves(&cube([1.0, 2.0, 3.0]));
        let adjacent = [
            (FaceTag::CubePosX, FaceTag::CubePosY),
            (FaceTag::CubePosX, FaceTag::CubeNegZ),
            (FaceTag::CubeNegY, FaceTag::CubeNegZ),
        ];
        for (fa, fb) in adjacent {
            let (geom, convexity) =
                resolve_edge(&csg, (0, fa), (0, fb)).expect("adjacent faces share an edge");
            assert_eq!(convexity, Convexity::Convex, "{fa:?}/{fb:?}");
            let EdgeGeom::Straight(e) = geom else {
                panic!("cube edge must be straight");
            };
            // Endpoints on both face planes and a full-length run.
            let expect_len = match (fa, fb) {
                (FaceTag::CubePosX, FaceTag::CubePosY) => 3.0,
                (FaceTag::CubePosX, FaceTag::CubeNegZ) => 2.0,
                _ => 1.0,
            };
            assert!(((e.end - e.start).length() - expect_len).abs() < 1e-9);
        }
        // Opposite faces are parallel: no edge.
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CubePosX), (0, FaceTag::CubeNegX)).unwrap_err(),
            EdgeReject::NotAnEdge
        );
    }

    #[test]
    fn cube_edge_directions_point_into_faces() {
        let csg = leaves(&cube([1.0, 1.0, 1.0]));
        let (geom, _) = resolve_edge(&csg, (0, FaceTag::CubePosX), (0, FaceTag::CubePosY)).unwrap();
        let EdgeGeom::Straight(e) = geom else {
            panic!("straight");
        };
        assert!((e.na - DVec3::X).length() < 1e-9);
        assert!((e.nb - DVec3::Y).length() < 1e-9);
        // Into face A (x = 1 plane): downward; into face B (y = 1): -x.
        assert!(e.da.dot(DVec3::Y) < -0.99);
        assert!(e.db.dot(DVec3::X) < -0.99);
        assert!((e.extent_a - 1.0).abs() < 1e-9);
        assert!((e.extent_b - 1.0).abs() < 1e-9);
    }

    #[test]
    fn same_leaf_cylinder_rims() {
        let csg = leaves(&cylinder(1.0, 2.0));
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("rim must be circular");
        };
        assert!((c.center - DVec3::new(0.0, 2.0, 0.0)).length() < 1e-9);
        assert!((c.axis - DVec3::Y).length() < 1e-9, "axis = cap outward");
        assert!((c.radius - 1.0).abs() < 1e-9);
        // Plane face setback runs inward; wall setback runs down.
        assert!(c.da2.x < -0.99);
        assert!(c.db2.y < -0.99);

        let (_, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylBottom)).unwrap();
        assert_eq!(convexity, Convexity::Convex);

        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CylTop), (0, FaceTag::CylBottom)).unwrap_err(),
            EdgeReject::NotAnEdge
        );
    }

    #[test]
    fn cone_base_rim_resolves() {
        let csg = leaves(&DynamicObject::Cone {
            radius: 1.0,
            height: 2.0,
            segments: 32,
        });
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::ConeSide), (0, FaceTag::ConeBottom)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.radius - 1.0).abs() < 1e-9);
        assert!(
            (c.axis + DVec3::Y).length() < 1e-9,
            "axis = base outward (-Y)"
        );
        // Wall setback climbs the slant: radially inward and axially away
        // from the base plane in meridian coords.
        assert!(c.db2.x < 0.0 && c.db2.y < 0.0);
    }

    #[test]
    fn hole_rim_is_convex_circle_with_outward_plane_dir() {
        let obj = plate_with_hole();
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CylSide)).unwrap();
        assert_eq!(convexity, Convexity::Convex, "hole rims are chamferable");
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.center - DVec3::new(2.0, 1.0, 2.0)).length() < 1e-9);
        assert!((c.radius - 0.5).abs() < 1e-9);
        assert!((c.axis - DVec3::Y).length() < 1e-9);
        // The plate surface survives *outside* the hole; the wall runs down.
        assert!(c.da2.x > 0.99);
        assert!(c.db2.y < -0.99);
    }

    #[test]
    fn boss_base_ring_is_concave() {
        let csg = leaves(&boss_on_plate());
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CylSide)).unwrap();
        assert_eq!(
            convexity,
            Convexity::Concave,
            "boss-base ring is a fillet case"
        );
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!(c.da2.x > 0.99, "plate survives outside the boss");
        assert!(c.db2.y > 0.99, "boss wall survives upward");
    }

    #[test]
    fn straight_concave_rejected() {
        // L-bracket: two cubes sharing a concave inner edge along Z at x=1,y=1.
        let obj = DynamicObject::Union {
            children: vec![
                cube([2.0, 1.0, 2.0]),
                translate([0.0, 1.0, 0.0], cube([1.0, 1.0, 2.0])),
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CubePosX)).unwrap_err(),
            EdgeReject::StraightConcave
        );
    }

    #[test]
    fn oblique_and_curved_pairs_rejected() {
        // Tilted cylinder through a cube top: oblique plane∩cylinder.
        let obj = DynamicObject::Union {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate(
                    [2.0, 0.0, 2.0],
                    DynamicObject::Rotate {
                        angles: [30.0, 0.0, 0.0],
                        child: Box::new(cylinder(0.5, 3.0)),
                    },
                ),
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CylSide)).unwrap_err(),
            EdgeReject::UnsupportedCurve
        );
        // Cylinder∩cylinder is never a supported curve.
        let obj = DynamicObject::Union {
            children: vec![
                cylinder(1.0, 2.0),
                translate([1.0, 0.0, 0.0], cylinder(1.0, 2.0)),
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CylSide), (1, FaceTag::CylSide)).unwrap_err(),
            EdgeReject::UnsupportedCurve
        );
    }

    #[test]
    fn mismatched_anisotropic_frames_reject() {
        // A boss cylinder standing on a squashed cylinder's cap: the cap's
        // elliptical extent needs leaf A's frame, the boss wall stays round
        // only in leaf B's — no single frame fits both.
        let obj = DynamicObject::Union {
            children: vec![
                DynamicObject::Scale {
                    factors: [2.0, 1.0, 1.0],
                    child: Box::new(cylinder(1.0, 2.0)),
                },
                translate([0.0, 2.0, 0.0], cylinder(0.25, 1.0)),
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CylTop), (1, FaceTag::CylSide)).unwrap_err(),
            EdgeReject::NonUniformScale
        );
        // Cube faces stay planes under any affine: edges still resolve.
        let obj = DynamicObject::Scale {
            factors: [2.0, 1.0, 1.0],
            child: Box::new(cube([1.0, 1.0, 1.0])),
        };
        let csg = leaves(&obj);
        assert!(resolve_edge(&csg, (0, FaceTag::CubePosX), (0, FaceTag::CubePosY)).is_ok());
    }

    #[test]
    fn elliptical_rim_resolves_in_leaf_frame() {
        // [2,1,1] squashes the cross-section: the rim is a world ellipse,
        // resolved as a unit circle in the leaf's frame.
        let obj = DynamicObject::Scale {
            factors: [2.0, 1.0, 1.0],
            child: Box::new(cylinder(1.0, 2.0)),
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert_ne!(c.frame, DAffine3::IDENTITY);
        assert!((c.radius - 1.0).abs() < 1e-9, "frame-local unit radius");
        // The world polyline traces the ellipse x²/4 + z² = 1 at y = 2.
        for p in geom.polyline(24) {
            assert!((p.y - 2.0).abs() < 1e-9);
            assert!((p.x * p.x / 4.0 + p.z * p.z - 1.0).abs() < 1e-9);
        }
        // World distances follow the ellipse (near the semi-major vertex).
        let d = geom.distance_to(DVec3::new(2.1, 2.0, 0.0));
        assert!((d - 0.1).abs() < 0.01, "distance to ellipse vertex: {d}");
    }

    #[test]
    fn spheroid_dome_rim_resolves() {
        // A [1,2,1]-scaled sphere sunk into a plate: an ellipsoid in world,
        // a unit sphere in its own frame — the dome ring resolves there.
        let obj = DynamicObject::Union {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate(
                    [2.0, 0.5, 2.0],
                    DynamicObject::Scale {
                        factors: [1.0, 2.0, 1.0],
                        child: Box::new(DynamicObject::Sphere {
                            radius: 1.0,
                            segments: 32,
                            stacks: 16,
                        }),
                    },
                ),
            ],
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::Sphere)).unwrap();
        assert_eq!(convexity, Convexity::Concave, "dome meets plate concavely");
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        // Plate top y = 1 is local y = 0.25 in the sphere frame: the ring
        // radius is √(1 − 0.25²), unscaled in world x/z.
        let expect = (1.0f64 - 0.0625).sqrt();
        assert!((c.radius - expect).abs() < 1e-9);
        for p in geom.polyline(16) {
            assert!((p.y - 1.0).abs() < 1e-9);
            let rho = DVec2::new(p.x - 2.0, p.z - 2.0).length();
            assert!((rho - expect).abs() < 1e-9);
        }
    }

    #[test]
    fn cap_chord_straight_edge_maps_from_frame() {
        // Squashed cylinder with its x > 0 half cut away: the cap meets the
        // cut wall in a straight chord — resolved in the leaf frame (the
        // cap's extent is elliptical in world) and mapped back exactly.
        let obj = DynamicObject::Difference {
            children: vec![
                DynamicObject::Scale {
                    factors: [2.0, 1.0, 1.0],
                    child: Box::new(cylinder(1.0, 2.0)),
                },
                translate([0.0, -1.0, -3.0], cube([3.0, 4.0, 6.0])),
            ],
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylTop), (1, FaceTag::CubeNegX)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Straight(e) = geom else {
            panic!("straight chord");
        };
        // Chord x = 0, y = 2, z ∈ [−1, 1] (the ellipse's minor vertices).
        assert!(e.start.x.abs() < 1e-6 && e.end.x.abs() < 1e-6);
        assert!((e.start.y - 2.0).abs() < 1e-6 && (e.end.y - 2.0).abs() < 1e-6);
        assert!(((e.end - e.start).length() - 2.0).abs() < 1e-3);
        assert!((e.na - DVec3::Y).length() < 1e-6, "cap normal up");
        assert!((e.nb - DVec3::X).length() < 1e-6, "cut wall normal +x");
        // In-face directions: away from the cut into the cap (−x), and up
        // isn't possible on the wall — material survives below? No: the wall
        // is the hole side; its surviving in-face direction runs down (−y).
        assert!(e.da.dot(DVec3::X) < -0.99);
        assert!(e.db.dot(DVec3::Y) < -0.99);
        assert!(e.extent_a > 0.5 && e.extent_b > 0.5);
    }

    #[test]
    fn axisymmetric_scale_keeps_rims_circular() {
        // Stretched taller ([1,3,1]): cross-section untouched, rim stays a
        // unit circle at the scaled height.
        let obj = DynamicObject::Scale {
            factors: [1.0, 3.0, 1.0],
            child: Box::new(cylinder(1.0, 2.0)),
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.center - DVec3::new(0.0, 6.0, 0.0)).length() < 1e-9);
        assert!((c.radius - 1.0).abs() < 1e-9);

        // Fattened uniformly across the cross-section ([2,1,2]).
        let obj = DynamicObject::Scale {
            factors: [2.0, 1.0, 2.0],
            child: Box::new(cylinder(1.0, 2.0)),
        };
        let csg = leaves(&obj);
        let (geom, _) = resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylBottom)).unwrap();
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!(c.center.length() < 1e-9);
        assert!((c.radius - 2.0).abs() < 1e-9);

        // Cone base rim under [2,1,2]: radius doubles, slant changes, still
        // a circle.
        let obj = DynamicObject::Scale {
            factors: [2.0, 1.0, 2.0],
            child: Box::new(DynamicObject::Cone {
                radius: 1.0,
                height: 2.0,
                segments: 32,
            }),
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::ConeSide), (0, FaceTag::ConeBottom)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.radius - 2.0).abs() < 1e-9);
    }

    #[test]
    fn axisymmetric_scale_composes_with_rotation() {
        // World Scale[2,2,1] over Rotate(x, 90°): the cylinder axis lands on
        // Z (unit world stretch) and the cross-section spans world X/Y
        // (scaled 2×) — axisymmetric about the *leaf* axis even though the
        // scale wrapper alone is not.
        let obj = DynamicObject::Scale {
            factors: [2.0, 2.0, 1.0],
            child: Box::new(DynamicObject::Rotate {
                angles: [90.0, 0.0, 0.0],
                child: Box::new(cylinder(1.0, 2.0)),
            }),
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.radius - 2.0).abs() < 1e-9);
        assert!(c.axis.cross(DVec3::Z).length() < 1e-9);
        assert!((c.center - DVec3::new(0.0, 0.0, 2.0)).length() < 1e-9);

        // The same wrapper pair with the scale squashing the cross-section
        // ([1,2,1] over the rotated cylinder → world X vs Y anisotropy):
        // resolves through the leaf frame as a world ellipse at z = 2.
        let obj = DynamicObject::Scale {
            factors: [1.0, 2.0, 1.0],
            child: Box::new(DynamicObject::Rotate {
                angles: [90.0, 0.0, 0.0],
                child: Box::new(cylinder(1.0, 2.0)),
            }),
        };
        let csg = leaves(&obj);
        let (geom, _) = resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert_ne!(c.frame, DAffine3::IDENTITY);
        for p in geom.polyline(16) {
            assert!((p.z - 2.0).abs() < 1e-9);
            assert!((p.x * p.x + p.y * p.y / 4.0 - 1.0).abs() < 1e-9);
        }
    }

    #[test]
    fn scaled_hole_rim_resolves() {
        // Plate with a hole whose cylinder is stretched deeper: the rim on
        // the plate top is still an exact circle (cross-leaf pair).
        let obj = DynamicObject::Difference {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate(
                    [2.0, -0.5, 2.0],
                    DynamicObject::Scale {
                        factors: [1.0, 2.0, 1.0],
                        child: Box::new(cylinder(0.5, 1.0)),
                    },
                ),
            ],
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CylSide)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        assert!((c.center - DVec3::new(2.0, 1.0, 2.0)).length() < 1e-9);
        assert!((c.radius - 0.5).abs() < 1e-9);
    }

    #[test]
    fn partial_circle_rejected_when_hole_leaves_plate() {
        // Hole centered on the plate boundary: half the rim is air.
        let obj = DynamicObject::Difference {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate([4.0, -0.5, 2.0], cylinder(0.5, 2.0)),
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::CylSide)).unwrap_err(),
            EdgeReject::PartialCircle
        );
    }

    #[test]
    fn consumed_edge_rejected_and_partial_edge_trimmed() {
        // A slab that shaves the top of the cube: the old top edges are gone.
        let consumed = DynamicObject::Difference {
            children: vec![
                cube([2.0, 2.0, 2.0]),
                translate([-1.0, 1.5, -1.0], cube([4.0, 2.0, 4.0])),
            ],
        };
        let csg = leaves(&consumed);
        assert!(
            resolve_edge(&csg, (0, FaceTag::CubePosY), (0, FaceTag::CubePosX)).is_err(),
            "an edge shaved off entirely must not resolve"
        );

        // A notch eating half the edge: the surviving run is trimmed.
        let notched = DynamicObject::Difference {
            children: vec![
                cube([2.0, 2.0, 2.0]),
                translate([1.0, 1.0, 1.0], cube([2.0, 2.0, 2.0])),
            ],
        };
        let csg = leaves(&notched);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosX), (0, FaceTag::CubePosY)).unwrap();
        assert_eq!(convexity, Convexity::Convex);
        let EdgeGeom::Straight(e) = geom else {
            panic!("straight");
        };
        // Original edge spans z ∈ [0, 2] at (2, 2, z); the notch removes z > 1.
        let len = (e.end - e.start).length();
        assert!(
            (len - 1.0).abs() < 0.05,
            "surviving run ≈ half the edge, got {len}"
        );
        let z_max = e.start.z.max(e.end.z);
        assert!(z_max < 1.05, "run must stay on the surviving side");
    }

    #[test]
    fn sphere_plane_circle_resolves() {
        // Sphere sunk halfway into a plate top.
        let obj = DynamicObject::Union {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate(
                    [2.0, 0.5, 2.0],
                    DynamicObject::Sphere {
                        radius: 1.0,
                        segments: 32,
                        stacks: 16,
                    },
                ),
            ],
        };
        let csg = leaves(&obj);
        let (geom, convexity) =
            resolve_edge(&csg, (0, FaceTag::CubePosY), (1, FaceTag::Sphere)).unwrap();
        assert_eq!(convexity, Convexity::Concave, "dome meets plate concavely");
        let EdgeGeom::Circular(c) = geom else {
            panic!("circular");
        };
        // Circle at plate-top height: radius √(1 − 0.5²).
        assert!((c.radius - (1.0f64 - 0.25).sqrt()).abs() < 1e-9);
        assert!((c.center - DVec3::new(2.0, 1.0, 2.0)).length() < 1e-9);
    }

    #[test]
    fn leaf_flipped_tracks_difference_nesting() {
        // a − (b − c): b flipped once, c twice (back to unflipped).
        let obj = DynamicObject::Difference {
            children: vec![
                cube([4.0, 4.0, 4.0]),
                DynamicObject::Difference {
                    children: vec![
                        translate([1.0, 1.0, 1.0], cube([2.0, 2.0, 2.0])),
                        translate([1.5, 1.5, 1.5], cube([1.0, 1.0, 1.0])),
                    ],
                },
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(leaf_flipped(&csg), vec![false, true, false]);
    }

    #[test]
    fn distance_and_polyline_helpers() {
        let csg = leaves(&cylinder(1.0, 2.0));
        let (geom, _) = resolve_edge(&csg, (0, FaceTag::CylSide), (0, FaceTag::CylTop)).unwrap();
        // A point just outside the rim.
        assert!((geom.distance_to(DVec3::new(1.1, 2.0, 0.0)) - 0.1).abs() < 1e-9);
        let poly = geom.polyline(16);
        assert_eq!(poly.len(), 16);
        for p in poly {
            assert!((p.y - 2.0).abs() < 1e-9);
            assert!((DVec3::new(p.x, 0.0, p.z).length() - 1.0).abs() < 1e-9);
        }
    }
}