ConformalLabpp/code/include/CGAL/Discrete_conformal_map.h

// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT
//
// Package: conformallab++ / Discrete_conformal_map (Phase 8 MVP, 2026-05-19)

/*!
\file CGAL/Discrete_conformal_map.h
\ingroup PkgConformalMapRef

User-facing entry point for the Discrete_conformal_map package.

This header provides a single function — `discrete_conformal_map_euclidean`
— that computes a Euclidean discrete-conformal flattening of an open or
closed triangle mesh.  Spherical and hyperbolic variants are scheduled
for Phase 8b.2 once the Euclidean pattern is validated by Phase 9a
(Inversive-Distance functional).

\section Example Simplest usage

\code{.cpp}
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Discrete_conformal_map.h>

using K = CGAL::Simple_cartesian<double>;
using Mesh = CGAL::Surface_mesh<K::Point_3>;

int main() {
    Mesh mesh = ...;  // load a triangle mesh
    auto result = CGAL::discrete_conformal_map_euclidean(mesh);
    if (!result.converged)
        return 1;
    // result.u_per_vertex[v] now holds the conformal scale factor at v.
}
\endcode

\section NamedParams Tuning via named parameters

\code{.cpp}
auto result = CGAL::discrete_conformal_map_euclidean(
    mesh,
    CGAL::parameters::gradient_tolerance(1e-12)
                     .max_iterations(500));
\endcode

\sa `CGAL::Default_conformal_map_traits`
\sa `CGAL::parameters::vertex_curvature_map`
*/

#ifndef CGAL_DISCRETE_CONFORMAL_MAP_H
#define CGAL_DISCRETE_CONFORMAL_MAP_H

#include <CGAL/Conformal_map_traits.h>
#include <CGAL/Conformal_map/internal/parameters.h>
#include <CGAL/Kernel_traits.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/property_map.h>

// Existing implementation headers (Layer 1 — unchanged).
#include "../euclidean_functional.hpp"
#include "../layout.hpp"
#include "../spherical_functional.hpp"
#include "../hyper_ideal_functional.hpp"
#include "../gauss_bonnet.hpp"
#include "../newton_solver.hpp"

#include <vector>
#include <unordered_map>

namespace CGAL {

// ════════════════════════════════════════════════════════════════════════════
// Result type
// ════════════════════════════════════════════════════════════════════════════

/*!
\ingroup PkgConformalMapRef

Result of `discrete_conformal_map_euclidean`.  Carries the converged
scale factors `u_v`, Newton diagnostics, and the convergence flag.
*/
template <typename FT = double>
struct Conformal_map_result
{
    /// Conformal scale factor `u_v` per vertex (indexed by raw vertex index).
    /// Length: `num_vertices(mesh)`.
    std::vector<FT> u_per_vertex;

    /// Number of Newton iterations performed.
    int iterations = 0;

    /// `‖G(u*)‖∞` at termination.
    FT gradient_norm = FT(0);

    /// `true` iff `gradient_norm < gradient_tolerance`.
    bool converged = false;

    /// `true` iff the linear solver used the SparseQR fallback at any
    /// Newton step (gauge mode on closed mesh without pinned vertex).
    bool sparse_qr_fallback_used = false;
};


// ════════════════════════════════════════════════════════════════════════════
// discrete_conformal_map_euclidean  —  user-facing entry
// ════════════════════════════════════════════════════════════════════════════

/*!
\ingroup PkgConformalMapRef

Compute the Euclidean discrete-conformal map of `mesh`.

This is the user-facing entry of Phase 8 MVP.  Internally it delegates
to the existing implementation in `code/include/euclidean_functional.hpp`
and `code/include/newton_solver.hpp` (Phase 1–7), so the algorithmic
behaviour is identical to the legacy API; this function only changes
the public façade.

\tparam TriangleMesh  A `CGAL::Surface_mesh<P>` for some point type `P`.
                      Other `FaceGraph` models are planned for Phase 8a.2.
\tparam NamedParameters  Optional CGAL named-parameter pack.

\param mesh  The input mesh (modified in place: property maps are attached).
\param np    Named parameters:
  \cgalParamNBegin{vertex_curvature_map}
    \cgalParamDescription{Property map `vertex → FT` of target cone angles Θᵥ.}
    \cgalParamDefault{2π at interior vertices, π at boundary vertices.}
  \cgalParamNEnd
  \cgalParamNBegin{gradient_tolerance}
    \cgalParamDescription{Newton stops when `‖G(u)‖∞ < tol`.}
    \cgalParamDefault{`1e-10`}
  \cgalParamNEnd
  \cgalParamNBegin{max_iterations}
    \cgalParamDescription{Hard limit on Newton steps.}
    \cgalParamDefault{`200`}
  \cgalParamNEnd
  \cgalParamNBegin{fixed_vertex_map}
    \cgalParamDescription{Property map `vertex → bool`; `true` ⇒ pinned.}
    \cgalParamDefault{The first vertex in `mesh.vertices()` is pinned.}
  \cgalParamNEnd

\returns A `Conformal_map_result<FT>` carrying `u_v` and Newton diagnostics.

\pre `mesh` is a triangle mesh.
\pre `mesh` satisfies the Gauss–Bonnet relation
     `Σ(2π − Θᵥ) = 2π·χ(mesh)` for the chosen target curvature map.
*/
template <typename TriangleMesh,
          typename CGAL_NP_TEMPLATE_PARAMETERS>
auto discrete_conformal_map_euclidean(
    TriangleMesh& mesh,
    const CGAL_NP_CLASS& np = parameters::default_values())
{
    // ── Type plumbing ──────────────────────────────────────────────────────
    //
    // Deduce the kernel from the mesh's Point_3 type rather than hard-coding
    // Simple_cartesian<double>.  This lets the wrapper work with any
    // Surface_mesh<P> whose P is a CGAL kernel point.  The user can override
    // the entire traits class via the `geom_traits(...)` named parameter
    // (Phase 8b.2 extension; default below covers the common case).
    using Point_type      = typename TriangleMesh::Point;
    using Default_kernel  = typename CGAL::Kernel_traits<Point_type>::Kernel;
    using Default_traits  = Default_conformal_map_traits<TriangleMesh, Default_kernel>;
    using Traits          = typename internal_np::Lookup_named_param_def<
                                internal_np::geom_traits_t,
                                CGAL_NP_CLASS,
                                Default_traits>::type;
    using FT              = typename Traits::FT;

    Conformal_map_result<FT> result;

    // ── 1. Set up property maps (legacy layer) ─────────────────────────────
    auto maps = ::conformallab::setup_euclidean_maps(mesh);
    ::conformallab::compute_euclidean_lambda0_from_mesh(mesh, maps);

    // ── 2. Target curvature: user-supplied or "natural-theta" default ─────
    auto theta_param = parameters::get_parameter(
        np, Conformal_map::internal_np::vertex_curvature_map);
    constexpr bool has_theta = !std::is_same_v<
        decltype(theta_param), internal_np::Param_not_found>;
    if constexpr (has_theta) {
        // User-provided Θ: copy into the property map and verify Gauss–Bonnet.
        // Throws std::runtime_error if the user-supplied Θ violates GB.
        for (auto v : mesh.vertices())
            maps.theta_v[v] = get(theta_param, v);
        ::conformallab::check_gauss_bonnet(mesh, maps);
    }
    // If no Θ is supplied, the default behaviour is "natural-theta": set Θ
    // such that x = 0 is the natural equilibrium (the actual angle sums at
    // x = 0 become the targets).  This matches the convention of the
    // existing test suite and guarantees that the default invocation
    // converges immediately for any well-formed triangle mesh.
    // The actual Θ adjustment is done after DOF assignment (step 5b below).

    // ── 3. Pin vertices: user map, or the first vertex by default ──────────
    //
    // The legacy v_idx property map has -1 as default (= pinned).  We must
    // first mark every vertex as "free" (any non-negative sentinel), then
    // pin the requested ones, then assign sequential DOF indices.
    constexpr int FREE = 0;
    for (auto v : mesh.vertices())
        maps.v_idx[v] = FREE;

    auto pin_param = parameters::get_parameter(
        np, Conformal_map::internal_np::fixed_vertex_map);
    constexpr bool has_pin = !std::is_same_v<
        decltype(pin_param), internal_np::Param_not_found>;

    bool any_pinned = false;
    if constexpr (has_pin) {
        for (auto v : mesh.vertices())
            if (get(pin_param, v)) {
                maps.v_idx[v] = -1;
                any_pinned = true;
            }
    }
    if (!any_pinned) {
        auto it = mesh.vertices().begin();
        if (it != mesh.vertices().end()) {
            maps.v_idx[*it] = -1;
            any_pinned = true;
        }
    }

    // ── 4. Assign DOF indices 0..n−1 to non-pinned vertices ─────────────────
    int idx = 0;
    for (auto v : mesh.vertices())
        if (maps.v_idx[v] != -1)
            maps.v_idx[v] = idx++;

    // ── 5. Read tolerances ─────────────────────────────────────────────────
    const FT  tol      = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance),
        FT(1e-10));
    const int max_iter = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::max_iterations),
        200);

    // ── 5b. Natural-theta default: shift Θ so that x = 0 is the equilibrium
    //
    // Only applied when the user did NOT supply a vertex_curvature_map.
    // The trick: evaluate G at x = 0, then subtract G_v from Θ_v.  After
    // this shift the new G(0) is identically zero, so Newton starts at the
    // optimum and immediately reports "converged".  This matches the
    // contract of the existing test suite ("natural-theta" pattern).
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    if constexpr (!has_theta) {
        auto G0 = ::conformallab::euclidean_gradient(mesh, x0, maps);
        for (auto v : mesh.vertices()) {
            const int j = maps.v_idx[v];
            if (j >= 0)
                maps.theta_v[v] -= G0[static_cast<std::size_t>(j)];
        }
    }

    // ── 6. Newton on x_0 = 0 ───────────────────────────────────────────────
    auto nr = ::conformallab::newton_euclidean(mesh, x0, maps, tol, max_iter);

    // ── 7. Pack result: u_v for every vertex, including pinned (u=0) ──────
    result.u_per_vertex.assign(num_vertices(mesh), FT(0));
    for (auto v : mesh.vertices()) {
        const int j = maps.v_idx[v];
        if (j >= 0)
            result.u_per_vertex[v.idx()] = nr.x[static_cast<std::size_t>(j)];
        // else: pinned ⇒ u_v stays 0
    }
    result.iterations    = nr.iterations;
    result.gradient_norm = nr.grad_inf_norm;
    result.converged     = nr.converged;

    // ── 8. Optional layout step (Phase 8b-Lite extension) ──────────────────
    //
    // If the caller supplied `output_uv_map(pmap)`, run the priority-BFS
    // trilateration on the converged x and write per-vertex `Point_2`
    // coordinates into `pmap`.  Optional `normalise_layout(true)` applies
    // the canonical PCA centroid + major-axis normalisation.
    auto uv_param = parameters::get_parameter(
        np, Conformal_map::internal_np::output_uv_map);
    constexpr bool has_uv = !std::is_same_v<
        decltype(uv_param), internal_np::Param_not_found>;
    if constexpr (has_uv) {
        if (nr.converged) {
            auto layout = ::conformallab::euclidean_layout(mesh, nr.x, maps);

            const bool do_norm = parameters::choose_parameter(
                parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout),
                false);
            if (do_norm) ::conformallab::normalise_euclidean(layout);

            for (auto v : mesh.vertices()) {
                const auto& uv = layout.uv[v.idx()];
                put(uv_param, v,
                    typename Traits::Kernel::Point_2(uv.x(), uv.y()));
            }
        }
    }

    return result;
}


// ════════════════════════════════════════════════════════════════════════════
// discrete_conformal_map_spherical  —  Phase 8b-Lite
// ════════════════════════════════════════════════════════════════════════════

/*!
\ingroup PkgConformalMapRef

Compute the spherical discrete-conformal map of a closed genus-0 mesh.

The spherical DCE energy is *concave*, so its Hessian is NSD at the
optimum and `newton_spherical()` factorises −H internally (handled by
the legacy implementation; no caller action required).  A gauge vertex
is pinned automatically to remove the rotational mode.

\tparam TriangleMesh  A `CGAL::Surface_mesh<P>` for some point type `P`.
\tparam NamedParameters  Optional CGAL named-parameter pack.

\param mesh  The input mesh (modified in place: property maps attached).
\param np    Same named parameters as `discrete_conformal_map_euclidean`.

\returns A `Conformal_map_result<FT>` carrying `u_v` per vertex and
         Newton diagnostics.

\pre `mesh` is a closed genus-0 triangle mesh.
\pre The user-supplied or natural-theta Θ satisfies the spherical
     Gauss–Bonnet relation `Σ(2π − Θᵥ) = 4π` (sphere).
*/
template <typename TriangleMesh,
          typename CGAL_NP_TEMPLATE_PARAMETERS>
auto discrete_conformal_map_spherical(
    TriangleMesh& mesh,
    const CGAL_NP_CLASS& np = parameters::default_values())
{
    using Point_type     = typename TriangleMesh::Point;
    using Default_kernel = typename CGAL::Kernel_traits<Point_type>::Kernel;
    using Default_traits = Default_conformal_map_traits<TriangleMesh, Default_kernel>;
    using Traits         = typename internal_np::Lookup_named_param_def<
                              internal_np::geom_traits_t,
                              CGAL_NP_CLASS,
                              Default_traits>::type;
    using FT             = typename Traits::FT;

    Conformal_map_result<FT> result;

    auto maps = ::conformallab::setup_spherical_maps(mesh);
    ::conformallab::compute_lambda0_from_mesh(mesh, maps);

    auto theta_param = parameters::get_parameter(
        np, Conformal_map::internal_np::vertex_curvature_map);
    constexpr bool has_theta = !std::is_same_v<
        decltype(theta_param), internal_np::Param_not_found>;
    if constexpr (has_theta) {
        for (auto v : mesh.vertices())
            maps.theta_v[v] = get(theta_param, v);
    }

    // Pin one vertex (gauge fix) — user-supplied or first vertex.
    constexpr int FREE = 0;
    for (auto v : mesh.vertices()) maps.v_idx[v] = FREE;

    auto pin_param = parameters::get_parameter(
        np, Conformal_map::internal_np::fixed_vertex_map);
    constexpr bool has_pin = !std::is_same_v<
        decltype(pin_param), internal_np::Param_not_found>;

    bool any_pinned = false;
    if constexpr (has_pin) {
        for (auto v : mesh.vertices())
            if (get(pin_param, v)) { maps.v_idx[v] = -1; any_pinned = true; }
    }
    if (!any_pinned) {
        auto it = mesh.vertices().begin();
        if (it != mesh.vertices().end()) { maps.v_idx[*it] = -1; any_pinned = true; }
    }

    int idx = 0;
    for (auto v : mesh.vertices())
        if (maps.v_idx[v] != -1) maps.v_idx[v] = idx++;

    const FT  tol      = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance),
        FT(1e-10));
    const int max_iter = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::max_iterations),
        200);

    // Natural-theta default for the spherical functional.
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    if constexpr (!has_theta) {
        auto G0 = ::conformallab::spherical_gradient(mesh, x0, maps);
        for (auto v : mesh.vertices()) {
            const int j = maps.v_idx[v];
            if (j >= 0) maps.theta_v[v] -= G0[static_cast<std::size_t>(j)];
        }
    }

    auto nr = ::conformallab::newton_spherical(mesh, x0, maps, tol, max_iter);

    result.u_per_vertex.assign(num_vertices(mesh), FT(0));
    for (auto v : mesh.vertices()) {
        const int j = maps.v_idx[v];
        if (j >= 0) result.u_per_vertex[v.idx()] = nr.x[static_cast<std::size_t>(j)];
    }
    result.iterations    = nr.iterations;
    result.gradient_norm = nr.grad_inf_norm;
    result.converged     = nr.converged;

    // Optional 3-D layout step (point on S²)
    auto uv_param = parameters::get_parameter(
        np, Conformal_map::internal_np::output_uv_map);
    constexpr bool has_uv = !std::is_same_v<
        decltype(uv_param), internal_np::Param_not_found>;
    if constexpr (has_uv) {
        if (nr.converged) {
            auto layout = ::conformallab::spherical_layout(mesh, nr.x, maps);
            const bool do_norm = parameters::choose_parameter(
                parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout),
                false);
            if (do_norm) ::conformallab::normalise_spherical(layout);
            for (auto v : mesh.vertices()) {
                const auto& p = layout.pos[v.idx()];
                put(uv_param, v,
                    typename Traits::Kernel::Point_3(p.x(), p.y(), p.z()));
            }
        }
    }

    return result;
}


// ════════════════════════════════════════════════════════════════════════════
// discrete_conformal_map_hyper_ideal  —  Phase 8b-Lite
// ════════════════════════════════════════════════════════════════════════════

/*!
\ingroup PkgConformalMapRef

Result of `discrete_conformal_map_hyper_ideal`.  Carries both vertex
DOFs `b_v` and edge DOFs `a_e` (hyper-ideal triangles in H³).
*/
template <typename FT = double>
struct Hyper_ideal_map_result
{
    /// Vertex DOFs `b_v` (length = num_vertices(mesh); pinned vertices = 0).
    std::vector<FT> b_per_vertex;
    /// Edge DOFs `a_e` (length = num_edges(mesh); pinned edges = 0).
    std::vector<FT> a_per_edge;

    /// Newton iterations actually performed (≤ `max_iterations`).
    int  iterations    = 0;
    /// Final infinity-norm of the gradient (Newton stopping criterion).
    FT   gradient_norm = FT(0);
    /// `true` iff `gradient_norm < gradient_tolerance` at exit.
    bool converged     = false;
};

/*!
\ingroup PkgConformalMapRef

Compute the hyper-ideal discrete-conformal map of a triangle mesh
(Springborn 2020 §4).

\note  Phase 8b-Lite scope: vertex DOFs `b_v` are assigned automatically
       to all vertices; edge DOFs `a_e` are similarly assigned.  The
       block-FD Hessian (Phase 9b) is used internally — see
       `newton_hyper_ideal` for the solver convention.
*/
template <typename TriangleMesh,
          typename CGAL_NP_TEMPLATE_PARAMETERS>
auto discrete_conformal_map_hyper_ideal(
    TriangleMesh& mesh,
    const CGAL_NP_CLASS& np = parameters::default_values())
{
    using Point_type     = typename TriangleMesh::Point;
    using Default_kernel = typename CGAL::Kernel_traits<Point_type>::Kernel;
    using Default_traits = Default_conformal_map_traits<TriangleMesh, Default_kernel>;
    using Traits         = typename internal_np::Lookup_named_param_def<
                              internal_np::geom_traits_t,
                              CGAL_NP_CLASS,
                              Default_traits>::type;
    using FT             = typename Traits::FT;

    Hyper_ideal_map_result<FT> result;

    auto maps = ::conformallab::setup_hyper_ideal_maps(mesh);
    // Hyper-ideal init does not derive from mesh geometry: the user's
    // Θ_v and θ_e are the model inputs.  Defaults from setup are
    // Θ_v = 2π, θ_e = π (orthogonal).

    auto theta_param = parameters::get_parameter(
        np, Conformal_map::internal_np::vertex_curvature_map);
    constexpr bool has_theta = !std::is_same_v<
        decltype(theta_param), internal_np::Param_not_found>;
    if constexpr (has_theta) {
        for (auto v : mesh.vertices())
            maps.theta_v[v] = get(theta_param, v);
    }

    const int n = ::conformallab::assign_all_dof_indices(mesh, maps);

    const FT  tol      = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::gradient_tolerance),
        FT(1e-8));
    const int max_iter = parameters::choose_parameter(
        parameters::get_parameter(np, Conformal_map::internal_np::max_iterations),
        200);

    // Initial point: b_v = 1.0 (positive log-scale), a_e = 0.5 (moderate).
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    for (auto v : mesh.vertices()) {
        int i = maps.v_idx[v];
        if (i >= 0) x0[static_cast<std::size_t>(i)] = 1.0;
    }
    for (auto e : mesh.edges()) {
        int i = maps.e_idx[e];
        if (i >= 0) x0[static_cast<std::size_t>(i)] = 0.5;
    }

    auto nr = ::conformallab::newton_hyper_ideal(mesh, x0, maps, tol, max_iter);

    result.b_per_vertex.assign(num_vertices(mesh), FT(0));
    result.a_per_edge   .assign(num_edges(mesh),    FT(0));
    for (auto v : mesh.vertices()) {
        int j = maps.v_idx[v];
        if (j >= 0) result.b_per_vertex[v.idx()] = nr.x[static_cast<std::size_t>(j)];
    }
    for (auto e : mesh.edges()) {
        int j = maps.e_idx[e];
        if (j >= 0) result.a_per_edge[e.idx()] = nr.x[static_cast<std::size_t>(j)];
    }
    result.iterations    = nr.iterations;
    result.gradient_norm = nr.grad_inf_norm;
    result.converged     = nr.converged;

    // Optional Poincaré-disk layout (2-D in the unit disk).
    auto uv_param = parameters::get_parameter(
        np, Conformal_map::internal_np::output_uv_map);
    constexpr bool has_uv = !std::is_same_v<
        decltype(uv_param), internal_np::Param_not_found>;
    if constexpr (has_uv) {
        if (nr.converged) {
            auto layout = ::conformallab::hyper_ideal_layout(mesh, nr.x, maps);
            const bool do_norm = parameters::choose_parameter(
                parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout),
                false);
            if (do_norm) ::conformallab::normalise_hyperbolic(layout);
            for (auto v : mesh.vertices()) {
                const auto& uv = layout.uv[v.idx()];
                put(uv_param, v,
                    typename Traits::Kernel::Point_2(uv.x(), uv.y()));
            }
        }
    }

    (void)n;   // already-counted by maps; silence unused-var warnings if any
    return result;
}

}  // namespace CGAL

#endif  // CGAL_DISCRETE_CONFORMAL_MAP_H