and shape prior information in a variational framework. The shape registration is considered the backbone of the approach where more general transformations ...
A Kidney Segmentation Approach from DCE-MRI Using Level Sets H. Abdelmunim Computational Biomedicine Lab Department of Computer Science, University of Houston Houston, TX, USA Aly A. Farag and W. Miller Computer Vision and Image Processing Lab Department of Electrical and Computer Engineering, University of Louisville Louisville, KY, USA Mohamed AboelGhar Medical School Mansoura University Mansoura, Egypy 1. Introduction
Abstract Acute rejection is the most common reason of graft failure after kidney transplantation, and early detection is crucial to survive the transplanted kidney function [1]. Automatic classification of normal and acute rejection transplants from Dynamic Contrast Enhanced Magnetic Resonance Imaging (DCEMRI), is of great importance. Kidney segmentation is the first step for such classification. The image intensity inside the kidney is used as an indication of failure/success. Differentiating between different cases cases is implemented by comparing subsequential kidney scans signals. So, this process is mainly dependent on segmentation. This paper introduces a new shape-based segmentation approach based on level sets. Training shapes are collected from different real data sets to represent the shape variations. Signed distance functions are used to represent these shapes. The methodology incorporates image and shape prior information in a variational framework. The shape registration is considered the backbone of the approach where more general transformations can be used to handle the process. We introduce a novel shape dissimilarity measure that enables the use of different (inhomogeneous) scales. The approach gives successful results compared with other techniques restricted to transformations with homogeneous scales. Results for segmenting kidney images will be illustrated and compared with other approaches to show the efficiency of the proposed technique.
In the United States, approximately 12000 renal transplants are performed annually [2], and considering the limited supply of donor organs, every effort is made to salvage the transplanted kidney [3]. Currently, the diagnosis of rejection is done via biopsy which has the downside effect of subjecting the patients to risks such as bleeding and infections. Moreover, the relatively small needle biopsies may lead to over or underestimation of the extent of inflammation in the entire graft [4]. Therefore, a noninvasive and repeatable technique is not only helpful but also needed in the diagnosis of acute renal rejection. In DCE-MRI, a contrast agent called Gd-DTPA is injected into the bloodstream, and as it perfuses into the organ, the kidneys are imaged rapidly and repeatedly. During the perfusion, Gd-DTPA causes a change in the relaxation times of the tissue and creates a contrast change in the images. As a result, the patterns of the contrast change gives functional information, while MRI provides good anatomical information which helps in distinguishing the diseases that affect different parts of the kidneys. However, even with an imaging technique like DCE-MRI, there are several problems such as, (i) the spatial resolution of the dynamic MR images is low due to fast scanning, (ii) the images suffer from the motion induced by the patient breathing which necessitates advanced registration techniques, (iii) the intensity of the kidney changes non-uniformly as the contrast agent perfuses into the cortex which complicates the segmentation procedures (see Fig. 1). To the best of our knowledge, there has been limited 1
978-1-4244-2340-8/08/$25.00 ©2008 IEEE
certain shape, we can define the following level set function φ to be the minimum Euclidean distance between the point X and the shape boundary. It is positive inside, negative outside, and zero on the shape boundary. Such representation can account for local deformations that are not visible for iso-contours that far away from the original shape and geometrical features of the shape can also be derived naturally from this representation.
3. Level Sets and Intensity Segmentation Within the level set formalism [8], the evolving curve is a propagating front embedded as the zero level of a 3D scalar function φ(x, y, t). The continuous change of φ can be described by the partial differential equation:
Figure 1. Different DCEMRI kidney images showing the challenges of the segmentation problem.
work on the dynamic MRI to overcome the problems of registration and segmentation. For the registration problem, Gerig et al. [5] proposed, using Hough transform, to register the edges in an image to the edges of a mask and Giele et al. [6] introduced a phase difference movement detection method to correct for kidney displacements. Both of these studies required building a mask manually by drawing the kidney contour on a 2D DCE-MRI image, followed by the registration of the time frames to this mask. Most of these efforts used healthy transplants in the image analysis, and edge detection algorithms were sufficient. However; in the case of acute rejection patients, the uptake of the contrast agent is decreased, so edge detection fails in giving connected contours. For this reason, we consider the combined usage of gray level and prior shape information to give better results.
2. Shape Modeling by Level Sets Shape representation is the main task of analysis of shapes. The selection of such representation is very important in several computer vision and medical applications such as registration and segmentation. There are several ways described in [7]. Although some of these ways are powerful enough to capture local deformations, they require a large number of parameters to deal with important shape deformations and have problems with changing topology. An emerging way to represent shapes can be derived using level sets [8]. This representation is invariant to translation and rotation. Given a curve V that represent boundaries of a
∂φ(X, t) + F |∇φ(X, t)| = 0, (1) ∂t where F is a scalar velocity function depending on the local geometric properties (local curvature) of the front and the external parameters related to the input data e.g, image gradient. The function φ deforms iteratively according to F , and the position of the 2D front is given at each iteration by solving the equation φ(X, t) = 0. The design of the velocity function F plays the major role in the evolutionary process. We have chosen the form F = ν − ²κ, where ν = 1 or −1 for the contracting or expanding front respectively, ² is a smoothing coefficient which is always small with respect to 1, and κ is the local curvature of the front. The latter parameter acts as a regularization term. The term (νg = ± 1) in the PDE specifies the direction of the front propagation either moving inward or outward (g is used to refer to the intensity gray level model associated with a function φg which will evolve according to the above PDE). The problem can be reformulated as classification of each point at the evolving front (points of the narrow band region). If the point belongs to the associated object, the front expands, otherwise (the background) it contracts. The point classification is based on the Bayesian decision [9] at point X. The term (νg ) for each point is replaced by the function νg (X) defined as follows: ½ −1 if πo po (I(X)) ≥ πb pb (I(X)) νg (X) = (2) +1 otherwise where π is the region prior probability and p(.) is the corresponding probability density function (pdf) for the object (o) and the background (b). We characterize each region by a Gaussian distribution with adaptive parameters. Expressions for the prior probability and the Gaussian parameters are found in [10].
4. Alignment of the Training Shapes Let the training set consist of a set of training shapes V1 , ..., VN , with signed distance level set functions defined
as above φ1 , ..., φN . The goal is to calculate the set of pose parameters used to jointly align these shapes, and hence remove any variations in shape due to pose differences. The objective is to find a set of transformations A1 , ..., AN that register an evolving mean shape represented by φM to the given training shapes respectively. We assume that each transformation has scaling components sx , sy , rotation angle θ (with a rotation matrix defined by R), and translations T = [Tx , Ty ]T . The scale matrix will be S = diag{sx , sy }. Global alignment details will come in the following sections.
4.1. Dissimilarity Measure and Energy Function Given two shapes α and β represented by signed distance functions φα and φβ , we need to calculate a transformation (A) defined as above to move the first shape to its target. The signed distance function is invariant to rotation and translation but unfortunately variant to scaling. Now we are going to formulate a dissimilarity measure to overcome this problem. Assume that the vector distance function can be expressed in terms of its projections in the coordinates directions as dα = [dx dy ]T at any point in the domain of the shape α. Applying a global transformation A on the given function φα results in a change of its 0 distance projections as dα = SRdα . This vector is di0 rectly having a magnitude: φα = kSdα k which implies 0 that φα ≤ max(sx , sy )φα . A dissimilarity measure can be directly formulated as r(X) = kSdα (X)k − φβ (A) to measure the difference between the transformed shape and its target representation. A natural energy function can be formulated by summing up the squared difference between the two representations as follows: Z 0 δε (φα , φβ )r2 dΩ (3) E1 = Ω
0
where δε is added to reduce the complexity of the problem as shown in [10] where ε is the width parameter of the band around the shape contour. It is straightforward to show that the given measure r satisfies the relation: r ≤ (sφα (X) − φβ (A)) where s = max(sx , sy ). Hence another energy can be obtained which is always less than or equal to E1 : Z 0 δ (φα , φβ )(sφα (X) − φβ (A))2 dΩ (4) E= Ω
The above energy function has a big advantage since it allows the use of inhomogeneous scales. Also adding the in0 dicator function δ makes the comparison limited to a band around the shapes contours. This is very important to make the perfect alignment with almost zero energy when ε → 0. Then we will obtain identical optimization parameters for both E1 and E. Gradient descent is used to handle the opti-
mization to estimate the transformation parameters. Unfortunately, the function s(sx , sy ) = max(sx , sy ) is not differentiable at the line (sx = sy ). So we use a smeared version of that function based on its original definition: s(sx , sy ) = max(sx , sy ) = sx Hε (sx −sy )+b(1−Hε (sx −sy )) (5) Based on the above formulation, the function returns sx if sx − sy ≥ 0, otherwise it results in sy . The smeared heaviside step function H (defined in [10]) is used to get a smooth transition around the line sx = sy and hence the function is differentiable everywhere. The function ∂s = derivatives will be directly calculated as follows ∂s x ∂s Hε (sx − sy ) + (sx − sy )δε (sx − sy ) and ∂sy = Hε (sy − sx ) + (sy − sx )δε (sy − sx ).
4.2. Mean Shape Estimation Registering the mean shape φM with the shapei represented by φi means si φM ≈ φi (Ai ) where i = 1..N and si = max(sx i , sy i ) as shown above. So the transformations parameters should minimize the following energy function: E(φM , φ1 , ..., φN ) =
N Z X i=1
0
Ω
δ (φM , φi )r2i dΩ
(6)
where ri = si φM − φi (Ai ), is the dissimilarity measure. Rotation and translation do not have any effect on the value of the signed distance map but scaling will result in increasing/decreasing projections in each direction which motivates the use of the smeared max function. The minimization of the energy with respect to the transformations parameters is done through the gradient descent. The mean level set function φM evolves according its calculus of variation using the following PDE: N Z N Z X X 0 ∂ 0 2 ∂ φM = −2 δi si ri dΩ − [ δ r ]dΩ. ∂t ∂φM i i i=1 Ω i=1 Ω (7)
5. The Shape Model Level Set Function The shape model is required to capture the variations in the training set. For this purpose, Each curve/surface is transformed to the domain of the mean function ΦM by its corresponding transformation. The model is considered to be a weighted sum of the transformed signed distance maps deviated from the mean as follows: φp = φM +
N X
wi (φti − φM )
(8)
i=1
where φti = S−1 i φi (Ai ) represents the transformed signed distance map marked by t. We define w = [w1 ...wN ]T
to be the weighting coefficient vector. By varying these weights, Φp can cover all values of the training distance maps and hence the shape model changes according to all the given shapes. In [11], principal component analysis is used to get eigen shapes and a limited number of them is used to build the model. The shape variability is restricted to the selected eigen shapes while in the proposed approach all the training shapes are taken into consideration to enhance the results.
6. Registering the Shape and Intensity Models The shape model is fitted to the image volume by means of registration using a transformation A (defined as before) registering the intensity model φg to the shape model φp . As stated in section 4, the registration is formulated as an energy minimization problem as follows: Z 0 δε (φg , φp )r2 dΩ. (9) E(φg , φp ) = Ω
where r = sφg − φp (A) is the dissimilarity measure. The transformation A and the weight vector w should minimize the objective function E: using the gradient descent flow. A closed form solution for the weight vector is used as shown in [12]. The shape parameters and the transformation give the steady state shape as the segmentation results.
7. Experimental Results A 39 data sets are used to build the kidney shape contours (Fig. 2.(1)a). The alignment and mean shape calculations are carried out to remove the differences between the training instances (Fig. 2.(1).b and (1).c). Image statistics are used to initialize the intensity segmentation level set function φg using the well known stochastic expectation maximization (SEM). Automatic seed initialization is used based on the region parameters estimated by the SEM. A similar approach to that in [10] but the contour seeds are decided based on the Bayesian decision rule. The contour evolves and reaches the steady state to mark the object region boundaries. This results in segmenting some parts of the background as kidney or at the same time missing kidney tissues. This motivates the incorporation of the shape model φp by estimating the weight and the global transformation parameters. Some examples of the shape segmentation results are illustrated in Fig. 2. The approach provides successful results since it can handle different scales, rotations, and translations associated with the registration process. A comparison between the proposed technique and that given in [12] is illustrated in 3. The homogeneous scales do not allow the shape model to correctly deform to mark the correct boundaries of the kidney while this problem is already solved by the proposed methodology.
After segmentation, the ultimate goal of the proposed algorithms is to successfully construct a renogram (mean intensity signal curves) from the DCE-MRI sequences, showing the behavior of the kidney as the contrast agent perfuses into the transplant. In acute rejection patients, the DCEMRI images show a delayed perfusion pattern and a reduced cortical enhancement. Three patients’ renograms are shown in Fig. 3d (1 and 2 are normal but 3 is acute rejection). The normal patient shows the expected abrupt increase to the higher signal intensities and the valley with a small slope. The acute rejection patients show a delay in reaching their peak signal intensities. According to [1], the relative peak signal intensity, time to peak signal intensity, the slope between the peak and the first minimum, and the slope between the peak and the signal measured from the last image in the sequence are the major four features in the renogram of the segmented kidney for classification.
8. Conclusions and Future Work We have presented an efficient approach for the shapebased segmentation problem using level sets. Prior shape information is gathered from real segmented structures for the object of interest. Shapes of the training set are aligned to form an active shape model which is embedded into the image domain by the means of registration. The process incorporate image information to represent a regional shape that is able to be registered with the active shape model. The shape registration process is handled in a new way. The approach allows using more general transformations by incorporating inhomogeneous scales by using a dissimilarity measure approximated by using a smeared version of the maximum function. The technique is used in segmenting the kidney from DCE-MRI images as the basic step for the automatic classification of normal and acute rejection transplants. From the results, the approach is successful when compared with other techniques that use less general transformations. Regarding future work, we are going to use the classification approach in [1] to differentiate between normal and acute rejection patients.
References [1] A. El-Baz, R. Fahmi, S. Esen Yuksel, A. A. Farag, W. Miller, M. Abou El-Ghar, T. Eldiasty. ”A New CAD System for the Evaluation of Kidney Diseases Using DCE-MRI,” Proc. of International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’06), Copenhagen, Denmark, October 1-6, 2006, pp. 446-453. 1, 4 [2] U.S. Department of Health and Human Services. Annual report of the U.S. scientific registry of transplant recipients and the organ procurement and transplantation network: transplant data: 1990-1999. Bureau of Health Resources Department, Richmond, VA; 2000. 1
(1)
(2)
(3)
(4) (a)
(b)
(c)
Figure 2. Training contours are given in 1.a, alignment in 1.b, and average shape is given in 1.c. Segmentation results(2–4): (a) level set initialization, (b) steady state result for the intensity model φg , (c) steady state results for the shape model φp . [3] M. Neimatallah, Q. Dong, S. Schoenberg, K. Cho, and M. Prince, Magnetic resonance imaging in renal transplantation, Journal of Magnetic Resonance Imaging, vol. 10(3), pp. 357368, Sep 1999. 1 [4] D. Yang, Q. Ye, M. Williams, Y. Sun, T. C. C. Hu, D. S. Williams, J. M. F. Moura, and C. Ho., USPIO-Enhanced Dy-
namic MRI: Evaluation of Normal and Transplanted Rat Kidneys, Magnetic Resonance in Medicine, vol. 46, 1152-1163, 2001. 1 [5] G. Gerig, R. Kikinis, W. Kuoni, G.K. van Schulthess, and O. Kubler, Semiautomated ROI analysis in dynamic MRI studies: Part I: image analysis tools for automatic correction of or-
(a)
(b)
70
Normalized Cortex Signal
60
50
40
30
1 2 3
20
10
0
10
20
30
40
50
60
70
Scan #
(c)
(d)
Figure 3. Comparison Segmentation results: (a) level set initialization, (b) results of the approach given in [12], (c) our segmentation results. Plot of different renograms are given in (d).
gan displacements, IEEE Transactions Image Processing, vol. 11:(2),pp. 221-232, 1992. 2 [6] E. Giele, Computer methods for semi-automatic MR renogram determination, Ph.D. dissertation, Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven, 2002. 2 [7] K. Siddiqi, A. Shokoufandeh, S. Dickinson, and S. Zucker. ”Shocks graphs and shape matching,”IEEE International Journal of Computer Vision, 35:1332, 1999. 2 [8] S. Osher and J. Sethian. ”Fronts Propagating with CurvatureDependent Speed: Algorithms Based on the Hamilton-Jacobi Formulation,”Journal of Computational Physics, 79:12–49, 1988. 2 [9] R. Duda, P. Hart, and D. Stork. ”Pattern Classification ,” John Wiley and Sons Inc., 2001. 2 [10] A. A. Farag and Hossam Hassan, “Adaptive Segmentation of Multi-modal 3D Data Using Robust Level Set Techniques“, in Proc International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’04), Saint Malo, France, pp. 143-150, September, 2004. 2, 3, 4 [11] A.Tsai, W. Wells, C. Tempany, E. Grimson, and A. Willsky. ”Mutual information in coupled multi-shape model for medical image segmentation,”In Medical Image Analysis, vol. 8, pp 429–445, December 2004. 4 [12] M. Rousson, N. Paragios and R. Deriche. ”Implicit Active Shape Models for 3D Segmentation in MRI Imaging,” Medical Image Computing and Computer Assisted Intervention
(MICCAI), Part 1, pp 209–216, Saint-Malo, France, September 26-29, 2004. 4, 6
