ConformalLabpp/doc/architecture/overall_pipeline.md

# conformallab++ — Architecture & Pipeline

## Positioning

conformallab++ is a **specialised research library for discrete conformal geometry** on
triangulated surfaces. It is not a general geometry-processing framework.

The library revives the algorithms of the original ConformalLab by Stefan Sechelmann
in a modern C++ setting. Its core concern is one precise mathematical question:

> Given a triangulated surface, find a conformally equivalent metric that satisfies
> prescribed curvature (angle-sum) constraints at each vertex.

Everything in the library serves this goal:

| What it is | What it is not |
|------------|----------------|
| Discrete conformal maps (Euclidean, spherical, hyperbolic) | General mesh processing |
| Newton solver for angle-sum energy functionals | Remeshing / boolean / smoothing |
| Priority-BFS layout into ℝ², S², Poincaré disk | Point-cloud or implicit-field processing |
| Holonomy, period matrices, fundamental domains | NURBS / parametric modelling |
| Research platform — experiment-first, CLI-ready | Production rendering engine |

**Relation to existing libraries.**  
conformallab++ sits *on top of* CGAL, not beside it.
`CGAL::Surface_mesh` is the internal mesh representation; there is no additional
abstraction layer. Eigen handles all linear algebra. libigl provides the optional
interactive viewer. The library adds the conformal-geometry layer that none of these
provide.

---

## The conformal geometry pipeline

The full pipeline runs in three stages. All stages operate on the same
`ConformalMesh` (= `CGAL::Surface_mesh<Point3>` with attached property maps).

```mermaid
graph TD
    IN["Input\n(OFF / OBJ / PLY / builder)"]
    ADAPTER["Input Adapter\nload_mesh() · make_triangle() · make_quad_strip() …"]

    subgraph PRE["① PREPROCESSING"]
        MAPS["Maps Setup\nsetup_*_maps() + compute_lambda0_from_mesh()"]
        DOF["DOF Assignment\nv_idx[v] · e_idx[e] · pin / free"]
        THETA["Target Angles\ntheta_v[v] — cone metric or natural equilibrium"]
        GB["Gauss–Bonnet Check\ncheck_gauss_bonnet() / enforce_gauss_bonnet()"]
        MAPS --> DOF --> THETA --> GB
    end

    subgraph CORE["② PROCESSING CORE"]
        NEWTON["Newton Solver\nnewton_euclidean/spherical/hyper_ideal()\n→ NewtonResult  x* ∈ ℝⁿ"]
        CG["Cut Graph  [closed surfaces]\ncompute_cut_graph()\n→ CutGraph  (2g seam edges)"]
        LAYOUT["Layout / Embedding\neuclidean/spherical/hyper_ideal_layout()\n→ Layout2D / Layout3D\n  · uv[v]  · halfedge_uv[h]\n  · HolonomyData"]
        NORM["Normalisation\nnormalise_euclidean / _hyperbolic / _spherical()"]
        PERIOD["Period Matrix  [genus 1]\ncompute_period_matrix()\n→ PeriodData  τ ∈ ℍ"]
        FD["Fundamental Domain\ncompute_fundamental_domain()\n→ FundamentalDomain + tiling"]
        NEWTON --> CG --> LAYOUT --> NORM --> PERIOD --> FD
    end

    subgraph POST["③ POSTPROCESSING"]
        SERIAL["Serialisation\nsave_result_json/xml()\nsave_layout_off()"]
        VIZ["Visualisation\nexample_viewer  (libigl / GLFW)"]
        CLI["CLI App\nconformallab_core  -i -g -o -j -x -s"]
    end

    IN --> ADAPTER --> PRE
    PRE --> CORE
    CORE --> POST

    style PRE  fill:#dbeafe,stroke:#3b82f6,stroke-width:2px,color:#000
    style CORE fill:#dcfce7,stroke:#16a34a,stroke-width:2px,color:#000
    style POST fill:#fef9c3,stroke:#ca8a04,stroke-width:2px,color:#000
```

---

## Stage ① — Preprocessing

### Input adapter

All mesh input flows through a single entry point:

```cpp
ConformalMesh mesh = load_mesh("input.off");   // OFF · OBJ · PLY
ConformalMesh mesh = make_quad_strip();         // built-in test meshes
```

**Precondition guarantee:** the adapter ensures the mesh is a valid, orientable,
triangulated surface with consistent halfedge structure (CGAL validity).
Downstream stages never check for degeneracies — they trust the adapter.

### Maps setup

Each geometry mode has its own maps struct that attaches property maps to the mesh:

| Geometry | Setup function | Key property maps |
|----------|---------------|-------------------|
| Euclidean (ℝ²) | `setup_euclidean_maps()` | `lambda0[e]`, `theta_v[v]`, `v_idx[v]` |
| Spherical (S²) | `setup_spherical_maps()` | `lambda0[e]`, `theta_v[v]`, `v_idx[v]` |
| Hyper-ideal (H²) | `setup_hyper_ideal_maps()` | `beta_v[v]`, `alpha_e[e]`, `v_idx[v]`, `e_idx[e]` |

`compute_*_lambda0_from_mesh()` initialises log-edge-lengths from the 3-D vertex
positions. After this step the solver works entirely in scale-factor space; the
original vertex positions are no longer needed.

### DOF assignment & target angles

```cpp
// Pin first vertex (gauge fix for open meshes)
maps.v_idx[*mesh.vertices().begin()] = -1;
int idx = 0;
for (auto v : rest_of_vertices) maps.v_idx[v] = idx++;

// Set target curvature (cone metric or natural equilibrium)
maps.theta_v[v] = 2 * M_PI;          // flat interior vertex
maps.theta_v[v] = M_PI / 3;          // 60° cone singularity
```

**Natural equilibrium shortcut:** evaluate the gradient at `x = 0`, subtract it from
`theta_v` — the solver then converges to `x* = 0` identically (useful for tests).

### Gauss–Bonnet check

Before solving, the prescribed angles must satisfy:

$$\sum_{v} (2\pi - \Theta_v) = 2\pi \cdot \chi(M)$$

```cpp
check_gauss_bonnet(mesh, maps);         // throws if violated
enforce_gauss_bonnet(mesh, maps);       // distributes residual uniformly
```

**This is the most common source of silent Newton non-convergence.**
Prescribing angles that violate Gauss–Bonnet means no conformal factor exists —
the solver will iterate without converging.

---

## Stage ② — Processing core

### The three geometry modes

conformallab++ implements discrete conformal geometry in three model spaces:

| Mode | Space | Curvature | Typical surfaces |
|------|-------|-----------|-----------------|
| **Euclidean** | ℝ² | K = 0 | Flat tori, developable surfaces, open patches |
| **Spherical** | S² | K = +1 | Genus-0 (sphere-like) surfaces |
| **Hyper-ideal** | H² (Poincaré disk) | K = −1 | Genus-g surfaces (g ≥ 1), hyperbolic structures |

All three share the same algorithmic structure — only the angle formula, the
trilateration geometry, and the Hessian sign differ.

### Newton solver

```
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** (angle-sum defect per vertex)
2. Evaluate Hessian **H** (analytical for Euclidean/Spherical; symmetric FD for HyperIdeal)
3. Solve **H·Δx = −G** — try `SimplicialLDLT`, fall back to `SparseQR` on rank deficiency
4. Backtracking line search (up to 20 halvings)

**Preconditions:** mesh triangulated + manifold · Gauss–Bonnet satisfied · DOFs assigned  
**Provides:** `NewtonResult.x` — converged scale factors; `converged`, `iterations`, `grad_inf_norm`

**SparseQR fallback** handles gauge modes on closed meshes without a pinned vertex.
The fallback is public API: `solve_linear_system(H, rhs, &used_fallback)`.

### Cut graph (closed surfaces)

For a closed genus-g surface, cutting along 2g independent homology cycles
turns it into a topological disk — a prerequisite for globally consistent BFS layout.

```cpp
CutGraph cg = compute_cut_graph(mesh);
// cg.cut_edge_flags[e.idx()] — true for seam edges
// cg.cut_edge_indices         — ordered list of 2g seam edges
// cg.genus                    — g
```

**Algorithm:** tree-cotree decomposition (Erickson–Whittlesey 2005).  
**Preconditions:** closed, orientable, triangulated mesh  
**Provides:** exactly `2g` seam edges whose removal makes the surface simply connected

### Layout / Embedding

BFS-trilateration unfolds the mesh into the target geometry.
The root face (largest 3-D area, 1.5× interior bonus) is placed first;
all other faces are placed in **priority order by BFS depth** (min-heap),
minimising trilateration error accumulation.

```cpp
HolonomyData hol;
Layout2D layout = euclidean_layout(mesh, result.x, maps, &cg, &hol, /*normalise=*/true);
Layout3D slayout = spherical_layout(mesh, result.x, smaps);
Layout2D hlayout = hyper_ideal_layout(mesh, result.x, hmaps, &cg, &hol);
```

**Key outputs:**

| Field | Content |
|-------|---------|
| `layout.uv[v.idx()]` | Primary UV — first / shallowest-BFS visit |
| `layout.halfedge_uv[h.idx()]` | UV of `source(h)` as seen from `face(h)` — seam-aware |
| `layout.has_seam` | True when a vertex was reached via two different paths |
| `hol.translations[i]` | Translation ω_i across cut edge i (Euclidean / spherical) |
| `hol.mobius_maps[i]` | Möbius isometry T_i ∈ SU(1,1) across cut edge i (hyperbolic) |

**`halfedge_uv` — texture atlas semantics:**  
At a seam edge the two opposite halfedges carry *different* UV values.
This gives each face its own UV copy of a seam vertex, enabling
a proper GPU texture atlas without vertex duplication.

**Trilateration — geometry by mode:**

| Mode | Method | Accuracy |
|------|--------|---------|
| Euclidean | Analytic formula in ℝ² | Exact |
| Spherical | Spherical law of cosines on S² | Exact |
| Hyper-ideal | Möbius + hyperbolic law of cosines in Poincaré disk | Exact |

**Preconditions:** `NewtonResult.converged` · mesh · maps · optional `CutGraph`  
**Provides:** `Layout2D/3D` with `uv`, `halfedge_uv` · `HolonomyData` with translations / Möbius maps

### Möbius maps

`MobiusMap` (T(z) = (az+b)/(cz+d)) is the central algebraic object for hyperbolic geometry:

```cpp
MobiusMap T = MobiusMap::from_three(z1,w1, z2,w2, z3,w3); // fit to 3 correspondences
MobiusMap S = T.inverse().compose(U);                       // group operations
Eigen::Vector2d p2 = T.apply(p);                           // apply to 2-D point
```

Used for: hyperbolic trilateration · holonomy tracking · normalisation centering.

### Normalisation

After layout, a canonical post-processing step brings the result into a standard position:

| Mode | Method | Effect |
|------|--------|--------|
| Euclidean | PCA — centroid → origin, major axis → x-axis | Translation + rotation |
| Hyperbolic | Weighted Möbius centering (Fréchet mean, 30 iterations) | Maps centroid to disk origin |
| Spherical | Rodrigues rotation | Maps centroid to north pole |

Both `uv` and `halfedge_uv` are transformed identically.

### Period matrix (genus 1)

From the two holonomy translations ω₁, ω₂ ∈ ℂ read off from the cut graph,
the conformal type of a flat torus is the SL(2,ℤ)-orbit of:

$$\tau = \omega_2 / \omega_1 \in \mathbb{H}$$

```cpp
PeriodData pd = compute_period_matrix(hol);
// pd.tau                   — complex period ratio
// pd.omega[i]              — lattice generators as complex numbers
// pd.in_fundamental_domain — after SL(2,ℤ) reduction
// pd.genus()               — g = omega.size() / 2
```

SL(2,ℤ) reduction alternates T: τ↦τ+1 and S: τ↦−1/τ steps until
τ ∈ F = {|τ| ≥ 1, −½ ≤ Re(τ) < ½}.

**Note:** The Siegel period matrix Ω ∈ H_g for genus g ≥ 2 requires integrating
holomorphic differentials — deferred to Phase 8.

### Fundamental domain

```cpp
FundamentalDomain fd = compute_fundamental_domain(hol);
// genus 1: CCW parallelogram {0, ω₁, ω₁+ω₂, ω₂}
// genus g > 1: empty — 4g-polygon boundary walk deferred to Phase 8

auto tiles = tiling_neighbourhood(layout, hol, /*m_max=*/2, /*n_max=*/2);
// returns (2·m_max+1)·(2·n_max+1) translated copies of the layout
// for visualising the universal cover
```

---

## Stage ③ — Postprocessing

### Serialisation

```cpp
// Layout as mesh file
save_layout_off("layout.off", mesh, layout);

// Full result: DOF vector + metadata + layout UVs
save_result_json("result.json", result, "euclidean", V, F, &layout);
save_result_xml ("result.xml",  result, "euclidean", V, F, &layout);

// Round-trip load
NewtonResult res2; std::string geom; Layout2D uv2;
load_result_json("result.json", &res2, &geom, &uv2);
```

### CLI app

```bash
conformallab_core \
  -i input.off        # input mesh
  -g euclidean        # geometry: euclidean | spherical | hyper_ideal
  -o layout.off       # layout output
  -j result.json      # JSON serialisation
  -x result.xml       # XML serialisation
  -s                  # show input in viewer
  -v                  # verbose solver output
```

### Interactive viewer

`example_viewer` (libigl / GLFW) shows the 3-D mesh and the 2-D layout
side-by-side. Built automatically with `-DWITH_CGAL=ON`.

---

## Processing unit contracts

Each stage has explicit preconditions and guarantees.
A pipeline is valid if every unit's preconditions are satisfied
by the outputs of all preceding units.

| Unit | Preconditions | Provides |
|------|--------------|---------|
| `load_mesh` | Valid file path, supported format | Manifold, oriented, triangulated `ConformalMesh` |
| `setup_*_maps` | Triangulated mesh | Initialised property maps; `lambda0` from 3-D positions |
| `check_gauss_bonnet` | `theta_v` set | Throws if Σ(2π−Θ_v) ≠ 2π·χ |
| `enforce_gauss_bonnet` | `theta_v` set | Σ(2π−Θ_v) = 2π·χ guaranteed |
| `newton_*` | GB satisfied · DOFs assigned | `NewtonResult.converged` · `x*` · gradient norm |
| `compute_cut_graph` | Closed, orientable mesh | `2g` seam edges · `CutGraph.genus` |
| `euclidean_layout` | `newton_euclidean` converged | `uv[v]` · `halfedge_uv[h]` · `HolonomyData` |
| `normalise_euclidean` | `layout.success == true` | Centroid at origin · major axis = x-axis |
| `compute_period_matrix` | `hol.translations.size() >= 2` | `τ ∈ ℍ` · optionally reduced to F |
| `compute_fundamental_domain` | `HolonomyData` (genus 1) | CCW parallelogram · edge identifications |

---

## The three geometry modes in detail

```
              Euclidean            Spherical            Hyper-ideal
              ─────────────────────────────────────────────────────
Space         ℝ²                   S²                   H² (Poincaré disk)
Curvature     K = 0                K = +1               K = −1
Genus         any (cone metrics)   0                    ≥ 1
Angle sum     Σα_v = Θ_v           Σα_v = Θ_v           Σβ_v = Θ_v
Hessian       PSD (cotangent-Lap.) NSD (sign-flip)      PSD (FD, strict conv.)
Holonomy      translations ω_i     rotations (2-D)      Möbius maps T_i ∈ SU(1,1)
Period        τ = ω₂/ω₁ ∈ ℍ        —                    axis of T_i
Normalise     PCA centring         Rodrigues to N pole  weighted Möbius centring
```

---

## Extension points

### Adding a new functional

1. Create `my_functional.hpp` with a `Maps` struct and `evaluate_my_functional()`.
2. The gradient must satisfy: `G_v = Σ(angle contributions) − theta_v[v]`.
3. Verify with a finite-difference gradient check (copy any `GradientCheck_*` test).
4. Pass to `solve_linear_system(H, -G)` or write a thin `newton_my` wrapper.

### Adding a new geometry mode

Implement `trilaterate_my()` with the correct local placement formula,
then follow the same BFS structure as `euclidean_layout` (see `layout.hpp`).
The priority BFS, holonomy tracking, and `halfedge_uv` population are geometry-agnostic.

### Adding a new processing unit

Declare its preconditions and capabilities explicitly (see table above).
The pipeline validation is currently manual (documented contracts); a
compile-time or runtime check is a natural Phase 8 extension.

---

## Declarative pipeline (target for Phase 8)

A lightweight YAML description for reproducible experiments:

```yaml
pipeline:
  name: flat_torus_period
  geometry: euclidean

  input:
    source: data/torus.off

  steps:
    - id: setup
      unit: setup_euclidean_maps
      provide: [maps_initialised]

    - id: gauss_bonnet
      unit: enforce_gauss_bonnet
      require: [maps_initialised]
      provide: [gauss_bonnet_satisfied]

    - id: solve
      unit: newton_euclidean
      require: [gauss_bonnet_satisfied]
      params:
        tol: 1.0e-10
        max_iter: 200
      provide: [x_converged]

    - id: cut
      unit: compute_cut_graph
      require: [mesh_closed]
      provide: [cut_graph]

    - id: layout
      unit: euclidean_layout
      require: [x_converged, cut_graph]
      params:
        normalise: true
      provide: [layout_uv, holonomy]

    - id: period
      unit: compute_period_matrix
      require: [holonomy]
      provide: [tau]

  output:
    layout:  out/torus_layout.off
    json:    out/torus_result.json
    tau:     out/torus_tau.txt
```

Each unit declares `require` (preconditions) and `provide` (capabilities).
A pipeline is valid if for every step, all `require` keys are satisfied by the
accumulated `provide` set of all preceding steps.
This mirrors exactly the contract table in this document.

---

## Recommended reading

| Source | Relevance in conformallab++ |
|--------|----------------------------|
| Springborn — *Ideal Hyperbolic Polyhedra and Discrete Uniformization* (2020) | HyperIdeal functional; ζ₁₃/ζ₁₄/ζ₁₅ in `hyper_ideal_geometry.hpp` |
| Pinkall, Polthier — *Computing Discrete Minimal Surfaces* (1993) | Cotangent-Laplace Hessian in `euclidean_hessian.hpp` |
| Bobenko, Springborn — *Variational Principles for Circle Patterns* (2004) | Angle-sum variational framework used throughout |
| Luo — *Combinatorial Yamabe Flow on Surfaces* (2004) | Inversive-distance functional (not yet ported) |
| Erickson, Whittlesey — *Greedy Optimal Homotopy Generators* (SODA 2005) | Tree-cotree algorithm in `cut_graph.hpp` |