Module ilsqr 

iLSQR: LSQR dipole inversion with streaking-artefact removal

ilsqr implements STI Suite 3.0’s QSM_iLSQR (Li et al. 2015), as called with QSM_iLSQR(phase, mask, 'TE', TE, 'B0', B0, 'H', H, 'padsize', [p p p], 'voxelsize', vs) with p = IlsqrParams::pad_mm. The defaults are STI_DEFAULT_NITER = 100 iterations, QSM_iLSQR’s own default, which applies when it is called as a function without 'niter' (as UK Biobank’s pipeline and STI Suite’s STI_Templates.m call it), and STI_GUI_PAD_MM = 12 mm of padding, what STI Suite 3.0’s QSM_GUI passes (see IlsqrParams). The GUI passes STI_GUI_NITER = 40 iterations. Called without 'padsize', QSM_iLSQR pads STI_DEFAULT_PAD_MM = 6; UK Biobank and STI Suite’s examples pass 64, STI_PAD_MM. STI Suite ships as protected P-code, so the algorithm was recovered by black-box probing (call tracing, capturing the operators and right-hand sides it hands to MATLAB’s lsqr, and fitting every stage until a MATLAB re-implementation reproduced STI to its own single-precision noise floor):

1. Crop field and mask to the mask’s bounding box (odd extents grown by one voxel to make them even), then zero-pad each axis to the next multiple of 16 strictly above crop + 2·pad_mm/voxel_size (default pad_mm = 12). Everything below runs on that grid.
2. Laplacian weights W: L = |F⁻¹(k²·F φ)| with k² in (cycles/mm)² and the DC term set to 10⁶ (STI’s value, which adds 10⁶·mean(φ) to every voxel, so the |·| decides which tail of the Laplacian is down-weighted). Ranked inside the mask: W = 1 below the 60th percentile, falling linearly to 0.2 at the 99.9th, 0 outside. The field itself is not re-masked.
3. Initial solution: LSQR on D W D χ = D W φ (D = dipole convolution), tolerance 0.01, at most max_iter iterations (MATLAB lsqr stopping rules).
4. FastQSM prior: sign(D)·Fφ, blended where |D| is small with a copy smoothed in k-space by a radius-3 ball closed with IPT’s imclose(·, [1 1 1]) (which, on this ball, only fills its centre), masked, blended again, rescaled to a TKD (threshold 1/8) solution by least squares inside the mask.
5. Gradient weights G_d (d = x, y, z): forward differences of the FastQSM map (no voxel-size scaling), pooled over the three axes inside the mask; G = 1 below the 50th percentile of |∇|, 0 above the 70th, linear between.
6. Streaking artefacts: LSQR for y in min ‖G ∇ F⁻¹(M_ic y) − G ∇(M·Re χ₁)‖ with M_ic = |D| < 0.1 (the ill-conditioned cone), tolerance tol (default 1e-3), at most max_iter iterations; artefact χ_sa = Re F⁻¹(M_ic y).
7. χ = (Re χ₁ − χ_sa) · mask, cropped back to the input grid.

Percentiles are STI’s: the sorted value at 1-based index round(p·N), not interpolated. The field may be in any linear unit (rad, Hz, ppm); χ comes out in the unit the field would have for χ = 1 (QSM.rs passes ppm, so χ is in ppm).

Performance and precision

Nearly all the time goes into the two LSQR solves, each costing two 3-D FFTs on the padded grid per operator call. The initial solve runs in k-space (2 FFTs per call instead of 4, the same iterates up to rounding), the transforms are pruned to the box holding the mask and carry the solves’ diagonal products in their first and last passes, and the artefact solve stores its vectors on the cone and the mask only. On a UK Biobank SWI field (352×400×96 padded) this takes ~15 s on 16 cores with a 1.9 GB peak. Against the image-space formulation of the initial solve, the k-space one gives the same stop flags and iteration counts in both solves there, and χ within a relative L2 of 5.8e-9 (the artefact solve amplifying rounding-level differences of the initial solution, which agree to 2.7e-15).

By default everything runs in double precision. STI Suite runs its two LSQR solves in single precision: their operators (single FFTs, single(D), single(W), single(G)), and MATLAB’s lsqr on single vectors and scalars, with the right-hand sides computed in double and rounded to single (probed: STI’s initial solve is reproduced to 2e-8 that way, against 7e-4 in double). IlsqrPrecision::Single does the same, for parity with STI; it runs the two solves on dense vectors over the padded grid (no pruned transforms), so it is slower and needs more memory. Double precision stays the default (QSM.rs keeps double where the reference uses single and it does no harm): on whole heads both agree with STI equally well (QSM-CI’s sim and KKI travelling head: 0.0024 % and 0.017 % NRMSE in either precision). LSQR stopped after a fixed number of iterations amplifies rounding, though, so on small noisy volumes only single precision follows STI (below: 0.17 % against 10.4 %).

Agreement with STI Suite

tests/ilsqr_sti_parity.rs (measured 2026-10-30; re-measured 2026-11-02 after the dipole kernel’s rounding, below):

• Against the re-implementation of QSM_iLSQR the algorithm was recovered with (tests/data/ilsqr_sti_clone.m, double precision, 5 iterations per solve): max |Δ| 5.7e-14 ppm on a 20×20×14 phantom (noise, phase outside the mask). This pins the algorithm.
• Against STI’s own output on that phantom (licensed MATLAB R2023b, 100 iterations): IlsqrPrecision::Single χ NRMSE 0.17 %, initial LSQR solution 0.010 %; double precision χ 10.4 % (r 0.9950), initial 0.36 %.
• Against the verified QSM-CI ilsqr-sti container (QSMxT/QSM-CI#304) at its call (100 iterations, padsize 12 and 8) on two noise-free 32×32×24 / 32×28×28 phantoms: double NRMSE 0.12 % and 0.071 % (r 0.9999993, 0.9999997); single 0.032 % and 0.0055 % (0.019 % before the dipole kernel’s rounding), at or below STI’s own spread between CPUs (its χ moves by 3.6e-4 and 1.2e-4 relative between the workstation that made the fixture and a GitHub-hosted runner). That fixture is regenerated in Actions by tests/data/ilsqr_sti_qsmci_small.py.
• Inside QSMART (tests/qsmart_qsmci_parity.rs, 2026-11-02: QSMART’s call, 50 iterations, STI’s default padding, field in rad/s, on the QSMART toolbox’s own SDF local fields of a 36×36×28 vessel phantom): 1.0e-4 % and 1.0e-3 % NRMSE for the two stages, once QSMART’s own γ mismatch is applied (it makes the field with γ = 2.675e8 rad/s/T and STI turns it back with 2π·42.575 MHz/T, so the original’s χ is 2.2e-5 low; QSM.rs keeps one γ and does not reproduce that: 2.5e-3 % without it). That is STI’s single-precision floor (its own χ moves by 1.3e-4 % / 1.3e-3 % between CPUs there).

Dipole kernel rounding (astewartau/QSM.rs#396, 2026-11-02). Until then the kernel’s k-space coordinates were [crate::fft::fftfreq]’s, i·(1/(n·vs)); STI computes i/(n·vs). The two round differently, and on the lattice points of the magic angle (k_x² + k_y² = 2k_z² for H = z) that decides whether D is 0 or −5.6e-17. FastQSM takes sign(D), so the FastQSM prior and from it the gradient weights moved (at ~1,600 voxels per axis on the QSMART phantom), while the initial solution agreed to 1e-10 %. Inside QSMART (48³ grid, 64 of its 266 magic-angle points differing in sign) χ was 0.15 % and 0.094 % NRMSE from STI’s for the two stages, in either precision; found by dumping the MATLAB re-implementation’s intermediates on QSMART’s call (which reproduces STI there to 2.6e-5 % / 5.7e-4 %; with fftfreq’s coordinates it gives the same 0.15 % / 0.094 %). Neither sign pattern is better than the other (both are rounding), so the kernel follows STI’s arithmetic (sti_dipole_kernel).

Until 2026-10-30 the FastQSM kernel was closed without IPT’s padding (−∞ beyond the array for the dilation, +∞ for the erosion, on the unpadded 7³ ball), which kept 32 dilated voxels at the ball’s ±y caps (179 non-zeros against imclose’s 147). FastQSM moved by 2 %, its gradient weights at ~800 voxels per axis, and χ by 3.4 % / 1.6 % NRMSE on the container phantoms, 0.99 % on QSM-CI’s sim and 1.78 % on the KKI travelling head; with the closing fixed, 0.0024 % and 0.017 % (QSMxT/QSM-CI parity runs 41019608126 → 41023156365, astewartau/QSM.rs#329). QSMART, which calls this iLSQR, moved from 16.71 % to 16.64 % on sim against the original QSMART (4.27 % on KKI). The licensed fixture above had been made with an imclose stand-in of the same kind, which hid it; it was regenerated with MATLAB’s own imclose.

Padding

IlsqrParams::pad_mm is STI’s 'padsize': millimetres of zero padding per side, the same on every axis. It sets the grid that both LSQR solves run on, so for small or thin volumes it dominates the cost: an 8³ volume at 1 mm is solved on 144³ with 64 mm, on 32³ with 6 mm. The default, STI_GUI_PAD_MM (12), is what STI Suite’s QSM_GUI passes (and what STI_Templates.m and the χ-separation toolbox’s Chisep_script.m use). STI_DEFAULT_PAD_MM (6) is what QSM_iLSQR uses when called without 'padsize' (as QSMART calls it). STI_PAD_MM (64) is what UK Biobank’s pipeline and STI Suite’s examples pass (set it to reproduce them exactly). 0 is allowed (the grid is still rounded up past the crop to a multiple of 16) but wraps thin volumes onto themselves. The padding changes χ because the dipole convolution is circular on the padded grid and the Laplacian and gradient weights are ranked over it.

On 150 UK Biobank subjects × 2 visits, the qsm-forward head phantom and thin slabs, comparing 0–128 mm (ukb-qsmxt docs/19), no test-retest, age or longitudinal metric differed from 64 mm. Above 32 mm the result only moves by the floor set by the capped LSQR iterations (32 and 128 mm are each ~2 ppb RMS from 64 mm). On a thin phantom slab 0–6 mm lose ~1.5 NRMSE points against the ground truth and 16 mm ~0.5; on whole heads every pad gave the same accuracy.

QSM.rs ≤ v0.38 shipped a different iLSQR (a port of QSM.m’s formulation: no padding, LSMR for the artefacts, other weights). It was removed: STI Suite’s algorithm is the reference iLSQR, and it is also what the original QSMART calls (see crate::pipeline::run_qsmart).

Reference pipeline

What STI Suite 3.0 runs around QSM_iLSQR. Its GUI and functions are P-code; these facts were probed in MATLAB R2023b by driving QSM_GUI headless with logging wrappers on the functions it calls (3-echo synthetic head, TE 4/10/16 ms, 1×1×1.5 mm):

• Chain (the GUI’s “3D GRE data” then “iLSQR” buttons): MRPhaseUnwrap on the raw, unmasked phase of every echo (default padsize, 12 voxels), the unwrapped echoes summed, V_SHARP(Σφ, mask, 'smvsize', 12) (12 is the GUI’s “Radius” field default), then QSM_iLSQR on V-SHARP’s eroded mask. See crate::unwrap::laplacian and crate::bgremove::vsharp for those stages.
• TE: the GUI passes TE = the sum of the echo times (its “sum_TE” field: 30 ms for 4 + 10 + 16), matching the summed phase. For iLSQR this only sets the unit of χ.
• Parameters the GUI passes (as a params struct): niter 40, padsize [12 12 12], tol_step1 0.01, tol_step2 0.001, Kthreshold 0.1. The struct call is bitwise equal to QSM_iLSQR(..., 'padsize', [12 12 12], 'niter', 40). Called without them, QSM_iLSQR uses niter 100 (bitwise equal to an explicit 100) and padsize 6 mm (see Padding). IlsqrParams::default takes the function’s 100 iterations (decided 2026-10-29: UK Biobank and STI_Templates.m l. 115 call QSM_iLSQR without 'niter', so 100 is what they run; 40 is only what the GUI passes) and the GUI’s 12 mm padding. (STI_Templates.m lists padsize = [12 12 12] as the default; the function’s own default is 6.) Until 2026-10-29 the default was the GUI’s 40; the QSM-CI parity pair against the verified ilsqr-sti container (which calls it without 'niter') differed by 17–18 % NRMSE at 40 and 1–2 % at 100 (astewartau/QSM.rs#329); that 1–2 % was the FastQSM closing (see Agreement with STI Suite).
• Masking: the GUI’s “Single button for QSM” path builds the mask by thresholding the magnitude (Support_Functions/CreateMask.m: rescale, threshold, largest connected component, fill holes; then the P-coded Getmasking/PolishMask), not with BET.
• B0 direction: the GUI’s H field defaults to [0 0 1] and is not filled from the data; the v3.0 manual says to compute it from the DICOM for oblique scans.
• Other callers: QSMART (QSMART.m l. 32–39), like Hongfu Sun’s qsm_r2s_prisma.m whose call it copies, passes the field in rad/s with 'niter', 50, 'TE', 1000 and no 'padsize'. QSMART used STI Suite v2.2 (paper §3.4); v2.2’s QSM_iLSQR agrees with v3.0’s to 8e-5 relative on two synthetic cases. UK Biobank passes 'padsize', [64 64 64].

References: Li, W., Wang, N., Yu, F., Han, H., Cao, W., Romero, R., Tantiwongkosi, B., Duong, T.Q., Liu, C. (2015). “A method for estimating and removing streaking artifacts in quantitative susceptibility mapping.” NeuroImage, 108:111-122. https://doi.org/10.1016/j.neuroimage.2014.12.043

Paige, C.C., Saunders, M.A. (1982). “LSQR: An algorithm for sparse linear equations and sparse least squares.” ACM Trans. Math. Softw. 8(1):43-71.

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): sti.unwrap-laplacian, sti.laplacian-solver, sti.echo-combination, sti.echo-weights, sti.bg-vsharp, sti.vsharp-no-deconvolve, sti.vsharp-radius, sti.vsharp-grid, sti.ilsqr-iterations, sti.ilsqr-pad.

Modules

timing 🔒

Wall-clock tracing of iLSQR for native benchmarks and parity checks. Not compiled for wasm32, which has no clock in std.

Structs

IlsqrParams

Parameters for the iLSQR algorithm.

IlsqrTrace

Where an ilsqr_traced run spent its time, and how its two LSQR solves ended.

LsqrOutcome

How an lsqr_matlab solve ended (MATLAB lsqr flags).

Enums

IlsqrPrecision

Arithmetic of iLSQR’s two LSQR solves (see the module docs, Performance and precision).

Constants

NONE 🔒

Sentinel of the compact-index maps in the artefact solve.

NORM_CHUNK 🔒

Chunk length of the LSQR norms: fixed, so every sum is independent of the thread count.

STI_DEFAULT_NITER

QSM_iLSQR’s own iteration cap, used when it is called without 'niter' (output bitwise equal to an explicit 100), as UK Biobank’s pipeline and STI Suite’s STI_Templates.m call it, and the default of IlsqrParams::max_iter.

STI_DEFAULT_PAD_MM

Zero padding STI’s QSM_iLSQR uses when called without 'padsize' (as QSMART calls it), in mm. Probed in STI Suite 3.0: output bitwise equal to an explicit padsize of 6 mm (and different at 5.5 and 6.5) on five crop / voxel-size combinations.

STI_GUI_NITER

Iteration cap (niter) STI Suite 3.0’s QSM_GUI passes to QSM_iLSQR. Not the default of IlsqrParams::max_iter (that is STI_DEFAULT_NITER); set it to reproduce the GUI.

STI_GUI_PAD_MM

Zero padding STI Suite 3.0’s QSM_GUI passes to QSM_iLSQR (padsize [12 12 12]), in mm, and the default of IlsqrParams::pad_mm. STI_Templates.m and the χ-separation toolbox’s Chisep_script.m use the same value.

STI_PAD_MM

Zero padding STI’s QSM_iLSQR is given by UK Biobank and in STI Suite’s examples, in mm.

Traits

Tracer 🔒

Instrumentation hooks for ilsqr_with. The library path (ilsqr) uses (), whose hooks do nothing, so it never reads a clock (std::time::Instant panics on wasm32-unknown-unknown). The hooks only observe; they cannot change the result.

Functions

artefact_solve 🔒

Step 6 of ilsqr: the streaking artefacts, by LSQR (tol) over the k-space cone |D| < 0.1 against the weighted gradient (weights gw, full grid, zero off the mask) of the initial solution x1. Returns the real part of F⁻¹ of the solution (the artefacts divided by √N) and how LSQR ended.

artefact_solve_single 🔒

Step 6 of ilsqr as STI computes it (IlsqrPrecision::Single): the right-hand side G ∇(M·Re χ₁) in double, rounded to single; the operator y ↦ G ∇(√N F⁻¹(M_ic y)) and its adjoint r ↦ M_ic F(∇ᴴ(G r)) / √N in single precision, on dense vectors (the unknown over the whole k-space grid, the residual over three whole image grids); MATLAB’s lsqr on single data. Returns Re F⁻¹(M_ic y) in double (the artefacts divided by √N), as artefact_solve does.

axpy 🔒

circ_next 🔒

Circular neighbour l + e_axis on the grid n.

cnorm32 🔒

‖x‖₂ of a single-precision vector, accumulated in f64 over dense_parts and rounded to f32 (the value MATLAB’s norm of a single vector rounds to, up to its BLAS’s order).

cnorm_parts 🔒

compact_parts 🔒

See dense_parts. pos are the (strictly increasing) positions in the long vector.

dense_parts 🔒

How an LSQR vector is split for its norms: lengths of consecutive parts, each summed sequentially, the partial sums then added in order.

div_inplace 🔒

grad_masked 🔒

Weighted circular forward differences (no voxel-size scaling), as STI’s iLSQR uses, at the mask voxels only: r[a·nm + i] = gw[a][i] · (val(l + e_a) − val(l)) with l = inside[i], nm = inside.len(), gw the weights compacted to the mask.

grad_masked_adj 🔒

Adjoint of grad_masked onto the box bx of the grid: acc[l] = Σ_a (q_a[l − e_a] − q_a[l]) for l in bx (the rest of acc is left as it is), summed in axis order from zero, with q_a = gw[a]·r_a on the mask and zero elsewhere (mask_of[l] = position of l in inside, or NONE). The result is zero outside the mask grown by one voxel along +e_a.

ilsqr

iLSQR as in STI Suite 3.0 (QSM_iLSQR), padded by params.pad_mm.

ilsqr_traced

super::ilsqr, also returning an IlsqrTrace (wall time per stage and per FFT, and how both LSQR solves ended). Native targets only; the result is identical to ilsqr’s.

ilsqr_with 🔒

ilsqr with instrumentation hooks.

initial_solve 🔒

Step 3 of ilsqr: LSQR (tol 0.01) on D W D χ = D W φ with D = F⁻¹ diag(dk) F. Returns χ and how LSQR ended.

initial_solve_single 🔒

Step 3 of ilsqr as STI computes it (IlsqrPrecision::Single): the right-hand side D W φ in double, rounded to single; the operator D W D evaluated in single precision (single 3-D FFTs, single(D), single(W)) and MATLAB’s lsqr on single data. Returns the real part of the solution in double, as STI’s double(x1).

lsqr

LSQR iterative solver for Ax = b

lsqr_matlab

LSQR (Paige & Saunders 1982) with MATLAB lsqr’s stopping rules, zero initial guess, no preconditioner: stops when ‖r‖ ≤ tol·‖b‖ or ‖Aᴴr‖ / (‖A‖_F·‖r‖) ≤ tol (both tested on the previous iterate, before it is updated), after three stagnant steps, or on breakdown. apply_a computes A x, apply_ah computes Aᴴ u.

lsqr_matlab_f32 🔒

MATLAB lsqr (zero initial guess, no preconditioner) on single data: the iteration and stopping rules of lsqr_matlab, with every vector and scalar in single precision, as MATLAB runs it when b and the operator’s results are single (a double scalar such as tol or eps meeting a single one is rounded to single first).

lsqr_matlab_into 🔒

lsqr_matlab with operators that write into a caller-owned buffer (out holds stale data on entry and must be overwritten), taking b by value as LSQR’s first u, so a solve allocates its six vectors once. parts = (range parts, domain parts): how the norms of u and v vectors are split (see dense_parts); the domain parts also give the length of the unknown. Performs exactly the same floating-point operations as the textbook loop (the norms of updated vectors are fused into the update passes, with the same parts).

norm 🔒

nsq 🔒

|z|² (exactly z.norm_sqr())

scale_inplace 🔒

sign0 🔒

split_parts 🔒

split_parts_mut 🔒

sti_bounding_box 🔒

STI’s MaskBoundingBox: tight bounding box (0-based, inclusive); an odd extent grows by one voxel at the high end, or at the low end when the high end is the array edge (and shrinks by one at the high end when both ends are edges).

sti_dipole_kernel 🔒

STI’s dipole kernel D = 1/3 − (k·H)² / |k|² (0 at DC) on the unshifted padded grid, in STI’s arithmetic: k from sti_freqs, k·H = k_x H_x + k_y H_y + k_z H_z and |k|² = k_x² + k_y² + k_z² in that order, and H used as given when it is a unit vector (normalised otherwise). See sti_freqs for why the rounding matters.

sti_fastqsm 🔒

STI’s FastQSM on the padded grid, given F φ (unshifted k-space). Returns the masked estimate, rescaled to the TKD (threshold 1/8) solution by least squares inside the mask.

sti_freqs 🔒

STI’s k-space coordinates (cycles/mm) on the unshifted layout of an even axis of n voxels: the signed index divided by n·vs, as STI computes (-n/2:n/2-1)/(n*vs).

sti_padded_dims 🔒

STI’s padded grid: per axis the next multiple of 16 strictly above c + 2·pad_mm/vs.

sti_percentile_index 🔒

0-based index of STI’s percentile: the sorted value at 1-based index round(p·N).

sti_percentiles 🔒

STI’s percentiles p_lo ≤ p_hi of v (which is reordered), by selection instead of a full sort: the k-th element in total_cmp order is unique, so this is exactly the sorted lookup.

transform 🔒

One of [Fft3dWorkspace]’s pruned / fused transforms (forward, or inverse if inverse).

update_and_norm 🔒

f applied elementwise to y (paired with x), then cnorm_parts(y), in one pass: the same values as the two passes, part by part.

update_d_x_and_norms 🔒

LSQR’s direction and solution updates d ← (v − θ d)/ρ, x ← x + φ d in one pass, returning (‖d‖, ‖x‖): the same values as two update_and_norm passes, part by part.