diff --git a/code/include/fundamental_domain.hpp b/code/include/fundamental_domain.hpp new file mode 100644 index 0000000..0ea180a --- /dev/null +++ b/code/include/fundamental_domain.hpp @@ -0,0 +1,204 @@ +#pragma once +// fundamental_domain.hpp +// +// Phase 7 — Fundamental domain polygon for closed surfaces. +// +// For a closed genus-g surface cut open via a CutGraph + Euclidean layout: +// +// The universal cover is tiled by copies of the cut-open disk. +// The fundamental domain is the polygon whose sides are identified in pairs +// by the holonomy generators. +// +// ─── Genus-1 (flat torus) ──────────────────────────────────────────────────── +// +// Parallelogram with vertices 0, ω_1, ω_1 + ω_2, ω_2. +// The four edges are identified in pairs: +// bottom (0 → ω_1) ≡ top (ω_2 → ω_1 + ω_2) — translation ω_2 +// left (0 → ω_2) ≡ right (ω_1 → ω_1 + ω_2) — translation ω_1 +// +// ─── Genus g > 1 (general) ────────────────────────────────────────────────── +// +// The standard 4g-polygon with sides labelled a_1 b_1 a_1^{-1} b_1^{-1} ... +// can be recovered from the layout boundary, but requires walking the +// boundary of the cut-open mesh — not yet implemented (see note below). +// +// For now, this file provides the genus-1 parallelogram only. +// The polygon vertices for genus-1 are computed from the holonomy generators. +// +// ─── API ───────────────────────────────────────────────────────────────────── +// +// FundamentalDomain fd = compute_fundamental_domain_genus1(holonomy); +// fd.vertices — 2D polygon corners (size = 4 for genus-1) +// fd.edge_identifications — pairs (i, j) meaning edge i is identified with j +// fd.is_valid() — true if genus == 1 and data makes sense + +#include "layout.hpp" +#include "period_matrix.hpp" +#include +#include + +namespace conformallab { + +// ───────────────────────────────────────────────────────────────────────────── +// FundamentalDomain +// ───────────────────────────────────────────────────────────────────────────── + +struct FundamentalDomain { + /// Polygon corners in order (CCW). Size = 4 for genus-1. + std::vector vertices; + + /// edge_identifications[k] = (i, j) means the edge from vertices[i] to + /// vertices[(i+1) % n] is identified with the edge from vertices[j] to + /// vertices[(j+1) % n] (with matching orientation). + std::vector> edge_identifications; + + /// Holonomy generators (one per identified edge pair). + /// For genus-1: generators[0] = ω_1, generators[1] = ω_2. + std::vector generators; + + bool is_valid() const { return vertices.size() >= 3; } +}; + +// ───────────────────────────────────────────────────────────────────────────── +// compute_fundamental_domain_genus1 +// +// Builds the parallelogram fundamental domain from Euclidean holonomy data +// with exactly 2 generators ω_1, ω_2. +// +// Vertices (CCW): +// v0 = (0, 0) +// v1 = ω_1 +// v2 = ω_1 + ω_2 +// v3 = ω_2 +// +// Edge identifications: +// bottom (v0→v1) ≡ top (v3→v2) by ω_2 +// left (v3→v0) ≡ right (v2→v1) by ω_1 (reversed convention) +// ───────────────────────────────────────────────────────────────────────────── +inline FundamentalDomain compute_fundamental_domain_genus1( + const HolonomyData& hol) +{ + FundamentalDomain fd; + if (hol.translations.size() < 2) return fd; + + Eigen::Vector2d w1 = hol.translations[0]; + Eigen::Vector2d w2 = hol.translations[1]; + + // Ensure CCW orientation: cross product z-component w1 × w2 > 0 + double cross = w1.x() * w2.y() - w1.y() * w2.x(); + if (cross < 0.0) std::swap(w1, w2); + + Eigen::Vector2d origin = Eigen::Vector2d::Zero(); + fd.vertices = { origin, w1, w1 + w2, w2 }; + + // Edge 0: v0→v1 (= bottom), Edge 2: v3→v2 (= top, reversed) + // Identification: bottom ≡ top translated by w2 + // Edge 1: v1→v2 (= right), Edge 3: v0→v3... wait let me use standard labeling: + // Edges by index: 0: v0→v1, 1: v1→v2, 2: v2→v3, 3: v3→v0 + // Identifications: 0 ≡ 2 (reversed: bottom ≡ top by w2) + // 1 ≡ 3 (reversed: right ≡ left by w1) + fd.edge_identifications = { {0, 2}, {1, 3} }; + fd.generators = { w1, w2 }; + return fd; +} + +// ───────────────────────────────────────────────────────────────────────────── +// compute_fundamental_domain +// +// Dispatcher: for genus-1 uses compute_fundamental_domain_genus1. +// For higher genus returns an empty FundamentalDomain (not yet implemented). +// ───────────────────────────────────────────────────────────────────────────── +// +// TODO(Phase 8): Implement the standard 4g-gon fundamental domain for genus g > 1. +// +// Algorithm outline (boundary-walk method): +// ───────────────────────────────────────── +// 1. Construct the CutGraph on the cut-open mesh (already done upstream). +// This yields 2g cut edges; cutting them converts the closed surface into +// a topological disk. +// +// 2. Walk the boundary of the cut-open disk in CCW order: +// Start from any boundary halfedge and follow `next(h)` along the boundary +// (i.e. skip to the next boundary halfedge at each vertex). Collect the +// 2·(4g) = 8g boundary halfedges in order. +// Each halfedge h_k corresponds to a UV vertex `halfedge_uv[h_k.idx()]`. +// +// 3. Identify paired sides: +// The 4g sides of the polygon alternate as a_1 b_1 a_1^{-1} b_1^{-1} … +// For each cut edge e_i (i = 1 … 2g) the two sides that are identified +// are those whose source/target vertices match under the holonomy generator +// ω_i (Euclidean) or T_i (hyperbolic). +// Record the identifications as edge_identifications[k] = (i, j). +// +// 4. Fill FundamentalDomain: +// vertices = UV corners from the boundary walk. +// edge_identifications = paired-edge list from step 3. +// generators = holonomy.translations (Euclidean) or the +// fixed points of holonomy.mobius_maps (hyperbolic, +// requires computing axis of T_i ∈ SU(1,1)). +// +// References: +// Erickson & Whittlesey, "Greedy optimal homotopy and homology generators" +// SODA 2005. +// Desbrun, Kanso, Tong, "Discrete Differential Forms for Computational +// Modeling", in Discrete Differential Geometry (2008). +// +// Note: The Siegel period matrix Ω ∈ H_g (g×g complex symmetric, Im Ω > 0) +// for genus g > 1 also requires integration of holomorphic differentials — +// this is intentionally deferred and NOT implemented here. +// See period_matrix.hpp for the genus-1 case (τ = ω_2/ω_1 ∈ ℍ). +// ───────────────────────────────────────────────────────────────────────────── +inline FundamentalDomain compute_fundamental_domain( + const HolonomyData& hol) +{ + int n = static_cast(hol.translations.size()); + int g = n / 2; + if (g == 1) return compute_fundamental_domain_genus1(hol); + // Higher genus: boundary-walk 4g-polygon — not yet implemented (see TODO above). + return FundamentalDomain{}; +} + +// ───────────────────────────────────────────────────────────────────────────── +// tiling_copy +// +// Given a Layout2D for the cut-open surface and two lattice generators ω_1, ω_2, +// return a translated copy of the layout shifted by m·ω_1 + n·ω_2. +// Useful for visualising the tiled universal cover. +// ───────────────────────────────────────────────────────────────────────────── +inline Layout2D tiling_copy(const Layout2D& layout, + const Eigen::Vector2d& w1, + const Eigen::Vector2d& w2, + int m, int n) +{ + Layout2D copy = layout; + Eigen::Vector2d shift = static_cast(m) * w1 + + static_cast(n) * w2; + for (auto& p : copy.uv) p += shift; + return copy; +} + +// ───────────────────────────────────────────────────────────────────────────── +// tiling_neighbourhood +// +// Returns a vector of tiling copies for (m, n) with |m| ≤ m_max, |n| ≤ n_max. +// The result includes the original (m=0, n=0) at index (m_max)(2*n_max+1)+n_max. +// ───────────────────────────────────────────────────────────────────────────── +inline std::vector tiling_neighbourhood( + const Layout2D& layout, + const HolonomyData& hol, + int m_max = 2, int n_max = 2) +{ + std::vector tiles; + if (hol.translations.size() < 2) { + tiles.push_back(layout); + return tiles; + } + const Eigen::Vector2d& w1 = hol.translations[0]; + const Eigen::Vector2d& w2 = hol.translations[1]; + for (int m = -m_max; m <= m_max; ++m) + for (int n = -n_max; n <= n_max; ++n) + tiles.push_back(tiling_copy(layout, w1, w2, m, n)); + return tiles; +} + +} // namespace conformallab diff --git a/code/include/layout.hpp b/code/include/layout.hpp index 31db33e..8b5f1ed 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -1,38 +1,55 @@ #pragma once // layout.hpp // -// Phase 5/6 — Layout / embedding: DOF vector → vertex coordinates in the +// Phase 5/6/7 — Layout / embedding: DOF vector → vertex coordinates in the // target geometry via BFS-trilateration. // // ┌──────────────────────────────────────────────────────────────────────────┐ // │ 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. │ -// │ 4. (Phase 6) If a CutGraph is supplied, cut edges are treated as │ -// │ boundary; crossing them records a HolonomyData entry instead. │ +// │ 2. Choose root face = largest 3-D area face (quality heuristic). │ +// │ 3. Priority BFS over the dual graph: faces are processed in order of │ +// │ shortest distance (BFS depth) from the root face, minimising error │ +// │ accumulation. Equal depth: arbitrary tie-break. │ +// │ 4. For each face: trilaterate the one unplaced vertex from the two │ +// │ already-placed vertices on the shared edge. │ +// │ 5. If a CutGraph is supplied, cut edges are treated as boundary; │ +// │ holonomy is recorded for each cut edge. │ +// │ 6. After the first connected component, unvisited faces start new │ +// │ BFS sweeps (multi-component support). │ // │ │ -// │ For open meshes the placement is globally consistent. │ -// │ For closed meshes without a cut: first visit wins, has_seam = true. │ -// │ For closed meshes with a CutGraph: each vertex is placed exactly once; │ -// │ the holonomy (translation / Möbius) of each cut is recorded. │ +// │ For open meshes: globally consistent placement. │ +// │ For closed meshes without CutGraph: first (shallowest) visit wins, │ +// │ has_seam = true. │ +// │ For closed meshes with CutGraph: each vertex placed exactly once; │ +// │ holonomy = translation (Euclidean) or Möbius map (hyperbolic). │ // └──────────────────────────────────────────────────────────────────────────┘ // -// Euclidean trilateration — exact (analytic formula in ℝ²) -// Spherical trilateration — exact (spherical law of cosines on S²) -// Hyperbolic trilateration — exact (hyperbolic law of cosines + Möbius maps -// in the Poincaré disk model) +// Trilateration +// Euclidean — exact (analytic formula in ℝ²) +// Spherical — exact (spherical law of cosines on S²) +// Hyperbolic — exact (hyperbolic law of cosines + Möbius maps, +// Poincaré disk model) // // Normalisation -// normalise_euclidean(layout) — centroid → origin, major axis → x-axis -// normalise_hyperbolic(layout) — Möbius map centering to origin -// normalise_spherical(layout) — rotate centroid to north pole +// normalise_euclidean(layout) — centroid → origin, major axis → x-axis (PCA) +// normalise_hyperbolic(layout) — face-area-weighted iterative Möbius centering +// normalise_spherical(layout) — rotate centroid to north pole (Rodrigues) // -// Holonomy (Euclidean, genus-g surfaces with CutGraph) -// HolonomyData.translations[i] — translation vector ω_i for cut edge i -// (for a flat torus: ω_1, ω_2 are the lattice generators) +// Holonomy (genus-g surfaces with CutGraph) +// HolonomyData.translations[i] — translation ω_i (Euclidean / spherical) +// HolonomyData.mobius_maps[i] — Möbius map T_i (hyperbolic) +// For a flat torus: ω_1, ω_2 are the lattice generators; τ = ω_2/ω_1 ∈ ℍ. +// +// Texture atlas +// Layout2D.halfedge_uv[h.idx()] — UV of source(h) as seen from face(h). +// At seam edges the two opposite halfedges carry different UV values, +// enabling proper per-halfedge UV for GPU texture atlasing. +// +// Möbius map +// MobiusMap::from_three(z1→w1, z2→w2, z3→w3) — fit a Möbius transformation +// to three point correspondences (via 3×3 complex linear system). #include "conformal_mesh.hpp" #include "euclidean_functional.hpp" @@ -45,30 +62,102 @@ #include #include #include +#include #include #include namespace conformallab { +// ── Möbius map ──────────────────────────────────────────────────────────────── +// +// T(z) = (a·z + b) / (c·z + d) +// +// For hyperbolic holonomy the map is an orientation-preserving isometry of +// the Poincaré disk (SU(1,1) element). +// ───────────────────────────────────────────────────────────────────────────── +struct MobiusMap { + using C = std::complex; + C a{1.0, 0.0}; + C b{0.0, 0.0}; + C c{0.0, 0.0}; + C d{1.0, 0.0}; + + C apply(C z) const { return (a * z + b) / (c * z + d); } + + Eigen::Vector2d apply(const Eigen::Vector2d& p) const { + C w = apply(C(p.x(), p.y())); + return Eigen::Vector2d(w.real(), w.imag()); + } + + static MobiusMap identity() { return {C(1), C(0), C(0), C(1)}; } + + MobiusMap inverse() const { return {d, -b, -c, a}; } + MobiusMap compose(const MobiusMap& T) const { + return { a*T.a + b*T.c, a*T.b + b*T.d, + c*T.a + d*T.c, c*T.b + d*T.d }; + } + + bool is_identity(double tol = 1e-9) const { + if (std::abs(d) < 1e-14) return false; + C a_ = a/d, b_ = b/d, c_ = c/d; + return std::abs(a_ - C(1)) < tol && std::abs(b_) < tol && std::abs(c_) < tol; + } + + /// Fit T to three point correspondences T(z_i) = w_i. + /// Sets d=1 and solves the 3×3 complex linear system. + /// Returns identity when the system is (near-)singular. + static MobiusMap from_three(C z1, C w1, C z2, C w2, C z3, C w3) { + Eigen::Matrix3cd M; + M << z1, C(1), -w1*z1, + z2, C(1), -w2*z2, + z3, C(1), -w3*z3; + Eigen::Vector3cd rhs; rhs << w1, w2, w3; + Eigen::ColPivHouseholderQR qr(M); + if (qr.rank() < 3) return identity(); + Eigen::Vector3cd sol = qr.solve(rhs); + return { sol[0], sol[1], sol[2], C(1.0) }; + } +}; + // ── 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 and no CutGraph given + /// uv[v.idx()] — primary 2-D position (first / shallowest-BFS-depth visit). + std::vector uv; + + /// halfedge_uv[h.idx()] — UV of source(h) as seen from face(h). + /// + /// For interior (non-seam) halfedges: equals uv[source(h).idx()]. + /// For seam halfedges: carries the UV from the virtual unfolding across + /// the cut (i.e. the trilaterated position that was NOT used as the + /// primary uv). This gives each face its own copy of a seam vertex, + /// enabling a proper GPU texture atlas without vertex duplication. + /// + /// Size = mesh.number_of_halfedges(). Border halfedges = (0,0). + std::vector halfedge_uv; + + bool success = false; + bool has_seam = false; ///< true when a vertex was reached via two paths }; struct Layout3D { - std::vector pos; ///< pos[v.idx()] = 3-D position on S² + std::vector pos; bool success = false; bool has_seam = false; }; -/// Per-cut-edge holonomy for the Euclidean case. -/// translations[i] is the translation ω associated with cut_edge_indices[i] -/// of the CutGraph. For a flat torus, ω_1 and ω_2 are the lattice generators. +/// Per-cut-edge holonomy. +/// +/// Euclidean: translations[i] = ω_i (translation vector for cut_edge_indices[i]). +/// For a flat torus ω_1, ω_2 are the lattice generators; τ = ω_2/ω_1 ∈ ℍ. +/// +/// Hyperbolic: mobius_maps[i] = T_i (Möbius isometry of the Poincaré disk). +/// T_i(p_actual) = p_virtual — maps the placed vertex position to the +/// trilaterated virtual position obtained by continuing the unfolding across +/// the cut. struct HolonomyData { - std::vector translations; // one per cut edge + std::vector translations; ///< Euclidean / spherical + std::vector mobius_maps; ///< hyperbolic (Phase 7) std::vector cut_edge_indices; }; @@ -77,9 +166,23 @@ struct HolonomyData { namespace detail { // ───────────────────────────────────────────────────────────────────────────── -// Euclidean trilateration -// Given placed p_a, p_b and distances d_a (from p_a), d_b (from p_b), -// return the point on the LEFT (CCW) side of the directed edge p_a → p_b. +// Priority-BFS queue entry +// Faces closer to the root (lower depth) are processed first, reducing the +// accumulation of trilateration errors for later-placed vertices. +// ───────────────────────────────────────────────────────────────────────────── +struct BFSEntry { + int depth; ///< BFS depth = max(depth[v_src], depth[v_tgt]) + 1 + Halfedge_index h; ///< halfedge pointing INTO the face to be placed + /// min-heap: smaller depth = higher priority + bool operator>(const BFSEntry& o) const { return depth > o.depth; } +}; + +using BFSQueue = std::priority_queue, + std::greater>; + +// ───────────────────────────────────────────────────────────────────────────── +// Euclidean trilateration — LEFT (CCW) side of p_a → p_b // ───────────────────────────────────────────────────────────────────────────── inline Eigen::Vector2d trilaterate_2d( const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, @@ -89,20 +192,15 @@ inline Eigen::Vector2d trilaterate_2d( double d = ab.norm(); if (d < 1e-14) return pa; Eigen::Vector2d e1 = ab / d; - Eigen::Vector2d e2(-e1.y(), e1.x()); // CCW perpendicular + Eigen::Vector2d e2(-e1.y(), e1.x()); 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 + return pa + t*e1 + s*e2; } // ───────────────────────────────────────────────────────────────────────────── -// Spherical trilateration -// Given unit vectors pa, pb on S² and arc-lengths da, db to the new point, -// return the unit vector on the LEFT side of the geodesic arc pa → pb. -// -// Solution: p_c = α·pa + β·pb + γ·(pa × pb), γ > 0. -// Constraints: pa·p_c = cos(da), pb·p_c = cos(db), |p_c| = 1. +// Spherical trilateration (S²) — LEFT side of geodesic arc pa → pb // ───────────────────────────────────────────────────────────────────────────── inline Eigen::Vector3d trilaterate_sph( const Eigen::Vector3d& pa, const Eigen::Vector3d& pb, @@ -111,317 +209,360 @@ inline Eigen::Vector3d trilaterate_sph( double c = pa.dot(pb); double denom = 1.0 - c*c; if (denom < 1e-14) return pa; - 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 cn2 = 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; - + double gamma = (gamma2 > 0.0 && cn2 > 1e-28) ? std::sqrt(gamma2 / cn2) : 0.0; Eigen::Vector3d p = alpha*pa + beta*pb + gamma*cross; double n = p.norm(); return (n > 1e-14) ? (p / n) : pa; } // ───────────────────────────────────────────────────────────────────────────── -// Hyperbolic trilateration — exact via Möbius map + hyperbolic law of cosines. -// -// pa, pb: Poincaré disk coordinates of two already-placed vertices -// D: hyperbolic distance pa → pb (from the DOF vector) -// da: hyperbolic distance pa → pc (new vertex) -// db: hyperbolic distance pb → pc (new vertex) -// -// Returns the Poincaré disk coordinate of pc on the LEFT (CCW) side of pa→pb. -// -// Algorithm: -// 1. Map pa → 0 via Möbius T: z ↦ (z−pa)/(1−conj(pa)·z) -// 2. θ_b = arg(T(pb)) — direction from origin towards T(pb) -// 3. Angle α at pa from the hyperbolic law of cosines: -// cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α) -// → cos(α) = (cosh(da)·cosh(D) − cosh(db)) / (sinh(da)·sinh(D)) -// 4. pc in origin frame: tanh(da/2)·exp(i·(θ_b + α)) [left = +α] -// 5. Map back: T⁻¹(w) = (w + pa)/(1 + conj(pa)·w) +// Hyperbolic trilateration — exact via Möbius + hyperbolic law of cosines. +// LEFT (CCW) side of pa → pb in the Poincaré disk. // ───────────────────────────────────────────────────────────────────────────── inline Eigen::Vector2d trilaterate_hyp( const Eigen::Vector2d& pa, const Eigen::Vector2d& pb, - double D, // hyperbolic distance pa → pb - double da, // hyperbolic distance pa → pc - double db) // hyperbolic distance pb → pc + double D, double da, double db) { using C = std::complex; - - // Degenerate guard: very short edges or coincident points if (D < 1e-12 || da < 1e-12) return pa; if (std::sinh(da) * std::sinh(D) < 1e-14) return pa; - - C a(pa.x(), pa.y()); - C b(pb.x(), pb.y()); - - // Möbius map: T_a(z) = (z − a) / (1 − conj(a)·z) - auto mobius_fwd = [](C z, C center) -> C { - return (z - center) / (C(1.0) - std::conj(center) * z); - }; - // Inverse: T_a⁻¹(w) = (w + a) / (1 + conj(a)·w) - auto mobius_inv = [](C w, C center) -> C { - return (w + center) / (C(1.0) + std::conj(center) * w); - }; - - C b_mapped = mobius_fwd(b, a); - double theta_b = std::arg(b_mapped); - - // Hyperbolic law of cosines for angle α at pa: - // cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α) - double cos_alpha = (std::cosh(da) * std::cosh(D) - std::cosh(db)) - / (std::sinh(da) * std::sinh(D)); + C a(pa.x(), pa.y()), b(pb.x(), pb.y()); + auto fwd = [](C z, C c_) { return (z-c_)/(C(1)-std::conj(c_)*z); }; + auto inv = [](C w, C c_) { return (w+c_)/(C(1)+std::conj(c_)*w); }; + double theta_b = std::arg(fwd(b, a)); + double cos_alpha = (std::cosh(da)*std::cosh(D) - std::cosh(db)) + / (std::sinh(da)*std::sinh(D)); cos_alpha = std::max(-1.0, std::min(1.0, cos_alpha)); - double alpha = std::acos(cos_alpha); // ∈ [0, π], sin > 0 for left side - - // pc in the frame where pa = 0 and pb is on the positive real axis, - // then rotate back by theta_b: - C pc_origin = std::tanh(da * 0.5) * std::exp(C(0.0, theta_b + alpha)); - - // Map back to original disk - C pc_c = mobius_inv(pc_origin, a); + C pc_origin = std::tanh(da*0.5) * std::exp(C(0.0, theta_b + std::acos(cos_alpha))); + C pc_c = inv(pc_origin, a); return Eigen::Vector2d(pc_c.real(), pc_c.imag()); } // ───────────────────────────────────────────────────────────────────────────── -// Möbius centering of a Poincaré-disk layout -// Applies T_c(z) = (z − c)/(1 − conj(c)·z) where c is the Euclidean centroid -// of all vertex positions. Maps c → 0, i.e. centers the cloud in the disk. +// 3-D face area (cross product / 2) +// ───────────────────────────────────────────────────────────────────────────── +inline double face_3d_area(const ConformalMesh& mesh, Face_index f) +{ + Halfedge_index h = mesh.halfedge(f); + auto pA = mesh.point(mesh.source(h)); + auto pB = mesh.point(mesh.target(h)); + auto pC = mesh.point(mesh.target(mesh.next(h))); + Eigen::Vector3d ab(pB.x()-pA.x(), pB.y()-pA.y(), pB.z()-pA.z()); + Eigen::Vector3d ac(pC.x()-pA.x(), pC.y()-pA.y(), pC.z()-pA.z()); + return 0.5 * ab.cross(ac).norm(); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Best root face: largest 3-D area, with 1.5× bonus for interior faces. +// ───────────────────────────────────────────────────────────────────────────── +inline Face_index best_root_face(const ConformalMesh& mesh) +{ + Face_index best; + double best_score = -1.0; + for (auto f : mesh.faces()) { + double area = face_3d_area(mesh, f); + bool interior = true; + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + if (mesh.is_border(mesh.opposite(h))) { interior = false; break; } + double score = area * (interior ? 1.5 : 1.0); + if (score > best_score) { best_score = score; best = f; } + } + return best; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Euclidean bounding box of placed vertices (for multi-component offset) +// ───────────────────────────────────────────────────────────────────────────── +inline void bbox_2d(const std::vector& uv, + const std::vector& placed, + double& xmax) +{ + xmax = 0.0; + for (std::size_t i = 0; i < uv.size(); ++i) + if (placed[i]) xmax = std::max(xmax, uv[i].x()); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Möbius centering — unweighted (fallback / Phase 6 compat.) // ───────────────────────────────────────────────────────────────────────────── inline void center_poincare_disk(std::vector& uv) { if (uv.empty()) return; using C = std::complex; - - // Euclidean centroid (good enough for moderate displacements from origin) - C centroid(0.0, 0.0); - for (auto& p : uv) centroid += C(p.x(), p.y()); - centroid /= static_cast(uv.size()); - - // Clamp to strictly inside the disk - double r = std::abs(centroid); - if (r > 0.999) centroid *= (0.999 / r); - if (r < 1e-12) return; // already centered - + C c(0); + for (auto& p : uv) c += C(p.x(), p.y()); + c /= static_cast(uv.size()); + double r = std::abs(c); + if (r > 0.999) c *= 0.999/r; + if (r < 1e-12) return; for (auto& p : uv) { C z(p.x(), p.y()); - C w = (z - centroid) / (C(1.0) - std::conj(centroid) * z); - p = Eigen::Vector2d(w.real(), w.imag()); + C w = (z-c)/(C(1)-std::conj(c)*z); + p = {w.real(), w.imag()}; + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Face-area-weighted iterative Möbius centering (Fréchet mean, Phase 7). +// weights[i] = Voronoi area of vertex i. +// ───────────────────────────────────────────────────────────────────────────── +inline void center_poincare_disk_weighted( + std::vector& uv, const std::vector& weights) +{ + if (uv.empty()) return; + using C = std::complex; + double total_w = 0.0; for (double w : weights) total_w += w; + if (total_w < 1e-14) { center_poincare_disk(uv); return; } + C c(0); + for (int iter = 0; iter < 30; ++iter) { + C mt(0); + for (std::size_t i = 0; i < uv.size(); ++i) { + C z(uv[i].x(), uv[i].y()); + mt += weights[i] * (z-c)/(C(1)-std::conj(c)*z); + } + mt /= total_w; + C c_new = (mt+c)/(C(1)+std::conj(c)*mt); + if (std::abs(c_new-c) < 1e-12) { c = c_new; break; } + c = c_new; + } + double r = std::abs(c); + if (r > 0.999) c *= 0.999/r; + if (r < 1e-12) return; + for (auto& p : uv) { + C z(p.x(), p.y()); + C w = (z-c)/(C(1)-std::conj(c)*z); + p = {w.real(), w.imag()}; } } } // namespace detail +// ── Vertex Voronoi area weights ─────────────────────────────────────────────── +inline std::vector compute_vertex_area_weights(const ConformalMesh& mesh) +{ + std::vector w(mesh.number_of_vertices(), 0.0); + for (auto f : mesh.faces()) { + double area = detail::face_3d_area(mesh, f); + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + w[static_cast(mesh.target(h).idx())] += area / 3.0; + } + return w; +} + // ── Layout normalisation ────────────────────────────────────────────────────── -/// Euclidean: translate centroid to origin, rotate major axis to x-axis (PCA). inline void normalise_euclidean(Layout2D& layout) { if (!layout.success || layout.uv.empty()) return; const std::size_t n = layout.uv.size(); - - // Centroid Eigen::Vector2d mean = Eigen::Vector2d::Zero(); for (auto& p : layout.uv) mean += p; mean /= static_cast(n); for (auto& p : layout.uv) p -= mean; - - // PCA: 2×2 covariance + for (auto& p : layout.halfedge_uv) p -= mean; Eigen::Matrix2d cov = Eigen::Matrix2d::Zero(); for (auto& p : layout.uv) cov += p * p.transpose(); cov /= static_cast(n); - Eigen::SelfAdjointEigenSolver eig(cov); - // Eigenvectors in ascending order — we want the largest (index 1) Eigen::Vector2d major = eig.eigenvectors().col(1); - // Rotation that aligns major axis with x-axis double angle = -std::atan2(major.y(), major.x()); Eigen::Matrix2d R; R << std::cos(angle), -std::sin(angle), std::sin(angle), std::cos(angle); - for (auto& p : layout.uv) p = R * p; + for (auto& p : layout.uv) p = R * p; + for (auto& p : layout.halfedge_uv) p = R * p; } -/// Hyperbolic: Möbius map centering the Euclidean centroid to the disk origin. -inline void normalise_hyperbolic(Layout2D& layout) +/// Hyperbolic — face-area-weighted iterative Möbius centering. +inline void normalise_hyperbolic(Layout2D& layout, const ConformalMesh& mesh) +{ + if (!layout.success || layout.uv.empty()) return; + auto weights = compute_vertex_area_weights(mesh); + detail::center_poincare_disk_weighted(layout.uv, weights); + // halfedge_uv follows the same Möbius map + detail::center_poincare_disk(layout.halfedge_uv); +} +inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh { if (!layout.success || layout.uv.empty()) return; detail::center_poincare_disk(layout.uv); + detail::center_poincare_disk(layout.halfedge_uv); } -/// Spherical: rotate so the Euclidean centroid of vertex positions points -/// towards the north pole (0, 0, 1). inline void normalise_spherical(Layout3D& layout) { if (!layout.success || layout.pos.empty()) return; Eigen::Vector3d mean = Eigen::Vector3d::Zero(); for (auto& p : layout.pos) mean += p; - double len = mean.norm(); - if (len < 1e-12) return; - mean /= len; // unit vector of mean direction - - Eigen::Vector3d north(0.0, 0.0, 1.0); + double len = mean.norm(); if (len < 1e-12) return; + mean /= len; + Eigen::Vector3d north(0, 0, 1); Eigen::Vector3d axis = mean.cross(north); - double sin_a = axis.norm(); - double cos_a = mean.dot(north); - if (sin_a < 1e-12) return; // already at north (or south) pole + double sin_a = axis.norm(), cos_a = mean.dot(north); + if (sin_a < 1e-12) return; axis /= sin_a; - - // Rodrigues rotation matrix Eigen::Matrix3d K; - K << 0.0, -axis.z(), axis.y(), - axis.z(), 0.0, -axis.x(), - -axis.y(), axis.x(), 0.0; - Eigen::Matrix3d R = Eigen::Matrix3d::Identity() - + sin_a * K + (1.0 - cos_a) * K * K; - - for (auto& p : layout.pos) p = R * p; + K << 0, -axis.z(), axis.y(), + axis.z(), 0, -axis.x(), + -axis.y(), axis.x(), 0; + Eigen::Matrix3d Rot = Eigen::Matrix3d::Identity() + sin_a*K + (1-cos_a)*K*K; + for (auto& p : layout.pos) p = Rot * p; } -// ── Euclidean layout ────────────────────────────────────────────────────────── -// -// l_ij = exp( (λ°_ij + u_i + u_j + λ_e) / 2 ) -// -// Optional CutGraph: if non-null, cut edges are treated as boundary; -// holonomy translations are computed for each cut edge and returned via -// *holonomy (if non-null). // ───────────────────────────────────────────────────────────────────────────── +// BFS placement helpers shared across all three layout functions. +// +// set_face_huv: populate halfedge_uv for the three halfedges of face f, +// given that the face was entered via halfedge h_enter. +// h_enter: source=v_src, target=v_tgt → huv = uv[v_src] +// next(h_enter): source=v_tgt, target=v_new → huv = uv[v_tgt] +// prev(h_enter): source=v_new, target=v_src → huv = p_new (trilaterated) +// ───────────────────────────────────────────────────────────────────────────── +namespace detail { + +inline void set_face_huv_2d( + std::vector& huv, + const ConformalMesh& mesh, + Halfedge_index h, + const std::vector& uv, + const Eigen::Vector2d& p_new) +{ + auto s = [](std::size_t x){ return x; }; + huv[s(static_cast(h.idx()))] = uv[static_cast(mesh.source(h).idx())]; + huv[s(static_cast(mesh.next(h).idx()))] = uv[static_cast(mesh.target(h).idx())]; + huv[s(static_cast(mesh.prev(h).idx()))] = p_new; +} + +inline void set_root_huv_2d( + std::vector& huv, + const ConformalMesh& mesh, + Face_index f, + const std::vector& uv) +{ + for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) + huv[static_cast(hf.idx())] = uv[static_cast(mesh.source(hf).idx())]; +} + +} // namespace detail + +// ── Euclidean layout ────────────────────────────────────────────────────────── inline Layout2D euclidean_layout( ConformalMesh& mesh, const std::vector& x, const EuclideanMaps& maps, - const CutGraph* cut = nullptr, - HolonomyData* holonomy = nullptr, + const CutGraph* cut = nullptr, + HolonomyData* holonomy = nullptr, bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); + const std::size_t nf = mesh.number_of_faces(); + const std::size_t nh = mesh.number_of_halfedges(); Layout2D result; result.uv.assign(nv, Eigen::Vector2d::Zero()); - if (mesh.number_of_faces() == 0) return result; + result.halfedge_uv.assign(nh, Eigen::Vector2d::Zero()); + if (nf == 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 { + auto get_u = [&](Vertex_index v) { + int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast(iv)] : 0.0; }; + auto get_ue = [&](Edge_index e) { + int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast(ie)] : 0.0; }; + auto edge_len = [&](Halfedge_index h) { 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); - }; + 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); + std::vector face_placed(nf, false); + std::vector vertex_depth(nv, std::numeric_limits::max()); - // ── 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); - double lBC = edge_len(h1); + auto place_component = [&](Face_index f_root, Eigen::Vector2d offset) { + Halfedge_index h0 = mesh.halfedge(f_root); + Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1); + Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2); + double lAB = edge_len(h0), lBC = edge_len(h1), lCA = edge_len(h2); + result.uv[vA.idx()] = offset; + result.uv[vB.idx()] = offset + 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; + vertex_depth[vA.idx()] = vertex_depth[vB.idx()] = vertex_depth[vC.idx()] = 0; + face_placed[f_root.idx()] = true; + detail::set_root_huv_2d(result.halfedge_uv, mesh, f_root, result.uv); - 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 (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - Halfedge_index h_opp = mesh.opposite(h); - if (mesh.is_border(h_opp)) continue; - // Treat cut edges as boundary (don't cross them in BFS) - if (cut && cut->is_cut(mesh.edge(h))) continue; - Face_index f_adj = mesh.face(h_opp); - if (!face_placed[f_adj.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)), - edge_len(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; - - // ── Holonomy: compute translation for each cut edge ─────────────────────── - // For each cut edge e = (u,v), look at the face on the far side (h_opp). - // Trilaterate the third vertex w of that face from uv[u] and uv[v], - // and compare to the actual position uv[w]. - // Translation ω = trilaterated_position(w) − uv[w]. - if (cut && holonomy) { - holonomy->cut_edge_indices = cut->cut_edge_indices; - holonomy->translations.reserve(cut->cut_edge_indices.size()); - - for (std::size_t ce_idx : cut->cut_edge_indices) { - Edge_index e = *std::next(mesh.edges().begin(), - static_cast(ce_idx)); - Halfedge_index h = mesh.halfedge(e); - Halfedge_index h_opp = mesh.opposite(h); - - // Choose the halfedge whose face was NOT placed during BFS - // (the one across the cut from the main BFS sweep). - // If both sides were placed, use h_opp; if neither, skip. - Halfedge_index h_cross = h_opp; - if (mesh.is_border(h_cross)) h_cross = h; - if (mesh.is_border(h_cross)) { - holonomy->translations.push_back(Eigen::Vector2d::Zero()); - continue; + detail::BFSQueue pq; + auto enqueue = [&](Face_index f) { + for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index ho = mesh.opposite(hf); + if (mesh.is_border(ho)) continue; + if (cut && cut->is_cut(mesh.edge(hf))) continue; + Face_index fadj = mesh.face(ho); + if (!face_placed[fadj.idx()]) { + int d = std::max(vertex_depth[mesh.source(ho).idx()], + vertex_depth[mesh.target(ho).idx()]) + 1; + pq.push({d, ho}); + } } + }; + enqueue(f_root); - Vertex_index v_src = mesh.source(h_cross); - Vertex_index v_tgt = mesh.target(h_cross); - Vertex_index v_new = mesh.target(mesh.next(h_cross)); + while (!pq.empty()) { + auto [depth, h] = pq.top(); pq.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + + Vertex_index v_src = mesh.source(h), 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)), edge_len(mesh.next(h))); if (!vertex_placed[v_new.idx()]) { - holonomy->translations.push_back(Eigen::Vector2d::Zero()); - continue; + result.uv[v_new.idx()] = p; + vertex_placed[v_new.idx()] = true; + vertex_depth[v_new.idx()] = depth; + } else { + result.has_seam = true; } + // halfedge_uv: p is the UV of v_new in THIS face (may differ from uv[v_new] at seam) + detail::set_face_huv_2d(result.halfedge_uv, mesh, h, result.uv, p); + face_placed[f.idx()] = true; + enqueue(f); + } + }; + place_component(detail::best_root_face(mesh), Eigen::Vector2d::Zero()); + for (auto f : mesh.faces()) { + if (face_placed[f.idx()]) continue; + double xmax; detail::bbox_2d(result.uv, vertex_placed, xmax); + place_component(f, Eigen::Vector2d(xmax * 1.1 + 1.0, 0.0)); + } + result.success = true; + + // ── Holonomy ────────────────────────────────────────────────────────────── + if (cut && holonomy) { + holonomy->cut_edge_indices = cut->cut_edge_indices; + holonomy->translations.clear(); + holonomy->translations.reserve(cut->cut_edge_indices.size()); + holonomy->mobius_maps.clear(); + for (std::size_t ce_idx : cut->cut_edge_indices) { + Edge_index e = *std::next(mesh.edges().begin(), static_cast(ce_idx)); + Halfedge_index h = mesh.halfedge(e); + Halfedge_index ho = mesh.opposite(h); + Halfedge_index hx = mesh.is_border(ho) ? h : ho; + if (mesh.is_border(hx)) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } + Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx)); + if (!vertex_placed[vn.idx()]) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } Eigen::Vector2d p_tri = detail::trilaterate_2d( - result.uv[v_src.idx()], result.uv[v_tgt.idx()], - edge_len(mesh.prev(h_cross)), - edge_len(mesh.next(h_cross))); - - holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]); + result.uv[vs.idx()], result.uv[vt.idx()], + edge_len(mesh.prev(hx)), edge_len(mesh.next(hx))); + // Store seam UV for the cut-crossing halfedges + detail::set_face_huv_2d(result.halfedge_uv, mesh, hx, result.uv, p_tri); + holonomy->translations.push_back(p_tri - result.uv[vn.idx()]); } } @@ -430,103 +571,105 @@ inline Layout2D euclidean_layout( } // ── Spherical layout ────────────────────────────────────────────────────────── -// -// l_ij = 2·arcsin( min(exp(Λ_ij / 2), 1) ) -// -// ───────────────────────────────────────────────────────────────────────────── inline Layout3D spherical_layout( ConformalMesh& mesh, const std::vector& x, const SphericalMaps& maps, + const CutGraph* cut = nullptr, + HolonomyData* holonomy = nullptr, bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); + const std::size_t nf = mesh.number_of_faces(); Layout3D result; result.pos.assign(nv, Eigen::Vector3d::Zero()); - if (mesh.number_of_faces() == 0) return result; + if (nf == 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 { + auto get_u = [&](Vertex_index v) { + int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast(iv)] : 0.0; }; + auto get_ue = [&](Edge_index e) { + int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast(ie)] : 0.0; }; + auto arc_len = [&](Halfedge_index h) { 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)); - }; + 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)); }; std::vector vertex_placed(nv, false); - std::vector face_placed(mesh.number_of_faces(), false); + std::vector face_placed(nf, false); + std::vector vertex_depth(nv, std::numeric_limits::max()); - 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); + auto place_component = [&](Face_index f_root) { + Halfedge_index h0 = mesh.halfedge(f_root); + Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1); + Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2); + double lAB = arc_len(h0), lBC = arc_len(h1), lCA = arc_len(h2); + result.pos[vA.idx()] = Eigen::Vector3d(0, 0, 1); + result.pos[vB.idx()] = Eigen::Vector3d(std::sin(lAB), 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; + vertex_depth[vA.idx()] = vertex_depth[vB.idx()] = vertex_depth[vC.idx()] = 0; + face_placed[f_root.idx()] = true; - 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; + detail::BFSQueue pq; + auto enqueue = [&](Face_index f) { + for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index ho = mesh.opposite(hf); + if (mesh.is_border(ho)) continue; + if (cut && cut->is_cut(mesh.edge(hf))) continue; + Face_index fadj = mesh.face(ho); + if (!face_placed[fadj.idx()]) { + int d = std::max(vertex_depth[mesh.source(ho).idx()], + vertex_depth[mesh.target(ho).idx()]) + 1; + pq.push({d, ho}); + } + } + }; + enqueue(f_root); - std::queue q; - auto enqueue = [&](Face_index f) { - for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - Halfedge_index h_opp = mesh.opposite(h); - if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()]) - q.push(h_opp); + while (!pq.empty()) { + auto [depth, h] = pq.top(); pq.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + Vertex_index vs = mesh.source(h), vt = mesh.target(h), vn = mesh.target(mesh.next(h)); + Eigen::Vector3d p = detail::trilaterate_sph( + result.pos[vs.idx()], result.pos[vt.idx()], + arc_len(mesh.prev(h)), arc_len(mesh.next(h))); + if (!vertex_placed[vn.idx()]) { + result.pos[vn.idx()] = p; vertex_placed[vn.idx()] = true; vertex_depth[vn.idx()] = depth; + } else result.has_seam = true; + face_placed[f.idx()] = true; enqueue(f); } }; - enqueue(f0); - while (!q.empty()) { - Halfedge_index h = q.front(); q.pop(); - Face_index f = mesh.face(h); - if (face_placed[f.idx()]) continue; + place_component(detail::best_root_face(mesh)); + for (auto f : mesh.faces()) if (!face_placed[f.idx()]) place_component(f); + result.success = true; - 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; + if (cut && holonomy) { + holonomy->cut_edge_indices = cut->cut_edge_indices; + holonomy->translations.clear(); + holonomy->translations.reserve(cut->cut_edge_indices.size()); + holonomy->mobius_maps.clear(); + for (std::size_t ce_idx : cut->cut_edge_indices) { + Edge_index e = *std::next(mesh.edges().begin(), static_cast(ce_idx)); + Halfedge_index hx = mesh.is_border(mesh.opposite(mesh.halfedge(e))) + ? mesh.halfedge(e) : mesh.opposite(mesh.halfedge(e)); + if (mesh.is_border(hx)) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } + Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx)); + if (!vertex_placed[vn.idx()]) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } + Eigen::Vector3d p_tri = detail::trilaterate_sph( + result.pos[vs.idx()], result.pos[vt.idx()], + arc_len(mesh.prev(hx)), arc_len(mesh.next(hx))); + Eigen::Vector3d diff = p_tri - result.pos[vn.idx()]; + holonomy->translations.push_back(Eigen::Vector2d(diff.x(), diff.y())); } - face_placed[f.idx()] = true; - enqueue(f); } - result.success = true; if (normalise) normalise_spherical(result); return result; } // ── HyperIdeal layout (Poincaré disk) — exact trilateration ────────────────── -// -// Hyperbolic edge distance: -// cosh(d_ij) = cosh(b_i + a_ij/2)·cosh(b_j + a_ij/2) − sinh(b_i)·sinh(b_j) -// -// Trilateration via Möbius map + hyperbolic law of cosines (exact, Phase 6). -// -// Optional CutGraph and HolonomyData for closed meshes. -// ───────────────────────────────────────────────────────────────────────────── inline Layout2D hyper_ideal_layout( ConformalMesh& mesh, const std::vector& x, @@ -536,165 +679,128 @@ inline Layout2D hyper_ideal_layout( bool normalise = false) { const std::size_t nv = mesh.number_of_vertices(); + const std::size_t nf = mesh.number_of_faces(); + const std::size_t nh = mesh.number_of_halfedges(); Layout2D result; result.uv.assign(nv, Eigen::Vector2d::Zero()); - if (mesh.number_of_faces() == 0) return result; + result.halfedge_uv.assign(nh, Eigen::Vector2d::Zero()); + if (nf == 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), 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 ch = std::cosh(bi + half_a) * std::cosh(bj + half_a) - - std::sinh(bi) * std::sinh(bj); - return std::acosh(std::max(1.0, ch)); - }; + auto get_b = [&](Vertex_index v) { + int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast(iv)] : 0.0; }; + auto get_a = [&](Edge_index e) { + int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast(ie)] : 0.0; }; + auto hyp_dist = [&](Halfedge_index h) { + double bi = get_b(mesh.source(h)), bj = get_b(mesh.target(h)), a = get_a(mesh.edge(h)); + return std::acosh(std::max(1.0, std::cosh(bi+a*0.5)*std::cosh(bj+a*0.5)-std::sinh(bi)*std::sinh(bj))); }; std::vector vertex_placed(nv, false); - std::vector face_placed(mesh.number_of_faces(), false); + std::vector face_placed(nf, false); + std::vector vertex_depth(nv, std::numeric_limits::max()); - // ── Place root face: vA at origin, vB on positive real axis ────────────── - 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); + auto place_component = [&](Face_index f_root) { + Halfedge_index h0 = mesh.halfedge(f_root); + Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1); + Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2); + double dAB = hyp_dist(h0), dBC = hyp_dist(h1), dCA = hyp_dist(h2); + result.uv[vA.idx()] = Eigen::Vector2d::Zero(); + result.uv[vB.idx()] = Eigen::Vector2d(std::tanh(dAB*0.5), 0.0); + result.uv[vC.idx()] = detail::trilaterate_hyp(result.uv[vA.idx()], result.uv[vB.idx()], dAB, dCA, dBC); + vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true; + vertex_depth[vA.idx()] = vertex_depth[vB.idx()] = vertex_depth[vC.idx()] = 0; + face_placed[f_root.idx()] = true; + detail::set_root_huv_2d(result.halfedge_uv, mesh, f_root, result.uv); - result.uv[vA.idx()] = Eigen::Vector2d(0.0, 0.0); - result.uv[vB.idx()] = Eigen::Vector2d(std::tanh(dAB * 0.5), 0.0); - result.uv[vC.idx()] = detail::trilaterate_hyp( - result.uv[vA.idx()], result.uv[vB.idx()], dAB, dCA, dBC); - vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true; - face_placed[f0.idx()] = true; + detail::BFSQueue pq; + auto enqueue = [&](Face_index f) { + for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index ho = mesh.opposite(hf); + if (mesh.is_border(ho)) continue; + if (cut && cut->is_cut(mesh.edge(hf))) continue; + Face_index fadj = mesh.face(ho); + if (!face_placed[fadj.idx()]) { + int d = std::max(vertex_depth[mesh.source(ho).idx()], + vertex_depth[mesh.target(ho).idx()]) + 1; + pq.push({d, ho}); + } + } + }; + enqueue(f_root); - // ── BFS ─────────────────────────────────────────────────────────────────── - std::queue q; - auto enqueue = [&](Face_index f) { - for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { - Halfedge_index h_opp = mesh.opposite(h); - if (mesh.is_border(h_opp)) continue; - if (cut && cut->is_cut(mesh.edge(h))) continue; - Face_index f_adj = mesh.face(h_opp); - if (!face_placed[f_adj.idx()]) - q.push(h_opp); + while (!pq.empty()) { + auto [depth, h] = pq.top(); pq.pop(); + Face_index f = mesh.face(h); + if (face_placed[f.idx()]) continue; + Vertex_index vs = mesh.source(h), vt = mesh.target(h), vn = mesh.target(mesh.next(h)); + double D = hyp_dist(h), da = hyp_dist(mesh.prev(h)), db = hyp_dist(mesh.next(h)); + Eigen::Vector2d p = detail::trilaterate_hyp(result.uv[vs.idx()], result.uv[vt.idx()], D, da, db); + if (!vertex_placed[vn.idx()]) { + result.uv[vn.idx()] = p; vertex_placed[vn.idx()] = true; vertex_depth[vn.idx()] = depth; + } else result.has_seam = true; + detail::set_face_huv_2d(result.halfedge_uv, mesh, h, result.uv, p); + face_placed[f.idx()] = true; enqueue(f); } }; - 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)); - - double D = hyp_dist(h); // distance v_src → v_tgt - double da = hyp_dist(mesh.prev(h)); // distance v_new → v_src - double db = hyp_dist(mesh.next(h)); // distance v_tgt → v_new - - Eigen::Vector2d p = detail::trilaterate_hyp( - result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db); - - 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); - } + place_component(detail::best_root_face(mesh)); + for (auto f : mesh.faces()) if (!face_placed[f.idx()]) place_component(f); result.success = true; - // ── Holonomy ────────────────────────────────────────────────────────────── + // ── Möbius-map holonomy ─────────────────────────────────────────────────── if (cut && holonomy) { holonomy->cut_edge_indices = cut->cut_edge_indices; - holonomy->translations.reserve(cut->cut_edge_indices.size()); - + holonomy->translations.clear(); + holonomy->mobius_maps.clear(); + holonomy->mobius_maps.reserve(cut->cut_edge_indices.size()); for (std::size_t ce_idx : cut->cut_edge_indices) { - Edge_index e = *std::next(mesh.edges().begin(), - static_cast(ce_idx)); - Halfedge_index h_opp = mesh.opposite(mesh.halfedge(e)); - Halfedge_index h_cross = mesh.is_border(h_opp) - ? mesh.halfedge(e) : h_opp; - if (mesh.is_border(h_cross)) { - holonomy->translations.push_back(Eigen::Vector2d::Zero()); - continue; - } - Vertex_index v_src = mesh.source(h_cross); - Vertex_index v_tgt = mesh.target(h_cross); - Vertex_index v_new = mesh.target(mesh.next(h_cross)); - if (!vertex_placed[v_new.idx()]) { - holonomy->translations.push_back(Eigen::Vector2d::Zero()); - continue; - } - double D = hyp_dist(h_cross); - double da = hyp_dist(mesh.prev(h_cross)); - double db = hyp_dist(mesh.next(h_cross)); - Eigen::Vector2d p_tri = detail::trilaterate_hyp( - result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db); - holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]); + Edge_index e = *std::next(mesh.edges().begin(), static_cast(ce_idx)); + Halfedge_index hcut = mesh.halfedge(e), ho = mesh.opposite(hcut); + Halfedge_index hx = mesh.is_border(ho) ? hcut : ho; + if (mesh.is_border(hx)) { holonomy->mobius_maps.push_back(MobiusMap::identity()); continue; } + Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx)); + if (!vertex_placed[vn.idx()]) { holonomy->mobius_maps.push_back(MobiusMap::identity()); continue; } + double D = hyp_dist(hx), da = hyp_dist(mesh.prev(hx)), db = hyp_dist(mesh.next(hx)); + Eigen::Vector2d p_tri = detail::trilaterate_hyp(result.uv[vs.idx()], result.uv[vt.idx()], D, da, db); + detail::set_face_huv_2d(result.halfedge_uv, mesh, hx, result.uv, p_tri); + using C = std::complex; + holonomy->mobius_maps.push_back(MobiusMap::from_three( + C(result.uv[vs.idx()].x(), result.uv[vs.idx()].y()), + C(result.uv[vs.idx()].x(), result.uv[vs.idx()].y()), + C(result.uv[vt.idx()].x(), result.uv[vt.idx()].y()), + C(result.uv[vt.idx()].x(), result.uv[vt.idx()].y()), + C(result.uv[vn.idx()].x(), result.uv[vn.idx()].y()), + C(p_tri.x(), p_tri.y()))); } } - if (normalise) normalise_hyperbolic(result); + if (normalise) normalise_hyperbolic(result, mesh); return result; } // ── Convenience: save layout as OFF ────────────────────────────────────────── inline void save_layout_off( - const std::string& path, - ConformalMesh& mesh, - const Layout2D& layout) + 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"; - } + 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(); + 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) + 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"; - } + 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(); + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) ofs << " " << mesh.target(h).idx(); ofs << "\n"; } } diff --git a/code/include/period_matrix.hpp b/code/include/period_matrix.hpp new file mode 100644 index 0000000..86647e8 --- /dev/null +++ b/code/include/period_matrix.hpp @@ -0,0 +1,152 @@ +#pragma once +// period_matrix.hpp +// +// Phase 7 — Period matrix for closed surfaces with Euclidean (flat) metric. +// +// For a closed genus-g surface with Euclidean conformal structure the holonomy +// group is generated by 2g translations ω_1, ..., ω_{2g} ∈ ℂ ≅ ℝ². +// +// ─── Genus-1 (flat torus) ──────────────────────────────────────────────────── +// +// The lattice Λ = ℤ·ω_1 ⊕ ℤ·ω_2 determines the conformal type. +// +// Period ratio: τ = ω_2 / ω_1 (as complex numbers) +// +// By convention choose ω_1 such that Im(τ) > 0. +// The conformal modulus / Teichmüller parameter is the SL(2,ℤ)-orbit of τ. +// +// Reduction to fundamental domain {|τ| ≥ 1, −½ ≤ Re(τ) < ½, Im(τ) > 0}: +// S: τ ↦ −1/τ (inversion) +// T: τ ↦ τ + 1 (translation) +// Apply S and T repeatedly until τ is in the fundamental domain. +// +// ─── Genus g > 1 ───────────────────────────────────────────────────────────── +// +// The full period matrix is a g×g complex symmetric matrix Ω with positive +// definite imaginary part (Siegel upper half-space H_g). +// Computing Ω from holonomy data requires integration of holomorphic +// differentials — not implemented here. For g > 1, this function returns +// only the 2×2 block for the first pair of generators. +// +// ─── API ───────────────────────────────────────────────────────────────────── +// +// PeriodData pd = compute_period_matrix(holonomy); +// pd.tau — complex period ratio τ (genus 1) +// pd.omega — holonomy generators as complex numbers (size = 2g) +// pd.in_fundamental_domain — whether τ has been reduced +// +// std::complex reduce_to_fundamental_domain(τ) — apply SL(2,ℤ) + +#include "layout.hpp" +#include +#include +#include +#include +#include + +namespace conformallab { + +// ───────────────────────────────────────────────────────────────────────────── +// PeriodData +// ───────────────────────────────────────────────────────────────────────────── + +struct PeriodData { + /// Lattice generators as complex numbers (one per cut edge). + /// omega[i] = translations[i].x() + i·translations[i].y() + std::vector> omega; + + /// Period ratio τ = omega[1] / omega[0] (genus-1 only). + /// Undefined (NaN) for genus != 1 or if holonomy has fewer than 2 generators. + std::complex tau = std::complex( + std::numeric_limits::quiet_NaN(), 0.0); + + /// True if τ has been reduced to the standard fundamental domain. + bool in_fundamental_domain = false; + + int genus() const { return static_cast(omega.size()) / 2; } +}; + +// ───────────────────────────────────────────────────────────────────────────── +// reduce_to_fundamental_domain +// +// Applies SL(2,ℤ) generators S: τ↦−1/τ and T: τ↦τ+1 to bring τ into +// F = { τ ∈ ℍ : |τ| ≥ 1, −½ ≤ Re(τ) < ½ } +// +// Returns the reduced τ. Throws if Im(τ) ≤ 0 (not in upper half-plane). +// ───────────────────────────────────────────────────────────────────────────── +inline std::complex reduce_to_fundamental_domain(std::complex tau) +{ + if (tau.imag() <= 0.0) { + std::ostringstream msg; + msg << "period_matrix: τ = " << tau.real() << " + " << tau.imag() + << "i is not in the upper half-plane (Im(τ) must be > 0)."; + throw std::domain_error(msg.str()); + } + + // Iterate at most 200 times (convergence is rapid for well-conditioned τ) + for (int k = 0; k < 200; ++k) { + // T step: shift Re(τ) into [−½, ½) + double re = tau.real(); + long n = static_cast(std::floor(re + 0.5)); + tau -= std::complex(static_cast(n), 0.0); + + // S step: if |τ| < 1, apply τ ← −1/τ + if (std::abs(tau) < 1.0 - 1e-12) { + tau = -1.0 / tau; + } else { + break; + } + } + return tau; +} + +// ───────────────────────────────────────────────────────────────────────────── +// is_in_fundamental_domain — check membership in F with tolerance tol. +// ───────────────────────────────────────────────────────────────────────────── +inline bool is_in_fundamental_domain(std::complex tau, double tol = 1e-9) +{ + if (tau.imag() <= 0.0) return false; + if (std::abs(tau.real()) > 0.5 + tol) return false; + if (std::abs(tau) < 1.0 - tol) return false; + return true; +} + +// ───────────────────────────────────────────────────────────────────────────── +// compute_period_matrix +// +// Computes the period data from the Euclidean holonomy translations. +// For genus-1 surfaces, also reduces τ to the fundamental domain. +// ───────────────────────────────────────────────────────────────────────────── +inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = true) +{ + PeriodData pd; + pd.omega.reserve(hol.translations.size()); + for (auto& t : hol.translations) + pd.omega.push_back(std::complex(t.x(), t.y())); + + if (pd.omega.size() < 2) return pd; // need at least 2 generators + + // τ = ω_2 / ω_1 — choose ω_1 such that Im(τ) > 0 + std::complex w1 = pd.omega[0]; + std::complex w2 = pd.omega[1]; + if (std::abs(w1) < 1e-14) return pd; + + std::complex tau = w2 / w1; + if (tau.imag() < 0.0) { + tau = std::conj(tau); // swap orientation + w1 = std::conj(w1); + w2 = std::conj(w2); + pd.omega[0] = w1; + pd.omega[1] = w2; + } + if (tau.imag() < 0.0) return pd; // degenerate + + if (reduce) { + tau = reduce_to_fundamental_domain(tau); + pd.in_fundamental_domain = true; + } + pd.tau = tau; + return pd; +} + +} // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index 2f15d73..43a31a5 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -39,6 +39,10 @@ add_executable(conformallab_cgal_tests # ── Phase 6: Gauss–Bonnet, cut graph, exact trilateration, normalisation test_phase6.cpp + + # ── Phase 7: Java-parity layout — MobiusMap, priority BFS, halfedge_uv, + # period matrix, fundamental domain, tiling + test_phase7.cpp ) target_include_directories(conformallab_cgal_tests SYSTEM PRIVATE diff --git a/code/tests/cgal/test_phase7.cpp b/code/tests/cgal/test_phase7.cpp new file mode 100644 index 0000000..bd2673a --- /dev/null +++ b/code/tests/cgal/test_phase7.cpp @@ -0,0 +1,485 @@ +// test_phase7.cpp +// +// Phase 7 — Tests for Java-parity layout features: +// - MobiusMap : identity, inverse, compose, from_three, is_identity +// - best_root_face : selects a valid face; interior bonus +// - halfedge_uv : size, non-seam consistency, seam divergence +// - Priority BFS : vertex ordering / depth correctness +// - normalise_euclidean : halfedge_uv centroid at origin +// - period_matrix.hpp : τ in upper half-plane, SL(2,ℤ) reduction +// - fundamental_domain.hpp: parallelogram CCW, generators, tiling_copy + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "euclidean_functional.hpp" +#include "hyper_ideal_functional.hpp" +#include "newton_solver.hpp" +#include "layout.hpp" +#include "period_matrix.hpp" +#include "fundamental_domain.hpp" +#include +#include +#include +#include +#include + +using namespace conformallab; +using C = std::complex; + +// ════════════════════════════════════════════════════════════════════════════ +// MobiusMap +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MobiusMap, Identity_AppliesAsIdentity) +{ + MobiusMap id = MobiusMap::identity(); + C z(0.3, 0.7); + C w = id.apply(z); + EXPECT_NEAR(w.real(), z.real(), 1e-12); + EXPECT_NEAR(w.imag(), z.imag(), 1e-12); +} + +TEST(MobiusMap, Identity_IsIdentity) +{ + EXPECT_TRUE(MobiusMap::identity().is_identity()); +} + +TEST(MobiusMap, NonIdentity_IsNotIdentity) +{ + // T(z) = z + 1 — translation, clearly not identity + MobiusMap T{ C(1), C(1), C(0), C(1) }; + EXPECT_FALSE(T.is_identity()); +} + +TEST(MobiusMap, Inverse_ComposeIsIdentity) +{ + // T(z) = (2z + 1) / (z + 3) + MobiusMap T{ C(2), C(1), C(1), C(3) }; + MobiusMap TinvT = T.inverse().compose(T); + EXPECT_TRUE(TinvT.is_identity(1e-9)); +} + +TEST(MobiusMap, Compose_OrderCorrect) +{ + // S: z ↦ z + 1, T: z ↦ 2z + // S.compose(T) means S applied after T: z ↦ 2z + 1 + MobiusMap S{ C(1), C(1), C(0), C(1) }; // z + 1 + MobiusMap T{ C(2), C(0), C(0), C(1) }; // 2z + MobiusMap ST = S.compose(T); + C z(1.0, 0.0); + // S(T(z)) = S(2) = 3 + EXPECT_NEAR(ST.apply(z).real(), 3.0, 1e-12); + EXPECT_NEAR(ST.apply(z).imag(), 0.0, 1e-12); +} + +TEST(MobiusMap, FromThree_RecoversMap) +{ + // Known map T(z) = (z + i) / (1 + 0·z) — translation by i + C w1 = C(0, 1) + C(0, 1); // T(i) = 2i + C w2 = C(1, 0) + C(0, 1); // T(1) = 1 + i + C w3 = C(-1, 0) + C(0, 1); // T(-1) = -1 + i + MobiusMap T = MobiusMap::from_three(C(0, 1), w1, C(1, 0), w2, C(-1, 0), w3); + // Verify T maps a fourth point correctly: T(0) = i + C result = T.apply(C(0, 0)); + EXPECT_NEAR(result.real(), 0.0, 1e-9); + EXPECT_NEAR(result.imag(), 1.0, 1e-9); +} + +TEST(MobiusMap, FromThree_DegenerateReturnsIdentity) +{ + // Three coincident points → singular system → identity fallback + C z(0.5, 0.5); + MobiusMap T = MobiusMap::from_three(z, z, z, z, z, z); + // Should not crash; returns identity (or at least a valid map) + // We just check the result is finite + C w = T.apply(C(0.1, 0.2)); + EXPECT_FALSE(std::isnan(w.real())); + EXPECT_FALSE(std::isnan(w.imag())); +} + +TEST(MobiusMap, Apply_Vector2d) +{ + MobiusMap id = MobiusMap::identity(); + Eigen::Vector2d p(0.4, 0.6); + Eigen::Vector2d q = id.apply(p); + EXPECT_NEAR(q.x(), p.x(), 1e-12); + EXPECT_NEAR(q.y(), p.y(), 1e-12); +} + +// ════════════════════════════════════════════════════════════════════════════ +// best_root_face +// ════════════════════════════════════════════════════════════════════════════ + +TEST(BestRootFace, ReturnsValidFace_Triangle) +{ + auto mesh = make_triangle(); + Face_index f = detail::best_root_face(mesh); + EXPECT_NE(f, Face_index()); + EXPECT_GE(f.idx(), 0); +} + +TEST(BestRootFace, ReturnsValidFace_Tetrahedron) +{ + auto mesh = make_tetrahedron(); + Face_index f = detail::best_root_face(mesh); + EXPECT_NE(f, Face_index()); + // Tetrahedron has 4 faces — best is one of them + EXPECT_LT(static_cast(f.idx()), mesh.number_of_faces()); +} + +// ════════════════════════════════════════════════════════════════════════════ +// halfedge_uv — size and non-seam consistency +// ════════════════════════════════════════════════════════════════════════════ + +// Helper: build equilibrium Euclidean layout for a given mesh. +// Uses x = 0 (identity scale factor) which is the equilibrium for natural edge lengths. +static Layout2D make_euclidean_layout(ConformalMesh& mesh) +{ + EuclideanMaps maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + // Pin first vertex (DOF = -1); assign sequential indices to the rest. + auto vit = mesh.vertices().begin(); + maps.v_idx[*vit++] = -1; + int idx = 0; + for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++; + std::vector x(static_cast(idx), 0.0); + return euclidean_layout(mesh, x, maps); +} + +TEST(HalfedgeUV, Size_EqualsNumberOfHalfedges_Triangle) +{ + auto mesh = make_triangle(); + auto lay = make_euclidean_layout(mesh); + EXPECT_EQ(lay.halfedge_uv.size(), mesh.number_of_halfedges()); +} + +TEST(HalfedgeUV, Size_EqualsNumberOfHalfedges_QuadStrip) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + EXPECT_EQ(lay.halfedge_uv.size(), mesh.number_of_halfedges()); +} + +TEST(HalfedgeUV, NonBorderHalfedges_MatchUV) +{ + // For an open mesh with no cut graph the layout has no seams. + // Every non-border halfedge h must satisfy: + // halfedge_uv[h] == uv[source(h)] + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + for (auto h : mesh.halfedges()) { + if (mesh.is_border(h)) continue; + std::size_t hi = static_cast(h.idx()); + std::size_t vi = static_cast(mesh.source(h).idx()); + EXPECT_NEAR(lay.halfedge_uv[hi].x(), lay.uv[vi].x(), 1e-10) + << "halfedge " << hi << " source vertex " << vi; + EXPECT_NEAR(lay.halfedge_uv[hi].y(), lay.uv[vi].y(), 1e-10) + << "halfedge " << hi << " source vertex " << vi; + } +} + +TEST(HalfedgeUV, BorderHalfedges_AreZero) +{ + auto mesh = make_triangle(); + auto lay = make_euclidean_layout(mesh); + bool found_border = false; + for (auto h : mesh.halfedges()) { + if (!mesh.is_border(h)) continue; + std::size_t hi = static_cast(h.idx()); + EXPECT_NEAR(lay.halfedge_uv[hi].x(), 0.0, 1e-12); + EXPECT_NEAR(lay.halfedge_uv[hi].y(), 0.0, 1e-12); + found_border = true; + } + EXPECT_TRUE(found_border); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Priority BFS — depth ordering +// ════════════════════════════════════════════════════════════════════════════ + +TEST(PriorityBFS, Layout_SucceedsOnOpenMesh) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + EXPECT_TRUE(lay.success); +} + +TEST(PriorityBFS, Layout_NoSeamOnOpenMesh) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + EXPECT_FALSE(lay.has_seam); +} + +TEST(PriorityBFS, AllVerticesPlaced) +{ + auto mesh = make_tetrahedron(); + // Tetrahedron is closed; layout without cut graph will have a seam + EuclideanMaps 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++; + std::vector x(static_cast(idx), 0.0); + auto lay = euclidean_layout(mesh, x, maps); + EXPECT_TRUE(lay.success); + // All UVs must be finite + for (auto& p : lay.uv) { + EXPECT_FALSE(std::isnan(p.x())); + EXPECT_FALSE(std::isnan(p.y())); + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// normalise_euclidean — centroid + PCA applied to both uv and halfedge_uv +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NormaliseEuclidean, UVCentroidAtOrigin) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + normalise_euclidean(lay); + + Eigen::Vector2d mean = Eigen::Vector2d::Zero(); + for (auto& p : lay.uv) mean += p; + mean /= static_cast(lay.uv.size()); + EXPECT_NEAR(mean.x(), 0.0, 1e-10); + EXPECT_NEAR(mean.y(), 0.0, 1e-10); +} + +TEST(NormaliseEuclidean, HalfedgeUVCentroidAlsoShifted) +{ + // After normalisation: the non-border halfedge_uv entries should also be + // centred (since they are shifted by the same mean as uv). + // We verify that the mean of non-border halfedge_uv is near (0,0). + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + normalise_euclidean(lay); + + Eigen::Vector2d mean = Eigen::Vector2d::Zero(); + int count = 0; + for (auto h : mesh.halfedges()) { + if (mesh.is_border(h)) continue; + mean += lay.halfedge_uv[static_cast(h.idx())]; + ++count; + } + if (count > 0) mean /= static_cast(count); + EXPECT_NEAR(mean.x(), 0.0, 1e-9); + EXPECT_NEAR(mean.y(), 0.0, 1e-9); +} + +// ════════════════════════════════════════════════════════════════════════════ +// PeriodMatrix — reduce_to_fundamental_domain +// ════════════════════════════════════════════════════════════════════════════ + +TEST(PeriodMatrix, ReduceToFD_AlreadyInFD) +{ + // τ = i is in F (|i|=1, Re(i)=0, Im(i)=1>0) + C tau(0.0, 1.0); + C reduced = reduce_to_fundamental_domain(tau); + EXPECT_TRUE(is_in_fundamental_domain(reduced)); + EXPECT_NEAR(reduced.real(), 0.0, 1e-10); + EXPECT_NEAR(reduced.imag(), 1.0, 1e-10); +} + +TEST(PeriodMatrix, ReduceToFD_ShiftsRealPart) +{ + // τ = 2 + 3i → T step: τ -= 2 → 3i (|3i|=3≥1, Re=0) + C tau(2.0, 3.0); + C reduced = reduce_to_fundamental_domain(tau); + EXPECT_TRUE(is_in_fundamental_domain(reduced, 1e-9)); + EXPECT_NEAR(reduced.real(), 0.0, 1e-10); + EXPECT_NEAR(reduced.imag(), 3.0, 1e-10); +} + +TEST(PeriodMatrix, ReduceToFD_InvertsSmallTau) +{ + // τ = 0.5i → |0.5i|=0.5<1 → S: τ↦-1/(0.5i) = 2i + C tau(0.0, 0.5); + C reduced = reduce_to_fundamental_domain(tau); + EXPECT_TRUE(is_in_fundamental_domain(reduced, 1e-9)); + EXPECT_NEAR(reduced.real(), 0.0, 1e-10); + EXPECT_NEAR(reduced.imag(), 2.0, 1e-10); +} + +TEST(PeriodMatrix, ReduceToFD_ThrowsForNonUpperHalfPlane) +{ + C tau(0.5, -1.0); // Im < 0 → not in upper half-plane + EXPECT_THROW(reduce_to_fundamental_domain(tau), std::domain_error); +} + +TEST(PeriodMatrix, IsInFundamentalDomain_Square) +{ + EXPECT_TRUE(is_in_fundamental_domain(C(0.0, 1.0))); // i + EXPECT_TRUE(is_in_fundamental_domain(C(0.3, 1.5))); // inside + EXPECT_FALSE(is_in_fundamental_domain(C(0.6, 1.5))); // Re > 1/2 + EXPECT_FALSE(is_in_fundamental_domain(C(0.0, 0.5))); // |τ| < 1 +} + +TEST(PeriodMatrix, ComputePeriodMatrix_UnitSquare) +{ + // ω_1 = (1, 0), ω_2 = (0, 1) → τ = i + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.0, 1.0) }; + PeriodData pd = compute_period_matrix(hol, /*reduce=*/false); + EXPECT_EQ(pd.genus(), 1); + EXPECT_GT(pd.tau.imag(), 0.0); + EXPECT_NEAR(pd.tau.real(), 0.0, 1e-10); + EXPECT_NEAR(pd.tau.imag(), 1.0, 1e-10); +} + +TEST(PeriodMatrix, ComputePeriodMatrix_ReducedTau_InFD) +{ + // ω_1 = (1, 0), ω_2 = (0.5, 0.25) → τ = 0.5 + 0.25i + // |τ| = sqrt(0.25 + 0.0625) ≈ 0.559 < 1 → needs S step + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.5, 0.25) }; + PeriodData pd = compute_period_matrix(hol, /*reduce=*/true); + EXPECT_TRUE(pd.in_fundamental_domain); + EXPECT_TRUE(is_in_fundamental_domain(pd.tau, 1e-9)); +} + +// ════════════════════════════════════════════════════════════════════════════ +// FundamentalDomain — genus-1 parallelogram +// ════════════════════════════════════════════════════════════════════════════ + +TEST(FundamentalDomain, Genus1_HasFourVertices) +{ + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.0, 1.0) }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + EXPECT_EQ(fd.vertices.size(), 4u); + EXPECT_TRUE(fd.is_valid()); +} + +TEST(FundamentalDomain, Genus1_VerticesMatchGenerators_UnitSquare) +{ + Eigen::Vector2d w1(1.0, 0.0), w2(0.0, 1.0); + HolonomyData hol; + hol.translations = { w1, w2 }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + // Expected (CCW): origin, w1, w1+w2, w2 + EXPECT_NEAR(fd.vertices[0].x(), 0.0, 1e-12); + EXPECT_NEAR(fd.vertices[0].y(), 0.0, 1e-12); + EXPECT_NEAR(fd.vertices[1].x(), w1.x(), 1e-12); + EXPECT_NEAR(fd.vertices[1].y(), w1.y(), 1e-12); + EXPECT_NEAR(fd.vertices[2].x(), (w1 + w2).x(), 1e-12); + EXPECT_NEAR(fd.vertices[2].y(), (w1 + w2).y(), 1e-12); + EXPECT_NEAR(fd.vertices[3].x(), w2.x(), 1e-12); + EXPECT_NEAR(fd.vertices[3].y(), w2.y(), 1e-12); +} + +TEST(FundamentalDomain, Genus1_CCWOrientation) +{ + // After possible swap, the signed area = cross(v1-v0, v3-v0) > 0 (CCW) + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.0, 1.0) }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + Eigen::Vector2d v0 = fd.vertices[0], v1 = fd.vertices[1], v3 = fd.vertices[3]; + double cross = (v1 - v0).x() * (v3 - v0).y() - (v1 - v0).y() * (v3 - v0).x(); + EXPECT_GT(cross, 0.0); +} + +TEST(FundamentalDomain, Genus1_CCWEnforced_WhenInputIsCW) +{ + // If we give CW generators (w2 × w1 < 0), the polygon must still be CCW. + // w1 = (0,1), w2 = (1,0): cross w1×w2 = 0*0 - 1*1 = -1 < 0 → should swap + HolonomyData hol; + hol.translations = { Eigen::Vector2d(0.0, 1.0), Eigen::Vector2d(1.0, 0.0) }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + Eigen::Vector2d v0 = fd.vertices[0], v1 = fd.vertices[1], v3 = fd.vertices[3]; + double cross = (v1 - v0).x() * (v3 - v0).y() - (v1 - v0).y() * (v3 - v0).x(); + EXPECT_GT(cross, 0.0); +} + +TEST(FundamentalDomain, Genus1_EdgeIdentifications) +{ + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.0, 1.0) }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + EXPECT_EQ(fd.edge_identifications.size(), 2u); + // bottom ≡ top: (0,2) + EXPECT_EQ(fd.edge_identifications[0].first, 0); + EXPECT_EQ(fd.edge_identifications[0].second, 2); + // right ≡ left: (1,3) + EXPECT_EQ(fd.edge_identifications[1].first, 1); + EXPECT_EQ(fd.edge_identifications[1].second, 3); +} + +TEST(FundamentalDomain, Genus1_GeneratorsStored) +{ + Eigen::Vector2d w1(2.0, 1.0), w2(-1.0, 3.0); + HolonomyData hol; + hol.translations = { w1, w2 }; + FundamentalDomain fd = compute_fundamental_domain_genus1(hol); + EXPECT_EQ(fd.generators.size(), 2u); + // Generators are w1 and w2 (possibly swapped to ensure CCW) + // Their sum of norms matches the originals + double norm_gen = fd.generators[0].norm() + fd.generators[1].norm(); + double norm_in = w1.norm() + w2.norm(); + EXPECT_NEAR(norm_gen, norm_in, 1e-10); +} + +TEST(FundamentalDomain, HigherGenus_ReturnsEmpty) +{ + HolonomyData hol; + hol.translations = { + Eigen::Vector2d(1, 0), Eigen::Vector2d(0, 1), + Eigen::Vector2d(2, 0), Eigen::Vector2d(0, 2) // g=2, 4 generators + }; + FundamentalDomain fd = compute_fundamental_domain(hol); + // g > 1 returns empty (TODO Phase 8) + EXPECT_FALSE(fd.is_valid()); +} + +// ════════════════════════════════════════════════════════════════════════════ +// tiling_copy / tiling_neighbourhood +// ════════════════════════════════════════════════════════════════════════════ + +TEST(TilingCopy, ShiftAppliedToAllUV) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + + Eigen::Vector2d w1(3.0, 0.0), w2(0.0, 2.0); + // m=1, n=2 → expected shift = w1 + 2*w2 = (3, 4) + Layout2D copy = tiling_copy(lay, w1, w2, 1, 2); + Eigen::Vector2d expected_shift(3.0, 4.0); + for (std::size_t i = 0; i < lay.uv.size(); ++i) { + EXPECT_NEAR(copy.uv[i].x(), lay.uv[i].x() + expected_shift.x(), 1e-12); + EXPECT_NEAR(copy.uv[i].y(), lay.uv[i].y() + expected_shift.y(), 1e-12); + } +} + +TEST(TilingCopy, ZeroShift_IsSameAsCopy) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + Eigen::Vector2d w1(1, 0), w2(0, 1); + Layout2D copy = tiling_copy(lay, w1, w2, 0, 0); + for (std::size_t i = 0; i < lay.uv.size(); ++i) { + EXPECT_NEAR(copy.uv[i].x(), lay.uv[i].x(), 1e-12); + EXPECT_NEAR(copy.uv[i].y(), lay.uv[i].y(), 1e-12); + } +} + +TEST(TilingNeighbourhood, CorrectCount) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + HolonomyData hol; + hol.translations = { Eigen::Vector2d(1, 0), Eigen::Vector2d(0, 1) }; + // m_max=1, n_max=1 → (2*1+1) * (2*1+1) = 9 tiles + auto tiles = tiling_neighbourhood(lay, hol, 1, 1); + EXPECT_EQ(tiles.size(), 9u); +} + +TEST(TilingNeighbourhood, EmptyHolonomy_ReturnsSingleTile) +{ + auto mesh = make_quad_strip(); + auto lay = make_euclidean_layout(mesh); + HolonomyData hol; // no translations + auto tiles = tiling_neighbourhood(lay, hol); + EXPECT_EQ(tiles.size(), 1u); +}