Module lbv 

Laplacian Boundary Value (LBV) background field removal

LBV removes the background field by solving the Laplace equation inside the mask: the local field fL satisfies ∇²fL = ∇²fT on the mask’s interior and is 0 on its boundary layer, so the background fT − fL is harmonic there and matches the total field fT on the boundary.

Port of LBV.m and mexMGv6

This module is a port of the authors’ code: functions/_LBV/LBV.m and its C++ solver mexMGv6.cpp (Dong Zhou, 2013) in the MEDI toolbox, release 2024.11.26. The solver is a full multigrid (FMG) with Gauss–Seidel relaxation, ported step for step, including its quirks:

• Grids. Level 0 is the input grid; level d+1 has (1 + n_d) / 2 voxels along each axis and samples the level-0 grid at 2^(d+1) × the index (MG::initializeMem). Its mask samples the (peeled) level-0 mask the same way (MG::initializeROI), and every level’s interior is its mask minus the voxels on the array edge and those with a 6-neighbour outside the mask (MG::createInterior). Only interior voxels are ever updated.
• Depth. depth is the coarsest level. Automatic (depth −1 in LBV.m, None here): ⌊log2 m⌋ − 3 for the smallest dimension m, at least 0 (MG::determineDepth). mexMGv6 holds at most 10 levels, so depth ≤ 8.
• Right-hand side. Each level’s source ρ_d is the 7-point Laplacian of the total field sampled onto that level, with the finest voxel size on every level (MG::initializeRho; the coarse grids are not rescaled by 2^d).
• Relaxation (MG::relax): lexicographic Gauss–Seidel in place, v ← (Σ_axes h⁻²(v₊ + v₋) − r) / (2 Σ h⁻²). On the coarsest level it ignores the sweep count and repeats until the summed absolute update over the interior is at most tol; the check runs before the first sweep with the update set to 1, so tol ≥ 1 skips the coarsest solve. tol is absolute, in the field’s units: LBV.m works on MEDI’s iFreq, radians accumulated over the echo spacing, so lbv converts the field from ppm to those units first (LbvParams::te) and back at the end.
• FMG (MG::FMG): solve the coarsest level, then on each finer level prolong the coarser solution (MG::C2F, no correction) and run one V-cycle on L v = ρ_d.
• V-cycle (MG::V_cycle): N1 relaxations; the residual f − L v on the interior, injected onto the next level (MG::updateRes, MG::F2C: the even-indexed samples, no full weighting); zero the coarse solution; recurse; add the prolonged coarse correction scaled by 4^−(d+1) for coarse level d+1 (the absolute level, MG::C2F’s prefac); N2 relaxations. On the coarser levels of a V-cycle the right-hand side is the residual buffer itself, which updateRes overwrites with the new residual before the recursion, so their post-relaxation runs on that new residual. Kept as in the original.
• Prolongation (MG::C2F): every coarse mask voxel is copied to its fine position, and the fine voxels between it and its +x/+y/+z neighbours (edges, face diagonals, the body diagonal) get the mean of the coarse voxels around them when those are all in the coarse mask. These writes reach the fine mask’s boundary layer and a few voxels just outside the mask, which the relaxation never resets, so the returned field is non-zero there too (in silico 1723 boundary voxels, 1.1 % of the norm inside the mask; in vivo 1192, and 64 outside the mask, 2.3 %). LBV.m returns them and MEDI’s pipeline passes that field to MEDI_L1 with the full mask, so they are kept; the returned mask is LBV.m’s interior.
• Finish (MG::postProcessing): N3 relaxations on level 0 (to convergence if depth is 0).
• Peel (MG::preProcessing): peel more layers are removed from the mask before the levels are built, so the interior is the mask minus peel + 1 layers.

The authors’ defaults are LbvParams::default(): tol 0.01, automatic depth, peel 0 and N1/N2/N3 = 30/100/100 (LBV.m:29-46). MEDI’s pipeline script has LBV as a commented alternative to PDF (README.m:129) with tol 0.005; Ashley’s decision (2026-10-30) was to keep the function’s own 0.01.

Precision

mexMGv6 computes in single precision (float), with the prolongation’s averages and the coarse-grid correction evaluated in double and rounded (MG::C2F). LbvPrecision::Single reproduces that arithmetic operation for operation. The default, LbvPrecision::Double, runs the same algorithm in double precision, which does no harm and avoids single precision’s rounding in the relaxation (QSM.rs keeps double where the reference uses single).

Parity

LbvPrecision::Single against the compiled mexMGv6 (MEDI 2024.11.26’s mexMGv6.cpp built with R2023b’s mex, as QSM-CI’s verified lbv-cornell container runs it), on the same iFreq: bit-identical in the local field and the interior, at LBV.m’s defaults, on QSM-CI’s in-silico head (164×205×205, 7 T) and an EMTP travelling-head scan (184×256×128, 3 T), and on tests/data/lbv_cornell_small.bin (tests/lbv_cornell_parity.rs; made with the container by tests/data/lbv_cornell_small.py), which also covers peel, explicit depths, another tol and other relaxation counts. The default, double precision, against the same mex: NRMSE 0.00077 % in silico and 0.0011 % in vivo (max |Δ| 5.9e-6 and 4.4e-6 rad). QSM-CI’s lbv pair (QSMxT, default double precision, run 41423039167): NRMSE 7.7e-4 % in silico and 0.0018 % on KKI (max |Δ| 3.9e-7 and 2.0e-6 ppm).

Not ported: with peel > 0 mexMGv6 also writes the peeled mask to mask_p<peel>.bin in the working directory (and crashes MATLAB where it cannot write there).

The previous solver

Until astewartau/QSM.rs#406 this module solved the boundary value problem by successive over-relaxation (ω 1.5, to a relative tolerance, at most min(3·max dim, 500) sweeps). That is a different solver, and it stopped elsewhere: against mexMGv6 it measured 7.72 % NRMSE in silico and 24.0 % on KKI (QSMxT/QSM-CI parity run 41384710780). It was not kept as an option, because it was not better: against the exact solution of the discrete problem (conjugate gradients to 1e-12), the multigrid at LBV.m’s defaults is off by 2.4 % in silico and 8.5 % in vivo (NRMSE over the interior), the SOR solver by 7.0 % and 16.7 %; against the in-silico head’s true local field (interior, mean removed) all three are at LBV’s own limit (exact solution 43.5 %, multigrid 43.4 %, SOR 42.4 %).

Reference: 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: MEDI toolbox 2024.11.26, functions/_LBV/LBV.m and mexMGv6.cpp (http://pre.weill.cornell.edu/mri/pages/qsm.html).

Structs

LbvParams

LBV parameters (LBV.m’s defaults)

Level 🔒

Level d of the hierarchy (MG’s per-depth arrays).

Mg 🔒

Enums

LbvPrecision

Arithmetic of the multigrid solver.

Rhs 🔒

Which buffer is the right-hand side of a relaxation or V-cycle on a level.

Constants

MAX_COARSEST_SWEEPS 🔒

Guard against a coarsest-level relaxation that never meets tol (e.g. tol 0 in single precision, where the update can cycle at the last bit). mexMGv6 has no limit; at its default tol the coarsest level converges in far fewer sweeps.

MAX_DEPTH 🔒

mexMGv6’s arrays hold 10 levels, and initializeMem writes the size of level depth + 1.

Traits

Real 🔒

The arithmetic mexMGv6’s steps need, in f32 (its own) or f64.

Functions

auto_depth 🔒

MG::determineDepth: ⌊log2⌋ of the smallest dimension, minus 3, at least 0.

c2f_put 🔒

One value MG::C2F writes at fine index t: the mean sum · w of the coarse values around it (w = 1/2, 1/4, 1/8), formed in double as the C++ does (/2.0 is a double literal). As a correction it is scaled by prefac and added to the fine value in double, then rounded.

create_interior 🔒

MG::createInterior: the mask minus voxels on the array edge or with a 6-neighbour outside.

lbv

LBV background field removal (LBV.m)

lbv_multigrid

mexMGv6 on a field already in the units tol is given in (MEDI’s iFreq), as [fL, interior] = LBV(iFreq, Mask, matrix_size, voxel_size, tol, depth, peel, N1, N2, N3). params.b0, echo_times and te are not used; the other parameters are not validated (see LbvParams::validate).