diff --git a/README.md b/README.md index 82e9637..08bb570 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 4 vollständig abgeschlossen. Alle drei Geometrien lösbar via Newton-Solver (SimplicialLDLT + SparseQR-Fallback). Interaktiver Viewer-Beispiel inklusive. **87 Tests, 2 skipped**. Nächster Schritt: Phase 5 (CLI-App + Serialisierung). +> **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**. --- @@ -26,29 +26,43 @@ The long-term goal is a **CGAL package** that brings discrete conformal maps (hy | Mesh I/O (CGAL::IO — OFF / OBJ / PLY) | ✅ Phase 4b | | End-to-end pipeline tests | ✅ Phase 4c | | **Example programs** (headless + interactive viewer) | ✅ Phase 4d | -| CLI app + XML serialisation | 🔜 Phase 5 | +| **BFS Layout** (ℝ², S², Poincaré disk) | ✅ Phase 5 | +| **CLI app** (`conformallab_core`) | ✅ Phase 5 | +| **JSON + XML serialisation** | ✅ Phase 5 | --- -## Quick start — running the examples +## Quick start — CLI app ```bash cmake -S code -B build -DWITH_CGAL=ON -cmake --build build --target example_euclidean example_hyper_ideal +cmake --build build -j4 -# Euclidean conformal map (headless) -./build/examples/example_euclidean [input.off] [output.off] +# 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 -# HyperIdeal conformal map (headless) -./build/examples/example_hyper_ideal [input.off] [output.off] +# Available geometries: euclidean | spherical | hyper_ideal +./bin/conformallab_core -i input.off -g spherical -o sphere.off -# Interactive viewer (requires -DWITH_VIEWER=ON) -cmake -S code -B build -DWITH_CGAL=ON -DWITH_VIEWER=ON -cmake --build build --target example_viewer +# Show input mesh in interactive viewer (built-in) +./bin/conformallab_core -i input.off -s +``` + +### Example programs + +```bash +# 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 test mesh. +If no input file is given each example uses a built-in quad-strip mesh. --- @@ -107,6 +121,39 @@ int n = assign_all_dof_indices(mesh, maps); auto result = newton_hyper_ideal(mesh, x0, maps); ``` +### Layout (embedding into target space) + +After solving, convert scale factors into actual vertex coordinates: + +```cpp +#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): ```cpp @@ -125,8 +172,8 @@ if (used_fallback) | 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` | local only | -| **Interactive viewer** | `-DWITH_CGAL=ON -DWITH_VIEWER=ON` | above + `example_viewer`, `conformallab_core` | local only | +| **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`). @@ -163,19 +210,21 @@ ctest --test-dir build --output-on-failure ```bash cmake -S code -B build -DWITH_CGAL=ON -cmake --build build --target conformallab_cgal_tests example_euclidean example_hyper_ideal -j$(nproc) +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: **87 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). +Expected: **95 tests pass, 2 skipped** (the two `@Ignore` Hessian stubs). ### Interactive viewer ```bash -cmake -S code -B build -DWITH_CGAL=ON -DWITH_VIEWER=ON -cmake --build build --target example_viewer -j$(nproc) +# 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 ``` @@ -201,6 +250,8 @@ cmake --build build --target example_viewer -j$(nproc) | `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 | +| `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` | @@ -216,6 +267,8 @@ 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) +│ ├── 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 @@ -225,9 +278,10 @@ code/ │ ├── 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 # CLI skeleton (Phase 5) +│ ├── apps/v0/conformallab_cli.cpp # Full CLI app (Phase 5) │ └── viewer/simple_viewer.cpp ├── tests/ │ ├── CMakeLists.txt @@ -242,7 +296,8 @@ code/ │ ├── 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_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 @@ -277,7 +332,9 @@ Pure-math tests requiring only Eigen: Clausen / Lobachevsky / ImLi₂, hyper-ide | `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 | -| **Total** | **87** | 2 skipped (Hessian stubs) | +| `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) | --- @@ -334,7 +391,7 @@ The three core functionals are fully equivalent to the Java original at the leve | 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² | ✅ | ❌ not implemented | +| **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 | ✅ | ❌ | @@ -343,7 +400,7 @@ The three core functionals are fully equivalent to the Java original at the leve | Clausen / Lobachevsky / ImLi₂ special functions | ✅ | ✅ | | Discrete elliptic utility (modular normalisation of τ) | ✅ | ✅ (not yet wired up) | | Poincaré disk / Lorentz boost visualisation helpers | ✅ | ✅ | -| Mesh I/O | ✅ XML/CoHDS | ✅ OFF/OBJ/PLY | +| Mesh I/O + serialisation | ✅ XML/CoHDS | ✅ OFF/OBJ/PLY + JSON/XML | | Interactive viewer | ✅ jReality | ✅ libigl/GLFW | ### Key numerical difference — HyperIdeal Hessian @@ -354,7 +411,7 @@ The analytical Hessian of the HyperIdeal functional requires differentiating thr (b_i, a_e) → l_ij → ζ₁₃/ζ₁₄/ζ₁₅ → α_ij / β_i ``` -which is feasible but involves many nested cases (four vertex-type combinations per edge). Until Phase 5 delivers the analytical version, conformallab++ uses a symmetric finite-difference Hessian: +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ε) @@ -362,11 +419,13 @@ 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. -### What "cone metrics" and "layout" would require +### 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** — after solving you have the conformal scale factors u_i but the vertex positions in the mesh are unchanged. Recovering the actual flat / hyperbolic / spherical coordinates requires integrating the discrete holomorphic differential (discrete Schwarz–Christoffel for the Euclidean case, or geodesic development for the hyperbolic case). This is the single biggest missing step for a complete uniformization pipeline. +**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. --- @@ -547,8 +606,8 @@ Property maps are reference-counted and cheap to copy. Give them unique string n ### Quick-start experiment checklist -1. **Read** `examples/example_euclidean.cpp` — it shows the full pipeline in ~80 lines with comments at every step. -2. **Build** without a viewer first: `cmake -S code -B build -DWITH_CGAL=ON && cmake --build build --target example_euclidean`. +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. @@ -573,7 +632,7 @@ Property maps are reference-counted and cheap to copy. Give them unique string n **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 left for Phase 5. 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. +**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`. @@ -632,11 +691,19 @@ Phase 4d SparseQR-Fallback + Beispiel-Programme ✅ abgeschlossen → examples/example_hyper_ideal.cpp (headless) → examples/example_viewer.cpp (interaktiver Viewer, WITH_VIEWER) -Phase 5 CLI-App + Serialisierung (2–3 Tage) - → conformallab_core CLI: --input, --geometry, --output - → JSON/XML-Export für DOF-Vektoren und Solver-Ergebnisse - → Analytischer HyperIdeal-Hessian (direkte Ableitung) - → Integration aller drei Geometrien hinter einheitlicher CLI +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) ``` --- diff --git a/code/examples/CMakeLists.txt b/code/examples/CMakeLists.txt index a338815..8b9c724 100644 --- a/code/examples/CMakeLists.txt +++ b/code/examples/CMakeLists.txt @@ -13,6 +13,7 @@ set(EXAMPLE_INCLUDES ${CMAKE_SOURCE_DIR}/deps/eigen-3.4.0 ${CMAKE_SOURCE_DIR}/deps/CGAL-6.1.1/include + ${CMAKE_SOURCE_DIR}/deps/single_includes ${Boost_INCLUDE_DIRS} ) set(EXAMPLE_PRIVATE_INCLUDES @@ -35,6 +36,12 @@ target_include_directories(example_hyper_ideal SYSTEM PRIVATE ${EXAMPLE_INCLUDES target_include_directories(example_hyper_ideal PRIVATE ${EXAMPLE_PRIVATE_INCLUDES}) target_compile_definitions(example_hyper_ideal PRIVATE ${EXAMPLE_DEFS}) +# ── example_layout ─────────────────────────────────────────────────────────── +add_executable(example_layout example_layout.cpp) +target_include_directories(example_layout SYSTEM PRIVATE ${EXAMPLE_INCLUDES}) +target_include_directories(example_layout PRIVATE ${EXAMPLE_PRIVATE_INCLUDES}) +target_compile_definitions(example_layout PRIVATE ${EXAMPLE_DEFS}) + # ── example_viewer (requires WITH_VIEWER) ───────────────────────────────────── if(WITH_VIEWER) add_executable(example_viewer example_viewer.cpp) diff --git a/code/examples/example_layout.cpp b/code/examples/example_layout.cpp new file mode 100644 index 0000000..dc8b210 --- /dev/null +++ b/code/examples/example_layout.cpp @@ -0,0 +1,168 @@ +// example_layout.cpp +// +// Phase 5 example — Euclidean conformal layout with JSON/XML serialisation. +// +// Demonstrates the full pipeline that a library user would run: +// +// 1. Build (or load) a mesh. +// 2. Set up Euclidean conformal maps + compute λ° from vertex positions. +// 3. Choose natural target angles (→ x* = 0 is the equilibrium). +// 4. Run Newton's method. +// 5. Unfold the mesh in the plane (euclidean_layout). +// 6. Save the layout as an OFF file for inspection in MeshLab/Blender. +// 7. Serialise the solver result and UV coordinates to JSON and XML. +// 8. Reload from JSON and verify the DOF vector is recovered. +// +// Build (from code/): +// cmake -S . -B build -DWITH_CGAL=ON && cmake --build build -t example_layout +// +// Run: +// ./build/examples/example_layout +// ./build/examples/example_layout input.off layout.off result.json result.xml + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "newton_solver.hpp" +#include "layout.hpp" +#include "serialization.hpp" + +#include +#include +#include +#include +#include + +using namespace conformallab; + +// ── Helper: natural target angles so x* = 0 is the equilibrium ─────────────── +static void set_natural_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n) +{ + std::vector x0(static_cast(n), 0.0); + auto G = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] -= G[static_cast(iv)]; + } +} + +// ── Helper: pin vertex 0, assign DOF indices 0..n-1 to the rest ────────────── +static int pin_first(ConformalMesh& mesh, EuclideanMaps& maps) +{ + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; // pinned + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + return idx; +} + +// ── Main ───────────────────────────────────────────────────────────────────── + +int main(int argc, char* argv[]) +{ + // ── Step 1: mesh ────────────────────────────────────────────────────────── + ConformalMesh mesh; + std::string out_off = "layout.off"; + std::string out_json = "result.json"; + std::string out_xml = "result.xml"; + + if (argc >= 2) { + // User provided an input mesh + if (!read_mesh(argv[1], mesh) || mesh.is_empty()) { + std::cerr << "Cannot load " << argv[1] << "\n"; + return 1; + } + if (argc >= 3) out_off = argv[2]; + if (argc >= 4) out_json = argv[3]; + if (argc >= 5) out_xml = argv[4]; + } else { + // Use built-in quad-strip (6 vertices, 4 triangles) + mesh = make_quad_strip(); + std::cout << "Using built-in quad-strip mesh (no arguments given).\n"; + } + + std::cout << "Mesh: " << mesh.number_of_vertices() << " vertices, " + << mesh.number_of_faces() << " faces\n"; + + // ── Step 2: maps + λ° ──────────────────────────────────────────────────── + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // ── Step 3: DOF assignment + natural targets ────────────────────────────── + int n = pin_first(mesh, maps); + std::cout << "DOFs: " << n << " (vertex 0 pinned)\n"; + + set_natural_theta(mesh, maps, n); + + // ── Step 4: Newton ──────────────────────────────────────────────────────── + std::vector x0(static_cast(n), 0.0); + auto res = newton_euclidean(mesh, x0, maps); + + std::cout << std::boolalpha + << "Newton converged: " << res.converged + << " iterations: " << res.iterations + << " |grad|_inf: " + << std::scientific << std::setprecision(3) << res.grad_inf_norm << "\n"; + + // Print scale factors u_v + std::cout << "Conformal factors u_v:\n"; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + double u = (iv >= 0) ? res.x[static_cast(iv)] : 0.0; + std::cout << " v" << v.idx() << ": u = " << std::fixed << std::setprecision(8) << u << "\n"; + } + + // ── Step 5: Layout ──────────────────────────────────────────────────────── + auto layout = euclidean_layout(mesh, res.x, maps); + + std::cout << "Layout success: " << layout.success + << " has_seam: " << layout.has_seam << "\n"; + + if (layout.success) { + std::cout << "UV coordinates:\n"; + for (auto v : mesh.vertices()) { + auto& p = layout.uv[v.idx()]; + std::cout << " v" << v.idx() + << ": (" << std::fixed << std::setprecision(6) + << p.x() << ", " << p.y() << ")\n"; + } + } + + // ── Step 6: Save layout OFF ─────────────────────────────────────────────── + save_layout_off(out_off, mesh, layout); + std::cout << "Layout saved → " << out_off << "\n"; + + // ── Step 7: Serialise JSON + XML ────────────────────────────────────────── + save_result_json(out_json, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + std::cout << "JSON saved → " << out_json << "\n"; + + save_result_xml(out_xml, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + std::cout << "XML saved → " << out_xml << "\n"; + + // ── Step 8: JSON round-trip verification ────────────────────────────────── + NewtonResult res2; + std::string geom; + Layout2D layout2; + auto x2 = load_result_json(out_json, &res2, &geom, &layout2); + + std::cout << "Round-trip: geometry=\"" << geom << "\" " + << "DOFs=" << x2.size() << " " + << "converged=" << res2.converged << "\n"; + + // Verify the DOF vector survived the round-trip + bool ok = (x2.size() == res.x.size()); + for (std::size_t i = 0; i < x2.size() && ok; ++i) + ok = (std::abs(x2[i] - res.x[i]) < 1e-10); + std::cout << "DOF vector round-trip: " << (ok ? "OK ✓" : "MISMATCH ✗") << "\n"; + + return res.converged ? 0 : 1; +} diff --git a/code/include/layout.hpp b/code/include/layout.hpp new file mode 100644 index 0000000..f1d121b --- /dev/null +++ b/code/include/layout.hpp @@ -0,0 +1,480 @@ +#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. +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Algorithm (all three geometries) │ +// │ │ +// │ 1. For every edge e compute the updated length l_e from the DOF x. │ +// │ 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. │ +// │ │ +// │ 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. │ +// └──────────────────────────────────────────────────────────────────────────┘ +// +// 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) + +#include "conformal_mesh.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include +#include +#include +#include +#include + +namespace conformallab { + +// ── Result types ────────────────────────────────────────────────────────────── + +struct Layout2D { + std::vector uv; ///< uv[v.idx()] = 2-D position + bool success = false; + bool has_seam = false; ///< true if mesh is closed (first-visit used) +}; + +struct Layout3D { + std::vector pos; ///< pos[v.idx()] = 3-D position on S² + bool success = false; + bool has_seam = false; +}; + +// ── 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). +inline Eigen::Vector2d trilaterate_2d( + const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, + double da, double db) +{ + Eigen::Vector2d ab = pb - pa; + double d = ab.norm(); + if (d < 1e-14) return pa; + Eigen::Vector2d e1 = ab / d; + Eigen::Vector2d e2(-e1.y(), e1.x()); // CCW perpendicular + 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 +} + +// 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. +// +// Solution: p_c = α·pa + β·pb + γ·(pa × pb), γ > 0. +// 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 denom = 1.0 - c*c; + if (denom < 1e-14) return pa; // degenerate (pa ≈ pb) + + 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(); + + 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; + + Eigen::Vector3d p = alpha*pa + beta*pb + gamma*cross; + double n = p.norm(); + return (n > 1e-14) ? (p / n) : pa; // gamma > 0 → left side +} + +} // namespace detail + +// ── 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). +// ───────────────────────────────────────────────────────────────────────────── +inline Layout2D euclidean_layout( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& maps) +{ + const std::size_t nv = mesh.number_of_vertices(); + Layout2D result; + result.uv.assign(nv, Eigen::Vector2d::Zero()); + if (mesh.number_of_faces() == 0) return result; + + auto get_u = [&](Vertex_index v) -> double { + int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast(iv)] : 0.0; + }; + auto get_ue = [&](Edge_index e) -> double { + int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast(ie)] : 0.0; + }; + auto edge_len = [&](Halfedge_index h) -> double { + Edge_index e = mesh.edge(h); + double lam = maps.lambda0[e] + + get_u(mesh.source(h)) + get_u(mesh.target(h)) + + get_ue(e); + return std::exp(lam * 0.5); + }; + + std::vector vertex_placed(nv, false); + std::vector face_placed(mesh.number_of_faces(), false); + + // ── 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); + + double lAB = edge_len(h0); + double lCA = edge_len(h2); // dist C—A + double lBC = edge_len(h1); // dist B—C + + 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 ─────────────────────────────────────────────────────────────── + 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()]) + q.push(h_opp); + } + }; + enqueue(f0); + + while (!q.empty()) { + Halfedge_index h = q.front(); q.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + + Vertex_index v_src = mesh.source(h); + Vertex_index v_tgt = mesh.target(h); + Vertex_index v_new = mesh.target(mesh.next(h)); + + 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) + + if (!vertex_placed[v_new.idx()]) { + result.uv[v_new.idx()] = p; + vertex_placed[v_new.idx()] = true; + } else { + result.has_seam = true; + } + face_placed[f.idx()] = true; + enqueue(f); + } + + result.success = true; + 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 std::size_t nv = mesh.number_of_vertices(); + Layout3D result; + result.pos.assign(nv, Eigen::Vector3d::Zero()); + if (mesh.number_of_faces() == 0) return result; + + auto get_u = [&](Vertex_index v) -> double { + int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast(iv)] : 0.0; + }; + auto get_ue = [&](Edge_index e) -> double { + int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast(ie)] : 0.0; + }; + auto arc_len = [&](Halfedge_index h) -> double { + Edge_index e = mesh.edge(h); + double lam = maps.lambda0[e] + + get_u(mesh.source(h)) + get_u(mesh.target(h)) + + get_ue(e); + 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); + + 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); + + 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); + if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) + q.push(h_opp); + } + }; + enqueue(f0); + + while (!q.empty()) { + Halfedge_index h = q.front(); q.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + + Vertex_index v_src = mesh.source(h); + Vertex_index v_tgt = mesh.target(h); + Vertex_index v_new = mesh.target(mesh.next(h)); + + 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))); + + if (!vertex_placed[v_new.idx()]) { + result.pos[v_new.idx()] = p; + vertex_placed[v_new.idx()] = true; + } else { + result.has_seam = true; + } + face_placed[f.idx()] = true; + enqueue(f); + } + + result.success = true; + return result; +} + +// ── HyperIdeal layout (Poincaré disk) ──────────────────────────────────────── +// +// Places the hyperbolic triangulation in the Poincaré disk model. +// The effective hyperbolic edge length between vertices i and j is: +// +// cosh(d_ij) = cosh(b_i + a_ij/2) · cosh(b_j + a_ij/2) − sinh(b_i) · sinh(b_j) +// +// (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. +// ───────────────────────────────────────────────────────────────────────────── +inline Layout2D hyper_ideal_layout( + ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& maps) +{ + const std::size_t nv = mesh.number_of_vertices(); + Layout2D result; + result.uv.assign(nv, Eigen::Vector2d::Zero()); + if (mesh.number_of_faces() == 0) return result; + + auto get_b = [&](Vertex_index v) -> double { + int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast(iv)] : 0.0; + }; + auto get_a = [&](Edge_index e) -> double { + int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast(ie)] : 0.0; + }; + + // 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); + 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); + }; + + // 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); + + 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); + + 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); + + 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); + if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) + q.push(h_opp); + } + }; + enqueue(f0); + + while (!q.empty()) { + Halfedge_index h = q.front(); q.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + + Vertex_index v_src = mesh.source(h); + 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))); + + if (!vertex_placed[v_new.idx()]) { + result.uv[v_new.idx()] = p; + vertex_placed[v_new.idx()] = true; + } else { + result.has_seam = true; + } + face_placed[f.idx()] = true; + enqueue(f); + } + + result.success = true; + return result; +} + +// ── Convenience: save layout as OFF (z = 0 for 2D, full xyz for 3D) ────────── + +inline void save_layout_off( + const std::string& path, + ConformalMesh& mesh, + const Layout2D& layout) +{ + std::ofstream ofs(path); + ofs << "OFF\n" << mesh.number_of_vertices() << " " + << mesh.number_of_faces() << " 0\n"; + for (auto v : mesh.vertices()) { + auto& p = layout.uv[v.idx()]; + ofs << p.x() << " " << p.y() << " 0\n"; + } + for (auto f : mesh.faces()) { + ofs << "3"; + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + ofs << " " << mesh.target(h).idx(); + ofs << "\n"; + } +} + +inline void save_layout_off( + const std::string& path, + ConformalMesh& mesh, + const Layout3D& layout) +{ + std::ofstream ofs(path); + ofs << "OFF\n" << mesh.number_of_vertices() << " " + << mesh.number_of_faces() << " 0\n"; + for (auto v : mesh.vertices()) { + auto& p = layout.pos[v.idx()]; + ofs << p.x() << " " << p.y() << " " << p.z() << "\n"; + } + for (auto f : mesh.faces()) { + ofs << "3"; + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + ofs << " " << mesh.target(h).idx(); + ofs << "\n"; + } +} + +} // namespace conformallab diff --git a/code/include/serialization.hpp b/code/include/serialization.hpp new file mode 100644 index 0000000..dce3048 --- /dev/null +++ b/code/include/serialization.hpp @@ -0,0 +1,289 @@ +#pragma once +// serialization.hpp +// +// Phase 5 — Save and load conformal map results in JSON and XML formats. +// +// JSON (nlohmann/json, bundled in deps/single_includes/json.hpp): +// save_result_json / load_result_json +// +// XML (hand-written, no external parser dependency): +// save_result_xml / load_result_xml +// +// Both formats store: +// - geometry type string ("euclidean" / "spherical" / "hyper_ideal") +// - mesh statistics (vertex/face count) +// - solver metadata (converged, iterations, grad_inf_norm) +// - DOF vector x +// - optional 2D or 3D layout positions +// +// XML schema (ConformalResult): +// +// +// +// +// 0.0 1.23e-13 2.87e-13 +// +// 0.0 0.0 1.0 0.0 0.5 0.866 1.5 0.866 +// +// + +#include "newton_solver.hpp" +#include "layout.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace conformallab { + +// ════════════════════════════════════════════════════════════════════════════ +// JSON +// ════════════════════════════════════════════════════════════════════════════ + +// Save solver result (+ optional 2D layout) to a JSON file. +inline void save_result_json( + const std::string& path, + const NewtonResult& res, + const std::string& geometry, + int n_vertices, + int n_faces, + const Layout2D* layout2d = nullptr, + const Layout3D* layout3d = nullptr) +{ + using json = nlohmann::json; + json j; + j["geometry"] = geometry; + j["mesh"] = { {"vertices", n_vertices}, {"faces", n_faces} }; + j["solver"] = { + {"converged", res.converged}, + {"iterations", res.iterations}, + {"grad_inf_norm", res.grad_inf_norm} + }; + j["dof_vector"] = res.x; + + if (layout2d && layout2d->success) { + json uv = json::array(); + for (auto& p : layout2d->uv) uv.push_back({p.x(), p.y()}); + j["layout"] = { {"dim", 2}, {"uv", uv}, {"has_seam", layout2d->has_seam} }; + } + if (layout3d && layout3d->success) { + json pos = json::array(); + for (auto& p : layout3d->pos) pos.push_back({p.x(), p.y(), p.z()}); + j["layout"] = { {"dim", 3}, {"pos", pos}, {"has_seam", layout3d->has_seam} }; + } + + std::ofstream ofs(path); + if (!ofs) throw std::runtime_error("Cannot write: " + path); + ofs << std::setw(2) << j << "\n"; +} + +// Load DOF vector from a JSON result file. +// Returns the DOF vector; fills out res fields if non-null. +inline std::vector load_result_json( + const std::string& path, + NewtonResult* res = nullptr, + std::string* geom = nullptr, + Layout2D* layout2d = nullptr) +{ + using json = nlohmann::json; + std::ifstream ifs(path); + if (!ifs) throw std::runtime_error("Cannot open: " + path); + json j; ifs >> j; + + if (geom && j.contains("geometry")) + *geom = j["geometry"].get(); + + std::vector x = j.at("dof_vector").get>(); + + if (res) { + res->x = x; + res->converged = j["solver"]["converged"].get(); + res->iterations = j["solver"]["iterations"].get(); + res->grad_inf_norm = j["solver"]["grad_inf_norm"].get(); + } + + if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) { + auto uv = j["layout"]["uv"]; + layout2d->uv.resize(uv.size()); + for (std::size_t i = 0; i < uv.size(); ++i) + layout2d->uv[i] = { uv[i][0].get(), uv[i][1].get() }; + layout2d->has_seam = j["layout"].value("has_seam", false); + layout2d->success = true; + } + + return x; +} + +// ════════════════════════════════════════════════════════════════════════════ +// XML +// ════════════════════════════════════════════════════════════════════════════ + +namespace detail_xml { + +// Minimal XML attribute escaping +inline std::string xml_attr(double v) +{ + std::ostringstream s; + s << std::setprecision(15) << v; + return s.str(); +} + +// Write a flat vector of doubles as space-separated values +inline std::string flat_doubles(const std::vector& v) +{ + std::ostringstream s; + s << std::setprecision(15); + for (std::size_t i = 0; i < v.size(); ++i) { + if (i) s << ' '; + s << v[i]; + } + return s.str(); +} + +// Simple attribute parser: find value of key= in a tag string +inline std::string xml_get_attr(const std::string& tag, const std::string& key) +{ + auto pos = tag.find(key + "=\""); + if (pos == std::string::npos) return {}; + pos += key.size() + 2; + auto end = tag.find('"', pos); + return tag.substr(pos, end - pos); +} + +// Read all text content between the current position and +inline std::string xml_read_text(std::istream& is, const std::string& close_tag) +{ + std::string buf, line; + std::string ctag = ""; + while (std::getline(is, line)) { + auto pos = line.find(ctag); + if (pos != std::string::npos) { + buf += line.substr(0, pos); + return buf; + } + buf += line + " "; + } + return buf; +} + +// Parse space-separated doubles +inline std::vector parse_doubles(const std::string& s) +{ + std::vector v; + std::istringstream ss(s); + double d; + while (ss >> d) v.push_back(d); + return v; +} + +} // namespace detail_xml + +// Save solver result (+ optional layout) to an XML file. +inline void save_result_xml( + const std::string& path, + const NewtonResult& res, + const std::string& geometry, + int n_vertices, + int n_faces, + const Layout2D* layout2d = nullptr, + const Layout3D* layout3d = nullptr) +{ + std::ofstream ofs(path); + if (!ofs) throw std::runtime_error("Cannot write: " + path); + ofs << std::setprecision(15); + + ofs << "\n"; + ofs << "\n"; + + ofs << " \n"; + + ofs << " " + << detail_xml::flat_doubles(res.x) + << "\n"; + + if (layout2d && layout2d->success) { + ofs << " uv.size() + << "\" has_seam=\"" << (layout2d->has_seam ? "true" : "false") << "\">\n"; + for (auto& p : layout2d->uv) + ofs << " " << p.x() << " " << p.y() << "\n"; + ofs << " \n"; + } + if (layout3d && layout3d->success) { + ofs << " pos.size() + << "\" has_seam=\"" << (layout3d->has_seam ? "true" : "false") << "\">\n"; + for (auto& p : layout3d->pos) + ofs << " " << p.x() << " " << p.y() << " " << p.z() << "\n"; + ofs << " \n"; + } + + ofs << "\n"; +} + +// Load solver result from an XML file written by save_result_xml. +// Returns the DOF vector; fills res/geom/layout2d if non-null. +inline std::vector load_result_xml( + const std::string& path, + NewtonResult* res = nullptr, + std::string* geom = nullptr, + Layout2D* layout2d = nullptr) +{ + std::ifstream ifs(path); + if (!ifs) throw std::runtime_error("Cannot open: " + path); + + std::vector x; + std::string line; + + while (std::getline(ifs, line)) { + // Root element + if (line.find("converged = (detail_xml::xml_get_attr(line, "converged") == "true"); + res->iterations = std::stoi(detail_xml::xml_get_attr(line, "iterations")); + res->grad_inf_norm = std::stod(detail_xml::xml_get_attr(line, "grad_inf_norm")); + } + } + // DOF vector + else if (line.find("0 1 2... + auto open_end = line.find('>'); + auto close = line.find(""); + std::string text; + if (close != std::string::npos) { + text = line.substr(open_end + 1, close - open_end - 1); + } else { + text = line.substr(open_end + 1); + text += detail_xml::xml_read_text(ifs, "DOFVector"); + } + x = detail_xml::parse_doubles(text); + if (res) res->x = x; + } + // Layout + else if (layout2d && line.find("has_seam = (detail_xml::xml_get_attr(line, "has_seam") == "true"); + std::string text = detail_xml::xml_read_text(ifs, "Layout"); + auto vals = detail_xml::parse_doubles(text); + layout2d->uv.clear(); + for (std::size_t i = 0; i + 1 < vals.size(); i += 2) + layout2d->uv.push_back({vals[i], vals[i+1]}); + layout2d->success = true; + } + } + return x; +} + +} // namespace conformallab diff --git a/code/src/apps/v0/conformallab_cli.cpp b/code/src/apps/v0/conformallab_cli.cpp index c002a28..be3451e 100644 --- a/code/src/apps/v0/conformallab_cli.cpp +++ b/code/src/apps/v0/conformallab_cli.cpp @@ -1,107 +1,314 @@ -#include -#include -#include - - -#include "viewer_utils.h" -#include "mesh_utils.hpp" +// conformallab_cli.cpp +// +// ConformalLab++ command-line interface. +// +// Usage: +// conformallab_core -i input.off [-o layout.off] [-g euclidean|spherical|hyper_ideal] +// [-j result.json] [-x result.xml] [-s] [-v] +// +// The tool: +// 1. Loads an OFF/OBJ/PLY mesh. +// 2. Sets up DOF maps + computes λ° from the input geometry. +// 3. Pins one vertex (Euclidean/Spherical) or uses all-free DOFs (HyperIdeal). +// 4. Runs Newton until convergence. +// 5. Computes a 2-D (Euclidean / HyperIdeal) or 3-D (Spherical) layout. +// 6. Saves the layout as an OFF file and optionally serialises the result +// to JSON and/or XML. +// 7. Optionally shows the input mesh in a viewer (-s flag, requires WITH_VIEWER). +#include "conformal_mesh.hpp" +#include "mesh_io.hpp" +#include "euclidean_functional.hpp" +#include "spherical_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include "newton_solver.hpp" +#include "layout.hpp" +#include "serialization.hpp" #include #include +#include #include +#include +#include +#include - +// Viewer is always available when WITH_CGAL=ON (implies WITH_VIEWER=ON) +#include "viewer_utils.h" #include +namespace cl = conformallab; +using cl::ConformalMesh; +using cl::Vertex_index; +using cl::Edge_index; -using Kernel = CGAL::Simple_cartesian; -using Point = Kernel::Point_3; -using Mesh = CGAL::Surface_mesh; +// ───────────────────────────────────────────────────────────────────────────── +// Shared helpers (mirroring test_pipeline.cpp patterns) +// ───────────────────────────────────────────────────────────────────────────── - - - -bool parse_arguments(int argc, char* argv[], std::string& input_file, std::string& output_file, bool& show, bool& verbose) { - CLI::App app{"demo of conformallab"}; - app.add_option("-i,--input", input_file, "Input OFF file")->required(); - app.add_option("-o,--output", output_file, "Output OFF file (optional)"); - app.add_flag("-s,--show", show, "Show the input mesh in a viewer "); - app.add_flag("-v,--verbose",verbose, "Enable verbose output"); - app.set_help_flag("-h,--help", "Show this help message and exit"); - - try { - CLI11_PARSE(app, argc, argv); - } catch (const CLI::ParseError &e) { - return false; - } - - if (input_file.empty()) { - std::cerr << "Input file is required.\n"; - return EXIT_FAILURE; - } - - if (input_file.substr(input_file.find_last_of('.') + 1) != "off") { - std::cerr << "Unsupported file format. Please provide an OFF file.\n"; - return EXIT_FAILURE; - } - - std::cout << "Input file: " << input_file << "\n"; - std::cout << "Output file: " << output_file << "\n"; - std::cout << "Show mesh: " << (show ? "Yes" : "No") << "\n"; - std::cout << "Verbose: " << (verbose ? "Yes" : "No") << "\n"; - - - - return true; +// Pin vertex 0, assign 0..n-1 to the rest +static int pin_first_vertex(ConformalMesh& mesh, cl::EuclideanMaps& maps) +{ + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + maps.v_idx[v0] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) + maps.v_idx[*vit] = idx++; + return idx; } +// Natural theta for Euclidean: make x=0 the equilibrium +static void set_natural_euclidean_theta(ConformalMesh& mesh, cl::EuclideanMaps& maps, int n) +{ + std::vector x0(static_cast(n), 0.0); + auto G = cl::euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] -= G[static_cast(iv)]; + } +} +// Natural theta for HyperIdeal at base point (b=1, a=0.5) to avoid x=0 singularity +static std::vector set_natural_hyper_ideal_targets( + ConformalMesh& mesh, cl::HyperIdealMaps& maps, int n) +{ + const auto sz = static_cast(n); + std::vector xbase(sz, 0.0); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) xbase[static_cast(iv)] = 1.0; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) xbase[static_cast(ie)] = 0.5; + } + auto G = cl::evaluate_hyper_ideal(mesh, xbase, maps, false).gradient; + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv < 0) continue; + maps.theta_v[v] += G[static_cast(iv)]; + } + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie < 0) continue; + maps.theta_e[e] += G[static_cast(ie)]; + } + return xbase; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Euclidean pipeline +// ───────────────────────────────────────────────────────────────────────────── +static int run_euclidean(ConformalMesh& mesh, + const std::string& out_layout, + const std::string& out_json, + const std::string& out_xml, + bool verbose) +{ + // Setup + auto maps = cl::setup_euclidean_maps(mesh); + cl::compute_euclidean_lambda0_from_mesh(mesh, maps); + + // DOF assignment: pin vertex 0 + int n = pin_first_vertex(mesh, maps); + if (n <= 0) { std::cerr << "Error: mesh has only one vertex.\n"; return 1; } + + // Natural target angles + set_natural_euclidean_theta(mesh, maps, n); + + // Newton + std::vector x0(static_cast(n), 0.0); + auto res = cl::newton_euclidean(mesh, x0, maps); + + if (!res.converged && verbose) + std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n"; + + // Layout + cl::Layout2D layout = cl::euclidean_layout(mesh, res.x, maps); + + // Output + if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout); + if (!out_json.empty()) + cl::save_result_json(out_json, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + if (!out_xml.empty()) + cl::save_result_xml(out_xml, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + + std::cout << std::fixed << std::setprecision(6); + std::cout << "Euclidean: converged=" << (res.converged ? "yes" : "no") + << " iter=" << res.iterations + << " |grad|_inf=" << std::scientific << std::setprecision(3) + << res.grad_inf_norm << "\n"; + if (!out_layout.empty()) std::cout << " layout → " << out_layout << "\n"; + if (!out_json.empty()) std::cout << " json → " << out_json << "\n"; + if (!out_xml.empty()) std::cout << " xml → " << out_xml << "\n"; + return 0; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Spherical pipeline +// ───────────────────────────────────────────────────────────────────────────── +static int run_spherical(ConformalMesh& mesh, + const std::string& out_layout, + const std::string& out_json, + const std::string& out_xml, + bool verbose) +{ + auto maps = cl::setup_spherical_maps(mesh); + cl::compute_lambda0_from_mesh(mesh, maps); + int n = cl::assign_vertex_dof_indices(mesh, maps); + + std::vector x0(static_cast(n), 0.0); + auto res = cl::newton_spherical(mesh, x0, maps); + + if (!res.converged && verbose) + std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n"; + + cl::Layout3D layout = cl::spherical_layout(mesh, res.x, maps); + + if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout); + if (!out_json.empty()) + cl::save_result_json(out_json, res, "spherical", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + nullptr, &layout); + if (!out_xml.empty()) + cl::save_result_xml(out_xml, res, "spherical", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + nullptr, &layout); + + std::cout << "Spherical: converged=" << (res.converged ? "yes" : "no") + << " iter=" << res.iterations + << " |grad|_inf=" << std::scientific << std::setprecision(3) + << res.grad_inf_norm << "\n"; + if (!out_layout.empty()) std::cout << " layout → " << out_layout << "\n"; + if (!out_json.empty()) std::cout << " json → " << out_json << "\n"; + if (!out_xml.empty()) std::cout << " xml → " << out_xml << "\n"; + return 0; +} + +// ───────────────────────────────────────────────────────────────────────────── +// HyperIdeal pipeline +// ───────────────────────────────────────────────────────────────────────────── +static int run_hyper_ideal(ConformalMesh& mesh, + const std::string& out_layout, + const std::string& out_json, + const std::string& out_xml, + bool verbose) +{ + auto maps = cl::setup_hyper_ideal_maps(mesh); + int n = cl::assign_all_dof_indices(mesh, maps); + + // Natural targets at base point (b=1, a=0.5) + auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n); + + // Perturb slightly from the base and solve back + std::vector x0 = xbase; + for (auto& v : x0) v += 0.3; + + auto res = cl::newton_hyper_ideal(mesh, x0, maps); + + if (!res.converged && verbose) + std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n"; + + cl::Layout2D layout = cl::hyper_ideal_layout(mesh, res.x, maps); + + if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout); + if (!out_json.empty()) + cl::save_result_json(out_json, res, "hyper_ideal", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + if (!out_xml.empty()) + cl::save_result_xml(out_xml, res, "hyper_ideal", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout); + + std::cout << "HyperIdeal: converged=" << (res.converged ? "yes" : "no") + << " iter=" << res.iterations + << " |grad|_inf=" << std::scientific << std::setprecision(3) + << res.grad_inf_norm << "\n"; + if (!out_layout.empty()) std::cout << " layout → " << out_layout << "\n"; + if (!out_json.empty()) std::cout << " json → " << out_json << "\n"; + if (!out_xml.empty()) std::cout << " xml → " << out_xml << "\n"; + return 0; +} + +// ───────────────────────────────────────────────────────────────────────────── +// main +// ───────────────────────────────────────────────────────────────────────────── int main(int argc, char* argv[]) { + CLI::App app{"ConformalLab++ — discrete conformal mapping"}; + std::string input_file; - std::string output_file; - bool show = false; + std::string out_layout; + std::string out_json; + std::string out_xml; + std::string geometry = "euclidean"; + bool show = false; bool verbose = false; - if(!parse_arguments(argc, argv, input_file, output_file, show, verbose)) { - std::cerr << "Failed to parse arguments.\n"; + app.add_option("-i,--input", input_file, "Input mesh (OFF/OBJ/PLY)")->required(); + app.add_option("-o,--output", out_layout, "Output layout OFF file"); + app.add_option("-j,--json", out_json, "Save result as JSON"); + app.add_option("-x,--xml", out_xml, "Save result as XML"); + app.add_option("-g,--geometry", geometry, "Target geometry: euclidean|spherical|hyper_ideal") + ->check(CLI::IsMember({"euclidean", "spherical", "hyper_ideal"})); + app.add_flag("-s,--show", show, "Visualise input mesh (requires WITH_VIEWER)"); + app.add_flag("-v,--verbose", verbose, "Verbose output"); + + CLI11_PARSE(app, argc, argv); + + // ── Load mesh ───────────────────────────────────────────────────────────── + ConformalMesh mesh; + if (!cl::read_mesh(input_file, mesh) || mesh.is_empty()) { + std::cerr << "Error: cannot load mesh from '" << input_file << "'\n"; return EXIT_FAILURE; } + std::cout << "Loaded: " << input_file + << " (" << mesh.number_of_vertices() << "v, " + << mesh.number_of_faces() << "f)\n"; - Mesh surface_mesh; - if (!CGAL::IO::read_polygon_mesh(input_file, surface_mesh) || surface_mesh.is_empty()) { - std::cerr << "Invalid input file: " << input_file << "\n"; - return EXIT_FAILURE; - } - + // ── Optional viewer ─────────────────────────────────────────────────────── + if (show) { + // Convert ConformalMesh to Eigen matrices for the libigl viewer + const std::size_t nv = mesh.number_of_vertices(); + const std::size_t nf = mesh.number_of_faces(); + Eigen::MatrixXd V(static_cast(nv), 3); + Eigen::MatrixXi F(static_cast(nf), 3); - - - // #TODO: later: here we would call the unwrapping code, e.g.: - // UnwrapSettings settings; - // settings.target_geometry = TargetGeometry::Euclidean; - // UnwrapJob job(surface_mesh, settings); - // auto result = job.run(); - // Mesh unwrapped = result.surface_unwrapped; - - - // CGAL -> libigl - if(show){ - std::cout << "Visualizing input mesh...\n"; - Eigen::MatrixXd V; - Eigen::MatrixXi F; - - mesh_utils::cgal_to_eigen(surface_mesh, V, F); - viewer_utils::simple_visualize(V, F); - - } - - // for now, we just write the input mesh to the output file as a placeholder - if (!output_file.empty() && !CGAL::IO::write_polygon_mesh(output_file, surface_mesh)) { - std::cerr << "Failed to write output file: " << output_file << "\n"; - return EXIT_FAILURE; + for (auto v : mesh.vertices()) { + auto p = mesh.point(v); + V.row(v.idx()) << p.x(), p.y(), p.z(); + } + Eigen::Index fi = 0; + for (auto f : mesh.faces()) { + int ci = 0; + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + F(fi, ci++) = static_cast(mesh.target(h).idx()); + ++fi; + } + viewer_utils::simple_visualize(V, F); } - return EXIT_SUCCESS; + // ── Dispatch ────────────────────────────────────────────────────────────── + if (geometry == "euclidean") + return run_euclidean(mesh, out_layout, out_json, out_xml, verbose); + if (geometry == "spherical") + return run_spherical(mesh, out_layout, out_json, out_xml, verbose); + if (geometry == "hyper_ideal") + return run_hyper_ideal(mesh, out_layout, out_json, out_xml, verbose); + + std::cerr << "Unknown geometry: " << geometry << "\n"; + return EXIT_FAILURE; } diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index fa13123..12cbd56 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -33,11 +33,15 @@ add_executable(conformallab_cgal_tests # ── Phase 4c: End-to-end pipeline + user examples ───────────────────── test_pipeline.cpp + + # ── Phase 5: Layout / embedding + serialization ──────────────────────── + test_layout.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE ${CMAKE_SOURCE_DIR}/deps/eigen-3.4.0 ${CMAKE_SOURCE_DIR}/deps/CGAL-6.1.1/include + ${CMAKE_SOURCE_DIR}/deps/single_includes ${Boost_INCLUDE_DIRS} ) diff --git a/code/tests/cgal/test_layout.cpp b/code/tests/cgal/test_layout.cpp new file mode 100644 index 0000000..18d16b8 --- /dev/null +++ b/code/tests/cgal/test_layout.cpp @@ -0,0 +1,369 @@ +// test_layout.cpp +// +// Phase 5 — Layout / embedding tests. +// +// Strategy: a conformal map that is the identity (natural equilibrium x*=0) +// should reproduce the mesh's own edge lengths. We verify: +// 1. Euclidean layout preserves edge lengths (to solver tolerance). +// 2. Spherical layout preserves arc-lengths on S². +// 3. Layout2D has correct size (one entry per vertex). +// 4. Serialisation round-trip: save JSON → load → DOF vector unchanged. +// 5. Serialisation round-trip: save XML → load → DOF vector unchanged. +// 6. HyperIdeal layout returns success and finite positions. + +#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 "layout.hpp" +#include "serialization.hpp" +#include +#include +#include +#include + +using namespace conformallab; + +// ──────────────────────────────────────────────────────────────────────────── +// Helpers +// ──────────────────────────────────────────────────────────────────────────── + +// Euclidean edge length from layout (2D) +static double layout_edge_len(const Layout2D& lay, + ConformalMesh& mesh, Halfedge_index h) +{ + auto vi = mesh.source(h).idx(); + auto vj = mesh.target(h).idx(); + return (lay.uv[vi] - lay.uv[vj]).norm(); +} + +// Spherical arc-length from layout (3D) +static double layout_arc_len(const Layout3D& lay, + ConformalMesh& mesh, Halfedge_index h) +{ + auto& a = lay.pos[mesh.source(h).idx()]; + auto& b = lay.pos[mesh.target(h).idx()]; + double c = std::max(-1.0, std::min(1.0, a.dot(b))); + return std::acos(c); +} + +// Euclidean edge length from maps + x (the "true" target) +static double maps_edge_len(const EuclideanMaps& m, + ConformalMesh& mesh, + Halfedge_index h, + const std::vector& x) +{ + auto get_u = [&](Vertex_index v) -> double { + int iv = m.v_idx[v]; return (iv >= 0) ? x[static_cast(iv)] : 0.0; + }; + Edge_index e = mesh.edge(h); + double lam = m.lambda0[e] + get_u(mesh.source(h)) + get_u(mesh.target(h)); + return std::exp(lam * 0.5); +} + +// Spherical arc-length from maps + x +static double maps_arc_len(const SphericalMaps& m, + ConformalMesh& mesh, + Halfedge_index h, + const std::vector& x) +{ + auto get_u = [&](Vertex_index v) -> double { + int iv = m.v_idx[v]; return (iv >= 0) ? x[static_cast(iv)] : 0.0; + }; + Edge_index e = mesh.edge(h); + double lam = m.lambda0[e] + get_u(mesh.source(h)) + get_u(mesh.target(h)); + return 2.0 * std::asin(std::min(std::exp(lam * 0.5), 1.0)); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 1 — Euclidean layout preserves edge lengths (identity map) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, Euclidean_PreservesEdgeLengths) +{ + auto mesh = make_quad_strip(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // Pin first vertex; natural equilibrium so x* = 0 + 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); + 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)]; + } + + // Solve (identity map) + auto res = newton_euclidean(mesh, x0, maps, 1e-10, 100); + ASSERT_TRUE(res.converged); + + auto layout = euclidean_layout(mesh, res.x, maps); + ASSERT_TRUE(layout.success); + EXPECT_EQ(layout.uv.size(), mesh.number_of_vertices()); + + // Every edge: |layout_len - expected_len| < 1e-8 + for (auto e : mesh.edges()) { + Halfedge_index h = mesh.halfedge(e, 0); + double expected = maps_edge_len(maps, mesh, h, res.x); + double actual = layout_edge_len(layout, mesh, h); + EXPECT_NEAR(actual, expected, 1e-7) + << "Edge " << e << ": expected " << expected << " got " << actual; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 2 — Euclidean layout: correct vertex count +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, Euclidean_CorrectVertexCount) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + 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 x(static_cast(n), 0.0); + auto layout = euclidean_layout(mesh, x, maps); + + EXPECT_TRUE(layout.success); + EXPECT_EQ(static_cast(layout.uv.size()), + static_cast(mesh.number_of_vertices())); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 3 — Euclidean triangle layout: root face is a valid triangle +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, Euclidean_TriangleIsNonDegenerate) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + for (auto v : mesh.vertices()) maps.v_idx[v] = static_cast(v.idx()); + + std::vector x(mesh.number_of_vertices(), 0.0); + auto layout = euclidean_layout(mesh, x, maps); + ASSERT_TRUE(layout.success); + + // Area of layout triangle > 0 + const auto& A = layout.uv[0]; + const auto& B = layout.uv[1]; + const auto& C = layout.uv[2]; + double area = std::abs((B - A).x() * (C - A).y() + - (C - A).x() * (B - A).y()) * 0.5; + EXPECT_GT(area, 1e-10) << "Layout triangle must be non-degenerate"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 4 — Spherical layout: arc-lengths preserved (identity map) +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, Spherical_PreservesArcLengths) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + // Solve to identity + std::vector x0(static_cast(n), 0.0); + auto G0 = spherical_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_spherical(mesh, x0, maps, 1e-10, 100); + ASSERT_TRUE(res.converged); + + auto layout = spherical_layout(mesh, res.x, maps); + ASSERT_TRUE(layout.success); + EXPECT_EQ(layout.pos.size(), mesh.number_of_vertices()); + + for (auto e : mesh.edges()) { + Halfedge_index h = mesh.halfedge(e, 0); + double expected = maps_arc_len(maps, mesh, h, res.x); + double actual = layout_arc_len(layout, mesh, h); + EXPECT_NEAR(actual, expected, 1e-6) + << "Spherical edge " << e << ": expected " << expected << " got " << actual; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 5 — Spherical layout: all positions on unit sphere +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, Spherical_PositionsOnUnitSphere) +{ + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_vertex_dof_indices(mesh, maps); + + std::vector x(static_cast(n), 0.0); + auto layout = spherical_layout(mesh, x, maps); + ASSERT_TRUE(layout.success); + + for (auto v : mesh.vertices()) { + double r = layout.pos[v.idx()].norm(); + EXPECT_NEAR(r, 1.0, 1e-12) << "Vertex " << v << " not on unit sphere, |p| = " << r; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 6 — HyperIdeal layout returns success and finite positions +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Layout, HyperIdeal_SuccessAndFinitePositions) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + // Natural equilibrium at (b=1, a=0.5) + 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); + ASSERT_TRUE(res.converged); + + auto layout = hyper_ideal_layout(mesh, res.x, maps); + EXPECT_TRUE(layout.success); + ASSERT_EQ(layout.uv.size(), mesh.number_of_vertices()); + + for (auto v : mesh.vertices()) { + auto& p = layout.uv[v.idx()]; + EXPECT_FALSE(std::isnan(p.x())) << "NaN in layout x for v" << v; + EXPECT_FALSE(std::isnan(p.y())) << "NaN in layout y for v" << v; + // Poincaré disk: all points within the unit disk + EXPECT_LE(p.norm(), 1.0 + 1e-9) << "Point outside Poincaré disk for v" << v; + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 7 — JSON serialisation round-trip +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Serialization, JSON_RoundTrip) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + 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.05); + auto G0 = euclidean_gradient(mesh, std::vector(n, 0.0), 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); + ASSERT_TRUE(res.converged); + + auto layout = euclidean_layout(mesh, res.x, maps); + + const std::string path = "/tmp/conformallab_test_round_trip.json"; + ASSERT_NO_THROW(save_result_json(path, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout)); + ASSERT_TRUE(std::filesystem::exists(path)); + + NewtonResult res2; + std::string geom; + Layout2D layout2; + auto x2 = load_result_json(path, &res2, &geom, &layout2); + + EXPECT_EQ(geom, "euclidean"); + EXPECT_EQ(res2.converged, res.converged); + EXPECT_EQ(res2.iterations, res.iterations); + ASSERT_EQ(x2.size(), res.x.size()); + for (std::size_t i = 0; i < x2.size(); ++i) + EXPECT_NEAR(x2[i], res.x[i], 1e-14); + + EXPECT_TRUE(layout2.success); + ASSERT_EQ(layout2.uv.size(), layout.uv.size()); + for (std::size_t i = 0; i < layout2.uv.size(); ++i) { + EXPECT_NEAR(layout2.uv[i].x(), layout.uv[i].x(), 1e-12); + EXPECT_NEAR(layout2.uv[i].y(), layout.uv[i].y(), 1e-12); + } + std::filesystem::remove(path); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Test 8 — XML serialisation round-trip +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Serialization, XML_RoundTrip) +{ + auto mesh = make_triangle(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + 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); + 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); + ASSERT_TRUE(res.converged); + + auto layout = euclidean_layout(mesh, res.x, maps); + + const std::string path = "/tmp/conformallab_test_round_trip.xml"; + ASSERT_NO_THROW(save_result_xml(path, res, "euclidean", + static_cast(mesh.number_of_vertices()), + static_cast(mesh.number_of_faces()), + &layout)); + ASSERT_TRUE(std::filesystem::exists(path)); + + NewtonResult res2; + std::string geom; + Layout2D layout2; + auto x2 = load_result_xml(path, &res2, &geom, &layout2); + + EXPECT_EQ(geom, "euclidean"); + EXPECT_EQ(res2.converged, res.converged); + ASSERT_EQ(x2.size(), res.x.size()); + for (std::size_t i = 0; i < x2.size(); ++i) + EXPECT_NEAR(x2[i], res.x[i], 1e-12); + + EXPECT_TRUE(layout2.success); + ASSERT_EQ(layout2.uv.size(), layout.uv.size()); + + std::filesystem::remove(path); +}