// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// conformal_quality.hpp
//
// Phase 9g.1 — Quantitative correctness metrics for computed conformal maps.
//
// Measures the quality and validity of a discrete conformal map layout:
// - IsothermicityMeasure: pointwise deviation from conformality (metric anisotropy).
// - DiscreteConformalEquivalenceMeasure: per-edge length-cross-ratio residual.
// - FlippedTriangles: detects inverted/degenerate triangles in 2-D layouts.
// - LengthCrossRatio: the discrete conformal invariant (per-edge).
// - ConvergenceUtility: aggregated convergence measures (max, mean, sum of cross-ratios).
//
// Mathematical references:
//   Springborn-Schröder-Pinkall 2008: discrete conformal invariant theory.
//   Bobenko-Springborn 2004: variational foundation.
//
// Java sources (ported from):
//   plugin/visualizer/IsothermicityMeasure.java
//   plugin/visualizer/DiscreteConformalEquivalencemMeasure.java
//   plugin/visualizer/FlippedTriangles.java
//   heds/adapter/types/LengthCrossRatio.java
//   convergence/ConvergenceUtility.java

#pragma once

#include "conformal_mesh.hpp"
#include "layout.hpp"
#include <Eigen/Dense>
#include <vector>
#include <cmath>
#include <algorithm>

namespace conformallab {

// ────────────────────────────────────────────────────────────────────────────
// LengthCrossRatio — the discrete conformal invariant
// ────────────────────────────────────────────────────────────────────────────

/// Compute the cross-ratio q = (a·c)/(b·d) of the four edges of a
/// quadrilateral formed by two adjacent triangles sharing an edge.
/// Input: edge lengths a, b, c, d in order around the quad.
/// Returns the cross-ratio q.
inline double length_cross_ratio(double a, double b, double c, double d)
{
    const double denom = b * d;
    if (denom < 1e-16) return 0.0;  // degenerate edge
    return (a * c) / denom;
}

// ────────────────────────────────────────────────────────────────────────────
// IsothermicityMeasure — pointwise metric anisotropy
// ────────────────────────────────────────────────────────────────────────────

/// Evaluate the isothermicity measure at a single vertex in a 2-D layout.
/// Isothermicity is the local conformality condition: the metric tensor
/// is a positive scalar multiple of the identity (no anisotropy).
/// Measure: pointwise deviation from a conformal map.
/// Returns the anisotropy ratio (1.0 = isotropic / conformal).
inline double isothermicity_measure_at_vertex(
    const ConformalMesh&          mesh,
    Vertex_index                  v,
    const Layout2D&               layout)
{
    // Collect all halfedges emanating from v.
    std::vector<Halfedge_index> hs;
    for (auto h : CGAL::halfedges_around_source(v, mesh))
        hs.push_back(h);

    if (hs.empty()) return 1.0;

    // Compute metric tensor components at v via edge pairs.
    // For a conformal map, the metric g = λ²I (λ > 0 scale factor, I identity).
    // Compute an empirical metric from the layout: edges adjacent to v
    // span the tangent space.
    double g11 = 0.0, g12 = 0.0, g22 = 0.0;
    int n_edges = 0;

    for (std::size_t i = 0; i < hs.size(); ++i) {
        auto h1 = hs[i];
        auto h2 = hs[(i + 1) % hs.size()];

        Vertex_index v2 = mesh.target(h1);  // = mesh.source(h2)
        Vertex_index v3 = mesh.target(h2);

        const auto& p1 = layout.uv[v.idx()];
        const auto& p2 = layout.uv[v2.idx()];
        const auto& p3 = layout.uv[v3.idx()];

        // Two edge vectors from v.
        double e1x = p2.x() - p1.x(), e1y = p2.y() - p1.y();
        double e2x = p3.x() - p1.x(), e2y = p3.y() - p1.y();

        // Metric tensor as outer product (unnormalised).
        g11 += e1x * e1x;
        g12 += e1x * e1y;
        g22 += e1y * e1y;

        // Also accumulate e2 contribution (for a rotationally averaged metric).
        g11 += e2x * e2x;
        g12 += e2x * e2y;
        g22 += e2y * e2y;

        n_edges += 2;
    }

    if (n_edges <= 0) return 1.0;

    g11 /= n_edges;
    g12 /= n_edges;
    g22 /= n_edges;

    // Eigenvalues of g: λ_± = (g11 + g22 ± √((g11-g22)² + 4g12²)) / 2.
    double trace = g11 + g22;
    double det = g11 * g22 - g12 * g12;

    if (trace < 1e-16 || det < 1e-16) return 1.0;  // degenerate

    double disc = (g11 - g22) * (g11 - g22) + 4.0 * g12 * g12;
    disc = std::sqrt(disc);

    double lambda_max = (trace + disc) / 2.0;
    double lambda_min = (trace - disc) / 2.0;

    if (lambda_min < 1e-16) return 1.0;  // degenerate

    // Anisotropy: λ_max / λ_min (conformal ⟺ ratio ≈ 1).
    return lambda_max / lambda_min;
}

/// Compute the isothermicity measure for the entire layout.
/// Returns a vector of anisotropy ratios, one per vertex.
inline std::vector<double> isothermicity_measure(
    const ConformalMesh&          mesh,
    const Layout2D&               layout)
{
    std::vector<double> result;
    result.reserve(mesh.number_of_vertices());
    for (auto v : mesh.vertices())
        result.push_back(isothermicity_measure_at_vertex(mesh, v, layout));
    return result;
}

// ────────────────────────────────────────────────────────────────────────────
// DiscreteConformalEquivalenceMeasure — length-cross-ratio residual
// ────────────────────────────────────────────────────────────────────────────

/// Evaluate the discrete conformal equivalence condition at a single edge.
/// For an edge e = (i,j), form the quad with the two adjacent triangles:
/// compute the cross-ratio q from the layout edge lengths.
/// The conformal condition is: q + 1/q = 2 (i.e. q = 1, isotropic scaling).
/// Measure: |q + 1/q - 2| (residual; 0 = conformal).
inline double discrete_conformal_equivalence_at_edge(
    const ConformalMesh&          mesh,
    Edge_index                    e,
    const Layout2D&               layout)
{
    // Find the two halfedges for this edge.
    auto h = mesh.halfedge(e);

    // Get the four vertices of the quad formed by the two adjacent triangles.
    Vertex_index v1 = mesh.source(h);
    Vertex_index v2 = mesh.target(h);
    Vertex_index v3 = mesh.source(mesh.next(h));
    Vertex_index v4 = mesh.source(mesh.next(mesh.opposite(h)));

    // Compute edge lengths from the layout.
    auto dist = [&layout](Vertex_index u1, Vertex_index u2) {
        const auto& p1 = layout.uv[u1.idx()];
        const auto& p2 = layout.uv[u2.idx()];
        double dx = p1.x() - p2.x();
        double dy = p1.y() - p2.y();
        return std::sqrt(dx * dx + dy * dy);
    };

    double a = dist(v1, v3);  // opposite to v4
    double b = dist(v1, v4);  // opposite to v3
    double c = dist(v2, v3);  // opposite to v4
    double d = dist(v2, v4);  // opposite to v3

    // Cross-ratio q = (a·c)/(b·d).
    double q = length_cross_ratio(a, b, c, d);

    // Conformal condition: q + 1/q = 2 (only satisfied when q = 1).
    if (q < 1e-16) return 1.0;  // degenerate
    double residual = q + 1.0 / q - 2.0;
    return std::abs(residual);
}

/// Compute the discrete conformal equivalence measure for all edges.
/// Returns a vector of residuals, one per edge.
inline std::vector<double> discrete_conformal_equivalence_measure(
    const ConformalMesh&          mesh,
    const Layout2D&               layout)
{
    std::vector<double> result;
    result.reserve(mesh.number_of_edges());
    for (auto e : mesh.edges())
        result.push_back(discrete_conformal_equivalence_at_edge(mesh, e, layout));
    return result;
}

// ────────────────────────────────────────────────────────────────────────────
// FlippedTriangles — embedded validity check
// ────────────────────────────────────────────────────────────────────────────

/// Check if a single triangle is flipped or degenerate in the 2-D layout.
/// A triangle is valid iff its signed area > 0 (positive orientation).
/// Degenerate: signed area ≈ 0 (collinear or nearly collinear vertices).
/// Returns true if the triangle is flipped or degenerate.
inline bool is_flipped_triangle(
    const ConformalMesh&          mesh,
    Face_index                    f,
    const Layout2D&               layout)
{
    // Extract the three vertices of the triangle.
    auto h = mesh.halfedge(f);
    Vertex_index v1 = mesh.source(h);
    Vertex_index v2 = mesh.source(mesh.next(h));
    Vertex_index v3 = mesh.source(mesh.next(mesh.next(h)));

    const auto& p1 = layout.uv[v1.idx()];
    const auto& p2 = layout.uv[v2.idx()];
    const auto& p3 = layout.uv[v3.idx()];

    // Signed area (× 2): (p2 - p1) × (p3 - p1) in ℝ².
    double signed_area_2x = (p2.x() - p1.x()) * (p3.y() - p1.y())
                          - (p2.y() - p1.y()) * (p3.x() - p1.x());

    // Positive area: valid orientation. Zero or negative: flipped/degenerate.
    return signed_area_2x <= 1e-14;
}

/// Count the number of flipped or degenerate triangles in the layout.
/// Returns the count (0 = valid layout).
inline int flipped_triangles(
    const ConformalMesh&          mesh,
    const Layout2D&               layout)
{
    int count = 0;
    for (auto f : mesh.faces())
        if (is_flipped_triangle(mesh, f, layout))
            count++;
    return count;
}

// ────────────────────────────────────────────────────────────────────────────
// ConvergenceUtility — aggregated convergence measures
// ────────────────────────────────────────────────────────────────────────────

/// Aggregated cross-ratio statistics for a layout.
struct CrossRatioStats {
    double max_cross_ratio;      ///< max of (q + 1/q) over all edges
    double mean_cross_ratio;     ///< mean of (q + 1/q)
    double sum_cross_ratio;      ///< sum of (q + 1/q)

    double max_multi_ratio;      ///< max per-face product of cross-ratios
    double mean_multi_ratio;     ///< mean per-face product
    double sum_multi_ratio;      ///< sum of per-face products

    double max_scale_invariant_circumradius;  ///< max of R/√A per face
    double mean_scale_invariant_circumradius; ///< mean of R/√A
    double sum_scale_invariant_circumradius;  ///< sum of R/√A
};

/// Compute convergence statistics for a layout.
/// - Cross-ratio (q + 1/q) per edge; aggregated max/mean/sum.
/// - Multi-ratio: per-face product ∏(q + 1/q) for the 3 edges of each face.
///   (Multi-ratio = 1 iff all edges are conformal.)
/// - Scale-invariant circumradius: R/√A per face (mesh quality metric).
inline CrossRatioStats convergence_utility(
    const ConformalMesh&          mesh,
    const Layout2D&               layout)
{
    CrossRatioStats stats = {};

    std::vector<double> cross_ratios;
    std::vector<double> multi_ratios;
    std::vector<double> scale_inv_circumradii;

    auto dist = [&layout](Vertex_index u1, Vertex_index u2) {
        const auto& p1 = layout.uv[u1.idx()];
        const auto& p2 = layout.uv[u2.idx()];
        double dx = p1.x() - p2.x();
        double dy = p1.y() - p2.y();
        return std::sqrt(dx * dx + dy * dy);
    };

    // Per-face metrics.
    for (auto f : mesh.faces()) {
        auto h = mesh.halfedge(f);
        Vertex_index v1 = mesh.source(h);
        Vertex_index v2 = mesh.source(mesh.next(h));
        Vertex_index v3 = mesh.source(mesh.next(mesh.next(h)));

        const auto& p1 = layout.uv[v1.idx()];
        const auto& p2 = layout.uv[v2.idx()];
        const auto& p3 = layout.uv[v3.idx()];

        // Signed area.
        double signed_area_2x = (p2.x() - p1.x()) * (p3.y() - p1.y())
                              - (p2.y() - p1.y()) * (p3.x() - p1.x());
        double area = std::abs(signed_area_2x) / 2.0;

        if (area < 1e-16) continue;  // degenerate

        // Three edge lengths of the triangle.
        double a = dist(v1, v2);
        double b = dist(v2, v3);
        double c = dist(v3, v1);

        // Circumradius R = abc / (4·Area).
        double circum_radius = (a * b * c) / (4.0 * area);

        // Scale-invariant: R / √A.
        double scale_inv_cr = circum_radius / std::sqrt(area);
        scale_inv_circumradii.push_back(scale_inv_cr);

        // Three cross-ratios (per edge/angle of the triangle).
        // For each edge, form the quad with the opposite vertex and its neighbors.
        double multi_product = 1.0;
        for (int ei = 0; ei < 3; ++ei) {
            auto he = mesh.halfedge(f);
            for (int k = 0; k < ei; ++k) he = mesh.next(he);

            Vertex_index eu1 = mesh.source(he);
            Vertex_index eu2 = mesh.target(he);
            Vertex_index eu3 = mesh.source(mesh.next(he));
            Vertex_index eu4 = mesh.source(mesh.next(mesh.opposite(he)));

            double ea = dist(eu1, eu3);
            double eb = dist(eu1, eu4);
            double ec = dist(eu2, eu3);
            double ed = dist(eu2, eu4);

            double q = length_cross_ratio(ea, eb, ec, ed);
            if (q > 1e-16) {
                double qf = q + 1.0 / q;
                cross_ratios.push_back(qf);
                multi_product *= qf;
            }
        }
        multi_ratios.push_back(multi_product);
    }

    // Aggregate statistics.
    if (!cross_ratios.empty()) {
        auto [min_it, max_it] = std::minmax_element(cross_ratios.begin(), cross_ratios.end());
        stats.max_cross_ratio = *max_it;
        stats.mean_cross_ratio = 0.0;
        for (double v : cross_ratios) stats.mean_cross_ratio += v;
        stats.mean_cross_ratio /= static_cast<double>(cross_ratios.size());
        stats.sum_cross_ratio = 0.0;
        for (double v : cross_ratios) stats.sum_cross_ratio += v;
    }

    if (!multi_ratios.empty()) {
        auto [min_it, max_it] = std::minmax_element(multi_ratios.begin(), multi_ratios.end());
        stats.max_multi_ratio = *max_it;
        stats.mean_multi_ratio = 0.0;
        for (double v : multi_ratios) stats.mean_multi_ratio += v;
        stats.mean_multi_ratio /= static_cast<double>(multi_ratios.size());
        stats.sum_multi_ratio = 0.0;
        for (double v : multi_ratios) stats.sum_multi_ratio += v;
    }

    if (!scale_inv_circumradii.empty()) {
        auto [min_it, max_it] = std::minmax_element(scale_inv_circumradii.begin(),
                                                     scale_inv_circumradii.end());
        stats.max_scale_invariant_circumradius = *max_it;
        stats.mean_scale_invariant_circumradius = 0.0;
        for (double v : scale_inv_circumradii)
            stats.mean_scale_invariant_circumradius += v;
        stats.mean_scale_invariant_circumradius /= static_cast<double>(scale_inv_circumradii.size());
        stats.sum_scale_invariant_circumradius = 0.0;
        for (double v : scale_inv_circumradii)
            stats.sum_scale_invariant_circumradius += v;
    }

    return stats;
}

}  // namespace conformallab