// 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, 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,
{
    /// Flip edges to improve triangle quality
    /// An edge is flipped if:
    /// 1. The two triangles are nearly coplanar (configurable threshold)
    /// 2. The new edge would be shorter than the current edge
    /// 3. The flip would not create degenerate triangles
    pub(crate) fn flip_edges_aspect_ratio(&mut self, aspect_quality_weight: S) -> 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;
            // 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_aspect_ratio(corner, aspect_quality_weight)
                })
                .collect();

            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, _edge)) = candidates.pop() {
                if c1 != self.corner_table.twin(c0) {
                    // c1 is no longer twin of c0. Could this ever happen?
                    continue;
                }
                // recheck if we should flip this edge, but ignore the aspect ratio quality this time
                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 new edge would be shorter than the current edge (TODO: should this be here?)
    /// 3. The flip would not create degenerate triangles
    /// 4. The valence quality is improved
    /// 5. The triangle quality is improved (better angles/aspect ratios)
    /// ```
    fn should_flip_edge_aspect_ratio(
        &self,
        c0: CornerIndex,
        quality_threshold: S,
    ) -> Option<(CornerIndex, CornerIndex, S, 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);

        // Calculate current edge length (V₀V₂)
        let current_edge_length_sq = (pos_v2 - pos_v0).magnitude_sq();

        // Calculate potential new edge length (V₃V₁)
        let new_edge_length_sq = (pos_v1 - pos_v3).magnitude_sq();

        // todo: should we really do this?
        // Only flip if new edge is shorter
        if new_edge_length_sq >= current_edge_length_sq * 1.5_f64.as_() {
            return None;
        }

        // 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_improvement = {
            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);

            // Hard constraint: don't worsen valence
            if valence_after > valence_before {
                return None;
            }

            valence_before - valence_after
        };

        // 6. Check triangle quality improvement

        let quality_before =
            self.min_triangle_quality_in_quad(pos_v0, pos_v1, pos_v2, pos_v3, false);
        let quality_after = self.min_triangle_quality_in_quad(pos_v0, pos_v1, pos_v2, pos_v3, true);
        let quality_improvement = quality_after - quality_before;

        // If valence improves, prioritize by valence but use quality as tiebreaker
        if valence_improvement > 0 {
            // Score: valence_improvement * 1000 + quality_improvement
            // Since quality is typically 0-1 range, this keeps valence dominant

            let score = (valence_improvement as f64).as_() * 1000.0_f64.as_() + quality_improvement;
            if Float::is_finite(score) {
                return Some((c0, c1, score, self.corner_table.canonical_edge_hash(c0)));
            } else {
                return None;
            }
        }

        // Valence neutral - only flip if quality improves enough
        if quality_improvement > quality_before * (quality_threshold - S::ONE) {
            return Some((
                c0,
                c1,
                quality_improvement,
                self.corner_table.canonical_edge_hash(c0),
            ));
        }

        None
    }

    fn min_triangle_quality_in_quad(
        &self,
        p0: S::Vec3Simd,
        p1: S::Vec3Simd,
        p2: S::Vec3Simd,
        p3: S::Vec3Simd,
        after_flip: bool,
    ) -> S {
        let (tri1_quality, tri2_quality) = if after_flip {
            // Triangles after flip: (v0, v2, v3) and (v1, v3, v2)
            (
                self.triangle_quality(p0, p2, p3),
                self.triangle_quality(p1, p3, p2),
            )
        } else {
            // Current triangles: (v0, v1, v2) and (v0, v3, v1)
            (
                self.triangle_quality(p0, p1, p2),
                self.triangle_quality(p0, p3, p1),
            )
        };

        Float::min(tri1_quality, tri2_quality)
    }

    fn triangle_quality(&self, p0: S::Vec3Simd, p1: S::Vec3Simd, p2: S::Vec3Simd) -> S {
        let min_angle = self.min_angle(p0, p1, p2);
        min_angle / 60.0_f64.to_radians().as_() // 60° is perfect
    }

    fn min_angle(&self, p0: S::Vec3Simd, p1: S::Vec3Simd, p2: S::Vec3Simd) -> S {
        let v0 = p1 - p0;
        let v0_magnitude = v0.magnitude();
        let v1 = p2 - p1;
        let v1_magnitude = v1.magnitude();
        let v2 = p0 - p2;
        let v2_magnitude = v2.magnitude();

        let angle0: S = (v0.dot(-v2) / (v0_magnitude * v2_magnitude)).acos();
        let angle1: S = (v1.dot(-v0) / (v1_magnitude * v0_magnitude)).acos();
        let angle2: S = (v2.dot(-v1) / (v2_magnitude * v1_magnitude)).acos();

        Float::min(Float::min(angle0, angle1), angle2)
    }
}
