Meshes¶
Routines that treat a triangle mesh as a graph over its vertices, with the edges taken from the faces.
Connected components¶
Labelling the connected components of a triangle mesh, i.e. finding which vertices are reachable from which through shared faces.
import navis_fastcore as fastcore
import numpy as np
# Two disconnected triangles
faces = np.array([[0, 1, 2], [3, 4, 5]], dtype=np.uint32)
fastcore.mesh_connected_components(faces, n_vertices=6)
# array([0, 0, 0, 3, 3, 3], dtype=uint32)
connectivity picks which of the three readings of "connected" you get, each strictly
finer than the one before it. The default, "vertex", joins two vertices whenever a face
names them both and labels every vertex; "face" joins two faces wherever they share an
edge and labels every face; "manifold" joins two faces only across an edge carrying
exactly two of them.
Each step drops a kind of junction. "vertex" → "face" drops the pinch points:
# Two triangles meeting at vertex 2 and nowhere else
faces = np.array([[0, 1, 2], [2, 3, 4]], dtype=np.uint32)
fastcore.mesh_connected_components(faces, n_vertices=5)
# array([0, 0, 0, 0, 0], dtype=uint32) -- one component: you can walk through the pinch
fastcore.mesh_connected_components(faces, connectivity="face")
# array([0, 1], dtype=uint32) -- two: you cannot step across it
"face" → "manifold" drops the seams — an edge three or more faces deep belongs to no
single surface:
# Three fins meeting along the spine (1, 2)
faces = np.array([[1, 2, 3], [1, 2, 4], [1, 2, 5]], dtype=np.uint32)
fastcore.mesh_connected_components(faces, connectivity="face")
# array([0, 0, 0], dtype=uint32) -- the spine is a shared edge like any other
fastcore.mesh_connected_components(faces, connectivity="manifold")
# array([0, 1, 2], dtype=uint32) -- three faces on it, so it joins nothing
Pick by what the components are for:
| joins across | use it for | |
|---|---|---|
"vertex" |
a shared corner | "can these vertices reach each other along mesh edges" — what geodesic_matrix_mesh answers with distances |
"face" |
any shared edge | splitting a mesh into the pieces you could walk across — trimesh's split(only_watertight=False) |
"manifold" |
an edge with exactly two faces | splitting it into pieces that are surfaces, each with a well-defined inside — before asking one for its volume or winding. Reproduces trimesh's face_adjacency |
Note there is no per-vertex form of the face answers: a pinch vertex belongs to several face components at once. A boundary edge — one face — joins nothing under any reading.
navis_fastcore.mesh_connected_components(faces, n_vertices=None, connectivity='vertex', threads=None)
¶
Find connected components of a triangle mesh.
Three readings of "connected", chosen with connectivity, each strictly finer
than the one before it:
"vertex"(default) joins two vertices whenever a face names them both, so a face joins its three corners and the answer is one label per vertex."face"joins two faces wherever they share an edge, and the answer is one label per face."manifold"joins two faces only across an edge carrying exactly two of them, and is likewise one label per face.
Each step drops a kind of junction. Going from "vertex" to "face" drops the
pinch points: two triangles meeting at a single corner are one component under the
first and two under the second, because there is no edge to step across. Going from
"face" to "manifold" drops the seams: three sheets meeting along one edge are
one component under "face" and three under "manifold", because an edge that
deep belongs to no single surface.
So pick by what the components are for. "vertex" answers "can these vertices
reach each other along mesh edges" — the question
:func:~navis_fastcore.geodesic_matrix_mesh answers with distances. "face"
splits a mesh into the pieces you could walk across, which is trimesh's
split(only_watertight=False). "manifold" splits it into pieces that are
surfaces, each with a well-defined inside — what you want before asking a piece for
its volume or its winding — and reproduces trimesh's face_adjacency exactly.
Face components cannot be reported per vertex, incidentally, and that is not an oversight in the interface: a pinch vertex belongs to several face components at once, so there is no per-vertex form of that answer.
All three are Union-Find (DSU) with path-halving. "vertex" builds no adjacency at
all — one serial sweep of the faces over a single integer array of length
n_vertices. The other two group the 3 * F edges the faces name first, which is
a parallel sort and the only reason they take threads; they then differ by one
test on how many faces each edge came back with, so "manifold" costs no more than
"face".
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
n_vertices
|
TYPE:
|
connectivity
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
components
|
For
TYPE:
|
Notes
n_vertices and threads each belong to particular connectivities, and passing
one where it does not apply is an error rather than something quietly dropped — the
rule :func:~navis_fastcore.smooth_mesh follows for its per-method parameters.
Sizing a call with n_vertices and then asking for connectivity="face" would
otherwise hand back an array of a different length than the caller expected, which is
the one failure that looks like success.
An edge with a single face on it — the boundary of an open mesh — joins nothing under
any reading; there is no second face to join it to. Self-loop edges from degenerate
faces are kept throughout, as in :func:~navis_fastcore.unique_edges.
"manifold" is trimesh.graph.face_adjacency, which builds its adjacency with
group_rows(edges_sorted, require_count=2) and so drops a deeper edge outright.
"face" is the reading the rest of this module uses, where a shared edge is a
shared edge however many faces are on it.
Examples:
Two triangles sharing an edge — one component either way:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
>>> fastcore.mesh_connected_components(faces, n_vertices=4)
array([0, 0, 0, 0], dtype=uint32)
>>> fastcore.mesh_connected_components(faces, connectivity="face")
array([0, 0], dtype=uint32)
Two disjoint triangles — two components, of three vertices or of one face:
>>> faces = np.array([[0, 1, 2], [3, 4, 5]], dtype=np.uint32)
>>> fastcore.mesh_connected_components(faces, n_vertices=6)
array([0, 0, 0, 3, 3, 3], dtype=uint32)
>>> fastcore.mesh_connected_components(faces, connectivity="face")
array([0, 1], dtype=uint32)
Where "vertex" and "face" part company — two triangles pinched together at
vertex 2. The vertex graph walks straight through it, the faces cannot:
>>> faces = np.array([[0, 1, 2], [2, 3, 4]], dtype=np.uint32)
>>> fastcore.mesh_connected_components(faces, n_vertices=5)
array([0, 0, 0, 0, 0], dtype=uint32)
>>> fastcore.mesh_connected_components(faces, connectivity="face")
array([0, 1], dtype=uint32)
And where "face" and "manifold" do — three fins meeting along the spine
(1, 2). That edge carries three faces, so it is a seam no one surface owns:
>>> faces = np.array([[1, 2, 3], [1, 2, 4], [1, 2, 5]], dtype=np.uint32)
>>> fastcore.mesh_connected_components(faces, connectivity="face")
array([0, 0, 0], dtype=uint32)
>>> fastcore.mesh_connected_components(faces, connectivity="manifold")
array([0, 1, 2], dtype=uint32)
Face labels are face indices, so they group the rows of faces directly:
Geodesic distances¶
The mesh counterpart to geodesic_matrix, which works on skeletons. A skeleton
is a tree, so distances there come from walking to the lowest common ancestor. A mesh is a
general cyclic graph, so this runs a Dijkstra per source instead — in parallel, which is
the whole point: scipy.sparse.csgraph.dijkstra holds the GIL, so you cannot get that
speedup from Python by threading it yourself.
import navis_fastcore as fastcore
import numpy as np
# Two triangles sharing the 1-2 edge, forming a unit square
faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
vertices = np.array([[0, 0, 0],
[1, 0, 0],
[0, 1, 0],
[1, 1, 0]], dtype=np.float64)
fastcore.geodesic_matrix_mesh(faces, vertices)
# array([[0. , 1. , 1. , 2. ],
# [1. , 0. , 1.4142135, 1. ],
# [1. , 1.4142135, 0. , 1. ],
# [2. , 1. , 1. , 0. ]], dtype=float32)
Mind the size of the output
A full V x V matrix is around 107 GB at V=164k, so for anything but a small mesh you
want sources and/or targets.
targets is worth calling out. scipy.sparse.csgraph.dijkstra has no notion of
targets: it always materialises all V columns and makes you slice afterwards. Passing
targets here means only those columns are ever allocated — for 200 sources and 100
targets on a 41k-vertex mesh that is 0.03 MB instead of 70 MB.
If you only need the nearest (or farthest) target, use
geodesic_nearest_mesh instead — its output is
O(sources) rather than O(sources x targets), and it is faster too, because the
search stops at the first target it settles.
This is the along-edge distance
Shortest paths are constrained to run along mesh edges, so on a coarse mesh they overshoot the true surface geodesic. This is the same approximation navis makes.
navis_fastcore.geodesic_matrix_mesh(faces, vertices=None, n_vertices=None, sources=None, targets=None, limit=None, threads=None, dtype=None)
¶
Calculate geodesic ("along-the-mesh-edge") distances on a triangle mesh.
This is the mesh counterpart to :func:~navis_fastcore.geodesic_matrix, which
works on skeletons. Where the skeleton version exploits the tree structure, a
mesh is a general cyclic graph, so this runs a parallel Dijkstra (or a BFS when
unweighted) over the vertex adjacency derived from faces.
Notes
This is the distance along mesh edges, not the exact surface geodesic: shortest paths are constrained to run along edges, so on a coarse mesh they overshoot the true surface distance.
Beware the size of the output. A full V x V matrix is ~107 GB at V=164k, so
for anything but a small mesh you want sources and/or targets. Unlike
scipy.sparse.csgraph.dijkstra — which has no notion of targets and always
materialises all V columns before you can slice them — targets here means
only those columns are ever allocated.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
n_vertices
|
TYPE:
|
sources
|
TYPE:
|
targets
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
matrix
|
Geodesic distances in the resolved
TYPE:
|
Examples:
Two triangles sharing the 1-2 edge, forming a unit square. Vertices 0 and 3 are the opposite corners, so they are not directly connected — the shortest path between them goes around, via 1 or 2:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
>>> vertices = np.array([[0, 0, 0],
... [1, 0, 0],
... [0, 1, 0],
... [1, 1, 0]], dtype=np.float64)
>>> fastcore.geodesic_matrix_mesh(faces, vertices)
array([[0. , 1. , 1. , 2. ],
[1. , 0. , 1.4142135, 1. ],
[1. , 1.4142135, 0. , 1. ],
[2. , 1. , 1. , 0. ]], dtype=float32)
Without vertices every edge has weight 1, so you get hop counts instead — the
shared diagonal 1-2 is now a single hop rather than sqrt(2):
navis_fastcore.geodesic_nearest_mesh(faces, vertices=None, n_vertices=None, sources=None, targets=None, limit=None, threads=None, dtype=None)
¶
For each source vertex, find the nearest target vertex on a mesh.
A memory-efficient alternative to :func:~navis_fastcore.geodesic_matrix_mesh:
it keeps only the nearest target and the distance to it, so the output is
O(len(sources)) rather than O(len(sources) * len(targets)). It is also
faster, because the search stops at the first target it settles instead of
exploring the whole connected component.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
n_vertices
|
TYPE:
|
sources
|
TYPE:
|
targets
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
distances
|
Distance from each source to its nearest target, in the resolved
TYPE:
|
nearest
|
Vertex index of that nearest target;
TYPE:
|
Notes
A source that is itself a target is matched to its nearest distinct target, never to itself (so the distance is never trivially 0). Ties break towards the lower vertex index.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
>>> vertices = np.array([[0, 0, 0],
... [1, 0, 0],
... [0, 1, 0],
... [1, 1, 0]], dtype=np.float64)
>>> dists, nearest = fastcore.geodesic_nearest_mesh(
... faces, vertices, sources=[0], targets=[2, 3]
... )
>>> dists
array([1.], dtype=float32)
>>> nearest
array([2], dtype=int32)
navis_fastcore.geodesic_farthest_mesh(faces, vertices=None, n_vertices=None, sources=None, targets=None, limit=None, threads=None, dtype=None)
¶
For each source vertex, find the farthest target vertex on a mesh.
The mirror image of :func:~navis_fastcore.geodesic_nearest_mesh, with the same
O(len(sources)) memory footprint. Unlike nearest, this cannot stop early
— it has to settle every target — but the farthest one then comes for free, since
the search settles vertices in increasing order of distance.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
n_vertices
|
TYPE:
|
sources
|
TYPE:
|
targets
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
distances
|
Distance from each source to its farthest target, in the resolved
TYPE:
|
farthest
|
Vertex index of that farthest target;
TYPE:
|
Notes
As with nearest, a source that is itself a target is matched to a distinct
target, never to itself.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
>>> vertices = np.array([[0, 0, 0],
... [1, 0, 0],
... [0, 1, 0],
... [1, 1, 0]], dtype=np.float64)
>>> dists, farthest = fastcore.geodesic_farthest_mesh(
... faces, vertices, sources=[0], targets=[2, 3]
... )
>>> farthest
array([3], dtype=int32)
Arbitrary graphs¶
The same kernel, over an explicit edge list rather than a mesh. Unlike
geodesic_matrix, this makes no tree assumption, so cycles are fine.
import navis_fastcore as fastcore
import numpy as np
# A triangle. The direct 0-2 edge has weight 5, so the shortest
# path between them goes the long way round via 1.
edges = np.array([[0, 1], [1, 2], [2, 0]], dtype=np.uint32)
weights = np.array([1, 1, 5], dtype=np.float32)
fastcore.geodesic_matrix_graph(edges, 3, weights=weights)
# array([[0., 1., 2.],
# [1., 0., 1.],
# [2., 1., 0.]], dtype=float32)
# The distances follow the weights' dtype: float64 in, float64 out.
fastcore.geodesic_matrix_graph(edges, 3, weights=weights.astype(np.float64))
# array([[0., 1., 2.],
# [1., 0., 1.],
# [2., 1., 0.]])
See Float return dtypes for when the wider width is
worth asking for, and for the dtype argument that overrides the default in either
direction.
navis_fastcore.geodesic_matrix_graph(edges, n_nodes, weights=None, directed=False, sources=None, targets=None, limit=None, threads=None, dtype=None)
¶
Calculate geodesic distances over an arbitrary graph.
The general form of :func:~navis_fastcore.geodesic_matrix_mesh. Unlike
:func:~navis_fastcore.geodesic_matrix, this makes no tree assumption — cycles
are fine.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
weights
|
TYPE:
|
directed
|
TYPE:
|
sources
|
TYPE:
|
targets
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
matrix
|
Geodesic distances in the resolved dtype;
TYPE:
|
Notes
Dijkstra sums one weight per hop, so a path of k hops carries up to k
roundings. float32 is right for mesh and skeleton work — a 24-bit mantissa
resolves a 100 mm neuron to ~6 nm, and the matrix is by far the largest thing
this allocates. float64 earns its keep when the accumulation is long rather
than the graph large (tens of thousands of hops), when weights span a wide
dynamic range, or when you are comparing against
scipy.sparse.csgraph, which works in float64 unconditionally.
Examples:
A triangle — a cycle, which the skeleton functions would reject. Note the direct 0-2 edge has weight 5, so the shortest path goes the long way round via 1:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 0]], dtype=np.uint32)
>>> weights = np.array([1, 1, 5], dtype=np.float32)
>>> fastcore.geodesic_matrix_graph(edges, 3, weights=weights)
array([[0., 1., 2.],
[1., 0., 1.],
[2., 1., 0.]], dtype=float32)
Hand it float64 weights and the distances come back float64:
>>> fastcore.geodesic_matrix_graph(edges, 3, weights=weights.astype(np.float64))
array([[0., 1., 2.],
[1., 0., 1.],
[2., 1., 0.]])
dtype overrides that in either direction — here asking for float64 from a
float32 input, which is what you want when the weights were measured coarsely
but the paths are long enough for the accumulation to matter:
Graph primitives¶
The handful of traversal operations that mesh algorithms actually need, taken straight off an edge list. These exist because reaching for a general-purpose graph library means paying to build a graph object first — on a 41k-vertex mesh that construction alone costs more than every query you then run against it.
import navis_fastcore as fastcore
import numpy as np
# A path 0-1-2, a lone edge 3-4, and an isolated node 5
edges = np.array([[0, 1], [1, 2], [3, 4]], dtype=np.uint32)
fastcore.connected_components_graph(edges, n_nodes=6)
# array([0, 0, 0, 3, 3, 5], dtype=uint32)
navis_fastcore.connected_components_graph(edges, n_nodes)
¶
Find connected components of a graph given as an edge list.
The edge-list counterpart of :func:~navis_fastcore.mesh_connected_components,
using the same Union-Find: a single integer array of length n_nodes, no
adjacency list. Use this when the graph is not a triangle mesh, or when you
already hold the deduplicated edges and would rather not walk the faces again.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
components
|
For each node, the smallest node index in its component.
TYPE:
|
Examples:
A path 0-1-2, a lone edge 3-4, and an isolated node 5:
Level sets¶
level_set_components is the one worth knowing
about. Given a label per node it finds the connected components of every label's
induced subgraph in one pass — which is the inner loop of wavefront-style mesh
skeletonization, where the label is a binned geodesic distance and each component is one
ring around the structure.
Done conventionally that loop costs one subgraph construction plus one component search per distinct level. Here it is a single sweep over the edges, unioning an edge only when its endpoints agree:
import navis_fastcore as fastcore
import numpy as np
faces = ... # your mesh
edges = fastcore.unique_edges(faces)
n = len(vertices)
# Cast a wave from vertex 0 and collapse each ring
dist = fastcore.geodesic_matrix_mesh(faces, n_vertices=n, sources=[0])[0]
rings, n_rings = fastcore.level_set_components(edges, n, dist.astype(np.int64))
# Ring ids are contiguous, so aggregating is a plain bincount
sizes = np.bincount(rings[rings >= 0], minlength=n_rings)
Note that dist is -1 where the search could not reach, and negative labels are
excluded rather than grouped — so an unreachable region does not become one bogus
level.
On a 41k-vertex mesh with ~200 levels this runs in ~0.3 ms against ~12 ms for the per-level-subgraph equivalent, on top of the ~28 ms of graph construction it avoids entirely.
navis_fastcore.level_set_components(edges, n_nodes, labels)
¶
Find the connected components of every level set at once.
Given a label per node, this finds the connected components of each subgraph induced by the nodes sharing a label — all labels in a single pass, by unioning an edge only when its two endpoints agree.
This is the primitive behind "which nodes were reached by the same wavefront
and are actually touching", where labels is a (binned) geodesic distance
and each component is one ring around the structure.
The point is that it replaces a loop. With a general-purpose graph library
the same result costs one induced-subgraph construction plus one component
search per distinct label, so a mesh with a thousand levels pays a thousand
graph builds; here it is one O(E) sweep over the edges, and the only
allocations are three n_nodes-sized integer arrays.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
labels
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ids
|
Component of each node in
TYPE:
|
n_components
|
Number of components found.
TYPE:
|
Examples:
A path 0-1-2-3-4 labelled 0, 0, 0, 1, 1: one run per label, so two
components.
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.uint32)
>>> ids, n = fastcore.level_set_components(edges, 5, [0, 0, 0, 1, 1])
>>> ids
array([0, 0, 0, 1, 1], dtype=int32)
Nodes sharing a label but not touching stay separate — here label 0 appears at both ends of the path:
>>> ids, n = fastcore.level_set_components(edges, 5, [0, 1, 1, 1, 0])
>>> ids
array([0, 1, 1, 1, 2], dtype=int32)
Aggregating per component is then a plain np.bincount:
navis_fastcore.contract_vertices(edges, mapping, threads=None)
¶
Contract nodes onto new ids and return the simplified edge list.
Both endpoints of every edge are pushed through mapping; edges that end up
with both ends on the same new node (self-loops) are dropped, and the rest are
deduplicated. This is igraph's contract_vertices() followed by
simplify(), fused — and, unlike igraph's version, it does not rewrite a
graph object in place, so contracting does not cost a copy of the graph.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
mapping
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
edges
|
The surviving edges as
TYPE:
|
Examples:
A square 0-1-2-3 with a diagonal, collapsing {0, 1} -> 0 and
{2, 3} -> 1. The 0-1 and 2-3 edges become self-loops and vanish; the
remaining three all become 0-1 and collapse to a single edge:
navis_fastcore.minimum_spanning_tree(edges, n_nodes, weights=None, maximize=False, threads=None)
¶
Find the minimum (or maximum) spanning forest of a graph.
Kruskal's algorithm on the same Union-Find as the component search: sort the
edges by weight, keep the ones that join two different components.
Disconnected input is fine — each component contributes its own tree, so this
is really a spanning forest, matching igraph's spanning_tree() and
scipy's minimum_spanning_tree.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
weights
|
TYPE:
|
maximize
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
indices
|
Row indices into
TYPE:
|
Examples:
A triangle with weights 1, 2, 3 — the spanning tree takes the two cheap edges and rejects the one that would close the cycle:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [0, 2]], dtype=np.uint32)
>>> weights = np.array([1, 2, 3], dtype=np.float32)
>>> keep = fastcore.minimum_spanning_tree(edges, 3, weights)
>>> edges[keep]
array([[0, 1],
[1, 2]], dtype=uint32)
Ask for the maximum instead and it takes the two expensive ones:
Turning an edge list into a tree¶
minimum_spanning_tree picks which edges
survive. parents_from_edges picks which way they point —
which is what turns a bag of undirected edges into something you can walk, root, or write
out as SWC. Cycles in the input are fine; each component contributes a spanning tree of
itself, so this doubles as the cycle-breaker networkx.bfs_tree is usually pressed into.
import navis_fastcore as fastcore
# Break cycles, orient, and re-index so parents come before their children
keep = fastcore.minimum_spanning_tree(edges, n, weights=lengths)
parents, order = fastcore.parents_from_edges(edges[keep], n)
new_ids = np.empty(n, dtype=np.int64)
new_ids[order] = np.arange(n)
swc_parents = np.where(parents < 0, -1, new_ids[parents])[order]
order is the second return for exactly that last step. A node always settles after its
parent, so relabelling by it guarantees parents get lower ids than their children — the SWC
requirement — and it comes free, since the search already visits nodes in that order.
One search, not one per component. The obvious construction is a shortest-path tree per
component, which is what geodesic_predecessors
gives you — and it costs O(components x n_nodes) in output alone. On a skeleton that
shatters into four thousand fragments that is a 2 GB array to answer a question whose answer
is one n_nodes-long column. Here the components are swept one after another into that
single column, so the cost is O(V + E) however finely the graph is fragmented:
| 100k-node graph | fastcore | igraph (BFS per component) | networkx (bfs_tree per component) |
|---|---|---|---|
| one arbor | 2.9 ms | 14 ms | 365 ms |
| ~4000 fragments | 2.7 ms | 4370 ms | 285 ms |
Weights are optional and change what you get: None gives the breadth-first tree, weights
give the shortest-path tree. Neither is the minimum spanning tree — for that, run
minimum_spanning_tree first and orient what it keeps, as above.
navis_fastcore.parents_from_edges(edges, n_nodes, weights=None, roots=None)
¶
Orient a graph into a rooted spanning forest — one parent per node.
The missing half of "I have an edge list and I want a tree".
:func:~navis_fastcore.minimum_spanning_tree picks which edges survive; this
picks which way they point, which is what turns a bag of undirected edges into
something you can walk, root, or write out as SWC. Cycles in the input are
fine — each component contributes a spanning tree of itself, so this doubles as
the cycle-breaker networkx.bfs_tree is usually pressed into.
One search covers the whole graph. The obvious construction — a shortest-path
tree per component — is what
:func:~navis_fastcore.geodesic_predecessors gives you, and it costs
O(components * n_nodes) in output alone: on a mesh that shatters into four
thousand specks that is a two-gigabyte array to answer a question whose answer
is one n_nodes-long column. Here the components are swept one after another
into that single column, so the cost is O(V + E) however finely the graph is
fragmented.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
weights
|
TYPE:
|
roots
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
parents
|
Parent of each node,
TYPE:
|
order
|
Every node in the order it settled. A node always settles after its
parent, so this is a topological order — relabel by it and parents are
guaranteed to have lower ids than their children, which is exactly the
SWC requirement. It comes free: the search already visits nodes in this
order, and deriving it afterwards from
TYPE:
|
Notes
Among equal-length routes the parent is whichever settled first, which is deterministic but otherwise arbitrary — as it is for any spanning tree of a graph with more than one.
Examples:
A path, written "backwards" to show the orientation comes from the search and not from the order the endpoints happen to be in:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[1, 0], [2, 1], [3, 2]], dtype=np.uint32)
>>> parents, order = fastcore.parents_from_edges(edges, 4)
>>> parents
array([-1, 0, 1, 2], dtype=int32)
Root it at the far end instead and every link reverses:
>>> parents, order = fastcore.parents_from_edges(edges, 4, roots=[3])
>>> parents
array([ 1, 2, 3, -1], dtype=int32)
>>> order
array([3, 2, 1, 0], dtype=uint32)
Cycles are broken; two components each get their own root:
>>> edges = np.array([[0, 1], [1, 2], [2, 0], [4, 5]], dtype=np.uint32)
>>> parents, order = fastcore.parents_from_edges(edges, 6)
>>> parents
array([-1, 0, 0, -1, -1, 4], dtype=int32)
order relabels a forest so parents come before their children:
Which edges are load-bearing¶
bridges is the counterpart to
minimum_spanning_tree rather than a variant of
it: the MST asks which edges to keep to stay connected, this asks which ones may not be
dropped. That is the question behind "prune this graph but do not shatter it", where you
have a set of edges you would like gone and need to know which of them are load-bearing.
# Drop the edges you don't want -- except the ones holding the graph together
unwanted = ... # bool mask over `edges`
safe = unwanted & ~fastcore.bridges(edges, n)
edges = edges[~safe]
Parallel edges are honoured: two nodes joined twice are joined by a cycle, so neither of those edges is a bridge. That is why this does not share the deduplicated adjacency the geodesic searches use — that would fuse a parallel pair into one edge and report a bridge that is not there. Self-loops are never bridges.
It is Tarjan's algorithm on an explicit stack, so a mesh strip tens of thousands of vertices
long does not overflow anything. Against igraph's Graph.bridges() on a 100k-node graph:
2.6 ms against 13.5 ms for one arbor, 2.2 ms against 207 ms once it fragments.
navis_fastcore.bridges(edges, n_nodes)
¶
Find the edges whose removal would disconnect their component.
Tarjan's algorithm: one depth-first sweep tracking, per node, the earliest
node reachable from its subtree by a single back edge. A tree edge (u, v)
is a bridge exactly when nothing under v can climb above it, i.e. there is
no second route around it.
The counterpart to :func:~navis_fastcore.minimum_spanning_tree rather than a
variant of it: the MST asks which edges to keep to stay connected, this asks
which ones may not be dropped. That is the question behind "prune this graph
but do not shatter it", where you have a set of edges you would like gone and
need to know which of them are load-bearing.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
mask
|
TYPE:
|
Notes
Parallel edges are honoured: two nodes joined twice are joined by a cycle, so
neither of those edges is a bridge. Self-loops are never bridges. This is why
bridges does not share the deduplicated adjacency the geodesic searches in
this module use — that would fuse a parallel pair into one edge and report a
bridge that is not there.
Examples:
Every edge of a tree is a bridge:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> path = np.array([[0, 1], [1, 2], [2, 3]], dtype=np.uint32)
>>> fastcore.bridges(path, 4)
array([ True, True, True])
Close it into a ring and none of them is:
>>> ring = np.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=np.uint32)
>>> fastcore.bridges(ring, 4)
array([False, False, False, False])
Two triangles joined by a single edge — only the link:
Spanning a subset by geodesic distance¶
Skeletonization ends up here: the mesh has been thinned to a scatter of surviving vertices
that must be rejoined along the surface rather than through space. That is a minimum
spanning tree over a k-node subset, weighted by geodesic distance in the mesh underneath.
The obvious route is to ask for the k x k geodesic matrix and hand it to a matrix MST.
That materialises k**2 distances to use k - 1 of them — 400 MB at k = 10_000, before
the O(k^2) MST itself — and needs k separate searches to fill.
geodesic_mst_mesh never forms the matrix:
edges, weights = fastcore.geodesic_mst_mesh(faces, keep, vertices)
# Rows index `keep`, so map them back yourself
skeleton_edges = keep[edges]
Following Mehlhorn's construction for the distance network, one multi-source search
partitions every vertex by which of keep is nearest, and then each mesh edge whose
endpoints fall in different cells offers one candidate: joining their two owners at
d(u) + w(u, v) + d(v). An MST over those candidates is an MST of the full distance
network, so one sweep and one Kruskal replace k searches and a dense matrix. The returned
weights come back exactly equal to the geodesic distances between the pairs they join, so
they are usable as lengths and not merely as an ordering.
The cost is flat in k, because it is one sweep whatever k is — which is the whole shape
of the table, on a 100k-node graph:
k |
fastcore | k x k matrix + MST |
matrix size |
|---|---|---|---|
| 250 | 12.7 ms | 187 ms | 0.3 MB |
| 1000 | 7.6 ms | 584 ms | 4 MB |
| 4000 | 8.3 ms | 7820 ms | 64 MB |
limit bounds how far apart two nodes may be and still be joined. The result is then the
MST of the graph on nodes keeping only pairs within limit, which is a forest when that
graph is disconnected — the same trade scipy.sparse.csgraph.dijkstra(limit=...) offers,
except that here it also prunes the sweep, so it buys time rather than merely discarding
results. Nodes in different components of the mesh are never joined either way.
navis_fastcore.geodesic_mst_mesh(faces, nodes, vertices=None, n_vertices=None, limit=None, threads=None, dtype=None)
¶
Minimum spanning tree over a subset of mesh vertices, by geodesic distance.
The tree that reconnects a scatter of surviving vertices through the mesh they were carved out of — the last step of a skeletonisation, where the mesh has been thinned to a few thousand vertices that must be rejoined along the surface rather than through space.
The obvious route is to ask for the k x k geodesic matrix and hand it to a
matrix MST. That materialises k**2 distances to use k - 1 of them —
400 MB at k = 10_000, before the O(k**2) MST itself — and it needs k
separate searches to fill. This never forms the matrix. Instead, following
Mehlhorn's construction for the distance network, one multi-source search
partitions every vertex by which of nodes is nearest, and then each mesh
edge whose endpoints fall in different cells offers one candidate: joining their
two owners at d(u) + w(u, v) + d(v). An MST over those candidates is an MST
of the full distance network, so one sweep and one Kruskal replace k searches
and a dense matrix.
The returned weights come back exactly equal to the geodesic distances between the pairs they join, so they are usable as lengths and not merely as an ordering.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
nodes
|
TYPE:
|
vertices
|
TYPE:
|
n_vertices
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
edges
|
Rows of positions in
TYPE:
|
weights
|
Geodesic distance across each of those edges, in the resolved
TYPE:
|
``M`` is ``len(nodes) - 1`` when every node can reach every other within
|
|
``limit``, and less when they cannot: vertices in different components of the
|
|
mesh are never joined.
|
|
Examples:
Two triangles sharing an edge, spanning three of the four vertices:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 2, 3]], dtype=np.uint32)
>>> verts = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]], dtype=float)
>>> edges, weights = fastcore.geodesic_mst_mesh(faces, [0, 1, 3], verts)
>>> edges
array([[0, 1],
[1, 2]])
>>> np.round(weights, 3)
array([1., 1.], dtype=float32)
The rows index nodes, so map them back to vertex ids yourself:
navis_fastcore.geodesic_mst_graph(edges, n_nodes, nodes, weights=None, limit=None, threads=None, dtype=None)
¶
Minimum spanning tree over a subset of graph nodes, by geodesic distance.
The edge-list form of :func:~navis_fastcore.geodesic_mst_mesh, which explains
why this never builds the k x k distance matrix the question seems to call
for. Always undirected — a minimum spanning tree of a directed graph is a
different problem (an arborescence) with a different algorithm.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
nodes
|
TYPE:
|
weights
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
edges
|
Rows of positions in
TYPE:
|
weights
|
Geodesic distance across each of those edges, in the resolved
TYPE:
|
Examples:
Two paths joined at their middle. Spanning the four endpoints costs three edges, each two hops long:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [1, 3], [3, 4]], dtype=np.uint32)
>>> mst, weights = fastcore.geodesic_mst_graph(edges, 5, nodes=[0, 2, 4])
>>> np.asarray([0, 2, 4])[mst]
array([[0, 2],
[0, 4]])
>>> weights
array([2., 3.], dtype=float32)
Nodes in different components are never joined, so the result is a forest:
Paths, not just distances¶
geodesic_matrix_graph answers how far;
geodesic_path and
geodesic_predecessors answer which way.
The motivating case is TEASAR-style skeletonization, which extracts a path, zeroes the edge weights along it so it is free to re-traverse, and searches again — so the graph changes between every call. That is why these take a bare edge list: there is no index to build, and nothing to invalidate when the weights move.
import navis_fastcore as fastcore
import numpy as np
edges, lengths = fastcore.unique_edges(faces, vertices)
weights = lengths.astype(np.float32)
# The route from the root to the farthest vertex...
dists, _ = fastcore.geodesic_predecessors(edges, n, weights, sources=[root])
farthest = int(np.argmax(dists[0]))
(path,) = fastcore.geodesic_path(edges, n, root, [farthest], weights=weights)
# ...and make it free to walk again
on_path = np.isin(edges, path).all(axis=1)
weights[on_path] = 0
Zero weights are explicitly supported. Among equal-length paths the route is picked deterministically, so repeated runs give the same skeleton.
navis_fastcore.geodesic_path(edges, n_nodes, source, targets, weights=None, directed=False)
¶
Node sequences of the shortest paths from source to each target.
The convenience form of :func:~navis_fastcore.geodesic_predecessors for the
common single-source case: one search, with the predecessor chains walked in
Rust rather than in Python. Because every target is known up front the search
also stops as soon as the last of them settles, so a short path in a large graph
costs a ball, not a sweep.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
source
|
TYPE:
|
targets
|
TYPE:
|
weights
|
TYPE:
|
directed
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
paths
|
One per target, ordered source-first / target-last (so
TYPE:
|
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 0]], dtype=np.uint32)
>>> weights = np.array([1, 1, 5], dtype=np.float32)
>>> fastcore.geodesic_path(edges, 3, 0, [2], weights=weights)
[array([0, 1, 2], dtype=uint32)]
An unreachable target gives an empty path:
navis_fastcore.geodesic_predecessors(edges, n_nodes, weights=None, directed=False, sources=None, limit=None, threads=None, dtype=None)
¶
Shortest path tree(s) - distances and the route to each node.
The predecessor-returning counterpart to
:func:~navis_fastcore.geodesic_matrix_graph. Use this when you need the path
itself; use geodesic_matrix_graph when the distance is enough, and
:func:~navis_fastcore.geodesic_path when you want the node sequences rather
than the raw chains.
Because this takes a bare edge list there is no index to build or invalidate between calls, which is what algorithms that re-weight the graph every iteration (TEASAR zeroes the edges along each path it extracts, then searches again) need.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
weights
|
TYPE:
|
directed
|
TYPE:
|
sources
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
dtype
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
distances
|
In the resolved dtype. As
TYPE:
|
predecessors
|
For each node, the node before it on the shortest path back to that
row's source.
TYPE:
|
Notes
Among equal-length paths the predecessor is the one reached first in the
search's own deterministic order, so results are reproducible run to run and do
not depend on threads.
Examples:
A triangle whose direct 0-2 edge is expensive, so the shortest path to 2 goes the long way round via 1:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 0]], dtype=np.uint32)
>>> weights = np.array([1, 1, 5], dtype=np.float32)
>>> dists, pred = fastcore.geodesic_predecessors(
... edges, 3, weights=weights, sources=[0]
... )
>>> dists
array([[0., 1., 2.]], dtype=float32)
>>> pred
array([[-1, 0, 1]], dtype=int32)
Clustering by geodesic radius¶
geodesic_clusters partitions a graph into
connected clusters of bounded radius: pick an unassigned seed, absorb everything within
max_dist of it that no earlier cluster claimed, repeat. Collapsing each cluster to its
centroid is a downsampling step — vertices come out spaced by roughly max_dist.
labels, n_clusters = fastcore.geodesic_clusters(edges, n, max_dist=2.0, weights=weights)
# Contiguous labels, so the centroids are one bincount away
centers = np.zeros((n_clusters, 3))
np.add.at(centers, labels, vertices)
centers /= np.bincount(labels, minlength=n_clusters)[:, None]
# ...and the coarse graph is just the contracted edge list
coarse = fastcore.contract_vertices(edges, labels.astype(np.uint32))
The radius is the true geodesic distance from the seed, not the length of the walk that reached it. The usual Python implementation of this is a recursive depth-first walk that accumulates distance along its own traversal path, which both gives worse clusters (a node close to a seed is dropped because the walk arrived the long way round) and recurses as deep as the cluster is large.
navis_fastcore.geodesic_clusters(edges, n_nodes, max_dist, weights=None, seeds=None)
¶
Greedily partition nodes into connected clusters of bounded radius.
Repeatedly takes an unassigned node as a seed and grows a cluster outwards from
it, absorbing any node reachable within max_dist that no earlier cluster has
already claimed. Collapsing each cluster to its centroid gives a coarser graph
whose nodes are spaced by roughly max_dist, which is what makes this useful
as mesh or skeleton downsampling.
The radius is the true geodesic distance from the seed, not the length of the walk that happened to reach it - so a node close to a seed is never excluded merely because a traversal arrived at it the long way round.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
max_dist
|
TYPE:
|
weights
|
TYPE:
|
seeds
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
labels
|
Cluster of each node, contiguous in
TYPE:
|
n_clusters
|
TYPE:
|
Notes
The greedy outer loop is inherently sequential - cluster n depends on
everything every earlier cluster claimed - so there is no threads argument.
Examples:
A path 0-1-...-5 with a radius of one hop:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5]], dtype=np.uint32)
>>> labels, n = fastcore.geodesic_clusters(edges, 6, 1)
>>> labels
array([0, 0, 1, 1, 2, 2], dtype=int32)
>>> n
3
Seeding from the middle instead:
Reusing a graph across many queries¶
Every function above takes an edge list and builds an adjacency index from it, answers one
question, and throws the index away. GeodesicGraph is the
same functionality with that build hoisted out. It carries the methods you would expect —
| method | equivalent to |
|---|---|
.distances() |
geodesic_matrix_graph |
.nearest() / .farthest() |
geodesic_nearest_mesh / geodesic_farthest_mesh |
.predecessors() / .path() |
geodesic_predecessors / geodesic_path |
.clusters() |
geodesic_clusters |
.components() |
connected_components_graph |
— each answering exactly what its counterpart does, so migrating is a mechanical change:
g = fastcore.GeodesicGraph(edges, n, weights=weights)
routes = [g.path(a, [b]) for a, b in pairs] # one index build, not one per pair
When this actually pays. The build is O(E) over the whole graph, so hoisting it out is
worth real time exactly when each query is small relative to the graph — many short paths,
a nearest with a tight limit, grow. On a 40k-vertex mesh, 500 short-path queries run
~100x faster as methods than as free-function calls. It buys nothing measurable when a
single query already sweeps the graph: one 50-source distance matrix on that same mesh takes
90 ms either way, against a 1 ms build. Reach for the class because you have a graph and
want to stop re-passing edges/weights/directed to everything — and take the speedup
where the query pattern happens to earn it.
Two further methods, grow and
farthest_seed, have no free-function
counterpart at all — they are the ones that are inherently called in a loop, and are
covered below.
.subset(nodes) carves out an induced subgraph without ever returning to your edge list —
masking and renumbering an edge list in numpy is both slower and easy to get subtly wrong:
labels = g.components()
biggest = g.subset(labels == np.bincount(labels).argmax())
biggest.parent_nodes # which original nodes these are
Note that distances inside a subset are not generally the parent's — a shortest path that left the subset is gone. Taking a whole connected component, as above, is the case where they agree.
Growing fixed-size regions¶
geodesic_clusters fixes each cluster's radius.
When what has to be fixed is its size — tiling a neuron into equal-length inputs for a
neural network, say — use GeodesicGraph instead. It
grows outwards from a seed and stops once it has gathered the requested number of
points, so each region is the geodesic ball that happens to hold exactly that many.
This one is a class rather than a function because the calling pattern is different: a tiling driver asks many small questions of one graph, and building the adjacency index per call — O(E) over the whole graph, against a query that explores a handful of nodes — would cost far more than the searches themselves. Build once, query in a loop:
g = fastcore.GeodesicGraph(edges, n, weights=weights)
claimed = np.zeros(n, dtype=bool)
fragments = []
while not claimed.all():
frag = g.grow(int(np.argmax(~claimed)), size=64, forbidden=claimed)
claimed[frag] = True
fragments.append(frag)
Feeding claimed back in as forbidden is what makes the fragments a genuine partition:
an already-claimed point is never collected again, and a node whose points are all
claimed becomes a wall, so a later fragment cannot tunnel through an earlier one and
come out somewhere unrelated. Every fragment is therefore a single connected piece.
Pass return_distances=True and each point's distance to the seed comes back alongside.
The search settles points in distance order, so it already holds the number and returning it
is free. It is what makes a patch non-uniform: thinning a grown region by radius — a dense
core with a sparse, far-reaching halo, so one point budget buys both local detail and
long-range context — needs to know each point's radius.
idx, dist = g.grow(seed, size=8 * 1024, return_distances=True) # oversized candidate pool
keep = thin_by_radius(dist) # your falloff of choice
patch = idx[keep]
Points sharing a node share a distance exactly, since a point's position is its node's, so a thinning keyed on these will not drift.
Points do not have to be the graph's nodes. Pass item_nodes to attach a cloud —
a resampled surface, say — to the graph, and growth counts and returns cloud points
while still travelling along the graph. Nodes carrying no point simply conduct, which is
what keeps a patch connected when the cloud is far sparser than the mesh beneath it:
# `source_id[i]` is the mesh vertex that cloud point `i` was sampled from
g = fastcore.GeodesicGraph(edges, n_vertices, weights=lengths, item_nodes=source_id)
patch = g.grow(seed, size=1024) # 1024 *cloud points*, geodesically connected
Spreading seeds evenly¶
Growing regions from arbitrary seeds clumps them. GeodesicGraph.farthest_seed
picks the point geodesically farthest from everything chosen so far — farthest-point
sampling — so a sequence of them tiles the graph evenly instead:
g = fastcore.GeodesicGraph(edges, n, weights=weights)
chosen = np.zeros(n, dtype=bool)
patches = []
for _ in range(k):
seed = g.farthest_seed(chosen)
if seed is None: # nothing left to seed
break
patch = g.grow(seed, size=64)
chosen[patch] = True # or just `chosen[seed] = True` to overlap less
patches.append(patch)
Only points reachable from something already chosen are candidates, and the search only jumps to a fresh component (largest first) once the reachable frontier is exhausted. That rule is what stops a mesh with a few hundred disconnected specks from seeding every speck before it returns to the main body.
The cost is what makes this usable at scale. The obvious implementation re-runs a
multi-source Dijkstra over the whole graph per seed, which is quadratic in the seed count —
placing 2560 seeds on a 160k-vertex mesh takes ~93 s through
scipy.sparse.csgraph.dijkstra(..., min_only=True). Here the distance field is updated
incrementally and the update is pruned against the running field, so each fold costs only
the region the new seed actually claims; the same 2560 seeds take ~0.35 s.
This assumes chosen only ever grows between calls, which is what the loop above does. It
is allowed to shrink — the field is rebuilt when that is detected, so the answer stays
correct — but that call pays a cold start.
Sweeping outwards from a set¶
GeodesicGraph.ball answers "what is within r of
any of these nodes, how far, and which one is nearest" — one multi-source search rather
than one search per source:
The query is scipy.sparse.csgraph.dijkstra(..., min_only=True, limit=r), and the
difference is what comes back: the ball itself, rather than three node-sized arrays with
the ball buried in them. For a radius that covers a fraction of the graph, allocating and
filling those arrays costs more than the search does — 275 µs against 5 µs on a 26k-vertex
mesh at a radius reaching ~170 nodes.
That gap is the point, because the callers that want this ask thousands of times: walking a
mesh invalidating a neighbourhood at a time, growing regions and marking what they covered.
Ties — a node equidistant from two sources — go to whichever settled first, deterministic
but otherwise arbitrary, exactly as min_only is elsewhere.
Changing the graph as you go¶
Some algorithms re-weight as they run. TEASAR is the standard example: it extracts the longest path, zeroes it so later paths may re-traverse it for free, and repeats — which is what makes it pick branches rather than re-walk the trunk.
set_weights edits those arcs in place:
g = fastcore.GeodesicGraph(edges, n, weights=lengths)
...
path = g.path(root, [farthest])[0]
g.set_weights(np.column_stack((path[:-1], path[1:])), 0.0) # free to re-traverse
Rebuilding the graph to change a few hundred edges costs O(E) and gives up the reason for holding a prepared graph at all — 728 µs against 2.6 µs for the edit, on that same mesh. It cannot add an edge, only re-weight one the graph already has, since growing the adjacency would mean rebuilding it; a pair that is not an edge raises rather than being ignored, an edge list out of step with the graph being the likely cause.
Distances change, so the incremental field behind farthest_seed is discarded on each edit
— interleaving the two costs a cold start per edit. Component labels survive.
navis_fastcore.GeodesicGraph
¶
A graph prepared once for many geodesic queries.
The module-level geodesic functions each build an adjacency index from your edge list, answer one question and throw it away. That is the right trade for a single all-pairs sweep, and the wrong one for algorithms that ask many small questions of the same graph - the index build is O(E) over the whole graph, so it dwarfs a query that only ever explores a small ball.
This class is for that second pattern. It builds the index once, keeps the search scratch space alive between calls, and lets each query cost only the ball it explores.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
n_nodes
|
TYPE:
|
weights
|
TYPE:
|
directed
|
TYPE:
|
item_nodes
|
TYPE:
|
| ATTRIBUTE | DESCRIPTION |
|---|---|
n_nodes |
TYPE:
|
n_items |
Equals
TYPE:
|
item_nodes |
The node each item sits on.
TYPE:
|
Notes
What's here. Most methods are the module-level functions with the index build
taken out, and answer exactly what their counterpart does: :meth:distances
(:func:~navis_fastcore.geodesic_matrix_graph), :meth:nearest, :meth:farthest,
:meth:predecessors, :meth:path, :meth:clusters and :meth:components. The two
with no counterpart - :meth:grow and :meth:farthest_seed - are the ones that only
make sense against a graph you keep, since they are called in a loop. :meth:subset
carves out an induced subgraph without going back to your edge list.
When the reuse pays. Hoisting the index build out of the call is worth real time
exactly when each query is small relative to the graph - many short paths, a
:meth:nearest with a tight limit, :meth:grow. On a 40k-vertex mesh, 500
short-path queries run ~100x faster through this class than through
:func:~navis_fastcore.geodesic_path. It buys nothing measurable when a single query
already sweeps the graph: one 50-source distance matrix on that mesh costs 90 ms
either way, against a 1 ms build. The other reason to reach for the class is simply
that you have a graph, and would rather not re-pass edges/weights/directed
to every call.
Width. float32 only, and the one place the module's "your dtype in, your dtype
out" rule does not apply: float64 weights are accepted but narrowed, and every
distance this class returns is float32. That is deliberate rather than an omission.
The class exists for "large graph, many small queries", which is exactly the case
where float32 is the right width and where doubling the several node-sized arrays it
holds resident for a whole run would be felt. If you need float64, use the module-level
functions, which rebuild the index per call and take a dtype.
Items. Optionally each node carries zero or more items - points of a cloud attached to the graph, one entry of a resampled surface say. By default the distinction vanishes entirely.
The rule for which index space a method speaks is short: :meth:grow,
:meth:farthest_seed and :meth:item_components count and return items;
everything else takes and returns graph nodes, exactly as the free function it
mirrors does. So growth follows the graph but is measured in cloud points - which is
what keeps a patch of a cloud far sparser than its mesh connected, since the empty
nodes in between conduct without contributing - while a distance matrix stays a
matrix over the graph.
Direction. With directed=True, :meth:grow gathers the out-reachable ball,
:meth:farthest_seed measures distance from the done set, and :meth:clusters
grows out-balls (so it differs from :func:~navis_fastcore.geodesic_clusters, which
is always undirected). :meth:components still reports weakly connected
components - a search has to start somewhere.
Threading. Queries share mutable scratch space, so concurrent calls on one
instance serialise. Build one instance per thread if you need real parallelism. The
threads argument on the matrix-style methods is unaffected: those parallelise
internally over sources.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5]], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 6)
>>> g
GeodesicGraph(n_nodes=6, n_items=6)
Grow a ball of three nodes around node 2 - nearest first:
Tile the whole graph into disjoint connected fragments by feeding what has already
been claimed back in as forbidden:
>>> claimed = np.zeros(6, dtype=bool)
>>> frags = []
>>> while not claimed.all():
... frag = g.grow(int(np.argmax(~claimed)), 2, forbidden=claimed)
... claimed[frag] = True
... frags.append(frag.tolist())
>>> frags
[[0, 1], [2, 3], [4, 5]]
item_nodes
property
¶
The node each item sits on, as a (n_items, ) uint32 array.
ball(sources, max_dist=None)
¶
Every node within max_dist of any of sources, and its nearest source.
One multi-source search, so this costs the ball it returns rather than one sweep
per source - and it returns the ball itself, not a (n_nodes, ) array with the
ball buried in it. Both halves matter for the pattern this exists for: sweeping a
graph a neighbourhood at a time, thousands of small radii against one big graph.
The nearest equivalent elsewhere,
scipy.sparse.csgraph.dijkstra(..., min_only=True, limit=...), allocates and
fills three node-sized arrays per call whatever the radius, which for a small one
costs more than the search does.
| PARAMETER | DESCRIPTION |
|---|---|
sources
|
TYPE:
|
max_dist
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
nodes
|
The nodes within reach, in increasing-distance order. Every source is
in here at distance 0. Nodes farther than
TYPE:
|
distances
|
Distance from each node to its nearest source.
TYPE:
|
sources
|
That source, as a node index. A source's own entry is itself. Ties are broken deterministically but arbitrarily.
TYPE:
|
Examples:
A path of 7 nodes, seeded from both ends - one hop out from each:
>>> edges = np.array([[i, i + 1] for i in range(6)], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 7)
>>> nodes, dist, src = g.ball([0, 6], 1)
>>> order = np.argsort(nodes) # settle order interleaves the two frontiers
>>> nodes[order]
array([0, 1, 5, 6], dtype=uint32)
>>> dist[order]
array([0., 1., 1., 0.], dtype=float32)
>>> src[order]
array([0, 0, 6, 6], dtype=uint32)
clusters(max_dist, seeds=None)
¶
Greedily partition nodes into connected clusters of bounded radius.
As :func:~navis_fastcore.geodesic_clusters, minus the adjacency build. This is
the radius-bounded sibling of :meth:grow's count-bounded ball: use this when
the cluster's physical extent is what must be fixed, :meth:grow when its size
is.
On a graph built with directed=True this grows out-balls, so it differs
from the always-undirected free function.
| PARAMETER | DESCRIPTION |
|---|---|
max_dist
|
TYPE:
|
seeds
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
labels
|
TYPE:
|
n_clusters
|
TYPE:
|
components()
¶
Component label of each node.
Labels are node indices - the smallest node index in the component, the same
convention :func:~navis_fastcore.connected_components_graph uses. See
:meth:item_components for the per-item view.
| RETURNS | DESCRIPTION |
|---|---|
labels
|
TYPE:
|
distances(sources=None, targets=None, limit=None, threads=None)
¶
Pairwise geodesic distances between sources and targets (nodes).
The same query as :func:~navis_fastcore.geodesic_matrix_graph, minus the
adjacency build. Use this when you slice the same graph repeatedly - a batch of
sources at a time, say - where the free function would rebuild its index on each
call.
| PARAMETER | DESCRIPTION |
|---|---|
sources
|
TYPE:
|
targets
|
TYPE:
|
limit
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
matrix
|
TYPE:
|
Examples:
farthest(sources=None, targets=None, limit=None, threads=None)
¶
For each source node, the distance to its farthest target and which one.
The mirror of :meth:nearest; see
:func:~navis_fastcore.geodesic_farthest_mesh. This one cannot stop early - it
has to settle every target - but the farthest is then free, since the kernels
settle in increasing distance order.
| RETURNS | DESCRIPTION |
|---|---|
distances
|
TYPE:
|
indices
|
TYPE:
|
farthest_seed(done)
¶
The undone item geodesically farthest from everything already done.
Calling this repeatedly with a growing done set is farthest-point sampling:
it spreads seeds evenly over the graph instead of letting them clump, which is
what you want when placing patches, landmarks or cluster centres. Pair it with
:meth:grow - seed, grow, mark what you covered, seed again.
Only items reachable from something in done are candidates. Unreachable
ones are infinitely far and would otherwise win every time, so a mesh with a few
hundred disconnected specks would seed every speck before returning to the main
body. Once the reachable frontier is exhausted - or when done is empty - this
jumps to a fresh component, largest first. Ties go to the lower item index.
| PARAMETER | DESCRIPTION |
|---|---|
done
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
seed
|
TYPE:
|
Notes
The distance field is maintained incrementally, and each update is pruned
against the running field, so a call costs the region the new sources actually
claim rather than a sweep of the whole graph. The usual implementation of this
(a fresh multi-source Dijkstra per seed, e.g. via
scipy.sparse.csgraph.dijkstra(..., min_only=True)) cannot prune that way and
so goes quadratic in the number of seeds: placing 2560 seeds on a 160k-vertex
mesh takes ~93 s that way against ~0.35 s here.
done is expected to only ever grow between calls. It may shrink - the field
is rebuilt from scratch when that is detected, so the answer stays correct - but
that call loses the incremental saving.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[i, i + 1] for i in range(8)], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 9)
>>> done = np.zeros(9, dtype=bool)
>>> done[0] = True
>>> picks = []
>>> for _ in range(4):
... s = g.farthest_seed(done)
... picks.append(s)
... done[s] = True
>>> picks # the far end, then repeated bisection
[8, 4, 2, 6]
grow(seed, size, forbidden=None, return_distances=False)
¶
Grow a connected region of up to size items outwards from seed.
Settles nodes in order of increasing geodesic distance from the seed's node,
collecting the items on each until size are gathered. The region is
therefore the geodesic ball around the seed that happens to hold size
items - not whatever a depth-first walk stumbled into - and it is always
connected, since every node reached bar the seed's own is reached through one
settled before it.
This is the count-bounded sibling of
:func:~navis_fastcore.geodesic_clusters, which bounds by radius instead. Use
this one when the fragment size is what must be fixed, e.g. tiling a neuron
into equal-length inputs for a neural network.
| PARAMETER | DESCRIPTION |
|---|---|
seed
|
TYPE:
|
size
|
TYPE:
|
forbidden
|
TYPE:
|
return_distances
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
region
|
Item indices, seed-first and in increasing-distance order.
TYPE:
|
distances
|
Only if
TYPE:
|
Examples:
A path of 7 nodes carrying a sparse cloud - two points at each end, nothing in between. The empty nodes conduct, so the patch spans them:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[i, i + 1] for i in range(6)], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 7, item_nodes=[0, 0, 6, 6])
>>> g.n_items
4
>>> g.grow(0, 4)
array([0, 1, 2, 3], dtype=uint32)
The distances come back alongside on request - here both points at each end share their node's distance:
item_components()
¶
Component label of each item.
Labels are node indices - specifically the smallest node index in the component,
the same convention
:func:~navis_fastcore.connected_components_graph uses - not a contiguous range.
:meth:farthest_seed deliberately does not offer a "random seed from the largest
component"; that would mean owning a random number generator, and a caller who
cares about reproducibility wants it to be theirs. This is the piece to build it
from:
import navis_fastcore as fastcore import numpy as np edges = np.array([[0, 1], [2, 3], [3, 4]], dtype=np.uint32) g = fastcore.GeodesicGraph(edges, 5) labels = g.item_components() labels array([0, 0, 2, 2, 2], dtype=uint32) pool = np.flatnonzero(labels == np.bincount(labels).argmax()) int(np.random.default_rng(0).choice(pool)) # a seed off the largest component 4
| RETURNS | DESCRIPTION |
|---|---|
labels
|
TYPE:
|
nearest(sources=None, targets=None, limit=None, threads=None)
¶
For each source node, the distance to its nearest target and which one.
As :func:~navis_fastcore.geodesic_nearest_mesh, minus the adjacency build.
O(sources) output instead of O(sources x targets), and faster than the
matrix too - each search stops at the first target it settles.
A source that is itself a target is matched to its nearest distinct target.
Sources with no reachable distinct target get -1 / -1.
| PARAMETER | DESCRIPTION |
|---|---|
sources
|
DEFAULT:
|
targets
|
DEFAULT:
|
limit
|
DEFAULT:
|
threads
|
DEFAULT:
|
| RETURNS | DESCRIPTION |
|---|---|
distances
|
TYPE:
|
indices
|
Index into the graph, not into
TYPE:
|
path(source, targets)
¶
Node sequences of the shortest paths from source to each of targets.
As :func:~navis_fastcore.geodesic_path, minus the adjacency build. One search,
stopped as soon as the last target settles.
| PARAMETER | DESCRIPTION |
|---|---|
source
|
TYPE:
|
targets
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
paths
|
One per target, source-first and target-last. An unreachable target
gives an empty array; a target equal to
TYPE:
|
Examples:
predecessors(sources=None, limit=None, threads=None)
¶
Shortest-path trees: distances and the route to every node.
As :func:~navis_fastcore.geodesic_predecessors, minus the adjacency build.
Use :meth:path when you want the node sequences rather than the raw chains.
| PARAMETER | DESCRIPTION |
|---|---|
sources
|
TYPE:
|
limit
|
DEFAULT:
|
threads
|
DEFAULT:
|
| RETURNS | DESCRIPTION |
|---|---|
distances
|
TYPE:
|
predecessors
|
The node before each node on its shortest path back to that row's
source;
TYPE:
|
set_weights(edges, weights)
¶
Re-weight edges in place, leaving the graph otherwise untouched.
For algorithms that change the graph as they run - TEASAR zeroing each path it extracts so later ones may re-traverse it for free is the case this was added for. Rebuilding the graph after each change costs O(E) against an edit of a few hundred edges, and gives up the reason for holding a prepared graph in the first place; this costs O(edits) and a binary search apiece.
| PARAMETER | DESCRIPTION |
|---|---|
edges
|
TYPE:
|
weights
|
TYPE:
|
Notes
Only available on a graph built with weights: there is no weight array on an unweighted one to write into, and materialising one would quietly turn every later search from a BFS into a Dijkstra.
Distances change, so the incremental field behind :meth:farthest_seed is
discarded here - a minimum folded under the old weights cannot be corrected under
the new ones. Interleaving farthest_seed with re-weighting therefore pays a cold
start after each edit; not interleaving them costs nothing. Component labels
survive, since which edges exist has not changed.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [1, 2], [0, 2]], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 3, weights=[1.0, 1.0, 5.0])
>>> g.distances(sources=[0], targets=[2])
array([[2.]], dtype=float32)
>>> g.set_weights([[0, 2]], [0.5])
>>> g.distances(sources=[0], targets=[2])
array([[0.5]], dtype=float32)
subset(nodes)
¶
The subgraph induced on nodes, as a graph in its own right.
New node i is old node nodes[i]. Edges with an endpoint outside the
subset are dropped, and items go wherever their node did - an item whose node was
dropped goes with it. The result carries parent_nodes and parent_items
so anything computed on it can be mapped back.
The subgraph is carved out of the adjacency already built rather than re-derived, so this never returns to your original edge list - which is the point, since masking and renumbering an edge list in numpy is both slower and easy to get wrong. Restricting to one connected component is the motivating case.
Distances within a subset are not generally the parent's: a shortest path that left the subset is gone. Taking a whole connected component is the case where they do agree.
| PARAMETER | DESCRIPTION |
|---|---|
nodes
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
subgraph
|
TYPE:
|
Examples:
Pull out the largest connected component:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> edges = np.array([[0, 1], [2, 3], [3, 4]], dtype=np.uint32)
>>> g = fastcore.GeodesicGraph(edges, 5)
>>> labels = g.components()
>>> sub = g.subset(labels == np.bincount(labels).argmax())
>>> sub.n_nodes
3
>>> sub.parent_nodes # which original nodes these are
array([2, 3, 4], dtype=uint32)
>>> sub.distances(sources=[0], targets=[2])
array([[2.]], dtype=float32)
Simplification¶
Every mesh simplifier gives you a smaller mesh. The problem is what happens to the data
you had attached to the old vertices — synapses, radii, compartment labels. Decimating
a neuron mesh with pyfqmr or meshopt orphans all of it, and re-attaching by nearest
neighbour afterwards is both slower and wrong: an edge collapse moves its survivor to the
quadric-optimal point, which is frequently nearer some other vertex than the one that
actually merged into it.
simplify_mesh returns that correspondence as a third
array. vertex_map[i] is the simplified vertex that original vertex i ended up in, or
-1 if it did not survive:
verts, faces, vmap = fastcore.simplify_mesh(faces, vertices, ratio=0.1)
# Push a per-vertex quantity onto the simplified mesh. `bincount` over the map is the
# whole operation -- no spatial query, no tolerance to pick.
live = vmap >= 0
counts = np.bincount(vmap[live], weights=synapses[live], minlength=len(verts))
The map is the forward direction — indexed by original vertex, valued in simplified
vertices — because that is the direction aggregation needs. It is int32 rather than
uint32 for the -1.
Pinning vertices¶
Positions move under decimation, so a vertex that carried a synapse is no longer exactly
where the synapse was. Where that matters, lock freezes it:
lock = np.zeros(len(vertices), dtype=bool)
lock[synapse_vertices] = True
verts, faces, vmap = fastcore.simplify_mesh(faces, vertices, ratio=0.1, lock=lock)
assert np.array_equal(verts[vmap[synapse_vertices]], vertices[synapse_vertices])
A locked vertex is never merged into another and never moved — the equality above is bitwise, not approximate. It may still absorb its neighbours, which is what keeps the face target reachable when the pinned set is large; freezing a vertex's whole one-ring instead would stall the sweep as soon as you pinned a few thousand synapses. The floor is that every locked vertex survives, so a target below the locked count cannot be met.
Lossless¶
simplify_mesh_lossless collapses only edges
whose quadric error is under epsilon and runs to a fixed point. It has no face budget:
it is for shedding over-tessellation — coplanar fans, duplicate vertices, degenerate
faces — rather than hitting a target.
"Lossless" is a claim about the surface, not the outline. A quadric measures distance
to the planes of the incident faces, and the plane of a flat patch says nothing about
where that patch ends, so on an open mesh a planar region will collapse its own boundary
inwards at zero measured cost. Pass preserve_border=True there.
Notes¶
Non-manifold input is fine. Nothing here checks for manifoldness, and each collapse guard skips what it cannot handle rather than failing. This is the reason the algorithm is implemented here rather than wrapped from a crate: everything built on a halfedge or corner table either refuses a mesh with an edge shared by three faces or silently drops the offending faces, and meshes out of EM segmentation are full of them.
It is the same algorithm pyfqmr runs — a port of Sven Forstmann's Simplify.h — so
expect comparable speed rather than a speed-up. On a clean mesh the two agree to the face
array, which is what the test suite checks. What is new is the vertex map, the pinning,
and that no C++ toolchain is involved, which is what lets the same code build for
pyodide and for the R source tarball.
Determinism. Same input, same output, every run and every machine — the sweep is
single-threaded and index-ordered. The result does depend on the order of faces, since
triangles are visited in input order; that is normal for this family of algorithm and not
a determinism failure.
navis_fastcore.simplify_mesh(faces, vertices, ratio=None, n_faces=None, aggressiveness=7.0, preserve_border=False, lock=None)
¶
Simplify a triangle mesh, tracking where every vertex went.
Iteratively contracts the edge whose collapse costs least, where the cost is the
Garland-Heckbert quadric error: the summed squared distance from the merged
vertex to the planes of every face that met at the two it replaces. This is the
algorithm pyfqmr runs — a port of the same MIT-licensed original — with the
one thing no implementation in this space returns, vertex_map.
Non-manifold input is fine. Meshes out of EM segmentation routinely have edges shared by three faces, bowtie vertices and zero-area triangles; nothing here checks for manifoldness, and each collapse guard skips what it cannot handle rather than failing.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
ratio
|
TYPE:
|
n_faces
|
TYPE:
|
aggressiveness
|
TYPE:
|
preserve_border
|
TYPE:
|
lock
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
vertices
|
Positions of the surviving vertices.
TYPE:
|
faces
|
Faces, indexing the returned
TYPE:
|
vertex_map
|
For each input vertex, the index of the output vertex it
ended up in;
TYPE:
|
Notes
When a vertex maps to -1. Being merged is not one of those cases —
that is what the map is for, and a collapsed vertex points at whatever it
merged into. Decimating a clean closed mesh to 5% of its faces yields no -1
at all. The rule is that vertex_map[i] is -1 exactly when the vertex
i ended up in is referenced by no surviving face, which happens in four
situations:
iis in no face to begin with, so it never takes part in a collapse.iappears only in zero-area faces. A face naming the same vertex twice carries no plane and no normal, so it is dropped on the way in, which reduces this to case 1.- The whole piece
ibelonged to was decimated away. A collapse deletes exactly the faces holding both its endpoints, so a survivor normally keeps faces — but the budget is global, with nothing reserved per component, so a small disconnected fragment is consumed entirely once the target is tight enough. Simplify per component if that matters. - The input is wholly degenerate, so the output mesh is empty.
Mask with vertex_map >= 0 before aggregating (see Examples).
Positions move. vertices_out[vertex_map[i]] is not
vertices_in[i]: a collapse moves its survivor to the quadric-optimal point,
so any vertex that took part in one has shifted. Use lock where the exact
position matters.
What lock guarantees is that a locked vertex is never merged into
another and never moved. It is not an absolute guarantee against -1: it
cannot conjure up a vertex that no surviving face references, so cases 1, 2 and
4 above still apply, and a locked vertex can in principle lose every face it sat
on if each of them happened to hold both endpoints of some other collapse.
Locking does set a floor on how far the mesh can shrink — every locked vertex
survives — so a target below the locked count is not reachable.
Determinism. The result depends on the order of faces, since the sweep
visits them in input order, but is otherwise deterministic: same input, same
output, every run. There is no threads argument because each collapse
invalidates its own neighbourhood, so the sweep cannot be parallelised without
changing the answer. The GIL is released for the duration, so simplifying
several meshes from a thread pool does scale.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> # A unit square as two triangles, subdivided once.
>>> faces = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.uint32)
>>> vertices = np.array([[0., 0., 0.], [1., 0., 0.],
... [0., 1., 0.], [1., 1., 0.]])
>>> v, f, vmap = fastcore.simplify_mesh(faces, vertices, n_faces=2)
>>> len(f)
2
>>> vmap
array([0, 1, 2, 3], dtype=int32)
Push a per-vertex quantity — synapse counts, say — onto the simplified mesh:
>>> syn = np.array([3, 0, 1, 2])
>>> live = vmap >= 0
>>> np.bincount(vmap[live], weights=syn[live], minlength=len(v))
array([3., 0., 1., 2.])
Pin the vertices that carry synapses so they keep their exact positions:
navis_fastcore.simplify_mesh_lossless(faces, vertices, epsilon=0.001, max_iterations=9999, preserve_border=False, lock=None)
¶
Simplify a triangle mesh without changing its shape.
Collapses only edges whose quadric error is below epsilon and repeats until
a whole pass changes nothing. There is no face budget: this is for shedding
over-tessellation — coplanar fans, duplicated vertices, degenerate faces —
rather than for hitting a target. Use :func:simplify_mesh for that.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
epsilon
|
TYPE:
|
max_iterations
|
TYPE:
|
preserve_border
|
TYPE:
|
lock
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
vertices
|
TYPE:
|
faces
|
TYPE:
|
vertex_map
|
As :func:
TYPE:
|
Notes
"Lossless" is a claim about the surface, not the outline. A quadric measures
distance to the planes of the incident faces, and the plane of a flat patch says
nothing about where that patch ends — so with preserve_border=False a planar
region will happily collapse its own boundary inwards at zero measured cost. Pass
preserve_border=True on open meshes.
The circumstances under which a vertex maps to -1 rather than to the vertex
it merged into, and what lock does and does not guarantee, are as
:func:simplify_mesh. Case 3 there — a whole piece consumed — arrives by a
different route here: there is no face budget, but epsilon is an absolute
error, so a component small enough that all of its edges fall under it collapses
away entirely. A lone triangle does this at any size; a tetrahedron 0.001 across
does it at the default epsilon, whether or not a larger mesh surrounds it.
Examples:
A flat 4x4 grid of vertices. The four interior ones are exactly coplanar with their neighbours, so removing them costs nothing; the rim is held in place.
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> vertices = np.array([[i, j, 0] for i in range(4) for j in range(4)],
... dtype=float)
>>> faces = np.array(
... [t for i in range(3) for j in range(3)
... for t in ([i * 4 + j, (i + 1) * 4 + j, (i + 1) * 4 + j + 1],
... [i * 4 + j, (i + 1) * 4 + j + 1, i * 4 + j + 1])],
... dtype=np.uint32,
... )
>>> v, f, vmap = fastcore.simplify_mesh_lossless(faces, vertices,
... preserve_border=True)
>>> len(faces), len(f)
(18, 12)
>>> bool(np.abs(v[:, 2]).max() < 1e-9) # still flat
True
The four interior vertices — 5, 6, 9 and 10 — merged into a single one:
Smoothing¶
The other half of mesh cleanup: moving vertices to take the noise out of a surface, without changing how many there are or which faces they form. The face array and the vertex order come back untouched, so anything you have indexed by vertex is still attached to the vertex it was attached to.
That default is Taubin's λ|μ filter, and it is the default because the obvious alternative is a trap.
Why not plain Laplacian¶
The plain Laplacian step — average each vertex with its neighbours, take a fraction of
the way — removes high frequencies quickly and low ones slowly. A closed surface's
enclosed volume is a low frequency, so it removes that too. At lamb=0.5 and five
iterations, which is what navis.smooth_mesh ships today, a neuron mesh comes out having
lost 88% of its volume.
Taubin alternates a shrinking λ pass with an inflating μ pass, tuned so the two cancel below a cut-off frequency and reinforce above it. Same fixture, twenty iterations: it holds its volume to within 5%.
# Explicitly, if you want the Laplacian anyway
smoothed = fastcore.smooth_mesh(faces, vertices, method="laplacian", lamb=0.5,
iterations=5)
One Taubin iteration is a full λ/μ pair
Two sweeps over the mesh, not one — trimesh.smoothing.filter_taubin counts
half-steps. Counting half-steps lets an odd iterations end on a λ pass, which is a
shrink that nothing undoes, and the whole point of the filter is that the passes come
in pairs. iterations=10 here equals iterations=20 there, and the two agree to
~1e-11.
method="humphrey" is the third option — the HC filter of Vollmer et al., which fights
shrinkage by pulling each vertex back towards where it started rather than towards a
lower frequency. It is the gentler of the two on fine detail worth keeping.
Weights¶
weights chooses how a vertex's one-ring is averaged. The default "uniform" counts
every neighbour equally, which also regularises the sampling: where the tessellation is
uneven it slides vertices along the surface towards even spacing. Sometimes that is what
you want; often it is drift you did not ask for.
"cotangent" is the discrete Laplace–Beltrami operator — each edge weighted by the
cotangents of the two angles opposite it. It is a function of the surface rather than of
the triangulation, so it moves vertices along the normal and leaves them alone within the
surface. On a UV sphere, whose rings crowd together at the poles, it drifts less than half
as far as the uniform umbrella for the same amount of smoothing.
Cotangents go negative on obtuse triangles, and a negative weight pushes a vertex away from its neighbour, so those contributions are clamped to zero — the usual remedy. A vertex whose weights all vanish that way falls back to the uniform umbrella. The cost is that on a surface of mostly-obtuse triangles cotangent weighting degrades towards uniform, which is the right way to fail.
Unlike trimesh, which builds its operator once from the input geometry and reuses it,
the geometry-dependent weightings here are recomputed from the current positions every
pass — the flow they are supposed to discretise rather than a snapshot taken before the
first step. It costs nothing, because the weights are never materialised at all: each is
derived as the one-ring is walked, which is about what reading it back from an array would
have cost anyway.
Boundaries and pinning¶
A boundary vertex's one-ring lies entirely to one side of it, so an open mesh's rim rolls
inwards under any of these filters. preserve_border=True pins it — a boundary vertex
being an endpoint of an edge used by exactly one face:
lock freezes an arbitrary set on top of that — the same name and the same meaning as
simplify_mesh's. A locked vertex comes back at bitwise
the same coordinates but still pulls on its neighbours, which is what makes it a boundary
condition rather than a hole:
lock = np.zeros(len(vertices), dtype=bool)
lock[synapse_vertices] = True
smoothed = fastcore.smooth_mesh(faces, vertices, lock=lock, preserve_border=True)
assert np.array_equal(smoothed[lock], vertices[lock])
Volume correction¶
volume_correction=True rescales the result about its centroid so the enclosed volume
matches the input's:
About the centroid is the one place this deliberately differs from
trimesh.smoothing.filter_laplacian, and the difference is not cosmetic. Upstream
rescales by (vol_before / vol_after) ** (1/3) about the origin, which is not a shape
operation:
- It translates the mesh. On the 722817260 test neuron at navis' own defaults, the constraint displaces the result by 41 µm. The mesh is 19–26 µm across.
- It is not translation invariant. The same mesh smoothed at two different offsets
comes out two different shapes; far enough from the origin the volume ratio goes
negative and the cube root returns
NaN. - It divides by the smoothed volume, so a mesh with a hole big enough to make that
zero is a
ZeroDivisionErrorrather than a diagnostic.
The correction here also runs once, at the end — which is not an approximation of
running it every iteration but exactly equal to it. Every filter is an affine combination
of a vertex and a normalised average of its neighbours, and those commute with a uniform
scaling, so scaling first and smoothing lands on the same vertices as smoothing and
scaling afterwards. Upstream pays a full pass over the faces and a (F, 3, 3) gather per
iteration — 40% of its runtime — for a result it could have had at the end for one pass.
When the volume is undefined
On a closed mesh the correction is exactly what it says. A mesh that is not closed still usually gets one, and deliberately: both measurements cone every face back to the same anchor, so their ratio stays a consistent measure of how much the surface shrank even where neither number is an enclosed volume on its own. That matters because meshes worth smoothing are almost never watertight — the 722817260 neuron is not.
What is left is the genuinely undecidable case: the ratio of the two signed volumes is
zero, infinite, NaN or negative, a flat sheet being the clean example. There the
vertices come back smoothed but unscaled and a RuntimeWarning says so. Consistently
inverted winding is not in that set — both volumes come out negative, the ratio is
positive, and the correction is as valid as ever.
Notes¶
Speed. On a 421k-vertex / 881k-face mesh, ten iterations with the volume correction:
trimesh.smoothing.filter_laplacian |
5.42 s |
fastcore.smooth_mesh, uniform |
0.03 s |
fastcore.smooth_mesh, cotangent |
0.06 s |
The arithmetic was never the cost. Upstream spends 57% of its time building the operator —
vertex_neighbors is a list of 421k Python lists, 636 MB of heap for 10 MB of vertices —
and another 40% in the volume constraint's per-iteration vertices[faces] gather. The
sparse matrix–vector product itself is 42 ms of the 5.4 s.
Non-manifold input is fine, as for simplification. An edge shared by three faces, a face naming the same vertex twice, a duplicated face and a vertex no face mentions are all merely data; nothing here reads more topology than "which vertices are adjacent to which". A vertex in no face never moves.
Determinism. Same input, same output, every run and at every threads setting. The
volume and centroid reductions are folded in fixed-size chunks and summed in order for
exactly that reason — floating-point addition is not associative, so a reduction tree that
depended on how rayon split the work would change the last bit of the scale factor between
runs.
navis_fastcore.smooth_mesh(faces, vertices, method='taubin', iterations=10, lamb=None, mu=None, alpha=None, beta=None, weights='uniform', preserve_border=False, lock=None, volume_correction=False, threads=None)
¶
Smooth a triangle mesh.
Moves vertices and touches nothing else: the face array, the vertex count and the vertex order all come back unchanged, so anything you have indexed by vertex — synapses, radii, labels — is still attached to the vertex it was attached to.
Three methods, chosen with method:
"taubin" (the default)
Alternating shrink and inflate passes, tuned so the two cancel below a
cut-off frequency. Removes noise without removing the shape, and is the
default for that reason.
"laplacian"
The plain diffusion step: simple, effective, and it shrinks. At
lamb=0.5 and five iterations — what navis.smooth_mesh ships — a
neuron mesh loses 88% of its enclosed volume. Reach for it when the mesh is
a means to an end rather than when its volume means something, or pair it
with volume_correction.
"humphrey"
The HC filter of Vollmer et al., which fights shrinkage by pulling each
vertex back towards where it started rather than towards a lower frequency.
The gentler of the two on fine detail worth keeping.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
vertices
|
TYPE:
|
method
|
TYPE:
|
iterations
|
TYPE:
|
lamb
|
TYPE:
|
mu
|
TYPE:
|
alpha
|
TYPE:
|
beta
|
TYPE:
|
weights
|
TYPE:
|
preserve_border
|
TYPE:
|
lock
|
TYPE:
|
volume_correction
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
vertices
|
New positions, in the same order as the input.
TYPE:
|
Notes
The volume correction scales about the centroid, not the origin. This is the
one place where the result deliberately differs from
trimesh.smoothing.filter_laplacian, which is what navis.smooth_mesh
calls today. Upstream rescales by (vol_before / vol_after) ** (1/3) about the
origin, which is not a shape operation: on the 722817260 test neuron at navis'
own defaults it displaces the mesh by 41 um, and the mesh is 19-26 um across. It
is also not translation invariant — the same mesh smoothed at two different
offsets comes out two different shapes, and far enough from the origin the volume
ratio goes negative and the cube root returns NaN. Scaling about the mesh's own
centroid is the same size change with none of that.
The correction also runs once, at the end, which is not an approximation of running it every iteration but exactly equal to it: every filter here is an affine combination of a vertex and a normalised average of its neighbours, and those commute with a uniform scaling. Upstream pays a full pass over the faces per iteration — 40% of its runtime — for a result it could have had at the end.
When the volume is undefined. On a closed mesh the correction is exactly what it says. A mesh that is not closed still usually gets one, and deliberately: both measurements cone every face back to the same anchor, so their ratio stays a consistent measure of how much the surface shrank even where neither number is an enclosed volume on its own. That matters because meshes worth smoothing are almost never watertight — the 722817260 test neuron is not — and refusing on that basis would refuse on nearly every mesh this exists for.
What is left is the genuinely undecidable case: the ratio of the two signed
volumes is zero, infinite, NaN or negative, so it has no cube root worth taking. A
flat sheet is the clean example, with both volumes exactly zero. There the
vertices come back smoothed but unscaled and a RuntimeWarning says so.
Consistently inverted winding is not in that set — both volumes come out
negative, the ratio is positive, and the correction is as valid as ever.
Non-manifold input is fine, as for :func:simplify_mesh. An edge shared by
three faces, a face naming the same vertex twice, a duplicated face and a vertex
no face mentions are all merely data; nothing here reads more topology than
"which vertices are adjacent to which". A vertex in no face never moves.
Examples:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> # A 5x5 grid with its middle vertex lifted out of the plane.
>>> faces = np.array([[i * 5 + j, (i + 1) * 5 + j, (i + 1) * 5 + j + 1]
... for i in range(4) for j in range(4)]
... + [[i * 5 + j, (i + 1) * 5 + j + 1, i * 5 + j + 1]
... for i in range(4) for j in range(4)], dtype=np.uint32)
>>> vertices = np.array([[i, j, 0.] for i in range(5) for j in range(5)])
>>> vertices[12, 2] = 1.0
>>> v = fastcore.smooth_mesh(faces, vertices, method="laplacian", lamb=1.0,
... iterations=1)
>>> float(v[12, 2]) # back in the plane its six neighbours span
0.0
Pin the rim so an open mesh does not roll inwards:
Capping holes¶
Subsetting a mesh drops every face that loses a corner, which leaves the cut cross-sections
standing open. These four functions find those openings and triangulate them shut. They are
separate rather than one fill_holes because the two ways in enter at different points:
you either have a mesh and want every hole in it closed, or you are about to cut one and
want only the holes the cut itself makes.
Only faces are ever added, never vertices. Every vertex index a caller already holds — in its own face array, in per-vertex data, in a connector table — still points at what it pointed at before, which is what lets the cap be applied after the subset rather than during it.
Every hole in a mesh¶
import navis_fastcore as fastcore
import numpy as np
halfedges = fastcore.boundary_halfedges(faces)
rings, offsets = fastcore.trace_loops(halfedges)
caps = fastcore.triangulate_rings(rings, offsets, vertices)
faces = np.vstack((faces, caps))
Only the holes a cut makes¶
Worked out on the original faces, before the subset, and applied after it — capping only adds faces, so every index handed out by the subset still stands.
exposed = fastcore.exposed_halfedges(faces, dropped)
# ... subset the mesh, then remap `exposed` onto the surviving vertices ...
renumber = np.full(len(vertices), -1, dtype=np.int64)
renumber[kept] = np.arange(len(kept))
exposed = renumber[exposed].astype(np.uint32)
rings, offsets = fastcore.trace_loops(exposed)
caps = fastcore.triangulate_rings(rings, offsets, new_vertices)
exposed_halfedges deliberately leaves out edges that were boundary already: those belong
to openings the mesh came with — a neurite truncated at the edge of the dataset, say — and
sealing those is boundary_halfedges' job. One consequence is worth knowing: where a cut
runs into an opening the mesh came with, what it exposes is an open chain rather than a
ring, and trace_loops abandons it. That hole stays open.
Why these are here¶
Almost all of it is boundary_halfedges. Grouping the 3F edges a face array names is the
whole cost of finding a boundary, and the obvious numpy spelling —
np.unique(keys, return_inverse=True, return_counts=True) — is a stable argsort: 75 ms of
an 84 ms call on a 578k-face mesh. That is not a formulation problem. The bare np.sort of
the same keys is already 51 ms, so no rearrangement in numpy can win. Sorting bare u64
keys in parallel and taking a second pass over the faces to recover each boundary edge's
direction brings the call to 8 ms.
On a 578k-face mesh with ~23k holes punched into it, against the equivalent numpy implementation:
| numpy | fastcore | ||
|---|---|---|---|
boundary_halfedges |
89 ms | 9.1 ms | 10x |
trace_loops |
38 ms | 0.67 ms | 56x |
triangulate_rings |
89 ms | 0.67 ms | 132x |
| end to end | 224 ms | 11 ms | 21x |
And the same mesh on the subset path, 400 twig cuts exposing 4.3k half-edges:
| numpy | fastcore | ||
|---|---|---|---|
exposed_halfedges |
6.4 ms | 0.80 ms | 8x |
trace_loops |
0.97 ms | 0.15 ms | 7x |
triangulate_rings |
2.9 ms | 0.18 ms | 16x |
| end to end | 10.7 ms | 0.87 ms | 12x |
One of those numbers deserves a caveat: triangulate_rings' 132x is not faster ear-clipping —
the C++ it replaces is the same algorithm — it is the disappearance of a 23,000-iteration
Python loop around it, and most of those rings are three vertices and never reach the
ear-clipper at all.
trace_loops is worth a word too, in the other direction. It is a sequential walk and does
not scale with cores at all; what makes it cheap is that it is proportional to the boundary
rather than to the mesh, and that its adjacency is a CSR keyed by vertex id rather than a hash
map of per-vertex lists. That is the difference between 0.67 ms and about 5.
Non-manifold boundaries, and why not a cycle basis¶
At a non-manifold boundary vertex several half-edges leave at once. trace_loops is greedy:
it takes whichever is still free, so every half-edge lands in exactly one ring and the whole
boundary is covered. A cycle basis — networkx.cycle_basis, which is what
trimesh.repair.fill_holes uses — quietly drops the edges that are not part of a simple
cycle, and those holes stay open.
Being greedy means the decomposition depends on the order the half-edges arrive in, which is
why boundary_halfedges and exposed_halfedges both return theirs in 3F edge-list order:
that is the one order that does not depend on how the parallel work happened to be split, so
the same mesh gives the same rings at every threads setting.
It also means a walk can leave a pinch vertex and come back to it, having gone right round
one of the other loops meeting there. What it traced is then a figure of eight rather than a
polygon, and nothing that triangulates a ring is defined on one of those. So a walk is cut
where it crosses itself and the pieces handed on separately: same half-edges, but every ring
trace_loops returns is simple, no vertex twice.
How a ring is closed¶
triangulate_rings ear-clips each ring, trying three things in order:
- The ring flattened through its area-weighted (Newell) normal. Cheaper than a best-fit plane and, on the rings a cut actually produces, it fails slightly less often too.
- The ring flattened through its best-fit plane, from the eigenvectors of its 3x3 scatter matrix.
- The ring as it stands, in three dimensions, clipping whichever ear is cheapest by
area + 0.05 * perimeter².
A ring only gets past step 1 if the flattening self-intersects, which is what makes
ear-clipping run out of ears part way through and yield fewer than the n - 2 triangles a
simple polygon always does. That is not the same as the ring being un-planar: a gently
curved ring can cast a crossed shadow and a folded one need not. Since steps 1 and 2 are
both projections they tend to fail together, which is what step 3 is for — it has no plane
to be defeated by. Over the 933 openings of an invaginated neuron mesh the split was 94.0%,
1.1% and 4.9%.
All three close the hole and wind it correctly, which is what callers depend on. Only step 3 is a heuristic about the shape of the cap: unlike a planar clip it cannot promise the cap does not fold over itself somewhere, which is why it goes last. What it buys is a cap that stays local — on a 1223-vertex opening around an invaginated soma, a median cap edge of 65 nm where a fan out of one vertex gives 2.7 µm.
The cap winds against its ring. The ring runs the way the faces it still has wind it, so a cap that agreed would have the two disagreeing about which side is out.
The ear-clipping is a Rust port of mapbox's earcut rather than a binding to it, so it needs
no extension module of its own. Against mapbox_earcut on the same rings it picks the same
triangles about 93% of the time and an equally valid alternative otherwise — same triangle
count, same total oriented area, same winding — so do not depend on the exact triangles, only
on the hole being closed the right way round.
Rings arriving from trace_loops are simple, and steps 1-3 assume it. A ring built by hand
that names the same vertex twice is a polygon touching itself; it still comes back closed and
correctly wound, but its shape is not specified. Such a ring is also the input that can send
mapbox_earcut — what navis reached for before this module existed — into an infinite
loop on its best-fit-plane retry.
navis_fastcore.boundary_halfedges(faces, threads=None)
¶
Find every edge of a mesh that has only one face on it.
An interior edge has two faces on it and a boundary edge has one, so this is
a grouping of the 3 * F edges the faces name. That grouping is the whole
cost, and it is why this is here: np.unique(keys, return_inverse=True,
return_counts=True) — the obvious way to write it — is a stable argsort,
75 ms of an 84 ms call on a 578k-face mesh, and numpy cannot be talked out of
it (the bare np.sort of the same keys is already 51 ms). Sorting the bare
keys in parallel and taking a second pass to recover the direction brings
that to 8 ms.
Use :func:~navis_fastcore.exposed_halfedges instead where you already know
which vertices are going away — it never looks at the mesh as a whole.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
halfedges
|
Directed half-edges, wound the way their one remaining face
winds them — which is what
:func:
TYPE:
|
Examples:
Two triangles sharing the edge (1, 2): that one is interior, the other
four are boundary.
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.uint32)
>>> fastcore.boundary_halfedges(faces)
array([[0, 1],
[2, 0],
[1, 3],
[3, 2]], dtype=uint32)
A closed mesh has no boundary at all:
navis_fastcore.exposed_halfedges(faces, dropped, threads=None)
¶
Find the edges a subset is about to expose.
Call this with the original faces, before subsetting.
A face survives only if all three of its corners do, so an edge ends up on a
new boundary exactly when it loses a face to the cut but keeps one. Both
halves of that test are local to the cut, so — unlike
:func:~navis_fastcore.boundary_halfedges — this never has to group the
edges of the whole mesh. Only a face losing exactly one corner can leave an
edge behind (lose two and there is no edge left with both ends standing),
which on a real prune is a percent or so of the faces that go.
Edges that were already boundary are left out: they belong to openings the mesh came with — a neurite truncated at the edge of the dataset, say — and sealing those is not this function's business.
| PARAMETER | DESCRIPTION |
|---|---|
faces
|
TYPE:
|
dropped
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
halfedges
|
Directed half-edges, wound the way the one face they have left winds them, with indices into the original vertices. Remap them onto the surviving vertices before capping.
TYPE:
|
Examples:
Two triangles sharing the edge (1, 2). Dropping vertex 0 kills the first
face and leaves (1, 2) newly open, wound the way the surviving face winds
it:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> faces = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.uint32)
>>> fastcore.exposed_halfedges(faces, np.array([True, False, False, False]))
array([[2, 1]], dtype=uint32)
Dropping vertex 3 instead costs the second face, and (1, 2) is exposed
the other way round:
navis_fastcore.trace_loops(halfedges)
¶
Walk directed half-edges into closed rings.
Greedy: at a non-manifold boundary vertex several half-edges leave at once,
so this takes whichever is still free. Every half-edge lands in exactly one
ring, which is what makes this cover the whole boundary — a cycle basis
(networkx.cycle_basis, which is what trimesh.repair.fill_holes uses)
quietly drops the edges that are not part of a simple cycle.
Every ring comes back simple — no vertex twice — which is what makes it a
polygon, and so the thing :func:~navis_fastcore.triangulate_rings is defined
on. A greedy walk does not give that on its own: at a pinch, where several
boundary edges meet at one point, it can leave and re-enter the same vertex,
and what it traced is then a figure of eight. A walk is therefore cut where it
crosses itself and the pieces handed on separately — the same half-edges,
grouped the way the caller can use.
A walk that runs into a dead end abandons what is still in hand, and so is a ring of fewer than three vertices. In both cases the half-edges it consumed stay consumed, so this always terminates — but it does mean the rings need not account for every half-edge handed in. Cycles already split off from an abandoned walk closed on their own account and are kept.
| PARAMETER | DESCRIPTION |
|---|---|
halfedges
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
rings
|
Every ring's vertices, end to end.
TYPE:
|
offsets
|
Where each ring starts and stops: ring
TYPE:
|
Examples:
The four boundary edges of two triangles trace into one ring:
navis_fastcore.triangulate_rings(rings, offsets, vertices, threads=None)
¶
Triangulate boundary rings, wound against the direction they run in.
Three attempts, in order: ear-clipping the ring flattened through its area-weighted (Newell) normal, then flattened through its best-fit plane, and failing both, ear-clipping it in three dimensions without flattening at all. All three are always closed and always correctly wound, which is what everything downstream depends on.
A ring only gets past the first attempt if the flattening self-intersects,
which is what makes ear-clipping run out of ears part way through. That is
not the same as the ring being un-planar — a gently curved ring can cast a
crossed shadow and a folded one need not — and since both of the first two
attempts are projections, they tend to fail together. The third has no plane
to be defeated by; it clips whichever ear is cheapest by
area + 0.05 * perimeter², which is a heuristic about the shape of the
cap rather than a guarantee about it, and is there because the alternative
is a fan reaching from one vertex of the opening to every other.
The cap winds against its ring, because the ring runs the way the faces it still has wind it — a cap that agreed would have the two disagreeing about which side is out.
Rings are independent and run one per worker.
| PARAMETER | DESCRIPTION |
|---|---|
rings
|
TYPE:
|
offsets
|
TYPE:
|
vertices
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
faces
|
New faces, indices into
TYPE:
|
Examples:
A square hole in the z = 0 plane, wound counter-clockwise seen from
+z — so the cap comes back wound the other way:
Projecting for a 2-D renderer¶
Drawing a mesh flat — what navis.plot2d does — takes four steps before a rasteriser
sees anything: project the vertices onto the view plane, drop the faces pointing away
from the viewer, sort what is left along the view axis so that painting it gives correct
occlusion, and lay the survivors out as polygons. project_mesh_2d does all four in one
parallel pass.
rings, bbox, ix, depth, normals = fastcore.project_mesh_2d(
vertices, faces, xy_ix=(0, 1), depth_ix=2, front=1
)
The view is axis-aligned and named by column: xy_ix are the two coordinate columns that
make up the picture, depth_ix is the remaining, into-the-screen one, and front says
which end of that axis the viewer is on — coordinates are never flipped, so a
right-to-left view is the caller's business, not the projection's.
Fusing the steps is the whole point. Each one written the obvious vectorised way in numpy, on an 8.4M-vertex, 16.9M-face neuron:
| step | cost |
|---|---|
project to (V, 2) |
76 ms |
| cull | 226 ms |
| gather the kept faces | 72 ms |
gather the kept corners into (K, 3, 2) |
191 ms |
| close each triangle into a ring | 173 ms |
| bounding box of the result | 534 ms |
| total | 1.27 s |
None of that is arithmetic-bound: they are single-threaded walks over arrays far larger than any cache, and the two gathers plus the ring layout write 900 MB between them to say something the mesh already said. Fused, it is 133 ms, and the bounding box comes out of the same pass that wrote the rings.
The cull is the part worth explaining. Whether a face points at the viewer is the sign of its normal's depth component, and that component is a 2x2 determinant of the two other columns of the edge vectors — which, for an axis-aligned view, are exactly the columns being projected onto. So it never forms the other two components of the cross product and never reads the depth column. It is the same test the full cross product applies, to the bit.
Faces come back as rings: four points each, the first repeated at the end. That is what
a path fill wants — a closed subpath, no separate close-path instruction — and rings[:, :3]
is a view of the plain triangles, not a copy. Emitting triangles and closing them afterwards
is the 173 ms row above plus a second buffer the size of the first.
Two things are optional because most callers do not read them. order=False skips the sort
and the depths: a mesh filled as one path in one colour is drawn under the nonzero winding
rule, which is blind to the order its subpaths arrive in, so the sort cannot change a pixel.
normals=False skips the per-face normals, which only a caller that is shading has any use
for.
Smooth shading is the one thing this does not do. Averaged vertex normals take a contribution from every face, back-facing ones included, so they cannot be computed from the survivors alone; a caller that wants them has to accumulate them itself.
navis_fastcore.project_mesh_2d(vertices, faces, xy_ix=(0, 1), depth_ix=2, front=1, order=True, normals=False, threads=None)
¶
Project a mesh into a view plane: cull, sort and lay out, in one pass.
The view is axis-aligned and named by column: xy_ix are the two coordinate
columns that make up the picture and depth_ix is the remaining,
into-the-screen one. Coordinates are never flipped - a right-to-left view is
the caller's business - which is why front is needed to say which end of
the depth axis the viewer is on.
Doing this in one pass is the point. Each step written the obvious vectorised way in numpy, on an 8.4M-vertex, 16.9M-face neuron:
========================================== ========
step cost
========================================== ========
project to (V, 2) 76 ms
cull 226 ms
gather the kept faces 72 ms
gather the kept corners into (K, 3, 2) 191 ms
close each triangle into a ring 173 ms
bounding box of the result 534 ms
========================================== ========
None of that is arithmetic-bound: they are single-threaded walks over arrays far larger than any cache, and the two gathers plus the ring layout write 900 MB between them to say something the mesh already said.
| PARAMETER | DESCRIPTION |
|---|---|
vertices
|
TYPE:
|
faces
|
TYPE:
|
xy_ix
|
TYPE:
|
depth_ix
|
TYPE:
|
front
|
TYPE:
|
order
|
TYPE:
|
normals
|
TYPE:
|
threads
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
rings
|
Each surviving face as a closed ring of projected corners, its
first repeated at the end - which is what a path fill wants. For
plain triangles take
TYPE:
|
bbox
|
TYPE:
|
ix
|
Index of each surviving face in
TYPE:
|
depth
|
Mean depth of each survivor along
TYPE:
|
normals
|
Unit face normals, zero for a degenerate face.
TYPE:
|
Examples:
Three triangles in the z = 5, z = 0 and z = 0 planes, the last of
them wound the other way round. Looking down +z, that last one is facing
away and goes:
>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> vertices = np.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 0.],
... [0., 0., 5.], [1., 0., 5.], [0., 1., 5.]])
>>> faces = np.array([[3, 4, 5], [0, 1, 2], [0, 2, 1]], dtype=np.uint32)
>>> rings, bbox, ix, depth, _ = fastcore.project_mesh_2d(vertices, faces)
>>> ix
array([1, 0])
The viewer is at +z, so the z = 0 triangle is the far one and comes
back first - painting them in that order is what gets the occlusion right:
Each ring is its triangle's projected corners, closed by a repeat of the first:
>>> rings[0]
array([[0., 0.],
[1., 0.],
[0., 1.],
[0., 0.]])
>>> bool(np.array_equal(rings[:, 3], rings[:, 0]))
True
The bounding box covers those rings, and rings[:, :3] is a view of the
plain triangles, taken without copying:
Flipping front turns the mesh around, and the face that was culled is the
only one left: