diff --git a/code/include/cp_euclidean_functional.hpp b/code/include/cp_euclidean_functional.hpp index bbc326e..82efb31 100644 --- a/code/include/cp_euclidean_functional.hpp +++ b/code/include/cp_euclidean_functional.hpp @@ -89,9 +89,25 @@ struct CPEuclideanMaps { CPFMapD phi_f; ///< target face-angle sum (default 2π) }; -// Create the property maps with sensible defaults. -// θ_e = π/2 produces an orthogonal circle packing (Koebe-Andreev-Thurston). -// φ_f = 2π is the natural target for a flat triangle. +/// Attach the three CP-Euclidean property maps to `mesh` with default +/// values and return their handles. +/// +/// Defaults: +/// * `theta_e[e] = π/2` for every edge — orthogonal circle packing +/// (Koebe-Andreev-Thurston). +/// * `phi_f[f] = 2π` for every face — flat target. +/// * `f_idx[f] = -1` for every face — all faces pinned initially; +/// call `assign_cp_euclidean_face_dof_indices()` next to assign +/// DOF indices to all faces except one gauge-pinned face. +/// +/// The maps are named with the `"cf:"` / `"ce:"` prefixes +/// (cf = circle-packing-face, ce = circle-packing-edge) so they do +/// not collide with the Euclidean / Spherical / HyperIdeal maps. +/// +/// \param mesh Input mesh. Modified in place: three property maps are +/// attached if not already present, otherwise the existing +/// maps are returned unchanged (CGAL property-map idempotence). +/// \returns A bundle of all three property maps for caller use. inline CPEuclideanMaps setup_cp_euclidean_maps(ConformalMesh& mesh) { CPEuclideanMaps m; @@ -101,8 +117,18 @@ inline CPEuclideanMaps setup_cp_euclidean_maps(ConformalMesh& mesh) return m; } -// Assign DOF indices 0..n-1 to all faces except `pinned`, which gets −1. -// Mirrors Java CPEuclideanFunctional's convention "skip face index 0". +/// Assign sequential DOF indices `0..n-1` to all faces except `pinned`, +/// which receives the sentinel `-1` (gauge-fixed face, `ρ_pinned = 0`). +/// +/// Mirrors the Java CPEuclideanFunctional's "skip face index 0" +/// convention from `evaluateEnergyAndGradient` (lines 184-185 of +/// CPEuclideanFunctional.java). The C++ port exposes the choice of +/// pinned face explicitly rather than hard-coding it. +/// +/// \param mesh The mesh. Read for face iteration only; not modified. +/// \param m Map bundle whose `f_idx` is overwritten. +/// \param pinned The face whose DOF is fixed at zero (the gauge). +/// \returns The number of free DOFs assigned (`num_faces(mesh) - 1`). inline int assign_cp_euclidean_face_dof_indices(ConformalMesh& mesh, CPEuclideanMaps& m, Face_index pinned) @@ -115,7 +141,9 @@ inline int assign_cp_euclidean_face_dof_indices(ConformalMesh& mesh, return idx; } -// Convenience: pin the first face in iteration order. +/// Convenience overload: pin the **first** face in `mesh.faces()` order. +/// Use this when any face works as the gauge (typically true for +/// closed mesh experiments). inline int assign_cp_euclidean_face_dof_indices(ConformalMesh& mesh, CPEuclideanMaps& m) { @@ -124,6 +152,8 @@ inline int assign_cp_euclidean_face_dof_indices(ConformalMesh& mesh, return assign_cp_euclidean_face_dof_indices(mesh, m, *it); } +/// Count the free DOFs (faces with `f_idx >= 0`). +/// Equivalent to `num_faces(mesh) - `. inline int cp_euclidean_dimension(const ConformalMesh& mesh, const CPEuclideanMaps& m) { diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index cadbcbf..1d9fc86 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -63,8 +63,17 @@ struct EuclideanMaps { EuclEMapD lambda0; ///< base log-length λ°_e (default 0.0) }; -// Create and attach property maps with sensible defaults. -// theta_v = 2π (flat vertex), phi_e = π (interior edge, flat surface). +/// Attach the five Euclidean property maps to `mesh` with sensible +/// defaults and return their handles. +/// +/// Defaults: +/// * `v_idx[v] = -1` (every vertex pinned; user must reassign before solving) +/// * `e_idx[e] = -1` (no edge DOFs by default; use `assign_euclidean_all_dof_indices` for cyclic functional) +/// * `theta_v[v] = 2π` (flat interior vertex target) +/// * `phi_e[e] = π` (interior edge turn angle target — flat surface) +/// * `lambda0[e] = 0` (placeholder; call `compute_euclidean_lambda0_from_mesh` next) +/// +/// Map name prefix: `"ev:"` (vertex) and `"ee:"` (edge). inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh) { EuclideanMaps m; @@ -76,7 +85,12 @@ inline EuclideanMaps setup_euclidean_maps(ConformalMesh& mesh) return m; } -// Assign DOF indices 0..n-1 for all vertices only (no edge DOFs). +/// Assign sequential DOF indices `0..n-1` to all vertices. +/// +/// **Note:** does NOT pin a gauge vertex. For closed meshes the caller +/// must set one `m.v_idx[v] = -1` either before or after this call to +/// remove the rotational mode (the Newton solver's SparseQR fallback +/// will otherwise pick a minimum-norm solution but at higher cost). inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m) { int idx = 0; @@ -84,7 +98,10 @@ inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMap return idx; } -// Assign DOF indices for all vertices AND all edges. +/// Assign DOF indices for all vertices AND all edges (vertex-DOFs first, +/// then edge-DOFs). Use this overload for the "cyclic" formulation that +/// includes per-edge log-length DOFs (`λ_e`) on top of per-vertex scale +/// factors (`u_v`). inline int assign_euclidean_all_dof_indices(ConformalMesh& mesh, EuclideanMaps& m) { int idx = 0; @@ -93,7 +110,7 @@ inline int assign_euclidean_all_dof_indices(ConformalMesh& mesh, EuclideanMaps& return idx; } -// Count variable DOFs (vertices + edges). +/// Count the free DOFs (vertices + edges with index `≥ 0`). inline int euclidean_dimension(const ConformalMesh& mesh, const EuclideanMaps& m) { int dim = 0; diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index 2c24e58..2a50868 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -53,8 +53,20 @@ struct HyperIdealMaps { EMapD theta_e; // target intersection angle θ_e }; -// Add all needed persistent property maps and return handles. -// Defaults: theta_v = 2π (regular cone), theta_e = π (orthogonal circles). +/// Attach the four HyperIdeal property maps to `mesh` and return their +/// handles. +/// +/// Defaults: +/// * `v_idx[v] = -1` (ideal vertex — i.e. the corresponding `b_v` is fixed at 0) +/// * `e_idx[e] = -1` (edge DOF fixed at 0) +/// * `theta_v[v] = 2π` (regular cone target) +/// * `theta_e[e] = π` (orthogonal-circle target) +/// +/// The map prefix `"v:"` / `"e:"` is intentionally generic for the +/// HyperIdeal functional — it is the canonical / Phase 3b model. +/// Other functionals use distinct prefixes (`"ev:"` Euclidean, `"sv:"` +/// Spherical, `"cf:"`/`"ce:"` CP-Euclidean, `"iv:"`/`"ie:"` +/// Inversive-Distance) so all models can coexist on the same mesh. inline HyperIdealMaps setup_hyper_ideal_maps(ConformalMesh& mesh) { HyperIdealMaps m; @@ -65,7 +77,7 @@ inline HyperIdealMaps setup_hyper_ideal_maps(ConformalMesh& mesh) return m; } -// Count variable DOFs: #variable_vertices + #variable_edges. +/// Count free DOFs: `#variable_vertices + #variable_edges`. inline int hyper_ideal_dimension(const ConformalMesh& mesh, const HyperIdealMaps& m) { int dim = 0; @@ -74,8 +86,15 @@ inline int hyper_ideal_dimension(const ConformalMesh& mesh, const HyperIdealMaps return dim; } -// Assign DOF indices 0..n-1: vertices first, then edges. -// All vertices and edges become variable. Returns total DOF count. +/// Make every vertex hyper-ideal and every edge variable, assigning +/// sequential DOF indices `0..n-1` (vertices first, edges after). +/// +/// This is the standard initialisation for the Springborn-2020 +/// hyper-ideal functional — gauge fixing is **not** needed because +/// the energy is strictly convex on the full DOF space (no +/// rotational mode for an all-hyper-ideal configuration). +/// +/// \returns total DOF count = `num_vertices(mesh) + num_edges(mesh)`. inline int assign_all_dof_indices(ConformalMesh& mesh, HyperIdealMaps& m) { int idx = 0; diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index 69e6aed..a932fa5 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -89,7 +89,19 @@ struct InversiveDistanceMaps { IDEMapD I_e; ///< inversive distance I_ij (per edge, constant) }; -// Create the property maps with sensible defaults. +/// Attach the four inversive-distance property maps to `mesh` and +/// return their handles. +/// +/// Defaults are intentionally trivial — every real use of this +/// functional must call `compute_inversive_distance_init_from_mesh()` +/// next to populate `r0` and `I_e` from the input geometry. +/// * `v_idx[v] = -1` (all vertices pinned initially) +/// * `theta_v[v] = 2π` (regular interior vertex) +/// * `r0[v] = 1.0` (placeholder) +/// * `I_e[e] = 1.0` (tangential default — overwritten by init step) +/// +/// The maps use the `"iv:"` / `"ie:"` prefix so they do not collide +/// with the Euclidean / Spherical / HyperIdeal / CP-Euclidean maps. inline InversiveDistanceMaps setup_inversive_distance_maps(ConformalMesh& mesh) { InversiveDistanceMaps m; @@ -100,8 +112,16 @@ inline InversiveDistanceMaps setup_inversive_distance_maps(ConformalMesh& mesh) return m; } -// Assign sequential DOF indices to all vertices (no gauge pinning here — -// the caller should set one v_idx to −1 before assigning). +/// Assign sequential DOF indices `0..n-1` to every vertex. +/// +/// **Note:** this overload does NOT pin a gauge vertex. The caller +/// is expected to either: +/// 1. set one `m.v_idx[v] = -1` *before* calling this function (then +/// the call is a no-op for that vertex) — OR — +/// 2. flip one assigned index back to `-1` *after* this function. +/// +/// For a closed mesh, exactly one pin is required to remove the +/// global rotational mode. inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& mesh, InversiveDistanceMaps& m) { @@ -110,7 +130,7 @@ inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh& m return idx; } -// Count free DOFs. +/// Count the free DOFs (vertices with `v_idx >= 0`). inline int inversive_distance_dimension(const ConformalMesh& mesh, const InversiveDistanceMaps& m) { @@ -119,18 +139,28 @@ inline int inversive_distance_dimension(const ConformalMesh& mesh, return dim; } -// ── Initialisation from initial mesh geometry ──────────────────────────────── -// -// Two-phase init mirroring "compute_lambda0" for euclidean_functional: -// 1. Choose r_i^(0). Simplest heuristic: r_i = (1/3)·(min adjacent ℓ). -// Other choices (max ℓ, mean ℓ, length-of-shortest-vertex-cycle) are -// possible; the user can override `m.r0[v]` between setup and init. -// 2. Compute I_ij = ( ℓ² − r_i² − r_j² ) / ( 2 r_i r_j ) for each edge. -// -// Note: a valid inversive-distance packing requires I_ij > −1 on every edge, -// and the triangle inequality must hold on every face under the resulting ℓ. -// The choice r_i = ⅓·min(ℓ_e adj v) keeps I_ij safely positive for most -// real meshes. +/// Two-phase initialisation from initial mesh geometry. Mirrors the +/// role of `compute_lambda0_from_mesh` in the Euclidean functional, but +/// adapted to Luo's vertex-based radius parametrisation. +/// +/// **Phase 1.** Pick a positive radius per vertex: +/// \code +/// r_i^(0) = (1/3) · min{ℓ_e : e adjacent to v_i} +/// \endcode +/// This is a heuristic — the user may override `m.r0[v]` for any +/// vertex between `setup_inversive_distance_maps()` and this call. +/// +/// **Phase 2.** Compute the per-edge inversive distance via the +/// Bowers-Stephenson 2004 identity: +/// \code +/// I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j ) +/// \endcode +/// +/// \pre Every edge has positive 3-D length. +/// \pre Radii produced in Phase 1 are positive (degenerate isolated +/// vertices fall back to `r_i = 1`). +/// \post Every `I_e[e] > -1` for a valid packing. The chosen +/// Phase-1 heuristic keeps `I_e > 0` for most real meshes. inline void compute_inversive_distance_init_from_mesh(ConformalMesh& mesh, InversiveDistanceMaps& m) { diff --git a/code/include/mesh_utils.hpp b/code/include/mesh_utils.hpp index bbb3b07..bef7ce7 100644 --- a/code/include/mesh_utils.hpp +++ b/code/include/mesh_utils.hpp @@ -1,14 +1,37 @@ #pragma once +// mesh_utils.hpp +// +// Conversions between CGAL::Surface_mesh and Eigen matrices. Used +// primarily by the viewer / example programs to bridge to libigl, which +// expects (V, F) matrix pairs rather than a halfedge data structure. +// +// All functions are templated on the kernel so the same code works +// with `Simple_cartesian` (production) and with any CGAL +// `Kernel_d::Point_3` (test scaffolding). + #include #include #include namespace mesh_utils { + +/// Copy `mesh` into an Eigen `(V, F)` pair (libigl convention). +/// +/// **Side effect:** `mesh` is triangulated in place via +/// `CGAL::Polygon_mesh_processing::triangulate_faces` so the output +/// `F` is guaranteed to be a 3-column matrix. If `mesh` is already a +/// triangle mesh this is a no-op. +/// +/// \param mesh Input surface mesh. **Modified in place** if any face +/// has more than 3 vertices. +/// \param V Output: `(num_vertices, 3)` matrix of vertex positions. +/// \param F Output: `(num_faces, 3)` matrix of vertex indices per +/// face (rows are individual triangles). template void cgal_to_eigen(CGAL::Surface_mesh& mesh, Eigen::MatrixXd& V, Eigen::MatrixXi& F) { - + CGAL::Polygon_mesh_processing::triangulate_faces(mesh); V.resize(mesh.num_vertices(), 3); @@ -30,13 +53,38 @@ void cgal_to_eigen(CGAL::Surface_mesh& mesh, face_idx++; } } + +/// Quick interactive visualisation via libigl + GLFW. +/// +/// **Requires** `WITH_VIEWER=ON` at CMake time (which is implied by +/// `WITH_CGAL=ON`). Blocks until the viewer window is closed. +/// Not suitable for CI / headless contexts. +/// +/// Typical use: +/// \code{.cpp} +/// Eigen::MatrixXd V; Eigen::MatrixXi F; +/// mesh_utils::cgal_to_eigen(mesh, V, F); +/// mesh_utils::simple_visualize_mesh(V, F); +/// \endcode template void simple_visualize_mesh(Eigen::MatrixXd& V, Eigen::MatrixXi& F) { igl::opengl::glfw::Viewer viewer; viewer.data().set_mesh(V, F); viewer.launch(); } -// Zero-copy map for V (optional) + +/// Zero-copy `Eigen::Map` view of `mesh`'s vertex positions. +/// +/// Returns a row-major `(N, 3)` `Eigen::Map` that aliases the +/// `mesh.points()` storage directly — no allocation, O(1). +/// +/// **Lifetime warning:** the returned `Map` references memory owned by +/// `mesh`. Adding or removing vertices may invalidate the underlying +/// storage; use the `Map` only as long as `mesh` is structurally stable. +/// +/// This is the read-write counterpart to `cgal_to_eigen` for cases +/// where the caller wants to *modify* vertex positions through Eigen +/// (e.g. apply a Möbius transformation) without an intermediate copy. template Eigen::Map> get_vertex_map(CGAL::Surface_mesh& mesh) { diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index 36ee0aa..7ae020b 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -57,8 +57,14 @@ struct SphericalMaps { // Defaults: theta_v = 2π, theta_e = π, lambda0 = 0. // lambda0 = 0 means exp(λ°/2)=1, i.e., l=π — degenerate unless u_i<0. -// For real meshes, set lambda0 from mesh geometry via -// compute_lambda0_from_mesh() below. +/// Attach the five spherical property maps to `mesh` and return their +/// handles. Mirrors `setup_euclidean_maps` but uses the `"sv:"` / +/// `"se:"` prefix so the two functionals can coexist on the same mesh +/// (useful for cross-validation tests). +/// +/// Defaults match the Euclidean defaults except that `lambda0 = 0` here +/// gives `l_e = π` which is degenerate on the unit sphere — always call +/// `compute_lambda0_from_mesh(mesh, m)` next on a real mesh. inline SphericalMaps setup_spherical_maps(ConformalMesh& mesh) { SphericalMaps m; @@ -70,7 +76,8 @@ inline SphericalMaps setup_spherical_maps(ConformalMesh& mesh) return m; } -// Assign DOF indices 0..n-1 for all vertices (only vertex DOFs). +/// Assign sequential DOF indices `0..n-1` to all vertices (no edge DOFs). +/// Caller is expected to pin one gauge vertex with `m.v_idx[v] = -1`. inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m) { int idx = 0; @@ -78,7 +85,9 @@ inline int assign_vertex_dof_indices(ConformalMesh& mesh, SphericalMaps& m) return idx; } -// Assign DOF indices for all vertices AND edges. +/// Assign DOF indices for all vertices AND all edges (vertex-DOFs first, +/// then edge-DOFs). Mirrors `assign_euclidean_all_dof_indices` for the +/// cyclic spherical formulation. inline int assign_all_spherical_dof_indices(ConformalMesh& mesh, SphericalMaps& m) { int idx = 0; @@ -87,7 +96,7 @@ inline int assign_all_spherical_dof_indices(ConformalMesh& mesh, SphericalMaps& return idx; } -// Count variable DOFs. +/// Count the free DOFs (vertices + edges with index `≥ 0`). inline int spherical_dimension(const ConformalMesh& mesh, const SphericalMaps& m) { int dim = 0; @@ -96,9 +105,14 @@ inline int spherical_dimension(const ConformalMesh& mesh, const SphericalMaps& m return dim; } -// Set lambda0 from mesh vertex positions (unit-sphere assumed): -// λ°_e = 2·log(sin(l_e / 2)) where l_e = arccos(p_i · p_j). -// Requires vertices to lie on the unit sphere. +/// Compute `λ°_e` for every edge from the input vertex positions, +/// assuming `mesh` has vertices on the unit sphere. +/// +/// Formula: `λ°_e = 2·log(sin(l_e / 2))` where `l_e = arccos(p_i · p_j)` +/// is the spherical arc length of edge `e`. +/// +/// \pre Every vertex `v` of `mesh` lies on the unit sphere (norm = 1). +/// \pre No edge is degenerate (`p_i ≠ p_j` and `p_i ≠ -p_j`). inline void compute_lambda0_from_mesh(ConformalMesh& mesh, SphericalMaps& m) { for (auto e : mesh.edges()) {