Module laplacian

Expand description

Laplacian-based phase unwrapping

The Laplacian of the wrapped phase equals the Laplacian of the true phase wherever neighbouring samples differ by less than π, so the true phase can be recovered by solving a Poisson equation. Path-independent and fast, unlike region-growing methods.

laplacian_unwrap unwraps only; the harmonic (background) component survives and the result is a total field. It has two solvers, named as in QSM.jl’s unwrap_laplacian:

`LaplacianSolver`	∇²φ estimate	Poisson solve	Reference
`Dct` (default)	finite differences of wrapped neighbour differences	DCT, Neumann boundary on the array	Ghiglia & Romero (1994); QSM.jl `:dct`
`Fft { pad }`	`cos φ·∇²(sin φ) − sin φ·∇²(cos φ)`, spectral ∇² (`	k	²`)

Dct is unweighted least-squares unwrapping: the Poisson equation for the wrapped phase differences, solved under a Neumann condition with a DCT (Ghiglia & Romero 1994). It reproduces QSM.jl’s :dct solver exactly (QSM.jl v0.5.4 on byte-identical input: r = 1.000000, rms difference 0.0). tests/laplacian_qsmjl_parity.rs keeps this checked against QSM.jl fde4b92 in the verified QSM-CI julia-qsm image (QSMxT/QSM-CI#299) on two wrapped phantoms, 1 mm and anisotropic: max |Δ| 8.9e-15 rad (measured 2026-10-30).

Through UK Biobank’s pipeline in the verified QSM-CI ukb-qsm container (QSMxT/QSM-CI#305), Fft with pad 64 per echo, combined with B0WeightType::AssumedDecay, gives the pipeline’s combined phase to 2.5e-14 rad (tests/ukb_qsm_parity.rs, two echoes at 1×1×1.5 mm, measured 2026-10-30).

Fft is the literal Schofield & Zhu method: the Laplacian of the true phase from the sin/cos identity with spectral derivatives, then an FFT Poisson solve. It reproduces STI Suite 3.0’s MRPhaseUnwrap(phase, 'voxelsize', vs, 'padsize', pad), which is this algorithm, to ≈1e-14 rad inside the mask (pad 64; in-vivo 3 T data and synthetic volumes at several voxel sizes). Steps: zero-pad by pad per side (an odd padded dimension gets one more zero plane at its end); K = |k|² on the padded FFT grid, k_d = m / (M_d·h_d), no 4π² factor; L = cos φ·F⁻¹[K·F(sin φ)] − sin φ·F⁻¹[K·F(cos φ)]; φ_u = F⁻¹[F(L)/K] with DC → 0; crop.

The two give nearly the same local field after background removal (simulated head: NRMSE 79.5 % for Fft against 81.3 % for Dct; on a 150-subject in-vivo test-retest cohort the same iron-region ICC, 0.924 against 0.925). They differ in the total field: zero padding puts a step wherever the phase at the array edge is not zero, which Fft turns into a harmonic error that Dct does not have. Hence Dct is the default.

Multi-echo field maps are built by crate::pipeline::run_field_mapping, whose laplacian_solver setting selects the solver: each echo is masked (unless laplacian_unmasked), unwrapped and combined with the configured B0WeightType. UK Biobank’s field map is Fft { pad: [64; 3] } with B0WeightType::AssumedDecay at 40 ms.

UnwrapMethod::Laplacian selects laplacian_unwrap with LaplacianSolver::Dct, since the pipeline removes background as a later stage.

§Reference pipeline (STI Suite)

Probed in MATLAB R2023b by driving STI Suite 3.0’s QSM_GUI headless with logging wrappers (3-echo synthetic head, TE 4/10/16 ms, 1×1×1.5 mm):

MRPhaseUnwrap’s default padsize is 12 voxels (bitwise equal to an explicit [12 12 12]; LaplacianSolver::FFT_DEFAULT_PAD), as STI_Templates.m says.
The GUI (“3D GRE data” and “Single button for QSM”) calls it once on all echoes, MRPhaseUnwrap(φ, 'voxelsize', vs) with the raw, unmasked wrapped phase (the input is bitwise the loaded phase, nonzero outside the mask), so the GUI does not do the mask · φ that run_field_mapping does by default and UK Biobank does (laplacian_unmasked is the GUI’s choice; rule sti.unwrap-unmasked). It then sums the unwrapped echoes (no per-echo offset removal, no weighting) and gives Σφ to V-SHARP, with TE = ΣTE for the inversion: in φ/TE form, the B0WeightType::TEs weighting (Σφₑ/ΣTEₑ).
STI Suite 2.01’s GUI also offered “Combine Laplacian” (the Laplacians of all echoes, combined) followed by iHARPERELLA (see crate::bgremove::iharperella).
QSMART does not use MRPhaseUnwrap but Hongfu Sun’s lapunwrap.m (mirror extension): LaplacianSolver::Mirror, per echo on the unmasked phase (unwrap_phase.m l. 5–10; laplacian_unmasked).

§Reference pipeline (other than STI Suite)

QSM.jl (kamesy/QSM.jl, local copy fde4b92): unwrap_laplacian’s default solver is :mgpcg, a Dirichlet condition on the ROI (uphas[!mask] = 0) that “combines phase unwrapping and harmonic background field removing”; :dct and :fft are the array-boundary options (src/unwrap/laplacian.jl:6,14-19,28-32). :mgpcg is not ported; the default here is :dct, which unwraps only. The README’s pipeline masks with FSL BET on the last echo (-f 0.5), erodes the mask by 5 voxels, unwraps every echo with unwrap_laplacian(phas, mask1, vsz) (so :mgpcg), scales echo e by 1/(B0·γ·TEₑ) and runs V-SHARP (README.md:34-49).
MriResearchTools.jl / CLEAR-SWI (korbinian90, 771eadf / 0949bc5): laplacianunwrap is the Schofield & Zhu sin/cos method with DCT transforms (Neumann, unpadded, no mask, no voxel size; src/laplacianunwrapping.jl:1-20), so it is neither of the solvers above exactly. CLEAR-SWI’s default unwrapping is laplacian (ext/ClearswiApp/argparse.jl:55-57), applied per echo to the unmasked phase, and it also offers laplacianslice, the same per 2D slice (src/phase_processing.jl:88-116); it does no MCPC-3D-S before combining the echoes with calculateB0_unwrapped (:phase_snr).
Nothing in the method assumes a field strength or B0 direction; Dct and Fft take the voxel size into the Laplacian, so isotropic voxels are not required.

§References

Ghiglia, D.C., Romero, L.A. (1994). “Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods.” Journal of the Optical Society of America A, 11(1):107-117. https://doi.org/10.1364/JOSAA.11.000107

Schofield, M.A., Zhu, Y. (2003). “Fast phase unwrapping algorithm for interferometric applications.” Optics Letters, 28(14):1194-1196. https://doi.org/10.1364/OL.28.001194

Li, W., Wu, B., Liu, C. (2011). “Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition.” NeuroImage, 55(4):1645-1656 (STI Suite’s use of the Schofield & Zhu method). https://doi.org/10.1016/j.neuroimage.2010.11.088

Reference implementation: https://github.com/kamesy/QSM.jl (unwrap_laplacian).

Enums§

LaplacianSolver: Poisson solver for laplacian_unwrap (names as in QSM.jl’s unwrap_laplacian).

Functions§

apply_k2 🔒: Multiply (or, with invert, divide) a spectrum by K = kx² + ky² + kz²; DC → 0 on divide.
laplacian_unwrap: Laplacian phase unwrapping, without background field removal.
laplacian_unwrap_fft_mirrored 🔒: Schofield & Zhu sin/cos unwrapping on the volume mirrored to 2n along each axis (lapunwrap.m). Returns the whole volume, unmasked.
laplacian_unwrap_fft_padded 🔒: Schofield & Zhu sin/cos unwrapping on the zero-padded volume (STI Suite’s MRPhaseUnwrap). Returns the whole volume, unmasked.
solve_poisson_dct 🔒: Solve ∇²u = f under a Neumann boundary condition on the array, in place, via DCT-II.
wrap 🔒: Wrap angle to [-π, π]
wrapped_laplacian_neumann 🔒: Wrapped Laplacian under a Neumann boundary: at each array face the missing neighbour is the sample itself, so the wrapped difference across the face is zero. This is exactly the half-sample even extension the DCT-II assumes, so pairing it with solve_poisson_dct reproduces the even-extended periodic solve without building the 2x-per-axis extension.