Fabien Pierre
email: fabien.pierre@univ-lorraine.fr

Mathieu Amendola

Clémence Bigeard

Timothé Ruel

Pierre-Frédéric Villard
email: pierre-frederic.villard@univ-lorraine.fr

Segmentation with Active Contours

Keywords: active contours, image segmentation, medical imaging

Active contours (also known as snakes) have shown their ability to introduce regularity on image segmentation. In contrast with level-set approaches, the active contours techniques based on a contour parameterization are able to maintain the initial topology of the area of interest. For this reason, it has been used in recent medical research for diaphragm segmentation. Most of the on-line codes for 2D/3D segmentation, as well as built-in Matlab toolboxes are based on level-set methods. Moreover, in the literature, the implementation details of active contours methods with meshes in three dimensions are tight, making tedious any reproduction of these techniques. In this paper, we give some details of the implementation of active contours in 2D/3D with meshes, especially about the choice of the use of a 2D/3D mesh and its refinement. We also explore the choice of the parameters with a quantitative study of their influence on the segmentation results. The 3D segmentation method has been applied to CT scan images of the lungs.

Introduction

Medical training is more and more done using computer-based simulators. They offer a wide variety of training procedures without being invasive to any real patient and provide useful metrics to evaluate trainee's performance. Companies like Mentice 2 or Simbionix 3 are well known actors on the market. However those commercial simulators use generic geometries for the anatomy, i.e. based on conceptual organ geometries not based on real patients, so they cannot be used for treatment planing. Recent research work is precisely focusing on filling this gap by building patient-based simulators [START_REF] Villard | Interventional radiology virtual simulator for liver biopsy[END_REF][START_REF] Kerrien | Blood vessel modeling for interactive simulation of interventional neuroradiology procedures[END_REF][START_REF] Mastmeyer | Population-based respiratory 4D motion atlas construction and its application for VR simulations of liver punctures[END_REF]. When other segmentation techniques fail and when deep-learning techniques lack of labeled data, segmenting 3D organs with low contrast is still an open problem. Various approaches exist and stateof-the-art techniques are atlas-based models [START_REF] Maggia | Traumatic Brain Lesion Quantification based on Mean Diffusivity Changes[END_REF] and machine-learning based methods [START_REF] Trullo | Segmentation of Organs at Risk in thoracic CT images using a SharpMask architecture and Conditional Random Fields[END_REF]. They rely on a consequent database that is not always available. In the following we will only focus on techniques efficient on a small amount of data.

Considering that the image is defined on a regular grid of pixels in 2D (resp. voxels in 3D), the aim of the image segmentation is the computation of N sets of pixels (resp. voxels) such that each set corresponds to an area of the image (these sets are called classes). These areas are computed considering contours, shapes, or statistics of their pixels (resp. voxels).

In the literature, two types of image segmentation methods can be distinguished, based on the discretization of the data. The first one is based on the natural regular grid definition of the image, in pixels or in voxels. The second one is based on the discretization of the boundary of the pixels sets. The first type has been widely used in the seminal segmentation approaches, especially with thresholding methods (see, e.g., [START_REF] Sezgin | Survey over image thresholding techniques and quantitative performance evaluation[END_REF] for a review of such techniques). One of the most known methods is Otsu's method, based on histograms [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]. Avoiding the dependency of the pixel discretization, the freely deformable models are based on a discretization of the boundary of the pixel classes. Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF] have introduced such models for image segmentation in a method called snakes, based on the energy minimization. The minimizer provides a trade-off between a regular shape and a data-driven one. The concept has been extended in 3 dimensions by Terzopoulos et al. [START_REF] Terzopoulos | Constraints on deformable models: Recovering 3D shape and nonrigid motion[END_REF].

The level-set approaches have been introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] and popularized for image analysis and computer vision by Malladi et al. [START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF]. In these methods, the boundary of the shape of interest is defined as the level-set of a function, which is considered as the minimizer of an energy. For these approaches, two variants exist, which are variational approaches (see, e.g., [START_REF] Chan | Active contours without edges[END_REF]) and a PDE approach [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]. These methods favor a regular level-set and are able to naturally tackle the issue of topology variations, whereas the free deformable models preserve the initial topology of the shape. While the preservation of the topology can be seen as an advantage, especially when the interest object is known (e.g., the lungs), the free-deformable approach has been extended to manage some changes of topology [START_REF] Mcinemey | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF]. In this paper, we will focus on the free deformable model with fixed topology.

In the following we will focus on the 2D formalization in Section 2, and the numerical computation of a solution in Section 3. It will be then extended to the 3D case in Section 5. The numerical results as well as the influence of the parameters will be studied in Section 4 (2D) and in Section 6 (3D).

Model for Active Contours in 2D

Active contour models aim to extract the boundary of an area contained in an image. To compute it, the main idea consists in the definition of a curve as the minimizer of a functional. This functional is composed of a sum of energies to perform a trade-off between the regularity of the shape and the contours defined by the gradient of the image. In our framework, the source code of the 2D implementation of active contours is in the release2D folder. It is composed of the main Matlab script MAIN EvolutionSnake2D.m and a functions folder, which contains all the utilities functions.

Model Definition

In the following, I denotes an image on a domain Ω ⊂ R 2 . The active contour model is defined as a functional applied to a curve in the image domain Ω. Formally, we aim to compute a path v

v : [0, 1] → Ω ⊂ R 2 s → v(s) = (x(s), y(s)), (1)
such that v(0) = v(1) and v ∈ C 2 . At the position s, the active contour is represented by v(s) = (x(s), y(s)) where x and y are snake point coordinate vectors, defined on Ω. Let us consider the following energy

E(s) = 1 0 E int (s) + E ext (s) ds, (2)
where E int denotes the internal energy and E ext the external energy. The internal energy penalizes non-smooth curves and its minimization favors the regularity of v. This is the a priori part of the functional. The external energy depends on the image and enables the curve to fit the area of interest. This last term can be seen as the data fidelity term of the functional.

Internal energy. To regularize the shape of the snake, the internal energy can be written as

E int = 1 2 α(s) v ′ (s) 2 + β(s) v ′′ (s) 2 , (3)
where α(s) and β(s) are weights, which control the snake regularity. In the following, α and β are considered constant to simplify the problem for a numerical application. The minimization of the term E int enforces the regularity of the curve.

External energy. The external energy represents the data-fidelity term. It is defined as a potential P on the image domain Ω

P : Ω → R (4) (x, y) → P (x, y), (5)
such that the minimum of P is the contour of interest in the image. In case of detection of a black contour in a white background, the potential can be the image itself. For smooth images, we consider that the external energy is only represented by the image gradient. Indeed, the active contour must be brought to high gradient areas (which represent image edges) to make the snake close to the contours. For this purpose, the image gradient has to be as high as possible. The external energy can therefore be represented by E ext = -P (v(s)), where

P (x, y) = ∂I ∂x (x, y), ∂I ∂y (x, y) 2 . ( 6
)
This potential has an analogous definition with the Canny edges detector [START_REF] Canny | A computational approach to edge detection[END_REF].

Global energy. Combining Equations ( 2), ( 3) and ( 6), our problem can be written as

min v 1 0 G(s, v, v ′ , v ′′ ) ds, (7) with
G(s, v, v ′ , v ′′ ) = 1 2 α v ′ (s) 2 + β v ′′ (s) 2 -P (v(s)). ( 8
)
The computation of the optimal condition for the minimization of this energy based on the calculus of the variations is detailed in the next section.

Optimal Condition

Let us now write the optimal condition to minimize the functional [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF].

The aim is to find the active contour which minimizes the previous equation. To solve it, we use Euler-Lagrange's theorem, which is given by Equation ( 9)

min v 1 0 G(s, v, v ′ , v ′′ ) ds ⇔ ∂G ∂v - ∂ ∂s ∂G ∂v ′ + ∂ 2 ∂s 2 ∂ 2 G ∂v ′′ = 0. ( 9
)
By analogy with our problem we write the following equalities (2) (s), and

∂G ∂v = -∇(P (v(s))), ∂ ∂s ∂G ∂v ′ = αv
∂ 2 ∂s 2 ∂ 2 G ∂v ′′ = βv (4) (s). ( 10
)
Summing up, we obtain the optimal condition ∀s ∈ [0, 1], αv (2) (s) -βv (4)

(s) -∇(P (v(s))) = 0. ( 11
)
Follow the numerical computation of the snakes, let us consider the discretization of the curve v in the next section.

Snake Discretization and Numerical Scheme

In the following, a discrete active contour, (as known as snake) V is a discretized version of a continuous path v, composed of N points (also called vertices) indexed from 1 to N , represented by the coordinates (v 1 , . . . , v N ) ∈ (R 2 ) N . The curve may have self-intersections. The snake V is a regular discretization given by v i = v((i -1)/N ). For all i, v i denotes the coordinates (x i , y i ) ∈ R 2 .

Discretization of the snake derivatives. In the case of a discrete active contour (i.e., when Equation ( 11) is discretized), we have to define v (2) and v (4) .

To discretize the derivatives, we use the Euler's method for the second order

∂ 2 v ∂s 2 ≃ v k+1 -2v k + v k-1 , (12)
as well as the fourth order

∂ 4 v ∂s 4 ≃ v k+2 -4v k+1 + 6v k -4v k-1 + v k+2 . (13)
Let us remark that the usual Euler's method is ∂v ∂s ≃ N •(v k+1 -v k ). Nevertheless, removing the factor N of all the computation makes the parameters of the model independent of the resolution of the discretization. Notice that to simplify the presentation and without loss of generality, Equations ( 12) and ( 13) are normalized to the case of constant arc length equal to 1.

These derivatives of discrete snakes can be seen as linear operator on the snakes. As boundary conditions, we expect the active contour to be closed, therefore we define v 0 := v N . With this cyclic boundary conditions, let A 2 denote the tridiagonal matrix (second derivative) and A 4 the pentagonal matrix (fourth derivative)

∂ 2 v ∂x 2 k ≃          -2 1 0 . . . 0 1 1 -2 1 . . . 0 0 0 1 -2 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 0 . . . -2 1 1 0 0 . . . 1 -2                 v 1 v 2 . . . v n-1 v n        = A 2 v, (14)
and

∂ 4 v ∂x 4 k ≃          6 -4 1 . . .
               v 1 v 2 . . . v n-1 v n        = A 4 v. (15)
To simplify the notations, let us denote A = -αA 2 + βA 4 . The computation of the matrix A can be performed by release2D/functions/computeA.m. To help the implementation, A 4 is approximated by A 2 2 (line 42 of computeA.m). The cubic spline refinement detailed in the paragraph "Adaptive discretization" aims at this reparameterization in order to avoid degeneration of the curve.

Temporal discretization for the numerical scheme. Let us denote by V the discretized path such as

(v 1 , . . . , v N ) ∈ (R 2 ) N and ∇(P (V )) ∈ (R 2 ) N , (16)
the gradient of the potential computed in the coordinates given in V . From Equation ( 11), the optimal condition can be now written as

AV -∇(P (V )) = 0. (17)
We aim to discretize temporally the equation to obtain an iterative algorithm. In this way, the snake is going to get closer, step-by-step, to the image contours.

Let V t denote the vector V at time t and let us consider the following fixed-point problem

AV -∇(P (V )) = 0 ⇐⇒ AV -∇(P (V )) + V = V. (18)
Let us denote g(V ) = AV -∇(P (V )) and γ the step size of Euler's temporal discretization. Using Euler implicit method for AV and explicit method for -∇(P (V )), we finally obtain the following numerical scheme

AV t -∇(P (V t-1 )) = - 1 γ (V t -V t-1 ). ( 19
)
By inverting the matrix A + γI d , we can obtain V t depending on V t-1

V t = (γA + I d ) -1 (V t-1 + γ∇(P (V t-1 ))) . ( 20
)
The implementation of Equation ( 20) is done in release2D/functions/EvolutionSnake2D.m at line 147, with a pre-computation of (γA + I d ) -1 at line 120. The properties of the functional [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF] does not enable us to prove the convergence of such scheme. Moreover, since some flatnesses may appear in the functional, some terms have to be artificially added to numerically ensure the evolution of the numerical scheme in some cases.

Balloon force in 2D. To avoid stationary effects during the minimization process due to constant parts of the image, we have added a balloon force as proposed in [START_REF] Cohen | Finite-element methods for active contour models and balloons for 2-D and 3-D images[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF].

The balloon force aims to naturally contract/expand the snake, even if the image potential does not attract it [START_REF] Cohen | On active contour models and balloons[END_REF]. This enables the snake to get close to the image when it could not have been achieved by the regularity forces. This balloon force is built from the curve normal vectors, computed in each point of the active contour from the previous and the following point. For instance, for the point v i , let us consider the tangential vector -→ t i , computed by averaging the previous and the next edges

- → t i = 1 2 ((v i+1 -v i ) + (v i -v i-1 )) = t x t y . (21)
For the balloon force, we naturally choose the orthogonal vector -→ n i to expand the active contour

- → n i = - → t i ⊺ t i = 1 t 2 y + t 2 x -t y t x . ( 22
)
By computing -→ n for each point of the snake, we obtain a force

F B = κ B - → n , (23)
with the sign of κ B controlling the expansion or contraction of the active contours. This force is integrated in the numerical scheme in the following way:

V t = (γA + I d ) -1 (V t-1 + γ [∇(P (V t-1 )) + F B ]) . ( 24
)
This force is computed by release2D/functions/balloonForce.m.

In order to obtain a balloon force oriented outside of the polygon, it is oriented clockwise thanks to the function polygonParity.m which uses a parity test to check the sense of the polygon once the user draws it. In our code, a negative coefficient κ B means that the balloon force is inner-oriented whereas a positive coefficient means an outer orientation of the force.

In the following we detail the discretization of ∇(P (V t-1 )).

Potential discretization and smoothing. To compute explicitly ∇(P (V t-1 )), the image gradient has to be discretized. The following formulas have been taken from [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF]. Let N denote the number of points and g a discrete function and define the gradient

∂g x,y ∂x =            g x+1,y -g x-1,y 2 if 1 < x < N, g 2,y -g 1,y 2 if x = 1, g N,y -g N -1,y 2 if x = N. (25)
Likewise in the y direction we have

∂g x,y ∂y =            g x,y+1 -g x,y-1 2 if 1 < y < N, g x,2 -g x,1 2 if y = 1, g x,N -g x,N -1 2 if y = N. ( 26
)
The computation of this gradient is implemented in release2D/functions/gradientCentered.m. To compute ∇(P (x, y)), we apply these formulas directly on the image with an interpolation available in release2D/functions/interpSnake2.m. We then compute pixel-wise the norm of the gradients. Finally, we apply these formulas again to compute P . To avoid some problems with the noise, the image is blurred with a Gaussian kernel with a standard deviation equal to σ. The parameter of this kernel as well as its influence on the result is discussed in Section 4.

Stopping criterion. The computation of the active contour being an iterative process, it has to be stopped. The model being highly non-convex the stopping criterion can only be heuristic. Our strategy is inspired by the convergent algorithm. When the difference between two iterations V t and V t+1 is small enough, the algorithm is stopped. To compute the difference between the two vectors, the standard L2-norm is used. In practice, when V t -V t+1 ≤ ε, the algorithm is stopped. For practical cases, we have experimentally used ε = 9 • 10 -2 , which is implemented with an interpolation available in release2D/functions/EvolutionSnake2D.m at lines 134 and 176. Since there is no convergence guaranty, the algorithm also stops when a maximum number of iterations is reached.

In Section 4, a discussion about the influence of the smoothness of the image and the time-step parameter is proposed.

Mickey Mouse Ears failure case. In this paragraph, we describe a failure case that may happen with some set of parameters with some specific initialization values. In the case of an initialization with a polygonal snake, the topology of the initial contour can be modified, and some blobs, looking as the ears of Mickey Mouse can appear. For instance in Figure 1, whereas the initial curve is a polygon with no multiple points, the iterative algorithm produces some double points (the curve crosses itself). To understand the origin of this problem, let us display locally the balloon forces around the snake, and focus on the angle of the snake (see, e.g., Figure 2).

One can see the balloon forces are crossed, that will produce a crossing of the points of the curves, dealing with a small blob. Moreover, after one crossing, the balloon force will continue to favor an increase of the blob.

In order to avoid this effect, the smoothing of the balloon force has been proposed. To this aim, the convolution of the balloon force with a normalized 1D Gaussian kernel centered in 0 and with a standard-deviation denoted by σ B is used to regularize this force. This convolution, considered with cyclic boundary conditions, is applied directly on each coordinates of the F B vector defined in Equation ( 23). The smoothing version of

F B = F B,x F B,y is defined as FB,i = (G σ B * F B,x ) i (G σ B * F B,y ) i . ( 27
)
The smoothing of the balloon force is implemented in release2D/functions/smoothForces.m. The computation of the convolution is done through FFT transformation in order to make the circular conditions easier to compute. The balloon force with and without regularization can be seen in Figure 2. With the Gaussian smoothing the crossing of the force is lower. In Figure 3 the smoothing of the balloon force shows some positive effect since it can remove the blob from the resulting snakes.

The computation of the contour following the snake evolution model including a smooth balloon force is described in Algorithm 1. A hack has been added in order to cancel the balloon force when the force P ′ (V ) deriving from the external energy is strong enough.

Adaptive discretization. In order to enable the snake to fit some small structures, a refinement of the discretization of the active contour is proposed. To this aim, the points of the snake are interpolated with cubic splines. These splines are defined as follows: the splines connect two adjacent points of the vector defining the active contours; the splines should be connected smoothly in the sense that the second order derivative of the spline has to be equal in the connecting points. The method and the implementation is widely inspired by [START_REF] Press | Numerical recipes in C++: the art of scientific computing[END_REF], Chapter 3. Let us define a discretized path in R 2 defined by the points

x i y i = x(t i ) y(t i ) , (28)
with i = 1, . . . , N . Let us assume the following conventions: t N +1 = t 1 , x N +1 = x 1 , and y N +1 = y 1 .

In the interval t ∈ [t i , t i+1 ], the linear interpolation can be written as

x(t) y(t) = A x i y i + B x i+1 y i+1 , (29)
with A = t i+1 -t t i+1 -t i and B = 1 -A. In order to get a smooth path, we aim to define the second derivative at the connection points

∂ 2 x ∂t 2 (t i ) = x ′′ i , (30)
∂ 2 y ∂t 2 (t i ) = y ′′ i , (31)
with x ′′ i and y ′′ i given values. As done before, we use the conventions x ′′ N +1 = x ′′ 1 , and y ′′ N +1 = y ′′ 1 . In order to constraint the smoothness of the curve, the interpolation is written as a third order polynomial x(t)

y(t) = A x i y i + B x i+1 y i+1 + C x ′′ i y ′′ i + D x ′′ i+1 y ′′ i+1 , (32)
Algorithm 1: Algorithm to return a polygon given by a snake evolution model.

Input : input image I (image domain: Ω ⊂ R 2 ), initial contour V init ∈ (R 2 ) n Output : final contour V ∈ (R 2 ) n V ← V init //
if FB,i ≤ P ′ (V ) i then C i = f alse if NOT C i then FB,i ← 0 V 1 ← (γA + I d ) -1 V + γ P ′ (V ) + FB,i // Numerical scheme Equation (24) // Cubic spline refinement if nb iteration mod 100 = 0 then V 1 ← refinement of V 1 ,
← nb iteration + 1 if V and V 1 have the same size then error ← V -V 1 2 2 V ← V 1 return V with C = 1 6 (A 3 -A)(t i+1 -t i ) 2 and D = 1 6 (B 3 -B)(t i+1 -t i ) 2 .
The second order derivative is thus a linear interpolation of the constraints when t ∈ [t i , t i+1 ]

∂ 2 x ∂t 2 (t) = Ax ′′ i + Bx ′′ i+1 , (33)
∂ 2 y ∂t 2 (t) = Ay ′′ i + By ′′ i+1 . (34)
We have now 2N unknowns x ′′ i and y ′′ i that we have to compute. To this aim, we assume the continuity of the first order derivative. With some rearrangements, for i = 1, . . . , N

x i+1 -x i t i+1 -t i - x i -x i-1 t i -t i-1 = t i -t i-1 6 x ′′ i-1 + t i+1 -t i-1 3 x ′′ i + t i+1 -t i 6 x ′′ i+1 (35) y i+1 -y i t i+1 -t i - y i -y i-1 t i -t i-1 = t i -t i-1 6 y ′′ i-1 + t i+1 -t i-1 3 y ′′ i + t i+1 -t i 6 y ′′ i+1 , (36)
with the convention t 0 = t N , x 0 = x N , y 0 = y N , x ′′ 0 = x N , and y ′′ 0 = y N . We then get 2N equations with 2N unknowns. Solving this linear system gives us the second order derivative constraints at each connection points. To solve it, the variable t is discretized with a step equal to 1 N . Equation (32) gives then the values of the spline between connecting points. Figure 4 shows an example. The set of points (x i , y i ) are the red circles, which are connected by smooth curves. Figure 4: The cubic splines connect the points (red circles) with smooth (polynomial) interpolations. This curve can be regularly sampled with respect to t (colored dots). Notice that, theoretically, this spline interpolation can change the curve topology, but it is not the case on our practical experiments.

Once the coefficients of Equation (32) have been computed, the splines are uniformly sampled. To this aim, a numerical computation of the length of each spline is done with a left Riemann sum to compute the following path length with a subdivision of the [0, 1] interval

1 0 ∂x ∂t (t), ∂y ∂t (t) t dt. ( 37
)
The left Riemann approximation can be written as

1 N N -1 i=0 ∂x ∂t (i/N ), ∂y ∂t (i/N ) t . ( 38
)
Then, the spline is homogeneously sampled: the complete path defined on [0, 1] is subdivided relatively to the lengths of each spline. To this aim, the interval [0, 1] is uniformly sampled, and for each time step in this interval it can be determined to which spline of the path it corresponds. The point of the path at this time step of the interval can be computed from Equation (32). All these computations are listed in Algorithm 2 and performed in release2D/functions/splinesInterpolation2D.m. The refinement strategy, as well as the influence of the splines refinement on a numerical example, are studied in Section 4.

Algorithm 2: Algorithm to return a polygon with a different number of points after refinement with cubic splines.

Input : Closed path V defined by points (X(i), Y (i)) with i ∈ [0, N ] Output : New path interpolated with cubic splines V 2 with N 2 points of coordinates (X 2 , Y 2 ). Parameters: N 2 the new number of wished points // 1-Compute second order derivatives of the interpolating function, both are computed from a system with 2N equations and 2N unknown Given X and Y , compute X" and Y " as defined in Equation ( 35) and (36). // 2-Compute the arc length of each spline M ← 100 // left Riemann sum resolution. for i = 1 to N , with i + 1 = 1 when i = N do Compute S ′ i with respect to t between (X(i), Y (i)) and (X(i + 1), Y (i + 1)) from the X" and Y " values (the derivative of the spline S i is formally computed from Equation ( 32)). Sample the function S ′ i at M points t j uniformly distributed between 0 and path length(i)

← j S ′ i (t j ) M •N
// Left Riemannian sum to approximate Equation (38).

// 3-Sample the splines to get the points of the new path Compute T 2 , a vector of N 2 uniformly distributed between 0 and N j=1 path length(j) for i = 1 to N , with i + 1 = 1 when i = N do for t ∈ T 2 such as being between i j=1 path length(j) and i+1 j=1 path length(j

) do A ← i+1 j=1 path length(j)-t path length(i+1) B ← 1 -A C ← 1/6(A 3 -A)(i/N -(i -1)/N ) 2 D ← 1/6(B 3 -B)(i/N -(i -1)/N ) 2 2 ← AX(i) + BX(i 1) + CX"(i) + DX"(i + 1) // Equation (32) Y 2 ← AY (i) + BY (i + 1) + CY "(i) + DY "(i + 1) // Equation (32) Add (X 2 , Y 2 ) to V 2 . return V 2

2D Segmentation Results

In this section, we propose to study the influence of the parameters of the snakes model in the two dimensions examples. Let us study the following parameters: the smoothing terms α, β, σ, the time-step γ, and the cubic spline refinement. For all the tested images, the initial contour V t=0 is given by coarsely defining manually various points outside of the target as illustrated by Figure 3(a).

Influence of α and β. In the internal energy defined in (3), the two parameters α and β tune the regularity of the snake shape. The α parameter rules the stretching between the points of the snake, whereas the β parameter rules the second-order derivative of the path.

The influence of the α parameter can be seen in Figure 5. The results have been provided with various α values. The higher α, the lower the perimeter of the curve. This result matches well with the theoretical aspect: when α is high, the differences between the points are penalized, thus they are closer to their neighbors. In case when α is too low, some Mickey Mouse Ears effect can occur as in Figure 3.

In order to see the influence of the β parameter, let us compare two results, one with a large

V t -V t+1 ≤ 9 • 10 -2
and the initialization snake is chosen as in Figure 1.

value and another with a small value. The results of the algorithm with β = 4 and with β = 40000 can be seen in Figure 6. With a small value of β, the contour has sharp edges, since the term based on the second-order derivative doesn't penalize much the curvature of the path. With a large value of β, large curvature values (the angular parts of the path) are penalized, thus the result is smooth. The stopping criterion is defined as V t -V t+1 ≤ 9 • 10 -2 and the initialization snake is chosen as in Figure 3.

Influence of σ and the step time Since there is no convergence guaranty for the numerical scheme, let us numerically study the behavior of the iterative algorithm.

The original problem [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF] is not convex, but the convolution of the potential with a Gaussian filter changes the properties of the energy. Since this filtering has a low-pass effect, the energy is smoother when σ is higher.

(a) γ = 0.0005 and σ = 1.6.

0 500 1000 1500 -12 -10 -8 -6 -4 -2 0 2 4 10 4
(b) γ = 0.0005 and σ = 0.6. Assume that f : R n → R is convex, and differentiable. If

∇f (x) -∇f (y) ≤ L x -y for any x, y then ∇f is Lipschitz continuous with constant L > 0.
The smoothing induced by the Gaussian filtering can be seen as a decrease of the Lipschitz constant of the gradient of the energy, thus it can improve the convergence properties of the numerical scheme (see, e.g. [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF], Theorem 2.2.4). Nevertheless, it modifies the data of the problem by blurring it, since the result may fit worse to the original image contours.

For the experiments of Figure 7, we have used the example of Figure 6. In Figure 7 we have drawn the value of the functional [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF] during the time evolution iterations of the numerical scheme [START_REF] Trullo | Segmentation of Organs at Risk in thoracic CT images using a SharpMask architecture and Conditional Random Fields[END_REF]. Two parameters have some influence into the convergence of the numerical scheme: σ and γ. The influence of the γ parameter is the standard behavior for the gradient descent time-step for nonconvex problems. When γ is small enough and the functional ( 7) is smooth enough (σ = 1.6), the energy of the numerical scheme decreases to an asymptotic value (see, e.g., Figure 7(a)). When σ becomes smaller (=0.6) the functional ( 7) is sharper, the numerical scheme has an energy that decreases and oscillates around a local minimum (see, e.g., Figure 7(b)). When γ is too large, the numerical scheme has a random behavior (see, e.g., Figure 7(c)).

Cubic spline path refinement. In order to refine the path, a sampling of the cubic spline interpolation can be used, as introduced in Section 3. Our strategy consists in the computation of the maximum distance between the points every 100 iterations. If it exceeds 2 pixels, the set of points is interpolated by cubic splines and the new path is sampled with 10% more points.

In Figure 8, the influence of the refinement is shown. The snake does not fit the contours in the concave parts with the original method. With the refinement, the snake is closer to the contour.

Practical example of the lungs. For non binary images, it is helpful to increase the contrast of the images to better fit the contours. To this aim, we have used a linear mapping that saturates 10% of the values of the pixels. With this transformation, the image range is between 0 and 255 for all examples. Since the active contour is sensitive to its parameters, this transformation reduces this sensitivity.

In the following, we propose to segment the lungs with active contours. Figure 9 shows the active contour method to segment lungs. The contour fits well a lung in its smooth part, whereas it is quite inaccurate in the interior part which is a complex contour. It is similar to the results existing in the literature [START_REF] Pinho | Keuhkot: a method for lung segmentation[END_REF] because such accuracy is enough for radiotherapy treatment planning or for training simulators.

1 presents various images with sets of suitable parameters. In the next section, an extension to the 3D case is proposed. Table 1: Various images with a set of valid parameters.

Extension in Three Dimensions

In this work, for the sake of simplicity, the 3D case has been considered by extending the 2D numerical scheme instead of redefining model [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF] with surface integrals as done in [START_REF] Cohen | Finite-element methods for active contour models and balloons for 2-D and 3-D images[END_REF]. In our framework, the source code relative of the 3D implementation of active contours is in the release3D folder. Similarly to the 2D case, it is composed of the main Matlab script MAIN EvolutionSnake3D.m and a functions folder.

3D extension of the 2D optimal conditions. In 2D, the optimal condition, after discretization, reads as Equation ( 17), i.e.

AV -∇(P (V )) = 0.

In this equation, AV represents a smoothing term which ensures the regularity of the snakes, whereas ∇(P (V )) is a data-fidelity term. According to Equations ( 25) and (26), the definition of this term is only based on the pixel representation of the original image. Thus, the extension of this term to the 3D case is obvious, just by extending the definition of the gradient.

The definition of A settles on the representation on the snake itself. In the following, we assume that the natural 3D extension of the snake is a triangular mesh. This mesh is defined as a set of 3D points (vertices) in combination with an adjacency matrix defining the edges of a graph. In the following, we assume that all the faces of this mesh are triangles, but this is not a restriction for our definition of the regularity matrix A.

Computation of the regularity matrix.

Let us define {v i with 1 ≤ i ≤ N } the set of the vertices of the given mesh. Precisely, we consider that v i is a point in a three dimensional space, thus Balloon force in 3D. To extend the 2D computation of the balloon force, the definition of the normal force of the mesh is required. This normal vector can be naturally computed on each face of the mesh. Let us consider a face composed of the vertices v 1 , v 2 and v 3 and let -→ N i denote the area of the i-th face, computed from the directed cross product -→

N i = (v 1 -v 2 ) ∧ (v 3 -v 1 ). (49)
Next, for each vertex, the normal vectors of the adjacent faces are averaged to give the normal vector at this vertex. The balloon force is obtained for every vertex with the same computation as Equation ( 23). This force is computed by release3D/functions/EvolutionSnake3D.m for each of the coordinates.

The triangles are given by an initial mesh from a sphere where the vertices are indexed such as the normals are all in the same direction. While the mesh is evolving, the topology does not change and the normals will still be computed from the same vertex indices.

Meshing and initializing. To initialize the active contour, a sphere is meshed from an octahedron subdivided two times such as new faces are created from the added vertex at the middle of each edge. Every vertex of the sphere is moved such as its distance to the sphere center is equal to the sphere radius. The code to compute this initial mesh is given by trisphere.m. Then we expand step by step the snake following the evolution equations. An illustration of the mesh evolution is given on Figure 10. Contrary to the 2D framework, the number of contour vertices does not evolve with the iterations, therefore the number of triangles remains the same and equal to the number constituting this initial sphere. The influence of the main parameters has been tested. The same behaviors as in the 2D case has been observed. Parameter α acts on the stretching between points to allow the smoothness of the The method has also been tested on other synthesised 3D images. The result is illustrated on Figure 12. They have been produced with the default parameters given at the beginning of this section. The first image (Figure 12 Similarly to the 2D case, our method has been carried out on the practical example of the lungs. Results are illustrated on Figure 13. The parameter set has been left unchanged. Similarly to the synthesized 3D images, the active contours fit well the image. The initial contour is also a sphere located inside the lung.

Conclusion

In this work, we have recalled the general principles of the active contours with meshes. The work has been performed both for 2D images as well as for 3 dimensions images. The discretization scheme and the numerical resolution are detailed. The various engineering hacks that are more or less explained in the literature are detailed in order to reproduce the method. The parameters are given and their effects on the results are discussed, still in 2D and 3D. We have also presented a concrete segmentation case with medical images and the extraction of lungs contours.

Finally, we can conclude that this approach preserves the initial topology of the contour, which can be an advantage when the shape of the object to segment is known. The drawback of this method is the sensitivity to the parameters but if the regularity of the shape does not vary too much, the tuning remains reliable.

  (a) Initialization of the algorithm. (b) Mickey Mouse Ears failure.

Figure 1 :

 1 Figure1: In some conditions, the evolution of the snake can produce some double points on the curve, even if the initialization does not contain one. This image has been produced with α = 1750, β = 700, κ B = -90, γ = 0.00005, σ = 1 and 10000 iterations.

  Field of balloon forces applied to a snake. The same field with a Gaussian smoothing.(c) Zoom on (a). Initial snake (red) and resulting snake (green).(d) Zoom on (b). Initial snake (red) and resulting snake (green).

Figure 2 :

 2 Figure 2: When the balloon force is large enough, and the snake (in red in Figure (c)) is modified with respect to this force, a loop occurs in the resulting snake (in green in Figure (c)). With the Gaussian smoothing, this phenomenon is avoided (see, e.g., (d)).

  (a) Initialization of the algorithm. No smoothing Smoothing parameter= 50 Smoothing parameter= 100 (b) Result after convergence.

Figure 3 :

 3 Figure 3: The smoothing of the balloon force has a valuable effect on the reduction of the blobs. With σ B = 100, there is no blob. Is the smoothing is low σ B = 50 or without smoothing, some blob appears in the bottom-left side. The others parameters are the same as in Figure 1.

Figure 5 :

 5 Figure 5: Influence of the parameter α on the segmentation. The results are produced with β = 500, κ B = -50, γ = 0.0001, σ = 1, σ B = 70. The initialization is displayed on the left-hand side. The stopping criterion is defined as V t -V t+1 ≤ 9 • 10 -2 and the initialization snake is chosen as in Figure 1.

(a) β = 4 .

 4 (b) β = 40000.

Figure 6 :

 6 Figure 6: Influence of the parameter β on the segmentation. The results are produced with α = 17.5, κ B = -50, γ = 0.0002, σ = 1, σ B = 70. The stopping criterion is defined as V t -V t+1 ≤ 9 • 10 -2 and the initialization snake is chosen as in Figure 3.

  γ = 0.05 and σ = 1.6.

Figure 7 :

 7 Figure 7: Value of the energy (7) with respect to the number of iterations in Algorithm 1. We can see the influence of the parameters σ and γ on the decreasing of the energy. The results are performed with α = 750, β = 1000, κ B = -30, σ B = 70 and 1500 iterations.

Figure 8 :

 8 Figure 8: Effect of the discretization refinement on the segmentation result with a shape containing concavities. The results are produced with α = 1500, β = 2000, κ B = -50, γ = 0.0001, σ = 1, σ B = 70 and 10000 iterations.

Figure 9 :

 9 Figure 9: Effect of the discretization refinement on the segmentation result with a shape containing concavities. The result is produced with α = 1500, β = 2000, κ B = -50, γ = 0.0001, σ = 1, σ B = 70 and 10000 iterations.

7 •

 7 10 3 3.7 • 10 3 5 • 10 -2 5 • 10 -4 -10 100 2000

Figure 10 :

 10 Figure 10: Evolution of the 3D contour to match a 3D image of a cube. The image voxels are illustrated on three slices cut in the middle of the image. The mesh of the contour is in red.

5 Figure 11 :

 511 Figure 11: Influence of parameters alpha and beta.

Figure 12 :

 12 Figure 12: Results on synthesised 3D images. The image voxels are illustrated with three slices and the mesh is in red.

  (a) and (b)) is an ellipsoid and the second image is composed of two perpendicular ellipsoids (Figure 12 (c) and (d)). Both the initial and the final steps of the contour evolution are displayed.

  (a) Initialization of the contour for the lung: the initial sphere is located inside the lung (b) Lung contour as a mesh display in the coronal direction (c) Lung contour display with Gouraud lighting in rough perspective

Figure 13 :

 13 Figure 13: Results on medical images of lungs.

  Contour initialization I ← G σ * I, with G σ a 2D Gaussian kernel with standard-deviation σ.

	// Computation of the regularization matrix
	Computing A 2 defined in Equation (14)
	A ← -αA 2 + βA 2 2 P ′ ← ∇ ∇I 2	// Equations (15) and (20) // Contribution of external energy with Equation (6)
	C ← (true) n ∈ ({true, f alse}) n	// Hack : enable the balloon force
	nb iteration ← 1 ; error ← 1		// Stopping conditions
	// Iterative process
	while nb iteration ≤ n max AND error > 9 • 10 -2 do
	FB computed from Equation (27)	// Computation of the balloon forces
	// Hack : disable balloon force when External energy is large enough

Image Credits

The images credits of the images in Table 1 are: shamrock, shield, cat, vase and coin, personal work. The medical image is a slice extracted of an anonimized 3D CT scan from the dir-lab database [3] labeled as "4DCT 1" on their webpage 4 .

. In this section, we assume that the adjacency matrix is defined as follows Λ = (α i,j ) 1≤i≤N, 1≤j≤N ,

with α i,j = 1 if v i neighbor of v j , 0 otherwise.

Let us remark that the number of neighbors of a vertex i is given by summing the i-th column of the matrix Λ.

The difference finite formulas on the curve (see Equation [START_REF] Mastmeyer | Population-based respiratory 4D motion atlas construction and its application for VR simulations of liver punctures[END_REF], and ( 13)) have to be adapted to the mesh. For this purpose, let n i denote the number of neighbors of vertex v i . Those vertices are given by an adjacent matrix within triangulation2adjacency.m. The second order (Equation ( 12)) derivative was defined as a linear combination of the current vertex with its direct neighbors. As an extension, we propose to define the second order gradient on the mesh by the sum of the neighbors minus vertex itself multiplied by the number of neighbors. Thus, the matrix A 2 for the second order is defined as

with a 2,i,j =

This formula is a direct extension of the second order derivative filter [1, -2, 1]. Notice that the computation using matrix A 2 does not take into account the spacing between the vertices. The matrix for the fourth order is simply given by

This matrix operation is easy to compute with Matlab, but in the C++ implementation it can be more easily computed by the following formula

with a 4,i,j =

and V = x y z ⊺ , we get the final numerical scheme, identical to the 2D case

The computation of matrix A is performed by release3D/functions/computeA3D.m. The implementation of Equation ( 48) is done in release3D/functions/EvolutionSnake3D.m for each of the coordinates, at lines 152 to 154, with a pre-computation of (γA + I d ) -1 at line 113.

∇(P (V t-1 )) is computed from the voxel image: the image gradient is discretized for each direction. As in the 2D case, the balloon force can help the convergence of the snake to an image contour.