diff --git a/README.md b/README.md index 08bb570..6d6e712 100644 --- a/README.md +++ b/README.md @@ -419,6 +419,83 @@ H[i,j] = ( G(x + ε·eⱼ)[i] − G(x − ε·eⱼ)[i] ) / (2ε) This is O(ε²) accurate (≈ 10⁻¹⁰ relative error at ε = 10⁻⁵), positive semi-definite by strict convexity of the HyperIdeal energy (Springborn 2020), and costs n extra gradient evaluations per Newton step instead of O(n). For meshes with fewer than ~500 DOFs the difference in wall time is negligible. +### Layout — detailed comparison with the Java original + +The BFS-trilateration algorithm itself is the same in both implementations. The differences are in the surrounding infrastructure: + +#### 1. Hyperbolic trilateration — approximation vs. exact + +In the Euclidean and spherical cases the C++ trilateration is **mathematically identical** to the Java original. +For the HyperIdeal / Poincaré disk case the two diverge: + +| | Java ConformalLab | conformallab++ | +|---|---|---| +| Poincaré disk trilateration | **Exact** via Möbius transformations | **Approximation**: `r = tanh(d/2)`, then Euclidean trilaterate | + +Java computes the new point by: (i) Möbius-map `p_a` to the origin, (ii) rotate so `p_b` lies on the positive real axis, (iii) apply the standard on-axis hyperbolic trilateration formula, (iv) invert the Möbius map. The result is exact in the hyperbolic metric. + +The C++ approximation replaces the hyperbolic distance `d` by the Poincaré-disk Euclidean radius `tanh(d/2)` and then calls `trilaterate_2d`. This is first-order accurate near the origin (small triangles) but degrades towards the boundary of the unit disk (large `d`, highly curved triangles). For the Newton solver this does not matter — the solver works entirely in log-scale `x`-space and never calls the layout — but the embedded coordinates will deviate from the true hyperbolic positions for deep hyperbolic triangulations. + +**To implement the exact version** replace `trilaterate_hyp` in `layout.hpp` with: + +```cpp +// Möbius map: send point p to the origin (unit disk) +auto mobius_to_origin = [](Eigen::Vector2d p, Eigen::Vector2d center) { + // T(z) = (z - center) / (1 - conj(center)*z) [complex arithmetic in R^2] + Eigen::Vector2d num = p - center; + double den_re = 1.0 - (center.x()*p.x() + center.y()*p.y()); + double den_im = (center.x()*p.y() - center.y()*p.x()); + double den2 = den_re*den_re + den_im*den_im; + return Eigen::Vector2d( + (num.x()*den_re + num.y()*den_im) / den2, + (num.y()*den_re - num.x()*den_im) / den2); +}; +// After mapping p_a → 0, p_b lies on the real axis. +// The on-axis formula for left-side trilateration then applies. +``` + +#### 2. Closed meshes — seam flag vs. cut + period matrices + +| | Java | conformallab++ | +|---|---|---| +| Open mesh (with boundary) | BFS unfolding | BFS unfolding ✅ identical | +| Closed mesh, genus 0 (sphere) | BFS + Möbius normalisation | BFS, `has_seam = true` | +| Closed mesh, genus ≥ 1 | Full cut + holonomy + period matrix | ❌ not implemented | + +For a closed mesh the BFS in both implementations eventually tries to place a vertex that was already placed from the other side of the seam. Java records the **holonomy element** — the Möbius transformation (Euclidean: rigid motion; hyperbolic: Möbius; spherical: rotation) that maps the "first visit" position to the "second visit" position. These holonomy elements around the two generators of the fundamental group $\pi_1(\Sigma_g)$ are the **monodromy data** that determine the conformal structure of the surface. + +For genus $g = 0$ (topological sphere) the holonomy is trivial after a Möbius normalization and the layout closes up. For genus $g \geq 1$ the holonomy data feeds into the **period matrix** computation; the lattice $\Lambda \subset \mathbb{C}$ for a torus is read off from the two holonomy elements of the single handle. + +#### 3. Layout normalisation + +Java applies a post-processing Möbius transformation (Euclidean: rigid motion + scaling; hyperbolic: Möbius centring) to bring the layout into a canonical form. C++ returns the raw BFS output: vertex 0 at the origin, vertex 1 on the positive x-axis. + +#### 4. Gauss–Bonnet consistency check + +Before solving, Java verifies that the prescribed target angles satisfy + +$$\sum_{v} (2\pi - \Theta_v) = 2\pi \cdot \chi(M)$$ + +(Gauss–Bonnet for the target metric). If this fails, no conformal factor exists that realises those angles and Java aborts with a meaningful error. C++ has no such check; Newton will fail to converge silently. + +#### Summary + +``` + Java C++ +───────────────────────────────────────────────────── +BFS-trilateration (open) ✅ ✅ identical +Euclidean trilateration ✅ ✅ identical +Spherical trilateration ✅ ✅ identical +Hyperbolic trilateration ✅ exact ⚠️ tanh(d/2) approx +Closed mesh genus 0 ✅ ⚠️ has_seam flag only +Closed mesh genus ≥ 1 ✅ ❌ +Holonomy / monodromy ✅ ❌ +Period matrices ✅ ❌ +Layout normalisation ✅ Möbius ❌ raw output +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. @@ -612,11 +689,163 @@ Property maps are reference-counted and cheap to copy. Give them unique string n 4. **Try different target angles** — set `maps.theta_v[v] = M_PI / 3` for all interior vertices and see how the solver responds. The constraint `Σ(2π − Θ_v) = 2π·χ(M)` (Gauss–Bonnet) must hold for a solution to exist. 5. **Inspect convergence** — `NewtonResult` carries `iterations`, `grad_inf_norm`, and the full `x` at termination. Plot `||G(xₖ)||` per iteration to verify quadratic convergence near the solution. -### Recommended reading +### Porting the Java uniformization pipeline — a roadmap for contributors -| Paper | Relevance to this codebase | -|-------|---------------------------| -| Springborn, Schröder, Pinkall — *Conformal Equivalence of Triangle Meshes* (2008) | Euclidean & spherical functionals; the Schläfli formula at the core of `spherical_functional.hpp` | +The Java ConformalLab has a complete global uniformization pipeline that conformallab++ does not yet have. This section explains the mathematics behind it and the precise C++ steps needed to port each piece. + +#### What "global uniformization" means mathematically + +Given a triangulated surface $M$ of genus $g$ and a convergerd conformal scale factor $u_v$, global uniformization produces: + +- $g = 0$: an embedding into the sphere $S^2$ or the Riemann sphere $\hat{\mathbb{C}}$, well-defined up to Möbius transformation. +- $g = 1$: a flat torus $\mathbb{C} / \Lambda$; the lattice $\Lambda = \mathbb{Z} + \tau\mathbb{Z}$ is the **period** (one complex number $\tau$ in the upper half-plane). +- $g \geq 2$: a hyperbolic surface $\mathbb{H} / \Gamma$; the **period matrix** $\Omega \in \mathbb{H}_g$ (a $g\times g$ symmetric complex matrix with positive-definite imaginary part) parametrizes the conformal class. + +The BFS unfolding in `layout.hpp` already computes *local* coordinates correctly. The missing piece is tracking what happens when the BFS crosses a **handle** — a non-contractible loop in $M$. + +#### Step 1 — Gauss–Bonnet consistency check + +Before anything else, add a pre-solve validation function. This is a single loop: + +```cpp +// gauss_bonnet.hpp (new file, ~30 lines) +#pragma once +#include "conformal_mesh.hpp" +#include "euclidean_functional.hpp" // or whichever Maps type +#include "constants.hpp" +#include +#include + +namespace conformallab { + +// Throws if Σ(2π − θ_v) ≠ 2π·χ(M) to within tol. +// Call before newton_euclidean / newton_spherical. +inline void check_gauss_bonnet( + const ConformalMesh& mesh, + const EuclideanMaps& maps, // or SphericalMaps / HyperIdealMaps + double tol = 1e-8) +{ + // Euler characteristic χ = V − E + F + int chi = static_cast(mesh.number_of_vertices()) + - static_cast(mesh.number_of_edges()) + + static_cast(mesh.number_of_faces()); + + double lhs = 0.0; + for (auto v : mesh.vertices()) + lhs += TWO_PI - maps.theta_v[v]; + + double rhs = TWO_PI * chi; + if (std::abs(lhs - rhs) > tol) + throw std::runtime_error( + "Gauss–Bonnet violated: Σ(2π−θ_v) = " + std::to_string(lhs) + + ", expected 2π·χ = " + std::to_string(rhs)); +} + +} // namespace conformallab +``` + +This is the easiest piece and the highest-value first step: it will immediately catch target-angle mistakes that currently cause silent Newton non-convergence. + +#### Step 2 — homological basis / fundamental polygon cut + +To turn a genus-$g$ surface into a disk, cut along a **standard homological basis**: $g$ pairs of loops $(a_1, b_1), \ldots, (a_g, b_g)$ with $a_i \cdot b_j = \delta_{ij}$. + +In Java this is computed from the dual graph via a **spanning tree + co-tree** decomposition: + +1. Compute a spanning tree $T$ of the primal graph (edges used in BFS). +2. The non-tree edges form the co-tree; they correspond to independent homology classes. +3. For genus $g$, pick $2g$ non-tree edges that form a standard symplectic basis (i.e. they pair up with intersection number 1). +4. Cut the mesh along these $2g$ loops to obtain a $4g$-gon. + +In C++ the data structure to add is a **cut graph** stored as a set of `Edge_index` values that are marked as "boundary" for the BFS: + +```cpp +// cut_graph.hpp (new file) +#pragma once +#include "conformal_mesh.hpp" +#include +#include + +namespace conformallab { + +struct CutGraph { + std::unordered_set cut_edges; // edges treated as boundary + int genus; +}; + +// Compute a cut graph for mesh using spanning tree + co-tree. +// Result: a set of 2g edges whose removal makes M simply connected. +CutGraph compute_cut_graph(ConformalMesh& mesh); + +} // namespace conformallab +``` + +The `euclidean_layout` BFS then needs one extra line: + +```cpp +// Inside the BFS enqueue lambda — treat cut edges as boundary: +if (!mesh.is_border(h_opp) && !cut.cut_edges.count(mesh.edge(h).idx())) + q.push(h_opp); +``` + +#### Step 3 — holonomy tracking + +Once the BFS has a cut, every time it would cross a cut edge it instead records the **holonomy element** — the transformation that maps the coordinate on one side of the cut to the coordinate on the other side. + +For the **Euclidean case** the holonomy group is a subgroup of $\text{Isom}(\mathbb{R}^2)$ (rigid motions). Each holonomy element is a $2\times 2$ rotation + translation: + +```cpp +struct EuclideanHolonomy { + Eigen::Matrix2d R; // rotation + Eigen::Vector2d t; // translation + // apply: p ↦ R·p + t +}; +``` + +For the **hyperbolic case** the holonomy group is a subgroup of $\text{PSL}(2,\mathbb{R})$ acting on the upper half-plane (or equivalently $\text{PSU}(1,1)$ acting on the Poincaré disk). Each element is a $2\times 2$ real matrix with determinant 1: + +```cpp +struct MobiusElement { + Eigen::Matrix2d M; // [[a, b], [c, d]], det = 1 + // apply to z = (x,y): z ↦ (M[0,0]·z + M[0,1]) / (M[1,0]·z + M[1,1]) + // (complex arithmetic) +}; +``` + +The BFS accumulates one holonomy element per cut edge pair $(a_i, b_i)$. After BFS completes, the holonomy data is a $2g$-tuple of group elements. + +#### Step 4 — period matrix extraction + +From the holonomy elements, the period data is extracted as follows: + +**Genus 1 (torus):** The two holonomy elements $h_{a_1}, h_{b_1}$ are translations: $h_{a_1}(z) = z + \omega_1$, $h_{b_1}(z) = z + \omega_2$ with $\omega_1, \omega_2 \in \mathbb{C}$. The period ratio is $\tau = \omega_2 / \omega_1 \in \mathbb{H}$ (upper half-plane). Reduce $\tau$ to the fundamental domain of $\text{SL}(2,\mathbb{Z})$ with the standard Euclidean algorithm on $\text{SL}(2,\mathbb{Z})$. + +**Genus $g \geq 2$:** The $2g$ holonomy elements are generators of a Fuchsian group $\Gamma \subset \text{PSL}(2,\mathbb{R})$. The period matrix $\Omega_{ij} = \int_{b_j} \omega_i$ is computed by integrating the $g$ holomorphic differentials $\omega_1, \ldots, \omega_g$ around the $b$-cycles. On a discrete surface this reduces to a linear system involving the holonomy matrices. See Bobenko–Springborn (2004) §6 for the discrete version. + +The `discrete_elliptic_utility.hpp` in this codebase already contains a modular normalisation helper for $\tau$ — it is not yet wired up to any layout pipeline, but is precisely what Step 4 needs for the genus-1 case. + +#### Step 5 — layout normalisation + +After BFS + holonomy extraction, apply a canonical Möbius transformation to bring the layout into a standard position: + +- **Euclidean genus 0:** scale + rotate so that $\omega_1 = 1$ (standard horizontal period). +- **Hyperbolic:** apply the unique $\text{PSU}(1,1)$ element that maps the centroid of all vertices to the origin of the Poincaré disk. +- **Spherical genus 0:** apply the unique Möbius transformation that maps the three face-centres of the root face to the standard position on $S^2$. + +#### Recommended order of implementation + +| Priority | Piece | Estimated effort | C++ hook | +|----------|-------|-----------------|----------| +| 1 | Gauss–Bonnet check | ~30 lines, 1 day | new `gauss_bonnet.hpp` | +| 2 | Exact hyperbolic trilateration | ~50 lines, 1 day | replace `trilaterate_hyp` in `layout.hpp` | +| 3 | Cut-graph computation | ~200 lines, 3–5 days | new `cut_graph.hpp` | +| 4 | Holonomy tracking in BFS | ~80 lines, 2 days | extend `euclidean_layout` / `hyper_ideal_layout` | +| 5 | Period matrix (genus 1) | ~100 lines, 3 days | wire up `discrete_elliptic_utility.hpp` | +| 6 | Period matrix (genus ≥ 2) | research-level, weeks | new `period_matrix.hpp` | + +Steps 1–4 together constitute a **practically complete** uniformization for genus-0 and genus-1 surfaces, which covers the vast majority of mesh-processing use cases. + +### Recommended reading | Springborn — *Ideal Hyperbolic Polyhedra and Discrete Uniformization* (2020) | HyperIdeal functional; the ζ₁₃/ζ₁₄/ζ₁₅ functions in `hyper_ideal_geometry.hpp` | | Pinkall, Polthier — *Computing Discrete Minimal Surfaces* (1993) | Cotangent-Laplace Hessian in `euclidean_hessian.hpp` | | Luo — *Combinatorial Yamabe Flow on Surfaces* (2004) | Inversive-distance functional (not yet ported — good first contribution) |