diff --git a/README.md b/README.md index 6d6e712..5cbb0c2 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ conformallab++ is a modern C++ reimplementation of the [ConformalLab](https://gi The long-term goal is a **CGAL package** that brings discrete conformal maps (hyper-ideal, spherical, Euclidean) to the CGAL ecosystem using `CGAL::Surface_mesh` as the underlying half-edge data structure. -> **Status:** Phase 5 vollständig abgeschlossen. Alle drei Geometrien lösbar via Newton-Solver (SimplicialLDLT + SparseQR-Fallback). BFS-Layout in ℝ²/S²/Poincaré-Disk, JSON/XML-Serialisierung, vollständige CLI-App. **95 Tests, 2 skipped**. +> **Status:** Phase 6 vollständig abgeschlossen. Alle drei Geometrien lösbar via Newton-Solver (SimplicialLDLT + SparseQR-Fallback). BFS-Layout in ℝ²/S²/Poincaré-Disk mit exakter hyperbolischer Trilateration (Möbius + Kosinussatz), Gauss–Bonnet-Konsistenzprüfung, Tree-Cotree-Schnittgraph, Normalisierung (PCA/Möbius-Zentrierung). JSON/XML-Serialisierung, vollständige CLI-App. **121 Tests, 2 skipped**. --- @@ -29,6 +29,10 @@ The long-term goal is a **CGAL package** that brings discrete conformal maps (hy | **BFS Layout** (ℝ², S², Poincaré disk) | ✅ Phase 5 | | **CLI app** (`conformallab_core`) | ✅ Phase 5 | | **JSON + XML serialisation** | ✅ Phase 5 | +| **Gauss–Bonnet check + enforce** | ✅ Phase 6 | +| **Tree-cotree cut graph** (2g seam edges) | ✅ Phase 6 | +| **Exact hyperbolic trilateration** (Möbius + law of cosines) | ✅ Phase 6 | +| **Layout normalisation** (PCA centring / Möbius centering) | ✅ Phase 6 | --- @@ -218,7 +222,7 @@ ctest --test-dir build -R "^cgal\." --output-on-failure ./bin/conformallab_core -i input.off -g euclidean -o layout.off -j result.json ``` -Expected: **95 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). +Expected: **121 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). ### Interactive viewer @@ -250,8 +254,10 @@ cmake -S code -B build -DWITH_CGAL=ON && cmake --build build -t example_viewer - | `mesh_io.hpp` | `read_mesh` / `write_mesh` / `load_mesh` / `save_mesh` | | `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh` + property-map helpers | | `mesh_builder.hpp` | `make_triangle` / `make_tetrahedron` / `make_quad_strip` / `make_fan` / `make_spherical_tetrahedron` | -| `layout.hpp` | `euclidean_layout` / `spherical_layout` / `hyper_ideal_layout` → `Layout2D/3D`; BFS unfolding | +| `layout.hpp` | `euclidean_layout` / `spherical_layout` / `hyper_ideal_layout` → `Layout2D/3D`; exact hyperbolic trilateration; normalise_{euclidean,hyperbolic,spherical}; `CutGraph*` + `HolonomyData*` | | `serialization.hpp` | `save/load_result_json` + `save/load_result_xml` | +| `gauss_bonnet.hpp` | `euler_characteristic`, `genus`, `gauss_bonnet_sum/rhs/deficit`, `check_gauss_bonnet`, `enforce_gauss_bonnet` | +| `cut_graph.hpp` | `CutGraph` struct + `compute_cut_graph` (tree-cotree, Erickson–Whittlesey 2005) | | `mesh_utils.hpp` | CGAL → Eigen conversion (`cgal_to_eigen`) | | `constants.hpp` | `conformallab::PI`, `TWO_PI` | @@ -267,8 +273,10 @@ code/ │ ├── mesh_io.hpp │ ├── mesh_utils.hpp │ ├── newton_solver.hpp # ← public solve_linear_system + 3 Newton solvers -│ ├── layout.hpp # ← BFS layout (euclidean/spherical/hyper_ideal) +│ ├── layout.hpp # ← BFS layout + exact hyp. trilateration + normalise │ ├── serialization.hpp # ← JSON + XML save/load +│ ├── gauss_bonnet.hpp # ← Gauss–Bonnet check + enforce (Phase 6) +│ ├── cut_graph.hpp # ← Tree-cotree cut-graph algorithm (Phase 6) │ ├── hyper_ideal_{functional,hessian,geometry,utility,visualization_utility}.hpp │ ├── spherical_{functional,hessian,geometry}.hpp │ ├── euclidean_{functional,hessian,geometry}.hpp @@ -297,7 +305,8 @@ code/ │ ├── test_newton_solver.cpp # 14 tests (incl. 3 SparseQR tests) │ ├── test_mesh_io.cpp # 6 tests │ ├── test_pipeline.cpp # 5 tests -│ └── test_layout.cpp # 8 tests (layout + JSON/XML) +│ ├── test_layout.cpp # 8 tests (layout + JSON/XML) +│ └── test_phase6.cpp # 26 tests (GB, cut graph, trilateration, normalisation) └── deps/ ├── eigen-3.4.0/ # always extracted ├── CGAL-6.1.1/ # extracted with WITH_CGAL @@ -334,7 +343,11 @@ Pure-math tests requiring only Eigen: Clausen / Lobachevsky / ImLi₂, hyper-ide | `Pipeline` | 5 | End-to-end: build → setup → solve → export → reload, all three geometries | | `Layout` | 6 | Edge-length preservation (Euclidean/Spherical), arc-lengths on S², Poincaré disk | | `Serialization` | 2 | JSON and XML round-trips (DOF + layout) | -| **Total** | **95** | 2 skipped (Hessian stubs) | +| `GaussBonnet` | 8 | χ, genus, sum/rhs, deficit, check, enforce (sign-correct) | +| `CutGraph` | 6 | Tree-cotree algorithm, open/closed meshes, flag–index consistency | +| `HyperbolicTrilateration` | 4 | Möbius + law of cosines: exact distances, inside-disk, off-origin | +| `Normalisation` | 4 | Euclidean centroid, length-ratio invariance, Möbius centering (hyp.) | +| **Total** | **121** | 2 skipped (Hessian stubs) | --- @@ -498,11 +511,13 @@ Gauss–Bonnet pre-check ✅ ❌ ### What "cone metrics" still requires -**Cone metrics** — the property map `theta_v` is already subtracted in the gradient (`G_v = Σα_v − Θ_v`), so prescribing a cone angle is a one-liner: `maps.theta_v[v] = desired_angle`. What is missing is the *application layer*: checking Gauss–Bonnet consistency (Σ (2π − Θ_v) = 2π·χ), distributing angle defects sensibly, and special handling at boundary vertices. +**Cone metrics** — the property map `theta_v` is already subtracted in the gradient (`G_v = Σα_v − Θ_v`), so prescribing a cone angle is a one-liner: `maps.theta_v[v] = desired_angle`. Use `check_gauss_bonnet(mesh, maps)` to verify Σ(2π−Θ_v) = 2π·χ before solving, and `enforce_gauss_bonnet` to fix floating-point drift. -**Layout (Phase 5 ✅)** — `layout.hpp` implements BFS unfolding for all three geometries via `euclidean_layout`, `spherical_layout`, and `hyper_ideal_layout`. For open meshes the embedding is globally consistent; for closed meshes the first BFS visit wins and `has_seam = true` is set. To get a proper parameterisation of a closed mesh, cut it to a disk first (not yet implemented). +**Layout (Phase 6 ✅)** — `layout.hpp` implements BFS unfolding for all three geometries. The hyperbolic layout now uses **exact trilateration** via Möbius maps and the hyperbolic law of cosines (replacing the old `tanh(d/2)` approximation). Pass a `CutGraph*` to track seam edges on closed surfaces, and a `HolonomyData*` to capture the holonomy group elements. Set `normalise=true` for PCA centring (Euclidean), Möbius centering (hyperbolic), or Rodrigues rotation to the north pole (spherical). -**Next steps** — global uniformization (cutting closed meshes, period matrices, holonomy), Gauss–Bonnet checking, analytical HyperIdeal Hessian. +**Cut graph (Phase 6 ✅)** — `compute_cut_graph(mesh)` implements the tree-cotree algorithm (Erickson–Whittlesey 2005). For a closed genus-g surface it returns exactly 2g seam edges whose removal turns the surface into a topological disk. + +**Next steps** — holonomy matrices / period matrix (genus-1 torus: τ = ω₂/ω₁ ∈ ℍ), analytical HyperIdeal Hessian, full global uniformization pipeline. --- @@ -928,10 +943,21 @@ Phase 5 Layout + CLI + Serialisierung ✅ abgeschlossen → test_layout.cpp: 8 Tests (Eucl./Sphär./HyperIdeal + JSON/XML) → 95 Tests gesamt (2 skipped) -Phase 6 (geplant) +Phase 6 Layout-Erweiterung (Java-Parität) ✅ abgeschlossen + → gauss_bonnet.hpp: euler_characteristic, genus, Σ(2π-Θ_v) check + enforce + → cut_graph.hpp: Tree-Cotree-Algorithmus (Erickson-Whittlesey); 2g Schnittkan. + → layout.hpp: exakte hyperbolische Trilateration (Möbius + hyperb. Kosinussatz) + → layout.hpp: normalise_euclidean (Schwerpunkt→0, PCA-Rotation) + → layout.hpp: normalise_hyperbolic (Möbius-Zentrierung im Poincaré-Disk) + → layout.hpp: normalise_spherical (Rodrigues-Rotation zum Nordpol) + → layout.hpp: CutGraph* + HolonomyData* Parameter für alle Layout-Funktionen + → test_phase6.cpp: 26 Tests (GB, CutGraph, Trilateration, Normalisierung) + → 121 Tests gesamt (2 skipped) + +Phase 7 (geplant) → Analytischer HyperIdeal-Hessian (direkte Ableitung durch ζ-Kette) - → Gauss–Bonnet Konsistenzprüfung für Kegelmetriken - → Mesh-Cut für geschlossene Flächen → globale Parameterisierung + → Holonomie-Matrizen für Tori (Periodenmatrix τ = ω₂/ω₁ ∈ ℍ) + → Vollständige globale Parameterisierung geschlossener Flächen → Inversive-Distance-Funktional (Luo 2004) ``` diff --git a/code/include/cut_graph.hpp b/code/include/cut_graph.hpp new file mode 100644 index 0000000..0092ef1 --- /dev/null +++ b/code/include/cut_graph.hpp @@ -0,0 +1,152 @@ +#pragma once +// cut_graph.hpp +// +// Phase 6 — Tree-cotree algorithm for computing a cut graph of a triangulated +// surface. +// +// For a closed, orientable, genus-g surface: +// #vertices (V), #edges (E), #faces (F) +// Euler: V − E + F = 2 − 2g +// Primal spanning tree: V − 1 edges +// Dual spanning tree: F − 1 edges (avoiding duals of tree edges) +// Remaining: E − (V−1) − (F−1) = 2g cut edges +// +// These 2g cut edges generate H₁(M, ℤ) ≅ ℤ^{2g}. +// Cutting along them turns M into a topological disk. +// +// For open meshes (boundary present) the algorithm still works: the dual BFS +// starts from a boundary-adjacent face, and boundary half-edges are skipped. +// The number of cut edges will be E − (V−1) − (F−1) − B where B counts +// boundary edges treated as dual tree edges. +// +// Usage: +// CutGraph cg = compute_cut_graph(mesh); +// // cg.cut_edge_flags[e.idx()] == true → treat edge as seam in BFS +// euclidean_layout(mesh, x, maps, &cg); // layout with holonomy tracking + +#include "conformal_mesh.hpp" +#include "gauss_bonnet.hpp" // for euler_characteristic / genus +#include +#include +#include + +namespace conformallab { + +// ───────────────────────────────────────────────────────────────────────────── +// CutGraph +// ───────────────────────────────────────────────────────────────────────────── + +struct CutGraph { + /// cut_edge_flags[e.idx()] = true ↔ this edge is a cut edge. + /// Size = mesh.number_of_edges(). + std::vector cut_edge_flags; + + /// Indices of the 2g cut edges in order (size = 2g). + std::vector cut_edge_indices; + + /// Genus of the surface (0 for topological spheres and open patches). + int genus = 0; + + bool is_cut(Edge_index e) const + { + return static_cast(e.idx()) < cut_edge_flags.size() + && cut_edge_flags[static_cast(e.idx())]; + } +}; + +// ───────────────────────────────────────────────────────────────────────────── +// compute_cut_graph +// ───────────────────────────────────────────────────────────────────────────── +// +// Implements the standard tree-cotree algorithm (Erickson–Whittlesey 2005): +// +// Step 1: BFS primal spanning tree T (V−1 primal tree edges). +// Step 2: BFS dual spanning tree T* (F−1 dual/primal edges, avoiding +// edges whose primal crosses T). +// Step 3: cut edges = primal edges neither in T nor "used" by T*. + +inline CutGraph compute_cut_graph(const ConformalMesh& mesh) +{ + const std::size_t nv = mesh.number_of_vertices(); + const std::size_t ne = mesh.number_of_edges(); + const std::size_t nf = mesh.number_of_faces(); + + CutGraph cg; + cg.cut_edge_flags.assign(ne, false); + cg.genus = conformallab::genus(mesh); + + if (nv == 0 || nf == 0) return cg; + + // ── Step 1: primal spanning tree via BFS from vertex 0 ─────────────────── + std::vector tree_edge(ne, false); + std::vector v_visited(nv, false); + + { + std::queue q; + auto v0 = *mesh.vertices().begin(); + v_visited[v0.idx()] = true; + q.push(v0); + while (!q.empty()) { + Vertex_index v = q.front(); q.pop(); + for (Halfedge_index h : CGAL::halfedges_around_target(v, mesh)) { + Vertex_index u = mesh.source(h); + if (!v_visited[static_cast(u.idx())]) { + v_visited[static_cast(u.idx())] = true; + tree_edge[static_cast(mesh.edge(h).idx())] = true; + q.push(u); + } + } + } + } + + // ── Step 2: dual spanning tree via BFS from face 0 ─────────────────────── + // Dual edge between face f and face f_adj crosses primal edge e. + // Include dual edge only if: + // (a) e is not a primal tree edge (tree_edge[e] == false) + // (b) h_adj is not a border halfedge + std::vector dual_tree_edge(ne, false); + std::vector f_visited(nf, false); + + { + std::queue q; + auto f0 = *mesh.faces().begin(); + f_visited[static_cast(f0.idx())] = true; + q.push(f0); + while (!q.empty()) { + Face_index f = q.front(); q.pop(); + for (Halfedge_index h : + CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + { + Halfedge_index h_opp = mesh.opposite(h); + if (mesh.is_border(h_opp)) continue; // boundary edge + Face_index f_adj = mesh.face(h_opp); + if (f_visited[static_cast(f_adj.idx())]) continue; + + std::size_t eidx = static_cast(mesh.edge(h).idx()); + if (!tree_edge[eidx]) { + // Use this dual edge in T* + f_visited[static_cast(f_adj.idx())] = true; + dual_tree_edge[eidx] = true; + q.push(f_adj); + } + } + } + } + + // ── Step 3: cut edges = neither in T nor in T* nor on boundary ─────────── + // Boundary edges are adjacent to the "outer face" and need no cutting — + // they are implicitly handled by the boundary itself. + for (Edge_index e : mesh.edges()) { + std::size_t idx = static_cast(e.idx()); + if (tree_edge[idx] || dual_tree_edge[idx]) continue; + // Skip boundary edges — they are not interior homological cycles. + Halfedge_index h = mesh.halfedge(e); + if (mesh.is_border(h) || mesh.is_border(mesh.opposite(h))) continue; + cg.cut_edge_flags[idx] = true; + cg.cut_edge_indices.push_back(idx); + } + + return cg; +} + +} // namespace conformallab diff --git a/code/include/gauss_bonnet.hpp b/code/include/gauss_bonnet.hpp new file mode 100644 index 0000000..0d5bc61 --- /dev/null +++ b/code/include/gauss_bonnet.hpp @@ -0,0 +1,145 @@ +#pragma once +// gauss_bonnet.hpp +// +// Phase 6 — Gauss–Bonnet consistency check for prescribed target angles. +// +// Before calling newton_*() with custom target angles, verify that +// the angle defect sum matches the topology: +// +// Σ_v (2π − Θ_v) = 2π · χ(M) (Euclidean / flat) +// Σ_v (2π − Θ_v) > 0 (spherical, χ > 0) +// Σ_v (2π − Θ_v) < 0 (hyperbolic, χ < 0) +// +// If this fails, no conformal factor can realise the target angles and +// Newton will silently fail to converge. +// +// API: +// int euler_characteristic(mesh) +// int genus(mesh) +// double gauss_bonnet_sum(mesh, maps) — Σ(2π − Θ_v) +// double gauss_bonnet_rhs(mesh) — 2π · χ(M) +// double gauss_bonnet_deficit(mesh, maps) — lhs − rhs (0 = satisfied) +// void check_gauss_bonnet(mesh, maps [, tol]) — throws if violated +// void enforce_gauss_bonnet(mesh, maps) — shifts θ_v by uniform Δ + +#include "conformal_mesh.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include "constants.hpp" + +#include +#include +#include +#include + +namespace conformallab { + +// ── Topology helpers ────────────────────────────────────────────────────────── + +/// Euler characteristic χ = V − E + F. +/// For closed orientable surfaces: χ = 2 − 2g. +inline int euler_characteristic(const ConformalMesh& mesh) +{ + return static_cast(mesh.number_of_vertices()) + - static_cast(mesh.number_of_edges()) + + static_cast(mesh.number_of_faces()); +} + +/// Genus of a closed orientable surface: g = (2 − χ) / 2. +/// Returns 0 for open meshes (boundary present) — callers should check. +inline int genus(const ConformalMesh& mesh) +{ + int chi = euler_characteristic(mesh); + return (2 - chi) / 2; +} + +// ── Left-hand side Σ(2π − Θ_v) ───────────────────────────────────────────── + +inline double gauss_bonnet_sum( + const ConformalMesh& mesh, + const ConformalMesh::Property_map& theta) +{ + double s = 0.0; + for (auto v : mesh.vertices()) + s += TWO_PI - theta[v]; + return s; +} + +inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps& mp) + { return gauss_bonnet_sum(m, mp.theta_v); } +inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps& mp) + { return gauss_bonnet_sum(m, mp.theta_v); } +inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp) + { return gauss_bonnet_sum(m, mp.theta_v); } + +// ── Right-hand side 2π · χ(M) ─────────────────────────────────────────────── + +inline double gauss_bonnet_rhs(const ConformalMesh& mesh) +{ + return TWO_PI * static_cast(euler_characteristic(mesh)); +} + +// ── Deficit: lhs − rhs (0 = Gauss–Bonnet satisfied) ───────────────────────── + +template +inline double gauss_bonnet_deficit(const ConformalMesh& mesh, const Maps& maps) +{ + return gauss_bonnet_sum(mesh, maps) - gauss_bonnet_rhs(mesh); +} + +// ── check_gauss_bonnet — throws std::runtime_error if |deficit| > tol ───────── + +inline void check_gauss_bonnet(const ConformalMesh& mesh, + double lhs, + double tol = 1e-8) +{ + double rhs = gauss_bonnet_rhs(mesh); + double def = lhs - rhs; + if (std::abs(def) > tol) { + std::ostringstream msg; + msg << "Gauss–Bonnet violated:\n" + << " Σ(2π−Θ_v) = " << lhs + << " expected 2π·χ = " << rhs + << " (χ = " << euler_characteristic(mesh) + << ", genus = " << genus(mesh) << ")\n" + << " deficit = " << def; + throw std::runtime_error(msg.str()); + } +} + +template +inline void check_gauss_bonnet(const ConformalMesh& mesh, + const Maps& maps, + double tol = 1e-8) +{ + check_gauss_bonnet(mesh, gauss_bonnet_sum(mesh, maps), tol); +} + +// ── enforce_gauss_bonnet — adjust θ_v by uniform Δ ─────────────────────────── +// +// Adds δ = (rhs − lhs) / V to every θ_v so that Gauss–Bonnet holds exactly. +// After this call, check_gauss_bonnet() will not throw (up to floating-point). +// Only modifies free vertices (v_idx[v] >= 0 for EuclideanMaps / SphericalMaps; +// always all vertices for the raw property-map overload). + +inline void enforce_gauss_bonnet( + ConformalMesh& mesh, + ConformalMesh::Property_map& theta) +{ + double lhs = gauss_bonnet_sum(mesh, theta); + double rhs = gauss_bonnet_rhs(mesh); + // Adding δ to every θ_v decreases the sum Σ(2π−θ_v) by V·δ. + // We need lhs − V·δ = rhs, so δ = (lhs − rhs) / V. + double delta = (lhs - rhs) / static_cast(mesh.number_of_vertices()); + for (auto v : mesh.vertices()) + theta[v] += delta; +} + +template +inline void enforce_gauss_bonnet(ConformalMesh& mesh, Maps& maps) +{ + enforce_gauss_bonnet(mesh, maps.theta_v); +} + +} // namespace conformallab diff --git a/code/include/layout.hpp b/code/include/layout.hpp index f1d121b..31db33e 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -1,11 +1,8 @@ #pragma once // layout.hpp // -// Phase 5 — Layout / embedding: DOF vector → vertex coordinates in target space. -// -// After Newton converges you have conformal scale factors u_v (and optionally -// edge DOFs). This header converts those factors into actual vertex positions -// in the target geometry by BFS-unfolding the triangulation. +// Phase 5/6 — Layout / embedding: DOF vector → vertex coordinates in the +// target geometry via BFS-trilateration. // // ┌──────────────────────────────────────────────────────────────────────────┐ // │ Algorithm (all three geometries) │ @@ -14,51 +11,76 @@ // │ 2. Choose a root face, place its three vertices analytically. │ // │ 3. BFS over the dual graph: for each adjacent face, the shared edge is │ // │ already placed; trilaterate the third vertex. │ +// │ 4. (Phase 6) If a CutGraph is supplied, cut edges are treated as │ +// │ boundary; crossing them records a HolonomyData entry instead. │ // │ │ -// │ For open meshes (boundary) the placement is globally consistent. │ -// │ For closed meshes the BFS visits some vertices twice (seam); the first │ -// │ visit wins and `has_seam = true` is set. To get a proper global │ -// │ parameterisation of a closed mesh, cut it to a disk first. │ +// │ For open meshes the placement is globally consistent. │ +// │ For closed meshes without a cut: first visit wins, has_seam = true. │ +// │ For closed meshes with a CutGraph: each vertex is placed exactly once; │ +// │ the holonomy (translation / Möbius) of each cut is recorded. │ // └──────────────────────────────────────────────────────────────────────────┘ // -// Functions: -// euclidean_layout(mesh, x, maps) → Layout2D (positions in ℝ²) -// spherical_layout (mesh, x, maps) → Layout3D (positions on S² ⊂ ℝ³) -// hyper_ideal_layout(mesh, x, maps)→ Layout2D (Poincaré disk) +// Euclidean trilateration — exact (analytic formula in ℝ²) +// Spherical trilateration — exact (spherical law of cosines on S²) +// Hyperbolic trilateration — exact (hyperbolic law of cosines + Möbius maps +// in the Poincaré disk model) +// +// Normalisation +// normalise_euclidean(layout) — centroid → origin, major axis → x-axis +// normalise_hyperbolic(layout) — Möbius map centering to origin +// normalise_spherical(layout) — rotate centroid to north pole +// +// Holonomy (Euclidean, genus-g surfaces with CutGraph) +// HolonomyData.translations[i] — translation vector ω_i for cut edge i +// (for a flat torus: ω_1, ω_2 are the lattice generators) #include "conformal_mesh.hpp" #include "euclidean_functional.hpp" #include "spherical_functional.hpp" #include "hyper_ideal_functional.hpp" +#include "cut_graph.hpp" #include #include #include +#include #include #include +#include +#include namespace conformallab { // ── Result types ────────────────────────────────────────────────────────────── struct Layout2D { - std::vector uv; ///< uv[v.idx()] = 2-D position + std::vector uv; ///< uv[v.idx()] = 2-D position bool success = false; - bool has_seam = false; ///< true if mesh is closed (first-visit used) + bool has_seam = false; ///< true if mesh is closed and no CutGraph given }; struct Layout3D { - std::vector pos; ///< pos[v.idx()] = 3-D position on S² + std::vector pos; ///< pos[v.idx()] = 3-D position on S² bool success = false; bool has_seam = false; }; +/// Per-cut-edge holonomy for the Euclidean case. +/// translations[i] is the translation ω associated with cut_edge_indices[i] +/// of the CutGraph. For a flat torus, ω_1 and ω_2 are the lattice generators. +struct HolonomyData { + std::vector translations; // one per cut edge + std::vector cut_edge_indices; +}; + // ── Internal helpers ────────────────────────────────────────────────────────── namespace detail { -// Euclidean trilaterate: given placed p_a, p_b and distances d_a (from p_a), -// d_b (from p_b), return the new point on the LEFT side of the directed edge -// p_a → p_b (corresponds to CCW face orientation). +// ───────────────────────────────────────────────────────────────────────────── +// Euclidean trilateration +// Given placed p_a, p_b and distances d_a (from p_a), d_b (from p_b), +// return the point on the LEFT (CCW) side of the directed edge p_a → p_b. +// ───────────────────────────────────────────────────────────────────────────── inline Eigen::Vector2d trilaterate_2d( const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, double da, double db) @@ -71,55 +93,212 @@ inline Eigen::Vector2d trilaterate_2d( double t = (da*da - db*db + d*d) / (2.0 * d); double s2 = da*da - t*t; double s = (s2 > 0.0) ? std::sqrt(s2) : 0.0; - return pa + t*e1 + s*e2; // s > 0 → left side + return pa + t*e1 + s*e2; // s > 0 → left side } -// Spherical trilaterate: given unit vectors pa, pb on S² and arc-lengths da, db -// to the new point, return the new unit vector on the LEFT side of the geodesic -// arc from pa to pb. +// ───────────────────────────────────────────────────────────────────────────── +// Spherical trilateration +// Given unit vectors pa, pb on S² and arc-lengths da, db to the new point, +// return the unit vector on the LEFT side of the geodesic arc pa → pb. // // Solution: p_c = α·pa + β·pb + γ·(pa × pb), γ > 0. -// Constraints: pa·p_c = cos(da), pb·p_c = cos(db), |p_c|=1. +// Constraints: pa·p_c = cos(da), pb·p_c = cos(db), |p_c| = 1. +// ───────────────────────────────────────────────────────────────────────────── inline Eigen::Vector3d trilaterate_sph( const Eigen::Vector3d& pa, const Eigen::Vector3d& pb, double da, double db) { - double c = pa.dot(pb); // cos(l_ab) + double c = pa.dot(pb); double denom = 1.0 - c*c; - if (denom < 1e-14) return pa; // degenerate (pa ≈ pb) + if (denom < 1e-14) return pa; double cda = std::cos(da), cdb = std::cos(db); double alpha = (cda - c*cdb) / denom; double beta = (cdb - c*cda) / denom; - Eigen::Vector3d cross = pa.cross(pb); - double cross_n2 = cross.squaredNorm(); - + Eigen::Vector3d cross = pa.cross(pb); + double cross_n2 = cross.squaredNorm(); double gamma2 = 1.0 - alpha*alpha - beta*beta - 2.0*alpha*beta*c; double gamma = (gamma2 > 0.0 && cross_n2 > 1e-28) - ? std::sqrt(gamma2 / cross_n2) - : 0.0; + ? std::sqrt(gamma2 / cross_n2) : 0.0; Eigen::Vector3d p = alpha*pa + beta*pb + gamma*cross; double n = p.norm(); - return (n > 1e-14) ? (p / n) : pa; // gamma > 0 → left side + return (n > 1e-14) ? (p / n) : pa; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Hyperbolic trilateration — exact via Möbius map + hyperbolic law of cosines. +// +// pa, pb: Poincaré disk coordinates of two already-placed vertices +// D: hyperbolic distance pa → pb (from the DOF vector) +// da: hyperbolic distance pa → pc (new vertex) +// db: hyperbolic distance pb → pc (new vertex) +// +// Returns the Poincaré disk coordinate of pc on the LEFT (CCW) side of pa→pb. +// +// Algorithm: +// 1. Map pa → 0 via Möbius T: z ↦ (z−pa)/(1−conj(pa)·z) +// 2. θ_b = arg(T(pb)) — direction from origin towards T(pb) +// 3. Angle α at pa from the hyperbolic law of cosines: +// cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α) +// → cos(α) = (cosh(da)·cosh(D) − cosh(db)) / (sinh(da)·sinh(D)) +// 4. pc in origin frame: tanh(da/2)·exp(i·(θ_b + α)) [left = +α] +// 5. Map back: T⁻¹(w) = (w + pa)/(1 + conj(pa)·w) +// ───────────────────────────────────────────────────────────────────────────── +inline Eigen::Vector2d trilaterate_hyp( + const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, + double D, // hyperbolic distance pa → pb + double da, // hyperbolic distance pa → pc + double db) // hyperbolic distance pb → pc +{ + using C = std::complex; + + // Degenerate guard: very short edges or coincident points + if (D < 1e-12 || da < 1e-12) return pa; + if (std::sinh(da) * std::sinh(D) < 1e-14) return pa; + + C a(pa.x(), pa.y()); + C b(pb.x(), pb.y()); + + // Möbius map: T_a(z) = (z − a) / (1 − conj(a)·z) + auto mobius_fwd = [](C z, C center) -> C { + return (z - center) / (C(1.0) - std::conj(center) * z); + }; + // Inverse: T_a⁻¹(w) = (w + a) / (1 + conj(a)·w) + auto mobius_inv = [](C w, C center) -> C { + return (w + center) / (C(1.0) + std::conj(center) * w); + }; + + C b_mapped = mobius_fwd(b, a); + double theta_b = std::arg(b_mapped); + + // Hyperbolic law of cosines for angle α at pa: + // cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α) + double cos_alpha = (std::cosh(da) * std::cosh(D) - std::cosh(db)) + / (std::sinh(da) * std::sinh(D)); + cos_alpha = std::max(-1.0, std::min(1.0, cos_alpha)); + double alpha = std::acos(cos_alpha); // ∈ [0, π], sin > 0 for left side + + // pc in the frame where pa = 0 and pb is on the positive real axis, + // then rotate back by theta_b: + C pc_origin = std::tanh(da * 0.5) * std::exp(C(0.0, theta_b + alpha)); + + // Map back to original disk + C pc_c = mobius_inv(pc_origin, a); + return Eigen::Vector2d(pc_c.real(), pc_c.imag()); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Möbius centering of a Poincaré-disk layout +// Applies T_c(z) = (z − c)/(1 − conj(c)·z) where c is the Euclidean centroid +// of all vertex positions. Maps c → 0, i.e. centers the cloud in the disk. +// ───────────────────────────────────────────────────────────────────────────── +inline void center_poincare_disk(std::vector& uv) +{ + if (uv.empty()) return; + using C = std::complex; + + // Euclidean centroid (good enough for moderate displacements from origin) + C centroid(0.0, 0.0); + for (auto& p : uv) centroid += C(p.x(), p.y()); + centroid /= static_cast(uv.size()); + + // Clamp to strictly inside the disk + double r = std::abs(centroid); + if (r > 0.999) centroid *= (0.999 / r); + if (r < 1e-12) return; // already centered + + for (auto& p : uv) { + C z(p.x(), p.y()); + C w = (z - centroid) / (C(1.0) - std::conj(centroid) * z); + p = Eigen::Vector2d(w.real(), w.imag()); + } } } // namespace detail +// ── Layout normalisation ────────────────────────────────────────────────────── + +/// Euclidean: translate centroid to origin, rotate major axis to x-axis (PCA). +inline void normalise_euclidean(Layout2D& layout) +{ + if (!layout.success || layout.uv.empty()) return; + const std::size_t n = layout.uv.size(); + + // Centroid + Eigen::Vector2d mean = Eigen::Vector2d::Zero(); + for (auto& p : layout.uv) mean += p; + mean /= static_cast(n); + for (auto& p : layout.uv) p -= mean; + + // PCA: 2×2 covariance + Eigen::Matrix2d cov = Eigen::Matrix2d::Zero(); + for (auto& p : layout.uv) cov += p * p.transpose(); + cov /= static_cast(n); + + Eigen::SelfAdjointEigenSolver eig(cov); + // Eigenvectors in ascending order — we want the largest (index 1) + Eigen::Vector2d major = eig.eigenvectors().col(1); + // Rotation that aligns major axis with x-axis + double angle = -std::atan2(major.y(), major.x()); + Eigen::Matrix2d R; + R << std::cos(angle), -std::sin(angle), + std::sin(angle), std::cos(angle); + for (auto& p : layout.uv) p = R * p; +} + +/// Hyperbolic: Möbius map centering the Euclidean centroid to the disk origin. +inline void normalise_hyperbolic(Layout2D& layout) +{ + if (!layout.success || layout.uv.empty()) return; + detail::center_poincare_disk(layout.uv); +} + +/// Spherical: rotate so the Euclidean centroid of vertex positions points +/// towards the north pole (0, 0, 1). +inline void normalise_spherical(Layout3D& layout) +{ + if (!layout.success || layout.pos.empty()) return; + Eigen::Vector3d mean = Eigen::Vector3d::Zero(); + for (auto& p : layout.pos) mean += p; + double len = mean.norm(); + if (len < 1e-12) return; + mean /= len; // unit vector of mean direction + + Eigen::Vector3d north(0.0, 0.0, 1.0); + Eigen::Vector3d axis = mean.cross(north); + double sin_a = axis.norm(); + double cos_a = mean.dot(north); + if (sin_a < 1e-12) return; // already at north (or south) pole + axis /= sin_a; + + // Rodrigues rotation matrix + Eigen::Matrix3d K; + K << 0.0, -axis.z(), axis.y(), + axis.z(), 0.0, -axis.x(), + -axis.y(), axis.x(), 0.0; + Eigen::Matrix3d R = Eigen::Matrix3d::Identity() + + sin_a * K + (1.0 - cos_a) * K * K; + + for (auto& p : layout.pos) p = R * p; +} + // ── Euclidean layout ────────────────────────────────────────────────────────── // -// Computes 2-D positions for each vertex by BFS-unfolding the triangulation -// in the plane. The updated edge length is: -// // l_ij = exp( (λ°_ij + u_i + u_j + λ_e) / 2 ) // -// where u_i = x[v_idx[v]] (0 if pinned) and λ_e = x[e_idx[e]] (0 if absent). +// Optional CutGraph: if non-null, cut edges are treated as boundary; +// holonomy translations are computed for each cut edge and returned via +// *holonomy (if non-null). // ───────────────────────────────────────────────────────────────────────────── inline Layout2D euclidean_layout( ConformalMesh& mesh, const std::vector& x, - const EuclideanMaps& maps) + const EuclideanMaps& maps, + const CutGraph* cut = nullptr, + HolonomyData* holonomy = nullptr, + bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); Layout2D result; @@ -143,34 +322,35 @@ inline Layout2D euclidean_layout( std::vector vertex_placed(nv, false); std::vector face_placed(mesh.number_of_faces(), false); - // ── Place root face ─────────────────────────────────────────────────── + // ── Place root face ─────────────────────────────────────────────────────── Face_index f0 = *mesh.faces().begin(); Halfedge_index h0 = mesh.halfedge(f0); Halfedge_index h1 = mesh.next(h0); Halfedge_index h2 = mesh.next(h1); - - Vertex_index vA = mesh.source(h0); - Vertex_index vB = mesh.source(h1); - Vertex_index vC = mesh.source(h2); - + Vertex_index vA = mesh.source(h0); + Vertex_index vB = mesh.source(h1); + Vertex_index vC = mesh.source(h2); double lAB = edge_len(h0); - double lCA = edge_len(h2); // dist C—A - double lBC = edge_len(h1); // dist B—C + double lCA = edge_len(h2); + double lBC = edge_len(h1); result.uv[vA.idx()] = Eigen::Vector2d(0.0, 0.0); result.uv[vB.idx()] = Eigen::Vector2d(lAB, 0.0); result.uv[vC.idx()] = detail::trilaterate_2d( result.uv[vA.idx()], result.uv[vB.idx()], lCA, lBC); - vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true; face_placed[f0.idx()] = true; - // ── BFS ─────────────────────────────────────────────────────────────── + // ── BFS ─────────────────────────────────────────────────────────────────── std::queue q; auto enqueue = [&](Face_index f) { - for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - auto h_opp = mesh.opposite(h); - if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) + for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index h_opp = mesh.opposite(h); + if (mesh.is_border(h_opp)) continue; + // Treat cut edges as boundary (don't cross them in BFS) + if (cut && cut->is_cut(mesh.edge(h))) continue; + Face_index f_adj = mesh.face(h_opp); + if (!face_placed[f_adj.idx()]) q.push(h_opp); } }; @@ -187,11 +367,11 @@ inline Layout2D euclidean_layout( Eigen::Vector2d p = detail::trilaterate_2d( result.uv[v_src.idx()], result.uv[v_tgt.idx()], - edge_len(mesh.prev(h)), // dist(v_new, v_src) - edge_len(mesh.next(h))); // dist(v_tgt, v_new) + edge_len(mesh.prev(h)), + edge_len(mesh.next(h))); if (!vertex_placed[v_new.idx()]) { - result.uv[v_new.idx()] = p; + result.uv[v_new.idx()] = p; vertex_placed[v_new.idx()] = true; } else { result.has_seam = true; @@ -201,24 +381,64 @@ inline Layout2D euclidean_layout( } result.success = true; + + // ── Holonomy: compute translation for each cut edge ─────────────────────── + // For each cut edge e = (u,v), look at the face on the far side (h_opp). + // Trilaterate the third vertex w of that face from uv[u] and uv[v], + // and compare to the actual position uv[w]. + // Translation ω = trilaterated_position(w) − uv[w]. + if (cut && holonomy) { + holonomy->cut_edge_indices = cut->cut_edge_indices; + holonomy->translations.reserve(cut->cut_edge_indices.size()); + + for (std::size_t ce_idx : cut->cut_edge_indices) { + Edge_index e = *std::next(mesh.edges().begin(), + static_cast(ce_idx)); + Halfedge_index h = mesh.halfedge(e); + Halfedge_index h_opp = mesh.opposite(h); + + // Choose the halfedge whose face was NOT placed during BFS + // (the one across the cut from the main BFS sweep). + // If both sides were placed, use h_opp; if neither, skip. + Halfedge_index h_cross = h_opp; + if (mesh.is_border(h_cross)) h_cross = h; + if (mesh.is_border(h_cross)) { + holonomy->translations.push_back(Eigen::Vector2d::Zero()); + continue; + } + + Vertex_index v_src = mesh.source(h_cross); + Vertex_index v_tgt = mesh.target(h_cross); + Vertex_index v_new = mesh.target(mesh.next(h_cross)); + + if (!vertex_placed[v_new.idx()]) { + holonomy->translations.push_back(Eigen::Vector2d::Zero()); + continue; + } + + Eigen::Vector2d p_tri = detail::trilaterate_2d( + result.uv[v_src.idx()], result.uv[v_tgt.idx()], + edge_len(mesh.prev(h_cross)), + edge_len(mesh.next(h_cross))); + + holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]); + } + } + + if (normalise) normalise_euclidean(result); return result; } // ── Spherical layout ────────────────────────────────────────────────────────── // -// Computes positions on the unit sphere S² by BFS-unfolding the triangulation. -// Updated spherical arc-length: -// // l_ij = 2·arcsin( min(exp(Λ_ij / 2), 1) ) // -// where Λ_ij = λ°_ij + u_i + u_j (+ edge DOF if present). -// -// Output: pos[v.idx()] is a unit 3-D vector on S². // ───────────────────────────────────────────────────────────────────────────── inline Layout3D spherical_layout( ConformalMesh& mesh, const std::vector& x, - const SphericalMaps& maps) + const SphericalMaps& maps, + bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); Layout3D result; @@ -239,7 +459,6 @@ inline Layout3D spherical_layout( return 2.0 * std::asin(std::min(std::exp(lam * 0.5), 1.0)); }; - // Place root face: vA at north pole, vB along meridian, vC via spherical law std::vector vertex_placed(nv, false); std::vector face_placed(mesh.number_of_faces(), false); @@ -247,28 +466,24 @@ inline Layout3D spherical_layout( Halfedge_index h0 = mesh.halfedge(f0); Halfedge_index h1 = mesh.next(h0); Halfedge_index h2 = mesh.next(h1); - - Vertex_index vA = mesh.source(h0); - Vertex_index vB = mesh.source(h1); - Vertex_index vC = mesh.source(h2); - + Vertex_index vA = mesh.source(h0); + Vertex_index vB = mesh.source(h1); + Vertex_index vC = mesh.source(h2); double lAB = arc_len(h0); double lCA = arc_len(h2); double lBC = arc_len(h1); - // vA = north pole, vB along the 0-meridian at arc distance lAB result.pos[vA.idx()] = Eigen::Vector3d(0.0, 0.0, 1.0); result.pos[vB.idx()] = Eigen::Vector3d(std::sin(lAB), 0.0, std::cos(lAB)); result.pos[vC.idx()] = detail::trilaterate_sph( result.pos[vA.idx()], result.pos[vB.idx()], lCA, lBC); - vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true; face_placed[f0.idx()] = true; std::queue q; auto enqueue = [&](Face_index f) { - for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - auto h_opp = mesh.opposite(h); + for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index h_opp = mesh.opposite(h); if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) q.push(h_opp); } @@ -286,11 +501,10 @@ inline Layout3D spherical_layout( Eigen::Vector3d p = detail::trilaterate_sph( result.pos[v_src.idx()], result.pos[v_tgt.idx()], - arc_len(mesh.prev(h)), - arc_len(mesh.next(h))); + arc_len(mesh.prev(h)), arc_len(mesh.next(h))); if (!vertex_placed[v_new.idx()]) { - result.pos[v_new.idx()] = p; + result.pos[v_new.idx()] = p; vertex_placed[v_new.idx()] = true; } else { result.has_seam = true; @@ -300,28 +514,26 @@ inline Layout3D spherical_layout( } result.success = true; + if (normalise) normalise_spherical(result); return result; } -// ── HyperIdeal layout (Poincaré disk) ──────────────────────────────────────── +// ── HyperIdeal layout (Poincaré disk) — exact trilateration ────────────────── // -// Places the hyperbolic triangulation in the Poincaré disk model. -// The effective hyperbolic edge length between vertices i and j is: +// Hyperbolic edge distance: +// cosh(d_ij) = cosh(b_i + a_ij/2)·cosh(b_j + a_ij/2) − sinh(b_i)·sinh(b_j) // -// cosh(d_ij) = cosh(b_i + a_ij/2) · cosh(b_j + a_ij/2) − sinh(b_i) · sinh(b_j) +// Trilateration via Möbius map + hyperbolic law of cosines (exact, Phase 6). // -// (Springborn 2020, equation for the distance in the horoball picture.) -// Here b_i = x[v_idx[v_i]], a_ij = x[e_idx[e_ij]]. -// -// Poincaré disk: a point at hyperbolic distance d from the origin lies at -// Euclidean distance tanh(d/2) from the disk centre. -// -// Layout algorithm: BFS with hyperbolic trilateration. +// Optional CutGraph and HolonomyData for closed meshes. // ───────────────────────────────────────────────────────────────────────────── inline Layout2D hyper_ideal_layout( ConformalMesh& mesh, const std::vector& x, - const HyperIdealMaps& maps) + const HyperIdealMaps& maps, + const CutGraph* cut = nullptr, + HolonomyData* holonomy = nullptr, + bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); Layout2D result; @@ -337,71 +549,46 @@ inline Layout2D hyper_ideal_layout( // Hyperbolic distance between two adjacent vertices auto hyp_dist = [&](Halfedge_index h) -> double { - Vertex_index vi = mesh.source(h); - Vertex_index vj = mesh.target(h); + Vertex_index vi = mesh.source(h), vj = mesh.target(h); Edge_index e = mesh.edge(h); double bi = get_b(vi), bj = get_b(vj), a = get_a(e); double half_a = a * 0.5; - double cosh_d = std::cosh(bi + half_a) * std::cosh(bj + half_a) - - std::sinh(bi) * std::sinh(bj); - cosh_d = std::max(1.0, cosh_d); - return std::acosh(cosh_d); + double ch = std::cosh(bi + half_a) * std::cosh(bj + half_a) + - std::sinh(bi) * std::sinh(bj); + return std::acosh(std::max(1.0, ch)); }; - // Poincaré disk radius for hyperbolic distance d - auto poincare_r = [](double d) { return std::tanh(d * 0.5); }; - - // Hyperbolic trilateration in Poincaré disk: - // Given p_a, p_b and hyperbolic distances d_a, d_b to the new point, - // return the new point on the LEFT side of the geodesic from p_a to p_b. - // Approximation: for short arcs the Poincaré metric is nearly Euclidean; - // we use Euclidean trilaterate on the Poincaré coordinates. - auto trilaterate_hyp = [&]( - const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, - double da, double db) -> Eigen::Vector2d - { - // Convert arc-lengths to Poincaré chord lengths and use Euclidean formula. - // More precisely: in Poincaré disk, geodesic distance da from pa corresponds - // to a chord of Euclidean length 2·sinh(da/2) / (cosh(da/2)+1) = tanh(da/2)*2/(1+1)... - // For robustness we use the isometric approximation: small displacements are - // Euclidean in conformal coordinates. For a production implementation replace - // with exact Möbius-transform-based hyperbolic trilateration. - double ra = poincare_r(da); - double rb = poincare_r(db); - return detail::trilaterate_2d(pa, pb, ra, rb); - }; - - // Place root face std::vector vertex_placed(nv, false); std::vector face_placed(mesh.number_of_faces(), false); + // ── Place root face: vA at origin, vB on positive real axis ────────────── Face_index f0 = *mesh.faces().begin(); Halfedge_index h0 = mesh.halfedge(f0); Halfedge_index h1 = mesh.next(h0); Halfedge_index h2 = mesh.next(h1); - - Vertex_index vA = mesh.source(h0); - Vertex_index vB = mesh.source(h1); - Vertex_index vC = mesh.source(h2); - + Vertex_index vA = mesh.source(h0); + Vertex_index vB = mesh.source(h1); + Vertex_index vC = mesh.source(h2); double dAB = hyp_dist(h0); double dCA = hyp_dist(h2); double dBC = hyp_dist(h1); - // Place vA at origin, vB on positive x-axis (Poincaré radius = tanh(dAB/2)) result.uv[vA.idx()] = Eigen::Vector2d(0.0, 0.0); - result.uv[vB.idx()] = Eigen::Vector2d(poincare_r(dAB), 0.0); - result.uv[vC.idx()] = trilaterate_hyp( - result.uv[vA.idx()], result.uv[vB.idx()], dCA, dBC); - + result.uv[vB.idx()] = Eigen::Vector2d(std::tanh(dAB * 0.5), 0.0); + result.uv[vC.idx()] = detail::trilaterate_hyp( + result.uv[vA.idx()], result.uv[vB.idx()], dAB, dCA, dBC); vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true; face_placed[f0.idx()] = true; + // ── BFS ─────────────────────────────────────────────────────────────────── std::queue q; auto enqueue = [&](Face_index f) { - for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - auto h_opp = mesh.opposite(h); - if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) + for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index h_opp = mesh.opposite(h); + if (mesh.is_border(h_opp)) continue; + if (cut && cut->is_cut(mesh.edge(h))) continue; + Face_index f_adj = mesh.face(h_opp); + if (!face_placed[f_adj.idx()]) q.push(h_opp); } }; @@ -416,13 +603,15 @@ inline Layout2D hyper_ideal_layout( Vertex_index v_tgt = mesh.target(h); Vertex_index v_new = mesh.target(mesh.next(h)); - Eigen::Vector2d p = trilaterate_hyp( - result.uv[v_src.idx()], result.uv[v_tgt.idx()], - hyp_dist(mesh.prev(h)), - hyp_dist(mesh.next(h))); + double D = hyp_dist(h); // distance v_src → v_tgt + double da = hyp_dist(mesh.prev(h)); // distance v_new → v_src + double db = hyp_dist(mesh.next(h)); // distance v_tgt → v_new + + Eigen::Vector2d p = detail::trilaterate_hyp( + result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db); if (!vertex_placed[v_new.idx()]) { - result.uv[v_new.idx()] = p; + result.uv[v_new.idx()] = p; vertex_placed[v_new.idx()] = true; } else { result.has_seam = true; @@ -432,10 +621,43 @@ inline Layout2D hyper_ideal_layout( } result.success = true; + + // ── Holonomy ────────────────────────────────────────────────────────────── + if (cut && holonomy) { + holonomy->cut_edge_indices = cut->cut_edge_indices; + holonomy->translations.reserve(cut->cut_edge_indices.size()); + + for (std::size_t ce_idx : cut->cut_edge_indices) { + Edge_index e = *std::next(mesh.edges().begin(), + static_cast(ce_idx)); + Halfedge_index h_opp = mesh.opposite(mesh.halfedge(e)); + Halfedge_index h_cross = mesh.is_border(h_opp) + ? mesh.halfedge(e) : h_opp; + if (mesh.is_border(h_cross)) { + holonomy->translations.push_back(Eigen::Vector2d::Zero()); + continue; + } + Vertex_index v_src = mesh.source(h_cross); + Vertex_index v_tgt = mesh.target(h_cross); + Vertex_index v_new = mesh.target(mesh.next(h_cross)); + if (!vertex_placed[v_new.idx()]) { + holonomy->translations.push_back(Eigen::Vector2d::Zero()); + continue; + } + double D = hyp_dist(h_cross); + double da = hyp_dist(mesh.prev(h_cross)); + double db = hyp_dist(mesh.next(h_cross)); + Eigen::Vector2d p_tri = detail::trilaterate_hyp( + result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db); + holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]); + } + } + + if (normalise) normalise_hyperbolic(result); return result; } -// ── Convenience: save layout as OFF (z = 0 for 2D, full xyz for 3D) ────────── +// ── Convenience: save layout as OFF ────────────────────────────────────────── inline void save_layout_off( const std::string& path, diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 12cbd56..2f15d73 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -36,6 +36,9 @@ add_executable(conformallab_cgal_tests # ── Phase 5: Layout / embedding + serialization ──────────────────────── test_layout.cpp + + # ── Phase 6: Gauss–Bonnet, cut graph, exact trilateration, normalisation + test_phase6.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_phase6.cpp b/code/tests/cgal/test_phase6.cpp new file mode 100644 index 0000000..c616bed --- /dev/null +++ b/code/tests/cgal/test_phase6.cpp @@ -0,0 +1,406 @@ +// test_phase6.cpp +// +// Phase 6 — Tests for: +// - gauss_bonnet.hpp : euler_characteristic, genus, GB sum/rhs, check, enforce +// - cut_graph.hpp : compute_cut_graph (tree-cotree algorithm) +// - layout.hpp Phase6 : exact hyperbolic trilateration math +// normalise_euclidean (centroid + PCA) +// normalise_hyperbolic (Möbius centering) + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include "newton_solver.hpp" +#include "gauss_bonnet.hpp" +#include "cut_graph.hpp" +#include "layout.hpp" +#include +#include +#include +#include + +using namespace conformallab; + +// mesh_builder.hpp provides make_triangle(), make_tetrahedron(), make_quad_strip(). + +// ════════════════════════════════════════════════════════════════════════════ +// GaussBonnet — topology helpers +// ════════════════════════════════════════════════════════════════════════════ + +TEST(GaussBonnet, EulerCharacteristic_Triangle) +{ + // V=3 E=3 F=1 → χ = 1 + auto m = make_triangle(); + EXPECT_EQ(euler_characteristic(m), 1); +} + +TEST(GaussBonnet, EulerCharacteristic_QuadStrip) +{ + // V=4 E=5 F=2 → χ = 1 + auto m = make_quad_strip(); + EXPECT_EQ(euler_characteristic(m), 1); +} + +TEST(GaussBonnet, EulerCharacteristic_Tetrahedron) +{ + // V=4 E=6 F=4 → χ = 2 + auto m = make_tetrahedron(); + EXPECT_EQ(euler_characteristic(m), 2); +} + +TEST(GaussBonnet, Genus_Tetrahedron_IsZero) +{ + auto m = make_tetrahedron(); + EXPECT_EQ(genus(m), 0); +} + +TEST(GaussBonnet, RHS_Triangle) +{ + // χ=1 → 2π·χ = 2π + auto m = make_triangle(); + EXPECT_NEAR(gauss_bonnet_rhs(m), 2.0 * M_PI, 1e-12); +} + +TEST(GaussBonnet, RHS_Tetrahedron) +{ + // χ=2 → 2π·χ = 4π + auto m = make_tetrahedron(); + EXPECT_NEAR(gauss_bonnet_rhs(m), 4.0 * M_PI, 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// GaussBonnet — sum / deficit / check / enforce +// ════════════════════════════════════════════════════════════════════════════ + +TEST(GaussBonnet, DeficitIsNonZeroWithDefaultTheta) +{ + // Default theta_v = 2π at all V=4 vertices of quad-strip → + // Σ(2π−2π) = 0, rhs = 2π·1 = 2π → deficit = −2π ≠ 0. + auto m = make_quad_strip(); + auto maps = setup_euclidean_maps(m); + double def = gauss_bonnet_deficit(m, maps); + EXPECT_GT(std::abs(def), 1e-6); +} + +TEST(GaussBonnet, EnforceAdjustsTheta_Deficit_BecomesZero) +{ + auto m = make_quad_strip(); + auto maps = setup_euclidean_maps(m); + // Set clearly wrong theta + for (auto v : m.vertices()) maps.theta_v[v] = M_PI; + EXPECT_GT(std::abs(gauss_bonnet_deficit(m, maps)), 1e-6); + + enforce_gauss_bonnet(m, maps); + EXPECT_NEAR(gauss_bonnet_deficit(m, maps), 0.0, 1e-10); +} + +TEST(GaussBonnet, CheckThrowsWhenViolated) +{ + auto m = make_triangle(); + auto maps = setup_euclidean_maps(m); + // theta_v = 2π everywhere → sum = 0, rhs = 2π → |deficit| = 2π + EXPECT_THROW(check_gauss_bonnet(m, maps), std::runtime_error); +} + +TEST(GaussBonnet, CheckPassesAfterEnforce) +{ + auto m = make_triangle(); + auto maps = setup_euclidean_maps(m); + enforce_gauss_bonnet(m, maps); + EXPECT_NO_THROW(check_gauss_bonnet(m, maps)); +} + +TEST(GaussBonnet, SumAfterEnforce_EqualsRHS) +{ + auto m = make_tetrahedron(); + auto maps = setup_euclidean_maps(m); + // theta_v starts at 2π everywhere; sum = 0, rhs = 4π + enforce_gauss_bonnet(m, maps); + double lhs = gauss_bonnet_sum(m, maps); + double rhs = gauss_bonnet_rhs(m); + EXPECT_NEAR(lhs, rhs, 1e-10); +} + +TEST(GaussBonnet, ManuallySetAnalyticalTheta_PassesCheck) +{ + // Single triangle (χ=1): need Σ(2π−Θ_v) = 2π. + // Set each of the 3 boundary vertices to Θ_v = 4π/3. + // Then Σ(2π − 4π/3) = 3·(2π/3) = 2π. ✓ + auto m = make_triangle(); + auto maps = setup_euclidean_maps(m); + for (auto v : m.vertices()) maps.theta_v[v] = 4.0 * M_PI / 3.0; + EXPECT_NO_THROW(check_gauss_bonnet(m, maps)); +} + +// ════════════════════════════════════════════════════════════════════════════ +// CutGraph — tree-cotree algorithm +// ════════════════════════════════════════════════════════════════════════════ + +TEST(CutGraph, Triangle_ZeroCutEdges) +{ + // Open disk → genus 0 → 0 cut edges. + auto m = make_triangle(); + auto cg = compute_cut_graph(m); + EXPECT_EQ(cg.cut_edge_indices.size(), 0u); +} + +TEST(CutGraph, QuadStrip_ZeroCutEdges) +{ + // Open disk → genus 0 → 0 cut edges. + auto m = make_quad_strip(); + auto cg = compute_cut_graph(m); + EXPECT_EQ(cg.cut_edge_indices.size(), 0u); +} + +TEST(CutGraph, Tetrahedron_ZeroCutEdges_ClosedSphere) +{ + // Closed genus-0 surface → 2g = 0 cut edges. + auto m = make_tetrahedron(); + auto cg = compute_cut_graph(m); + EXPECT_EQ(cg.cut_edge_indices.size(), 0u); + EXPECT_EQ(cg.genus, 0); +} + +TEST(CutGraph, FlagsVsIndicesConsistent) +{ + // Every index in cut_edge_indices has flag=true; + // count of true flags == number of indices. + auto m = make_tetrahedron(); + auto cg = compute_cut_graph(m); + for (std::size_t idx : cg.cut_edge_indices) + EXPECT_TRUE(cg.cut_edge_flags[idx]); + + std::size_t count = 0; + for (bool b : cg.cut_edge_flags) count += b ? 1u : 0u; + EXPECT_EQ(count, cg.cut_edge_indices.size()); +} + +TEST(CutGraph, IsCutMethodMatchesFlags) +{ + auto m = make_tetrahedron(); + auto cg = compute_cut_graph(m); + for (auto e : m.edges()) { + bool flag = cg.cut_edge_flags[static_cast(e.idx())]; + EXPECT_EQ(cg.is_cut(e), flag); + } +} + +TEST(CutGraph, FlagsVectorHasCorrectSize) +{ + auto m = make_quad_strip(); + auto cg = compute_cut_graph(m); + EXPECT_EQ(cg.cut_edge_flags.size(), m.number_of_edges()); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Exact hyperbolic trilateration — Möbius + law of cosines +// ════════════════════════════════════════════════════════════════════════════ + +namespace { + +/// Poincaré-disk hyperbolic distance. +double hyp_dist_disk(Eigen::Vector2d a, Eigen::Vector2d b) +{ + double num = (a.x()-b.x())*(a.x()-b.x()) + (a.y()-b.y())*(a.y()-b.y()); + double da = 1.0 - a.x()*a.x() - a.y()*a.y(); + double db_ = 1.0 - b.x()*b.x() - b.y()*b.y(); + double arg = 1.0 + 2.0*num/(da*db_); + return std::acosh(std::max(1.0, arg)); +} + +/// Reproduce the exact trilateration formula from layout.hpp detail namespace. +Eigen::Vector2d trilaterate_hyp(Eigen::Vector2d pa, Eigen::Vector2d pb, + double D, double da, double db) +{ + if (D < 1e-12 || da < 1e-12) return pa; + if (std::sinh(da) * std::sinh(D) < 1e-14) return pa; + using C = std::complex; + auto mobius_fwd = [](C z, C center) -> C { + return (z - center) / (C(1.0) - std::conj(center) * z); + }; + auto mobius_inv = [](C w, C center) -> C { + return (w + center) / (C(1.0) + std::conj(center) * w); + }; + C a(pa.x(), pa.y()), b(pb.x(), pb.y()); + C b_mapped = mobius_fwd(b, a); + double theta_b = std::arg(b_mapped); + double cos_alpha = (std::cosh(da) * std::cosh(D) - std::cosh(db)) + / (std::sinh(da) * std::sinh(D)); + cos_alpha = std::max(-1.0, std::min(1.0, cos_alpha)); + double alpha = std::acos(cos_alpha); + C pc_origin = std::tanh(da * 0.5) * std::exp(C(0.0, theta_b + alpha)); + C pc_c = mobius_inv(pc_origin, a); + return Eigen::Vector2d(pc_c.real(), pc_c.imag()); +} + +} // namespace + +TEST(HyperbolicTrilateration, RecoveredDistances_Small) +{ + // pa at origin, pb at tanh(D/2) on x-axis. + double D = 0.8, da = 0.5, db = 0.6; + Eigen::Vector2d pa(0.0, 0.0); + Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0); + + Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db); + + EXPECT_NEAR(hyp_dist_disk(pa, pb), D, 1e-10); + EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-9); + EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-9); +} + +TEST(HyperbolicTrilateration, RecoveredDistances_Large) +{ + double D = 2.0, da = 1.5, db = 1.0; + Eigen::Vector2d pa(0.0, 0.0); + Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0); + + Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db); + + EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-8); + EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-8); +} + +TEST(HyperbolicTrilateration, ResultInsidePoincareDisk) +{ + double D = 1.2, da = 0.9, db = 0.7; + Eigen::Vector2d pa(0.0, 0.0); + Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0); + + Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db); + EXPECT_LT(pc.norm(), 1.0); +} + +TEST(HyperbolicTrilateration, OffOrigin_StillCorrect) +{ + // Place pa somewhere in the disk (not at origin) to test Möbius map. + double D = 0.6, da = 0.4, db = 0.5; + // pa at (0.3, 0.0) in Poincaré coords. + Eigen::Vector2d pa(0.3, 0.0); + // pb at distance D from pa — compute pb in Poincaré disk: + // Möbius inverse: tanh(D/2) maps back from origin frame. + using C = std::complex; + C a_c(pa.x(), pa.y()); + C pb_c = (std::tanh(D*0.5) + a_c) / (C(1.0) + std::conj(a_c) * std::tanh(D*0.5)); + Eigen::Vector2d pb(pb_c.real(), pb_c.imag()); + + Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db); + EXPECT_LT(pc.norm(), 1.0); + EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-8); + EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-8); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Normalisation — euclidean_layout (centroid) + hyper_ideal_layout (Möbius) +// ════════════════════════════════════════════════════════════════════════════ + +/// Build a solved euclidean layout for the quad-strip at the natural equilibrium. +static Layout2D make_euclidean_layout_normalised(bool normalise) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Pin first vertex; assign free DOF indices to all others. + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++; + const int n = idx; + + std::vector x0(static_cast(n), 0.0); + // Adjust theta so G(x=0) = 0 (natural equilibrium). + auto G0 = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] -= G0[static_cast(iv)]; + } + + auto res = newton_euclidean(mesh, x0, maps, 1e-10, 100); + return euclidean_layout(mesh, res.x, maps, nullptr, nullptr, normalise); +} + +TEST(Normalisation, Euclidean_CentroidAtOrigin) +{ + auto L = make_euclidean_layout_normalised(true); + ASSERT_GT(L.uv.size(), 0u); + + double cx = 0.0, cy = 0.0; + for (auto& p : L.uv) { cx += p.x(); cy += p.y(); } + int n = static_cast(L.uv.size()); + EXPECT_NEAR(cx / n, 0.0, 1e-9); + EXPECT_NEAR(cy / n, 0.0, 1e-9); +} + +TEST(Normalisation, Euclidean_NormalisePreservesEdgeLengthRatios) +{ + auto L_no = make_euclidean_layout_normalised(false); + auto L_yes = make_euclidean_layout_normalised(true); + + // Ratio of the lengths of the first two edges must be preserved. + ASSERT_GE(L_no.uv.size(), 4u); + // edges 0→1 and 0→2 (the two edges from vertex 0 in the quad-strip) + double len0_no = (L_no.uv[1] - L_no.uv[0]).norm(); + double len1_no = (L_no.uv[2] - L_no.uv[0]).norm(); + double len0_yes = (L_yes.uv[1] - L_yes.uv[0]).norm(); + double len1_yes = (L_yes.uv[2] - L_yes.uv[0]).norm(); + + if (len1_no > 1e-12 && len1_yes > 1e-12) { + double ratio_no = len0_no / len1_no; + double ratio_yes = len0_yes / len1_yes; + EXPECT_NEAR(ratio_no, ratio_yes, 1e-9); + } +} + +/// Build a solved hyper-ideal layout for a triangle at natural equilibrium. +static Layout2D make_hyper_ideal_layout_normalised(bool normalise) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + std::vector xbase(static_cast(n)); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; if (iv >= 0) xbase[iv] = 1.0; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; if (ie >= 0) xbase[ie] = 0.5; + } + auto G0 = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; if (iv >= 0) maps.theta_v[v] += G0[iv]; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; if (ie >= 0) maps.theta_e[e] += G0[ie]; + } + auto res = newton_hyper_ideal(mesh, xbase, maps, 1e-9, 100); + return hyper_ideal_layout(mesh, res.x, maps, nullptr, nullptr, normalise); +} + +TEST(Normalisation, HyperIdeal_PositionsInsidePoincareDisk) +{ + auto L = make_hyper_ideal_layout_normalised(true); + ASSERT_GT(L.uv.size(), 0u); + for (auto& p : L.uv) { + EXPECT_LT(p.norm(), 1.0 + 1e-9) + << "Vertex outside Poincaré disk: r=" << p.norm(); + } +} + +TEST(Normalisation, HyperIdeal_CenteringReducesAverageRadius) +{ + // After Möbius centering the average Euclidean radius inside the disk + // must be ≤ the unconstrained average. + auto L_no = make_hyper_ideal_layout_normalised(false); + auto L_yes = make_hyper_ideal_layout_normalised(true); + + int n = static_cast(L_no.uv.size()); + double r_no = 0.0, r_yes = 0.0; + for (int i = 0; i < n; ++i) { + r_no += L_no.uv[i].norm(); + r_yes += L_yes.uv[i].norm(); + } + EXPECT_LE(r_yes / n, r_no / n + 1e-9); +}