Susceptibility distortion in EPI

In this tour you correct the susceptibility distortion of a measured EPI pair. In an echo-planar image, the phase-encoding direction is sampled once per echo spacing rather than once per dwell time, so its bandwidth per pixel is only a few tens of hertz. A spin off resonance by \(\Delta f\) is displaced along phase encoding by \(\Delta f\) divided by that bandwidth, which near the frontal sinus and the petrous bone comes to several millimetres at 3 T. The orbitofrontal cortex and the temporal poles are compressed or stretched, and signal piles up where neighbouring voxels are displaced onto the same place.

The displacement changes sign with the direction in which you traverse k-space. Two acquisitions with opposite phase-encoding polarity, blip-up and blip-down, are distorted in opposite directions, and the displacement field that brings them into register is the correction [1]. You correct such a pair, measured at 3 T and published on OpenNeuro [3], with bartorch.tools.correct_susceptibility(), which runs PyHySCO [2], and you check the estimated displacement against a gradient-echo field map of the same subject.

Prerequisites: lesson 1, 1. Your data as tensors.

Learning objectives

• Read a BIDS EPI series and its sidecar, and compute the bandwidth per pixel along phase encoding and the displacement per hertz of off-resonance.

• Tell the direction of the displacement in anatomical terms from the phase-encoding direction and the image orientation.

• Recognise the compression, stretching and pile-up of susceptibility distortion, and their reversal between the two polarities.

• Estimate the displacement field from a reversed phase-encoding pair and correct both images, intensity included.

• Check the estimate against an independent field map, and the correction by the agreement of the two corrected images.

The data

The dataset holds gradient-echo EPI series of one subject acquired on a Siemens Prisma at 3 T: a 64 x 64 matrix, 44 axial slices, 3.75 x 3.75 x 4 mm, with the same protocol and shim and opposite phase-encoding polarity, and a dual-echo gradient-echo field map on the same grid. Each EPI series has five volumes, and you average them to raise the signal-to-noise ratio. read_nifti() returns the volumes as (volumes, z, y, x) together with the affine from voxel indices (x, y, z) to RAS millimetres. The JSON sidecars hold the timings of the readout. The download into CACHE is hidden.

import json

import torch

import bartorch.tools as bt

from bartorch.io import read_nifti
from nibabel.orientations import aff2axcodes

blip_up_file = CACHE / "sub-1/fmap/sub-1_dir-PA_epi"
blip_down_file = CACHE / "sub-1/func/sub-1_task-rest_acq-AP_bold"
blip_up = read_nifti(f"{blip_up_file}.nii.gz")
blip_down = read_nifti(f"{blip_down_file}.nii.gz")
sidecar_up = json.loads(Path(f"{blip_up_file}.json").read_text())
sidecar_down = json.loads(Path(f"{blip_down_file}.json").read_text())

affine = blip_up.affine
spacing = affine[:3, :3].norm(dim=0).flip(0)
voxel_mm = tuple(float(v) for v in spacing)  # (z, y, x)

up = blip_up.image.mean(0)
down = blip_down.image.mean(0)

(5, 44, 64, 64) volumes, voxel (z, y, x) (4.0, 3.75, 3.75) mm
voxel axes x, y, z point to LAS
blip-up   BOLD_3mm_matched_0.55esp_topup: j
blip-down BOLD_3mm_matched_0.55esp: j-

The voxel axis x points to the subject’s left and y to anterior, so PhaseEncodingDirection j is posterior to anterior and j- anterior to posterior. Phase encoding is along the anterior-posterior axis of the head, and the readout is left-right. The figures show the slices with anterior at the top and the subject’s left on the right.

Displacement along phase encoding

You acquire \(N\) phase-encoding lines at the effective echo spacing \(\Delta t_{esp}\), the echo spacing divided by any parallel-imaging acceleration. Off resonance, the phase you accrue is linear in \(k_y\), and a linear phase in k-space is a shift in the image:

\[d = \frac{\Delta f}{\mathrm{BW}_{PE}}\ \text{voxels}, \qquad \mathrm{BW}_{PE} = \frac{1}{N\, \Delta t_{esp}},\]

with \(\mathrm{BW}_{PE}\) the bandwidth per pixel along phase encoding. BIDS writes \((N - 1)\,\Delta t_{esp}\) as TotalReadoutTime. For a positive \(\Delta f\) the displacement is along PhaseEncodingDirection, and reversing the blips reverses it. The readout, at a bandwidth per pixel of PixelBandwidth, is displaced by a fraction of a voxel, and you neglect it.

lines = sidecar_up["AcquisitionMatrixPE"]
echo_spacing = sidecar_up["EffectiveEchoSpacing"]
bandwidth_pe = 1 / (lines * echo_spacing)

64 lines at an effective echo spacing of 0.55 ms
bandwidth per pixel: 28.4 Hz along phase encoding (sidecar 28.409 Hz), 2055 Hz along the readout
total readout time 34.65 ms
100 Hz off resonance: 3.5 voxels = 13.2 mm along phase encoding

The bandwidth you compute agrees with the one the sidecar states. Along the readout the bandwidth per pixel is about seventy times larger, which is why only phase encoding is distorted.

The reference field map

To check the correction you need an independent measurement of the field. The gradient-echo field map is the phase difference \(\Delta\phi\) between two echoes, stored by the scanner as integers from 0 to 4095 for \(-\pi\) to \(\pi\), and the off-resonance is \(\Delta f = \Delta\phi / (2\pi\, \Delta\mathrm{TE})\). You smooth the phase over one voxel (as a complex exponential, hidden below). The displacement it predicts for the blip-up image is \(\Delta f / \mathrm{BW}_{PE}\) voxels towards anterior.

phase_file = CACHE / "sub-1/fmap/sub-1_acq-v4_phasediff"
magnitude_file = CACHE / "sub-1/fmap/sub-1_acq-v4_magnitude1"
sidecar_phase = json.loads(Path(f"{phase_file}.json").read_text())
stored = read_nifti(f"{phase_file}.nii.gz").image[0]
phase = (stored - 2048) / 2048 * torch.pi
magnitude = read_nifti(f"{magnitude_file}.nii.gz").image[0]
delta_te = sidecar_phase["EchoTime2"] - sidecar_phase["EchoTime1"]


off_resonance = phase / (2 * torch.pi * delta_te)
predicted_mm = off_resonance / bandwidth_pe * voxel_mm[1]

echo times 4.00 and 6.46 ms: unambiguous within +/-203 Hz

Reversed-gradient correction

PyHySCO looks for the displacement field \(b\) for which the blip-up image sampled at \(y + b\) and the blip-down image sampled at \(y - b\), each multiplied by the Jacobian determinant of its transformation, \(1 \pm \partial b / \partial y\), agree. The Jacobian factor restores the intensity of voxels compressed into a pile-up or stretched over several voxels. A smoothness penalty on \(b\), and a constraint that keeps both transformations invertible, regularise the problem.

The estimation is three-dimensional: you pass the volumes with their voxel size and the phase-encoding axis y. You get \(b\) back in millimetres along the voxel axis y, that is towards anterior, on the faces between voxels, and average neighbouring faces to put it at the voxel centres.

result = bt.correct_susceptibility(
    up,
    down,
    voxel_size=voxel_mm,
    phase_encoding_axis=1,
)
faces = result.field_map
estimated_mm = (0.5 * (faces[:, 1:] + faces[:, :-1])).float()
corrected_up = result.blip_up.float()
corrected_down = result.blip_down.float()

The corrected pair

If the correction is right, the two corrected images are the same image. You measure their agreement over the brain, where the corrected EPI and the field map’s magnitude both have signal, as the correlation coefficient and as the root-mean-square difference relative to the mean image.

• 
• 

blip-up vs blip-down, acquired : correlation 0.871, RMS difference 21 % of the mean image
blip-up vs blip-down, corrected: correlation 0.989, RMS difference 6 % of the mean image

Slice 21 passes through the orbitofrontal cortex and slice 14 through the temporal poles and the cerebellum. The lower row of each figure enlarges the anterior part of the slice, in the box, with its own grey scale.

Above the frontal sinus the off-resonance is positive. The blip-up image, phase-encoded from posterior to anterior, displaces the orbitofrontal cortex anteriorly and compresses its anterior edge into a bright rim; the blip-down image displaces it posteriorly and flattens it. At the temporal poles, above the petrous bone and the mastoid air cells, the off-resonance is negative and the directions swap: the blip-down image stretches the poles anteriorly into streaks, and the blip-up image compresses them into a few dark voxels. The corrected image, the mean of the two corrected images, puts the cortex between the two and restores its intensity.

In these overlays the blip-up image is magenta and the blip-down image green, so tissue where the two agree is grey. Before correction the frontal, temporal and occipital edges carry a magenta fringe on one side and a green fringe on the other; after correction they are grey. What is left is green at the temporal poles, where the blip-up image has lost the signal that the blip-down image still has: signal lost within a voxel in one acquisition cannot be brought back by a displacement.

Comparison with the field map

Now you compare the estimate with the displacement the gradient-echo field map predicts. You do it over the whole brain, and over the brain voxels where the field map predicts a displacement of more than one voxel, \(|\Delta f| > \mathrm{BW}_{PE}\): that is where there is a distortion to measure. Over the rest of the brain the field is close to zero, and a correlation there measures mostly noise. The field map was acquired with its own shim, so the two fields may also differ by a smooth field of first and second order.

distorted = brain & (off_resonance.abs() > bandwidth_pe)

• 
• 

shim, EPI       [1746, -3889, -287, 575, 111, -87, -1, -28]
shim, field map [1747, -3883, -285, 583, 125, -23, -45, -40]
field over the brain: -45 to +47 Hz (1st to 99th percentile)
displaced by more than one voxel: 7 % of the brain voxels
estimate vs field map over the brain: correlation 0.38, slope 0.31
estimate vs field map over the displaced voxels: correlation 0.61, slope 0.30

The estimate follows the pattern of the field-map prediction: anterior displacement above the frontal sinus and in the cerebellum, posterior displacement at the temporal poles. It is smaller in amplitude: the fitted slope is about 0.3, and in the joint histograms the voxels lie flatter than the diagonal. Two likely causes, neither established here, are the smoothness penalty on the displacement field, which at 3.75 mm voxels flattens a displacement that changes within a few voxels, and the difference between the shims of the two acquisitions. Neither EPI image has signal in the voids above the sinus, where the field map predicts the largest displacement, so the estimate there is an extrapolation by the smoothness penalty.

In practice you acquire the pair as two short series of a few volumes, estimate the displacement field once, and apply it to every volume of the functional or diffusion series acquired with one of the two polarities.

As a spec

What this tour built, stated the way you would ask an agent for it:

Correct the susceptibility distortion of a reversed phase-encoding EPI
pair (BIDS, 3 T, phase encoding along y) with bartorch. Read both
series with bartorch.io.read_nifti and average their volumes; compute
the bandwidth per pixel along phase encoding as
1 / (AcquisitionMatrixPE * EffectiveEchoSpacing) from the sidecar.
Run bt.correct_susceptibility(up, down, voxel_size=(z, y, x) in mm,
phase_encoding_axis=1) and move the returned displacement (mm, on the
faces between voxels) to the voxel centres. Report the correlation and
RMS difference of the two images over the brain before and after, and
compare the displacement with the one a dual-echo GRE field map
predicts (correlation and slope).

References

[1]

Andersson JLR, Skare S, Ashburner J. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage 20(2):870-888 (2003). https://doi.org/10.1016/S1053-8119(03)00336-7

[2]

Julian A, Ruthotto L. PyHySCO: GPU-enabled susceptibility artifact distortion correction in seconds. Front Neurosci (2024).

[3]

Cieslak M, Elliott M, Satterthwaite T. Example Fieldmaps. OpenNeuro, accession ds001600. https://openneuro.org/datasets/ds001600. Licensed under CC-BY-SA.