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.

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

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.