gridflood

Multi-source flood fill on regular grids, plus the algorithms that are the same traversal with a different queue: connected components, nearest-seed allocation (raster Voronoi with obstacles), cost distance and cost allocation. Dependency-free, with a tiled driver that reproduces the global result exactly, so it runs out of core and in parallel.

What it does

One kernel, three entry points:

function	what	queue
fill(.., Cost::Unit)	multi-source flood fill; nearest seed by hop count	FIFO (BFS)
fill(.., Cost::Euclidean) / Cost::PerCell(&cost)	cost distance + cost allocation	heap (Dijkstra)
components(..)	connected-component labelling	FIFO, lazily seeded
relax(..)	continue any fill from new frontier cells	heap

Every result is a label and a dist per cell. The defining rule is that each reachable cell gets the lexicographically smallest (dist, label) over all seeds. Ties in distance go to the lower label. That single rule makes the answer independent of seed order, traversal order and tiling, which is what lets tile::fill_tiled and tile::components_tiled give bit-identical output to the global kernels.

Layout

Row-major, cell = row * ncol + col, 0-based. Inputs are slices: a &[bool] fillable mask and, for Cost::PerCell, a &[f64] cost. Labels are u32 with NO_LABEL (u32::MAX) for unreached cells; distances are f64 with INFINITY for unreached.

use gridflood::{fill, Connectivity, Cost, Dims, Seed};

let dims = Dims::new(20, 30);
let fillable = vec![true; dims.len()];
let seeds = [Seed::new(dims.cell(0, 0), 1), Seed::new(dims.cell(19, 29), 2)];
let f = fill(dims, &fillable, &seeds, Connectivity::Four, Cost::Unit);
// f.label: which seed won each cell; f.dist: Manhattan hops to it

Tiling

use gridflood::tile::{fill_tiled, CostModel, SliceSource, TileScheme};

let src = SliceSource::new(dims, &fillable);           // or your own TileSource
let out = fill_tiled(dims, TileScheme::new(512, 512), &src, &seeds,
                     Connectivity::Eight, CostModel::Unit);
let f = out.assemble();                                 // full-grid Fill

Each tile is filled from its own seeds, then the driver walks every seam: where a cell could be improved by its neighbour across the seam it is stamped and the tile is re-relaxed from those cells. Distances only decrease, so this converges, typically in as many rounds as the longest wavefront crosses tiles (5-6 on a 4000 x 4000 grid with 200 seeds). No halo is needed. Components are merged with union-find over seam adjacencies and renumbered in scan order.

TileSource is the out-of-core hook: implement fillable(&TileRect) (and cost) over a GDAL dataset, a COG, a Zarr store, and the driver never holds more than the tiles. SliceSource is the in-memory version.

Enable features = ["parallel"] for rayon across tiles.

Semantics

• 4-connectivity with Cost::Unit gives Manhattan distance; 8-connectivity gives Chebyshev. Cost::Euclidean charges sqrt(2) for diagonal steps.
• Cost::PerCell charges step_length * (cost_from + cost_to) / 2, the usual raster cost-distance convention. Non-finite or negative cost cells are treated as not fillable.
• Seeds on non-fillable cells are ignored. Two seeds on one cell: lower label wins.
• relax lets you add seeds to a finished fill, or push a fill into a region you have just unblocked, without starting over.

Binding from R

The kernels take plain slices, so an extendr wrapper is a few lines.

gridflood's row/col meaning is the same as ximage's: an R matrix m with m[1, 1] at the top-left, rows running down the image and columns across, exactly as the matrix prints. The only difference is storage order. R holds that matrix column-major (cell = col * nrow + row, 0-based); gridflood reads row-major (cell = row * ncol + col). Because every algorithm here is symmetric in the two axes, the cheapest bridge is to swap the dimension names rather than the data:

# m is an nrow x ncol matrix in ximage orientation
out <- gridflood_fill(as.vector(m), dims = c(ncol(m), nrow(m)), ...)
matrix(out, nrow(m), ncol(m))   # back in ximage orientation, no transpose

On the Rust side that is Dims::new(ncol, nrow): gridflood's "row" is R's column, each of length nrow. Seed cells given as R matrix cell numbers ((col - 1) * nrow + (row - 1), or which(m == x) - 1) pass straight through unchanged, and 4/8 connectivity, distances and labels are unaffected by the swap.

Two cautions. Tile sizes swap too: TileScheme::new(tile_ncol, tile_nrow) in that call. And cell numbers from terra (cellFromRowCol) are row-major, so they do not match the swapped layout; convert them with rowColFromCell and recompute, or transpose the data instead of swapping the dims.

NO_LABEL maps to NA_integer_ and labels are 0-based; add 1 at the boundary.

Performance

4000 x 4000 grid, 200 seeds, 30 percent blocked, 8-connectivity, single thread on a modest container:

	ms
global fill, unit cost (BFS)	1100
global fill, per-cell cost (Dijkstra)	4900
tiled fill 1024 x 1024, unit	1400 (880 with parallel, 2 cores)
global components	800
tiled components 512 x 512	880 (590 with parallel)

cargo run --release --example bench reproduces this.

Not here (yet)

• Priority-flood depression filling. Same heap, different update rule; a natural addition.
• Jump flooding for GPU/wasm targets.
• A tolerance-based fillable predicate (magic wand). Compute the mask yourself for now; a closure-based fillable is the obvious extension.

License

MIT OR Apache-2.0