Off-resonance correction of spiral imaging

In this tour you correct the blurring that \(B_0\) off-resonance causes in a spiral image. A spiral readout reaches each k-space radius at its own time, so a spin off resonance by \(f\) picks up a phase \(2\pi f\, t(k)\) that grows from the centre of k-space outwards. In a Cartesian scan the same phase is linear along the readout and only shifts the spin; on a spiral it is radial, and it spreads the spin into a ring. With a readout of tens of milliseconds, the field near the sinuses and the ear canals is enough to smear the orbitofrontal and temporal cortex over several voxels.

You simulate an axial spiral scan of the head at 3 T, in a field computed from the susceptibility of the head, and you correct it with the known field map in two ways: by multifrequency interpolation (MFI) on the gridded image, which is what Gadgetron’s spiral deblurring gadget does, and by a model-based reconstruction whose encoding operator includes the field map by time segmentation.

Prerequisites: lessons 7 and 8, 7. Basics of non-Cartesian imaging and 8. Radial SENSE.

Learning objectives

• Relate the blurring of a spiral image to the off-resonance frequency and the readout duration.

• Describe a spiral readout by its time map \(t(|k|)\) and factorise its off-resonance transfer with fit_transfer().

• Deblur the gridded image with deblur(), and choose the number of demodulation frequencies.

• Set up a time-segmented reconstruction with FieldCorrected(), and see where it does better than conjugate-phase methods such as MFI.

The object and its field map

Your object is an axial slice of the BrainWeb T1-weighted head [1] through the orbitofrontal cortex and the temporal lobes, 128 x 128 over a 220 mm field of view (1.7 mm in-plane). Its field map is the \(B_0\) offset at 3 T that the susceptibility difference between air and tissue, \(\Delta\chi = 9.4\) ppm, produces: computed in 3D with the dipole kernel, less a second-order shim fitted over the brain, rounded to 1 Hz, limited to \(\pm 150\) Hz, and set to zero outside the head, where a measured field map would have no signal to be estimated from. How both are built is not the subject here, so that code is hidden.

import math

import torch

import bartorch
import bartorch.tools as bt

from bartorch import linop, optim

SIZE, FOV_MM = 128, 220.0

field over the head: -109 to +150 Hz
field over the brain: -41 to +44 Hz

The largest offsets are in the scalp. Inside the brain, outlined in grey, the field reaches about \(-40\) Hz in the lateral temporal lobes and \(+40\) Hz in the orbitofrontal cortex. Below you will see that this is about one cycle of phase over the readout.

The spiral readout

You acquire four interleaves of an Archimedean spiral that reach \(k_{max}\) at the Nyquist edge of the 128 matrix, each in a 24 ms readout of 6000 samples (4 µs dwell time). The interleaves are rotations of one arm by \(2\pi/4\), and the arm makes 16 turns, so the rings of the four interleaves together are one grid unit apart.

Near the centre of k-space the slew rate limits how fast the spiral can turn, and further out the gradient amplitude limits how fast it can travel. The readout time \(t/T = (s^2 + 0.2\, s) / 1.2\) along the arm coordinate \(s = |k|/k_{max}\) models both. The trajectory is in grid units.

INTERLEAVES, SAMPLES, READOUT_S = 4, 6000, 24e-3
TURNS = SIZE / 2 / INTERLEAVES

sample_time = torch.linspace(0.0, READOUT_S, SAMPLES, dtype=torch.float64)
arm = torch.sqrt(0.01 + 1.2 * sample_time / READOUT_S) - 0.1  # s(t)
rotations = torch.arange(INTERLEAVES)[:, None] / INTERLEAVES
angle = 2 * math.pi * (TURNS * arm + rotations)
radius = SIZE / 2 * arm
trajectory = torch.stack(
    [
        radius * torch.cos(angle),
        radius * torch.sin(angle),
        torch.zeros_like(angle),
    ],
    dim=-1,
).float()
sample_time = sample_time.float()

The acquisition

A voxel at off-resonance frequency \(f\) contributes \(x(r)\, e^{-2\pi i k \cdot r}\, e^{2\pi i f\, t(k)}\) to the sample at \(k\). You can simulate that exactly for this field map, because it takes only integer frequencies: transform each frequency’s part of the object with bartorch.nufft(), give it the phase it accrues at each sample time, and add up the parts. Your reference is the same acquisition on resonance.

def acquire(field):
    samples = torch.zeros(INTERLEAVES, SAMPLES, 1, dtype=torch.complex64)
    for frequency in torch.unique(field[head]):
        part = image * ((field == frequency) & head)
        phase = 2 * math.pi * float(frequency) * sample_time
        accrued = torch.polar(torch.ones(SAMPLES), phase)
        samples += bartorch.nufft(part, trajectory) * accrued[:, None]
    return samples


on_resonance = bartorch.nufft(image * head, trajectory)
off_resonance = acquire(field_map)

Gridding

You reconstruct both by the density-compensated adjoint NUFFT, with Pipe-Menon weights from bartorch.estimate_density(). A voxel off resonance by \(f\) is spread over a ring whose size grows with the phase \(2\pi f T\) it has accrued by the end of the readout: one full cycle at \(f = 1/T \approx 42\) Hz.

kxy = trajectory[..., :2].reshape(-1, 2)
density = bartorch.estimate_density(kxy, (SIZE, SIZE))
density = density.reshape(INTERLEAVES, SAMPLES, 1)


def grid(samples):
    return bartorch.nufft_adjoint(
        samples * density,
        trajectory,
        image_shape=(SIZE, SIZE),
    )


reference = grid(on_resonance)
blurred = grid(off_resonance)

Off resonance, the scalp, where the field is largest, is spread into a halo that overlaps the frontal and temporal cortex, and the cortex of the temporal lobes is smeared over a few voxels. The two boxes mark the regions you will enlarge below.

Multifrequency interpolation

Conjugate-phase reconstruction [2] undoes the phase voxel by voxel,

\[\hat x(r) = \sum_k w_k\, y_k\, e^{2\pi i k \cdot r}\, e^{-2\pi i f(r)\, t_k},\]

which costs one transform per voxel. MFI [3] costs one transform per demodulation frequency instead. It approximates the transfer over the readout as \(e^{-2\pi i f t} \approx \sum_m a_m(f)\, e^{-2\pi i f_m t}\), demodulates the gridded image at each \(f_m\) in k-space, and combines the demodulated images voxel by voxel with the weights \(a_m(f(r))\).

Since every interleaf is a rotation of the first, the readout time depends on \(|k|\) alone. ReadoutTiming tabulates it from one interleaf. fit_transfer() then places the demodulation frequencies uniformly over the band of the field map and, unless you give it a number, takes the fewest that approximate the transfer to 1 % RMS, starting from \(\lceil 2.5\, f_{max}\, T \rceil\), the number Gadgetron’s MFIOperator uses.

timing = bt.ReadoutTiming.from_trajectory(
    trajectory[0, :, :2],
    duration=READOUT_S,
)
band = float(field_map[head].abs().max())
transfer = bt.fit_transfer(timing, band=band)
deblurred = bt.deblur(blurred, field_map, transfer)

band +-150 Hz: 12 demodulation frequencies, RMS error of the transfer 6.2e-03

Time-segmented reconstruction

The other way is to put the field map into the encoding operator,

\[y = \sum_l \operatorname{diag}(b_l)\, E\, \operatorname{diag}(c_l)\, x,\]

with \(E\) the NUFFT, and temporal and spatial coefficients \(b_l(t)\) and \(c_l(r)\) fitted to \(e^{2\pi i f(r) t}\) over the histogram of the field map [4]. Then you solve that model like any other. FieldCorrected() builds the operator from the field map and the sample times, and CG solves the normal equations. Unlike a conjugate-phase method, this does not assume the field to be constant over the size of the blurring.

As references, you also run the same CG without the field map, on the data off and on resonance.

coils = torch.ones(1, SIZE, SIZE, dtype=torch.complex64)
E = linop.NoncartesianSense(coils, (SIZE, SIZE), traj=trajectory)
A = linop.FieldCorrected(
    E,
    field_map,
    readout_time=sample_time,
    mask=head,
    segments=transfer.terms,
)
solve = optim.CG(maxiter=20)
model_based = solve(off_resonance[None, ..., 0], A)
model_uncorrected = solve(off_resonance[None, ..., 0], E)
model_reference = solve(on_resonance[None, ..., 0], E)

Results

The spiral samples a disc of radius \(k_{max}\), and truncating the object’s spectrum at its edge leaves ringing around the scalp in every reconstruction, on resonance too. So the best you can get with this readout is the on-resonance reconstruction of the same kind, and that is what you compare each image with, as the normalised root-mean-square error (NRMSE) over the brain. The two on-resonance reconstructions themselves are compared with the object after a least-squares scaling.

• 
• 
• 
• 

gridding, on resonance  0.040 against the object
CG, on resonance        0.016 against the object
gridding, uncorrected   0.079
gridding, MFI           0.046
CG, uncorrected         0.088
CG, time-segmented      0.028

Both corrections restore the position of the cortex, in the overview and in the two enlarged regions, and both bring the NRMSE over the brain well below that of their uncorrected counterpart. What MFI leaves, in the last figure, sits mostly in the scalp, where the field changes by tens of hertz across the ring a voxel is blurred into. Like every conjugate-phase method, MFI assumes the field is constant over that ring, and where it is not, the demodulation leaves an intensity error. The time-segmented reconstruction models the phase of each voxel up to the segmentation error and removes most of that residual too. The cost is an iterative solve, one NUFFT pair per segment and iteration, against one FFT pair per demodulation frequency for MFI.

Exercise: demodulation frequencies

The MFI approximation is poor until the demodulation frequencies are about \(1/T\) apart. Ask fit_transfer() for a fixed number of them and watch three things: the RMS error of the transfer, the largest \(\sum_m |a_m(f)|\) (amplification, the factor by which noise and residual error are amplified), and the NRMSE of the corrected image over the brain.

rows = []
for terms in (9, 13, 17, 21, 25):
    trial = bt.fit_transfer(timing, band=band, terms=terms)
    corrected = bt.deblur(blurred, field_map, trial)
    rows.append((terms, trial, corrected))

terms  transfer  sum|a|  brain NRMSE
    9   1.6e-01     2.3        0.177
   13   1.9e-03    14.8        0.046
   17   8.7e-06   191.6        0.046
   21   1.7e-08  2853.1        0.046
   25   3.0e-10  1777.4        0.046

With too few frequencies the transfer is poorly approximated and the image is off. Once the transfer error is small, more frequencies leave the image as it is, and the amplification grows by orders of magnitude. Gadgetron’s gpuSpiralDeblurGadget sets the band from the echo spacing of its field map, \(f_{max} = 1.2 / (2\,\Delta TE)\), estimates the field map from a low-pass filtered two-echo spiral, and applies MFI with \(\lceil 2.5\, f_{max}\, T \rceil\) frequencies. Here you knew the field map exactly; in practice its own error and smoothing limit both corrections.

As a spec

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

Correct the off-resonance blurring of a 2D spiral acquisition (4
interleaves of 6000 samples, 24 ms readout, 128 matrix) with a known
field map in Hz, with bartorch. Grid with the density-compensated
adjoint NUFFT (bartorch.estimate_density weights), then deblur by
multifrequency interpolation: bt.ReadoutTiming.from_trajectory on one
interleaf, bt.fit_transfer over the band of the field map, bt.deblur.
Also reconstruct with linop.FieldCorrected around a single-coil
linop.NoncartesianSense (as many segments as MFI frequencies, masked
to the head) and optim.CG, 20 iterations. Report the brain NRMSE of
each against the on-resonance reconstruction of the same kind, and
tabulate the MFI error for 9 to 25 frequencies.

References

[1]

Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imaging 17(3):463-468 (1998). https://doi.org/10.1109/42.712135

[2]

Noll DC, Meyer CH, Pauly JM, Nishimura DG, Macovski A. A homogeneity correction method for magnetic resonance imaging with time-varying gradients. IEEE Trans Med Imaging 10(4):629-637 (1991). https://doi.org/10.1109/42.108599

[3]

Man LC, Pauly JM, Macovski A. Multifrequency interpolation for fast off-resonance correction. Magn Reson Med 37(5):785-792 (1997). https://doi.org/10.1002/mrm.1910370523

[4]

Sutton BP, Noll DC, Fessler JA. Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities. IEEE Trans Med Imaging 22(2):178-188 (2003). https://doi.org/10.1109/TMI.2002.808360