Module tikhonov 

Tikhonov regularization for QSM

Tikhonov regularization adds an L2 penalty term to stabilize the inversion:

chi = argmin_x ||Dx - f||_2^2 + lambda||Gamma*x||_2^2

This has a closed-form solution in k-space: chi_hat = D* * f_hat / (|D|^2 + lambda|Gamma|^2)

Reference: Bilgic, B., et al. (2014). “Fast image reconstruction with L2-regularization.” Journal of Magnetic Resonance Imaging, 40(1):181-191. https://doi.org/10.1002/jmri.24365

Reference implementation: https://github.com/kamesy/QSM.jl

Reference pipeline

QSM.jl (commit fde4b92) tikh (src/inversion/direct.jl L100-155): lambda 1e-2 and reg = :identity in the code (L149-150); its docstring’s signature line says :gradient (L109) but its keyword list says :identity (L131). Same README pipeline as super::rts.

QSM.jl is this method’s reference: its authors released no code, and 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, src/inversion/direct.jl:182-184), and so does this port by default (TikhonovParams::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: identical (NRMSE 0) at QSM.jl’s defaults and with padding off on both. Before this port padded, its defaults differed from QSM.jl’s by NRMSE 1.9 %. tests/inversion_qsmjl_parity.rs checks the identity and gradient regularisers against QSM.jl’s output on small phantoms.

Structs

TikhonovParams

Tikhonov algorithm parameters

Enums

Regularization

Regularization type for Tikhonov

Functions

tikhonov

Tikhonov regularization for dipole inversion

tikhonov_unpadded 🔒

tikhonov on the already padded grid.