Module laplacian 

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.

This module provides two algorithms that are not interchangeable, because the boundary condition used to solve the Poisson equation decides whether the harmonic (background) component of the field survives:

Function	Boundary condition	Category
laplacian_unwrap	Neumann, on the array	Phase unwrapping
laplacian_unwrap_bfr (deprecated)	∇² masked to the ROI	Phase unwrapping + background field removal

laplacian_unwrap_bfr zeroes ∇²φ outside the mask, which discards every field source outside the ROI. Background fields are harmonic inside the ROI, and ∇²(harmonic) = 0 carries no information about them, so they cannot be recovered afterwards — the function returns a partially background-removed field, not a total field. That is the same combination HARPERELLA and iHARPERELLA perform, and is why it is categorised with them rather than with ROMEO.

Pair laplacian_unwrap_bfr with a separate background-removal stage only deliberately: doing so removes background twice, by an amount that is not controlled.

The background removal works by zeroing ∇² outside the mask — which deletes the exterior sources that generate the background, since a field produced outside the ROI is harmonic inside it — and then solving under a homogeneous Dirichlet condition on the ROI. On the project’s test data it reaches r = 0.887 against the ground-truth local field, against r = 0.909 for [crate::bgremove::lbv] on the same field and 0.879 for V-SHARP.

UnwrapMethod::Laplacian selects laplacian_unwrap, since the pipeline removes background as a later stage; reach for laplacian_unwrap_bfr when you want the two together.

References

Laplacian unwrapping: 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

The background-removal half of laplacian_unwrap_bfr (solving the Laplacian as a boundary value problem on the ROI): Zhou, D., Liu, T., Spincemaille, P., Wang, Y. (2014). “Background field removal by solving the Laplacian boundary value problem.” NMR in Biomedicine, 27(3):312-319. https://doi.org/10.1002/nbm.3064

Reference implementation: https://github.com/kamesy/QSM.jl — its unwrap_laplacian exposes the same split through its solver keyword (:dct/:fft impose the boundary condition on the array and unwrap only; :mgpcg imposes it on the ROI and also removes the harmonic background).

Both functions here were cross-checked against it by running QSM.jl v0.5.4 on byte-identical input (a wrapped harmonic ramp plus a non-harmonic blob, 64³): laplacian_unwrap reproduces :dct exactly (r = 1.000000, rms difference 0.0), and laplacian_unwrap_bfr matches :mgpcg to r = 0.999983 — the residual being Gauss-Seidel against their multigrid-preconditioned CG on the same equation. Both implementations return the harmonic component as zero and the non-harmonic component at r > 0.9999.

Functions

laplacian_unwrap

Laplacian phase unwrapping, without background field removal.

laplacian_unwrap_bfrDeprecated

solve_poisson_dct 🔒

Solve ∇²u = f under a Neumann boundary condition on the array, in place, via DCT-II.

solve_poisson_dirichlet_roi 🔒

Laplacian phase unwrapping combined with background field removal.

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.

wrapped_laplacian_periodic 🔒

Compute wrapped Laplacian of phase with periodic boundary conditions