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 "../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;

    return result;
}

}  // namespace CGAL

#endif  // CGAL_DISCRETE_CONFORMAL_MAP_H