Region competition based active contour for medical object extraction

Available online at www.sciencedirect.com

Computerized Medical Imaging and Graphics 32 (2008) 109–117

Region competition based active contour for medical object extraction Yanfeng Shang a,b,∗ , Xin Yang a , Lei Zhu a , Rudi Deklerck b , Edgard Nyssen b a

Institute of Image Processing & Pattern Recognition, Shanghai Jiaotong University, Shanghai 200240, PR China b Department of Electronics and Informatics, Vrije Universiteit Brussel, 1050 Brussel, Belgium Received 1 March 2006; received in revised form 6 September 2007; accepted 15 October 2007

Abstract In this paper, a probabilistic and level set model for three-dimensional medical object extraction is proposed, which is called region competition based active contour. The algorithms are derived by minimizing a region based probabilistic energy function and implemented in a level set framework. An additional speed-controlling term makes the active contour quickly convergent to the actual contour on strong edges, whereas a probabilistic model makes the active contour performing well for weak edges. Prior knowledge about the initial contour and the probabilistic distribution contributes to more efficient extraction. The developed model has been applied to a variety of medical images, from CTA and MRA of the coronary to rotationally scanned and real-time three-dimensional echocardiography images of the mitral valve. As the results show, the algorithm is fast, convergent, adapted to a broad range of medical objects and produces satisfactory results. © 2007 Elsevier Ltd. All rights reserved. Keywords: Medical image; Object extraction; Active contour; Region competition; Level set; Snake model; Prior knowledge

1. Introduction In clinical applications, automatic medical object extraction is often a necessary preprocessing step for three-dimensional (3D) image reconstruction, quantitative analysis and computeraided diagnosis. Due to the wide variety of shapes, the complexity of the topology, the presence of noise in a complicated background, and the diversity of imaging techniques, extracting medical objects automatically and accurately is a very difficult task. Active contours are among the most successful image segmentation techniques used in a variety of applications, such as computer vision, pattern recognition, and biomedical image processing. Depending on the involved image features, the active contour based extraction method can be categorized as: edge based [1–4], region based [5–10], higher order tensor based [11] and other higher level prior knowledge based methods [12–15]. Edge based object extraction is historically the first and therefore also the most widespread instance of the active contour model. This model derives external forces from the magnitude of the gradient of the image, which drives the active contour

∗

Corresponding author. Tel.: +86 21 34201324. E-mail address: [email protected] (Y. Shang).

0895-6111/$ – see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2007.10.004

towards the object boundaries characterized by maximal values of the magnitude of the image gradient. Examples of such models are balloon models [3,4], geodesic active contours [1], geometric active contours [2] and other extensions of these models. A good review can be found in [16,17]. The disadvantage of these models is in general that the active contour may leak out of the ideal contour when the edges are weak. Region based object extracting models show more advantages than edge based models on weak edges. The motion of the active contour depends – in addition to the gradient – on region and pixel features whereas the speed is in general determined by a statistical distance metric expressing the difference between local pixel, object and background features. A number of applications show the validity of these models [5–10]. However, a common disadvantage of these models is that the active contour will move to and fro with a high speed on strong edges because of the large deviation of the features between neighboring pixels. This may lead to imprecise segmentation results. An attractive approach is the region competition model [6], which uses generalized Bayes and minimum description length (MDL) criteria to make contours converge to a local optimum within a Snake model. One drawback is that it cannot deal well in a natural way with topological changes: additional region growing algorithms have to be involved to merge the

110

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

contours. In addition, since it is hard to split an active contour, it requires to initialize a large number of contours and to evolve them in turn which increases the computing complexity. Another appealing region based model is the active contours without edges approach [5]. It uses a single feature, i.e. image intensity, to control the movement of the active contour within a level set framework. This model has however a few disadvantages. There may be some imprecision on weak edges in medical images and the active contour will oscillate about strong edges. The aforementioned models may be further complemented with higher order gray level derivatives and higher level prior knowledge for enabling the extraction of complex medical objects. Higher order gray level tensors offer more information on local geometric properties than image intensity and gradient alone, e.g. the Hessian matrix, a second-order tensor, can be used for vascular object extraction [11]. Prior knowledge about the shape of the specific object can enhance the extracting process in terms of robustness and efficiency [12–15], but it is at the cost of generality. For example, the a prior shape of the left ventricle may be a powerful constraint for tracking its contour [18], but it may fail when the cardiac anatomy deviates strongly from the normal heart. Motivated by the above problems, we propose a novel object extracting model, called region competition based active contour model, the RCAC. This model combines edge and region features in a level set evolution equation and has been extended to 3D to extract medical objects. Since it is based on a level set algorithm, it can deal with topological changes in a natural way. A damping term reduces oscillations when the active contour reaches a sharp edge, whereas a prior distribution of the object gray level intensities makes the extraction process more robust in a cluttered background. The remainder of this paper is organized as follows. In Section 2, concepts and mathematics related to edge based and region based active contour models are presented. In Section 3, the RCAC model is formulated in detail and an RCAC based 3D medical object extraction approach is presented. Medical applications and experiments are discussed in Section 4. Finally, a conclusion is formulated in Section 5. 2. Background The edge based active contour finds an optimal object boundary by minimizing an energy function over all closed curves [1–3]. The energy function is defined as a curve integral:  1   2 C (q) dq E(C) = α 0



+β 0

1

 C (q)2 dq − γ



1

most used methods to find the optimal solution. Geodesic active contour (GAC) [1] minimizes the energy function by steepest descent, which leads to a level set formula. The evolution equation is: ∂φ = g(|∇I|)(k + v0 )|∇φ| + ∇g(|∇I|) · ∇φ ∂t

where k is the Euclidean curvature, v0 is a constant which makes the active contour expanding (v0 > 0) or shrinking (v0 < 0), and g: [0,∞] → + is a strictly decreasing function, such as g(|I|) = (1 + |I|)−1 . The region based active contour extracts the objects by minimizing energy function over regions [6,7], others than curves. The energy function is always defined as a region integral of some statistical features over all objects. In a statistical and variational approach, region competition [6] uses a generalized Bayes/MDL criterion, which is also a generalized Mumford–Shah function [7], E[Γi , {αi }] ⎧ ⎫ ⎪ ⎪    k ⎨ ⎬  μ = ds − λ log P{I(x, y) |αi }dxdy + γ ⎪ 2 Γi ⎪ ⎭ i=1 ⎩ Ωi

(3) where Γi and αi are the boundary and the statistical distribution of region Ωi (i = 1, . . ., k), respectively. The first term in the right hand side is the length of the boundary curve for the region Ωi . The second term is the sum of the cost for coding the intensity of every pixel (x,y) inside region Ωi according to a distribution P{I(x,y) : (x,y)∈Ri |αi }. γ is assumed common to all regions. Region competition minimizes above criteria function in a Snake model framework. The model initializes a large number of Snakes in different homogeneous regions, and evolves these Snakes in turn. The neighboring Snakes will compete between neighboring regions. When they are homogeneous, they would be merged by a region-merging algorithm which is based on squared Fisher distance. Active contour without edge [5], which is called C-V model, after the initials of its authors, is also a region based active contour. The C-V model minimizes a Mumford–Shah energy function in a level set framework. The curve evolution equation is: ∂φ = δε (φ){μ · div(∇φ/∇φ) − v − λ1 (I(x, y) − uo )2 ∂t +λ2 (I(x, y) − ub )2 }

|∇I(C(q))| dq

(2)

(4)

(1)

0

where I: [0,a] × [0,b] → + is a given target image, C: [0,1] → 2 is a parameterized closed curve, and α, β and γ are real positive constants. The Snake model and level set are the

where uo and ub are the mean intensity of object and background, respectively. δε (φ) is a one-dimensional Dirac function, v is a constant expanding or shrinking term. μ, λ1 and λ2 are constant parameters.

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

3. Model and system description 3.1. Region competition based active contour model For an n-dimensional image, the objects extracting problem can be considered as a classification of object and background. The energy function is defined as:  E[Γ, {αo , αb }] = −λ{ log P(I(x)|αo )dx  +

Ωo

Ωb

Γ

ds

Ωo = {(x)|φ(x) > 0} = H(φ(x)) Backgound : Ωb = {(x)|φ(x) < 0} = 1 − H(φ(x))



H(φ) =

(6)

Γ = {(x)|φ(x) = 0} = δ(φ(x)) 0 φ(x) ≤ 0 1 φ(x) > 0

δ(φ(x)) = H  (φ(x))

and δ(φ) is a Dirac function which is zero on Ωo and Ωb and nonzero on Γ .

ds is the integral of the object’s surface. The surface Γ integral of a function f(x) over the boundary Γ is defined as

f (x)δ(φ(x))|∇φ(x)|dx [19]. Then, the third item in formula Ω (5) could be written as:  δ(φ(x))|∇φ(x)|dx (7) Ω

Substituting (6) and (7) to (5), yields,  E[φ(x), ∇φ(x), {αo , αb }] = δ(φ(x))|∇φ(x)|dx −λ{  +

Ω

Ω

log P(I(x)|αo )H(φ(x))dx

log P(I(x)αb )(1 − H(φ(x)))dx}

(9)

The energy E will be minimized when the Euler-Lagrange differential equation is satisfied [20]: ∂F ∂ ∂F ∂ ∂F ∂ ∂F − − =0 − ∂φ ∂x ∂φx ∂y ∂φy ∂z ∂φz

(8)

(10)

This equation has been written for the three dimensional case, but it can be readily extended to n-dimensions. The first term in Eq. (10) is ∂F dδ(φ) = −λ[log P(I|αo ) − log P(I|αb )] · δ(φ) + μ |∇φ| ∂φ dφ (11) The remaining terms in Eq. (10) are ∂ ∂F ∂ ∂F ∂ ∂F + + ∂x ∂φx ∂y ∂φy ∂z ∂φz   ∂F ∂F ∂F = div , , ∂φx ∂φy ∂φz   ∂|∇φ| ∂|∇φ| ∂|∇φ| = div μδ(φ) , μδ(φ) , μδ(φ) ∂φx ∂φy ∂φz   φx φy φz = div μδ(φ) , μδ(φ) , μδ(φ) |∇φ| |∇φ| |∇φ|   ∇φ =μ div δ(φ) |∇φ|   ∇φ ∇φ = μδ(φ) div +μ∇δ(φ) |∇φ| |∇φ|   ∇φ dδ(φ) = μδ(φ) div +μ |∇φ| |∇φ| dφ

(12)

It then follows that    ∇φ δ(φ) −λ[log P(I|αo ) − log P(I|αb )] − μdiv |∇φ| =0

Ω



− λlog P(I|αb )(1 − H(φ)) + μδ(φ)|∇φ|}dx  = F (x, φ, ∇φ)dx

(5)

Object :

where

Ω

Ω

where Ωo and Ωb are the object and the background regions, respectively, Ωo ∪Ωb = Ω. αo and αb are statistical distribution of object and background, respectively. Γ is the closed boundary of the objects x = [x1 , x2 , . . . ,xn ]T . The two regions Ωo and Ωb will compete on the boundary Γ according to region distribution α to minimize function E, and then the active contour evolves according to the local distribution, the local smoothness, and the global distribution of the object and background. The region competition procedure is implemented in a level set framework which can deal with region splitting and merging naturally. The level set embeds the contour curves or surfaces as a zero level set in a higher dimensional function φ(x), the boundary Γ is obtained for φ(x) = 0. φ(x) is always defined as a signed distance function (SDF) [19], |φ(x)| = 1. The object and background region can be formulated by a Heaviside function H(φ):

Contour :

For a fixed image, the distribution α of object and background depends on the contour Γ . Γ is determined by φ(x). Therefore, the energy E is a function of φ(x) and φ(x). It can be rewritten as:  {−λlog P(I|αo )H(φ) E[φ, ∇φ] =



log P(I(x)αb )dx} + μ

111

(13)

In order to find the solution to the static equation, a standard approach consists of making it dynamic, so that an initial level set function φ, different from the solution to the static equation, can gradually evolve towards it. Note that when φ has converged (i.e. no longer changes in time), ∂φ/∂t = 0 and the static equation is

112

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

indeed satisfied    ∇φ ∂φ = δ(φ) λ[log P(I|αo ) − log P(I|αb )] + μdiv ∂t |∇φ| (14)

The coefficient δ(φ) limits the equation to be active at the object contours which leads to a narrow band in the numeric implementation. A higher evolving speed can make the active contour run out of the narrow band which results in an unexpected error. Replacing δ(φ) by |φ|, enables the zero level set to move in a broader region. The intensity of image I(x) varies strongly on object edges, which leads to substantial fluctuations of [log P(I|αb ) − log P(I|αo )] for slight positional changes of the boundaries near the edges. Hence, this will cause the zero level set oscillating rather intensively on sharp edges and leading to unacceptable segmentation errors. To achieve convergence at image edges, it proved to be necessary to add a speed-controlling term 1/(1 + |I × φ|). Only when the magnitude of I and φ are high and they have approximately the same direction, this speed-controlling term will be quite small, and will greatly slow down the progression of the active contour, causing it to become more convergent on sharp edges. Preprocessing with a narrow Gaussian filter will further improve the steadiness of the speed-controlling term near the edges. Hence, a new evolution equation is formulated as ∂φ |∇(φ)| = ∂t 1 + |∇G(I) · ∇φ|    ∇φ × λ[log P(I|αo ) − log P(I|αb )] + μdiv |∇φ| (15) The active contour evolves as a result of competition between object and background regions. It expands when P(I|αo ) > P(I|αb ) and shrinks when P(I|αo ) < P(I|αb ). So we call it region competition based active contour (RCAC). Formula (15) is a generalized iterative equation for object extraction. We can embed object and background distributions or any other classifiers in it. Supposing that the object and the background are Gaussian distributed, where the distribution parameters are αo ∼ (μo ,σ o ), αb ∼ (μb ,σ b ). μ is the mean intensity and σ is the standard deviation. It follows: |∇(φ)| ∂φ = ∂t 1 + |∇G(I) · ∇φ|      (I−μo )2 (I−μb )2 ∇φ × −λ +μdiv +v − 2σo2 |∇φ| 2σb2 (16) where v is a parameter related to σ o and σ b , and it is set to a constant value in our tests. When parameters (μ, σ) are unknown

Fig. 1. Flow chart of medical object extraction based on RCAC.

in advance, they can be estimated at each iterating step. That is   I(x)H(φ)dx (I(x) − μo )2 H(φ)dx 2 Ω Ω  μo = , σo = H(φ)dx H(φ)dx Ω  Ω I(x)(1 − H(φ))dx μb = Ω , (17) (1 − H(φ))dx  Ω (I(x) − μb )2 (1 − H(φ))dx σb2 = Ω  (1 − H(φ))dx Ω

3.2. System descriptions We developed a medical object extracting system based on RCAC. Fig. 1 shows the flow chart of medical object extraction. There are three inputs for the RCAC model: original images, initial zero level set, and prior distribution. The RCAC model yields object contours as output, which can be used as an initial zero level set for the next slice or volume, in case there is sufficient similarity between neighboring images in the same volume or neighboring volumes of a time sequence. A region of interest (ROI) can also be obtained from it in order to calculate the prior distribution of object and background. In this object extracting system, the original images, which are pre-processed, could be 2D or 3D, CT, MR or Ultrasound. The initial zero level set could consist of several mouse clicks, then a large number of iterations would be required to make them evolve to the whole organ. It could also be a prior organ mask and a smaller number of iterations would be necessary to obtain a satisfactory result. Surely, the latter will be more efficient. Taking into account a prior Gaussian probability structure, the following equation is obtained ∂φ |∇(φ)| = ∂t 1 + |∇G(I) · ∇φ|



×



−λ

(I − μo )2 2σ  2o

−

(I − μb )2 2σ  2b

 + μdiv

 ∇φ  |∇φ|

 +v

(18)

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

where (μo , σo ) and (μb , σb ) are the statistical distribution parameters of object and background. When an exact prior distribution is known, μ = μprior and σ  = σ prior are taken. When the prior distribution is only roughly known, μ = (μ + μprior )/2 and σ  2 = 2 )/2 are taken. (σ 2 + σprior In general the background in medical images can be quite complex and difficult to be modeled as one Gaussian distribution. However, it can be treated as one Gaussian distribution in a small ribbon/ROI around the object. Therefore, we first try to obtain an approximate contour of the object, which is then further expanded into the background by using a simplified RCAC model     ∂φ |∇(φ)| ∇φ = μdiv +v (19) ∂t 1 + |∇G(I) · ∇φ| |∇φ| where v is a positive constant which makes the active contour expand. In order to make the estimators more robust, a morphological erosion is performed both on the approximate object region and on the background ribbon. Preprocessing is important because noise is inherent in some medical images. A smoothing filter could improve the accuracy of classification and speed up the evolution of the active contour in noise stained images. The modified curvature diffusion equation (MCDE) [21,22] is a good choice, since it preserves the edge and reduces the noise. The equation is given as ut (x) = c(|∇u(x)|) u(x) + ∇c(|∇u(x)|)∇u(x) where u0 is the original image and c(|∇u|) = e K a constant parameter.

(20) −(|∇u|/K)2

, with

4. Experiments and evaluations The designed model has been applied to two problems: the coronary artery extraction in CTA and MRA and the mitral valve extraction from 3D ultrasound data. In the former case, the performance of the RCAC scheme has been compared to that of the GAC and C-V algorithms. In the second case, the results of automatic mitral valve segmentation are compared with ones achieved manually.

by isolating the terms characterising the region properties in (4) and (16): −λ[(I − μo )2 /(2σo2 ) − (I − μb )2 /(2σb2 )] and −λ1 (I − μo )2 + λ2 (I − μb )2 for the RCAC and the C-V models, respectively. It can be noted that the C-V algorithm assumes equal variances. Given λ1 = λ2 = λ, the model behaves like a linear discriminant classifier with as decision boundary (μo + μb )/2 under the supplementary assumption that the prior probabilities for the object (o) and the background (b) class are equal: i.e. P(o) = P(b). On the other hand, the RCAC based classifier takes into account the variances of the classes, and behaves therefore like a quadratic discriminant classifier. Under the assumption of equal prior probabilities P(o) = P(b), it will separate the object and the background at a decision boundary corresponding to the intersections of the two Gaussian intensity probability distributions, which is more optimal than the decision boundary obtained with the linear discriminant classifier related to the C-V-model (see Fig. 2). Results achieved by the C-V and RCAC algorithms on the extraction of the coronary artery in one slice of a 3D CT study, are presented in Fig. 3. As the arrows indicate, the C-V model loses some details due to the presence of weak edges and this results in a narrower object as σ o > σ b . To facilitate the comparison of the accuracy between the two models, we adapted the C-V model by introducing RCAC’s convergence speed coefficient: 1/(1 + |∇G(I) · ∇φ|). When tracking the coronary artery, the previously mentioned ROI strategy is adopted to calculate the statistical distribution of object and background, since we are dealing, in this case, with a complicated background. Results of coronary artery extraction by the GAC algorithm are shown in Fig. 4. The GAC model is more sensitive to the initial contour than the RCAC model. The initial contour in Fig. 3(a) leads to segmentation error when applying the GAC algorithm. Two more accurate initial contours have to be chosen which are shown in Fig. 4(a). The green one leads to (b) which is slightly larger than the actual object and which is inferior compared to RCAC’s result in Fig. 3(c). We set the number of iterations to 300 and v = −2, since the initial contour should shrink. The red initial contours lead to (c) which is unacceptable. Here, we set

4.1. Coronary artery extraction in CTA and MRA Coronary artery disease is considered as the most frequent of all diseases that affect the heart. The coronary arteries get narrowed or blocked due to deposition of fat, cholesterol and other substances. A precise extraction is a necessary step to provide the information allowing for an aided diagnosis. As the coronary artery is only a few pixels wide, a slight segmentation error may lead to misdiagnosis of stenosis and hence a robust extracting model is mandatory. The RCAC model offers a solution for these medical objects. A comparison of coronary artery extraction is made between RCAC, C-V and GAC models as follows. Two Bayesian two-class classifiers – derived from a Normal distribution model and taking positive values for pixels belonging to the object and negative values for the background pixels – are derived from the RCAC and C-V models,

113

Fig. 2. Results of classification of RCAC and C-V model.

114

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

Fig. 3. Coronary artery extracted in a CT slice by C-V and RCAC models. (a) Original CT image and initial contour, (b) coronary artery extracted by the C-V model and (c) coronary artery extracted by the RCAC model.

Fig. 4. Coronary artery extracted in a CT slice by the GAC model. (a) Initial contours for the GAC model, (b) coronary artery extracted by the GAC model from the green initial contour outside the object and (c) coronary artery extracted by the GAC model from the red initial contours inside the object.

the number of iterations to 300 and v = 1, since we want the initial contours to expand. The RCAC, C-V and GAC algorithms have been compared in terms of their convergence and their impact on the results. Fig. 5(b) presents the change of the object area over subsequent iterations. The plot shows that the C-V model’s output keeps on oscillating around the target area, while the progress towards convergence of the RCAC algorithm is steady and the

final border is reached after only 9 iterations. In fact, a high value of log P(I|αo ) − log P(I|αb ) in Eq. (15) on strong edges will lead to a high evolving speed and may drive the active contour to pass over the actual edge. Yet, the speed coefficient 1/(1 + |G(I)·φ|) will act as a counterbalance and will slow down the active contour when it approaches the edge of the target object, so that an optimal contour extraction can be achieved. As to the GAC model, it reaches the actual area after 25 iterations,

Fig. 5. Comparison of convergence on edge between RCAC, C-V and GAC models. (a) Coronary extracted in 2D MRA image, (b) time–area curve of RCAC, C-V and GAC of (a).

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

Fig. 6. Coronary reconstructed with the C-V (a) and RCAC (b) models.

but then it further shrinks (in the case of the initial green contour) so that the object area continuously decreases. The GAC model is always problematic in the presence of weak edges. A visual comparison (Fig. 6) shows that the results of the extraction while using the RCAC model manifests higher smoothness, when both algorithms are applied in a slice-by-slice approach. The bumpiness of the surface present in the C-V output is introduced by the oscillations around the object edges as described before. We have not reconstructed the coronary artery

115

by the GAC model because of the aforementioned shrinking problem, experienced in certain slices. In addition, a 3D implementation of the RCAC model has been developed and applied on the 3D CTA and MRA data sets in consecutive chunks consisting of 15 slices (Fig. 7). After convergence, the result obtained for the current chunk is used as the initial surface for the next chunk of 15 slices. For each chunk, the RCAC-algorithm will estimate the statistical parameters of the object and background intensities, since they are quite varying over the full volume. The arrow on the CTA case indicates a case of stenosis, and shows how easily a stenosis can be identified, when a quite accurate extraction is provided. The resolutions of two 3D images are 512 × 512 × 45 and 512 × 512 × 80, respectively. We set the number of iterations to 45, the time step in a way that the active contour maximally moves 2.5 pixels, the smoothing coefficient μ to 20 and the expanding or shrinking constant ν to zero. The total calculation time is 24 and 44 s on a 3.4 GHz Pentium 4 computer, respectively. 4.2. Mitral valve extraction The mitral valve disorder is any defect in the valve which lies between the left atrium and the left ventricle. It is one of the

Fig. 7. Coronary reconstructed using the 3D implementation of the RCAC model from (a) CTA and (b) MRA.

Fig. 8. Mitral valve reconstructed by the RCAC model from images obtained with the HP Sonos 5500 ultrasound scanner. (a) Original ultrasound image; (b) mitral valve extracted by the RCAC model and (c) surface reconstruction of the extracted mitral valves.

116

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

Fig. 9. Mitral valve reconstructed by the RCAC model from a real-time 3D ultrasound volume. (a) A slice of the original ultrasound volume; (b) volume rendering of the mitral valve; (c) the mitral valve extracted by the RCAC model and (d) surface reconstruction of the extracted mitral valve.

most common heart-valve disorders among people of all ages. The primary diagnostic method of valvular defect is to observe its shape and movement by echocardiography images. Thus, the efficient extraction of the mitral valve can provide an important added value to the diagnosis. We applied the RCAC model to an application involving mitral valve extraction. For the purpose of our experiments, two medical databases have been acquired, from a rotational ultrasound scanner HP Sonos 5500 with TTO probe and a real-time 3D ultrasound scanner Philips Sonos 7500. Examples of the scan data are depicted in Figs. 8 and 9. In the first case, the data has been acquired at 3◦ angular slice spacing over a span of 180◦ at a frame rate of 24 frames per second. Thus, one patient scan consists of a series of 2D images obtained at different angles around a pivot point. Indeed, a volume time series is generated where each volume represents the heart at different time instances during one cardiac cycle. Typically, one patient scan consists of approximately 1000 images. The resolution of an image is 240 × 256. The patient database consists of 12 scans (8 male and 4 female patients) in the age range between 6 months and 14 years. For the extraction of the mitral valve by the RCAC model, we set the number of iterations to 100, the time step in a way that the active contour maximally moves 2.5 pixels, the smoothing coefficient μ to 15 and the expanding or shrinking constant ν to zero. The total calculation time is about 10 s for each volume on a 3.4 GHz Pentium 4 PC. In the second case, a single data set consists of 9–18 volumes captured within one cardiac cycle. The resolution of each volume is 208 × 144 × 160. The database consists of 40 data sets, from which 50% are pathological cases (17 mitral regurgitation, 1 mitral stenosis and 2 mitral valve prolapse). In this application, we set the number of iterations to 80, the time step in a way that the active contour maximally moves 2.5 pixels, μ and ν again to

Table 1 The mean, maximum and standard deviation of the distance between RCAC results and the manual delineations of the 2D rotational echocardiography data (1 mm = 2.43 pixels)

Patient 1 Patient 2

Mean distance (mm)

Maximum distance (mm)

Standard deviation (mm)

0.7012 0.6188

2.7673 1.9646

0.5107 0.4715

respectively 15 and zero. The total calculating time is about 8 s for each volume on a 3.4 GHz Pentium 4 PC. To give a quantitative and objective evaluation of the results of our method, a validation has been performed by comparing the RCAC segmentation results to the results of the manual delineation performed by two experienced cardiologists. Due to the large amount of manual work, only a limited number of patient data has been manually outlined. The mean, maximum and standard deviation of the distance between the two segmentation results are calculated. The measurements (presented in Tables 1 and 2) together with the visual assessment of the results demonstrate that the RCAC model based extraction yields accurate object delineations with regard to the manual outlines. Table 2 The mean, maximum and standard deviation of the distance between RCAC results and the manual delineations of the real-time 3D echocardiography data (1 mm = 1.68 pixels)

Patient 1 Patient 2 Patient 3

Mean distance (mm)

Maximum distance (mm)

Standard deviation (mm)

0.7148 0.8339 0.6552

1.7643 2.6685 1.5029

0.5199 0.6040 0.4726

Y. Shang et al. / Computerized Medical Imaging and Graphics 32 (2008) 109–117

5. Conclusion In this paper, we presented the Region Competition based Active Contour model, which is a novel 3D object extracting approach based on a probabilistic level set model. It expresses the minimization of a region based probabilistic energy function, implemented in a level set framework. An additional speed constraint enables the active contour to extract the objects accurately at an improved speed and without oscillation at the object boundaries. We successfully applied the model to the medical applications of coronary artery and mitral valve extraction from CTA/MRA and ultrasound images, respectively. The experiments show that the algorithm is fast, convergent, adaptable to a broad range of medical objects and producing satisfactory results, in comparison to the C-V model of [5], the Geodesic active contours of [1] and manual segmentations. Acknowledgements This work was partially supported by the National Natural Science Foundation of China (No. 60572154) and the National Basic Research Program of China (973 Project) (No. 2003CB716104). References [1] Caselles V, Kimmel R, Shapiro G. Geodesic active contours. Int J Comput Vis 1997;22:61–79. [2] Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propagation. IEEE Trans Pattern Anal Mach Intell 1995;17(2):158–75. [3] Cohen L. On active contour models and balloons. Computer vision, graphics, and image processing. Image Understanding 1991;53(2):211–8. [4] Cohen L, Cohen I. Finite-element methods for active contour models and balloons for 2-d and 3-d images. IEEE Trans Pattern Anal Mach Intell 1993;15(11):1131–46. [5] Chan TF, Vese LA. Active contours without edges. IEEE Trans Image Process 2001;10(2):266–77. [6] Zhu SC, Yuille A. Region competition: unifying snakes, region grouping, and Bayes/MDL for multiband image segmentation. IEEE Trans Pattern Anal Mach Intell 1996;18(9):884–900. [7] Mumford D, Shah J. Optimal approximation by piece-wise smooth functions and associated variational problems. Communication Pure Appl Math 1989;42(5):577–685. [8] Frederic P, Michel B, Thierry B, et al. Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm. IEEE Trans Image Process 2005;14(7):910–24. [9] Mukherjee DP, Ray N, Acton ST. Level set analysis for leukocyte detection and tracking. IEEE Trans Image Process 2004;13(4):562–72. [10] Freedman D, Zhang T. Active contours for tracking distributions. IEEE Trans Image Process 2004;13(4):518–26. [11] Lorigo L, Faugeras O, Grimson WEL, et al. Curves: curve evolution for vessel segmentation. Med Image Anal 2001;5:195–206. [12] Daniel C, Florian T, Joachim W, et al. Diffusion snakes: introducing statistical shape knowledge into the Mumford–Shah function. Int J Comput Vis 2002;50:295–313.

117

[13] Michael EL, Grimson WEL, Olivier F. Statistical shape influence in geodesic active contours. Comput Vis Image Und 2000;1:316–23. [14] Chen Y, Hemant D, Tagare ST, et al. Using prior shapes in geometric active contours in a variational framework. Int J Comput Vis 2002;50:315–28. [15] Shang YF, Yang X, Zhu M, et al. Prior based cardiac valve segmentation in echocardiographic sequences: geodesic active contour guided by region and shape prior. IbPRIA 2005. Proceedings, Part II (LNCS Vol. 3523): 447–54. [16] Jasjit SS, Kecheng L, Sameer S, et al. Shape recovery algorithms using level sets in 2-D/3-D medical imagery: a state-of-the-art review. IEEE Trans Inf Technol Biomed 2002;6(1):8–28. [17] McInerney T, Terzopoulos D. Deformable models in medical image analysis: a survey. Med Image Anal 1996;1(2):91– 108. [18] Nikos P. A level set approach for shape-driven segmentation and tracking of the left ventricle. IEEE Trans Med Imaging 2003;22(6):773–6. [19] Sethian James A. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge, UK: Cambridge University Press; 1999. [20] Arfken G. Mathematical methods for physicists. 3rd ed. Orlando, FL: Academic Press; 1985. [21] Pietro P, Jalhandra M. Scale-space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell 1990;12:629–39. [22] Whitaker R, Xue X. Variable-conductance, level-set curvature for image denoising. ICIP 2001:142–5. Yanfeng Shang received the MS degree in water conservancy and electric power engineering in 2003 from North China University of Water Conservancy and Electric Power. He has been a PhD student in the Institute of Image Processing and Pattern recognition at Shanghai Jiaotong University since 2003. His research interests include medical image analysis and pattern recognition. Xin Yang holds the MS degree in control engineering from Northwestern Polytechnic University, Xi’an, China in 1982 and the doctor of applied science degree in electronic engineering from the Vrije Universiteit Brussel in 1995. Since 1997, he has been a professor in the Institute of Image Processing and Pattern Recognition, Shanghai Jiao Tong University, Shanghai, China. His current research activities are in the area of medical image analysis and partial differential equations in image procession. Prof. Yang is the author of more than 70 papers in journal and refereed conference proceedings. Lei Zhu received the MS degree in software engineering in 2004 from Shanghai Jiaotong University. He has been a PhD student in the Institute of Image Processing and Pattern recognition at Shanghai Jiaotong University from 2004. His research interests include medical image visualization and analysis. Rudi Deklerck received the PhD degree in applied sciences in 1999 from Vrije Universiteit Brussel (VUB), Belgium. He is a professor at the Department of Electronics and Informatics (ETRO), VUB and is active within the Interdisciplinary institute for BroadBand Technology (IBBT) since April 2004. His research is mainly situated in the domain of medical image processing and covers a number of topics: i.e. 3D surface extraction, visualization of medical images, medical image analysis and computer aided detection. Edgard Nyssen is a professor at the Applied Sciences Faculty of the Vrije Universiteit Brussel. His teaching topics include digital circuit design, biomedical statistics, pattern recognition, and program development with C, C++, and Java. His research interests include the application of pattern recognition techniques and statistical modelling to image analysis, as well as basic research on pattern classifiers – more specifically, evaluation techniques and new pattern classifier models.