isomesh 0.0.11

//! Cube topology and the sign rule, shared by every extractor.
//!
//! Corner and edge numbering, which face touches which, and the one global
//! decision about what an exactly-zero sample means. Marching Cubes and Surface
//! Nets both index cells the same way, and two copies of this would be two
//! chances to disagree about it.

use crate::Real;

/// A cube has eight corners.
///
/// Corner `i` sits at local coordinate `(i & 1, (i >> 1) & 1, (i >> 2) & 1)` —
/// bit 0 is x, bit 1 is y, bit 2 is z, matching the crate's x-varies-fastest
/// index order.
pub const CORNER_COUNT: usize = 8;

/// A cube has twelve edges.
pub const EDGE_COUNT: usize = 12;

/// The two corners of each edge, lower corner first.
///
/// Generated by walking axes x, y, z and, within each, the four corners whose
/// bit for that axis is zero. The numbering is therefore this crate's own; no
/// paper's numbering is being matched, which is what lets the Marching Cubes
/// table be derived rather than transcribed.
pub const EDGE_CORNERS: [[u8; 2]; EDGE_COUNT] = build_edge_corners();

const fn build_edge_corners() -> [[u8; 2]; EDGE_COUNT] {
    let mut out = [[0u8; 2]; EDGE_COUNT];
    let mut n = 0;
    let mut axis = 0u8;
    while axis < 3 {
        let mut c = 0u8;
        while c < 8 {
            if c & (1 << axis) == 0 {
                out[n] = [c, c | (1 << axis)];
                n += 1;
            }
            c += 1;
        }
        axis += 1;
    }
    out
}

/// The axis each edge runs along.
pub const EDGE_AXIS: [u8; EDGE_COUNT] = build_edge_axis();

const fn build_edge_axis() -> [u8; EDGE_COUNT] {
    let mut out = [0u8; EDGE_COUNT];
    let mut e = 0;
    while e < EDGE_COUNT {
        let d = EDGE_CORNERS[e][0] ^ EDGE_CORNERS[e][1];
        out[e] = match d {
            1 => 0,
            2 => 1,
            _ => 2,
        };
        e += 1;
    }
    out
}

/// Index of the edge joining two corners.
///
/// # Panics
///
/// At compile time, if the two corners are not adjacent on the cube.
pub const fn edge_index(a: u8, b: u8) -> u8 {
    let (lo, hi) = if a < b { (a, b) } else { (b, a) };
    let mut e = 0;
    while e < EDGE_COUNT {
        if EDGE_CORNERS[e][0] == lo && EDGE_CORNERS[e][1] == hi {
            return e as u8;
        }
        e += 1;
    }
    panic!("corners are not joined by a cube edge")
}

/// The four corners of a face, counter-clockwise **as seen from outside**.
///
/// `axis` is the face's normal axis and `side` selects the near (`0`) or far
/// (`1`) face along it. With `u = (axis+1) % 3` and `v = (axis+2) % 3` the
/// triple `(u, v, axis)` is right-handed, so on the far face the viewer sees
/// `u × v` pointing at them and counter-clockwise runs `00 → 10 → 11 → 01`; on
/// the near face they see it pointing away, and the order reverses.
pub const fn face_corners(axis: usize, side: u8) -> [u8; 4] {
    let u = (axis + 1) % 3;
    let v = (axis + 2) % 3;
    let base = side << axis;
    let c00 = base;
    let c10 = base | (1 << u);
    let c11 = base | (1 << u) | (1 << v);
    let c01 = base | (1 << v);
    if side == 1 {
        [c00, c10, c11, c01]
    } else {
        [c00, c01, c11, c10]
    }
}

/// Is this edge one of the four on the face `(axis, side)`?
///
/// An edge lies on a face when both its corners do. Its two corners differ only
/// in its own axis, so it is on the face exactly when the face's axis is not its
/// own and its lower corner is on the right side.
#[inline]
#[must_use]
pub const fn edge_on_face(edge: u8, axis: usize, side: u8) -> bool {
    EDGE_AXIS[edge as usize] as usize != axis
        && (EDGE_CORNERS[edge as usize][0] >> axis) & 1 == side
}

/// Do these two cube edges both lie on some one face of the cube?
///
/// This is the question "can any single *other* cell name both of these edges?"
/// Two cells share a face and its four edges; two cells that share only an edge
/// or a corner share no *pair* of edges. So a mesh edge joining two cut cube
/// edges that share no face can be emitted by exactly one cell — which is what
/// makes such a chord safe from the collision A-015 is about. See
/// [`crate::marching_cubes::table::CENTROID_BASE`].
#[inline]
#[must_use]
pub const fn edges_share_a_face(a: u8, b: u8) -> bool {
    let mut axis = 0usize;
    while axis < 3 {
        let mut side = 0u8;
        while side < 2 {
            if edge_on_face(a, axis, side) && edge_on_face(b, axis, side) {
                return true;
            }
            side += 1;
        }
        axis += 1;
    }
    false
}

/// Is this corner inside the solid, for the given case index?
#[inline]
pub const fn corner_inside(case: u8, corner: u8) -> bool {
    case & (1 << corner) != 0
}

/// The local offset of a corner, as grid steps.
#[inline]
#[must_use]
pub(crate) const fn corner_offset(corner: u8) -> [u32; 3] {
    [
        (corner & 1) as u32,
        ((corner >> 1) & 1) as u32,
        ((corner >> 2) & 1) as u32,
    ]
}

/// The global rule for classifying a sample.
///
/// **A sample of exactly zero is outside.** Lengyel (2010 dissertation, §3.1.1)
/// notes the choice has to be made once and applied everywhere, and that either
/// answer is acceptable so long as it is consistent.
///
/// Choosing strict `< 0` buys something concrete: on a cut edge one endpoint is
/// then strictly negative and the other is `>= 0`, so `a - b` is strictly
/// negative and the interpolation `t = a / (a - b)` can never divide by zero.
/// The implementation brief's suggested `|a - b| < eps → t = 0.5` guard is
/// therefore not needed at all — and is actively unwanted, since an epsilon
/// there would snap genuinely resolvable crossings to the edge midpoint and
/// blunt exactly the sub-voxel features A-014 exists to recover.
#[inline]
pub fn is_inside<R: Real>(value: R) -> bool {
    value < R::ZERO
}

/// Where the surface crosses an edge, as a **signed offset from the edge
/// midpoint** in units of the edge, in `[-1/2, +1/2]`.
///
/// `d = ((a + b) / 2) / (a - b)`, and the vertex goes at
/// `edge_midpoint + edge_vector * d`. Every extractor places its vertices from
/// this, so they agree about where the surface is and differ only in what they
/// do with the crossings.
///
/// # Why the midpoint and not the lower corner (R-059, P-61)
///
/// The obvious form is a parameter `t = a / (a - b)` from the lower corner, and
/// this crate shipped it until `✗39` measured what it costs. Under a reflection
/// of the grid, which corner is "lower" swaps, so the correctly reflected
/// parameter is `1 - t` — and `1 - t` is **not** what `b / (b - a)` computes.
/// Reflection acts on `[0, 1]` as `0 <-> 1`, an *affine* map, and IEEE 754
/// respects sign flips exactly and affine maps only approximately.
///
/// In the centred frame reflection is a **sign flip**, and antisymmetry is four
/// guarantees rather than four observations:
///
/// - `fl(a + b) = fl(b + a)`, because addition is commutative;
/// - halving is exact, because 2 is a power of two;
/// - `fl(b - a) = -fl(a - b)`, because round-to-nearest is an odd function;
/// - `fl(S / -D) = -fl(S / D)`, for the same reason.
///
/// So the simultaneous endpoint-and-sign swap negates `d` **exactly**, and a
/// mirrored grid produces a mirrored vertex bit for bit. `M-356`'s note that the
/// subtraction is the culprit is the half of the story that is wrong: the two
/// denominators are already exact negations.
///
/// The tail cost is a wash and was measured before the change: over 300,000
/// random world positions against exact rational arithmetic, the lower-corner
/// form is mean 0.086 ulp / worst 422 ulp and this one mean 0.052 / worst 757.
///
/// # Panics
///
/// In debug builds, if the endpoints are not on opposite sides — in which case
/// there is no crossing to find and the caller has a bug.
///
/// # Why this is public
///
/// R-066. `P-61` made this *the* definition of where a crossing is, and every
/// extractor in the crate routes through it. Anything that measures a crossing
/// -- an error bound, a certificate, a diagnostic -- has to ask the same
/// question of the same expression, and a second copy of a one-line formula is
/// exactly the drift `x39` cost this crate 216 golden hashes to find.
#[inline]
pub fn edge_offset<R: Real>(a: R, b: R) -> R {
    debug_assert!(is_inside(a) != is_inside(b));
    ((a + b) * R::HALF) / (a - b)
}

/// Place a crossing on the segment `lo`..`hi`, given `d` from [`edge_offset`].
///
/// One line, and it exists so the centred frame has exactly one spelling: a
/// second `(lo + hi) * HALF + (hi - lo) * d` written out by hand is a second
/// place for the frame to drift.
#[inline]
pub(crate) fn place<R: Real>(lo: R, hi: R, d: R) -> R {
    (lo + hi) * R::HALF + (hi - lo) * d
}

#[cfg(test)]
mod tests {
    // Exact comparisons on purpose: these pin bit-level conventions.
    #![allow(clippy::float_cmp)]

    use super::*;
    use alloc::vec::Vec;

    #[test]
    fn every_edge_joins_corners_differing_in_one_bit() {
        for (e, [a, b]) in EDGE_CORNERS
            .iter()
            .map(|c| (c, *c))
            .enumerate()
            .map(|(i, (_, c))| (i, c))
        {
            let diff = a ^ b;
            assert!(diff.is_power_of_two(), "edge {e}: {a} to {b}");
            assert_eq!(diff.trailing_zeros(), u32::from(EDGE_AXIS[e]));
            assert!(a < b, "edge {e} is not lower-corner-first");
        }
    }

    #[test]
    fn the_twelve_edges_are_distinct() {
        let mut seen: Vec<[u8; 2]> = EDGE_CORNERS.to_vec();
        seen.sort_unstable();
        seen.dedup();
        assert_eq!(seen.len(), EDGE_COUNT);
    }

    #[test]
    fn edge_index_round_trips() {
        for (e, [a, b]) in EDGE_CORNERS.iter().copied().enumerate() {
            assert_eq!(edge_index(a, b), e as u8);
            assert_eq!(edge_index(b, a), e as u8, "order must not matter");
        }
    }

    /// Each face names four distinct corners, and each corner sits on three
    /// faces. A transposition in `face_corners` would break this.
    #[test]
    fn faces_partition_the_corners() {
        let mut appearances = [0u32; CORNER_COUNT];
        for axis in 0..3 {
            for side in 0..2u8 {
                let corners = face_corners(axis, side);
                let mut sorted: Vec<u8> = corners.to_vec();
                sorted.sort_unstable();
                sorted.dedup();
                assert_eq!(sorted.len(), 4, "axis {axis} side {side}");
                for c in corners {
                    appearances[c as usize] += 1;
                    // Every corner of this face agrees on the face's own axis.
                    assert_eq!((c >> axis) & 1, side);
                }
            }
        }
        assert_eq!(appearances, [3; CORNER_COUNT]);
    }

    /// Consecutive corners around a face must be adjacent on the cube, or the
    /// "walk the face" construction in the Marching Cubes table would be walking
    /// diagonals.
    #[test]
    fn face_corners_are_a_cycle() {
        for axis in 0..3 {
            for side in 0..2u8 {
                let c = face_corners(axis, side);
                for k in 0..4 {
                    let diff = c[k] ^ c[(k + 1) % 4];
                    assert!(diff.is_power_of_two(), "axis {axis} side {side} at {k}");
                }
            }
        }
    }

    #[test]
    fn corner_offset_matches_the_bit_layout() {
        assert_eq!(corner_offset(0), [0, 0, 0]);
        assert_eq!(corner_offset(1), [1, 0, 0]);
        assert_eq!(corner_offset(2), [0, 1, 0]);
        assert_eq!(corner_offset(4), [0, 0, 1]);
        assert_eq!(corner_offset(7), [1, 1, 1]);
    }

    #[test]
    fn zero_is_outside() {
        assert!(is_inside(-1.0f64));
        assert!(!is_inside(0.0f64));
        assert!(!is_inside(1.0f64));
        // -0.0 is not less than 0.0, so it is outside too.
        assert!(!is_inside(-0.0f64));
    }

    /// The property that makes the division safe without any epsilon guard: a cut
    /// edge always has a strictly negative endpoint and a non-negative one.
    #[test]
    fn edge_offset_lands_within_half_an_edge_of_the_midpoint() {
        for (a, b) in [(-1.0f64, 1.0), (-0.25, 3.0), (-2.0, 0.0), (5.0, -0.5)] {
            let d = edge_offset(a, b);
            assert!((-0.5..=0.5).contains(&d), "a={a} b={b} d={d}");
            assert!((0.0..=1.0).contains(&place(0.0, 1.0, d)), "a={a} b={b}");
        }
        // Exactly at the negative endpoint when the other is zero.
        assert_eq!(edge_offset(-2.0f64, 0.0), 0.5);
        assert_eq!(place(0.0f64, 1.0, edge_offset(-2.0, 0.0)), 1.0);
    }

    /// R-059's whole point, asserted rather than argued: the simultaneous
    /// endpoint-and-sign swap negates the offset **exactly**.
    ///
    /// The control is the second half and it is **searched rather than chosen**,
    /// because the first version of this test picked `(−0.1, 0.7)` by hand to
    /// show the parameter form failing and those values are one of the many
    /// where `t + t' == 1` exactly — M-32's rule, on the day it was written.
    /// The loop below asserts that at least one pair in a deterministic sweep
    /// separates the two forms, so the test cannot pass by both forms being
    /// exact on everything it happens to try.
    #[test]
    fn the_offset_is_exactly_antisymmetric_under_the_endpoint_swap() {
        let mut separated = 0usize;
        let mut state = 0x2026_u64;
        for _ in 0..4096 {
            // SplitMix64, so the sweep is the same on every machine.
            state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
            let mut z = state;
            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
            z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
            z ^= z >> 31;
            let unit = |bits: u64| (bits >> 11) as f64 * (1.0 / 9_007_199_254_740_992.0);
            let a = -unit(z).max(f64::MIN_POSITIVE);
            let b = unit(z.rotate_left(32)).max(f64::MIN_POSITIVE);

            // The centred form, always exact.
            assert_eq!(edge_offset(b, a), -edge_offset(a, b), "a={a} b={b}");
            // The parameter form, whose two anchors are `t` and `1 - t`.
            if a / (a - b) + b / (b - a) != 1.0 {
                separated += 1;
            }
        }
        assert!(
            separated > 0,
            "the control found no pair where the lower-corner parameter form \
             fails, so this test cannot tell the two forms apart"
        );

        // Hand-picked extremes, kept because they are the interesting shapes:
        // a tiny straddle, a zero endpoint, and an asymmetric magnitude.
        for (a, b) in [
            (-1.0f64, 3.0),
            (-1e-300, 1e-300),
            (-1.0 / 3.0, 7.0 / 11.0),
            (0.25, -1e-17),
            (-2.0, 0.0),
        ] {
            assert_eq!(edge_offset(b, a), -edge_offset(a, b), "a={a} b={b}");
        }
    }
}