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