Module romeo 

ROMEO weight calculation for phase unwrapping

Calculates edge weights for region-growing phase unwrapping based on:

• Phase coherence
• Phase gradient coherence (multi-echo)
• Magnitude coherence
• Magnitude weights

Reference: Dymerska, B., Eckstein, K., Bachrata, B., Siow, B., Trattnig, S., Shmueli, K., Robinson, S.D. (2021). “Phase unwrapping with a rapid opensource minimum spanning tree algorithm (ROMEO).” Magnetic Resonance in Medicine, 85(4):2294-2308. https://doi.org/10.1002/mrm.28563

Reference implementation: https://github.com/korbinian90/MriResearchTools.jl (romeo/unwrap re-exported from https://github.com/korbinian90/ROMEO.jl).

Reference pipeline

Facts from ROMEO.jl (commit 08776d4, 2026-09-02) and MriResearchTools.jl (771eadf), the authors’ code; file:line refer to those trees.

• Weights. :romeo (phase coherence, phase gradient coherence, phase linearity, magnitude coherence) is the API’s default (src/weights.jl:56,75-76), which RomeoParams::default matches; :romeo4 sets the same four flags (src/weights.jl:81-82). The CLI’s -w romeo default becomes romeo3 (no phase linearity) whenever a magnitude is given and romeo4 otherwise (ext/RomeoApp/caller.jl:49-55): RomeoParams::romeo_cli.
• Weight bins. Each weight in [0, 1] goes to the bin max(round((1 − w)·255), 1) (rescale, src/weights.jl:139-146); the bins order the growth, so they are ROMEO.jl’s here (stored as 256 − bin). An edge is weighted when its first voxel is in the mask (mask[I] && checkbounds(J), src/weights.jl:99-100), with the magnitude masked, so the growth can step one voxel out of the mask along +x, +y or +z (and through such voxels).
• Multi-echo mode. Temporal (template) unwrapping is the default: the template echo (template=1, i.e. index 0 here) is unwrapped spatially, with the echo after it (or before it, for a later template) as phase2 for the gradient-coherence weight, and the others temporally outward from it (src/unwrapping.jl:68-93). Individual unwrapping is opt-in, and the CLI says it “might be necessary if there is large movement (timeseries) or phase-offset-correction is not applicable” (argparse.jl:95-100). Both match RomeoParams::default (template: 0, individual: false).
• temporal_uncertain_unwrapping. 0.5 in the API (src/unwrapping.jl:52,70), which RomeoParams::default and CLEAR-SWI’s romeo_combine use; the CLI’s default is 0 and it marks the option EXPERIMENTAL (argparse.jl:165-172); RomeoParams::romeo_cli has 0.
• correct_global is off by default in both the API and the CLI (src/unwrapping.jl:43, argparse.jl:129-135). (QSMxT’s Python releases ran the CLI with --individual-unwrapping --correct-global, qsmxt/interfaces/nipype_interface_romeo.py:66; that is QSMxT’s choice, not ROMEO’s default.) Where it applies, as in ROMEO.jl’s 3D unwrap!: to a single volume, to the template echo before the others are unwrapped temporally from it, and in individual mode to each echo and then across echoes (correct_multi_echo_wraps); see unwrap_individual for the two places (threads and even medians) where QSM.rs is deterministic and ROMEO.jl is not.
• Undefined phase (NaN). ROMEO.jl takes NaN phase (e.g. MCPC-3D-S’s output where its offset is undefined, crate::utils::phase_offset_removal) as “no data”: any edge touching it gets a NaN weight, which rescale makes invalid, so the growth never reaches it and it stays NaN; the global corrections skip non-finite voxels. QSM.rs does the same, including Julia’s NaN-propagating max(0, x) in the gradient-coherence and phase-linearity weights (phase linearity is 0.5 next to a NaN, as in ROMEO.jl; before 2026-11 it was 0 here).
• Mask. The CLI unwraps inside robustmask of the template echo’s magnitude by default (argparse.jl:38-42, caller.jl:201-204), a magnitude threshold with hole filling (MriResearchTools src/masking.jl:7-29), not a brain mask; with no magnitude it unwraps unmasked. CLEAR-SWI calls ROMEO with no mask at all (CLEARSWI.jl src/phase_processing.jl:122). The QSM.rs pipeline passes its brain mask.
• Phase offsets. When the CLI computes a B0 map from multi-echo data it turns MCPC-3D-S phase-offset correction on (monopolar) unless told otherwise (caller.jl:80-91), with smoothing sigma 7 mm converted to voxels (argparse.jl:79-91, caller.jl:137-141); bipolar correction is opt-in and needs ≥ 3 echoes (argparse.jl:71-76). MCPC-3D-S itself is in crate::utils::phase_offset_removal; its defaults are in the module docs of utils::multi_echo. calculateB0_unwrapped states that “the phase offsets have to be removed prior” (MriResearchTools src/romeofunctions.jl:3-7).
• Seed. ROMEO.jl seeds at the voxel whose three forward edges have the best summed bin (ties to the highest linear index) and corrects the seed’s 2π count from phase2 with multi-echo input, or wraps it into (−π, π] otherwise (src/seed.jl). So does every path here ([crate::region_grow::grow_region_unwrap_full]). Before, the default single-seed path seeded at the mask’s centre of mass (find_seed_point, no longer used here), which could change the global 2π offset and the growth order.
• Resolution, field strength, orientation. Nothing in ROMEO depends on them: the weights use phase differences between 6-connected neighbours and echo-time ratios only.

Agreement with the reference

In double precision QSM.rs reproduces ROMEO.jl’s unwrap exactly (the same number of turns at every voxel) in template mode with the API’s and the CLI’s settings and in individual mode with global correction, and MriResearchTools’ mcpc3ds to 9e-16 rad (tests/romeo_romeojl_parity.rs, a noisy three-echo phantom; ROMEO.jl 08776d4 and MriResearchTools 771eadf in the QSM-CI julia-qsm image sha256:b3c048ac…). Given the CLI’s own MCPC-3D-S output on the QSM-CI in-vivo travelling-head scan (Prisma 2.89 T, 1 mm, five echoes), QSM.rs’s template unwrapping gives the CLI’s unwrapped echoes exactly.

tests/romeo_validity_romeojl_parity.rs (ROMEO.jl 1.6.0, MriResearchTools 3.8.0, same image, one Julia thread) adds a phantom whose brain mask reaches past MCPC-3D-S’s mask: the NaN validity, individual mode with and without global correction, and template mode with and without it, on the corrected and on the raw phase, all give ROMEO.jl’s turns at every voxel and its B0 to 1e-14 (measured 2026-11-04). Known gap: template mode’s re-unwrapping of uncertain voxels (temporal_uncertain_unwrapping > 0, the API default) is not yet ROMEO.jl’s. ROMEO.jl rates the voxels with voxelquality (Float32 weights before binning, no mask, the sum over six edges divided by six, NaN for an edge touching a NaN) and grows from every voxel it marks visited, those outside the mask included (temporal_uncertain_unwrapping!); temporal_uncertain_rewrap uses voxel_quality_romeo (binned weights, in-mask neighbours only) and grows from the certain voxels inside the mask. On that phantom 14 of 3,000 voxels of echoes 2 and 3 differ after MCPC-3D-S and 850 of 3,392 on the raw phase; the fixture’s template cases therefore set it to 0 (the CLI’s default). QSMxT’s default (individual mode) does not use it.

The romeo CLI itself reads the NIfTI phase and magnitude as Float32 and computes in single precision; QSM.rs computes in f64 (kept: double precision throughout, Ashley 2026-10-28). That alone leaves a residual on the in-vivo scan. The CLI’s defaults (romeo -B with a brain mask) against the pipeline with RomeoParams::romeo_cli: B0 NRMSE 2.9 % (corr 0.9995), from 2π differences in ~1,250 voxels of the MCPC-3D-S HIP unwrap that propagate through the offset smoothing (fed the CLI’s own single-precision HIP, QSM.rs differs in 36 voxels; fed its corrected phase, in none). Before the seed, weight bins, edge masking, mcpc3ds’s own mask and smoothing were ROMEO.jl’s, the same comparison gave NRMSE 25 %; with individual unwrapping, global correction, romeo4 and a 4 mm offset smoothing on both sides it went from 6.4 % to 1.1 %.

Structs

RomeoParams

Parameters for ROMEO phase unwrapping.

Enums

RomeoWeightType

Weight calculation scheme for ROMEO unwrapping.

Constants

TWO_PI 🔒

Functions

calculate_weights_bestpath

Calculate BestPath edge weights.

calculate_weights_romeo

Calculate ROMEO edge weights for phase unwrapping.

calculate_weights_romeo_configurable

Calculate ROMEO edge weights with configurable weight components.

calculate_weights_romeo_with_flags

Calculate ROMEO edge weights with full 6-component flag control.

calculate_weights_single_echo

Simplified weight calculation for single-echo data (no phase2)

compute_weights_from_params 🔒

Compute weights from params (dispatches BestPath vs ROMEO).

correct_global_offset 🔒

Correct global 2π offset for a single 3D volume.

correct_multi_echo_wraps

Correct global 2π offsets between echoes using median wrap counting.

correct_regions 🔒

Correct each region’s median to nearest 0 by subtracting n·2π.

find_seed_point

Find a seed point for region growing: the centre of mass of the mask, or - when that falls outside the mask, as it does for a disconnected or concave one - the nearest voxel that is inside it.

finish_3d_unwrap 🔒

What ROMEO.jl’s 3D unwrap! does after the growth (src/unwrapping.jl, src/algorithm.jl): region merging and correction, then, with correct_global, wrapped .-= 2π·median(round(φ[mask] / 2π)). Its multi-echo modes call that same 3D unwrap! for the template echo and for every echo in individual mode, so all three paths finish here.

gamma 🔒

ROMEO.jl’s γ (src/utility.jl): one wrap only, x ± 2π when x is outside [−π, π]. For the difference of two wrapped phases it is the full wrap; for an unwrapped phase2 (the previous echo in individual mode) it is not, and ROMEO.jl’s weights are this one’s.

idx3d 🔒

Index into a 3D array in Fortran order (column-major, matches NIfTI)

julia_max0 🔒

Julia’s max(0, x): NaN stays NaN (Rust’s f64::max(0.0) would return 0).

merge_regions_post 🔒

Merge adjacent regions by computing inter-region phase offsets.

percentile_95 🔒

Compute 95th percentile of magnitude within mask (matching ROMEO.jl).

phase_linearity_edge 🔒

Phase linearity for an edge between two voxels.

phase_linearity_triplet 🔒

Phase linearity of three consecutive phase values.

rescale_weight 🔒

Rescale a weight in [0, 1] (1 best) to the queue priority, as ROMEO.jl bins it.

temporal_uncertain_rewrap 🔒

Re-unwrap voxels that are uncertain after temporal unwrapping.

unwrap_individual 🔒

Individual multi-echo unwrapping: ROMEO.jl’s unwrap_individual! (src/unwrapping.jl, v1.6.0).

unwrap_romeo

Unwrap a single 3D phase volume using ROMEO.

unwrap_romeo_multi_echo

Unwrap 4D multi-echo phase data using ROMEO.

unwrap_template 🔒

Template-based multi-echo unwrapping.

unwrap_voxel 🔒

Temporal unwrap: remove 2π wraps by comparing to a reference value.

voxel_quality_romeo

Calculate per-voxel quality map from ROMEO edge weights

wrap_angle 🔒

Wrap angle to [-π, π]