Module tv 

Total Variation (TV) regularized dipole inversion using ADMM

Solves the L1-regularized inverse problem: min_x ||Dx - f||_2^2 + lambda||grad(x)||_1

using Alternating Direction Method of Multipliers (ADMM).

Reference: Bilgic, B., Fan, A.P., Polimeni, J.R., et al. (2014). “Fast quantitative susceptibility mapping with L1-regularization and automatic parameter selection.” Magnetic Resonance in Medicine, 72(5):1444-1459. https://doi.org/10.1002/mrm.25029

Reference implementation: https://github.com/kamesy/QSM.jl (tv), whose defaults these are: lambda = 1e-3, rho = 100*lambda, tol = 1e-3, 250 iterations. Like QSM.jl, the field is in ppm and the first x-update starts from the scaled split variable set to the gradient of the field. QSM.jl’s optional data and TV weights (W, Wtv) are not ported; their default (nothing) is what this function computes.

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).

QSM.jl pads the field to its fast FFT size first (pad = (0, 0, 0), padfastfft), and so does this port by default (TvParams::pad).

Half spectrum (tilted B0)

QSM.jl applies its k-space multipliers on FFTW’s real-to-complex half spectrum and inverts with irfft. With B0 tilted from the grid axes the dipole kernel is not even on the Nyquist planes of an even axis (D(k) ≠ D(−k) there, since the index of −½ maps to itself), and irfft resolves that differently from a full complex FFT and its real part. This port uses the multiplier irfft implies ([crate::fft::c2r_kernel]/[crate::fft::c2r_spectrum], astewartau/QSM.rs#329), which reproduces QSM.jl there too. Before (2026-10-29) the full-FFT port differed from QSM.jl on the oblique KKI travelling head (Philips Achieva, B0 tilted) by 0.5–1.2 % NRMSE for TKD, TSVD, Tikhonov and RTS (QSMxT/QSM-CI parity workflow); tests/inversion_qsmjl_parity.rs now has even-sized oblique cases.

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: NRMSE 3.4e-14 % at QSM.jl’s defaults. Before this port padded, its defaults differed from QSM.jl’s by NRMSE 6.9 %. tests/inversion_qsmjl_parity.rs checks it against QSM.jl’s output on a small phantom (relative L2 ~1e-15). As in super::rts, the DC term of the Laplacian is exactly 0 here; on anisotropic voxels QSM.jl inverts its rounding, which shifts χ by an arbitrary constant.

Reference pipeline

QSM.jl’s README pipeline (commit fde4b92, L35-52: BET on the last echo, 5-voxel erosion, Laplacian unwrapping, V-SHARP), as for super::rts.

Structs

TvParams

TV-ADMM algorithm parameters

Functions

tv_admm

TV-ADMM dipole inversion

tv_admm_unpadded 🔒

tv_admm on the already padded grid.