Module hdqsm 

HD-QSM: Hybrid data-fidelity two-stage linear dipole inversion.

A linear dipole-inversion method that runs in two stages. Stage 1 solves an L1 data-fidelity problem (robust to phase/model errors) and derives a spatially-varying discrepancy weighting. Stage 2 solves an L2 data-fidelity problem reweighted by that discrepancy map. Both stages use ADMM with an L1 total-variation regularizer.

Data weight

w is the magnitude normalised as in HDQSM.m’s usage example, mask.*mag/max(mag), or the mask when there is no magnitude (HDQSM.m: “Not used if not specified”). The magnitude is the paper’s Eq. 2 echo combination, sum(TE_i*mag_i^2) / sum(TE_i*mag_i) (hdqsm_magnitude), which the QSM-CI original (HD-QSM c456a94) uses too; the pipeline makes it from every echo’s magnitude and requires them (QSM.rs#403: QSMxT gave the port the root sum of squares, part of a 9.1 % / 7.1 % gap on sim / KKI).

Input

The field is masked first, as in HDQSM.m’s usage example (params.input = mask.*phase_use/phase_scale) and the QSM-CI original. It matters where a local field is non-zero outside the mask (stage 1’s L1 term and the stage-2 discrepancy see it).

Dipole kernel

HDQSM.m takes the kernel as an input and ships none; it is “based on” FANSI, and the original is run with FANSI’s dipole_kernel_angulated (mode 1, discrete). That is the default here (HdQsmParams::kernel, FansiKernel has the decision and the measurements); the continuous kernel is an option.

Finite differences

As in HDQSM.m, the gradients are unit-spacing forward differences whatever the voxel size; only the dipole kernel uses the voxel size.

Units

Like the other inversions, HD-QSM takes the local field in ppm and returns ppm. HDQSM.m’s usage example passes phase/phase_scale, and the authors’ scripts on the 2019 challenge data divide the frequency map by gamma*B0, i.e. they work in ppm, so no conversion is made. Stage 2 is not scale invariant, so alpha_l2 is a value for fields in ppm.

Reference: Lambert, M., Tejos, C., Langkammer, C., Milovic, C. (2022). “Hybrid data fidelity term approach for quantitative susceptibility mapping.” Magnetic Resonance in Medicine, 88(2):962-972. https://doi.org/10.1002/mrm.29218

Reference implementation: HDQSM.m, https://github.com/mglambert/HD-QSM (MIT).

Reference pipeline

Lambert et al. 2022, Section 3.3 (in vivo): Siemens 3T, 6 echoes, TE1 = spacing = 4.92 ms, 0.8x0.8x2 mm; SEGUE unwrapping, PDF background removal, a magnitude-weighted least-squares echo fit for the local field, and the residual background removed with WH-QSM’s harmonic estimate. For the 2019-challenge simulations the field was the magnitude-weighted least-squares fit, zero-padded to 256^3, with 300 iterations (Section 3.2); this port does not pad. HDQSM.m’s header says tol_update defaults to 1.0, its code to 1e-2 (L68-72); the code value is used here.

Parity

Against the verified original (QSM-CI pair hdqsm; FANSI-toolbox 93044a3 / HD-QSM c456a94, 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.0015 %, KKI 0.00086 %, from 9.1 % / 7.1 % (sim / KKI) before. Three causes: the kernel, the Eq. 2 magnitude (QSMxT gave the root sum of squares) and the masked input. What remains is single (MATLAB) vs double (here) precision (HD-QSM works in ppm, so no γ). 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): hdqsm.bg-pdf.

Structs

HdQsmParams

HD-QSM algorithm parameters.

Functions

apply_dipole 🔒

Forward dipole convolution: Dx = real(ifft(kernel .* fft(x))).

hdqsm

HD-QSM dipole inversion.

hdqsm_magnitude

HD-QSM’s magnitude, the paper’s Eq. 2 echo combination (Lambert et al. 2022), as the QSM-CI original computes it: mag_use = sum(TE.*mag.^2, 4) ./ sum(TE.*mag, 4), with non-finite values (no signal in any echo) set to 0.

norm2 🔒

L2 norm of a slice.

percent_update 🔒

Percent update 100 * ||x - x_prev|| / ||x||.