axiolid-mesh-boolean-boolmesh 0.3.2

//--- Copyright (C) 2025 Saki Komikado <komietty@gmail.com>,
//--- This Source Code Form is subject to the terms of the Mozilla Public License v.2.0.
#![allow(clippy::needless_range_loop)]

use crate::csg::{Real, Vec2u, Vec3, Vec3u};
#[cfg(feature = "parallel")]
use rayon::prelude::*;
use std::f64::consts::PI;

/// Hmesh preserves the order of pos and idx in any cases.
/// Edges are ordered so as the edge is forward (tail idx < head idx)
#[derive(Debug, Clone)]
pub(in crate::csg::manifold) struct Hmesh {
    pub nv: usize,
    pub nf: usize,
    pub nh: usize,
    pub twin: Vec<usize>,
    pub head: Vec<usize>,
    pub tail: Vec<usize>,
    pub vns: Vec<Vec3>,
    pub fns: Vec<Vec3>,
}

fn edge_topology(
    pos: &[Vec3],
    idx: &[Vec3u],
    e2v: &mut Vec<Vec2u>,
    e2f: &mut Vec<Vec2u>,
    f2e: &mut Vec<Vec3u>,
) -> Result<(), String> {
    if pos.is_empty() {
        return Err("empty pos matrix".into());
    }
    if idx.is_empty() {
        return Err("empty idx matrix".into());
    }

    // Exactly three entries per face, known up front: reserving avoids the
    // ~log2(3n) reallocations and copies a push-grown Vec pays.
    let mut ett: Vec<[usize; 4]> = Vec::with_capacity(idx.len() * 3);

    for (i, idx_) in idx.iter().enumerate() {
        for j in 0..3 {
            let mut v1 = idx_[j];
            let mut v2 = idx_[(j + 1) % 3];
            if v1 > v2 {
                std::mem::swap(&mut v1, &mut v2);
            }
            ett.push([v1, v2, i, j]);
        }
    }
    // Entries end with (face, corner), which is unique per element, so no
    // two rows compare equal and stability is vacuous.
    ett.sort_unstable();

    // Walk the sorted table in runs sharing an (v1, v2) key. One run is
    // one edge: a single row is a border edge, two rows an interior one.
    //
    // Three or more rows means three or more faces meet on that edge,
    // which this structure cannot represent -- `e2f` holds exactly two
    // face slots. The previous count collapsed any run to one edge while
    // the fill loop below emits ceil(k/2) of them, so a non-manifold edge
    // overran the allocation and aborted the process (kernel#102).
    //
    // Refusing is deliberate rather than widening the count to match:
    // making the arithmetic agree would keep the first two faces, drop
    // the rest, and return a plausible-looking mesh that silently lost
    // geometry. A caller cannot detect that; it can handle an error.
    let mut ne = 0;
    let mut i = 0;
    while i < ett.len() {
        let mut j = i + 1;
        while j < ett.len() && ett[j][0] == ett[i][0] && ett[j][1] == ett[i][1] {
            j += 1;
        }
        if j - i > 2 {
            return Err(format!(
                "edge ({}, {}) has {} incident faces; a half-edge mesh admits at most two",
                ett[i][0],
                ett[i][1],
                j - i
            ));
        }
        ne += 1;
        i = j;
    }

    e2v.resize(ne, Vec2u::MAX);
    e2f.resize(ne, Vec2u::MAX);
    f2e.resize(idx.len(), Vec3u::MAX);
    ne = 0;

    let mut i = 0;
    while i < ett.len() {
        if i == ett.len() - 1 || !((ett[i][0] == ett[i + 1][0]) && (ett[i][1] == ett[i + 1][1])) {
            // Border edge
            let [v1, v2, i, j] = ett[i];
            e2v[ne][0] = v1;
            e2v[ne][1] = v2;
            e2f[ne][0] = i;
            f2e[i][j] = ne;
        } else {
            let r1 = ett[i];
            let r2 = ett[i + 1];
            e2v[ne][0] = r1[0];
            e2v[ne][1] = r1[1];
            e2f[ne][0] = r1[2];
            e2f[ne][1] = r2[2];
            f2e[r1[2]][r1[3]] = ne;
            f2e[r2[2]][r2[3]] = ne;
            i += 1; // skip the next one
        }
        ne += 1;
        i += 1;
    }

    for i in 0..e2f.len() {
        let fid = e2f[i][0];
        let mut flip = true;
        for j in 0..3 {
            if idx[fid][j] == e2v[i][0] && idx[fid][(j + 1) % 3] == e2v[i][1] {
                flip = false;
            }
        }

        if flip {
            let tmp = e2f[i][0];
            e2f[i][0] = e2f[i][1];
            e2f[i][1] = tmp;
        }
    }
    Ok(())
}

impl Hmesh {
    pub fn new(pos: &[Vec3], idx: &[Vec3u]) -> Result<Self, String> {
        let mut e2v = Default::default();
        let mut e2f = Default::default();
        let mut f2e = Default::default();
        edge_topology(pos, idx, &mut e2v, &mut e2f, &mut f2e)?;

        let nv = pos.len();
        let nf = idx.len();
        let ne = e2v.len();
        let nh = e2v.len() * 2;
        let np = 3;
        let mut v2h = vec![usize::MAX; nv];
        let mut e2h = vec![usize::MAX; ne];
        let mut f2h = vec![usize::MAX; nf];
        let mut next = vec![usize::MAX; nh];
        let mut prev = vec![usize::MAX; nh];
        let mut twin = vec![usize::MAX; nh];
        let mut head = vec![usize::MAX; nh];
        let mut tail = vec![usize::MAX; nh];
        let mut edge = vec![usize::MAX; nh];
        let mut face = vec![usize::MAX; nh];

        for it in 0..nf {
            for ip in 0..np {
                let ih_bgn = it * np;
                let iv = idx[it][ip];
                let ie = f2e[it][ip];
                let ih = ih_bgn + ip;
                next[ih] = ih_bgn + (ip + 1) % np;
                prev[ih] = ih_bgn + (ip + np - 1) % np;
                head[ih] = idx[it][(ip + 1) % np];
                tail[ih] = iv;
                edge[ih] = ie;
                face[ih] = it;
                if f2h[it] == usize::MAX {
                    f2h[it] = ih;
                }
                if v2h[iv] == usize::MAX {
                    v2h[iv] = ih;
                }
                if e2h[ie] == usize::MAX {
                    e2h[ie] = ih;
                } else {
                    twin[ih] = e2h[ie];
                    twin[e2h[ie]] = ih;
                }
            }
        }

        if twin.iter().any(|v| v == &usize::MAX) {
            return Err("Input mesh must not contain boundary edges.".into());
        }

        // `half` used to be materialised here as an identity permutation
        // over 0..nh, then immediately walked in that same order by the
        // caller. Removed: the caller now just uses 0..nh directly.
        let mut vns = vec![Vec3::ZERO; nv];
        let mut fns = vec![Vec3::ZERO; nf];

        #[cfg(feature = "parallel")]
        fns.par_iter_mut().enumerate().for_each(|(i, n)| {
            let ih = f2h[i];
            let p2 = pos[head[ih]];
            let p1 = pos[tail[ih]];
            let p0 = pos[tail[prev[ih]]];
            let x = p2 - p1;
            let t = (p1 - p0) * -1.;
            *n = x.cross(t).normalize();
        });

        #[cfg(not(feature = "parallel"))]
        for i in 0..nf {
            let ih = f2h[i];
            let p2 = pos[head[ih]];
            let p1 = pos[tail[ih]];
            let p0 = pos[tail[prev[ih]]];
            let x = p2 - p1;
            let t = (p1 - p0) * -1.;
            fns[i] = x.cross(t).normalize();
        }

        for i in 0..nf {
            for j in 0..3 {
                let i_curr = idx[i][j];
                let v_prev = pos[idx[i][(j + 2) % 3]];
                let v_curr = pos[i_curr];
                let v_next = pos[idx[i][(j + 1) % 3]];
                let e_curr = (v_next - v_curr).normalize();
                let e_prev = (v_curr - v_prev).normalize();
                if e_curr.is_nan() || e_prev.is_nan() {
                    continue;
                }
                let dot = -e_prev.dot(e_curr);
                let phi = if dot >= 1. {
                    0.
                } else if dot <= -1. {
                    PI as Real
                } else {
                    dot.acos()
                };
                vns[i_curr] += fns[i] * phi;
            }
        }

        #[cfg(feature = "parallel")]
        vns.par_iter_mut().for_each(|n| *n = n.normalize_or_zero());

        #[cfg(not(feature = "parallel"))]
        for n in &mut vns {
            *n = n.normalize_or_zero();
        }

        Ok(Hmesh {
            nv,
            nf,
            nh,
            twin,
            head,
            tail,
            vns,
            fns,
        })
    }
}

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

    /// kernel#102: an edge shared by three faces must refuse.
    ///
    /// Three triangles fanned around the shared edge (0, 1), so that
    /// edge has three incident faces. The half-edge table has two
    /// face slots per edge, so this is unrepresentable and must be
    /// rejected -- it used to overrun the allocation and abort.
    #[test]
    fn an_edge_with_three_faces_is_refused() {
        let pos = vec![
            Vec3::new(0., 0., 0.),
            Vec3::new(1., 0., 0.),
            Vec3::new(0., 1., 0.),
            Vec3::new(0., -1., 0.),
            Vec3::new(0., 0., 1.),
        ];
        let idx = vec![
            Vec3u::new(0, 1, 2),
            Vec3u::new(0, 1, 3),
            Vec3u::new(0, 1, 4),
        ];
        let (mut e2v, mut e2f, mut f2e) = (vec![], vec![], vec![]);
        let r = edge_topology(&pos, &idx, &mut e2v, &mut e2f, &mut f2e);
        let err = r.expect_err("three faces on one edge must refuse");
        assert!(err.contains("incident faces"), "unexpected: {err}");
    }

    /// The two-face case still builds, so the guard did not
    /// over-reject ordinary interior edges.
    #[test]
    fn an_edge_with_two_faces_still_builds() {
        let pos = vec![
            Vec3::new(0., 0., 0.),
            Vec3::new(1., 0., 0.),
            Vec3::new(0., 1., 0.),
            Vec3::new(0., -1., 0.),
        ];
        let idx = vec![Vec3u::new(0, 1, 2), Vec3u::new(0, 1, 3)];
        let (mut e2v, mut e2f, mut f2e) = (vec![], vec![], vec![]);
        assert!(edge_topology(&pos, &idx, &mut e2v, &mut e2f, &mut f2e).is_ok());
    }
}