rscad_core/feature/

geom.rs

//! Analytic ray/membership evaluation of a [`DynamicObject`] tree.
//!
//! A tree is flattened into supported primitive leaves ([`GeomLeaf`], each
//! carrying its accumulated world transform) plus a boolean expression over
//! them ([`CsgExpr`]). Rays are intersected per leaf in local space and the
//! resulting `t`-intervals are combined through the booleans (classic ray-CSG
//! with normal flipping under `Difference`), so a hit identifies the *final
//! solid's* surface: leaf, face tag, and outward normal — with no meshing.
//!
//! ## Conventions
//!
//! - Primitive anchors follow the backend contract pinned in
//!   `rscad-csg-traits/src/dynamic.rs`: cube corner-anchored spanning
//!   `[0, size]`, cylinder/cone base on `y = 0` extending along `+Y`, sphere
//!   centered at the origin.
//! - `Rotate` accumulates as **Rz·Ry·Rx** (OpenSCAD order), matching every
//!   renderer (`DynamicObject::extract_top_transforms`, the SDF compiler,
//!   Manifold). Note `SceneGraph::world_transform` currently uses XYZ order —
//!   a pre-existing divergence; this module must match the pixels.
//! - A zero-length `Mirror` axis yields empty geometry (Manifold behavior).
//! - Bare 2D leaves in 3D context evaluate to empty (Manifold semantics).
//! - Rays keep their (unnormalized) world parameterization in local space, so
//!   `t` values are directly comparable across leaves.

use crate::feature::{FaceTag, LeafPath};
use crate::prelude_::*;

/// `t`-interval tie/degeneracy tolerance, in world ray-parameter units.
const T_EPS: f64 = 1.0e-9;
/// Segments used to tessellate circular face outlines (exact under any
/// affine transform, unlike transforming a circle's radius).
const OUTLINE_SEGMENTS: usize = 48;

/// Callback supplying a conservative **local-space** AABB for a subtree this
/// module cannot evaluate analytically (sketch/extrude/revolve/offset based).
/// Returning `None` drops the subtree from the model (it stops occluding).
pub type AabbHint<'a> = &'a dyn Fn(&DynamicObject) -> Option<(DVec3, DVec3)>;

/// A supported primitive in its local frame.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum LeafPrim {
    /// Corner-anchored box spanning `[0, size]`.
    Cube { size: DVec3 },
    /// Base at `y = 0`, axis `+Y`. `segments` is the render tessellation —
    /// edge-feature masks match it so boolean chords coincide.
    Cylinder {
        radius: f64,
        height: f64,
        segments: usize,
    },
    /// Base disk at `y = 0`, apex at `y = height`.
    Cone {
        radius: f64,
        height: f64,
        segments: usize,
    },
    /// Centered at the origin.
    Sphere { radius: f64 },
    /// Conservative stand-in for an unsupported subtree.
    Aabb { min: DVec3, max: DVec3 },
}

/// A primitive leaf with its accumulated world transform.
#[derive(Clone, Copy, Debug)]
pub struct GeomLeaf {
    pub prim: LeafPrim,
    /// Local → world (renderer rotation convention, see module docs).
    pub world: DAffine3,
    /// World → local (cached inverse of `world`).
    pub world_inv: DAffine3,
}

/// Boolean structure over leaf indices into [`CsgLeaves::leaves`].
#[derive(Clone, Debug, PartialEq)]
pub enum CsgExpr {
    Leaf(usize),
    Union(Vec<CsgExpr>),
    Inter(Vec<CsgExpr>),
    /// First operand minus the rest.
    Diff(Vec<CsgExpr>),
}

/// Flattened analytic model of one tree: primitive leaves + boolean structure.
#[derive(Clone, Debug)]
pub struct CsgLeaves {
    pub leaves: Vec<GeomLeaf>,
    pub expr: CsgExpr,
    /// Branch-path of each leaf (parallel to `leaves`), relative to the walk
    /// root passed to [`CsgLeaves::from_dynamic`].
    pub paths: Vec<LeafPath>,
}

/// One surface crossing along a ray.
#[derive(Clone, Copy, Debug)]
pub struct HitBoundary {
    /// Ray parameter in world units (ray direction's world length preserved).
    pub t: f64,
    /// Index into [`CsgLeaves::leaves`].
    pub leaf: usize,
    pub face: FaceTag,
    /// Surface belongs to a `Difference` subtrahend: its outward normal is
    /// the leaf's inward normal.
    pub flipped: bool,
}

/// A solid stretch of the ray: `enter.t < exit.t`.
#[derive(Clone, Copy, Debug)]
pub struct HitSpan {
    pub enter: HitBoundary,
    pub exit: HitBoundary,
}

impl CsgLeaves {
    /// Flatten `obj` into leaves + boolean structure. `world` is the
    /// transform above the walk root; `aabb_hint` supplies proxy bounds for
    /// unsupported subtrees (see [`AabbHint`]).
    pub fn from_dynamic(obj: &DynamicObject, world: DAffine3, aabb_hint: Option<AabbHint>) -> Self {
        Self::build(obj, world, aabb_hint, false)
    }

    /// Like [`Self::from_dynamic`], but edge-feature masks are omitted: the
    /// solid is the tree as if no chamfer/fillet were applied, with leaf
    /// paths keeping the same shape (a feature's child region stays `[0]`).
    /// Feature lowering resolves its edges against this view so a sibling
    /// feature's cut cannot trim a stored edge short — masks then span the
    /// full edge and overlap at shared corners instead of leaving a sliver
    /// of the original sharp edge standing (mixed chamfer/fillet corners).
    pub fn from_dynamic_unfeatured(
        obj: &DynamicObject,
        world: DAffine3,
        aabb_hint: Option<AabbHint>,
    ) -> Self {
        Self::build(obj, world, aabb_hint, true)
    }

    fn build(
        obj: &DynamicObject,
        world: DAffine3,
        aabb_hint: Option<AabbHint>,
        strip_features: bool,
    ) -> Self {
        let mut csg = CsgLeaves {
            leaves: Vec::new(),
            expr: CsgExpr::Union(Vec::new()),
            paths: Vec::new(),
        };
        let mut path = Vec::new();
        csg.expr = csg.walk(obj, world, &mut path, aabb_hint, strip_features);
        csg
    }

    /// The empty solid.
    pub fn is_empty(&self) -> bool {
        self.leaves.is_empty()
    }

    fn push_leaf(&mut self, prim: LeafPrim, world: DAffine3, path: &[u16]) -> CsgExpr {
        let idx = self.leaves.len();
        self.leaves.push(GeomLeaf {
            prim,
            world,
            world_inv: world.inverse(),
        });
        self.paths.push(LeafPath(path.to_vec()));
        CsgExpr::Leaf(idx)
    }

    fn walk_children(
        &mut self,
        children: &[DynamicObject],
        world: DAffine3,
        path: &mut Vec<u16>,
        hint: Option<AabbHint>,
        strip: bool,
    ) -> Vec<CsgExpr> {
        let branching = children.len() > 1;
        children
            .iter()
            .enumerate()
            .map(|(i, child)| {
                if branching {
                    path.push(i as u16);
                }
                let expr = self.walk(child, world, path, hint, strip);
                if branching {
                    path.pop();
                }
                expr
            })
            .collect()
    }

    fn walk(
        &mut self,
        obj: &DynamicObject,
        world: DAffine3,
        path: &mut Vec<u16>,
        hint: Option<AabbHint>,
        strip: bool,
    ) -> CsgExpr {
        let empty = || CsgExpr::Union(Vec::new());
        match obj {
            DynamicObject::Cube { size } => self.push_leaf(
                LeafPrim::Cube {
                    size: DVec3::from_array(*size),
                },
                world,
                path,
            ),
            DynamicObject::Sphere { radius, .. } => {
                self.push_leaf(LeafPrim::Sphere { radius: *radius }, world, path)
            }
            DynamicObject::Cylinder {
                radius,
                height,
                segments,
            } => self.push_leaf(
                LeafPrim::Cylinder {
                    radius: *radius,
                    height: *height,
                    segments: *segments,
                },
                world,
                path,
            ),
            DynamicObject::Cone {
                radius,
                height,
                segments,
            } => self.push_leaf(
                LeafPrim::Cone {
                    radius: *radius,
                    height: *height,
                    segments: *segments,
                },
                world,
                path,
            ),

            // Bare 2D geometry in 3D context is empty (Manifold semantics).
            DynamicObject::Circle { .. }
            | DynamicObject::Square { .. }
            | DynamicObject::Polygon { .. }
            | DynamicObject::Sketch { .. }
            | DynamicObject::Empty => empty(),

            // Unsupported solid subtrees become conservative AABB proxies so
            // they keep occluding picks; without a hint they vanish.
            DynamicObject::LinearExtrude { .. }
            | DynamicObject::Revolve { .. }
            | DynamicObject::Offset { .. } => match hint.and_then(|f| f(obj)) {
                Some((min, max)) => self.push_leaf(LeafPrim::Aabb { min, max }, world, path),
                None => empty(),
            },

            // Edge features evaluate through their lowered masks, with fixed
            // path indices (0 = child, 1 = subtractive masks, 2 = additive
            // masks) so graph-side path replay stays aligned whether or not
            // any mask resolved. Mask leaves all map to the feature node
            // itself. Masks walk without an AABB hint: Revolve-based fillet
            // rings must vanish here (like the SDF sentinel), not become
            // proxy boxes that would poison picking around the rim. In the
            // unfeatured view (`strip`) the masks are omitted entirely — the
            // child region keeps its `[0]` path so stored paths still match.
            DynamicObject::EdgeFeature { data, child } => {
                if strip {
                    path.push(0);
                    let child_expr = self.walk(child, world, path, hint, strip);
                    path.pop();
                    return child_expr;
                }
                let (masks, _) = crate::feature::lower::lower_masks(data, child);
                path.push(0);
                let child_expr = self.walk(child, world, path, hint, strip);
                path.pop();
                path.push(1);
                let sub_exprs = self.walk_children(&masks.sub, world, path, None, strip);
                path.pop();
                path.push(2);
                let add_exprs = self.walk_children(&masks.add, world, path, None, strip);
                path.pop();
                let cut = CsgExpr::Diff(vec![child_expr, CsgExpr::Union(sub_exprs)]);
                CsgExpr::Union(core::iter::once(cut).chain(add_exprs).collect())
            }

            DynamicObject::Union { children } | DynamicObject::Group { children } => {
                CsgExpr::Union(self.walk_children(children, world, path, hint, strip))
            }
            DynamicObject::Intersection { children } => {
                CsgExpr::Inter(self.walk_children(children, world, path, hint, strip))
            }
            DynamicObject::Difference { children } => {
                CsgExpr::Diff(self.walk_children(children, world, path, hint, strip))
            }

            DynamicObject::Translate { offset, child } => {
                let world = world * DAffine3::from_translation(DVec3::from_array(*offset));
                self.walk(child, world, path, hint, strip)
            }
            DynamicObject::Rotate { angles, child } => {
                // Rz·Ry·Rx — OpenSCAD order, matching every renderer.
                let q = DQuat::from_euler(
                    EulerRot::ZYX,
                    angles[2].to_radians(),
                    angles[1].to_radians(),
                    angles[0].to_radians(),
                );
                self.walk(child, world * DAffine3::from_quat(q), path, hint, strip)
            }
            DynamicObject::Scale { factors, child } => {
                let world = world * DAffine3::from_scale(DVec3::from_array(*factors));
                self.walk(child, world, path, hint, strip)
            }
            DynamicObject::Mirror { axis, child } => {
                let axis = DVec3::from_array(*axis);
                let len_sq = axis.length_squared();
                if len_sq < 1.0e-24 {
                    return empty(); // Manifold: zero mirror normal → empty
                }
                let n = axis / len_sq.sqrt();
                let m = DMat3::IDENTITY - 2.0 * DMat3::from_cols(n * n.x, n * n.y, n * n.z);
                self.walk(child, world * DAffine3::from_mat3(m), path, hint, strip)
            }
            DynamicObject::Color { child, .. } => self.walk(child, world, path, hint, strip),
            DynamicObject::Material { child, .. } => self.walk(child, world, path, hint, strip),
        }
    }

    /// All solid stretches of the ray through the final solid, sorted by `t`.
    /// `rd` need not be normalized — `t` stays in `ro + t·rd` units.
    pub fn spans(&self, ro: DVec3, rd: DVec3) -> Vec<HitSpan> {
        self.expr_spans(&self.expr, ro, rd)
    }

    /// First surface crossing with `t > t_min` (typically `0.0`).
    pub fn first_hit(&self, ro: DVec3, rd: DVec3, t_min: f64) -> Option<HitBoundary> {
        self.spans(ro, rd)
            .into_iter()
            .flat_map(|s| [s.enter, s.exit])
            .find(|b| b.t > t_min)
    }

    fn expr_spans(&self, expr: &CsgExpr, ro: DVec3, rd: DVec3) -> Vec<HitSpan> {
        match expr {
            CsgExpr::Leaf(idx) => leaf_spans(*idx, &self.leaves[*idx], ro, rd),
            CsgExpr::Union(children) => children
                .iter()
                .map(|c| self.expr_spans(c, ro, rd))
                .reduce(|a, b| merge_spans(a, b, BoolOp::Union))
                .unwrap_or_default(),
            CsgExpr::Inter(children) => {
                let mut iter = children.iter().map(|c| self.expr_spans(c, ro, rd));
                let Some(first) = iter.next() else {
                    return Vec::new();
                };
                iter.fold(first, |a, b| merge_spans(a, b, BoolOp::Inter))
            }
            CsgExpr::Diff(children) => {
                let mut iter = children.iter().map(|c| self.expr_spans(c, ro, rd));
                let Some(first) = iter.next() else {
                    return Vec::new();
                };
                iter.fold(first, |a, b| merge_spans(a, b, BoolOp::Diff))
            }
        }
    }

    /// Whether the final solid contains `p` (closed: boundary counts as inside).
    pub fn contains(&self, p: DVec3) -> bool {
        self.expr_contains(&self.expr, p)
    }

    fn expr_contains(&self, expr: &CsgExpr, p: DVec3) -> bool {
        match expr {
            CsgExpr::Leaf(idx) => leaf_contains(&self.leaves[*idx], p),
            CsgExpr::Union(children) => children.iter().any(|c| self.expr_contains(c, p)),
            CsgExpr::Inter(children) => {
                !children.is_empty() && children.iter().all(|c| self.expr_contains(c, p))
            }
            CsgExpr::Diff(children) => match children.split_first() {
                None => false,
                Some((first, rest)) => {
                    self.expr_contains(first, p) && !rest.iter().any(|c| self.expr_contains(c, p))
                }
            },
        }
    }

    /// Outward normal of the **final solid** at a boundary returned by
    /// [`Self::spans`]/[`Self::first_hit`] (unit length, world space).
    pub fn boundary_normal(&self, boundary: &HitBoundary, ro: DVec3, rd: DVec3) -> DVec3 {
        let leaf = &self.leaves[boundary.leaf];
        let p_local = leaf.world_inv.transform_point3(ro + boundary.t * rd);
        let n_local = local_face_normal(&leaf.prim, boundary.face, p_local);
        // Gradient chain rule: normals map by the inverse-transpose.
        let n_world = (leaf.world_inv.matrix3.transpose() * n_local).normalize_or_zero();
        if boundary.flipped { -n_world } else { n_world }
    }
}

// ── Per-leaf ray intersection ────────────────────────────────────

/// Intersect the ray with one leaf in its local frame. Returns 0 or 1 span
/// (every supported primitive is convex).
fn leaf_spans(idx: usize, leaf: &GeomLeaf, ro: DVec3, rd: DVec3) -> Vec<HitSpan> {
    let ro_l = leaf.world_inv.transform_point3(ro);
    let rd_l = leaf.world_inv.matrix3 * rd;
    let boundary = |t: f64, face: FaceTag| HitBoundary {
        t,
        leaf: idx,
        face,
        flipped: false,
    };
    let interval = match leaf.prim {
        LeafPrim::Cube { size } => slab_interval(ro_l, rd_l, DVec3::ZERO, size, false),
        LeafPrim::Aabb { min, max } => slab_interval(ro_l, rd_l, min, max, true),
        LeafPrim::Sphere { radius } => sphere_interval(ro_l, rd_l, radius),
        LeafPrim::Cylinder { radius, height, .. } => cylinder_interval(ro_l, rd_l, radius, height),
        LeafPrim::Cone { radius, height, .. } => cone_interval(ro_l, rd_l, radius, height),
    };
    match interval {
        Some(((t0, f0), (t1, f1))) if t1 - t0 > T_EPS => vec![HitSpan {
            enter: boundary(t0, f0),
            exit: boundary(t1, f1),
        }],
        _ => Vec::new(),
    }
}

type FaceInterval = Option<((f64, FaceTag), (f64, FaceTag))>;

/// Slab test against `[min, max]`, tagging each end with the binding axis
/// face (or `Unknown` for AABB proxies).
fn slab_interval(ro: DVec3, rd: DVec3, min: DVec3, max: DVec3, proxy: bool) -> FaceInterval {
    let mut t_in = (f64::NEG_INFINITY, FaceTag::Unknown);
    let mut t_out = (f64::INFINITY, FaceTag::Unknown);
    for axis in 0..3 {
        let (o, d) = (ro[axis], rd[axis]);
        if d == 0.0 {
            if o < min[axis] || o > max[axis] {
                return None;
            }
            continue;
        }
        let t_min_face = (min[axis] - o) / d;
        let t_max_face = (max[axis] - o) / d;
        // Moving along +d crosses the min (negative) face first; along −d the
        // max (positive) face first.
        let (near, far) = if d > 0.0 {
            (
                (t_min_face, tag_for(axis, false, proxy)),
                (t_max_face, tag_for(axis, true, proxy)),
            )
        } else {
            (
                (t_max_face, tag_for(axis, true, proxy)),
                (t_min_face, tag_for(axis, false, proxy)),
            )
        };
        if near.0 > t_in.0 {
            t_in = near;
        }
        if far.0 < t_out.0 {
            t_out = far;
        }
    }
    (t_in.0 < t_out.0).then_some((t_in, t_out))
}

fn tag_for(axis: usize, positive: bool, proxy: bool) -> FaceTag {
    if proxy {
        FaceTag::Unknown
    } else {
        FaceTag::cube(axis, positive)
    }
}

fn sphere_interval(ro: DVec3, rd: DVec3, radius: f64) -> FaceInterval {
    let a = rd.dot(rd);
    if a == 0.0 {
        return None;
    }
    let b = 2.0 * ro.dot(rd);
    let c = ro.dot(ro) - radius * radius;
    let disc = b * b - 4.0 * a * c;
    if disc <= 0.0 {
        return None;
    }
    let sq = disc.sqrt();
    let (t0, t1) = ((-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a));
    Some(((t0.min(t1), FaceTag::Sphere), (t0.max(t1), FaceTag::Sphere)))
}

fn cylinder_interval(ro: DVec3, rd: DVec3, radius: f64, height: f64) -> FaceInterval {
    // Radial constraint in the XZ plane.
    let a = rd.x * rd.x + rd.z * rd.z;
    let radial: FaceInterval = if a == 0.0 {
        (ro.x * ro.x + ro.z * ro.z <= radius * radius).then_some((
            (f64::NEG_INFINITY, FaceTag::CylSide),
            (f64::INFINITY, FaceTag::CylSide),
        ))
    } else {
        let b = 2.0 * (ro.x * rd.x + ro.z * rd.z);
        let c = ro.x * ro.x + ro.z * ro.z - radius * radius;
        let disc = b * b - 4.0 * a * c;
        (disc > 0.0).then(|| {
            let sq = disc.sqrt();
            let (t0, t1) = ((-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a));
            (
                (t0.min(t1), FaceTag::CylSide),
                (t0.max(t1), FaceTag::CylSide),
            )
        })
    };
    let radial = radial?;
    // Axial slab [0, height]: entering upward crosses the bottom cap.
    let axial: FaceInterval = if rd.y == 0.0 {
        (0.0..=height).contains(&ro.y).then_some((
            (f64::NEG_INFINITY, FaceTag::CylBottom),
            (f64::INFINITY, FaceTag::CylTop),
        ))
    } else {
        let ta = (0.0 - ro.y) / rd.y;
        let tb = (height - ro.y) / rd.y;
        Some(if rd.y > 0.0 {
            ((ta, FaceTag::CylBottom), (tb, FaceTag::CylTop))
        } else {
            ((tb, FaceTag::CylTop), (ta, FaceTag::CylBottom))
        })
    };
    let axial = axial?;
    let t_in = if radial.0.0 > axial.0.0 {
        radial.0
    } else {
        axial.0
    };
    let t_out = if radial.1.0 < axial.1.0 {
        radial.1
    } else {
        axial.1
    };
    (t_in.0 < t_out.0).then_some((t_in, t_out))
}

/// Cone solid `{ √(x²+z²) ≤ k·(h−y), 0 ≤ y ≤ h }`, `k = r/h`. Convex, so the
/// interval is found from boundary-crossing candidates + midpoint membership
/// (robust against the double-cone sheet and apex degeneracies).
fn cone_interval(ro: DVec3, rd: DVec3, radius: f64, height: f64) -> FaceInterval {
    if height <= 0.0 || radius <= 0.0 {
        return None;
    }
    let k = radius / height;
    let inside = |t: f64| -> bool {
        let p = ro + t * rd;
        if p.y < 0.0 || p.y > height {
            return false;
        }
        (p.x * p.x + p.z * p.z).sqrt() <= k * (height - p.y)
    };

    let mut candidates: Vec<(f64, FaceTag)> = Vec::with_capacity(4);
    // Sheet crossings: x² + z² = k²(h − y)².
    let a = rd.x * rd.x + rd.z * rd.z - k * k * rd.y * rd.y;
    let b = 2.0 * (ro.x * rd.x + ro.z * rd.z + k * k * rd.y * (height - ro.y));
    let c = ro.x * ro.x + ro.z * ro.z - k * k * (height - ro.y) * (height - ro.y);
    let mut sheet = |t: f64| {
        let y = ro.y + t * rd.y;
        if (-T_EPS..=height + T_EPS).contains(&y) {
            candidates.push((t, FaceTag::ConeSide));
        }
    };
    if a.abs() > 1.0e-14 {
        let disc = b * b - 4.0 * a * c;
        if disc > 0.0 {
            let sq = disc.sqrt();
            sheet((-b - sq) / (2.0 * a));
            sheet((-b + sq) / (2.0 * a));
        }
    } else if b.abs() > 1.0e-14 {
        sheet(-c / b);
    }
    // Base-plane crossing y = 0, inside the base disk.
    if rd.y != 0.0 {
        let t = -ro.y / rd.y;
        let p = ro + t * rd;
        if p.x * p.x + p.z * p.z <= radius * radius + T_EPS {
            candidates.push((t, FaceTag::ConeBottom));
        }
    }
    candidates.sort_by(|x, y| x.0.total_cmp(&y.0));
    candidates.dedup_by(|x, y| (x.0 - y.0).abs() < T_EPS);
    // Convex solid: the inside stretch is bounded by consecutive candidates.
    candidates
        .windows(2)
        .find(|w| inside(0.5 * (w[0].0 + w[1].0)))
        .map(|w| (w[0], w[1]))
}

// ── Per-leaf membership ──────────────────────────────────────────

fn leaf_contains(leaf: &GeomLeaf, p: DVec3) -> bool {
    let p = leaf.world_inv.transform_point3(p);
    match leaf.prim {
        LeafPrim::Cube { size } => {
            (0.0..=size.x).contains(&p.x)
                && (0.0..=size.y).contains(&p.y)
                && (0.0..=size.z).contains(&p.z)
        }
        LeafPrim::Aabb { min, max } => {
            (min.x..=max.x).contains(&p.x)
                && (min.y..=max.y).contains(&p.y)
                && (min.z..=max.z).contains(&p.z)
        }
        LeafPrim::Sphere { radius } => p.length_squared() <= radius * radius,
        LeafPrim::Cylinder { radius, height, .. } => {
            (0.0..=height).contains(&p.y) && p.x * p.x + p.z * p.z <= radius * radius
        }
        LeafPrim::Cone { radius, height, .. } => {
            if height <= 0.0 {
                return false;
            }
            let k = radius / height;
            (0.0..=height).contains(&p.y) && (p.x * p.x + p.z * p.z).sqrt() <= k * (height - p.y)
        }
    }
}

// ── Normals ──────────────────────────────────────────────────────

/// Leaf-local outward normal of `face` at local point `p` (on that face).
fn local_face_normal(prim: &LeafPrim, face: FaceTag, p: DVec3) -> DVec3 {
    match face {
        FaceTag::CubePosX => DVec3::X,
        FaceTag::CubeNegX => -DVec3::X,
        FaceTag::CubePosY | FaceTag::CylTop => DVec3::Y,
        FaceTag::CubeNegY | FaceTag::CylBottom | FaceTag::ConeBottom => -DVec3::Y,
        FaceTag::CubePosZ => DVec3::Z,
        FaceTag::CubeNegZ => -DVec3::Z,
        FaceTag::CylSide => DVec3::new(p.x, 0.0, p.z).normalize_or_zero(),
        FaceTag::ConeSide => {
            let LeafPrim::Cone { radius, height, .. } = prim else {
                return DVec3::Y;
            };
            let rho = (p.x * p.x + p.z * p.z).sqrt();
            if rho < 1.0e-12 {
                return DVec3::Y; // apex: any outward direction; pick the axis
            }
            // Gradient of √(x²+z²) − k(h−y):  (x/ρ, k, z/ρ).
            DVec3::new(p.x / rho, radius / height, p.z / rho).normalize()
        }
        FaceTag::Sphere => p.normalize_or_zero(),
        FaceTag::Unknown => aabb_proxy_normal(prim, p),
    }
}

/// Nearest-plane normal for AABB proxy faces (best effort — proxies are
/// never selectable, this only keeps hover math finite).
fn aabb_proxy_normal(prim: &LeafPrim, p: DVec3) -> DVec3 {
    let LeafPrim::Aabb { min, max } = prim else {
        return DVec3::Y;
    };
    let mut best = (f64::INFINITY, DVec3::Y);
    for axis in 0..3 {
        let mut n = DVec3::ZERO;
        n[axis] = -1.0;
        if (p[axis] - min[axis]).abs() < best.0 {
            best = ((p[axis] - min[axis]).abs(), n);
        }
        if (max[axis] - p[axis]).abs() < best.0 {
            best = ((max[axis] - p[axis]).abs(), -n);
        }
    }
    best.1
}

// ── Interval booleans ────────────────────────────────────────────

#[derive(Clone, Copy, PartialEq)]
enum BoolOp {
    Union,
    Inter,
    Diff,
}

/// Combine two disjoint sorted span lists. Difference events from `b` flip
/// their boundary orientation (subtrahend exit becomes result entry).
fn merge_spans(a: Vec<HitSpan>, b: Vec<HitSpan>, op: BoolOp) -> Vec<HitSpan> {
    // Fast paths.
    match op {
        BoolOp::Union if b.is_empty() => return a,
        BoolOp::Union if a.is_empty() => return b,
        BoolOp::Inter if a.is_empty() || b.is_empty() => return Vec::new(),
        BoolOp::Diff if a.is_empty() || b.is_empty() => return a,
        _ => {}
    }

    struct Event {
        t: f64,
        from_b: bool,
        is_enter: bool,
        boundary: HitBoundary,
    }
    let mut events: Vec<Event> = Vec::with_capacity(2 * (a.len() + b.len()));
    for (spans, from_b) in [(&a, false), (&b, true)] {
        for s in spans.iter() {
            events.push(Event {
                t: s.enter.t,
                from_b,
                is_enter: true,
                boundary: s.enter,
            });
            events.push(Event {
                t: s.exit.t,
                from_b,
                is_enter: false,
                boundary: s.exit,
            });
        }
    }
    events.sort_by(|x, y| x.t.total_cmp(&y.t));

    let combined = |ia: bool, ib: bool| match op {
        BoolOp::Union => ia || ib,
        BoolOp::Inter => ia && ib,
        BoolOp::Diff => ia && !ib,
    };
    let mut inside_a = false;
    let mut inside_b = false;
    let mut inside = false;
    let mut open: Option<HitBoundary> = None;
    let mut out: Vec<HitSpan> = Vec::new();
    for ev in events {
        if ev.from_b {
            inside_b = ev.is_enter;
        } else {
            inside_a = ev.is_enter;
        }
        let now = combined(inside_a, inside_b);
        if now == inside {
            continue;
        }
        inside = now;
        let mut boundary = ev.boundary;
        if op == BoolOp::Diff && ev.from_b {
            boundary.flipped = !boundary.flipped;
        }
        if now {
            open = Some(boundary);
        } else if let Some(enter) = open.take() {
            out.push(HitSpan {
                enter,
                exit: boundary,
            });
        }
    }
    // Normalize: drop degenerate slivers, join touching spans (their shared
    // boundary is interior to the combined solid, not a real surface).
    let mut norm: Vec<HitSpan> = Vec::with_capacity(out.len());
    for span in out {
        if span.exit.t - span.enter.t <= T_EPS {
            continue;
        }
        match norm.last_mut() {
            Some(prev) if span.enter.t - prev.exit.t <= T_EPS => prev.exit = span.exit,
            _ => norm.push(span),
        }
    }
    norm
}

// ── Face outlines (for viewport highlights) ──────────────────────

impl GeomLeaf {
    /// Closed world-space polylines outlining `face`. Circles are tessellated
    /// locally and transformed, so non-uniform scale renders the true ellipse.
    /// Spheres and proxies return no outline.
    pub fn face_outline(&self, face: FaceTag) -> Vec<Vec<DVec3>> {
        let circle = |y: f64, r: f64| -> Vec<DVec3> {
            (0..OUTLINE_SEGMENTS)
                .map(|i| {
                    let a = core::f64::consts::TAU * i as f64 / OUTLINE_SEGMENTS as f64;
                    self.world
                        .transform_point3(DVec3::new(r * a.cos(), y, r * a.sin()))
                })
                .collect()
        };
        match (&self.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 Vec::new(),
                };
                let (u, v) = ((axis + 1) % 3, (axis + 2) % 3);
                let corners = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)];
                let quad = corners
                    .iter()
                    .map(|&(cu, cv)| {
                        let mut p = DVec3::ZERO;
                        p[axis] = if positive { size[axis] } else { 0.0 };
                        p[u] = cu * size[u];
                        p[v] = cv * size[v];
                        self.world.transform_point3(p)
                    })
                    .collect();
                vec![quad]
            }
            (LeafPrim::Cylinder { radius, height, .. }, FaceTag::CylTop) => {
                vec![circle(*height, *radius)]
            }
            (LeafPrim::Cylinder { radius, .. }, FaceTag::CylBottom) => vec![circle(0.0, *radius)],
            (LeafPrim::Cylinder { radius, height, .. }, FaceTag::CylSide) => {
                vec![circle(0.0, *radius), circle(*height, *radius)]
            }
            (LeafPrim::Cone { radius, .. }, FaceTag::ConeBottom | FaceTag::ConeSide) => {
                vec![circle(0.0, *radius)]
            }
            _ => Vec::new(),
        }
    }
}

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

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

    fn unit_ray_x() -> (DVec3, DVec3) {
        (DVec3::new(-10.0, 0.5, 0.5), DVec3::X)
    }

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

    #[test]
    fn cube_axis_ray_hits_neg_then_pos_x() {
        let csg = leaves(&cube([1.0, 1.0, 1.0]));
        let spans = csg.spans(unit_ray_x().0, unit_ray_x().1);
        assert_eq!(spans.len(), 1);
        let s = spans[0];
        assert!((s.enter.t - 10.0).abs() < 1e-12);
        assert!((s.exit.t - 11.0).abs() < 1e-12);
        assert_eq!(s.enter.face, FaceTag::CubeNegX);
        assert_eq!(s.exit.face, FaceTag::CubePosX);
        let n = csg.boundary_normal(&s.enter, unit_ray_x().0, unit_ray_x().1);
        assert!((n - (-DVec3::X)).length() < 1e-12);
    }

    #[test]
    fn cube_miss_and_axis_parallel_guard() {
        let csg = leaves(&cube([1.0, 1.0, 1.0]));
        // Parallel to X at y outside the slab.
        assert!(csg.spans(DVec3::new(-10.0, 2.0, 0.5), DVec3::X).is_empty());
        // Ray with a zero component inside the slab still hits.
        assert_eq!(csg.spans(DVec3::new(0.5, 0.5, -5.0), DVec3::Z).len(), 1);
    }

    #[test]
    fn translated_rotated_cube_matches_renderer_convention() {
        // rotate([0,0,90]) then translate([2,0,0]): +X face maps to +Y side.
        let obj = DynamicObject::Translate {
            offset: [2.0, 0.0, 0.0],
            child: Box::new(DynamicObject::Rotate {
                angles: [0.0, 0.0, 90.0],
                child: Box::new(cube([1.0, 2.0, 3.0])),
            }),
        };
        let csg = leaves(&obj);
        // Rz(90°): local +X → world +Y. Cube occupies x∈[0,-2]? No:
        // corner-anchored [0,1]×[0,2]×[0,3] rotated by Rz90 → x∈[-2,0], y∈[0,1],
        // then translated → x∈[0,2] shifted: x∈[2-2,2]=[0,2], y∈[0,1], z∈[0,3].
        let ro = DVec3::new(1.0, 10.0, 1.5);
        let spans = csg.spans(ro, -DVec3::Y);
        assert_eq!(spans.len(), 1);
        // Entering downward crosses the face whose world normal is +Y: the
        // rotated local +X face.
        assert_eq!(spans[0].enter.face, FaceTag::CubePosX);
        let n = csg.boundary_normal(&spans[0].enter, ro, -DVec3::Y);
        assert!((n - DVec3::Y).length() < 1e-9);
    }

    #[test]
    fn non_uniform_scale_keeps_world_t() {
        let obj = DynamicObject::Scale {
            factors: [2.0, 3.0, 4.0],
            child: Box::new(cube([1.0, 1.0, 1.0])),
        };
        let csg = leaves(&obj);
        let (ro, rd) = (DVec3::new(-10.0, 1.0, 1.0), DVec3::X);
        let spans = csg.spans(ro, rd);
        assert_eq!(spans.len(), 1);
        assert!((spans[0].enter.t - 10.0).abs() < 1e-12);
        assert!(
            (spans[0].exit.t - 12.0).abs() < 1e-12,
            "scaled span is 2 long"
        );
    }

    #[test]
    fn mirror_flips_geometry_not_normals_outwardness() {
        let obj = DynamicObject::Mirror {
            axis: [1.0, 0.0, 0.0],
            child: Box::new(cube([1.0, 1.0, 1.0])),
        };
        let csg = leaves(&obj);
        // Mirrored cube occupies x ∈ [-1, 0].
        let (ro, rd) = (DVec3::new(-10.0, 0.5, 0.5), DVec3::X);
        let spans = csg.spans(ro, rd);
        assert_eq!(spans.len(), 1);
        assert!((spans[0].enter.t - 9.0).abs() < 1e-12);
        let n = csg.boundary_normal(&spans[0].enter, ro, rd);
        assert!(
            (n - (-DVec3::X)).length() < 1e-12,
            "outward normal faces the ray"
        );
    }

    #[test]
    fn zero_mirror_axis_is_empty() {
        let obj = DynamicObject::Mirror {
            axis: [0.0, 0.0, 0.0],
            child: Box::new(cube([1.0, 1.0, 1.0])),
        };
        assert!(leaves(&obj).is_empty());
    }

    #[test]
    fn sphere_and_cylinder_intervals() {
        let csg = leaves(&DynamicObject::Sphere {
            radius: 2.0,
            segments: 8,
            stacks: 8,
        });
        let spans = csg.spans(DVec3::new(-10.0, 0.0, 0.0), DVec3::X);
        assert_eq!(spans.len(), 1);
        assert!((spans[0].enter.t - 8.0).abs() < 1e-12);
        assert!((spans[0].exit.t - 12.0).abs() < 1e-12);

        let cyl = leaves(&DynamicObject::Cylinder {
            radius: 1.0,
            height: 4.0,
            segments: 8,
        });
        // Down the axis: enters the top cap, exits the bottom.
        let spans = cyl.spans(DVec3::new(0.0, 10.0, 0.0), -DVec3::Y);
        assert_eq!(spans.len(), 1);
        assert_eq!(spans[0].enter.face, FaceTag::CylTop);
        assert_eq!(spans[0].exit.face, FaceTag::CylBottom);
        assert!((spans[0].enter.t - 6.0).abs() < 1e-12);
        // Sideways through the wall at mid-height.
        let spans = cyl.spans(DVec3::new(-5.0, 2.0, 0.0), DVec3::X);
        assert_eq!(spans[0].enter.face, FaceTag::CylSide);
        assert!((spans[0].enter.t - 4.0).abs() < 1e-12);
        let n = cyl.boundary_normal(&spans[0].enter, DVec3::new(-5.0, 2.0, 0.0), DVec3::X);
        assert!((n - (-DVec3::X)).length() < 1e-12);
    }

    #[test]
    fn cone_side_and_base() {
        let cone = leaves(&DynamicObject::Cone {
            radius: 2.0,
            height: 4.0,
            segments: 8,
        });
        // Horizontal ray at half height: local radius there is 1.
        let (ro, rd) = (DVec3::new(-5.0, 2.0, 0.0), DVec3::X);
        let spans = cone.spans(ro, rd);
        assert_eq!(spans.len(), 1);
        assert_eq!(spans[0].enter.face, FaceTag::ConeSide);
        assert!((spans[0].enter.t - 4.0).abs() < 1e-9);
        assert!((spans[0].exit.t - 6.0).abs() < 1e-9);
        let n = cone.boundary_normal(&spans[0].enter, ro, rd);
        // Outward, tilted upward by the slope: n ∝ (−1, k, 0), k = r/h = 0.5.
        let expect = DVec3::new(-1.0, 0.5, 0.0).normalize();
        assert!((n - expect).length() < 1e-9);
        // Vertical ray through the base.
        let spans = cone.spans(DVec3::new(0.5, -5.0, 0.0), DVec3::Y);
        assert_eq!(spans[0].enter.face, FaceTag::ConeBottom);
        assert_eq!(spans[0].exit.face, FaceTag::ConeSide);
    }

    #[test]
    fn union_merges_touching_and_overlapping() {
        let obj = DynamicObject::Union {
            children: vec![
                cube([1.0, 1.0, 1.0]),
                DynamicObject::Translate {
                    offset: [1.0, 0.0, 0.0],
                    child: Box::new(cube([1.0, 1.0, 1.0])),
                },
            ],
        };
        let csg = leaves(&obj);
        let spans = csg.spans(unit_ray_x().0, unit_ray_x().1);
        assert_eq!(spans.len(), 1, "touching cubes merge into one span");
        assert!((spans[0].exit.t - 12.0).abs() < 1e-12);
    }

    #[test]
    fn difference_flips_subtrahend_normal() {
        // 2×1×1 plate minus a unit cube in its middle-left: ray along X sees
        // plate wall, then the hole wall (subtracted cube's face, flipped).
        let obj = DynamicObject::Difference {
            children: vec![
                cube([3.0, 1.0, 1.0]),
                DynamicObject::Translate {
                    offset: [1.0, -0.5, -0.5],
                    child: Box::new(cube([1.0, 2.0, 2.0])),
                },
            ],
        };
        let csg = leaves(&obj);
        let (ro, rd) = unit_ray_x();
        let spans = csg.spans(ro, rd);
        assert_eq!(spans.len(), 2, "plate is split by the through-slot");
        // First span: [x=0, x=1], exit face is the slot wall = subtracted
        // cube's NegX face, flipped.
        let exit = spans[0].exit;
        assert_eq!(exit.face, FaceTag::CubeNegX);
        assert!(exit.flipped);
        let n = csg.boundary_normal(&exit, ro, rd);
        assert!(
            (n - DVec3::X).length() < 1e-12,
            "flipped normal points into the slot (out of the material)"
        );
        // contains: inside plate material, not inside the slot.
        assert!(csg.contains(DVec3::new(0.5, 0.5, 0.5)));
        assert!(!csg.contains(DVec3::new(1.5, 0.5, 0.5)));
    }

    #[test]
    fn intersection_binds_to_tighter_solid() {
        let obj = DynamicObject::Intersection {
            children: vec![
                cube([2.0, 2.0, 2.0]),
                DynamicObject::Translate {
                    offset: [1.0, 0.0, 0.0],
                    child: Box::new(cube([2.0, 2.0, 2.0])),
                },
            ],
        };
        let csg = leaves(&obj);
        let (ro, rd) = (DVec3::new(-10.0, 1.0, 1.0), DVec3::X);
        let spans = csg.spans(ro, rd);
        assert_eq!(spans.len(), 1);
        assert!(
            (spans[0].enter.t - 11.0).abs() < 1e-12,
            "second cube binds entry"
        );
        assert!(
            (spans[0].exit.t - 12.0).abs() < 1e-12,
            "first cube binds exit"
        );
    }

    #[test]
    fn empty_leaves_and_2d_are_invisible() {
        assert!(leaves(&DynamicObject::Empty).is_empty());
        assert!(leaves(&DynamicObject::Square { size: [1.0, 1.0] }).is_empty());
        let diff = DynamicObject::Difference { children: vec![] };
        assert!(
            leaves(&diff)
                .spans(unit_ray_x().0, unit_ray_x().1)
                .is_empty()
        );
    }

    #[test]
    fn unsupported_subtree_uses_aabb_hint() {
        let extrude = DynamicObject::LinearExtrude {
            height: 2.0,
            twist: 0.0,
            scale: 1.0,
            slices: 1,
            child: Box::new(DynamicObject::Square { size: [1.0, 1.0] }),
        };
        // Without a hint the subtree vanishes.
        assert!(leaves(&extrude).is_empty());
        // With a hint it occludes as a box, tagged Unknown.
        let hint = |obj: &DynamicObject| match obj {
            DynamicObject::LinearExtrude { height, .. } => {
                Some((DVec3::ZERO, DVec3::new(1.0, *height, 1.0)))
            }
            _ => None,
        };
        let csg = CsgLeaves::from_dynamic(&extrude, DAffine3::IDENTITY, Some(&hint));
        let spans = csg.spans(DVec3::new(-5.0, 0.5, 0.5), DVec3::X);
        assert_eq!(spans.len(), 1);
        assert_eq!(spans[0].enter.face, FaceTag::Unknown);
        assert!(!spans[0].enter.face.is_selectable());
    }

    #[test]
    fn branch_paths_skip_single_child_wrappers() {
        // Union[cube, Translate(Rotate(cube))] — second leaf's path is just [1].
        let obj = DynamicObject::Union {
            children: vec![
                cube([1.0; 3]),
                DynamicObject::Translate {
                    offset: [5.0, 0.0, 0.0],
                    child: Box::new(DynamicObject::Rotate {
                        angles: [0.0, 45.0, 0.0],
                        child: Box::new(cube([1.0; 3])),
                    }),
                },
            ],
        };
        let csg = leaves(&obj);
        assert_eq!(csg.paths[0], LeafPath(vec![0]));
        assert_eq!(csg.paths[1], LeafPath(vec![1]));
        // A single-child Difference consumes no index.
        let single = DynamicObject::Difference {
            children: vec![cube([1.0; 3])],
        };
        let csg = leaves(&single);
        assert_eq!(csg.paths[0], LeafPath(Vec::new()));
    }

    #[test]
    fn first_hit_skips_behind_camera() {
        let csg = leaves(&cube([1.0, 1.0, 1.0]));
        // Origin inside the cube: first hit in front is the exit face.
        let hit = csg.first_hit(DVec3::splat(0.5), DVec3::X, 0.0).unwrap();
        assert_eq!(hit.face, FaceTag::CubePosX);
    }

    #[test]
    fn face_outline_quad_and_circles() {
        let csg = leaves(&cube([1.0, 2.0, 3.0]));
        let outline = csg.leaves[0].face_outline(FaceTag::CubePosX);
        assert_eq!(outline.len(), 1);
        assert_eq!(outline[0].len(), 4);
        assert!(outline[0].iter().all(|p| (p.x - 1.0).abs() < 1e-12));

        let cyl = leaves(&DynamicObject::Cylinder {
            radius: 1.0,
            height: 2.0,
            segments: 8,
        });
        assert_eq!(cyl.leaves[0].face_outline(FaceTag::CylSide).len(), 2);
        assert_eq!(cyl.leaves[0].face_outline(FaceTag::CylTop).len(), 1);

        let sph = leaves(&DynamicObject::Sphere {
            radius: 1.0,
            segments: 8,
            stacks: 8,
        });
        assert!(sph.leaves[0].face_outline(FaceTag::Sphere).is_empty());
    }
}