rscad_core/feature/

lower.rs

//! Lower an [`EdgeFeatureData`] into pure-CSG masks against its child.
//!
//! Chamfers and fillets are not primitives anywhere in the pipeline — an
//! `EdgeFeature` node *lowers* to booleans over ordinary
//! `Cube`/`Cylinder`/`Cone`/`Revolve` solids, so every backend (Manifold,
//! BSP, SDF ray-march, projection) renders them through variants it already
//! supports. Subtractive masks cut as `Difference(child, Union(sub))`;
//! additive masks (concave ring fillets) union on top of that.
//!
//! Edges are **re-resolved against the current child on every lowering**
//! ([`resolve_edge`]), so features follow parameter edits; selections that no
//! longer resolve are skipped (callers surface the count from
//! [`LowerReport`]). Masks are built in the child's local frame — callers
//! apply outer transforms above the feature node as usual.
//!
//! ## Chamfer masks
//!
//! - **Straight convex edge**: the cut plane passes through the two setback
//!   lines `P_F = edge + s·d̂_F`; the mask is an oriented box with one face on
//!   that plane, deep enough to swallow the edge, emitted as
//!   `Translate ∘ Rotate ∘ Cube`. Extending past the setback lines is safe —
//!   beyond them the box sits on the air side of both faces (convexity) — and
//!   the tiny ε overhang along the edge keeps corners watertight (adjacent
//!   chamfers stay un-mitered in v1).
//! - **Circular convex rim**: the meridian setback points revolve into a cut
//!   cone (a cylinder when the setbacks share a radius). The mask is
//!   `Difference(slab, cone)` when the rim's material lies radially inside
//!   (outer rims) or `Intersection(slab, cone)` when it lies outside
//!   (hole-rim countersinks). The slab is clamped axially to the setback
//!   band, so blind holes are never deepened and walls below the chamfer are
//!   untouched; radial/axial ε only ever extends into air.
//!
//! ## Fillet masks
//!
//! For faces meeting at interior angle `γ` (`cos γ = d̂_A·d̂_B`), a radius-`r`
//! arc tangent to both faces has its center a distance `h = r/sin(γ/2)` down
//! the bisector and touches each face at setback `s_T = r/tan(γ/2)` from the
//! edge. The clamp below applies to `s_T` (that is what must fit on the
//! faces), and `r` rescales with it.
//!
//! - **Straight convex edge**: the chamfer box through the `s_T` setbacks,
//!   minus the tangent cylinder along the edge — the cylinder wall becomes
//!   the fillet surface. Fully SDF-visible.
//! - **Circular rim** (convex *or* concave): the meridian wedge between the
//!   setbacks minus the tangent circle, sampled into a polygon
//!   (`data.segments` arc chords) and revolved around the rim axis. Convex
//!   rims subtract the ring; concave rims (boss-base junctions) *add* the
//!   corner-fill ring. `Revolve` masks are invisible in the SDF preview
//!   (existing sentinel) and ignored by the analytic picking walk — mesh
//!   backends render them exactly.
//!
//! Feature size is clamped per edge to `0.45 ×` the smaller adjacent face
//! extent, so a slider can't push a mask through the far side of a face.
//!
//! Rim masks tessellate at the **curved leaf's own segment count and
//! angular phase** ([`rim_segments`]/[`place_on_rim`]): where a mask surface
//! runs along the leaf's wall (countersink tangency, fillet tangent bands),
//! coincident chords make the boolean junction a clean polygon ring instead
//! of an interleaved sliver band (wobbly seams, streaky shading).
//!
//! ## Frame-mapped (elliptical) rims
//!
//! Rims of anisotropically scaled leaves resolve as circles in the leaf's
//! own frame ([`CircularEdge::frame`](crate::feature::edge::CircularEdge)).
//! Their masks are built in that frame by the same meridian constructions
//! and wrapped with the frame transform (`Translate∘Rotate∘Scale∘Rotate`
//! from an SVD split), so the cut is the exact affine image of the circular
//! one: the feature stretches with the shape, and `size`/extents are
//! frame-local for those edges. A constant-width cut along an ellipse is
//! not expressible with the primitive vocabulary (its offset curves are
//! not conics), so stretching is the exact-CSG semantic choice.

use std::collections::BTreeMap;

use crate::feature::edge::{Convexity, EdgeGeom, resolve_edge};
use crate::feature::geom::CsgLeaves;
use crate::feature::{EdgeFeatureData, EdgeFeatureKind, LeafPath};
use crate::prelude_::*;

/// Watertightness margin as a fraction of the (clamped) feature size.
const EPS_FRAC: f64 = 1.0e-3;
/// Fraction of the smaller adjacent face extent a feature may consume.
const MAX_EXTENT_FRAC: f64 = 0.45;
/// Tessellation of revolved/round mask surfaces that are *not* user-tunable
/// (cut cones, slabs, ring revolution) — matches typical body tessellation.
const MASK_SEGMENTS: usize = 64;

/// What a lowering pass did — callers log/report the skip count.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct LowerReport {
    pub applied: usize,
    pub skipped: usize,
}

/// Masks for one lowering pass, in the child's frame: `sub` cuts material
/// (chamfers, convex fillets), `add` fills it in (concave ring fillets).
#[derive(Clone, Debug, Default)]
pub struct LoweredMasks {
    pub sub: Vec<DynamicObject>,
    pub add: Vec<DynamicObject>,
}

/// A single edge's mask with its boolean role.
enum Mask {
    Sub(DynamicObject),
    Add(DynamicObject),
}

/// Lower `data` against `child`: subtract the cut masks, union the fill
/// masks, or pass the child through unchanged when nothing resolves.
pub fn lower(data: &EdgeFeatureData, child: &DynamicObject) -> (DynamicObject, LowerReport) {
    let (masks, report) = lower_masks(data, child);
    let mut out = child.clone();
    if !masks.sub.is_empty() {
        out = DynamicObject::Difference {
            children: vec![
                out,
                DynamicObject::Union {
                    children: masks.sub,
                },
            ],
        };
    }
    if !masks.add.is_empty() {
        let children = core::iter::once(out).chain(masks.add).collect();
        out = DynamicObject::Union { children };
    }
    (out, report)
}

/// The masks for `data` against `child`, in the child's frame.
/// Unresolvable/unsupported edges are skipped and counted.
///
/// Edges resolve against the **unfeatured** child (inner chamfers/fillets
/// stripped): a sibling feature's cut must not trim a stored edge's
/// surviving run, or the mask stops short of shared corners and a sliver of
/// the original sharp edge pokes through where mixed-kind features meet.
/// Full-length masks simply overlap at corners, matching the same-kind
/// merge behavior.
pub fn lower_masks(data: &EdgeFeatureData, child: &DynamicObject) -> (LoweredMasks, LowerReport) {
    let csg = CsgLeaves::from_dynamic_unfeatured(child, DAffine3::IDENTITY, None);
    let index: BTreeMap<&LeafPath, usize> =
        csg.paths.iter().enumerate().map(|(i, p)| (p, i)).collect();
    let mut masks = LoweredMasks::default();
    let mut report = LowerReport::default();
    for sel in &data.edges {
        let resolved = index
            .get(&sel.a.path)
            .zip(index.get(&sel.b.path))
            .ok_or(crate::feature::edge::EdgeReject::Unsupported)
            .and_then(|(&la, &lb)| resolve_edge(&csg, (la, sel.a.face), (lb, sel.b.face)));
        let mask = match resolved {
            Ok((geom, convexity)) => edge_mask(data, &geom, convexity),
            Err(_) => None,
        };
        match mask {
            Some(Mask::Sub(m)) => {
                masks.sub.push(m);
                report.applied += 1;
            }
            Some(Mask::Add(m)) => {
                masks.add.push(m);
                report.applied += 1;
            }
            None => report.skipped += 1,
        }
    }
    (masks, report)
}

/// Mask for one resolved edge, or `None` when this kind/shape combination is
/// not lowerable (concave chamfers are undefined; straight concave creases
/// are rejected at resolve time and never reach here).
fn edge_mask(data: &EdgeFeatureData, geom: &EdgeGeom, convexity: Convexity) -> Option<Mask> {
    match (data.kind, convexity, geom) {
        (EdgeFeatureKind::Chamfer, Convexity::Convex, EdgeGeom::Straight(e)) => {
            let s = clamp_size(data.size, e.extent_a, e.extent_b)?;
            Some(Mask::Sub(straight_chamfer_mask(e, s)))
        }
        (EdgeFeatureKind::Chamfer, Convexity::Convex, EdgeGeom::Circular(c)) => {
            let s = clamp_size(data.size, c.extent_a, c.extent_b)?;
            circular_chamfer_mask(c, s).map(Mask::Sub)
        }
        (EdgeFeatureKind::Fillet, Convexity::Convex, EdgeGeom::Straight(e)) => {
            straight_fillet_mask(e, data.size, data.segments).map(Mask::Sub)
        }
        (EdgeFeatureKind::Fillet, convexity, EdgeGeom::Circular(c)) => {
            circular_fillet_mask(c, data.size, data.segments, convexity).map(match convexity {
                Convexity::Convex => Mask::Sub,
                Convexity::Concave => Mask::Add,
            })
        }
        _ => None,
    }
}

/// Per-edge size clamp: at most [`MAX_EXTENT_FRAC`] of the smaller adjacent
/// face extent. `None` when there is no room at all.
fn clamp_size(size: f64, extent_a: f64, extent_b: f64) -> Option<f64> {
    let cap = MAX_EXTENT_FRAC * extent_a.min(extent_b);
    let s = size.min(cap);
    (s > 1.0e-9).then_some(s)
}

// ── Straight chamfer ─────────────────────────────────────────────

/// Oriented box mask cutting the corner between the setback lines.
fn straight_chamfer_mask(e: &crate::feature::edge::StraightEdge, s: f64) -> DynamicObject {
    let eps = EPS_FRAC * s;
    let u_hat = e.dir();
    let len = (e.end - e.start).length();
    let pa = e.start + s * e.da;
    let pb = e.start + s * e.db;

    // Frame: x = edge, y = across the chamfer strip, z = out of the cut
    // plane toward the removed corner (air side, `z·(na+nb) > 0`).
    let mut y_hat = (pb - pa).normalize();
    let mut z_hat = u_hat.cross(y_hat);
    if z_hat.dot(e.na + e.nb) < 0.0 {
        y_hat = -y_hat;
        z_hat = -z_hat;
    }
    let (from, width) = if y_hat.dot(pb - pa) >= 0.0 {
        (pa, (pb - pa).length())
    } else {
        (pb, (pb - pa).length())
    };
    let depth = (e.start - pa).dot(z_hat);

    let corner = from - eps * u_hat - eps * y_hat;
    oriented_cube(
        corner,
        DMat3::from_cols(u_hat, y_hat, z_hat),
        DVec3::new(len + 2.0 * eps, width + 2.0 * eps, depth + eps),
    )
}

/// `Translate ∘ Rotate ∘ Cube` for a corner-anchored box in the frame given
/// by rotation matrix `r` (orthonormal, right-handed).
fn oriented_cube(corner: DVec3, r: DMat3, size: DVec3) -> DynamicObject {
    DynamicObject::Translate {
        offset: corner.to_array(),
        child: Box::new(DynamicObject::Rotate {
            angles: euler_zyx_degrees(r),
            child: Box::new(DynamicObject::Cube {
                size: size.to_array(),
            }),
        }),
    }
}

/// Euler angles (degrees, `[x, y, z]` storage) reproducing `r` under the
/// renderer's `Rz·Ry·Rx` convention.
fn euler_zyx_degrees(r: DMat3) -> [f64; 3] {
    let (z, y, x) = DQuat::from_mat3(&r).to_euler(EulerRot::ZYX);
    [x.to_degrees(), y.to_degrees(), z.to_degrees()]
}

// ── Circular chamfer ─────────────────────────────────────────────

/// Revolved cut mask for a convex rim. Meridian frame: x = outward radial,
/// y = rim axis (the planar face's outward normal). Built in the rim's
/// [`frame`](crate::feature::edge::CircularEdge::frame) and wrapped with it,
/// so elliptical rims (anisotropic leaves) get the frame image of the
/// circular cut — the feature stretches with the shape, and `s` is
/// frame-local there.
fn circular_chamfer_mask(c: &crate::feature::edge::CircularEdge, s: f64) -> Option<DynamicObject> {
    let eps = EPS_FRAC * s;
    // Setback points in the meridian plane; the rim sits at (radius, 0).
    let pa = DVec2::new(c.radius, 0.0) + s * c.da2;
    let pb = DVec2::new(c.radius, 0.0) + s * c.db2;
    if pa.x <= eps || pb.x <= eps {
        return None; // setback crosses the axis — no ring left to cut
    }

    // The removed wedge sits radially inside the cut for hole rims
    // (plane-face material outside the rim) and outside it for outer rims.
    let removed_inside = c.da2.x > 0.0;
    let (y_lo, y_hi) = (pa.y.min(pb.y), pa.y.max(pb.y));
    let slab_r = pa.x.max(pb.x).max(c.radius) + eps;
    // ε extends only into air: above the planar face (+y) and radially; the
    // low end stays exactly at the setback so walls/holes below are safe.
    let slab = axial_cylinder(c, slab_r, y_lo, y_hi + eps);
    let cut = cut_cone(c, pa, pb, y_lo, y_hi + eps)?;

    let children = vec![slab, cut];
    let mask = if removed_inside {
        DynamicObject::Intersection { children }
    } else {
        DynamicObject::Difference { children }
    };
    Some(wrap_frame(&c.frame, mask))
}

/// Angular resolution for a rim's masks: the curved leaf's own tessellation
/// when known (chords then coincide in the mesh booleans — no sliver fins),
/// else [`MASK_SEGMENTS`].
fn rim_segments(c: &crate::feature::edge::CircularEdge) -> usize {
    c.leaf_segments.unwrap_or(MASK_SEGMENTS).max(3)
}

/// `Translate ∘ Rotate` for rim masks: local `+Y` along `up_sign · axis`,
/// local `+X` along the leaf's tessellation origin (`ref_dir`), placed at
/// `center + y_offset · axis` — mask polygon vertices land on the leaf's.
fn place_on_rim(
    c: &crate::feature::edge::CircularEdge,
    y_offset: f64,
    up_sign: f64,
    obj: DynamicObject,
) -> DynamicObject {
    let y = up_sign * c.axis;
    let x = c.ref_dir;
    let r = DMat3::from_cols(x, y, x.cross(y));
    DynamicObject::Translate {
        offset: (c.center + y_offset * c.axis).to_array(),
        child: Box::new(DynamicObject::Rotate {
            angles: euler_zyx_degrees(r),
            child: Box::new(obj),
        }),
    }
}

/// Solid of revolution whose meridian boundary passes through `pa` and `pb`
/// (a cone, or a cylinder when the radii match), covering at least the
/// axial slab band `[y_lo, y_hi]` so the boolean against the slab is exact.
fn cut_cone(
    c: &crate::feature::edge::CircularEdge,
    pa: DVec2,
    pb: DVec2,
    y_lo: f64,
    y_hi: f64,
) -> Option<DynamicObject> {
    let d = pb - pa;
    if d.x.abs() < 1.0e-12 {
        // Vertical cut: a plain cylinder through both points.
        let pad = y_hi - y_lo;
        return Some(axial_cylinder(c, pa.x, y_lo - pad, y_hi + pad));
    }
    if d.y.abs() < 1.0e-12 {
        return None; // horizontal cut can't bound a corner wedge
    }
    let slope = d.y / d.x; // dy per dρ, nonzero both ways
    let y_apex = pa.y - pa.x * slope;
    // The setback segment never crosses the axis, so the apex lies outside
    // its own y-band; but it must also clear the (ε-padded) slab band for
    // the cone to bound the cut across the whole slab.
    let sigma = if y_apex >= y_hi {
        1.0 // apex above: cone narrows upward, build apex-up directly
    } else if y_apex <= y_lo {
        -1.0 // apex below: build apex-up in the flipped frame ỹ = −y
    } else {
        return None; // ε-sliver degeneracy
    };
    // Build coords ỹ = σ·y: apex at apex_t, base beyond the far slab end.
    let apex_t = sigma * y_apex;
    let far_t = sigma * if sigma > 0.0 { y_lo } else { y_hi };
    let base_t = far_t - (y_hi - y_lo).max(1.0e-9);
    let height = apex_t - base_t;
    if height <= 0.0 {
        return None;
    }
    // Radius at the build-frame base, from the meridian line through pa.
    let radius = pa.x * height / (apex_t - sigma * pa.y);
    if !radius.is_finite() || radius <= 0.0 {
        return None;
    }
    // Base sits at rim-frame y = σ·base_t, apex points along σ·axis.
    Some(place_on_rim(
        c,
        sigma * base_t,
        sigma,
        DynamicObject::Cone {
            radius,
            height,
            segments: rim_segments(c),
        },
    ))
}

/// Cylinder around the rim axis spanning rim-frame `y ∈ [y0, y1]`.
fn axial_cylinder(
    c: &crate::feature::edge::CircularEdge,
    radius: f64,
    y0: f64,
    y1: f64,
) -> DynamicObject {
    place_on_rim(
        c,
        y0,
        1.0,
        DynamicObject::Cylinder {
            radius,
            height: (y1 - y0).max(1.0e-9),
            segments: rim_segments(c),
        },
    )
}

// ── Fillets ──────────────────────────────────────────────────────

/// Fillet frame for faces meeting at interior angle `γ` (`cos γ = d̂_A·d̂_B`):
/// the tangency setback `s_T = r/tan(γ/2)` (clamped to the face extents, with
/// the radius rescaled to match) and the arc-center distance `h = r/sin(γ/2)`
/// along the bisector. `None` for flush faces, knife edges, or no room.
fn fillet_frame(
    size: f64,
    cos_gamma: f64,
    extent_a: f64,
    extent_b: f64,
) -> Option<(f64, f64, f64)> {
    let sin_half = ((1.0 - cos_gamma) / 2.0).max(0.0).sqrt();
    let cos_half = ((1.0 + cos_gamma) / 2.0).max(0.0).sqrt();
    if sin_half < 1.0e-6 || cos_half < 1.0e-6 {
        return None;
    }
    let setback = clamp_size(size * cos_half / sin_half, extent_a, extent_b)?;
    let radius = setback * sin_half / cos_half;
    Some((radius, setback, radius / sin_half))
}

/// Box-minus-cylinder mask rounding a straight convex edge: the chamfer box
/// through the tangency setbacks, minus the tangent cylinder whose wall
/// becomes the fillet surface.
fn straight_fillet_mask(
    e: &crate::feature::edge::StraightEdge,
    size: f64,
    segments: usize,
) -> Option<DynamicObject> {
    let (r, setback, h) = fillet_frame(size, e.da.dot(e.db), e.extent_a, e.extent_b)?;
    let cut_box = straight_chamfer_mask(e, setback);
    let eps = EPS_FRAC * setback;
    let u_hat = e.dir();
    let len = (e.end - e.start).length();
    // Tangent cylinder along the edge, overshooting the box on both ends so
    // the subtraction is watertight.
    let axis = e.start + h * (e.da + e.db).normalize();
    let round = place_on_axis(
        axis - 2.0 * eps * u_hat,
        u_hat,
        DynamicObject::Cylinder {
            radius: r,
            height: len + 4.0 * eps,
            segments: segments.max(3),
        },
    );
    Some(DynamicObject::Difference {
        children: vec![cut_box, round],
    })
}

/// Revolved fillet mask for a rim: the meridian wedge between the tangency
/// setbacks minus the tangent circle, swept around the rim axis. Convex rims
/// subtract the ring (rounding the corner off); concave rims (boss-base
/// junctions) return the corner-fill ring for the caller to union in.
fn circular_fillet_mask(
    c: &crate::feature::edge::CircularEdge,
    size: f64,
    segments: usize,
    convexity: Convexity,
) -> Option<DynamicObject> {
    let (r, setback, h) = fillet_frame(size, c.da2.dot(c.db2), c.extent_a, c.extent_b)?;
    let eps = EPS_FRAC * setback;
    let e2 = DVec2::new(c.radius, 0.0);
    let pa = e2 + setback * c.da2;
    let pb = e2 + setback * c.db2;
    if pa.x <= eps || pb.x <= eps {
        return None; // setback crosses the axis — no ring left
    }
    let center = e2 + h * (c.da2 + c.db2).normalize();

    // ε pushes the wall segments across the faces: into air for a convex
    // corner cut, into material for a concave corner fill. The planar face's
    // meridian outward normal is +y by construction.
    let out = match convexity {
        Convexity::Convex => 1.0,
        Convexity::Concave => -1.0,
    };
    let na = out * DVec2::Y;
    let nb = out * c.nb2;

    let mut points = vec![pa + eps * na, e2 + eps * (na + nb), pb + eps * nb, pb];
    // The short arc from pb back to pa — the one bulging toward the edge.
    let a0 = (pb - center).to_angle();
    let mut sweep = (pa - center).to_angle() - a0;
    if sweep > core::f64::consts::PI {
        sweep -= core::f64::consts::TAU;
    } else if sweep < -core::f64::consts::PI {
        sweep += core::f64::consts::TAU;
    }
    let chords = segments.max(3);
    points.extend(
        (1..chords).map(|i| center + r * DVec2::from_angle(a0 + sweep * i as f64 / chords as f64)),
    );
    points.push(pa);
    ensure_ccw(&mut points);

    Some(wrap_frame(
        &c.frame,
        place_on_rim(
            c,
            0.0,
            1.0,
            DynamicObject::Revolve {
                angle_degrees: 360.0,
                segments: rim_segments(c),
                child: Box::new(DynamicObject::Polygon {
                    points: points.iter().map(|p| p.to_array()).collect(),
                }),
            },
        ),
    ))
}

/// Reverse `points` in place if the loop winds clockwise (backends fill CCW).
fn ensure_ccw(points: &mut [DVec2]) {
    let area2: f64 = (0..points.len())
        .map(|i| {
            let p = points[i];
            let q = points[(i + 1) % points.len()];
            p.x * q.y - q.x * p.y
        })
        .sum();
    if area2 < 0.0 {
        points.reverse();
    }
}

/// `Translate ∘ Rotate` aligning local +Y with `up` (unit) at `origin`.
fn place_on_axis(origin: DVec3, up: DVec3, obj: DynamicObject) -> DynamicObject {
    let rot = DQuat::from_rotation_arc(DVec3::Y, up);
    let (z, y, x) = rot.to_euler(EulerRot::ZYX);
    DynamicObject::Translate {
        offset: origin.to_array(),
        child: Box::new(DynamicObject::Rotate {
            angles: [x.to_degrees(), y.to_degrees(), z.to_degrees()],
            child: Box::new(obj),
        }),
    }
}

/// Wrap `obj` in `Translate ∘ Rotate ∘ Scale ∘ Rotate` wrappers reproducing
/// the affine `frame` (identity passes through). The linear part splits as
/// `M = R₁·Σ·R₂` via the eigendecomposition of `MᵀM` (an SVD); a reflection
/// folds into a negated scale factor so both rotations stay proper.
fn wrap_frame(frame: &DAffine3, obj: DynamicObject) -> DynamicObject {
    if *frame == DAffine3::IDENTITY {
        return obj;
    }
    let m = frame.matrix3;
    let (q, lam) = jacobi_eigen(m.transpose() * m);
    let sigma = DVec3::new(
        lam.x.max(0.0).sqrt(),
        lam.y.max(0.0).sqrt(),
        lam.z.max(0.0).sqrt(),
    );
    if sigma.min_element() < 1.0e-12 * sigma.max_element().max(1.0e-300) {
        return obj; // degenerate (flattened) frame — nothing sane to emit
    }
    // N = M·Q·Σ⁻¹ is orthogonal; improper N moves its reflection into Σ.
    let n = m * q * DMat3::from_diagonal(sigma.recip());
    let (r1, sigma) = if n.determinant() < 0.0 {
        let flip = DVec3::new(1.0, 1.0, -1.0);
        (n * DMat3::from_diagonal(flip), sigma * flip)
    } else {
        (n, sigma)
    };
    let inner = DynamicObject::Rotate {
        angles: euler_zyx_degrees(q.transpose()),
        child: Box::new(obj),
    };
    let scaled = DynamicObject::Scale {
        factors: sigma.to_array(),
        child: Box::new(inner),
    };
    let rotated = DynamicObject::Rotate {
        angles: euler_zyx_degrees(r1),
        child: Box::new(scaled),
    };
    DynamicObject::Translate {
        offset: frame.translation.to_array(),
        child: Box::new(rotated),
    }
}

/// Eigendecomposition `m = Q·diag(λ)·Qᵀ` of a symmetric 3×3 matrix by cyclic
/// Jacobi rotations; `Q` is returned proper (det +1).
fn jacobi_eigen(m: DMat3) -> (DMat3, DVec3) {
    // a[r][c], symmetric.
    let mut a = [
        [m.x_axis.x, m.y_axis.x, m.z_axis.x],
        [m.x_axis.y, m.y_axis.y, m.z_axis.y],
        [m.x_axis.z, m.y_axis.z, m.z_axis.z],
    ];
    let mut q = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
    for _sweep in 0..12 {
        let off = a[0][1].abs() + a[0][2].abs() + a[1][2].abs();
        let diag = a[0][0].abs() + a[1][1].abs() + a[2][2].abs();
        if off <= 1.0e-15 * diag.max(1.0e-300) {
            break;
        }
        for (p, r) in [(0usize, 1usize), (0, 2), (1, 2)] {
            let apr = a[p][r];
            if apr.abs() <= 1.0e-30 {
                continue;
            }
            // Classic Jacobi rotation zeroing a[p][r].
            let theta = 0.5 * (a[r][r] - a[p][p]) / apr;
            let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt());
            let c = 1.0 / (t * t + 1.0).sqrt();
            let s = t * c;
            let rot_cols = |row: &mut [f64; 3]| {
                let (kp, kr) = (row[p], row[r]);
                row[p] = c * kp - s * kr;
                row[r] = s * kp + c * kr;
            };
            a.iter_mut().for_each(rot_cols);
            // Rows p and r (p < r in every pair above).
            let (head, tail) = a.split_at_mut(r);
            for (apk, ark) in head[p].iter_mut().zip(tail[0].iter_mut()) {
                let (x, y) = (*apk, *ark);
                *apk = c * x - s * y;
                *ark = s * x + c * y;
            }
            q.iter_mut().for_each(rot_cols);
        }
    }
    let mut qm = DMat3::from_cols(
        DVec3::new(q[0][0], q[1][0], q[2][0]),
        DVec3::new(q[0][1], q[1][1], q[2][1]),
        DVec3::new(q[0][2], q[1][2], q[2][2]),
    );
    if qm.determinant() < 0.0 {
        qm.z_axis = -qm.z_axis;
    }
    (qm, DVec3::new(a[0][0], a[1][1], a[2][2]))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::feature::{EdgeSel, FaceRef, FaceTag};

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

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

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

    fn face(path: &[u16], tag: FaceTag) -> FaceRef {
        FaceRef {
            path: LeafPath(path.to_vec()),
            face: tag,
        }
    }

    fn chamfer(size: f64, edges: Vec<EdgeSel>) -> EdgeFeatureData {
        EdgeFeatureData {
            kind: EdgeFeatureKind::Chamfer,
            size,
            segments: 16,
            edges,
        }
    }

    fn contains(obj: &DynamicObject, p: DVec3) -> bool {
        CsgLeaves::from_dynamic(obj, DAffine3::IDENTITY, None).contains(p)
    }

    #[test]
    fn straight_chamfer_cuts_the_corner_and_nothing_else() {
        let child = cube([2.0, 2.0, 2.0]);
        // Edge between +X and +Y faces (along Z at x=2, y=2), chamfer 0.4.
        let data = chamfer(
            0.4,
            vec![EdgeSel::new(
                face(&[], FaceTag::CubePosX),
                face(&[], FaceTag::CubePosY),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!(report.applied, 1);
        assert_eq!(report.skipped, 0);

        // The old edge line is gone along its whole length…
        for z in [0.05, 1.0, 1.95] {
            assert!(
                !contains(&lowered, DVec3::new(1.95, 1.95, z)),
                "corner material at z={z} must be cut"
            );
        }
        // …the cut plane sits at the setback (x + y = 4 − 0.4): just under
        // it is material, just over it is air.
        assert!(contains(&lowered, DVec3::new(1.75, 1.80, 1.0)));
        assert!(!contains(&lowered, DVec3::new(1.85, 1.85, 1.0)));
        // Untouched regions stay: the opposite edge and the faces away from
        // the setback band.
        assert!(contains(&lowered, DVec3::new(0.05, 0.05, 1.0)));
        assert!(contains(&lowered, DVec3::new(1.95, 0.5, 1.0)));
        assert!(contains(&lowered, DVec3::new(0.5, 1.95, 1.0)));
    }

    #[test]
    fn rim_chamfer_cuts_ring_without_touching_wall_below() {
        let child = cylinder(1.0, 2.0);
        let data = chamfer(
            0.2,
            vec![EdgeSel::new(
                face(&[], FaceTag::CylSide),
                face(&[], FaceTag::CylTop),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));

        // Rim corner removed all around.
        for theta in [0.0f64, 1.0, 2.5, 4.0, 5.5] {
            let (s, c) = theta.sin_cos();
            assert!(
                !contains(&lowered, DVec3::new(0.97 * c, 1.97, 0.97 * s)),
                "rim corner at θ={theta} must be cut"
            );
        }
        // The 45° cut surface: material just inside, air just outside.
        assert!(contains(&lowered, DVec3::new(0.85, 1.90, 0.0)));
        assert!(!contains(&lowered, DVec3::new(0.93, 1.95, 0.0)));
        // Wall below the setback and the cap center are untouched.
        assert!(contains(&lowered, DVec3::new(0.99, 1.75, 0.0)));
        assert!(contains(&lowered, DVec3::new(0.0, 1.99, 0.0)));
        assert!(contains(&lowered, DVec3::new(0.99, 0.5, 0.0)));
    }

    #[test]
    fn hole_rim_countersink_clips_to_the_funnel() {
        // Plate with a through-hole at (2, ·, 2), r = 0.5; countersink 0.15.
        let child = DynamicObject::Difference {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate([2.0, -0.5, 2.0], cylinder(0.5, 2.0)),
            ],
        };
        let data = chamfer(
            0.15,
            vec![EdgeSel::new(
                face(&[0], FaceTag::CubePosY),
                face(&[1], FaceTag::CylSide),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));

        // Funnel ring at the top: what was plate material at the rim is gone.
        assert!(!contains(&lowered, DVec3::new(2.55, 0.97, 2.0)));
        // Below the countersink depth the wall is intact…
        assert!(contains(&lowered, DVec3::new(2.55, 0.5, 2.0)));
        // …and the hole is not deepened (it is a through hole: the bore
        // itself stays air, the material right next to it stays material).
        assert!(!contains(&lowered, DVec3::new(2.0, 0.5, 2.0)));
        assert!(contains(&lowered, DVec3::new(2.7, 0.97, 2.0)));
    }

    fn fillet(size: f64, edges: Vec<EdgeSel>) -> EdgeFeatureData {
        EdgeFeatureData {
            kind: EdgeFeatureKind::Fillet,
            ..chamfer(size, edges)
        }
    }

    fn boss_on_plate() -> (DynamicObject, EdgeSel) {
        let child = DynamicObject::Union {
            children: vec![
                cube([4.0, 1.0, 4.0]),
                translate([2.0, 1.0, 2.0], cylinder(0.5, 1.0)),
            ],
        };
        let sel = EdgeSel::new(face(&[0], FaceTag::CubePosY), face(&[1], FaceTag::CylSide));
        (child, sel)
    }

    #[test]
    fn concave_chamfer_skips() {
        // Boss on a plate: the base ring is concave — chamfer must skip it.
        let (child, sel) = boss_on_plate();
        let data = chamfer(0.2, vec![sel]);
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (0, 1));
        assert_eq!(lowered, child, "nothing to lower — child passes through");
    }

    #[test]
    fn straight_fillet_keeps_the_arc_interior() {
        // Fillet 0.4 on the +X/+Y edge of a 2³ cube: arc center (1.6, 1.6),
        // r = 0.4 (90° dihedral: setback = radius).
        let child = cube([2.0, 2.0, 2.0]);
        let data = fillet(
            0.4,
            vec![EdgeSel::new(
                face(&[], FaceTag::CubePosX),
                face(&[], FaceTag::CubePosY),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));

        // The sharp corner is rounded off along the whole edge…
        for z in [0.05, 1.0, 1.95] {
            assert!(
                !contains(&lowered, DVec3::new(1.95, 1.95, z)),
                "corner material at z={z} must be cut"
            );
        }
        // …but unlike a chamfer, points beyond the chord and inside the
        // tangent circle survive: (1.85, 1.85) has x+y > 3.6 yet lies within
        // r of the arc center.
        assert!(contains(&lowered, DVec3::new(1.85, 1.85, 1.0)));
        // Faces outside the setback band are untouched.
        assert!(contains(&lowered, DVec3::new(1.95, 0.5, 1.0)));
        assert!(contains(&lowered, DVec3::new(0.5, 1.95, 1.0)));
    }

    #[test]
    fn convex_rim_fillet_is_a_subtractive_ring() {
        let data = fillet(
            0.2,
            vec![EdgeSel::new(
                face(&[], FaceTag::CylSide),
                face(&[], FaceTag::CylTop),
            )],
        );
        let (masks, report) = lower_masks(&data, &cylinder(1.0, 2.0));
        assert_eq!((report.applied, report.skipped), (1, 0));
        assert_eq!((masks.sub.len(), masks.add.len()), (1, 0));
        // The ring is a revolved meridian polygon around the rim axis.
        assert_polygon_ring(&masks.sub[0], 0.2, DVec2::new(0.8, -0.2));
    }

    #[test]
    fn boss_base_ring_fillet_is_an_additive_ring() {
        let (child, sel) = boss_on_plate();
        let data = fillet(0.2, vec![sel]);
        let (masks, report) = lower_masks(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));
        assert_eq!((masks.sub.len(), masks.add.len()), (0, 1));
        // Corner-fill arc center in boss meridian coords: radially outside
        // the boss wall, above the plate top, by (r, r) for the 90° corner.
        assert_polygon_ring(&masks.add[0], 0.2, DVec2::new(0.7, 0.2));

        // The lowered tree unions the ring on top of the child.
        let (lowered, _) = lower(&data, &child);
        match &lowered {
            DynamicObject::Union { children } => {
                assert_eq!(children.len(), 2);
                assert_eq!(children[0], child);
            }
            other => panic!("expected Union(child, ring), got {:?}", other.name()),
        }
    }

    /// Ring mask sanity: a revolved polygon whose arc chords sit at radius
    /// `r` from `center` (meridian coords) and whose points stay off-axis.
    fn assert_polygon_ring(mask: &DynamicObject, r: f64, center: DVec2) {
        let mut node = mask;
        loop {
            match node {
                DynamicObject::Translate { child, .. }
                | DynamicObject::Rotate { child, .. }
                | DynamicObject::Revolve { child, .. } => node = child,
                DynamicObject::Polygon { points } => {
                    assert!(points.iter().all(|p| p[0] > 0.0), "profile crosses axis");
                    let on_arc = points
                        .iter()
                        .filter(|p| ((DVec2::from_array(**p) - center).length() - r).abs() < 1.0e-9)
                        .count();
                    assert!(
                        on_arc >= 3,
                        "expected arc samples at radius {r} of {center}"
                    );
                    return;
                }
                other => panic!("unexpected ring mask node {:?}", other.name()),
            }
        }
    }

    #[test]
    fn mixed_chamfer_fillet_corner_leaves_no_horn() {
        // Chamfer on the PosX/PosY edge, fillet stacked on the PosX/PosZ
        // edge — they share the (2,2,2) corner. The fillet must resolve its
        // edge against the unfeatured cube: if the chamfer's cut trims it,
        // the fillet stops short and a sliver of the original sharp edge
        // survives at the corner (the "horn").
        let inner = DynamicObject::EdgeFeature {
            data: chamfer(
                0.3,
                vec![EdgeSel::new(
                    face(&[], FaceTag::CubePosX),
                    face(&[], FaceTag::CubePosY),
                )],
            ),
            child: Box::new(cube([2.0, 2.0, 2.0])),
        };
        let stacked = DynamicObject::EdgeFeature {
            data: EdgeFeatureData {
                kind: EdgeFeatureKind::Fillet,
                ..chamfer(
                    0.3,
                    vec![EdgeSel::new(
                        face(&[0], FaceTag::CubePosX),
                        face(&[0], FaceTag::CubePosZ),
                    )],
                )
            },
            child: Box::new(inner),
        };
        // Horn region: on the old sharp PosX/PosZ edge, inside the chamfer's
        // trim zone (y > 1.7) but below its cut plane (x + y < 3.7), outside
        // the fillet's tangent cylinder — must be removed by the full-length
        // fillet mask.
        assert!(!contains(&stacked, DVec3::new(1.98, 1.71, 1.98)));
        // The fillet's kept round survives away from the corner…
        assert!(contains(&stacked, DVec3::new(1.8, 1.0, 1.8)));
        // …the chamfer still cuts its own edge…
        assert!(!contains(&stacked, DVec3::new(1.9, 1.9, 1.0)));
        // …and plain faces stay.
        assert!(contains(&stacked, DVec3::new(1.98, 0.5, 0.5)));
    }

    #[test]
    fn rim_masks_adopt_leaf_tessellation_and_phase() {
        // 20-segment hole: mask cones/cylinders must tessellate at 20 with
        // the hole's angular origin (local +X), so their chords coincide
        // with the bore's in mesh booleans; the fillet ring revolve too.
        let hole = translate(
            [2.0, -0.5, 2.0],
            DynamicObject::Cylinder {
                radius: 0.5,
                height: 2.0,
                segments: 20,
            },
        );
        let child = DynamicObject::Difference {
            children: vec![cube([4.0, 1.0, 4.0]), hole],
        };
        let sel = || {
            vec![EdgeSel::new(
                face(&[0], FaceTag::CubePosY),
                face(&[1], FaceTag::CylSide),
            )]
        };

        let (masks, report) = lower_masks(&chamfer(0.15, sel()), &child);
        assert_eq!((report.applied, report.skipped), (1, 0));
        let csg = CsgLeaves::from_dynamic(&masks.sub[0], DAffine3::IDENTITY, None);
        assert!(!csg.leaves.is_empty());
        for leaf in &csg.leaves {
            let segments = match leaf.prim {
                LeafPrim::Cylinder { segments, .. } | LeafPrim::Cone { segments, .. } => segments,
                _ => continue,
            };
            assert_eq!(segments, 20, "mask must match the hole tessellation");
            // Phase: the mask's local +X maps onto the hole's angular
            // origin (world +X here — the hole leaf carries no rotation).
            let x_img = (leaf.world.matrix3 * DVec3::X).normalize();
            assert!(
                (x_img - DVec3::X).length() < 1e-9,
                "mask phase must anchor at the leaf origin, got {x_img:?}"
            );
        }

        // Fillet ring: the revolve sweep adopts the hole's 20 segments.
        let fillet_data = EdgeFeatureData {
            kind: EdgeFeatureKind::Fillet,
            ..chamfer(0.15, sel())
        };
        let (masks, report) = lower_masks(&fillet_data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));
        fn find_revolve_segments(obj: &DynamicObject) -> Option<usize> {
            match obj {
                DynamicObject::Revolve { segments, .. } => Some(*segments),
                DynamicObject::Translate { child, .. }
                | DynamicObject::Rotate { child, .. }
                | DynamicObject::Scale { child, .. } => find_revolve_segments(child),
                _ => None,
            }
        }
        assert_eq!(find_revolve_segments(&masks.sub[0]), Some(20));
    }

    #[test]
    fn wrap_frame_reproduces_affine() {
        let cases = [
            // Rotation + anisotropic scale + rotation + translation.
            DAffine3::from_translation(DVec3::new(1.0, 2.0, 3.0))
                * DAffine3::from_quat(DQuat::from_rotation_z(0.5))
                * DAffine3::from_scale(DVec3::new(2.0, 1.0, 0.5))
                * DAffine3::from_quat(DQuat::from_rotation_x(0.3)),
            // Reflection (negative determinant) folds into a scale factor.
            DAffine3::from_scale(DVec3::new(-2.0, 1.0, 0.5))
                * DAffine3::from_quat(DQuat::from_rotation_y(0.7)),
        ];
        for f in cases {
            let wrapped = wrap_frame(&f, cube([1.0, 1.0, 1.0]));
            let csg = CsgLeaves::from_dynamic(&wrapped, DAffine3::IDENTITY, None);
            let got = csg.leaves[0].world;
            assert!(
                (got.matrix3.x_axis - f.matrix3.x_axis).length() < 1e-9
                    && (got.matrix3.y_axis - f.matrix3.y_axis).length() < 1e-9
                    && (got.matrix3.z_axis - f.matrix3.z_axis).length() < 1e-9
                    && (got.translation - f.translation).length() < 1e-9,
                "wrapper chain must reproduce the frame: {got:?} vs {f:?}"
            );
        }
    }

    #[test]
    fn elliptical_rim_chamfer_cuts_in_world() {
        // Anisotropic [2,1,1] cylinder: the rim resolves in the leaf frame
        // and the mask is the frame image of the circular cut — the corner
        // ring is removed at both ellipse azimuths, walls and cap survive.
        let child = DynamicObject::Scale {
            factors: [2.0, 1.0, 1.0],
            child: Box::new(cylinder(1.0, 2.0)),
        };
        let data = chamfer(
            0.2,
            vec![EdgeSel::new(
                face(&[], FaceTag::CylSide),
                face(&[], FaceTag::CylTop),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));
        // Rim corner cut on the major (x) and minor (z) vertices.
        assert!(!contains(&lowered, DVec3::new(1.94, 1.97, 0.0)));
        assert!(!contains(&lowered, DVec3::new(0.0, 1.97, 0.97)));
        // Wall below the setback band and the cap center survive.
        assert!(contains(&lowered, DVec3::new(1.98, 1.5, 0.0)));
        assert!(contains(&lowered, DVec3::new(0.0, 1.5, 0.99)));
        assert!(contains(&lowered, DVec3::new(0.0, 1.99, 0.0)));
    }

    #[test]
    fn stale_paths_skip_and_valid_ones_apply() {
        let child = cube([2.0, 2.0, 2.0]);
        let data = chamfer(
            0.2,
            vec![
                EdgeSel::new(face(&[], FaceTag::CubePosX), face(&[], FaceTag::CubePosY)),
                // A path that doesn't exist in this child.
                EdgeSel::new(face(&[7], FaceTag::CubePosX), face(&[7], FaceTag::CubePosY)),
            ],
        );
        let (_, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 1));
    }

    #[test]
    fn nested_feature_child_paths_resolve() {
        // An outer feature whose child is itself an EdgeFeature addresses
        // the inner child region with a leading 0 — the path the engine
        // stores when stacking mixed kinds. Must resolve, not skip.
        let inner = DynamicObject::EdgeFeature {
            data: chamfer(
                0.2,
                vec![EdgeSel::new(
                    face(&[], FaceTag::CubePosX),
                    face(&[], FaceTag::CubePosY),
                )],
            ),
            child: Box::new(cube([2.0, 2.0, 2.0])),
        };
        let data = chamfer(
            0.2,
            vec![EdgeSel::new(
                face(&[0], FaceTag::CubeNegX),
                face(&[0], FaceTag::CubeNegY),
            )],
        );
        let (lowered, report) = lower(&data, &inner);
        assert_eq!((report.applied, report.skipped), (1, 0));
        // Both the inner and the outer corner are cut.
        assert!(!contains(&lowered, DVec3::new(1.95, 1.95, 1.0)));
        assert!(!contains(&lowered, DVec3::new(0.05, 0.05, 1.0)));
    }

    #[test]
    fn chamfer_follows_parameter_edits() {
        let edges = vec![EdgeSel::new(
            face(&[], FaceTag::CubePosX),
            face(&[], FaceTag::CubePosY),
        )];
        let data = chamfer(0.3, edges);
        // Same feature, two child sizes: the cut follows the moving edge.
        let (small, _) = lower(&data, &cube([2.0, 2.0, 2.0]));
        let (large, _) = lower(&data, &cube([3.0, 3.0, 3.0]));
        assert!(!contains(&small, DVec3::new(1.95, 1.95, 1.0)));
        assert!(contains(&large, DVec3::new(1.95, 1.95, 1.0)));
        assert!(!contains(&large, DVec3::new(2.95, 2.95, 1.0)));
    }

    #[test]
    fn size_clamps_to_adjacent_extent() {
        // Chamfer "1000" on a 2×2×2 cube: clamped to 0.45·2 = 0.9 per face,
        // so the body must survive.
        let child = cube([2.0, 2.0, 2.0]);
        let data = chamfer(
            1000.0,
            vec![EdgeSel::new(
                face(&[], FaceTag::CubePosX),
                face(&[], FaceTag::CubePosY),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!(report.applied, 1);
        assert!(contains(&lowered, DVec3::new(1.0, 0.2, 1.0)));
        assert!(
            contains(&lowered, DVec3::new(1.9, 1.05, 1.0)),
            "clamped cut leaves the face below the setback"
        );
        assert!(!contains(&lowered, DVec3::new(1.9, 1.9, 1.0)));
    }

    #[test]
    fn transformed_leaf_masks_follow_the_leaf_frame() {
        // The chamfered cube sits rotated + translated inside the child.
        let child = translate(
            [5.0, 0.0, 0.0],
            DynamicObject::Rotate {
                angles: [0.0, 0.0, 45.0],
                child: Box::new(cube([2.0, 2.0, 2.0])),
            },
        );
        let data = chamfer(
            0.4,
            vec![EdgeSel::new(
                face(&[], FaceTag::CubePosX),
                face(&[], FaceTag::CubePosY),
            )],
        );
        let (lowered, report) = lower(&data, &child);
        assert_eq!((report.applied, report.skipped), (1, 0));
        // Probe in the leaf's frame mapped to child coords: the +X/+Y edge
        // midpoint of the rotated cube.
        let world = DAffine3::from_translation(DVec3::new(5.0, 0.0, 0.0))
            * DAffine3::from_quat(DQuat::from_rotation_z(45.0_f64.to_radians()));
        let old_edge = world.transform_point3(DVec3::new(1.95, 1.95, 1.0));
        let deep = world.transform_point3(DVec3::new(1.0, 1.0, 1.0));
        assert!(!contains(&lowered, old_edge));
        assert!(contains(&lowered, deep));
    }
}