Module chisepnet

Expand description

χ-sepnet deep-learning χ-separation (onnx feature).

Two entry points. chisepnet_toolbox is the authors’ route, the χ-separation toolbox v1.1.3’s chi_sepnet_general (see The authors’ route below), and the pipeline’s default (ChiSepNetRoute::Toolbox). chisepnet is the network alone on the acquired grid with the caller’s QSM (ChiSepNetRoute::Native), described first.

χ-sepnet (SNU-LIST) is a 3D U-Net that maps a 3-channel patch — [QSM (χ_total, ppm), local field (ppm), R2′/Dr] each z-scored by training statistics — to paramagnetic (χ+) and diamagnetic (χ−) source magnitudes (z-scored). The network is a fixed 192×192×128 patch; we run it as an overlapping sliding window over the whole volume and average the overlaps, then de-normalize.

Recipe (mirrors the SNU-LIST / QSM-CI recon.py): z-score each channel, zero outside the mask, end-pad each dim up to the patch size, tile with a 0.75 stride plus a final flush patch, average overlaps, de-normalize, crop, and mask. Dr (114 Hz/ppm, the network’s COSMOS-referenced relaxivity) scales R2′ into the network’s ppm-equivalent input channel. Negative de-normalised magnitudes are set to 0, as the authors’ test.py does (“zero-truncation”). We return χ− as a signed (≤ 0) value to match the crate’s separation convention (chi_pos ≥ 0, chi_neg ≤ 0, chi_total).

Array orientation. The network takes no B0 direction, so whatever orientation the caller hands it is the one it reconstructs in. Kim et al. (Hum Brain Mapp 2025, doi:10.1002/hbm.70136) say what that should be, and the two halves are not equally strict. The training scans were acquired so that “the imaging slab was axial so that the z-axis was oriented along the B0 field”, and the patches were then augmented by rotating “only … within the plane perpendicular to B0 (degree: −90° to +90°) to avoid complication from B0 orientation dependent R2* in white matter”. So the slice axis carries a real requirement, while the in-plane order is deliberately covered by augmentation across the full quarter turn.

The complete order is (A-P, L-R, I-S). The in-plane half of that is not in the paper: it comes from the brain mask shipped in Data/test_file.mat — a dummy file, but the mask is real, and axis 1 is the mirror-symmetry axis in all six COSMOS orientations (Dice 0.87-0.97 against 0.78-0.88) while axis 0 is the longest — corroborated by the authors’ own axis naming in train_data_patch.py. The left/right sign stays undetermined, and the numbers below say it does not matter. Nothing in the reference code enforces any of this: neither the toolbox’s Chisep_script.m, nor the authors’ test.py, nor the QSM-CI recon.py ever permutes or flips an axis, and neither does this function — it is handed a Grid (dims and voxel sizes), not an affine, so the caller owns the orientation.

That turns out to cost very little. Measured on the QSM Reconstruction Challenge 2.0 head phantom (QSM-CI’s χ-separation dataset, which builds χ+/χ− on that head model): correlation against the ground-truth source maps, one 192³ pass, the volume relabelled and the result relabelled back.

array orientation	χ+ r	χ− r
RAS — what QSMxT feeds today	0.9480	0.8072
all 8 in-plane relabellings (slice axis still along B0)	0.9477–0.9483	0.8053–0.8078
slice axis moved to A-P (coronal-like)	0.9471	0.7914
slice axis moved to L-R (sagittal-like)	0.9485	0.7855

which is what the ±90° in-plane augmentation predicts: the in-plane order is immaterial, and even putting the slice axis across B0 costs only about 0.02 of χ− correlation and nothing of χ+, less than the training recipe alone would lead you to expect. Two of the three input channels determine χ± algebraically and carry no orientation at all (χ+ + χ− = χ_total, χ+ − χ− ≈ R2′/Dr), which is most of why. The outputs are not identical across orientations — two in-plane orientations of the same in-vivo volume differ by ~12% NRMSE, against a ~3% floor from a pure 16-voxel translation — but neither is nearer the truth, so there is no orientation to prefer. Running on an oblique acquisition is a different question, and still wants resample_to_axial: a relabelling can only fix axes that are already cardinal.

Patch size. The published graph declares that patch as a fixed input shape, but it is fully convolutional (Conv/Relu/MaxPool/ConvTranspose/Concat only), so the hosted chi-sepnet.onnx has its spatial axes re-declared as dynamic — bit-identical at the authors’ patch, and able to run smaller ones. ChiSepNetParams::patch defaults to the authors’ 192×192×128; a 32-bit host (WASM) cannot afford it — one 64-channel activation at that size is 1.2 GB — and uses WASM_PATCH instead. Patch dimensions must be multiples of 16 (four pooling levels) or the skip-connection concatenations misalign.

Weights are not bundled; the caller passes the exported chi-sepnet.onnx bytes.

§Parity

Measured against QSM-CI’s chisepnet entry (the SNU-LIST toolbox’s 240904_xsepnet.onnx and its normalisation .mat, run by QSM-CI’s recon.py) on a 36×32×28 in-vivo box of prisma-bridge-run1, at the authors’ 192×192×128 patch (QSM.rs#330):

The hosted chi-sepnet.onnx (spatial axes re-declared dynamic) gives exactly the toolbox graph’s outputs on the same patch: max |Δ| 0.
The normalisation constants above are the .mat’s values.
The only difference is zero-truncation. This function sets negative de-normalised magnitudes to 0, as the authors’ test.py does. QSM-CI’s recon.py keeps them: on this box 2 487 of 19 996 χ+ voxels and 5 χ− voxels are negative. Against the entry as it runs, that is 1.73 % NRMSE for χ+ and 0.045 % for χ−. Against the entry’s output truncated the same way it is 0.0000 % (max |Δ| 1.4e-8 ppm). QSMxT/QSM-CI#333 adds the truncation to the entry.
The SNU MATLAB toolbox truncates too. Its chi_sepnet_general.p → chi_sepnet.p (v1.1.3) was run as a black box on Bunya, MATLAB R2023b, 2026-10-30. The ONNX import was replaced by a one-layer network with constant z-scored outputs. χ+ = −2 and χ− = −3 came back as 0 and 0. χ+ = +2 came back as 2·0.03526 + 0.02591 = 0.0964. χ− = +0.5 came back as 0.0407, with χ_total = χ+ − |χ−|. So negative de-normalised magnitudes are set to 0, and positive ones pass through unchanged.
The same run shows that the toolbox’s χ-sepnet path works differently from this function. It takes the local field in Hz and R2′ (not a QSM). It makes the QSM channel itself with its own 240904_QSMnet.onnx, after rotating to B0 (CalculateRotMat) and resampling to 1 mm (interp3): chisepnet_toolbox, The authors’ route below.

§Reference pipeline

Kim et al. 2025 (arXiv:2409.14030, §2.1.1 and §2.2):

Training data: 3 T (Siemens Tim Trio), 1 × 1 × 1 mm, 6 echoes, TE 7.70–32.85 ms (ΔTE ≈ 5 ms), six head orientations, with 64³ training patches.
Field chain: magnitude by sum of squares, phase by MCPC-3D. Mask by FSL BET on the combined magnitude. Each echo unwrapped with ROMEO, background removal by V-SHARP with “a spherical mean value size of 20” per echo, then weighted echo averaging (Wu 2012) into the local field. The toolbox driver instead runs V-SHARP smvsize = 12 after echo averaging (Chisep_script.m:113-139); the two disagree (open question).
QSM channel: QSMnetME, a QSMnet trained on multi-echo data, applied to the local field (§2.2.2): [crate::inversion::qsmnetme], which chisepnet_toolbox runs.
R2’ channel: R2* by ARLO; R2 by StimFit dictionary matching of a 2D multi-echo SE with stimulated-echo correction (QSM.rs: r2_epg); R2’ = R2* − R2 with negatives set to 0. With GRE only, R2PRIMEnet (crate::relaxometry) predicts R2’.
Resolution: trained at 1 mm isotropic. The toolbox either interpolates the inputs to 1 mm and returns 1 mm maps, or runs a resolution-generalisation pipeline for QSMnetME (Chisep_script.m:188-195, chi_sepnet_general vs ..._wResolGen). The README says the networks handle resolutions above 0.6 mm. chisepnet_toolbox runs the first (the script’s resgen = false; despite the script’s comment it returns its maps on the acquired grid); chisepnet runs on the native grid with neither. The resolution-generalisation branch (chi_sepnet_general_new_wResolGen) is not ported.
Field strength: 3 T only. At 7 T the authors keep this 3 T network and feed it the 7 T local field, a QSMnetME QSM and a 3 T-equivalent R2’ from R2PRIMEnet7T (Kim J. et al., arXiv:2410.12239, §2.2.1). They name the choice of D_r at 7 T as unresolved (§4). Feeding a 7 T measured R2’ here is off-distribution.
Orientation: “For non-axial input data, the orientation was rotated such that the B0 orientation is along the z-direction” (§2.1.2), as covered above.

§Agreement with the original

tests/models_onnx/separation.rs runs it on a 36×32×28 box of in-vivo data (1 mm, 2.89 T, across the brain’s edge) against its original’s verified QSM-CI entry, SNU-LIST’s 240904_xsepnet.onnx in QSM-CI’s sliding window (chisepnet) (QSMxT/QSM-CI#316; the fixture is regenerated in Actions by tests/data/dl_separation_qsmci_small.py, from the entry with the authors’ zero-truncation, QSMxT/QSM-CI#333): at the authors’ 192×192×128 patch, χ+ NRMSE 5.8e-5 %, χ− 3.4e-5 % (max |Δ| 2.4e-7 ppm; measured 2026-11-02), float32 rounding. The 1.7 % / 0.045 % first measured was the entry’s missing truncation (above; astewartau/QSM.rs#330). The patch matters for parity: a volume smaller than the patch is zero-padded up to it, and the network’s activations in that padding differ from the convolutions’ own zero padding, so a smaller ChiSepNetParams::patch changes the result near the volume’s far edges.

§The authors’ route

chisepnet_toolbox is chi_sepnet_general of the toolbox v1.1.3, which Chisep_script.m calls with resgen = false (l. 194). The P-code was probed as a black box in licensed MATLAB R2023b (QSM.rs#377, 2026-11-03; method in [crate::utils::chisep_frame]): the inputs go into the training frame (B0 along the third axis, 1 mm, centred, padded to multiples of 16), QSMnetME makes the QSM channel there ([crate::inversion::qsmnetme]), R2′ is the measured map or, for R2*, R2PRIMEnet’s output, χ-sepnet runs once over the whole frame, and χ± come back onto the acquired grid. The steps and what was measured for each are listed on chisepnet_toolbox. Findings that settle earlier open points:

Output grid. The toolbox returns its maps on the acquired grid (interp3 'spline' at R⁻¹·x), not at 1 mm as the script’s comment says, so no resampling of our own is involved.
R2 is not clipped* before R2PRIMEnet, and R2PRIMEnet’s output reaches χ-sepnet unclipped (the toolbox clips its de-normalised copy, which it does not use).
Whole volume, no tiling: the toolbox replaces each network’s input layer to fit the frame. ChiSepNetParams::frame_patch tiles instead, for 32-bit hosts (not the authors’).
B0 direction: the toolbox does not normalise it; a non-unit vector gives it a slightly non-orthogonal rotation. QSM.rs normalises first (identical for unit vectors). For B0 = −ẑ the toolbox’s rotation is NaN; QSM.rs takes the half turn about the first axis.
7 T: the authors feed χ-sepnet a 3 T-equivalent R2′ from R2PRIMEnet7T (toolbox v1.2’s R2PNET_7T.onnx, only through SNU’s form, not in v1.1.3 and not available to us). So the R2* route refuses fields outside 3 T ± B0_TOLERANCE_T, unless ChiSepNetParams::ignore_field_mismatch; a measured R2′ is taken at any field.
Not part of the route: Chisep_script.m Tukey-windows the complex multi-echo data in k-space before everything (tukey_fac 0.4 for Siemens and GE, 0 for Philips, l. 69-73), so its mask, R2* fit and phase see windowed data. That happens upstream of the field and R2* maps this function takes, so it is the caller’s: [crate::utils::tukey_window] ports it, and the rules chi-sepnet.tukey-siemens-ge and chi-sepnet.tukey-philips suggest the authors’ value by manufacturer on this route, as the authors’ setting (Kind::Reference: shown, never preselected).

Agreement with the toolbox. On the probe volumes (64×60×48 at 0.8×0.9×1.1 mm with B0 along (0.1, 0.2, 0.97)/|·| and R2*, and at 1.1×0.7×0.9 mm with B0 along (−0.15, 0.05, 0.95)/|·| and R2′; MATLAB R2023b with the networks computed by ONNX Runtime), chisepnet_toolbox matches the toolbox’s χ+ to 5.0e-5 % / 5.2e-5 % NRMSE and χ− to 5.0e-5 % / 5.1e-5 % (max |Δ| 1.8e-7 ppm): float32 rounding, tract against ONNX Runtime. QSM-CI’s Python reproduction (algorithms/chisepnet/chisep_route.py) matches the toolbox to 1.4e-5 to 2.1e-5 %. tests/models_onnx/snu_route.rs checks chisepnet_toolbox (R2′ and R2*) and [crate::inversion::qsmnetme] against that reproduction on the in-vivo box, as acquired and with a tilted, anisotropic geometry (tests/data/dl_snu_route_small.py).

§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-sepnet.min-resolution, chi-sepnet.tukey-siemens-ge, chi-sepnet.tukey-philips.

Structs§

ChiSepNetNorm: Training z-score constants for χ-sepnet (xsepnet_train_patch_norm_factor…mat): each field is (mean, std). Inputs are normalized (x-mean)/std; outputs de-normalized y*std + mean. dr scales R2′ (Hz) into the ppm-equivalent channel (r2prime/dr - mean)/std.
ChiSepNetParams: Inference parameters for chisepnet and chisepnet_toolbox.
ChiSepNetToolboxModels: The weights chisepnet_toolbox runs: χ-sepnet, QSMnetME, and R2PRIMEnet for ChiSepNetR2::R2Star.

Enums§

ChiSepNetR2: χ-sepnet’s R2 input for chisepnet_toolbox.
ChiSepNetRoute: Which route χ-sepnet’s inputs take (ChiSepNetParams::route).

Constants§

AUTHORS_PATCH: The authors’ patch: (D, H, W) = (x, y, z).
POOL_MULTIPLE: Four pooling levels, so every patch dimension must be a multiple of this. Public (and built without onnx) so a host can check a patch, or state the rule, without restating the number.
WASM_PATCH: A patch that fits a 32-bit (WASM) heap. Peak memory, measured natively (single-threaded) on a 184×256×128 volume: 1.7 GB at this patch, 4.5 GB at AUTHORS_PATCH. On the QSM-CI χ-separation phantom it agrees with the authors’ patch to r ≥ 0.996 on both maps on ChiSepNetRoute::Native; as ChiSepNetParams::frame_patch on the default ChiSepNetRoute::Toolbox (QSMnetME and χ-sepnet both windowed) it agrees with the whole-frame passes to r = 0.990 (χ+) / 0.983 (χ−) (tests/chisep_xsepnet.rs).

Functions§

chisepnet: Run χ-sepnet χ-separation.
chisepnet_toolbox: χ-sepnet by the authors’ route: the χ-separation toolbox v1.1.3’s chi_sepnet_general (Chisep_script.m:177-194), recovered by black-box probing in licensed MATLAB R2023b (QSM.rs#377; method in [crate::utils::chisep_frame]). Unlike chisepnet it takes no QSM: