Module decompose 

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.

Like the QSM-CI reference, the phase is synthesized from the provided conventional QSM (χ_total) rather than re-derived from raw multi-echo phase — this isolates DECOMPOSE’s separation step. The observed complex data per echo is y = |mag_norm| · exp(−i·(2/3)·χ_total·γ·B0·TE) (magnitude normalised by its global maximum).

The per-voxel fit is a 3-stage alternating bounded least-squares, 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. 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 chi_bound, R2*₀ and the amplitudes at zero) and returns χ+ = 0 exactly, 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 original method, published at 3 T, does not meet this limit in practice — the threshold there is ~0.2 ppm.

Unwrapping is only unambiguous while the phase step between echoes stays under π, so the χ search is additionally capped at π / ((2/3)·γ·B0·ΔTE_max) — 0.31 ppm at 7 T with 8 ms spacing, inside the 0.5 ppm default DecomposeParams::chi_bound. Past that the model matches a given χ on more than one branch, the stage-3 landscape grows spurious minima, and scattered voxels fall into them; capping the search is what removes the speckle of dropouts that unwrapping alone left in the recovered pallidum.

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.

Output-mapping note. The QSM-CI reference comment claims a para/dia output swap (χ+ = |DSC|); on the qsm-forward phantom that swap anti-correlates with the ground truth, while the physically-consistent mapping used here (χ+ = |PSC|, from the paramagnetic sub-model) correlates strongly. We therefore do not replicate the reference’s swap.

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:118735. https://doi.org/10.1016/j.neuroimage.2021.118735 Reference implementation (QSM-CI port of Tim Ho’s open MATLAB DECOMPOSE-QSM).

Structs

DecomposeParams

Parameters for decompose.

LmScratch 🔒

Scratch buffers for lm_bounded, sized once for a residual of length m and reused across every stage, LM iteration and voxel. Contents are always written before they are read.

Scratch 🔒

Reusable scratch for one voxel’s fit. Every buffer is fully overwritten before it is read, so a Scratch shared across voxels (or across a rayon worker’s whole share of them) cannot carry state between them.

Constants

GAMMA 🔒

Gyromagnetic constant used by the DECOMPOSE reference: 42.58 · 2π.

LM_MAX_N 🔒

Largest parameter count any DECOMPOSE stage fits (stage 1: [C+, C−, C₀]). Keeps the LM state on the stack.

TAU 🔒

2π, the period of the phase branch cut.

Functions

cexp_decay 🔒

amp · exp(−(decay + i·freq)·t) as (re, im).

clog 🔒

Complex natural log: log(z) = log|z| + i·angle(z).

decompose

DECOMPOSE source separation from a provided QSM and multi-echo magnitude.

dot 🔒

dsc_model 🔒

Diamagnetic-only reconstruction sub-model (dscModel.m). Note the flipped imaginary sign (2/3)·(−χ−)·γ·B0 relative to signal_model.

fit_voxel 🔒

The three-stage alternating fit for a single voxel. Returns (C+, C−, C₀, R2*₀, χ+, χ−).

lm_bounded 🔒

Minimise ‖resid(x)‖² over the box [lb, ub] by Levenberg–Marquardt with a forward finite-difference Jacobian and per-step projection onto the box. Sized for the 1–3 parameter DECOMPOSE stages.

psc_model 🔒

Paramagnetic-only reconstruction sub-model (pscModel.m).

recon_phase 🔒

Reconstruct a source value: −Σ angle(model(TE)) / den, with the phase unwrapped along the echo axis so a strong source is not folded back to zero (see unwrap_near).

signal_model 🔒

Complex DECOMPOSE signal model at echo time t (seconds). Returns (re, im).

solve_small 🔒

Solve A x = b for small n (≤ LM_MAX_N) by Gaussian elimination with partial pivoting. a is row-major n×n. Returns None if singular. Only the first n entries of the inputs and the result are used.

unwrap_near 🔒

Put ang on the 2π-branch nearest prev.