Expand description
Weak-Harmonic QSM (WH-QSM) dipole inversion
Jointly estimates the susceptibility map x AND a residual harmonic
background field phi_h, so that any harmonic (Laplacian-null) field
remaining after background-field removal is absorbed into phi_h instead
of corrupting the susceptibility estimate. This makes the reconstruction
robust to imperfect background-field removal.
The optimization solves a nonlinear total-variation problem via ADMM,
with an additional weak-harmonic regularization term coupling x and the
harmonic field. The data-fidelity term is nonlinear (sine of the phase),
and is solved with an inner Newton iteration.
Reference: Milovic, C., Bilgic, B., Zhao, B., Acosta-Cabronero, J., Tejos, C. (2019). “Weak-harmonic regularization for quantitative susceptibility mapping.” Magnetic Resonance in Medicine, 81(2):1399-1411. https://doi.org/10.1002/mrm.27483
Ported from FANSI’s WH_nlTV.m (https://gitlab.com/cmilovic/FANSI-toolbox).
§Units
Like the other inversions, WH-QSM takes the local field in ppm and returns ppm. WH_nlTV.m
works on the phase in radians, so the field is converted with 2*pi*gamma*b0*te first and the
result converted back at the end; left in ppm, sin(x) ~ x and the data term would not be
nonlinear. b0 and echo_times are scan parameters that the pipeline fills from metadata, and
WhQsmParams::te picks the echo time the conversion uses (by default the reference’s choice;
see there).
§Data weight
WH_nlTV.m squares params.weight (recommended: the magnitude). FANSI’s WH-QSM example
(script_DICOM_Philips.m, the only script that calls WH_nlTV) sets
weight = m0.*mask./cnoise, divided by its maximum, where m0 is the proton density
from its nonlinear T2* fit (nlT2sTV) and cnoise the complex-domain noise from its
nonlinear multi-echo field fit (nlme_L2gradient). The dipole stage has neither map. The
default (WhQsmWeight::Flat) is WH_nlTV.m’s own default when no weight is given: 1 over
the whole volume (W = ones(N), lines 70-75), with the mask given separately (it confines
x and the harmonic penalty beta), as the QSM-CI original runs it. The data term outside
the mask is part of the method: an earlier port weighted only the mask, which with the
continuous kernel left it 21.4 % / 50.7 % (sim / KKI χ NRMSE) from the verified original;
a numpy transcription of WH_nlTV.m on the original’s own inputs with the mask as the
weight diverges on sim, and with ones(N) and FANSI’s kernel it is 0.0006 % from it
(QSM.rs#402). The magnitude weights are options:
WhQsmWeight::MaxNormalised is the example’s formula with the magnitude for m0 and no
noise term, mask.*mag / max(mask.*mag), and WhQsmWeight::InVivo FANSI’s
script_invivo.m form, min(2*mask*mag/max(mag), 1), as NDI uses. How the weight and the
echo time (WhQsmParams::te) interact on the phantom is noted there.
§Newton safeguard
The inner Newton step for z2 (WH_nlTV.m l. 258-268) minimises, voxel by voxel,
W*(1 - cos(z2 - phase)) + mu2/2*(z2 - rhs/mu2)^2, whose minimiser lies within W/mu2 of
rhs/mu2 (where the derivative W*sin(z2 - phase) + mu2*(z2 - rhs/mu2) vanishes). With
the flat weight, a field that is non-zero outside the mask (QSM-CI’s simulated local
field, up to ~6 rad) puts some voxels where W*cos + mu2 nears 0, and the unguarded step
bounces: on that sim the original’s Newton happened to settle after ~24 inner steps, while
the port’s (same formula, FFT rounding apart) ran to its 50-step cap at 1e16 and the solve
went NaN. The port clamps each Newton iterate to that interval. This changes nothing where
the step stays inside it (it does wherever the original converges normally), and the fixed
point is the same minimiser; it is not in WH_nlTV.m. With it the port is 0.0006 % (sim)
and 0.0009 % (KKI, no clamping needed) from the verified original (χ in radians, QSM.rs#402).
§Finite differences
As in WH_nlTV.m, the gradient and Laplacian are unit-spacing differences whatever the
voxel size; only the dipole kernel uses it.
§Reference pipeline
FANSI script_DICOM_Philips.m (commit 93044a3) runs WH_nlTV on the same local field and
TE scaling as nlTV (see super::fansi: SEGUE, multi-echo fit, MSMV background removal,
deltaTE with “any individual echo time” allowed, L364), with mask = mask2, beta = 200,
muh = 5 and the advice that beta is “usually … 20 to 50 times larger than muh” and that
alpha1 may be slightly lower than FANSI’s (L425-436). This port uses WH_nlTV.m’s own
defaults for beta (150) and muh (5), not the script’s beta of 200. Without params.mask,
WH_nlTV.m takes the mask as W > 0; here the mask is always given.
§Parity
Against the verified original (QSM-CI pair whqsm; FANSI-toolbox 93044a3, identical to
licensed MATLAB R2023b), final-map χ NRMSE at both sides’ defaults on QSM-CI’s simulated head
and the KKI travelling head (QSM-CI parity run 41403629314, QSM.rs#399–#403): sim 0.0018 %, KKI 0.0038 %, from
21.4 % / 50.7 % (sim / KKI) before. Two causes: the kernel, and WH_nlTV.m’s flat weight over the whole volume (see Data weight); the Newton safeguard keeps the sim run from diverging (see Newton safeguard). What remains is
single (MATLAB) vs double (here) precision and FANSI’s γ, 42.58 against the crate’s 42.576
(the wrappers keep theirs). Fixture: tests/fansi_qsmci_parity.rs (relative L2 below 1e-5
on an oblique, anisotropic phantom).
§Pipeline rules
The facts above that a configuration can depart from are published as rules in
crate::pipeline::rules::RULES (each with its strength, basis, reason and source): fansi.no-laplacian-multi-echo, whqsm.flat-weight-te-spacing, whqsm.magnitude-weight-te-spacing.
Structs§
- WhQsm
Params - WH-QSM parameters.
Enums§
- WhQsm
Weight - The data weight (see the module docs).