dijkstra3d-sparse

Dijkstra shortest paths, distance fields and connected components over sparse 3D voxel sets given as an (N, 3) integer coordinate array — a sparse analogue of seung-lab/dijkstra3d, which operates on dense 3D arrays. Rust core, Python/NumPy frontend.

Why

dijkstra3d is fast because it never builds an explicit graph: it walks an implicit rectangular grid where a voxel's neighbours are generated by coordinate offset. The only thing making it "dense" — and the reason it needs memory proportional to the bounding-box volume W·H·D — is that its coordinate → payload lookup is a dense array sized to the full box.

For sparse objects (a thin structure inside a large box, N ≪ W·H·D) that is wasteful. This library keeps the implicit-grid walk and swaps that one dense component for a sparse hash coordinate → compact index [0, N). Everything else — binary heap, edge relaxation, parent tracking, path reconstruction — is unchanged. No adjacency list is ever materialized (this is not a CSR/explicit-graph Dijkstra), and all working memory is O(N), independent of the bounding box.

From benchmarks/RESULTS.md: a 1.5M-voxel helical tube in a 16,267 × 4,005 × 4,006 bounding box solves in ~0.4 s within ~470 MiB peak RSS — where a dense field over the same box would need 4 TiB. Going coordinates → distance field through scipy.sparse.csgraph.dijkstra instead takes ~3.5 s and 1.6 GiB peak on the same workload: SciPy's solver itself is fast, but it first needs the ~30M-edge CSR graph materialized — exactly the step the implicit-grid walk skips.

One caveat: if you already hold a CSR graph and only solve on it repeatedly, SciPy's solver alone is competitive (~0.14 s on this workload once the graph exists). The advantage here is going from raw coordinates to a field — the typical starting point for voxel data — without ever paying the time and memory to build an edge list.

Install

pip install dijkstra3d-sparse

Pre-built wheels cover Linux / macOS / Windows, Python 3.9+. Building from source needs a Rust toolchain (pip invokes it automatically via maturin).

Quickstart

import numpy as np
import dijkstra3d_sparse as ds

# a sparse voxel set: (N, 3) integer coordinates, any origin, unsorted OK
voxels = np.argwhere(volume > 0).astype(np.int32)   # e.g. from a dense mask
# ... or coordinates that never lived in a dense array at all

# distance + predecessor field from voxel row 0
dist, pred = ds.dijkstra_field(voxels, sources=0, connectivity=26,
                               anisotropy=(16.0, 16.0, 40.0))

# shortest path to the voxel farthest from the source
target = int(np.argmax(np.where(np.isfinite(dist), dist, -1)))
coords = ds.path(voxels, pred, target, dist=dist)    # (M, 3), source → target

# connected components over the same implicit grid
n_components, labels = ds.connected_components(voxels, connectivity=26)

# hold coordinates instead of row indices? map them first
src = ds.index_of(voxels, [[10, 4, 2], [0, 0, 0]])
dist, pred = ds.dijkstra_field(voxels, src)          # multi-source: dist to nearest

dist/pred are 1-D arrays aligned 1:1 with the rows of voxels (the key difference from dijkstra3d, whose field is a dense 3D array). Unreached voxels get dist = +inf, pred = -1; -1 matches SciPy's "no predecessor" sentinel, so (dist, pred) is a drop-in for scipy.sparse.csgraph.dijkstra(..., return_predecessors=True) on the equivalent explicit graph.

API

dijkstra_field(voxels, sources, *, node_cost=None, connectivity=26,
               anisotropy=(1.0, 1.0, 1.0), cost_mode="vertex",
               free_mask=None, free_eps=1e-6, min_only=True,
               stop_mask=None, stop_count=1,
               index_kind="hash") -> (dist, pred)

shortest_path(voxels, source, target, **kw) -> (path, cost)  # early exit

shortest_path_to_set(voxels, source, stop_mask, **kw) -> (path, hit, cost)

path(voxels, pred, target, *, dist=None) -> (M, 3) int32   # source → target

connected_components(voxels, *, group=None,
                     connectivity=26) -> (n_components, labels)

label_adjacency(voxels, labels, *, connectivity=26) -> (K, 2) int64

exposed_faces(voxels, *, index_kind="hash") -> (N,) uint8   # surface-face mask

factorize(voxels, *, return_index=False,                    # dedup coords -> labels
          index_kind="hash") -> (n_labels, labels[, reps])

ray_exits(voxels, origins, directions, *, max_dist,         # boundary crossings
          max_crossings=1, index_kind="hash") -> (t, n_hits)

index_of(voxels, coords, *, strict=True) -> int | (M,) int64

Graph(voxels, *, index_kind="hash")   # reusable handle, methods below

Reusable Graph handle

Every free function above rebuilds the coordinate → row spatial index — the one O(N) setup cost — on each call. For repeated queries over the same voxel set, build a Graph once; it holds the index and exposes the same operations as methods, minus the voxels/index_kind arguments:

g = ds.Graph(voxels, index_kind="hash")   # O(N) index build happens here, once

dist, pred = g.dijkstra_field(0, cost_mode="geometric")     # reuses the index
dist2, _   = g.dijkstra_field([3, 7], node_cost=penalty,    # different cost model,
                              cost_mode="additive")         # same handle
coords, hit, cost = g.shortest_path_to_set(q, anchors)      # grafting primitive
n_comp, labels = g.connected_components()
n_rings, rings = g.connected_components(group=level)        # components per group
ring_edges = g.label_adjacency(rings)                       # which rings touch
face_mask = g.exposed_faces()                               # surface-face mask
t, n_hits = g.ray_exits(origins, dirs, max_dist=caps)       # ray boundary crossings
rows = g.index_of(coords)
g.n, g.voxels, g.index_kind                                 # introspection

Only voxels and index_kind are fixed at construction — connectivity, anisotropy, cost_mode, node_cost and the masks stay per-call, so one handle serves queries with different cost models. Results are identical to the free functions (same code runs; only where the index is built moves), duplicate coordinates are rejected at construction, and the handle keeps its own copy of the coordinates, so it is unaffected by later mutation of the input array. The payoff scales with call count — grafting loops that issue one shortest_path_to_set per path are the motivating case (see benchmarks/RESULTS.md).

Edge-cost model

Step lengths are precomputed per offset from anisotropy = (wx, wy, wz), matching dijkstra3d exactly: axis moves cost wx/wy/wz, face diagonals sqrt(wa² + wb²), corner diagonals sqrt(wa² + wb² + wc²). The cost of the directed edge cur → nbr is then:

With node_cost=None every mode reduces to the pure geometric step length. Costs must be finite and non-negative (Dijkstra invariant; validated at the boundary).

free_mask: edges into masked voxels cost free_eps (small, strictly positive) in total. This supports incremental path extraction where later paths should ride an already-selected node set for ~free before diverging.

min_only=False runs one Dijkstra per source and returns (S, N) arrays, mirroring SciPy; the default True returns a single (N,) field of distances to the nearest source.

Early termination & search-to-a-set

Dijkstra settles nodes in non-decreasing distance order, so the moment a node is popped its distance and path are final. stop_mask exploits this: the search stops as soon as stop_count masked voxels have been settled (default 1 — i.e. at the nearest member of the set), returning a partial field that is exact on everything it touched and +inf/-1 beyond. SciPy's limit distance cutoff cannot express "stop when you reach node X / this set". Two wrappers make this ergonomic:

# point → point, terminating the instant the target settles
coords, cost = ds.shortest_path(voxels, source, target)

# point → nearest member of an anchor set
coords, hit, cost = ds.shortest_path_to_set(voxels, source, anchor_mask)
# hit = row index of the anchor reached (-1 + empty path if unreachable)

This is the primitive for incremental tree construction (grafting — e.g. centerline/skeleton extraction): repeatedly connect a query voxel to a growing anchor set, where each query only explores the local catchment between the query and the nearest anchor instead of the full voxel set:

anchors = np.zeros(len(voxels), dtype=bool)
anchors[seed] = True
for query in queries:
    coords, hit, cost = ds.shortest_path_to_set(voxels, query, anchors)
    anchors[ds.index_of(voxels, coords)] = True   # graft the spur

On the benchmark tube (1.5M voxels), 60 such grafts run in ~1.4 s total, with per-query touched voxels falling from ~10% of N (sparse anchors) to ~0.3% (dense anchors) — versus 100% of N per query for repeated full fields. stop_mask composes with everything else: with multiple sources and min_only=True it means "grow a field from all sources until it first touches the anchor set". It is also the recommended replacement for free_mask-based grafting tricks — cleaner (no cost distortion) and cheaper (early exit); if both are given they stay independent (free_mask changes edge costs, stop_mask only changes termination).

Graph contraction without an edge list

Two primitives run over the same implicit-grid probe as everything else, so a caller can contract the voxel graph — collapse voxels into groups, then ask which groups touch — without ever materializing adjacencies:

• connected_components(voxels, group=values) connects two voxels only when they are connectivity-adjacent and group[u] == group[v], i.e. the components of the sub-graph induced by each group value. Grouping only constrains unions: the same value in two spatially separate places stays two components, and a voxel whose group differs from all its neighbours' becomes a singleton. group=None is the plain component labelling.
• label_adjacency(voxels, labels) returns the distinct pairs of different labels that touch, as a sorted (K, 2) array — the edges of the quotient graph. It deduplicates during the probe, so the (typically enormous) intermediate adjacency count never exists. labels need not be dense or non-negative.

Together they express level-set / Reeb-graph constructions such as wavefront skeletonization in three passes over one Graph:

g = ds.Graph(voxels)
n_comp, comp = g.connected_components()                    # one wave per component
seeds = [int(np.flatnonzero(comp == c)[0]) for c in range(n_comp)]
dist, _ = g.dijkstra_field(seeds, cost_mode="geometric")

level = np.floor(dist / step_size).astype(np.int64)        # geodesic level sets
n_rings, rings = g.connected_components(group=level)       # rings = level components
skeleton_edges = g.label_adjacency(rings)                  # contract onto rings

On the benchmark tube (1.5M voxels) that pipeline runs in 1.1 s at 514 MiB peak RSS, versus 5.7 s at 3.6 GiB for the same result via an explicit edge list plus SciPy — 30.5M adjacencies materialized to yield 2,429 distinct ring pairs (see benchmarks/RESULTS.md).

Surface faces

exposed_faces(voxels) answers, for every voxel in one pass, which of its six face-neighbours are absent from the set — the first stage of any voxel mesher / surface extraction. It returns an (N,) uint8 mask, bit k set iff the neighbour across face k is missing, in the order +x, -x, +y, -y, +z, -z:

mask = ds.exposed_faces(voxels)             # (N,) uint8; 0 = interior, 63 = isolated
right = voxels[(mask & (1 << 0)) != 0]      # voxels whose +x face is exposed

Unlike the other free functions it does not route through a Graph — and builds no spatial index at all. It sorts one packed 16-byte key per voxel and sweeps the three positive face offsets across it as three linear merges (a hit clears the far voxel's opposite bit, so the other three offsets are free), then drops that array inside the one call: 16 B/voxel of working set, and the surface pass leaves nothing behind for a later stage's peak to stack on. Already-sorted input — what np.argwhere and np.unique(..., axis=0) hand over — makes the sort a single scan. index_kind is accepted for signature parity only, as with factorize.

(Graph.exposed_faces() exists too, but reuses — and keeps — the handle's index, so it gives up that transient-memory win. It is also the one query where the backends diverge: a "sorted" handle sweeps its key array with no lookups; a "hash" one must probe.)

Deduplicating coordinates

factorize(voxels) is the sparse np.unique(coords, axis=0, return_inverse=True): it assigns every row a dense label, equal exactly when the coordinates are equal, in one pass over the rows instead of a sort. It is the only primitive here that accepts duplicate coordinates — collapsing them is the whole point (the others reject repeats).

n, labels = ds.factorize(cells)                      # (E, 3) with repeats -> labels
n, labels, reps = ds.factorize(cells, return_index=True)
unique = cells[reps]                                 # one row per label...
assert np.array_equal(unique[labels], cells)         # ...and labels index back

Labels are 0 .. n-1 in order of first appearance by row; reps[k] is the first row carrying label k. This is the dedup a voxel mesher runs on its per-quad corners, and the coarse-cell assignment (fine // scale) a downsampler runs on its nodes — the same operation, so it lives here rather than as a hand-rolled argsort in each caller. It groups by exact coordinate equality, which is a different question from connected_components (spatial adjacency, unique input) despite the shared (n, labels) return.

Both of those inputs are derived from a voxel grid, so their bounding box is compact — and a compact box is resolved through a direct-address table instead of a hash map: no hashing, no collisions, and less memory than the map would have reserved. Coordinates scattered across a wide box fall back to hashing. A cost decision measured from the input's bounding box; the labels are identical either way, and index_kind is signature parity only.

Ray exits

ray_exits(voxels, origins, directions, max_dist=…) walks each ray through the voxel set and reports where it crosses the object's boundary. Voxel centres sit on integer coordinates, so cell c occupies [c - 0.5, c + 0.5); each ray p(t) = origin + t · direction (t ≥ 0) is stepped with a 3-D DDA (Amanatides & Woo) that visits exactly the cells it passes through — no sampling, no interpolation — and every t at which occupancy flips is reported:

g = ds.Graph(voxels)                        # build the index once...
t, n_hits = g.ray_exits(origins, dirs,      # ...then cast in chunks
                        max_dist=caps, max_crossings=2)

radius    = np.where(n_hits > 0, t[:, 0], caps)   # first exit = the radius
escaped   = n_hits == 0                           # never left within max_dist
reentered = n_hits > 1                            # not star-shaped this way

t[r, 0] is the first exit — the cross-section radius along that direction — and later entries strictly alternate re-entry / exit, so n_hits > 1 says the object is not star-shaped about the origin along that ray. Padding beyond n_hits[r] is +inf.

• directions are index-space and need not be unit. Pass a physically-unit direction divided by the voxel spacing and t comes back as a physical distance — which is why there is no anisotropy parameter here: spacing is the caller's metric, not the library's.
• max_dist bounds t, not the cell count, and is per-ray (scalar or (R,)). A ray reaching it without crossing gets n_hits = 0 — it escaped.
• The origin cell is assumed occupied and never itself reported; a ray starting in an empty cell has nothing to exit and returns n_hits = 0.
• Unlike exposed_faces/factorize, index_kind is not inert here: a ray walk is a stream of unpredictable point probes with no exploitable ordering, which is what "hash" (the default) is for. "sorted" returns identical results, more slowly.

Prefer Graph.ray_exits over the free function: rays are normally fired in chunks against one voxel set, and rebuilding the index per chunk would dominate the walk. The per-ray state is 13 scalars in registers and nothing is allocated in the loop — the point being that a vectorized DDA must instead keep the live ray set in (R, 3) arrays and touch all of it to advance any single ray by one cell, which costs an order of magnitude more than the index probes it wraps.

Notes

• Multiple sources: seed them all — one pass computes distance-to-nearest-source and predecessors pointing back to each voxel's nearest source.
• Output is deterministic: heap ties break on row index, so identical inputs give identical fields across runs and platforms.
• Duplicate coordinates in voxels raise ValueError.
• Coordinates may be negative and use the full int32 range; there is no bounding-box extent limit.
• index_kind selects the spatial-index backend ("hash" FxHashMap probes, default; "sorted" binary search over sorted keys, slightly lower memory). Results are identical.

Development

uv venv && source .venv/bin/activate
uv pip install numpy scipy pytest maturin
maturin develop --release --uv   # build the Rust extension into the venv
pytest                           # Python test suite (SciPy parity + properties)
cargo test                       # Rust unit tests
python benchmarks/bench.py      # benchmark + O(N) memory gate

The test suite asserts parity with scipy.sparse.csgraph on the equivalent explicit CSR graph for all cost modes, connectivities and anisotropies, plus structural invariants (source distance 0, triangle inequality along edges, path adjacency/cost).

License

GPL-3.0-or-later, like dijkstra3d.