Existence and uniqueness results for the Gradient Vector Flow and geodesic active contours mixed model

This article deals with the so called GVF (Gradient Vector Flow) introduced by C. Xu, J.L. Prince . We give existence and uniqueness results for the front propagation flow for boundary extraction that was initiated by Paragios, Mellina-Gottardo et Ralmesh . The model combines the geodesic active contour flow and the GVF to determine the geometric flow. The motion equation is considered within a level set formulation to result an Hamilton-Jacobi equation.


Introduction
We consider an image segmentation model : I is a given image and we want to detect boundaries without connexity or convexity assumptions on contours. Therefore we are interested in the Gradient Vector Flow (GVF) as a front propagation flow model. This model builds a class of vector fields derived from images and has been introduced by Chenyang Xu and Jerry L. Prince in [14]. The GVF can be viewed as external forces for active contour models: it allows to solve problems where classical methods convergence fail to deal with boundary concavities. On the other hand, a new front boundary-based geometric flow for boundary extraction was proposed by Paragios, Mellina-Gottardo and Ralmesh in [11,12]. In this model the GVF is used to revise the geodesic active contour model of V. Caselles, R. Kimmel, G. Sapiro [2] resulting on a bidirectional geometric flow. The classical parametric active contour model was proposed by D.Terzopoulos, A.Witkin et M.Kass [13]. It derives from an energy functional minimization. Curves are drawn toward the object boundaries under potential external forces action which can be written as the negative gradient of a scalar potential function derived from images. Other forces (as pressure forces) may be added. However, the performance of methods is limited by unstable initialization process and poor convergence when the boundary is concave. Chenyang Xu and Jerry L. Prince [14] set a new external force, dealing with these limitations.
Most of time snakes model provide a local minimum of the functional cost. So, C.Xu et JL.Prince generalized this model from the balance equation at the equilibrium between internal and external forces (Euler equation of the minimization problem). They replace the standard external force F (p) ext = −∇P x (where P is a potential edges detector), by a more general external force.
This new external force field is called the Gradient Vector Flow (GVF). It includes a divergence-free component and a curl-free component [7,16]. Therefore, this new active contour model cannot be formulated as an energy minimization problem. The external force, denoted V below is introduced via the balance equation that can be written : The resulting parametrized curve solving the equation (1.1) is called " GVF-snake ". The final configuration of a GVF-snake satisfies an equilibrium equation which is not a variational problem Euler equation, since V (x, y) is not an irrotational field.
In next section, we introduce the GVF and present the advantages of such a field in view of an active contour model and we present and give a precise definition of the Gradient Vector Flow. In section 3. we present this model which is inspired by the geodesic active contour model combined with the GVF. Finally, the level set strategy leads to an Hamilton-Jacobi equation that has not been studied yet (to our knowledge): in the last section we give existence and uniqueness results.

The Gradient Vector Flow (GVF)
The GVF V is a 2-dimensional vector field that should minimize the following objective function [14,11].
where Ω is an open, bounded subset of R 2 , u x , u y , v x , v y denotes the spatial derivatives (with respect to x and y) of the field V = (u, v), µ > 0 and f : Ω → R a continuous edge detector.
There are many choices for f : in [15], C.Xu and JL.Prince consider where G σ is a Gaussian kernel, so that the image is filtered. Here we choose detector proposed by R. Deriche and O. Faugeras in [3] f Here | · | denotes the R 2 -euclidean norm. The first term of the functional E is a regularization term whereas the second term is a data-driven component. If f |∇f | 2 is small, the energy is dominated by the first term and we get a slowly varying field. On the hand, when f |∇f | 2 is large, the second term forces V to decrease to ∇f . Therefore V is close to the gradient of the edge map when it is large (this is classical external force for snakes) and does not evolve quickly in homogeneous regions (which increase the snake capture area). The parameter µ govern the tradeoff between the two integrands of the cost functional and should be set according to the image noise level (more noise increases µ) see [15].

On GVF existence
First we precise the notations : Ω is an open, bounded subset of R 2 with C ∞ boundary Γ: it is the image domain. Let V = (u, v) and W = (ξ, χ) in H(Ω) := H 1 (Ω) × H 1 (Ω); then DV (x) = (∇u(x), ∇v(x)) ∈ L 2 (Ω) 4 . The inner product in L 2 (Ω) 4 is A first "definition" of the GVF could be the following: The Gradient Vector Flow field is defined as a solution of the following optimization problem (P) : where the energy functional E(V ) is defined by (2.1). Unfortunately, such a definition is not correct since the minimum is not necessarily attained. The study of E will provide another definition of the GVF as the solution to a decoupled system parabolic partial differential equations. From now and in the sequel we assume that the edge function f verifies (H 1 ) : (Ω) and f ≥ 0 (Ω) and f |∇f | 2 ∈ L ∞ (Ω) which are the "minimal" assumptions in a first step.
The bilinear form π is continuous on H(Ω) for the where C(µ, f ) is a constant that only depends on data. L is obviously linear and continuous on H(Ω). We deduce the continuity of E. ✷

Theorem 2.1
The functional E is Gâteaux-differentiable on H(Ω) and for every V = (u, v) and W = (ξ, χ) in H(Ω) we get Proof -The Gâteaux -differentiability is clear. In addition for every V = (u, v) and W = (ξ, χ) in H(Ω) we have then π is coercive on H(Ω). This implies the coercivity and the strict-convexity of E on H(Ω) and the problem (P) has a unique solution.
Proof -π has been define above and we get therefore π is H 1 (Ω) 2 -coercive on H(Ω) and strictly convex. Therefore E is coercive and strictly convex on H(Ω) as well. It follows that if (2.5) is verified then (P) admits a unique solution. ✷ Let us formally write the Euler-equations of problem (P): assume V * = (u * , v * ) is a solution to (P). By convexity, it is a stationary point and ∇E(V * ) = 0. Let W = (ξ, χ) ∈ H(Ω). Integrating by parts expression (2.4 ) gives where ν(x) is the outer unit normal of ∂Ω = Γ at x. We first suppose that ξ and χ belong to D(Ω) (the space of C ∞ (Ω) functions with compact support in Ω). Then where f x , f y are spatial derivatives of f . With (2.6), we obtain : Finally, if a solution V * = (u * , v * ) to problem (P) exists it must verify : (2.7) The above equations are equilibrium equations if V * realizes the minimum of energy E. However, we cannot ensure the existence of such a minimum. Indeed, assumption (2.5) is not realistic : the same gray level for image I on a significant area implies ∇f = 0. In this case the functional E is a priori non longer coercive and we do not know if (P) has a solution. So, instead of computing the "exact" GVF, minimizer of the functional E, we approach it by a minimizing sequence and we consider it is the equilibrium state of a time evolving vectors fields. The stationary problem becomes a dynamic one : A simple way to impose a motion to the vectors field is to impose the velocity This leads to parabolic partial differential equations. The gradient vector flow is initialized as the gradient of the edge detector f : We obtain the following dynamic formulation which is the GVF suitable definition : The gradient vector flow V = (u, v) is defined by the following decoupled equations respectively verified by each of its coordinates u and v :

GVF properties.
In this subsection, we give regularity properties of solutions to (2.10) and (2.11). Of course, it is sufficient to consider equation (2.10). First, the existence of a unique solution is given by a classical theorem (see for example [9,4,5]). The bilinear form a associated to equation (2.10) is .It satisfies : Note that here a does not depend on t. So we may assert that (2.10) has a unique solution Theorem 2.3 Let T > 0, and assume (H 1 ). Then the GVF (u, v) (solution of (2.10) and Proof -We use a generic regularity result ( [9]). We prove the result for the component u.
So the solution u of ( (2.10) belongs to C 1 ([ε, T ] × Ω) for all ε > 0. Moreover ∇f ∈ C 0 (Ω) according to (H 1 ) and compatibility conditions are satisfied with respect to boundary and initial data. So we may conclude. ✷ The GVF is built as a spatial diffusion of the edge detector f gradient. This is equivalent to a progressive construction of the gradient vector flow starting from the object boundaries and moving toward the flat background. In [15], the GVF is normalized to obtain a more efficient propagation. It is denotedV ( to give the new external force of the geometric flow called GVF-snake model. The velocity of the contour C is given by the equation : C. Xu and J.L. Prince have shown that a such flow is not dependent of initial conditions and deal with concave regions. However, it depends on curve parametrization, cannot manage topological changes, and involves second and fourth order derivatives that are difficult to estimate. The objective of N. Paragios, O. Mellina-Gottardo and V. Ralmesh in [11] is to eliminate these shortcomings by integrating the GVF with the geodesic active contour and implementing it using the level set method [12]. Our aim is to study the Hamilton-Jacobi equation derived from this model. ) and directly defined by the contour evolution velocity. It is based on the remark that the Gradient Vector Flow field after the rescaling refers to the direction that has to be followed to locally deform the contour and to reach the closest object boundaries.
On the other hand, given the fact that the propagation of a contour often occurs along the normal direction, the propagation will be optimal whenV and the unit outward normal ν are colinear. So we choose to project the normalized gradient vector flow onto the outward normal. Then we multiply the velocity by an edge detector function g (that may be different from f ), which represents the contour information. The contour evolution velocity is then given by the equation : where I σ is the filtered image and g where h is defined by (2.2) When there is no boundary information ( |∇I σ | 2 ≪ 1), the contour evolution is driven by the inner product between the Normalized Gradient Vector Flow (NGVF) and the normal direction : it is adapted to deal with concave regions. When the curve reaches the object boundaries neighbourhood ( |∇I σ | 2 ≃ +∞ ) then g ≈ 0 that is, the flow becomes inactive and the equilibrium state is reached.
It is classical to impose a regularity condition on the contour propagation adding a curvature term and a "balloon force" H. The evolution equation becomes where β > 0.

Level set implementation
Problems of topologic changes can be solved using the level set method [10]. The moving 2Dcurve is viewed as the zero level set of a 3D surface which equation is z − Φ(x, y) = 0. We denote Φ t the partial derivative ∂Φ ∂t of Φ towards t.
Theorem 3.1 The evolution of the 3D-surface Φ is described by : ± d the distance between x and Γ ; the positive sign (resp. negative) is chosen if the point x is outside (resp. inside) Γ. Proof -The curve C is the zero set level of Φ : Φ(t, C(x)) = 0. Deriving formally with respect to t gives for almost every x ∈ Ω: The final flow can be decomposed in • (a) : a term that provides propagation regularity aims and shrinks the curve toward the object boundaries, • (b) : a bidirectional flow that moves the curve toward the internal and external objects boundaries, • (c) : an adaptative balloon force that drives the propagation of the curve when the boundary term becomes inactive [12].

Choice of the edge detector f
As mentionned before, we focus on the Gaussian edge detector proposed by R. Deriche et O. Faugeras in [3]: where h is defined by (2.2). We must ensure that assumption (H 1 ) is satisfied for the edge detector f . Therefore, we have to set additional hypothesis the image intensity function that we have denoted I. I is supposed to have compact support included in Ω that is image frame: we decide that I is equal 0 out of the "true image". I ∈ L ∞ (Ω) ⊂ L 2 (Ω) ⊂ L 1 (Ω). On the other hand I / ∈ H 1 (Ω). Indeed, the image gradient norm becomes infinite at objects contours. Furthermore I / ∈ C 1 (Ω). We correct this lack of regularity using a filtering process that makes the filtered image C ∞ (and of course C 1 ).
The filtered image is supposed to be C ∞ , with compact support. So, we cannot choose G σ * I since the resulting image has no compact support. So we consider a fixed compact subset X of Ω and Y a compact subset of Ω containing X. Let us consider G X σ a C ∞ (Ω) projected of the Gaussian kernel G σ such as : So the regularized (filtered) image I σ = G X σ * I verifies I σ ∈ C ∞ c (Ω) and ∀x ∈ Ω, ∇I σ (x) = ∇(G X σ * I)(x) = ∇G X σ * I(x). Of course the support of I σ , contained in supp(G X σ ) + supp(I), is not necessarily included in Ω, but even if we must extend the frame Ω of the image, we consider that the filtered image has a compact support in Ω. More precisely : is C ∞ (Ω), with compact support in Ω and so bounded on Ω.
∇I σ ∈ L ∞ (Ω) ∩ C ∞ c (Ω) From now and in the sequel, we denote I σ by I and we make the following hypothesis on I: (Ω) Note that I may be extended by 0 to Ω. Now we can give f regularity properties : Proof -It is clear since I is C ∞ (Ω), with compact support in Ω and h is obvioulsy C ∞ . In addition f is nonnegative and bounded by 1 on Ω. ✷ Let us denote g the function We have seen in Proposition 3.2 that g = 1 − f ∈ C ∞ (Ω).
In the sequel, we need the following lemma : Lemma 3.1 The function g belongs to C 1 (Ω) : it is Lipschitz continuous on Ω with constant K 1 .

Propagation equation study
Equation (3.15) has been obtained quite formally with the level-set method. Now, we give existence and uniqueness of a viscosity solution ([1] for example) . Of course, we assume that the image I verifies (H I ) so that f satisfies (H 1 ). The outward unit normal ν is a C 1,1 Lipschitz vector field.

Viscosity theory framework
Let us precise the notations: for p = p 1 p 2 ∈ R 2 , the matrix 1 |p| 2 We choose β ≥ 0 and assume for simplicity that the curvature term κ is the standard mean curvature : where I is the identity matrix. Setting and F (x, p, X) := The Hamiltonian F : R 2 × R 2 × S 2 → R, is independent of t and Φ. Here S N is N × N the symmetric matrices space endowed with th classical order. First, we have to verify we may use the viscosity theory framework.
Definition 4.1 Let F be a function : F is proper if : -it is a non-decreasing function with respect to the first variable t : -it is a decreasing function with respect to the last variable X where σ i stands for the ith column of σ. Let X, Y ∈ S 2 be given. Let us assume Y ≥ X, then : The constant β is positive and the function g as well, so we deduce that : We conclude that the function F is proper. ✷ The general viscosity theory framework is thus well posed.

Ishii and Sato theorem
In what follows we use a theorem by Ishii and Sato [6] that gives existence of viscosity solution to singular degenerate parabolic (Hamilton-Jacobi) equations with nonlinear oblique derivative boundary conditions.
In this subsection we recall this theorem. In the sequel we denote ρ the function defined from R 2 in R by ρ(p, q) = min |p − q| min(|p|, |q|) , 1 .
Here ν(z) is the unit outer normal vector of Ω at z ∈ ∂Ω.

Existence and uniqueness of the solution of the evolution problem
In the sequel we assume that the balloon force H verifies (H 2 ) : H is Lipschitz continuous on Ω.
In order to use Ishii-Sato theorem we have to verify every hypothesis. The Hamiltonian F(t, x, r, p, X) = F (x, p, X) is defined by (4.23): Let us define the symmetric, semi-definite positive matrix A(x, p) : In this case F does not depend neither on t nor λ.
We choose a Neumann-type boundary condition : In addition we assumed that the balloon function H is continuous on Ω . Therefore, the Hamiltonian F is continuous on Ω × R 2 \ {0 R 2 } × S 2 [H2] Let us show there exists a constant γ ∈ R such that for each (x, p, X) ∈ ×Ω × (R 2 \ {0 R 2 }) × S 2 , the function λ → F (x, p, X) − γλ is non-decreasing on R. Since F does not explicitly depend on λ, any negative constant γ is suitable.
For the three last hypothesis B1, B2 and B3 on the boundary condition, the proof is the same as in C. Le Guyader [8].
[B1] This hypothesis consists in showing that the function B which defines the boundary condition on ]0, +∞[×∂Ω is C(R N × R N ) ∩ C 1,1 (R N × (R N \ {0 R N })). We have chosen a Neumann-type condition, by denoting ν(x) the outward unit normal vector of ∂Ω at the point x, our boundary condition is written : The point [B1] will be satisfied if ν is a C 1,1 vector field which is the case if Ω has a C 2 boundary.

Conclusion
We have recalled Gradient Vector Flow model that we hahe justified and we have given regularity properties. Then we proved existence and uniqueness fo viscosity solutions of the Hamilton-Jacobi equation derived fron the GVF-geodesic active contour model.
Next step is to perform the numerical realization of this GVF-geodesic active contour process solving he Hamilton-Jacobi equation (4.35). We shall use his method to the perform "tuffeau" tomographic images segmentation. This material is has been used during past centuries to build monuments as castles and churches in the Val-de-Loire area. These images allow to get information on the structure of the damaged material : we have to identify different phases as calcite (light grey), silice (dark grey) and porosity (black).
We shall combine the GVF-geodesic model with a region segmentaion approach to identify the three constituents of tuffeau . The segmentation of these images is a step of pretreatment which aims at reconstructing the stone porosity domain.