isomesh 0.0.9

//! 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 fraction from `a` to `b`.
///
/// Linear interpolation of the field along the edge: `t = a / (a - b)`. Both
/// extractors place their vertices from this, so they agree about where the
/// surface is and differ only in what they do with the crossings.
///
/// # 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.
#[inline]
pub(crate) fn edge_crossing<R: Real>(a: R, b: R) -> R {
    debug_assert!(is_inside(a) != is_inside(b));
    a / (a - b)
}

#[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 interpolation division safe without any
    /// epsilon guard: a cut edge always has a strictly negative endpoint and a
    /// non-negative one.
    #[test]
    fn edge_crossing_lands_between_the_endpoints() {
        for (a, b) in [(-1.0f64, 1.0), (-0.25, 3.0), (-2.0, 0.0), (5.0, -0.5)] {
            let t = edge_crossing(a, b);
            assert!((0.0..=1.0).contains(&t), "a={a} b={b} t={t}");
        }
        // Exactly at the negative endpoint when the other is zero.
        assert_eq!(edge_crossing(-2.0f64, 0.0), 1.0);
    }
}