EPI Nyquist ghost and ramp sampling

Note

Go to the end to download the full example code.

EPI Nyquist ghost and ramp sampling#

Open in Colab

An echo-planar readout acquires k-space in a train of lines of alternating readout gradient polarity, and two corrections are applied before the lines form a Cartesian k-space. A timing error between the readout gradient and the ADC, and eddy currents, displace the echoes of the reversed lines relative to the forward ones; the resulting odd/even phase produces a Nyquist ghost, a copy of the object displaced by half the field of view along the phase-encoding direction. Sampling during the ramps of the readout gradient shortens the echo spacing, and the samples, uniform in time, are not uniform in \(k_x\).

This example simulates an eight-channel EPI acquisition of a Shepp-Logan phantom with a known ADC delay and constant phase error, estimates the odd/even phase from a three-line navigator, and compares the ghosted and the corrected image with the delay-free image. It then resamples a ramp-sampled readout onto a uniform \(k_x\) grid and compares the result with linear interpolation.

Learning objectives

  • Relate an ADC delay to a linear phase in hybrid space and to the Nyquist ghost at half the field of view.

  • Estimate the odd/even phase from a three-line navigator with estimate_epi_phase(), and apply it to the reversed lines with correct_lines().

  • Quantify the ghost by the ghost-to-signal ratio.

  • Resample a ramp-sampled readout with epi_ramp_operator(), and check the sampling condition under which the resampling is exact.

import math

import numpy as np
import torch

import bartorch
import bartorch.tools as bt

The Nyquist ghost#

A delay \(\delta\) of the ADC relative to the readout gradient shifts the echo of every line by the same time: towards positive \(k_x\) on a forward line and towards negative \(k_x\) on a reversed one. In hybrid space – after the inverse Fourier transform along the readout – a shift of \(\delta\) dwell times is a linear phase \(\pi \delta u\) over the readout coordinate \(u \in [-1, 1]\). A constant phase \(\phi_0\), from eddy currents or a \(B_0\) offset during the readout, adds to it with the same alternating sign. The phase difference between odd and even lines modulates k-space at the Nyquist frequency of the phase-encoding direction, and the image acquires a ghost at half the field of view.

Forward lines carry \(+(\pi \delta u + \phi_0)\) and reversed lines \(-(\pi \delta u + \phi_0)\). A reversed line is stored in the order it was digitised, that is flipped along the readout.

SIZE = 128
DELAY = 0.2  # dwell times
PHASE_0 = 0.25  # rad
SLOPE = math.pi * DELAY * (SIZE - 1) / SIZE  # rad over u in [-1, 1]

kspace = bt.phantom(SIZE, kspace=True, coils=8)[:, 0]  # (coils, ky, kx)
hybrid = bartorch.fft(kspace, axes=(-1,), inverse=True, unitary=True)
u = torch.linspace(-1.0, 1.0, SIZE)


def acquire(row, polarity):
    """One readout of polarity +1 or -1, digitised in the order it was played."""
    delayed = row * torch.polar(torch.ones(SIZE), polarity * (SLOPE * u + PHASE_0))
    line = bartorch.fft(delayed, axes=(-1,), unitary=True)
    return torch.flip(line, [-1]) if polarity < 0 else line


train = [(acquire(hybrid[:, ky], 1 - 2 * (ky % 2)), ky % 2 == 1) for ky in range(SIZE)]

The navigator#

Three lines of alternating polarity acquired without phase-encoding blips sample the same line of k-space, so the phase between the reversed line and the mean of its two neighbours is the odd/even phase alone. estimate_epi_phase() fits a polynomial to it, weighted by the signal magnitude and summed over coils. The reversed navigator line is passed already flipped into readout order.

centre = hybrid[:, SIZE // 2]
navigator = [acquire(centre, 1), torch.flip(acquire(centre, -1), [-1]), acquire(centre, 1)]

fit = bt.estimate_epi_phase(navigator)
print(f"fitted     constant {float(fit[0]):+.3f}  linear {float(fit[1]):+.3f} rad")
print(f"impressed  constant {2 * PHASE_0:+.3f}  linear {2 * SLOPE:+.3f} rad")

fitted     constant +0.500  linear +1.247 rad
impressed  constant +0.500  linear +1.247 rad

The fitted phase is twice the impressed one, because the navigator measures the difference between a forward and a reversed line. The correction is applied to the reversed lines only: it brings them into phase with the forward lines, whose remaining phase is common to all lines and does not change the magnitude image.

The ghost is measured by the ghost-to-signal ratio: the mean signal outside the object divided by the mean signal inside it.

def reconstruct(lines):
    coil_images = bartorch.fft(torch.stack(lines, dim=1), axes=(-2, -1), inverse=True, unitary=True)
    return bartorch.rss(coil_images, axes=(0,)).abs()


ideal = reconstruct([kspace[:, ky] for ky in range(SIZE)])
ghosted = reconstruct(bt.correct_lines(train))
corrected = reconstruct(bt.correct_lines(train, fit))

inside = ideal > 0.05 * ideal.max()
outside = ideal < 0.01 * ideal.max()


def nrmse(estimate, target):
    return float((estimate - target).norm() / target.norm())


def ghost_to_signal(image):
    return float(image[outside].mean() / image[inside].mean())


for name, image in (("delay-free", ideal), ("flipped only", ghosted), ("corrected", corrected)):
    print(
        f"{name:>12}: ghost-to-signal {100 * ghost_to_signal(image):5.2f} %, "
        f"NRMSE {nrmse(image, ideal):.1e}"
    )

  delay-free: ghost-to-signal  1.38 %, NRMSE 0.0e+00
flipped only: ghost-to-signal 18.22 %, NRMSE 3.2e-01
   corrected: ghost-to-signal  1.38 %, NRMSE 1.4e-07
  • delay-free, flipped only, corrected
  • navigator phase, central column

The images are displayed from zero to 30 % of the peak, with the phase-encoding direction vertical. Without the phase correction the ghost appears at the top and bottom of the field of view and overlaps the object where it wraps. With the first-order fit the ghost-to-signal ratio returns to that of the delay-free image, whose signal outside the object is the truncation ringing of the phantom, since the simulated phase error is exactly first order. On measured data, higher orders of the phase, and phase errors that differ between lines of the same polarity, leave a residual ghost.

Ramp sampling#

Sampling during the ramps of the trapezoidal readout gradient places the samples at \(k_x(t) = \gamma \int_0^t G_x(\tau)\, d\tau\), which is denser on the ramps than on the plateau. The readout is the Fourier transform of an object that spans support pixels, so samples at any positions determine it, provided no two neighbouring samples are further apart than one over the support, the Nyquist spacing of that object. epi_ramp_operator() [1] is the regularized least-squares inverse of the transform at the sampled positions followed by the transform at the uniform ones.

The readout gradient is a trapezoid whose ramps each take 30 % of the ADC window, sampled with 160 samples, for a one-dimensional object of 64 pixels; positions are in cycles per pixel.

SUPPORT, SAMPLES = 64, 160
profile = torch.zeros(SUPPORT, dtype=torch.complex128)
profile[16:48] = torch.linspace(0.3, 1.0, 32) * torch.exp(1j * torch.linspace(0.0, 2.0, 32))
pixels = (torch.arange(SUPPORT) - SUPPORT // 2).double()

time = torch.linspace(0.0, 1.0, SAMPLES, dtype=torch.float64)
RAMP = 0.3
gradient = torch.clamp(torch.minimum(time / RAMP, (1 - time) / RAMP), 0, 1)
sampled_at = torch.cumsum(gradient, 0)
sampled_at = sampled_at - sampled_at.mean()
sampled_at = 0.5 * sampled_at / sampled_at.abs().max()
uniform_at = torch.arange(SAMPLES, dtype=torch.float64) / SAMPLES - 0.5


def encode(positions):
    return torch.exp(-2j * math.pi * torch.outer(positions, pixels)) @ profile


measured, truth = encode(sampled_at), encode(uniform_at)

operator = bt.epi_ramp_operator(sampled_at, uniform_at, SUPPORT)
resampled = (measured.to(torch.complex64)[None] @ operator.T)[0].to(torch.complex128)

Linear interpolation between neighbouring samples is the comparison. Both are assessed on the image profile, the inverse transform of the uniform samples.

linear = torch.complex(
    torch.from_numpy(np.interp(uniform_at, sampled_at, measured.real)),
    torch.from_numpy(np.interp(uniform_at, sampled_at, measured.imag)),
)


def image_profile(samples):
    return torch.fft.fftshift(torch.fft.ifft(torch.fft.ifftshift(samples))).abs()


step = float(torch.diff(sampled_at).max()) * SUPPORT
print(f"largest step x support: {step:.2f}  (the samples determine the object below 1)")
for name, estimate in (("band-limited resampling", resampled), ("linear interpolation", linear)):
    print(f"{name:>24}: image NRMSE {nrmse(image_profile(estimate), image_profile(truth)):.1e}")

largest step x support: 0.57  (the samples determine the object below 1)
 band-limited resampling: image NRMSE 1.4e-06
    linear interpolation: image NRMSE 3.2e-02
  • trapezoidal readout gradient
  • image profile, |error| / peak

The band-limited resampling reproduces the profile of a uniformly sampled readout to the precision of the operator, which is returned in single precision. Linear interpolation errs most on the plateau, where the samples are furthest apart; in the image its error is spread over the whole field of view, inside and outside the object, at up to a few per cent of the peak. The ratio printed above is the condition to check on a measured trajectory: where a step exceeds one over the support, the readout is aliased and no resampling recovers it.

References#

Total running time of the script: (0 minutes 1.761 seconds)

Gallery generated by Sphinx-Gallery