manifold-rust 0.9.2

// Phase 2: Core Types — ported from include/manifold/common.h, include/manifold/polygon.h,
// include/manifold/manifold.h (MeshGLP/Error), src/shared.h (Halfedge, TriRef, Barycentric, TmpEdge)

use std::collections::BTreeMap;
use crate::linalg::{Vec2, Vec3, Vec4, Mat3, Mat3x4};
use crate::math;

// ---------------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------------

pub const K_PI: f64 = std::f64::consts::PI;
pub const K_TWO_PI: f64 = std::f64::consts::TAU;
pub const K_HALF_PI: f64 = std::f64::consts::FRAC_PI_2;
/// Precision used for epsilon calculations relative to bounding-box scale.
pub const K_PRECISION: f64 = 1e-12;

pub const DEFAULT_SEGMENTS: i32 = 0;
pub const DEFAULT_ANGLE: f64 = 10.0;
pub const DEFAULT_LENGTH: f64 = 1.0;

// ---------------------------------------------------------------------------
// Scalar utilities
// ---------------------------------------------------------------------------

#[inline]
pub fn radians(a: f64) -> f64 {
    a * K_PI / 180.0
}

#[inline]
pub fn degrees(a: f64) -> f64 {
    a * 180.0 / K_PI
}

/// Smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1.
#[inline]
pub fn smoothstep(edge0: f64, edge1: f64, a: f64) -> f64 {
    let x = ((a - edge0) / (edge1 - edge0)).clamp(0.0, 1.0);
    x * x * (3.0 - 2.0 * x)
}

/// Sine function where multiples of 90 degrees come out exact.
///
/// Matches C++ `sind` (common.h), which reduces the argument with
/// `std::remquo(x, 90.0, &quo)` — round-to-nearest (ties to even), remainder
/// in [-45, 45]. A floor-based reduction (remainder in [0, 90)) is
/// mathematically equal but differs by ~1 ULP for reduced arguments in
/// (45, 90), which breaks bit-exactness with the C++ reference (e.g. cylinder
/// circle vertices).
pub fn sind(x: f64) -> f64 {
    if !x.is_finite() {
        return f64::NAN;
    }
    if x < 0.0 {
        return -sind(-x);
    }
    // Reconstruct std::remquo(x, 90.0, &quo): quo = nearest integer to the
    // exact x/90 (ties to even), remainder computed exactly. Round the
    // computed quotient, then fix up the rare off-by-one where the rounded
    // double quotient disagrees with the exact one.
    let mut quo = (x / 90.0).round_ties_even() as i64;
    let mut r = x - quo as f64 * 90.0;
    if r > 45.0 {
        quo += 1;
        r -= 90.0;
    } else if r < -45.0 {
        quo -= 1;
        r += 90.0;
    } else if r == 45.0 && quo % 2 != 0 {
        // Exact tie: remquo rounds the quotient to even.
        quo += 1;
        r = -45.0;
    } else if r == -45.0 && quo % 2 != 0 {
        quo -= 1;
        r = 45.0;
    }
    match ((quo % 4) + 4) % 4 {
        0 => math::sin(radians(r)),
        1 => math::cos(radians(r)),
        2 => -math::sin(radians(r)),
        3 => -math::cos(radians(r)),
        _ => 0.0,
    }
}

/// Cosine function where multiples of 90 degrees come out exact.
#[inline]
pub fn cosd(x: f64) -> f64 {
    sind(x + 90.0)
}

// ---------------------------------------------------------------------------
// Polygon types
// ---------------------------------------------------------------------------

/// Single polygon contour, wound CCW. First and last point are implicitly connected.
pub type SimplePolygon = Vec<Vec2>;

/// Set of polygons with holes (arbitrary nesting).
pub type Polygons = Vec<SimplePolygon>;

/// Polygon vertex with index.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct PolyVert {
    pub pos: Vec2,
    pub idx: i32,
}

/// Single indexed polygon contour, wound CCW.
pub type SimplePolygonIdx = Vec<PolyVert>;

/// Set of indexed polygons with holes.
pub type PolygonsIdx = Vec<SimplePolygonIdx>;

// ---------------------------------------------------------------------------
// Box (3D axis-aligned bounding box)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Box {
    pub min: Vec3,
    pub max: Vec3,
}

impl Default for Box {
    fn default() -> Self {
        Box {
            min: Vec3::splat(f64::INFINITY),
            max: Vec3::splat(f64::NEG_INFINITY),
        }
    }
}

impl Box {
    /// Default is an infinite box containing all space.
    pub fn new() -> Self {
        Self::default()
    }

    /// Box containing the two given points.
    pub fn from_points(p1: Vec3, p2: Vec3) -> Self {
        Box {
            min: Vec3::new(p1.x.min(p2.x), p1.y.min(p2.y), p1.z.min(p2.z)),
            max: Vec3::new(p1.x.max(p2.x), p1.y.max(p2.y), p1.z.max(p2.z)),
        }
    }

    /// A box containing a single point.
    pub fn from_point(p: Vec3) -> Self {
        Box { min: p, max: p }
    }

    /// True when the box has no volume (min > max on any axis).
    pub fn is_empty(&self) -> bool {
        self.min.x > self.max.x || self.min.y > self.max.y || self.min.z > self.max.z
    }

    pub fn size(&self) -> Vec3 {
        self.max - self.min
    }

    pub fn center(&self) -> Vec3 {
        (self.max + self.min) * 0.5
    }

    /// Absolute-largest coordinate value.
    pub fn scale(&self) -> f64 {
        let abs_min = Vec3::new(self.min.x.abs(), self.min.y.abs(), self.min.z.abs());
        let abs_max = Vec3::new(self.max.x.abs(), self.max.y.abs(), self.max.z.abs());
        let m = Vec3::new(
            abs_min.x.max(abs_max.x),
            abs_min.y.max(abs_max.y),
            abs_min.z.max(abs_max.z),
        );
        m.x.max(m.y).max(m.z)
    }

    pub fn contains_point(&self, p: Vec3) -> bool {
        p.x >= self.min.x && p.x <= self.max.x
            && p.y >= self.min.y && p.y <= self.max.y
            && p.z >= self.min.z && p.z <= self.max.z
    }

    pub fn contains_box(&self, other: &Box) -> bool {
        other.min.x >= self.min.x && other.max.x <= self.max.x
            && other.min.y >= self.min.y && other.max.y <= self.max.y
            && other.min.z >= self.min.z && other.max.z <= self.max.z
    }

    /// Expand in-place to include the given point.
    pub fn union_point(&mut self, p: Vec3) {
        self.min.x = self.min.x.min(p.x);
        self.min.y = self.min.y.min(p.y);
        self.min.z = self.min.z.min(p.z);
        self.max.x = self.max.x.max(p.x);
        self.max.y = self.max.y.max(p.y);
        self.max.z = self.max.z.max(p.z);
    }

    /// Return the union of this box with another.
    pub fn union_box(&self, other: &Box) -> Box {
        Box {
            min: Vec3::new(
                self.min.x.min(other.min.x),
                self.min.y.min(other.min.y),
                self.min.z.min(other.min.z),
            ),
            max: Vec3::new(
                self.max.x.max(other.max.x),
                self.max.y.max(other.max.y),
                self.max.z.max(other.max.z),
            ),
        }
    }

    /// Transform by axis-aligned affine transform (Mat3x4 * vec4(pt, 1)).
    pub fn transform(&self, t: &Mat3x4) -> Box {
        use crate::linalg::Vec4 as V4;
        let min_t = *t * V4::new(self.min.x, self.min.y, self.min.z, 1.0);
        let max_t = *t * V4::new(self.max.x, self.max.y, self.max.z, 1.0);
        Box {
            min: Vec3::new(min_t.x.min(max_t.x), min_t.y.min(max_t.y), min_t.z.min(max_t.z)),
            max: Vec3::new(min_t.x.max(max_t.x), min_t.y.max(max_t.y), min_t.z.max(max_t.z)),
        }
    }

    pub fn does_overlap_box(&self, other: &Box) -> bool {
        self.min.x <= other.max.x && self.min.y <= other.max.y && self.min.z <= other.max.z
            && self.max.x >= other.min.x && self.max.y >= other.min.y && self.max.z >= other.min.z
    }

    /// Does the given point project within the XY extent (including equality)?
    pub fn does_overlap_point_xy(&self, p: Vec3) -> bool {
        p.x >= self.min.x && p.x <= self.max.x && p.y >= self.min.y && p.y <= self.max.y
    }

    pub fn is_finite(&self) -> bool {
        self.min.x.is_finite() && self.min.y.is_finite() && self.min.z.is_finite()
            && self.max.x.is_finite() && self.max.y.is_finite() && self.max.z.is_finite()
    }
}

impl std::ops::Add<Vec3> for Box {
    type Output = Box;
    fn add(self, shift: Vec3) -> Box {
        Box { min: self.min + shift, max: self.max + shift }
    }
}
impl std::ops::AddAssign<Vec3> for Box {
    fn add_assign(&mut self, shift: Vec3) {
        self.min = self.min + shift;
        self.max = self.max + shift;
    }
}
impl std::ops::Mul<Vec3> for Box {
    type Output = Box;
    fn mul(self, scale: Vec3) -> Box {
        Box { min: self.min * scale, max: self.max * scale }
    }
}
impl std::ops::MulAssign<Vec3> for Box {
    fn mul_assign(&mut self, scale: Vec3) {
        self.min = self.min * scale;
        self.max = self.max * scale;
    }
}

// ---------------------------------------------------------------------------
// Rect (2D axis-aligned bounding box)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Rect {
    pub min: Vec2,
    pub max: Vec2,
}

impl Default for Rect {
    fn default() -> Self {
        Rect {
            min: Vec2::splat(f64::INFINITY),
            max: Vec2::splat(f64::NEG_INFINITY),
        }
    }
}

impl Rect {
    pub fn new() -> Self {
        Self::default()
    }

    pub fn from_points(a: Vec2, b: Vec2) -> Self {
        Rect {
            min: Vec2::new(a.x.min(b.x), a.y.min(b.y)),
            max: Vec2::new(a.x.max(b.x), a.y.max(b.y)),
        }
    }

    pub fn size(&self) -> Vec2 {
        self.max - self.min
    }

    pub fn area(&self) -> f64 {
        let sz = self.size();
        sz.x * sz.y
    }

    pub fn scale(&self) -> f64 {
        let abs_min = Vec2::new(self.min.x.abs(), self.min.y.abs());
        let abs_max = Vec2::new(self.max.x.abs(), self.max.y.abs());
        let m = Vec2::new(abs_min.x.max(abs_max.x), abs_min.y.max(abs_max.y));
        m.x.max(m.y)
    }

    pub fn center(&self) -> Vec2 {
        (self.max + self.min) * 0.5
    }

    pub fn contains_point(&self, p: Vec2) -> bool {
        p.x >= self.min.x && p.x <= self.max.x && p.y >= self.min.y && p.y <= self.max.y
    }

    pub fn contains_rect(&self, other: &Rect) -> bool {
        other.min.x >= self.min.x && other.max.x <= self.max.x
            && other.min.y >= self.min.y && other.max.y <= self.max.y
    }

    pub fn does_overlap(&self, other: &Rect) -> bool {
        self.min.x <= other.max.x && self.min.y <= other.max.y
            && self.max.x >= other.min.x && self.max.y >= other.min.y
    }

    pub fn is_empty(&self) -> bool {
        self.max.y <= self.min.y || self.max.x <= self.min.x
    }

    pub fn is_finite(&self) -> bool {
        self.min.x.is_finite() && self.min.y.is_finite()
            && self.max.x.is_finite() && self.max.y.is_finite()
    }

    pub fn union_point(&mut self, p: Vec2) {
        self.min.x = self.min.x.min(p.x);
        self.min.y = self.min.y.min(p.y);
        self.max.x = self.max.x.max(p.x);
        self.max.y = self.max.y.max(p.y);
    }

    pub fn union_rect(&self, other: &Rect) -> Rect {
        Rect {
            min: Vec2::new(self.min.x.min(other.min.x), self.min.y.min(other.min.y)),
            max: Vec2::new(self.max.x.max(other.max.x), self.max.y.max(other.max.y)),
        }
    }
}

impl std::ops::Add<Vec2> for Rect {
    type Output = Rect;
    fn add(self, shift: Vec2) -> Rect {
        Rect { min: self.min + shift, max: self.max + shift }
    }
}
impl std::ops::AddAssign<Vec2> for Rect {
    fn add_assign(&mut self, shift: Vec2) {
        self.min = self.min + shift;
        self.max = self.max + shift;
    }
}
impl std::ops::Mul<Vec2> for Rect {
    type Output = Rect;
    fn mul(self, scale: Vec2) -> Rect {
        Rect { min: self.min * scale, max: self.max * scale }
    }
}
impl std::ops::MulAssign<Vec2> for Rect {
    fn mul_assign(&mut self, scale: Vec2) {
        self.min = self.min * scale;
        self.max = self.max * scale;
    }
}

// ---------------------------------------------------------------------------
// OpType
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum OpType {
    Add,
    Subtract,
    Intersect,
}

// ---------------------------------------------------------------------------
// Error
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Error {
    NoError,
    NonFiniteVertex,
    NotManifold,
    VertexOutOfBounds,
    PropertiesWrongLength,
    MissingPositionProperties,
    MergeVectorsDifferentLengths,
    MergeIndexOutOfBounds,
    TransformWrongLength,
    RunIndexWrongLength,
    FaceIdWrongLength,
    InvalidConstruction,
    ResultTooLarge,
    InvalidTangents,
}

impl Error {
    pub fn to_str(self) -> &'static str {
        match self {
            Error::NoError => "No Error",
            Error::NonFiniteVertex => "Non-Finite Vertex",
            Error::NotManifold => "Not Manifold",
            Error::VertexOutOfBounds => "Vertex Out of Bounds",
            Error::PropertiesWrongLength => "Properties Wrong Length",
            Error::MissingPositionProperties => "Missing Position Properties",
            Error::MergeVectorsDifferentLengths => "Merge Vectors Different Lengths",
            Error::MergeIndexOutOfBounds => "Merge Index Out of Bounds",
            Error::TransformWrongLength => "Transform Wrong Length",
            Error::RunIndexWrongLength => "Run Index Wrong Length",
            Error::FaceIdWrongLength => "Face ID Wrong Length",
            Error::InvalidConstruction => "Invalid Construction",
            Error::ResultTooLarge => "Result Too Large",
            Error::InvalidTangents => "Invalid Tangents",
        }
    }
}

impl std::fmt::Display for Error {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.to_str())
    }
}

// ---------------------------------------------------------------------------
// Quality (static global for circle quantization)
// ---------------------------------------------------------------------------

use std::sync::OnceLock;
use std::sync::Mutex;

struct QualityState {
    min_circular_angle: f64,
    min_circular_edge_length: f64,
    circular_segments: i32,
}

static QUALITY_STATE: OnceLock<Mutex<QualityState>> = OnceLock::new();

fn quality_state() -> &'static Mutex<QualityState> {
    QUALITY_STATE.get_or_init(|| {
        Mutex::new(QualityState {
            min_circular_angle: DEFAULT_ANGLE,
            min_circular_edge_length: DEFAULT_LENGTH,
            circular_segments: DEFAULT_SEGMENTS,
        })
    })
}

pub struct Quality;

impl Quality {
    pub fn set_min_circular_angle(angle: f64) {
        quality_state().lock().unwrap().min_circular_angle = angle;
    }

    pub fn set_min_circular_edge_length(length: f64) {
        quality_state().lock().unwrap().min_circular_edge_length = length;
    }

    pub fn set_circular_segments(n: i32) {
        quality_state().lock().unwrap().circular_segments = n;
    }

    pub fn get_circular_segments(radius: f64) -> i32 {
        let q = quality_state().lock().unwrap();
        if q.circular_segments > 0 {
            return q.circular_segments;
        }
        // Match C++ exactly: int truncation (not ceil), fmin (not fmax), round down to multiple of 4
        let n_seg_a = (360.0 / q.min_circular_angle) as i32;
        let n_seg_l = (2.0 * radius.abs() * K_PI / q.min_circular_edge_length) as i32;
        let mut n_seg = n_seg_a.min(n_seg_l) + 3;
        n_seg -= n_seg % 4;
        n_seg.max(4)
    }

    pub fn reset_to_defaults() {
        let mut q = quality_state().lock().unwrap();
        q.min_circular_angle = DEFAULT_ANGLE;
        q.min_circular_edge_length = DEFAULT_LENGTH;
        q.circular_segments = DEFAULT_SEGMENTS;
    }
}

// ---------------------------------------------------------------------------
// ExecutionParams
// ---------------------------------------------------------------------------

#[derive(Clone, Debug)]
pub struct ExecutionParams {
    pub intermediate_checks: bool,
    pub self_intersection_checks: bool,
    pub process_overlaps: bool,
    pub suppress_errors: bool,
    pub cleanup_triangles: bool,
    pub verbose: i32,
}

impl Default for ExecutionParams {
    fn default() -> Self {
        ExecutionParams {
            intermediate_checks: false,
            self_intersection_checks: false,
            process_overlaps: true,
            suppress_errors: false,
            cleanup_triangles: true,
            verbose: 0,
        }
    }
}

// ---------------------------------------------------------------------------
// Smoothness
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Smoothness {
    /// The halfedge index = 3 * tri + i
    pub halfedge: usize,
    /// 0 = sharp, 1 = smooth
    pub smoothness: f64,
}

// ---------------------------------------------------------------------------
// MeshGLP / MeshGL / MeshGL64
// ---------------------------------------------------------------------------

/// GL-style mesh representation. Generic over precision (f32/f64) and index type (u32/u64).
#[derive(Clone, Debug, Default)]
pub struct MeshGLP<P: Copy + Default, I: Copy + Default = u32> {
    /// Number of properties per vertex, always >= 3.
    pub num_prop: I,
    /// Flat interleaved vertex properties: [x, y, z, ...] × num_verts.
    pub vert_properties: Vec<P>,
    /// Triangle vertex indices, 3 per triangle (CCW from outside).
    pub tri_verts: Vec<I>,
    /// Optional: merge-from vertex indices.
    pub merge_from_vert: Vec<I>,
    /// Optional: merge-to vertex indices.
    pub merge_to_vert: Vec<I>,
    /// Optional: run start indices into triVerts.
    pub run_index: Vec<I>,
    /// Optional: original mesh ID per run.
    pub run_original_id: Vec<u32>,
    /// Optional: 3×4 column-major transform per run (12 elements each).
    pub run_transform: Vec<P>,
    /// Optional: source face ID per triangle.
    pub face_id: Vec<I>,
    /// Optional: halfedge tangent vectors (4 per halfedge).
    pub halfedge_tangent: Vec<P>,
    /// Optional: per-run flags; 1 = backside (normals need flipping).
    pub run_flags: Vec<u8>,
    /// Tolerance for mesh simplification.
    pub tolerance: P,
}

impl MeshGLP<f32, u32> {
    pub fn num_vert(&self) -> usize {
        if self.num_prop == 0 { 0 } else { self.vert_properties.len() / self.num_prop as usize }
    }

    pub fn num_tri(&self) -> usize {
        self.tri_verts.len() / 3
    }

    pub fn get_vert_pos(&self, v: usize) -> [f32; 3] {
        let offset = v * self.num_prop as usize;
        [self.vert_properties[offset], self.vert_properties[offset + 1], self.vert_properties[offset + 2]]
    }

    pub fn get_tri_verts(&self, t: usize) -> [u32; 3] {
        let offset = 3 * t;
        [self.tri_verts[offset], self.tri_verts[offset + 1], self.tri_verts[offset + 2]]
    }

    pub fn get_tangent(&self, h: usize) -> [f32; 4] {
        let offset = 4 * h;
        [
            self.halfedge_tangent[offset],
            self.halfedge_tangent[offset + 1],
            self.halfedge_tangent[offset + 2],
            self.halfedge_tangent[offset + 3],
        ]
    }

    /// Merges coincident vertices based on position within tolerance.
    /// Uses BVH collision detection to find open edges, then groups
    /// coincident vertices via union-find. Returns true if new merges
    /// were found, false if the mesh was already fully merged.
    pub fn merge(&mut self) -> bool {
        use crate::collider::Collider;
        use crate::disjoint_sets::DisjointSets;
        use crate::sort::morton_code;
        use std::collections::BTreeSet;

        let num_vert = self.num_vert();
        let num_tri = self.num_tri();

        // Build initial merge map from existing merge vectors
        let mut merge_map: Vec<usize> = (0..num_vert).collect();
        for i in 0..self.merge_from_vert.len() {
            merge_map[self.merge_from_vert[i] as usize] = self.merge_to_vert[i] as usize;
        }

        // Find open (non-manifold) edges
        let next = [1usize, 2, 0];
        let mut open_edges: BTreeSet<(usize, usize)> = BTreeSet::new();
        for tri in 0..num_tri {
            for i in 0..3 {
                let a = merge_map[self.tri_verts[3 * tri + next[i]] as usize];
                let b = merge_map[self.tri_verts[3 * tri + i] as usize];
                let edge = (a, b);
                // Look for the reverse edge
                let rev = (b, a);
                if open_edges.contains(&rev) {
                    open_edges.remove(&rev);
                } else {
                    open_edges.insert(edge);
                }
            }
        }

        if open_edges.is_empty() {
            return false;
        }

        // Collect unique open vertices — only the START vertex of each open
        // halfedge, matching C++ which stores (start,end) and takes edge.first=start.
        // Our BTreeSet stores (end,start) so we take edge.1 (= b = start vertex).
        let open_verts: Vec<usize> = {
            let mut vset = std::collections::BTreeSet::new();
            for (_a, b) in &open_edges {
                vset.insert(*b);
            }
            vset.into_iter().collect()
        };
        let num_open = open_verts.len();

        // Compute bounding box
        let mut bbox = Box::default();
        for v in 0..num_vert {
            let pos = self.get_vert_pos(v);
            let p = Vec3::new(pos[0] as f64, pos[1] as f64, pos[2] as f64);
            bbox.union_point(p);
        }

        let tolerance = f64::max(
            self.tolerance as f64,
            f32::EPSILON as f64 * bbox.scale(),
        );

        // Build BVH boxes and morton codes for open vertices
        let mut vert_box: Vec<Box> = Vec::with_capacity(num_open);
        let mut vert_morton: Vec<u32> = Vec::with_capacity(num_open);
        for &v in &open_verts {
            let pos = self.get_vert_pos(v);
            let center = Vec3::new(pos[0] as f64, pos[1] as f64, pos[2] as f64);
            let half_tol = tolerance / 2.0;
            let bx = Box::from_points(
                center - Vec3::new(half_tol, half_tol, half_tol),
                center + Vec3::new(half_tol, half_tol, half_tol),
            );
            vert_box.push(bx);
            vert_morton.push(morton_code(center, &bbox));
        }

        // Sort by morton code
        let mut order: Vec<usize> = (0..num_open).collect();
        order.sort_by_key(|&i| vert_morton[i]);

        let sorted_box: Vec<Box> = order.iter().map(|&i| vert_box[i]).collect();
        let sorted_morton: Vec<u32> = order.iter().map(|&i| vert_morton[i]).collect();
        let sorted_verts: Vec<usize> = order.iter().map(|&i| open_verts[i]).collect();

        // Build collider and find coincident vertex pairs
        let collider = Collider::new(sorted_box.clone(), sorted_morton);
        let uf = DisjointSets::new(num_vert as u32);

        collider.collisions_with_boxes(&sorted_box, false, |a, b| {
            uf.unite(sorted_verts[a] as u32, sorted_verts[b] as u32);
        });

        // Also merge from existing merge vectors
        for i in 0..self.merge_from_vert.len() {
            uf.unite(self.merge_from_vert[i], self.merge_to_vert[i]);
        }

        // Rebuild merge vectors
        self.merge_from_vert.clear();
        self.merge_to_vert.clear();
        for v in 0..num_vert {
            let merge_to = uf.find(v as u32) as usize;
            if merge_to != v {
                self.merge_from_vert.push(v as u32);
                self.merge_to_vert.push(merge_to as u32);
            }
        }

        true
    }

    /// True if triangle run `run` is on the backside (e.g. from a subtraction).
    /// run_flags is a bitmask (#1718): bit 0 = backside. Informational only —
    /// the framework already orients stored normals on the standard flow.
    pub fn backside(&self, run: usize) -> bool {
        run < self.run_flags.len() && (self.run_flags[run] & 1) != 0
    }

    /// True if the first three extra-property channels (slots 3, 4, 5) of run
    /// `run` carry world-frame vertex normals (set by `CalculateNormals(0)`,
    /// round-tripped via run_flags bit 1, #1718). Consumers should treat the
    /// slot as normals and skip re-applying run_transform to it.
    pub fn has_normals(&self, run: usize) -> bool {
        run < self.run_flags.len() && (self.run_flags[run] & 2) != 0
    }

    /// Applies run transforms to normals stored at `normal_idx` in each vertex's properties,
    /// then clears run_transform and run_flags. Matches C++ MeshGL::UpdateNormals(normalIdx).
    ///
    /// The normal transform is the inverse-transpose of the 3×3 rotation part of the run
    /// transform. For backside runs (run_flags bit 0 set), normals are additionally negated.
    pub fn update_normals(&mut self, normal_idx: usize) {
        if normal_idx < 3 || normal_idx + 3 > self.num_prop as usize {
            return;
        }
        let num_vert = self.num_vert();
        let num_run = self.run_original_id.len();
        let np = self.num_prop as usize;
        let mut vert_updated = vec![false; num_vert];

        for run in 0..num_run {
            // Build the 3x3 normal transform from the column-major 3x4 run transform
            let offset = 12 * run;
            let has_transform = offset + 12 <= self.run_transform.len();

            // Extract mat3 (upper-left 3x3 of the 3x4 transform)
            let (m00, m01, m02,
                 m10, m11, m12,
                 m20, m21, m22) = if has_transform {
                let t = &self.run_transform[offset..offset + 12];
                // Column-major: col0=[t0,t1,t2], col1=[t3,t4,t5], col2=[t6,t7,t8]
                (t[0] as f64, t[3] as f64, t[6] as f64,
                 t[1] as f64, t[4] as f64, t[7] as f64,
                 t[2] as f64, t[5] as f64, t[8] as f64)
            } else {
                (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
            };

            // Normal transform = inverse(transpose(M)) = (M^T)^{-1}
            // For a rotation matrix R: (R^T)^{-1} = R itself.
            // For a general transform with scale s: det = s^3, inv_trans = M / s^2.
            // We compute full adjugate/determinant to match C++ la::inverse(la::transpose(M)).
            let det = m00*(m11*m22 - m12*m21) - m01*(m10*m22 - m12*m20) + m02*(m10*m21 - m11*m20);
            let (n00, n01, n02, n10, n11, n12, n20, n21, n22) = if det.abs() < 1e-30 {
                (1.0f64, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
            } else {
                let inv = 1.0 / det;
                // Adjugate of transpose(M) = transpose of adjugate(M)
                let a00 = (m11*m22 - m12*m21) * inv;
                let a01 = (m02*m21 - m01*m22) * inv;
                let a02 = (m01*m12 - m02*m11) * inv;
                let a10 = (m12*m20 - m10*m22) * inv;
                let a11 = (m00*m22 - m02*m20) * inv;
                let a12 = (m02*m10 - m00*m12) * inv;
                let a20 = (m10*m21 - m11*m20) * inv;
                let a21 = (m01*m20 - m00*m21) * inv;
                let a22 = (m00*m11 - m01*m10) * inv;
                (a00, a01, a02, a10, a11, a12, a20, a21, a22)
            };

            let sign = if self.backside(run) { -1.0f64 } else { 1.0 };

            // Determine run's vertex range
            let start = if run < self.run_index.len() { self.run_index[run] as usize } else { 0 };
            let end = if run + 1 < self.run_index.len() { self.run_index[run + 1] as usize } else { self.tri_verts.len() };

            for idx in (start..end).step_by(1) {
                let vert = self.tri_verts[idx] as usize;
                if vert >= num_vert || vert_updated[vert] { continue; }
                vert_updated[vert] = true;
                let prop_start = vert * np + normal_idx;
                let nx = self.vert_properties[prop_start] as f64;
                let ny = self.vert_properties[prop_start + 1] as f64;
                let nz = self.vert_properties[prop_start + 2] as f64;
                // Apply normal transform
                let tx = n00*nx + n01*ny + n02*nz;
                let ty = n10*nx + n11*ny + n12*nz;
                let tz = n20*nx + n21*ny + n22*nz;
                // SafeNormalize
                let len = (tx*tx + ty*ty + tz*tz).sqrt();
                let (tx, ty, tz) = if len > 0.0 {
                    (sign * tx / len, sign * ty / len, sign * tz / len)
                } else {
                    (0.0, 0.0, 0.0)
                };
                self.vert_properties[prop_start] = tx as f32;
                self.vert_properties[prop_start + 1] = ty as f32;
                self.vert_properties[prop_start + 2] = tz as f32;
            }
        }
        self.run_transform.clear();
        self.run_flags.clear();
    }
}

impl MeshGLP<f64, u64> {
    pub fn num_vert(&self) -> usize {
        if self.num_prop == 0 { 0 } else { self.vert_properties.len() / self.num_prop as usize }
    }

    pub fn num_tri(&self) -> usize {
        self.tri_verts.len() / 3
    }

    pub fn get_vert_pos(&self, v: usize) -> [f64; 3] {
        let offset = v * self.num_prop as usize;
        [self.vert_properties[offset], self.vert_properties[offset + 1], self.vert_properties[offset + 2]]
    }

    pub fn get_tri_verts(&self, t: usize) -> [u64; 3] {
        let offset = 3 * t;
        [self.tri_verts[offset], self.tri_verts[offset + 1], self.tri_verts[offset + 2]]
    }

    pub fn get_tangent(&self, h: usize) -> [f64; 4] {
        let offset = 4 * h;
        [
            self.halfedge_tangent[offset],
            self.halfedge_tangent[offset + 1],
            self.halfedge_tangent[offset + 2],
            self.halfedge_tangent[offset + 3],
        ]
    }
}

/// Single-precision mesh (standard for graphics).
pub type MeshGL = MeshGLP<f32, u32>;

/// Double-precision, 64-bit index mesh (for huge meshes).
pub type MeshGL64 = MeshGLP<f64, u64>;

// ---------------------------------------------------------------------------
// RayHit (from include/manifold/common.h)
// ---------------------------------------------------------------------------

/// Result of a RayCast query: a single triangle-ray intersection.
#[derive(Clone, Debug, Default)]
pub struct RayHit {
    /// Triangle index that was hit.
    pub face_id: u64,
    /// Parametric distance along the ray segment in [0, 1].
    /// 0 = origin, 1 = endpoint.
    pub distance: f64,
    /// 3D position of the hit point.
    pub position: crate::linalg::Vec3,
    /// Geometric face normal at the hit.
    pub normal: crate::linalg::Vec3,
}

// ---------------------------------------------------------------------------
// Halfedge (from src/shared.h)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct Halfedge {
    pub start_vert: i32,
    pub end_vert: i32,
    pub paired_halfedge: i32,
    pub prop_vert: i32,
}

impl Halfedge {
    pub fn is_forward(&self) -> bool {
        self.start_vert < self.end_vert
    }
}

impl PartialOrd for Halfedge {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for Halfedge {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        if self.start_vert == other.start_vert {
            self.end_vert.cmp(&other.end_vert)
        } else {
            self.start_vert.cmp(&other.start_vert)
        }
    }
}

// ---------------------------------------------------------------------------
// Barycentric (from src/shared.h)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Barycentric {
    pub tri: i32,
    pub uvw: Vec4,
}

// ---------------------------------------------------------------------------
// TriRef (from src/shared.h)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct TriRef {
    /// Unique ID of the mesh instance of this triangle.
    pub mesh_id: i32,
    /// OriginalID of the mesh this triangle came from.
    pub original_id: i32,
    /// Source face ID.
    pub face_id: i32,
    /// Triangles with same coplanar_id are coplanar.
    pub coplanar_id: i32,
}

impl TriRef {
    pub fn same_face(&self, other: &TriRef) -> bool {
        self.mesh_id == other.mesh_id
            && self.coplanar_id == other.coplanar_id
            && self.face_id == other.face_id
    }
}

// ---------------------------------------------------------------------------
// TmpEdge (from src/shared.h)
// ---------------------------------------------------------------------------

#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct TmpEdge {
    pub first: i32,
    pub second: i32,
    pub halfedge_idx: i32,
}

impl TmpEdge {
    pub fn new(start: i32, end: i32, idx: i32) -> Self {
        TmpEdge {
            first: start.min(end),
            second: start.max(end),
            halfedge_idx: idx,
        }
    }
}

impl PartialOrd for TmpEdge {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for TmpEdge {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        if self.first == other.first {
            self.second.cmp(&other.second)
        } else {
            self.first.cmp(&other.first)
        }
    }
}

// ---------------------------------------------------------------------------
// MeshRelationD (from src/impl.h)
// ---------------------------------------------------------------------------

/// Transform relation between meshes.
#[derive(Clone, Debug)]
pub struct Relation {
    pub original_id: i32,
    pub transform: Mat3x4,
    pub back_side: bool,
    /// True when this meshID's contribution to `properties_` slots 0..2 holds
    /// world-frame vertex normals (set by `CalculateNormals` at slot 0).
    /// Carries through Transforms and Booleans; exported as run_flags bit 1.
    /// Per C++ #1718.
    pub has_normals: bool,
}

impl Default for Relation {
    fn default() -> Self {
        Relation {
            original_id: -1,
            transform: Mat3x4::identity(),
            back_side: false,
            has_normals: false,
        }
    }
}

impl Relation {
    /// Normal transform: inverse-transpose of the 3×3 linear part.
    /// Multiply stored-property normals by this to get world-space normals.
    /// Matches C++ Relation::GetNormalTransform()
    pub fn get_normal_transform(&self) -> Mat3 {
        let sign = if self.back_side { -1.0 } else { 1.0 };
        // NormalTransform(M) = inverse(transpose(M)) = (M^T)^{-1}
        self.transform.rotation().transpose().inverse() * sign
    }

    /// Inverse normal transform: transpose of the 3×3 linear part.
    /// Multiply world-space normals by this before storing in properties.
    /// Matches C++ Relation::GetInverseNormalTransform()
    pub fn get_inverse_normal_transform(&self) -> Mat3 {
        let sign = if self.back_side { -1.0 } else { 1.0 };
        // InverseNormalTransform(M) = M^T
        self.transform.rotation().transpose() * sign
    }
}

/// Mesh relation table stored on ManifoldImpl.
#[derive(Clone, Debug, Default)]
pub struct MeshRelationD {
    /// originalID of this Manifold if it is an original; -1 otherwise.
    pub original_id: i32,
    // C++ uses std::map (ordered by meshID); several sites iterate this map
    // and feed the order into output runs and fresh-ID assignment, so an
    // unordered map here breaks determinism and C++ parity.
    pub mesh_id_transform: BTreeMap<i32, Relation>,
    pub tri_ref: Vec<TriRef>,
}

impl MeshRelationD {
    pub fn new() -> Self {
        MeshRelationD {
            original_id: -1,
            mesh_id_transform: BTreeMap::new(),
            tri_ref: Vec::new(),
        }
    }
}

// ---------------------------------------------------------------------------
// Inline utility from shared.h
// ---------------------------------------------------------------------------

/// Return next halfedge index within the same triangle (wraps 0→1→2→0).
#[inline]
pub fn next_halfedge(current: i32) -> i32 {
    let n = current + 1;
    if n % 3 == 0 { n - 3 } else { n }
}

/// Returns the previous halfedge index within the same triangle.
/// For triangle t: PrevHalfedge(3t+i) = 3t + (i+2)%3
pub fn prev_halfedge(current: i32) -> i32 {
    let base = current - (current % 3);
    let pos = (current % 3 + 2) % 3;
    base + pos
}

/// Return next index within 0..3 (wraps 0→1→2→0).
#[inline]
pub fn next3(i: i32) -> i32 {
    let n = i + 1;
    if n == 3 { 0 } else { n }
}

// ---------------------------------------------------------------------------
#[cfg(test)]
#[path = "types_tests.rs"]
mod tests;