"""Surface algorithms on a ``TriMesh`` as free functions. All indices 0-based."""
from __future__ import annotations
from collections import deque
from collections.abc import Callable, Iterator, Sequence
from typing import Any
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from scipy.spatial import cKDTree
from .._typing import BoolArray, FloatArray, IntArray, Point, PointArray
from ..core import conform
from ..linemesh import LineMesh
from ..linemesh.assemble import loft as line_loft
from ..linemesh.assemble import loft_fn as line_loft_fn
from ..linemesh.assemble import merge as line_merge
from ..linemesh.morph import reverse as line_reverse
from ..linemesh.shape import arclength_fractions
from .trimesh import TriMesh
# -- Laplace ------------------------------------------------------------
[docs]
def cotan_laplacian(surface: TriMesh) -> sp.csr_matrix:
"""Linear (P1) cotangent-weight stiffness matrix (recomputed each call)."""
xyz, tri = surface.points, surface.tris
nv, nt = xyz.shape[0], tri.shape[0]
cot_clamp = 1e3
p = xyz[tri] # (nt,3,3)
# Build triples in the original loop order so duplicate summation is bit-identical.
I = np.empty((nt, 3, 4))
J = np.empty((nt, 3, 4))
S = np.empty((nt, 3, 4))
valid = np.zeros((nt, 3), dtype=bool)
for c in range(3):
a, b = (c + 1) % 3, (c + 2) % 3
e1 = p[:, a, :] - p[:, c, :]
e2 = p[:, b, :] - p[:, c, :]
area2 = np.linalg.norm(np.cross(e1, e2), axis=1)
ok = area2 >= 1e-14
cotc = np.divide(np.sum(e1 * e2, axis=1), area2,
out=np.zeros(nt), where=ok)
cotc = np.clip(cotc, -cot_clamp, cot_clamp)
wgt = 0.5 * cotc
ia, ib = tri[:, a], tri[:, b]
I[:, c, 0], I[:, c, 1], I[:, c, 2], I[:, c, 3] = ia, ib, ia, ib
J[:, c, 0], J[:, c, 1], J[:, c, 2], J[:, c, 3] = ib, ia, ia, ib
S[:, c, 0], S[:, c, 1], S[:, c, 2], S[:, c, 3] = -wgt, -wgt, wgt, wgt
valid[:, c] = ok
mask = np.repeat(valid[:, :, None], 4, axis=2).ravel()
If, Jf, Sf = I.ravel()[mask], J.ravel()[mask], S.ravel()[mask]
return sp.coo_matrix((Sf, (If, Jf)), shape=(nv, nv)).tocsr()
[docs]
def solve_dirichlet(surface: TriMesh, dpoints: IntArray, dvals: FloatArray) -> FloatArray:
"""Solve the Laplace system with Dirichlet values at ``dpoints``; returns an ``(nv,)`` field."""
L = cotan_laplacian(surface)
nv = surface.n_points
dpoints = np.asarray(dpoints, dtype=np.int64).ravel()
dvals = np.asarray(dvals, dtype=float).ravel()
u = np.zeros(nv)
u[dpoints] = dvals
is_d: BoolArray = np.zeros(nv, dtype=bool)
is_d[dpoints] = True
fpoints = np.flatnonzero(~is_d)
Lc = L.tocsr()
rhs = -Lc[fpoints, :][:, dpoints] @ u[dpoints]
A = Lc[fpoints, :][:, fpoints].tocsc()
u[fpoints] = spla.spsolve(A, np.asarray(rhs).ravel())
return u
# -- boundary loops -----------------------------------------------------
def _boundary_edges(surface: TriMesh) -> IntArray:
tri = surface.tris
E = np.vstack([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]])
E = np.sort(E, axis=1)
# ``conform.unique_rows``, not ``np.unique(axis=0)``: the latter views each row as a
# void scalar and argsorts that, which measured 137 ms on a 48k-vertex surface and
# was the entire cost of ``boundary_loops`` -- the walk underneath it is free.
# Packing the two ids into one int64 gives the same lexicographic order, so the
# result is identical rather than merely equivalent.
Eu, _ic, cnt = conform.unique_rows(E, return_counts=True)
return np.asarray(Eu[cnt == 1, :], dtype=np.int64)
[docs]
def boundary_loops(surface: TriMesh) -> list[IntArray]:
"""Group open-boundary edges into connected loops; returns a list of vertex-index arrays."""
bnd_edges = _boundary_edges(surface)
bnd_verts: IntArray = np.unique(bnd_edges.ravel())
if bnd_verts.size == 0:
return []
# Index the adjacency by *boundary* vertex, not by every vertex on the surface.
# Allocating one list per surface point cost 48k empty lists to hold the 876 useful
# ones on a scan-sized patch -- the walk itself was never the expensive part, so this
# keeps the BFS, and with it the discovery order the loop ordering depends on,
# exactly as it was.
local: IntArray = np.searchsorted(bnd_verts, bnd_edges)
adj: list[list[int]] = [[] for _ in range(bnd_verts.size)]
for i, j in local:
adj[i].append(j)
adj[j].append(i)
visited: BoolArray = np.zeros(bnd_verts.size, dtype=bool)
loops = []
for v0 in range(bnd_verts.size):
if visited[v0]:
continue
comp = [v0]
visited[v0] = True
queue = deque([v0])
while queue:
point = queue.popleft()
for w in adj[point]:
if not visited[w]:
visited[w] = True
comp.append(w)
queue.append(w)
loops.append(bnd_verts[np.array(comp, dtype=np.int64)])
return loops
[docs]
def order_boundary_loop(surface: TriMesh, lv: IntArray) -> IntArray:
"""Return the vertices of one boundary loop ``lv`` in cyclic order."""
nv = surface.n_points
be = _boundary_edges(surface)
inlv: BoolArray = np.zeros(nv, dtype=bool)
inlv[lv] = True
be = be[inlv[be[:, 0]] & inlv[be[:, 1]], :]
adj: dict[int, list[int]] = {}
for a, b in be:
adj.setdefault(a, []).append(b)
adj.setdefault(b, []).append(a)
start = int(lv[0])
ordv = [start]
prev = -1
cur = start
while True:
nb = [x for x in adj[cur] if x != prev]
if not nb or nb[0] == start:
break
prev = cur
cur = nb[0]
ordv.append(cur)
return np.array(ordv, dtype=np.int64)
# -- trifurcation splitting ----------------------------------------------
#
# A trifurcation (three boundary loops on one open surface) is cut into three legs by
# solving three Laplace fields and reading their sign pattern -- shared, verbatim, by
# every mesher in this toolkit that starts from a trifurcated wall (``examples/carotid.py``,
# ``examples/femoral.py``). Identifying *which* physical loop is leg 1/2/3 is left to the
# caller: that identification is a property of one specific junction's layout, not of the
# splitting algorithm below it -- see each example's own ``order_openings``.
[docs]
def conduction_field(
surf: TriMesh, gloops: Sequence[IntArray], neumann: int, dvals: Sequence[float],
) -> FloatArray:
"""Laplace with Neumann on loop ``neumann``, Dirichlet ``dvals`` on the other two,
shifted to zero mean on the free loop."""
dpoints, dv = [], []
for k in range(3):
if k == neumann:
continue
g = np.asarray(gloops[k]).ravel()
dpoints.append(g)
dv.append(dvals[k] * np.ones(g.size))
u = solve_dirichlet(surf, np.concatenate(dpoints), np.concatenate(dv))
return u - np.mean(u[gloops[neumann]])
[docs]
def seam_fields(surf: TriMesh, gloops: Sequence[IntArray]) -> FloatArray:
"""The three conduction seam fields ``U`` ``(nv, 3)``: column ``k`` is 0 on opening
``k``, 1 on the next opening round, and Neumann (free) on the third -- so its sign
is what :func:`cut_surface_from_fields` reads to place seam ``k``."""
nan = np.nan
dvals = [[nan, 0, 1], [1, nan, 0], [0, 1, nan]]
U = np.zeros((surf.n_points, 3))
for k in range(3):
U[:, k] = conduction_field(surf, gloops, k, dvals[k])
return U
[docs]
def leg_label(F: FloatArray) -> IntArray:
"""Which of the three legs (1, 2 or 3) each row of seam-field values ``F``
``(n, 3)`` falls in, by the sign pattern :func:`seam_fields` gives each leg."""
a, b, c = F[:, 0], F[:, 1], F[:, 2]
lab: IntArray = np.zeros(a.shape[0], dtype=np.int64)
lab[(b > 0) & (c < 0)] = 1
lab[(a < 0) & (c > 0)] = 2
lab[(a > 0) & (b < 0)] = 3
return lab
[docs]
def cut_surface_from_fields(
surf: TriMesh, F: FloatArray,
) -> tuple[PointArray, list[IntArray]]:
"""Cut ``surf`` into three legs defined by seam fields ``F`` ``(nv, 3)``,
retriangulating triangles a seam crosses. Returns ``(V, faces)``: ``V`` is the
surface's own points plus every new seam-crossing and triple point, and ``faces``
is a 3-element list of ``(T,3)`` arrays -- one leg's triangles each -- both
indexing into ``V``."""
xyz = surf.points
tri = surf.tris
nv = xyz.shape[0]
V = [xyz[i, :].copy() for i in range(nv)]
faces: dict[int, list[list[int]]] = {1: [], 2: [], 3: []}
lab = leg_label(F)
ecache: dict[tuple[int, int, int], int] = {}
def edge_pt(vi: int, vj: int, fi: int) -> int:
key = (min(vi, vj), max(vi, vj), fi)
if key in ecache:
return ecache[key]
fi_i = F[vi, fi]
fi_j = F[vj, fi]
t = fi_i / (fi_i - fi_j)
V.append(xyz[vi, :] + t * (xyz[vj, :] - xyz[vi, :]))
idv = len(V) - 1
ecache[key] = idv
return idv
def triple_pt(v: IntArray) -> int:
Aeq = np.array([[F[v[0], 0], F[v[1], 0], F[v[2], 0]],
[F[v[0], 1], F[v[1], 1], F[v[2], 1]],
[1.0, 1.0, 1.0]])
lam = np.linalg.solve(Aeq, np.array([0.0, 0.0, 1.0]))
V.append(lam @ xyz[v, :])
idv = len(V) - 1
return idv
for e in range(tri.shape[0]):
v = tri[e, :]
l = lab[v]
ul = np.unique(l)
if ul.size == 1:
faces[int(ul[0])].append([v[0], v[1], v[2]])
elif ul.size == 2:
cnt = np.array([np.sum(l == u) for u in ul])
lone = int(ul[cnt == 1][0])
pair = int(ul[cnt == 2][0])
p = int(v[l == lone][0])
qr = v[l == pair]
q = int(qr[0])
r = int(qr[1])
fi = (6 - lone - pair) - 1
e1 = edge_pt(p, q, fi)
e2 = edge_pt(p, r, fi)
faces[lone].append([p, e1, e2])
faces[pair].append([q, r, e2])
faces[pair].append([q, e2, e1])
else:
T = triple_pt(v)
l0, l1, l2 = int(l[0]), int(l[1]), int(l[2])
m12 = edge_pt(int(v[0]), int(v[1]), (6 - l0 - l1) - 1)
m23 = edge_pt(int(v[1]), int(v[2]), (6 - l1 - l2) - 1)
m31 = edge_pt(int(v[2]), int(v[0]), (6 - l2 - l0) - 1)
faces[l0].append([int(v[0]), m12, T])
faces[l0].append([int(v[0]), T, m31])
faces[l1].append([int(v[1]), m23, T])
faces[l1].append([int(v[1]), T, m12])
faces[l2].append([int(v[2]), m31, T])
faces[l2].append([int(v[2]), T, m23])
Varr: PointArray = np.array(V, dtype=float)
faces_list = [np.array(faces[k], dtype=np.int64).reshape(-1, 3) for k in (1, 2, 3)]
return Varr, faces_list
[docs]
def leg_field(
V: PointArray, faces: list[IntArray], leg: int, gloops: Sequence[IntArray],
) -> tuple[TriMesh, FloatArray, IntArray, IntArray, IntArray]:
"""One leg as a sub-mesh with Laplace solved (0 on opening, 1 on seam).
Returns ``(sub, us, opening, seam, vids)``."""
sub, vids = TriMesh.from_faces(V, faces[leg])
sloops = [c for c in boundary_loops(sub) if c.size >= 3]
gset = set(int(x) for x in gloops[leg])
opencnt = np.array([np.sum([1 for x in vids[c] if int(x) in gset]) for c in sloops])
oi = int(np.argmax(opencnt))
rest = [i for i in range(len(sloops)) if i != oi]
si = rest[int(np.argmax([sloops[i].size for i in rest]))]
opening = sloops[oi]
seam = sloops[si]
us = solve_dirichlet(
sub,
np.concatenate([opening, seam]),
np.concatenate([np.zeros(opening.size), np.ones(seam.size)]))
return sub, us, opening, seam, vids
def _split_two_arcs(
ordv: IntArray, iA1: int, iA2: int, segX: IntArray, segY: IntArray,
) -> tuple[IntArray, IntArray]:
ordv = np.asarray(ordv).ravel()
k1 = int(np.flatnonzero(ordv == iA1)[0])
ordv = np.roll(ordv, -k1)
k2 = int(np.flatnonzero(ordv == iA2)[0])
a = ordv[0:k2 + 1]
b = np.concatenate([ordv[k2:], ordv[:1]])[::-1]
interiorA = a[1:-1]
segXset = set(int(x) for x in segX)
if all(int(x) in segXset for x in interiorA):
return a, b
return b, a
[docs]
def arc_curve(
V: PointArray, arcverts: IntArray, iA1: int, n: int, keep: float = 0.5,
) -> Callable[[FloatArray], PointArray]:
"""Analytic ``A1 -> A2`` seam arc: the scanned polyline's deviation from its own
chord, refit as a truncated **sine** series in the normalized arc-length parameter
``s``, and returned as a callable ``p(s)`` on ``[0, 1]``.
This is :func:`fourier_ring`'s open-curve sibling, and the basis is chosen for the
one property the seam needs: every ``sin(k*pi*s)`` vanishes at both ends, so ``A1``
and ``A2`` -- the triple points where all three arcs meet -- stay **bit-exact**
however many modes are kept. That is what lets the three legs keep welding: each
arc is fitted once, globally, and the *same* mesh is handed to both legs that share
it, so no leg ever sees a privately-refitted seam (which is why the seam ring
could not simply go through :func:`fourier_wall` -- refitting per leg would mix the
two arcs that leg happens to see).
Dropping the upper modes removes the STL facet noise for the same reason it does on
a station ring, and ``keep`` is applied to the modes the *mesh* can resolve
(``n`` elements -> ``n-1`` interior nodes -> modes ``1..n-1``), mirroring
:func:`fourier_ring`'s ``keep * (M // 2 + 1)``."""
arcverts = np.asarray(arcverts).ravel()
if arcverts[0] != iA1:
arcverts = arcverts[::-1]
N = 4 * n # dense uniform-arc-length fit samples
s = np.linspace(0.0, 1.0, N + 1)
Q = resample_polyline(V[arcverts, :], s)
a1, a2 = Q[0, :], Q[-1, :]
dev = Q - (a1 + s[:, None] * (a2 - a1)) # vanishes at both ends by construction
K = max(1, int(keep * (n - 1)))
kk = np.arange(1, K + 1)
# discrete sine transform (type I) of the interior samples, truncated to K modes
S = np.sin(np.pi * np.outer(np.arange(1, N), kk) / N) # (N-1, K)
B = (2.0 / N) * (S.T @ dev[1:N, :]) # (K,3) sine coefficients
def p(t: FloatArray) -> PointArray:
t = np.atleast_1d(np.asarray(t, dtype=float))
return (a1 + t[:, None] * (a2 - a1)
+ np.sin(np.pi * np.outer(t, kk)) @ B)
return p
[docs]
def arc_mesh(
p: Callable[[FloatArray], PointArray], n: int, order: int, element_tag: str,
) -> LineMesh:
"""Mesh an :func:`arc_curve` into ``n`` high-order elements of the given ``order``,
evenly spaced by arc length. ``LineMesh.loft_fn`` evaluates ``p`` at every node --
corners *and* the private GLL interiors -- so no node lands on a chord."""
return line_loft_fn(p, arclength_fractions(p, n), order=order, element_tags=element_tag)
def _ring(p: LineMesh, q: LineMesh) -> LineMesh:
"""Close two shared-endpoint ``A1 -> A2`` arcs into one loop by welding them at
``A1`` and ``A2`` (:func:`linemesh.merge <nekmeshpy.linemesh.assemble.merge>`); ``q`` is
reversed so the traversal runs ``A1 -> A2`` down ``p`` then ``A2 -> A1`` back up
``q`` without crossing.
``reverse`` carries ``q``'s high-order nodes with it; re-lofting its points
would straight-subdivide them and lose the curve at ``order > 1``."""
return line_merge([p, line_reverse(q)])
[docs]
def seam_rings(
V: PointArray, faces: list[IntArray], gloops: Sequence[IntArray], n_half: int, *,
order: int, keep: float = 0.5, element_tag: str = "wall",
) -> tuple[list[LineMesh], Point, Point, LineMesh]:
"""Build the three conformal seam rings + spine.
Returns ``(rings, A1, A2, spine)``."""
segV: list[Any] = [None, None, None]
ordV: list[Any] = [None, None, None]
for leg in range(3):
sub, _, _, seam, vids = leg_field(V, faces, leg, gloops)
segV[leg] = vids[seam]
ordV[leg] = vids[order_boundary_loop(sub, seam)]
common = np.intersect1d(np.intersect1d(segV[0], segV[1]), segV[2])
Pc = V[common, :]
diff = Pc[:, None, :] - Pc[None, :, :]
Dc = np.sum(diff ** 2, axis=2)
# MATLAB max(Dc(:)) scans column-major; replicate its tie-break.
mx = int(np.argmax(Dc.ravel(order="F")))
ia, ib = np.unravel_index(mx, Dc.shape, order="F")
iA1 = int(common[ia])
iA2 = int(common[ib])
A1: Point = V[iA1, :]
A2: Point = V[iA2, :]
arcAB, arcAC = _split_two_arcs(ordV[0], iA1, iA2, segV[1], segV[2])
bc1, bc2 = _split_two_arcs(ordV[1], iA1, iA2, segV[2], segV[0])
segAset = set(int(x) for x in segV[0])
arcBC = bc1 if all(int(x) in segAset for x in bc2) else bc2
# Each shared arc is refit and meshed exactly once here, so the two legs that
# share it are handed bit-identical geometry and the blocks still weld.
abP = arc_mesh(arc_curve(V, arcAB, iA1, n_half, keep), n_half, order, element_tag)
acP = arc_mesh(arc_curve(V, arcAC, iA1, n_half, keep), n_half, order, element_tag)
bcP = arc_mesh(arc_curve(V, arcBC, iA1, n_half, keep), n_half, order, element_tag)
rings = [_ring(abP, acP), _ring(abP, bcP), _ring(acP, bcP)]
spine = line_loft((abP.points + acP.points + bcP.points) / 3.0, order=order)
return rings, A1, A2, spine
# -- analytic wall from a scanned ring -----------------------------------
[docs]
def fourier_ring(P: PointArray, keep: float = 0.5) -> Callable[[FloatArray], PointArray]:
"""Refit a closed scanned ring as a truncated Fourier series ``p(t)``, periodic
on ``[0, 2*pi)`` with sample ``j`` at ``t = 2*pi*j/M``.
``x``, ``y`` and ``z`` are transformed independently (``np.fft.rfft`` against the
uniform parameter) and only the lowest ``keep`` fraction of the modes is retained;
the rest -- facet-scale noise from the STL projection, which a high-order wall would
otherwise resolve faithfully -- is dropped. The result is a genuine closed-form
curve, so ``LineMesh.loft_fn`` can place every GLL node *on* it instead of on the
chords between the samples (see the straight-GLL-subdivision trap in CLAUDE.md).
"""
P = np.asarray(P, dtype=float)
M = P.shape[0]
nk = max(2, int(keep * (M // 2 + 1))) # modes 0..nk-1; never the Nyquist mode
C = np.fft.rfft(P, axis=0)[:nk, :] / M
C[1:, :] *= 2.0 # each retained mode is a conjugate pair
k = np.arange(nk)
def p(t: FloatArray) -> PointArray:
t = np.atleast_1d(np.asarray(t, dtype=float))
return np.real(np.exp(1j * t[:, None] * k[None, :]) @ C)
return p
[docs]
def fourier_wall(P: PointArray, order: int, element_tag: str, keep: float = 0.5) -> LineMesh:
"""Closed high-order wall ``LineMesh`` through :func:`fourier_ring`, sampled at the
same ``M`` parameters as ``P`` so index 0 / ``M//2`` still land on the spine rails.
``LineMesh.loft_fn`` only ever makes an open chain, so weld the two ends into a loop
the way the toolkit expects (:func:`linemesh.merge <nekmeshpy.linemesh.assemble.merge>`)."""
M = np.asarray(P).shape[0]
return line_merge([line_loft_fn(
fourier_ring(P, keep=keep), np.linspace(0.0, 2.0 * np.pi, M + 1),
order=order, element_tags=element_tag)])
# -- isocontours --------------------------------------------------------
def _chain_segments(segs: FloatArray) -> LineMesh | None:
"""Chain unordered 3-D segments ``(S,6)`` into a single closed ordered loop (the
largest connected component), or ``None`` if they form no loop."""
if segs is None or len(segs) == 0:
return None
segs = np.asarray(segs, dtype=float)
ns = segs.shape[0]
pts_raw = np.vstack([segs[:, 0:3], segs[:, 3:6]])
# Weld coincident endpoints by a scale-relative *radius*, not a grid. Snapping to
# a lattice missed any two endpoints that straddled a bin boundary however close
# they were, and a single missed weld here does not raise -- it breaks the chain, so
# the walk below returns a shorter loop, or ``None``.
scl = conform.bbox_scale(pts_raw)
tol = 1e-6 * scl if scl > 0 else 1.0
_, ic = np.unique(conform.coincident_clusters(pts_raw, tol), return_inverse=True)
ic = ic.ravel()
npts = int(ic.max()) + 1
coord = np.zeros((npts, 3))
cnt = np.zeros(npts)
for i in range(2 * ns):
coord[ic[i], :] += pts_raw[i, :]
cnt[ic[i]] += 1
coord = coord / cnt[:, None]
# node -> incident neighbours and each node pair -> its segment indices,
# so the walk finds the next segment by incidence.
seg_ids = np.column_stack([ic[:ns], ic[ns:]])
adj: list[list[int]] = [[] for _ in range(npts)]
pair_segs: dict[tuple[int, int], list[int]] = {}
for s in range(ns):
i, j = int(seg_ids[s, 0]), int(seg_ids[s, 1])
adj[i].append(j)
adj[j].append(i)
pair_segs.setdefault((i, j) if i <= j else (j, i), []).append(s)
visited_seg = np.zeros(ns, dtype=bool)
loops = []
for s in range(ns):
if visited_seg[s]:
continue
start = int(seg_ids[s, 0])
cur = int(seg_ids[s, 1])
visited_seg[s] = True
order = [start, cur]
while cur != start:
found = False
for nb in adj[cur]:
for s2 in pair_segs[(cur, nb) if cur <= nb else (nb, cur)]:
if visited_seg[s2]:
continue
visited_seg[s2] = True
cur = nb
order.append(cur)
found = True
break
if found:
break
if not found:
break
loops.append(line_loft(coord[order, :], loop=True))
if not loops:
return None
return max(loops, key=lambda ring: ring.n_points)
# -- polyline / ring conformalization -----------------------------------
def _lerp_along(P: PointArray, arclen: FloatArray, targets: FloatArray) -> PointArray:
"""Piecewise-linear resample of polyline ``P`` at each query arc length in
``targets`` (``arclen`` = ``P``'s cumulative arc length); callers pass already
clamped ``targets``. Returns ``(len(targets), 3)``."""
K = P.shape[0]
# ``searchsorted(side="right") - 1`` is the last index with ``arclen[idx] <= s`` --
# exactly what the old per-target ``flatnonzero(arclen <= s)[-1]`` scan returned,
# ties included, but O(log K) instead of O(K) and without the Python loop.
idx: IntArray = np.clip(
np.searchsorted(arclen, targets, side="right") - 1, 0, K - 2)
span: FloatArray = arclen[idx + 1] - arclen[idx]
t: FloatArray = np.where(span > 0.0,
(targets - arclen[idx]) / np.where(span > 0.0, span, 1.0),
0.0)
out: PointArray = P[idx, :] + t[:, None] * (P[idx + 1, :] - P[idx, :])
return out
[docs]
def resample_polyline(
points: PointArray, fractions: FloatArray, *, closed: bool = False,
) -> PointArray:
"""Arc-length resample of a polyline at normalized ``fractions`` in ``[0, 1]``. With
``closed=True`` the polyline wraps (fraction ``1`` returns to the start). Returns
``(len(fractions), 3)``."""
P: PointArray = np.asarray(points, dtype=float).reshape(-1, 3)
fr: FloatArray = np.atleast_1d(np.asarray(fractions, dtype=float))
path: PointArray = np.vstack([P, P[0:1, :]]) if closed else P
arclen: FloatArray = np.concatenate(
[[0.0], np.cumsum(np.sqrt(np.sum(np.diff(path, axis=0) ** 2, axis=1)))])
targets: FloatArray = np.clip(fr, 0.0, 1.0) * arclen[-1]
return _lerp_along(path, arclen, targets)
def _resample_loop(pts: PointArray, n: int) -> PointArray:
"""Arc-length resample of a closed loop to ``n`` points (wrapping, endpoint
excluded)."""
return resample_polyline(pts, np.linspace(0.0, 1.0, n, endpoint=False), closed=True)
def _align_loop(A: PointArray, B: PointArray) -> PointArray:
"""Cyclically shift (and possibly flip) closed loop ``A`` to best match ``B`` in
least squares; returns the reordered points."""
M = A.shape[0]
best = np.inf
bestA: PointArray = A
for f in (0, 1):
Af: PointArray = A[::-1, :].copy() if f else A
for s in range(M):
As: PointArray = np.roll(Af, s, axis=0)
d = float(np.sum((As - B) ** 2))
if d < best:
best = d
bestA = As
return bestA
def _split_ring(pts: PointArray, f: float, nh: int) -> PointArray:
"""Resample closed loop ``pts`` into ``2*nh`` points so index 0 is ``pts[0]`` and
index ``nh`` sits at arc-length fraction ``f``. Returns ``(2*nh, 3)``."""
R: PointArray = np.asarray(pts, dtype=float).reshape(-1, 3)
Rc: PointArray = np.vstack([R, R[0:1, :]])
al: FloatArray = np.concatenate(
[[0.0], np.cumsum(np.sqrt(np.sum(np.diff(Rc, axis=0) ** 2, axis=1)))])
tot = al[-1]
p1 = _lerp_along(Rc, al, np.clip(np.linspace(0.0, f * tot, nh + 1), 0.0, tot))
p2 = _lerp_along(Rc, al, np.clip(np.linspace(f * tot, tot, nh + 1), 0.0, tot))
return np.vstack([p1[0:nh, :], p2[0:nh, :]])
# -- projection ---------------------------------------------------------
def _vertex_fans(T: IntArray, nv: int) -> list[IntArray]:
"""``[triangles touching vertex v]`` for every ``v``, assembled as a sparse matrix
rather than by appending to ``nv`` Python lists once per triangle corner."""
tri = np.asarray(T, dtype=np.int64).reshape(-1, 3)
rows = tri.ravel()
cols: IntArray = np.repeat(np.arange(tri.shape[0], dtype=np.int64), 3)
m = sp.coo_matrix((np.ones(rows.shape[0], dtype=np.int8), (rows, cols)),
shape=(nv, max(tri.shape[0], 1))).tocsr()
ind, ptr = m.indices, m.indptr
return [np.asarray(ind[ptr[v]:ptr[v + 1]], dtype=np.int64) for v in range(nv)]
[docs]
def project_points(surface: TriMesh, P: PointArray) -> PointArray:
"""Snap points onto the nearest vertex's triangle fan (points assumed near the surface)."""
P = np.atleast_2d(np.asarray(P, dtype=float))
Vx, T = surface.points, surface.tris
nv = Vx.shape[0]
VT: list[IntArray] = _vertex_fans(T, nv)
# one batched tree query instead of a full O(V) scan per point. Only the *search*
# changes: the fan chosen and the arithmetic on it are untouched. Exact ties are
# the one place this can differ -- ``argmin`` took the lowest vertex index, a tree
# takes whichever it reaches first -- but a tie means two vertices equidistant from
# the query, and their fans overlap on the triangles that matter.
nearest: IntArray = cKDTree(Vx).query(P)[1]
Q = P.copy()
for i in range(P.shape[0]):
p = P[i, :]
vs = int(nearest[i])
best = np.inf
bq = Vx[vs, :]
for e in VT[vs]:
q = _closest_on_tri(p, Vx[T[e, 0], :], Vx[T[e, 1], :], Vx[T[e, 2], :])
d = np.sum((q - p) ** 2)
if d < best:
best = d
bq = q
Q[i, :] = bq
return Q
[docs]
def project_to_surface(
surface: TriMesh, P: PointArray, faces: IntArray | None = None,
) -> PointArray:
"""Closest-point projection over triangles; ``faces`` defaults to all triangles."""
P = np.atleast_2d(np.asarray(P, dtype=float))
Vx = surface.points
T = surface.tris if faces is None else np.asarray(faces, dtype=np.int64).reshape(-1, 3)
m = P.shape[0]
Q = P.copy()
best = np.full(m, np.inf)
for e, pts in _projection_candidates(Vx, T, P):
A, B, C = Vx[T[e, 0], :], Vx[T[e, 1], :], Vx[T[e, 2], :]
q, d2 = _closest_on_tri_vec(P[pts, :], A, B, C)
upd = d2 < best[pts]
if np.any(upd):
sel = pts[upd]
best[sel] = d2[upd]
Q[sel, :] = q[upd, :]
return Q
def _projection_candidates(
Vx: PointArray, T: IntArray, P: PointArray,
) -> Iterator[tuple[int, IntArray]]:
"""``(triangle, the query points it could possibly be closest to)``, ascending by
triangle -- the broad phase that keeps :func:`project_to_surface` off an O(T*P) scan.
Exactness, not approximation. ``r0`` is the distance from a query point to the
nearest *vertex* of the surface, so the true closest point is at most ``r0`` away;
a triangle holding a point that close has its centroid within ``r0 + reach``, where
``reach`` is the largest centroid-to-vertex distance of any triangle. Every triangle
that could win is therefore in the ball, and the ones dropped provably could not.
Yielding in ascending triangle order matters as much as the bound: the caller keeps
a strict ``d2 < best``, so on an exact tie the lowest triangle index wins, exactly as
the full scan did. Same kernel, same arithmetic, same reduction order -- the result
is bit-identical to testing every triangle."""
if T.shape[0] == 0 or P.shape[0] == 0:
return
cen: PointArray = Vx[T].mean(axis=1)
reach = float(np.max(np.linalg.norm(Vx[T] - cen[:, None, :], axis=2)))
# the bound must be taken over *this* triangle set's own vertices, not the whole
# surface's: with a ``faces`` subset the nearest vertex overall can be nearer than
# anything in the subset, which would shrink the ball below the true answer and drop
# the winning triangle.
r0: FloatArray = cKDTree(Vx[np.unique(T)]).query(P)[0]
hits = cKDTree(cen).query_ball_point(P, r=r0 + reach)
counts: IntArray = np.fromiter((len(h) for h in hits), dtype=np.int64,
count=len(hits))
if not counts.sum():
return
tri: IntArray = np.concatenate([np.asarray(h, dtype=np.int64)
for h in hits if len(h)])
pt: IntArray = np.repeat(np.arange(P.shape[0], dtype=np.int64), counts)
# group the flat (triangle, point) pairs by triangle, ascending
order = np.argsort(tri, kind="stable")
tri, pt = tri[order], pt[order]
edges: IntArray = np.flatnonzero(np.diff(tri)) + 1
for lo, hi in zip(np.concatenate([[0], edges]),
np.concatenate([edges, [tri.shape[0]]])):
yield int(tri[lo]), pt[lo:hi]
# -- closest-point helpers ----------------------------------------------
def _closest_on_tri(p: Point, a: Point, b: Point, c: Point) -> Point:
ab = b - a
ac = c - a
ap = p - a
d1 = np.dot(ab, ap)
d2 = np.dot(ac, ap)
if d1 <= 0 and d2 <= 0:
return a
bp = p - b
d3 = np.dot(ab, bp)
d4 = np.dot(ac, bp)
if d3 >= 0 and d4 <= d3:
return b
vc = d1 * d4 - d3 * d2
if vc <= 0 and d1 >= 0 and d3 <= 0:
return a + (d1 / (d1 - d3)) * ab
cp = p - c
d5 = np.dot(ab, cp)
d6 = np.dot(ac, cp)
if d6 >= 0 and d5 <= d6:
return c
vb = d5 * d2 - d1 * d6
if vb <= 0 and d2 >= 0 and d6 <= 0:
return a + (d2 / (d2 - d6)) * ac
va = d3 * d6 - d5 * d4
if va <= 0 and (d4 - d3) >= 0 and (d5 - d6) >= 0:
w = (d4 - d3) / ((d4 - d3) + (d5 - d6))
return b + w * (c - b)
den = 1.0 / (va + vb + vc)
v = vb * den
w = vc * den
return a + ab * v + ac * w
def _closest_on_seg_vec(
P: PointArray, U: Point, V: Point,
) -> tuple[PointArray, FloatArray]:
uv = V - U
L2 = np.dot(uv, uv)
if L2 == 0:
q = np.tile(U, (P.shape[0], 1))
else:
t = np.clip(((P - U) @ uv) / L2, 0.0, 1.0)
q = U + t[:, None] * uv
return q, np.sum((q - P) ** 2, axis=1)
def _closest_on_tri_vec(
P: PointArray, A: Point, B: Point, C: Point,
) -> tuple[PointArray, FloatArray]:
m = P.shape[0]
INF = np.full(m, np.inf)
qab, dab = _closest_on_seg_vec(P, A, B)
qbc, dbc = _closest_on_seg_vec(P, B, C)
qca, dca = _closest_on_seg_vec(P, C, A)
ab = B - A
ac = C - A
n = np.cross(ab, ac)
nn = np.dot(n, n)
if nn > 0:
ap = P - A
t = (ap @ n) / nn
qf = P - t[:, None] * n
v0, v1, v2 = ab, ac, qf - A
d00 = np.dot(v0, v0)
d01 = np.dot(v0, v1)
d11 = np.dot(v1, v1)
d20 = v2 @ v0
d21 = v2 @ v1
den = d00 * d11 - d01 * d01
v = (d11 * d20 - d01 * d21) / den
w = (d00 * d21 - d01 * d20) / den
u = 1 - v - w
inside = (u >= 0) & (v >= 0) & (w >= 0)
df = np.where(inside, np.sum((qf - P) ** 2, axis=1), INF)
else:
qf, df = P.copy(), INF
q = qab.copy()
d2 = dab.copy()
for qc, dc in ((qbc, dbc), (qca, dca), (qf, df)):
upd = dc < d2
q[upd, :] = qc[upd, :]
d2 = np.minimum(d2, dc)
return q, d2