remesh 0.0.5

Isotropic remeshing library
Documentation
// SPDX-License-Identifier: MIT OR Apache-2.0
// Copyright (c) 2025 lacklustr@protonmail.com https://github.com/eadf

use crate::common::macros::{integrity_assert_eq, integrity_assert_unique, integrity_println};
use crate::corner_table::CornerIndex;
use crate::isotropic_remesh::IsotropicRemeshAlgo;
use rayon::iter::{IntoParallelIterator, ParallelIterator};
use std::fmt::Debug;
use vector_traits::num_traits::AsPrimitive;
use vector_traits::prelude::SimdUpgradable;

impl<S, V, const ENABLE_UNSAFE: bool> IsotropicRemeshAlgo<S, V, ENABLE_UNSAFE>
where
    S: crate::common::sealed::ScalarType,
    f64: AsPrimitive<S>,
    V: Debug + Copy + From<[S; 3]> + Into<[S; 3]> + Sync + 'static,
{
    /// Flip edges to improve overall valence
    /// An edge is flipped if:
    /// 1. The two triangles are nearly coplanar (configurable threshold)
    /// 2. The flip would not create degenerate triangles
    /// 3. The valence quality is improved
    pub(crate) fn flip_edges_dihedral(&mut self) -> bool {
        let num_corners = self.corner_table.vertex_of_corners().len() as u32;
        let mut total_changes = false;

        // Usually converges in 2-3 passes
        let mut pass = 1;
        while pass < 6 {
            let mut changes_this_pass = false;

            #[cfg(any(feature = "integrity_check", debug_assertions))]
            let mut last_quality = 1000_i16;

            // Iterate through all edges
            let mut candidates: Vec<_> = (0..num_corners)
                .into_par_iter()
                .map(CornerIndex)
                .filter_map(|corner| {
                    if self.corner_table.is_corner_deleted(corner) {
                        return None;
                    }

                    if corner > self.corner_table.twin(corner) {
                        return None; // Will process from twin instead
                    }

                    self.should_flip_edge_dihedral(corner)
                })
                .collect();

            // sort candidates
            candidates.sort_unstable_by(|a, b| {
                (a.2, a.3)
                    .partial_cmp(&(b.2, b.3))
                    .unwrap_or(std::cmp::Ordering::Equal)
            });

            while let Some((c0, c1, _quality, _)) = candidates.pop() {
                if c1 != self.corner_table.twin(c0) {
                    // c1 is no longer twin of c0. Could this ever happen?
                    continue;
                }
                #[cfg(any(feature = "integrity_check", debug_assertions))]
                {
                    assert!(last_quality >= _quality, "{last_quality} > {_quality}");
                    last_quality = _quality;
                }

                // recheck if we should flip this edge
                if let Some((c0, _c1, _q, _)) = self.should_flip_edge_dihedral(c0) {
                    self.flip_manifold_edge(c0);
                    changes_this_pass = true;

                    #[cfg(feature = "integrity_check")]
                    self.check_mesh_integrity(
                        format!(
                            "after flip edge {}-{} quality:{_q}",
                            self.corner_table.data.dbg_corner(c0),
                            self.corner_table.data.dbg_corner(_c1)
                        )
                        .as_str(),
                    )
                    .unwrap();
                }
            }

            if !changes_this_pass {
                break;
            }
            total_changes |= changes_this_pass;
            pass += 1;
        }
        integrity_println!("flip_edges converged after {} passes", pass);
        total_changes
    }

    /// Determines if an edge should be flipped
    /// ```text
    ///
    ///              c₀o
    ///      V₃ ─────────────── V₂        V₃ ─────────────── V₂
    ///      │c₀p         c₀n / │         │ \ c₁n        c₁  │
    ///      │              /   │         │   \              │
    ///      │     t₀     / c₁  │         │ c₀p\     t₁      │
    ///      │          /       │    =>   │      \           │
    /// c₀no │        /         │ c₁no    │        \         │
    ///      │      /           │         │          \       │
    ///      │c₀  /      t₁     │         │    t₀      \ c₁p │
    ///      │  /               │         │              \   │
    ///      │/ c₁n         c₁p │         │ c₀         c₀n \ │
    ///      V₀ ────────────── V₁         V₀ ────────────── V₁
    ///              c₁o
    ///
    /// Input: c₀: The corner index at V₀ in the V₀V₂V₃ triangle.
    /// An edge should be flipped if:
    /// 1. The two triangles are nearly coplanar (configurable threshold)
    /// 2. The flip would not create degenerate triangles
    /// 3. The valence quality is improved (lowering V₀ and V₂ valence while keeping valence of V₁ and V₃ acceptable)
    /// ```
    /// Returns `Option<(corner, corner.twin, valence improvement)>`
    pub(crate) fn should_flip_edge_dihedral(
        &self,
        c0: CornerIndex,
    ) -> Option<(CornerIndex, CornerIndex, i16, u64)> {
        // Get the four vertices of the quad formed by the two triangles
        // c0 is in triangle ABC, c1 is in triangle BAD
        let (c0n, c0p) = self.corner_table.next_prev(c0); // corner V₃, V₂ in t₀
        let c1p = self.corner_table.opposite(c0p); // corner at V₁ in t₁
        let c1 = self.corner_table.next(c1p);
        integrity_assert_eq!(c1, self.corner_table.twin(c0));

        let v0 = self.corner_table.vertex(c0);
        let v2 = self.corner_table.vertex(c0n);
        let v3 = self.corner_table.vertex(c0p);
        let v1 = self.corner_table.vertex(c1p);
        integrity_assert_eq!(v0, self.corner_table.vertex(self.corner_table.prev(c1p)));

        // Check if all vertices are valid and distinct
        integrity_assert_unique!(v0, v1, v2, v3);

        // Get actual vertex positions
        let pos_v0 = self.vertex(v0).to_simd();
        let pos_v1 = self.vertex(v1).to_simd();
        let pos_v2 = self.vertex(v2).to_simd();
        let pos_v3 = self.vertex(v3).to_simd();

        integrity_assert_unique!(pos_v0, pos_v1, pos_v2, pos_v3);

        // Check if triangles are nearly coplanar and the new triangles are ok
        if !Self::quad_is_convex_coplanar(
            pos_v0,
            pos_v1,
            pos_v2,
            pos_v3,
            self.params.coplanar_threshold_sq,
        ) {
            return None;
        }

        let valence_v0 = self.corner_table.valence(c0);
        let valence_v1 = self.corner_table.valence(c1p);
        let valence_v2 = self.corner_table.valence(c0n);
        let valence_v3 = self.corner_table.valence(c0p);

        // 5. Check valence improvement
        let valence_before =
            self.quad_valence_deviation(valence_v0, valence_v1, valence_v2, valence_v3, false);
        let valence_after =
            self.quad_valence_deviation(valence_v0, valence_v1, valence_v2, valence_v3, true);

        // Only flip if it IMPROVES (reduces deviation)
        if valence_after < valence_before {
            return Some((
                c0,
                c1,
                valence_before - valence_after,
                self.corner_table.canonical_edge_hash(c0),
            ));
        }
        None // Don't flip if equal or worse
    }

    pub(crate) fn quad_valence_deviation(
        &self,
        v0: i16,
        v1: i16,
        v2: i16,
        v3: i16,
        after_flip: bool,
    ) -> i16 {
        let target_valence = 6; // For interior vertices

        let (val0, val1, val2, val3) = if after_flip {
            // After flip: V0 and V2 lose the edge between them
            //             V1 and V3 gain the new edge between them
            (v0 - 1, v1 + 1, v2 - 1, v3 + 1)
        } else {
            (v0, v1, v2, v3)
        };

        // Sum of squared deviations from ideal valence
        (val0 - target_valence).pow(2)
            + (val1 - target_valence).pow(2)
            + (val2 - target_valence).pow(2)
            + (val3 - target_valence).pow(2)
    }
}