957f506679
Follow-up to S2: surface the new NewtonResult diagnostics through the public CGAL result types so callers of the high-level API see them too. - Add `CGAL::Newton_status` (alias of conformallab::NewtonStatus) and three fields — `status`, `sparse_qr_fallback_used`, `min_ldlt_pivot` — to Conformal_map_result, Hyper_ideal_map_result and Circle_packing_result. (sparse_qr_fallback_used already existed on Conformal_map_result but was never populated; it is now wired through.) - Map them from NewtonResult at all five entry points (euclidean / spherical / hyper_ideal / inversive_distance / cp_euclidean). - Test: SingleTriangleConverges now asserts status == Converged and the diagnostics propagate. Additive only — existing `converged`/`iterations`/`gradient_norm` semantics unchanged. 301/301 CGAL tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
599 lines
25 KiB
C++
599 lines
25 KiB
C++
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
//
|
||
// Package: conformallab++ / Discrete_conformal_map (Phase 8 MVP, 2026-05-19)
|
||
|
||
/*!
|
||
\file CGAL/Discrete_conformal_map.h
|
||
\ingroup PkgConformalMapRef
|
||
|
||
User-facing entry points for the Discrete_conformal_map package (Phase 8b-Lite).
|
||
|
||
This header provides three functions covering all three DCE geometries:
|
||
|
||
- `CGAL::discrete_conformal_map_euclidean` — flat conformal map (ℝ²), open or closed mesh
|
||
- `CGAL::discrete_conformal_map_spherical` — spherical uniformisation (S²), genus-0 mesh
|
||
- `CGAL::discrete_conformal_map_hyper_ideal` — hyperbolic conformal map (H²), genus ≥ 1
|
||
|
||
For circle-packing models see the companion headers:
|
||
- `<CGAL/Discrete_circle_packing.h>` — `discrete_circle_packing_euclidean` (CP-Euclidean)
|
||
- `<CGAL/Discrete_inversive_distance.h>` — `discrete_inversive_distance_map` (Luo 2004)
|
||
|
||
\section Example Simplest usage
|
||
|
||
\code{.cpp}
|
||
#include <CGAL/Simple_cartesian.h>
|
||
#include <CGAL/Surface_mesh.h>
|
||
#include <CGAL/Discrete_conformal_map.h>
|
||
|
||
using K = CGAL::Simple_cartesian<double>;
|
||
using Mesh = CGAL::Surface_mesh<K::Point_3>;
|
||
|
||
int main() {
|
||
Mesh mesh = ...; // load a triangle mesh
|
||
auto result = CGAL::discrete_conformal_map_euclidean(mesh);
|
||
if (!result.converged)
|
||
return 1;
|
||
// result.u_per_vertex[v] now holds the conformal scale factor at v.
|
||
}
|
||
\endcode
|
||
|
||
\section NamedParams Tuning via named parameters
|
||
|
||
\code{.cpp}
|
||
auto result = CGAL::discrete_conformal_map_euclidean(
|
||
mesh,
|
||
CGAL::parameters::gradient_tolerance(1e-12)
|
||
.max_iterations(500));
|
||
\endcode
|
||
|
||
\sa `CGAL::Default_conformal_map_traits`
|
||
\sa `CGAL::parameters::vertex_curvature_map`
|
||
*/
|
||
|
||
#ifndef CGAL_DISCRETE_CONFORMAL_MAP_H
|
||
#define CGAL_DISCRETE_CONFORMAL_MAP_H
|
||
|
||
#include <CGAL/Conformal_map_traits.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/property_map.h>
|
||
|
||
// Existing implementation headers (Layer 1 — unchanged).
|
||
#include "../euclidean_functional.hpp"
|
||
#include "../layout.hpp"
|
||
#include "../spherical_functional.hpp"
|
||
#include "../hyper_ideal_functional.hpp"
|
||
#include "../gauss_bonnet.hpp"
|
||
#include "../newton_solver.hpp"
|
||
|
||
#include <vector>
|
||
#include <unordered_map>
|
||
|
||
namespace CGAL {
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Result type
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
/// Why the Newton solve terminated (re-exported from the solver; I1 audit).
|
||
/// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`.
|
||
using Newton_status = ::conformallab::NewtonStatus;
|
||
|
||
/*!
|
||
\ingroup PkgConformalMapRef
|
||
|
||
Result of `discrete_conformal_map_euclidean`. Carries the converged
|
||
scale factors `u_v`, Newton diagnostics, and the convergence flag.
|
||
*/
|
||
template <typename FT = double>
|
||
struct Conformal_map_result
|
||
{
|
||
/// Conformal scale factor `u_v` per vertex (indexed by raw vertex index).
|
||
/// Length: `num_vertices(mesh)`.
|
||
std::vector<FT> u_per_vertex;
|
||
|
||
/// Number of Newton iterations performed.
|
||
int iterations = 0;
|
||
|
||
/// `‖G(u*)‖∞` at termination.
|
||
FT gradient_norm = FT(0);
|
||
|
||
/// `true` iff `gradient_norm < gradient_tolerance`.
|
||
bool converged = false;
|
||
|
||
/// Why the solve stopped (more informative than `converged` alone): one of
|
||
/// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`.
|
||
Newton_status status = Newton_status::MaxIterations;
|
||
|
||
/// `true` iff the linear solver used the SparseQR fallback at any
|
||
/// Newton step (gauge mode on closed mesh without pinned vertex).
|
||
bool sparse_qr_fallback_used = false;
|
||
|
||
/// Smallest `|Dᵢᵢ|` of the last LDLT factorisation — a cheap
|
||
/// near-singularity proxy (small ⇒ ill-conditioned solve). `0` if LDLT
|
||
/// never succeeded.
|
||
FT min_ldlt_pivot = FT(0);
|
||
};
|
||
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// discrete_conformal_map_euclidean — user-facing entry
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
/*!
|
||
\ingroup PkgConformalMapRef
|
||
|
||
Compute the Euclidean discrete-conformal map of `mesh`.
|
||
|
||
This is the user-facing entry of Phase 8 MVP. Internally it delegates
|
||
to the existing implementation in `code/include/euclidean_functional.hpp`
|
||
and `code/include/newton_solver.hpp` (Phase 1–7), so the algorithmic
|
||
behaviour is identical to the legacy API; this function only changes
|
||
the public façade.
|
||
|
||
\tparam TriangleMesh A `CGAL::Surface_mesh<P>` for some point type `P`.
|
||
Other `FaceGraph` models are planned for Phase 8a.2.
|
||
\tparam NamedParameters Optional CGAL named-parameter pack.
|
||
|
||
\param mesh The input mesh (modified in place: property maps are attached).
|
||
\param np Named parameters:
|
||
\cgalParamNBegin{vertex_curvature_map}
|
||
\cgalParamDescription{Property map `vertex → FT` of target cone angles Θᵥ.}
|
||
\cgalParamDefault{2π at interior vertices, π at boundary vertices.}
|
||
\cgalParamNEnd
|
||
\cgalParamNBegin{gradient_tolerance}
|
||
\cgalParamDescription{Newton stops when `‖G(u)‖∞ < tol`.}
|
||
\cgalParamDefault{`1e-10`}
|
||
\cgalParamNEnd
|
||
\cgalParamNBegin{max_iterations}
|
||
\cgalParamDescription{Hard limit on Newton steps.}
|
||
\cgalParamDefault{`200`}
|
||
\cgalParamNEnd
|
||
\cgalParamNBegin{fixed_vertex_map}
|
||
\cgalParamDescription{Property map `vertex → bool`; `true` ⇒ pinned.}
|
||
\cgalParamDefault{The first vertex in `mesh.vertices()` is pinned.}
|
||
\cgalParamNEnd
|
||
|
||
\returns A `Conformal_map_result<FT>` carrying `u_v` and Newton diagnostics.
|
||
|
||
\pre `mesh` is a triangle mesh.
|
||
\pre `mesh` satisfies the Gauss–Bonnet relation
|
||
`Σ(2π − Θᵥ) = 2π·χ(mesh)` for the chosen target curvature map.
|
||
*/
|
||
template <typename TriangleMesh,
|
||
typename CGAL_NP_TEMPLATE_PARAMETERS>
|
||
auto discrete_conformal_map_euclidean(
|
||
TriangleMesh& mesh,
|
||
const CGAL_NP_CLASS& np = parameters::default_values())
|
||
{
|
||
// ── Type plumbing ──────────────────────────────────────────────────────
|
||
//
|
||
// Deduce the kernel from the mesh's Point_3 type rather than hard-coding
|
||
// Simple_cartesian<double>. This lets the wrapper work with any
|
||
// Surface_mesh<P> whose P is a CGAL kernel point. The user can override
|
||
// the entire traits class via the `geom_traits(...)` named parameter
|
||
// (Phase 8b.2 extension; default below covers the common case).
|
||
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;
|
||
|
||
// ── 1. Set up property maps (legacy layer) ─────────────────────────────
|
||
auto maps = ::conformallab::setup_euclidean_maps(mesh);
|
||
::conformallab::compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// ── 2. Target curvature: user-supplied or "natural-theta" default ─────
|
||
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) {
|
||
// User-provided Θ: copy into the property map and verify Gauss–Bonnet.
|
||
// Throws std::runtime_error if the user-supplied Θ violates GB.
|
||
for (auto v : mesh.vertices())
|
||
maps.theta_v[v] = get(theta_param, v);
|
||
::conformallab::check_gauss_bonnet(mesh, maps);
|
||
}
|
||
// If no Θ is supplied, the default behaviour is "natural-theta": set Θ
|
||
// such that x = 0 is the natural equilibrium (the actual angle sums at
|
||
// x = 0 become the targets). This matches the convention of the
|
||
// existing test suite and guarantees that the default invocation
|
||
// converges immediately for any well-formed triangle mesh.
|
||
// The actual Θ adjustment is done after DOF assignment (step 5b below).
|
||
|
||
// ── 3. Pin vertices: user map, or the first vertex by default ──────────
|
||
//
|
||
// The legacy v_idx property map has -1 as default (= pinned). We must
|
||
// first mark every vertex as "free" (any non-negative sentinel), then
|
||
// pin the requested ones, then assign sequential DOF indices.
|
||
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;
|
||
}
|
||
}
|
||
|
||
// ── 4. Assign DOF indices 0..n−1 to non-pinned vertices ─────────────────
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices())
|
||
if (maps.v_idx[v] != -1)
|
||
maps.v_idx[v] = idx++;
|
||
|
||
// ── 5. Read tolerances ─────────────────────────────────────────────────
|
||
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);
|
||
|
||
// ── 5b. Natural-theta default: shift Θ so that x = 0 is the equilibrium
|
||
//
|
||
// Only applied when the user did NOT supply a vertex_curvature_map.
|
||
// The trick: evaluate G at x = 0, then subtract G_v from Θ_v. After
|
||
// this shift the new G(0) is identically zero, so Newton starts at the
|
||
// optimum and immediately reports "converged". This matches the
|
||
// contract of the existing test suite ("natural-theta" pattern).
|
||
std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
|
||
if constexpr (!has_theta) {
|
||
auto G0 = ::conformallab::euclidean_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)];
|
||
}
|
||
}
|
||
|
||
// ── 6. Newton on x_0 = 0 ───────────────────────────────────────────────
|
||
auto nr = ::conformallab::newton_euclidean(mesh, x0, maps, tol, max_iter);
|
||
|
||
// ── 7. Pack result: u_v for every vertex, including pinned (u=0) ──────
|
||
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)];
|
||
// else: pinned ⇒ u_v stays 0
|
||
}
|
||
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);
|
||
|
||
// ── 8. Optional layout step (Phase 8b-Lite extension) ──────────────────
|
||
//
|
||
// If the caller supplied `output_uv_map(pmap)`, run the priority-BFS
|
||
// trilateration on the converged x and write per-vertex `Point_2`
|
||
// coordinates into `pmap`. Optional `normalise_layout(true)` applies
|
||
// the canonical PCA centroid + major-axis normalisation.
|
||
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 layout = ::conformallab::euclidean_layout(mesh, nr.x, maps);
|
||
|
||
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;
|
||
}
|
||
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// 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
|
||
Gauss–Bonnet 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;
|
||
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);
|
||
|
||
// Optional 3-D layout step (point on S²)
|
||
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 layout = ::conformallab::spherical_layout(mesh, nr.x, maps);
|
||
const bool do_norm = parameters::choose_parameter(
|
||
parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout),
|
||
false);
|
||
if (do_norm) ::conformallab::normalise_spherical(layout);
|
||
for (auto v : mesh.vertices()) {
|
||
const auto& p = layout.pos[v.idx()];
|
||
put(uv_param, v,
|
||
typename Traits::Kernel::Point_3(p.x(), p.y(), p.z()));
|
||
}
|
||
}
|
||
}
|
||
|
||
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;
|
||
|
||
/// 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);
|
||
};
|
||
|
||
/*!
|
||
\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;
|
||
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);
|
||
|
||
// Optional Poincaré-disk layout (2-D in the unit disk).
|
||
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 layout = ::conformallab::hyper_ideal_layout(mesh, nr.x, maps);
|
||
const bool do_norm = parameters::choose_parameter(
|
||
parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout),
|
||
false);
|
||
if (do_norm) ::conformallab::normalise_hyperbolic(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()));
|
||
}
|
||
}
|
||
}
|
||
|
||
(void)n; // already-counted by maps; silence unused-var warnings if any
|
||
return result;
|
||
}
|
||
|
||
} // namespace CGAL
|
||
|
||
#endif // CGAL_DISCRETE_CONFORMAL_MAP_H
|