Module star 

STAR-QSM: streaking artifact reduction for QSM

star_qsm implements STI Suite 3.0’s QSM_star (Wei et al. 2015), as called with QSM_star(phase, mask, 'TE', TE, 'B0', B0, 'H', H, 'padsize', pad, 'voxelsize', vs) with pad = StarParams::pad. Its default, STAR_GUI_PAD = [8, 8, 20] voxels, is what STI Suite’s QSM_GUI passes (see StarParams::pad); QSM_star called without 'padsize' pads QSM_STAR_DEFAULT_PAD = [0, 0, 0]. STI Suite ships as protected P-code, so the algorithm was recovered by black-box probing (call tracing through logging wrappers around the functions it calls, capturing the arrays it hands them, then a MATLAB re-implementation that reproduces QSM_star to within 4e-16 relative on phantoms and an in-vivo field):

1. Zero-pad the phase φ (radians at the echo time) and the mask M by pad[a] voxels on both sides of axis a. STI’s dipole kernel errors on odd extents, so QSM_star does too; here an odd padded extent grows by one more voxel at the high end instead. Everything below runs on that grid.
2. Dipole kernel D = 1/3 − (k·H)²/|k|² (D = 0 at k = 0) with k in cycles/mm (computed as STI does, to the bit: see sti_kernel), and the k-space mask K = |D| ≥ 0.065. STI builds both on a centred k-grid and transforms with centred FFTs (fftshift ∘ fftn ∘ fftshift, unitary); on an even grid every operator below is the same as the plain FFT convolution used here.
3. Initial guess: z = M·F⁻¹(D⁺ F(Mφ)) with D⁺ = 1/D on K and ±1/0.065 off it (the sign of D, + where D = 0); then twice z ← M·F⁻¹(K·F z); then x₀ = z / max|z|. All of it stays complex.
4. Solve min ½‖A x − Mφ‖² + τ‖x‖₁ (A x = F⁻¹(D F(M x)), τ = 10⁻⁶, complex x) from x₀ with SpaRSA 2.0 (Wright, Nowak & Figueiredo) at the options STI passes: 'AT' = A itself (not its adjoint M F⁻¹ D F), stopping criterion 0 (stop once fewer than 1 % of the components change between zero and non-zero, after at least 5 and at most 1000 iterations), and SpaRSA’s defaults otherwise: Barzilai–Borwein step α = ‖A Δx‖²/‖Δx‖² clamped to [10⁻³⁰, 10³⁰], α = 1 at the start, no monotone or safeguard test (every step is accepted), no continuation, no debiasing, and SpaRSA’s complex soft threshold x ↦ max(|x| − τ/α, 0)/(max(|x| − τ/α, 0) + τ/α) · x. STI also passes a TV operator as 'PHI'; SpaRSA ignores a 'PHI' given without 'PSI'. If τ ≥ max|A φ| SpaRSA returns zero.
5. χ = 1.1 · Re(x) / (2π·γ·B0·TE) with γ = 42.575 MHz/T (STI’s ScalingFactor), cropped back to the input grid and masked. The factor 1.1 is a constant in STI’s code.

SpaRSA takes 5–12 iterations on the phantoms below and 5 on the in-vivo field.

Units

Like the other inversions, star_qsm takes the field in ppm and returns χ in ppm, but unlike them the result depends on the field’s scale: x₀ is normalised to a peak of one (in radians), τ is fixed and the solve stops after a handful of iterations, so the same field in ppm gives a different χ at another TE or B0. StarParams::te and StarParams::b0 are therefore part of the method: the field is converted to the phase φ = 2π·42.575·B0·TE · field it would have at the echo time te picks from StarParams::echo_times, which is what STI is given. By default that is the sum of the echo times: STI Suite’s QSM_GUI gives QSM_star the V-SHARP of the echoes’ summed unwrapped phases with 'TE' their summed echo times (see StarParams::te). With that conversion star_qsm(field) is QSM_star(φ) exactly (to rounding), 1.1 included.

Agreement with STI Suite

On six 48×56×40 phantoms (3 T and 7 T, TE 5–20 ms, axial and oblique B0, voxels 0.8×0.9×1.5 mm, pad 6, 8 and 12 on every axis) and an in-vivo 7 T local field (164×206×206), χ agrees with STI Suite’s output to rounding: relative L2 difference ≤ 1.1e-14, max |Δχ| ≤ 2.6e-15 ppm, with the same iteration counts. Through STI Suite’s GUI (QSM_GUI: load, “3D GRE data”, “star QSM”, with QSM_star’s inputs and output captured) on 64×64×48 at 1×1×1.5 mm, 80×72×40 at 0.8×0.9×2 mm and 64×64×64 at 1 mm, the default STAR_GUI_PAD reproduces the GUI’s χ to rounding as well. tests/star_sti_parity.rs checks three of the phantoms and the GUI padding on the anisotropic one (licensed MATLAB R2023b, tests/data/star_sti_small.m). It also checks three 32×36×28 phantoms (3 T axial, anisotropic voxels, oblique 7 T with pad 8) against the verified QSM-CI star-sti container (QSMxT/QSM-CI#304), which writes χ as float32: QSM.rs is within half a float32 ulp of every voxel (relative L2 2.4–3.0e-8, the float32 rounding; measured 2026-10-30). That fixture is regenerated in Actions by tests/data/star_sti_qsmci_small.py.

Performance

Each SpaRSA iteration costs four 3-D FFTs on the padded grid (A twice); the forward transforms skip the zero rows outside the mask’s bounding box and multiply by D in their last pass. The initial guess and setup cost ten more. Memory is three complex and one real array on the padded grid. On the in-vivo field above (188×230×230 padded, 5 iterations) it takes 1.2–1.7 s on 16 cores; iLSQR takes 9 s on the same field and cores.

Reference pipeline

What STI Suite 3.0’s QSM_GUI runs around QSM_star, probed headless in MATLAB R2023b with logging wrappers (3-echo synthetic head, TE 4/10/16 ms, 1×1×1.5 mm). The “3D GRE data”

• “star QSM” buttons and the “Single button for QSM” path call the same functions:

• Chain: MRPhaseUnwrap on the raw, unmasked phase of every echo (default padsize, 12 voxels); the unwrapped echoes summed (V-SHARP’s input is bitwise sum(φ, 4)); V_SHARP(Σφ, mask, 'smvsize', 12) (no deconvolution; see crate::bgremove::vsharp); then QSM_star on V-SHARP’s eroded mask.
• TE: QSM_star gets 'TE' = the sum of the echo times (the GUI’s “sum_TE” field: 30 ms for 4 + 10 + 16; the GUI’s “Frequency” output is TissuePhase/(2π·ΣTE)), so the phase it inverts is the summed phase. STAR’s result depends on the phase scale (see Units), so this matters; StarParams::te defaults to EchoTime::Sum for this.
• Padding: [8, 8, 20] from the GUI (STAR_GUI_PAD). Called without 'padsize', QSM_star pads 0: its output is bitwise equal to 'padsize', [0 0 0], and differs from [12 12 12] by 8–15 % relative, on three grids (64×64×48 at 1×1×1.5 mm and at 1 mm, 54×60×40 at 0.8×0.9×2 mm). [12, 12, 12] (STI_TEMPLATES_STAR_PAD) is the default that STI_Templates.m lists, not what the function does.
• Masking: the GUI thresholds the magnitude (Support_Functions/CreateMask.m: rescale, threshold, largest connected component, fill holes; then the P-coded Getmasking/PolishMask), not BET.
• Data: the v3.0 manual’s test data are 3T multi-echo human brain GRE (0.9×0.9×2 mm, 3 echoes), single-echo human and mouse brain, and 2D GRE-EPI.

Reference: Wei, H., Dibb, R., Zhou, Y., Sun, Y., Xu, J., Wang, N., Liu, C. (2015). “Streaking artifact reduction for quantitative susceptibility mapping of sources with large dynamic range.” NMR in Biomedicine, 28(10):1294-1303. https://doi.org/10.1002/nbm.3383

Wright, S.J., Nowak, R.D., Figueiredo, M.A.T. (2009). “Sparse reconstruction by separable approximation.” IEEE Trans. Signal Process. 57(7):2479-2493. https://doi.org/10.1109/TSP.2009.2016892

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, star.padding, star.te-sum.

Structs

StarParams

Parameters for star_qsm. The defaults are STI Suite’s: the GUI’s padding, and QSM_star’s defaults otherwise.

Constants

ALPHA_MAX 🔒

ALPHA_MIN 🔒

SpaRSA’s default AlphaMin / AlphaMax.

CHUNK 🔒

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

KSPACE_THRESHOLD 🔒

|D| below which the initial guess clamps 1/D, and outside which K keeps k-space.

MAX_ITER 🔒

SpaRSA’s MaxiterA as STI passes it.

MIN_ITER 🔒

SpaRSA’s default MiniterA.

QSM_STAR_DEFAULT_PAD

Zero padding QSM_star uses when called without 'padsize', in voxels per side: none. Probed in MATLAB R2023b on STI Suite 3.0: QSM_star without 'padsize' is bitwise equal to 'padsize', [0 0 0] and differs from [12 12 12] by 8-15 % relative, on 64×64×48 at 1×1×1.5 mm and at 1 mm and on 54×60×40 at 0.8×0.9×2 mm.

STAR_GUI_PAD

Zero padding STI Suite 3.0’s QSM_GUI gives QSM_star, in voxels per side on each axis, and the default of StarParams::pad.

STI_CHI_FACTOR 🔒

The constant STI multiplies the solution by.

STI_GAMMA_MHZ_PER_T 🔒

γ/2π in MHz/T as STI Suite’s ScalingFactor uses it: the reference implementation’s value, not crate::constants. Keep it.

STI_TEMPLATES_STAR_PAD

Zero padding STI_Templates.m lists as QSM_star’s default (padsize = [12 12 12]), in voxels per side. Not what the function does without 'padsize' (QSM_STAR_DEFAULT_PAD).

TAU 🔒

SpaRSA’s L1 weight τ as STI passes it.

TOL_A 🔒

SpaRSA’s default ToleranceA (stopping criterion 0: fraction of components changing).

Functions

dinv 🔒

The initial guess’s inverse kernel: 1/D where |D| ≥ 0.065, else ±1/0.065 by the sign of D (+ at D = 0).

initial_guess 🔒

Step 3: x₀ = z / max|z| with z = M F⁻¹(K F(M F⁻¹(K F(M F⁻¹(D⁺ F y)))))). y is zero outside inside, which lies in support. None if z is zero.

mask_inplace 🔒

nsq 🔒

padded_dims 🔒

Padded grid: pad[a] voxels on both sides of axis a, plus one at the high end where that is odd.

soft 🔒

SpaRSA’s complex soft threshold (soft.m): max(|v| − t, 0)/(max(|v| − t, 0) + t) · v.

sparsa 🔒

SpaRSA 2.0 at STI’s options (module docs, step 4) for A x = op(M x), Aᵀ = A, with y zero outside inside. x holds the initial guess on entry and the solution on return. apply_core applies op in place. Returns the number of iterations (0 if τ ≥ max|A y|, when x is set to zero as SpaRSA does).

star_qsm

STAR-QSM as in STI Suite 3.0 (QSM_star).

sti_kernel 🔒

STI’s dipole kernel (calcD2Matrix) in unshifted order: D = 1/3 − (k·H)²/|k|², 0 at k = 0, with k_d = f/(n_d·vs_d) for the signed frequency index f (−n/2 at Nyquist) and H normalised.