diff --git a/code/include/CGAL/Conformal_layout.h b/code/include/CGAL/Conformal_layout.h new file mode 100644 index 0000000..8dc4311 --- /dev/null +++ b/code/include/CGAL/Conformal_layout.h @@ -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 +#include +#include + +#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 +Layout2D euclidean_layout( + TriangleMesh& mesh, + const std::vector& 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 +Layout3D spherical_layout( + TriangleMesh& mesh, + const std::vector& 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 +Layout2D hyper_ideal_layout( + TriangleMesh& mesh, + const std::vector& 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 diff --git a/code/include/CGAL/Discrete_circle_packing.h b/code/include/CGAL/Discrete_circle_packing.h new file mode 100644 index 0000000..175975d --- /dev/null +++ b/code/include/CGAL/Discrete_circle_packing.h @@ -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 +#include +#include +#include +#include +#include +#include + +#include "../cp_euclidean_functional.hpp" +#include "../newton_solver.hpp" + +namespace CGAL { + +// ── Default traits for CP-Euclidean ─────────────────────────────────────────── + +template > +struct Default_cp_euclidean_traits; + +template +struct Default_cp_euclidean_traits, K> +{ + using Kernel = K; + using FT = typename K::FT; + using Point_3 = typename K::Point_3; + using Triangle_mesh = CGAL::Surface_mesh; + + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using Edge_descriptor = typename boost::graph_traits::edge_descriptor; + using Face_descriptor = typename boost::graph_traits::face_descriptor; + + // CP-Euclidean property maps — note the *face* DOF index map. + using Face_index_pmap = typename Triangle_mesh::template Property_map; + using Theta_e_pmap = typename Triangle_mesh::template Property_map; + using Phi_f_pmap = typename Triangle_mesh::template Property_map; +}; + +// ── 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 +struct Circle_packing_result +{ + /// Face DOFs `ρ_f = log R_f` (length = num_faces(mesh); pinned face = 0). + std::vector 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

`. +\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` 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 +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::Kernel; + using Default_traits = Default_cp_euclidean_traits; + 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 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 x0(static_cast(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(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(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 diff --git a/code/include/CGAL/Discrete_conformal_map.h b/code/include/CGAL/Discrete_conformal_map.h index 64de271..3a42b16 100644 --- a/code/include/CGAL/Discrete_conformal_map.h +++ b/code/include/CGAL/Discrete_conformal_map.h @@ -59,6 +59,8 @@ auto result = CGAL::discrete_conformal_map_euclidean( // Existing implementation headers (Layer 1 — unchanged). #include "../euclidean_functional.hpp" +#include "../spherical_functional.hpp" +#include "../hyper_ideal_functional.hpp" #include "../gauss_bonnet.hpp" #include "../newton_solver.hpp" @@ -269,6 +271,220 @@ auto discrete_conformal_map_euclidean( 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

` 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` 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 + Gauss–Bonnet relation `Σ(2π − Θᵥ) = 4π` (sphere). +*/ +template +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::Kernel; + using Default_traits = Default_conformal_map_traits; + 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 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 x0(static_cast(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(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(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 +struct Hyper_ideal_map_result +{ + /// Vertex DOFs `b_v` (length = num_vertices(mesh); pinned vertices = 0). + std::vector b_per_vertex; + /// Edge DOFs `a_e` (length = num_edges(mesh); pinned edges = 0). + std::vector 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 +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::Kernel; + using Default_traits = Default_conformal_map_traits; + 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 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 x0(static_cast(n), 0.0); + for (auto v : mesh.vertices()) { + int i = maps.v_idx[v]; + if (i >= 0) x0[static_cast(i)] = 1.0; + } + for (auto e : mesh.edges()) { + int i = maps.e_idx[e]; + if (i >= 0) x0[static_cast(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(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(j)]; + } + 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 diff --git a/code/include/CGAL/Discrete_inversive_distance.h b/code/include/CGAL/Discrete_inversive_distance.h new file mode 100644 index 0000000..341b346 --- /dev/null +++ b/code/include/CGAL/Discrete_inversive_distance.h @@ -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 +#include +#include +#include +#include +#include +#include + +#include // for Conformal_map_result + +#include "../inversive_distance_functional.hpp" +#include "../newton_solver.hpp" + +namespace CGAL { + +// ── Default traits for Inversive-Distance ──────────────────────────────────── + +template > +struct Default_inversive_distance_traits; + +template +struct Default_inversive_distance_traits, K> +{ + using Kernel = K; + using FT = typename K::FT; + using Point_3 = typename K::Point_3; + using Triangle_mesh = CGAL::Surface_mesh; + + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Edge_descriptor = typename boost::graph_traits::edge_descriptor; + + // Inversive-distance specific property maps. + using Vertex_index_pmap = typename Triangle_mesh::template Property_map; + using Theta_v_pmap = typename Triangle_mesh::template Property_map; + using R0_pmap = typename Triangle_mesh::template Property_map; + using I_e_pmap = typename Triangle_mesh::template Property_map; +}; + +// ── 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

`. +\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` 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 +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::Kernel; + using Default_traits = Default_inversive_distance_traits; + 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 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 x0(static_cast(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(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(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 diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 65479b6..b26795c 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -84,6 +84,11 @@ add_executable(conformallab_cgal_tests # Convergence tests for newton_cp_euclidean (analytic Hessian) and # newton_inversive_distance (FD Hessian). test_newton_phase9a.cpp + + # ── Phase 8b-Lite: CGAL entry wrappers for the 4 non-Euclidean modes ───── + # Spherical, HyperIdeal, CircleP-Euclidean, Inversive-Distance via + # public API + Conformal_layout.h wrapper. + test_cgal_phase8b_lite.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_cgal_phase8b_lite.cpp b/code/tests/cgal/test_cgal_phase8b_lite.cpp new file mode 100644 index 0000000..f23e2be --- /dev/null +++ b/code/tests/cgal/test_cgal_phase8b_lite.cpp @@ -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 +#include +#include +#include + +#include "mesh_builder.hpp" +#include "conformal_mesh.hpp" + +#include +#include + +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 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())); + } +}