Module bestpath 

Best-path phase unwrapping (3D-SRNCP).

Sorting by reliability following a non-continuous path, after H. Abdul-Rahman, M. Gdeisat, D. Burton and M. Lalor, “Fast three-dimensional phase-unwrapping algorithm based on sorting by reliability following a non-continuous path”, Proc. SPIE 5856 (2005) 32–40, and the journal version, H.S. Abdul-Rahman, M.A. Gdeisat, D.R. Burton, M.J. Lalor, F. Lilley, C.J. Moore, “Fast and robust three-dimensional best path phase unwrapping algorithm”, Applied Optics 46(26) (2007) 6623–6635, https://doi.org/10.1364/AO.46.006623.

Reference pipeline

The authors’ C code is distributed by scikit-image (skimage.restoration.unwrap_phase, see below); it takes a wrapped volume and an optional mask and has no parameters besides the periodic-boundary flags, so there are no reference defaults to match beyond those BestPathParams documents. It assumes nothing about voxel size, field strength or orientation (reliabilities are unscaled second differences). It is not selectable in the QSM pipeline’s field_mapping.unwrapping_algorithm (ROMEO, Laplacian or PRELUDE); it is available as UnwrapMethod::BestPath to the slice-wise unwrapper and to MCPC-3D-S’s HIP unwrapping, and its reliabilities as ROMEO’s :bestpath weights.

The method is a path-following unwrapper that never actually follows a path through the volume. Instead:

1. Every voxel gets a reliability — the summed squared second difference of the wrapped phase over the 13 directions of the 26-neighbourhood. A low value means the local phase is smooth and the voxel is trustworthy.
2. Every 6-connected pair of in-mask voxels becomes an edge, whose reliability is the sum of its two endpoints’. Each edge also records the integer wrap count implied by the phase difference across it.
3. The edges are sorted by reliability, most reliable first, and walked in that order. Each edge joins the two voxel groups at its ends, offsetting one group by the wraps the edge implies. Edges whose endpoints already share a group are skipped.

So the unwrapping order is global rather than local: the most trustworthy phase differences in the whole volume are committed first, and noisy regions are reached last, through whatever reliable route survives. That is what makes it robust to residues compared with a seeded flood fill.

Relationship to ROMEO

calculate_weights_bestpath also implements Abdul-Rahman reliabilities, but as edge weights fed into ROMEO’s minimum-spanning-tree growing (matching ROMEO.jl’s calculateweights_bestpath, including its inversion to 1/sqrt(D²) and quantisation to u8). That is ROMEO with different weights; this module is the original algorithm, which merges groups over a globally sorted edge list rather than growing a region from a seed.

Deviations from the reference implementation

The algorithm here is a reimplementation, informed by the C reference that scikit-image distributes as skimage/restoration/unwrap_3d_ljmu.c under the BSD-3-Clause licence (see THIRD_PARTY_NOTICES.md). Three deliberate departures:

• Group merging uses a union-find with integer potentials (UnionFind) rather than the reference’s explicit per-group linked lists. The merge decisions and the resulting wrap counts are the same; only the bookkeeping cost differs (O(α) per edge instead of O(size of the smaller group)).
• Masked voxels keep their original wrapped phase on output. The reference overwrites them with the minimum of the unwrapped phase, which is a display convenience; leaving them alone matches unwrap_romeo and the rest of this crate.
• Wrap-around (periodic) boundary options are not implemented. They exist in the reference for interferometry; MRI volumes are not periodic.

An upstream quirk, deliberately preserved

The reference has an extend_mask step that reads as though it should exclude mask-adjacent voxels from the reliability calculation — it tests whether all 26 neighbours are unmasked. It does not: extended_mask is calloc’d to NOMASK, and the function only ever assigns NOMASK, never MASK, so the buffer is uniformly NOMASK and the guard that reads it is always true. The mask therefore plays no part in reliability at all; it only decides which edges exist.

This module reproduces the effective behaviour, not the apparent intent, and that is a deliberate choice rather than bug-compatibility for its own sake. Implementing the intent — stand-in reliabilities for mask-adjacent voxels — measurably unwraps worse near the ROI boundary on noisy data, because it throws away a graded quality signal and replaces it with a shuffle. A boundary voxel whose neighbours lie in noise already earns a poor reliability on its own and sorts late for the right reason.

The absolute 2π offset of each connected component is arbitrary — it is fixed only up to the choice of which voxel ends up as the component’s representative. Set BestPathParams::correct_global for a defined choice.

Verification

Checked against skimage.restoration.unwrap_phase (scikit-image 0.26.0) on 48×48×32 synthetic fields — smooth, noisy, and masked to an irregular ROI with holes. Every voxel matched bit-for-bit, across several of scikit-image’s own rng seeds. The two diverge only once the added noise approaches ±π, where the wrap counts are genuinely ambiguous: there they disagree on 0.07% of voxels, which is the same margin by which scikit-image disagrees with itself between rng seeds (0.06%).

tests/bestpath_skimage_parity.rs keeps this checked, against scikit-image 0.26.0 in the verified QSM-CI py-originals image (QSMxT/QSM-CI#302; tests/data/bestpath_skimage_small.py regenerates it in Actions): on two small phantoms, one noisy, every voxel has scikit-image’s wrap count, and the values agree to 1.8e-15 rad (measured 2026-10-30). In the noisy one the whole volume sits one 2π above scikit-image’s: the component’s absolute offset (above). On a travelling head, every voxel took scikit-image’s wrap count too.

It was this comparison that turned up the extend_mask quirk above. Following the apparent intent cost an order of magnitude in accuracy near the ROI boundary at high noise — 0.51% of wrap counts wrong against ground truth, versus 0.06% for the reference’s actual behaviour, which this module now reproduces exactly.

Structs

BestPathParams

Parameters for unwrap_bestpath.

Edge 🔒

One edge of the voxel graph, awaiting its turn in the sorted merge.

Constants

DIRECTIONS 🔒

The 13 direction vectors of the 26-neighbourhood, one per antipodal pair.

PI 🔒

TWO_PI 🔒

Functions

build_edges 🔒

Build the 6-connected edge list, skipping any pair with a masked endpoint.

compute_reliability 🔒

Per-voxel reliability: the summed squared second difference of the wrapped phase over the 13 directions of the 26-neighbourhood. Lower is better.

correct_global_offset 🔒

Subtract the median wrap count over the mask, so the result stays centred near the input range. Mirrors ROMEO’s correct_global_offset.

find_wrap 🔒

Number of 2π that separates two wrapped phases, as the reference’s find_wrap.

sort_edges 🔒

Sort edges most-reliable (lowest value) first.

splitmix64 🔒

SplitMix64, used to reproduce the reference’s rand() initialisation deterministically and independently of iteration order.

stand_in_reliability 🔒

A pseudo-random stand-in reliability, uniform over [0, 2^31).

unwrap_bestpath

Unwrap a 3D phase volume with the best-path (3D-SRNCP) algorithm.

wrap_pi 🔒

Wrap a value already known to lie within ±3π into [-π, π].