axiolid-mesh-boolean-boolmesh 0.3.2

//! The `MeshBoolean` implementation.
//!
//! # Whose fault an error is
//!
//! Input faults are the caller's: an operand that fails the orientation
//! gate (`InvalidInput`, `Degenerate`, naming which argument) or that the
//! kernel cannot build as a manifold (`NotManifold`). A result that fails
//! [`check_result`] is `BackendContractViolation`: the inputs were already
//! validated, so a bad result is this provider's defect, never blamed on
//! the caller. A failure inside the kernel's own solve (for example an odd
//! edge-point count in `pair_up`, #101) is `BackendContractViolation` for
//! the same reason: the operands passed every input gate, so the solve
//! failing on them is this provider's defect.

use crate::csg::{compute_boolean, OpType};
use axiolid_contracts::{
    Backend, BackendDescriptor, BackendId, CancellationGranularity, Determinism, ExecutionOptions,
    ExecutionTarget, GeomError, GeomResult, ScratchRequirement,
};
use axiolid_core::BooleanOperator;
use axiolid_mesh::{AttributeChannel, AttributeFate, Blend, DropReason, TriMesh};
use axiolid_mesh_boolean_contract::{
    merge_fates, symmetric_difference_via_composition, BooleanEvidence, BooleanOutcome, MeshBoolean,
};

use crate::attributes::{carry, sources_by_plane, FaceSource};
use crate::convert::{from_boolean_mesh, six_signed_volume, to_manifold};

/// Mesh boolean built on the algorithm absorbed from `boolmesh` (pure Rust,
/// `glam`-only, MPL-2.0).
///
/// Adopted in `docs/adr/0014` and absorbed into this crate's private `csg`
/// module by `docs/adr/0047`. This type owns conversion and contract
/// enforcement; no `csg` type is part of the public API, so the kernel can
/// be replaced without a consumer noticing.
#[derive(Debug, Clone, Copy, Default)]
pub struct BoolmeshBoolean;

impl BoolmeshBoolean {
    /// Stable identifier used in errors and explicit provider selection.
    pub const ID: BackendId = BackendId::new("boolmesh");

    /// Construct the provider.
    pub const fn new() -> Self {
        Self
    }

    /// Subtract axis-aligned box cutters from an axis-aligned box subject using
    /// an exact closed-form construction, bypassing the general solver.
    ///
    /// Returns `Ok(None)` when the operands are outside this path's competence,
    /// so the caller falls back to [`MeshBoolean::subtract_many`] rather than
    /// receiving an approximate answer. Declining cases:
    ///
    /// * subject or any tool is not an axis-aligned box (recognised
    ///   structurally -- triangle count, corner lattice, and face planes -- not
    ///   by bounding box, since every mesh has one of those)
    /// * no tools, or no tool overlapping the subject
    /// * the induced grid would exceed `max_cells`
    ///
    /// # Why opt-in rather than automatic
    ///
    /// Dispatching on shape automatically would make output topology depend on
    /// input in a way the caller cannot predict: the same wall would produce
    /// different triangle counts depending on whether its openings happened to
    /// be axis-aligned. Callers that want the speed ask for it and handle the
    /// `None`; callers that want one predictable code path never see this.
    ///
    /// # Guarantees
    ///
    /// The result is watertight by construction (adjacent cells address shared
    /// grid vertices by integer index) and deterministic across processes
    /// (vertex identity is an ordered map, not a randomly-seeded hash map).
    ///
    /// The returned evidence has [`BooleanEvidence::analytic_path`] set, so a
    /// caller can tell this result apart from a general-solver result.
    pub fn subtract_boxes_analytic(
        &self,
        subject: &TriMesh,
        tools: &[TriMesh],
        options: &ExecutionOptions,
        max_cells: usize,
    ) -> GeomResult<Option<BooleanOutcome>> {
        options.check_cancelled()?;

        // Tolerance is scale-relative: an absolute epsilon would reject a
        // building-scale box and accept a millimetre-scale non-box.
        let d = subject.bounds().diagonal();
        let eps = d.x.max(d.y).max(d.z) * 1e-9;

        let Some(host) = crate::box_detect::recognise(subject, eps) else {
            return Ok(None);
        };
        let mut cutters = Vec::with_capacity(tools.len());
        for tool in tools {
            let Some(b) = crate::box_detect::recognise(tool, eps) else {
                return Ok(None);
            };
            cutters.push(b);
        }

        let Some(mut result) = crate::cellular::subtract_boxes(&host, &cutters, max_cells) else {
            return Ok(None);
        };

        // The analytic path is exact, but it is still a provider result: hold it
        // to the same orientation contract as the general path rather than
        // trusting the construction.
        check_result(&result, BooleanOperator::Difference)?;

        let tool_refs: Vec<&TriMesh> = tools.iter().collect();
        let fates = carry_analytic(subject, tools, &mut result);
        let evidence = evidence_for(subject, &tool_refs, &result, 1)
            .with_analytic_path(true)
            .with_attribute_fates(fates);
        Ok(Some(BooleanOutcome::new(result, evidence)))
    }
}

impl Backend for BoolmeshBoolean {
    fn descriptor(&self) -> BackendDescriptor {
        BackendDescriptor::new(Self::ID, ExecutionTarget::PortableCpu)
    }
}

/// Verify the result against the contract before handing it back.
///
/// `boolmesh` reports its own manifoldness; we additionally check orientation,
/// because a returned inside-out solid would poison every subsequent operation
/// in a subtract chain. Blamed on the backend, not the caller: the inputs were
/// already validated on the way in.
fn check_result(result: &TriMesh, operation: BooleanOperator) -> GeomResult<()> {
    // An empty result is legitimate: subtracting a tool that fully contains the
    // subject leaves nothing. Orientation is undefined for it, so stop here.
    //
    // Removing this early return is behaviour-preserving (an empty mesh sums to
    // zero, not a negative), so no test can distinguish it. It is kept because
    // it states the intent at the point of the decision.
    if result.indices.is_empty() {
        return Ok(());
    }
    let six_volume = six_signed_volume(&result.positions, &result.indices);
    if !six_volume.is_finite() {
        return Err(GeomError::BackendContractViolation {
            backend: BoolmeshBoolean::ID,
            detail: format!("{operation:?} returned non-finite signed volume"),
        });
    }
    if six_volume < 0.0 {
        return Err(GeomError::BackendContractViolation {
            backend: BoolmeshBoolean::ID,
            detail: format!(
                "{operation:?} returned an inside-out solid (signed volume {:.6})",
                six_volume / 6.0
            ),
        });
    }
    Ok(())
}

impl MeshBoolean for BoolmeshBoolean {
    /// Honest about the general path, which is not byte-reproducible.
    ///
    /// Upstream `boolmesh` dedups vertices through a randomly-seeded
    /// `HashMap`, so its output ordering varies between processes even
    /// single-threaded. That is the same root cause documented in
    /// [`Self::subtract_boxes_analytic`], which avoids it with a `BTreeMap`. The general
    /// path therefore cannot claim [`Determinism::Bitwise`] at any thread
    /// count, and enabling `parallel` does not make it weaker than it
    /// already is.
    ///
    /// [`Determinism::Topological`] is the honest ceiling: the result's
    /// connectivity is stable, its vertex ordering is not. Callers needing
    /// byte reproducibility use [`Self::subtract_boxes_analytic`],
    /// whose ordered vertex identity makes it reproducible across processes.
    fn determinism(&self) -> Determinism {
        Determinism::Topological
    }

    /// Measured, not guessed.
    ///
    /// `boolmesh` builds a Morton collider and intersection tables sized by the
    /// combined input. It exposes no bound, so this was measured directly with
    /// a counting global allocator (`src/bin/scratch_probe/`) across all four
    /// operations at 24 to 1,536 input triangles, after one discarded warmup
    /// call so the first operation measured is not charged for process
    /// startup (#110):
    ///
    /// ```text
    ///  triangles   peak bytes   bytes/triangle   worst operation
    ///         24        52424             2184   SymmetricDifference
    ///         96        98920             1030   SymmetricDifference
    ///        384       343768              895   SymmetricDifference
    ///       1536      1349668              878   SymmetricDifference
    /// ```
    ///
    /// Consumption is linear in input size, converging to under 1 KiB per
    /// triangle; the higher ratio at small inputs is fixed setup cost being
    /// divided by few triangles. `SymmetricDifference` is worst because it is
    /// composed from three passes and holds intermediates alive.
    ///
    /// The declared bound is 4 KiB per triangle: about 1.9x the worst
    /// observed ratio, as headroom for allocator variance and future operand
    /// shapes. `tests/scratch_bound.rs` measures the same workloads and fails
    /// if any peak exceeds it. A declared bound that is occasionally too low
    /// is worse than `Unbounded`, because it makes a budget look enforced
    /// when it is not.
    fn scratch_requirement(&self) -> ScratchRequirement {
        ScratchRequirement::PerElement {
            bytes_per_element: 4096,
        }
    }

    /// `boolmesh` takes no cancellation handle, so nothing can interrupt a
    /// single boolean once it starts. Declared honestly: the batch override
    /// polls between groups, which is the only real poll point available.
    fn cancellation_granularity(&self) -> CancellationGranularity {
        CancellationGranularity::BetweenOperations
    }

    /// Group mutually disjoint cutters and remove each group with one
    /// boolean, instead of one boolean per cutter.
    ///
    /// Rests on `(S \ A) \ B == S \ (A union B)` and on a concatenation of
    /// disjoint solids being their union. Bounding-box grouping over-separates
    /// but never wrongly fuses, so the result is identical to the sequential
    /// default -- gated by volume equality in `tests/batch.rs`.
    fn subtract_many(
        &self,
        subject: &TriMesh,
        tools: &[TriMesh],
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        self.subtract_grouped(subject, tools, options)
    }

    fn union_many(
        &self,
        solids: &[TriMesh],
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        self.union_tree(solids, options)
    }

    fn boolean(
        &self,
        subject: &TriMesh,
        tool: &TriMesh,
        operation: BooleanOperator,
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        self.boolean_impl(subject, tool, operation, options, false)
    }
}

impl BoolmeshBoolean {
    /// Alternative to [`MeshBoolean::boolean`] using a cheaper winding-number
    /// classification for the same result.
    ///
    /// See `csg::boolean03::kernel03::winding03_fast`'s doc comment for the
    /// guarantee: not an approximation, but a bug in edge-break detection
    /// would mislabel a whole connected component instead of one vertex.
    /// Opt-in rather than the default until this has run against a wider
    /// correctness corpus than the differential tests already covering it.
    ///
    /// `SymmetricDifference` is refused rather than silently composed from
    /// three slow-path calls: a caller asking for the fast path should get
    /// it or a clear refusal, not a surprise fallback.
    pub fn boolean_fast(
        &self,
        subject: &TriMesh,
        tool: &TriMesh,
        operation: BooleanOperator,
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        if operation == BooleanOperator::SymmetricDifference {
            return Err(GeomError::Unsupported {
                backend: BoolmeshBoolean::ID,
                operation: axiolid_contracts::Operation::MeshBoolean,
            });
        }
        self.boolean_impl(subject, tool, operation, options, true)
    }

    fn boolean_impl(
        &self,
        subject: &TriMesh,
        tool: &TriMesh,
        operation: BooleanOperator,
        options: &ExecutionOptions,
        fast_winding: bool,
    ) -> GeomResult<BooleanOutcome> {
        // `SymmetricDifference` has no `boolmesh` counterpart. Compose it from
        // the three primitives rather than pretending the backend's operation
        // set is the contract's; `sub_operations` in the evidence records that
        // the caller got a composed result.
        if operation == BooleanOperator::SymmetricDifference {
            return symmetric_difference_via_composition(self, subject, tool, options);
        }

        let op = match operation {
            BooleanOperator::Union => OpType::Add,
            BooleanOperator::Intersection => OpType::Intersect,
            BooleanOperator::Difference => OpType::Subtract,
            BooleanOperator::SymmetricDifference => unreachable!("composed above"),
            // A future operand added to the contract. Refusing is honest and
            // lets the registry fall back to a provider that supports it.
            _ => {
                return Err(GeomError::Unsupported {
                    backend: BoolmeshBoolean::ID,
                    operation: axiolid_contracts::Operation::MeshBoolean,
                })
            }
        };
        options.check_cancelled()?;
        let subject_manifold = to_manifold(subject, "subject")?;
        let tool_manifold = to_manifold(tool, "tool")?;

        let output = match compute_boolean(&subject_manifold, &tool_manifold, op, fast_winding) {
            Ok(output) => output,
            // An empty result is a legitimate value under the contract: the
            // intersection of disjoint solids is nothing. `boolmesh` signals
            // that by failing to build a mesh from an empty vertex matrix, so
            // translate it back into the empty solid the contract specifies
            // rather than surfacing a backend quirk as a geometry error.
            Err(reason) if is_empty_result(&reason) => {
                let mut empty = TriMesh::new(Vec::new(), Vec::new());
                // Zero faces carry every channel vacuously: nothing lost,
                // nothing derived. Keeps a union/fold over it consistent.
                let fates = carry(subject, tool, &mut empty, &[]);
                let evidence =
                    evidence_for(subject, &[tool], &empty, 1).with_attribute_fates(fates);
                return Ok(BooleanOutcome::new(empty, evidence));
            }
            Err(reason) => {
                return Err(GeomError::BackendContractViolation {
                    backend: BoolmeshBoolean::ID,
                    detail: format!("{operation:?} failed inside the solve: {reason}"),
                })
            }
        };

        let mut result = from_boolean_mesh(&output);
        check_result(&result, operation)?;
        let sources: Vec<FaceSource> = output
            .src
            .iter()
            .map(|&(operand, triangle)| FaceSource { operand, triangle })
            .collect();
        let fates = carry(subject, tool, &mut result, &sources);
        let evidence = evidence_for(subject, &[tool], &result, 1).with_attribute_fates(fates);
        Ok(BooleanOutcome::new(result, evidence))
    }
}

/// Whether a `boolmesh` failure actually means "the result is empty".
///
/// Matched on the message because the backend does not expose a typed empty
/// case. Narrow on purpose: any other failure stays a real error rather than
/// being silently converted into an empty solid, which would turn a genuine
/// fault into a plausible-looking zero-volume answer.
fn is_empty_result(reason: &impl std::fmt::Display) -> bool {
    let text = reason.to_string().to_lowercase();
    text.contains("empty pos matrix") || text.contains("empty mesh")
}

/// Build evidence for one completed operation.
///
/// `output_components` counts connected components by union-find over shared
/// vertex indices. A difference that severs a wall reports two components, so
/// the caller learns about the split here rather than in a later quantity.
fn evidence_for(
    subject: &TriMesh,
    tools: &[&TriMesh],
    result: &TriMesh,
    sub_operations: usize,
) -> BooleanEvidence {
    let disjoint = tools
        .iter()
        .filter(|tool| !subject.bounds().intersects(&tool.bounds()))
        .count();
    let evidence = BooleanEvidence::record(
        subject.triangle_count(),
        tools.iter().map(|t| t.triangle_count()).sum(),
        result.triangle_count(),
        axiolid_mesh::component_count(result),
    )
    .with_disjoint_tools(disjoint)
    .with_sub_operations(sub_operations)
    // `boolmesh` does not report coincident-face encounters, so claiming
    // detection would be a lie. Left false until a provider can answer.
    .with_coincident_faces(false)
    // Fallback fate, overwritten by every path that carries channels (all
    // of them since #116). What remains here is the answer when no
    // provenance exists: the `carry` refusal path, and nothing else.
    .with_attribute_fates(
        subject
            .attributes
            .iter()
            .map(|channel| {
                let reason = match channel.blend {
                    axiolid_mesh::Blend::None => axiolid_mesh::DropReason::NotBlendable,
                    _ => axiolid_mesh::DropReason::ProviderLimitation,
                };
                (channel.name.clone(), AttributeFate::Dropped(reason))
            })
            .collect(),
    );

    match relative_overlap(subject, tools) {
        Some(ratio) => evidence.with_relative_overlap(ratio),
        None => evidence,
    }
}

/// Carry channels through the analytic box path.
///
/// The cellular result is rebuilt, not cut, so it has no `Tref`s. Every face
/// lies in a host face (as is) or a cutter face (reversed); plane lookup
/// recovers the source, and the pairwise sampler does the rest. The cutters
/// are fused so the sampler sees one tool, as in the general path.
fn carry_analytic(
    subject: &TriMesh,
    tools: &[TriMesh],
    result: &mut TriMesh,
) -> Vec<(String, AttributeFate)> {
    let members: Vec<&TriMesh> = tools.iter().collect();
    let tool = crate::grouping::fuse_with_channels(&members, &subject.attributes);
    match sources_by_plane(result, &[(subject, 1.0), (&tool, -1.0)]) {
        Some(sources) => carry(subject, &tool, result, &sources),
        // A face with no source means the construction is not what the
        // lookup assumes: attach nothing rather than a wrong value.
        None => dropped_fates(subject, DropReason::ProviderLimitation),
    }
}

/// Every subject channel as dropped, with `reason` (`NotBlendable` wins for
/// a `Blend::None` channel regardless: that is the data's own answer).
fn dropped_fates(subject: &TriMesh, reason: DropReason) -> Vec<(String, AttributeFate)> {
    subject
        .attributes
        .iter()
        .map(|c| {
            let r = match c.blend {
                Blend::None => DropReason::NotBlendable,
                _ => reason,
            };
            (c.name.clone(), AttributeFate::Dropped(r))
        })
        .collect()
}

/// Per-channel fates, in the subject's channel order.
type Fates = Vec<(String, AttributeFate)>;

/// Every subject channel `Preserved`: the starting point a composed path
/// merges each step's fates onto.
fn seed_fates(subject: &TriMesh) -> Fates {
    subject
        .attributes
        .iter()
        .map(|c| (c.name.clone(), AttributeFate::Preserved))
        .collect()
}

/// `mesh` carrying exactly `template`'s channels: its own where one matches
/// (name, width, blend), an all-`UNMAPPED` channel where none does.
///
/// Union is symmetric but channel reporting is keyed on the first operand,
/// so a tree node's left child must carry every channel the batch reports
/// on -- whichever solid it came from. Unmapped keeps the absence honest.
fn conform(mesh: &TriMesh, template: &[AttributeChannel]) -> TriMesh {
    let mut out = mesh.clone();
    out.attributes = crate::grouping::fuse_with_channels(&[mesh], template).attributes;
    out
}

/// Remove every channel whose composed fate is `Dropped`.
fn strip_dropped(mesh: &mut TriMesh, fates: &[(String, AttributeFate)]) {
    mesh.attributes.retain(|c| {
        !fates
            .iter()
            .any(|(n, f)| n == &c.name && matches!(f, AttributeFate::Dropped(_)))
    });
}
/// Smallest relative overlap between the subject and any tool.
///
/// Measured from operand bounds, not from the result: the result cannot show
/// how thin the sliver that produced it was. For each intersecting tool, the
/// overlap box's shortest side is divided by the operand size, giving a
/// scale-free number (ADR 0045).
///
/// Disjoint tools are skipped — they contribute no intersection to condition.
/// `None` when nothing intersects or the operands are degenerate, which is
/// honest: no intersection was constructed, so none was ill conditioned.
fn relative_overlap(subject: &TriMesh, tools: &[&TriMesh]) -> Option<f64> {
    let subject_bounds = subject.bounds();
    let mut worst: Option<f64> = None;

    for tool in tools {
        let tool_bounds = tool.bounds();
        if !subject_bounds.intersects(&tool_bounds) {
            continue;
        }

        // The overlap box: componentwise max of mins, min of maxes.
        let lo = subject_bounds.min.max(tool_bounds.min);
        let hi = subject_bounds.max.min(tool_bounds.max);
        let overlap = hi - lo;

        // Scale: the larger operand's diagonal, so the ratio is relative to
        // the model rather than to whichever operand happens to be smaller.
        let scale = subject_bounds
            .diagonal()
            .length()
            .max(tool_bounds.diagonal().length());
        if !scale.is_finite() || scale <= 0.0 {
            continue;
        }

        // The shortest overlap side governs: a wide, thin sliver is exactly
        // as ill conditioned as its thin dimension makes it.
        let thinnest = overlap.x.min(overlap.y).min(overlap.z);
        if !thinnest.is_finite() {
            continue;
        }
        let ratio = (thinnest / scale).max(0.0);
        worst = Some(worst.map_or(ratio, |w: f64| w.min(ratio)));
    }

    worst
}

/// Batch difference: fuse mutually disjoint cutters, then subtract per group.
///
/// The sequential default runs N booleans, each against a subject that has
/// already been cut N-1 times. This override runs one boolean per GROUP of
/// mutually disjoint cutters, which for the dominant layout (a wall with
/// non-overlapping openings) collapses to a single boolean.
///
/// Correctness rests on `(S \ A) \ B == S \ (A union B)`, plus the fact that a
/// concatenation of disjoint solids IS their union. Grouping uses bounding
/// boxes, which over-separate but never wrongly fuse, so the result is
/// identical to the sequential path -- asserted by the volume gates.
impl BoolmeshBoolean {
    /// Subtract every tool, grouping disjoint ones into single operations.
    pub(crate) fn subtract_grouped(
        &self,
        subject: &TriMesh,
        tools: &[TriMesh],
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        let bounds: Vec<_> = tools.iter().map(TriMesh::bounds).collect();
        let groups = crate::grouping::disjoint_groups(&bounds);

        let mut current = subject.clone();
        let mut sub_operations = 0;
        // Seeded Preserved: the subject enters untouched, and each step's
        // fates compose onto this (`merge_fates`).
        let mut fates = seed_fates(subject);
        for group in &groups {
            // The only real poll point: between groups. Cancelling here returns
            // no mesh at all rather than a partially cut one.
            options.check_cancelled()?;
            sub_operations += 1;
            // A single-member group gains nothing from fusing, so skip the
            // copy and subtract the tool directly.
            let step = if let [only] = group.as_slice() {
                self.difference(&current, &tools[*only], options)?
            } else {
                let members: Vec<&TriMesh> = group.iter().map(|&i| &tools[i]).collect();
                let fused = crate::grouping::fuse_with_channels(&members, &subject.attributes);
                self.difference(&current, &fused, options)?
            };
            fates = merge_fates(&fates, step.evidence.attribute_fates);
            current = step.mesh;
        }
        let borrowed: Vec<&TriMesh> = tools.iter().collect();
        let evidence =
            evidence_for(subject, &borrowed, &current, sub_operations).with_attribute_fates(fates);
        Ok(BooleanOutcome::new(current, evidence))
    }

    /// Union every solid by balanced pairwise reduction.
    ///
    /// The sequential fold the trait default performs makes step `i` union
    /// an accumulator already holding `i` solids against one more, so the
    /// subject is re-walked on every step and total work is quadratic in the
    /// operand count. Reducing in pairs -- union neighbours, then pairs of
    /// those, until one remains -- issues the SAME number of booleans, but
    /// the operands stay small until the final levels.
    ///
    /// Measured on a k^3 grid of icospheres (see the `benchmarks` sibling
    /// repo, `sphere_grid`): 1.4x at 8 solids, 7.8x at 125, 28.9x at 512
    /// (44.4s to 1.5s). The gap widens with n because the difference is
    /// complexity, not a constant factor.
    ///
    /// # Known cost: intermediates are rebuilt at every level
    ///
    /// Each level hands the next one a [`TriMesh`], so every intermediate is
    /// converted back from the internal half-edge form and rebuilt by
    /// `to_manifold` on the level above. Instrumenting `to_manifold` against
    /// `compute_boolean` on a disjoint sphere grid:
    ///
    /// ```text
    ///   n   wall_ms  build_ms  build%  triangles rebuilt vs input
    ///   8       8.9       3.9   43.8%   3.00x
    ///  27      51.2      22.2   43.4%   4.85x
    ///  64     145.7      65.2   44.7%   6.00x
    /// 125     344.1     147.5   42.9%   6.98x
    /// ```
    ///
    /// Construction is ~43% of wall time and the redundancy grows with `n`:
    /// at 125 solids the same triangles are rebuilt seven times. Threading
    /// the built structure through the levels instead would have a ceiling
    /// of about **1.6x** for this path (perfect reuse, zero conversion
    /// cost), so the real figure would be lower.
    ///
    /// It has NOT been done, deliberately. It requires a new internal type
    /// boundary carrying a collider and planar grid between levels, which
    /// makes intermediates heavier in memory, and it touches the code path
    /// with the recorded ordering nondeterminism -- while this path is
    /// currently STABLE at 1000 solids. That is a real property to risk for
    /// a bounded win.
    ///
    /// Note the ceiling applies to BATCH paths only. A single
    /// [`MeshBoolean::boolean`] call has no intermediates -- both operands
    /// are leaves -- so construction there is irreducible and reuse saves
    /// exactly nothing. It is not a lead on the single-boolean gap against
    /// upstream Manifold, which profiling showed to be flat rather than
    /// hotspot-shaped.
    ///
    /// Union is associative and commutative, so any reduction order yields
    /// the same solid; unlike `subtract_grouped` this needs no disjointness
    /// precondition and no fusing, and therefore has no correctness cliff.
    /// Floating-point differences remain -- a differently ordered triangle
    /// list sums to a marginally different volume -- so the gates compare
    /// volumes to a relative tolerance, not bitwise.
    pub(crate) fn union_tree(
        &self,
        solids: &[TriMesh],
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        // The union of nothing is nothing. Returning an empty solid keeps
        // this total rather than making every caller special-case it.
        if solids.is_empty() {
            let borrowed: Vec<&TriMesh> = Vec::new();
            let empty = TriMesh::default();
            let evidence = evidence_for(&empty, &borrowed, &empty, 0);
            return Ok(BooleanOutcome::new(empty, evidence));
        }

        // Every solid carries solids[0]'s channel set, so a channel the
        // result keeps is defined (or explicitly unmapped) on every piece.
        // Each node carries the fates of the path that built it.
        let template = &solids[0].attributes;
        let mut level: Vec<(TriMesh, Fates)> = solids
            .iter()
            .map(|s| (conform(s, template), seed_fates(&solids[0])))
            .collect();
        let mut sub_operations = 0;
        while level.len() > 1 {
            // Every pair at one level is independent: no pair reads another
            // pair's output, so the level is a pure map. That is the whole
            // reason this axis can be threaded while intra-solve threading
            // cannot pay for itself -- there is no coordination inside the
            // level, only a join at the end of it.
            let mut pairs = level.chunks_exact(2);
            let step = |pair: &[(TriMesh, Fates)]| -> GeomResult<(TriMesh, Fates)> {
                let out = self.union(&pair[0].0, &pair[1].0, options)?;
                // Both halves' histories, then this union.
                let history = merge_fates(&pair[0].1, pair[1].1.clone());
                Ok((
                    out.mesh,
                    merge_fates(&history, out.evidence.attribute_fates),
                ))
            };

            // Cancellation is polled ONCE per level rather than per pair.
            // Under rayon the pairs are not ordered, so a per-pair poll
            // would abort at a schedule-dependent point; per level, the
            // cut point is deterministic and the partial work discarded is
            // the same either way.
            options.check_cancelled()?;

            #[cfg(feature = "parallel-batch")]
            let mut next: Vec<(TriMesh, Fates)> = {
                use rayon::prelude::*;
                // `par_iter` over an indexed slice, NOT `par_bridge`: the
                // bridge does not preserve order, which would permute the
                // level and silently change the tree's shape from run to
                // run. An indexed parallel iterator collects positionally,
                // so the output is identical to the sequential path.
                let chunks: Vec<&[(TriMesh, Fates)]> = pairs.clone().collect();
                let merged: Result<Vec<(TriMesh, Fates)>, GeomError> =
                    chunks.par_iter().map(|pair| step(pair)).collect();
                merged?
            };

            #[cfg(not(feature = "parallel-batch"))]
            let mut next: Vec<(TriMesh, Fates)> = {
                let mut acc = Vec::with_capacity(level.len().div_ceil(2));
                for pair in pairs.clone() {
                    acc.push(step(pair)?);
                }
                acc
            };

            sub_operations += next.len();
            // `chunks_exact` was cloned for the merge above, so advance the
            // original to reach its remainder.
            for _ in pairs.by_ref() {}
            // An odd trailing solid rides to the next level uncombined
            // rather than being unioned into an already-merged neighbour,
            // which would unbalance the tree this exists to keep balanced.
            if let Some(last) = pairs.remainder().first() {
                next.push(last.clone());
            }
            level = next;
        }

        let (mut result, fates) = level.into_iter().next().expect("non-empty input");
        // A channel some step dropped is not on the result, or only partly:
        // strip it so the mesh never carries a channel its fate calls lost.
        strip_dropped(&mut result, &fates);
        let borrowed: Vec<&TriMesh> = solids.iter().collect();
        // Evidence names the first operand as the subject and the rest as
        // tools, matching how a caller reads a batch: one solid grown by the
        // others. The reduction order is an implementation detail.
        let (subject, tools) = borrowed.split_first().expect("non-empty input");
        let evidence =
            evidence_for(subject, tools, &result, sub_operations).with_attribute_fates(fates);
        Ok(BooleanOutcome::new(result, evidence))
    }

    /// One union, routed through the same validation as `boolean`.
    fn union(
        &self,
        subject: &TriMesh,
        tool: &TriMesh,
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        self.boolean(subject, tool, BooleanOperator::Union, options)
    }

    /// One difference, routed through the same validation as `boolean`.
    fn difference(
        &self,
        subject: &TriMesh,
        tool: &TriMesh,
        options: &ExecutionOptions,
    ) -> GeomResult<BooleanOutcome> {
        self.boolean(subject, tool, BooleanOperator::Difference, options)
    }
}

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

    /// A tetrahedron with reversed winding: the shape a faulty backend would
    /// return if it inverted its output.
    fn inside_out_tetrahedron() -> TriMesh {
        let positions = vec![
            Point3::new(0.0, 0.0, 0.0),
            Point3::new(1.0, 0.0, 0.0),
            Point3::new(0.0, 1.0, 0.0),
            Point3::new(0.0, 0.0, 1.0),
        ];
        // Verified: signed volume -1/6, i.e. inward-facing normals.
        let indices = vec![0, 1, 2, 0, 3, 1, 0, 2, 3, 1, 3, 2];
        TriMesh::new(positions, indices)
    }

    fn positive_infinite_volume() -> TriMesh {
        let large = 1.0e308;
        let small = 1.0e-308;
        TriMesh::new(
            vec![
                Point3::new(small, small, 1.0),
                Point3::new(large, 1.0, small),
                Point3::new(small, large, 1.0),
            ],
            vec![0, 1, 2],
        )
    }

    /// `check_result` guards against an upstream regression that inverts its
    /// output. With validated inputs the current `boolmesh` release never does
    /// this -- verified by instrumenting the branch across the whole suite and
    /// observing zero hits -- so the guard is exercised directly here rather
    /// than left as untested defensive code.
    #[test]
    fn an_inside_out_result_is_blamed_on_the_backend() {
        let error = check_result(&inside_out_tetrahedron(), BooleanOperator::Difference)
            .expect_err("an inside-out result must be rejected");
        match error {
            GeomError::BackendContractViolation { backend, detail } => {
                assert_eq!(backend, BoolmeshBoolean::ID);
                assert!(detail.contains("inside-out"), "{detail}");
            }
            other => panic!("must blame the backend, not the caller: {other:?}"),
        }
    }

    #[test]
    fn a_non_finite_result_is_blamed_on_the_backend() {
        let error = check_result(&positive_infinite_volume(), BooleanOperator::Union)
            .expect_err("a non-finite result volume must be rejected");
        assert!(
            matches!(error, GeomError::BackendContractViolation { backend, ref detail }
                if backend == BoolmeshBoolean::ID && detail.contains("non-finite signed volume")),
            "must blame the backend, got {error:?}"
        );
    }

    /// An empty result is legitimate (tool fully contains subject) and must not
    /// be mistaken for an orientation fault.
    #[test]
    fn an_empty_result_is_accepted() {
        assert!(check_result(&TriMesh::default(), BooleanOperator::Difference).is_ok());
    }

    /// A correctly oriented result passes.
    #[test]
    fn an_outward_result_is_accepted() {
        let mut mesh = inside_out_tetrahedron();
        for corner in mesh.indices.chunks_exact_mut(3) {
            corner.swap(1, 2);
        }
        assert!(check_result(&mesh, BooleanOperator::Difference).is_ok());
    }
}