Module decompose

Expand description

DECOMPOSE-QSM: signal-domain paramagnetic/diamagnetic source separation.

DECOMPOSE (Chen et al., NeuroImage 2021) separates paramagnetic (χ+, iron) and diamagnetic (χ−, myelin·calcium) susceptibility by fitting a three-compartment complex multi-echo gradient-echo signal per voxel. Each voxel’s signal is

  S(t) = C+ · exp(−( a·χ+  + R2*₀ + i·(2/3)·χ+·γ·B0)·t)
       + C− · exp(−(−a·χ−  + R2*₀ + i·(2/3)·χ−·γ·B0)·t)
       + C₀ · exp(−R2*₀·t)

with γ = 42.58·2π, a = 2π·γ·B0 / (9√3) the static-dephasing broadening coefficient, χ+ ≥ 0, χ− ≤ 0 (ppm-scale). C+, C−, C₀ are the paramagnetic, diamagnetic and neutral compartment amplitudes and R2*₀ the baseline decay.

As in the authors’ DECOMPOSE_analysis.m, the phase is synthesized from a QSM of each echo (theirs: STAR-QSM per echo, QSM_multi): the complex data per echo is y(TE) = |mag_norm(TE)| · exp(−i·(2/3)·χ(TE)·γ·B0·TE) (magnitude normalised by its global maximum). decompose takes that per-echo QSM; producing it is the caller’s. A voxel whose Σ y is zero is skipped, as theirs is.

The per-voxel fit is a 3-stage alternating bounded least-squares (lsqcurvefit), repeated n_inner times: (1) amplitudes [C+,C−,C₀] against the linear signal, then (2) R2*₀ and (3) [χ+,χ−] against log(signal) (complex log). The residual is the complex difference packed as [Re; −Im], matching the reference.

Each source is reconstructed from the fitted parameters via a per-compartment phase accumulation, −Σ angle(model) / ((2/3)·γ·B0·ΣTE): the paramagnetic sub-model (pscModel) gives χ+ and the diamagnetic sub-model (dscModel) gives |χ−|. χ+ is returned ≥ 0, χ− ≤ 0.

Phase unwrapping (off by default: DecomposeParams::unwrap_phase). The authors’ code takes every phase on its principal branch (MATLAB’s log and angle), and so does this port by default. With the option on: every phase the fit and the reconstruction see arrives through atan2, on (−π, π]. A compartment accumulates (2/3)·χ·γ·B0·TE radians, so above |χ_total| = π / ((2/3)·γ·B0·TE_max) the late echoes fold back onto that branch — 0.09 ppm at 7 T with a 28 ms last echo, which is ordinary globus pallidus. Left folded, the alternation converges to a degenerate fixed point (χ− pinned at the χ bound, R2*₀ and the amplitudes at zero) and returns χ+ ≈ 0, so the most paramagnetic structure in the volume reads as a hole. Both the stage-2/3 residuals and recon_phase therefore walk the echoes in order and keep each phase on the branch nearest its predecessor (unwrap_near); the synthesized data phase is known in closed form and is used unwrapped rather than round-tripped through atan2. The option does only this; it does not change the χ search (that is DecomposeParams::chi_bound). The original method, published at 3 T, does not meet this limit in practice — the threshold there is ~0.2 ppm.

The χ bound (DecomposeParams::chi_bound, a tuning value). The authors search 0 ≤ χ+ ≤ 0.5, −0.5 ≤ χ− ≤ 0 ppm, a bound “roughly calculated using TE = 25 ms and B0 = 3T to ensure the phase does not exceed 2π” (Chen et al. 2021, §2.2); that is the default (DecomposeChiBound::Ppm(0.5)). DecomposeChiBound::EchoSpacing instead takes the bound from the scan: π / ((2/3)·γ·B0·ΔTE_max), with ΔTE_max the largest phase step the unwrapping walks (TE₁ from zero, or any echo spacing). Unwrapping along the echoes is unambiguous only while the phase step between echoes stays under π, and this is the largest |χ| for which it does: 0.31 ppm (0.3145) at 7 T with 8 ms spacing, 0.73 ppm at 3 T. Past it the model matches a given χ on more than one branch and the stage-3 landscape grows spurious minima.

At 7 T. The authors designed the fit for 3 T, and at 7 T the phase wraps between echoes inside ordinary deep grey matter. Measured 2026-10-30 on the QSM-CI 7 T chisep phantom (164×205×205, 1 mm, TE 4/12/20/28 ms; its provided χ given at every echo) with this code at astewartau/QSM.rs@add3d6e (the lsqcurvefit port, FDLIBM), n_inner 10:

QSMCI_CHISEP=<phantom> DECOMPOSE_UNWRAP=<0|1> DECOMPOSE_CHI_BOUND=<ppm|echo-spacing> \
  cargo test --release --features parallel --test chisep_relaxometry_qsmci \
  test_decompose_qsmci -- --ignored --nocapture

“Pallidum” is the phantom’s globus pallidus (dseg label 2, 3683 voxels, true mean χ+ 0.137 ppm); “collapsed” counts its voxels whose χ+ is under a quarter of the truth.

`unwrap_phase`	`chi_bound`	χ+ r / XSIM / NRMSE	χ− r / XSIM / NRMSE	pallidum mean χ+	collapsed
off (authors)	0.5 (authors)	0.484 / 0.216 / 10.9 %	0.745 / 0.325 / 4.8 %	0.002	3541
on	0.5	0.658 / 0.250 / 10.0 %	0.751 / 0.324 / 4.7 %	0.027	2932
off	echo spacing (0.3145)	0.587 / 0.367 / 9.7 %	0.846 / 0.471 / 4.1 %	0.005	3372
on	echo spacing (0.3145)	0.895 / 0.489 / 7.1 %	0.865 / 0.477 / 4.0 %	0.111	17

So the two help in different ways. The bound alone improves χ− (r 0.75 → 0.85) and χ+‘s XSIM, but leaves the pallidum collapsed; unwrapping alone recovers little of it; together they recover it. (The rows were first measured with the bound hard-wired into the unwrapping option, as min(0.5, π / ((2/3)·γ·B0·ΔTE_max)); that option, and a plain chi_bound of 0.3145 with and without unwrapping, gave the same metrics to the digits shown, and so does DecomposeChiBound::EchoSpacing. An older note here, from the bounded-LM port, said the bound removed “a speckle of dropouts” that unwrapping alone left; with this solver, unwrapping alone leaves most of the pallidum collapsed, and voxels at exactly χ+ = 0 number 1–8 of 1.31 M in every row.) The defaults stay the authors’ (off, 0.5); the pipeline rules decompose.unwrap-7t and decompose.chi-bound-7t suggest both at 7 T (from 6.5 T, so a Siemens 7 T’s 6.98 T counts). An earlier figure, χ+ r 0.59 for the authors’ default against 0.87 (QSM.rs#333’s CI run 40786263426), compared ports that also differed in the solver, so it isolated neither change.

Note that −Σ angle / den remains a compressive estimator of χ: its gain rises with χ itself (≈0.03 at white-matter χ+, ≈0.87 in the pallidum on the QSM-CI phantom), so it flattens contrast between weak sources. That is inherent to the reference’s reconstruction, not to this port — the fitted compartment amplitudes set the gain, and they are only weakly identifiable from data that carries no sub-voxel compartment signature.

Even on its own model the fit underestimates. Data generated noise-free from the signal model above was given as the per-echo QSM, χ(TE) = −∠S(TE) / ((2/3)·γ·B0·TE) (Eq. 4). That is 400 voxels with C± 0.05–0.4, χ± ±(0.01–0.15) ppm, R2*₀ 15–25 s⁻¹, at 3 T with 16 echoes from 2.7 to 46.2 ms; no phase wraps. Measured 2026-11-01, n_inner 10:

                  r vs Eq. 6/7 at the truth   median out/truth   slope through 0   exactly 0
  χ+ (vs PCS)     0.85                        0.71               0.86              44 / 400
  χ− (vs DCS)     0.83                        0.70               0.85              48 / 400

The fitted C± and χ± correlate with the truth at only r 0.15–0.40. The 10-pass alternation stops at a fixed point short of the global minimum from the reference’s start values. On such data that minimum is exact: a global fit of all six parameters recovers the truth. MATLAB R2023b, running the reference loop (the QSM-CI decompose-qsm code) on the same voxels, gives the same table. This port is within 0.6 ppb of it in every voxel (NRMSE 0.47 % χ+, 0.29 % χ−). So the underestimate belongs to the method, not to the port.

Output mapping. The port fits what the authors fit: their ydata = squeeze(y_data(i,j,k,:))' is MATLAB’s conjugate transpose, so the fit is to conj(y), in their parameter order and with their bounds and initial values. Conjugating the data mirrors the fit (conj(S(C+, C−, χ+, χ−)) = S(C−, C+, −χ−, −χ+)), which is why their script writes DSC to its paramagnetic file (XPSC) and PSC to the diamagnetic one; the port writes χ+ = |DSC| and χ− = −|PSC| the same way. (Fitting y itself and writing PSC as χ+ gives the same maps mathematically, as earlier versions of this port did, but not the same rounding, and the solver’s path through a near-degenerate voxel follows rounding.)

Solver. Every stage is MATLAB lsqcurvefit with its default options, which the authors use: the trust-region-reflective method, ported in [crate::solvers::lsq_trr] and verified against MATLAB R2023b call by call on this fit (its module docs, Matching MATLAB: 92.6 % of the fixture’s 11340 stage calls take MATLAB’s exact path, all stage-2 calls among them). The solver of each stage is selectable (DecomposeParams::solvers), for the paper’s choice. Before that this port used a bounded Levenberg–Marquardt, which left χ± 2.9–3.3 % (NRMSE) from the authors’ code on the KKI crop below and 11.1 % / 9.0 % on the full chisep-mc phantom (astewartau/QSM.rs#329).

§Agreement with the authors’ code

Against the QSM-CI decompose-qsm container (the authors’ utils/ and fit loop compiled from licensed MATLAB R2023b, verified identical to it, QSMxT/QSM-CI#309), given the same per-echo QSM:

port	data	χ+ NRMSE (r)	χ− NRMSE (r)
unwrapping and χ cap on, LM 30 iterations, stop at 1e-10	KKI 24×24×16 crop (3 T, 5 echoes)	92.0 % (0.587)	34.9 % (0.933)
principal-branch phases, LM 30 iterations	KKI crop	8.2 % (0.9956)	10.1 % (0.9951)
LM to 400 iterations, stop at 1e-12	KKI crop	2.9 % (0.9994)	3.3 % (0.9994)
LM to 400 iterations, stop at 1e-12	`tests/data/decompose_qsm_small` (378 voxels)	1.1 % (0.99985)	11.6 % (0.9915)
`lsqcurvefit` trust-region-reflective (default)	`tests/data/decompose_qsm_small`	0.24 % (0.99999)	0.90 % (0.99995)

What remains is concentrated in voxels where one compartment collapses (an amplitude or χ within rounding of 0): there the Jacobian column of the vanished parameter is rounding noise, and MATLAB’s and Rust’s elementary functions round differently (see lsq_trr).

§Threads

Voxels are fitted in parallel (rayon, feature parallel), and the result is bit-identical at any thread count. The fit allocates nothing per voxel (see the loop). On a 6-core i5-11500, 12000 synthetic 6-echo voxels, static musl build: 1.37 s on 1 thread, 0.70 s on 2, 0.37 s on 4, 0.26 s on 6, 0.23 s on 12. (With bounded LM to 400 iterations: 28.0 s, 8.0 s, 6.4 s on 1, 4, 6 threads; lsqcurvefit’s 1e-6 tolerances end each stage far sooner than LM’s 1e-12.)

§Reference pipeline

Chen et al. 2021 (NeuroImage 242:118477, §3.6 “Image preparation and QSM reconstruction”; PMC8720043):

Data: in vivo at 3 T (GE MR750w), 0.86 × 0.86 × 1.0 mm, 16 echoes, TE1/ΔTE/TE16 = 2.7/2.9/46.2 ms. The phantom was 1 mm iso, 16 echoes 1.47–25.9 ms; the ex vivo data were 7 T.
Mask: FSL BET. Unwrapping: Laplacian. Background removal: V-SHARP.
QSM input: STAR-QSM computed for each echo (“Note that QSM of each echo, χ(t), is used instead of the average QSM across all echoes”, §2.2; “it is necessary to use TE-dependent QSM as the input”, §4.1). A host gets it from a per-echo reconstruction (crate::pipeline::per_echo). The signal model (Eq. 4) is y(t) = M(t)·exp(−i·(2/3)·χ(t)·γ·B0·t) with the per-echo χ(t), not one echo-combined χ. decompose takes that per-echo χ. (The authors’ script runs STAR-QSM per echo at padsize 12 after MRPhaseUnwrap and V_SHARP(smvsize 25) of that echo.)
Referencing: χ referenced to the whole-brain mean.
Constants: 0.5 > |χ+|, |χ−| > 0 (§2.2, DecomposeParams::chi_bound’s default), and a = 323.5 Hz/ppm at 3 T (§2.1), which the formula above reproduces and which scales with B0 by construction.

DecomposeParams’s b0 defaults to 7 T. The pipeline dispatcher overrides it from the scan metadata.

Reference: Chen, J., et al. (2021). “Decompose quantitative susceptibility mapping (QSM) to sub-voxel diamagnetic and paramagnetic components based on gradient-echo MRI data.” NeuroImage 242:118477. https://doi.org/10.1016/j.neuroimage.2021.118477 Reference implementation (QSM-CI port of Tim Ho’s open MATLAB DECOMPOSE-QSM).

§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): decompose.mask-bet, decompose.unwrap-laplacian, decompose.bg-vsharp, decompose.qsm-star, decompose.qsm-per-echo (Recommended: stated by the authors, §2.2 and §4.1), decompose.reference-mean, and at 7 T decompose.unwrap-7t and decompose.chi-bound-7t (Recommended: measured, At 7 T above).

Structs§

DecomposeParams: Parameters for decompose.
DecomposeStageSolvers: Which solver fits each stage of the alternating fit.
Fit 🔒: One voxel’s fitted parameters, in the authors’ order and sign convention (fitted to conj(y)).
Scratch 🔒: Reusable scratch for one voxel’s fit. Every buffer is fully overwritten before it is read, so a Scratch shared across voxels cannot carry state between them.
StageSolvers 🔒: lsqcurvefit workspaces for the three stages (3, 1 and 2 parameters).

Enums§

DecomposeChiBound: The bound on |χ| in DECOMPOSE’s stage-3 search ([0, −bound] to [bound, 0]), see the module docs, The χ bound.
DecomposeStageSolver: A stage solver for decompose.

Constants§

GAMMA 🔒: Gyromagnetic constant used by the DECOMPOSE reference: 42.58 · 2π — the reference implementation’s value (signalModel.m), not crate::constants. Keep it.
TAU 🔒: 2π, the period of the phase branch cut.

Functions§

broadening 🔒: a = 2π·γ·B0 / (9√3) (signalModel.m).
compartment 🔒: C · exp(−(decay + i·freq)·t) evaluated as MATLAB does it: exp of the complex exponent (−decay·t, −freq·t) is e^{re}·(cos im, sin im), then scaled by C.
decompose: DECOMPOSE source separation from per-echo QSM and the multi-echo magnitude.
dsc_model 🔒: The authors’ dscModel.m: the C_minus compartment with its frequency negated ((2/3)·(−χ−)·γ·B0), the rest at the baseline decay.
fit_voxel 🔒: The authors’ three-stage alternating fit for one voxel (DECOMPOSE_analysis.m):
matlab_log 🔒: MATLAB’s complex log(z) as (log|z|, angle z): FDLIBM’s atan2 and MATLAB’s rule for the real part ([fdlibm::log_abs]), bit for bit.
psc_model 🔒: The authors’ pscModel.m: the C_plus compartment, the rest at the baseline decay.
recon_phase 🔒: Reconstruct a source value: −Σ angle(model(TE)) / den, with (if unwrap) the phase unwrapped along the echo axis so a strong source is not folded back to zero (see unwrap_near).
signal_model 🔒: The authors’ signalModel.m at echo time t (seconds). Returns (re, im).
unwrap_near 🔒: Put ang on the 2π-branch nearest prev (only with DecomposeParams::unwrap_phase).