Bi-planar image segmentation based on variational geometrical active contours with shape priors

Medical Image Analysis 17 (2013) 165–181

Contents lists available at SciVerse ScienceDirect

Medical Image Analysis journal homepage: www.elsevier.com/locate/media

Bi-planar image segmentation based on variational geometrical active contours with shape priors El Hadji S. Diop ⇑, Valérie Burdin Image and Information Department, LaTIM-INSERM U 650, Telecom Bretagne, Technopôle Brest Iroise, CS 83818, 29238 Brest Cedex 3, France

a r t i c l e

i n f o

Article history: Received 12 January 2012 Received in revised form 18 September 2012 Accepted 28 September 2012 Available online 26 October 2012 Keywords: Image segmentation Variational formulation Euler–Lagrange equations Partial differential equations Level sets

a b s t r a c t This work proposes an image segmentation model based on active contours. For a better handling of regions where anatomical structures are poorly contrasted and/or missing, we propose to incorporate a priori shape information in a variational formulation. Based on a level set approach, the proposed functional is composed of four terms. The first one makes the level set keep the important signed distance function property, which is necessary to guarantee the good level set evolution. Doing so results in avoiding the classical re-initialization process, contrary to most existing works where a partial differential equation is used instead. The second energy term contains the a priori information about admissible shapes of the target object, the latter being integrated in the level set evolution. An energy that drives rapidly the level set towards objects of interest is defined in the third term. A last term is defined on prior shapes thanks to a complete and modified Mumford–Shah model. The segmentation model is derived by solving the Euler–Lagrange equations associated to the functional minimization. Efficiency and robustness of our segmentation model are validated on synthetic images, digitally reconstructed images, and real image radiographs. Quantitative evaluations of segmentation results are also provided, which also show the importance of prior shapes in the context of image segmentation. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Image segmentation is one of the most important tasks in the areas of image processing and computer vision. That is one of the reasons why new segmentation techniques are always being developed, and more elaborated models are still carried out. In medical imaging, despite it constitutes a crucial part, in most cases it is only a first step towards a real target application, which could be a study, an analysis or a therapeutic evaluation of segmented medical parts and structures. One of the most known and easiest segmentation methods is the image thresholding. It is known that in the case of Computed Tomography (CT) images, that, most of the times, such an approach could be satisfactory to segment anatomical structures. Unfortunately, for most of other image acquisition techniques, an image thresholding and other low level segmentation methods are not efficient1 enough, and most of times, at all, in order to get the objects and/or regions of interests. That were the major reasons of why ⇑ Corresponding author. Address: Center of Mathematical Morphology, Department of Mathematics and Systems, MINES ParisTech, 35 rue Saint-Honoré, 77305 Fontainebleau Cedex, France. E-mail addresses: [email protected] (E.H.S. Diop), valerie. [email protected] (V. Burdin). 1 Here, ‘‘efficiency of the segmentation’’ refers to the segmented objects of interests themselves and the accuracy of the segmentations. 1361-8415/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.media.2012.09.006

much more elaborated and much more sophisticated techniques have been developed since then; particularly, the active contours or snakes methods Kass et al. (1987). Local contour deformations are obtained by minimizing a defined functional (or energy) such that its minimum is reached on the edges of the target object (or objects when dealing with many objects of interests) one would like to segment. Depending on the way to formulate and implement the active contours, different approaches were proposed. Indeed, the active contours methods can basically be grouped in two major classes: the parametric active contours (PACs) and the geometric active contours (GACs). PAC is defined in a Lagrangian domain, and explicitly formulated with a parametrized curve to be deformed afterwards. On the other hand, GAC are defined in an Eulerian domain, and represented implicitly with level sets of a bidimensional function. In addition to better robustness regarding the initialization of the active contour, GAC really brought noticeable improvements, in comparison to the PAC. Furthermore, besides more efficient and more robust numerical schemes, due mainly to the level set formulation, topological changes were allowed during GAC evolution. However, whatever the formulation and implementation that are used, one should keep in mind that (P/G) AC strongly rely on image gradients and other local features. And then, it appears clearly that if the image one would like to segment has low gradients, the active contours will not become stationary at object boundaries, and such approaches will be definitely inappropriate and inefficient to get good segmentations.

166

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

Unfortunately, what we just described is truly the case in most medical images. Indeed, maybe except for CT images, the majority of medical imaging modalities; particularly, X-rays modality for which the bi-planar EOS system and images we are dealing with are also issued, presents generally many anatomical structures having sensibly the same gray level values (e.g., it is the case for bones). Also, images have weak contrast, and very often, many parts of anatomical structures boundaries are not salient or do not appear at all; this is also amplified by occlusions phenomena in bones joints. Furthermore, due to their complicated forms in many cases (the femur and tibia bones for instance), another issue regarding anatomical structures is about the fact that it is very difficult to describe or model them just by using parametric models. If that was possible, then, one could use Bardinet et al. (1996) approach where superquadrics were proposed, followed by some refinement formulated as an inverse problem. Therefore, it becomes clear that using only the local image features is definitely insignificant for characterizing an object (or objects) of interests in this kind of images. One known way to resolve this problem and obtain better accuracy in the segmentation of the desired structures, is to incorporate a priori information about the boundary of the target object shape. 1.1. Motivation Images we are currently dealing with are obtained with the EOS acquisition system which provides very low dose X-rays radiographs, simultaneously in the frontal and profile positions, along with 3D bones reconstructions gained only from the already obtained 2D frontal and profile radiographs (see Laporte et al., 2003). However, in other positions different from the standing one, the 3D reconstruction of bones is impossible, and even if one can obtain a 3D reconstruction, and in an extremely difficult way, the latter one is very far away from the exact 3D bone. A 2D/3D image registration method, as proposed in Zosso et al. (2008) and Jerbi et al. (2009), is a way to overcome that issue. Indeed, the registration provides a way to find the positions of the structures reconstructed in the bi-planar radiographs of other positions, from the standing position. However, that approach works if and only if the only structures present in the bi-planar radiographs, are the structures of interests; this is in fact not the case at all for real radiographs. The algorithm proposed in Jerbi et al. (2009) was successfully applied on digitally reconstructed radiographs (DRRs) images. The case of real radiographs is more complicated though, and results are not actually satisfactory. This is mainly due to the presence of various structures in real radiographs, in addition to problems related to gray level image intensities between structures in the image, shadings, gaps in object edges, contours saliency, occlusions phenomena in bones joints, etc. Segmentation results will be used to make the 2D/3D registration more robust and more accurate, in order to perform better analyses of the hip and knee joint pathologies. Therefore, the accuracy and computation time are essential factors to be taken into account. 1.2. Main contributions We propose a new segmentation model formulated in a variational framework with a functional composed of four terms. For each of the energy terms, the explanation of its potential advantages on the active contour evolution is provided. The first term allows the level set to keep during its evolution the essential property of SDF. This way, instead of solving periodically a PDE as it is done classically, this first term eliminates the re-initialization procedure. Shape prior information is defined in the second term and incorporated in the level set evolution. We define in the third term an energy that drives quickly the level set

towards the boundaries of the target object. The last term is a complete and modified Mumford–Shah functional. Integrating all the four terms together allows us to provide a theoretical explanation of the good behaviors of the segmentation method, which is shown to be robust in the presence of occlusions and missing parts. 1.3. Paper organization The outline of the paper is as follows. In Section 2, we present a brief review of some shape prior-based image segmentation methods. Our proposed model based on geometric active contours is detailed in Section 3. Qualitative and quantitative results are presented in Section 4, for both synthetic, DRRs and real radiographs images. Concluding remarks and perspectives of this work are proposed in Section 5. 2. Brief review on shape prior-based segmentation In this section, we propose a review of some segmentation approaches based on active contours and a priori shapes information. The active contours are formulated either in a parametric way as for PAC or in geometric form as for GAC. In either case, prior shapes are incorporated in the active contour evolution model, or are performed in an independent stage of the segmentation process. We first quote the work of Leventon et al. (2000) who proposed to integrate a set of deformable shapes during the level set evolution. Prior shapes were obtained by defining a probability distribution on variances of the training set elements. The proposed shape model was based on a Principal Components Analysis (PCA), which was then applied on SDFs built with the training set composed of target objects contours. Then, at each step of the proposed segmentation algorithm, the level set evolved locally thanks to the intrinsic image features, such as the gradients and curvature, but the evolution was also guided by the estimations of the maximum a posteriori of the shape prior position and shape. However, as pointed out in Bresson et al. (2006), due to the fact that a posteriori probability is maximized at each iteration of the algorithm through an independent optimization process, then, the final evolution equation is no more a PDE. In fact, two independent steps are mandatory to evolve the level set. Proposed evolution equation in Leventon et al. (2000) is defined as follow:

uðt þ 1Þ ¼ uðtÞ þ k1 ðgðc þ jÞjruðtÞj þ ruðtÞ  prgÞ þ k2 ðuH ðtÞ  uðtÞÞ: ð1Þ Chen et al. (2002) proposed a different method to integrate the shape prior information in the active contours evolution. In fact, contrary to the work of Leventon et al. (2000) where the shape prior was obtained in a probabilistic manner, Chen et al. proposed to construct it as the mean of the training contours, which were primarily rigidly aligned. Also, the proposed evolution equations are truly PDEs, and authors also proposed a theoretical proof of the existence of the minimum of their functional; the proof was based on bounded variations functions. The following functional was proposed:

EðC; l; R; TÞ ¼

Z 0

1



 k 2 gðjrIjðCðpÞÞÞ þ d ðlRCðpÞ þ TÞ jC 0 ðpÞjdp: 2 ð2Þ

Following the same principle, Thiruvenkadam et al. (2008) proposed to resolve the segmentation issue in the case where occlusions appeared specifically in edges. The problem was then formulated as a segmentation with depth; the objectives being to find objects contours, object intensities and spatial order. The following functional was proposed:

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

  12 N C B Z B p N N C X X X C B EðU; C obj ; C int Þ ¼ cp Hð/p Þ  ð1Þp1 cp;k Sp;k  ~c e S b C dx BI 0  C B X p¼1 p¼2 k¼1 A @ 0

þk

Z X Z X N N b Sð/k Þdx: jrHð/k Þjdx þ b X k¼1

X k¼1

ð3Þ

We end this section with the Bayesian approach combined with prior shape information, which was proposed by Chang et al. (2009). With their method, it is in fact possible to include several prior shapes and then, choose the most suitable one. A major benefit of that work is the feasibility of the segmentation of more than one object of interest in the image. It is also important to notify that with their proposed segmentation algorithm, prior shapes are constructed interactively by hand. We do not wish to give more practical details, but let us just indicate that the Bayesian optimization was mainly performed with the graph cuts algorithm proposed by Boykov and Kolmogorov (2004). Next, we introduce our proposed segmentation model. 3. Proposed variational image segmentation model Our segmentation algorithm is presented in details in this section. In fact, a variational approach of the segmentation problem is proposed and formulated with level sets combined with prior shapes of contours of the ought-to-be segmented object. As discussed earlier, adding a priori shape information in the segmentation formulation, guides the active towards expected zones, i.e. target objects edges. Indeed, the construction of prior shapes is crucial for real radiograph images, mainly because of the kind of anatomical structures depicted in medical images, the almost same gray level image intensities between objects of interests (femur or tibia for instance) and other structures (soft tissues, pelvis, femur, tibia, fibula, tarsus, metatarsus, etc.), problems of edges saliences, occlusions phenomena in bones joints, etc. In this section, we first start with the prior shapes construction. Our proposed energy functional is presented in Section 3.3, while the variational formulation is introduced in Section 3.4. 3.1. Prior shapes designing Let us first recall that many methods have been proposed in Leventon et al. (2000), Cremers et al. (2002), Rousson and Paragios (2002), Chen et al. (2002), Chan and Zhu (2003), Thiruvenkadam et al. (2008), and Chang et al. (2009), to construct a priori objects shapes. In all approaches, the main goal was to statistically take the maximum of advantages on the possible object deformations regarding the training dataset, in order to constrain, during the segmentation process, the level set to evolve only on the basis of the predefined set of admissible (deformable or not) shapes. Let us just make a short review of some existing approaches. Based on a set of a priori contours, Cremers et al. (2002) proposed a nonlinear approach to construct the prior shapes. Their method relied on an estimation of the probability density of the explicit data provided by the training contours, and was based on a Gaussian kernel assumption modeling, and also on the Mercer kernels (see Courant and Hilbert (1953) for more details). Authors proposed a sort of extension of the PCA in a probabilistic domain. A different approach was proposed by Rousson and Paragios (2002) who chose to include the local deformations of the shapes, i.e. their variances, built from the training contours data. In fact, using a variational approach and under the assumptions that each grid of the prior shape could be modeled with a Gaussian probability density, authors pro-

167

posed to estimate the ‘‘best’’ shape by looking for the maximum likelihood of local densities of the training elements with respect to the shape image and variances of the shapes deformations. The idea followed by authors is in fact the same as the one proposed by Wang and Staib (1998), even if the application context was different. Indeed, the prior shape modeling was designed for a non rigid image registration that was based on an elastic model combined with statistical prior shapes. In the work proposed by Chen et al. (2002), authors provided a simpler and more intuitive approach to design a priori shape, which was then incorporated in the level set evolution. Indeed, the shape prior was constructed statistically as being equal to the mean of contours that were comprising the training dataset. Let us notify that all the training contours were in advance rigidly aligned to one of them. Finally, the shape prior was also built in Chan and Zhu (2003), Thiruvenkadam et al. (2008) and Chang et al. (2009) as a single shape. But, in addition, in order to take into account the possible shape deformations, other degrees of freedom were allowed in the shape model, on one hand, by the means of rigid transformations as in Chan and Zhu (2003), in the other hand, with the use of static and dynamic labeling functions to make active regions where the shape prior should precisely take place in Cremers et al. (2003). In Chang et al. (2009), a priori shape was selected interactively by hand, while its effectiveness was incorporated dynamically in the segmentation algorithm, and in a more subtle way, by making the information provided by the spatial order estimation to enforce the shape prior constraints only on occluded edges in Thiruvenkadam et al. (2008). For many reasons, particularly because of the robustness and suitability to be incorporated in our energy functional presented in Section 3.3, we choose in this work to follow Leventon et al. (2000) proposed approach to construct the set of admissible and deformable prior shapes. Let us start by considering N aligned contours having the same shape with the target object, and let us denote them by {Cn}16n6N. For each element n of the dataset, let us define the contour Cn in an implicit way, i.e. by embedding it as a signed distance function (SDF) denoted by /n (see e.g. Osher and Fedkiw, 2003; Osher and Paragios, 2003 for more details about SDFs).  be the mean of Let {/n}16n6N be our training dataset, and let /  is given by /  ¼ 1 PN / . the SDF. / n n¼1 N To learn the local deformations of our training dataset, we follow the main idea proposed in Leventon et al. (2000), and compute the variances of the /n by performing a PCA on {/n}16n6N, instead of doing it on the contours {Cn}16n6N. Although such an approach is different and presents undoubtedly many advantages, the philosophy behind this way of learning the shape deformations is quite similar to the works of Turk and Pentland (1991) and Wang and Staib (1998). For instance, Turk and Pentland (1991) used the so-introduced notions of eigenfaces consisted in applying the PCA on each face image in order to build a recognition model. With our approach, SFDs are defined implicitly, and then, there is no need to handle curve parametrization terms. Furthermore, as also pointed out in Leventon et al. (2000), this way of proceeding increases in fact the robustness regarding a slightly misalignment of the dataset contours, which yields to an inaccurate objects matching, which induces then a wrong learning of the shapes deformations. The main reason of that relies on the fact that when one uses the parametrized contours, then, a small misalignment compares the image intensity (could be gray level values) variance between independent objects in Turk and Pentland (1991), while correlations in SDFs of insignificant misaligned pixels are very strong. Another benefit in considering SDFs is the fact that it avoids the tedious problem of pixel-by-pixel contours matching, and the bunch in errors prone resulting from such a registration approach. In a theoretical point of view, when applied to {/n}16n6N, the PCA method looks for the best orthonormal basis, we denote by

168

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

{uk}16k6K with of course K < N, such that the projection of /n on that basis has a minimal distance in the sens of the L2 norm.  for Let M be the matrix containing terms of the form uk  /, k = 1,2, . . . , N. Using a single value decomposition, the covariance matrix given by N1 MMt could be decomposed as:

V RV t ¼

1 MM t ; N

ð4Þ

where the columns of the matrix V hold the orthogonal modes of variations of the target object shape, and R represents a diagonal matrix whose elements are the singular values. Let U be the matrix composed of the p first columns of V. Then, newly built shapes, denoted by /, can be obtained in fact by performing a projection on the space of eigenvectors. Thus, one has:

 þ Uk; /¼/

ð5Þ

where U = [ui]16i6p designates the matrix of eigenvectors, and k ¼ ½ki t16i6p represents the vector of weights or shape parameters. Optimal values of k will be determined by solving the PDE described in Section 3.4, and obtained by minimizing our proposed segmentation functional. Thus, depending on the values of ki, a priori shape  will be deformed from the mean SDF /. Next, different shape parameters are used for examples to show some newly created shapes that are obtained thanks to Eq. (5). 3.2. Shape deformation examples In this study, we mainly consider three datatypes. Synthetic images are first considered. Afterwards, we deal with digitally reconstructed radiographs (DRRs), and end with real radiograph images. The first example copes with synthetic images. Our training dataset is composed of N = 30 shapes embedded as SDFs. The PCA yields two principal modes with a fitting percentage of 98%. New created shapes can then be obtained by using principal modes and the following values ki = ±3, i = 1, 2, in Eq. (5). Results are displayed in Fig. 1. A last real radiograph example is provided where new shapes are obtained by using a dataset of N = 21 right aligned contours of femurs. Performing PCA on embedded corresponding SDFs gives five principal modes of variations, with a fitting probability to f/n gNn¼1 of about 95%. Fig. 2 illustrates created shapes based on Eq. (5) and for ki = ±1.3, i = 1, 2. 3.3. Energy functional for the segmentation model This section is dedicated to the proposed functional from which the segmentation model presented in Section 3.4 will be derived from. For many reasons, previously discussed, we choose a level set formulation in our segmentation energy functional. Let I be the image, and let us consider the open and bounded set X  R2 be the image domain. We assume that the boundary of X denoted by oX is Lipschitz. Based on local image features, and using the set of admissible prior shapes information whose construction was preliminary detailed in Section 3.1, as in Diop et al. (2010), we define our segmentation functional as:

Fðu; k; V; Iin ; Iout Þ ¼

a

F 1 ðuÞ þ F 2 ðu; k; VÞ þ bF 3 ðuÞ 2 þ mF 4 ðk; V; Iin ; Iout Þ;

ð6Þ

where a, b and m are positive parameters to counterbalance the different energies by giving weights to energy functionals F1, F2, F3 and F4, respectively defined as follows:

Z

ðjruj  1Þ2 dx; Z h i c ng þ /2 ðk; hðxÞÞ dðuÞjrujdx; F 2 ðu; k; VÞ ¼ 2 X Z gHðuÞdx; F 3 ðuÞ ¼

F 1 ðuÞ ¼

ð7Þ

X

ð8Þ ð9Þ

X

and last, we have:

F 4 ðk; V; Iin ; Iout Þ ¼

Z

jI  Iin j2 þ ljrIin j2 Hð/ðk; hðxÞÞÞdx þ

Z

X

jI

X 2

2

 Iout j þ ljrIout j ð1  Hð/ðk; hðxÞÞÞÞdx Z þ f dH1 ðCðk; hðxÞÞÞ:

ð10Þ

X

In the notations we used, u represents the level set function; n and c are positive parameters; d() is the unidimensional Dirac distribution; g is an edge detector function; Iin and Iout designate the Mumford–Shah functional terms; m is a positive parameter. H() is the Heaviside function, that is: H : R ! f0; 1g, s.t. H = 1 in Rþ ; H ¼ 0 0 in R H , and d = H in the sense of distributions. Also, herein, / represents the shape prior; k is the vector of shape parameters; h is a function that accounts the spatial rigid transformations achieved through a vector V which holds the scale parameter factor, the rotation matrix and the translation vector; C(k, h(x)) = {x 2 X, s.t. /(k, h(x)) = 0}; and last, H1 is the one dimensional Hausdorff measure. Next, we come back on each term of the overall proposed functional in order to give both the expected contributions and the interpretation of each energy in the segmentation process. 3.3.1. Functional energy term F1 Let us first recall that level sets are implicitly defined in a higher dimension space than contours, and are embedded as SDFs. In their classical formulation (see e.g., Osher and Sethian, 1988; Malladi et al., 1995; Caselles et al., 1997), it is very important to make them keep, at least around the zero level set neighborhood, the SDF property, so that they are prevented from being either too flat or very sharp, which induces then, in either case, a wrong evolution of the level set. To avoid that major level set drawback, the most common way is to apply the re-initialization procedure, which consists of solving periodically after some iterations on the level set function, as in Sussman et al. (1994), Zhao et al. (1996), and Chan and Vese (2001), the following PDE:

8 < @u ¼ signðu0 Þð1  jrujÞ; @t : uðx; 0Þ ¼ u0 ðxÞ;

ð11Þ

where u0 refers to the level set to be re-initialized. As stated by Chan and Vese (2001), this step could be interpreted as a rescaling and regularization process. However, as mentioned in Osher and Fedkiw (2003), the re-initialization process should not be too much important, because doing so keeps interior contours from growing. Based on variants of the PDE (11), different methods for maintaining the SDF property were proposed (see Peng et al., 1999; Sussman and Fatemi, 1999; Osher and Fedkiw, 2003). It is important to notice that the efficiency of all the proposed re-initialization techniques is only guaranteed if and only if the function to be re-initialized is close to a SDF, prior to the application of such a procedure. In fact, the re-initialization process definitely suffers against a major side effect, which is the fundamental difference between the theory and the practice of level set methods. Indeed, as pointed out in Gomes and Faugeras (2000), the solution of the PDE (11) is not a SDF, and consequently, this makes arise many questionable points for which the most crucial ones could be addressed as at what step it should be applied and also how to do it efficiently. Actually, there

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

169

(a)

Fig. 1. Examples of newly built shapes obtained with the shapes deformations Eq. (5) and shape parameters ki = ±3, i = 1, 2. Aligned contours of the training basis are  is plotted in red. Black curves from second line represent the zero levels of resulted / functions. (For interpretation of the references to displayed in (a), and the mean SDF / color in this figure legend, the reader is referred to the web version of this article.)

are yet no accurate responses about that. Let us also mention the fact that most of the time, the re-initialization process is used as for refinements of segmentation results, by tuning the number of its applications and steps in order to drive the active contour in the expected target object edges. Gomes and Faugeras (2000) proposed a different way to overcome the re-initialization step, and considered the more general Hamilton–Jacobi equation:

8 @u > > < @t ¼ bjruj; jruj ¼ 1; > > : uðx; 0Þ ¼ u0 ðxÞ;

ð12Þ

F 1 ðuÞ ¼

ðjruj  1Þ2 dx:

F 2 ðu; k; VÞ ¼

Z h

ð13Þ

ð14Þ

X

This functional avoids the re-initialization process and guarantees the SDF property. In fact, when keeping jruj bounded, one ensures

i c ng þ /2 ðk; hðxÞÞ dðuÞjrujdx; 2

ð15Þ

where n, c > 0; g is an edge detector, which is a positive and strictly non increasing function. In this work, it is taken as:

g¼

In this study, we also account that constraint in a different manner, by integrating it in the overall segmentation energy. Our approach is similar to Li et al. (2005) method, which was proposed in a context of non shape prior based segmentation. The functional F1 is defined as follow:

Z

3.3.2. Functional energy term F2 The functional F2 is defined as follow:

X

where authors did impose in fact the constraint jruj = 1. Indeed, it is known that a SDF u satisfies jruj = 1; but conversely, any implicit function u satisfying jruj = 1 is the SDF to a surface plus a constant (the constant equals to 0 on the surface) (see Arnold, 1983). Thus, by taking into account that constraint, authors showed that the resulted PDE, which is in fact not an Hamilton–Jacobi equation, contrary to Eq. (12), should write:

8 < @u ¼ bðx  uruÞ; @t : uðx; 0Þ ¼ u0 ðxÞ:

the correct computations of u derivatives without any need to reinitialize the function u.

1 1 þ gjrGr HIj2

;

ð16Þ

where w stands for the convolution operator, g > 0, Gr is a Gaussian kernel with a variance r. To take into account the pose of the target object in the image, a rigid spatial transformation function is also incorporated in the segmentation model, by introducing the function h defined as follow:

h :X ! X; x#hðxÞ ¼ sRh x þ T;

ð17Þ

where s > 0 represents the scale parameter, h is the rotation angle and T = [Tx, Ty] designates the translation vector. The rotation matrix Rh is thus given by:

 Rh ¼

 cos h sin h :  sin h cos h

ð18Þ

We denote by V ¼ ½s; h; T 2 XV  Rþ  ½0; 2p  R2þ , the vector holding the rigid transformation parameters, where XV refers to the domain of V. The energy F2 makes evolve the active contour towards high gradients areas and also, most of all, towards regions that are similar to the prior shape, until the stationary case is reached at the expected edges of the target object. Indeed, the minimization

170

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

(a)

Fig. 2. Examples of newly built shapes obtained with the shapes deformations Eq. (5) and shape parameters ki = ±1.3, i = 1, 5. Aligned contours of the training basis are  is plotted in red. Black curves in second and third lines represent the zero levels of resulted / functions. (For interpretation of the displayed in (a), and the mean SDF / references to color in this figure legend, the reader is referred to the web version of this article.)

of F2 increases the similarity between the active contour and the shape prior. In the expression of the functional F2,/2(k, h(x)) is in fact a similarity measure between a point x 2 X of the active contour (i.e. zero level set of u), and newly constructed shapes gained from Eq. (5) and thanks to the shape parameters k and the eigenvectors U. Let us draw the remark that a comparable energy was proposed by in Chen et al. (2002) (cf. Eq. (2)), with major differences. Indeed, the segmentation energy to be minimized was only that one, and furthermore, the shape prior was constructed, as recalled in Section 2, as a unique shape taken as the statistical mean of the training contours comprising the dataset, which were aligned beforehand, on the basis of a rigid transformation. 3.3.3. Functional energy term F3 The energy F3 is defined as:

F 3 ðuÞ ¼

Z

gHðuÞdx;

3.3.4. Functional energy term F4 The last functional F4 is defined as a modified version of the complete Mumford–Shah functional energy (see Mumford and Shah, 1989). Such a model is the most suitable one for our bones (femur, tibia) segmentation problem, in the sense that those anatomic objects are represented in bi-planar images as regions having almost homogeneous gray level image intensities. Let us denote K(in,out) = jI  I(in,out)j2 + ljr I(in,out)j2. The energy term F4 is then defined as follow:

F 4 ðk; V; Iin ; Iout Þ ¼

Z

Kin Hð/ðk; hðxÞÞÞ þ Kout ð1 Z  Hð/ðk; hðxÞÞÞÞdx þ f dH1 ðCðk; hðxÞÞÞ; X

ð20Þ

X

ð19Þ

X

where H() is the Heaviside function. This functional could be seen as a weighted area term of the target object of interests. In fact, in the case where g is taken as equal to 1,F3 is then the area of the region {x 2 X, s.t. u(x) < 0}, which is simply the area of the object of interests. Thus, the minimization of F3 yields another force that will allow the active contour to move quickly towards object edges.

where C(k, h(x)) = {x 2 X, s.t. /(k, h(x)) = 0}, and H1 represents the one dimensional Hausdorff measure. In fact, the first term of that energy was proposed in Bresson et al. (2006) in the context of ventricles segmentations in Magnetic Resonance Imaging. The Hausdorff measure term was not added, and authors justified their choice in doing so by the fact that newly built shapes by the means of the PCA were smooth enough. Our study context is different though. Indeed, bones structures, for instance femora, tibias, etc., have more complications than ventricles. Therefore, it is very important to add another constraint in order to have more accurate

171

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

segmentations, which is in fact resolved here by constraining the length of the object in the second term of F4. We show in Section 4 dedicated to the numerical results, the usefulness of the such added term in the segmentations. In the next part, we present the variational formulation of our segmentation problem, and related evolution equations that constitute our proposed segmentation model.

8 @I H out > < @t ¼ I  Iout þ lMIout in Xout  Rþ ; @Iout ¼ 0 on ð@ Xout [ fx 2 X s:t: /ðk; hðxÞÞ ¼ 0gÞ  RH þ; > : @nout Iout ðx; 0Þ ¼ IðxÞ; x 2 Xout ;

ð25Þ

where nout is the outward pointing unit normal vector of the boundary of Xout denoted by oXout which is assumed to be Liptschitz.

3.4. Variational formulation and PDE-based segmentation model In preceding sections, we detailed the construction of the set of admissible deformable shapes (cf. Section 3.1), and presented the proposed segmentation energy functional (cf. Section 3.3). Based on that, we present now the variational formulation of the segmentation problem, along with the segmentation model described through the resolutions of related PDEs, which will be obtained by solving the Euler–Lagrange equations associated with the minimization of the proposed functional. The level set formulation of the variational problem is thus given by:

u; k; V; Iin ; Iout ¼ arg minu;k;V;Iin ;Iout Fðu; k; V; Iin ; Iout Þ a ¼ arg minu;k;V;Iin ;Iout F 1 ðuÞ þ F 2 ðu; k; VÞ 2 þ bF 3 ðuÞ þ mF 4 ðk; V; Iin ; Iout Þ;

Mumford–Shah functional defined in the energy F4 (cf. Section 3.3.4). Euler–Lagrange equations associated to Iin and Iout in the minimization of Eq. (21), are thus given by:

ð21Þ

where F is in fact given by Eq. (6), and F1, F2, F3 and F4 are given respectively by Eqs. (7)–(10), as presented in details in Section 3.3. By the means of calculus of variations (see e.g., Evans (1998) and Aubert and Kornprobst (2002) for more information about that notion), we solve the Euler–Lagrange equations associated to the variational formulation (Eq. (21)) to finally obtain the following evolution equations:

8    h i @u u u > ¼ a Mu  div jr þ dðuÞdiv ðng þ c/ðk; hðxÞÞÞ jr > > @t r u j r u j > > < bgdðuÞ; in X  RH þ; > @u H > > @n ¼ 0 on @ X  Rþ ; > > : uðx; 0Þ ¼ u0 ðxÞ; x 2 X; ð22Þ where n is the outward pointing unit normal vector of the boundary of X denoted by oX which is assumed to be Liptschitz.

R R 8 @k ¼ 2c X ðrk /ðk; hðxÞÞÞ/ðk; hðxÞÞdðuÞjrujdx  f X rk /ðk;hðxÞÞ > @t > > 0 < d ð/ðk;hðxÞÞÞdx; R >  m X ðKin  Kout Þrk /ðk;hðxÞÞdð/ðk;hðxÞÞÞdx; in Xk  RHþ ; > > : kðx; 0Þ ¼ k0 ; x 2 Xk ;

8 @Iin H > < @t ¼ I  Iin þ lMIin ; in Xin  Rþ ; @Iin ¼ 0 on ð@ Xin [ fx 2 X t:q: /ðk; hðxÞÞ ¼ 0gÞ  RH þ; > : @nin Iin ðx; 0Þ ¼ IðxÞ; x 2 Xin ;

where nin is the outward pointing unit normal vector of the boundary of Xin denoted by oXin which is assumed to be Liptschitz. Our segmentation model is now constituted by Eqs. (22)–(25). 4. Numerical results In this section, we show the segmentation results obtained on synthetic images (Section 4.1), on digitally reconstructed radiographs (DRRs) (Section 4.2), and finally, on real image radiographs (Section 4.3). Our segmentation model is basically composed of Eqs. (22)–(25), which are implemented in that order, i.e. first, Eq. (25); next, Eq. (24); afterwards, Eq. (23); and finally, Eq. (22). Segmentation results obtained with our model are compared to ones gained from both a non shape prior-based segmentation approach as it was proposed in Li et al. (2005), and another segmentation method incorporating a priori shapes information as in Bresson et al. (2006). Furthermore, quantitative results are also provided in Section 4.4 for a segmentation evaluation purpose, in order to complete obtained qualitative results, which will be presented in the next three subsections. From now on, we will refer to the non shape prior-based segmentation as classical active contours method. Some numerical aspects should be clarified before going further. In fact, in the numerical implementations, the Heaviside H() and Dirac d() functions are approximated respectively by smoothed C2 and C1 functions, as it was the case in Sussman et al. (1994) and Zhao et al. (1996). Approximated functions used for H() and d() are then respectively given by:

8 > < 1; H ðxÞ ¼ 0; > :1

2

(

ð23Þ where Xk is the domain of shape parameters which are used to build new shapes as in (5).

R 8 @V ¼ 2c X ðrV hðxÞÞðr/ðk; hðxÞÞÞ/ðk; hðxÞÞdðuÞjrujdx; > > > @t R > > f X ðrV hðxÞÞðr/ðk; hðxÞÞÞd0 ð/ðk; hðxÞÞÞ þ hrV ðr/ðk; hðxÞÞÞ; > > < E r/ðk;hðxÞÞ dð/Þdx; j r /ðk;hðxÞÞj > > R > > > m X ðKin  Kout Þrk /ðk; hðxÞÞdð/ðk; hðxÞÞÞdx; in XV  RH > þ > : Vðx; 0Þ ¼ V 0 ; x 2 XV ; ð24Þ where XV represents the domain of V ¼ ½s; h; T 2 XV  Rþ  ½0; 2p  R2þ , which is the vector holding the rigid transformation parameters (17) with a scale s, a rotation matrix Rh (18) of angle h and a translation vector T = [Tx, Ty]. Let us denote Xin = {x 2 X, s.t. /(k, h(x)) > 0} and Xout = {x 2 X, s.t. /(k, h(x)) < 0}. The terms Iin and Iout are related to the

ð26Þ

d ðxÞ ¼

0; 1 2

if x > ; if x < ;   1 þ x þ p1 sin px ; if jxj 6 ;



if jxj > ;   1 þ cos px ; if jxj 6 :

In our experiments, the parameter

ð27Þ

ð28Þ

 is set as  = 1.5.

4.1. Synthetic images We first start by illustrating the proposed segmentation method on synthetic images. The training basis we use and a priori shapes were already built and designed previously in Section 3.1, and were displayed in Fig. 1. The robustness of the approach is evaluated on its capabilities to segment the object of interests, and even in the case where occlusions appear in edges or some missing parts appear in the image. Our approach is then primarily illustrated on an occluded image, for which occlusions appear in the nearby of image edges. Occlusions have first the same gray level intensities as the true object (cf. Fig. 3a). Fig. 3 shows noticeable benefits that could be

172

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

Fig. 3. (a) Original synthetic image. (b) Test image: artificial occlusions created on image contours. (c) Active contour initialization. (d) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (e) Our approach. (f) Proposed method in Bresson et al. (2006). (g) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

gained from shape priors-based segmentation approach, in comparison to classical active contours (cf. Fig. 3g). In fact, in the latter case, while taking into account only the image local features, such as gray level image intensities, gradients, and curvature, the classical active contours method considers the target object as the overall image constituting not only of the true image, but also, of the occluded regions that are artificially created in zones around image edges. This is the principal reason of why the final classical active contours based segmentation (Fig. 3g) is very far away from the ground truth image (Fig. 3a); especially, in zones where the artificial occlusions take place. On the other hand, although one can also notice the neat improvements of the segmentation gained from the proposed approach in Bresson et al. (2006) (Fig. 3f), result is much more better with our proposed model (Fig. 3e), due mainly to the good behaviors of the functional F4 in Eq. (10) (see Section 3.3.4) in the overall segmentation energy. Indeed, its minimization takes into account the length of the zero level of the shape priors, which allows then a better handling of contours. Concretely, the evolution of the active contour based on our approach is first achieved through the effects of the energy F4, which then makes it move towards the target object areas; then, based on the forces generated by the energy F2 (Eq. (8)), it evolves towards the edges of the object of interests. The progression of the active contour into the directions of the zones of the object interests, is then performed quickly thanks to the functional F3 (Eq. (9)), while the good evolution of the corresponding embedded level set function is guaranteed by the functional F1 (Eq. (7)), by maintaining its SDF property. Parameters

we use in the simulations are: a = 3, n = 1,c = 1/3, b = 0.1, m = 10, l = 50, f = 1 and dt = 0.4. We also apply the proposed model on a different kind of occluded image (Fig. 4a); gray level image intensities in occlusions are different from the object of interest. Segmentations results are shown in Fig. 4. The same remarks previously drawn for the first synthetic image, also hold true for this last synthetic image. The proposed segmentation model behaviors are also illustrated on a synthetic image having some missing parts created artificially (Fig. 5a). Segmentation results obtained with our method and proposed approaches in Li et al. (2005) and Bresson et al. (2006) are displayed in Fig. 5. Results show again the benefits of our segmentation method on the accuracy of the segmentations (Fig. 5d), in comparison, on one hand and with no surprise, to the classical active contours as in Li et al. (2005) (Fig. 5f), but also to the proposed approach in Bresson et al. (2006) (Fig. 5e). Let us just mention that the latter segmentation method is much more better than the classical active contours, for the same reasons previously discussed. 4.2. Digitally reconstructed radiographs Contrary to the case of synthetic images, a primarily and important step consists in the alignment of the training dataset contours. It is really a crucial stage, and the PCA results will strongly rely on it. In fact, it is clear that if one aims at analyzing the shapes deformations of a set of data objects, a minimum requirement in order to correctly do that will be to make sure of the shapes similarity, at

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

173

Fig. 4. (a) Synthetic occluded image test: gray level based occlusions on edges artificially created from the original image displayed in Fig. 3a. (b) Zero level set initialization. (c) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (d) Our approach. (e) Proposed method in Bresson et al. (2006). (f) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. (a) Synthetic image test: missing parts are taken from the original image displayed in Fig. 3a. (b) Zero level set initialization. (c) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (d) Our approach. (e) Proposed method in Bresson et al. (2006). (f) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

least in regions where one would like to learn the shapes deformations. Therefore, if some objects (without loss of generality let us consider two objects O1 and O2) were not aligned at all or not well aligned, a direct consequence would be to get some counterparts of the object O1 be covered or/and masked by other structures of the object O2, or vice versa, and then, the study of the shape deformations would be inaccurate in regions where structures were contiguous. To be more concrete, let us take the example of two femora, which are not aligned at all or not well aligned. In this case, bones could overlap each other, and then, a region of the femoral head could mask the small or greater trochanter, which could result in an incorrect shape analysis in those anatomical structures. We

overcome this issue by performing in advance an alignment procedure of the contours training dataset {Cn}16n6N. Without loss of generality, the alignment is achieved as regards as the first contour C1. It is well known that shape alignment and point-wise correspondence problems are difficult to solve, and become hectic issues; especially, if shapes are represented explicitly by (parametrized) contours. Hopefully, finding correspondences had been intensively studied (see e.g., Zhang, 1994; Veltkamp and Hagedoorn, 1999; Sebastien et al., 2001; Mori et al., 2005; Huang et al., 2006). In this study, we follow Veltkamp and Hagedoorn (1999) and Diop et al. (2010) to propose an alignment process which consists in finding the optimal parameters for a rigid transformation

174

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

between objects to be aligned (Cn, n – 1) and the reference object, which is C1 here. Thus, it comes to find for all n 2 s2,hNt thei scale factor sn, rotation angle hn and translation vector T n ¼ T nx ; T ny , such that the aligned object denoted by C al n could write as:

C al n ¼ s n Rhn C n þ T n :

ð29Þ

The similarity measure is then the minimum of the following quantity:

    al m C 1 ; C al ¼ area of A1 [ Aal n  A1 \ An ; n

ð30Þ

where An represents the area of the surface enclosed by the contour Cn. Let us mention that this measure was also used in Chen et al. (2002) and Bresson et al. (2006). Once the alignment process of the contours training basis {Cn} is performed, the next step is the construction of the SDFs basis {/n}, for which the /n are the embedded functions of corresponding aligned contours C al n . And then, we apply, as for the case of synthetic images, the proposed segmentation model; the training basis we use and a priori shapes are built and designed in the way presented in Section 3.1. In Fig. 6, we show segmentation results obtained on a DRR image with our approach, proposed one in Bresson et al. (2006) and the classical active contours method of Li et al. (2005). As expected, shape prior-based segmentation method gives better results (Fig. 6d and e) than the classical active contours (Fig. 6f). Indeed, in the latter approach, the active contour, only on the basis of local image features, is rather attracted by the lateral epicondyle; this region having stronger image gradients. A more accurate segmentation is provided with our approach, though; especially, in regions nearby the lateral epicondyle and the lesser trochanter. Parameters used in the simulations are as follows: a = 2.6, n = 1, c = 5, b = 0.3, m = 10, l = 40, f = 1 and dt = 0.07. It is important to point out the fact that the DRR image (Fig. 6a) is truly a projection in a 2D space of a 3D volume, which yields then bi-planar images (here we have only the profile view). This explains the variations of gray level image intensities in the DRR of the femur bone, which are in fact occlusions phenomena. For a better illustration of that, we display in Fig. 7 the edge detector function g (Eq. (16)) of the DRR. 4.3. Real radiographic bi-planar images As in the case of DRRs images (Section 4.2), the first step for the segmentation of real bi-planar radiographs consists in aligning the training contours. Our approach is the same as in the preceding section, i.e. it is based on Eqs. (29) and (30). Once the contours alignment procedure is done, we build the SDFs {/n} basis components, and then apply the segmentation algorithm like before.

In order to better show the robustness of the proposed method, we first apply it on a real image radiograph containing the femur bone, for which all other anatomical structures are removed by using a mask. For that example, other test images and of course for real radiographic images for which the target anatomical structure is the femur bone, the training basis we use and a priori shapes were already built and designed previously in Section 3.1, and were displayed in Fig. 2. The masked image is first binarized, before we use the real femur gray level image intensities. This first test is to ensure the expected good behaviors of the segmentation model. Results are displayed in Fig. 8, and as expected, segmentation results are similar to ones carried out with Li et al. (2005) classical active contours approach. Parameters used to obtain the segmentation are the followings: a = 0.2, n = 22, c = 0.85, b = 0.5, m = 60, l = 40, f = 1 and dt = 0.4. Secondly, we use the true gray level values of the femur, instead of the binary intensities, as in the earlier example. Results are displayed in Fig. 9, and show same segmentation results for the different approaches. Fig. 9d illustrates how the prior shape (green curve) guides the active contour (red curve), but also, how the latter one takes the maximum of advantage on the local image properties to achieve the femur segmentation, and in a better way than the prior shape (the final position of the zero level of the prior is a little deviated from the femur contours). Let us also draw the remark that our approach needed only 36 iterations, while Bresson et al. (2006) method required 60 iterations to reach the stationary case, and furthermore we put much more weights (m = 100) in the Mumford–Shah energy functional without the term related to the Hausdorff measure (cf. Eq. (10) without second term), compared to m = 60 with our approach. Without these changes, the greater trochanter could not be segmented though. On the other hand, that issue is avoided with the proposed method, by adding the Hausdorff measure as for another constraint, which is then taken into account in the overall segmentation. Parameters we used in this case are the followings: a = 1, n = 21, c = 0.85, b = 0.5, m = 60, l = 40, f = 1 and dt = 0.4. Next, we show the robustness and efficiency of our approach by applying it on femur images with some missing parts. For doing so, we deal with three different kinds of missing anatomical structures in the femur image: first, we consider some missing parts in the linea aspera (Fig. 10a); afterwards, missing parts appear in the medial epicondyle (Fig. 11a); and last, we have in the same femur image, missing structures in both the femoral head, the linea aspera and the medial epicondyle (Fig. 12a). Results displayed respectively in Figs. 10–12 definitely show the importance of a priori information about the target object shape in the segmentation process. On one hand, incorporating prior shapes in the segmenta-

Fig. 6. (a) DRR of femur bone. (b) Zero level set initialization. (c) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (d) Our approach. (e) Proposed method in Bresson et al. (2006). (f) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

175

Fig. 7. Edge detector function g.

Fig. 8. (a) Real image radiograph on which a binary mask is applied. (b) Binary masked image. (c) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (d) Our approach. (e) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tion helps fill up gaps both in the linea aspera (Fig. 10b) and in the medial epicondyle (Fig. 11b), so that the overall segmentations are the whole femur bone. Nonetheless, in the case where many gaps appear in the image (Fig. 12a), despite that a prior information, the femur could not be fully recovered by Bresson et al. (2006) approach; especially, the femoral head (Fig. 12c), in comparison to the proposed segmentation model (Fig. 12b). Contrary to latter approaches, without that a priori information, a classical active contours method clearly fails to correctly segment the whole femur in all cases (Figs. 10c, 11c and 12d). Notice the topological changes in the classical active contours results due to the level set formulation whose benefits were discussed earlier in the paper. In the following last examples, we deal with real EOS bi-planar image radiographs without any masking techniques: in the first image, the object of interests is the femur bone (Fig. 13a), and in the last one, the target object is the tibia bone (Fig. 16a). For both

cases, it is clear that a classical active contours method (e.g., Li et al., 2005) will be indubitably inefficient (Figs. 13d and 16d) as for a segmentation method, even if the initialization is very close to the femur contours, as we did here. With no surprise gives Bresson et al. (2006) shape prior-based approach better results (Figs. 13c and 16c) than a classical active contour method (Figs. 13d and 16d). The proposed segmentation model gives better segmentations and outperforms the classical active contour approach, as shown by the results (Figs. 13b and 16b). In fact, in the femur segmentation displayed in Fig. 13c, one can see some left over parts; particularly, in regions nearby the lateral epicondyle, the medial condyle and the greater trochanter, in comparison to the proposed segmentation model (Fig. 13b). For instance, we zoom on those regions we just discuss to highlight areas where our proposed method (see Fig. 15) does better than Bresson et al. (2006) approach (see Fig. 14). Parameters values we used to get the segmentation, are

176

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

Fig. 9. (a) Real image radiograph on which a mask is applied. (b) Image of femur bone after masking. (c) Red: active contour evolution with our approach. Green: shape prior evolution. Segmentation results obtained with: (d) Our approach. (e) Classical active contour method Li et al. (2005). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. (a) Real masked image radiograph with missing parts on the linea aspera. Segmentation results obtained with: (b) Our approach. (c) Classical active contour method Li et al. (2005).

the followings: a = 0.2, n = 19, c = 0.85, b = 0.5, m = 250, l = 40, f = 1 and dt = 0.4. The segmentation of the tibia bone is in fact more complicated due, in addition to its complicated shape, to the superposition of the fibula and the tibia in some zones (see Fig. 16a). Let us mention that the image from which the tibia image is cropped is different from the previous femur image.

Fig. 11. (a) Real masked image radiograph with missing parts on the medial epicondyle. Segmentation results obtained with: (b) Our approach. (c) Classical active contour method Li et al. (2005).

Our training dataset is composed of N = 18 shapes embedded as SDFs. The PCA yields five principal modes with a fitting percentage of 98%. The tibia segmentations are shown in Fig. 16. One can notice in the Bresson et al. (2006) segmentation result (Fig. 13c) that there are few structures left over in the segmentation; mainly, in the distal tibiofibular joint and the lateral condyle. Those regions

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

177

Fig. 12. (a) Real masked image radiograph with missing parts on femoral head, linea aspera and medial epicondyle. Segmentation results obtained with: (b) Our approach. (c) Proposed method in Bresson et al. (2006). (d) Classical active contour method Li et al. (2005).

Fig. 14. (a) Segmentation result of femur bone (Fig. 13a) with Bresson et al. (2006) method. (a) Zooms respectively on greater trochanter (yellow), lesser trochanter (lime) and lateral epicondyle (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 13. (a) Real profile radiograph containing the target femur bone. Segmentation results obtained with: (b) Our approach. (c) Proposed method in Bresson et al. (2006). (d) Classical active contour method Li et al. (2005).

are better handled with the proposed approach (Fig. 16b). For instance, we also zoom on those areas to better show the efficiency of our approach (see Fig. 18), in comparison to Bresson et al. (2006) method (see Fig. 17). Parameters values used to get the tibia segmentation are the followings: a = 0.2, n = 9, c = 0.8, b = 0.4, m = 350, l = 40, f = 1 and dt = 0.4. 4.4. Segmentations evaluation This last section is dedicated to the quantitative evaluation of segmentation results obtained with our approach and methods

proposed in Bresson et al. (2006) and Li et al. (2005). This task is really difficult, because it is not obvious to find well suitable metrics or measures in order to compare contours represented as a set of points. Most of current metrics are not very robust to noise and outliers, and sometimes, one can get small errors while segmentation results are far away from the ground truth. And conversely, one can also get a huge error just for a tiny deviation of the segmentation result from the gold star. Therefore, this quantitative evaluation section should not be taken alone, but must be completed by the qualitative results obtained in preceding sections. For instance, let us quote the study proposed in Chabrier et al. (2008), where authors compared some well-known supervised evaluation criteria used to quantify the accuracy of contours carried out with some edge detector algorithms. Other scalable discrepancy measures were also proposed in Odet et al. (2002), where authors took into account both under and over detected points, along with their relative position. Different metrics are used in this work in order to achieve the segmentation evaluations. The first ones are the usual L2 norm and the Minimum Squared Error (MSE). We also work with the commonly used Hausdorff (referred to as CH) metric (see e.g., Baddeley, 1992). As pointed out in Baddeley (1992), Lp norms are inappropriate to express human vision, and on the other hand, Hausdorff metric is very sensitive to noise and even to changes in a single pixel. Therefore, we also use the modified Hausdorff (MH) distance as proposed in Dubuisson and Jain (1994). In order to better evaluate segmentation results, we also use the well

178

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

Fig. 17. (a) Segmentation result of tibia bone (Fig. 16a) with method proposed in Bresson et al. (2006). (a) Zooms respectively on lateral condyle (yellow) and tibiofibular joint (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 15. (a) Segmentation result of femur bone (Fig. 13a) using our approach. (a) Zooms respectively on greater trochanter (yellow), lesser trochanter (lime) and lateral epicondyle (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16. (a) Real profile radiograph image containing the tibia bone. Segmentation results obtained with: (b) Our approach. (c) Proposed method in Bresson et al. (2006). (d) Classical active contour method Li et al. (2005).

known Kullback–Leibler (KL) divergence (see Kullback and Leibler, 1951; Kullback, 1959). To better account both the ground truth and the segmented image contours, we use the modified KL (MKL) divergence obtained by symmetrizing the KL divergence (see Rao

Fig. 18. (a) Segmentation result of tibia bone (Fig. 16a) with our approach. (a) Zooms respectively on lateral condyle (yellow) and tibiofibular joint (orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

179

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181 Table 1 Segmentation evaluations with methods based on the proposed one (referred to as New), ones presented in Bresson et al. (2006) (BAL) and Li et al. (2005) (LAL). Methods

L2

MSE (103)

CH

MKL

NSJ (103)

SJ (102)

FOM

Fig. 3 New BAL LAL

16.93 18.23 250.03

1.03 1.11 15.3

2 2 18

0.72 0.58 5.35

2.51 3.34 588.5

4.65 5.43 624.7

1.85 2.19 690.7

0.9 0.92 0.57

Fig. 5 New BAL LAL

26.87 24.96 247.28

1.64 1.53 15.1

4 5 71

1.5 0.95 25.3

6.21 6.42 574.0

4.14 7.78 422.04

2.41 1.98 140.7

0.77 0.86 0.64

Fig. 4 New BAL LAL

22.97 17.70 245.85

1.4 1.1 15.0

7 7 41

4.25 4.87 9.37

4.43 3.04 568.0

3.4 3.0 611.6

1.8 1.0 459.1

0.84 0.91 0.51

Fig. 6 New BAL LAL

25.81 26.47 253.0

0.2 0.21 2.02

11 11 78

4.82 4.88 29.33

6.26 6.59 601.6

11.6 12.3 623.0

2.88 3.05 221.7

0.92 0.91 0.69

Fig. 10 New BAL LAL

30.75 29.9 251.4

0.36 0.35 2.9

38 38 32

12.6 12.7 14.7

8.89 8.40 594.0

8.54 9.41 502.8

2.94 2.79 193.6

0.80 0.83 0.79

Fig. 11 New BAL LAL

29.28 29.36 252.8

0.33 0.34 2.93

38 38 54

12.65 12.66 14.94

8.06 8.11 600.5

8.54 8.56 541.8

2.67 2.69 228.8

0.83 0.83 0.84

Fig. 12 New BAL LAL

32.75 35.65 249.2

0.37 0.41 2.88

38 38 59

12.6 12.7 16.29

10.1 12.0 583.7

9.6 14.3 446.7

3.34 4.05 158.9

0.78 0.79 0.63

Fig. 13 New BAL LAL

25.05 28.10 252.7

0.29 0.32 2.93

35 35 21

11.29 11.41 2.39

5.90 7.42 600.36

5.77 8.34 575.5

1.93 2.47 285.3

0.85 0.84 0.75

Fig. 16 New BAL LAL

26.39 28.52 250.7

0.38 0.42 3.7

30 30 17

9.88 9.88 3.34

6.55 7.65 590.67

10.5 12.6 602.35

2.78 3.27 226.9

0.87 0.85 0.61

MH

and Nayak, 1985). In addition to KL, we also work with two other measures brought from information theory: the symmetric Jensen-like (SJ) and non-symmetric Jensen-like (NSJ) divergences (see Michel et al. (1994)), which are based on the class of Rényi entropies (see Rényi, 1960). The last metric we use for segmentations evaluation, is the most popular and widely used Abdou and Pratt (1979) Figure of Merit (FOM). FOM is not symmetric and 0 6 FOM(I, J) 6 1. Error types of that measure are studied and well illustrated in Baddeley (1992). Using above presented measures, we report in Table 1 quantitative results obtained with the different segmentation methods, which confirm the performances of shape prior based approaches, and also the good behaviors of the proposed algorithm. For instance, the most adequate measures that really express obtained qualitative results are the L2 norm, the MSE and metrics based on information theory (errors are expressed in logarithm, so they could be important!); especially, MKL, SJ and NSJ; the lower the error, the better the segmentation in the sense of accuracy. Other measures give sometimes very good results, but fail in some cases to express the qualitative results, e.g., we obtain with CH (see Baddeley, 1992) the lowest errors by using the Li et al. (2005) classical active contour method (Table 1, for Figs. 10, 13 and 16), which is in fact definitely not true! 5. Discussion and perspectives We have proposed here a new image segmentation method that is based on active contours and formulated with a level set ap-

proach. In fact, despite noticeable developments in imaging techniques, most of medical images are suffering from poor contrast image intensities and occlusions phenomena; especially, in bones joints. Therefore, it is very often in medical images that some anatomical structures are missing, or in other cases, the latter ones are not well displayed or do not appear at all. In addition to occlusions phenomena and missing parts in images, one has to deal with the complicated forms of anatomical structures, which make hard their segmentations. As illustrated in the numerical results in Section 4, it is clear that classical active contours methods definitely fail to segment desired structures. In this study we follow the known way to resolve segmentation issues, that is, to incorporate a priori information about the form of the target object in such a way that the zero level set is also guided, in addition to intrinsic image features, by a set of a priori admissible and deformable shapes, which are then taken into account in the variational formulation of the segmentation problem. The functional we proposed in Section 3 is composed of four terms. The first one (cf. Eq. (7)) is designed so that the level set keeps the important signed distance function property, which is necessary to guarantee its good evolution. In doing so, we avoid the re-initialization process which is in fact a crucial step in classical level set-based works (see e.g. Sussman et al., 1994; Zhao et al., 1996; Chan and Vese, 2001), and is actually solved independently by the resolution of the classical PDE (11). Issues of such an approach were discussed in Section 3.3.1. The second functional (cf. Eq. (8)) holds a priori information about the admissible deformable shapes of the target object, and is integrated in the level set

180

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181

evolution. An energy that drives rapidly the level set towards the edges of the object of interests, is proposed in the third term (cf. Eq. (9)), which can be interpreted as a weighted area term of the target object. A last term applied on deformable prior shapes is defined (cf. Eq. (9)) on the basis of a complete and modified Mumford–Shah energy functional. Our segmentation model, obtained by solving related Euler–Lagrange equations associated to the minimization of our functional (Eq. (6)), is presented in Section 3.4. The robustness and efficiency of the approach are shown in Section 4 on synthetic images (Section 4.1), digitally reconstructed radiographs (Section 4.2) and real radiographic images (Section 4.3), along with quantitative evaluations of segmentation results, which were reported in Table 1. Both qualitative and quantitative results confirm the good behaviors of the proposed algorithm. However, despite the obtained good results on the presented different kinds of images, the approach should be improved again for more real images, as depicted by Fig. 16. The segmentation of the tibia bone is in fact really tedious, because of gray level intensity issues, its complicated shape, the superposition of the fibula and the tibia in some regions. Nonetheless, the segmentation carried out with our method is better than compared approaches, and is really good, but not perfect in comparison to the ground truth. Indeed, most of occlusion problems are resolved, except in a tiny zone where the level set is attracted a little bit by the fibula (Fig. 16). This could be resolved by adding other constraints in the overall segmentation energy, for example with Thiruvenkadam et al. (2008) spatial order estimation approach, in order to better handle intersection regions where occlusions appear or are supposed to take place. This is very challenging, because, as remarked in Thiruvenkadam et al. (2008), the application of the proposed work was in fact restrictive to many real images, due to some assumptions in the model. In this work also, the study of the minimum of the proposed functional (Eq. (6)) is not done, as it was the case in some works (see e.g., Chen et al., 2002; Bresson et al., 2006) by proving the existence of a minimizer in the space of Bounded Variations functions. Actually, the theoretical study of our functional, along with to add the spatial order constraints are left for future works. Acknowledgements Authors are really grateful to Dr. Julien Leboucher from Telecom Bretagne/LaTIM for kindly providing the femur DRR image, training datasets for DRRs, real femur and tibia images, and also for providing the femur and tibia ground truth segmentations, which were used in the segmentation evaluations. We also thank Dr. Xavier Bresson from City University of Hong Kong for his help on numerical implementations. References Abdou, I.E., Pratt, W.K., 1979. Quantitative design and evaluation of enhancement/ thresholding edge detectors. Proceedings of the IEEE 67, 753–763. Arnold, V.I., 1983. Geometrical Methods in the Theory of Ordinary Differential Equations. Springer-Verlag. Aubert, G., Kornprobst, P., 2002. Mathematical problems in image processing: partial differential equations and the calculus of variations. Applied Mathematical Sciences, vol. 147. Springer-Verlag, Inc., New York, NY. Baddeley, A.J., 1992. Errors in binary images and an Lp version of the Hausdorff metric. Nieuw Archief voor Wiskunde 10, 157–183. Bardinet, E., Cohen, L.D., Ayache, N., 1996. Tracking and motion analysis of the left ventricle with deformable superquadrics. Medical Image Analysis 1, 129–149. Boykov, Y., Kolmogorov, V., 2004. An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision. IEEE Transactions on Pattern Analysis and Machine Intelligence 26, 1124–1137. Bresson, X., Vandergheynst, P., Thiran, J.P., 2006. A variational model for object segmentation using boundary information and shape prior driven by the Mumford–Shah functional. International Journal of Computer Vision 68, 145– 162.

Caselles, V., Kimmel, R., Sapiro, G., 1997. Geodesic active contours. International Journal of Computer Vision 22, 61–79. Chabrier, S., Laurent, H., Rosenberger, C., Emile, B., 2008. Comparative study of contour detection evaluation criteria based on dissimilarity measures. EURASIP Journal on Image and Video Processing, 13. Chan, T., Zhu, W., 2003. Level Set Based Shape Prior Segmentation. Technical Report. Department of Mathematics, UCLA. Chan, T.F., Vese, L.A., 2001. Active contours without edges. IEEE Transactions on Image Processing 10, 266–277. Chang, H., Yang, Q., Parvin, B., 2009. A Bayesian Approach for Image Segmentation with Shape Priors. Technical Report. Lawrence Berkeley National Laboratory. Chen, Y., Tagare, H.D., Thiruvenkadam, S., Huang, F., Wilson, D., Gopinath, K.S., Briggs, R.W., Geiser, E.A., 2002. Using prior shapes in geometric active contours in a variational framework. International Journal of Computer Vision 50, 315– 328. Courant, R., Hilbert, D., 1953. Methods of Mathematical Physics, vol. 1. Interscience Publishers, Inc., New York. Cremers, D., Kohlberger, T., Schnörr, C., 2002. Nonlinear shape statistics in Mumford–Shah based segmentation. In: ECCV. Springer, Copenhagen, Denmark, pp. 93–108. Cremers, D., Sochen, N., Schnörr, C., 2003. Towards recognition-based variational segmentation using shape priors and dynamic labeling. In: Scale-Space Methods in Computer Vision. Springer, Isle of Skye, pp. 388–400. Diop, E.H.S., Ba, S.O., Jerbi, T., Burdin, V., 2010. Variational and shape prior-based level set model for image segmentation. In: International Conference of Numerical Analysis and Applied Mathematics (ICNAAM), Rhodes, Greece, pp. 2139–2142. Dubuisson, M.P., Jain, A.K., 1994. A modified Hausdorff distance for object matching. In: International Conference on Pattern Recognition, Jerusalem, Israel, pp. 566– 568. Evans, L.C., 1998. Partial differential equations. Graduate Studies in Mathematics, vol. 19. American Mathematical Society, Providence, Rhode Island. Gomes, J., Faugeras, O., 2000. Reconciling distance functions and level sets. Journal of Visual Communication and Image Representation 11, 209–223. Huang, X., Paragios, N., Metaxas, D.N., 2006. Shape registration in implicit spaces using information theory and free form deformations. IEEE Transactions on Pattern Analysis and Machine Intelligence 28, 1303–1318. Jerbi, T., Burdin, V., Stindel, E., Roux, C., 2009. Registration of low dose bi-planar acquisitions for motion analysis. In: IEEE EMBS, Minnesota, USA, pp. 90–93. Kass, M., Witkin, A., Terzopoulos, D., 1987. Snakes: active contours models. International Journal of Computer Vision 1988, 321–331. Kullback, S., 1959. Information Theory and Statistics. John Wiley and Sons, New York. Kullback, S., Leibler, R.A., 1951. On Information and sufficiency. The Annals of Mathematical Statistics 22, 79–86. Laporte, S., Skalli, W., Guise, J.A.D., Lavaste, F., Mitton, D., 2003. A biplanar reconstruction method based on 2D and 3D contours: application to the distal femur. Computer Methods in Biomechanics and Biomedical Engineering 6, 1–6. Leventon, M.E., Grimson, W.E.L., Faugeras, O., 2000. Statistical shape influence in geodesic active contours. In: CVPR, South Carolina, USA, pp. 316–323. Li, C., Xu, C., Gui, C., Fox, M.D., 2005. Level set evolution without re-initialization: a new variational formulation. In: IEEE CVPR, Washington DC, USA, pp. 430–436. Malladi, R., Sethian, J.A., Vemuri, B.C., 1995. Shape modeling with front propagation: a level set approach. IEEE Transactions on Pattern Analysis and Machine Intelligence 17, 158–175. Michel, O., Baraniuk, R.G., Flandrin, P., 1994. Time-frequency based distance and divergence measures. In: IEEE International Symposium on Time-Frequency and Time-Scale Analysis, Philadelphia, USA, pp. 64–67. Mori, G., Belongie, S., Malik, J., 2005. Efficient shape matching using shape contexts. IEEE Transactions on Pattern Analysis and Machine Intelligence 27, 1832–1837. Mumford, D., Shah, J., 1989. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics 42, 577–685. Odet, C., Belaroussi, B., Benoit-Cattin, H., 2002. Scalable discrepancy measures for segmentation evaluation. In: IEEE International Conference on Image Processing (ICIP), New York, USA, pp. 785–788. Osher, S., Fedkiw, R., 2003. Level Set Methods and Dynamic Implicit Surfaces. Applied Mathematical Sciences. Springer. Osher, S., Paragios, N., 2003. Geometric Level Set Methods in Imaging, Vision, and Graphics. Springer-Verlag, Inc., New York. Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics 79, 12–49. Peng, D., Merriman, B., Osher, S., Zhao, H., Kang, M., 1999. A PDE-based fast local level set method. Journal of Computational Physics 155, 410–438. Rao, C.R., Nayak, T.K., 1985. Cross entropy, dissimilarity measures, and characterizations of quadratic entropy. IEEE Transactions on Information Theory IT-31, 589–593. Rényi, A., 1960. On measures of entropy and information. In: Berkeley Symposium on Mathematical Statistics and Probability, Berkeley, California, pp. 547–561. Rousson, M., Paragios, N., 2002. Shape priors for level set representations. In: ECCV, Copenhagen, pp. 78–92. Sebastien, T.B., Klein, P.N., Kimia, B.B., 2001. Alignment-based recognition of shape outlines. In: Arcelli, C., Cordella, L.P., di Baja, G.S. (Eds.), International Workshop on Visual Form. Springer, Capri, Italy, pp. 606–618.

E.H.S. Diop, V. Burdin / Medical Image Analysis 17 (2013) 165–181 Sussman, M., Fatemi, E., 1999. An efficient interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow. SIAM Journal on Scientific Computing 20, 1165–1191. Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics 114, 146–159. Thiruvenkadam, S.R., Chan, T.F., Hong, B.W., 2008. Segmentation under occlusions using selective shape prior. SIAM Journal on Imaging Sciences 1, 115–142. Turk, M., Pentland, A., 1991. Eigenfaces for recognition. Journal of Cognitive Neuroscience 3, 71–86. Veltkamp, R., Hagedoorn, M., 1999. State-of-the-Art in Shape Matching. Technical Report UU-CS-1999-27, Utrecht University.

181

Wang, Y., Staib, L.H., 1998. Elastic model based non-rigid registration incorporating statistical shape information. In: Medical Image Computing and ComputerAssisted Intervention (MICCAI), Cambridge, MA, USA. pp. 1162–1173. Zhang, Z., 1994. Iterative point matching for registration of free-form curves and surfaces. International Journal of Computer Vision 13, 119–152. Zhao, H.K., Chan, T., Merriman, B., Osher, S., 1996. A variational level set approach to multiphase motion. Journal of Computational Physics 127, 179–195. Zosso, D., Callennec, B.L., Cuadra, M.B., Aminian, K., Jolles, B.M., Thiran, J.P., 2008. Biplanar 2D-to-3D Registration in Fourier Domain for Stereoscopic X-Ray Motion Tracking. In: Reinhardt, J.M., Pluim, J.P. (Eds.), Medical Imaging 2008: Image Processing.