Module nltv 

Nonlinear total variation (NLTV) dipole inversion

TV-regularised dipole inversion with a nonlinear data term on the wrapped phase:

min_x || W (exp(i D x) - exp(i f)) ||_2^2 + lambda || grad x ||_1

solved with ADMM. The data term gets its own auxiliary variable z = D x, whose subproblem W^2 sin(z - f) + mu (z - (D x + u)) = 0 is solved per voxel with Newton’s method.

Reference: Milovic, C., Bilgic, B., Zhao, B., Acosta-Cabronero, J., Tejos, C. (2018). “Fast nonlinear susceptibility inversion with variational regularization.” Magnetic Resonance in Medicine, 80(2):814-821. https://doi.org/10.1002/mrm.27073

Ported line for line from QSM.jl’s nltv (src/inversion/nltv.jl, https://github.com/kamesy/QSM.jl), including its defaults, its zero padding to a fast FFT size, and its starting point (z - u = grad f for the TV term).

Two places where this port has to choose, because QSM.jl’s code is not well defined:

• QSM.jl’s W is optional, but its unweighted path (W = nothing) never refreshes the FFT of x before forming D x for the Newton step, and with W = 1 and mu = 1 the Newton step’s derivative W^2/mu*cos(z - f) + 1 reaches zero: QSM.jl itself diverges on fields of more than a few radians (the qsm-forward 7T phantom). This port therefore requires the magnitude.
• QSM.jl takes the data weight W from the caller and gives no rule for making it. Here W = min(2*mask*mag/max(mag), 1), the normalisation of FANSI’s in vivo example by the same authors as the paper (and of [super::ndi]); QSM.jl squares it.
• QSM.jl allocates the data-term multiplier without initialising it; here it starts at zero.

QSM.jl uses real FFTs, this port complex FFTs. When B0 is tilted from the grid axes the dipole kernel is not symmetric on the Nyquist planes, and the two would treat those bins differently (about 1e-4 ppm RMS on a 3T scan tilted by 6 degrees); before each inverse FFT this port therefore replaces the spectrum by the Hermitian one QSM.jl’s irfft reads from its stored half ([crate::fft::c2r_spectrum], astewartau/QSM.rs#329, from 2026-10-29). Not fixture-tested (QSM.jl’s nltv is not reproducible on small arrays; see below).

QSM.jl is this method’s reference: QSM.jl, Kames’s maintained successor to QSM.m, is the implementation designated for this family (QSM.rs’s reference rule: the authors’ code where it exists, otherwise one designated implementation per family, the others’ values recorded).

Agreement with QSM.jl

On the QSM-CI in-vivo travelling-head scan (Siemens Prisma 2.89 T, 1 mm, 184×256×128, V-SHARP local field and brain mask), against QSM.jl in the QSM-CI julia-qsm image (sha256:b3c048ac…, QSM.jl fde4b92), both outputs saved as float32: with QSM.jl given this port’s data weight W and the field in radians at the first echo time, NRMSE 1.8e-10 %, max |Δ| 1.2e-10 ppm (195 iterations here). QSM.jl’s nltv cannot be run on small arrays for a fixture test: its data-term multiplier ue is allocated with similar and never initialised, and on small arrays (which reuse freed memory) its output is NaN and differs between runs. The fresh, zeroed allocations of the in-vivo arrays start it at zero, as this port does.

Units

Like the other inversions, NLTV takes the local field in ppm and returns ppm. The nonlinear data term needs the phase in radians, so the field is converted with 2*pi*gamma*b0*te first and the result converted back at the end. b0 and echo_times are scan parameters that the pipeline fills from metadata, and NltvParams::te picks the echo time the conversion uses (by default the reference’s choice; see there).

Reference pipeline

QSM.jl’s README example (commit fde4b92, L35-52), the pipeline its author ships for all its inversions: brain mask by FSL BET on the last echo’s magnitude (-f 0.5), eroded by 5 voxels; Laplacian unwrapping of each echo (which also removes harmonic background), each echo scaled to ppm by its own TE; V-SHARP with QSM.jl’s defaults (radii 18*min(vsz) down to 2*max(vsz) mm in steps of 2*max(vsz), threshold 0.05; src/bgremove/sharp.jl L136-137). QSM.jl’s inversions take 3D or 4D (per-echo) fields; the pipeline here inverts one combined field.

Structs

NltvParams

NLTV algorithm parameters (QSM.jl nltv defaults)

Functions

nltv

NLTV dipole inversion (QSM.jl nltv)

shrink 🔒