Augmented Generalized SENSE Reconstruction to Correct for Rigid Body Motion

THEORY

For a better understanding of how motion affects k-space and thus the reconstructed image, we will first describe the formation of an MRI data sample at a particular point in time, t_i (i = 1…κ), using discrete sampling in the time and spatial domain, r_u (u = 1…N²):

Therefore, the entire acquisition of κ k-space points can be stored in a measurement vector, m = [m(t₁) m(t₂)…m(t_κ)]^T, so that the signal formation in matrix form can be expressed as:

where v = [ρ(r₁)ρ(r₂)…ρ(r_N²)]^T is the spatial signal distribution in vector form, k_i = [k_x(t_i)k_y(t_i)k_z(t_i)]^T|_i=1…κ is an arbitrary k-space sampling trajectory, and r_u = [r_x,ur_y,ur_z,u]^T|_u=1…N² is the position vector. One can then obtain the solution vector v for the image signal distribution by computing the Moore-Penrose inverse of E:

The encoding matrix E can be replaced by the matrix product GF so that Eq. [2] becomes

Here, instead of populating the matrix E with N² × κ complex exponentials and performing a direct computation of Eq. [2], an FFT F (C_N²×_N²) is usually performed followed by a resampling of the Cartesian k-space data along the arbitrary k-space waveform using the matrix G (i.e., the inverse gridding matrix). For optimal sinc interpolation, one element in the inverse gridding matrix becomes G_l,i = sinc(k̃_l − k_i), where k̃_l and k_i represent the k-space grid locations. These positions, k_i, are defined by the arbitrary sampling trajectory. To boost computation efficiency, more compact resampling kernels (e.g., Kaiser-Bessel window) are normally used instead of sinc interpolation (19). The formalism described in Eq. [4] is also known as the inverse gridding operation (20). From a Cartesian image, inverse gridding is able to generate k-space sample points along an arbitrary trajectory defined by k_i, and therefore it represents the inverse operation to conventional gridding reconstruction (19).

Most frequently, MR data that are sampled along non-Cartesian k-space trajectories are reconstructed using gridding reconstruction (19); however, Eq. [3] can also be used to form an image from these arbitrarily sampled k-space data. Because of computational complexity, the direct inversion of Eq. [3] to find v is prohibitively inefficient; therefore, a least-squares solution for v,

is normally found in an iterative fashion by means of the conjugate-gradient method (21). This approach has recently gained popularity because it provides the foundation for iterative SENSE reconstruction (22) and nonlinear phase correction in diffusion-weighted imaging (DWI) (23).

In the next step, we demonstrate how motion impacts the acquired k-space data. Specifically, we assume that motion occurs between interleaves but not during interleaves. This assumption is legitimate for all MR applications for which the positional changes during one interleave are negligible. With FF^H = I, where I is the identity matrix, let Ω be a warping matrix that operates on an unperturbed image v₀ that yields the image v, and let Λ be the transformation rule in k-space that corresponds to Ω in the image domain, so that

Consequently, in the presence of motion, the signal equation for the acquired MR data changes to

Equation [7] clearly demonstrates that the solution of Eq. [4] yields a solution image v that is corrupted by motion. To obtain the true image v₀, the spatial warping Ω must be included in the encoding matrix E. An augmented version of Eq. [3] that considers spatial transformation is therefore

The term in parentheses on the left-hand side of Eq. [8] reflects the gridding of the acquired data, correction for motion, and inverse FFT. Equation [8] can be solved for v₀ following the same scheme as outlined in Eqs. [3] and [5]. Here, D(D_(γκ),_(γκ) = k_density(k_i)⁻¹) is introduced in Eq. [8] as a sampling density correction matrix exactly as in Pruessmann et al.’s (22) original approach. Also, the rotation of k-space data will affect the sampling density in k-space. Therefore, in this study the true sampling density was determined by a Voronoi tessellation approach (20) after motion correction was applied to the trajectories. Note, however, that accurate knowledge of the sampling density distribution is not so critical for this estimation problem because the iterative conjugate-gradient approach works even without prior knowledge of the sampling density information (22). Nevertheless, the exclusion of the density matrix causes much slower convergence rates because the density correction acts as a preconditioner. Additional preconditioning can be achieved by adding intensity correction, as suggested by Eq. [24] in Ref. 22.

Despite density correction, severe motion can lead to inconsistencies of k-space data (i.e., some areas in k-space are less densely sampled than required by the Nyquist sampling criterion). In other words, depending on the severity of motion, the linear dependency in Eq. [8] increases and therefore the conditioning of the system matrix E′ worsens, yielding considerable reconstruction artifacts. We hypothesize that parallel imaging helps to ameliorate the conditioning of this matrix system by adding additional encoding information using spatially varying coil sensitivities. Following Pruessmann et al.’s (22) seminal work, the k-space data for nc coils can be written as

where $E_{(i,\gamma )u}^{\prime }=c^{\gamma }({\mathbf{r}}_u)·exp(-\text{j}{\mathbf{k}}_i·{\mathbf{r}}_u)$, c^γ(r_u) is the complex coil sensitivity of the γ^-th coil at position r_u, and m₁ … m_nc are the k-space sampling data vectors for each of the nc coils. For each coil, one can combine E with c^γ on a per-pixel basis using E · diag(c^γ), and v can be estimated based on the aforementioned iterative conjugate-gradient method (22). For increased processing speed and without loss of generality, E is again replaced by GF. Instead of the conventional gridding and inverse gridding, we use in our implementation the convolution approach introduced by Wajer and Pruessmann (24). Alternatively, the “fastest Fourier transform in the West” (FFTW) algorithm (25) could be used for non-power-of-2 matrices and non-integer oversampling factors to speed up the gridding and inverse gridding steps during the conjugate-gradient iteration (26).

Assuming that there is an inverse warping operation (i.e., there exists a matrix Ω_inv so that Ω Ω_inv s₀ = s₀ (27)), our motion-correction approach in Eq. [7] can be expanded to incorporate the coil sensitivity information by rewriting the signal equation of the γ^th coil as

so that when all coils are combined, Eq. [10] can be expanded to

In this context, it is important to understand that in the case of patient motion the conventional frame of reference defined by the imaging gradients and the patient frame of reference are not necessarily aligned. Here the matrices G and F and the vector c^γ are all defined in the gradient frame of reference, whereas the final image v₀ is given in the patient frame of reference. Specifically, the matrix Ω transforms voxels from the patient frame into the gradient frame of reference. Notice that from the perspective of the patient frame of reference for each shot, the coil sensitivities do not remain static and instead “move around the object” (Fig. 2). To transfer the “coil motion” into the patient frame, the warping Ω_inv c^γ is used. In this study we assumed that the coil elements are static (because we used a rigid eight-channel head array coil). If needed, any true coil motion in the gradient frame of reference (theoretically, even without subject motion) can be modeled by including the motion pattern into the warping matrix Ω_inv. For example, coil motion relative to the patient might occur in cardiac or abdominal imaging during respiratory excursions. However, for self-calibrated data, such as variable-density spirals (Fig. 3a) or a preceding spiral-IN navigator (Fig. 3b), from which the coil sensitivity information can be extracted for each shot, the spatial warping of coil sensitivities is not required since these sensitivity maps are inherently registered to the trajectories to form the desired image.

In this study we assumed rigid body motion only, which can be modeled by translations and rotations. Specifically, one can apply rotations most efficiently during gridding by applying the appropriate rotation matrix to the initially prescribed k-space trajectory vector, k′ = R_φ,ι,ξ,n k, for each profile/interleave, n, where

is a 3D rotation around the x, y, and z axes using the angles φ_n, ϕ_n, and ξ_n, respectively. For simple rotation, the expression GΛ can be replaced by the gridding matrix G′, which uses k′ instead of k to compute the elements of G′. The direct rotation of the trajectory avoids interpolation errors that would otherwise occur using image (in our case k-space) rotation functions. In a similar fashion, one can achieve subvoxel translational motion, δx = (δx δy δz)^T, most efficiently by adding corresponding phase ramps in k-space during gridding by phasing the sampled k-space data:

Theoretically, the order in which rotational and translational motion is applied is irrelevant as long as it follows the order in which these parameters are acquired during the registration process (see next section). For implementation reasons, in this work we applied the rotation first, followed by the translation.

MATERIALS AND METHODS

RESULTS

DISCUSSION

Patient motion can significantly reduce the diagnostic quality of MR studies. Often the condition of the patient precludes him or her from being cooperative or remaining still for the duration of the examination. This is especially true for longer examination periods, and the likelihood that the patient will become uncomfortable and begin to move increases as the scan progresses. This in turn can require repeated examinations and lead to overly prolonged studies. In the pediatric setting, often sedation or anesthesia is needed to avoid motion artifacts, which leads to increased logistic problems and higher examination costs. It is therefore of great relevance to reduce the degree of motion artifacts in MRI scans. In the past several years, several attempts have been made to minimize artifacts from patient motion. These correction approaches have been carried out either prospectively or retrospectively (6,7,42,43).

In this study we propose a variant of parallel imaging reconstruction for arbitrary k-space acquisitions that allows for efficient retrospective in-plane motion correction. An excellent review of the potential of parallel imaging to correct for motion was presented recently by Larkman et al. (1). The current work is somewhat similar to the studies of Atkinson and Hill (4), Batchelor et al. (27), and Pipe et al. (2). The novelty of our approach is that it provides 1) a framework for motion correction that can be applied to arbitrarily sampled k-space data; 2) the capacity to overcome local k-space undersampling and different sampling densities by improving the conditioning of the inverse problem by means of additional spatial encoding information from coil sensitivities; 3) a rapid translation correction by the proper phasing of each k-space sample point and rotation correction by applying a rotation matrix to each k-space vector k(t_i) that is part of the gridding step; and 4) the extraction of motion information from either a navigator scheme built into the trajectory itself or an additional navigator image. This new method relies heavily on the pioneering work of Pruessmann et al. (22) and Kannengiesser et al. (44) regarding the fast reconstruction of undersampled, non-Cartesian k-space data. A related iterative approach for motion correction was recently proposed by Batchelor et al. (27), although they did not capitalize on the additional encoding capabilities from coil sensitivities to improve the conditioning of the inverse problem. Potential future developments include finding alternative ways to reconstruct the data by synthesizing missing k-space parts using information from adjacent points in k-space, such as GRAPPA-based methods (17,45).

A comparative evaluation of image quality for motion-corrupted data before and after correction for rigid body motion was performed both in vivo and using computer simulations. The proposed rigid body motion correction (i.e., correcting k-space data for rotation and translation) with conventional gridding reconstruction was able to reduce artifacts to some extent; however, considerable residual artifacts were still apparent, due primarily to local violation of the Nyquist sampling criterion in k-space. With the use of augmented iterative parallel imaging reconstruction, the level of residual artifacts was significantly reduced. The amount of motion simulated in this study was fairly large (±30° rotation and ±15 mm translation), but was selected to cover the typical amounts of motion encountered in clinical studies.

Exact knowledge of the rotation and translation components is crucial for any motion-correction scheme to be able to fully correct image artifacts. Like many other correction schemes, the efficacy of our method strongly depends on the ability of the navigator image and the coregistration technique to pick up on small amounts of motion in vivo. However, for the navigated interleaved spiral trajectories that we tested in vivo, the small residual motion errors did not cause any significant residual artifacts. If the residual motion error increases, mismatch will first be noticeable as a slight blurring in the corrected image and/or by “swirl artifacts.” Conversely, the great advantage of low-resolution navigator images over higher resolution images is that the coregistration is much more robust (i.e., the parameter search is not so easily trapped in the local optima of the cost function, and the smaller image size allows very rapid registration). In assessing the influence of navigator resolution on the registration quality, we found that the fluctuations resulting from the registration were in the submillimeter range. A potential difficulty for the registration step is the variable coil sensitivity modulation that affects the image intensity. These additional signal intensity modulations may reduce the effectiveness of traditional similarity measures, and this effect needs to be investigated further. To minimize this confounder, we first intensity-corrected the navigators. In future studies, similarity measures that are known to be more insensitive to differences in image contrast, such as (normalized) mutual information (30,31,37), may be more suitable than the Pearson correlation coefficient. However, in this study we did not find that the choice of the similarity measure had any influence on the registration, regardless of whether Pearson correlation or normalized mutual information was used.

The motion-correction scheme presented here can also be used in conjunction with parallel imaging for scan acceleration. There is essentially no difference between undersampled k-space resulting from conventional parallel imaging and that resulting from local sampling density variations caused by reorganizing k-space for motion correction. The degree to which aliasing can be avoided is determined by the amount of net regional undersampling in k-space. In other words, with increasing data reduction as a result of parallel imaging, the method’s ability to correct for motion diminishes. If more emphasis is to be placed on reducing severe motion, the parallel imaging reduction factor has to be derated accordingly. The ability of the reconstruction ultimately depends on the number of coils used, the geometry factor, and the general conditioning of the design matrix.

One of the shortcomings of the proposed method is its inability to correct for through-plane motion. We found that spin history effects can reduce the quality of the reconstruction, although the final outcome is still significantly better than that achieved without any correction. Ultimately, spin history effects will always challenge 2D acquisition techniques that are exposed to potential movement with higher degrees of motional freedom. Nevertheless, in this work our goal was to investigate the possibility of correcting motion and spatiotemporal coil sensitivity variations using an augmented generalized parallel imaging reconstruction method. This method, in its general form, is expressed as a 3D problem with six degrees of freedom. The application of this technique to 3D acquisitions should also improve the quality of 3D scans, and further work on retrieving 3D motion information in a rapid fashion is warranted. It is important to note that even with prospective motion correction, and regardless of whether Cartesian or arbitrary k-space acquisitions are used, modifications to the original sensitivity information (i.e., translation and rotation) might be required because the patient will be exposed to different coil sensitivities during the course of motion. Therefore, there is a chance that conventional Cartesian parallel imaging algorithms or standard iterative SENSE algorithms do not work properly even when the scan plane can follow the movement of the object.

Although we initially demonstrated the navigator image technique in this study by using a SE approach, we also applied the navigator scheme to regular GRE scans. The limitation is that the minimum TE for the pulse sequence is delayed by a few milliseconds to accommodate the short spiral-IN navigator image. If necessary, and if T₂* decay permits, the navigator signal can also be acquired at the end of the spiral-OUT train (i.e., after the interleaf is acquired that contributes to the high-resolution image). For SEs the spiral-IN navigator is particularly useful because it takes advantage of the high SNR available as the signal increases prior to the SE formation. As shown in Fig. 3a and b, a variable-density spiral-out scan can be used as well. The image quality of variable-density spiral scans is affected somewhat more by off-resonance effects and signal decay compared to conventional, constant-pitch spirals, and it depends on the pitch factor and overall readout length. Nevertheless, this trajectory has great potential because of its inherent navigation capabilities, and is a viable alternative that has been used successfully for self-navigated DTI (23,29). Because a spiral readout is used, the overall image quality of our sequence is slightly impaired by the T₂* decay during the readout, but one could easily improve this by increasing the number of interleaves. The spiral-IN navigator and the parallel imaging-based correction scheme can be also combined with conventional FSE readouts.

CONCLUSIONS

We have proposed a new motion-correction scheme and pulse sequence that takes advantage of the new parallel imaging techniques. Our motion-correction scheme capitalizes on the ability of parallel imaging methods to correct for missing data in k-space. The amount of motion-induced image displacement and rotation is determined by means of registration algorithms applied to navigator images, and can be included in the iterative parallel imaging reconstruction. An appropriate navigator image can be produced with either a variable-density spiral waveform or an additional low-resolution, single-shot spiral scan prior to regular image formation. This new correction approach may avoid the well-known shortcomings of regular Cartesian encoding in the presence of patient motion, and may become a helpful asset for correcting such artifacts, especially in studies of vulnerable subjects. A variant of this method could also be used to allow parallel imaging in prospectively motion-corrected studies.