ConformalLabpp/README.md

conformallab++

conformallab++ is a modern C++ reimplementation of the ConformalLab software by Stefan Sechelmann for experiments in discrete conformal geometry and related mesh transformations.

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.

---

Features

Area	Status
Clausen / Lobachevsky / ImLi₂ functions	✅ Phase 1
Hyper-ideal geometry (ζ, lᵢⱼ, αᵢⱼ, σᵢ, σᵢⱼ)	✅ Phase 2
CGAL Surface_mesh infrastructure + mesh builders	✅ Phase 3a
Hyper-ideal functional (energy + gradient)	✅ Phase 3b
Spherical functional (energy + gradient + gauge-fix)	✅ Phase 3c/3e
Euclidean functional (energy + gradient)	✅ Phase 3d
Euclidean Hessian (cotangent-Laplace, Pinkall–Polthier)	✅ Phase 3f
Spherical Hessian (∂α/∂u from law of cosines)	✅ Phase 3f
Hyper-ideal Hessian (numerical FD, symmetrised)	✅ Phase 4a
Newton solver — all three geometries	✅ Phase 4a
SparseQR fallback for rank-deficient H (gauge modes)	✅ Phase 4a
Mesh I/O (CGAL::IO — OFF / OBJ / PLY)	✅ Phase 4b
End-to-end pipeline tests	✅ Phase 4c
Example programs (headless + interactive viewer)	✅ Phase 4d
BFS Layout (ℝ², S², Poincaré disk)	✅ Phase 5
CLI app (conformallab_core)	✅ Phase 5
JSON + XML serialisation	✅ Phase 5

---

Quick start — CLI app

cmake -S code -B build -DWITH_CGAL=ON
cmake --build build -j4

# Run the conformal map CLI on any OFF/OBJ/PLY mesh
./bin/conformallab_core -i input.off -g euclidean -o layout.off -j result.json -x result.xml

# Available geometries: euclidean | spherical | hyper_ideal
./bin/conformallab_core -i input.off -g spherical -o sphere.off

# Show input mesh in interactive viewer (built-in)
./bin/conformallab_core -i input.off -s

Example programs

# Layout + JSON/XML round-trip demo
./build/examples/example_layout  [input.off]  [layout.off]  [result.json]  [result.xml]

# Headless pipelines
./build/examples/example_euclidean   [input.off]  [output.off]
./build/examples/example_hyper_ideal [input.off]  [output.off]

# Interactive viewer (requires -DWITH_CGAL=ON, viewer is built automatically)
./build/examples/example_viewer  [input.off]

If no input file is given each example uses a built-in quad-strip mesh.

---

Library usage — minimal Euclidean pipeline

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "newton_solver.hpp"

using namespace conformallab;

int main() {
    // 1. Load mesh
    ConformalMesh mesh = load_mesh("input.off");

    // 2. Set up functional maps
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // 3. Assign DOFs — pin first vertex (gauge fix)
    auto vit = mesh.vertices().begin();
    maps.v_idx[*vit++] = -1;   // pinned: u[v0] = 0
    int idx = 0;
    for (; vit != mesh.vertices().end(); ++vit)
        maps.v_idx[*vit] = idx++;
    const int n = idx;

    // 4. Set target angles (natural equilibrium: x* = 0)
    std::vector<double> x0(n, 0.0);
    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[iv];
    }

    // 5. Solve
    auto result = newton_euclidean(mesh, std::vector<double>(n, -0.1), maps);

    // 6. Save result
    if (result.converged)
        save_mesh("output.off", mesh);

    return result.converged ? 0 : 1;
}

For HyperIdeal geometry:

auto maps = setup_hyper_ideal_maps(mesh);
int n = assign_all_dof_indices(mesh, maps);
// … set theta_v / theta_e targets …
auto result = newton_hyper_ideal(mesh, x0, maps);

Layout (embedding into target space)

After solving, convert scale factors into actual vertex coordinates:

#include "layout.hpp"
#include "serialization.hpp"

// Euclidean: flat 2-D positions in ℝ²
Layout2D layout = euclidean_layout(mesh, result.x, maps);
// layout.uv[v.idx()] = Eigen::Vector2d

// Spherical: unit vectors on S² ⊂ ℝ³
Layout3D slayout = spherical_layout(mesh, result.x, smaps);
// slayout.pos[v.idx()] = Eigen::Vector3d

// HyperIdeal: Poincaré disk coordinates in ℝ²
Layout2D hlayout = hyper_ideal_layout(mesh, result.x, hmaps);

// Save layout as OFF
save_layout_off("layout.off", mesh, layout);

// Serialise to JSON / XML
save_result_json("result.json", result, "euclidean",
    mesh.number_of_vertices(), mesh.number_of_faces(), &layout);
save_result_xml("result.xml", result, "euclidean",
    mesh.number_of_vertices(), mesh.number_of_faces(), &layout);

// Load back
NewtonResult res2; std::string geom; Layout2D uv2;
auto x = load_result_json("result.json", &res2, &geom, &uv2);

Using solve_linear_system directly (with fallback detection):

#include "newton_solver.hpp"

bool used_fallback = false;
auto dx = conformallab::solve_linear_system(H, rhs, &used_fallback);
if (used_fallback)
    std::cout << "SparseQR was used (H is rank-deficient)\n";

---

Build modes

Mode	CMake flags	What gets built	CI
Tests only (default)	(none)	conformallab_tests — Eigen + GTest only	✅ automatic
CGAL tests + examples	-DWITH_CGAL=ON	conformallab_cgal_tests, example_euclidean, example_hyper_ideal, example_layout, conformallab_core (CLI)	local only
Interactive viewer	-DWITH_CGAL=ON (implied)	above + example_viewer + viewer linked into CLI	local only

External dependencies are bundled as tarballs in code/deps/tarballs/ and extracted lazily at CMake configure time (GTest is fetched via FetchContent).

Boost is required only with -DWITH_CGAL=ON (header-only use by CGAL 6.x).

---

Prerequisites

Tool	Minimum
C++ compiler (GCC or Clang)	C++17
CMake	3.20
Boost headers	1.70 (only with -DWITH_CGAL=ON)

---

Getting started

git clone https://codeberg.org/TMoussa/ConformalLabpp
cd ConformalLabpp

Tests only (CI default — no system deps)

cmake -S code -B build
cmake --build build --target conformallab_tests -j$(nproc)
ctest --test-dir build --output-on-failure

CGAL tests + headless examples (needs system Boost)

cmake -S code -B build -DWITH_CGAL=ON
cmake --build build -j$(nproc)
ctest --test-dir build -R "^cgal\." --output-on-failure
./build/examples/example_euclidean
./build/examples/example_hyper_ideal
./build/examples/example_layout
./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).

Interactive viewer

# WITH_CGAL=ON already implies viewer; example_viewer is built automatically
cmake -S code -B build -DWITH_CGAL=ON && cmake --build build -t example_viewer -j$(nproc)
./build/examples/example_viewer data/off/example.off

---

Public headers (code/include/)

Header	Description
clausen.hpp	Clausen Cl₂, Lobachevsky Л, ImLi₂
hyper_ideal_geometry.hpp	ζ functions, lᵢⱼ, αᵢⱼ, σᵢ, σᵢⱼ
hyper_ideal_utility.hpp	Tetrahedron volume (Meyerhoff / Kolpakov–Mednykh)
hyper_ideal_functional.hpp	HyperIdeal energy + gradient on ConformalMesh
hyper_ideal_hessian.hpp	HyperIdeal Hessian (numerical FD, symmetrised)
hyper_ideal_visualization_utility.hpp	Poincaré disk projection, circumcircle helpers
spherical_geometry.hpp	Spherical arc length, half-angle formula
spherical_functional.hpp	Spherical energy + gradient + gauge-fix
spherical_hessian.hpp	Spherical Hessian (∂α/∂u, law of cosines)
euclidean_geometry.hpp	Euclidean corner-angle (t-value / atan2)
euclidean_functional.hpp	Euclidean energy + gradient
euclidean_hessian.hpp	Cotangent-Laplace Hessian (Pinkall–Polthier)
newton_solver.hpp	newton_euclidean / newton_spherical / newton_hyper_ideal + public solve_linear_system
mesh_io.hpp	read_mesh / write_mesh / load_mesh / save_mesh
conformal_mesh.hpp	ConformalMesh = CGAL::Surface_mesh<Point3> + 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
serialization.hpp	save/load_result_json + save/load_result_xml
mesh_utils.hpp	CGAL → Eigen conversion (cgal_to_eigen)
constants.hpp	conformallab::PI, TWO_PI

---

Project structure

code/
├── include/                         # All public headers (header-only library)
│   ├── conformal_mesh.hpp
│   ├── mesh_builder.hpp
│   ├── mesh_io.hpp
│   ├── mesh_utils.hpp
│   ├── newton_solver.hpp            # ← public solve_linear_system + 3 Newton solvers
│   ├── layout.hpp                   # ← BFS layout (euclidean/spherical/hyper_ideal)
│   ├── serialization.hpp            # ← JSON + XML save/load
│   ├── hyper_ideal_{functional,hessian,geometry,utility,visualization_utility}.hpp
│   ├── spherical_{functional,hessian,geometry}.hpp
│   ├── euclidean_{functional,hessian,geometry}.hpp
│   ├── clausen.hpp
│   └── constants.hpp
├── examples/                        # Standalone example programs
│   ├── CMakeLists.txt
│   ├── example_euclidean.cpp        # Headless Euclidean pipeline
│   ├── example_hyper_ideal.cpp      # Headless HyperIdeal pipeline
│   ├── example_layout.cpp           # Solve → layout → OFF/JSON/XML round-trip
│   └── example_viewer.cpp           # Interactive libigl viewer (WITH_VIEWER)
├── src/
│   ├── apps/v0/conformallab_cli.cpp # Full CLI app (Phase 5)
│   └── viewer/simple_viewer.cpp
├── tests/
│   ├── CMakeLists.txt
│   ├── *.cpp                        # conformallab_tests (no CGAL)
│   └── cgal/
│       ├── CMakeLists.txt
│       ├── test_conformal_mesh.cpp          # 14 tests
│       ├── test_hyper_ideal_functional.cpp  #  7 tests (1 skipped)
│       ├── test_spherical_functional.cpp    # 11 tests (1 skipped)
│       ├── test_euclidean_functional.cpp    # 11 tests
│       ├── test_euclidean_hessian.cpp       #  8 tests
│       ├── test_spherical_hessian.cpp       #  8 tests
│       ├── 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)
└── deps/
    ├── eigen-3.4.0/        # always extracted
    ├── CGAL-6.1.1/         # extracted with WITH_CGAL
    ├── libigl-2.6.0/       # extracted with WITH_VIEWER
    ├── glfw-3.4/           # extracted with WITH_VIEWER
    ├── libigl-glad/
    └── single_includes/    # CLI11, json.hpp

---

Test suites

conformallab_tests (CI — always built)

Pure-math tests requiring only Eigen: Clausen / Lobachevsky / ImLi₂, hyper-ideal geometry, tetrahedron volumes.

conformallab_cgal_tests (local — -DWITH_CGAL=ON)

Suite	Tests	What it checks
ConformalMeshTopology	4	Euler characteristic, vertex/edge/face counts
ConformalMeshTraversal	4	Halfedge iteration, valence, opposite
ConformalMeshProperties	5	Property maps (λ, θ, idx, α, geometry type)
ConformalMeshValidity	1	CGAL validity for all factory meshes
HyperIdealFunctional	7	FD gradient checks + Hessian symmetry
SphericalFunctional	11	Angle formula + gradient + gauge-fix
EuclideanFunctional	11	Angle formula + gradient
EuclideanHessian	8	Cotangent-Laplace structure, FD agreement, PSD, null space
SphericalHessian	8	Derivative correctness, NSD at equilibrium
NewtonSolver	11	Convergence (Euclidean ×3, Spherical ×4, HyperIdeal ×4)
SparseQRFallback	3	Full-rank LDLT path · singular matrix triggers QR · closed-mesh gauge-mode
MeshIO	6	OFF/OBJ round-trips, error handling
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)

---

Newton solver & SparseQR fallback

newton_solver.hpp exposes three solvers with a unified interface:

NewtonResult newton_euclidean (mesh, x0, maps [, tol, max_iter])
NewtonResult newton_spherical (mesh, x0, maps [, tol, max_iter])
NewtonResult newton_hyper_ideal(mesh, x0, maps [, tol, max_iter, hess_eps])

Each iteration:

1. Evaluate gradient G
2. Evaluate Hessian H (analytical for Euclidean / Spherical; numerical FD for HyperIdeal)
3. Solve H·Δx = −G — try Eigen::SimplicialLDLT, fall back to Eigen::SparseQR on failure
4. Backtracking line search (up to 20 halvings)

Gradient sign conventions:

Geometry	G	H sign
Euclidean	Θ_v − Σα_v	PSD → LDLT on H
Spherical	Θ_v − Σα_v	NSD → LDLT on −H
HyperIdeal	Σβ_v − Θ_v	PSD → LDLT on H

SparseQR fallback (solve_linear_system):
When SimplicialLDLT fails (singular/rank-deficient H, e.g. gauge modes on closed meshes without a pinned vertex), SparseQR finds the minimum-norm Newton step orthogonal to the null space. Because the gradient always lies in the row space of H, the solver converges correctly without requiring the caller to pin a vertex.

The fallback is a public API:

bool used_fallback = false;
auto dx = conformallab::solve_linear_system(H, rhs, &used_fallback);

---

Mathematical scope — C++ vs. Java original

The three core functionals are fully equivalent to the Java original at the level of energy, gradient, and Hessian formulas. The table below shows where parity holds, where there is a numerical difference, and what is not yet ported.

Mathematical layer	Java ConformalLab	conformallab++
Euclidean functional — energy, gradient	✅	✅
Spherical functional — energy, gradient, gauge-fix	✅	✅
HyperIdeal functional — energy, gradient	✅	✅
Inversive-distance functional (Luo 2004, Bowers–Stephenson)	✅	❌ not ported
Euclidean Hessian — cotangent-Laplace (Pinkall–Polthier 1993)	✅ analytical	✅ analytical
Spherical Hessian — ∂α/∂u from law of cosines	✅ analytical	✅ analytical
HyperIdeal Hessian — through ζ → lᵢⱼ → β/α chain	✅ analytical	⚠️ symmetric FD
Newton solver	✅	✅
SparseQR fallback for gauge-mode null spaces	?	✅
Cone metrics — prescribed Θ_v ≠ 2π	✅ full pipeline	⚠️ data structure only
Boundary conditions — Dirichlet u=f, Neumann, free boundary	✅	⚠️ pin-only
Layout / embedding — DOF vector → vertex coordinates in ℝ² / H² / S²	✅	✅ BFS unfolding (all three geometries)
Gauss–Bonnet consistency check on target angles	✅	❌
Global uniformization for genus g ≥ 1	✅	❌
Period matrices — Teichmüller parameters for g ≥ 2	✅	❌
Holonomy / monodromy	✅	❌
HyperIdeal generator — constructing geometrically valid meshes	✅	❌ only test meshes
Clausen / Lobachevsky / ImLi₂ special functions	✅	✅
Discrete elliptic utility (modular normalisation of τ)	✅	✅ (not yet wired up)
Poincaré disk / Lorentz boost visualisation helpers	✅	✅
Mesh I/O + serialisation	✅ XML/CoHDS	✅ OFF/OBJ/PLY + JSON/XML
Interactive viewer	✅ jReality	✅ libigl/GLFW

Key numerical difference — HyperIdeal Hessian

The analytical Hessian of the HyperIdeal functional requires differentiating through the chain

(b_i, a_e) → l_ij  →  ζ₁₃/ζ₁₄/ζ₁₅  →  α_ij / β_i

which is feasible but involves many nested cases (four vertex-type combinations per edge). Until Phase 6 delivers the analytical version, conformallab++ uses a symmetric finite-difference Hessian:

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:

// 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.

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).

Next steps — global uniformization (cutting closed meshes, period matrices, holonomy), Gauss–Bonnet checking, analytical HyperIdeal Hessian.

---

For mathematicians — extending the library

This section explains how to add new functionals, test conjectures numerically, and hook into the existing solver infrastructure, with no assumed prior knowledge of the codebase.

Mental model

The library is built around one central idea: a discrete conformal functional E(x) whose critical points are the conformally equivalent metrics. Everything else is infrastructure for evaluating E, its gradient G = ∂E/∂x, and its Hessian H = ∂²E/∂x².

ConformalMesh           — half-edge mesh (CGAL::Surface_mesh)
    + property maps     — per-vertex / per-edge data (λ, θ, α, DOF index, …)

Maps struct             — collects all property maps for one functional
    theta_v[v]          — target angle at vertex v  (your input)
    v_idx[v]            — DOF index, or −1 if pinned
    e_idx[e]            — DOF index for edge DOFs (HyperIdeal only)

x ∈ ℝⁿ                 — the DOF vector the solver optimises

evaluate_*(mesh, x, maps)  → { energy, gradient, … }
newton_*(mesh, x0, maps)   → { x*, iterations, converged, … }

The mesh geometry (vertex positions) is only used to initialise the log edge-lengths λ°. From then on the solver works entirely in the x-space.

Adding a new functional — step-by-step

Copy euclidean_functional.hpp as a template (it is the simplest of the three). You need to provide:

1. A Maps struct that holds the property maps your functional needs:

// my_functional.hpp
#pragma once
#include "conformal_mesh.hpp"

namespace conformallab {

struct MyMaps {
    // property maps attached to the mesh
    ConformalMesh::Property_map<Vertex_index, double> lambda;   // log edge-lengths
    ConformalMesh::Property_map<Vertex_index, double> theta_v;  // target angles
    ConformalMesh::Property_map<Vertex_index, int>    v_idx;    // DOF indices

    // any extra parameters your functional needs
    double my_parameter = 1.0;
};

inline MyMaps setup_my_maps(ConformalMesh& mesh) { … }

2. An energy + gradient function:

struct MyResult {
    double              energy;
    std::vector<double> gradient;
};

inline MyResult evaluate_my_functional(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const MyMaps&              m,
    bool                       compute_energy = true)
{
    MyResult res;
    res.gradient.assign(x.size(), 0.0);

    for (auto f : mesh.faces()) {
        // iterate halfedges around face
        // compute your per-face contribution to E and G
        // accumulate: res.gradient[m.v_idx[v]] += …
    }

    // subtract target-angle term
    for (auto v : mesh.vertices()) {
        int iv = m.v_idx[v];
        if (iv < 0) continue;
        res.gradient[iv] -= m.theta_v[v];   // G_v = actual - target
    }

    return res;
}

3. A gradient check — before trusting your formula, verify it numerically. There is a ready-made helper in hyper_ideal_functional.hpp you can call directly, or write your own:

// Finite-difference gradient check for any functional
bool my_gradient_check(ConformalMesh& mesh,
                       const std::vector<double>& x,
                       const MyMaps& m,
                       double eps = 1e-6, double tol = 1e-5)
{
    auto r0 = evaluate_my_functional(mesh, x, m, false);
    const int n = static_cast<int>(x.size());
    for (int i = 0; i < n; ++i) {
        auto xp = x; xp[i] += eps;
        auto xm = x; xm[i] -= eps;
        double fd = (evaluate_my_functional(mesh, xp, m).energy
                   - evaluate_my_functional(mesh, xm, m).energy) / (2*eps);
        if (std::abs(fd - r0.gradient[i]) > tol * (1 + std::abs(fd)))
            return false;
    }
    return true;
}

Add a TEST(MyFunctional, GradientCheck_Triangle) in tests/cgal/ and it will be picked up automatically by CTest.

4. Hook into the Newton solver. Once your gradient and Hessian are correct, plug in solve_linear_system or write a thin wrapper in the style of newton_euclidean:

// Use a numerical Hessian first (safe starting point)
#include "newton_solver.hpp"
#include <Eigen/Sparse>

// Build H by FD of your gradient, then:
bool ok = false;
auto dx = detail::solve_with_fallback(H, -G, ok);

Or supply an analytical Hessian as a sparse matrix and pass it directly.

Where the key mathematical objects live

Object	File	What to look for
Corner angle formula (Euclidean)	euclidean_geometry.hpp	euclidean_corner_angle() — inputs are log half-edge lengths
Spherical angle formula	spherical_geometry.hpp	spherical_corner_angle() — uses spherical law of cosines
HyperIdeal angle (ζ₁₃/ζ₁₄/ζ₁₅)	hyper_ideal_geometry.hpp	zeta13/14/15(), alpha_ij() — the four vertex-type cases
Per-face energy term	*_functional.hpp	the inner loop over mesh.faces()
Gradient accumulation	*_functional.hpp	grad[v_idx[v]] += … after the face loop
Cotangent-Laplace structure	euclidean_hessian.hpp	euclidean_hessian() — shows the sparse-triplet pattern
Spherical Hessian derivation	spherical_hessian.hpp	comments give the ∂α/∂u formula step by step
Special functions	clausen.hpp	Cl2(), lobachevsky(), imLi2() — all take a double angle

How to navigate the half-edge mesh

for (auto f : mesh.faces()) {
    // The three halfedges of face f:
    auto h0 = mesh.halfedge(f);
    auto h1 = mesh.next(h0);
    auto h2 = mesh.next(h1);

    // Vertices opposite to each halfedge (the vertex NOT on h):
    Vertex_index v0 = mesh.target(h2);  // opposite to edge h0-h1
    Vertex_index v1 = mesh.target(h0);  // opposite to edge h1-h2
    Vertex_index v2 = mesh.target(h1);  // opposite to edge h0-h2 (= h2 target)

    // Access DOF index (−1 = pinned):
    int i0 = maps.v_idx[v0];

    // The opposite halfedge (for the adjacent face, if not on boundary):
    auto h_opp = mesh.opposite(h0);
    bool is_boundary = mesh.is_border(h_opp);
}

Attaching new data to a mesh

// Add a per-vertex curvature field (survives mesh copy):
auto [curv, created] = mesh.add_property_map<Vertex_index, double>("v:my_curv", 0.0);

// Write and read:
curv[v] = 1.234;
double k = curv[v];

// Pass it through your Maps struct so functions can access it.

Property maps are reference-counted and cheap to copy. Give them unique string names to avoid collision.

Quick-start experiment checklist

1. Read examples/example_layout.cpp — it shows the full pipeline (load → setup → solve → layout → JSON/XML save → reload) in ~120 lines with comments at every step.
2. Build with cmake -S code -B build -DWITH_CGAL=ON && cmake --build build --target example_layout and run ./build/examples/example_layout.
3. Add a gradient check test in tests/cgal/ — copy any TEST(…, GradientCheck_…) block and swap out the functional. Run with ctest -R your_test_name.
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.

Porting the Java uniformization pipeline — a roadmap for contributors

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:

// 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 <stdexcept>
#include <cmath>

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<int>(mesh.number_of_vertices())
            - static_cast<int>(mesh.number_of_edges())
            + static_cast<int>(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:

// cut_graph.hpp (new file)
#pragma once
#include "conformal_mesh.hpp"
#include <unordered_set>
#include <vector>

namespace conformallab {

struct CutGraph {
    std::unordered_set<std::size_t> 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:

// 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:

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:

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) | | Bobenko, Springborn — Variational Principles for Circle Patterns (2004) | Background for the angle-sum variational framework used throughout |

---

Key design decisions

CGAL as CoHDS replacement. CGAL::Surface_mesh<Point3> replaces the Java CoHDS half-edge data structure. Vertex/edge/face/halfedge descriptors are typed integers — no raw handles, no RTTI.

Property maps. mesh.add_property_map<Vertex_index, double>("v:lambda", 0.0) replaces the Java adapter/decorator pattern. Multiple maps attach to one mesh without subclassing.

DOF vector convention. All functionals use x indexed by v_idx[v] / e_idx[e] (−1 = pinned). This matches the Java FunctionalTest gradient-check convention and is uniform across all three geometries.

HyperIdeal Hessian via FD. The analytical Hessian through ζ13/14/15 → lij → β/α is deferred to Phase 6. A symmetric FD Hessian H[i,j] = (G(x+ε·eⱼ)[i] − G(x−ε·eⱼ)[i]) / (2ε) is O(ε²) accurate, PSD by strict convexity, and sufficient for Newton on < 500 DOFs.

Spherical Hessian sign. The spherical energy is concave (not convex) — the Hessian H is NSD at equilibrium. Newton solves (−H)·Δx = G, so the sign flip is handled transparently inside newton_spherical.

Natural theta trick. Tests set theta_v = Σα_v(x=x_base) to make x_base the known equilibrium, avoiding any need to manufacture reference solutions. For HyperIdeal x_base = (b=1.0, a=0.5) is used (x=0 is degenerate in log-space).

---

CI

Tests run automatically on push to main, dev, and claude/** branches via a self-hosted Gitea Actions runner (eulernest, ARM64). The CI image contains cmake, g++, git, and Node.js. Only conformallab_tests runs in CI (no Boost/CGAL dependency there).

# Rebuild and push the CI image when the Dockerfile changes
docker buildx build \
  --platform linux/arm64 \
  -f .gitea/docker/Dockerfile.ci-cpp \
  -t git.eulernest.eu/conformallab/ci-cpp:latest \
  --push \
  .gitea/docker/

---

Roadmap

Phase 1   Clausen / Lobachevsky / ImLi₂                ✅ abgeschlossen

Phase 2   Hyper-ideal Geometrie (ζ, lᵢⱼ, αᵢⱼ, σᵢ)    ✅ abgeschlossen

Phase 3a  CGAL Surface_mesh Infrastruktur              ✅ abgeschlossen
Phase 3b  HyperIdealFunctional                         ✅ abgeschlossen
Phase 3c  SphericalFunctional                          ✅ abgeschlossen
Phase 3d  EuclideanCyclicFunctional                    ✅ abgeschlossen
Phase 3e  Gauge-Fix für SphericalFunctional            ✅ abgeschlossen
Phase 3f  Analytische Hessians (Eucl. + Sphär.)        ✅ abgeschlossen
Phase 3g  PI-Konstante konsolidieren                   ✅ abgeschlossen

Phase 4a  Newton-Solver (alle drei Geometrien)         ✅ abgeschlossen
          → newton_euclidean / newton_spherical / newton_hyper_ideal
          → detail::solve_with_fallback → public solve_linear_system
          → Backtracking-Line-Search
          → hyper_ideal_hessian.hpp (numerischer FD-Hessian)

Phase 4b  CGAL::IO Mesh-Import/Export                  ✅ abgeschlossen
          → mesh_io.hpp: read/write/load/save
          → Format-Erkennung aus Dateiendung (OFF, OBJ, PLY)

Phase 4c  End-to-End-Pipeline Tests                    ✅ abgeschlossen
          → test_pipeline.cpp: 5 Tests (alle 3 Geometrien, I/O, full loop)

Phase 4d  SparseQR-Fallback + Beispiel-Programme       ✅ abgeschlossen
          → solve_linear_system als öffentliche API mit fallback_used-Flag
          → 3 dedizierte SparseQR-Tests (full-rank, singular, closed mesh)
          → examples/example_euclidean.cpp   (headless)
          → examples/example_hyper_ideal.cpp (headless)
          → examples/example_viewer.cpp      (interaktiver Viewer, WITH_VIEWER)

Phase 5   Layout + CLI + Serialisierung                ✅ abgeschlossen
          → layout.hpp: BFS-Einbettung in ℝ² (Euclidean/HyperIdeal) und S² (Spherical)
          → serialization.hpp: JSON (nlohmann/json) + XML (hand-written) save/load
          → conformallab_core CLI: -i/-o/-g/-j/-x/-s/-v Flags, alle drei Geometrien
          → example_layout.cpp: Solve → Layout → OFF/JSON/XML + Round-Trip-Check
          → test_layout.cpp: 8 Tests (Eucl./Sphär./HyperIdeal + JSON/XML)
          → 95 Tests gesamt (2 skipped)

Phase 6   (geplant)
          → Analytischer HyperIdeal-Hessian (direkte Ableitung durch ζ-Kette)
          → Gauss–Bonnet Konsistenzprüfung für Kegelmetriken
          → Mesh-Cut für geschlossene Flächen → globale Parameterisierung
          → Inversive-Distance-Funktional (Luo 2004)

---

License

conformallab++ is released under the MIT License (see LICENSE).