ConformalLabpp/code/include/CGAL/Discrete_circle_packing.h

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

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

User-facing entry for the **face-based** circle-packing functional of
Bobenko-Pinkall-Springborn 2015.  See `cp_euclidean_functional.hpp`
for the underlying algorithm and `doc/architecture/phase-9a-validation.md`
for the line-by-line mapping to the Java original
`CPEuclideanFunctional.java`.

This functional has a fundamentally different DOF structure to the
classical Euclidean / Spherical / HyperIdeal modes — one log-radius
`ρ_f` per **face** rather than one log-scale `u_v` per vertex.  We
therefore expose it via a dedicated header with its own default-trait
class (Strategy C of the Phase 8b architecture audit).
*/

#ifndef CGAL_DISCRETE_CIRCLE_PACKING_H
#define CGAL_DISCRETE_CIRCLE_PACKING_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/Surface_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <boost/graph/graph_traits.hpp>

#include "../cp_euclidean_functional.hpp"
#include "../newton_solver.hpp"

#include <stdexcept>

namespace CGAL {

// ── Default traits for CP-Euclidean ───────────────────────────────────────────

/*!
\ingroup PkgConformalMapConcepts
\brief Traits class for `discrete_circle_packing_euclidean()` —
declares the kernel, mesh and property-map types used by the
BPS-2015 face-based circle-packing functional.

Primary template; specialise it for non-`Surface_mesh` triangle meshes.
*/
template <typename TriangleMesh,
          typename Kernel_ = CGAL::Simple_cartesian<double>>
struct Default_cp_euclidean_traits;

/*!
\ingroup PkgConformalMapConcepts
\brief Specialisation for `CGAL::Surface_mesh<P>`; the only one shipped
in Phase 8b-Lite.
*/
template <typename K>
struct Default_cp_euclidean_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
{
    /// CGAL kernel parameter (defaults to `Simple_cartesian<double>`).
    using Kernel        = K;
    /// Scalar field type used for all CP-Euclidean DOFs (`ρ_f`, `θ_e`, `φ_f`).
    using FT            = typename K::FT;
    /// 3-D point type (vertex coordinates).
    using Point_3       = typename K::Point_3;
    /// Triangle-mesh type this specialisation targets.
    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;

    /// Boost-graph vertex descriptor for `Triangle_mesh`.
    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
    /// Boost-graph half-edge descriptor for `Triangle_mesh`.
    using Halfedge_descriptor = typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor;
    /// Boost-graph edge descriptor for `Triangle_mesh`.
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
    /// Boost-graph face descriptor for `Triangle_mesh`.
    using Face_descriptor     = typename boost::graph_traits<Triangle_mesh>::face_descriptor;

    // CP-Euclidean property maps — note the *face* DOF index map.

    /// Property map face → contiguous integer DOF index (legacy `cf:idx`).
    using Face_index_pmap = typename Triangle_mesh::template Property_map<Face_descriptor,   int>;
    /// Property map edge → intersection angle θₑ (legacy `ce:theta`).
    using Theta_e_pmap    = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
    /// Property map face → target angle sum φ_f (legacy `cf:phi`).
    using Phi_f_pmap      = typename Triangle_mesh::template Property_map<Face_descriptor,   FT>;
};

// ── Result type ───────────────────────────────────────────────────────────────

/*!
\ingroup PkgConformalMapRef

Result of `discrete_circle_packing_euclidean`.  Carries face DOFs
`ρ_f = log R_f` rather than the vertex DOFs of the classical modes.
*/
/// Why the Newton solve terminated (re-exported from the solver; I1 audit).
using Newton_status = ::conformallab::NewtonStatus;

template <typename FT = double>
struct Circle_packing_result
{
    /// Face DOFs `ρ_f = log R_f` (length = num_faces(mesh); pinned face = 0).
    std::vector<FT> rho_per_face;

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

    /// Why the solve stopped: `Converged` / `MaxIterations` /
    /// `LinearSolverFailed` / `LineSearchStalled`.
    Newton_status status = Newton_status::MaxIterations;
    /// `true` iff any Newton step fell back to SparseQR (gauge singularity).
    bool sparse_qr_fallback_used = false;
    /// Smallest `|Dᵢᵢ|` of the last LDLT (near-singularity proxy; 0 if none).
    FT   min_ldlt_pivot = FT(0);
};

// ── Entry function ────────────────────────────────────────────────────────────

/*!
\ingroup PkgConformalMapRef

Compute the BPS-2015 face-based circle-packing of `mesh`.

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

\param mesh  Input triangle mesh.
\param np    Named parameters (subset of those documented on
             `discrete_conformal_map_euclidean`; the curvature-map
             parameter `vertex_curvature_map` is **not** used in this
             face-based mode — instead the per-face target angle sum
             `φ_f` and per-edge intersection angle `θ_e` are set via
             the property maps on `mesh` before this call, or left at
             their defaults `φ_f = 2π`, `θ_e = π/2`).

\returns A `Circle_packing_result<FT>` with `ρ_f` per face.

\pre `mesh` is a triangle mesh.
\pre `φ_f` and `θ_e` satisfy the BPS-2015 admissibility conditions
     (Σ_f φ_f = 2π·χ + Σ_e (π − θ_e), see paper §6).
*/
template <typename TriangleMesh,
          typename CGAL_NP_TEMPLATE_PARAMETERS>
auto discrete_circle_packing_euclidean(
    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_cp_euclidean_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;

    Circle_packing_result<FT> result;

    auto maps = ::conformallab::setup_cp_euclidean_maps(mesh);

    // Pin first face by default; `fixed_vertex_map` is reused here as the
    // "fixed face" override hook (the parameter tag is generic enough).
    // For a richer API, a dedicated `fixed_face_map` tag could be added.
    auto it = mesh.faces().begin();
    if (it == mesh.faces().end()) {
        return result;   // empty mesh; trivial
    }
    const int n = ::conformallab::assign_cp_euclidean_face_dof_indices(mesh, maps, *it);

    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-phi default: shift φ_f so the gradient at ρ = 0 is zero.
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto G0 = ::conformallab::cp_euclidean_gradient(mesh, x0, maps);
    for (auto f : mesh.faces()) {
        int i = maps.f_idx[f];
        if (i >= 0) maps.phi_f[f] -= G0[static_cast<std::size_t>(i)];
    }

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

    result.rho_per_face.assign(num_faces(mesh), FT(0));
    for (auto f : mesh.faces()) {
        int j = maps.f_idx[f];
        if (j >= 0) result.rho_per_face[f.idx()] = nr.x[static_cast<std::size_t>(j)];
    }
    result.iterations    = nr.iterations;
    result.gradient_norm = nr.grad_inf_norm;
    result.converged     = nr.converged;
    result.status                  = nr.status;
    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);

    // ── output_uv_map (Phase 8b-Lite extension) ────────────────────────────
    //
    // The CP-Euclidean functional carries one DOF per *face* (the log of the
    // face-circle radius `ρ_f = log R_f`), not per vertex.  A faithful
    // layout therefore produces a circle packing in ℝ² — each face f is
    // mapped to a circle of radius `R_f` at some centre `c_f`, with
    // adjacent circles meeting at the prescribed intersection angle `θ_e`.
    // That is a per-face output, not the per-vertex Point_2 that
    // `output_uv_map` is typed for.
    //
    // For Phase 8b-Lite we deliberately don't fake it.  If the caller
    // supplies `output_uv_map(pmap)` we throw `std::runtime_error` with a
    // clear pointer to Phase 9c (BPS-2015 §6 face-based circle-packing
    // layout, ~150 lines, on the porting roadmap).  Failing loudly is
    // better than silently writing zeros.
    //
    // Users who want a UV-like coordinate today can:
    //   1. Solve a Euclidean DCE on the same mesh (vertex DOFs),
    //   2. Use `discrete_inversive_distance_map(... output_uv_map(pmap))`,
    //   3. Or compute face-centre positions by hand from `result.rho_per_face`
    //      + the per-edge `θ_e` values, plus a priority-BFS of their own.
    {
        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) {
            throw std::runtime_error(
                "CGAL::discrete_circle_packing_euclidean: the "
                "`output_uv_map(...)` named parameter is not yet supported "
                "for face-based CP-Euclidean.  The faithful output is a "
                "circle packing in the plane (per-face), not per-vertex "
                "UVs.  Tracked as Phase 9c; "
                "see doc/architecture/locked-vs-flexible.md and "
                "doc/tutorials/add-output-uv-map.md.");
        }
    }

    return result;
}

}  // namespace CGAL

#endif  // CGAL_DISCRETE_CIRCLE_PACKING_H