Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) in the liver enables quantification of hepatic perfusion for assessment of tumor response to radiation therapy, as well as normal tissue damage (1–5). One of the challenges in liver DCE-MRI is respiratory motion during acquisition, which can cause blurring and artifacts in the images and intensity variations in the dynamic curves, and subsequently, errors in perfusion quantification using pharmacokinetic (PK) models. To reduce respiratory motion-induced degradation of liver images during DCE-MRI acquisition, a breathing control paradigm has been proposed, during which the patient is instructed to hold their breath multiple times (taking breaths in between) (5). The dynamic image volumes acquired during the 2–3 minutes of contrast uptake are still misaligned because of poor reproducibility of the liver position during multiple breath holds, loss of control during breath holding, and occasionally gross movement of the patient. Therefore, a voxel cannot be assumed to cover the same piece of tissue throughout the whole examination, but, instead, it covers different pieces at different points in time. Therefore, if not corrected for, breathing and gastrointestinal motions can induce artifacts, such as intensity spikes and wide bumps, in the contrast uptake curves. To estimate the arterial and portal venous perfusion in the liver, a sufficient number of reliable samples of the uptake curve must be collected to fully describe the contrast uptake dynamics. The shape of the precontrast plateau and the beginning of the intensity ramp following contrast agent administration can have a large effect on the estimated perfusion parameters, in particular on hepatic arterial perfusion. Unwanted artifacts introduced into the uptake curve by patient motion reduce the number of reliable samples, increasing the risk of inaccurate PK model parameter estimation.
To correct for image misalignment in the DCE series, image registration has been used to align the magnetic resonance (MR) images at different points in time or phases (1, 6–8). However, registration of the DCE-MRI phases has some unique challenges, for example, a large number of dynamic image volumes, low spatial resolution, and low contrast-to-noise ratio that are tradeoffs for the high temporal resolution. In addition, the image contrast is changed completely over the whole course from precontrast to rapid contrast uptake and to full contrast enhancement. Combined with all these effects, the dynamic contrast agent-induced intensity changes pose a difficult problem even to registration methods using multimodality image-to-image metrics such as mutual information (9–11).
A number of methods have been developed to tackle the DCE-MRI registration problem. Some methods attempt to register each phase to a reference phase selected among all phases. Such methods use either region-limited rigid registration (12) or deformable registration (9, 10). Deformable registration may yield incorrect volume changes around contrast-enhanced lesions unless the method is adapted to explicitly counteract such changes (11). As an alternative to reference phase-based registration, image intensity changes can be modeled using a PK model and the model incorporated into the registration algorithm (13–16). Such PK-based methods can be slow because of the need for repeated PK model parameter estimation. For a faster procedure, principal component analysis (PCA) can be used to generate reference images for registration (17). However, PCA may not be able to adequately describe the spatially varying time–intensity curves. Therefore, the phases could alternatively be decomposed into a low-rank and a sparse image component in a procedure called robust PCA that can be combined with image registration (18, 19). Robust PCA can also be used to combine registration with reconstruction, but this requires access to raw data from the MRI scanner (20). However, when comparing the 3 representative registration methods (21), 1 PK model-based, 1 PCA-based, and 1 sequentially registering each phase to its preceding phase, the sequential registration method performed better (average error 14.7% in estimated parameters) than the more complex iterative, PCA-based and PK model-based methods (average errors: 39.5% and 39.2%, respectively). Nevertheless, all image registration methods reduced errors in the estimated parameters in the tumors in comparison with without image registration methods.
The performance of the sequential method can be understood by considering that the changes between neighboring phases are small. Therefore, registration has a greater chance of success. However, misregistration errors are accumulated and propagated from the earlier phases to the later phases. If registration at the early phase fails, all following phases will be incorrectly aligned with the first phase.
Inspired by the findings in Rajaraman S et al.'s study (21), and taking into account the aforementioned drawback of the sequential registration procedure, we introduce an overdetermined system for achieving robust registration for the DCE image series. In this procedure, rather than registering each phase to a reference phase or to a preceding phase, each phase is registered to every other phase. This produces an overdetermined system of transform equations that can be solved using robust statistical methods to reject outliers corresponding to failed registrations. This new method is evaluated on a set of 100 liver DCE-MRI examinations, and the results are compared with those of a conventional registration method where all phases are registered to 1 postcontrast reference phase. A robust registration method for liver DCE-MRI would greatly simplify the workflow of clinical liver perfusion studies with a large number of participating patients by eliminating the need for manual intervention following faulty registration.
Materials and Methods
Under institutional review board approval, 100 DCE-MRI examinations of the liver from 35 patients were included in this study (women, 8; men, 27; age at examination, 51–83 years; number of examinations per patient, 1–4). The patients were imaged for about 3 minutes using a repeated breath-hold paradigm (5). For this paradigm, the patients were instructed to initially hold their breath for as long as they could and then hold their breath repeatedly with a single deep inhalation in between each breath hold. Images acquired during deep inhalations showed large liver movement and severe motion blur. Therefore, these were excluded from further analysis. Dynamic MR imaging was started at the beginning of the first breath hold and a gadobenate-dimeglumine-based contrast agent (MultiHance, Bracco S.p.A., Italy) was administered intravenously.
A 3 T MRI scanner (Magnetom Skyra, Siemens Medical Systems, Erlangen, Germany) was used for imaging. A 3-dimensional gradient echo sequence that speeds up acquisition using stochastic sampling with view sharing (time-resolved angiography with stochastic trajectories; TWIST) was used for 10 patients, and a sequence that enhances speed through partial Fourier reconstruction (volume interpolated breath hold examination; VIBE) was used for 25 patients. Sequence parameters varied among the patients (Table 1).
During a breath hold, the liver should ideally be still. Unfortunately, partial inhalation during a breath hold, imperfect reproducibility of the exhale state, and changes to patient posture result in liver movement between images. These movements, however, are smaller than those observed during free breathing. For the spatial resolution, signal-to-noise ratio, and artifact level seen in this study, we, therefore, chose to use rigid registration to correct for the observed range of motions.
Before registration, a number of preprocessing steps were performed. These were the same for the robust registration procedure (section Robust Registration) and the conventional procedure (section Conventional Registration) to which it was compared.
Reference Phase Selection.
The conventional procedure involved alignment of all phases to 1 reference phase. For each examination, this reference phase was selected among all phases as the phase with the smallest sum of mean absolute distances to the other phases in intensity, which is calculated as follows:
Region-Limited Image Registration.
Because we wish to correct for the motion of the liver with a rigid transform and the rest of the body inside the field of view (FOV) does not move in the same way as the liver, region-limited registration was used. For this purpose, a region of interest (ROI) was drawn around the liver on the reference phase. This ROI was drawn snugly along the edge of the liver. For the region-limited registration, it is advantageous to include a margin outside the liver ROI. This margin covers the contrast-rich interface between the liver and the surrounding tissue, and it helps to drive the registration algorithm. The margin was set to 15 mm for all registrations.
Bias Field Correction.
MR images commonly contain a spatial intensity variation that is not the result of intrinsic tissue properties but of the position of the tissue with respect to the transmission and receiver radiofrequency coils. This “bias field” can cause misregistration if it is the dominant feature in an image. Fortunately, the bias field typically varies slowly with respect to the spatial position. The effect of the bias field on the images used for registration was reduced by dividing each phase by a filtered version of itself. The modified phases used for registration were then given by the following equation:
The robust registration procedure for each DCE time series was divided into 5 stages—3 translational registration stages and 2 rigid registration stages with both translation and rotation. For the robust procedure, the workflow for each stage is illustrated in Figure 1. Either each possible pair of phases or each pair in a subset of all possible pairs was registered. Thus, each registered pair gave rise to a rigid transform describing the relative position of the 2 phases of the pair. For the registration of each pair, the first phase of the pair was used as a fixed image and the second as a moving image. The dilated liver ROI was applied as a mask to the fixed image to restrict the domain for which the registration metric was evaluated. Because of this asymmetry, the registration of pair , where is the index of the fixed image and is the index of the moving image, did not produce the same transform as the registration of pair .
After the registration of all pairs of 1 of the steps, the resulting transforms were used to produce a set of consistent transforms as described in the next section. These transforms were then applied to the phases before proceeding to the next stage.
The registration of the 2 phases of a pair was performed using the open source registration software Plastimatch, version 1.6.0-beta (http://plastimatch.org). The settings of the translational and rigid registration steps are shown in Table 2. For the robust procedure, the first 2 translational steps only registered a subset with pairs selected randomly among the possible pairs. These 2 steps helped roughly align the liver ROI for all the phases in preparation for the full registrations performed by the third through the fifth stage.
The registration of a pair of phases and results in a rigid transform described by a rotation matrix and a translation vector . Thus, given a structure found in phase at spatial coordinates , the corresponding structure, as indicated by the registration, can be found at position in phase j.
Let us assume that there is a transform given by rotation matrix and translation vector , which transforms the coordinates of any structure in the reference phase into the coordinates of the same structure in phase . If all of the coordinates for the voxels inside the liver ROI are inserted into the columns of a matrix , the corresponding points in phase are given by , where . However, if the registration of phase pair is successful, the same coordinates are also given by the application of the true transform from the reference phase to phase followed by the registration transform from to such that: such that:(1), 1 for each registered phase pair, and up to are registered for a given stage. However, there are only true transforms given by and . Therefore, the system in equation (1) is overdetermined if . Because the registration of a phase pair may not produce the true transform but only an estimate, equation (1) is not valid for all registered pairs in general. Therefore, we must choose in what sense to solve the system in equation (1). If we define equation (2) as follows: (3) would produce a good solution for and . However, because some of the phase pair registrations could fail and produce severely erroneous transforms, and , the sum in equation (3) could be dominated by a few large terms that would thwart the accurate estimation of and . Therefore, to lessen the impact of outlier registration transforms, and can be robustly estimated by minimizing the following equation: (5) is iterated until convergence. To avoid division by 0 in equation (6), a nonzero regularization term is needed. As is the root mean squared difference of the liver voxel coordinates produced by the left-hand side (LHS) and right-hand side (RHS) of equation (1), should be smaller than the best expected registration accuracy. For this work, was selected as 0.01 mm.
Each iteration of equation (5) amounts to solving a constrained, weighted-least-squares problem, where the constraint is that all must be rotation matrices. Because the minimization problem in equation (5) is a rigid-body problem, it can be broken up into 2 independent problems, 1 problem for the translation of the liver center of mass and 1 problem for the rotation of the liver about its center of mass. This can be seen by noting that:(5) can be written as the following 2 problems: (14) is a problem to find the rotation matrices only, whereas equation (15) is a problem to find the transformed translation vectors only and . The problem in equation (15) is a system of linear equations of the following form: (14) is used to fix to a constant by setting The problem in equation (14) can similarly be written as a system of matrix equations as follows: (18) can be decomposed into its singular value decomposition to produce a set of the following iteration equations: (20) is selected to make the determinant of positive and equation (21) is needed to ensure that . The stopping criteria for equations (19) to (21) are as follows:
Empirically, the outer recursive sequence given by equation (5) was observed to converge within 20 iterations for a few test cases. Therefore, 20 iterations of equation (5) were used for all examinations.
A conventional registration procedure was used as a benchmark for the robust registration. This conventional procedure registered each phase to the reference phase using the same 5 registration stages as for the robust procedure. However, phases were not registered to anything but the reference phase, and no processing of the resulting transforms was performed except the combination of the transforms from the 5 stages for each phase. The same ROI was used to limit the registration region, and the settings in Table 2 were used for alignment.
The effect of registration on the mean signal intensity versus time for a 1-cm3 cubic ROI placed 1 cm below the dome of the liver for 1 patient is illustrated in Figure 2. The spikes seen in the intensity for the nonregistered time series are caused by tissue moving in and out of the ROI due to breathing.
Given the appearance of the signal variations seen in Figure 2, a suitable metric for the quality of the DCE-MRI time series should reflect how smooth the time–intensity curves are for the voxels inside the liver ROI. The temporal resolution of the time series as compared with the rise time of the time–intensity curve should also be reflected by the alignment quality metric. If only a few samples are collected during the rise time of a curve, the impact of 1 erroneous sample will be greater than that if a larger number of samples are collected. The variation of temporal sampling density across examinations can be taken into account by transforming the time of a given phase into a normalized time as follows:
Furthermore, errors in the baseline will be more troublesome for parameter estimation than had they appeared in the postcontrast plateau. In light of these considerations, we have chosen to use an approximation of the local spike area divided by the mean temporally local signal intensity to calculate the metric (see right panel of Figure 2 for illustration).
If the normalized temporal spacing between phases is , an approximation of the time-normalized, intensity-normalized spike area (TISA) metric, , at phase and position is then given by the following equation:(25) for a phase, equation (25) cannot be evaluated for phases adjacent to others omitted because of deep inhalations during the examination.
After the consistent registration solution has been found from the overdetermined system of transform equations by minimizing equation (4), the individual norms show how well a particular transform agrees with the overall solution. This is illustrated in Figure 3 for the last registration stage of 1 examination. The small dark square in the lower left corner of the figure represents the registration errors among the precontrast phases. The larger dark square in the upper right corner represent the errors among the postcontrast phases. These 2 groups have relatively small deviations from the consistent solution compared with the errors found for registrations between pre- and postcontrast phases as represented by the brighter rectangles in the upper left and lower right corners. This pattern was seen for a majority of the examinations registered and indicates that registrations between pre- and postcontrast phases are more likely to produce transforms that do not agree with the consistent solution.
By sampling the intensity along a 1×1-cm-thick column extending 20 cm above and 20 cm below the liver dome, the effect of registration on the position of the liver can be illustrated. Figure 4 shows the intensity along such a column for the phases of 1 examination along with the time–intensity curves for a 1-cm3 cubic ROI placed 1 cm below the liver dome. The inferosuperior intensity profile exhibits motion before registration. After conventional registration, the motion is reduced for most phases, but it fails for a few precontrast phases. The intensity profile produced by the robust registration method has no such obviously failed registrations. Failed registrations like the one seen in Figure 4C–D were not found among any of the registrations produced by the robust method.
The mean TISA, , for all examinations is shown with respect to the normalized time in Figure 5B. As a reference to help interpret the normalized time, , the mean of all liver time–intensity curves is shown in Figure 5A. This figure further illustrates where a particular value of corresponds to the precontrast region (), the intensity ramp (, or the postcontrast region ). The TISA curves are consistently higher in the precontrast region compared with those in the postcontrast region. This is a result of the lower mean intensity of the phases in the precontrast region. In addition, at the beginning of the intensity ramp (around , the conventional registration procedure performs almost as poorly as no registration, as seen by the close proximity of their TISA curves (green and red). Confidence intervals for the paired TISA differences between the robust method, the conventional method, and no registration are given in Table 3. The robust method is significantly better than no registration for the whole time interval (paired t test, P < .001). It is also significantly better than the conventional method for the precontrast and early contrast regions (P < .001).
In this study, we have implemented an overdetermined registration procedure for liver DCE-MRI examinations. We found that our procedure significantly reduced the image registration errors in the DCE time curves, particularly in the precontrast and early contrast uptake phases, compared with a conventional reference-phase-based registration method. Registration errors in the precontrast baseline could profoundly affect the calculation of signal intensity changes and the conversion to contrast concentration, whereas errors in the early contrast uptake curve propagate into derived parameters that are sensitive to the early curve dynamics, for example, hepatic arterial perfusion in liver DCE-MRI. The improved DCE curves describe the contrast agent uptake dynamics more correctly, and thereby, provide a better foundation for PK model parameter estimation.
DCE images acquired at different parts of the dynamic curves have different sensitivities to image registration errors. We found that registering precontrast phases to postcontrast phases is more likely to result in registration errors than registering within the precontrast phases or within the postcontrast phases. The conventional reference-based registration method likely results in a higher misregistration rate for the precontrast phases, as well as the early contrast uptake phases, if the postcontrast phase is chosen as reference. As an effect of this phenomenon, spikes can appear in or around the baseline plateau in the time–intensity curves. The parameters of a PK model may still be possible to estimate with such spikes present by using robust regression methods. However, each registration failure reduces the number of reliable samples. Even worse, if sample points are lost on the intensity ramp by removing the misregistration points, the remaining samples may not be able to correctly describe the dynamics of the curve, resulting in faulty PK parameter estimates. Some parameters are highly sensitive to the shape of the initial ramp of the uptake curve, for instance, hepatic arterial perfusion in the liver. The robust registration procedure presented reduces the number of failed registrations in the precontrast region by rejecting them as outliers in an overdetermined system of transform equations and thereby provides a greater number of reliable samples for parameter estimation.
For the postcontrast plateau, the robust registration method and the conventional registration method perform almost equally well. This is a result of the higher signal-to-noise ratio and more stable contrast for the phases of this region. Both registration methods produce smoother curves than no registration. The robust method is seen to have a slightly lower TISA than the conventional method. This is because the robust method takes the mutual information of all pairs of phases into account when estimating the transform. Therefore, the robust method acts as a more efficient estimator for the true transforms than does the conventional method.
In this work, we performed each registration separately and then combined the resulting transforms using a method that rejects outliers. Another possible approach would be to maximize the sum of the mutual information metrics of all pairs of phases using strategies similar to those reported by others (22–24). Such an approach might improve the registration of DCE-MRI further. In addition, although the method presented in this paper is only applicable to rigid registration, a method based on the maximization of pairwise mutual information could aid in improving the robustness of deformable registration for DCE-MRI. Robust deformable registration could allow free breathing examinations, which would relieve subjects from the burden of repeatedly holding their breath.
The robust method produced a consistent registration quality across all examinations, without failed alignments such as those observed for the conventional method. It therefore reduces the need for manual intervention following registration and helps to simplify the workflow of clinical studies with many subjects.
A method has been presented that registers every possible pair of phases from a DCE-MRI examination and derives a final set of transforms from them using minimization. The method is robust, insofar as to be able to reject such failed registrations that can appear when registering precontrast to postcontrast phases. The robust registration method improves the smoothness of the resulting time–intensity curves of voxels in the liver by eliminating spurious intensity spikes induced by motion or failed registration. Therefore, a baseline with reduced bias and a larger number of reliable samples in the intensity ramp are made available for the PK model parameter estimation. The consistent registration quality produced for all examinations studied shows that the method could simplify the workflow of clinical studies with many patients by eliminating the need for manual intervention following registration.