Module wavesep 

WaveSep: wavelet-based susceptibility source separation.

WaveSep splits net susceptibility (χ_total / QSM) into paramagnetic (χ+, iron) and diamagnetic (χ−, myelin·calcium) sources using an R2’ map to break the para/dia degeneracy. It solves two voxel-wise data-fidelity terms under a wavelet-domain L1 (soft-thresholding) sparsity prior by proximal gradient (ISTA):

  χ+ + χ−  ≈ χ_total          (net susceptibility)
  χ+ − χ−  ≈ R2' / Dr         (static-dephasing R2', single relaxivity kernel)

with the sign convention χ+ ≥ 0, χ− ≤ 0. Dr is the static-dephasing relaxivity (Hz/ppm); the qsm-forward phantom’s single kernel is 137. WaveSep’s QSM path uses no B0 direction (unlike an STI path), so single-orientation data needs no reorientation.

Each ISTA iteration is: a gradient step on the two quadratic fidelities, a wavelet-L1 proximal step (forward db4 periodic transform → soft-threshold all coefficients by alpha·lambda → inverse), then a sign projection (χ+ = max(χ+,0), χ− = min(χ−,0)) restricted to the mask. It stops when the relative change falls below tol.

The volume is zero-padded so every axis is a multiple of 2^L (L = the periodic max decomposition level, [dwt_max_level]), matching PyWavelets’ periodization round-trip, then cropped back.

The wavelet step is PyWavelets’ own transform: wavedecn/waverecn with mode="periodization" ([Periodization::PyWavelets]: the filter centred at 2k + lf/2, and PyWavelets’ db4 taps). Until 2026-10-30 WaveSep used MATLAB’s 'per' alignment (2k + 1, the convention AMP-PE needs) and MATLAB’s db4 taps, which is a differently aligned basis: the proximal step was 3.4 % from PyWavelets’ on a smooth 32³ field, and the separation 10.2 % (χ+) / 10.0 % (χ−) NRMSE from WaveSep f3a192e on the 24³ case of tests/wavesep_parity.rs (0.11 % / 1.08 % on QSMxT/QSM-CI#302’s 164×205×205 case). With PyWavelets’ transform the 24³ case agrees to 2.8e-8 / 2.5e-8 NRMSE (the fixture’s float32 rounding), with the same 13 iterations (astewartau/QSM.rs#327).

Reference pipeline

• Data: in-vivo 3 T, 1 mm isotropic, six head orientations per subject (Fang et al. 2023, in-vivo experiments). The reference χ+/χ− were χ-separation from all six orientations.
• QSM input: an LPCNN reconstruction. The repo’s example feeds LPCNN_data/test_results/3ori/..._Ori001-Ori002-Ori003.nii.gz (data/yml/example_qsm.yml:27), and the paper tests with 1 to 3 orientations. No single-orientation, conventional-QSM input appears in the authors’ pipeline.
• R2’: one R2’ map per orientation (R2p_fn_list, data/yml/template_qsm.yml:5-8). The R2’ gradient is averaged over them (wavesep/utils/solver_wavesep_qsm.py:202-204).
• D_r: read per dataset from params.json (Dr_pos, Dr_neg; wavesep/utils/run_wavesep_qsm.py:18). The repo ships no value, and Dr_pos ≠ Dr_neg has been supported only since the 2025-04-23 update (README), with a warning. The QSM.rs default DR_SHIN_2021 is not taken from WaveSep (open question).
• Defaults: alpha = 0.2, wavelet = db4, level = None (maximum), Lambda = 0.02 (wavesep/qsm_sep.py:14-18), maxit = 100 (run_wavesep_qsm.py:32). These match WaveSepParams.
• Orientation: “All images should be in the LPS orientation” (data/yml/template_qsm.yml:19). The QSM path uses no B0 direction, so this only fixes the wavelet axes.

Reference: Fang, Z., Shin, H.-G., van Zijl, P., Li, X., Sulam, J. (2023). “WaveSep: A flexible wavelet-based approach for source separation in susceptibility imaging.” Machine Learning in Clinical Neuroimaging (MLCN), MICCAI 2023, Springer LNCS. https://doi.org/10.1007/978-3-031-44858-4_6

Reference implementation: https://github.com/ZhenghanFang/WaveSep

Agreement with WaveSep

Against WaveSep f3a192e in the verified QSM-CI wavesep image (QSMxT/QSM-CI#302), on the 24³ case of tests/wavesep_parity.rs: χ+ 2.8e-8 and χ− 2.5e-8 NRMSE (the fixture’s float32), the same 13 iterations (measured 2026-10-30). The wavelet step’s history and the earlier gap (10.2 % / 10.0 % here; 0.11 % / 1.08 % on a 164×205×205 case) are above, under the padding.

Structs

WaveSepParams

Parameters for wavesep.

Functions

pad3d 🔒

Zero-pad a column-major 3D array from from to to (extra voxels appended at the high end of each axis), matching numpy.pad(a, [(0, p-d), ...]).

pad_spec 🔒

Padded dimensions and decomposition level, matching WaveSep’s pad_spec: grow each axis to a common multiple P (starting at 16) until every padded dim is divisible by 2^L, where L = dwt_max_level(min_padded, dec_len).

prox_wavelet_l1 🔒

Wavelet-L1 proximal operator for an orthonormal transform: soft-threshold all coefficients (approximation + details) by th, in place.

soft_threshold 🔒

unpad3d 🔒

Crop a column-major 3D array from from back to to (inverse of pad3d).

wavesep

WaveSep source separation from a QSM and an R2’ map.