Skip to main content

manifold_rust/
impl_mesh.rs

1// Phase 4: Mesh Data Structure — ported from src/impl.h, src/impl.cpp, src/properties.cpp
2//
3// This module implements the core ManifoldImpl struct: halfedge mesh representation,
4// bounding box, epsilon, manifold checks, and shape constructors.
5//
6// Phases 5-9 will fill in SortGeometry, CleanupTopology, SetNormalsAndCoplanar, etc.
7
8use std::sync::atomic::{AtomicU32, Ordering};
9use crate::linalg::{Vec3, Vec4, Mat3x4, IVec3, cross, dot, normalize, length2};
10use crate::types::{
11    Box as BBox, Error, Halfedge, MeshRelationD, Relation, TriRef, K_PRECISION,
12};
13
14// ---------------------------------------------------------------------------
15// Global mesh ID counter (mirrors Manifold::Impl::meshIDCounter_)
16// ---------------------------------------------------------------------------
17
18static MESH_ID_COUNTER: AtomicU32 = AtomicU32::new(1);
19
20pub fn reserve_ids(n: u32) -> u32 {
21    MESH_ID_COUNTER.fetch_add(n, Ordering::Relaxed)
22}
23
24// ---------------------------------------------------------------------------
25// Helpers
26// ---------------------------------------------------------------------------
27
28pub const K_REMOVED_HALFEDGE: i32 = -2;
29
30/// Next halfedge within the same triangle: 0→1→2→0.
31#[inline]
32pub fn next_halfedge(current: i32) -> i32 {
33    let n = current + 1;
34    if n % 3 == 0 { n - 3 } else { n }
35}
36
37#[inline]
38pub fn next3(i: usize) -> usize {
39    if i == 2 { 0 } else { i + 1 }
40}
41
42/// Safe normalize: returns zero vector if input is zero or non-finite.
43fn safe_normalize(v: Vec3) -> Vec3 {
44    let n = normalize(v);
45    if n.x.is_finite() { n } else { Vec3::new(0.0, 0.0, 0.0) }
46}
47
48fn max_epsilon(min_epsilon: f64, bbox: &BBox) -> f64 {
49    let epsilon = min_epsilon.max(K_PRECISION * bbox.scale());
50    if epsilon.is_finite() { epsilon } else { -1.0 }
51}
52
53// ---------------------------------------------------------------------------
54// ManifoldImpl — the core mesh representation
55// ---------------------------------------------------------------------------
56
57/// Internal halfedge mesh representation, mirroring `Manifold::Impl` in C++.
58#[derive(Clone)]
59pub struct ManifoldImpl {
60    pub bbox: BBox,
61    pub epsilon: f64,
62    pub tolerance: f64,
63    pub num_prop: usize,
64    pub status: Error,
65    pub vert_pos: Vec<Vec3>,
66    pub halfedge: Vec<Halfedge>,
67    pub properties: Vec<f64>,
68    pub vert_normal: Vec<Vec3>,
69    pub face_normal: Vec<Vec3>,
70    pub halfedge_tangent: Vec<Vec4>,
71    pub mesh_relation: MeshRelationD,
72    // Collider will be added in Phase 10
73}
74
75impl Default for ManifoldImpl {
76    fn default() -> Self {
77        ManifoldImpl {
78            bbox: BBox::new(),
79            epsilon: -1.0,
80            tolerance: -1.0,
81            num_prop: 0,
82            status: Error::NoError,
83            vert_pos: Vec::new(),
84            halfedge: Vec::new(),
85            properties: Vec::new(),
86            vert_normal: Vec::new(),
87            face_normal: Vec::new(),
88            halfedge_tangent: Vec::new(),
89            mesh_relation: MeshRelationD::new(),
90        }
91    }
92}
93
94impl ManifoldImpl {
95    pub fn new() -> Self {
96        Self::default()
97    }
98
99    // -----------------------------------------------------------------------
100    // Basic accessors
101    // -----------------------------------------------------------------------
102
103    pub fn num_vert(&self) -> usize {
104        self.vert_pos.len()
105    }
106
107    pub fn num_halfedge(&self) -> usize {
108        self.halfedge.len()
109    }
110
111    pub fn num_edge(&self) -> usize {
112        self.halfedge.len() / 2
113    }
114
115    pub fn num_tri(&self) -> usize {
116        self.halfedge.len() / 3
117    }
118
119    pub fn num_prop_vert(&self) -> usize {
120        if self.num_prop == 0 {
121            self.num_vert()
122        } else {
123            self.properties.len() / self.num_prop
124        }
125    }
126
127    pub fn is_empty(&self) -> bool {
128        self.num_tri() == 0
129    }
130
131    // -----------------------------------------------------------------------
132    // MakeEmpty
133    // -----------------------------------------------------------------------
134
135    pub fn make_empty(&mut self, status: Error) {
136        self.bbox = BBox::new();
137        self.vert_pos.clear();
138        self.halfedge.clear();
139        self.vert_normal.clear();
140        self.face_normal.clear();
141        self.halfedge_tangent.clear();
142        self.mesh_relation = MeshRelationD::new();
143        self.status = status;
144    }
145
146    // -----------------------------------------------------------------------
147    // ForVert — iterate halfedges around a vertex
148    // -----------------------------------------------------------------------
149
150    /// Apply `func` to each halfedge index around the vertex starting from `halfedge_idx`.
151    pub fn for_vert<F: FnMut(usize)>(&self, halfedge_idx: usize, mut func: F) {
152        let mut current = halfedge_idx;
153        loop {
154            current = next_halfedge(self.halfedge[current].paired_halfedge) as usize;
155            func(current);
156            if current == halfedge_idx {
157                break;
158            }
159        }
160    }
161
162    // -----------------------------------------------------------------------
163    // CalculateBBox
164    // -----------------------------------------------------------------------
165
166    pub fn calculate_bbox(&mut self) {
167        let mut bbox = BBox::new();
168        for v in &self.vert_pos {
169            if !v.x.is_nan() {
170                bbox.union_point(*v);
171            }
172        }
173        self.bbox = bbox;
174        if !self.bbox.is_finite() {
175            self.make_empty(Error::NoError);
176        }
177    }
178
179    // -----------------------------------------------------------------------
180    // SetEpsilon
181    // -----------------------------------------------------------------------
182
183    pub fn set_epsilon(&mut self, min_epsilon: f64, use_single: bool) {
184        self.epsilon = max_epsilon(min_epsilon, &self.bbox);
185        let mut min_tol = self.epsilon;
186        if use_single {
187            let float_eps = (f32::EPSILON as f64) * self.bbox.scale();
188            min_tol = min_tol.max(float_eps);
189        }
190        self.tolerance = self.tolerance.max(min_tol);
191    }
192
193    // -----------------------------------------------------------------------
194    // IsFinite
195    // -----------------------------------------------------------------------
196
197    pub fn is_finite(&self) -> bool {
198        self.vert_pos.iter().all(|v| v.x.is_finite() && v.y.is_finite() && v.z.is_finite())
199    }
200
201    // -----------------------------------------------------------------------
202    // IsManifold / Is2Manifold
203    // -----------------------------------------------------------------------
204
205    /// Check that the halfedge data structure is consistent (oriented even manifold).
206    pub fn is_manifold(&self) -> bool {
207        if self.halfedge.is_empty() {
208            return true;
209        }
210        if self.halfedge.len() % 3 != 0 {
211            return false;
212        }
213        for (edge, h) in self.halfedge.iter().enumerate() {
214            // Valid removed halfedge
215            if h.start_vert == -1 && h.end_vert == -1 && h.paired_halfedge == -1 {
216                continue;
217            }
218            // Neighbors in same triangle must not be removed
219            let n1 = next_halfedge(edge as i32) as usize;
220            let n2 = next_halfedge(n1 as i32) as usize;
221            if self.halfedge[n1].start_vert == -1 || self.halfedge[n2].start_vert == -1 {
222                return false;
223            }
224            if h.paired_halfedge == -1 {
225                return false;
226            }
227            let paired_idx = h.paired_halfedge as usize;
228            let paired = &self.halfedge[paired_idx];
229            if paired.paired_halfedge != edge as i32 {
230                return false;
231            }
232            if h.start_vert == h.end_vert {
233                return false;
234            }
235            if h.start_vert != paired.end_vert {
236                return false;
237            }
238            if h.end_vert != paired.start_vert {
239                return false;
240            }
241        }
242        true
243    }
244
245    /// Check that the mesh is a 2-manifold (no duplicate edges).
246    pub fn is_2_manifold(&self) -> bool {
247        if self.halfedge.is_empty() {
248            return true;
249        }
250        if !self.is_manifold() {
251            return false;
252        }
253        // Sort halfedges and check for duplicates
254        let mut sorted = self.halfedge.clone();
255        sorted.sort_unstable();
256        for i in 0..sorted.len().saturating_sub(1) {
257            let h = &sorted[i];
258            let h1 = &sorted[i + 1];
259            // Skip removed halfedges
260            if h.start_vert == -1 && h.end_vert == -1 && h.paired_halfedge == -1 {
261                continue;
262            }
263            if h.start_vert == h1.start_vert && h.end_vert == h1.end_vert {
264                return false; // Duplicate edge
265            }
266        }
267        true
268    }
269
270    // -----------------------------------------------------------------------
271    // CreateHalfedges
272    // -----------------------------------------------------------------------
273
274    /// Build the halfedge data structure from triangle lists.
275    ///
276    /// - `tri_prop`: property vertex indices per triangle (also geometry if `tri_vert` is empty)
277    /// - `tri_vert`: geometry vertex indices per triangle (may be empty)
278    ///
279    /// When `tri_vert` is empty, `tri_prop` is used for both geometry and properties.
280    /// When `tri_vert` is present, `tri_prop[i][j]` = `propVert`, `tri_vert[i][j]` = `startVert`.
281    pub fn create_halfedges(&mut self, tri_prop: &[IVec3], tri_vert: &[IVec3]) {
282        let num_tri = tri_prop.len();
283        if num_tri == 0 {
284            self.halfedge.clear();
285            return;
286        }
287        let num_halfedge = 3 * num_tri;
288        let num_edge = num_halfedge / 2;
289
290        self.halfedge.clear();
291        self.halfedge.resize(num_halfedge, Halfedge {
292            start_vert: -1,
293            end_vert: -1,
294            paired_halfedge: -1,
295            prop_vert: -1,
296        });
297
298        let use_prop = tri_vert.is_empty();
299
300        // Build halfedges and compute edge sort key
301        // key = [forward_bit:1][min_vert:31][max_vert:32]
302        // forward: v0 < v1 → bit=1; backward: v0 > v1 → bit=0
303        // After sorting: backward halfedges first, then forward, both by (min,max)
304        let mut edge_keys = vec![0u64; num_halfedge];
305
306        for tri in 0..num_tri {
307            for i in 0usize..3 {
308                let j = next3(i);
309                let e = 3 * tri + i;
310                let v0 = if use_prop { tri_prop[tri][i] } else { tri_vert[tri][i] };
311                let v1 = if use_prop { tri_prop[tri][j] } else { tri_vert[tri][j] };
312                self.halfedge[e] = Halfedge {
313                    start_vert: v0,
314                    end_vert: v1,
315                    paired_halfedge: -1,
316                    prop_vert: tri_prop[tri][i],
317                };
318                let fwd = if v0 < v1 { 1u64 } else { 0u64 };
319                let min_v = v0.min(v1) as u64;
320                let max_v = v0.max(v1) as u64;
321                edge_keys[e] = (fwd << 63) | (min_v << 32) | max_v;
322            }
323        }
324
325        // Sort halfedge indices by edge key. C++ CreateHalfedges uses a
326        // STABLE sort here (impl.cpp), and the #1687 fix ensures its parallel
327        // stable_sort matches std::stable_sort. When two halfedges share an
328        // edge key (duplicate directed edges in degenerate/intermediate
329        // meshes) the tie must break on original halfedge-index order, so we
330        // use a stable sort to stay bit-identical to C++.
331        let mut ids: Vec<usize> = (0..num_halfedge).collect();
332        ids.sort_by_key(|&i| edge_keys[i]);
333
334        // ids[0..num_edge] = backward halfedges (startVert > endVert), sorted by (min,max)
335        // ids[num_edge..] = forward halfedges (startVert < endVert), sorted by (min,max)
336
337        // Sequential pairing with opposed-triangle detection
338        let segment_end = num_edge;
339        let mut consecutive_start = 0usize;
340
341        for i in 0..num_edge {
342            let pair0 = ids[i];
343            let h0_sv = self.halfedge[pair0].start_vert;
344            let h0_ev = self.halfedge[pair0].end_vert;
345
346            let mut k = consecutive_start + num_edge;
347            'inner: loop {
348                if k >= segment_end + num_edge {
349                    break 'inner;
350                }
351                let pair1 = ids[k];
352                let h1_sv = self.halfedge[pair1].start_vert;
353                let h1_ev = self.halfedge[pair1].end_vert;
354
355                if h0_sv != h1_ev || h0_ev != h1_sv {
356                    break 'inner; // Different edge direction — no match
357                }
358
359                if self.halfedge[pair1].paired_halfedge != K_REMOVED_HALFEDGE {
360                    // Check for opposed triangle: same undirected edge, same third vertex
361                    let next0 = next_halfedge(pair0 as i32) as usize;
362                    let next1 = next_halfedge(pair1 as i32) as usize;
363                    if self.halfedge[next0].end_vert == self.halfedge[next1].end_vert {
364                        // Opposed triangles: mark both for removal.
365                        // Reorder ids so the remaining valid forward halfedge (at i+num_edge)
366                        // moves to position k, and pair1 (the opposed one) goes to i+num_edge.
367                        // This matches C++ which does: ids[k] = ids[i+numEdge]; ids[i+numEdge] = pair1;
368                        self.halfedge[pair0].paired_halfedge = K_REMOVED_HALFEDGE;
369                        self.halfedge[pair1].paired_halfedge = K_REMOVED_HALFEDGE;
370                        if i + num_edge != k {
371                            ids.swap(k, i + num_edge);
372                        }
373                        break 'inner;
374                    }
375                }
376
377                k += 1;
378            }
379
380            // Update consecutive_start for next iteration
381            if i + 1 < segment_end {
382                let next_sv = self.halfedge[ids[i + 1]].start_vert;
383                let next_ev = self.halfedge[ids[i + 1]].end_vert;
384                if next_sv != h0_sv || next_ev != h0_ev {
385                    consecutive_start = i + 1;
386                }
387            }
388        }
389
390        // Final pairing pass
391        for i in 0..num_edge {
392            let pair0 = ids[i];
393            let pair1 = ids[i + num_edge];
394            if self.halfedge[pair0].paired_halfedge != K_REMOVED_HALFEDGE {
395                self.halfedge[pair0].paired_halfedge = pair1 as i32;
396                self.halfedge[pair1].paired_halfedge = pair0 as i32;
397            } else {
398                // Invalidate both (opposed triangles removed)
399                self.halfedge[pair0] = Halfedge { start_vert: -1, end_vert: -1, paired_halfedge: -1, prop_vert: 0 };
400                self.halfedge[pair1] = Halfedge { start_vert: -1, end_vert: -1, paired_halfedge: -1, prop_vert: 0 };
401            }
402        }
403    }
404
405    // -----------------------------------------------------------------------
406    // InitializeOriginal
407    // -----------------------------------------------------------------------
408
409    /// Set up the mesh relation for a newly created original mesh.
410    pub fn initialize_original(&mut self) {
411        // Per C++ #1718: preserve the AND-across-old-Relations hasNormals state
412        // so AsOriginal keeps the recording when it builds a fresh Relation.
413        // Primitives start with an empty map → all_have_normals() is false.
414        let had_normals = self.all_have_normals();
415        let mesh_id = reserve_ids(1) as i32;
416        self.mesh_relation.original_id = mesh_id;
417        let num_tri = self.num_tri();
418        self.mesh_relation.tri_ref.resize(num_tri, TriRef::default());
419        for (tri, tri_ref) in self.mesh_relation.tri_ref.iter_mut().enumerate() {
420            tri_ref.mesh_id = mesh_id;
421            tri_ref.original_id = mesh_id;
422            tri_ref.face_id = -1;
423            tri_ref.coplanar_id = tri as i32;
424        }
425        self.mesh_relation.mesh_id_transform.clear();
426        self.mesh_relation.mesh_id_transform.insert(mesh_id, Relation {
427            original_id: mesh_id,
428            transform: Mat3x4::identity(),
429            back_side: false,
430            has_normals: had_normals,
431        });
432    }
433
434    /// True only when every meshID carries normals at slot 0..2 — the condition
435    /// under which `get_mesh_gl(-1)` can safely auto-substitute that slot. A
436    /// mixed Boolean output (some meshIDs with normals, some without) returns
437    /// false; the output MeshGL's per-run bit 1 still marks the with-normals
438    /// runs individually. AND semantics across meshIDs. Per C++ #1718.
439    pub fn all_have_normals(&self) -> bool {
440        let map = &self.mesh_relation.mesh_id_transform;
441        !map.is_empty() && map.values().all(|m| m.has_normals)
442    }
443
444    /// True iff the meshID owning `tri` has hasNormals set. False when the
445    /// meshID isn't in mesh_id_transform (treat as no-normals). Per C++ #1718.
446    pub fn tri_has_normals(&self, tri: usize) -> bool {
447        let mesh_id = self.mesh_relation.tri_ref[tri].mesh_id;
448        self.mesh_relation
449            .mesh_id_transform
450            .get(&mesh_id)
451            .map(|m| m.has_normals)
452            .unwrap_or(false)
453    }
454
455    /// Eager-transform slot 0..2 of `properties` for propVerts whose meshID
456    /// carries hasNormals. Used by both `transform` and `compose` so world-frame
457    /// normals stay in sync with vert_pos / face_normal across any sequence of
458    /// transforms (including mixed-input Boolean/Compose outputs where some
459    /// meshIDs carry normals and others don't). Per C++ #1718.
460    ///
461    /// `properties` is laid out as `properties[(offset + prop) * stride + i]`,
462    /// so callers can target an in-place properties_ vector (offset=0) or a
463    /// per-node slice of a combined array (offset=propVertIndices,
464    /// stride=numPropOut). Re-normalizes as it transforms so non-orthogonal
465    /// transforms (scale) and upstream barycentric interpolation don't leave
466    /// non-unit values that compound downstream.
467    pub fn eager_transform_prop_normals(
468        halfedge: &[crate::types::Halfedge],
469        mesh_relation: &crate::types::MeshRelationD,
470        normal_transform: crate::linalg::Mat3,
471        properties: &mut [f64],
472        num_prop_vert: usize,
473        stride: usize,
474        offset: usize,
475    ) {
476        // OR semantics (any meshID has normals), unlike all_have_normals():
477        // mixed inputs still need the per-meshID iteration below to rotate the
478        // with-normals subset.
479        if !mesh_relation.mesh_id_transform.values().any(|m| m.has_normals) {
480            return;
481        }
482        let tri_has_normals = |tri: usize| -> bool {
483            let mid = mesh_relation.tri_ref[tri].mesh_id;
484            mesh_relation
485                .mesh_id_transform
486                .get(&mid)
487                .map(|m| m.has_normals)
488                .unwrap_or(false)
489        };
490        let mut visited = vec![false; num_prop_vert];
491        for e in 0..halfedge.len() {
492            if !tri_has_normals(e / 3) {
493                continue;
494            }
495            let prop = halfedge[e].prop_vert;
496            if prop < 0 {
497                continue;
498            }
499            let prop = prop as usize;
500            if visited[prop] {
501                continue;
502            }
503            visited[prop] = true;
504            let base = (offset + prop) * stride;
505            let n = Vec3::new(properties[base], properties[base + 1], properties[base + 2]);
506            let nt = safe_normalize(normal_transform * n);
507            properties[base] = nt.x;
508            properties[base + 1] = nt.y;
509            properties[base + 2] = nt.z;
510        }
511    }
512
513    // -----------------------------------------------------------------------
514    // IncrementMeshIDs — port of C++ Manifold::Impl::IncrementMeshIDs()
515    // -----------------------------------------------------------------------
516
517    /// Allocates fresh unique mesh IDs and remaps all triRef.meshID values.
518    /// This ensures boolean results don't collide with source mesh IDs.
519    pub fn increment_mesh_ids(&mut self) {
520        use std::collections::{BTreeMap, HashMap};
521
522        // Build old -> new ID mapping. Iteration order determines which old
523        // ID gets which fresh ID, so it must be sorted like C++ std::map.
524        let old_transforms: BTreeMap<i32, Relation> =
525            std::mem::take(&mut self.mesh_relation.mesh_id_transform);
526        let num_mesh_ids = old_transforms.len() as u32;
527        if num_mesh_ids == 0 { return; }
528        let mut next_mesh_id = reserve_ids(num_mesh_ids) as i32;
529        let mut old2new: HashMap<i32, i32> = HashMap::new();
530        for (old_id, relation) in old_transforms {
531            old2new.insert(old_id, next_mesh_id);
532            self.mesh_relation.mesh_id_transform.insert(next_mesh_id, relation);
533            next_mesh_id += 1;
534        }
535
536        // Update all triRef.meshID
537        for tri_ref in &mut self.mesh_relation.tri_ref {
538            if let Some(&new_id) = old2new.get(&tri_ref.mesh_id) {
539                tri_ref.mesh_id = new_id;
540            }
541        }
542    }
543
544    // -----------------------------------------------------------------------
545    // DedupePropVerts — port of C++ Manifold::Impl::DedupePropVerts()
546    // -----------------------------------------------------------------------
547
548    /// Deduplicates property vertices that share identical property values
549    /// across paired halfedges within the same mesh.
550    pub fn dedupe_prop_verts(&mut self) {
551        let num_prop = self.num_prop;
552        if num_prop == 0 { return; }
553
554        let n_edges = self.halfedge.len();
555        // Collect (prop0, prop1) pairs for edges where properties match
556        let mut vert2vert: Vec<(i32, i32)> = vec![(-1, -1); n_edges];
557        for edge_idx in 0..n_edges {
558            let edge = self.halfedge[edge_idx];
559            if edge.paired_halfedge < 0 { continue; }
560            let edge_face = edge_idx / 3;
561            let pair_face = edge.paired_halfedge as usize / 3;
562
563            if self.mesh_relation.tri_ref[edge_face].mesh_id
564                != self.mesh_relation.tri_ref[pair_face].mesh_id
565            {
566                continue;
567            }
568
569            let prop0 = self.halfedge[edge_idx].prop_vert;
570            let prop1 = self.halfedge[next_halfedge(edge.paired_halfedge) as usize].prop_vert;
571            if prop0 < 0 || prop1 < 0 { continue; }
572
573            let mut prop_equal = true;
574            for p in 0..num_prop {
575                let idx0 = num_prop * prop0 as usize + p;
576                let idx1 = num_prop * prop1 as usize + p;
577                if idx0 >= self.properties.len() || idx1 >= self.properties.len() {
578                    prop_equal = false;
579                    break;
580                }
581                if self.properties[idx0] != self.properties[idx1] {
582                    prop_equal = false;
583                    break;
584                }
585            }
586            if prop_equal {
587                vert2vert[edge_idx] = (prop0, prop1);
588            }
589        }
590
591        // Use union-find to merge equivalent property vertices
592        let num_prop_vert = self.num_prop_vert();
593        let mut ds = crate::disjoint_sets::DisjointSets::new(num_prop_vert as u32);
594        for &(a, b) in &vert2vert {
595            if a >= 0 && b >= 0 {
596                ds.unite(a as u32, b as u32);
597            }
598        }
599        let mut vert_labels = Vec::new();
600        let num_labels = ds.connected_components(&mut vert_labels);
601
602        // Build label -> canonical vert mapping
603        let mut label2vert = vec![0i32; num_labels as usize];
604        for v in 0..num_prop_vert {
605            label2vert[vert_labels[v] as usize] = v as i32;
606        }
607
608        // Remap all prop_vert indices
609        for edge in &mut self.halfedge {
610            if edge.prop_vert >= 0 && (edge.prop_vert as usize) < num_prop_vert {
611                edge.prop_vert = label2vert[vert_labels[edge.prop_vert as usize] as usize];
612            }
613        }
614    }
615
616    // -----------------------------------------------------------------------
617    // RemoveUnreferencedVerts
618    // -----------------------------------------------------------------------
619
620    /// Mark unreferenced vertices as NaN (to be cleaned up by later passes).
621    pub fn remove_unreferenced_verts(&mut self) {
622        let num_vert = self.num_vert();
623        let mut keep = vec![false; num_vert];
624        for h in &self.halfedge {
625            if h.start_vert >= 0 {
626                keep[h.start_vert as usize] = true;
627            }
628        }
629        for (i, k) in keep.iter().enumerate() {
630            if !k {
631                self.vert_pos[i] = Vec3::new(f64::NAN, f64::NAN, f64::NAN);
632            }
633        }
634    }
635
636    // -----------------------------------------------------------------------
637    // SetNormalsAndCoplanar (stub — implemented in Phase 9)
638    // -----------------------------------------------------------------------
639
640    /// Compute face normals, assign coplanar IDs, and calculate vertex normals.
641    pub fn set_normals_and_coplanar(&mut self) {
642        crate::face_op::set_normals_and_coplanar(self);
643    }
644
645    // -----------------------------------------------------------------------
646    // SortGeometry (stub — implemented in Phase 5)
647    // -----------------------------------------------------------------------
648
649    /// Reorder mesh geometry for cache efficiency using Morton codes.
650    pub fn sort_geometry(&mut self) {
651        crate::sort::sort_geometry(self);
652    }
653
654    // -----------------------------------------------------------------------
655    // Shape constructors
656    // -----------------------------------------------------------------------
657
658    pub fn tetrahedron(transform: &Mat3x4) -> Self {
659        let vert_pos_raw: Vec<[f64; 3]> = vec![
660            [-1.0, -1.0,  1.0],
661            [-1.0,  1.0, -1.0],
662            [ 1.0, -1.0, -1.0],
663            [ 1.0,  1.0,  1.0],
664        ];
665        let tri_verts: Vec<IVec3> = vec![
666            IVec3::new(2, 0, 1),
667            IVec3::new(0, 3, 1),
668            IVec3::new(2, 3, 0),
669            IVec3::new(3, 2, 1),
670        ];
671        Self::from_shape(vert_pos_raw, tri_verts, transform)
672    }
673
674    pub fn cube(transform: &Mat3x4) -> Self {
675        let vert_pos_raw: Vec<[f64; 3]> = vec![
676            [0.0, 0.0, 0.0],
677            [0.0, 0.0, 1.0],
678            [0.0, 1.0, 0.0],
679            [0.0, 1.0, 1.0],
680            [1.0, 0.0, 0.0],
681            [1.0, 0.0, 1.0],
682            [1.0, 1.0, 0.0],
683            [1.0, 1.0, 1.0],
684        ];
685        let tri_verts: Vec<IVec3> = vec![
686            IVec3::new(1, 0, 4), IVec3::new(2, 4, 0),
687            IVec3::new(1, 3, 0), IVec3::new(3, 1, 5),
688            IVec3::new(3, 2, 0), IVec3::new(3, 7, 2),
689            IVec3::new(5, 4, 6), IVec3::new(5, 1, 4),
690            IVec3::new(6, 4, 2), IVec3::new(7, 6, 2),
691            IVec3::new(7, 3, 5), IVec3::new(7, 5, 6),
692        ];
693        Self::from_shape(vert_pos_raw, tri_verts, transform)
694    }
695
696    pub fn octahedron(transform: &Mat3x4) -> Self {
697        let vert_pos_raw: Vec<[f64; 3]> = vec![
698            [ 1.0,  0.0,  0.0],
699            [-1.0,  0.0,  0.0],
700            [ 0.0,  1.0,  0.0],
701            [ 0.0, -1.0,  0.0],
702            [ 0.0,  0.0,  1.0],
703            [ 0.0,  0.0, -1.0],
704        ];
705        let tri_verts: Vec<IVec3> = vec![
706            IVec3::new(0, 2, 4), IVec3::new(1, 5, 3),
707            IVec3::new(2, 1, 4), IVec3::new(3, 5, 0),
708            IVec3::new(1, 3, 4), IVec3::new(0, 5, 2),
709            IVec3::new(3, 0, 4), IVec3::new(2, 5, 1),
710        ];
711        Self::from_shape(vert_pos_raw, tri_verts, transform)
712    }
713
714    fn from_shape(vert_pos_raw: Vec<[f64; 3]>, tri_verts: Vec<IVec3>, transform: &Mat3x4) -> Self {
715        use crate::linalg::Vec4;
716        let mut m = Self::new();
717        m.vert_pos = vert_pos_raw
718            .iter()
719            .map(|v| {
720                let p = Vec3::new(v[0], v[1], v[2]);
721                // Apply transform: m * vec4(p, 1)
722                *transform * Vec4::new(p.x, p.y, p.z, 1.0)
723            })
724            .collect();
725
726        m.create_halfedges(&tri_verts, &[]);
727        m.initialize_original();
728        m.calculate_bbox();
729        m.set_epsilon(-1.0, false);
730        m.sort_geometry();
731        m.set_normals_and_coplanar();
732        m
733    }
734
735    // -----------------------------------------------------------------------
736    // Transform
737    // -----------------------------------------------------------------------
738
739    /// Apply affine transform, returning a new ManifoldImpl.
740    pub fn transform(&self, t: &Mat3x4) -> Self {
741        use crate::linalg::{Vec4, Mat3};
742        let identity = Mat3x4::identity();
743        if t == &identity {
744            // Clone self — this is a simplified version (full version uses Collider)
745            return self.shallow_clone();
746        }
747
748        let mut result = Self::new();
749        if self.status != Error::NoError {
750            result.status = self.status;
751            return result;
752        }
753
754        result.mesh_relation = self.mesh_relation.clone();
755        // Scale epsilon by spectral norm of transform, matching C++:
756        // result.epsilon_ *= SpectralNorm(mat3(transform_));
757        let m3_for_norm = Mat3::from_cols(
758            Vec3::new(t.x.x, t.x.y, t.x.z),
759            Vec3::new(t.y.x, t.y.y, t.y.z),
760            Vec3::new(t.z.x, t.z.y, t.z.z),
761        );
762        result.epsilon = self.epsilon * crate::svd::spectral_norm(m3_for_norm);
763        result.tolerance = self.tolerance;
764        result.num_prop = self.num_prop;
765        result.properties = self.properties.clone();
766        result.bbox = self.bbox;
767        result.halfedge = self.halfedge.clone();
768        result.mesh_relation.original_id = -1;
769
770        // Update mesh transforms
771        for (_, rel) in result.mesh_relation.mesh_id_transform.iter_mut() {
772            // rel.transform = t * Mat4(rel.transform) — combine transforms
773            rel.transform = mat3x4_mul_mat3x4(t, &rel.transform);
774        }
775
776        // Transform vertex positions
777        result.vert_pos = self.vert_pos.iter().map(|&v| {
778            *t * Vec4::new(v.x, v.y, v.z, 1.0)
779        }).collect();
780
781        // Transform normals (using inverse-transpose of 3x3 part)
782        let m3 = Mat3::from_cols(
783            Vec3::new(t.x.x, t.x.y, t.x.z),
784            Vec3::new(t.y.x, t.y.y, t.y.z),
785            Vec3::new(t.z.x, t.z.y, t.z.z),
786        );
787        let normal_t = m3.inverse().transpose();
788
789        result.face_normal = self.face_normal.iter().map(|&n| {
790            safe_normalize(normal_t * n)
791        }).collect();
792        result.vert_normal = self.vert_normal.iter().map(|&n| {
793            safe_normalize(normal_t * n)
794        }).collect();
795
796        // Per #1718: the properties clone above doesn't go through the vertPos /
797        // faceNormal transform, so eager-transform slot 0..2 per-meshID to keep
798        // recorded world-frame normals in sync. tri_ref / hasNormals flags are
799        // identical in self and result; iterate by prop vert (winding flip
800        // below only reorders halfedges, not prop assignments).
801        if self.num_prop >= 3 {
802            Self::eager_transform_prop_normals(
803                &self.halfedge,
804                &self.mesh_relation,
805                normal_t,
806                &mut result.properties,
807                self.num_prop_vert(),
808                self.num_prop,
809                0,
810            );
811        }
812
813        let invert = m3.determinant() < 0.0;
814
815        // Transform tangents — C++ TransformTangents
816        // Must happen BEFORE FlipTris (matches C++ order)
817        if !self.halfedge_tangent.is_empty() {
818            result.halfedge_tangent = vec![Vec4::new(0.0, 0.0, 0.0, 0.0); self.halfedge_tangent.len()];
819            for edge_out in 0..self.halfedge_tangent.len() {
820                let edge_in = if invert {
821                    let tri = edge_out / 3;
822                    let vert = 2 - (edge_out - 3 * tri);
823                    let flipped = 3 * tri + vert;
824                    self.halfedge[flipped].paired_halfedge as usize
825                } else {
826                    edge_out
827                };
828                let old_t = self.halfedge_tangent[edge_in];
829                let xyz = m3 * Vec3::new(old_t.x, old_t.y, old_t.z);
830                result.halfedge_tangent[edge_out] = Vec4::new(xyz.x, xyz.y, xyz.z, old_t.w);
831            }
832        }
833
834        if invert {
835            // Flip triangle winding — matches C++ FlipTris
836            for tri in 0..result.num_tri() {
837                // Swap first and third halfedge within tri
838                result.halfedge.swap(3 * tri, 3 * tri + 2);
839                // For each halfedge: swap startVert/endVert and remap pairedHalfedge
840                for i in 0..3 {
841                    let idx = 3 * tri + i;
842                    let h = &mut result.halfedge[idx];
843                    std::mem::swap(&mut h.start_vert, &mut h.end_vert);
844                    // FlipHalfedge: within the paired tri, mirror the edge index
845                    let paired = h.paired_halfedge;
846                    if paired >= 0 {
847                        let p = paired as usize;
848                        let p_tri = p / 3;
849                        let p_vert = 2 - (p - 3 * p_tri);
850                        h.paired_halfedge = (3 * p_tri + p_vert) as i32;
851                    }
852                }
853            }
854        }
855
856        result.calculate_bbox();
857        result.set_epsilon(result.epsilon, false);
858        result
859    }
860
861    /// Clone without collider (shallow copy for transform operations)
862    fn shallow_clone(&self) -> Self {
863        ManifoldImpl {
864            bbox: self.bbox,
865            epsilon: self.epsilon,
866            tolerance: self.tolerance,
867            num_prop: self.num_prop,
868            status: self.status,
869            vert_pos: self.vert_pos.clone(),
870            halfedge: self.halfedge.clone(),
871            properties: self.properties.clone(),
872            vert_normal: self.vert_normal.clone(),
873            face_normal: self.face_normal.clone(),
874            halfedge_tangent: self.halfedge_tangent.clone(),
875            mesh_relation: self.mesh_relation.clone(),
876        }
877    }
878}
879
880// ---------------------------------------------------------------------------
881// Transform helpers
882// ---------------------------------------------------------------------------
883
884/// Multiply two Mat3x4 transforms as affine matrices (t1 * t2 = (t1 * to_mat4(t2)).
885/// Result is t1 applied after t2.
886fn mat3x4_mul_mat3x4(t1: &Mat3x4, t2: &Mat3x4) -> Mat3x4 {
887    use crate::linalg::Vec4;
888    // Column vectors of t2 (as Vec4 with w=0 for rotation cols, w=1 for translation)
889    let c0 = *t1 * Vec4::new(t2.x.x, t2.x.y, t2.x.z, 0.0);
890    let c1 = *t1 * Vec4::new(t2.y.x, t2.y.y, t2.y.z, 0.0);
891    let c2 = *t1 * Vec4::new(t2.z.x, t2.z.y, t2.z.z, 0.0);
892    let c3 = *t1 * Vec4::new(t2.w.x, t2.w.y, t2.w.z, 1.0);
893    Mat3x4 { x: c0, y: c1, z: c2, w: c3 }
894}
895
896// ---------------------------------------------------------------------------
897// Mat3 for normal transform (we need inverse + transpose from linalg)
898// ---------------------------------------------------------------------------
899
900// These are already in linalg.rs but we need to use them here.
901// The Mat3 methods: inverse(), transpose(), determinant()
902
903// ---------------------------------------------------------------------------
904#[cfg(test)]
905#[path = "impl_mesh_tests.rs"]
906mod tests;