+1(978)310-4246 credencewriters@gmail.com
  

Magnetic Resonance in Medicine 50:343–353 (2003)
Nonlinear Phase Correction for Navigated
Diffusion Imaging
Karla L. Miller* and John M. Pauly
Motion during diffusion-weighted imaging (DWI) introduces
phase errors that can cause significant artifacts in brain images. One method of correcting these errors uses additional
navigator data to measure the phase corruptions. Standard
navigator methods correct for rigid-body motion but cannot
correct for nonrigid deformations of the brain related to the
cardiac cycle. This work derives a generalized reconstruction
that corrects for nonrigid motion based on a least-squares
formulation. Since this reconstruction has the disadvantage of
being computationally expensive, an approximation is presented, called a refocusing reconstruction. The refocusing reconstruction is both efficient and straightforward. Each readout
is multiplied in image space by the phase conjugate of the
navigator image, and these rephased readouts are then
summed. The conditions under which the refocusing reconstruction is sufficient are considered and methods to improve
the quality of refocused images are discussed. In particular,
synchronization of the acquisition to the cardiac cycle can
provide data that is well-conditioned to the refocusing reconstruction without incurring the large time penalty traditionally
associated with cardiac gating. These methods are applied to
steady-state DWI, a promising pulse sequence that is particularly sensitive to motion-induced phase artifacts. The refocusing reconstruction is shown to significantly improve SS-DWI
over standard rigid-body corrections. Magn Reson Med 50:
343–353, 2003. © 2003 Wiley-Liss, Inc.
Key words: diffusion MRI; motion correction; artifact reduction;
navigator; reconstruction
Diffusion-weighted imaging (DWI) is a powerful tool for
probing tissue microstructure that uses large bipolar gradient pairs to impart phase to moving particles. Diffusing
water spins within a voxel experience phase dispersion,
causing signal attenuation that is dependent on the
amount of diffusion. This attenuation results in diffusion
contrast, which has been shown to be useful in brain
imaging for detecting tissue disease states, as well as revealing structural information. A key problem in DWI is
nondiffusive bulk motion, which causes phase shifts in the
voxel signal (1). When data are gathered in a multishot
acquisition, phase shifts from different excitations can destructively interfere to cause serious image artifacts.
A number of methods have been proposed to deal with
motion artifacts in DWI. One simple way to avoid motioninduced artifacts is to gather the image data in a single
acquisition, in which case the image magnitude is uncor-
Magnetic Resonance Systems Research Laboratory, Department of Electrical
Engineering, Stanford University, Stanford, California.
Grant sponsor: NIH; Grant numbers: R01 CA50948; R01 AR46904.
*Correspondence to: Karla L. Miller, Packard Electrical Engineering Building,
Room 210, Stanford University, Stanford, CA 94305-9510. E-mail:
bison@stanford.edu
Received 29 August 2002; revised 28 March 2003; accepted 2 April 2003.
DOI 10.1002/mrm.10531
Published online in Wiley InterScience (www.interscience.wiley.com).
© 2003 Wiley-Liss, Inc.
rupted. However, these single-shot acquisitions have limited resolution and are sensitive to susceptibility and
eddy-current effects. Another approach is to correct motion-induced phase artifacts before combining the data
from multiple acquisitions. Phase correction imposes no
intrinsic limits on resolution, but does require the acquisition of additional navigator data that resolve the phase
errors in image space. Early works acquired 1D projection
navigators (1–3) or pairs of orthogonal projections (4), and
subsequent studies have extended this approach to a variety of pulse sequences and trajectories (5–11). An important disadvantage of 1D navigation is the inability to correct for phase shifts orthogonal to the direction of navigation. Recent studies have addressed this problem with 2D
navigation (12–14).
Navigated DWI methods typically assume rigid-body
motion, the effects of which have been well characterized
(1,15). However, a significant component of brain motion
is nonrigid (16 –18). During systole, the brain deforms
about the ventricles and displaces inferiorly through the
foramen magnum as though being pulled by the spinal
cord. The more superior and lateral regions of the brain are
nearly motionless, whereas the more inferior and medial
regions tend to experience large displacements. Even during the relative quiescence of diastole, motion induced by
the cardiac cycle follows this nonlinear spatial pattern as
the brain slowly returns to its resting position.
In the presence of diffusion-weighting gradients, these
deformations cause phase corruptions that vary nonlinearly across the object magnetization. The effect of brain
motion on the DWI signal is shown in Fig. 1. The signal
phase and amplitude are correlated with the cardiac cycle
and exhibit nonlinear spatial variation. The loss in signal
amplitude during peak systole is due to intravoxel dephasing caused by intense motion. The highly focal nature of
this dephasing highlights the spatial nonlinearity of phase
corruptions.
These phase errors are problematic for multishot acquisitions because each readout experiences a different image-space phase corruption (␾ j (r) for the jth readout). This
causes phase interference when readouts are combined to
form a high-resolution image. Navigated techniques correct phase errors with low-resolution navigator data that
are acquired along with each high-resolution data frame
(2). Provided the navigator data 1) samples the same phase
error ␾ j (r) as the high-resolution data, and 2) covers kspace sufficiently to resolve ␾ j (r), the navigator can be
used to correct phase errors in the high-resolution data
before the interleaves are combined. Based on the assumption of rigid-body motion, standard navigator methods correct for zeroth- and first-order phase terms (1). In the
presence of the nonrigid motions described above, a rigidbody correction is a first-order approximation to the de-
343
344
Miller and Pauly
function, and can be described conveniently using matrix
formalism. For an arbitrary trajectory that covers k-space
in C interleaves of length R, corresponding to an N â«» N
image, the acquired data are:
d ⫽ GFPm ⫽ RPm
FIG. 1. Illustration of the effects of motion on DW magnetization.
The magnitude (top) and phase (bottom) of a series of low-resolution
axial single-shot DWI images (2D navigators) are shown over a
single cardiac cycle. The magnitude and phase of a representative
voxel are plotted in the middle over four cardiac cycles. The spatial
variation in the image magnitude and phase is nonlinear and
strongly correlated to the cardiac cycle.
sired high-order correction. Rigid-body corrections leave
residual errors that are strongly correlated to the cardiac
cycle (11). Given that a 2D navigator contains high-order
phase information, image reconstruction should be improved by correcting for high-order phase terms.
This work presents a generalized reconstruction method
that corrects for nonrigid motion based on least-squares
estimation. The full least-squares estimate requires prohibitive computation, but suggests an efficient approximation
called the “refocusing reconstruction.” We consider the
conditions for this approximation, and propose efficient
algorithms for the refocusing reconstruction. We also present
a time-efficient approach to cardiac gating that improves the
performance of the refocusing reconstruction. These methods are applied to steady-state DWI (SS-DWI) (19), a potentially powerful DWI method that is particularly sensitive to
motion-induced phase corruptions.
THEORY
In this section, we examine the effect of nonlinear phase
corruption on DWI data and derive a nonlinear phase
correction based on least-squares estimation. Computational issues require that the least-squares reconstruction
be approximated, and an approximation called the “refocusing reconstruction” is presented. This reconstruction
can be performed in either image-space or k-space.
Least-Squares Reconstruction
Phase corruption can be expressed at each point in image
space as a multiplication of the magnetization by a phase
[1]
where d is a CR â«» 1 vector containing the k-space data, m
is a N 2 â«» 1 vector containing the object magnetization in
Cartesian coordinates, and P (CN 2 â«» N 2 ), F (CN 2 â«» CN 2 ),
and G (CR â«» CN 2 ) are matrices respectively representing
image-space phase corruption, discrete Fourier transform,
and resampling from a Cartesian grid onto the k-space
trajectory. The phase corruption matrix P is a stack of
diagonal matrices containing the phase corruption for each
readout. This phase corruption is measured by the navigator phase, which must be interpolated to match the
high-resolution readout (since the n â«» n navigator is of
lower resolution than the high-resolution readout). The
matrix R ⫽ GF describes the basis set used in the k-space
trajectory. Expressions for these vectors and matrices are
given in Appendix A.
The magnetization can be estimated from the data in Eq.
[1] using least-squares estimation. Assuming an accurate
measurement of the phase errors P, the least-squares reconstruction is:
m̂ ⫽ 共P*R*RP兲⫺1 P*R*d
⫽ M ⫺1 P*R*d
[2]
where * denotes the conjugate transpose. It is helpful to
consider this reconstruction as occurring in three steps: 1)
reconstruction R* of the individual data frames in imagespace, which we call the “interleave reconstruction”; 2)
phase correction P*, which we call the “refocusing operator”; and 3) removal of remaining artifacts by M⫺1 ⫽
(P*R*RP)⫺1, which we call the “unmixing operator.”
These steps are diagrammed in Fig. 2 using the notation
described in this section and in Appendix A.
The data are first partially reconstructed using the interleave reconstruction R* ⫽ F*G*, which separately resamples each readout onto a Cartesian grid at the highresolution readout size, and inverse Fourier transforms the
resampled data without combining the interleaves. The
resulting vector a ⫽ R*d contains the image-space data
corresponding to each readout of k-space data. Unless
each readout covers a full field of view (FOV), the images
in the vector a are aliased.
In the next step, the refocusing operator P* is applied to
the interleave images to yield the refocused image m̃ ⫽ P*a.
The refocusing operator performs the essential task of
rephasing the interleaves based on the phase corruption
measured by the navigator. Each interleave image is multiplied by the phase conjugate of the navigator image to
rephase the unaliased component. The sum of these
rephased interleaves is the refocused image. Although the
unaliased components will add coherently in this summation, the aliased components will not necessarily cancel,
since the orthogonality of the Fourier encoding has been
disturbed. The refocused image is thus the desired image
plus some aliased artifact. To the extent that the phase
Nonlinear Phase Correction for DWI
345
FIG. 2. Diagram of the least-squares matrix reconstruction given in Eq. [2]. Each raw data interleave dj is gridded and inverse Fouriertransformed by the interleave reconstruction to form an aliased image aj. These images are phase-corrected and summed by the refocusing
operator to form the refocused image m̃. Remaining aliased energy is removed by the unmixing operator M⫺1 to give the least-squares
optimal image m̂. If the aliased energy in the refocused image is sufficiently low, the reconstruction can stop at m̃, avoiding the
computationally expensive inverse calculation M⫺1. In this case, we call the algorithm a “refocusing reconstruction.”
errors are random, this aliased energy will be incoherent
and will tend to phase-cancel. The characterization and
minimization of the aliased signal are discussed later.
The final step of the least-squares reconstruction is the
application of the unmixing operator M⫺1 to the refocused
image. The unmixing operator removes two artifacts from
the refocused image: nonuniformities in the sampling trajectory, and motion-induced aliasing (see Appendix B).
The compensation for trajectory nonuniformities, a standard correction in the least-squares reconstruction of nonCartesian data, will be neglected here. More important to
our discussion is the removal of the aliased energy discussed above. The unmixing operator uses a linear combination of voxel signals to remove this coupling between
voxels. If the refocused image contains little aliasing, the
mixing matrix is close to the identity matrix.
Refocusing Reconstruction
Calculation of the unmixing operator M⫺1 is problematic
due to its computational complexity. Not only is the N 2 â«»
N 2 inversion computationally expensive, it also cannot be
precomputed because it depends on the phase corruptions. Further, the calculation of the mixing matrix M itself
requires CRN 2 ⬇ N 4 operations for each element of the
N 2 â«» N 2 matrix (see Eq. [12]). These computational requirements make direct calculation of the unmixing matrix (and hence the full least-squares reconstruction) intractable.
If the sources of coupling are sufficiently low, the mixing matrix has little off-diagonal energy and is close to the
identity matrix. Approximating the mixing matrix as the
identity yields the “refocusing reconstruction”:
m̃ ⫽ P*R*d.
[3]
Provided the unmixing matrix is nearly diagonal, the refocusing reconstruction is a good approximation to the
least-squares estimate that is much faster to compute. The
refocusing reconstruction is also conceptually elegant: before the readouts are combined, each readout is refocused
in image space by multiplication with the navigator phase
conjugate. This interpretation is not only computationally
useful, but also suggests a more efficient algorithm than
the matrix computation in Eq. [3] for refocusing the data.
This algorithm will be discussed later.
Some important properties of the refocusing reconstruction can be understood by considering the problem in
k-space. Corruption by a spatial phase function (e i␾(r)) is
equivalent to convolution by the transform of that phase
function (Ᏺ{e i␾(r)}) in the spectral domain. If all of k-space
were convolved with the same phase corruption kernel,
we could restore the original spectrum by deconvolving
with this kernel. However, in motion-corrupted DWI data,
each readout samples a spectrum that has been convolved
with a different kernel (Ᏺ{e i␾ j(r)} for the jth readout).
Interpreted in k-space, the refocusing operator deconvolves each interleave separately using the phase-corruption kernel from that excitation. Refocusing uses the correct deconvolution kernel for points from the same excitation. For data from different excitations, the
deconvolution will work well if the kernels are similar for
adjacent positions in k-space, but will be unable to remove
artifacts if the kernels are dissimilar. Hence, it is desirable
for the kernel to change slowly and smoothly over k-space
so that the deconvolution is locally consistent. An additional concern is periodicity in the phase corruptions that
causes a periodic weighting of the data in k-space. Any
residual periodicity left by the refocusing operator will
cause coherent ghosts, which are a more disruptive artifact
than the diffuse aliasing associated with discontinuities.
The refocusing reconstruction will produce high-fidelity
346
Miller and Pauly
(f) Multiply (e) by the phase conjugate of (c) (refocusing step).
(g) Add (f) to the sum accumulating in step 1.
3. The refocused image is the cumulative sum from (g).
For every interleave, this reconstruction requires two
Fourier transforms and a matrix multiplication at the highresolution matrix size (N â«» N). Each Fourier transform
requires time N 2 log N 2 , and the matrix multiplication
requires time N 3 for a total complexity of O(CN 2 log N 2 ⫹
CN 3 ) ⬇ O(CN 3 ).
The k-space refocusing algorithm (Fig. 3b) is:
FIG. 3. Block diagram of two versions of the refocusing reconstruction in which refocusing is performed in (a) image space and (b)
k-space. The boxed expressions indicate the operations of gridding
(Ᏻj grids the jth readout trajectory), inverse Fourier transform (Ᏺ⫺1),
and refocusing (multiplication by ᏼj*(r) ⫽ e⫺i␾j(r) in image-space or
convolution by Ᏺ{ᏼj* } in k-space).
images when phase corruptions are smooth and nonperiodic in k-space.
MATERIALS AND METHODS
Efficient Refocusing Reconstruction
Although the matrix formulation given in the previous
section is a useful way to derive the refocusing algorithm,
it is an inefficient method for calculating the refocused
image. Directly implementing the matrix formulation has
prohibitive memory requirements and is unlikely to take
full advantage of the highly structured nature of these
matrices. It is much more efficient to use operations equivalent to those given in the matrix formulation, since these
operations often have highly optimized algorithms (such
as the Fast Fourier Transform (FFT)). The implementation
of image-space refocusing shown in Fig. 3a can be calculated more efficiently than the matrix calculation of m̃
shown in Fig. 2. Refocusing can also be implemented as a
k-space deconvolution, as shown in Fig. 3b.
The image-space refocusing algorithm in Fig. 3a involves the following steps:
1. Initialize a zero-filled N â«» N matrix (where we will
accumulate the refocused image).
2. For each interleave:
(a) Reconstruct the navigator data in k-space (n â«»
n).
(b) Zero-pad (a) to the final image resolution (N â«»
N).
(c) Inverse Fourier transform (b) (to get the navigator
image).
(d) Reconstruct the high-resolution data frame in kspace at the final image resolution (N â«» N).
(e) Inverse Fourier transform (d) (to get the interleave
image).
1. Initialize a zero-filled N â«» N matrix (where we will
accumulate the refocused spectrum).
2. For each interleave:
(a) Reconstruct the navigator data in k-space (n â«»
n).
(b) Zero-pad (a) by a factor of 2 (2n â«» 2n).
(c) Inverse Fourier transform (b).
(d) Calculate the phase conjugate of (c).
(e) Fourier transform (d) (to get the refocusing kernel).
(f) Reconstruct the high-resolution data frame in kspace at the final image resolution (N â«» N).
(g) Convolve (f) with (e) (refocusing step).
(h) Add (g) to the sum accumulating in step 1.
3. Inverse Fourier transform the cumulative sum from
(h) to get the refocused image.
The refocusing kernel must be calculated at twice the
navigator bandwidth (step 2b) because normalizing in image-space (step 2d) extends the kernel in k-space. Calculating the kernel at twice the navigator bandwidth avoids
Gibbs ringing in the refocused image due to this extension.
For every interleave, the k-space algorithm calculates two
2n â«» 2n Fourier transforms and one convolution of an
N â«» N matrix with a 2n â«» 2n kernel. This algorithm runs
in approximately O(CN 2 n 2 ⫹ Cn 2 log n 2 ) ⬇ O(CN 2 n 2 )
time (where constant multipliers have been dropped, as is
standard in O-notation).
In most imaging situations, the k-space algorithm will
perform better than the image-space algorithm, since we
will usually have n 2 ⬍ N. However, this is not always the
case, as will be discussed later.
Cardiac Synchronization
If minor subject restraints are used to reduce bulk motion,
the primary source of motion artifacts in the brain is deformation related to the cardiac cycle (see Fig. 1). Cardiac
triggering information can be used to synchronize interleave ordering to the cardiac cycle to reduce discontinuities and periodicities in the k-space data (20,21). For example, in a spin-warp trajectory, phase encodes are normally gathered incrementally, resulting in a periodic
weighting over k-space with large discontinuities (see Fig.
4a). These effects can be lessened by choosing a phaseencode increment that traverses the extent of k-space during a single cardiac cycle. By interleaving the phase encodes over multiple cardiac cycles in this manner, adjacent lines in k-space are collected at the same portion of
the cardiac cycle, and the resultant k-space weighting is
smooth and nonperiodic (Fig. 4b). This method for cardiac
Nonlinear Phase Correction for DWI
347
a total scan time of 6:15 per image. To date, 10 subjects
have been scanned with variants of this sequence, five of
whom were scanned with the prescription specified above.
In some experiments, acquisition was synchronized to
the cardiac cycle using plethysmographic triggering. To
maintain the steady state, dummy cycles without acquisition were gathered while waiting for cardiac triggers. Following a systolic trigger, data were collected in an interleaved order such that the phase encodes acquired during
a single cardiac cycle were evenly distributed across k y space. This process was repeated for successive cardiac
cycles until all phase encodes had been acquired. Additionally, two to three dummy cycles (80 –120 ms) were
inserted following detection of a trigger to avoid the largest
motions associated with systole. Cardiac synchronization
usually increased scan times by 1–2 min.
Images were reconstructed using three methods: no correction, the standard linear correction, and the k-space
refocusing correction proposed in this work. Reconstruction of the data was implemented off-line in C code using
the FFTW library for Fourier transforms (22). Both the
image-space and k-space versions were implemented;
however, all refocused images shown here were calculated
with the k-space implementation.
FIG. 4. The effect of 2DFT phase-encode ordering on k-space. The
phase history of a representative voxel over time is shown in the middle
plot. The gray swathes represent cardiac-gated imaging windows. The
k-space modulation imparted by this phase history is shown for (a)
standard incremental phase-encode ordering and (b) an alternate
phase-encode ordering. The latter ordering interleaves the data from
different cardiac cycles to achieve a smoother modulation.
synchronization of phase encoding should not require
large increases in scan time, since data can be acquired
throughout a large portion of the cardiac cycle. Cardiac
synchronization in this manner reduces the types of kspace modulations that are most problematic for the refocusing reconstruction— discontinuity and periodicity—
without incurring the large increases in scan time that
usually accompany cardiac gating.
RESULTS
Refocusing Reconstruction
The different navigator corrections are compared on untriggered data in Fig. 6. Motion-induced phase artifacts
interfere to cause signal attenuation in areas with large
motions, and little or no artifact outside the object. Motion
artifacts are manifested as signal loss rather than the more
familiar ghosting artifacts due to the large number of excitations in the acquisition. All of the uncorrected images
(top row) are corrupted, with the most serious artifacts in
the S/I-weighted image. Arrows in the top row indicate
artifacted areas in the R/L- and A/P-weighted images. The
linear correction (middle row) restores signal in some corrupted regions (e.g., posterior areas in the S/I-weighted
image), but the images remain attenuated in the medial
Experiments
Brain images were acquired on healthy volunteers with the
navigated steady-state diffusion-weighted imaging (SSDWI) pulse sequence shown in Fig. 5. To minimize large
motions, padding was placed at the subject’s temples. The
readout consisted of a spiral navigator (FOV ⫽ 24 cm,
matrix ⫽ 16 ⫻ 16, in-plane resolution ⫽ 15 ⫻ 15 mm2,
receive bandwidth (BW) ⫽ 62.5 kHz) followed by a spinwarp high-resolution acquisition (FOV ⫽ 24 cm, matrix ⫽
192 ⫻ 192, in-plane resolution ⫽ 1.25 ⫻ 1.25 mm2, receive
BW ⫽ 7.8125 kHz). Experiments were performed on a
1.5 T GE Signa CV/i research scanner (General Electric Co.,
Milwaukee, WI) with 40 mT/m gradients switchable at
150 mT/m/ms and a transmit/receive quadrature birdcage
head coil. DW images were acquired using diffusion gradients with G ⫽ 40 mT/m and ␶ ⫽ 6.5 ms (b eff ⫽ 980 s/
mm2 (19)), which could be applied along any direction.
Some studies also acquired an image with minimal diffusion weighting (G ⫽ 0.3 mT/m and ␶ ⫽ 6.5 ms). Other
imaging parameters were: TR ⫽ 40 ms, ␣ ⫽ 30°, and
thickness ⫽ 6 mm. The sequence was repeated 48 times for
FIG. 5. Navigated steady-state DWI pulse sequence used in this
study. The sequence is fully refocused except for the diffusion
gradient (with amplitude G and duration ␶), which rephases DW
echoes over multiple TRs.
348
Miller and Pauly
FIG. 6. Untriggered spiral-navigated
2DFT images with varying levels of correction. Rows (from top to bottom) have
no correction, the standard linear correction, and the refocusing correction.
Columns (from left to right) show images
with diffusion weighting along the A/P,
R/L, and S/I directions. Arrows in the
A/P- and R/L-weighted uncorrupted images (top row) indicate examples of
dephasing due to motion. Arrows in the
middle row indicate artifacts introduced
by the linear correction. Arrows in the
bottom row indicate structures with increased conspicuity after the refocusing
correction.
areas where motion is more severe. Additionally, the linear correction introduces attenuation to the frontal lobes in
the R/L- and A/P-weighted images due to a poor linear fit
to the nonlinear phase artifacts (arrows in middle row of
Fig. 6). The refocused images (bottom row) have restored
signal and minimal artifacts. In particular, the medial regions show structures that were partially or completely
suppressed in the uncorrected images (e.g., the external
and internal capsule, indicated by arrows). The images
weighted along the primary axis of motion (S/I) have a
lower signal-to-noise ratio (SNR) due to the large motions
in the S/I direction. In peak systole, these motions cause
intravoxel dephasing, which cannot be corrected by postprocessing.
FIG. 7. Cardiac-synchronized, spiral-navigated 2DFT images acquired
with the phase-encode interleaving
scheme shown in Fig. 4. Diffusion
weighting is along the (a) A/P, (b) R/L,
and (c) S/I directions.
Cardiac Triggering
A set of refocused images acquired with cardiac synchronization is shown in Fig. 7. These images show strong
diffusion contrast with clear delineation of small whitematter structures. The S/I-weighted image does not suffer
from the degraded SNR seen in untriggered images. The
A/P- and R/L-weighted images are also improved with
triggering, and show fine detail (such as cortical white
matter) that is less apparent in the untriggered images.
Diffusion Tensor Imaging
A full set of diffusion tensor images acquired with the
standard directions for DTI (23) is shown in Fig. 8. These
Nonlinear Phase Correction for DWI
349
FIG. 8. Navigated SS-DWI images
acquired with diffusion weighting
along the standard directions for diffusion tensor. Images are cardiacsynchronized and were reconstructed using the refocusing reconstruction.
images were acquired with cardiac-synchronized phaseencode ordering, and reconstructed using the refocusing
reconstruction. Intravoxel dephasing is minor in these images due to the systolic delay (120 ms) and the fact that
diffusion weighting was never applied along the principal
axis of brain motion (S/I). These images show strong diffusion weighting with little visible artifact. Figure 9 shows
the trace image and fractional anisotropy (FA) map for
these data before and after correction. The correction improves both images, although the most dramatic improvement is found in the FA map.
DISCUSSION
We have presented a method for correction of nonlinear
phase corruption. This method does not require additional
FIG. 9. Trace and FA maps for the data in Fig.
8 before and after the refocusing correction. Images
before and after correction are identically windowed.
The trace image shows a moderate improvement in
contrast. The FA map shows a more marked improvement. The uncorrected FA map shows apparent anisotropy in gray-matter regions due to motioninduced signal loss in a subset of the images. The
refocusing correction restores this signal to give a
more accurate map of anisotropy.
350
Miller and Pauly
navigator data beyond those used for 2D linear corrections.
Rather, the refocusing correction represents a more complete use of the information in the navigator. The refocusing reconstruction algorithm is straightforward to implement: refocusing is simply a k-space deconvolution that is
performed before individual readouts are combined. Alternatively, it can be thought of as image-space multiplication of individual data frames by the navigator phase conjugate.
Generality of Reconstruction
The least-squares reconstruction is general in the sense
that its derivation does not make assumptions about the
trajectory or pulse sequence used to gather the data. Although the refocusing reconstruction will not always be a
good approximation of the least-squares solution, it is
expected to perform well provided the k-space modulation
induced by phase corruptions has minimal discontinuity
and periodicity. In addition to the 2DFT data shown here,
we previously used the refocusing reconstruction for selfnavigated spiral data with similar results (24). Pipe and
colleagues (14) have also had good success using a refocusing reconstruction with an FSE PROPELLER sequence
(however, in their case, the refocusing reconstruction was
the same as the least-squares reconstruction, since each
readout was unaliased).
While full least-squares reconstruction is generally intractable due to the inversion of the mixing matrix, it may
be manageable in special cases or with efficient approximations. For example, the mixing matrix M will be block
diagonal and sparse for trajectories in which the signal
from a voxel can only alias onto a subset of the other
voxels (see Appendix B). In this case, the matrix inversion
of M can be calculated by individually inverting the
smaller block submatrices. In our experience, this calculation is still very time-consuming. Alternatively, the mixing matrix could be approximated with an iterative
scheme similar to that proposed for multicoil reconstructions (25).
Fidelity of the Refocusing Reconstruction
The impulse response of the mixing operator is a useful
way to characterize the fidelity of the refocused image. The
off-peak power in the impulse response indicates the
amount of aliased energy in the refocused image, and
therefore the relative benefit of applying the unmixing
operator. Of particular interest are coherences in the voxel
impulse responses, which indicate coherent ghosts in the
refocused image. The impulse response of the mixing operator is shown in Fig. 10a for several representative voxels in an ungated data set. Although there is significant
energy off-center, most of the aliased signal has random
phase and will add incoherently. The impulse responses
also show widening of the central peak and coherent offcenter peaks, which will lead to blurring and ghosting,
respectively. These effects in the impulse response must
be reduced in order for the refocusing operator to provide
accurate reconstructions.
The quality of the refocusing reconstruction is dependent on the number of excitations that make up an image.
Aliased energy in the refocused image decreases as the
FIG. 10. Impulse responses of the mixing operator for a 2DFT
acquisition with (a) no averaging and (b) four averages. Each subfigure overlays the impulse responses for the same 10 voxels taken
from a representative navigator data set. The arrow in a indicates
coherences in the mixing impulse responses that lead to coherent
ghosts. Averaging in b suppresses much of the off-peak energy.
number of excitations increases, since the incoherent
sources of aliasing are sampled more often. As shown in
Fig. 10b, averaging suppresses the incoherent off-peak energy found in unaveraged data. In this case, averaging also
suppresses the coherent peaks since data were acquired
asynchronous to the cardiac cycle. Some averaging is always necessary to correct strong phase corruptions, which
warp the k-space trajectory. Acquiring multiple averages
ensures that the perturbed encoding functions will cover
k-space more fully. One side effect of this need to average
is the improvement of the refocusing reconstruction.
Cardiac Synchronization
The need for cardiac gating has been a common criticism
of navigator methods. Since linear corrections assume
Nonlinear Phase Correction for DWI
351
Computational Efficiency
FIG. 11. Simulated impulse response magnitude of the mixing operator for a 2DFT acquisition with and without cardiac synchronization (nav ⫽ 1). a: The impulse responses for standard incremental
order have strong coherences (arrows) that will lead to ghosting. b:
These coherences are eliminated when cardiac-synchronized
phase-encode ordering is used.
rigid-body motion, cardiac gating is required to limit
imaging to the more stable portions of the cardiac cycle,
in which motions are better approximated as linear.
With the nonlinear corrections introduced here, postsystolic delays only have to be long enough to avoid
intravoxel dephasing during peak systole, which cannot
be corrected in postprocessing. For the 2DFT SS-DWI
sequence used here, cardiac synchronization requires
only moderate increases in scan time (25–33%). This
imaging efficiency is possible because cardiac information is being used to synchronize the acquisition to the
cardiac cycle, rather than simply limiting the portion of
the cardiac cycle in which data are gathered. Our current implementation approximates the heart rate for the
entire scan. Monitoring changes in heart rate during
scanning would further reduce the overhead time for
cardiac synchronization.
Improvements to the refocused image due to cardiac
synchronization can be seen by examining the impulse
response of the mixing operator. The impulse response for
the standard phase-encode ordering scheme (Fig. 11a) has
significant off-peak energy, as well as a series of coherent
peaks (arrows) due to the periodic weighting of k-space
introduced by the cardiac cycle (several harmonics are
clearly visible). Cardiac-synchronized ordering (Fig. 11b)
eliminates off-center coherent peaks and produces a more
well-behaved impulse response. The attenuation band
around the impulse is due to the suppression of low frequencies in the k-space modulation.
Even though they are functionally equivalent, the imagespace and k-space refocusing reconstructions differ in
terms of computational efficiency. The choice of algorithm
depends on the relative sizes of the navigator and highresolution acquisitions, as well as the details of the highresolution trajectory. However, in almost all practical imaging scenarios, the k-space implementation is expected to
outperform the image-space algorithm. As presented
above, the k-space algorithm will perform better if n 2 ⬍ N,
because the k-space convolution will be faster than the
large matrix multiply.
An additional consideration in the comparison of the
two methods is the specific k-space trajectory used for the
high-resolution acquisition. Our previous analysis of the
k-space algorithm assumed that the convolution must be
performed over the entire N â«» N matrix, which is not true
for all trajectories. For some trajectories, the convolution
time can be significantly reduced by performing a local
convolution. For example, with a spinwarp trajectory, the
convolution need only include the area around each
phase-encode line, and can therefore be performed in
O(Nn 2 ) time. For this trajectory, the k-space algorithm
will in general perform better when n ⬍ N, which should
always be the case.
Thus far we have only considered 2D acquisitions. For
3D acquisitions, all 2D FFTs, matrix multiplications, and
convolutions are replaced with 3D operations. The k-space
algorithm is expected to scale more favorably than the
image-space algorithm for 3D acquisitions, since a local 3D
convolution can be performed significantly faster than
large 3D Fourier transforms and matrix multiplications
can be calculated. Early work with 3D SS-DWI knee imaging used a k-space algorithm with good results (26). For
this application, the image-space algorithm was found to
be prohibitively slow.
For the experiments performed here (a spinwarp trajectory with n ⫽ 16, N ⫽ 192), the k-space algorithm would
be expected to drastically outperform the image-space algorithm. Our implementation computed the k-space reconstruction in 2.6 s per average, and the image-space
reconstruction in 6.9 s per average. The modest improvement of the k-space algorithm is likely due to the use of
highly-optimized FFT software and unoptimized convolution calculations.
Weighted Refocusing Reconstruction
Previous work on navigated DWI usually included a
weighting or thresholding of data frames that are suspected to contain uncorrectable phase artifacts. The
method proposed here does not include such a correction,
primarily because it has not proven necessary under the
current imaging circumstances. However, imaging methods that acquire fewer data frames are more impacted by
artifacted data frames, and might benefit from a weighting
scheme.
Examining Fig. 1, we see that during periods of the most
intense motion (such as systole), navigator images exhibit
signal loss due to intravoxel dephasing. Navigator magnitude can thus be used as an indicator of uncorrectable
phase corruption and can be incorporated into the recon-
352
Miller and Pauly
struction. If we weight each image-space sample by its
navigator magnitude prior to correction, the least-squares
estimate will be biased to more reliable data. This can be
accomplished by solving the preconditioned problem:
Qd ⫽ QRPm
[4]
where Q is a preconditioning matrix that accomplishes the
desired image-space weighting. Note that preconditioning
is applied to the problem in k-space.
The least-squares solution to Eq. [4] is:
m̂ Q ⫽ 共P*R*Q*QRP兲⫺1 P*R*Q*Qd.
[5]
The estimate in Eq. [5] filters the k-space data by Q*Q prior
to refocusing. The magnitude effects of the filter are then
removed by the unmixing matrix (P*R*Q*QRP)⫺1. To
weight each image-space voxel by the navigator magnitude, Q*Q must convolve each high-resolution interleave
with the Fourier transform of the navigator magnitude. In
other words, R*Q*Qd should contain the interleave images weighted by the navigator magnitude.
A weighted refocusing reconstruction can then be defined by approximating the mixing matrix P*R*Q*QRP as
diagonal by zeroing the off-diagonal elements. The approximated unmixing matrix removes the weighting caused by
the preconditioning, but does not attempt to remove any
aliasing. Regularization can be used to improve the stability of the inversion in noise regions with low signal.
single image at roughly the same portion of the cardiac
cycle. The cardiac cycle could then be filled out by acquiring more slices (in multislice), more z-phase encodes (for
3D), or more averages. All of these methods use cardiac
information to achieve well-behaved modulations of the
k-space data without requiring large increases in scan
time.
CONCLUSIONS
We have presented a method for correcting nonlinear
phase errors in diffusion-weighted imaging using 2D navigation. The refocusing reconstruction, which is an approximation to the computationally-expensive leastsquares estimate, is a straightforward algorithm that makes
complete use of the information in a 2D navigator. The
heart of the reconstruction is the refocusing operator,
which rephases the unaliased component of each acquisition via multiplication with the navigator phase conjugate.
This reconstruction is a good approximation to the full
least-squares estimate provided the phase errors cause
smooth, nonperiodic modulations of k-space, and a sufficient number of excitations are acquired. Cardiac triggering information can be used to achieve desirable modulations through ordering of data interleaves. The refocusing
reconstruction is expected to be applicable in a wide range
of imaging conditions.
APPENDIX A
Application to Other Sequences and Trajectories
Structure of Reconstruction Matrices
The refocusing reconstruction will provide a good reconstruction of the data provided 1) a large number of excitations are acquired, and 2) k-space modulations induced by
motion are smooth and nonperiodic.
The requirement for many excitations minimizes aliasing energy by allowing the incoherent phase errors to
cancel. We have already discussed the improvement in the
refocused image from averaging; other methods of increasing the number of excitations will have the same effect. For
example, a 3D acquisition should be better-conditioned to
a refocusing reconstruction than a 2D multislice acquisition. Simply increasing the number of interleaves in an
acquisition will also improve the refocused image. The
need for a large number of excitations is a disadvantage for
some sequences, such as spin echo, which require long
echo times. These methods may benefit from a weighted
refocusing reconstruction, which will tend to suppress
these errors, particularly if regularization is used.
The requirement for smooth and nonperiodic k-space
modulations ensures that the refocusing operator will be
able to remove phase corruptions. This requirement can be
addressed in most trajectories with cardiac synchronization. Some trajectories cannot simultaneously achieve
smoothness and nonperiodicity, in which case we prefer
to meet the smoothness constraint since the refocusing
reconstruction does a poor job of removing discontinuities. Spiral and radial scans can achieve a smooth k-space
weighting by distributing the rotated interleaves from 0° to
180° over each cardiac cycle (this scheme would result in
azimuthal periodicity in k-space). Segmented EPI scans
can also be made smooth if we acquire each segment of a
The motion-corrupted magnetization during the cth readout is denoted as m c (r) ⫽ 兩m c (r)兩e i␾ c(r). The phase corruption during the cth excitation is expressed by the N 2 ⫻ N 2
matrix Pc:
Pc ⫽
冢
ei␾ c共r1兲
0
···
0
0
·
·
·
0
ei␾ c共r2兲
·
·
·
0
···
0
·
·
·
···
冣
.
ei␾ c共rN 兲
2
[6]
Pcm is the phase-corrupted magnetization that is the
source of signal during the cth readout. Each of the C
interleaves has a different phase corruption, represented
by stacking the Pc’s to form P:
冢冣
P1
P⫽
P2
·
·
·
PC
.
[7]
P splits the magnetization vector m into C separate phasecorrupted magnetization vectors. These vectors are
stacked to form the CN 2 â«» 1 vector Pm representing the
phase-corrupted magnetization that is sampled over the
course of the experiment.
Each of the phase-corrupted magnetization vectors is
separately Fourier transformed by the matrix F (CN 2 â«»
CN 2 ):
Nonlinear Phase Correction for DWI
F⫽
冢
FN
0
···
0
·
·
·
0
N
···
F
·
·
·
0
353
0
0
·
·
·
· · · FN
冣
REFERENCES
.
[8]
F is composed of C 2DFT matrices FN (N 2 â«» N 2 ) consisting of the Fourier kernels:
N
F lm
⫽ ei2␲kl䡠rm.
[9]
The product FPm is the CN 2 â«» 1 vector representing the C
phase-corrupted spectra, sampled in Cartesian coordinates.
Finally, the spectra given by FPm are resampled onto the
k-space sampling trajectories by the resampling matrix G
(CR â«» CN 2 ):
G⫽
冢
G1
0
0
·
·
·
0
G2 · · · 0
·
·
·
·
·
·
0 · · · GC
···
0
冣
.
[10]
Each submatrix Gc (R â«» N 2 ) interpolates the Cartesian
k-space points kl onto the sampling points for the cth
c
readout sm
:
c
⫽ sinc共kl ⫺ scm 兲.
G lm
[11]
Note that for trajectories that collect data on Cartesian
coordinates (such as spinwarp or interleaved EPI), the
elements of G are discrete delta functions.
APPENDIX B
Mixing Matrix
Using the notation of Appendix A, the elements of the
mixing matrix M for an arbitrary trajectory are:
冘
C
M lm ⫽
c⫽1
i共␾ c 共rm 兲⫺␾ c 共rl 兲兲
e
冘冋冘
R
N2
r⫽1
n⫽1
sinc共kn ⫺ src 兲e⫹i2␲kn䡠rm
冋冘
N2
â«»
n⫽1
册
册
sinc共kn ⫺ src 兲e⫺i2␲kn䡠rl . [12]
The summation over each readout r describes non-idealities in the k-space trajectory that causes coupling between
voxels within a readout. The phase difference describes
the aliasing caused by motion-induced phase errors. These
phase offsets are manifested as perturbations of the encoding functions used in the trajectory.
ACKNOWLEDGMENTS
K.M. thanks Dr. Brian Hargreaves for useful discussions
regarding the method and manuscript.
1. Anderson AW, Gore JC. Analysis and correction of motion artifacts in
diffusion weighted imaging. Magn Reson Med 1994;32:379 –387.
2. Ordidge RJ, Helpern JA, Qing ZX, Knight RA, Nagesh V. Correction of
motion artifacts in diffusion-weighted MR images using navigator echoes. Magn Reson Imaging 1994;12:455– 460.
3. deCrespigny AJ, Marks MP, Enzmann DR, Moseley ME. Navigated
diffusion imaging of normal and ischemic human brain. Magn Reson
Med 1995;33:720 –728.
4. Butts K, Crespigny A, Pauly JM, Moseley M. Diffusion-weighted interleaved echo-planar imaging with a pair of orthogonal navigator echoes.
Magn Reson Med 1996;35:763–770.
5. Mori S, van Zijl PCM. A motion correction scheme by twin-echo navigation for diffusion-weighted magnetic resonance imaging with multiple RF echo acquisition. Magn Reson Med 1998;40:511–516.
6. Williams CFM, Redpath TW, Norris DG. A novel fast split-echo multishot diffusion-weighted MRI method using navigator echoes. Magn
Reson Med 1999;41:734 –742.
7. Brockstedt S, Moore JR, Thomsen C, Holtas S, Stahlberg F. Highresolution diffusion imaging using phase-corrected segmented echoplanar imaging. Magn Reson Imaging 2000;18:649 – 657.
8. Clark CA, Barker GJ, Tofts PS. Improved reduction of motion artifacts in
diffusion imaging using navigator echoes and velocity compensation. J
Magn Reson 2000;142:358 –363.
9. Norris DG, Driesel W. Online motion correction for diffusion-weighted
imaging using navigator echoes: application to RARE imaging without
sensitivity loss. Magn Reson Med 2001;45:729 –733.
10. Bosak E, Harvey PR. Navigator motion correction of diffusion weighted
3D SSFP imaging. MAGMA 2001;12:167–176.
11. Jiang H, Golay X, van Zijl PCM, Mori S. Origin and minimization of
residual motion-related artifacts in navigator-corrected segmented diffusion-weighted EPI of the human brain. Magn Reson Med 2002;47:
818 – 822.
12. Butts K, Pauly J, Crespigny A, Moseley M. Isotropic diffusion-weighted
and spiral-navigated interleaved EPI for routine imaging of acute
stroke. Magn Reson Med 1997;38:741–749.
13. Atkinson D, Porter DA, Hill DLG, Calamante F, Connelly A. Sampling
and reconstruction effects due to motion in diffusion-weighted interleaved echo planar imaging. Magn Reson Med 2000;44:101–109.
14. Pipe JG, Farthing VG, Forbes KP. Multishot diffusion-weighted FSE
using PROPELLER MRI. Magn Reson Med 2002;47:42–52.
15. Norris DG. Implications of bulk motion for diffusion-weighted imaging
experiments: effects, mechanisms and solutions. J Magn Reson Imaging
2001;13:486 – 495.
16. Enzmann DR, Pelc NJ. Brain motion: measurement with phase-contrast
MR imaging. Radiology 1992;185:653– 660.
17. Greitz D, Wirestam R, Franck A, Nordell B, Thomsen C, Stahlberg F.
Pulsatile brain movement and associated hydrodynamics studied by
magnetic resonance phase imaging. Neuroradiology 1992;34:370 –380.
18. Poncelet BP, Wedeen VJ, Weisskoff RM, Cohen MS. Brain parenchyma
motion: measurement with cine echo-planar MR imaging. Radiology
1992;185:645– 651.
19. Buxton RB. The diffusion sensitivity of fast steady-state free precession
imaging. Magn Reson Med 1993;29:235–243.
20. Bailes DR, Gilderdale DJ, Bydder GM, Collins AG, Firmin DN. Respiratory ordered phase encoding (ROPE): a method for reducing respiratory motion artifacts in MR imaging. J Comput Assist Tomogr 1985;9:
835.
21. Haacke EM, Patrick JL. Reducing motion artifacts in two-dimensional
Fourier transform imaging. Magn Reson Imaging 1986;4:359 –376.
22. Frigo M, Johnson SG. FFTW: an adaptive software architecture for the
FFT. Proc ICASSP 1998;3:1381–1384.
23. Basser PJ, Pierpaoli C. A simplified method to measure the diffusion
tensor from seven MR images. Magn Reson Med 1998;39:928 –934.
24. Miller KL, Meyer CH, Pauly JM. Self-navigated spirals for high resolution steady-state diffusion imaging. In: Proceedings of the 10th Annual
Meeting of ISMRM, Honolulu, 2002.
25. Preussmann KP, Weiger M, Bornert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med
2001;46:638 – 651.
26. Miller KL, Hargreaves BA, Gold GE, Pauly JM. Navigated 3D steadystate diffusion imaging of knee cartilage. In: Proceedings of the 11th
Annual Meeting of ISMRM, Toronto, 2003.
Magnetic Resonance in Medicine 50:343–353 (2003)
Nonlinear Phase Correction for Navigated
Diffusion Imaging
Karla L. Miller* and John M. Pauly
Motion during diffusion-weighted imaging (DWI) introduces
phase errors that can cause significant artifacts in brain images. One method of correcting these errors uses additional
navigator data to measure the phase corruptions. Standard
navigator methods correct for rigid-body motion but cannot
correct for nonrigid deformations of the brain related to the
cardiac cycle. This work derives a generalized reconstruction
that corrects for nonrigid motion based on a least-squares
formulation. Since this reconstruction has the disadvantage of
being computationally expensive, an approximation is presented, called a refocusing reconstruction. The refocusing reconstruction is both efficient and straightforward. Each readout
is multiplied in image space by the phase conjugate of the
navigator image, and these rephased readouts are then
summed. The conditions under which the refocusing reconstruction is sufficient are considered and methods to improve
the quality of refocused images are discussed. In particular,
synchronization of the acquisition to the cardiac cycle can
provide data that is well-conditioned to the refocusing reconstruction without incurring the large time penalty traditionally
associated with cardiac gating. These methods are applied to
steady-state DWI, a promising pulse sequence that is particularly sensitive to motion-induced phase artifacts. The refocusing reconstruction is shown to significantly improve SS-DWI
over standard rigid-body corrections. Magn Reson Med 50:
343–353, 2003. © 2003 Wiley-Liss, Inc.
Key words: diffusion MRI; motion correction; artifact reduction;
navigator; reconstruction
Diffusion-weighted imaging (DWI) is a powerful tool for
probing tissue microstructure that uses large bipolar gradient pairs to impart phase to moving particles. Diffusing
water spins within a voxel experience phase dispersion,
causing signal attenuation that is dependent on the
amount of diffusion. This attenuation results in diffusion
contrast, which has been shown to be useful in brain
imaging for detecting tissue disease states, as well as revealing structural information. A key problem in DWI is
nondiffusive bulk motion, which causes phase shifts in the
voxel signal (1). When data are gathered in a multishot
acquisition, phase shifts from different excitations can destructively interfere to cause serious image artifacts.
A number of methods have been proposed to deal with
motion artifacts in DWI. One simple way to avoid motioninduced artifacts is to gather the image data in a single
acquisition, in which case the image magnitude is uncor-
Magnetic Resonance Systems Research Laboratory, Department of Electrical
Engineering, Stanford University, Stanford, California.
Grant sponsor: NIH; Grant numbers: R01 CA50948; R01 AR46904.
*Correspondence to: Karla L. Miller, Packard Electrical Engineering Building,
Room 210, Stanford University, Stanford, CA 94305-9510. E-mail:
bison@stanford.edu
Received 29 August 2002; revised 28 March 2003; accepted 2 April 2003.
DOI 10.1002/mrm.10531
Published online in Wiley InterScience (www.interscience.wiley.com).
© 2003 Wiley-Liss, Inc.
rupted. However, these single-shot acquisitions have limited resolution and are sensitive to susceptibility and
eddy-current effects. Another approach is to correct motion-induced phase artifacts before combining the data
from multiple acquisitions. Phase correction imposes no
intrinsic limits on resolution, but does require the acquisition of additional navigator data that resolve the phase
errors in image space. Early works acquired 1D projection
navigators (1–3) or pairs of orthogonal projections (4), and
subsequent studies have extended this approach to a variety of pulse sequences and trajectories (5–11). An important disadvantage of 1D navigation is the inability to correct for phase shifts orthogonal to the direction of navigation. Recent studies have addressed this problem with 2D
navigation (12–14).
Navigated DWI methods typically assume rigid-body
motion, the effects of which have been well characterized
(1,15). However, a significant component of brain motion
is nonrigid (16 –18). During systole, the brain deforms
about the ventricles and displaces inferiorly through the
foramen magnum as though being pulled by the spinal
cord. The more superior and lateral regions of the brain are
nearly motionless, whereas the more inferior and medial
regions tend to experience large displacements. Even during the relative quiescence of diastole, motion induced by
the cardiac cycle follows this nonlinear spatial pattern as
the brain slowly returns to its resting position.
In the presence of diffusion-weighting gradients, these
deformations cause phase corruptions that vary nonlinearly across the object magnetization. The effect of brain
motion on the DWI signal is shown in Fig. 1. The signal
phase and amplitude are correlated with the cardiac cycle
and exhibit nonlinear spatial variation. The loss in signal
amplitude during peak systole is due to intravoxel dephasing caused by intense motion. The highly focal nature of
this dephasing highlights the spatial nonlinearity of phase
corruptions.
These phase errors are problematic for multishot acquisitions because each readout experiences a different image-space phase corruption (␾ j (r) for the jth readout). This
causes phase interference when readouts are combined to
form a high-resolution image. Navigated techniques correct phase errors with low-resolution navigator data that
are acquired along with each high-resolution data frame
(2). Provided the navigator data 1) samples the same phase
error ␾ j (r) as the high-resolution data, and 2) covers kspace sufficiently to resolve ␾ j (r), the navigator can be
used to correct phase errors in the high-resolution data
before the interleaves are combined. Based on the assumption of rigid-body motion, standard navigator methods correct for zeroth- and first-order phase terms (1). In the
presence of the nonrigid motions described above, a rigidbody correction is a first-order approximation to the de-
343
344
Miller and Pauly
function, and can be described conveniently using matrix
formalism. For an arbitrary trajectory that covers k-space
in C interleaves of length R, corresponding to an N â«» N
image, the acquired data are:
d ⫽ GFPm ⫽ RPm
FIG. 1. Illustration of the effects of motion on DW magnetization.
The magnitude (top) and phase (bottom) of a series of low-resolution
axial single-shot DWI images (2D navigators) are shown over a
single cardiac cycle. The magnitude and phase of a representative
voxel are plotted in the middle over four cardiac cycles. The spatial
variation in the image magnitude and phase is nonlinear and
strongly correlated to the cardiac cycle.
sired high-order correction. Rigid-body corrections leave
residual errors that are strongly correlated to the cardiac
cycle (11). Given that a 2D navigator contains high-order
phase information, image reconstruction should be improved by correcting for high-order phase terms.
This work presents a generalized reconstruction method
that corrects for nonrigid motion based on least-squares
estimation. The full least-squares estimate requires prohibitive computation, but suggests an efficient approximation
called the “refocusing reconstruction.” We consider the
conditions for this approximation, and propose efficient
algorithms for the refocusing reconstruction. We also present
a time-efficient approach to cardiac gating that improves the
performance of the refocusing reconstruction. These methods are applied to steady-state DWI (SS-DWI) (19), a potentially powerful DWI method that is particularly sensitive to
motion-induced phase corruptions.
THEORY
In this section, we examine the effect of nonlinear phase
corruption on DWI data and derive a nonlinear phase
correction based on least-squares estimation. Computational issues require that the least-squares reconstruction
be approximated, and an approximation called the “refocusing reconstruction” is presented. This reconstruction
can be performed in either image-space or k-space.
Least-Squares Reconstruction
Phase corruption can be expressed at each point in image
space as a multiplication of the magnetization by a phase
[1]
where d is a CR â«» 1 vector containing the k-space data, m
is a N 2 â«» 1 vector containing the object magnetization in
Cartesian coordinates, and P (CN 2 â«» N 2 ), F (CN 2 â«» CN 2 ),
and G (CR â«» CN 2 ) are matrices respectively representing
image-space phase corruption, discrete Fourier transform,
and resampling from a Cartesian grid onto the k-space
trajectory. The phase corruption matrix P is a stack of
diagonal matrices containing the phase corruption for each
readout. This phase corruption is measured by the navigator phase, which must be interpolated to match the
high-resolution readout (since the n â«» n navigator is of
lower resolution than the high-resolution readout). The
matrix R ⫽ GF describes the basis set used in the k-space
trajectory. Expressions for these vectors and matrices are
given in Appendix A.
The magnetization can be estimated from the data in Eq.
[1] using least-squares estimation. Assuming an accurate
measurement of the phase errors P, the least-squares reconstruction is:
m̂ ⫽ 共P*R*RP兲⫺1 P*R*d
⫽ M ⫺1 P*R*d
[2]
where * denotes the conjugate transpose. It is helpful to
consider this reconstruction as occurring in three steps: 1)
reconstruction R* of the individual data frames in imagespace, which we call the “interleave reconstruction”; 2)
phase correction P*, which we call the “refocusing operator”; and 3) removal of remaining artifacts by M⫺1 ⫽
(P*R*RP)⫺1, which we call the “unmixing operator.”
These steps are diagrammed in Fig. 2 using the notation
described in this section and in Appendix A.
The data are first partially reconstructed using the interleave reconstruction R* ⫽ F*G*, which separately resamples each readout onto a Cartesian grid at the highresolution readout size, and inverse Fourier transforms the
resampled data without combining the interleaves. The
resulting vector a ⫽ R*d contains the image-space data
corresponding to each readout of k-space data. Unless
each readout covers a full field of view (FOV), the images
in the vector a are aliased.
In the next step, the refocusing operator P* is applied to
the interleave images to yield the refocused image m̃ ⫽ P*a.
The refocusing operator performs the essential task of
rephasing the interleaves based on the phase corruption
measured by the navigator. Each interleave image is multiplied by the phase conjugate of the navigator image to
rephase the unaliased component. The sum of these
rephased interleaves is the refocused image. Although the
unaliased components will add coherently in this summation, the aliased components will not necessarily cancel,
since the orthogonality of the Fourier encoding has been
disturbed. The refocused image is thus the desired image
plus some aliased artifact. To the extent that the phase
Nonlinear Phase Correction for DWI
345
FIG. 2. Diagram of the least-squares matrix reconstruction given in Eq. [2]. Each raw data interleave dj is gridded and inverse Fouriertransformed by the interleave reconstruction to form an aliased image aj. These images are phase-corrected and summed by the refocusing
operator to form the refocused image m̃. Remaining aliased energy is removed by the unmixing operator M⫺1 to give the least-squares
optimal image m̂. If the aliased energy in the refocused image is sufficiently low, the reconstruction can stop at m̃, avoiding the
computationally expensive inverse calculation M⫺1. In this case, we call the algorithm a “refocusing reconstruction.”
errors are random, this aliased energy will be incoherent
and will tend to phase-cancel. The characterization and
minimization of the aliased signal are discussed later.
The final step of the least-squares reconstruction is the
application of the unmixing operator M⫺1 to the refocused
image. The unmixing operator removes two artifacts from
the refocused image: nonuniformities in the sampling trajectory, and motion-induced aliasing (see Appendix B).
The compensation for trajectory nonuniformities, a standard correction in the least-squares reconstruction of nonCartesian data, will be neglected here. More important to
our discussion is the removal of the aliased energy discussed above. The unmixing operator uses a linear combination of voxel signals to remove this coupling between
voxels. If the refocused image contains little aliasing, the
mixing matrix is close to the identity matrix.
Refocusing Reconstruction
Calculation of the unmixing operator M⫺1 is problematic
due to its computational complexity. Not only is the N 2 â«»
N 2 inversion computationally expensive, it also cannot be
precomputed because it depends on the phase corruptions. Further, the calculation of the mixing matrix M itself
requires CRN 2 ⬇ N 4 operations for each element of the
N 2 â«» N 2 matrix (see Eq. [12]). These computational requirements make direct calculation of the unmixing matrix (and hence the full least-squares reconstruction) intractable.
If the sources of coupling are sufficiently low, the mixing matrix has little off-diagonal energy and is close to the
identity matrix. Approximating the mixing matrix as the
identity yields the “refocusing reconstruction”:
m̃ ⫽ P*R*d.
[3]
Provided the unmixing matrix is nearly diagonal, the refocusing reconstruction is a good approximation to the
least-squares estimate that is much faster to compute. The
refocusing reconstruction is also conceptually elegant: before the readouts are combined, each readout is refocused
in image space by multiplication with the navigator phase
conjugate. This interpretation is not only computationally
useful, but also suggests a more efficient algorithm than
the matrix computation in Eq. [3] for refocusing the data.
This algorithm will be discussed later.
Some important properties of the refocusing reconstruction can be understood by considering the problem in
k-space. Corruption by a spatial phase function (e i␾(r)) is
equivalent to convolution by the transform of that phase
function (Ᏺ{e i␾(r)}) in the spectral domain. If all of k-space
were convolved with the same phase corruption kernel,
we could restore the original spectrum by deconvolving
with this kernel. However, in motion-corrupted DWI data,
each readout samples a spectrum that has been convolved
with a different kernel (Ᏺ{e i␾ j(r)} for the jth readout).
Interpreted in k-space, the refocusing operator deconvolves each interleave separately using the phase-corruption kernel from that excitation. Refocusing uses the correct deconvolution kernel for points from the same excitation. For data from different excitations, the
deconvolution will work well if the kernels are similar for
adjacent positions in k-space, but will be unable to remove
artifacts if the kernels are dissimilar. Hence, it is desirable
for the kernel to change slowly and smoothly over k-space
so that the deconvolution is locally consistent. An additional concern is periodicity in the phase corruptions that
causes a periodic weighting of the data in k-space. Any
residual periodicity left by the refocusing operator will
cause coherent ghosts, which are a more disruptive artifact
than the diffuse aliasing associated with discontinuities.
The refocusing reconstruction will produce high-fidelity
346
Miller and Pauly
(f) Multiply (e) by the phase conjugate of (c) (refocusing step).
(g) Add (f) to the sum accumulating in step 1.
3. The refocused image is the cumulative sum from (g).
For every interleave, this reconstruction requires two
Fourier transforms and a matrix multiplication at the highresolution matrix size (N â«» N). Each Fourier transform
requires time N 2 log N 2 , and the matrix multiplication
requires time N 3 for a total complexity of O(CN 2 log N 2 ⫹
CN 3 ) ⬇ O(CN 3 ).
The k-space refocusing algorithm (Fig. 3b) is:
FIG. 3. Block diagram of two versions of the refocusing reconstruction in which refocusing is performed in (a) image space and (b)
k-space. The boxed expressions indicate the operations of gridding
(Ᏻj grids the jth readout trajectory), inverse Fourier transform (Ᏺ⫺1),
and refocusing (multiplication by ᏼj*(r) ⫽ e⫺i␾j(r) in image-space or
convolution by Ᏺ{ᏼj* } in k-space).
images when phase corruptions are smooth and nonperiodic in k-space.
MATERIALS AND METHODS
Efficient Refocusing Reconstruction
Although the matrix formulation given in the previous
section is a useful way to derive the refocusing algorithm,
it is an inefficient method for calculating the refocused
image. Directly implementing the matrix formulation has
prohibitive memory requirements and is unlikely to take
full advantage of the highly structured nature of these
matrices. It is much more efficient to use operations equivalent to those given in the matrix formulation, since these
operations often have highly optimized algorithms (such
as the Fast Fourier Transform (FFT)). The implementation
of image-space refocusing shown in Fig. 3a can be calculated more efficiently than the matrix calculation of m̃
shown in Fig. 2. Refocusing can also be implemented as a
k-space deconvolution, as shown in Fig. 3b.
The image-space refocusing algorithm in Fig. 3a involves the following steps:
1. Initialize a zero-filled N â«» N matrix (where we will
accumulate the refocused image).
2. For each interleave:
(a) Reconstruct the navigator data in k-space (n â«»
n).
(b) Zero-pad (a) to the final image resolution (N â«»
N).
(c) Inverse Fourier transform (b) (to get the navigator
image).
(d) Reconstruct the high-resolution data frame in kspace at the final image resolution (N â«» N).
(e) Inverse Fourier transform (d) (to get the interleave
image).
1. Initialize a zero-filled N â«» N matrix (where we will
accumulate the refocused spectrum).
2. For each interleave:
(a) Reconstruct the navigator data in k-space (n â«»
n).
(b) Zero-pad (a) by a factor of 2 (2n â«» 2n).
(c) Inverse Fourier transform (b).
(d) Calculate the phase conjugate of (c).
(e) Fourier transform (d) (to get the refocusing kernel).
(f) Reconstruct the high-resolution data frame in kspace at the final image resolution (N â«» N).
(g) Convolve (f) with (e) (refocusing step).
(h) Add (g) to the sum accumulating in step 1.
3. Inverse Fourier transform the cumulative sum from
(h) to get the refocused image.
The refocusing kernel must be calculated at twice the
navigator bandwidth (step 2b) because normalizing in image-space (step 2d) extends the kernel in k-space. Calculating the kernel at twice the navigator bandwidth avoids
Gibbs ringing in the refocused image due to this extension.
For every interleave, the k-space algorithm calculates two
2n â«» 2n Fourier transforms and one convolution of an
N â«» N matrix with a 2n â«» 2n kernel. This algorithm runs
in approximately O(CN 2 n 2 ⫹ Cn 2 log n 2 ) ⬇ O(CN 2 n 2 )
time (where constant multipliers have been dropped, as is
standard in O-notation).
In most imaging situations, the k-space algorithm will
perform better than the image-space algorithm, since we
will usually have n 2 ⬍ N. However, this is not always the
case, as will be discussed later.
Cardiac Synchronization
If minor subject restraints are used to reduce bulk motion,
the primary source of motion artifacts in the brain is deformation related to the cardiac cycle (see Fig. 1). Cardiac
triggering information can be used to synchronize interleave ordering to the cardiac cycle to reduce discontinuities and periodicities in the k-space data (20,21). For example, in a spin-warp trajectory, phase encodes are normally gathered incrementally, resulting in a periodic
weighting over k-space with large discontinuities (see Fig.
4a). These effects can be lessened by choosing a phaseencode increment that traverses the extent of k-space during a single cardiac cycle. By interleaving the phase encodes over multiple cardiac cycles in this manner, adjacent lines in k-space are collected at the same portion of
the cardiac cycle, and the resultant k-space weighting is
smooth and nonperiodic (Fig. 4b). This method for cardiac
Nonlinear Phase Correction for DWI
347
a total scan time of 6:15 per image. To date, 10 subjects
have been scanned with variants of this sequence, five of
whom were scanned with the prescription specified above.
In some experiments, acquisition was synchronized to
the cardiac cycle using plethysmographic triggering. To
maintain the steady state, dummy cycles without acquisition were gathered while waiting for cardiac triggers. Following a systolic trigger, data were collected in an interleaved order such that the phase encodes acquired during
a single cardiac cycle were evenly distributed across k y space. This process was repeated for successive cardiac
cycles until all phase encodes had been acquired. Additionally, two to three dummy cycles (80 –120 ms) were
inserted following detection of a trigger to avoid the largest
motions associated with systole. Cardiac synchronization
usually increased scan times by 1–2 min.
Images were reconstructed using three methods: no correction, the standard linear correction, and the k-space
refocusing correction proposed in this work. Reconstruction of the data was implemented off-line in C code using
the FFTW library for Fourier transforms (22). Both the
image-space and k-space versions were implemented;
however, all refocused images shown here were calculated
with the k-space implementation.
FIG. 4. The effect of 2DFT phase-encode ordering on k-space. The
phase history of a representative voxel over time is shown in the middle
plot. The gray swathes represent cardiac-gated imaging windows. The
k-space modulation imparted by this phase history is shown for (a)
standard incremental phase-encode ordering and (b) an alternate
phase-encode ordering. The latter ordering interleaves the data from
different cardiac cycles to achieve a smoother modulation.
synchronization of phase encoding should not require
large increases in scan time, since data can be acquired
throughout a large portion of the cardiac cycle. Cardiac
synchronization in this manner reduces the types of kspace modulations that are most problematic for the refocusing reconstruction— discontinuity and periodicity—
without incurring the large increases in scan time that
usually accompany cardiac gating.
RESULTS
Refocusing Reconstruction
The different navigator corrections are compared on untriggered data in Fig. 6. Motion-induced phase artifacts
interfere to cause signal attenuation in areas with large
motions, and little or no artifact outside the object. Motion
artifacts are manifested as signal loss rather than the more
familiar ghosting artifacts due to the large number of excitations in the acquisition. All of the uncorrected images
(top row) are corrupted, with the most serious artifacts in
the S/I-weighted image. Arrows in the top row indicate
artifacted areas in the R/L- and A/P-weighted images. The
linear correction (middle row) restores signal in some corrupted regions (e.g., posterior areas in the S/I-weighted
image), but the images remain attenuated in the medial
Experiments
Brain images were acquired on healthy volunteers with the
navigated steady-state diffusion-weighted imaging (SSDWI) pulse sequence shown in Fig. 5. To minimize large
motions, padding was placed at the subject’s temples. The
readout consisted of a spiral navigator (FOV ⫽ 24 cm,
matrix ⫽ 16 ⫻ 16, in-plane resolution ⫽ 15 ⫻ 15 mm2,
receive bandwidth (BW) ⫽ 62.5 kHz) followed by a spinwarp high-resolution acquisition (FOV ⫽ 24 cm, matrix ⫽
192 ⫻ 192, in-plane resolution ⫽ 1.25 ⫻ 1.25 mm2, receive
BW ⫽ 7.8125 kHz). Experiments were performed on a
1.5 T GE Signa CV/i research scanner (General Electric Co.,
Milwaukee, WI) with 40 mT/m gradients switchable at
150 mT/m/ms and a transmit/receive quadrature birdcage
head coil. DW images were acquired using diffusion gradients with G ⫽ 40 mT/m and ␶ ⫽ 6.5 ms (b eff ⫽ 980 s/
mm2 (19)), which could be applied along any direction.
Some studies also acquired an image with minimal diffusion weighting (G ⫽ 0.3 mT/m and ␶ ⫽ 6.5 ms). Other
imaging parameters were: TR ⫽ 40 ms, ␣ ⫽ 30°, and
thickness ⫽ 6 mm. The sequence was repeated 48 times for
FIG. 5. Navigated steady-state DWI pulse sequence used in this
study. The sequence is fully refocused except for the diffusion
gradient (with amplitude G and duration ␶), which rephases DW
echoes over multiple TRs.
348
Miller and Pauly
FIG. 6. Untriggered spiral-navigated
2DFT images with varying levels of correction. Rows (from top to bottom) have
no correction, the standard linear correction, and the refocusing correction.
Columns (from left to right) show images
with diffusion weighting along the A/P,
R/L, and S/I directions. Arrows in the
A/P- and R/L-weighted uncorrupted images (top row) indicate examples of
dephasing due to motion. Arrows in the
middle row indicate artifacts introduced
by the linear correction. Arrows in the
bottom row indicate structures with increased conspicuity after the refocusing
correction.
areas where motion is more severe. Additionally, the linear correction introduces attenuation to the frontal lobes in
the R/L- and A/P-weighted images due to a poor linear fit
to the nonlinear phase artifacts (arrows in middle row of
Fig. 6). The refocused images (bottom row) have restored
signal and minimal artifacts. In particular, the medial regions show structures that were partially or completely
suppressed in the uncorrected images (e.g., the external
and internal capsule, indicated by arrows). The images
weighted along the primary axis of motion (S/I) have a
lower signal-to-noise ratio (SNR) due to the large motions
in the S/I direction. In peak systole, these motions cause
intravoxel dephasing, which cannot be corrected by postprocessing.
FIG. 7. Cardiac-synchronized, spiral-navigated 2DFT images acquired
with the phase-encode interleaving
scheme shown in Fig. 4. Diffusion
weighting is along the (a) A/P, (b) R/L,
and (c) S/I directions.
Cardiac Triggering
A set of refocused images acquired with cardiac synchronization is shown in Fig. 7. These images show strong
diffusion contrast with clear delineation of small whitematter structures. The S/I-weighted image does not suffer
from the degraded SNR seen in untriggered images. The
A/P- and R/L-weighted images are also improved with
triggering, and show fine detail (such as cortical white
matter) that is less apparent in the untriggered images.
Diffusion Tensor Imaging
A full set of diffusion tensor images acquired with the
standard directions for DTI (23) is shown in Fig. 8. These
Nonlinear Phase Correction for DWI
349
FIG. 8. Navigated SS-DWI images
acquired with diffusion weighting
along the standard directions for diffusion tensor. Images are cardiacsynchronized and were reconstructed using the refocusing reconstruction.
images were acquired with cardiac-synchronized phaseencode ordering, and reconstructed using the refocusing
reconstruction. Intravoxel dephasing is minor in these images due to the systolic delay (120 ms) and the fact that
diffusion weighting was never applied along the principal
axis of brain motion (S/I). These images show strong diffusion weighting with little visible artifact. Figure 9 shows
the trace image and fractional anisotropy (FA) map for
these data before and after correction. The correction improves both images, although the most dramatic improvement is found in the FA map.
DISCUSSION
We have presented a method for correction of nonlinear
phase corruption. This method does not require additional
FIG. 9. Trace and FA maps for the data in Fig.
8 before and after the refocusing correction. Images
before and after correction are identically windowed.
The trace image shows a moderate improvement in
contrast. The FA map shows a more marked improvement. The uncorrected FA map shows apparent anisotropy in gray-matter regions due to motioninduced signal loss in a subset of the images. The
refocusing correction restores this signal to give a
more accurate map of anisotropy.
350
Miller and Pauly
navigator data beyond those used for 2D linear corrections.
Rather, the refocusing correction represents a more complete use of the information in the navigator. The refocusing reconstruction algorithm is straightforward to implement: refocusing is simply a k-space deconvolution that is
performed before individual readouts are combined. Alternatively, it can be thought of as image-space multiplication of individual data frames by the navigator phase conjugate.
Generality of Reconstruction
The least-squares reconstruction is general in the sense
that its derivation does not make assumptions about the
trajectory or pulse sequence used to gather the data. Although the refocusing reconstruction will not always be a
good approximation of the least-squares solution, it is
expected to perform well provided the k-space modulation
induced by phase corruptions has minimal discontinuity
and periodicity. In addition to the 2DFT data shown here,
we previously used the refocusing reconstruction for selfnavigated spiral data with similar results (24). Pipe and
colleagues (14) have also had good success using a refocusing reconstruction with an FSE PROPELLER sequence
(however, in their case, the refocusing reconstruction was
the same as the least-squares reconstruction, since each
readout was unaliased).
While full least-squares reconstruction is generally intractable due to the inversion of the mixing matrix, it may
be manageable in special cases or with efficient approximations. For example, the mixing matrix M will be block
diagonal and sparse for trajectories in which the signal
from a voxel can only alias onto a subset of the other
voxels (see Appendix B). In this case, the matrix inversion
of M can be calculated by individually inverting the
smaller block submatrices. In our experience, this calculation is still very time-consuming. Alternatively, the mixing matrix could be approximated with an iterative
scheme similar to that proposed for multicoil reconstructions (25).
Fidelity of the Refocusing Reconstruction
The impulse response of the mixing operator is a useful
way to characterize the fidelity of the refocused image. The
off-peak power in the impulse response indicates the
amount of aliased energy in the refocused image, and
therefore the relative benefit of applying the unmixing
operator. Of particular interest are coherences in the voxel
impulse responses, which indicate coherent ghosts in the
refocused image. The impulse response of the mixing operator is shown in Fig. 10a for several representative voxels in an ungated data set. Although there is significant
energy off-center, most of the aliased signal has random
phase and will add incoherently. The impulse responses
also show widening of the central peak and coherent offcenter peaks, which will lead to blurring and ghosting,
respectively. These effects in the impulse response must
be reduced in order for the refocusing operator to provide
accurate reconstructions.
The quality of the refocusing reconstruction is dependent on the number of excitations that make up an image.
Aliased energy in the refocused image decreases as the
FIG. 10. Impulse responses of the mixing operator for a 2DFT
acquisition with (a) no averaging and (b) four averages. Each subfigure overlays the impulse responses for the same 10 voxels taken
from a representative navigator data set. The arrow in a indicates
coherences in the mixing impulse responses that lead to coherent
ghosts. Averaging in b suppresses much of the off-peak energy.
number of excitations increases, since the incoherent
sources of aliasing are sampled more often. As shown in
Fig. 10b, averaging suppresses the incoherent off-peak energy found in unaveraged data. In this case, averaging also
suppresses the coherent peaks since data were acquired
asynchronous to the cardiac cycle. Some averaging is always necessary to correct strong phase corruptions, which
warp the k-space trajectory. Acquiring multiple averages
ensures that the perturbed encoding functions will cover
k-space more fully. One side effect of this need to average
is the improvement of the refocusing reconstruction.
Cardiac Synchronization
The need for cardiac gating has been a common criticism
of navigator methods. Since linear corrections assume
Nonlinear Phase Correction for DWI
351
Computational Efficiency
FIG. 11. Simulated impulse response magnitude of the mixing operator for a 2DFT acquisition with and without cardiac synchronization (nav ⫽ 1). a: The impulse responses for standard incremental
order have strong coherences (arrows) that will lead to ghosting. b:
These coherences are eliminated when cardiac-synchronized
phase-encode ordering is used.
rigid-body motion, cardiac gating is required to limit
imaging to the more stable portions of the cardiac cycle,
in which motions are better approximated as linear.
With the nonlinear corrections introduced here, postsystolic delays only have to be long enough to avoid
intravoxel dephasing during peak systole, which cannot
be corrected in postprocessing. For the 2DFT SS-DWI
sequence used here, cardiac synchronization requires
only moderate increases in scan time (25–33%). This
imaging efficiency is possible because cardiac information is being used to synchronize the acquisition to the
cardiac cycle, rather than simply limiting the portion of
the cardiac cycle in which data are gathered. Our current implementation approximates the heart rate for the
entire scan. Monitoring changes in heart rate during
scanning would further reduce the overhead time for
cardiac synchronization.
Improvements to the refocused image due to cardiac
synchronization can be seen by examining the impulse
response of the mixing operator. The impulse response for
the standard phase-encode ordering scheme (Fig. 11a) has
significant off-peak energy, as well as a series of coherent
peaks (arrows) due to the periodic weighting of k-space
introduced by the cardiac cycle (several harmonics are
clearly visible). Cardiac-synchronized ordering (Fig. 11b)
eliminates off-center coherent peaks and produces a more
well-behaved impulse response. The attenuation band
around the impulse is due to the suppression of low frequencies in the k-space modulation.
Even though they are functionally equivalent, the imagespace and k-space refocusing reconstructions differ in
terms of computational efficiency. The choice of algorithm
depends on the relative sizes of the navigator and highresolution acquisitions, as well as the details of the highresolution trajectory. However, in almost all practical imaging scenarios, the k-space implementation is expected to
outperform the image-space algorithm. As presented
above, the k-space algorithm will perform better if n 2 ⬍ N,
because the k-space convolution will be faster than the
large matrix multiply.
An additional consideration in the comparison of the
two methods is the specific k-space trajectory used for the
high-resolution acquisition. Our previous analysis of the
k-space algorithm assumed that the convolution must be
performed over the entire N â«» N matrix, which is not true
for all trajectories. For some trajectories, the convolution
time can be significantly reduced by performing a local
convolution. For example, with a spinwarp trajectory, the
convolution need only include the area around each
phase-encode line, and can therefore be performed in
O(Nn 2 ) time. For this trajectory, the k-space algorithm
will in general perform better when n ⬍ N, which should
always be the case.
Thus far we have only considered 2D acquisitions. For
3D acquisitions, all 2D FFTs, matrix multiplications, and
convolutions are replaced with 3D operations. The k-space
algorithm is expected to scale more favorably than the
image-space algorithm for 3D acquisitions, since a local 3D
convolution can be performed significantly faster than
large 3D Fourier transforms and matrix multiplications
can be calculated. Early work with 3D SS-DWI knee imaging used a k-space algorithm with good results (26). For
this application, the image-space algorithm was found to
be prohibitively slow.
For the experiments performed here (a spinwarp trajectory with n ⫽ 16, N ⫽ 192), the k-space algorithm would
be expected to drastically outperform the image-space algorithm. Our implementation computed the k-space reconstruction in 2.6 s per average, and the image-space
reconstruction in 6.9 s per average. The modest improvement of the k-space algorithm is likely due to the use of
highly-optimized FFT software and unoptimized convolution calculations.
Weighted Refocusing Reconstruction
Previous work on navigated DWI usually included a
weighting or thresholding of data frames that are suspected to contain uncorrectable phase artifacts. The
method proposed here does not include such a correction,
primarily because it has not proven necessary under the
current imaging circumstances. However, imaging methods that acquire fewer data frames are more impacted by
artifacted data frames, and might benefit from a weighting
scheme.
Examining Fig. 1, we see that during periods of the most
intense motion (such as systole), navigator images exhibit
signal loss due to intravoxel dephasing. Navigator magnitude can thus be used as an indicator of uncorrectable
phase corruption and can be incorporated into the recon-
352
Miller and Pauly
struction. If we weight each image-space sample by its
navigator magnitude prior to correction, the least-squares
estimate will be biased to more reliable data. This can be
accomplished by solving the preconditioned problem:
Qd ⫽ QRPm
[4]
where Q is a preconditioning matrix that accomplishes the
desired image-space weighting. Note that preconditioning
is applied to the problem in k-space.
The least-squares solution to Eq. [4] is:
m̂ Q ⫽ 共P*R*Q*QRP兲⫺1 P*R*Q*Qd.
[5]
The estimate in Eq. [5] filters the k-space data by Q*Q prior
to refocusing. The magnitude effects of the filter are then
removed by the unmixing matrix (P*R*Q*QRP)⫺1. To
weight each image-space voxel by the navigator magnitude, Q*Q must convolve each high-resolution interleave
with the Fourier transform of the navigator magnitude. In
other words, R*Q*Qd should contain the interleave images weighted by the navigator magnitude.
A weighted refocusing reconstruction can then be defined by approximating the mixing matrix P*R*Q*QRP as
diagonal by zeroing the off-diagonal elements. The approximated unmixing matrix removes the weighting caused by
the preconditioning, but does not attempt to remove any
aliasing. Regularization can be used to improve the stability of the inversion in noise regions with low signal.
single image at roughly the same portion of the cardiac
cycle. The cardiac cycle could then be filled out by acquiring more slices (in multislice), more z-phase encodes (for
3D), or more averages. All of these methods use cardiac
information to achieve well-behaved modulations of the
k-space data without requiring large increases in scan
time.
CONCLUSIONS
We have presented a method for correcting nonlinear
phase errors in diffusion-weighted imaging using 2D navigation. The refocusing reconstruction, which is an approximation to the computationally-expensive leastsquares estimate, is a straightforward algorithm that makes
complete use of the information in a 2D navigator. The
heart of the reconstruction is the refocusing operator,
which rephases the unaliased component of each acquisition via multiplication with the navigator phase conjugate.
This reconstruction is a good approximation to the full
least-squares estimate provided the phase errors cause
smooth, nonperiodic modulations of k-space, and a sufficient number of excitations are acquired. Cardiac triggering information can be used to achieve desirable modulations through ordering of data interleaves. The refocusing
reconstruction is expected to be applicable in a wide range
of imaging conditions.
APPENDIX A
Application to Other Sequences and Trajectories
Structure of Reconstruction Matrices
The refocusing reconstruction will provide a good reconstruction of the data provided 1) a large number of excitations are acquired, and 2) k-space modulations induced by
motion are smooth and nonperiodic.
The requirement for many excitations minimizes aliasing energy by allowing the incoherent phase errors to
cancel. We have already discussed the improvement in the
refocused image from averaging; other methods of increasing the number of excitations will have the same effect. For
example, a 3D acquisition should be better-conditioned to
a refocusing reconstruction than a 2D multislice acquisition. Simply increasing the number of interleaves in an
acquisition will also improve the refocused image. The
need for a large number of excitations is a disadvantage for
some sequences, such as spin echo, which require long
echo times. These methods may benefit from a weighted
refocusing reconstruction, which will tend to suppress
these errors, particularly if regularization is used.
The requirement for smooth and nonperiodic k-space
modulations ensures that the refocusing operator will be
able to remove phase corruptions. This requirement can be
addressed in most trajectories with cardiac synchronization. Some trajectories cannot simultaneously achieve
smoothness and nonperiodicity, in which case we prefer
to meet the smoothness constraint since the refocusing
reconstruction does a poor job of removing discontinuities. Spiral and radial scans can achieve a smooth k-space
weighting by distributing the rotated interleaves from 0° to
180° over each cardiac cycle (this scheme would result in
azimuthal periodicity in k-space). Segmented EPI scans
can also be made smooth if we acquire each segment of a
The motion-corrupted magnetization during the cth readout is denoted as m c (r) ⫽ 兩m c (r)兩e i␾ c(r). The phase corruption during the cth excitation is expressed by the N 2 ⫻ N 2
matrix Pc:
Pc ⫽
冢
ei␾ c共r1兲
0
···
0
0
·
·
·
0
ei␾ c共r2兲
·
·
·
0
···
0
·
·
·
···
冣
.
ei␾ c共rN 兲
2
[6]
Pcm is the phase-corrupted magnetization that is the
source of signal during the cth readout. Each of the C
interleaves has a different phase corruption, represented
by stacking the Pc’s to form P:
冢冣
P1
P⫽
P2
·
·
·
PC
.
[7]
P splits the magnetization vector m into C separate phasecorrupted magnetization vectors. These vectors are
stacked to form the CN 2 â«» 1 vector Pm representing the
phase-corrupted magnetization that is sampled over the
course of the experiment.
Each of the phase-corrupted magnetization vectors is
separately Fourier transformed by the matrix F (CN 2 â«»
CN 2 ):
Nonlinear Phase Correction for DWI
F⫽
冢
FN
0
···
0
·
·
·
0
N
···
F
·
·
·
0
353
0
0
·
·
·
· · · FN
冣
REFERENCES
.
[8]
F is composed of C 2DFT matrices FN (N 2 â«» N 2 ) consisting of the Fourier kernels:
N
F lm
⫽ ei2␲kl䡠rm.
[9]
The product FPm is the CN 2 â«» 1 vector representing the C
phase-corrupted spectra, sampled in Cartesian coordinates.
Finally, the spectra given by FPm are resampled onto the
k-space sampling trajectories by the resampling matrix G
(CR â«» CN 2 ):
G⫽
冢
G1
0
0
·
·
·
0
G2 · · · 0
·
·
·
·
·
·
0 · · · GC
···
0
冣
.
[10]
Each submatrix Gc (R â«» N 2 ) interpolates the Cartesian
k-space points kl onto the sampling points for the cth
c
readout sm
:
c
⫽ sinc共kl ⫺ scm 兲.
G lm
[11]
Note that for trajectories that collect data on Cartesian
coordinates (such as spinwarp or interleaved EPI), the
elements of G are discrete delta functions.
APPENDIX B
Mixing Matrix
Using the notation of Appendix A, the elements of the
mixing matrix M for an arbitrary trajectory are:
冘
C
M lm ⫽
c⫽1
i共␾ c 共rm 兲⫺␾ c 共rl 兲兲
e
冘冋冘
R
N2
r⫽1
n⫽1
sinc共kn ⫺ src 兲e⫹i2␲kn䡠rm
冋冘
N2
â«»
n⫽1
册
册
sinc共kn ⫺ src 兲e⫺i2␲kn䡠rl . [12]
The summation over each readout r describes non-idealities in the k-space trajectory that causes coupling between
voxels within a readout. The phase difference describes
the aliasing caused by motion-induced phase errors. These
phase offsets are manifested as perturbations of the encoding functions used in the trajectory.
ACKNOWLEDGMENTS
K.M. thanks Dr. Brian Hargreaves for useful discussions
regarding the method and manuscript.
1. Anderson AW, Gore JC. Analysis and correction of motion artifacts in
diffusion weighted imaging. Magn Reson Med 1994;32:379 –387.
2. Ordidge RJ, Helpern JA, Qing ZX, Knight RA, Nagesh V. Correction of
motion artifacts in diffusion-weighted MR images using navigator echoes. Magn Reson Imaging 1994;12:455– 460.
3. deCrespigny AJ, Marks MP, Enzmann DR, Moseley ME. Navigated
diffusion imaging of normal and ischemic human brain. Magn Reson
Med 1995;33:720 –728.
4. Butts K, Crespigny A, Pauly JM, Moseley M. Diffusion-weighted interleaved echo-planar imaging with a pair of orthogonal navigator echoes.
Magn Reson Med 1996;35:763–770.
5. Mori S, van Zijl PCM. A motion correction scheme by twin-echo navigation for diffusion-weighted magnetic resonance imaging with multiple RF echo acquisition. Magn Reson Med 1998;40:511–516.
6. Williams CFM, Redpath TW, Norris DG. A novel fast split-echo multishot diffusion-weighted MRI method using navigator echoes. Magn
Reson Med 1999;41:734 –742.
7. Brockstedt S, Moore JR, Thomsen C, Holtas S, Stahlberg F. Highresolution diffusion imaging using phase-corrected segmented echoplanar imaging. Magn Reson Imaging 2000;18:649 – 657.
8. Clark CA, Barker GJ, Tofts PS. Improved reduction of motion artifacts in
diffusion imaging using navigator echoes and velocity compensation. J
Magn Reson 2000;142:358 –363.
9. Norris DG, Driesel W. Online motion correction for diffusion-weighted
imaging using navigator echoes: application to RARE imaging without
sensitivity loss. Magn Reson Med 2001;45:729 –733.
10. Bosak E, Harvey PR. Navigator motion correction of diffusion weighted
3D SSFP imaging. MAGMA 2001;12:167–176.
11. Jiang H, Golay X, van Zijl PCM, Mori S. Origin and minimization of
residual motion-related artifacts in navigator-corrected segmented diffusion-weighted EPI of the human brain. Magn Reson Med 2002;47:
818 – 822.
12. Butts K, Pauly J, Crespigny A, Moseley M. Isotropic diffusion-weighted
and spiral-navigated interleaved EPI for routine imaging of acute
stroke. Magn Reson Med 1997;38:741–749.
13. Atkinson D, Porter DA, Hill DLG, Calamante F, Connelly A. Sampling
and reconstruction effects due to motion in diffusion-weighted interleaved echo planar imaging. Magn Reson Med 2000;44:101–109.
14. Pipe JG, Farthing VG, Forbes KP. Multishot diffusion-weighted FSE
using PROPELLER MRI. Magn Reson Med 2002;47:42–52.
15. Norris DG. Implications of bulk motion for diffusion-weighted imaging
experiments: effects, mechanisms and solutions. J Magn Reson Imaging
2001;13:486 – 495.
16. Enzmann DR, Pelc NJ. Brain motion: measurement with phase-contrast
MR imaging. Radiology 1992;185:653– 660.
17. Greitz D, Wirestam R, Franck A, Nordell B, Thomsen C, Stahlberg F.
Pulsatile brain movement and associated hydrodynamics studied by
magnetic resonance phase imaging. Neuroradiology 1992;34:370 –380.
18. Poncelet BP, Wedeen VJ, Weisskoff RM, Cohen MS. Brain parenchyma
motion: measurement with cine echo-planar MR imaging. Radiology
1992;185:645– 651.
19. Buxton RB. The diffusion sensitivity of fast steady-state free precession
imaging. Magn Reson Med 1993;29:235–243.
20. Bailes DR, Gilderdale DJ, Bydder GM, Collins AG, Firmin DN. Respiratory ordered phase encoding (ROPE): a method for reducing respiratory motion artifacts in MR imaging. J Comput Assist Tomogr 1985;9:
835.
21. Haacke EM, Patrick JL. Reducing motion artifacts in two-dimensional
Fourier transform imaging. Magn Reson Imaging 1986;4:359 –376.
22. Frigo M, Johnson SG. FFTW: an adaptive software architecture for the
FFT. Proc ICASSP 1998;3:1381–1384.
23. Basser PJ, Pierpaoli C. A simplified method to measure the diffusion
tensor from seven MR images. Magn Reson Med 1998;39:928 –934.
24. Miller KL, Meyer CH, Pauly JM. Self-navigated spirals for high resolution steady-state diffusion imaging. In: Proceedings of the 10th Annual
Meeting of ISMRM, Honolulu, 2002.
25. Preussmann KP, Weiger M, Bornert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med
2001;46:638 – 651.
26. Miller KL, Hargreaves BA, Gold GE, Pauly JM. Navigated 3D steadystate diffusion imaging of knee cartilage. In: Proceedings of the 11th
Annual Meeting of ISMRM, Toronto, 2003.

Purchase answer to see full
attachment

  
error: Content is protected !!