Expand description
L1-QSM: nonlinear L1 data-fidelity dipole inversion with TV regularization
Solves the nonlinear QSM inverse problem with an L1 data-fidelity term (robust to phase outliers / noise) and a total-variation prior:
min_x || W .* ( exp(i * D x) - exp(i * phi) ) ||_1 + alpha1 * ||grad(x)||_1
where D is the dipole operator and phi the local field phase in radians.
The nonlinear fidelity is linearized around the current auxiliary phase and
solved with an inner complex-argument Newton iteration, all embedded in an
ADMM splitting.
Reference: Milovic C, Lambert M, Langkammer C, Bredies K, Irarrazaval P, Tejos C. “Streaking artifact suppression of quantitative susceptibility mapping reconstructions via L1-norm data fidelity optimization (L1-QSM).” Magnetic Resonance in Medicine, 2022;87:457-473. https://doi.org/10.1002/mrm.28957
Ported from the FANSI toolbox nlL1TV.m
(https://gitlab.com/cmilovic/FANSI-toolbox).
§Units and weights
Like the other inversions, L1-QSM takes the local field in ppm and returns ppm.
nlL1TV.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. The
scale matters here: the L1 proximal step thresholds |exp(i*Dx) - exp(i*phi)|
at W/mu3, so in ppm every residual falls below it and the data are fitted
exactly. b0 and echo_times are scan parameters
that the pipeline fills from metadata, and L1QsmParams::te picks the echo time the
conversion uses (by default the reference’s choice; see there).
W = lambda * min(2*mask*mag/max(mag), 1): the magnitude in [0, 1] as in
FANSI’s in vivo example, scaled by lambda (nlL1TV.m: “lambdamagnitude_data
or lambdamask”). The magnitude is required: with a flat (mask) weight the
nonlinear data term can diverge where the phase nears +-pi.
The finite differences are FANSI’s, on unit spacing whatever the voxel size.
§Reference pipeline
FANSI script_DICOM_Philips.m (commit 93044a3) runs nlL1TV on the same local field and
TE scaling as nlTV (see super::fansi), with weight lambda*m0.*mask./cnoise / max,
lambda = 0.8 (L443-450), and says the regularisation weight must be re-tuned for L1-QSM,
not reused from FANSI or WH-QSM (L442-443). nlL1TV.m defaults: mu1 = 100*alpha1, mu2 = mu3
= 1, 50 iterations, tolUpdate 1, weight ones (L48-112).
§Parity
Against the verified original (QSM-CI pair l1qsm; 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.0022 %, KKI 0.0026 %, from
14.4 % / 15.5 % (sim / KKI) before. The kernel was the whole gap (a numpy nlL1TV.m with FANSI’s kernel is 0.0004 % from the original on sim; the isPrecond start and the 0.05 floor already matched). 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, l1qsm.te-first.
Structs§
- L1Qsm
Params - L1-QSM algorithm parameters (FANSI
nlL1TV.mdefaults where it has them).