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.

pub mod bounds;
pub mod collider;
pub mod hmesh;

use super::hmesh::Hmesh;
use crate::csg::collider::{morton_code, MortonCollider, PlanarGrid, K_NO_CODE};
use crate::csg::{next_of, Half, Real, Tref, Vec3, Vec3u, K_PRECISION};
use bounds::BBox;
#[cfg(feature = "parallel")]
use rayon::prelude::*;
use std::cmp::Ordering;
use std::collections::HashMap;

#[derive(Clone, Debug)]
pub struct Manifold {
    pub ps: Vec<Vec3>,            // positions
    pub hs: Vec<Half>,            // halfedges
    pub nv: usize,                // number of vertices
    pub nf: usize,                // number of faces
    pub nh: usize,                // number of halfedges
    pub eps: Real,                // epsilon
    pub face_normals: Vec<Vec3>,  //
    pub vert_normals: Vec<Vec3>,  //
    pub collider: MortonCollider, //
    /// A second broad phase over x/y alone, O(n) to build.
    ///
    /// `winding03`'s point-in-polygon test reads only x and y (see
    /// `BPos::overlaps_node`), and is run exactly once per operand --
    /// there is no second call to amortize a tree's O(n log n) build
    /// against. A uniform grid is O(n) to build and, for a subdivided
    /// mesh where faces are close to uniform in size, resolves the
    /// query about as well as a tree does.
    pub planar_grid: PlanarGrid,
    pub coplanar: Vec<i32>, // indices of coplanar faces
    /// The caller's triangle each face came from.
    ///
    /// `new` welds vertices, drops triangles the weld collapsed, and
    /// `new_impl` Morton-sorts the faces, so face `f` here is generally NOT
    /// input triangle `f`. `Tref::fid` indexes faces of THIS structure;
    /// composing it with `face_src` is how a result triangle is traced back
    /// to the triangle -- and therefore the attribute values -- the caller
    /// supplied (#116).
    pub face_src: Vec<usize>,
}

impl Manifold {
    pub fn new(pos: &[f64], idx: &[usize]) -> Result<Self, String> {
        if pos.len() % 3 != 0 {
            return Err("pos must be a multiple of 3".into());
        }
        if idx.len() % 3 != 0 {
            return Err("idx must be a multiple of 3".into());
        }

        // dedup vertices
        let mut hash = HashMap::with_capacity(pos.len() / 3);
        let mut weld = Vec::with_capacity(pos.len() / 3);
        // One entry per VERTEX, not per coordinate: `i` below counts
        // `pos.chunks(3)`. Upstream allocated `pos.len()`, three times what
        // the loop can address.
        let mut rmap = vec![0; pos.len() / 3];

        for (i, p) in pos.chunks(3).enumerate() {
            let v = Vec3::new(p[0] as Real, p[1] as Real, p[2] as Real);
            let k = (v.x.to_bits(), v.y.to_bits(), v.z.to_bits());
            if let Some(&w) = hash.get(&k) {
                rmap[i] = w;
            } else {
                let n = weld.len();
                weld.push(v);
                hash.insert(k, n);
                rmap[i] = n;
            }
        }

        // remove collapsed triangles, remembering which input triangle each
        // survivor was
        let (idx, kept): (Vec<_>, Vec<_>) = idx
            .chunks(3)
            .map(|i| Vec3u::new(rmap[i[0]], rmap[i[1]], rmap[i[2]]))
            .enumerate()
            .filter(|(_, is)| is.x != is.y && is.y != is.z && is.z != is.x)
            .map(|(t, is)| (is, t))
            .unzip();

        let mut mfd = Self::new_impl(weld, idx)?;
        // `new_impl` recorded positions in ITS input; lift them to ours.
        for src in &mut mfd.face_src {
            *src = kept[*src];
        }
        Ok(mfd)
    }

    pub fn new_impl(ps: Vec<Vec3>, idx: Vec<Vec3u>) -> Result<Self, String> {
        let bb = BBox::new(None, &ps);
        let (mut f_bb, mut f_mt) = compute_face_morton(&ps, &idx, &bb);
        let (hm, face_src) = sort_faces(&ps, &idx, &mut f_bb, &mut f_mt)?;
        let hs = (0..hm.nh)
            .map(|i| Half::new(hm.tail[i], hm.head[i], hm.twin[i]))
            .collect::<Vec<_>>();

        let mut e = K_PRECISION * bb.scale();
        e = if e.is_finite() { e } else { -1. };
        let eps = e;
        let collider = MortonCollider::new(&f_bb, &f_mt);
        let planar_grid = PlanarGrid::new(&f_bb, &bb);
        let coplanar = compute_coplanar_idx(&ps, &hm.fns, &hs, eps);

        let mfd = Manifold {
            nv: hm.nv,
            nf: hm.nf,
            nh: hm.nh,
            ps,
            hs,
            vert_normals: hm.vns,
            face_normals: hm.fns,
            eps,
            collider,
            planar_grid,
            coplanar,
            face_src,
        };

        if !mfd.is_manifold() {
            return Err("The input mesh is not manifold".into());
        }
        Ok(mfd)
    }

    pub fn is_manifold(&self) -> bool {
        self.hs.iter().enumerate().all(|(i, h)| {
            if h.tail().is_none() || h.head().is_none() {
                return true;
            }
            match h.pair() {
                None => false,
                Some(pair) => {
                    let mut good = true;
                    good &= self.hs[pair].pair() == Some(i);
                    good &= h.tail != h.head;
                    good &= h.tail == self.hs[pair].head;
                    good &= h.head == self.hs[pair].tail;
                    good
                }
            }
        })
    }
}

fn compute_face_morton(pos: &[Vec3], idx: &[Vec3u], bb: &BBox) -> (Vec<BBox>, Vec<u32>) {
    let n = idx.len();
    let mut bbs = vec![BBox::default(); n];
    let mut mts = vec![0; n];

    #[cfg(feature = "parallel")]
    {
        bbs.par_iter_mut()
            .zip(mts.par_iter_mut())
            .zip(idx.par_iter())
            .for_each(|((bb_, mt_), f)| {
                let p0 = pos[f.x];
                let p1 = pos[f.y];
                let p2 = pos[f.z];
                bb_.union(&p0);
                bb_.union(&p1);
                bb_.union(&p2);
                *mt_ = morton_code(&((p0 + p1 + p2) / 3.), bb);
            });
    }

    #[cfg(not(feature = "parallel"))]
    {
        for (i, f) in idx.iter().enumerate() {
            let p0 = pos[f.x];
            let p1 = pos[f.y];
            let p2 = pos[f.z];
            bbs[i].union(&p0);
            bbs[i].union(&p1);
            bbs[i].union(&p2);
            mts[i] = morton_code(&((p0 + p1 + p2) / 3.), bb);
        }
    }

    (bbs, mts)
}

fn sort_faces(
    pos: &[Vec3],
    idx: &[Vec3u],
    face_bboxes: &mut Vec<BBox>,
    face_morton: &mut Vec<u32>,
) -> Result<(Hmesh, Vec<usize>), String> {
    let mut map = (0..face_morton.len()).collect::<Vec<_>>();
    // Morton codes are u32 keys and the permutation is rebuilt from scratch,
    // so equal-key order is not observable: the unstable sort is free here.
    map.sort_unstable_by_key(|&i| face_morton[i]);
    *face_bboxes = map.iter().map(|&i| face_bboxes[i]).collect::<Vec<_>>();
    *face_morton = map.iter().map(|&i| face_morton[i]).collect::<Vec<_>>();

    // `Hmesh::new` keeps face order, so sorted face `f` is `idx[map[f]]`:
    // `map` IS the face-to-input record, and costs nothing extra to keep.
    let hm = Hmesh::new(pos, &map.iter().map(|&i| idx[i]).collect::<Vec<_>>())?;
    Ok((hm, map))
}

fn compute_coplanar_idx(ps: &[Vec3], ns: &[Vec3], hs: &[Half], tol: Real) -> Vec<i32> {
    let nt = hs.len() / 3;
    let mut priority = vec![];
    let mut res = vec![-1; nt];

    for t in 0..nt {
        let i = t * 3;
        let area = if hs[i].tail().is_none() {
            0.
        } else {
            let p0 = ps[hs[i].tail];
            let p1 = ps[hs[i].head];
            let p2 = ps[hs[i + 1].head];
            (p1 - p0).cross(p2 - p0).length_squared()
        };
        priority.push((area, t));
    }

    priority.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(Ordering::Equal));

    let mut interior = vec![];
    for (_, t) in priority.iter() {
        if res[*t] != -1 {
            continue;
        }
        res[*t] = *t as i32;

        let i = t * 3;
        let p = ps[hs[i].tail];
        let n = ns[*t];

        interior.clear();
        interior.extend_from_slice(&[i, i + 1, i + 2]);

        while let Some(hi) = interior.pop() {
            let h1 = next_of(hs[hi].pair);
            let t1 = h1 / 3;

            if res[t1] != -1 {
                continue;
            }

            if (ps[hs[h1].head] - p).dot(n).abs() < tol {
                res[t1] = *t as i32;
                if interior.last().copied() == Some(hs[h1].pair) {
                    interior.pop();
                } else {
                    interior.push(h1);
                }
                interior.push(next_of(h1));
            }
        }
    }
    res
}

pub fn cleanup_unused_verts(ps: &mut Vec<Vec3>, hs: &mut Vec<Half>, rs: &mut Vec<Tref>) {
    let bb = BBox::new(None, ps);
    let mt = ps.iter().map(|p| morton_code(p, &bb)).collect::<Vec<_>>();

    let mut new2old = (0..ps.len()).collect::<Vec<_>>();
    let mut old2new = vec![0; ps.len()];
    new2old.sort_by_key(|&i| mt[i]);
    for (new, &old) in new2old.iter().enumerate() {
        old2new[old] = new;
    }

    // reindex verts
    for h in hs.iter_mut() {
        if h.pair().is_none() {
            continue;
        }
        h.tail = old2new[h.tail];
        h.head = old2new[h.head];
    }

    // truncate pos container
    let nv = new2old
        .iter()
        // `K_NO_CODE` is `u32::MAX`, so `>=` could only ever mean `==`.
        // Written explicitly: same behaviour, and it no longer reads as if a
        // greater-than case existed.
        .position(|&v| mt[v] == K_NO_CODE)
        .unwrap_or(new2old.len());

    new2old.truncate(nv);

    *ps = new2old.iter().map(|&i| ps[i]).collect();
    // A removed face is three default halfedges, so faces are kept or
    // dropped whole; its ref goes with it, keeping `rs` face-parallel.
    if rs.len() * 3 == hs.len() {
        *rs = hs
            .chunks(3)
            .zip(rs.iter())
            .filter(|(face, _)| face[0].pair().is_some())
            .map(|(_, r)| *r)
            .collect();
    }
    *hs = hs.iter().filter(|h| h.pair().is_some()).cloned().collect();
}

/// Whether a half-edge array describes a closed two-manifold.
///
/// Extracted verbatim from `Manifold::is_manifold` so the result of a
/// boolean can be validated without building a `Manifold` around it.
/// Same predicate, same acceptance of unset tail/head as vacuously fine.
pub(crate) fn halfedges_are_two_manifold(hs: &[Half]) -> bool {
    hs.iter().enumerate().all(|(i, h)| {
        if h.tail().is_none() || h.head().is_none() {
            return true;
        }
        match h.pair() {
            None => false,
            Some(pair) => {
                let mut good = true;
                good &= hs[pair].pair() == Some(i);
                good &= h.tail != h.head;
                good &= h.tail == hs[pair].head;
                good &= h.head == hs[pair].tail;
                good
            }
        }
    })
}