Knowledge-based Registration & Segmentation of the Left Ventricle: A Level Set Approach Nikos Paragios





Mikael Rousson

Siemens Corporate Research 755 College Road East Princeton, NJ 08540, USA



Visvanathan Ramesh






I.N.R.I.A. BP. 93, 2004 Route des Lucioles 06902 Sophia Antipolis Cedex, France

e-mail: [email protected]

Abstract

segmentation by fitting a boundary template to the image [10] while non-parametric are based on the local propagation of regular curves (in most of the cases using level set methods [12]) under the influence of visual information [9] and shape constraints. Robustness is the main strength of parametric methods while not dealing with important local deformations their limitation. Non-parametric methods can deal with local variations of the structure of interest and complex topologies. Sensitivity to noise and their inability to account for occlusions is the main limitation of such techniques. Registration is also critical component of medical imaging [8]. One can consider registering the shapes of medical structures, their intensity characteristics or combining shape and intensity properties in a joint registration/segmentation framework. In the latest case, registration is performed between the actual segmentation map and either a shape model or a previous segmentation result. In this paper we are interested in detecting and segmenting the endocardium and the epicardium over various stages of the cardiac cycle [13]. The left ventricle (Figure (1.a)) is of particular interest to the physicians because it pumps oxygenated blood to distant tissue in the entire body. To this end, we propose a variational framework that translates high level application constraints to low level segmentation modules, able to deal with the lack of visual information. Methods aim at solving the same problem can be found in [4], while visual components similar to the ones used in this paper were considered in [20]. Constraints driven from the human anatomy are also part of the literature [3, 18] while the use of prior shape models in level set methods was also considered [7] to perform registration and segmentation [2, 17]. Magnetic resonance imaging is a high quality image modality. In this particular application one can claim that visual information is physically corrupted. Papillary mus-

In this paper, we propose a level set formulation to deal with the segmentation and registration of the left ventricle in Magnetic Resonance (MR) images. Our approach is based on the integration of visual information, anatomical constraints and a flexible shape-driven cardiac model. The visual information is expressed through an intensity-based grouping module. The anatomical constraint accounts for the relative positions of the structures of interest. Global shape consistency is introduced by seeking for the lowest potential of the distance between the solution and the prior model. Registration is obtained using the same criterion where the transformation that aligns the latest segmentation map to either the shape model or to the previous segmentation result (temporal domain) is to be recovered.

1 Introduction Segmentation and registration are of great interest in medical imaging. Boundary-based methods for visual grouping are based on the generation of a strength image and the extraction of prominent edges [5]. They can deal with important local deformations and provide accurate segmentation maps. Region-driven techniques [20] explore the assumption of homogeneous intensities properties of the structures of interest with the advantage of being robust to the presence of noise. Precision is a strict requirement during the automatic extraction of structures of interest in the medical field. To this end, curve evolution techniques based on boundary or global regional information were employed to tackle the segmentation problem. These methods can be either parametric or non-parametric. Parametric techniques perform



The authors would like to thank Marie-Piere Jolly and Gareth FunkaLea for fruitful discussions.

1

The snake model [5] can be used to define an objective function that aims at recovering the position of the myocardium by the propagation of two curves towards the desired image properties (strong edges). Lagrangian methods can be used to implement such flows. (a) (b) (c) Figure 1. (a) Input Image, (b) Expected Segmentation Map (manually determined), (c) Part of the Problematic Areas. cles present different intensity properties with respect to the endocardium even if they are considered to be part of it as shown in (Figure (1).b). Furthermore, boundaries between the epicardium and other heart elements cannot be determined reliably using image processing techniques. Consequently, we have to be cautious enough when deriving automatic solutions according to the observed visual information. Our approach is based on three basic principles; the segmentation map has to (i) be supported from the observed visual information, (ii) respect the physiology of the heart, (iii) follow some global shape consistency by component with respect to a generic cardiac model. Towards this end, we propose a level set framework that integrates an intensity-based region component, an anatomical constraint and a term that accounts for global shape-driven consistency. The remainder of this paper is organized as follows. In Section 2 we introduce the application context and the level set representations, while in Section 3 the visual information module is presented. In Section 4, we integrate the anatomical constraint and the shape prior model. Experimental results and discussion appear in Section 5.

2 Application Context The extraction of the myocardium, can be viewed as a three-modal image segmentation problem, seeking the separation between the endocardium, the epicardium and other cardiac components. The structure of interest is the muscle between the endocardium and the epicardium. We consider a common technique in computer vision to recover such solution; the propagation of curves (interfaces);

 Let   be the epicardium interface as shown in (Figure (1).b) (bright contour),

 Let



be the endocardium interface as shown in (Figure (1).b) (dark contour),

A three modal image partition can be defined  using these contours as shown in (Figure (1).b): (i) endocardium,    

  (ii) myocardium, (iii)  background. We recall that the region of interest is the myocardium.

2.1 Level Set Representations Level set representations [11, 12] is an alternative to the Lagrangian approach for tracking moving interfaces. They are implicit, intrinsic, parameter free and can deal with changes of topology. The central idea behind these techniques is to represent the evolving interface a the zero-level set of a structure defined in a higher dimension. In order to formally introduce these representations, we            will consider an evolving (2D) interface          "  # !  . One can assume that the propagation   $&%'! takes place in the direction of the inward normal ( according to a scalar function of the curvature ) ;

*

*

%  2



 

+!-,/.



(   ) +  !'! "!

1 Let 0 be a Lipschitz function that refers to a level set representation: 0

43

8 9 :A@"=CBEDGFIHKJL>NMCD OP >C>A@"=CBED&=CHKJL>NMQDCDGR;S=?>A@"=CBEDGFIJL>NMCD T P >C>A@"=CBED&=CHKJL>NMQDCDGU;S=?>A@"=CBEDGFWV X T JL>NMQDZY 

'5+6'%'!7,

where the evolving interface [ iso-surface

0

] is represented as the zero

[ 

%'!&!\,]   '5_! and ^ distance 43 %'!'! is the minimum Euclidean   ' % ! between the pixel '5_! and the interface .  Then, in order to track using the zero iso-surface of 0 , one can prove [12] that 0 has to evolve using the following PDE: *      * 0 +!`, . ) +!'!a bL0 +!a % '[3

These representations can refer to complex topologies and provide a straight-forward way to estimate the geometrical properties of the evolving interface. A step further is to use these representations within an optimization framework. To this end, one can define the approximations of D IRAC and H EAVISIDE distributions [19]:

cedfZg ƒI„ rg

!…,

!\,ih

;jo = k lmkjRn prqtsCu Owvyxez s {}€ ~ |  =‚k lmkjUn

8 9: u = l†R ;‡o = O l†U o Cz ˆ?‰ s{e~ | p sCu q | O    = / k l$kjU {

nT n

n

These functions can be used to define contour-based as well as region-based energetic modules for the evolving interface

muscles and (iii) the air-filled lungs (Figure (1)). Their intensity properties can be discriminated fairly well. To this end, the observed distribution (histogram) of the epicardium region is a mixture model with three components:

  [Ç

Ž[Ç ZÇ !7,ÉÈ  ! È  ! ȬÊL"Ê ! Ã



à  where È ½È &ÈGʲ6e È È ÈGÊ1,Ë are the a priori probabilities for the endocardium, à à the myocardium and the 

Figure 2. Extraction of the visual properties of the structure of interest and the background. in the level set space [19]:

>AŠZDŒ‹m‹Ž –

q

>[‘…>A@=CBEDCDr’E>A“K>A@"=ZBEDCDr”•@f”•B —&˜ ™ šQ›rœ&?ž'Ÿ~ ¢¡GžQ£&¤e ¥›

q ‹ ‹ …¦ [> ‘…>A@"=CBEDCDr§}>A“K>A@"=ZBEDCDyk ¨©‘…>A@"=#BEDyk ”•@f”•B >AŠ4ŠrD m – —&˜ ™ ª ž'¤Ÿe£'~šQ«$¡¬ž'£&¤ ¥› where ­ and ® are region and boundary positive monotonically decreasing data-driven functions. The first term [i] is a grouping component that accounts for some regional properties (modulo the definition of ­ ) of the area defined by the evolving interface. The second term [ii] is a combination of a boundary attraction term (modulo the definition of ¯ ) and a smoothness component [1, 6]. The extraction of the myocardium involves  for

two level set representations; one the endocardium [ 0 ] and one for the epicardium [ 0 ].

3 Visual Grouping Segmentation can be potentially addressed using boundary information, a powerful cue in image processing. Such technique is a compromise between attraction to the desired image characteristics and internal geometric smoothness constraints [5]. Papillary muscles can provide nonoptimal boundary information (Section 2) within the considered application and therefore an energy component that accounts for smoothness is considered:

* °²±   0 0 !7,³´³fµ – ¦ q [> ‘·¶—&Dy˜ k ¨©‘·¶j™ k – ¦ q >[‘´ÅG—&Dy˜ k ¨©‘´Å\™ k  ¸C¹‡º»½¼ d¾ ºy¿ÁÀ¢ÂÄà Z¸ Æ¿N¼ d¾ ºy¿ÁÀ¢Â This component aims at minimizing the length of the evolving interface and preserves regularity. Regional/global information was considered to improve performance of boundary-based segmentation techniques known to be sensitive to the initial conditions. Regiondriven methods use the evolving interface to define an image partition that is optimal with respect to some grouping criterion. In MR sequences of the left ventricle we can assume the existence of three populations [4]: (i) the blood , (ii) the

[Ç

background hypotheses. The unknown parameters of this model are estimated using the expectation-maximization principle. An example of this analysis is shown in (Figure (2)). Then, according to [14], the image partition that best accounts for the expected image characteristics of the different hypotheses is the one that maximizes the posterior segmentation probability. Such partition can be obtained by minimizing the following cost function:

ÌÎ́ÏAÐmÑrÒCÐ$Ó_ÔÕÖØ×× AÏ Ú ÔCÏ[Û¬Ö q ÏAÚ Ô[Ô log Ï Ü ÏÞÝeÔ[Ô …Ù q Å Ù ¶ Å â ß à#á ¡G«½žQãr~šC£&Á¤¡ Ö©×× ÏNÚ Ô log Ï Ü ÏÞÝeÔ[Ô ÖØ×× Ï[Û$Ö q ÏAÚ Ô[Ô log Ï Ü‡ä´ÏÞÝÔZÔ q  Ù Ù ¶ â ß  ß à#¶ á à#á Å â ª ›#Ÿe£'žQãr~šC£&Á¤¡ ~yã#å&œ&šQž'¤Ÿe£ where the (-log) can be replaced with any positive monotonically decreasing function to avoid stability problems. For the interpretation of this term we recall  the followƒæ„ that  0 !çèj , (ii) ing conditions hold: (i) the endocardium ƒI„   ƒI„ 

 0 … ! é ç _   0 … ! é ê •   the myocardium and the backƒæ„ 

ƒI„   0 !…êë 0 !…êé• . ground 

 will proThe calculus of variations with respect to 0 0 vide adaptive balloon forces [20] that aim at shrinking or expanding the evolving interfaces to provide optimal grouping;

” ‘·¶ì ¦ q [> ‘·¶D log í ¶ >A“‡D k ¨©‘·¶jk ”¢M í Å7>A“‡D q¦ ” Å Å´>A“‡D ‘ ì >[‘ Å D log í ä k ¨©‘ Å k ”¢M í >A“‡D In order to interpret these flows, we consider a given interface point; then the two conflicting hypotheses are evaluated. If the hypothesis corresponding to the interior region is stronger than the one of the exterior, then the interface will expand in the normal direction. In the opposite case it will shrink. The endocardium and the myocardium prob

ability laws are used for 0 , while a comparison between the myocardium and the background density functions de termine the evolution of 0 . Visual-based flows can provide sub-optimal segmentation results. As explained earlier, the intensity properties of the papillary muscles do not follow the endocardium distribution and the distinction between myocardium and other cardiac muscles is problematic. The use of shape-driven constraints can be considered to deal with these limitations.

4 Shape & Physiological Constraints Shape-driven constraints were considered within the propagation of curves in various ways. To this end, one has first to select an appropriate shape representation when introducing such constraints. Moreover, the extraction of an optimal set of parameters able to describe these constraints is to be done given a set of training examples. Deformable models/templates, B-splines, active shape and appearance models, Fourier descriptors are representations that were used to enforce shape consistency. Level set representations is an alternative to these methods with numerous potentials as well as certain limitations. We consider a pixel-wise stochastic level set representation [16] where prior knowledge is represented using a level set function;

"îðï ñ'ò‡ï óeô

rg

õ öj÷mø

!,

î



43

 

'5_!‡ùEúû?üeýEþ•ÿØû
ÿ û

    

defined in the pixel level that consists of two unknown variables;

 The shape image 0 î ,

 The local degrees of variability image ø î . Distance transforms are used as embedding function (Figure (3.c)) in the definition of 0 î . Such prior model consists of a variability image that describes the confidence of the prior model. In areas where important local deformations are plausible high variability estimates are present. This stochastic level set representation is parameter free, implicit, intrinsic and can describe objects with multiple non-connected components. Variational principles according to the maximum likelihood criterion between the model and a training set are usedø to determine the function 0 î and the variability estimates î [16]. This model can be used within the segmentation process to enforce global shape consistency. Let 0 be a level set representation to which we would like to introduce a global rigid-invariant shape constraint according to the model 0 î . We assume that 0 is part of the family of shapes that consists of all possible rigidtransformations of the model. Introducing such constraint can be done by updating locally the evolving representation to meet the model properties; optimal local match. Correspondence is determined through a rigid registration. 0 , we assume the existence Thus, given the current state  of an ideal transformation between the evolving representation and the shape model. In order to better account for the nature of the structure of interest, we assume that the optimal registration corresponds to the maximum likelihood between the representation and the model; 

«

 

í

 

«  "> ! #

>A@"=CBED >A@"=#BED q >A@"=rBŽDCD%$ &Î>A@=#BED('*) +> # >A@"=#BEDCD(,;

where - is the scale factor of the registration model. Level Set Representations with distance transforms as embedding function are invariant to translation and rotation but not to scale variations. One can easily predict that their values are scaled accordingly. Consequently scale appears as a multiplicative factor in the matching process. Solving segmentation/registration now is equivalent with finding a represen tation 0 and a global registration model ;

° 

0Œ



!\, ³7³fµ

ƒI„ /.

0"1243

!

ø î

+5

!'! Ã



   6 }0 ö‡ø 0 î '! !  6  ! î -

7

This functional consists of two unknown variables; (i) a level set representation that is optimal when it becomes a rigid transformation of the prior model, (ii) a transformation (registration) between the evolving current representation and the model. This term is defined in a qualitative manner; model parts with low variability are more significant than the ones that undergo important local deformations. Within the considered application, we can assume two ø shape prior models 8 ; 9 . an endocardium ø  9 .

different ñ ï îLô‡ ñ ï îLô;: and an epicardium ñ ï îLô‡ ñ ï îLô;: one. Clinical expertice can be used to derive these models and introduce shape-drive knowledge-based segmentation components that aim to (i) register the evolving representations to the prior models, (ii) force these representations to respect the prior models:

Ì= Ò+> ÔÕ ¶ Å ¶ Å GBP Ð ÖðÐ K ÏJN Ô O ×× ÏAÐ Ô@BADCFEHGJIK  ÏBN Ô"O
¶ ¶ p ¶  ML ¶ q I J Ï N Ô Q …Ù K ¶ ¶ ML ¶  ¶ RL ¶ ß à#á ›#Ÿ£&žQãZ~šC£½?¤e¡UT V~W›XWyš' ž'š GYP Ð ÖðÐ=K  ÏBN ÏAÐ  Ô @BADCFE G I K  JÏ N Ô O
×× Å Å Q I Kp Å ÏJRN L Ô  Ù q Å Å ML Å  Å ML Å à#á ß › W ãr~yšQ£&?¤eU 6 ¡ TZVy%~ WyX› WšQ ž'š 5

p

Å

ÔJS
â ÔO p

ÔJS â

where is the registration (rigid) transformation between the ø endocardium and its prior shape model 9 .

ñ ï îLô‡ ñ ï îLô : . Segmentation maps are constrained to be rigid transformations of the prior models. An element that has not been taken into account up to now is that within the considered application the evolving interfaces are used to describe the behavior of heart components. Therefore, their positions should respect some given (anatomical) constraints. The relative distances between the endocardium and the epicardium are constrained to be within a given range [3, 13, 18]. To be more specific, one should expect that globally the closest distance between an endocardium point and the epicardium has to be respect the physiology of the

o

A compact representation to reduce the method complexity that combines endocardium and epicardium shapes to a single model can also be considered.

(a)

(b)

(c)

Figure 3. (a) Shape Registration, (b) Level Set Shape Prior Model, (c) Coupling function that is used to introduce the anatomical constraints. heart. This condition is to be met for all pixels of the endocardium. Level set representations using distance transforms as embedding functions provide a natural way to impose this constraint. One can consider a coupling function [ as shown in Figure (3.c). Then, following [3] we can introduce an anatomydriven constraint that aims at preserving the distance between the endocardium and the epicardium within a given range

°]\   0 0 \ ! ,  c „   

c „   

³ ³ µ 7 ³7³ µ 0 %! [ a 0 a !ea bL0 a 0 ! [ a 0 a !ea bL0 a à with the following interpretation; If the distance between the evolving interfaces is within the acceptable range  ^ `_  , then this component becomes inactive. On the other hand, when the constraint is not satisfied, then it tends to propagate the endocardium and the epicardium towards a direction that satisfies the anatomical constraint.

5 The Complete Criterion Integration of visual information, anatomical constraints and the prior shape knowledge can be done towards the automatic segmentation of the myocardium. The optimization criterion consists of four different components (two sets of unknown variables). The evolving level set representations (epicardium, endocardium) and their rigid registration parameters with respect to the prior shape-driven models: a





aHb

>4‘·¶¢='‘´Å´= ¶}= Å$D$ì n – Oed –

agf

[> ‘·—'˜ ¶}=Q‘´ÅG™ D šQ›rœ&?ž'Ÿ

OHh –

>[‘·—'˜ ¶¢='‘´Å¬™ D šQ›rœ&¤e  ~šQ
cN«

a 

[> ‘·—&¶}˜ =&‘´Å¬™ D ~Ÿ~Fcޞ'¡·«

O

¦ –

a
[‘·< ¶¢='‘´—&Å7˜ = ¶}= Å$™ D Vy% ~ Wyj› išQ ž&š

The minimization of the objective function can be done using a gradient descent method and the calculus of variations. Due to the lack of space the corresponding PDEs are not presented in the paper but are available in [15]. As far the energy parameters are concerned, according to our experiments the anatomical constraint holds a minor role. Furthermore, the regularity condition is somehow overwritten by the use of the shape prior term that privileges

Figure 4. Left Ventricle segmentation using Visual information, anatomical constraints and prior shape knowledge. Segmentation is obtained progressively by the propagation of curves using level set methods (presented in a raster scan format). a certain topology. Therefore, we mainly consider two components, the region module and the shape prior term. The performance of the method for cases that do not respect the prior knowledge is questionable. The proposed framework aims to recover a solution that satisfies two independent constraints; (i) respect the prior knowledge (ii) follow the data. Therefore, when the visual information supports strongly the result, the method will be able to deviate from the shape constraint.

6 Conclusions, Experimental Results In this paper we have proposed a variational framework for segmentation and registration of the left ventricle in Magnetic Resonance images. The foundation of our approach is a powerful mathematical tool for evolving interfaces: the level set representations. The registration parameters are estimated reliably due to the information space that is used, the space of signed distance transforms. This model can be either a generic model obtained through a learning approach or some previous segmentation result. Encouraging experimental results were obtained by the application of the proposed framework as shown in (Figures (4,5)). Regarding the registration performance of our approach, a detailed validation on medical shape examples can be found in [15]. The proposed framework has certain limitations. The proper integration of the different segmentation modules ök is an issue. Furthermore, our approach is based on visual information despite the fact that medical information is k acquired in the l space. Therefore, the extension of the method has to be done. Another promising direction is to consider the problem using a single level set representation that corresponds to the myocardium. Last but not least the same framework can be used for segmentation, registration and tracking when replacing the prior model with the seg-

Figure 5. Segmentation and Registration of the Left Ventricle using Visual information, anatomical constraints and prior shape knowledge. mentation map of the previous frame as shown in (Figure (6)).

References [1] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic Active Contours. In IEEE ICCV, pages 694–699, Boston, USA, 1995. [2] Y. Chen, H. Thiruvenkadam, H. Tagare, F. Huang, and D. Wilson. On the Incorporation of Shape Priors int Geometric Active Contours. In IEEE VLSM, pages 145–152, 2001. [3] R. Goldenberg, R. Kimmel, E. Rivlin, and M. Rudzsky. Cortex Segmentation - A Fast Variational Geometric Approach. In IEEE VLSM, pages 127–135, 2001. [4] M.P. Jolly. Combining Edge, Region, and Shape Information to Segment the Left Ventricle in Cardiac MR Images. In MICCAI, pages 482–490, October 2001. [5] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active Contour Models. In IEEE ICCV, pages 261–268, 1987. [6] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi. Gradient flows and geometric active contour models. In IEEE ICCV, pages 810–815, Boston, USA, 1995. [7] M. Leventon, E. Grimson, and O. Faugeras. Statistical Shape Influence in Geodesic Active Controus. In IEEE CVPR, pages I:316–322, 2000. [8] J. Maintz and M. Viergever. A Survey for Medical Image Registration. Medical Image Analysis, 2:1–36, 1998. [9] R. Malladi and J. Sethian. A Real-Time Algorithm for Medical Shape Recovery. In IEEE ICCV, pages 304–310, Bombay, India, 1998. [10] T. McInerney and D. Terzopoulos. Deformable Models in Medical Image Analysis: A Survey. IEEE TMI, 1(2):91–108, 1996. [11] S. Osher and N. Paragios. Geometric Level Set Methods in Imaging, Vision and Graphics. Springer Verlag, 2002. [12] 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.

(a)

(b)

(c)

Figure 6. Cardiac Segmentation/Registration for Magnetic Resonance Images presented in raster-scan format. (a) segmented % frame, (c) segmented %  frame, (d) motion estimation (flow) for the structure à of interest between frame % and frame %  up-scaled four times for demonstration purposes. à [13] N. Paragios. A Variational Approach for the Segmentation of the Left Ventricle in Cardiac Images. In IEEE VLSM, pages 153–160, Vancouver, Canada, 2001. [14] N. Paragios and R. Deriche. Geodesic Active regions for Supervised Texture Segmentation. In IEEE ICCV, pages 926–932, Corfu, Greece, 1999. Previous: INRIA Research Report, RR 3440, June 1998, http://www.inria.fr/RRRT/RR-3440.html. [15] N. Paragios, M. Rousson, and V. Ramesh. Matching Distance Functions: A Shape-to-Area Variational Approach for Global-to-Local Registration. In ECCV, pages II:775–790, Copenhangen, Denmark, 2002. [16] M. Rousson and N. Paragios. Shape Priors for Level Set Representations. In ECCV, pages II:78–93, Copenhangen, Denmark, 2002. [17] A. Yezzi, L. Zollei, and T. Kapur. A Variational Framework for Joint Segmentation and Registration. In IEEE MMBIA, pages xx–yy, 2001. [18] X. Zeng, L. Staib, R. Schukz, and J. Duncan. Volumetric Layer Segmentation Using Coupled Surfaces Propagation. In IEEE CVPR, pages 708–715, Santa Barbara, USA, 1998. [19] H-K. Zhao, T. Chan, B. Merriman, and S. Osher. A variational Level Set Approach to Multiphase Motion. Journal of Computational Physics, 127:179–195, 1996. [20] S. Zhu and A. Yuille. Region Competition: Unifying Snakes, Region Growing, and Bayes/MDL for Multiband Image Segmentation. IEEE PAMI, 18:884–900, 1996.