fn solve_poisson_dirichlet_roi(
f: &[f64],
mask: &[u8],
grid: &Grid,
tol: f64,
max_iter: usize,
) -> Vec<f64>Expand description
Laplacian phase unwrapping combined with background field removal.
Solves the Poisson equation with the Laplacian zeroed outside mask. That discards the
field sources outside the ROI, and a field generated outside the ROI is harmonic inside
it — so the background component is removed along with the wraps.
The result is not a total field. It is unwrapped and partially background-removed, by an amount that depends on the mask and the field geometry. Following this with a separate background-removal stage (V-SHARP, PDF, …) removes background twice.
Use laplacian_unwrap to unwrap without removing background.
Because ∇²(harmonic) = 0, the discarded component leaves no trace in the input to the Poisson solve and cannot be restored afterwards.
§Intended input
One wrapped phase volume. On the project’s test data it then matches
[crate::bgremove::lbv] on the same field (r = 0.887 against 0.909, identical residual
smooth content). It takes wrapped phase, so in a multi-echo pipeline the only way to
apply it is to each echo before combining; that usage is not what the algorithm
describes, has not been validated, and leaves visibly more background than either
unwrapping then a field-map background removal or this function on a single volume.
For multi-echo data use laplacian_unwrap or ROMEO, combine, then a background
removal from crate::bgremove.
§Echo time
Accuracy falls off with the amount of phase to unwrap. On the 7 T test data, against
the ground-truth local field: r = 0.84 at TE = 4 ms, 0.70 at 8 ms, 0.50 at 12 ms —
where unwrapping then [crate::bgremove::lbv] gives 0.86, 0.84, 0.69 and
[crate::bgremove::vsharp] holds near 0.82 throughout. Prefer the earliest echo.
§Arguments
phase- Wrapped phase (nx * ny * nz)mask- Binary mask (nx * ny * nz), 1 = inside ROIgrid- Volume grid (dimensions and voxel sizes)
§Returns
Unwrapped, partially background-removed phase, zero outside mask.
§References
Schofield & Zhu (2003) for the unwrapping; Zhou et al. (2014) for the
boundary-value formulation of the background removal. See the module docs.
Solve ∇²u = f inside mask with u = 0 outside it (homogeneous Dirichlet on the ROI),
by Gauss-Seidel with successive over-relaxation.
Masking the source term is only half of the ROI formulation: the solution has to be constrained at the ROI boundary too, or it picks up an arbitrary harmonic component. Solving the masked source over the whole volume with a periodic FFT does not constrain it, and that component is large — which is what this replaces.
Mirrors the solver in [crate::bgremove::lbv], which solves the homogeneous case
(f = 0) with boundary values taken from the field.