Module tkd 

Truncated k-space division (TKD) / TSVD for QSM

TKD is the simplest dipole inversion method. It directly divides the field in k-space by the dipole kernel, with truncation to avoid division by small values near the magic angle.

Reference: Shmueli, K., de Zwart, J.A., van Gelderen, P., Li, T.Q., Dodd, S.J., Duyn, J.H. (2009). “Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data.” Magnetic Resonance in Medicine, 62(6):1510-1522. https://doi.org/10.1002/mrm.22135

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

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 (TkdParams::pad).

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) for tkd and tsvd at QSM.jl’s defaults and with padding off on both. Before this port padded, its defaults differed from QSM.jl’s by NRMSE 7.5 % (TKD) and 15.6 % (TSVD). tests/inversion_qsmjl_parity.rs checks both against QSM.jl’s output on small phantoms, including tilted B0 on even sizes.

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.

Reference pipeline

QSM.jl (commit fde4b92) tkd/tsvd (src/inversion/direct.jl L1-96), threshold 0.15 for both: TKD replaces 1/D by sign(D)/thr where |D| <= thr, TSVD by 0. QSM.jl’s README pipeline (L35-52: BET on the last echo, 5-voxel erosion, Laplacian unwrapping, V-SHARP) is the one its inversions are shown with; see super::rts.

Structs

TkdParams

TKD algorithm parameters

Functions

tkd

Truncated k-space division (TKD) for dipole inversion

tkd_unpadded 🔒

tkd on the already padded grid.

tsvd

Truncated singular value decomposition (TSVD) variant

tsvd_unpadded 🔒

tsvd on the already padded grid.