Remeshing: Implement ACVD (uniform, adaptive, and metric-dependent)

[csparker247]

Background

Implement ACVD (Approximated Centroidal Voronoi Diagrams) as a mesh remeshing/coarsening
utility in libcore. Companion to educelab/OpenABF#62, which tracks the same feature for the
flattening pipeline — see Where this lives below.

The three reference papers describe one clustering engine, not three algorithms. What
differs between them is the weight attached to each mesh element and where the output vertex
lands:

Layer	Contribution
Valette & Chassery 2004	The engine. Cluster mesh elements into n connected regions minimizing a weighted compactness energy; build the output as the dual of that clustering. Weight = element area.
Valette, Chassery & Prost 2008	Generalizes the scalar weight to a per-element SPD metric tensor. Yields curvature-adaptive gradation (scalar density) and anisotropic elements (full tensor), plus quadric-based vertex placement and feature/boundary constraints.
Audette et al. 2011	An application of the 2008 machinery: density comes from a spatial field (proximity to a surgical approach corridor) rather than from curvature. Collapses to a caller-supplied callable.

Where this lives

The engine goes here; OpenABF consumes it through a thin adapter (educelab/OpenABF#62).

• ACVD needs no sparse linear algebra. Dense 3×3 is the largest system it ever solves. The
  sole exception is the symmetric eigensolve behind anisotropic metrics — a self-contained
  Jacobi routine (~80 lines). So libcore can host all three papers without gaining a single
  dependency. OpenABF requires Eigen, so hosting the engine there would make Eigen the price
  of admission for any consumer that wants remeshing and nothing else.
• Mesh is a sufficient substrate. ACVD never mutates topology — it reads adjacency,
  writes one integer label per item, and emits a brand-new mesh. No edge collapse, no splicing,
  no incremental connectivity update. A half-edge structure is not required; four flat CSR
  arrays built in a single O(V+F) pass are.
• Dependency direction is OpenABF → libcore (heavy → light), never the reverse.

Proposal

Layered tracks, in dependency order. A + B is the useful milestone and stands alone.

• A — Topology. CSR vertex/vertex, vertex/face and edge adjacency; boundary loops; crease
  detection by dihedral angle; a small index-view type (C++17 has no std::span). Read-only,
  ~200 lines. → types/MeshTopology.hpp
• B — Uniform ACVD (2004). Item weights, seeding, energy minimizer, cluster repair, dual
  construction, trait-aware attribute transfer, quality statistics. →
  utils/Remeshing.hpp, utils/MeshQuality.hpp
• C — Adaptive + placement (2008, scalar). Curvature estimate, gradation weights, quadric
  (QEM) vertex placement, feature/boundary constraints. → utils/Curvature.hpp,
  utils/LinearAlgebra.hpp additions
• D — Anisotropic (2008, full). Curvature tensors, symmetric 3×3 eigensolve, tensor
  accumulators and the metric energy path. Highest risk, least load-bearing for our meshes —
  defer until something needs anisotropy.
• E — Approach-guided (2011). Field-callable overload, ROI distance helpers, size↔density
  calibration, quality report.
• F — Optional. Subdivision pre-pass for inputs too coarse for the requested cluster count;
  spatial index; UV/attribute resampling onto the new mesh.
• OpenABF wrapper. Build the CSR from HalfEdgeMesh (vertices already carry idx, so it's
  one pass over faces()), call the engine, return through
  insert_vertices/insert_faces. Tracked in [Feature] Implement ACVD OpenABF#62.

Proposed API surface

// include/educelab/core/utils/Remeshing.hpp
enum class ClusterDomain { Vertices, Faces };
enum class SitePlacement { Centroid, Projected, Quadric, QuadricRegularized };

template <typename T = float>
struct RemeshOptions {
    std::size_t   clusters{0};        // == desired output vertex count
    ClusterDomain domain{ClusterDomain::Vertices};
    SitePlacement placement{SitePlacement::QuadricRegularized};
    T             gradation{0};       // 0 = uniform (2004); >0 = curvature-adaptive (2008)
    T             featureAngle{to_radians<T>(60)};
    bool          preserveBoundary{true};
    std::size_t   maxSweeps{200};
    T             tolerance{1e-6};
    std::uint64_t seed{0};            // explicit: educelab::random() is not seedable
    Signal<std::size_t, double>* progress{nullptr};  // (sweep, energy)
};

template <typename T, std::size_t Dims, typename VTraits>
struct RemeshResult {
    Mesh<T, Dims, VTraits>   mesh;    // the remeshed output
    std::vector<std::size_t> labels;  // cluster id per input item
    double                   energy{};
    std::size_t              sweeps{};
    std::size_t              repairs{};
};

// 2004 uniform, and 2008 curvature-adaptive via opts.gradation
template <typename T, std::size_t Dims, typename VTraits>
[[nodiscard]] auto acvd_remesh(const Mesh<T, Dims, VTraits>& mesh,
                               const RemeshOptions<T>& opts)
    -> RemeshResult<T, Dims, VTraits>;

// 2008 generic / 2011 approach-guided: caller supplies the field.
//   WeightFn: T(const Vec<T,Dims>&)          -> scalar density
//   WeightFn: Mat<3,3,T>(const Vec<T,Dims>&) -> SPD metric tensor
template <typename T, std::size_t Dims, typename VTraits, typename WeightFn>
[[nodiscard]] auto acvd_remesh(const Mesh<T, Dims, VTraits>&,
                               const RemeshOptions<T>&, WeightFn&&)
    -> RemeshResult<T, Dims, VTraits>;

// Asymptotic CVD size relation: rho = h^-4
template <typename T>
[[nodiscard]] constexpr auto density_for_edge_length(T h) -> T;

Pipeline

1. Topology & preconditions — build CSR adjacency; filter zero-area faces; assert
   clusters <= item count, warn below ~8 items per cluster (below that the dual degenerates
   and the input needs subdividing first).
2. Item weights — lumped area (one third of the summed area of incident triangles), times
   curvature^gradation for adaptive or a caller-supplied spatial density for approach-guided.
3. Per-item metric tensors — anisotropic path only.
4. Seeding — n density-weighted seeds, then simultaneous multi-source BFS so every
   initial cluster is connected by construction. Every connected component must get ≥1 seed or
   it vanishes from the output.
5. Energy minimization — sweep the boundary-item queue applying the O(1) test below; on an
   accepted move, re-push the item's neighbours.
6. Cluster repair — reseed empty clusters by splitting the highest-energy cluster; split
   clusters that have become disconnected; re-minimize.
7. Site placement — centroid, surface-projected centroid, or quadric minimization
   regularized toward the centroid. Boundary clusters take line quadrics from their boundary
   edges; crease corners get pinned.
8. Dual construction & validity — one output triangle per input face carrying three
   distinct labels (winding carries over, so orientation is free), then repair what the dual
   can't express: clusters with <3 neighbours, duplicate triangles, edges with >2 incident
   faces.
9. Attribute transfer & report — average normals/colors per cluster behind if constexpr
   on the existing traits::has_normal / has_color. Return labels, final energy, quality
   statistics.

Implementation notes

Derived from the papers and verified numerically; recording here so they don't have to be
rediscovered.

The reassignment test. With cluster mass m and centroid γ, moving item s (weight ρ,
position p) from cluster a to b lowers the energy iff

m_a‖γ_a − p‖² / (m_a − ρ_s)  >  m_b‖γ_b − p‖² / (m_b + ρ_s)

Exact, translation-invariant, O(1). Verified against brute-force energy recomputation over 400
random reassignments: max relative error 1.8e-13 (double). Guard with m_a − ρ_s > 0 (never
empty a cluster — that silently drops an output vertex) and a strict improvement margin (an
epsilon-free comparison cycles between equal-energy states).

Do not evaluate the energy as Q − Σ‖c‖²/m. It is algebraically correct but both terms are
origin-dependent and nearly equal. Measured in float32 on 5,000 items offset to
(8000, 4000, 12000) with 200 clusters: true energy 57,442.31, subtractive form returned
−262,144.00 — a negative energy, which will send a sign-driven minimizer wandering. Mean-center
positions once up front and accumulate m/c in double regardless of the mesh's T.

Density ↔ target edge length. Asymptotic optimal site density on a surface goes as ρ^(1/2),
and site density is 1/h², so h ∝ ρ^(-1/4) ⇔ ρ = h^-4. Treat as a calibration and verify the
achieved size ratio numerically — a two-region size field is an easy test and the natural place
to catch an off-by-one in the exponent.

Metric generalization. A_i = ΣM_s, b_i = ΣM_s x_s, k_i = Σx_sᵀM_s x_s, giving
γ_i = A_i⁻¹b_i and E_i = k_i − b_iᵀA_i⁻¹b_i. Ten doubles per cluster instead of four;
M = ρI recovers the scalar case exactly, so one implementation covers both — but specialize
the scalar path, since the general one costs two 3×3 SPD solves per candidate move instead of
two squared distances. A is a sum of SPD tensors, so Cholesky is safe.

Gotchas in existing code

• normalize() has no zero guard, so Mesh::face_normal returns NaN for a zero-area triangle,
  and those NaNs propagate into quadrics and curvature tensors. Filter degenerate faces in
  stage 1 and treat their weight as zero.
• interior_angle() does not clamp its acos argument — dot products can land marginally
  outside [-1, 1] in float. vertex_normal is built on it.
• norm() accumulates in the element type. Fine for one vector, wrong for cluster
  accumulators.
• Mesh::insert_vertex clears the whole face-normal cache (O(F) per insertion), so interleaving
  vertex and face insertion during output construction is quadratic. Insert all vertices first.
• Mesh::vertex_faces returns a vector<vector> — ~48 B/vertex of overhead and poor locality
  at scale, hence the CSR structure in track A.
• docs/citations.bib is wired into Doxygen (CITE_BIB_FILES) and empty. Add the three papers
  there and reference them with @cite from the new headers.

Acceptance criteria

Track A

• MeshTopology provides vertex→vertex, vertex→face and face→face adjacency in CSR form,
  an edge table with incident faces, boundary-loop extraction, and crease flags.
• Built in a single O(V+F) pass; no allocation per adjacency query.
• Tests in tests/src/TestMeshTopology.cpp, registered in tests/CMakeLists.txt.

Track B

• acvd_remesh produces exactly opts.clusters output vertices for valid inputs.
• Incremental cluster accumulators match a from-scratch recompute after N random
  reassignments.
• The O(1) move test matches a full energy recomputation (regression test for the formula
  above).
• Energy is non-increasing across sweeps; convergence within maxSweeps on all fixtures.
• Identity case: clusters == item count reproduces the input mesh.
• Uniformity: on a regular grid patch and a subdivided icosphere, edge-length coefficient of
  variation is below a documented bound; sphere output vertices lie on the sphere within
  tolerance.
• Topology: every output edge has ≤2 incident faces, orientation is consistent, Euler
  characteristic is preserved for closed genus-0 input, boundary loop count is preserved for
  an open patch.
• Determinism: same seed → identical output, in-process and across Debug/Release.
• Degenerate inputs have defined, tested behavior: empty mesh, clusters == 0,
  clusters > item count, zero-area faces, duplicate vertices, disconnected components,
  non-manifold edges.
• Per-vertex normals and colors transfer when the vertex traits carry them, and are absent
  when they don't.

Tracks C–E

• Gradation > 0 concentrates vertices in high-curvature regions (measured, not eyeballed).
• Quadric placement preserves a crease better than centroid placement on a fixture with a
  known sharp edge.
• Open-mesh boundaries do not shrink; boundary vertex positions stay on the input boundary
  curve within tolerance.
• A two-region density field achieves the predicted mean edge-length ratio, validating the
  ρ = h^-4 exponent.
• Anisotropic metrics produce measurably elongated elements in the prescribed direction.

Docs & build

• Three papers added to docs/citations.bib, referenced via @cite from the new headers.
• New headers listed in public_hdrs in the root CMakeLists.txt and in the header-only
  list in README.md.
• Doxygen builds cleanly for the new public API.

Open questions

To settle against the paper PDFs before the minimizer hardens around a choice:

1. Vertices or faces as the clustered items? This plan assumes vertex clustering, which makes
   the requested cluster count equal the output vertex count and yields triangles directly by
   duality. Face clustering changes the dual rule and the boundary handling.
2. Boundary treatment for open meshes. A separate 1D clustering along each boundary loop
   (so the output boundary is a polyline through boundary clusters) is the natural approach, and
   matters more for our meshes than for the papers' test cases — most EduceLab meshes are open
   patches, where a shrinking or ragged boundary is disqualifying.
3. The exact gradation formula — the curvature→density map and the clamping that keeps it
   stable in flat regions.
4. Which curvature estimator the 2008 paper uses; standard choices differ noticeably on
   noisy meshes.
5. How much topology repair is enough? The papers acknowledge that clusters can become
   disconnected and that the dual can be invalid; a complete repair strategy is likely sketched
   in the paper and load-bearing in practice. Budget for iteration.
6. UV maps. Remeshing invalidates a UVMap outright. Dropping it is honest; resampling
   requires nearest-face queries and a spatial index. Decide the contract explicitly rather than
   leaving callers to discover it.

Items 1 and 2 are the decisions the rest of the code shapes itself around — cheap to settle now,
expensive to revisit.

Out of scope

• Deviation-metric-driven iterative coarsening (coarsen until surface deviation exceeds a
  tolerance). Scoped and deliberately deferred: it needs a spatial index plus two-sided sampled
  Hausdorff measurement, and the ACVD energy is a compactness metric, not a deviation metric, so
  error control has to be layered on rather than read off.
• QEM edge-collapse decimation. Better matched to "fewest vertices under a hard error bound",
  but a different algorithm, and the one case where a mutable half-edge structure would genuinely
  earn its keep.
• Parallelism. The boundary sweep is inherently sequential; single-threaded is the target.

Caveats

There is a reference implementation of this algorithm available online. This must be a pure
clean-room implementation from the reference papers — we MUST NOT examine that reference
library. Everything in this issue was derived from the papers' described methods and verified
numerically, with no reference implementation consulted.

References

Valette, Sébastien, and Jean-Marc Chassery. "Approximated centroidal Voronoi diagrams for
uniform polygonal mesh coarsening." Computer Graphics Forum 23, no. 3 (2004): 381–389.
[PDF]

Valette, Sébastien, Jean-Marc Chassery, and Rémy Prost. "Generic remeshing of 3D triangular
meshes with metric-dependent discrete Voronoi diagrams." IEEE Transactions on Visualization and
Computer Graphics 14, no. 2 (2008): 369–381. doi:10.1109/TVCG.2007.70430
[PDF]

Audette, Michel, Denis Rivière, Matthew Ewend, Andinet Enquobahrie, and Sébastien Valette.
"Approach-guided controlled resolution brain meshing for FE-based interactive neurosurgery
simulation." Workshop on Mesh Processing in Medical Image Analysis, MICCAI 2011: 176–186.

• Companion issue: [Feature] Implement ACVD OpenABF#62 (ACVD for the flattening pipeline; will consume this
  engine via an adapter)