Joint Reconstruction of Multiple Images and Motion in MRI Application to Free-Breathing Myocardial T2 Quantification
Abstract
Introduction
MRI examinations consist of acquiring multiple images which often contain a large amount of redundant information. This is because the anatomy is preserved provided the patient does not move significantly between the different imaging sequences. Exploiting these redundancies can be formalized as the joint reconstruction of these images using dedicated regularization schemes. Constraints that can take advantage of the anatomical redundancy can be based, for instance, on the co-occurrence of the spatial gradients (i.e. the gradients should be located in the same place in all images) or on joint statistics (i.e. joint entropy, mutual information etc…). The principle of joint image reconstruction has also been applied to multi-modal imaging such as the joint reconstruction of PET and MRI images
The application of this concept is particularly challenging in cardiovascular and body MRI applications. All images should be perfectly aligned in order to benefit best from the redundancy constraints. Organs may be misaligned due to involuntary motion, due to physiological changes (respiratory, cardiac, peristaltic motion etc…) or due to geometric distortions caused by hardware imperfections (e.g. in echo planar imaging sequences). In addition to correcting misalignments from image to image, it is desirable to correct for ghosting/blurring artifacts caused by motion during the k-space sampling (i.e. between the acquisition of different phase encoding steps), which can affect each image individually. These two types of motion will be termed inter-image and intra-image motion in the remainder. In the case of a single image acquisition, only intra-image motion correction needs to be considered. This can be formalized as the joint optimization of the image and a nonrigid motion model. Such an approach has also been applied to CT reconstruction
Myocardial T2 quantification is an example application that might benefit from a joint reconstruction framework. Myocardial T2 is altered in a number of cardiopathies including heart transplant rejection. Patient motion is commonly addressed by using a single-shot acquisition technique consisting of a breath-held T2-prepared balanced steady-state free precession (T2-prep bSSFP) sequence with three different echo times (TE). Such acquisitions are fast but have limited spatial resolution. Furthermore they are sensitive to artifacts because only three data points (three values of TE) are available for the two parameter exponential fit of the form aexp(−(TE/b)). The latter point can be addressed by covering a wider range of TE values. However the acquisition becomes significantly longer and has to be performed during free breathing. Such an approach, combined with nonrigid image registration, has been shown to improve the accuracy of T2 measurements . Finally in order to increase the spatial resolution, single-shot techniques have to be replaced by multi-shot ones. A black-blood fast spin echo sequence can be used in that case as . This sequence provides good contrast between blood and myocardium due to the blood signal cancellation.
In this paper a framework is proposed for the joint reconstruction of multiple images and motion, which is designed for a wide range of thoracic and abdominal MRI applications. An application to free-breathing T2quantification is presented with the data acquisition consisting of a series of 10 fast spin echo sequences with varying effective times of echo (TEeff), resulting in 10 images with varying contrast. The objectives of the study are (i) to evaluate the benefit of the proposed framework compared to the independent reconstruction of each image in terms of intra-image motion correction (ghosting/blurring artifact reduction) and inter-image motion correction (alignment) and (ii) to validate the method clinically by comparing the quantitative myocardial T2 values obtained from free-breathing to the conventional breath-held protocol in a population of heart-transplant patients.
Methods
MRI Experiments
Nine heart transplant patients were recruited for this study which was approved by the local ethics committee. Written informed consent was obtained from all patients. These patients were scanned with a 1.5T SignaHDxt system (General Electric Healthcare, Milwaukee, USA). A myocardial T2 mapping protocol was used which has been validated for detecting and predicting allograft rejection. The sequence was an ECG-triggered black-blood fast spin echo (FSE) and was repeated 10 times with different effective times of echo (TEeff): 10, 15, 20, 30, 40, 50, 55, 60, 70 and 80 ms. One short-axis slice was chosen in the middle of the left ventricle and each of the 10 images was acquired during a breath-hold period (gold standard, i.e. 10 separate breath-hold periods were necessary to obtain a T2 dataset) and during free breathing (for MI-GRICS reconstructions). The matrix size was 256 × 160 for breath-held sequences and 256 × 256 for free-breathing sequences. Images were eventually reconstructed to 512 × 512. Other parameters were: 10 mm slice thickness, typical field of view of 320 mm, pixel size of 1.25×2 mm2 (interpolated to 0.625×0.625 mm2), echo train length=16, number of excitations Nex=1, trigger time set for acquisition in end-diastole (depending on the heart rate), repetition time TR=2 heart beats. The MR system’s body coil was used for transmission and reception in breath-held acquisitions. For the free-breathing scans an 8-element cardiac coil array was used for signal reception. The body coil is used for the clinical diagnostic sequence in our center because the intensity weighting obtained with surface coils depends on the vendor’s calibration and coil combination algorithms which are not always sufficiently documented. For the free-breathing sequences, we used our own calculation of coil sensitivity maps from the raw data of the vendor’s calibration sequence (i.e. dividing each surface coil image by the body coil image). The coil combination was embedded in our reconstruction framework through the use of these sensitivity maps in the acquisition model.
Physiological signals were collected using an in-house real-time system. Respiration signals were provided by two pneumatic belts placed on the patient’s chest and abdomen. The system also recorded acquisition windows so that it was possible to relate each acquired k-space line to the corresponding motion signal samples. The two breathing signals as well as their first-order temporal derivatives were used as the sk(t) signals (i.e. Ns=4 signals) needed by GRICS and MI-GRICS reconstructions.
Free-breathing data were reconstructed using the conventional GRICS following (independent reconstruction of each image and motion model). Here conventional GRICS is expected to yield suboptimal results since we used Nex=1 (severe ill-conditioning) as explained in the Theory section. They were also reconstructed using multi-image GRICS following : (i) firstly we turned off the inter-image registration (i.e. T(u~i)=Identity) and the joint image regularization (i.e. R=Identity), i.e. the only difference with conventional GRICS was the use of a shared motion model; (ii) then we only turned off the inter-image registration (i.e. using T(u~i=Identity)); (iii) finally we used the full version providing the complete inter- and intra-motion correction. The latter three reconstructions will be termed MI-GRICS1, MI-GRICS2 and MI-GRICS3 respectively. For the nonrigid registration step required by MI-GRICS3 we used a technique which was extensively validated for T2 mapping in heart transplant patients]. It consists of a pre-processing by histogram matching in order to reduce intensity differences in the image series; then a sum-of-squared-differences similarity criterion is optimized using a multi-resolution Gauss-Newton scheme.
Implementation of the Reconstruction
Reconstructions were implemented in C++ language. The reconstruction code was run on a cluster of workstations comprised of 16 machines on an InfiniBand network. Each machine was equipped with Intel Xeon Quad Core CPU (X550, 2.67 GHz) and 24 GB RAM. Our reconstruction library supports both multi-node parallelization (MPI standard) and multicore parallelization (OpenMP standard). The most computationally intensive part lies in the conjugate gradient linear solvers . The application of the underlying operators (EHE and JHJ) comprises an outer loop on motion states and an inner loop on coil receivers, which were distributed respectively on the 16 nodes (MPI parallelization) and on 8 cores (OpenMP parallelization).
The regularization tuning parameters were set as follows: λ=10−8, μ=0.005, v=0.01 (λ,μ,v were set relative to the norm of the right hand-side of the linear system in which they are involved). They were optimized empirically on one dataset and then the same values were used for all other datasets.
Assessment of the Commutation Assumption Used To Derive MI-GRICS
In order to derive and the subsequent Algorithm (2), the assumption has been made that, in our application, two transformation operators can be commuted with little error ). The impact of this assumption was quantified retrospectively. Images and motion fields generated by MI-GRICS3 reconstruction were used for testing the commutation assumption1. For each patient dataset the displacement fields with maximal amplitude were extracted from both the motion model and the inter-image motion fields. Those displacement fields, U(1)max and U(2)max respectively, were then applied to each of the 10 images (corresponding to the 10 echo times) reconstructed by MI-GRICS3. After application of TU(1)maxTU(2)maxor TU(2)maxTU(1)max, the normalized root mean squared errors (NRMSE) between the resulting images were computed. Among the 10 images, only the worst case corresponding to maximal NRMSE was recorded. We also tested the commutation assumption when the amplitude of those maximal displacements was divided by 2 and by 4.
Assessment of Intra-Image Motion Correction (blurring/ghosting artifacts)
The quality of intra-image motion correction was assessed using the gradient entropy as an image quality metric, as expressed by:
Here ρ=[ρ1…ρNi]Tis the entire reconstructed image series, ∇ is the spatial gradient and H is the Shannon entropy based on the probability P of pixel value occurrences. The entropy of an information channel is increased when more uncertainty is associated with the information source and H(X)is maximal when all intensity values in X are equally probable. Motion artifacts (blurring and ghosting) come from a dispersion of the signal in the image. This results in a smoothing in the image histogram or in the gradient image histogram (the gradient image highlights the details/features), and thereby an increased entropy. We computed Qintra for the different reconstructions (GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3) and tested for statistical differences with a Wilcoxon signed rank test, setting the significance level to 5%.
Assessment of Inter-Image Motion Correction (registration)
Inter-image motion correction was first assessed qualitatively by visualizing a movie of the reconstructed image series. An average image of the spatial gradients (i.e. averaged over the 10 echo times) was also used to represent in a single static image the amount of residual misalignment in the reconstructed series.
An objective metric Qinterwas also computed to quantify the quality of inter-image motion correction, based on the normalized mutual information (NMI) between the first image of the reconstructed series and all other images:
The NMI was chosen because it is one of the most popular similarity criterions used in medical image registration and is fairly robust to contrast changes that occur in multi-contrast or multi-modal imaging. It is fair to use this metric for validation in this study because our optimization criterions are not based on NMI (our motion model optimizes the fidelity of the acquisition model to the data; our registration optimizes the sum-of-squared-differences of the images after histogram matching). Wilcoxon signed rank tests were performed again to compare the different reconstructions.
Quantitative Myocardial T2 Analysis
Two of the nine patients were excluded from the quantitative T2 analysis. The trigger delay for these two patients was improperly set (with echo trains falling in the systolic contraction) leading to unreliable signal intensity.
T2 quantification was performed in the remaining seven patients. The breath-held image series was considered to be the reference standard and was processed for T2 quantification in 3 steps: (i) nonrigid registration of the images using the same technique as described earlier, which was previously validated for myocardial T2 quantification; (ii) manual segmentation of endocardium and epicardium (J.-M. E., 30 years of experience in experimental and clinical T2 quantification) and separation into 6 segments according to American Heart Association (AHA) guidelines; (iii) segment-wise calculation of T2 values by an exponential fitting of the averaged segmental signals using a Levenberg-Marquardt algorithm. For the motion-corrected images (GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3 reconstructions) the post-processing only consisted of steps (ii) and (iii).
Segment-wise T2 values obtained from free-breathing data (uncorrected, GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3 reconstructions) were compared with breath-held data. We first compared T2 values from all myocardial segments of our patient population (N=42 segments). Then we also compared T2 values from myocardial segments where the exponential fit was good, as defined by the coefficient of determination R2>0.97. This criterion has been suggested in previous studies as a quality control index for excluding low quality datasets. The data were tested for statistical significance with Wilcoxon signed rank tests again.
Finally we compared the diagnostic accuracy of the proposed reconstruction (MI-GRICS3) with the reference technique (breath-hold): in our institution, suspicion of allograft rejection is reported if the average T2 of the 6 myocardial segments is above 62 ms.
Discussion
Joint Reconstruction of Multiple Images and Motion in MRI Application to Free-Breathing Myocardial T2 Quantification,A framework has been proposed for the joint reconstruction of multiple images and motion and has been applied to free-breathing myocardial T2quantification in heart transplant patients. Our results show that the joint reconstruction of multiple images and motion (MI-GRICS1, MI-GRICS2 and MI-GRICS3) has a benefit over the joint reconstruction of a single image and motion (GRICS). Specifically, the analysis of the intra- and inter-image motion correction metrics revealed the following significant improvements: (i) the sharing of the motion model improved inter-image correction over conventional GRICS (improved Qinter) which might be explained by a more stable fit of the motion model due to dimensionality reduction; (ii) the joint image regularization improved intra-image correction (improved Qintra) which might be explained by the more stable image reconstruction step; (iii) the added inter-image registration improved inter-image correction further (improved Qinter). The latter improvement was implemented with a limited extra computational cost of 13% compared to MI-GRICS2, including 9% for the registration steps alone. This was possible due to the commutation assumption described and we have checked that this assumption did not introduce large errors in our acquisition model, even in the worst case motion amplitudes. It should be noted that MI-GRICS3 is different from applying MI-GRICS2 followed by registration. This is because the inter-image motion fields ensure that the images are optimally aligned when applying the multi-image regularizer. In our results, although this did not translate into a statistically significant improvement of the gradient entropy metric, a reduction of the artifacts was visible in some cases .
The results also show the feasibility of using the technique for myocardial T2quantification in heart transplant patients. Motion-compensation strategies have been proposed previously for such applications but have been restricted to the registration of relatively low-resolution single-shot images. The proposed method should be particularly useful in patients who have difficulties holding their breath repeatedly and when lengthy high-resolution (multi-shot) images are desired. Another benefit is that the images generated by the reconstruction are already registered and are readily available for segmentation. Importantly segment-wise myocardial T2values were found to be in good agreement with the reference breath-held images. The bias (62.5 ms for MI-GRICS3 against 62.2 ms for the reference) was small compared to the inter-observer variability of T2 measurement with this technique which was evaluated to be 2.7 ms in Ref. The diagnosis with MI-GRICS3 was in agreement with that of the reference protocol in all seven patients. The acquisition time was nearly the same for the free-breathing and for the breath-held protocols although the k-space size was increased in the free-breathing experiments (from 160 to 256 in the phase encoding dimension). This was possible because the breath-held protocol comprises many idle periods so that the patient can rest between successive apneas, unlike the free-breathing protocol. Finally the reconstruction time for a T2 dataset was 9 min which is close to the scan time, making the technique useable clinically.
In this study we used a generally applicable constraint to take advantage of anatomical redundancies in the images based on the spatial gradients (i.e. gradient co-occurrence constraint). Other redundancy constraints would be of interest, in particular those based on joint statistics of the images such as the joint entropy or the mutual information. As an alternative in specific applications it might be possible to impose an analytical MR signal model (exponential models for T2 or T1 mapping, pharmacokinetic models for perfusion imaging etc…), mathematical models such as B-splines], low-rank signal evolution models or sparsity constraints in other domains such as the x-f space or x-p space with p the varying acquisition parameter space. Other models based on the singular value decomposition have also been proposed recently for accelerated parameter mappings. Here the multi-image reconstruction framework was applied to multi-contrast MRI datasets. It might also be adapted to dynamic imaging (with or without significant contrast changes) such as perfusion imaging or cardiac cine imaging and is therefore conceptually close to the recent work. Another possible use of the framework would be in the context of multi-modal imaging. Finally other motion sensors might be used as surrogates of (or in combination with) the external sensors used here. In particular, partial MRI data (navigators or self-gated approaches) might give signals better correlated with the motion of organs and improve the results further.
The joint reconstruction of image and motion, as defined in (1), (4), (7), is a non-convex optimization procedure. The proposed approach relies on a multi-scale implementation and a linearization with respect to the motion parameters . The latter step involved the optic flow equation which is known to be ill-posed (one equation with two unknowns, known as the aperture problem). Occlusions are also known to render motion estimation difficult. Here occlusions might result from large motion as structures move in and out of the imaging plane. This might cause the algorithms to be trapped in a local minimum. When applied in 2D, the method may fail when the slice thickness is small compared to the through-plane motion. In practice the clinical slice thicknesses of 8–10 mm used in cardiac MRI generally give good results. Otherwise the method might be applied in 3D to resolve through-plane motion.
One limitation of the method is the choice of the registration technique which might not be suitable for all applications. The use of a mutual information based similarity criterion might improve the registration step slightly. More sophisticated techniques might still be necessary in case of complex contrast changes from one image to another, as can occur for instance in first-pass perfusion imaging. Another limitation of this study is the small patient population.
In conclusion we have presented a joint optimization technique that allows reconstructing N images together with both intra-image and inter-image motion models. The method allows taking advantage of anatomical redundancies in the image dataset despite the complex organ deformations during the acquisition. The feasibility of the technique was demonstrated with free-breathing myocardial T2 quantification in heart transplant patients providing equivalent diagnostic outcome compared to the reference technique.
