remesh 0.0.5

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

use crate::common::VertexIndex;
use crate::common::macros::{
    integrity_assert, integrity_assert_eq, integrity_assert_unique, integrity_println,
};
use crate::{
    corner_table::{CornerIndex, TriangleIndex},
    isotropic_remesh::IsotropicRemeshAlgo,
};
use std::fmt::Debug;
use vector_traits::num_traits::{AsPrimitive, Float};
use vector_traits::prelude::{GenericVector3, 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,
{
    const SUBSCRIPT_CHARS: [char; 10] = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'];

    pub(crate) fn subscript(n: usize) -> String {
        n.to_string()
            .chars()
            .map(|c| Self::SUBSCRIPT_CHARS[c.to_digit(10).unwrap() as usize])
            .collect()
    }

    #[inline(always)]
    pub(crate) fn triangle_raw_normal(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
    ) -> S::Vec3Simd {
        assert!(a.is_finite());
        assert!(b.is_finite());
        assert!(c.is_finite());
        (b - a).cross(c - a)
    }

    pub(crate) fn iter_valid_vertices(&self) -> impl Iterator<Item = VertexIndex> {
        (0..self.vertices.len() as u32)
            .map(VertexIndex)
            .filter(|&v| self.corner_table.corner(v).is_valid())
    }

    fn check_integrity(&self) -> Result<(), String> {
        self.corner_table
            .check_integrity(self.vertices.len() as u32)?;

        for v in (0..self.vertices.len() as u32).map(VertexIndex) {
            if !self.corner_table.corner(v).is_valid() {
                assert!(!self.vertex_pool.is_used(v));
            } else {
                assert!(self.vertex_pool.is_used(v));
            }
        }

        for v in self.iter_valid_vertices() {
            self.check_vertex_fan(v)?;
        }

        // Checks for triangles with opposing normals that could create zero-volume regions
        let triangle_count = self.corner_table.vertex_of_corners().len() as u32 / 3;

        for triangle in (0..triangle_count).map(TriangleIndex) {
            {
                let [c, cn, cp] = triangle.corners();
                if self.corner_table.is_triangle_deleted(triangle) {
                    self.corner_table.check_deleted_triangle(triangle);
                    continue;
                }

                let vic = self.corner_table.vertex(c);
                let vc = self.vertex_of_corner(c).to_simd();
                let vn = self.vertex_of_corner(cn).to_simd();
                let vin = self.corner_table.vertex(cn);
                let vp = self.vertex_of_corner(cp).to_simd();
                let vip = self.corner_table.vertex(cp);
                let l = Self::triangle_raw_normal(vc, vn, vp).magnitude();
                if (!Float::is_finite(l)) || Float::abs(l) < S::EPSILON {
                    return Err(format!(
                        "Zero area triangle: {c:?}{vic:?},{cn:?}{vin:?},{cp:?}{vip:?} normal length:{l}",
                    ));
                }
            }

            for i in 0..3 {
                // ```text
                //           vp
                //          / │
                //        /   │
                //       vc   │
                //     /   \  │
                //   /      \ │
                // vo ─────── vn
                // ```

                let c = CornerIndex(triangle.0 * 3 + i);
                let (cn, cp) = self.corner_table.next_prev(c); // next
                let vic = self.corner_table.vertex(c);
                let vin = self.corner_table.vertex(cn);
                let co = self.corner_table.opposite(cp);
                debug_assert!(co.is_valid());
                debug_assert!(vic.is_valid());
                debug_assert!(vin.is_valid());

                let vio = self.corner_table.vertex(co);
                debug_assert!(vio.is_valid());
                let vip = self.corner_table.vertex(cp);
                debug_assert!(vip.is_valid());
                integrity_assert_unique!(vic, vio, vin, vip);

                // Only process edge in one direction to avoid duplicates
                if vic.0 < vin.0 {
                    // tests abc:(vc,vn,vp) vs dba:(vo,vn,vc)
                    let coplanarity = self.coplanarity(c);
                    assert!(Float::is_finite(coplanarity));

                    // If normals are opposing (dot product close to -1), report it
                    if coplanarity < self.params.inversion_validation_threshold {
                        let co = self.corner_table.opposite(cp);
                        let (con, cop) = self.corner_table.next_prev(co);
                        return Err(format!(
                            "adjacent triangles with opposite normals: {},{},{} & {},{},{} cop:{coplanarity}",
                            self.corner_table.data.dbg_corner(c),
                            self.corner_table.data.dbg_corner(cn),
                            self.corner_table.data.dbg_corner(cp),
                            self.corner_table.data.dbg_corner(con),
                            self.corner_table.data.dbg_corner(cop),
                            self.corner_table.data.dbg_corner(co),
                        ));
                    }
                }
            }
        }
        Ok(())
    }

    /// Debug method to check mesh integrity at any point
    pub fn check_mesh_integrity(&self, context: &str) -> Result<(), String> {
        match self.check_integrity() {
            Ok(()) => {
                integrity_println!("✓ Mesh integrity OK at: {context}");
                Ok(())
            }
            Err(e) => {
                eprintln!("✗ Mesh integrity FAILED at {context}: {e}");
                self.print_debug_stats();
                Err(e)
            }
        }
    }

    /// Quick integrity check for performance-critical paths
    pub fn check_mesh_integrity_fast(&self, context: &str) -> Result<(), String> {
        match self
            .corner_table
            .check_basic_integrity(self.vertices.len() as u32)
        {
            Ok(()) => Ok(()),
            Err(e) => {
                eprintln!("✗ Basic mesh integrity FAILED at {context}: {e}");
                Err(e)
            }
        }
    }

    pub(crate) fn dbg_edge(&self, c: CornerIndex) -> String {
        let cn = self.corner_table.next(c);
        //let vic = self.corner_table.get_vertex_of_corner(c);
        //let vicn = self.corner_table.get_vertex_of_corner(cn);
        //let vc = self.get_vertex_of_corner(c).into();
        //let vnc = self.get_vertex_of_corner(cn).into();
        format!(
            "{}->{}",
            self.corner_table.data.dbg_corner(c),
            self.corner_table.data.dbg_corner(cn)
        )
    }

    pub(crate) fn dbg_corner(&self, c: CornerIndex) -> String {
        self.corner_table.data.dbg_corner(c)
    }

    pub(crate) fn dbg_vertex(&self, vi: VertexIndex) -> String {
        let v = self.vertex(vi);
        let c = self.corner_table.corner(vi);
        let fan = self.corner_table.ccw_vertex_fan(c);
        format!("{vi:?}{:?} fan:{:?}", v.into(), fan)
    }

    pub(crate) fn dbg_vertex_fan(&self, fan: &[CornerIndex]) -> String {
        let parts: Vec<String> = fan
            .iter()
            .enumerate()
            .map(|(i, &c)| {
                let (n, p) = self.corner_table.next_prev(c);
                let o = self.corner_table.opposite(c);
                format!(
                    "t{}({}(O:{}),{},{})",
                    Self::subscript(i),
                    self.corner_table.data.dbg_corner(c),
                    self.corner_table.data.dbg_corner(o),
                    self.corner_table.data.dbg_corner(n),
                    self.corner_table.data.dbg_corner(p)
                )
            })
            .collect();
        format!("{}: valence:{}", parts.join(","), fan.len())
    }

    pub(crate) fn dbg_corner_table(&self, count: usize) {
        let num_corners = self.corner_table.len() as u32;
        let num_triangles = num_corners / 3;

        // Print first few triangles for manual inspection
        println!(
            "First {count} triangles (of {num_triangles}): [vertex_id:corner_id,opposite vertex_id:corner_id+1,opposite vertex_id:corner_id+2, opposite]"
        );

        for triangle in (0..num_triangles)
            .map(TriangleIndex)
            .filter(|t| !self.corner_table.is_triangle_deleted(*t))
            .take(count)
        {
            let c0 = CornerIndex(triangle.0 * 3);
            let c1 = CornerIndex(triangle.0 * 3 + 1);
            let c2 = CornerIndex(triangle.0 * 3 + 2);

            if c2.0 < num_corners {
                let vi0 = self.corner_table.vertex(c0);
                let vi1 = self.corner_table.vertex(c1);
                let vi2 = self.corner_table.vertex(c2);

                if !self.corner_table.is_triangle_deleted(triangle) {
                    let v0 = self.vertex(vi0).to_simd();
                    let v1 = self.vertex(vi1).to_simd();
                    let v2 = self.vertex(vi2).to_simd();

                    let n = Self::triangle_raw_normal(v0, v1, v2);
                    let v0: [S; 3] = v0.into();
                    let v1: [S; 3] = v1.into();
                    let v2: [S; 3] = v2.into();
                    if n.magnitude() < S::EPSILON {
                        println!(
                            "  {triangle:?}: ({}, {}, {}) ({v0:?},{v1:?},{v2:?}) zero area triangle",
                            self.corner_table.data.dbg_corner(c0),
                            self.corner_table.data.dbg_corner(c0),
                            self.corner_table.data.dbg_corner(c0)
                        );
                    } else {
                        println!(
                            "  {triangle:?}: ({}, {}, {}) ({v0:?},{v1:?},{v2:?})",
                            self.corner_table.data.dbg_corner(c0),
                            self.corner_table.data.dbg_corner(c0),
                            self.corner_table.data.dbg_corner(c0)
                        );
                    }
                } else {
                    println!(
                        "  T{triangle:?}: (C{} C{} C{}) deleted",
                        triangle.0 * 3,
                        triangle.0 * 3 + 1,
                        triangle.0 * 3 + 2
                    );
                }
            }
        }
    }

    /// Print detailed corner table statistics. Warning: O(n^2) or more
    pub fn print_debug_stats(&self) {
        let num_vertices = self.vertices.len() as u32;
        let num_corners = self.corner_table.len() as u32;
        let num_triangles = num_corners / 3;

        println!("=== Corner Table Debug Stats ===");
        println!("Vertices: {num_vertices}, Triangles: {num_triangles}, Corners: {num_corners}");

        println!(
            " free_triangles: {}",
            self.corner_table.deleted_triangle_count()
        );

        // Count boundary edges
        let boundary_edges = self
            .corner_table
            .opposite_corners()
            .iter()
            .enumerate()
            .map(|(c, &o)| (CornerIndex(c as u32), o))
            .filter(|(c, o)| {
                !self.corner_table.is_triangle_deleted((*c).into())
                    && self.corner_table.is_triangle_deleted((*o).into())
            })
            .count();
        println!("Boundary edges: {boundary_edges}");

        // Count connected components (simplified check)
        let mut vertex_degrees = vec![0_u32; num_vertices as usize];
        for &vertex in self.corner_table.vertex_of_corners().iter() {
            if vertex.0 < num_vertices {
                vertex_degrees[vertex.0 as usize] += 1;
            }
        }

        let isolated_vertices = vertex_degrees
            .iter()
            .enumerate()
            .filter(|(v, degree)| !self.is_vertex_deleted(VertexIndex(*v as u32)) && **degree == 0)
            .count();
        println!("Isolated vertices: {isolated_vertices}");

        if false {
            self.dbg_corner_table(20);
            println!();

            print!("  ");
            for v in self
                .vertices
                .iter()
                .enumerate()
                .map(|(v, _)| VertexIndex(v as u32))
            {
                if self.corner_table.corner_of_vertices()[v.0 as usize].is_valid() {
                    print!(
                        "{v:?}:{} ",
                        self.corner_table
                            .data
                            .dbg_corner(self.corner_table.corner_of_vertices()[v.0 as usize])
                    );
                } else {
                    print!("{v:?}:deleted ");
                }
            }
            println!();
        }
        println!("==============================");
    }

    /// ```text
    ///           vo
    ///          / │
    ///        /   │
    ///       vp   │
    ///     /   \  │
    ///   /      \ │
    /// vc ─────── vn
    /// ```
    /// returns the dot product between the normals of the triangles [c,c.next(),c.prev()]
    /// and [c.prev(),c.next(),c.opposite()]
    pub(crate) fn coplanarity(&self, c: CornerIndex) -> S {
        let (cn, cp) = self.corner_table.next_prev(c);

        let vic = self.corner_table.vertex(c);
        integrity_assert!(vic.is_valid());

        let vicn = self.corner_table.vertex(cn);
        integrity_assert!(vicn.is_valid());

        let vicp = self.corner_table.vertex(cp);
        integrity_assert!(vicp.is_valid());

        let co = self.corner_table.opposite(c);
        let (con, cop) = self.corner_table.next_prev(co);
        integrity_assert_eq!(self.corner_table.vertex(con), vicp);
        integrity_assert_eq!(self.corner_table.vertex(cop), vicn);

        let vico = self.corner_table.vertex(co);
        integrity_assert!(vico.is_valid());
        integrity_assert_unique!(vic, vicn, vicp, vico);

        // Get actual vertex positions
        let vc = self.vertex(vic).to_simd();
        let vn = self.vertex(vicn).to_simd();
        let vp = self.vertex(vicp).to_simd();
        let vo = self.vertex(vico).to_simd();

        // Check coplanarity between the two triangles
        assert_eq!(
            Self::normal_similarity_crease_angle(
                vc,
                vn,
                vp,
                vo,
                self.params.crease_limit_threshold_sq
            ),
            Some(true),
            "{:?},{:?},{:?},{:?} were too opposite",
            vc.into(),
            vn.into(),
            vp.into(),
            vo.into()
        );
        Self::normal_similarity_opposite_winding_slow(vc, vn, vp, vo)
    }

    #[inline(always)]
    /// ```text
    ///           c
    ///          /│
    ///        /  │
    ///       a   │
    ///     /  \  │
    ///   /     \ │
    /// d ─────── b
    /// ```
    /// This function returns the cosine of the dihedral angle between the two triangles
    /// (via the dot product of their normals).
    /// Triangles abc and dba share vertices a and b
    pub(crate) fn dihedral_angle_cosine_check(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
        threshold: S,
    ) -> bool {
        Self::normal_similarity_opposite_winding_slow_wrong_order(a, b, c, d) >= threshold
    }

    /// make sure the vertex fan of v is ok
    pub(crate) fn check_vertex_fan(&self, v: VertexIndex) -> Result<(), String> {
        if !v.is_valid() {
            return Err(format!("Vertex {v:?} was not valid"));
        }
        if self.is_vertex_deleted(v) {
            return Err(format!("Vertex {v:?} was valid but deleted"));
        }
        let c0 = self.corner_table.corner(v);
        if !c0.is_valid() {
            return Err(format!("Corner of Vertex {v:?} was invalid"));
        }
        if self.corner_table.is_triangle_deleted(c0.into()) {
            return Err(format!("Corner {c0:?} was already deleted"));
        }
        let c0_fan = self.corner_table.ccw_vertex_fan(c0);

        for &c in &c0_fan {
            if self.corner_table.is_triangle_deleted(c.into()) {
                return Err(format!("Corner {c:?} was already deleted"));
            }

            let o = self.corner_table.opposite(c);
            if self.corner_table.is_triangle_deleted(o.into()) {
                return Err(format!("opposite corner {o:?} was already deleted"));
            }
            let (cn, cp) = self.corner_table.next_prev(c);

            let vio = self.corner_table.vertex(o);
            let vic = self.corner_table.vertex(c);
            let vicn = self.corner_table.vertex(cn);
            let vicp = self.corner_table.vertex(cp);

            integrity_assert_unique!(vio, vic, vicn, vicp);

            let vo = self.vertex_of_corner(o).to_simd();
            let vc = self.vertex_of_corner(c).to_simd();
            let vcn = self.vertex_of_corner(cn).to_simd();
            let vcp = self.vertex_of_corner(cp).to_simd();

            match Self::normal_similarity_crease_angle(
                vc,
                vcn,
                vcp,
                vo,
                self.params.inversion_validation_threshold_sq,
            ) {
                Some(true) => {}
                v => {
                    integrity_println!("ci triangle:{}", self.corner_table.dbg_triangle(c));
                    integrity_println!("c.fan.prev: {}", self.dbg_fan_prev(&c0_fan));
                    return Err(format!(
                        "fail: dihedral_angle_cosine_test_negative_threshold()->{:?} {:?},{:?},{:?},{:?} {:?},{:?},{:?},{:?} {}",
                        v,
                        vic,
                        vicn,
                        vicp,
                        vio,
                        vc.into(),
                        vcn.into(),
                        vcp.into(),
                        vo.into(),
                        self.params.inversion_validation_threshold
                    ));
                }
            }
        }
        Ok(())
    }

    /// make sure the vertex fan of v is ok
    pub(crate) fn check_vertex_fan_dbg(&self, v: VertexIndex) -> Result<(), String> {
        let c0 = self.corner_table.corner(v);
        if self.corner_table.is_triangle_deleted(c0.into()) {
            return Err(format!("Corner {c0:?} was already deleted"));
        }
        for c in self.corner_table.iter_ccw_swing(c0) {
            if self.corner_table.is_triangle_deleted(c.into()) {
                return Err(format!("Corner {c:?} was already deleted"));
            }

            let o = self.corner_table.opposite(c);
            if self.corner_table.is_triangle_deleted(o.into()) {
                return Err(format!("opposite corner {o:?} was already deleted"));
            }
            let (cn, cp) = self.corner_table.next_prev(c);

            let vio = self.corner_table.vertex(o);
            let vic = self.corner_table.vertex(c);
            let vicn = self.corner_table.vertex(cn);
            let vicp = self.corner_table.vertex(cp);

            integrity_assert_unique!(vio, vic, vicn, vicp);

            let vo = self.vertex_of_corner(o).to_simd();
            let vc = self.vertex_of_corner(c).to_simd();
            let vcn = self.vertex_of_corner(cn).to_simd();
            let vcp = self.vertex_of_corner(cp).to_simd();

            match Self::normal_similarity_crease_angle(
                vc,
                vcn,
                vcp,
                vo,
                self.params.inversion_validation_threshold_sq,
            ) {
                Some(true) => {}
                Some(false) => {
                    return Err(format!(
                        "fail: dihedral_angle_cosine_test_negative_threshold()->Some(false) {:?},{:?},{:?},{:?} {:?},{:?},{:?},{:?} {}",
                        vic,
                        vicn,
                        vicp,
                        vio,
                        vc.into(),
                        vcn.into(),
                        vcp.into(),
                        vo.into(),
                        self.params.inversion_validation_threshold,
                    ));
                }
                None => {
                    return Err(format!(
                        "fail:dihedral_angle_cosine_test_negative_threshold()->None {:?},{:?},{:?},{:?} {:?},{:?},{:?},{:?} {}",
                        vic,
                        vicn,
                        vicp,
                        vio,
                        vc.into(),
                        vcn.into(),
                        vcp.into(),
                        vo.into(),
                        self.params.inversion_validation_threshold,
                    ));
                }
            }
            /*integrity_println!(
                "passed dihedral_angle_cosine_test({:?},{:?},{:?},{:?})",
                vic,
                vicn,
                vicp,
                vio
            );*/
        }
        Ok(())
    }

    pub(crate) fn dbg_fan_prev(&self, v_fan: &[CornerIndex]) -> String {
        v_fan
            .iter()
            .enumerate()
            .map(|(i, c)| {
                let p = self.corner_table.prev(*c);
                let vertex = self.corner_table.vertex(p);
                format!("V{}:{:?}", Self::subscript(i), vertex)
            })
            .collect::<Vec<_>>()
            .join(",")
    }

    pub(crate) fn dbg_fan_next(&self, v_fan: &[CornerIndex]) -> String {
        v_fan
            .iter()
            .enumerate()
            .map(|(i, c)| {
                let n = self.corner_table.next(*c);
                let vertex = self.corner_table.vertex(n);
                format!("V{}:{:?}", Self::subscript(i), vertex)
            })
            .collect::<Vec<_>>()
            .join(",")
    }
}