|Home | About | Journals | Submit | Contact Us | Français|
Dynamic enhancement causes serious problems for registration of contrast enhanced breast MRI, due to variable uptakes of agent on different tissues or even same tissues in the breast. We present an iterative optimization algorithm to de-enhance the dynamic contrast-enhanced breast MRI and then register them for avoiding the effects of enhancement on image registration. In particular, the spatially varying enhancements are modeled by a Markov Random Field, and estimated by a locally smooth function with boundaries using a graph cut algorithm. The de-enhanced images are then registered by conventional B-spline based registration algorithm. These two steps benefit from each other and are repeated until the results converge. Experimental results show that our two-step registration algorithm performs much better than conventional mutual information based registration algorithm. Also, the effects of tumor shrinking in the conventional registration algorithms can be effectively avoided by our registration algorithm.
Image registration is a fundamental problem in medical imaging. Classical image registration techniques assume corresponding pixels in two images have consistent intensities or at least are strongly correlated . For instance, it is commonly assumed that if two pixels have the same intensity in one image, they will have the same intensity in the other image . However, this assumption is not always valid in many applications. For example, in dynamic contrast enhanced (DCE) magnetic resonance (MR) images, tumor and non-tumor tissues may have similar intensities at a pre-contrast image, but then appear very dissimilar at post-contrast images due to variable signal enhancements [3,2]. Moreover, even same tissue (such as tumor tissue)  can appear different at various contrast enhancement stages, due to different uptakes of the contrast agent. These spatially varying enhancements lead to serious problems for the conventional registration algorithms which assume consistent intensities or strongly correlated intensities in the two images under registration . One typical problem is that a tumor region can be significantly deformed during registration of DCE-MR images in order to maximize the mutual information as reported in .
Many research efforts have been proposed to design good correlation metrics in order to cancel out intensity inconsistencies across the images. Normalized mutual information (NMI) [7,3] has been used to account for the change of image intensity and contrast in registering DCE-MR images. However, NMI only works well in the regions with similar enhancement characteristics  and poorly handles the images with spatially varying enhancement like the DCE-MR breast images . Model-based approaches, such as the pharmacokinetic (PK) model , have also been used to characterize pixel-wise enhancement-time signal. In theory, one can combine the PK model with conventional registration algorithms to get better performance, provided that the temporal enhanced signals comply well with the PK model. In practice, the PK model alone is often insufficient to model the intensity variations. Although the PK model with addition of a noise term can potentially better represent temporal enhanced signals , the noise term alone is insufficient to model inconsistency of intensity enhancements between the two images.
In this paper, we present a different approach to solve the intensity inconsistency of DCE-MR images. Instead of finding the physical model that represent this inconsistency, we set out to find the intensity transformation between the images, i.e., estimating enhancements. Specifically, we calculate the ratio of the intensity for each tissue point over the two images. The resulting ratio image is called the enhancement map. For perfectly registered images, the enhancement map is simply the ratio between the reference and the target images. For time-lapse DCE-MR images, computing the enhancement map requires registration between the target and the reference image due to possible patient motion.
We model the enhancement map by a Markov Random Field, and estimate it by a locally smooth function with boundaries using a graph cut algorithm . Then, we use the estimated enhancement map to eliminate the enhancements between two images under registration, and further use a B-spline based registration method to register those de-enhanced images. These two steps are incorporated into an iterative optimization algorithm, and are repeated iteratively until convergence, similar to . Experimental results show that the elimination of spatially variable enhancements in the DCE-MR images can significantly improve the registration accuracy, compared to the conventional mutual information based registration algorithms .
Given two images I and J, our goal is to find the correspondence mapping T that warps every pixel from I to J. We use notion T (I) to represent the resulting image after applying T on I. An ideal registration should perfectly align T (I) with J. However, T (I) and J can be different due to intensity changes η between I and J, e.g., enhancement variations in DCE-MR breast images. In addition, both I and J can be corrupted by noise ε. Therefore, registration methods that directly measure similarity of pixels intensities between the images might fail.
Notice, if we know η in prior, we can simply cancel out the enhancement between I and J and then apply existing registration algorithms to recover T. Similarly, if we know T in prior, we can then solve η using T (I) and J. However, neither η nor T is known to us. Therefore, our goal is to simultaneously find the optimal T and η by maximizing a posteriori (MAP) given I and J:
We propose an iterative optimization algorithm to approximate the MAP of T and η. We start with a rough registration using existing algorithms. At each iteration, we first estimate the optimal η given the current estimation of T. We then re-compute the registration by cancelling out η in J. The iteration continues until the results converge.
Assume we have estimated the registration T from I to J, our goal is to compute the enhancement map η between I′ and J, where I′ = T (I). The enhancement is defined as
where s represent the pixels inside image region Ω.
Enhancement map has been widely used to analyze DCE-MR images in breast cancer diagnosis [12,4]. We assume the enhancements change spatially in a smooth way, with discontinuities appeared between different kinds of enhancements, e.g., tumor boundaries.
If the registration T is exact, then η can be directly computed using equation (2). In practice, T can be inaccurate and both images are corrupted by noise ε. Notice, the spatial smoothness and tumor boundary discontinuity of the enhancement map are very similar to those seen in disparity maps in computer vision. Recent research  have shown that fields satisfying such properties can be modeled as a Markov Random Field. Methods such as belief propagation or graph-cut are very efficient to solve the MAP of Markov Random Fields.
We start with an initial estimation of the enhancement map by computing C′ using equation (2). We then discretize the range of the enhancements into a finite number (e.g., 15) of levels, which are initially obtained by clustering with k-means algorithm. Our goal is to assign every pixel in C′ with an enhancement label and find the MAP of the enhancement MRF η. In our implementation, we modify the graph cut algorithm  to compute η that maximizes the labeling consistency and overall smoothness. The energy function we choose to minimize is defined as
E1 computes the difference between the assigned level value and enhancement computed by equation (2) for all pixels, which is defined as:
E2 corresponds to the smoothness term in standard graph cut algorithms. It measures the consistency of the enhancement levels between the neighboring pixels in :
where δ is a Kronecker delta function.
Finally, we introduce a new term E3 that further constrains the enhancement variations between pixels, where
The use of E3 guarantees that even if the neighboring pixels have different enhancement levels, their differences should be constrained within a small range. In our experiments, we find the use of E3 is critical to maintain spatial smoothness of the recovered enhancement map.
In equation (3), λ1 and λ2 are used for adjusting the relative importance of the three energy terms. They are both set to 1 in our experiments. We also find 15 levels clustered by k-means are sufficient to accurately approximate an enhancement map for a typical pair of DCE-MR breast images. This can be shown from the results of enhancement map estimated by our graph cut on a pair of typical enhanced images in Fig. 1.
Once we recover the enhancement map η between I′ and J, we can then “de-enhance” J to the same level of I′ using η. Notice Jde−enhanced should be approximately the same as I′. Therefore, we can reuse equation (2) to compute Jde−enhanced as
In Fig. 1, we use the estimated η to de-enhance the post-contrast MR breast image. The resulting image maintains similar spatial characteristics as the unenhanced one while preserving important features such as tumor boundaries.
We have developed a two-step iterative optimization algorithm to simultaneously recover the MAP registration T * and enhancement map η*. We start with an initial registration using the conventional B-spline based non-rigid registration on the pre-contrast and the post-contrast MR images . We then:
In step (a) of each iteration, we de-enhance the target image using equation (7), to obtain Jde−enhanced. We then estimate T′ by finding the optimal T′ = arg maxT (T |I, Jde−enhanced). This optimization can be easily computed using many existing registration algorithm. We use the B-spline based registration algorithm . In step (b), we use the estimated T′ to register I with J. We then use the graph-cut algorithm to find the optimal enhancement map η′.
Fig. 2 illustrates the intermediate results using our method on MR breast images for registration and enhancement map estimation. From the results, we can see the steps of enhancement estimation and registration help each other for achieving better results.
The performance of our enhancement estimation and image registration algorithm is extensively evaluated on both simulated data and real data as described by separate sections next.
For generating the simulated dataset, we first select a pre-contrast image of a real patient, such as the one shown in Fig. 1 (a). Second, we use the enhancements computed between this pre-contrast image and one post-contrast image of the same patient by equation (2) as a simulated enhancement map, to enhance our selected pre-contrast image and obtain a simulated post-contrast image. Third, we further use simulated motions, represented by B-Splines, to deform this simulated post-contrast image. Finally, we add Gaussian noise (with zero mean and standard deviation of 20) to this simulated post-contrast image.
The pre-contrast image and its final simulated post-contrast image are shown in Fig. 3 (a) and (b), respectively. To be clear, we only show the region around tumor. Notice that both simulated enhancement map and motions are pre-known.
We have compared our algorithm with the NMI-based method on this simulated dataset. Notice that both algorithms use the B-spline-based deformation representation  with 1 degree and 13 control points for each axis. To compensate for the intensity inconsistencies between the pre- and post-contrast images, the NMI-based approach estimates different intensity distributions in the two images while our algorithm first computes and removes the enhancement map and then directly uses the intensity difference metric. The results are shown in Fig. 3.
These results demonstrate at least three pieces of better performance of our algorithm. First, after removing the estimated enhancements (Fig. 3(c)), the simulated post-contrast image looks very similar to the original pre-contrast image (Fig. 3(a)). This shows the effectiveness of our graph cut algorithm in estimating enhancements. Second, the conventional registration algorithm produced errors especially at tumor as reflected by an enlarged tumor (Fig. 3(d)), while our algorithm still worked very well (Fig. 3(e)). Finally, the better registration result by our algorithm can also be observed as more similarity between the estimated (Fig. 3(h)) and the simulated enhancement maps (Fig. 3(f)), compared to that produced by the conventional registration algorithm (Fig. 3(g)). This is further demonstrated by the difference maps between the simulated and the estimated enhancements for our registration algorithm (Fig. 3(j)) and conventional registration algorithm (Fig. 3(i)), respectively.
Quantitatively, our registration algorithm is also better than the conventional registration algorithm. For example, by comparing the simulated motions with the estimated motions, our registration algorithm produced mean error of 0.44 pixel, standard deviation of 0.38 pixel, and maximum error of 2.4 pixels.
In contrast, the conventional registration algorithm produced mean error of 0.82 pixel, standard deviation of 0.57 pixel, and maximum error of 10.8 pixels.
Data acquired from eight patients are also used to further compare the performances of our method and conventional registration algorithm. For each DCE-MR image sequence, the pre-contrast and the last post-contrast images are selected for evaluation of both registration algorithms. Typical results on three subjects are shown in Fig. 4. Each row corresponds to one individual subject. The first column shows the pre-contrast images, while the second column shows the post-contrast images. The de-enhanced images by our method are shown in the third column. It can be observed that the de-enhanced images are very similar to the pre-contrast images in the first column, which indicates the effectiveness of our enhancement estimation method.
The difference maps between the warped pre-contrast images and the deenhanced images are shown in the fourth and the fifth columns for conventional registration method and our registration method, respectively. It can be observed that our method generally produces smaller difference maps, indicating better registration results.
Our algorithm completes the registration within 15 seconds on a PC of 1.2 GHz Intel Pentium M CPU and 2GB memory for images of a resolution 512 × 512 while the conventional algorithm takes about 3 seconds. For all the experiments demonstrated in this paper, our two-step optimization converges after 3 iterations.
We have presented an iterative two-step optimization algorithm to simultaneously estimate dynamic enhancements in DCE-MR images and eliminate them for robust registration. The estimation of dynamic enhancement map is particularly formulated by a MRF model, and completed by a graph cut algorithm. The experimental results on both real and simulated datasets show the relative accuracy of registration after de-enhancing the DCE-MR images by our method. Also, by comparing with conventional registration methods using NMI, our registration algorithm can effectively register the tumor regions across different images, thus potentially avoiding the effects of tumor shrinking or enlarging that commonly happen in conventional registration methods. In the future, we will test our algorithm on more datasets available in our institute, adapt it to 3D case, and re-formulate it for group-wise registration of multiple contrast-enhanced breast images.