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 2010.  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"

namespace CGAL {

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

template <typename TriangleMesh,
          typename Kernel_ = CGAL::Simple_cartesian<double>>
struct Default_cp_euclidean_traits;

template <typename K>
struct Default_cp_euclidean_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
{
    using Kernel        = K;
    using FT            = typename K::FT;
    using Point_3       = typename K::Point_3;
    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;

    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
    using Halfedge_descriptor = typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor;
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
    using Face_descriptor     = typename boost::graph_traits<Triangle_mesh>::face_descriptor;

    // CP-Euclidean property maps — note the *face* DOF index map.
    using Face_index_pmap = typename Triangle_mesh::template Property_map<Face_descriptor,   int>;
    using Theta_e_pmap    = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
    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.
*/
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;

    int  iterations    = 0;
    FT   gradient_norm = FT(0);
    bool converged     = false;
};

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

/*!
\ingroup PkgConformalMapRef

Compute the BPS-2010 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-2010 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;
    return result;
}

}  // namespace CGAL

#endif  // CGAL_DISCRETE_CIRCLE_PACKING_H