Module chi_sep_ilsqr

Expand description

χ-separation (iLSQR): the SNU χ-separation toolbox’s chi_sep_iLSQR

chi_sep_ilsqr implements chi_sep_iLSQR of the χ-separation toolbox v1.1.3 (SNU-LIST, Jongho Lee’s lab; functions/corefunc/chi_sep_iLSQR.p and its chi_sep_iLSQR_step1_func.p, step2_func.p, step3_func.p), as its script Chisep_script.m calls it (l. 167-172): chi_sep_iLSQR(mag, local_field_hz, r2prime, mask, params, option_data) with params.b0_dir, params.CF, params.voxel_size, option_data.qsm and option_data.N_std. The toolbox ships as protected P-code; with the authors’ permission the algorithm was recovered by black-box probing in licensed MATLAB R2023b, as for [super::chi_sep_medi]: the profiler (which lists the bodies of the P-code’s anonymous functions verbatim), logging wrappers on the functions it calls that save the caller’s whole workspace (evalin('caller', 'save ...')), replaying the function handles captured in those workspaces, and step-by-step comparison of each step function’s inputs and outputs. A MATLAB re-implementation built from that reproduces chi_sep_iLSQR to rounding with the same iteration counts (see Agreement with the toolbox).

Earlier QSM.rs versions ported Shin et al. 2021’s Supplementary Methods instead: an edge-masked-TV Gauss-Newton solve with an SNR weight from the magnitude. The toolbox runs a different algorithm, an iLSQR-style three-step solve (Li et al. 2015) on the two sources; the old port’s output differed from the toolbox’s by relative L2 0.51 (χ+) / 0.54 (|χ−|) in vivo and 1.06 / 1.30 on QSM-CI’s chisep-mc phantom (the cases below).

§The algorithm

The unknowns are x+ = χ+ and x− = |χ−| in ppm. Everything is converted to ppm of the central frequency first: the local field f, R2′ as r = max(R2′/CF·10⁶, 0) (negative R2′ is clamped, after the conversion) and d± = D_r,±/CF·10⁶. With Dc(x) = real(ifftn(D.*fftn(x))):

Defaults (chi_sep_iLSQR): B0 = CF/42.58·10⁶; D_r,pos = D_r,neg = 137·(B0/2.9) (toolbox_default_dr); pad_size [16 16 16]; cg_max_iter_LSQR 50, cg_tol_LSQR 0.02, cg_max_iter_SA 50, cg_tol_SA 0.02, D2_thr 0.11, prc_W_I_min/max 60/99.8, WG_prc_min/max 50/70 (see ChiSepIlsqrParams); not parameters: the NNLS tolerance 0.01 and iter_max 10 (NNLS_TOL, NNLS_ITER_MAX).
Padding: every volume is zero-padded by pad_size voxels at both ends of every axis; the steps below run on that grid and crop their output back.
Dipole kernel: MEDI’s dipole_kernel(N, vs, B0_dir, 'kspace'), the same as chi_sep_MEDI’s, including its half-bin offset on odd axes.
Step 1, LSQR (chi_sep_iLSQR_step1_func). Weights: W_I = clamp((hi − L)/(hi − lo), 0, 1) with L = |del2(f, vs(1), vs(2), vs(3))| and lo, hi its prctiles 60 and 99.8 inside the mask. (del2 takes its first spacing argument for the second array axis, so on a grid with vs(1) ≠ vs(2) the first two spacings are swapped; kept.) w_mag = W_I·m, w_phs = 2π·W_I·m (m the mask). With Ar(x, d) = Dc(w_mag.*(d.*x)), ArH(y, d) = d.*(w_mag.*Dc(y)) and Af = AfH = Dc(w_phs.*Dc(·)), MATLAB’s pcg solves [ArH(Ar(x+,d+) + Ar(x−,d−), d+) + AfH(Af(x+ − x−)); ArH(…, d−) − AfH(Af(x+ − x−))] = [ArH(Dc(w_mag.*r), d+) + AfH(Dc(w_phs.*f)); ArH(Dc(w_mag.*r), d−) − AfH(Dc(w_phs.*f))] for the concatenated [x+ x−] (tolerance cg_tol_LSQR, at most cg_max_iter_LSQR iterations). The start is the voxelwise solution from the QSM χ: x+ = max((r + d−·χ)/(d+ + d−), 0), x− = max((r − d+·χ)/(d+ + d−), 0), masked. (Without a QSM the toolbox uses χ = 0.)
Step 2, streaking artefacts (chi_sep_iLSQR_step2_func), per source:
- a fast QSM χ_fs of the field: W = clamp((|D|^0.001 − a)/(b − a), 0, 1) with a, b the prctiles 1 and 30 of |D|^0.001 over the whole grid; K the k-space smoothing by a 2.5 mm partial-volume ball (pv_ball below); k1 = ifftshift(fftn(f).*sign(D)), then χ2 = real(ifftn(fftshift(W.*k1 + (1 − W).*K(k1)))). That multiplies the weight, laid out for an unshifted spectrum, with the shifted spectrum (a quirk of the toolbox, kept). Then k = fftn(χ2.*m), χ3 = real(ifftn(W.*k + (1 − W).*K(k))).*m, and χ_fs = c₁ + c₂·χ3 with (c₁, c₂) the least-squares fit of χ3 to the TKD solution (threshold 1/8, masked) over the whole grid;
- R2′ smoothed in k-space by a Tukey window (tukeywin(N, 0.5) per axis, applied as fftshift(ifftshift(fftn(r)).*w));
- gradient weights W_G± = clamp((hi − |g±|)/(hi − lo), 0, 1)·m per axis, lo/hi the prctiles 50/70 of |g±| inside the mask, g+ = (∇r_s + d−·∇χ_fs)/(d+ + d−), g− = (∇r_s − d+·∇χ_fs)/(d+ + d−) (∇: MEDI’s fgrad with the voxel size, on the masked volumes);
- the streaking component, by MATLAB’s pcg on a complex k-space unknown y restricted to the cone M_ic = |D| < D2_thr: M_ic.*fftn(∇ᵀ(W_G²·∇ real(ifftn(M_ic.*y)))) = M_ic.*fftn(∇ᵀ(W_G²·∇x)) (∇ here with unit spacing, ∇ᵀ its exact adjoint; tolerance cg_tol_SA, at most cg_max_iter_SA iterations, zero start), and the output is (x − real(ifftn(M_ic.*y)).*m).*m.
Step 3, NNLS (chi_sep_iLSQR_step3_func): repeat step 1 with one pcg iteration starting from the current iterate (not clamped), then step 2, then project both sources onto x ≥ 0. res = (‖Δx+‖ + ‖Δx−‖)/(‖x+‖ + ‖x−‖) of the projected iterates; the loop runs while res > 0.01 and at most iter_max + 1 = 11 times (its counter starts at 1 and is tested with ≤ iter_max + 1). The output is the iterate with the smallest res.
Output: x_para = x+, x_dia = x− (≥ 0), x_tot = x+ − x−, zero outside the mask. QSM.rs returns (χ+, χ− = −x−, χ+ + χ−).

MATLAB pcg. The solver is R2023b’s pcg.m without a preconditioner, line for line: the residual recomputed when the recurrence says converged or stagnates (up to min(n/50, 5, n − maxit) extra steps), and, when it does not converge, the iterate with the smallest residual returned instead of the last one (so a one-iteration solve that does not lower the residual returns its start). For the complex unknown of step 2 it uses the conjugate inner product and compares p'·q ≤ 0 by its real part, as MATLAB does.

§Inputs

Chisep_script.m passes the magnitude, the local field in Hz, R2′, the mask and option_data.qsm (STI Suite QSM_iLSQR) and option_data.N_std. chi_sep_iLSQR never uses mag or N_std: neither reaches its step functions, and the re-implementation that ignores them matches the toolbox (and its output is unchanged when they change, see the fixture generator). ChiSepIlsqrInputs therefore takes the local field, R2′, the QSM and the mask. The QSM is required, as the script passes it. The toolbox’s GUI (WrapperChisepiLSQR.p, probed with a stand-in app object) passes no QSM (option_data holds only N_std and zero mask_Unrel_Mag/mask_Unrel_Phs), so its runs start step 1 from zero; QSM.rs does not offer that path.

§Units

The toolbox takes the local field in Hz and converts it to ppm; chi_sep_ilsqr takes it in ppm directly. R2′ is in Hz, the QSM and the output in ppm. ChiSepIlsqrParams::cf sets the ppm scale of R2′ and D_r and the default D_r.

§Agreement with the toolbox

Compared with chi_sep_iLSQR itself (toolbox v1.1.3 in licensed MATLAB R2023b; QSM.rs in double precision, as the toolbox runs):

Phantoms (tests/chi_sep_ilsqr_snu_parity.rs, fixture from tests/data/chi_sep_ilsqr_snu_small.m): 40×36×28 at 3 T; 33×29×21 at 7 T, 0.8×0.9×1.5 mm, oblique B0; negative R2′ with D_r 120/150. Same NNLS iteration counts (2, 4, 4); χ+ and |χ−| to relative L2 ≤ 4·10⁻¹⁰ (3·10⁻¹⁴ on the second). The toolbox itself moves by up to 1.5·10⁻⁹ on these phantoms between 8 threads and 1.
In vivo (the EMTP travelling-head scan of QSM-CI’s fidelity check, 3 T, 184×256×128; its local field, R2′ and χ, 1.39 M mask voxels): 4 NNLS iterations each; χ+ relative L2 3.1·10⁻¹¹ (max |Δ| 2.5·10⁻¹⁰ ppm), |χ−| 3.5·10⁻¹¹ (3.2·10⁻¹⁰ ppm).
QSM-CI’s χ-separation phantom (chisep-mc, 7 T, 164×205×205, same inputs): 4 NNLS iterations each; χ+ relative L2 2.5·10⁻¹² (max |Δ| 8.8·10⁻¹² ppm), |χ−| 3.2·10⁻¹² (7.2·10⁻¹² ppm). The toolbox run with 4 instead of 8 threads moves by 1.1·10⁻¹² / 1.5·10⁻¹² here and 2.6·10⁻¹² / 2.9·10⁻¹² in vivo.
The verified QSM-CI chi-sep-ilsqr container (QSMxT/QSM-CI#306; its own STI Suite QSM, tests/chi_sep_ilsqr_qsmci_parity.rs, regenerated in Actions by tests/data/chi_sep_ilsqr_qsmci_small.py): χ+ and |χ−| to relative L2 3·10⁻⁸ on a 24×22×18 phantom, the container’s float32 output (measured 2026-10-30). On chisep-mc and the in-vivo scan above, given the container’s own QSM (its intermediate) and its CF (42.5774 MHz/T·B0): χ+ and |χ−| to relative L2 2.6·10⁻⁸ and 2.5·10⁻⁸, max |Δ| 7·10⁻⁹ and 3·10⁻⁸ ppm, its float32 rounding (examples/chi_sep_file.rs, measured 2026-11-01). QSM-CI’s chisep-ilsqr / chisep-ilsqr-qsmrs parity pair checks this on chisep-mc.

The residual is rounding inside the iterations. QSM.rs evaluates the same expressions with its own FFT and, where that saves work, in a different order (step 1’s operator applies Dc∘Dc once to the sum of the two sources’ R2′ terms instead of to each; the fast QSM’s fit uses centred sums instead of a QR solve), so every operation rounds differently, where the toolbox’s thread count changes only its reductions; in vivo that puts QSM.rs at about ten times the toolbox’s thread-count spread. The solve does not amplify noise in its input: a ±10⁻¹⁵ relative perturbation of the field moves QSM.rs’s in-vivo result by 4·10⁻¹⁵.

§Performance

Each step-1 pcg iteration costs six dipole convolutions (12 complex 3-D FFTs) on the padded grid; each step-2 iteration two FFTs; an NNLS iteration is one step-1 iteration (three operator applications) and two step-2 solves. The fixed parts (weights, fast QSM, smoothed R2′) are computed once, where the toolbox recomputes them in every NNLS iteration. In vivo (216×288×160 padded, 4 NNLS iterations) the toolbox takes 433 s on 8 threads and on chisep-mc (196×237×237) 1113 s (743 s on 4); QSM.rs took 280 s on 6 threads and 398 s with the memory layout before the one below, which is faster again (155 s → 91 s on 12 threads on a phantom of the in-vivo size).

Memory. The arithmetic is the toolbox’s, in the order above; what is stored is cut to what each step needs, without changing a single operation:

the dipole kernel is evaluated where it is used, from per-axis tables, by the same expressions as chi_sep_medi’s kernel (bit for bit, tested);
step 1’s operator chains its convolutions through one complex buffer (each convolution’s real part is the next one’s input, as before) and keeps its R2′ term in its output;
W_I and the NNLS iterates are stored at the mask’s voxels only (they are zero elsewhere); step 2’s weights W_G² are recomputed where used from the masked fast QSM and smoothed R2′ (kept on the mask’s bounding box) and their percentiles, by the same expressions; the gradient of step 2’s regulariser is fused into its adjoint;
step 2’s complex pcg stores only the cone’s voxels (19 % of k-space at the default D2_thr): its vectors are zero outside the cone, and its sums add the same terms in the same chunks;
pcg copies its best iterate (xmin) only once it differs from the current one, and the NNLS iterations’ start is rebuilt from the iterates rather than stored;
the k-space shifts run in place.

The peak is step 1’s first solve: pcg’s x, xmin, q, r, p (10 f64 per padded voxel), the right-hand side (2) and the FFT buffer (2), plus what the NNLS iterations keep (the mask, W_I, the two volumes W_G² is computed from, the cone): about 14.25 f64 per padded voxel + 1.5 per mask voxel + 2 per voxel of the mask’s bounding box. Later, with the current and best iterates held (4 per mask voxel), it is 12.25 + 5.5 per mask voxel + 2 per box voxel, which is more only if the mask fills over half the padded grid; a bound for any mask is 14.25 per padded voxel + 7.5 per input voxel. A larger D2_thr widens the cone and step 2’s share (it exceeds step 1’s above about 0.25). Measured on the phantom of examples/chi_sep_ilsqr_memory.rs at the in-vivo size (184×256×128, 16 % of the padded grid in the mask): a 1.20 GB heap peak, 15.1 f64 per padded voxel (1.34 GB resident), where the earlier layout needed 3.43 GB (43 per padded voxel). The results are the same numbers: χ+, χ− and χ_total equal the earlier layout’s bit for bit inside the mask (on the parity fixtures and that phantom, with and without threads); outside it both are zero, now +0 throughout where the earlier one left some −0.

§Deviations

None in the algorithm. Options of chi_sep_iLSQR that Chisep_script.m does not use are not exposed: option_data.mask_Unrel_Mag / mask_Unrel_Phs (masks of unreliable magnitude and phase voxels whose data weights drop; both default to zero), and params.flag_strong_x_correction / strong_x_thresh (off by default). The start without a QSM (the GUI’s) is not offered.

§Reference pipeline

What Chisep_script.m (toolbox v1.1.3) runs before chi_sep_iLSQR is what it runs before chi_sep_MEDI (see [super::chi_sep_medi]’s Reference pipeline): Tukey windowing, MEDI’s BET, ROMEO with weighted echo averaging (W = TE·exp(−TE/40 ms)), V-SHARP smvsize 12, the local field in Hz (l. 139), and the QSM from STI Suite QSM_iLSQR(..., 'padsize', [12 12 12]) (l. 145, chisep_toolbox_ilsqr) passed as option_data.qsm (l. 170). R2′ is the user’s. N_std (l. 76) is passed (l. 171) but not used.

Reference: Shin, H.-G., Lee, J., Yun, Y. H., Yoo, S. H., Jang, J., Oh, S.-H., Nam, Y., Jung, S., Kim, S., Fukunaga, M., Kim, W., Choi, H. J., Lee, J. (2021). “χ-separation: Magnetic susceptibility source separation toward iron and myelin mapping in the brain.” NeuroImage, 240:118371. https://doi.org/10.1016/j.neuroimage.2021.118371. The three steps follow Li, W., Wang, N., Yu, F., et al. (2015). “A method for estimating and removing streaking artifacts in quantitative susceptibility mapping.” NeuroImage, 108:111-122.

§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): chi-sep.mask-bet, chi-sep.unwrap-romeo, chi-sep.echo-combination, chi-sep.echo-weights, chi-sep.bg-vsharp, chi-sep-ilsqr.qsm-ilsqr, and by manufacturer chi-sep.tukey-siemens-ge and chi-sep.tukey-philips (Chisep_script.m Tukey-windows the echoes before this method as before χ-sepnet; [crate::utils::tukey_window]). All are the authors’ setting (Kind::Reference).

Structs§

BoxVol 🔒: A volume stored on a box of the padded grid, zero outside it.
ChiSepIlsqrInputs: The volumes chi_sep_ilsqr takes. chi_sep_iLSQR(mag, local_field_hz, r2prime, mask, params, option_data) also takes the magnitude and option_data.N_std, but never uses them (see the module docs). All are on the grid, column-major.
ChiSepIlsqrParams: Parameters for chi_sep_ilsqr. The defaults are the toolbox’s (see the module docs); the field names after the toolbox’s params fields are given with each.
Kernel 🔒: The toolbox’s dipole kernel (chi_sep_medi::toolbox_dipole_kernel) evaluated where it is used, from its per-axis k coordinates: the same expressions on the same values, so the same numbers, without a volume of it.
PcgOut 🔒: What pcg returns: the solution, MATLAB’s flag and iter.
Problem 🔒: The arrays of one run, on the padded grid. What is fixed across the NNLS iterations is kept in the least memory that gives the same numbers (see Performance in the module docs).
Support 🔒: A set of voxels of the padded grid, for vectors stored on those voxels only (the voxels of the full vector outside it are zero): their indices, ascending, and for each CHUNK of the full grid where its entries start. Sums over such a vector (Self::sum) add the same terms in the same chunks as det_sum over the full vector, less the zeros, so they are equal.
Workspace 🔒: FFT buffers reused across operator applications.

Enums§

Start 🔒: Where pcg starts.
XMin 🔒: Where pcg’s best iterate so far (MATLAB’s xmin) is: the start (rewritten by Start::Fill when needed), the current iterate, or its own copy.

Constants§

BALL_RADIUS_MM 🔒: Radius, in mm, of the ball the fast QSM smooths k-space with.
CHUNK 🔒: Chunk length of det_sum’s partial sums.
NNLS_ITER_MAX: iter_max of chi_sep_iLSQR_step3_func (not a parameter). The NNLS loop runs at most NNLS_ITER_MAX + 1 = 11 times: its counter starts at 1 and is tested with ≤ iter_max + 1.
NNLS_TOL: NNLS stopping tolerance of chi_sep_iLSQR_step3_func (its tol_norm_ratio, not a parameter).
TKD_THRESHOLD 🔒: TKD threshold of the fast QSM’s scaling reference.
TUKEY_RATIO 🔒: Tukey window ratio of the R2′ smoothing.

Traits§

PcgElem 🔒: An element of a pcg vector: real (step 1) or complex (step 2).

Functions§

chi_sep_ilsqr: χ-separation by the SNU χ-separation toolbox’s chi_sep_iLSQR (see the module docs).
chisep_toolbox_ilsqr: The iLSQR call of the χ-separation toolbox (Chisep_script.m l. 142-145) for the QSM that chi_sep_ilsqr and [super::chi_sep_medi] start from: STI_DEFAULT_NITER iterations (QSM_iLSQR’s own default, as the toolbox calls it without 'niter'), STI_GUI_PAD_MM padding, STI’s tolerances.
chunked_sum 🔒: f summed over each of chunks ranges in order, then the partial sums in order
circshift 🔒: fftshift (forward, shift by floor(n/2)) or ifftshift (shift by ceil(n/2)) per axis, in place: a rotation of each x-row, of each z-plane by whole rows, and of the volume by whole planes
clamp01 🔒: MATLAB’s min(max(v, 0), 1), which maps NaN to 0 (clamp would keep it)
del2 🔒: MATLAB’s del2(f, hx, hy, hz) with the spacings already in array-axis order (h[0] along the first axis): per axis, half the centred second difference inside and linear extrapolation at the ends, summed over the axes and divided by 3. Coordinates are h·(1:n) and their differences are taken as MATLAB takes them.
det_sum 🔒: A sum over 0..n in fixed chunks, so the result does not depend on the thread count
dp_mul 🔒: d * v, named so the operation order of the toolbox’s d.*(w.*·) stays visible
par_chunked 🔒: f(c, chunk) for each CHUNK of a (over threads with the parallel feature)
par_chunked2 🔒: par_chunked over two volumes at once
par_rows 🔒: f(j, k, l0, row) for each x-row of a volume (l0 its first index; over threads with the parallel feature)
pcg 🔒: MATLAB R2023b’s pcg(A, b, tol, maxit, [], [], x0) (no preconditioner), line for line.
prc_on 🔒: (lo, hi), the prctiles p of v(l) over the mask’s voxels l, for the weights min(max((hi − v)/(hi − lo), 0), 1)
prctile 🔒: MATLAB’s prctile(x, p) (R2023b, 'exact', no NaNs): with r = (p/100)·n and k = floor(r + 0.5), (0.5 + r − k)·x(k+1) + (0.5 − r + k)·x(k) of the sorted values, clamped to the first and last and copied where the two neighbours are equal. x is reordered.
psum 🔒: det_sum of f over a vector of n, or over its voxels in on (the rest zero)
pv_ball 🔒: The fast QSM’s smoothing ball: radius r mm on the grid (−n/2 : n/2−1)·vs per axis (half integers on an odd axis, as the dipole kernel), 1 where the voxel centre is within r − h (h the half diagonal), and partial-volume weighted by 20³ sub-samples at offsets (−0.475:0.05:0.475)·vs where it is within r + h; normalised to sum 1. Centred.
sign0 🔒: MATLAB’s sign (0 at 0)
tukeywin 🔒: MATLAB’s tukeywin(n, ratio) for 0 < ratio < 1