Echo planar imaging

A gradient echo acquires one k-space line per excitation. Echo planar imaging (EPI) acquires the whole matrix after a single excitation. This Tour first follows the readout gradient with further readout gradients of alternating polarity, each forming an echo, and then adds a phase-encode blip between the echoes so that each echo acquires a different line.

You then measure the two effects that limit such an acquisition: the signal decay over an echo train tens of milliseconds long, which blurs the image along the phase-encode direction, and the sensitivity of a train of alternating readouts to a delay between the gradient waveform and the acquisition, which produces the half-field-of-view ghost.

Prerequisites: lessons 1 and 2 of the course.

Readout train

A readout gradient traverses one k-space line from one end to the other. A second gradient of the opposite polarity traverses it back, so a train of them needs no rewinder between the echoes and forms one echo per gradient. The gradients are played back to back, so the echo spacing is the duration of one readout gradient, ramps included. Every second echo is acquired in the opposite direction, and its samples are in the reverse order of the odd echoes’.

import numpy as np

import pypulseqpp as pp

system = pp.Opts(
    max_grad=32.0,
    grad_unit="mT/m",
    max_slew=130.0,
    slew_unit="T/m/s",
    rf_dead_time=100e-6,
    rf_ringdown_time=20e-6,
    adc_dead_time=10e-6,
)

FOV = 220e-3
MATRIX = 64
THICKNESS = 5e-3
FLIP_ANGLE_DEG = 90.0
DWELL = 4e-6

rf, gz, gz_reph = pp.make_sinc_pulse(
    flip_angle=np.deg2rad(FLIP_ANGLE_DEG),
    duration=2e-3,
    slice_thickness=THICKNESS,
    apodization=0.5,
    time_bw_product=4.0,
    delay=system.rf_dead_time,
    system=system,
    use="excitation",
    return_gz=True,
)

acquisition = MATRIX * DWELL
raster = system.grad_raster_time
gx = pp.make_trapezoid(
    channel="x",
    amplitude=MATRIX / FOV / acquisition,
    flat_time=raster * np.ceil(acquisition / raster),
    system=system,
)
adc = pp.make_adc(num_samples=MATRIX, dwell=DWELL, delay=gx.rise_time, system=system)
gx_pre = pp.make_trapezoid(
    channel="x",
    area=-(gx.amplitude * gx.rise_time / 2 + (MATRIX / 2 + 0.5) / FOV),
    duration=5e-4,
    system=system,
)

ECHOES = 6

train = pp.Sequence(system=system)
train.add_block(rf, gz)
train.add_block(gx_pre, gz_reph)
for echo in range(ECHOES):
    train.add_block(pp.scale_grad(gx, (-1.0) ** echo), adc)

The analysis gives the k-space location of every sample of every acquisition window. Along the readout axis the train is a triangle wave between the two ends of the line, and an echo is where it crosses zero.

k_train, _, t_train_excitation, _, t_train = train.calculate_kspace()
kx_train = k_train[0] * FOV / MATRIX * 2
echo_times = np.array(
    [
        t_train[
            echo * MATRIX + int(np.argmin(np.abs(kx_train[echo * MATRIX :][:MATRIX])))
        ]
        - t_train_excitation[0]
        for echo in range(ECHOES)
    ]
)
print(
    "echo times (ms): "
    + ", ".join(f"{1e3 * time:.2f}" for time in echo_times)
    + f"\nspacing {1e3 * np.diff(echo_times).mean():.3f} ms, "
    f"readout gradient {1e3 * pp.calc_duration(gx):.3f} ms"
)

echo times (ms): 1.99, 2.69, 3.39, 4.09, 4.79, 5.49
spacing 0.700 ms, readout gradient 0.700 ms

The samples the analysis reports are not one line acquired six times: the even echoes run from \(+k_\mathrm{max}\) to \(-k_\mathrm{max}\), so a reconstruction has to reverse them before they are lines of the same matrix. The last section measures what a delay between the gradient and the acquisition does to such a train.

Single-shot acquisition

A phase-encode blip between the echoes advances the acquisition by one line, and the train runs the length of the matrix. The prewinder on the phase-encode axis starts the train at one edge of k-space.

gy_pre = pp.make_trapezoid(
    channel="y", area=-MATRIX / (2 * FOV), duration=5e-4, system=system
)
blip = pp.make_phase_blip(channel="y", fov=FOV, steps=1, system=system)

seq = pp.Sequence(system=system)
seq.add_block(rf, gz)
seq.add_block(gx_pre, gy_pre, gz_reph)
for echo in range(MATRIX):
    readout = pp.scale_grad(gx, (-1.0) ** echo)
    if echo == 0:
        seq.add_block(readout, adc)
    else:
        seq.add_block(readout, blip, adc)

ok, errors = seq.check_timing()
echo_spacing = pp.calc_duration(gx)
print(
    f"timing {ok}, {seq.num_blocks} blocks, "
    f"echo spacing {1e3 * echo_spacing:.3f} ms, "
    f"train {1e3 * MATRIX * echo_spacing:.1f} ms, "
    f"blip {1e6 * pp.calc_duration(blip):.0f} us"
)

seq.paper_plot()

timing True, 66 blocks, echo spacing 0.700 ms, train 44.8 ms, blip 80 us

k-space trajectory

One shot, coloured by the rank of each echo in the train: the acquisition starts at one corner of k-space and works across it, reversing direction at every line.

pp.plot.plot_kspace(seq, color_by="echo", plane="xy")

Signal decay along the train

The last line is acquired tens of milliseconds after the first, and the transverse signal has decayed over that interval. Along the phase-encode direction the acquired data are therefore the true k-space multiplied by a decaying envelope, and the image is convolved with that envelope’s transform.

The echo times come from the analysis, and the line each echo lands on from the k-space it reports, so the envelope follows the sequence’s own ordering.

k_adc, _, t_excitation, _, t_adc = seq.calculate_kspace()
echo_time = t_adc.reshape(MATRIX, MATRIX).mean(axis=1) - t_excitation[0]
line = np.round(k_adc[1].reshape(MATRIX, MATRIX)[:, 0] * FOV).astype(int)
order = np.argsort(line)

T2_STARS = (20e-3, 40e-3, 60e-3, 100e-3)


def point_spread(t2_star):
    """The phase-encode point-spread function the decay produces."""
    envelope = np.exp(-echo_time[order] / t2_star)
    spread = np.abs(np.fft.fftshift(np.fft.fft(np.fft.ifftshift(envelope))))
    return spread / spread.max()


def width(spread):
    """Full width at half maximum, in pixels, by linear interpolation."""
    above = np.flatnonzero(spread >= 0.5)
    first, last = above[0], above[-1]
    left = np.interp(0.5, spread[first - 1 : first + 1], [first - 1, first])
    right = np.interp(0.5, spread[last : last + 2][::-1], [last + 1, last])
    return right - left


widths = {t2_star: width(point_spread(t2_star)) for t2_star in T2_STARS}

   T2*    decay over the train    point spread
 20 ms                   0.100         1.51 px
 40 ms                   0.316         1.21 px
 60 ms                   0.464         1.13 px
100 ms                   0.631         1.08 px

The envelope is not centred on the middle of k-space: the train runs from one edge to the other, so the decay is monotonic across the lines rather than symmetric about the line the echo is on. The result is a point-spread function that is both wider than one pixel and asymmetric: the blurring along the phase-encode direction of an echo planar image. At the shortest \(T_2^*\) here the signal at the last line is a tenth of the first and the point spread is half again as wide as a pixel; at the longest it is within a tenth of a pixel of the unblurred width.

The remedies are a shorter echo spacing, fewer lines per shot (a segmented acquisition), and acquiring fewer lines outright, by partial Fourier or by parallel imaging.

Gradient delay and odd-echo shift

A delay between the gradient waveform and the acquisition displaces every sample along the readout direction by the distance k-space travels in that delay. The readouts alternate in polarity, so the displacement alternates in sign: the odd and the even lines of the matrix are shifted in opposite directions.

calculate_kspace() takes the delay as a parameter, so the displacement is measured from the trajectory the analysis reports rather than computed beside it.

DELAYS = np.array([0.0, 1e-6, 2e-6, 4e-6, 8e-6])

displacement = []
for delay in DELAYS:
    delayed = seq.calculate_kspace(trajectory_delay=delay)[0]
    kx = delayed[0].reshape(MATRIX, MATRIX)
    # The echo of each line, where the ideal trajectory crosses zero.
    centre = kx[:, MATRIX // 2] * FOV
    displacement.append(
        {
            "delay": delay,
            "samples": float(np.abs(centre[::2].mean() - centre[1::2].mean()) / 2),
        }
    )

 delay     odd-even displacement
0.0 us                 0.000 samples
1.0 us                 0.250 samples
2.0 us                 0.500 samples
4.0 us                 1.000 samples
8.0 us                 2.000 samples

The displacement is the delay divided by the dwell time, and it is the same for every line, so the odd and the even lines differ only in its sign. A quantity that alternates from one line to the next along the phase-encode direction is, after the transform, an image displaced by half the field of view, which is the ghost an uncorrected echo planar acquisition shows. Measuring the delay and correcting for it are reconstruction steps, outside the scope of this package.

Specification

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

Write a pypulseqpp single-shot EPI for a scanner with 32 mT/m and
130 T/m/s, 100 us RF dead time, 20 us ringdown and 10 us ADC dead
time. Protocol: 220 mm FOV, 64x64 matrix, 5 mm slice, 90 degree
sinc excitation (2 ms, time-bandwidth product 4), 4 us dwell. First
play a train of 6 readouts of alternating polarity and locate the
echoes with calculate_kspace. Then play 64 readouts with a one-line
phase-encode blip between them, check the timing, and plot the
k-space trajectory coloured by echo. Compute the phase-encode
point-spread function for T2* of 20, 40, 60 and 100 ms from the
echo times, and the odd-even k-space displacement for gradient
delays of 0 to 8 us with calculate_kspace(trajectory_delay=...).