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::macros::integrity_assert;
use crate::isotropic_remesh::IsotropicRemeshAlgo;
use std::fmt::Debug;
use vector_traits::num_traits::AsPrimitive;
use vector_traits::prelude::GenericVector3;

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,
{
    #[cfg(any(feature = "integrity_check", debug_assertions))]
    #[inline(always)]
    /// 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 normal_similarity_opposite_winding_slow_wrong_order(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
    ) -> S {
        Self::normal_similarity_opposite_winding_slow(c, a, b, d)
    }

    #[cfg(any(feature = "integrity_check", debug_assertions))]
    #[inline(always)]
    /// ```text
    ///           C
    ///          /│\
    ///        /  │  \
    ///      /    │    \
    ///    /      │      \
    ///  A ────── B ────── D
    /// ```
    /// This function returns the cosine of the dihedral angle between the two triangles
    /// (via the dot product of their normals).
    /// Triangles ABC and DCB share edge B-C  traversed in the *opposite* direction
    pub(crate) fn normal_similarity_opposite_winding_slow(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
    ) -> S {
        // Calculate normals of both triangles
        // Triangle ABC: vertices A, B, C in CCW order
        let normal_abc = Self::triangle_raw_normal(a, b, c);

        // Triangle DCB: vertices D, C, B in CCW order (sharing edge B-C with ABC)
        let normal_dcb = Self::triangle_raw_normal(d, c, b);
        assert!(normal_abc.magnitude_sq() > S::EPSILON_SQ);

        normal_abc.normalize().dot(normal_dcb.normalize())
    }

    #[cfg(any(feature = "integrity_check", debug_assertions))]
    #[inline(always)]
    /// ```text
    ///   D ───── C
    ///    \     /│
    ///      \ /  │
    ///      / \  │
    ///    /     \│
    ///   A ───── B
    /// ```
    /// This function returns the cosine of the dihedral angle between the two triangles
    /// (via the dot product of their normals).
    /// Triangles ABC and DBC share edge B-C traversed in the *same* direction
    pub(crate) fn normal_similarity_same_winding_slow(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
    ) -> S {
        -Self::normal_similarity_opposite_winding_slow(a, b, c, d)
    }

    #[inline(always)]
    /// Helper that computes the scaled dot product for triangle normal comparison. ABC vs DCB
    /// Always computes as if same winding.
    /// Returns `None` if triangle ABC or BCD is degenerate, otherwise `Some((dot, len_product_sq))`
    fn compute_scaled_normal_dot(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
    ) -> Option<(S, S)> {
        let bc = b - c;
        let ac = a - c;

        // normal(ABC) computed as normal(CAB) = (a-c).cross(b-c)
        let normal_abc = ac.cross(bc);
        let len_sq_abc = normal_abc.magnitude_sq();

        if !Self::validate_triangle_normal(len_sq_abc) {
            return None; // ABC degenerate
        }

        // normal(DBC) computed as normal(CDB) = (d-c).cross(b-c)
        // Reuses bc from above!
        let dc = d - c;
        let normal_dbc = dc.cross(bc);
        let len_sq_dbc = normal_dbc.magnitude_sq();
        if !Self::validate_triangle_normal(len_sq_dbc) {
            return None; // BCD degenerate
        }

        let dot = normal_abc.dot(normal_dbc);
        let len_product_sq = len_sq_abc * len_sq_dbc;

        Some((dot, len_product_sq))
    }

    #[inline(always)]
    /// ```text
    ///           C
    ///          /│\
    ///        /  │  \
    ///      /    │    \
    ///    /      │      \
    ///  A ────── B ────── D
    /// ```
    /// This function returns the result of a cosine of the dihedral angle test between the two
    /// triangles ABC & DCB (using the dot product of their normals).
    /// Triangles ABC and DCB share vertices B and C
    /// The DCB triangle is not checked for zero area, only ABC (only A is actually moving)
    ///
    /// Returns:
    /// - `Some(true)`: triangles are compatible (not inverted), collapse is safe
    /// - `Some(false)`: triangle ABC is degenerate (zero area).
    /// - `None`: triangles have opposite winding (normals point in strongly opposite directions),
    pub(crate) fn normal_similarity_opposite_winding(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
        cos_threshold_sq: S,
    ) -> Option<bool> {
        debug_assert!(cos_threshold_sq <= S::ONE);
        debug_assert!(cos_threshold_sq >= S::ZERO);

        let (dot, len_product_sq) = match Self::compute_scaled_normal_dot(a, b, c, d) {
            None => return Some(false), // degenerate
            Some((d, l)) => (-d, l),    // negate dot for opposite winding
        };

        if dot < S::ZERO {
            return None;
        }

        #[cfg(any(feature = "integrity_check", debug_assertions))]
        let slow_dihedral = Self::normal_similarity_opposite_winding_slow(a, b, c, d);
        if dot * dot >= cos_threshold_sq * len_product_sq {
            integrity_assert!(slow_dihedral >= cos_threshold_sq.sqrt());
            Some(true)
        } else {
            integrity_assert!(slow_dihedral < cos_threshold_sq.sqrt());
            None
        }
    }

    #[inline(always)]
    /// ```text
    ///   D ───── C
    ///    \     /│
    ///      \ /  │
    ///      / \  │
    ///    /     \│
    ///   A ───── B
    /// ```
    /// tests if the squared dot product between triangle normal of ABC vs the normal of triangle DBC
    /// are larger or equal than `cos_normal_sq_threshold`
    /// Returns:
    /// - `Some(true)`: triangles are OK (normals similar enough)
    /// - `Some(false)`: resulting triangle ABC would be degenerate (zero area)
    /// - `None`: triangles rejected (normals too different and/or opposite)
    pub(crate) fn normal_similarity_same_winding(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
        cos_threshold_sq: S,
    ) -> Option<bool> {
        #[cfg(any(feature = "integrity_check", debug_assertions))]
        use vector_traits::num_traits::Float;

        debug_assert!(cos_threshold_sq <= S::ONE);
        debug_assert!(cos_threshold_sq >= S::ZERO);

        let (dot, len_product_sq) = match Self::compute_scaled_normal_dot(a, b, c, d) {
            None => return Some(false), // degenerate
            Some(vals) => vals,         // use as-is for same winding
        };

        if dot < S::ZERO {
            return None;
        }
        #[cfg(any(feature = "integrity_check", debug_assertions))]
        let slow_dihedral = Self::normal_similarity_same_winding_slow(a, b, c, d);

        if dot * dot >= cos_threshold_sq * len_product_sq {
            #[cfg(any(feature = "integrity_check", debug_assertions))]
            if Float::is_finite(slow_dihedral) {
                integrity_assert!(
                    slow_dihedral >= cos_threshold_sq.sqrt(),
                    "{slow_dihedral} !>= {}",
                    cos_threshold_sq.sqrt()
                );
            }
            Some(true)
        } else {
            #[cfg(any(feature = "integrity_check", debug_assertions))]
            if Float::is_finite(slow_dihedral) {
                integrity_assert!(slow_dihedral < cos_threshold_sq.sqrt());
            }
            None
        }
    }

    #[inline(always)]
    /// ```text
    ///           C
    ///          /│\
    ///        /  │  \
    ///      /    │    \
    ///    /      │      \
    ///  A ────── B ────── D
    /// ```
    /// This function checks if two triangles ABC & DCB have normals that are NOT too opposite.
    /// Triangles ABC and DCB share vertices B and C (in opposite direction).
    ///
    /// The dihedral angle is measured between the two triangle normals. This function checks
    /// if the normals are pointing in dangerously opposite directions (indicating a non-convex
    /// or folded configuration).
    ///
    /// # Angle Interpretation
    /// - 0° = normals perfectly aligned (coplanar, same direction)
    /// - -90° = normals perpendicular
    /// - -180° = normals perfectly opposite (completely folded/inverted)
    ///
    /// # Threshold Parameter
    /// it is expected that the `cos_threshold` parameter, that `cos_threshold_sq` derives from, is negative.
    /// `cos_threshold_sq` is the **squared** cosine of the maximum allowed angle.
    /// - For angle = -154°: cos(-154°) ≈ -0.9, so cos_threshold_sq = 0.81
    /// - For angle = -165°: cos(-165°) ≈ -0.966, so cos_threshold_sq = 0.933
    ///
    /// A larger `cos_threshold_sq` value means stricter checking (rejecting angles closer to -180°).
    ///
    /// # Returns
    /// - `Some(true)`: normals are NOT too opposite (angle < threshold), mesh is "ok".
    /// - `Some(false)`: triangle ABC is degenerate (zero area)
    /// - `None`: normals are too opposite (angle >= threshold), indicating dangerous folding
    pub(crate) fn normal_similarity_crease_angle(
        a: S::Vec3Simd,
        b: S::Vec3Simd,
        c: S::Vec3Simd,
        d: S::Vec3Simd,
        cos_max_crease_angle_sq: S,
    ) -> Option<bool> {
        debug_assert!(cos_max_crease_angle_sq <= S::ONE);
        debug_assert!(cos_max_crease_angle_sq >= S::ZERO);

        let (dot, len_product_sq) = match Self::compute_scaled_normal_dot(a, b, c, d) {
            // degenerate
            None => return Some(false),
            // negate dot for opposite winding
            Some((d, l)) => (-d, l),
        };

        // For negative threshold check, dot product must be negative
        if dot >= S::ZERO {
            return Some(true); // Not opposite, passes the test
        }

        // Check if angle is LESS opposite than threshold (closer to perpendicular)
        // dot is negative, we want: |dot| < cos_threshold * sqrt(len_product_sq)
        // Which means: dot² < cos_threshold² * len_product²
        if dot * dot < cos_max_crease_angle_sq * len_product_sq {
            Some(true) // Less opposite than threshold, passes
        } else {
            None // Too opposite, fails
        }
    }

    /// Validates two consecutive triangles in a fan using pre-computed normals.
    ///
    /// This is the specialized version of `normal_similarity_crease_angle` for vertex fans,
    /// where triangles share the center vertex and consecutive edge vertices.
    ///
    /// # Arguments
    /// * `normal_abc` - Pre-computed normal of the first triangle (unnormalized)
    /// * `len_sq_abc` - Pre-computed squared length of `normal_abc`
    /// * `normal_bcd` - Pre-computed normal of the second triangle (unnormalized)
    /// * `len_sq_bcd` - Pre-computed squared length of `normal_bcd`
    /// * `cos_max_crease_angle_sq` - Squared cosine of maximum allowed crease angle
    ///
    /// # Returns
    /// - `Some(true)`: normals are NOT too opposite (angle < threshold), mesh is "ok"
    /// - `Some(false)`: one of the triangles is degenerate (zero area)
    /// - `None`: normals are too opposite (angle >= threshold), indicating dangerous folding
    ///
    /// # Note
    /// Unlike `normal_similarity_crease_angle`, this method expects normals computed with
    /// SAME winding (both normals point in the same general direction for coplanar triangles).
    /// Therefore, we check for NEGATIVE dot products to detect opposite-facing normals.
    pub(crate) fn validate_fan_triangle_pair(
        normal_abc: S::Vec3Simd,
        len_sq_abc: S,
        normal_bcd: S::Vec3Simd,
        len_sq_bcd: S,
        cos_max_crease_angle_sq: S,
    ) -> Option<bool> {
        debug_assert!(cos_max_crease_angle_sq <= S::ONE);
        debug_assert!(cos_max_crease_angle_sq >= S::ZERO);

        // Check for degenerate triangles
        if !Self::validate_triangle_normal(len_sq_abc) {
            return Some(false);
        }
        if !Self::validate_triangle_normal(len_sq_bcd) {
            return Some(false);
        }

        let dot = normal_abc.dot(normal_bcd);

        // For fan triangles with same winding, dot product should be positive
        // if they're facing roughly the same direction
        if dot >= S::ZERO {
            return Some(true); // Not opposite, passes the test
        }

        // Check if angle is LESS opposite than threshold
        // dot is negative, we want: |dot| < cos_threshold * sqrt(len_product_sq)
        // Which means: dot² < cos_threshold² * len_product²
        let len_product_sq = len_sq_abc * len_sq_bcd;

        if dot * dot < cos_max_crease_angle_sq * len_product_sq {
            Some(true) // Less opposite than threshold, passes
        } else {
            None // Too opposite, fails
        }
    }

    #[inline(always)]
    /// Validates that a single triangle normal is not degenerate.
    ///
    /// # Returns
    /// - `true`: triangle has sufficient area (not degenerate)
    /// - `false`: triangle is degenerate (zero or near-zero area)
    pub(crate) fn validate_triangle_normal(len_sq: S) -> bool {
        len_sq >= S::EPSILON_SQ
    }

    /// Determine if an edge could be flipped by testing:
    /// ```text
    ///   1. V₀V₁V₂ is coplanar with V₀V₂V₃
    ///   2. V₀V₁V₂V₃ is strictly convex (V₂ and V₀ on opposite sides of V₁V₃)
    ///      with sufficient margin to ensure positive area after flip
    ///
    ///     V₃───── V₂     V₃ ─── V₂
    ///     │      /│      │ \    │
    ///     │    /  │  =>  │  \   │
    ///     │  /    │      │   \  │
    ///     │/      │      │    \ │
    ///     V₀ ──── V₁     V₀ ─── V₁
    ///
    /// We know V₀V₁V₂ and V₀V₂V₃ has area > 0
    /// Everything is CCW
    /// ```
    pub(crate) fn quad_is_convex_coplanar(
        v0: S::Vec3Simd,
        v1: S::Vec3Simd,
        v2: S::Vec3Simd,
        v3: S::Vec3Simd,
        cos_threshold_sq: S,
    ) -> bool {
        let e01 = v1 - v0;
        let e02 = v2 - v0;
        let e03 = v3 - v0;

        // Triangles along diagonal v0-v2
        let n_012 = e01.cross(e02); // Triangle (v0, v1, v2)
        let n_023 = e02.cross(e03); // Triangle (v0, v2, v3)

        // Coplanarity: same-side normals should align
        let dot1 = n_012.dot(n_023);
        if dot1 < S::ZERO {
            return false;
        }

        let len_sq_012 = n_012.magnitude_sq();
        let len_sq_023 = n_023.magnitude_sq();

        if dot1 * dot1 < cos_threshold_sq * len_sq_012 * len_sq_023 {
            return false;
        }

        // Triangles along diagonal v1-v3
        let e12 = e02 - e01; // v2 - v1
        let e13 = e03 - e01; // v3 - v1
        let n_123 = e12.cross(e13); // Triangle (v1, v2, v3)
        let neg_n_130 = e13.cross(e01); // Triangle (v1, v3, v0) - inverted sign because we use -e10==e01 instead of e10

        // Convexity: opposite-side normals should point opposite directions
        // Now expecting negative dot since we removed the negation
        let dot2 = n_123.dot(neg_n_130);
        if dot2 >= S::ZERO {
            return false;
        }

        let len_sq_123 = n_123.magnitude_sq();
        let len_sq_130 = neg_n_130.magnitude_sq();

        // Check for degenerate triangles after flip
        if len_sq_123 < S::EPSILON_SQ || len_sq_130 < S::EPSILON_SQ {
            return false;
        }

        // Strict convexity: ensure |dot2|² > threshold² * |n_123|² * |n_130|²
        // Since dot2 < 0, we check: dot2² > threshold² * len_sq_123 * len_sq_130
        dot2 * dot2 > cos_threshold_sq * len_sq_123 * len_sq_130
    }

    /*
    #[allow(dead_code)]
    /// Validate consecutive triangles in the fan using a rolling window optimization.
    ///
    /// This method efficiently checks consecutive triangle pairs by computing each
    /// triangle normal only once and reusing it for the dihedral angle check.
    ///
    /// Important: This method requires the leading vertex of each triangle as the vertex fan
    /// is rotated CCW. This is the reason by `skip_leading` must  be >= 1.
    /// I.e. `get_vertex(corner_index.prev())` for a CCW mesh.
    ///
    /// # Arguments
    /// * `fan_vertices` an iterator over the vertices of the vertex fan. Must be `V(corner.prev)`
    /// * `center` - The center vertex position of the fan (V₀)
    /// * `skip_leading` - Number of triangles to skip at the beginning of the fan
    /// * `skip_trailing` - Number of triangles to skip at the end of the fan
    /// * `validator` - Closure that receives two consecutive triangle normals and their
    ///   squared lengths, returns true if the dihedral angle is acceptable
    ///
    /// # Returns
    /// `true` if all checked triangle pairs pass validation, `false` otherwise
    pub(crate) fn validate_consecutive_triangles_itertools<
        I: ExactSizeIterator<Item = S::Vec3Simd>,
        F: FnMut(S::Vec3Simd, S, S::Vec3Simd, S) -> bool,
    >(
        fan_vertices: I,
        center: S::Vec3Simd,
        skip_leading: usize,
        skip_trailing: usize,
        mut validator: F,
    ) -> bool {
        // Note: This implementation does NOT handle wraparound validation
        // (checking the pair between last and first triangle). It's designed
        // for validating a contiguous window of triangles, not the full fan.
        // Some leading triangles must be skipped.
        debug_assert!(
            skip_leading > 0,
            "Full fan validation with wraparound not supported"
        );

        let total_vertices = fan_vertices.len();
        let total_triangles = total_vertices;
        let triangles_to_check = total_triangles
            .saturating_sub(skip_leading)
            .saturating_sub(skip_trailing);
        let pairs_to_check = triangles_to_check.saturating_sub(1);
        let vertices_needed = pairs_to_check + 2;

        // We need vertices starting from (skip_leading - 1) because:
        // - Triangle at index skip_leading uses vertices [skip_leading-1, skip_leading]
        let start_vertex = skip_leading.saturating_sub(1);
        let end_vertex = (start_vertex + vertices_needed).min(total_vertices);

        println!(
            "validate_consecutive_triangles2 vertices:[{:?}..{:?}], triangles:[{:?}..{:?}]",
            start_vertex,
            end_vertex,
            skip_leading,
            total_triangles.saturating_sub(skip_trailing)
        );

        if end_vertex <= start_vertex + 1 {
            integrity_println!(
                "validate_consecutive_triangles2: nothing to compare. {} <= {}",
                end_vertex,
                start_vertex + 1
            );
            return true;
        }

        fan_vertices
            .skip(start_vertex)
            .take(end_vertex - start_vertex)
            .inspect(|&v| println!("#got vertex:{:?}", v.into()))
            // vertices -> edges from center
            .map(|v| v - center)
            .inspect(|&e| println!("#created edge:{:?}", e.into()))
            // edges -> window(2) -> triangle normals with lengths
            .tuple_windows::<(_, _)>()
            .map(|(e1, e2)| {
                let normal = e2.cross(e1);
                println!(
                    "#created normal:{:?} from {:?} and {:?}",
                    normal.into(),
                    e1.into(),
                    e2.into()
                );
                (normal, normal.magnitude_sq())
            })
            // triangle normals -> window(2) -> validate pairs
            .tuple_windows::<(_, _)>()
            .all(|((prev_n, prev_len), (curr_n, curr_len))| {
                println!(
                    "#comparing normal {:?} & {:?}",
                    prev_n.into(),
                    curr_n.into()
                );
                validator(prev_n, prev_len, curr_n, curr_len)
            })
    }*/

    /// Validate consecutive triangles in the fan using a rolling window optimization.
    ///
    /// This method efficiently checks consecutive triangle pairs by computing each
    /// triangle normal only once and reusing it for the dihedral angle check.
    ///
    /// Important: This method requires the leading vertex of each triangle as the vertex fan
    /// is rotated CCW. This is the reason by `skip_leading` must  be >= 1.
    /// I.e. `get_vertex(corner_index.prev())` for a CCW mesh.
    ///
    /// # Arguments
    /// * `leading_fan_vertices` an iterator over the leading vertices of the vertex fan. Must be `V(corner.prev())`
    /// * `center` - The center vertex position of the fan (V₀)
    /// * `skip_leading` - Number of triangles to skip at the beginning of the fan
    /// * `skip_trailing` - Number of triangles to skip at the end of the fan
    /// * `validator` - Closure that receives two consecutive triangle normals and their
    ///   squared lengths, returns true if the dihedral angle is acceptable
    ///
    /// # Returns
    /// `true` if all checked triangle pairs pass validation, `false` otherwise
    pub(crate) fn validate_consecutive_triangles<
        I: ExactSizeIterator<Item = S::Vec3Simd>,
        F: FnMut(S::Vec3Simd, S, S::Vec3Simd, S) -> bool,
    >(
        leading_fan_vertices: I,
        center: S::Vec3Simd,
        skip_leading: usize,
        skip_trailing: usize,
        mut validator: F,
    ) -> bool {
        debug_assert!(
            skip_leading > 0,
            "Full fan validation with wraparound not supported"
        );

        let total_vertices = leading_fan_vertices.len();
        let total_triangles = total_vertices;

        let triangles_to_check = total_triangles
            .saturating_sub(skip_leading)
            .saturating_sub(skip_trailing);
        let pairs_to_check = triangles_to_check.saturating_sub(1);
        let vertices_needed = pairs_to_check + 2;

        let start_vertex = skip_leading.saturating_sub(1);
        let end_vertex = (start_vertex + vertices_needed).min(total_vertices);

        if end_vertex <= start_vertex + 2 {
            return true;
        }

        let mut vertices = leading_fan_vertices
            .skip(start_vertex)
            .take(end_vertex - start_vertex);

        // Get first three vertices to bootstrap
        let v0 = match vertices.next() {
            Some(v) => v - center,
            None => return true,
        };
        let v1 = match vertices.next() {
            Some(v) => v - center,
            None => return true,
        };
        let v2 = match vertices.next() {
            Some(v) => v - center,
            None => return true,
        };

        // Compute first two triangle normals
        let mut prev_normal = v1.cross(v0);
        let mut prev_len = prev_normal.magnitude_sq();

        let mut curr_normal = v2.cross(v1);
        let mut curr_len = curr_normal.magnitude_sq();

        // Check first pair
        if !validator(prev_normal, prev_len, curr_normal, curr_len) {
            return false;
        }

        // Rolling window: for each new vertex, compute new triangle and check pair
        let mut prev_edge = v2;
        for vertex in vertices {
            let edge = vertex - center;

            // Shift window
            prev_normal = curr_normal;
            prev_len = curr_len;

            // Compute new triangle
            curr_normal = edge.cross(prev_edge);
            curr_len = curr_normal.magnitude_sq();

            // Validate pair
            if !validator(prev_normal, prev_len, curr_normal, curr_len) {
                return false;
            }

            prev_edge = edge;
        }

        true
    }
}