Fix multi-source heat-method geodesic distances (Dirichlet-pinned Poisson solve)

[jf---]

Problem

heat_geodesic_distances / HeatGeodesicSolver (and the two isoline functions built on them) return badly wrong distances whenever more than one source vertex is given — the more sources, the worse. The flagship use case "distance from the boundary" (seed every boundary vertex of an open mesh) reads ~0.5x the true distance near the boundary, with most near-boundary vertices reading below the euclidean straight-line lower bound, which is impossible for a distance. Single-source results are fine, which is why the examples (1–8 sources) never showed it.

Root cause

CGAL's Heat_method_3::Surface_mesh_geodesic_distances_3 recovers the Poisson potential phi only up to one additive constant, and phi is not constant across a multi-vertex source set (the normalized gradient field is not exactly integrable). Its value_at_source_set() then computes

d(i) = min over sources s of |phi(i) - phi(s)|

— the 1-D distance from the value phi(i) to the set of source values. Every vertex whose phi lands inside the source-value spread is folded toward 0, and the far field is shifted by max_s phi(s) instead of a single constant. With one source this degenerates to the correct shift; with many sources it collapses the near field. The fold is many-to-one, so no post-processing of the returned distances can repair it.

Fix

Replace the backend with an own Eigen implementation of the heat method (Crane, Weischedel & Wardetzky 2017) that solves the Poisson step Dirichlet-constrained: source vertices are removed from the unknowns and pinned to 0, so the spread cannot exist in the first place. Same discretization ingredients as before (cotan Laplacian with the same 1e-8 diagonal regularization, barycentric mass, t = mean_edge_length^2), one heat factorization reused across HeatGeodesicSolver.solve() calls.

Contract improvements that come with the pin:

• distances are exactly 0 at every source (previously only approximately, via the fold),
• distances are non-negative,
• correct for arbitrary source sets, including dense bands and solid source regions,
• empty source sets raise ValueError; meshes with a connected component containing no source raise RuntimeError (both previously returned meaningless values),
• the docstring claim that the intrinsic-Delaunay variant was used is dropped (the binding instantiated the class default, Direct).

Validation (reviewers: one command)

pytest tests/test_geodesics.py -q

The new test_heat_geodesic_multisource_boundary_sources uses the new committed fixture data/flat_star_disk_irregular_rim.off: a flat star-shaped disk (z = 0) with irregular, variable-density boundary sampling, all 508 boundary vertices seeded as sources. Because the mesh is flat, the exact geodesic distance to the source set in the near band equals the euclidean distance to the nearest source (verified to ~1e-15 against exact polyhedral MMP geodesics when the fixture was generated), and the euclidean distance is a hard lower bound everywhere — so the test needs no external oracle. It fails on main and passes here:

	band mean d / d_true	band vertices below the euclidean LOWER bound
main (CGAL Heat_method_3)	0.544	83.4 %
this PR	0.979	1.1 % (marginal, sub-edge-length)

Cross-checks against exact polyhedral geodesics (pygeodesic MMP, not a test dependency):

case	main	this PR
12.9k-vertex open slice, 434 boundary sources, near-rim band	0.408	0.966
elephant.off, single source (near-source band)	0.995	0.995
elephant.off, 5 spread sources	0.951	1.002

All 69 existing tests pass unchanged (test_heat_geodesic_distances_multiple_sources already pinned d(source) == 0 and d >= 0 — contracts the fold satisfied by construction and the pin satisfies structurally).

Upstream

The defect is CGAL's, not this binding's; a minimal upstream fix (mean-over-sources shift in value_at_source_set, validated on the same fixture: 0.54 → 1.13, lower-bound violations 83 % → 2.4 %) is filed as CGAL/cgal#9565, with the same fixture and an equivalent C++ test. This PR removes the dependency on the affected class either way and makes the distances exact at the sources, which the mean-shift cannot.

[jf---]

Fixture provenance — data/flat_star_disk_irregular_rim.off is fully synthetic (CC0), generated with the script below (numpy + scipy.spatial.Delaunay, fixed seed, parameters n_rim=500, n_int=2500, amp=0.35, dens_pow=1.5, seed=7), then frozen as a file so the test is deterministic across scipy/Qhull versions. The flat-truth claim (euclid-to-nearest-source == exact geodesic in the 2h..6h band) was verified against exact polyhedral MMP geodesics (pygeodesic, cross-validated bit-identical with tvb-gdist): max disagreement in the band 6.9e-16.

generator

import numpy as np
from scipy.spatial import Delaunay

def build(n_rim=500, n_int=2500, amp=0.35, dens_pow=1.5, seed=7):
    """Flat star disk; rim sampled with variable density (smooth warp of the parameter)."""
    rng = np.random.RandomState(seed)
    u = np.sort(rng.uniform(0, 1, n_rim))
    t = 2 * np.pi * (u + 0.35 / (2 * np.pi * 3) * np.sin(2 * np.pi * 3 * u) * dens_pow)
    r = 1.0 + amp * np.sin(7 * t)
    rim = np.c_[r * np.cos(t), r * np.sin(t)]
    pts = []
    while len(pts) < n_int:
        p = rng.uniform(-1.4, 1.4, size=(n_int * 4, 2))
        ang = np.arctan2(p[:, 1], p[:, 0])
        rad = np.hypot(p[:, 0], p[:, 1])
        keep = rad < (1.0 + amp * np.sin(7 * ang)) * 0.985
        pts.extend(p[keep].tolist())
    P = np.vstack([rim, np.array(pts[:n_int])])
    F = Delaunay(P).simplices.astype(np.int32)
    c = P[F].mean(axis=1)                      # drop triangles outside the concave region
    ang = np.arctan2(c[:, 1], c[:, 0])
    F = F[np.hypot(c[:, 0], c[:, 1]) < 1.0 + amp * np.sin(7 * ang)]
    V = np.c_[P, np.zeros(len(P))]
    used = np.unique(F)                        # drop unused vertices, remap
    remap = -np.ones(len(V), dtype=np.int64)
    remap[used] = np.arange(len(used))
    return np.ascontiguousarray(V[used]), np.ascontiguousarray(remap[F].astype(np.int32))

Why this shape: the defect scales with the spread of the Poisson potential over the source set, and the spread is driven by irregular boundary sampling density (a regular grid barely shows it: ratio 0.96). This mesh maximizes the failure (0.54) while both correct normalizations stay within a few percent.

[petrasvestartas]

Please use input as a compas mesh and internally in the wrapper function convert it to numpy arrays.
The wrappers must be user friendly and very minimal python compas style functions.

[petrasvestartas]

For file imports use:

import pathlib

Define input file path

IFILE = pathlib.Path(file).parent.parent / "data" / "shell_final.json"

[petrasvestartas]

Add in docs example with an image

[petrasvestartas]

For example: