ConformalLabpp/code/include/CGAL/Discrete_inversive_distance.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_inversive_distance.h
\ingroup PkgConformalMapRef

User-facing entry for the **vertex-based** inversive-distance circle-
packing functional of Luo (2004), with the Bowers-Stephenson (2004)
initialisation.  See `inversive_distance_functional.hpp` for the
underlying algorithm and `doc/roadmap/research-track.md` (item 9a.2)
for the research-track classification — this functional has **no Java
original** (verified empirically), it is from-the-literature research.

DOF structure
─────────────
* Per-vertex `u_i = log r_i` (compatible with the classical Euclidean
  trait).
* Per-edge constant `I_ij` computed once by Bowers-Stephenson from the
  input mesh geometry (handled internally by
  `compute_inversive_distance_init_from_mesh`).

Because the per-edge constant has a different meaning from the
Euclidean `λ°_e`, this entry has its own default-trait class
`Default_inversive_distance_traits`.
*/

#ifndef CGAL_DISCRETE_INVERSIVE_DISTANCE_H
#define CGAL_DISCRETE_INVERSIVE_DISTANCE_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 <CGAL/Discrete_conformal_map.h>   // for Conformal_map_result<FT>

#include "../inversive_distance_functional.hpp"
#include "../newton_solver.hpp"

namespace CGAL {

// ── Default traits for Inversive-Distance ────────────────────────────────────

/*!
\ingroup PkgConformalMapConcepts
\brief Traits class for `discrete_inversive_distance_map()` — declares
the kernel, mesh and property-map types used by Luo's 2004 vertex-based
inversive-distance circle packing.

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

/*!
\ingroup PkgConformalMapConcepts
\brief Specialisation for `CGAL::Surface_mesh<P>`; the only one shipped
in Phase 8b-Lite.
*/
template <typename K>
struct Default_inversive_distance_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 inversive-distance DOFs.
    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 edge descriptor for `Triangle_mesh`.
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;

    // Inversive-distance specific property maps.

    /// Property map vertex → contiguous integer DOF index (legacy `iv:idx`).
    using Vertex_index_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, int>;
    /// Property map vertex → target cone angle Θᵥ in radians (legacy `iv:theta`).
    using Theta_v_pmap      = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
    /// Property map vertex → initial radius r⁰ᵥ (legacy `iv:r0`).
    using R0_pmap           = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
    /// Property map edge → inversive distance Iᵢⱼ (legacy `ie:I`).
    using I_e_pmap          = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
};

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

/*!
\ingroup PkgConformalMapRef

Compute the Luo-2004 vertex-based inversive-distance circle packing of `mesh`.

The per-edge constant `I_ij` is computed once at the start from the input
3-D geometry via the Bowers-Stephenson identity
    `I_ij = (ℓ_ij² − r_i² − r_j²) / (2 r_i r_j)`,
with `r_i^(0) = (1/3) min{ℓ_e : e adj v_i}` as the default initial radii.
The user can override the initial radii by writing into the `r0`
property map before calling this function.

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

\param mesh  Input triangle mesh.
\param np    Named parameters:
  - `vertex_curvature_map(pmap)` — per-vertex Θ_v target.
  - `fixed_vertex_map(pmap)` — pinning override.
  - `gradient_tolerance(ε)` — Newton stop.
  - `max_iterations(n)` — Newton iteration cap.

\returns A `Conformal_map_result<FT>` with `u_per_vertex[v] = log r_v`
         (the converged log-radius at each vertex).

\pre `mesh` is a triangle mesh with positive edge lengths.
\pre The user-supplied or natural-theta Θ satisfies Gauss–Bonnet.

\note  Convergence is sensitive to the initial point and to extreme
       `I_ij` values.  For testing purposes the natural-theta default
       (Θ_v shifted so that u = 0 is the equilibrium) always converges
       in zero iterations.
*/
template <typename TriangleMesh,
          typename CGAL_NP_TEMPLATE_PARAMETERS>
auto discrete_inversive_distance_map(
    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_inversive_distance_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_inversive_distance_maps(mesh);
    ::conformallab::compute_inversive_distance_init_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 first vertex by default; user can override with fixed_vertex_map.
    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.
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    if constexpr (!has_theta) {
        auto G0 = ::conformallab::inversive_distance_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_inversive_distance(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 layout step (Phase 8b-Lite extension) ─────────────────────
    //
    // If the caller supplied `output_uv_map(pmap)`, lay out the converged
    // packing in ℝ² and write per-vertex `Point_2` coordinates into `pmap`.
    //
    // Method: the converged Inversive-Distance radii `r_i = exp(u_i)`
    // together with the fixed per-edge `I_ij` constants determine effective
    // Euclidean edge lengths via the Bowers-Stephenson identity
    //   ℓᵢⱼ² = rᵢ² + rⱼ² + 2·Iᵢⱼ·rᵢ·rⱼ
    // so we can populate a temporary `EuclideanMaps` whose `lambda0` carries
    // `log(ℓᵢⱼ²)` per edge and then reuse `euclidean_layout(mesh, 0, eucl)`
    // — the existing priority-BFS trilateration on the resulting triangle
    // metric.  All vertex/edge DOF indices stay at −1 (pinned), so the empty
    // DOF vector `0` produces lengths driven purely by `lambda0`.
    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 eucl = ::conformallab::setup_euclidean_maps(mesh);
            for (auto e : mesh.edges()) {
                auto h = mesh.halfedge(e);
                const double u_i = result.u_per_vertex[mesh.source(h).idx()];
                const double u_j = result.u_per_vertex[mesh.target(h).idx()];
                const double I   = maps.I_e[e];
                const double l2  = ::conformallab::id_detail::edge_length_squared(u_i, u_j, I);
                eucl.lambda0[e]  = (l2 > 0.0) ? std::log(l2) : -30.0;
            }
            // Empty DOF vector: every vertex is pinned (idx=-1), so the
            // layout depends purely on the lambda0 we just computed.
            std::vector<double> zero;
            auto layout = ::conformallab::euclidean_layout(mesh, zero, eucl);

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

}  // namespace CGAL

#endif  // CGAL_DISCRETE_INVERSIVE_DISTANCE_H