59a26123c8
MINOR-1 (spherical_functional.hpp:426)
Wrong comment said "second derivative < 0 for a convex functional".
The spherical energy is *concave* (NSD Hessian); the monotone-f argument
applies to both convex and concave functionals equally. Comment rewritten
to explain the actual physics: increasing scale increases all angles and
thus reduces Σ G_v.
MINOR-2 (spherical_functional.hpp:494-495)
Forward finite difference O(ε) → central finite difference O(ε²):
old: dft = (sum_Gv(t + fd_eps) - ft) / fd_eps
new: dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2*fd_eps)
Same cost when the extra sum_Gv(t - fd_eps) replaces the cached ft.
MINOR-3 (euclidean_functional.hpp, spherical_functional.hpp,
inversive_distance_functional.hpp)
New header gauss_legendre.hpp centralises the 10-point Gauss-Legendre
nodes and weights (gl10_nodes() / gl10_weights()). The three energy
functions now use the shared accessors instead of duplicated local
static arrays.
MINOR-4 (euclidean_functional.hpp, spherical_functional.hpp,
hyper_ideal_functional.hpp, inversive_distance_functional.hpp)
halfedge_to_index() centralised in conformal_mesh.hpp. All four local
aliases (eucl_hidx, spher_hidx, hidx, id_detail::hidx) now delegate to
it as one-line wrappers; the aliases are kept for now to avoid a larger
call-site churn, clearly documented as thin wrappers.
MINOR-5 (clausen.hpp:33-38)
Added a comment above inits() explaining the intentional off-by-one
return value and how it interacts with csevl() — matching the Java
Clausen.inits() / csevl() contract.
277/277 CGAL + 26/26 pure-math tests pass, 0 failed.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
392 lines
18 KiB
C++
392 lines
18 KiB
C++
#pragma once
|
||
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// inversive_distance_functional.hpp
|
||
//
|
||
// Phase 9a.2 — Inversive-distance circle-packing functional (Luo 2004).
|
||
//
|
||
// VERTEX-based circle packing. Each vertex carries a circle of radius
|
||
// r_i = exp(u_i). The inversive distance I_ij between two adjacent
|
||
// circles is a constant of the edge, derived once from the initial
|
||
// geometry via Bowers-Stephenson 2004.
|
||
//
|
||
// This is the FACE-DUAL of CPEuclideanFunctional (Phase 9a.1). The
|
||
// correspondence is I_ij = cos θ_e (Glickenstein 2011 §5).
|
||
//
|
||
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||
// │ Mathematical model │
|
||
// │ ────────────────── │
|
||
// │ │
|
||
// │ Variables: u_i = log r_i (per vertex; r_i is the radius) │
|
||
// │ Constants: I_ij (per edge; inversive distance) │
|
||
// │ Θ_v (per vertex; target cone angle) │
|
||
// │ │
|
||
// │ Bowers-Stephenson (init from initial geometry): │
|
||
// │ I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j ) │
|
||
// │ │
|
||
// │ Edge length (Luo 2004 §3, Glickenstein 2011 eq. 2.1): │
|
||
// │ ℓ_ij(u)² = exp(2 u_i) + exp(2 u_j) + 2 I_ij exp(u_i + u_j) │
|
||
// │ = r_i² + r_j² + 2 I_ij r_i r_j │
|
||
// │ │
|
||
// │ Triangle angles: same half-tangent law of cosines as the │
|
||
// │ Euclidean functional (numerically stable). │
|
||
// │ │
|
||
// │ Gradient (Luo 2004 Lemma 3.1): │
|
||
// │ ∂E/∂u_v = Θ_v − Σ_{T ∋ v} α_v(T) │
|
||
// │ │
|
||
// │ Energy: path integral E(u) = ∫₀¹ ⟨G(tu), u⟩ dt │
|
||
// │ (Luo's 1-form is closed; we use 10-point Gauss-Legendre │
|
||
// │ quadrature, identical to euclidean_functional.hpp) │
|
||
// │ │
|
||
// │ Hessian: finite-difference for the MVP port; an analytic form is │
|
||
// │ given in Glickenstein 2011 eq. (4.6) and may be added │
|
||
// │ later for performance. │
|
||
// └──────────────────────────────────────────────────────────────────────────┘
|
||
//
|
||
// Relation to euclidean_functional.hpp
|
||
// ────────────────────────────────────
|
||
// The two are structurally identical in:
|
||
// • DOF layout (per vertex), DOF index sentinel (−1 = pinned)
|
||
// • Gradient pattern (Θ − Σ α)
|
||
// • Energy via path integral (same Gauss-Legendre constants)
|
||
// • Halfedge convention (h0/h1/h2, source pattern, α opposite-edge)
|
||
//
|
||
// They differ ONLY in:
|
||
// • Per-edge constant: λ°_ij (log²-length) vs I_ij (inversive distance)
|
||
// • Edge-length formula:
|
||
// Euclidean: ℓ_ij = exp((λ°_ij + u_i + u_j) / 2)
|
||
// Inversive distance: ℓ_ij² = exp(2u_i) + exp(2u_j)
|
||
// + 2 I_ij exp(u_i + u_j)
|
||
//
|
||
// In particular at the tangential limit I_ij = 1 the inversive-distance length
|
||
// reduces to (exp(u_i) + exp(u_j))² ⇒ ℓ_ij = r_i + r_j (tangential circles),
|
||
// which is *different* from the Euclidean-conformal length even at the same
|
||
// initial geometry. The two functionals describe distinct geometric objects.
|
||
//
|
||
// Property-map name prefix: "iv:" (vertex) and "ie:" (edge).
|
||
|
||
#include "conformal_mesh.hpp"
|
||
#include "constants.hpp"
|
||
#include "gauss_legendre.hpp"
|
||
#include "euclidean_geometry.hpp" // euclidean_angles(λ12, λ23, λ31)
|
||
#include <CGAL/boost/graph/iterator.h>
|
||
#include <vector>
|
||
#include <cmath>
|
||
#include <cstdint>
|
||
#include <iostream>
|
||
|
||
namespace conformallab {
|
||
|
||
// ── Property-map type aliases ────────────────────────────────────────────────
|
||
|
||
/// Property map vertex → `int` for the Inversive-Distance functional.
|
||
using IDVMapI = ConformalMesh::Property_map<Vertex_index, int>;
|
||
/// Property map vertex → `double` for the Inversive-Distance functional.
|
||
using IDVMapD = ConformalMesh::Property_map<Vertex_index, double>;
|
||
/// Property map edge → `double` for the Inversive-Distance functional.
|
||
using IDEMapD = ConformalMesh::Property_map<Edge_index, double>;
|
||
|
||
// ── Persistent map bundle ─────────────────────────────────────────────────────
|
||
|
||
/// Bundle of the four property maps consumed by the Inversive-Distance
|
||
/// circle-packing functional (Luo 2004 / Bowers-Stephenson 2004).
|
||
struct InversiveDistanceMaps {
|
||
IDVMapI v_idx; ///< DOF index per vertex (−1 = pinned / u_v = 0)
|
||
IDVMapD theta_v; ///< target cone angle Θ_v (default 2π)
|
||
IDVMapD r0; ///< initial radius r_i^(0) (default 1)
|
||
IDEMapD I_e; ///< inversive distance I_ij (per edge, constant)
|
||
};
|
||
|
||
/// Attach the four inversive-distance property maps to `mesh` and
|
||
/// return their handles.
|
||
///
|
||
/// Defaults are intentionally trivial — every real use of this
|
||
/// functional must call `compute_inversive_distance_init_from_mesh()`
|
||
/// next to populate `r0` and `I_e` from the input geometry.
|
||
/// * `v_idx[v] = -1` (all vertices pinned initially)
|
||
/// * `theta_v[v] = 2π` (regular interior vertex)
|
||
/// * `r0[v] = 1.0` (placeholder)
|
||
/// * `I_e[e] = 1.0` (tangential default — overwritten by init step)
|
||
///
|
||
/// The maps use the `"iv:"` / `"ie:"` prefix so they do not collide
|
||
/// with the Euclidean / Spherical / HyperIdeal / CP-Euclidean maps.
|
||
inline InversiveDistanceMaps setup_inversive_distance_maps(ConformalMesh& mesh)
|
||
{
|
||
InversiveDistanceMaps m;
|
||
m.v_idx = mesh.add_property_map<Vertex_index, int> ("iv:idx", -1 ).first;
|
||
m.theta_v = mesh.add_property_map<Vertex_index, double>("iv:theta", TWO_PI ).first;
|
||
m.r0 = mesh.add_property_map<Vertex_index, double>("iv:r0", 1.0 ).first;
|
||
m.I_e = mesh.add_property_map<Edge_index, double>("ie:I", 1.0 ).first;
|
||
return m;
|
||
}
|
||
|
||
/// Assign sequential DOF indices `0..n-1` to every vertex.
|
||
///
|
||
/// **Note:** this overload assigns indices to ALL vertices unconditionally.
|
||
/// Any `v_idx` set before the call is overwritten. To pin a gauge vertex,
|
||
/// either use the two-argument overload below, or set `m.v_idx[v] = -1`
|
||
/// **after** this call. Pinning before this call has no effect.
|
||
///
|
||
/// For a closed mesh, exactly one pin is required to remove the
|
||
/// global rotational mode.
|
||
inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& mesh,
|
||
InversiveDistanceMaps& m)
|
||
{
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
|
||
return idx;
|
||
}
|
||
|
||
/// Assign sequential DOF indices to all vertices, pinning `gauge`
|
||
/// (`m.v_idx[gauge] = -1`). Use this overload on closed meshes to fix
|
||
/// the rotational gauge mode in a single call.
|
||
///
|
||
/// \returns The number of free DOFs assigned (`num_vertices − 1`).
|
||
inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& mesh,
|
||
InversiveDistanceMaps& m,
|
||
Vertex_index gauge)
|
||
{
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices())
|
||
m.v_idx[v] = (v == gauge) ? -1 : idx++;
|
||
return idx;
|
||
}
|
||
|
||
/// Count the free DOFs (vertices with `v_idx >= 0`).
|
||
inline int inversive_distance_dimension(const ConformalMesh& mesh,
|
||
const InversiveDistanceMaps& m)
|
||
{
|
||
int dim = 0;
|
||
for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim;
|
||
return dim;
|
||
}
|
||
|
||
/// Two-phase initialisation from initial mesh geometry. Mirrors the
|
||
/// role of `compute_lambda0_from_mesh` in the Euclidean functional, but
|
||
/// adapted to Luo's vertex-based radius parametrisation.
|
||
///
|
||
/// **Phase 1.** Pick a positive radius per vertex:
|
||
/// \code
|
||
/// r_i^(0) = (1/3) · min{ℓ_e : e adjacent to v_i}
|
||
/// \endcode
|
||
/// This is a heuristic — the user may override `m.r0[v]` for any
|
||
/// vertex between `setup_inversive_distance_maps()` and this call.
|
||
///
|
||
/// **Phase 2.** Compute the per-edge inversive distance via the
|
||
/// Bowers-Stephenson 2004 identity:
|
||
/// \code
|
||
/// I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j )
|
||
/// \endcode
|
||
///
|
||
/// \pre Every edge has positive 3-D length.
|
||
/// \pre Radii produced in Phase 1 are positive (degenerate isolated
|
||
/// vertices fall back to `r_i = 1`).
|
||
/// \post Every `I_e[e] > -1` for a valid packing. The chosen
|
||
/// Phase-1 heuristic keeps `I_e > 0` for most real meshes.
|
||
inline void compute_inversive_distance_init_from_mesh(ConformalMesh& mesh,
|
||
InversiveDistanceMaps& m)
|
||
{
|
||
// Phase 1: r_i = (1/3) · min adjacent edge length.
|
||
for (auto v : mesh.vertices()) {
|
||
double min_len = std::numeric_limits<double>::infinity();
|
||
for (auto h : CGAL::halfedges_around_target(v, mesh)) {
|
||
auto p1 = mesh.point(mesh.source(h));
|
||
auto p2 = mesh.point(mesh.target(h));
|
||
double dx = p1.x() - p2.x();
|
||
double dy = p1.y() - p2.y();
|
||
double dz = p1.z() - p2.z();
|
||
double len = std::sqrt(dx*dx + dy*dy + dz*dz);
|
||
if (len < min_len) min_len = len;
|
||
}
|
||
m.r0[v] = (std::isfinite(min_len) && min_len > 1e-15)
|
||
? min_len / 3.0
|
||
: 1.0;
|
||
}
|
||
|
||
// Phase 2: I_ij from initial geometry.
|
||
for (auto e : mesh.edges()) {
|
||
auto h = mesh.halfedge(e);
|
||
auto vi = mesh.source(h);
|
||
auto vj = mesh.target(h);
|
||
auto p1 = mesh.point(vi);
|
||
auto p2 = mesh.point(vj);
|
||
double dx = p1.x() - p2.x();
|
||
double dy = p1.y() - p2.y();
|
||
double dz = p1.z() - p2.z();
|
||
double l2 = dx*dx + dy*dy + dz*dz;
|
||
double ri = m.r0[vi];
|
||
double rj = m.r0[vj];
|
||
m.I_e[e] = (l2 - ri*ri - rj*rj) / (2.0 * ri * rj);
|
||
}
|
||
}
|
||
|
||
// ── Internal helpers ──────────────────────────────────────────────────────────
|
||
|
||
namespace id_detail {
|
||
|
||
inline double dof_val(int idx, const std::vector<double>& x) noexcept
|
||
{
|
||
return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
|
||
}
|
||
|
||
// halfedge_to_index is defined in conformal_mesh.hpp.
|
||
inline std::size_t hidx(Halfedge_index h) noexcept { return halfedge_to_index(h); }
|
||
|
||
// Inversive-distance edge length squared: ℓ² = exp(2u_i) + exp(2u_j) + 2 I r_i r_j
|
||
// where r_i = exp(u_i), so: ℓ² = r_i² + r_j² + 2 I r_i r_j.
|
||
// Returns -1 if the result is non-positive (degenerate; the caller skips the face).
|
||
inline double edge_length_squared(double u_i, double u_j, double I_ij) noexcept
|
||
{
|
||
double ri = std::exp(u_i);
|
||
double rj = std::exp(u_j);
|
||
double l2 = ri*ri + rj*rj + 2.0 * I_ij * ri * rj;
|
||
return l2 > 0.0 ? l2 : -1.0;
|
||
}
|
||
|
||
} // namespace id_detail
|
||
|
||
/// Inversive-Distance gradient `G_v = Θ_v − Σ_faces α_v(face)`. Same
|
||
/// half-edge corner-angle storage convention as `euclidean_gradient`.
|
||
inline std::vector<double> inversive_distance_gradient(
|
||
const ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const InversiveDistanceMaps& m)
|
||
{
|
||
const int n = inversive_distance_dimension(mesh, m);
|
||
std::vector<double> G(static_cast<std::size_t>(n), 0.0);
|
||
|
||
const std::size_t nh = mesh.number_of_halfedges();
|
||
std::vector<double> h_alpha(nh, 0.0);
|
||
|
||
// Pass 1 — per face, compute corner angles via the law of cosines.
|
||
// We reuse euclidean_angles(λ12, λ23, λ31) which takes 2·log(ℓ) per edge.
|
||
for (auto f : mesh.faces()) {
|
||
Halfedge_index h0 = mesh.halfedge(f);
|
||
Halfedge_index h1 = mesh.next(h0);
|
||
Halfedge_index h2 = mesh.next(h1);
|
||
|
||
Vertex_index v1 = mesh.source(h0);
|
||
Vertex_index v2 = mesh.source(h1);
|
||
Vertex_index v3 = mesh.source(h2);
|
||
|
||
Edge_index e12 = mesh.edge(h0);
|
||
Edge_index e23 = mesh.edge(h1);
|
||
Edge_index e31 = mesh.edge(h2);
|
||
|
||
double u1 = id_detail::dof_val(m.v_idx[v1], x);
|
||
double u2 = id_detail::dof_val(m.v_idx[v2], x);
|
||
double u3 = id_detail::dof_val(m.v_idx[v3], x);
|
||
|
||
double l12sq = id_detail::edge_length_squared(u1, u2, m.I_e[e12]);
|
||
double l23sq = id_detail::edge_length_squared(u2, u3, m.I_e[e23]);
|
||
double l31sq = id_detail::edge_length_squared(u3, u1, m.I_e[e31]);
|
||
|
||
// A non-real circle configuration (ℓ² ≤ 0) has no limiting angle — skip it.
|
||
if (l12sq <= 0 || l23sq <= 0 || l31sq <= 0) continue;
|
||
|
||
// euclidean_angles expects 2·log(ℓ) per edge — feed log(ℓ²).
|
||
// For a triangle-inequality-violating face euclidean_angles returns the
|
||
// *limiting* angles (π opposite the over-long edge, 0/0 otherwise) with
|
||
// valid=false. We deliberately do NOT skip on !fa.valid: using those
|
||
// limiting angles is the convex C¹ extension onto the infeasible region,
|
||
// so Newton can pass through a flip instead of stalling (Finding 9 —
|
||
// mirrors the Euclidean/Spherical fix in Finding 1). The angles come
|
||
// from genuine Euclidean side lengths, so the extension is geometric.
|
||
auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq));
|
||
|
||
h_alpha[id_detail::hidx(h0)] = fa.alpha3;
|
||
h_alpha[id_detail::hidx(h1)] = fa.alpha1;
|
||
h_alpha[id_detail::hidx(h2)] = fa.alpha2;
|
||
}
|
||
|
||
// Pass 2 — accumulate vertex gradient.
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = m.v_idx[v];
|
||
if (iv < 0) continue;
|
||
double sum_alpha = 0.0;
|
||
for (auto h : CGAL::halfedges_around_target(v, mesh)) {
|
||
if (mesh.is_border(h)) continue;
|
||
sum_alpha += h_alpha[id_detail::hidx(mesh.prev(h))];
|
||
}
|
||
G[static_cast<std::size_t>(iv)] = m.theta_v[v] - sum_alpha;
|
||
}
|
||
|
||
return G;
|
||
}
|
||
|
||
/// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated
|
||
/// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`).
|
||
inline double inversive_distance_energy(
|
||
const ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const InversiveDistanceMaps& m)
|
||
{
|
||
const double* gl_s = gl10_nodes();
|
||
const double* gl_w = gl10_weights();
|
||
|
||
const std::size_t n = x.size();
|
||
double E = 0.0;
|
||
std::vector<double> tx(n);
|
||
for (int k = 0; k < 10; ++k) {
|
||
double t = (1.0 + gl_s[k]) * 0.5;
|
||
double wt = gl_w[k] * 0.5;
|
||
for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i];
|
||
auto G = inversive_distance_gradient(mesh, tx, m);
|
||
double dot = 0.0;
|
||
for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i];
|
||
E += wt * dot;
|
||
}
|
||
return E;
|
||
}
|
||
|
||
/// FD gradient check for the Inversive-Distance functional (central diff).
|
||
/// Uses the same **relative** error criterion as every other gradient check:
|
||
/// `|analytic − fd| / max(1, |analytic|) < tol`.
|
||
inline bool gradient_check_inversive_distance(
|
||
const ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const InversiveDistanceMaps& m,
|
||
double eps = 1e-5,
|
||
double tol = 1e-4)
|
||
{
|
||
auto G = inversive_distance_gradient(mesh, x, m);
|
||
const std::size_t n = G.size();
|
||
bool ok = true;
|
||
|
||
for (std::size_t i = 0; i < n; ++i) {
|
||
std::vector<double> xp = x, xm = x;
|
||
xp[i] += eps;
|
||
xm[i] -= eps;
|
||
double Ep = inversive_distance_energy(mesh, xp, m);
|
||
double Em = inversive_distance_energy(mesh, xm, m);
|
||
double fd = (Ep - Em) / (2.0 * eps);
|
||
double err = std::abs(G[i] - fd);
|
||
double scale = std::max(1.0, std::abs(G[i]));
|
||
if (err / scale > tol) {
|
||
std::cerr << "[inversive-distance] FD gradient mismatch at DOF " << i
|
||
<< ": analytic=" << G[i]
|
||
<< " FD=" << fd
|
||
<< " rel-err=" << (err / scale) << "\n";
|
||
ok = false;
|
||
}
|
||
}
|
||
return ok;
|
||
}
|
||
|
||
/// Newton equilibrium check: returns `true` iff the gradient at `x`
|
||
/// is below `tol` in infinity norm (Σ adj-face angles equal Θ_v).
|
||
inline bool is_inversive_distance_equilibrium(
|
||
const ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const InversiveDistanceMaps& m,
|
||
double tol = 1e-8)
|
||
{
|
||
auto G = inversive_distance_gradient(mesh, x, m);
|
||
for (double g : G)
|
||
if (std::abs(g) > tol) return false;
|
||
return true;
|
||
}
|
||
|
||
} // namespace conformallab
|