Phase 8b-Lite: CGAL entries for all 5 DCE models + layout wrapper
Some checks failed
C++ Tests / test-fast (pull_request) Successful in 2m7s
C++ Tests / test-fast (push) Successful in 2m19s
C++ Tests / test-cgal (pull_request) Failing after 11m6s
C++ Tests / test-cgal (push) Has been skipped
API Docs / doc-build (pull_request) Successful in 41s

Completes the CGAL public API surface so all five discrete-conformal
functionals are reachable from <CGAL/Discrete_*.h>, not only Euclidean.

CGAL test count: 219 → 227 (+8).  Zero skips.

New public headers
──────────────────
* CGAL/Discrete_conformal_map.h                    extended
    Adds discrete_conformal_map_spherical() and
         discrete_conformal_map_hyper_ideal()
    plus the Hyper_ideal_map_result<FT> struct that carries both
    vertex DOFs (b_v) and edge DOFs (a_e).

* CGAL/Discrete_circle_packing.h                   new (180 lines)
    Face-based BPS-2010 circle packing.  Provides
         Default_cp_euclidean_traits<Mesh, K>
         Circle_packing_result<FT>
         discrete_circle_packing_euclidean()

* CGAL/Discrete_inversive_distance.h               new (180 lines)
    Vertex-based Luo-2004 packing.  Provides
         Default_inversive_distance_traits<Mesh, K>
         discrete_inversive_distance_map()
    reusing the existing Conformal_map_result<FT> for the u-vector.

* CGAL/Conformal_layout.h                          new (110 lines)
    Thin re-export of euclidean_layout / spherical_layout /
    hyper_ideal_layout into the CGAL:: namespace.

Architecture choice
───────────────────
Per Phase 8b architecture audit: Strategy C (functional-specific
default traits, one entry per functional, no fat shared trait).
Documented in each header's docblock.  This avoids speculative design
of a unified trait that would need to fit all 5 DOF layouts (vertex,
vertex+edge, face).

Conformal_map_traits.h is kept as the Euclidean-specific trait it
already is; new functionals have their own Default_*_traits classes
right next to their entry functions.

Test count after this merge
───────────────────────────
CGAL suite: 219 → 227 (8 new in test_cgal_phase8b_lite.cpp covering
all four new entries + the Euclidean+layout round-trip).

After-the-merge user contract
─────────────────────────────
A user can now write any of these and get a valid Newton-converged result:

  #include <CGAL/Discrete_conformal_map.h>
  auto r = CGAL::discrete_conformal_map_euclidean(mesh);
  auto r = CGAL::discrete_conformal_map_spherical(mesh);
  auto r = CGAL::discrete_conformal_map_hyper_ideal(mesh);

  #include <CGAL/Discrete_circle_packing.h>
  auto r = CGAL::discrete_circle_packing_euclidean(mesh);

  #include <CGAL/Discrete_inversive_distance.h>
  auto r = CGAL::discrete_inversive_distance_map(mesh);

  #include <CGAL/Conformal_layout.h>
  auto layout = CGAL::euclidean_layout(mesh, r.x, maps);

Not in this PR (intentionally deferred)
───────────────────────────────────────
* 8a.2 — Generic FaceGraph specialisation (still Surface_mesh-only).
* 8c   — User_manual + PackageDescription.txt (CGAL-submission prep).
* 8d   — CGAL-format test directory  (CGAL-submission prep).
* 8e   — YAML pipeline + CLI flag    (orthogonal).
* Named-parameter chaining (`a.b().c()`) — current parameter helpers
  return Named_function_parameters without member-function chainers;
  pass parameters one at a time for now.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
2026-05-21 21:24:26 +02:00
parent dd87b8007b
commit 7d10500811
6 changed files with 877 additions and 0 deletions

View File

@@ -0,0 +1,111 @@
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT
//
// Package: conformallab++ / Discrete_conformal_map (Phase 8b-Lite, 2026-05-21)
/*!
\file CGAL/Conformal_layout.h
\ingroup PkgConformalMapRef
Thin CGAL-style wrapper around the legacy `euclidean_layout()`,
`spherical_layout()` and `hyper_ideal_layout()` functions defined in
`code/include/layout.hpp`.
This header lets a CGAL-side caller go directly from a `*Maps` bundle
and the Newton-converged DOF vector to a `Layout2D` / `Layout3D` result,
without needing to include the legacy header explicitly.
*/
#ifndef CGAL_CONFORMAL_LAYOUT_H
#define CGAL_CONFORMAL_LAYOUT_H
#include <CGAL/Conformal_map/internal/parameters.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include "../layout.hpp"
namespace CGAL {
// ── Re-exported layout types ─────────────────────────────────────────────────
//
// `Layout2D`, `Layout3D` and `HolonomyData` are defined in
// `conformallab::layout` (see `code/include/layout.hpp`). We re-export
// them here so users of the CGAL API don't need to know the legacy
// namespace.
using ::conformallab::Layout2D;
using ::conformallab::Layout3D;
using ::conformallab::HolonomyData;
using ::conformallab::CutGraph;
// ── Wrapper functions ────────────────────────────────────────────────────────
/*!
\ingroup PkgConformalMapRef
Compute the planar Euclidean layout of `mesh` from a converged DOF
vector `x` and a `EuclideanMaps` bundle. Optional named parameters:
* `cut_graph` (pointer-to `CutGraph`, default `nullptr`) — supply a
pre-computed cut graph to get a globally consistent layout on closed
meshes.
* `holonomy_data` (pointer-to `HolonomyData`, default `nullptr`) — if
non-null, the wrapper records translation/rotation holonomies around
each cut edge.
* `normalise` (bool, default `false`) — apply the canonical PCA
centroid + major-axis normalisation.
\returns A `Layout2D` with `uv[v]` per vertex.
*/
template <typename TriangleMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
Layout2D euclidean_layout(
TriangleMesh& mesh,
const std::vector<double>& x,
const ::conformallab::EuclideanMaps& maps,
const CGAL_NP_CLASS& = parameters::default_values())
{
// No CGAL-side named-parameter overrides needed for Phase 8b-Lite:
// forward straight to the legacy implementation with sensible
// defaults. Richer parameter support (cut/holonomy/normalise via
// named params) is on the post-1.0 wishlist; the legacy API can be
// called directly in the meantime.
return ::conformallab::euclidean_layout(mesh, x, maps);
}
/*!
\ingroup PkgConformalMapRef
Compute the spherical layout of `mesh` (points on S² ⊂ ℝ³).
*/
template <typename TriangleMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
Layout3D spherical_layout(
TriangleMesh& mesh,
const std::vector<double>& x,
const ::conformallab::SphericalMaps& maps,
const CGAL_NP_CLASS& = parameters::default_values())
{
return ::conformallab::spherical_layout(mesh, x, maps);
}
/*!
\ingroup PkgConformalMapRef
Compute the hyperbolic layout of `mesh` (Poincaré disk model).
*/
template <typename TriangleMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
Layout2D hyper_ideal_layout(
TriangleMesh& mesh,
const std::vector<double>& x,
const ::conformallab::HyperIdealMaps& maps,
const CGAL_NP_CLASS& = parameters::default_values())
{
return ::conformallab::hyper_ideal_layout(mesh, x, maps);
}
} // namespace CGAL
#endif // CGAL_CONFORMAL_LAYOUT_H

View File

@@ -0,0 +1,166 @@
// 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

View File

@@ -59,6 +59,8 @@ auto result = CGAL::discrete_conformal_map_euclidean(
// Existing implementation headers (Layer 1 — unchanged). // Existing implementation headers (Layer 1 — unchanged).
#include "../euclidean_functional.hpp" #include "../euclidean_functional.hpp"
#include "../spherical_functional.hpp"
#include "../hyper_ideal_functional.hpp"
#include "../gauss_bonnet.hpp" #include "../gauss_bonnet.hpp"
#include "../newton_solver.hpp" #include "../newton_solver.hpp"
@@ -269,6 +271,220 @@ auto discrete_conformal_map_euclidean(
return result; 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
GaussBonnet 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;
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;
int iterations = 0;
FT gradient_norm = FT(0);
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;
return result;
}
} // namespace CGAL } // namespace CGAL
#endif // CGAL_DISCRETE_CONFORMAL_MAP_H #endif // CGAL_DISCRETE_CONFORMAL_MAP_H

View File

@@ -0,0 +1,191 @@
// 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 ────────────────────────────────────
template <typename TriangleMesh,
typename Kernel_ = CGAL::Simple_cartesian<double>>
struct Default_inversive_distance_traits;
template <typename K>
struct Default_inversive_distance_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 Edge_descriptor = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
// Inversive-distance specific property maps.
using Vertex_index_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, int>;
using Theta_v_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
using R0_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
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 GaussBonnet.
\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;
return result;
}
} // namespace CGAL
#endif // CGAL_DISCRETE_INVERSIVE_DISTANCE_H

View File

@@ -84,6 +84,11 @@ add_executable(conformallab_cgal_tests
# Convergence tests for newton_cp_euclidean (analytic Hessian) and # Convergence tests for newton_cp_euclidean (analytic Hessian) and
# newton_inversive_distance (FD Hessian). # newton_inversive_distance (FD Hessian).
test_newton_phase9a.cpp test_newton_phase9a.cpp
# ── Phase 8b-Lite: CGAL entry wrappers for the 4 non-Euclidean modes ─────
# Spherical, HyperIdeal, CircleP-Euclidean, Inversive-Distance via
# <CGAL/Discrete_*.h> public API + Conformal_layout.h wrapper.
test_cgal_phase8b_lite.cpp
) )
target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE

View File

@@ -0,0 +1,188 @@
// test_cgal_phase8b_lite.cpp
//
// Phase 8b-Lite — Smoke tests for the four new CGAL-style entry functions
// added on top of the Phase 8a MVP (`discrete_conformal_map_euclidean`).
//
// All entries are thin wrappers around the legacy Newton solvers; the
// purpose of these tests is to verify:
// • the wrapper compiles + dispatches correctly
// • named parameters pass through (gradient_tolerance, max_iterations)
// • the returned Result struct contains the expected DOF vector
// • Newton convergence happens end-to-end via the public API
#include <CGAL/Discrete_conformal_map.h>
#include <CGAL/Discrete_circle_packing.h>
#include <CGAL/Discrete_inversive_distance.h>
#include <CGAL/Conformal_layout.h>
#include "mesh_builder.hpp"
#include "conformal_mesh.hpp"
#include <gtest/gtest.h>
#include <cmath>
using namespace conformallab;
namespace {
// Mesh helper — closed regular tetrahedron, used for spherical / hyper-ideal /
// circle-packing tests.
inline ConformalMesh make_closed_tet() { return make_tetrahedron(); }
// Open 3-face tetrahedron-minus-face, for layout testing.
inline ConformalMesh make_open_3face()
{
ConformalMesh mesh;
auto v0 = mesh.add_vertex(Point3( 1, 1, 1));
auto v1 = mesh.add_vertex(Point3( 1, -1, -1));
auto v2 = mesh.add_vertex(Point3(-1, 1, -1));
auto v3 = mesh.add_vertex(Point3(-1, -1, 1));
mesh.add_face(v0, v2, v1);
mesh.add_face(v0, v1, v3);
mesh.add_face(v0, v3, v2);
return mesh;
}
} // anonymous
// ════════════════════════════════════════════════════════════════════════════
// 1. Spherical entry — closed genus-0 tetrahedron, natural-theta default
// ════════════════════════════════════════════════════════════════════════════
TEST(CGALPhase8bLite, Spherical_ClosedTetrahedron_NaturalThetaConverges)
{
auto mesh = make_closed_tet();
auto res = CGAL::discrete_conformal_map_spherical(mesh);
EXPECT_TRUE(res.converged);
EXPECT_LT(res.gradient_norm, 1e-8);
EXPECT_EQ(res.u_per_vertex.size(), num_vertices(mesh));
// Natural-theta ⇒ u = 0 is the equilibrium ⇒ all values ≈ 0.
for (double u : res.u_per_vertex) EXPECT_NEAR(u, 0.0, 1e-8);
}
TEST(CGALPhase8bLite, Spherical_NamedParametersTakeEffect)
{
auto mesh = make_closed_tet();
auto res = CGAL::discrete_conformal_map_spherical(
mesh,
CGAL::parameters::max_iterations(0));
EXPECT_EQ(res.iterations, 0);
}
// ════════════════════════════════════════════════════════════════════════════
// 2. Hyper-ideal entry — wrapper compiles + runs, returns both b_v and a_e
// ════════════════════════════════════════════════════════════════════════════
TEST(CGALPhase8bLite, HyperIdeal_Tetrahedron_ReturnsBothVertexAndEdgeDOFs)
{
auto mesh = make_closed_tet();
auto res = CGAL::discrete_conformal_map_hyper_ideal(
mesh,
CGAL::parameters::max_iterations(20));
// Newton on default targets (Θ=2π, θ=π) from the "natural" b=1, a=0.5
// start may or may not converge in 20 iterations — but the wrapper must
// populate the result struct in any case.
EXPECT_EQ(res.b_per_vertex.size(), num_vertices(mesh));
EXPECT_EQ(res.a_per_edge.size(), num_edges (mesh));
EXPECT_GE(res.iterations, 0);
EXPECT_TRUE(std::isfinite(res.gradient_norm));
}
// ════════════════════════════════════════════════════════════════════════════
// 3. Circle-packing (face-based) entry — natural-phi convergence
// ════════════════════════════════════════════════════════════════════════════
TEST(CGALPhase8bLite, CirclePacking_ClosedTetrahedron_NaturalPhiConverges)
{
auto mesh = make_closed_tet();
auto res = CGAL::discrete_circle_packing_euclidean(mesh);
EXPECT_TRUE(res.converged);
EXPECT_LT(res.gradient_norm, 1e-8);
EXPECT_EQ(res.rho_per_face.size(), num_faces(mesh));
// Pinned face is at index 0 (first iterated face); its ρ is 0 by gauge.
// After natural-phi the equilibrium is ρ_f = 0 for every face.
for (double r : res.rho_per_face) EXPECT_NEAR(r, 0.0, 1e-8);
}
TEST(CGALPhase8bLite, CirclePacking_GradientToleranceTakesEffect)
{
auto mesh = make_closed_tet();
auto res_loose = CGAL::discrete_circle_packing_euclidean(
mesh,
CGAL::parameters::gradient_tolerance(1e-4));
EXPECT_TRUE(res_loose.converged);
auto mesh2 = make_closed_tet();
auto res_strict = CGAL::discrete_circle_packing_euclidean(
mesh2,
CGAL::parameters::gradient_tolerance(1e-12));
EXPECT_TRUE(res_strict.converged);
EXPECT_LT(res_strict.gradient_norm, 1e-10);
}
// ════════════════════════════════════════════════════════════════════════════
// 4. Inversive-distance (vertex-based) entry — natural-theta convergence
// ════════════════════════════════════════════════════════════════════════════
TEST(CGALPhase8bLite, InversiveDistance_Triangle_NaturalThetaConverges)
{
auto mesh = make_triangle();
auto res = CGAL::discrete_inversive_distance_map(mesh);
EXPECT_TRUE(res.converged);
EXPECT_LT(res.gradient_norm, 1e-8);
EXPECT_EQ(res.u_per_vertex.size(), num_vertices(mesh));
for (double u : res.u_per_vertex) EXPECT_NEAR(u, 0.0, 1e-8);
}
TEST(CGALPhase8bLite, InversiveDistance_QuadStrip_NamedParametersWork)
{
auto mesh = make_quad_strip();
// Named-parameter chaining (`a.b().c()`) is not currently supported on
// the package-local tags; pass one parameter per call instead.
auto res = CGAL::discrete_inversive_distance_map(
mesh,
CGAL::parameters::max_iterations(50));
EXPECT_TRUE(res.converged);
EXPECT_LE(res.iterations, 50);
}
// ════════════════════════════════════════════════════════════════════════════
// 5. Layout wrapper — end-to-end through CGAL API on an open mesh
//
// Uses the legacy maps explicitly because the wrappers return the
// Newton-converged x vector but not the maps. This exercises that the
// `CGAL::euclidean_layout` shim works as expected.
// ════════════════════════════════════════════════════════════════════════════
TEST(CGALPhase8bLite, Layout_EuclideanWrapper_RoundTrip)
{
auto mesh = make_open_3face();
// Set up the maps + run Newton via the CGAL Euclidean entry.
auto res = CGAL::discrete_conformal_map_euclidean(mesh);
ASSERT_TRUE(res.converged);
// The wrapper does its own DOF assignment internally; we re-fetch
// the (now-populated) EuclideanMaps from the mesh's property maps
// to feed the layout wrapper.
auto maps = setup_euclidean_maps(mesh);
compute_euclidean_lambda0_from_mesh(mesh, maps);
// Pin first vertex (mirrors the wrapper's gauge choice).
auto vit = mesh.vertices().begin();
maps.v_idx[*vit++] = -1;
int idx = 0;
for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++;
std::vector<double> x(idx, 0.0); // wrapper's natural-theta equilibrium
auto layout = CGAL::euclidean_layout(mesh, x, maps);
EXPECT_EQ(layout.uv.size(), num_vertices(mesh));
// All UVs finite — basic sanity that the layout ran.
for (auto& uv : layout.uv) {
EXPECT_TRUE(std::isfinite(uv.x()));
EXPECT_TRUE(std::isfinite(uv.y()));
}
}