Module rts 

Rapid Two-Step (RTS) dipole inversion

Two steps:

1. Well-conditioned k-space: a few LSMR iterations on D χ̂ = F f (all of k-space, from zero), then the result is masked.
2. Ill-conditioned k-space: TV-regularised ADMM that holds the well-conditioned coefficients (|D| > delta) of step 1 with weight mu, starting from the gradient of the step-1 map, and raising the penalty rho by gamma whenever the primal residual falls by less than a factor tau.

Reference: Kames, C., Wiggermann, V., Rauscher, A. (2018). “Rapid two-step dipole inversion for susceptibility mapping with sparsity priors.” NeuroImage, 167:276-283. https://doi.org/10.1016/j.neuroimage.2017.11.018

Reference implementation: QSM.jl’s rts (https://github.com/kamesy/QSM.jl, by the method’s first author), ported line for line with its defaults. The field is in ppm, as in QSM.jl. QSM.jl works on real-to-complex FFTs; here the full complex FFT is used, and the LSMR inner products are taken over the half spectrum QSM.jl stores (kx index 0..=nx/2) so its scalars are the same.

QSM.jl is this method’s reference: RTS is Kames’s method and QSM.jl is his code, so this is the authors’ code (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/rts.jl:109-111), and so does this port by default (RtsParams::pad); the padding changes the result, since the filters act on the periodic extension of the array.

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 1.0e-8 % at QSM.jl’s defaults, 1.3e-9 % with padding off on both. Before this port padded, its defaults differed from QSM.jl’s by NRMSE 0.45 %. tests/inversion_qsmjl_parity.rs checks it against QSM.jl’s output on small phantoms (relative L2 ~1e-15), axial and oblique.

One deviation, where QSM.jl’s result is not well defined:

• DC of the Laplacian. On anisotropic voxels the FFT of the Laplacian stencil is not exactly 0 at DC (4.4e-16 at 0.9 × 1.0 × 1.3 mm), and QSM.jl inverts that rounding, which puts an arbitrary constant into χ (it moved χ by a uniform 1e-3 ppm on a small phantom, and changed the iteration count). Here the DC term is exactly 0, as it is on isotropic voxels.

Tilted B0 on even sizes is reproduced too: before each inverse FFT the spectrum is replaced by the Hermitian spectrum QSM.jl’s irfft reads from its stored half ([crate::fft::c2r_spectrum]; see Half spectrum below).

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’s README example (commit fde4b92, L35-52) ends in rts: FSL BET on the last echo’s magnitude (-f 0.5), eroded by 5 voxels; Laplacian unwrapping of each echo, each 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); then rts on the (per-echo) local field, with pad = (0, 0, 0), padding to QSM.jl’s fast FFT size only.

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): rts.bg-vsharp.

Structs

RtsParams

RTS algorithm parameters (QSM.jl rts keywords and defaults)

Functions

half_norm 🔒

Norm of a Hermitian-symmetric spectrum as QSM.jl sees it: over the kx = 0..=nx/2 half that its real-to-complex FFT stores.

lsmr_diag 🔒

LSMR for the diagonal operator diag(d) on a spectrum, from x = 0, with lambda = 0, atol = btol = 0 and conlim = 1e8, as QSM.jl’s rts calls its copy of IterativeSolvers.jl’s lsmr! (v0.9.2).

rts

RTS dipole inversion (QSM.jl rts)

rts_unpadded 🔒

rts on the already padded grid.