Module chi_sep_medi

Expand description

χ-separation (MEDI): the SNU χ-separation toolbox’s chi_sep_MEDI

chi_sep_medi implements chi_sep_MEDI of the χ-separation toolbox v1.1.3 (SNU-LIST, Jongho Lee’s lab; functions/corefunc/chi_sep_MEDI.p and chi_sep_MEDI_func.p), as its script Chisep_script.m calls it (l. 153-161): chi_sep_MEDI(mag, local_field_hz, r2prime, N_std, mask, params, option_data) with params.lambda = 1, params.lambda_CSF = 1, option_data.qsm, option_data.mask_CSF 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: 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 ...')), and replaying the operators captured in those workspaces (their function handles stay callable after load). A MATLAB re-implementation built from that reproduces chi_sep_MEDI to 3·10⁻¹⁶ relative on six phantoms (see Agreement with the toolbox).

Earlier QSM.rs versions ported Shin et al. 2021’s Supplementary Methods instead, with their own weights, λ and no QSM, CSF mask or noise map; the toolbox’s algorithm below differs in every step after the forward model.

§The algorithm

The unknowns are x+ = χ+ ≥ 0 and x− = |χ−| ≥ 0 in ppm. With D the dipole convolution in Hz ((CF/10⁶)·F⁻¹ D F) and D_r,± the relaxometric constants (Hz/ppm):

Defaults (chi_sep_MEDI): B0 = CF/42.58·10⁶; D_r,pos = D_r,neg = 137·(B0/2.9) (toolbox_default_dr); pad_size 5, tol_norm_ratio 0.03, max_iter 30, cg_max_iter 100, cg_tol 0.05 (see ChiSepMediParams); and, not parameters, at least CHI_SEP_MEDI_MIN_ITER = 6 iterations, ε = 10⁻⁵ and wei_mag = 1/(2π).
Edge masks on the input grid: wG = gradient_mask(mag, mask, 0.9) (MEDI’s gradient_mask.m: |∇(mag·mask)| thresholded so that 90 % of the mask’s voxel count of gradient components lie above the level; 1 = not an edge). The R2′ edge mask wG_r2p is the same magnitude mask.
Padding: every volume is zero-padded by pad_size voxels at both ends of the third axis only; wG, wG_r2p are padded with 1. Everything below runs on that grid. Exactly zero R2′ voxels (including the padding) become eps (2.2·10⁻¹⁶).
Weights and targets: m = dataterm_mask(N_std, mask) (MEDI’s: mask/N_std, 0 where not finite, divided by its mean over the mask); R2′ residual weight wei_mag·m, field residual weight m. Targets: R2′·mask where R2′ > 0 (0 elsewhere) and the local field in Hz.
Dipole kernel: MEDI’s dipole_kernel(N, vs, B0_dir, 'kspace'): the k-grid (−N/2 : N/2−1)/(N·vs) per axis, D = 1/3 − (k·B0)²/|k|² (0 at k = 0), fftshifted. On an axis with an odd number of voxels that grid steps through half-integers and misses k = 0, so the kernel is sampled half a frequency bin off and the convolution has an imaginary part; it is kept as the toolbox keeps it (complex arithmetic throughout, real part taken where the toolbox takes it).
Initial guess from the QSM χ (option_data.qsm): the voxelwise solution of D_r,pos·x+ + D_r,neg·x− = R2′, x+ − x− = χ, each clamped at 0: x+ = max((R2′ + D_r,neg·χ)/(D_r,pos + D_r,neg), 0), x− = max((R2′ − D_r,pos·χ)/(…), 0) (with the raw R2′, negative values included). The iterate is that, masked. A copy with the CSF voxels set to 0 is what the first iteration’s L1 weights see (a quirk of the toolbox: the iterate itself keeps the CSF values).
Gauss-Newton while (res > tol or iter < 6) and iter < max_iter:
- L1 weights V_s = 1/√(|wG·∇x_s|² + ε) for s = +, −, and V_t = 1/√(|wG·∇(x+ − x−)|² + ε) (∇: MEDI’s fgrad, forward differences scaled by the voxel size, 0 at the last voxel).
- Gradients b_s of `‖wei_mag·m·(D_r,pos·x+ + D_r,neg·x− − R2′)‖² + ‖m·(D x+ − D x− − f)‖²
  - λ Σ_s ‖wG V_s^½ ∇x_s‖² + 2λ ‖wG V_t^½ ∇(x+ − x−)‖² + λ_CSF Σ_s ‖LT(|x_s|)‖²(as the toolbox evaluates them:2·D_rᵀ M²·residual, ∇ᵀ(wG V wG ∇·)with∇ᵀthe exact adjoint offgrad), where LT(v) = M_CSF·(v − mean(v over M_CSF))`.
- One CG solve per source (cgsolve.m, MEDI’s: zero start, last iterate, residual recomputed every 50 iterations): A_s dx_s = −b_s, where A_s is the same expression’s linearisation in x_s with the other source’s step set to zero, CSF term included (λ_CSF·LT(LT(|dx|)), nonlinear through the abs; the toolbox rebinds its reg_p/reg_n to include it before building the CG operators).
- x_s ← max(x_s + dx_s·mask, 0); res = ‖Δx‖/‖x‖ over both sources inside the mask eroded by a sphere (imerode(mask, strel('sphere', erode_size)), voxels outside the array counted as inside), with Δx the projected step. erode_size = min(round(min(effective_mask_size)/10), 5) voxels, effective_mask_size being the number of slices along each axis that hold any mask voxel (5 on a head); with erode_size 0 the whole padded grid counts. The iterate with the smallest res so far is kept.
Output: that iterate, cropped back to the input grid: x_para = x+, x_dia = x− (≥ 0), x_tot = x+ − x−. QSM.rs returns (χ+, χ− = −x−, χ+ + χ−).

§Units

The toolbox takes the local field in Hz; chi_sep_medi takes it in ppm and multiplies by cf/10⁶ (so the result depends on ChiSepMediParams::cf only through D_r’s default and the relative weight of the two residuals). R2′ is in Hz, the QSM and the output in ppm.

§Agreement with the toolbox

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

Phantoms (tests/chi_sep_medi_snu_parity.rs, fixture from tests/data/chi_sep_medi_snu_small.m): 40×36×28 at 3 T with a CSF mask; 33×29×21 at 7 T, 0.8×0.9×1.5 mm, oblique B0, with a CSF mask; negative R2′, D_r 120/150, λ = 2, no CSF mask. Same Gauss-Newton iteration counts (7, 6, 6), χ+ and |χ−| to relative L2 ≤ 1.1·10⁻¹⁵, max |Δ| ≤ 8·10⁻¹⁷ ppm.
The verified QSM-CI chi-sep-medi container (v4, identical to licensed MATLAB; QSMxT/QSM-CI#351): option_data.qsm = χ_total, N_std = Fit_ppm_complex‘s dp1 from the echoes’ magnitudes, mask_CSF = [], CF = 42.5774 MHz/T·B0. QSM.rs’s N_std (crate::utils::fit_ppm_complex_noise) equals the container’s bit for bit in every case.
- a 24×22×18 phantom whose voxels decay at different rates (tests/chi_sep_medi_qsmci_parity.rs, regenerated in Actions by tests/data/chi_sep_medi_qsmci_small.py): χ+ and |χ−| to relative L2 2.6·10⁻⁸, the container’s float32 output;
- in vivo (the EMTP travelling-head scan of QSM-CI’s fidelity check, 3 T, 184×256×128, 1.39 M mask voxels; its local field, R2′ and χ): 6 iterations each, the same res at every iteration (to the 4 digits the toolbox prints); χ+ and |χ−| to relative L2 2.5·10⁻⁸, max |Δ| 3·10⁻⁸ / 9·10⁻⁸ ppm, the float32 output (measured 2026-11-01, examples/chi_sep_file.rs);
- QSM-CI’s χ-separation phantom (chisep-mc, 7 T, 164×205×205, same inputs): 6 iterations each, res equal at the first; χ+ relative L2 1.7·10⁻³ (max |Δ| 8.3·10⁻³ ppm), |χ−| 9.3·10⁻³ (7.9·10⁻³ ppm), r ≥ 0.9999. That residual is rounding, amplified: on this phantom every inner CG solve runs 42-50 iterations (5-10 in vivo), and the result moves as much under changes at the last bit. The toolbox itself, run with 4 instead of 8 threads, moves by 3.7·10⁻³ / 1.13·10⁻² (more than QSM.rs differs from it), and QSM.rs moves by 2.1·10⁻³ / 4.5·10⁻³ when CF changes by one ulp. QSM-CI’s chisep-medi / chisep-medi-qsmrs parity pair checks this phantom.
Before QSM.rs#382 these cases ran with N_std = 1/|mag| (in vivo relative L2 8·10⁻¹³ against licensed MATLAB with that input), and the convergence test eroded the mask by 3 voxels instead of the toolbox’s erode_size (5 on these scans; see step 7), which changed res by about 5 % but neither the iteration counts nor the results here.

§Performance

Each CG iteration costs four 3-D FFTs (two dipole convolutions, complex) and two fused gradient-adjoint passes; a Gauss-Newton iteration is two CG solves. In vivo (184×256×138 padded, 6 iterations) chi_sep_medi takes 106-133 s on 6 threads with a 2.7 GB peak (the toolbox: 286 s on 8); on chisep-mc (164×205×215, CG solves of ~45 iterations) 282-292 s (the toolbox: 931 s). Memory is about 45 f64 per padded voxel.

§Deviations

An all-zero CSF mask means no CSF term. The toolbox applies the CSF term whenever option_data.mask_CSF is not empty; with a mask that has no voxels, the mean over it is NaN, NaN·λ_CSF is NaN even for λ_CSF = 0, and every step is lost (the solver returns its initial guess). Its own GUI passes [] when it has no CSF mask, which turns the term off. The CSF mask is a required input here, as in Chisep_script.m; passing an all-zero mask reproduces the GUI’s no-CSF behaviour (mask_CSF = []).
No wall-clock limit. The toolbox also stops after params.time_limit (default 1800 s) of iterations. That makes its output depend on the machine; it is not reproduced.
Options of chi_sep_MEDI that Chisep_script.m does not use are not exposed: option_data.wG/wG_r2p (precomputed edge masks), mask_FastRelax/mask_SlowRelax (voxels whose R2′ residual weight drops to 0.1; both default to empty), and the TKD initialisation it uses without option_data.qsm (params.init_option).

§Reference pipeline

What Chisep_script.m (toolbox v1.1.3) runs before chi_sep_MEDI, with its line numbers:

Tukey k-space windowing of the complex multi-echo data, factor 0.4 for Siemens/GE and 0 for Philips (l. 70-73); mag = root-sum-of-squares over echoes (l. 75).
N_std from Preprocessing4Phase(mag, phase) (l. 76), which is MEDI’s Fit_ppm_complex noise map (probed: equal to 3·10⁻¹⁸; QSM.rs crate::utils::fit_ppm_complex_noise).
Mask: MEDI’s BET on mag (l. 82); R2* by ARLO with more than three equally spaced echoes, NNLS otherwise (l. 87-93); mask_CSF = extract_CSF(r2star, mask, voxel_size) (MEDI’s extract_CSF.m, l. 96; QSM.rs crate::inversion::extract_csf).
ROMEO with phase-offset correction and weighted echo averaging, W = TE·exp(−TE/40 ms) (l. 99-123); V-SHARP smvsize 12 (l. 136-139); the local field in Hz (l. 139).
QSM: STI Suite QSM_iLSQR(..., 'padsize', [12 12 12]) (l. 145), passed as option_data.qsm (l. 158).
R2′ is the user’s (the script’s input; its README offers a pseudo-R2 map otherwise).

The toolbox’s GUI (WrapperChisepMEDI.p, probed with a stand-in app object) differs: it passes mask_CSF = [] (no CSF term), the RSS of the Tukey-windowed magnitude, λ = λ_CSF = 1, and option_data.qsm only when its “QSM initialisation” option is on. Its N_std is the script’s: Fit_ppm_complex’s dp1 from the multi-echo magnitude, zero outside the mask (QSM.rs#382; probed in licensed R2023b on phantoms whose voxels decay at different rates, where the GUI’s and Preprocessing4Phase’s N_std both equal dp1 to 1.1·10⁻¹⁵). An earlier reading of the GUI’s as c/|mag| came from a phantom with uniform decay, where dp1 is proportional to 1/RSS; with per-voxel decay dp1·RSS varies (0.72-0.87 there).

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

§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-medi.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§

ChiSepMediInputs: The volumes chi_sep_medi takes, as chi_sep_MEDI(mag, local_field_hz, r2prime, N_std, mask, params, option_data) takes them. All are on the grid, column-major.
ChiSepMediParams: Parameters for chi_sep_medi. The defaults are the toolbox’s (see the module docs).
Problem 🔒: The arrays of one run, on the padded grid.
Workspace 🔒: Buffers reused across operator applications.

Constants§

CHI_SEP_MEDI_MIN_ITER: Gauss-Newton iterations chi_sep_MEDI always runs (its min_iter, not a parameter), even once the relative change is below ChiSepMediParams::tol.
EDGE_PERCENTAGE 🔒: MEDI gradient_mask’s percentage as the toolbox passes it: the number of gradient components marked as edges, as a fraction of the mask’s voxel count.
EPS_L1 🔒: ε in the L1 weights 1/sqrt(|wG·∇χ|² + ε) (chi_sep_MEDI_func’s eps__).
ERODE_MAX 🔒: The largest radius, in voxels, of chi_sep_MEDI’s convergence-mask erosion.
TOOLBOX_DR_REFERENCE_TESLA: The field strength, in tesla, at which the toolbox’s default relaxometric constant is super::DR_SHIN_2021 (137 Hz/ppm): chi_sep_MEDI and chi_sep_iLSQR default to D_r = 137·(B0/2.9) for both sources, B0 = CF/42.58e6, when params.Dr_pos/Dr_neg are not given (probed: 141.7241379310345 at CF = 127.74 MHz, bit for bit). Shin et al. 2021 measured 137 on a nominal 3 T (2.89 T) scanner; the toolbox scales it linearly with the field.
TOOLBOX_GYRO_HZ_PER_T: The gyromagnetic ratio the toolbox’s chi_sep_MEDI and chi_sep_iLSQR use to get the field strength from params.CF (their gyro_), in Hz/T. It differs from crate::constants::GAMMA_BAR_HZ_PER_T (42.576 MHz/T); kept so the field-dependent defaults match the toolbox.
WEI_MAG 🔒: Weight of the R2′ residual against the field residual, both in Hz (wei_mag = 1/(2π)).

Functions§

chi_sep_medi: χ-separation by the SNU χ-separation toolbox’s chi_sep_MEDI (see the module docs).
convergence_mask 🔒: chi_sep_MEDI’s convergence mask (Mask_erode): the mask eroded by a sphere of radius erode_size = min(round(min(effective_mask_size)/10), 5) voxels, where effective_mask_size counts, along each axis, the slices that hold any mask voxel (probed in licensed R2023b: two 10³ cubes far apart give [20 20 20]; voxel sizes and pad_size play no part; a head gives 5). A radius of 0 (an axis with fewer than five such slices) makes it the whole padded grid, as in the toolbox. round is MATLAB’s, halves away from zero (15 slices give 2, 25 give 3), as Rust’s.
erode_sphere 🔒: imerode(mask, strel('sphere', r)): a voxel stays if every voxel of the ball of radius r voxels around it is in the mask, counting voxels outside the array as in the mask (MATLAB’s border rule for erosion).
toolbox_default_dr: The toolbox’s default D_r in Hz/ppm at central frequency cf (Hz): 137·((cf/42.58e6)/2.9), in that order of operations.
toolbox_dipole_kernel 🔒: MEDI’s dipole_kernel(matrix_size, voxel_size, B0_dir, 'kspace'), as the toolbox builds it: the k-grid (-N/2 : N/2-1)/(N·vs) per axis (half-integer steps on an odd axis), D = 1/3 − (k·B0)²/|k|², NaN (k = 0) set to 0, then fftshifted (circular shift by floor(N/2)).