indicatrix 0.5.0

//! Enumeration and deduplication of a plane arrangement's vertices: the
//! shared primitive [`measure_solid`](super::measure_solid) and
//! [`build_solid_mesh`](super::build_solid_mesh) both build on.

use glam::DVec3;

use super::{BLANK_HALF_EXTENT, EPS_FEAS, MIN_TRIPLE_DET, VERTEX_DEDUP};

/// One deduplicated vertex of the solid.
pub(super) struct SolidVertex {
    pub(super) v: DVec3,
}

/// Accepted vertices, in first-seen order, plus an index over them sorted by
/// `x` so a new candidate's duplicate check only has to scan the vertices
/// that could possibly be within [`VERTEX_DEDUP`] of it on every axis instead
/// of the full accepted set.
///
/// `verts`' order is exactly the push order of [`insert_if_new`](Self::insert_if_new)
/// calls (the divergence-theorem sum in [`measure_solid`](super::measure_solid)
/// depends on it); `by_x` is a lookup structure only, never observed by
/// callers.
#[derive(Default)]
struct VertexAccumulator {
    verts: Vec<SolidVertex>,
    /// Indices into `verts`, kept sorted ascending by `verts[i].v.x`.
    by_x: Vec<usize>,
}

impl VertexAccumulator {
    /// Inserts `v` unless some already-accepted vertex is within
    /// [`VERTEX_DEDUP`] of it on every axis -- identical to a linear scan
    /// testing `(s.v - v).abs().max_element() < VERTEX_DEDUP` against every
    /// prior vertex, just restricted up front to the `x`-sorted window that
    /// could possibly match (any vertex outside `[v.x - VERTEX_DEDUP, v.x +
    /// VERTEX_DEDUP]` fails the `x`-axis check alone, so narrowing to that
    /// window changes no accept/reject decision).
    fn insert_if_new(&mut self, v: DVec3) {
        let verts = &self.verts;
        let lo = self
            .by_x
            .partition_point(|&i| verts[i].v.x < v.x - VERTEX_DEDUP);
        let hi = self
            .by_x
            .partition_point(|&i| verts[i].v.x <= v.x + VERTEX_DEDUP);
        for &i in &self.by_x[lo..hi] {
            if (self.verts[i].v - v).abs().max_element() < VERTEX_DEDUP {
                return;
            }
        }
        let idx = self.verts.len();
        self.verts.push(SolidVertex { v });
        let pos = self.by_x.partition_point(|&i| self.verts[i].v.x < v.x);
        self.by_x.insert(pos, idx);
    }
}

/// Drops duplicate planes (same normal and offset within tight tolerance) so a
/// tier that lists the same index twice can't double-count its face's area.
pub(super) fn dedup_planes(planes: &[(DVec3, f64)]) -> Vec<(DVec3, f64)> {
    let mut out: Vec<(DVec3, f64)> = Vec::with_capacity(planes.len());
    for &(n, m) in planes {
        let dup = out
            .iter()
            .any(|&(n2, m2)| n.dot(n2) > 1.0 - 1e-12 && (m - m2).abs() < 1e-9);
        if !dup {
            out.push((n, m));
        }
    }
    out
}

/// Drains one full (or final partial) [`crate::simd::TripleBatch`] solve into
/// `verts`, in ascending lane order. Returns `None` the moment a vertex
/// escapes to the blank box, propagated by the caller via `?` -- matching the
/// original loop's immediate `return None`. Shared by
/// [`feasible_vertices`]'s batching loop.
fn flush_solid_batch(
    batch: &crate::simd::TripleBatch,
    soa: &crate::simd::PlanesSoA64,
    acc: &mut VertexAccumulator,
) -> Option<()> {
    let sol = crate::simd::solve_triple_batch(batch);
    for lane in 0..batch.len {
        if sol.det[lane].abs() < MIN_TRIPLE_DET {
            continue;
        }
        let v = DVec3::new(sol.vx[lane], sol.vy[lane], sol.vz[lane]);
        if v.abs().max_element() > BLANK_HALF_EXTENT + 1.0 {
            continue;
        }
        if crate::simd::any_violation(soa, v, EPS_FEAS) {
            continue;
        }
        // A feasible vertex at the blank box means the real planes never
        // closed the solid up -- there is no finite stone to measure.
        if v.abs().max_element() > BLANK_HALF_EXTENT - 1.0 {
            return None;
        }
        acc.insert_if_new(v);
    }
    Some(())
}

/// Enumerates the solid's distinct vertices: every well-conditioned plane triple
/// whose intersection satisfies all half-spaces (within [`EPS_FEAS`]), then
/// deduplicated by position. Returns `None` when any vertex reaches the bounding
/// blank (the real planes don't bound a finite solid).
///
/// Batched through `crate::simd`, matching
/// `meet_solver::enumerate_candidate_vertices`: one `PlanesSoA64` built up
/// front (owner is irrelevant to this owner-free scan, so every plane is
/// pushed with owner 0), triples solved via `solve_triple_batch` via
/// [`flush_solid_batch`], and the `any()` feasibility scan replaced by
/// `any_violation` -- bit-identical per lane to the `glam` `DMat3` sequence
/// and scalar scan they replace (see `src/simd.rs`'s determinism contract).
/// Lanes are drained in ascending order and triples are still generated by
/// the same nested loops in the same order, so vertex order and every
/// decision here (determinant check, bounds check, feasibility,
/// blank-escape) match the unbatched scalar version exactly.
pub(super) fn feasible_vertices(planes: &[(DVec3, f64)]) -> Option<Vec<SolidVertex>> {
    let mut all: Vec<(DVec3, f64)> = planes.to_vec();
    for n in [
        DVec3::X,
        DVec3::NEG_X,
        DVec3::Y,
        DVec3::NEG_Y,
        DVec3::Z,
        DVec3::NEG_Z,
    ] {
        all.push((n, BLANK_HALF_EXTENT));
    }

    let mut soa = crate::simd::PlanesSoA64::with_capacity(all.len());
    for &(n, m) in &all {
        soa.push(n, m, 0);
    }

    let p = all.len();
    let mut acc = VertexAccumulator::default();
    let mut batch = crate::simd::TripleBatch::default();
    for a in 0..p {
        for b in (a + 1)..p {
            for c in (b + 1)..p {
                let (pa, pb, pc) = (all[a], all[b], all[c]);
                if batch.push((pa.0, pa.1), (pb.0, pb.1), (pc.0, pc.1)) {
                    flush_solid_batch(&batch, &soa, &mut acc)?;
                    batch = crate::simd::TripleBatch::default();
                }
            }
        }
    }
    if batch.len > 0 {
        flush_solid_batch(&batch, &soa, &mut acc)?;
    }
    Some(acc.verts)
}

/// Maps each plane of `deduped` (in order) back to its index in `original`.
///
/// Valid because [`dedup_planes`] is a pure filter: it copies each retained
/// plane through unchanged (no numeric transform) and never reorders
/// anything, only drops exact repeats -- so `deduped` is, element for
/// element, an order-preserving subsequence of `original`. A single lockstep
/// scan (advancing `original`'s cursor past every plane it consumes, whether
/// kept or skipped) recovers the mapping without re-implementing
/// `dedup_planes`'s own equality test, and without ever revisiting an
/// `original` index already assigned to an earlier `deduped` entry (so two
/// bit-identical original planes are told apart by position, not silently
/// both mapped to the first).
pub(super) fn dedup_origin_indices(
    original: &[(DVec3, f64)],
    deduped: &[(DVec3, f64)],
) -> Vec<usize> {
    let same = |a: (DVec3, f64), b: (DVec3, f64)| {
        a.0.x.to_bits() == b.0.x.to_bits()
            && a.0.y.to_bits() == b.0.y.to_bits()
            && a.0.z.to_bits() == b.0.z.to_bits()
            && a.1.to_bits() == b.1.to_bits()
    };
    let mut out = Vec::with_capacity(deduped.len());
    let mut cursor = 0usize;
    for &d in deduped {
        while cursor < original.len() && !same(original[cursor], d) {
            cursor += 1;
        }
        // `cursor == original.len()` here would mean `deduped` contains a
        // plane `dedup_planes` could not have produced from `original` --
        // an invariant violation, not a real runtime case. Clamping instead
        // of panicking keeps this diagnostic helper infallible even if that
        // invariant is ever broken by a future edit to `dedup_planes`.
        out.push(cursor.min(original.len().saturating_sub(1)));
        cursor += 1;
    }
    out
}

/// Diagnostic re-scan of `planes` (already deduped), used only when
/// [`feasible_vertices`] reports the arrangement unbounded: replays the same
/// augmented-triple enumeration and the same escape test
/// (`flush_solid_batch`'s `v.abs().max_element() > BLANK_HALF_EXTENT - 1.0`),
/// but instead of stopping at the first escaping vertex, visits every triple
/// and collects which of the REAL (non-blank) plane indices participate in
/// at least one.
///
/// Not SIMD-batched, unlike [`feasible_vertices`]: this only runs on an
/// editor's invalid intermediate states, never on the success path measured
/// by any perf budget, so a plain `glam` solve per triple (the same
/// arithmetic `feasible_vertices`'s batches compute, just one triple at a
/// time) is the right tradeoff here -- obviously correct beats fast.
pub(super) fn escaping_plane_indices(planes: &[(DVec3, f64)]) -> Vec<usize> {
    let real_count = planes.len();
    let mut all: Vec<(DVec3, f64)> = planes.to_vec();
    for n in [
        DVec3::X,
        DVec3::NEG_X,
        DVec3::Y,
        DVec3::NEG_Y,
        DVec3::Z,
        DVec3::NEG_Z,
    ] {
        all.push((n, BLANK_HALF_EXTENT));
    }

    let p = all.len();
    let mut escaping: Vec<usize> = Vec::new();
    for a in 0..p {
        for b in (a + 1)..p {
            for c in (b + 1)..p {
                let (na, ma) = all[a];
                let (nb, mb) = all[b];
                let (nc, mc) = all[c];
                let mat = glam::DMat3::from_cols(na, nb, nc).transpose();
                let det = mat.determinant();
                if det.abs() < MIN_TRIPLE_DET {
                    continue;
                }
                let v = mat.inverse() * DVec3::new(ma, mb, mc);
                if v.abs().max_element() > BLANK_HALF_EXTENT + 1.0 {
                    continue;
                }
                if all.iter().any(|&(n, m)| n.dot(v) - m > EPS_FEAS) {
                    continue;
                }
                if v.abs().max_element() > BLANK_HALF_EXTENT - 1.0 {
                    for idx in [a, b, c] {
                        if idx < real_count && !escaping.contains(&idx) {
                            escaping.push(idx);
                        }
                    }
                }
            }
        }
    }
    escaping.sort_unstable();
    escaping
}