Magnetic Resonance in Medicine 59:1401–1411 (2008)
Generalized MRI Reconstruction Including Elastic Physiological Motion and Coil Sensitivity Encoding Freddy Odille,1,2 Nicolae Cîndea,1,3 Damien Mandry,1,2 Cédric Pasquier,1,2 Pierre-André Vuissoz,1,2 and Jacques Felblinger1,2
This article describes a general framework for multiple coil MRI reconstruction in the presence of elastic physiological motion. On the assumption that motion is known or can be predicted, it is shown that the reconstruction problem is equivalent to solving an integral equation—known in the literature as a Fredholm equation of the first kind—with a generalized kernel comprising Fourier and coil sensitivity encoding, modified by physiological motion information. Numerical solutions are found using an iterative linear system solver. The different steps in the numerical resolution are discussed, in particular it is shown how over-determination can be used to improve the conditioning of the generalized encoding operator. Practical implementation requires prior knowledge of displacement fields, so a model of patient motion is described which allows elastic displacements to be predicted from various input signals (e.g., respiratory belts, ECG, navigator echoes), after a free-breathing calibration scan. Practical implementation was demonstrated with a moving phantom setup and in two free-breathing healthy subjects, with images from the thoracic-abdominal region. Results show that the method effectively suppresses the motion blurring/ghosting artifacts, and that scan repetitions can be used as a source of over-determination to improve the reconstruction. Magn Reson Med 59:1401–1411, 2008. © 2008 Wiley-Liss, Inc. Key words: Fredholm equation; magnetic resonance; motion correction; parallel imaging; reconstruction
INTRODUCTION The duration of MR acquisition is constrained, mainly by the physics underlying spin relaxation, and cannot be shortened beyond certain limits. The choice of MR sequence parameters is the result of a compromise between the desired spatial and/or temporal resolution, the signalto-noise ratio (SNR), the type of contrast, and the acquisition duration. In the context of cardiac or abdominal imaging, the latter is generally the limiting parameter, as patients cannot hold their breath much longer than 20 sec. Physiological motion is therefore a major impediment to the acquisition of high quality images. Imperfect breath 1 Imagerie Adaptative Diagnostique et Interventionnelle, Nancy University, Nancy, France 2 Inserm, ERI13, Nancy, France 3 Institut Elie Cartan de Nancy Ministère Nancy University, Nancy, France Grant sponsor: Ministére de l’Industrie et de Recherche, Region of Lorraine, France; Grant number: RNTS2003; Grant sponsors: Insern, Region of Zorraine *Correspondence to: Jacques Felblinger, IADI, Inserm ERI13, CHU de Nancy Brabois, 54511 Vandoeuvre-lès-Nancy, France. E-mail:
[email protected] Received 15 June 2007; revised 2 October 2007; accepted 20 November 2007. DOI 10.1002/mrm.21520 Published online 17 April 2008 in Wiley InterScience (www.interscience.wiley. com).
© 2008 Wiley-Liss, Inc.
holds result in motion artifacts that, depending on the kspace acquisition trajectory, can take the form of blurring or aliasing of the moving structures onto different parts of the image—also called “ghost” artifacts (1). Even with cooperative patients, the position of organs is likely to drift before the scan is complete (2). In addition to respiratory motion, cardiac contraction is another problem that requires using specific techniques, such as retrospective or prospective gating. Although these techniques are well known and routinely used in clinical practice, they have the drawback of increasing acquisition duration. As a consequence, combined cardiac/respiratory gating results in longer acquisitions, which are often not compatible with clinical protocols. Numerous methods have been proposed for handling physiological motion better, with several means of correcting motion-induced artifacts being put forward. In (3), the authors proposed a fully automatic method, based on minimization of a cost function quantifying the reconstructed image quality, with respect to the motion parameters describing motion. The only information needed by the method is a description of motion using a model with few parameters (e.g., an affine model), so that the optimization process is feasible in practice. The choice of the cost function remains a critical point, because classical image quality criteria like those based on Shannon entropy result in nonlinear optimization problems, which are difficult to solve, especially in high dimensional parameter search spaces. Other methods have been proposed that capitalize on prior knowledge acquired from a model of patient motion. These methods involve either adapting the reconstruction algorithm (4–6) or adapting the pulse sequence in real-time (7–9), known as prospective motion correction. Indeed the authors show that, assuming that motion is predictable or can be described by a model (a parametric model or information concerning complete displacement fields), it is possible to invert the process of spatial encoding of moving structures that underlies the production of artifacts. In all these methods, implementation was restricted to rigid or affine motion. Yet, in theory, the reconstruction proposed in (5) is not limited to such a model. While simple translations and rotations, or affine models, are generally sufficient for describing head motion, they are less accurate for assessing complex elastic deformation such as can occur in thoracic or abdominal organs. It would be an advantage to have a general framework for handling arbitrary physiological motion for these kinds of imaging protocols. The difficulty here lies in determining an appropriate model
1401
1402
Odille et al.
that could take into account such elastic motion. Moreover, for the framework to be complete, it would be desirable to include the parallel imaging theory, as in (6). This article proposes a general formulation of the MR image reconstruction process in the presence of elastic physiological motion and in the context of multiple coil acquisition. Although our formulation is slightly different, the method can be seen as a combination of the generalized SENSE parallel imaging algorithm (10), and the general motion compensation technique presented in (5). To be applicable in practice, the method requires prior knowledge of patient motion. A predictive model is therefore proposed to estimate the displacement fields needed by the reconstruction algorithm. The motion model involves predicting the displacement fields at certain control points defined by a Cartesian grid, using linear combinations of various input signals derived from respiratory belts and ECG data. The whole framework was validated with real data from a moving phantom, achieving periodic translational motion, as well as in healthy volunteers who were subjected to free-breathing ECG-triggered pulse sequences. The solution of the generalized reconstruction process was compared to the standard Fourier reconstruction, and to a reference breath hold (i.e. artifact-free) acquisition. THEORY Generalized Reconstruction Framework The proposed formulation relates the k-space data s(k), with k = [kx , ky , kz , γ ]T ∈ , collected by the different coils (indexed by γ = 1 · · · Nγ ), to the proton density ρ. In a classic reconstruction approach (Fourier transform or parallel imaging), the imaged subject is assumed to be static during the whole acquisition process. To account for motion, the position r = [x, y , z]T ∈ should be considered as a variable of the acquisition time, i.e., r = f (r0 , t), and thus as a variable of the k-space position, i.e., t = t(k)—since the MR acquisition process, for each individual coil, is sequential. We also introduce the downsampling operator d (d(k) = 0 or 1 for a Cartesian sampling) that can be used to accelerate k-space acquisition, as well as coil sensitivity σ (γ , r) – or equally σ (k, r). Neglecting the relaxation terms, the MR signal in a given point of the coil k-space array is given by: s(k) = d(k)ρ(r, t(k))σ (k, r)e−ik.r dr [1]
with the notation k · r = kx x + ky y + kz z. It should be noted that, for simplicity, we do not introduce explicit dependency of coil sensitivity over time [e.g. we do not write σ (k, r, t)]. This amounts to considering the coils as static. In addition, when the patient breathes, the coil loading changes slightly, so we are also disregarding the variation in sensitivity induced by these changes in loading. However, the sensitivity changes seen by the voxels, due to motion, are taken into account because of the dependency of σ on r. We assume that the motion function f is known or can be estimated by a predictive model. Furthermore, we assume that, for all acquisition times t(k), φk : r0 → r = f (r0 , t(k)) is “nearly bijective”, i.e. that it is injective, differentiable, has continuous partial derivatives, and has nonzero Jacobian
Jφk for all r0 ∈ . Under the latter hypotheses, a substitution rule can be applied in the integral. For simplicity, we also assume that φk () = , which means that the border ∂ of the field of view (FOV) is not deformed during the acquisition. The substitution between r and r0 yields the following equation: K (k, r0 )ρ0 (r0 )dr0 [2] s(k) =
with K (k, r0 ) = d(k)σ (k, φk (r0 ))e−ik·φk (r0 ) |Jφk (r0 )|, and ρ0 denoting the proton density at the reference time t0 from which the displacement fields are computed. Using operator notations, Eq. [2] reads: s = Eρ0
[3]
where E is a generalized encoding linear operator. Eqs. [2] or [3] relates the MR signal s to the proton density at the reference position ρ0 , via an integral equation of kernel K, comprising the downsampling and the coil sensitivity encoding operators used in parallel imaging, and a Fourier encoding term, modified with physiological motion information. The generalized reconstruction problem amounts to inverting Eq. [3]. Equation [2] is known in the literature as a Fredholm equation of the first kind. Generally the assumption is made that K is compact. This assumption can be made here because, in practice, MR experiments consist in finite discrete sampling of the coil k-spaces . Such integral equations are known to be ill-posed in Hadamard’s sense, that is, existence, uniqueness, and stability of the solution cannot be proven in the general case. However, the problem can be discretized to obtain numerical solutions, and Tikhonov regularization (11) may be used to impose additional smoothness constraints on the solution. A practical formulation of Eq. [3] can be given, using the adjoint operator E such that: K (r, k)u(k)dk = K (k, r)u(k)dk [4] (E u)(r) =
Nc kxmax kymax kzmax
where we noted · = [3] is then equivalent to:
γ =1 0
b = Aρ0
0
0
·. Equation [5]
with b = E s and A = E E. Numerical Resolution On the basis of formulation in [5] the discretization leads to a large linear algebraic system. Without loss of generality, in the following we will consider the problem of 2D spatial encoding, 3D would be treated similarly; it should be noted that in the case of 2D, motion is assumed to be in-plane. If the image to reconstruct is of size N ×N , ρ0 and b are N 2 ×1 vectors, and A is an N 2 × N 2 matrix. A crucial point in understanding generalized reconstruction algorithms of this kind—i.e. involving an N 2 × N 2 linear system—is that although the operator is large and dense, explicit knowledge of its N 4 coefficients is not required to solve the inverse problem. Instead, an iterative strategy can be adopted that only requires explicit
Generalized MRI Reconstruction
knowledge of the function x → Ax. Examples of such iterative linear system solvers include the conjugate gradient (CG), the generalized minimum residual (GMRES), and the biconjugate gradient (BiCG). In the present problem, Ax can be rewritten in matrix formalism, by separating the different encoding steps, as is done in (5). With Dt representing the downsampling matrices corresponding to the tth acquisition (i.e., the operator corresponding to the tth line(s) selection), σγ the discrete coil sensitivity operators, F the 2D Fourier operator, Tt the spatial transformation matrices, in discrete version, A reads : Nγ T TtH [6] A= σγH F H DtH Dt F σγ Tt t=1
γ =1
with the subscript X H denoting the Hermitian transpose. All operators involved in this sum are N 2 ×N 2 matrices. A is advantageously described in this form because in Eq. [6], F is the only operator that is dense, but can be computed efficiently via the fast Fourier transform (FFT) algorithm rather than by matrix multiplication. Indeed the coil weightings σ γ are diagonal, and the Tt are sparse. x → Tt x can be computed by means of an interpolation algorithm, a point that will be discussed later (the number of nonzero elements in Tt depends on the interpolation basis that is used). Conditioning of the Generalized Encoding Operator Here A has the advantageous property of being square and symmetric by construction. Moreover, if A is fullrank, it is positive definite, and thereby invertible. The above-mentioned iterative algorithms have good convergence properties when the latter conditions are fulfilled. However, in practice, A comprises coil sensitivity maps and displacement field predictions, which are both determined on the basis of experimental data. With regard to parallel imaging, over-determination is introduced to address this problem, by ensuring that the ratio of the number of coils over the acceleration factor is greater than 1. This stratagem amounts to forcing E to have more rows than columns, and hence increases the rank of A, or, equally, the number of significant eigenvalues of A. In addition to this, another source of over-determination should be used in the present problem to take into account bad conditioning induced by physiological motion encoding. The reconstruction will indeed be sensitive to motion prediction errors. It should also be noted that, even if no error were introduced by the prediction itself, the k-space trajectory deformation induced by motion is expected to yield certain k-space regions that are oversampled and others that are undersampled. Another source of over-determination directly linked to the motion encoding part of A should be introduced in order to overcome this problem. This can be easily done by repeating the k-space acquisition, i.e., by setting up more than one number of excitation (NEX). Unlike the conventional approach in which data from different excitations are averaged prior to reconstructing, here data from each excitation are handled as independent shots, because they are not acquired in the same cardiac/respiratory configuration. If A is still not well conditioned for numerical resolution, due to motion modeling errors for example, regularization can be employed. Tikhonov regularization (11) is a
1403
well-known technique which, in practice, is equivalent to replacing A by A + λId in [5], with Id being the identity operator, and λ a (small) tuning parameter. Interpolation Issue In Eq. [6] the downsampling, Fourier, and sensitivity weighting operators have images that remain on a Cartesian grid, which is not the case for the spatial transformation applied at each acquisition time. Similarly to image registration problems, we avoid creating holes in the transformed image by using the inverse transformation. Thus, data from the transformed image are collected by interpolating intensity values in the original image. Bilinear interpolation is a common choice in image processing problems, which offers a good compromise between regularity of the interpolated image and time needed for the computation. Other bases can be used, such as cubic splines or sinc based kernels, depending on the kind of regularity we want to impose on the reconstructed image. Computational Cost In certain cases the algorithm defined by Eq. [6] can be optimized. In particular, if the acquisition scheme is a line by line (or set of lines by set of lines) covering of k-space, the raw data experience, at each shot, a downsampling in the phase-encoding direction(s) only. As a result, aliasing at the origin of ghost artifacts also appears in the phase-encoding direction(s) only. Then, in a similar way to that used in the SPACE-RIP parallel imaging algorithm (12), the computation x → F H DtH Dt Fx can be reduced by one dimension, by applying all operations in the phase-encoding dimension(s), and keeping the frequency dimension unchanged. In addition to this optimization, the Fourier operator can be computed very efficiently by the FFT or FFTW algorithm. Hence the resolution will be constrained mainly by the interpolation, especially if a more sophisticated basis is desired. The Motion Model As mentioned previously, generalized reconstruction assumes that displacement fields can be predicted by a model. Several models which are driven by various input signals have already been proposed. Examples of signals that can serve as inputs include navigator echoes (6,13–15), respiratory belts (16), or combinations of various external sensor signals (belts, R-wave amplitude of the ECG) (17). Figure 1 illustrates how motion is predicted. The model consists of mapping these input signals to true motion. To do this, a calibration scan is performed covering several cardio-respiratory cycles. Then a regression method is employed to determine linear combinations of inputs that best estimate true displacement fields, or a few parameters describing motion (e.g., rigid or affine transformation). Thus the linear model coefficients can be derived from auto-regressive or adaptive models (16), principal component analysis (PCA) (13,14,18) or partial least squares regression (PLSR) (15). The latter regression methods have been proposed in order to handle the problem of possible collinear signals within the model inputs. We proposed an
1404
Odille et al.
represents K maps of N × N coefficients for motion prediction in the x direction, and K maps of N × N coefficients for the y direction. Once the model coefficients αk (r) are determined, displacement fields can be predicted in any subsequent MR scan, using the formula: ˆ t) = α(r)S(t) = u(r,
K
αk (r)Sk (t)
[8]
k=1
FIG. 1. Schematic diagram illustrating the predictive motion modeling method. From a fast MR imaging series, simultaneous acquisition of image and sensor data (pneumatic belts, ECG data, ...) allows to construct, in a calibration step based on a regression method, a linear model described by the coefficients αk (r ). Using newly acquired sensor data, the model is able to estimate the displacement fields, in the plane/volume of interest, needed by the generalized reconstruction method.
alternative regression method in (17), based on a variational method, that helps handle arbitrary elastic motion, and has the advantage over the other methods mentioned of using an additional local constraint, based on geometry, to solve the regression problem. Elastic motion is represented by displacement values at certain control points, defined by a Cartesian grid of size N × N . The basic assumption here is that, if the displacement of a particular point in the field of view (FOV) is followed, its position can be approximated using a linear combination of model inputs (predictability assumption). The second assumption is that, if two particular points within the FOV are close to each other, their displacements are highly correlated, and hence their linear relationship to the model inputs should be similar (local consistency assumption). We use a variational approach to find linear combinations that satisfy these two hypotheses. We call S(t) = [S1 (t) . . . SK (t)]T the model input signals and u(r, t) the displacement field at time t ∈ [0, T ], where T is the duration of the dynamic calibration series acquisition. The objective is to find a map of linear combinations of the Sk (t) that estimate the true displacement fields with minimal error. To do this we define αk (r) (k = 1 · · · K ) as being the projections of the r point displacement on the (Sk (t))k=1···K basis. Hence we are looking for the vector α = [α1 , . . . , αK ]T , minimizing the following functional: F=
0
T
||S(t)α(r) − u(r, t)||2 dt dr + µ R(∇α)dr
[7]
with the local consistency term (or regularization term), being R = ||∇α||2 (Tikhonov regularization), or R = ||∇α|| (TV regularization). The above functional is composed of two terms, corresponding to the predictability and the local consistency assumption, respectively (the parameter µ controls the trade-off between the two terms). The TV constraint is also proposed because the Tikhonov method is known, from image processing applications of variational methods, for possibly producing excessive smoothing. Resolution of the variational problem is achieved by writing the Euler–Lagrange equations. Their discretization leads to a sparse algebraic system. Its solution vector α
It should be noted that, in the limit case where µ = 0, the solution decomposition is simply the projection of each individual voxel displacement onto span({Sk }k=1···K ), which is equivalent to the PCA method in (13) or (14) when no eigenvalue is truncated. It should be possible to use either model. Specifically, if µ is small or if almost all eigenvalues are kept, the variational or PCA models will yield similar results. The rational underlying the variational method is that when sources of indetermination exist in the model (collinear inputs or image areas, which are poorly correlated with inputs), the model coefficients are chosen based on the coefficients of the neighboring areas, according to the local geometric constraint. METHODS Imaging and Instrumentation MR experiments were performed using a 1.5 T GE Signa Excite HDX MR system (General Electric, Milwaukee, WI). Physiological signals were collected using a modified version of the Maglife (Schiller Médical, Wissembourg, France) patient monitoring system. These signals were converted into optical signals, transmitted outside the MR bore through optical fibers, and then recorded on the signal analyzer and event controller (SAEC) computer and electronics system, presented in (19). Two pneumatic belts (one each on the thorax and abdomen) were used to monitor respiration. One ECG sensor was used for cardiac triggering and information about the cardiac phase. The SAEC system allows synchronous acquisition and recording of physiological signals (respiration and ECG) and MR signals (image acquisition windows and MR gradients), real-time processing of ECG signals (correction of MR gradient switching interference (19,20) and R-wave detection). MRI data were acquired with an ECG-triggered black blood RARE pulse sequence (black blood FSE, TE = 35.9 msec, ETL = 16, TI = 650 msec, 256 × 256, NEX = 1, 1.4 mm square pixels, 10-mm slice thickness), using an eight element cardiac coil. This pulse sequence was chosen because it is used routinely in clinical practice for high resolution cardiac imaging, and requires, with the conventional method, a breath hold of approximately 20 sec. Full sampling of the k-space was used in all experiments, to assess the amount of artifact reduction due to motion only (as opposed to parallel imaging artifacts due to imperfect unfolding). Prior to the pulse sequence of interest, a free-breathing calibration scan was performed to determine the linear motion model coefficients (2D FIESTA, TE = 0.9 msec, TR = 2.6 msec, 128 × 128 matrix, 2.8 mm square pixels, 10-mm slice thickness), resulting in a temporal resolution of 3.6 frames per second. The displacement field
Generalized MRI Reconstruction
estimates needed for the calibration were determined by pyramidal implementation of the Lucas and Kanade optical flow method (21) from the Intel open source computer vision library (22). Other methods could have been used to find these displacement data, such as those presented in (23,24) for example. The motion model coefficients were found using the proposed variational method (grid size of 128 × 128, TV regularization, µ = 0.01). Once determined, the linear model coefficients were used in combination with sensor measurements from newly acquired data to predict elastic displacement fields in the entire FOV.
Moving Phantom Setup A moving platform setup was developed to validate the present method in phantom experiments. The setup comprised both static and moving phantoms, resulting in piecewise translations in the FOV. The moving part achieved 1D periodic translation, with a pneumatic command generated by a ventilator simulating respiratory motion. The cardiac synchronization signal was generated by an ECG simulator device, in order that the configuration was the same as in our experiments with volunteer subjects. The imaging plane was chosen so that the moving phantom remains in plane. The setup produced displacements of 30 mm in amplitude in the up/down direction, and 5 mm in the right/left direction (see images in Fig. 2), and the period for translational motion was set to 5 sec. These settings are typical of adult respiratory motion. The high resolution scan was performed four times (NEX = 4) while the phantom was moving, and a reference scan was also acquired in static mode for purposes of comparison. The static scan was also used to compute the coil sensitivity maps, by extracting the 32 central k-space lines. In the motion model, because of its simplicity and reproducibility, only two input signals were considered: the first one was provided by a standard pneumatic belt placed on the moving phantom setup, and the second one was the derivative of the pneumatic belt signal. Using the derivative signal allows to take into account possible delays or nonlinearities between the pneumatic signal and true motion, and provides information similar to the precursory navigators proposed in (14).
Free-Breathing Subject Data Two healthy volunteers were subjected to our protocol for experiments with real elastic displacements from the thoraco-abdominal region. Two particular slice locations were prescribed: one in the right liver and lung (first subject), the other in the heart and left liver (second subject). In both subjects, the imaging plane was chosen sagittal to minimize through-plane motion effects. These imaging planes were not targeted for a specific clinical application, our primary concern was rather to put forward examples of elastic (nonrigid or affine) in-plane motion. For each experiment, the high resolution free-breathing scan was performed three times, leading to an over-determination factor of 3 in the motion encoding operator. The high resolution scan was also acquired once in breath hold for comparison purposes. The breath hold scan was also used to compute the
1405
coil sensitivity maps, by extracting the 32 central k-space lines. Motion model inputs included four input signals intended for monitoring respiratory motion—two pneumatic belts (one each on the thorax and abdomen), and their derivatives (as in the phantom experiment). Two additional inputs were used to monitor cardiac contraction, based on ECG data. With ti being the R-wave peak time of the ith ECG cycle, and RRi = ti+1 − ti , the two 2π additional cardiac phase signals were defined by cos( RR t) i 2π t), ∀t ∈ [t , t ]. This is equivalent to approxiand sin( RR i i+1 i mating the displacements due to cardiac contraction by a pseudo sinusoid, with frequency adjustable to heart rate, and with tunable amplitude and phase. The latter two signals were of use essentially for the model calibration scan, because it was not ECG-triggered, and as a result displacement fields are dependent on both the cardiac and respiratory phase. However, during the high resolution scans, these two signals were approximately constant due to the cardiac synchronization, so the predicted displacement fields have been computed for the cardiac phase of interest (in diastole here, at t = 650 msec in the cardiac cycle). Validation During the calibration phase, all displacement fields were computed relative to a reference image. In order to assess the reconstruction quantitatively, we chose as reference the image in the calibration series which was closest to the image acquired in breath hold, as assessed by the normalized mutual information (NMI). Hence the reconstructed image was approximately in the same respiratory phase as the breath hold image, because data from each shot were “registered” relative to that reference. However, since the calibration series has lower spatial resolution, a slight misalignment can be expected. Before comparing the proposed reconstruction with the breath hold scan, the slight misalignment was compensated by an elastic registration, using the same registration scheme as for the motion model calibration (optical-flow-based method). The interpolation involved in this registration process was performed with a windowed sinc kernel. The quality of the reconstructed image was assessed in terms of various criteria for similarity with the reference breath hold image (after the minor realignment). We chose the following criteria, which are commonly used in image processing: mean absolute error (MAE), correlation coefficient (CC), joint entropy (JE), and normalized mutual information (NMI). Moreover, we assessed image quality of the standard Fourier reconstruction and of the proposed reconstruction using the two following methods: image entropy was computed to quantify signal dispersion (entropy increases with noise or blurring) and gradient entropy was computed to quantify singularity (i.e., the amount of fine details) of the image. Implementation Details The reconstruction algorithm was implemented in the Matlab (MathWorks, Natick, MA) development environment.
1406
Odille et al.
FIG. 2. Reconstruction from the moving phantom, achieving in-plane piecewise translation: (a) standard Fourier reconstructions from NEX = 1 to NEX = 4, (b) generalized reconstructions with the proposed method from NEX = 1 to NEX = 4, (c) the static scan image for reference, and (d) the difference images between the generalized reconstructions and the reference from NEX = 1 to NEX = 4. The moving part of the setup (bottle, balls, and resolution phantom on the right of the image) achieved periodic translation with a maximum amplitude of 29.5 mm in the SI direction (up/down) and 5.1 mm in the AP direction (right/left).
Since the pulse sequence used in our experiment was a multishot scan, each term of the sum in [6] corresponds to one single echo train. This amounts to considering that the period needed to acquire this echo train—or equally, this set of k-space lines—was negligible with respect to physiological motion. Moreover, in [6], each term in the sum (a sum over echo trains here) is independent, making the algorithm highly parallelizable. The results presented in the following section were computed with custom parallelization software, based on Matlab, on four processors for a relatively fast reconstruction (less than 5 min). Bilinear interpolation was chosen for implementation of the spatial transformation operator involved in [6], and the Fourier operators employed the FFTW algorithm. The Tikhonov regularization parameter was empirically set to λ = 0.1 (a more systematic method could be employed, e.g., based on the study of L-curves (25)). The algorithm performed iterations until convergence was reached, as assessed by a tolerance that was fixed at = 10−3 . The GMRES was used as the iterative method to solve [5], because in our experiments, it was found to converge slightly faster to the desired tolerance than the conjugate gradient. With these settings, the time needed for the reconstruction ranged from 2 to 5 min according to the NEX and conditioning of the system.
RESULTS Moving Phantom Results Results from the moving phantom reconstructions are shown in Fig. 2. Images reconstructed by the standard Fourier transform show typical ghost artifacts produced by the moving part of the setup in the phase encoding direction. When employing averaging, the Fourier reconstruction converges to a blurred solution. Indeed it was shown in (26) that the effect of averaging is a convolution with a kernel whose width is approximately the range of motion. The generalized reconstruction, using displacement fields predicted by the motion model proposed in the second section, was able to remove these artifacts very efficiently. The static phantom reconstruction is also shown for reference, as well as the difference images between the generalized reconstruction and the static reference. Quantitative assessment is given in Table 1. Results clearly show that, for the various similarity criteria used for comparison with the static reference image and for image quality assessment, the generalized reconstruction was better than the standard Fourier reconstruction, and improved when the NEX increased. In this example, increasing the NEX up to 2 or 3 produced the most significant improvement.
Generalized MRI Reconstruction
1407
Table 1 Comparison of the Generalized Reconstruction with the Standard Fourier Reconstruction in the Moving Phantom Experiment NEX Criteria Similarity with static reference
MAE (less is better) CC (more is better) JE (less is better) NMI (more is better)
Image quality
Entropy (less is better)
Gradient entropy (more is better)
Reconstruction
1
2
3
4
Fourier Generalized Fourier Generalized Fourier Generalized Fourier Generalized Reference Fourier Generalized Reference Fourier Generalized
0.078 0.057 0.789 0.924 8.536 7.988 1.214 1.221
0.073 0.031 0.815 0.973 8.436 8.033 1.228 1.269
0.068 0.024 0.831 0.983 8.398 8.011 1.234 1.285
0.067 0.024 0.843 0.982 8.391 7.998 1.236 1.285
5.249 4.913
5.247 5.138
5.247 5.192
5.252 5.187
2.142 2.529
2.061 2.373
1.989 2.348
2.009 2.369
Free-Breathing Subject Data Displacement fields predicted by the motion model are represented in Figs. 3 and 4 for two example subjects, with imaging planes located in the liver and in the heart respectively. In these figures, the fields are represented in one typical respiratory cycle taken from the free-breathing scan. The maximum displacement (i.e., the amplitude between inspiration and expiration), as predicted by the model, was found to be 21.3 mm in the superior–inferior direction (SI) and 6.2 mm in the anterior-posterior direction (AP) for the first subject (slice in the liver), and 12.6 mm in the SI and 6.9 mm in the AP direction for the second subject (slice in the heart). The color maps in these examples show how the motion model was used to approximate elastic displacements from the thoraco-abdominal region, whereas rigid or affine models generally fail in these situations, due to the limited number of degrees of freedom. Reconstructed images from the two subjects are presented in Figs. 5 and 6. Visual comparison of these images with the reference breath hold scan clearly shows artifact reduction in the generalized reconstruction, compared with the Fourier reconstruction, in terms of blurring/ghosting artifacts. Anatomical details have been annotated. In the liver image, the vessels in the liver and in the lung are visible in the proposed reconstruction, whereas they are blurred in the standard Fourier reconstruction. The diaphragm also appears better. In the cardiac image, various anatomical details are improved significantly, such as the pericardium, the pulmonary veins, and the left anterior descending (LAD) coronary artery. Blurring in the liver is also reduced. Quantitative results for these two subjects are grouped in Table 2. Improvement from NEX = 1 to NEX = 3 is put in perspective by the different similarity criteria (MAE, CC, NMI, JE). Comparison of the reference with the reconstruction with (generalized) and without (Fourier) motion correction is more difficult to assess, because not all images have the same characteristics in terms of signal-to-noise ratio (SNR). Moreover, since the breath hold can hardly be perfect, residual motion always exists, and the reference image cannot be totally free of artifacts. However, the CC
5.116
2.310
and entropy based criteria (NMI and JE) show that the generalized reconstruction was always better than the Fourier reconstruction for NEX = 2 or NEX = 3. Image quality, as assessed by signal dispersion (entropy), also shows improvement with the proposed method.
DISCUSSIONS AND CONCLUSIONS We have presented a general formulation of the MR reconstruction problem in the presence of elastic physiological motion. On the basis of the work in (5,6,10), we extended the method for general motion compensation to the use of multiple coils. We wrote the reconstruction problem in continuous formulation, which explains the requirements of generalized reconstruction more clearly, in particular the hypotheses that are made concerning displacement fields. Theoretically they should be “nearly bijective” (i.e., they should have continuous partial derivatives and nonzero Jacobian) and be null at the borders of the FOV. We have discussed the numerical implementation and, importantly, we have proposed a simple means of adding over-determination to the problem by repeating the high resolution free-breathing scan several times. Experimental results, obtained with a specially designed moving phantom setup as well as from healthy subject data, confirmed that the reconstruction clearly improves when the NEX increases, and that artifacts from these elastic displacements can be removed efficiently. Unlike previous methods, the present reconstruction was not restricted, in either theory or practice, to rigid or affine models, so that it can be applied to body regions that experience complex elastic displacements such as are observed in cardiac or abdominal imaging. This was possible thanks to a motion model allowing such displacements. The model needed a free-breathing calibration scan, prior to the scan of interest, to identify the model coefficients. The model was able to predict elastic displacement fields in the FOV by linear combinations of newly acquired sensor data (derived from respiratory belts and ECG data). The accuracy of the predicted fields, in the results presented here, was
1408
Odille et al.
FIG. 3. Displacement fields from a sagittal slice located in the liver, predicted by the motion model, using 2 input signals (1 pneumatic belt signal, located on the abdominal wall, and its derivative). Fields in the SI and AP directions are shown in 5 positions from the free-breathing scan, taken in one typical respiratory cycle. The maximum amplitude of the predicted displacements (difference between expiration and inspiration) was 21.3 mm in SI and 6.2 mm in AP.
FIG. 4. Displacement fields from a sagittal slice located in the heart, predicted by the motion model, using 4 input signals (2 pneumatic belt signals, one each on the thorax and abdomen, and their first derivatives). Fields in the SI and AP directions are shown in 5 positions from the free-breathing scan, taken in one typical respiratory cycle. The maximum amplitude of the predicted displacements (difference between expiration and inspiration) was 12.6 mm in SI and 6.9 mm in AP.
Generalized MRI Reconstruction
1409
FIG. 5. Reconstruction from a free-breathing healthy subject (sagittal slice located in the right liver, NEX = 3): (a) standard Fourier reconstruction, (b) generalized reconstruction with the proposed method, (c) the breath hold image for reference. The predicted displacement fields used by the generalized reconstruction are shown in Fig. 3.
FIG. 6. Reconstruction from a free-breathing healthy subject (sagittal slice located in the heart, NEX = 3): (a) standard Fourier reconstruction, (b) generalized reconstruction with the proposed method, (c) the breath hold image for reference. Anatomical structures, such as LAD, myocardium, RV trabeculae, are finely depicted on the generalized reconstruction, at least as well as on the reference image. The predicted displacement fields used by the generalized reconstruction are shown in Fig. 4. Table 2 Comparison of the Generalized Reconstruction with the Standard Fourier Reconstruction in Experiments with Two Healthy Subjects Subject 1 (liver) NEX Criteria Similarity with static reference
MAE (less is better) CC (more is better) JE (less is better) NMI (more is better)
Image quality
Entropy (less is better)
Gradient entropy (more is better)
Subject 2 (heart) NEX
Reconstruction
1
2
3
1
Fourier Generalized Fourier Generalized Fourier Generalized Fourier Generalized Reference Fourier Generalized Reference Fourier Generalized
0.027 0.043 0.986 0.980 7.482 7.461 1.248 1.241
0.028 0.034 0.982 0.984 7.506 7.406 1.244 1.257 4.668 4.672 4.638 2.117 1.953 2.007
0.029 0.031 0.979 0.985 7.546 7.402 1.239 1.260
0.017 0.020 0.983 0.981 7.188 7.139 1.217 1.220
4.678 4.656
4.427 4.372
1.951 2.026
1.880 1.877
4.670 4.591 1.995 1.968
2 0.015 0.015 0.987 0.989 7.112 6.994 1.231 1.247 4.337 4.418 4.385 1.969 1.861 1.859
3 0.014 0.013 0.987 0.990 7.080 6.982 1.237 1.251 4.419 4.398 1.856 1.828
1410
sufficient to put forward significant improvements of the reconstructed image quality. The work presented here was in 2D, with (approximately) in-plane motion, in order to limit the reconstruction time (less than 5 min for a 256 × 256 slice, with an 8 element coil array), but the theory is applicable in 3D too. It should be noted that the computation time depends on how the sum in [6] is computed. Lines whose acquisition experiences similar displacements should be grouped as much as possible, to reduce the number of spatial transform operations. In this work, all lines in one echo train were grouped. The theory includes parallel imaging: when the displacement fields are null the algorithm is equivalent to the generalized SENSE algorithm (10) (this work was restricted to the Cartesian case, but noncartesian sampling would be treated similarly). Hence, it is possible to accelerate the scan by undersampling the k-space. Keeping in mind that parallel imaging reconstructions are penalized by an SNR loss, this property could be used to compensate, in part, the time spent on scan repetitions. The method has several limitations that should be mentioned. The implicit assumption was made that motion artifacts result from spatial encoding errors, i.e., that motion occurs between k-space line acquisitions (interview motion) and not during k-space line acquisition (intraview motion). Pulse sequences with long TR should be treated carefully because motion between a preparation pulse and spatial encoding pulses is likely to induce signal modifications (signal loss or incomplete blood or fat suppression for example). In the pulse sequence we used here, the black blood inversion occurred TI = 650 msec before spatial encoding and effective MR acquisition. However no signal modification was observed, because through-plane motion was negligible, due to the slice orientation being sagittal, and to the relatively thick slice (10 mm). More generally, the implicit assumption is that the mass of protons is conserved within the imaged slice/volume. If this hypothesis is violated, such as in the case of through-plane motion or intensity changes due to the arrival of a bolus, the method will not be applicable. Experiments carried out here did not require any pulse sequence modification, as opposed to navigator echo-based methods. This is an important asset because all existing MR pulse sequences could potentially take advantage of the technique. The proposed approach relies rather on a dedicated system to monitor various items of information from external sensors. This information provides indirect representation of motion, whereas navigator pulses provide direct image data but with less temporal accuracy, because navigator pulses are not acquired simultaneously with the image and may require significant acquisition time. When navigator echoes are available though, they can be added to the model inputs in order to improve prediction accuracy. The method also relies on how well the motion model is able to predict displacement fields. The assumption was made that the displacement of any voxel in the FOV can be approximated by linear combinations of the chosen sensor signals. This assumption may be wrong especially when large displacements occur. In this case, a more complex model should be used to account for the non-linearities, by
Odille et al.
adding new linearly independent inputs to the model (e.g. other motion sensors).
ACKNOWLEDGMENTS The authors thank Schiller Médical for providing an important part of the instrumentation required.
REFERENCES 1. Wood ML, Henkelman RM. MR image artifacts from periodic motion. Med Phys 1985;12:143–151. 2. Plathow C, Ley S, Zaporozhan J, Schöbinger M, Gruenig E, Puderbach M, Eichinger M, Meinzer HP, Zuna I, Kauczor HU. Assessment of reproducibility and stability of different breath-hold maneuvres by dynamic MRI: Comparison between healthy adults and patients with pulmonary hypertension. Eur Radiol 2006;16:173–179. 3. Atkinson D, Hill DLG, Stoyle PN, Summers PE, Clare S, Bowtell R, Keevil SF. Automatic compensation of motion artifacts in MRI. Magn Reson Med 1999;41:163–170. 4. Pipe JG. Motion correction with PROPELLER MRI: Application to head motion and free-breathing cardiac imaging. Magn Reson Med 1999;42:963–969. 5. Batchelor PG, Atkinson D, Irarrazaval P, Hill DLG, Hajnal J, Larkman D. Matrix description of general motion correction applied to multishot images. Magn Reson Med 2005;54:1273–1280. 6. Bammer R, Aksoy M, Liu C. Augmented generalized SENSE reconstruction to correct for rigid body motion. Magn Reson Med 2007;57:90–102. 7. Ehman RL, Felmlee JP. Adaptive technique for high-definition MR imaging of moving structures. Radiology 1989;173:255–263. 8. Thesen S, Heid O, Mueller E, Schad LR. Prospective acquisition correction for head motion with image-based tracking for real-time fMRI. Magn Reson Med 2000;44:457–465. 9. Nehrke K, Börnert P. Prospective correction of affine motion for arbitrary MR sequences on a clinical scanner. Magn Reson Med 2005;54:1130– 1138. 10. Pruessmann KP, Weiger M, Börnert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med 2001;46:638–651. 11. Tikhonov AN, Arsenin VY. Solution of Ill-posed problems. Washington: Winston; 1977. 12. Kyriakos WE, Panych LP, Kacher DF, Westin CF, Bao SM, Mulkern RV, Jolesz FA. Sensitivity profiles from an array of coils for encoding and reconstruction in parallel (SPACE RIP). Magn Reson Med 2000;44:301– 308. 13. McLeish K, Hill D, Atkinson D, Blackall J, Razavi R. A study of the motion and deformation of the heart due to respiration. IEEE Trans Med Imag 2002;21:1142–1150. 14. Manke D, Nehrke K, Börnert P. Novel prospective respiratory motion correction approach for free-breathing coronary MR angiography using a patient-adapted affine motion model. Magn Reson Med 2003;50:122– 131. 15. Ablitt NA, Gao J, Keegan J, Stegger L, Firmin DN, Yang GZ. Predictive cardiac motion modeling and correction with partial least squares regression. IEEE Trans Med Imag 2004;23:1315–1324. 16. Priatna A, Paschal CB, Shiavi RG. Evaluation of linear diaphragmchest expansion models for magnetic resonance imaging motion artifact correction. Comput Biol Med 1999;29:111–127. 17. Odille F, Pasquier C, Vuissoz PA, Felblinger J. Linear predictive modeling of patient motion using external sensors. In: Proceedings 15th Scientific Meeting, International Society for Magnetic Resonance in Medicine, Berlin, 2007. 18. Pasquier C, Odille F, Abächerli R, Vuissoz PA, Felblinger J. Modeling of organ’s motion using external sensors. In: Proceedings 14th Scientific Meeting, International Society for Magnetic Resonance in Medicine, Seattle, 2006. 19. Odille F, Pasquier C, Abächerli R, Vuissoz PA, Zientara GP, Felblinger J. Noise cancellation signal processing method and computer system for improved real-time electrocardiogram artifact correction during MRI data acquisition. IEEE Trans Biomed Eng 2007;54:630–640. 20. Abächerli R, Pasquier C, Odille F, Kraemer M, Schmid JJ, Felblinger J. Suppression of MR gradient artefacts on electrophysiological signals
Generalized MRI Reconstruction based on an adaptive real-time filter with LMS coefficient updates. Magn Reson Mater Phys 2005;18:41–50. 21. Lucas BD, Kanade T. An iterative image registration technique with an application to stereo vision. In: Proceedings of imaging understanding workshop, Washington, D.C. 1981. pp 121–130. 22. Intel. OpenCV: Open Source Computer Vision Library. 23. Christensen G, Rabbitt R, Miller M. Deformable templates using large deformation kinematics. IEEE Trans Med Imag 1996;10:1435–1447.
1411 24. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes DJ. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Trans Med Imag 1999;18:712–721. 25. Hansen C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992;34:561–580. 26. Wang Y, Riederer SJ, Ehman RL. Respiratory motion of the heart: Kinematics and the implications for the spatial resolution in coronary imaging. Magn Reson Med 1995;33:713–719.
