Module ndi 

Nonlinear Dipole Inversion (NDI) for QSM

Solves the dipole inversion problem with a nonlinear data-consistency term, modelling the wrapped phase directly via sin(Dx - phase) and minimizing it with a simple L2-regularized gradient descent:

min_x || W .* (exp(iDx) - exp(iphase)) ||^2 + alpha ||x||^2

which yields the update

x <- x - tau * D^H(W^2 .* sin(Dx - phase)) - taualphax

W is the data weight: the magnitude, normalised as in FANSI’s in vivo example (min(2*mask*mag/max(mag), 1)). The magnitude is required: with a uniform weight the update can diverge where the phase nears ±π.

Because the dipole kernel D is real, conj(D) = D and the adjoint is just another forward dipole application.

Reference: Polak, D., Chatnuntawech, I., Yoon, J., Iyer, S.S., Milovic, C., Lee, J., Bachert, P., Adalsteinsson, E., Setsompop, K., Bilgic, B. (2020). “Nonlinear dipole inversion (NDI) enables robust quantitative susceptibility mapping (QSM).” NMR in Biomedicine, 33(12):e4271. https://doi.org/10.1002/nbm.4271

Ported from FANSI’s ndi.m, and ndi_auto.m for the automatic stopping criterion (Milovic et al., ISMRM 2021, p. 3982) (https://gitlab.com/cmilovic/FANSI-toolbox).

Units

Like the other inversions, NDI takes the local field in ppm and returns ppm. FANSI’s solver 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 NdiParams::te picks the echo time the conversion uses (by default the reference’s choice; see there).

Reference pipeline

• Paper (Polak et al. 2020, Methods): 3T data at 1 mm isotropic, single TE 25 ms, and 7T Wave-CAIPI at 0.5 mm isotropic, TE 19 ms; BET brain masking, Laplacian unwrapping and SMV background filtering; a small Tikhonov term (“lambda = 0.1%”) to stop over-fitting.
• FANSI script_DICOM_Philips.m (commit 93044a3) runs ndi at the end of the same chain as nlTV (see super::fansi for its masking, SEGUE unwrapping, multi-echo fit and MSMV background removal): input local*deltaTE and weight m0.*mask./cnoise, normalised (L457-461), 300 iterations (L463). ndi.m itself defaults the weight to ones (L51-56).
• The defaults here are ndi.m’s (tau 2, alpha 1e-5, 100 iterations).

Parity

Against the verified original (QSM-CI pair ndi; 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.0011 %, KKI 0.00009 %, from 10.6 % / 9.5 % (sim / KKI) before. The kernel was the whole gap. 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, ndi.te-first.

Structs

NdiParams

NDI algorithm parameters

Functions

apply_dipole 🔒

Apply the forward dipole operator D * x in-place buffers.

auto_stop_masks 🔒

The two k-space bands ndi_auto.m compares (FANSI create_freqmasks_auto.m): a radial shell of high frequencies, split by |D| into a band close to the cone (0.15 < |D| <= 0.2) and one further off it (0.225 < |D| <= 0.275). k is the kernel in FFT order.

ndi

Nonlinear Dipole Inversion (NDI)