Introduction

WBFR : a wavelet based approach to background field removal QSM is a reproducible computational method for the extraction of tissue magnetization field maps within the spinal chord developed as part of the 2026 of Brainhack School.

WBFR aims to offer the field of quantitative susceptibility mapping an alternative to the more established background field removal algorithms. Broadly speaking, the traditional BFR algorithms (PDF, LBV, and the SHARP family of algorithms) are all based on the assumption of a harmonic background field within the cord ROI.

WBFR is truly novel in the sense that it does not make that assumption but leverages the multi-scale analysis capabilities of the discrete wavelet transform.

Background

Magnetic susceptibility is a property of materials that describes their reaction to being submerged in a magnetic field. The magnetic susceptibility distribution within some volume of tissue can reveal important information about its chemical composition. Myelin — a fatty material that surrounds and protects neurons — has negative magnetic susceptibility, whereas iron has positive magnetic susceptibility. Mapping the concentrations of myelin and iron within the spinal cord and the brain has proven immensely helpful in the diagnosis and prognosis of some degenerative diseases — most importantly, multiple sclerosis (MS).

Quantitative Susceptibility Mapping

QSM is a method of processing MRI data in a way that creates a map of the magnetic susceptibility distribution within a given region of interest (ROI). The steps below describe how QSM takes the phase and magnitude data produced by MRI scanners to output a susceptibility field map within the ROI (a local field map or chi-map).

Acquire phase and magnitude data $\implies$ Estimate the total field $\implies$ BFR i.e. remove contributions from regions outside the ROI $\implies$ Perform a dipole inversion to extract a local field map.

BFR is arguably the most important step in the QSM pipeline. It is where most of the information about the local field is either gained (by the successful removal of background artefacts) or lost.

Modelling the Background Field Contributions Within the ROI

A classical result of electromagnetism tells us that for any bounded domain $\Omega$ with a sufficiently regular boundary $\partial \Omega$, a magnetic field $\mathbf{B}(\mathbf{r})$ with sources outside $\Omega$ obeys the Laplace equation

$$\nabla \mathbf{B}(\mathbf{r}) = 0, \quad \forall \mathbf{r} \in \Omega$$

with $\mathbf{B}|_{\partial \Omega} = f$ for some sufficiently smooth function $f$ (this second equation is called the Dirichlet boundary condition). All of the more established BFR algorithms use some variation of this theoretical basis. More precisely, they all leverage a property of harmonic functions (i.e. any twice continuously function $\mathbf{B} : \mathbb{R}^3 \to \mathbb{R}$ such that $\nabla B = 0$) called the Spherical-Mean Property (SVM).

For any point $x$ in the domain $\Omega$ (where $\mathbf{B}$ is supposed to be harmonic) and a sphere $S(x,r)$ of radius $r$ and centered at $x$, if the sphere $S$ is still entirely contained in $\Omega$, then the value of $B$ at $x$ can be inferred by computing the mean value of $B$ on the sphere (here we assume $S$ to be a hollow sphere i.e. just a surface).

$$B(x) = \frac{1}{\text{Volume enclosed by }S} \int_{S} B(s),dS$$

The reasoning behind using the SVM for BFR is mathematically sound, but it clashes with the practical, physical limitations of MRI imagery. Just as the smallest unit of precision of a 2D image is the pixel, the smallest unit of precision of an MRI is the voxel (volumetric pixel) — a cube of side length approx. $0.4 \text{ mm}$ on most machines. A sphere, in the word of discrete signal processing is a type of 3D Swiss flag. Its radius cannot be smaller that one voxel side length. As opposed to the generic sphere $S$ in our statement of the SVM which can be arbitrarily small. Here is a visual representation of how this limitation can interfere with estimation of the background within our ROI.

As mentioned in our statement of the SVM, we need $S$ to be entirely contained in $\Omega$ for the property to hold. $S$ cannot come into contact with the boundary region (dashed red line in the figure). Hence a region one pixel wide on either side of the boundary represents a kind of "danger zone" for harmonicity-based BFR.

The easiest, most drastic way of dealing with what happens in the boundary region is simply to erode the mask — removing the boundary region from the ROI entirely. For a large ROI such as the brain, eroding the mask can be done without sacrificing too much valuable information. The cord, on the other hand, is $1 \text{ cm}$ in diameter. It needs all the voxels it can get !

What tells most classical BFR algorithms from one another is the assumptions made on the behavior of the background within the danger zone that is the boundary region.

Assumptions	Acronym	Algorithms
NO magnetic field Sources on the boundary of the ROI	NOS	SHARP (V-SHARP)
Assumption of No Harmonic internal and boundary i.e. that the value of the background field in the boundary region is equal to that of the total field.	NOHA	LBV
Minimization of an Objective function Involving a Norm: Takes a family of approximations for the background field and selects the one with the smallest $L^2$ norm.	MOIN	PDF

Research at Neuropoly has concluded that PDF is the best-performing algorithm for BFR of the spinal cord.

The Wavelet Transform: A Change in Paradigm

Another difference between QSM of the brain and QSM of the spinal cord is the greater diversity of sources for the background field of the cord. The lungs, the trachea, tissue interface between the cord and CSF, the vertebrae all create contributions that vary with sometimes indiscernible frequencies.

To differentiate between low frequency contributions like the spine, intermediate frequency like the trachea, and high frequency like tissue interface we look toward the wavelet transform. A mathematical tool that decomposes signals into different spatial and frequency scales.

A wavelet takes the form $\psi_{mn}(x)= 2^{-m/2} \psi(2^{-m}x-n)$ where $\psi$ is called a mother wavelet. Loosely speaking good $\psi$ functions are continuous, real functions, with compact support (with continuity, this implies that they are bell-shaped), and have enough sharp variations that they can't be described by low order polynomials (see vanishing moment property for a more rigorous caracterization). $m$ scales the mother wavelet in frequency and $n$ shifts it in time.

An image can be decomposed into an infinite sum of terms of the form $a_{mn} \cdot \psi_{mn}$ where $a_{mn}$ is a real (wavelet) coefficient à la Fourier transform.

Here an image of a man holding a camera is decomposed using the Discrete Wavelet-Transform (DWT). The output of the transform is always of the form $(cA,(cH_{1},cV_{1},cD_{1}), \dots ,(cH_{L},cV_{L}, c_D{L}))$ where the approximation coefficient $cA$ captures the main features of the image (high spatial frequency and low time frequency) coefficient tuple $(cH_{i},cV_{i},cD_{i})$ are respectively horizontal, vertical and diagonal detail coefficients at decomposition level $i$ in an $L$-level decomposition. The figure above is an example of a $1$ level wavelet transform using the Haar wavelet.

Methodology

To mesure the effectiveness of the wavelet transform as a means of removing background field in the spinal cord, we need a total fieldmap on which to apply the DWT and local field that we know is accurate, a ground truth, with which we'll compare the output of the transform.

An accurate tissue field is impossible to guarantee in in-vivo data. That is why we rely on the whole-body phantom : a numerical "crash-test dummy" for MRI methods developed by the Neuropoly lab. In a nutshell, a 3D model of the human body was made along with a map of all the magnetic properties of tissue relevant for MRI. The two put together can be used to accurately simulate MRI acquisition.

Creating the Total Field

Using the tissue-to-MR module in Neuropoly's shimming toolbox, and passing the whole-body phantom as input we create maps for the main MRI properties of tissue $T_{1},T_{2}, T^*_{2}, \text{PD}, \chi$ where $PD$ is the proton density and $\chi$, the quantity most important for our inquiry, the magnetic susceptibility. The module susceptibility-to-fieldmap-FFT takes a susceptibility map as input and outputs its associated field map.

$$\text{whole-body phantom} \implies \chi \text{ map} \implies B_{0}(\chi \ast D) \implies \text{total field}$$

Creating the Local Field

If we isolate the spinal cord by applying a mask that only preserves the white matter and grey matter of the phantom and apply the steps to create the total field map but to this restricted region of the phantom, we get our local field ground truth. $$\text{whole-body phantom} \implies \text{spinal cord phantom} \implies \chi \text{ map} \implies B_{0}(\chi \ast D) \implies \text{local field}$$

Optimization Loop

STEP 0 (Setup): Divide the local field into slices along the $z-$axis. Choose a wavelet and a level of decomposition. For a fixed level $L$ of decomposition, there are $3L+1$ coefficients (bands) in the 2D DWT : $1$ approximation coefficient $cA$;

$L$ vertical detail coefficients $V_1, \dots V_{L}$;

$L$ horizontal detail coefficients $H_1, \dots,H_{L}$;

$L$ diagonal detail coefficients $D_1, \dots, D_{L}$.

Below an example with $L=2$ with a mask of the cord dilated by 20 voxels Initialize the list candidate_bands to be the entire set of $3L$ detail coefficients of the DWT, initialize an empty list of selected_bands, an empty list of tested_bands, and a number baseline_score to +Inf.

STEP 1 : Select a candidate band out of the candidate_bands that were not tested.

STEP 2 : Compute the inverse DWT using trial_bands = selected_bands + candidate.

STEP 3 : Compute the quality loss function on this trial image.

STEP 4 : If trial_score < baseline : move the candidate from the candidate_bands to the selected bands. Set baseline = trial score. else : move candidate from candidate_bands to tested_bands. Repeat steps 1 through 4 until there is no noticeable improvement in the baseline score.

Repeating this optimization loop for all 106 discrete wavelets available in the PyWavelets library and every decomposition level compatible with our image size (between 2 and 5 depending on the choice of wavelet) proved to be too computationally demanding for a single desktop computer. That's we ran a multithreaded optimizer on the Rorqual HPC cluster provided by Alliance Canada. Here is an anonymized version the sbatch-compatible script that was used to run our optimizer on Rorqual :

#!/bin/bash

#SBATCH --account=[grp-name]
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=16
#SBATCH --mem-per-cpu=4000M
#SBATCH --time=02:00:00
#SBATCH --job-name=wvt_experiments_opt
module load python/3.11 nibabel/5.2.0 scipy-stack
source ~/links/projects/grp-name/username/wvt_experiments/env/bin/activate
python ~/links/projects/grp-name/username/wvt_experiments/src/wavelet_band_selector_gt_parallel.py

submit.sh

Our choice of optimization loop is by no means efficient, nor is its solution guaranteed to be optimal, but similar greedy algorithms have proven to be useful heuristics in a number of real-world setting. We thought it would be a good place to start, and it did lead to some interesting results.

Preliminary Results

TLDR : It works wonderfully, and I have no earthly idea why !

What to Look For: A Crash-Course in Spine and Spinal Cord Anatomy

The spine is divided into four regions: cervical, thoracic, lumbar and sacral. Slices of vertebrae C3, T2, and L1 provide a fair sample of the diversity of background field sources present in the spinal cord. The portion of the cord located in the cervical spine experiences background field contributions from the trachea and the lungs, the thoracic portion experiences contributions from the lungs and from the geometry of the spine, the lumbar portion mainly experiences contributions from the geometry of the spine.

The cord is composed of white matter and grey matter. Within a slice of cord, we expect to see white matter surrounding the characteristic "grey-matter butterfly".

Cervical Spine Slice

This is the region where the expect the greatest amount of perturbation, as it is close to the edge of our mask. Notice that while the background field isn't entirely removed in the WBFR result, the image has considerably more contrast than the corresponding slice in the PDF output.

Thoracic Spine Slice

We've struck gold ! WBFR preserves the most relevant structures of the cord —namely, the grey-matter butterfly —with again, more contrast than PDF.

Lumbar Spine Slice

We're beginning to notice a pattern here : Although the background isn't entirely removed, WBFR preserves contrast and structure. It also preserves more detail around the boundary of the cord than PDF does.

Metrics: Lies, Damned Lies, and Statistics

Using three standard metrics of image processing — the fields of compression and denoising in particular — we try to quantity the differences between the PDF output and that of WBFR. Mean-Squared-Error (MSE)

$$\text{MSE} = |\text{ground truth} - \text{BFR output}|^2$$

Correlation

$$r_{xy}=\frac{\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}};\sqrt{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}}$$

Where $x_i$ is a voxel value of the ground truth and $y_i$ the corresponding voxel value in one of our estimates.

Structural-Similarity Index Measure (SSIM)

$$\text{SSIM}(x,y)=\frac{(2\mu_{x}\mu_{y}+c_{1})(2\sigma_{xy}+c_{2})}{(\mu_{x}^{2}+\mu_{y}^{2}+c_{1})(\sigma_{x}^{2}+\sigma_{y}^{2}+c_{2})}$$

Where $c_1 = (0.01L)^2$ and $c_2=(0.03L)^2$, with $L$ being the dynamic range of the pixel values. Typically, $L = 2^{\text{number of bits per pixel}}-1$.

Comparison with the Ground Truth and PDF

It is quite surprising to find that WBFR doesn't perform particularly well in any of the standard metrics.

And yet, a qualitative analysis of WBFR by Neuropoly experts indicates that it keeps more of the clinically relevant details. It seems WBFR has a structure preserving property that we haven't been able to quantify. This warrants further investigation - which is, unfortunately, outside of the scope of this project.

Reproducibility

• This repo follows the project structure guidelines set in The Good Data-science Handbook
• A virtual environment wvt_env_2 can be used to manage dependencies.
• A copy of the whole-body phantom can acquired by requesting access to the Neuropoly through osf
• Processing of the whole-body phantom was done using the open software in Neuropoly's Spinal-Cord Toolbox.

Quick Start

Run WBFR on one sample subject.

1. Install dependencies

python3 -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt

2. Run a minimal example (subject 02)

mkdir -p output
PYTHONPATH=src python - <<'PY'
import nibabel as nib
from dwt_utils import load_nifti_as_array, dwt_bfr_volume

field_path = "data/sub-02/sub-02_desc-b0_fieldmap.nii.gz"
mask_path = "data/sub-02/sub-02_desc-SC_mask.nii.gz"
out_path = "output/sub-02_B_loc_dwt_quickstart.nii.gz"

field = load_nifti_as_array(field_path)
mask = load_nifti_as_array(mask_path) > 0
affine = nib.load(field_path).affine

B_loc = dwt_bfr_volume(
	field,
	mask,
	wavelet="rbio1.3",
	level=3,
	keep_detail_levels=(3, 2),
	sig_ext_mode="zero",
)

nib.save(nib.Nifti1Image(B_loc, affine=affine), out_path)
print("Saved:", out_path)
PY

3. Expected Result

A local-field estimate is saved at output/sub-02_B_loc_dwt_quickstart.nii.gz.

Sources

1. Lee, G. R., Gommers, R., Waselewski, F., Wohlfahrt, K., & O'Leary, A. (2019). PyWavelets: A Python package for wavelet analysis. Journal of Open Source Software, 4(36), 1237. https://doi.org/10.21105/joss.01237

2. van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., Boulogne, F., Warner, J. D., Yager, N., Gouillart, E., Yu, T., & the scikit-image contributors (2014). scikit-image: image processing in Python. PeerJ, 2, e453. https://doi.org/10.7717/peerj.453

3. De Leener, B., Lévy, S., Dupont, S. M., Fonov, V. S., Stikov, N., Collins, D. L., Callot, V., & Cohen-Adad, J. (2017). SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data. NeuroImage, 145, 24–43. https://doi.org/10.1016/j.neuroimage.2016.10.009

4. Bollmann, S., Rasmussen, K. G. B., Kristensen, M., Blendal, R. G., Østergaard, L. R., Plocharski, M., O'Brien, K., Langkammer, C., Janke, A., & Barth, M. (2017). DeepMRI: quantitative susceptibility mapping using deep learning. NMR in Biomedicine. https://doi.org/10.1002/nbm.3604

5. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. (2020). Array programming with NumPy. Nature, 585, 357–362. https://doi.org/10.1038/s41586-020-2649-2

6. Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in Science & Engineering, 9(3), 90–95. https://doi.org/10.1109/MCSE.2007.55

Licence

• This project falls under the GNU licence.