ConformalLabpp/doc/architecture/design-decisions.md

Key Design Decisions

Rationale for the architectural choices that distinguish conformallab++ from the Java original and from generic geometry-processing frameworks.

---

CGAL Surface_mesh as the halfedge data structure

The Java library uses CoHDS — a custom intrusive halfedge data structure with CoVertex, CoEdge, CoFace types that carry domain-specific data directly as fields.

conformallab++ replaces this with CGAL::Surface_mesh<Point3> and attaches data via named property maps:

// Java: vertex.getLambda()           → C++: maps.lambda0[v]
// Java: edge.getAlpha()              → C++: maps.e_alpha[e]
// Java: vertex.getSolverIndex()      → C++: maps.v_idx[v]

auto [lambda0, ok] = mesh.add_property_map<Edge_index, double>("e:lambda0", 0.0);
lambda0[e] = 1.234;

This decouples the mesh topology from the algorithm data, makes it straightforward to attach multiple independent data sets to the same mesh, and will enable the Phase 8 traits-class design to work with any CGAL-conforming mesh type.

---

DOF vector convention

All three functionals use the same indexing scheme: a flat std::vector<double> x indexed by v_idx[v] (vertices) and e_idx[e] (edges, HyperIdeal only). Index value -1 means "pinned" — the DOF is fixed at zero and excluded from the Newton system.

x[maps.v_idx[v]]  = uᵥ   (conformal scale factor, Euclidean/Spherical)
x[maps.e_idx[e]]  = λₑ   (edge log-length variable, HyperIdeal only)
-1                        pinned — u_v = 0 / λ_e = 0

This is consistent across all three geometry modes, enabling the same Newton solver and linear system infrastructure to serve all three without branching.

---

Scalar type: double, with one localized exception

The kernel is CGAL::Simple_cartesian<double> and the whole numerical core (functionals, Hessians, Newton, Eigen sparse solvers) is double. Conformal flattening minimises a smooth energy over transcendental (floating-point) lengths and angles, so exact/extended arithmetic buys nothing there — the Java original is double here too. The CGAL wrapper (include/CGAL/*.h) is templated on FT for API genericity, but delegates to the double core.

The one exception (deferred, Phase 9c/10): hyperbolic uniformization (FundamentalPolygonUtility, CanonicalFormUtility) composes products of isometry generators whose entries grow exponentially; double then fails to verify the group relation ∏gᵢ = Id. The Java original handles this with RnBig/PnBig/P2Big at MathContext(50). When ported, this needs a localized high-precision substrate (boost::multiprecision::cpp_dec_float_50 or MPFR mpreal) inside the uniformization module only — never the core or the Eigen solver. See CLAUDE.md Phase-9 note and doc/roadmap/java-parity.md (*Big exception).

---

Priority-BFS layout

A naive BFS layout places faces in arbitrary order; trilateration errors accumulate along the BFS frontier. conformallab++ uses a priority min-heap on BFS depth:

depth(face) = max(depth[v_src], depth[v_tgt]) + 1   for each new face

Faces with smaller depth (closer to the root) are placed first. This means each face's trilateration uses the two most accurately-placed adjacent vertices, minimising error propagation across the mesh.

Root face selection: largest 3-D area face, with an additional 1.5× bonus for interior faces over boundary faces. This heuristic places the root where metric distortion is lowest.

---

halfedge_uv semantics

layout.uv[v.idx()] gives the primary UV coordinate of vertex v — the position from the shallowest BFS visit. At seam edges this is insufficient for GPU rendering: two faces sharing a seam vertex need different UV values for that vertex.

layout.halfedge_uv[h.idx()] stores the UV of source(h) as seen from face(h):

halfedge h  →  face(h)  →  source(h) has UV = halfedge_uv[h.idx()]
opposite(h) →  face(h') →  source(h) has UV = halfedge_uv[opposite(h).idx()]
                            (different value at a seam)

At seam halfedges the two opposite halfedges carry different UV values — each face gets its own copy of the seam vertex. This enables a proper GPU texture atlas without vertex duplication in the index buffer.

---

Spherical Hessian sign convention

The spherical energy functional is concave (negative semidefinite Hessian). Standard Newton would require solving H·Δx = −G with NSD H, which Cholesky cannot handle.

newton_spherical() solves (−H)·Δx = G instead — algebraically identical, but −H is PSD and SimplicialLDLT works correctly. This sign flip is handled transparently inside newton_spherical(); callers need not be aware of it.

The gradient sign in spherical mode is also flipped vs. Euclidean:

• Euclidean: G_v = actual_sum − Θᵥ
• Spherical: G_v = Θᵥ − actual_sum

Both conventions drive the same equilibrium condition G = 0.

---

HyperIdeal Hessian via finite differences

The analytic HyperIdeal Hessian requires differentiating through the chain (bᵢ, aₑ) → lᵢⱼ → ζ₁₃/ζ₁₄/ζ₁₅ → αᵢⱼ/βᵢ with four vertex-type combinations per edge — substantial implementation complexity.

conformallab++ uses a symmetric finite-difference Hessian instead:

H[i,j] = (G(x + ε·eⱼ)[i] − G(x − ε·eⱼ)[i]) / (2ε),   ε = 1e-5

Properties:

• O(ε²) accuracy — relative error ≈ 10⁻¹⁰ at ε = 10⁻⁵
• PSD guaranteed by strict convexity of the HyperIdeal energy (Springborn 2020)
• Symmetrised automatically: H = (H + Hᵀ) / 2
• Cost: n extra gradient evaluations per Newton step (acceptable for < 500 DOFs)

The analytic Hessian is deferred to Phase 9b. See roadmap/java-parity.md.