Christian Gout
email: christiangout<chrisgout@cal.berkeley.edu>

Carole Le Guyader
email: carole.le-guyader@insa-rennes.fr

Segmentation of complex geophysical structures with well data

Keywords: Geodesic Active Contour, Level Set Method, viscosity solution, interpolation of points. AMS classification: 49L25, 74G65, 68U10

Segmentation of complex geophysical structures with well data

Christian Gout, Carole Le

1I n t r o d u c t i o n

From geophysical data, the general purpose of this paper is to design segmentation models integrating geometrical constraints while satisfying the classical criteria of detection, with the regularity needed on the contour (see Le Guyader [START_REF] Guyader | Imagerie Mathématique: Segmentation sous contraintes géométriques[END_REF] for more details). The segmentation process has to locate layers and/or faults, but it also needs to approximate interpolation conditions which correspond to well data (see Figures 1 and2).

Figure 1. In geophysics (oil reservoir), it is necessary to determine layers accurately for post-processing, moreover, some interpolation conditions can be integrated in the model (well data). (Image: USGS, www.usgs.gov/)

Figure 2. A 2D view of a velocity seismic model (without fault). The white points correspond to interpolation conditions. The horizontal segment is the initial guess. At the end of the segmentation process, one wants to both segment the image (i.e. find the interface between 2 regions) and interpolate the well data (curve).

This kind of problems is of crucial interest in Geosciences as it consitutes an important exploration gap for reservoir characterization. Morevoer, it is well known that interpretation of faults in seismic data is today a time consuming manual task... and reducing time from exploration to production of an oil field has great economical benefits. In the exploration phase, the geological interpretation of seismic data is one of the most time consuming tasks. This is usually done manually by interpreters, and much time could be saved by an automatization of these tasks using a segmentation process (see for instance Monsen et al. [START_REF] Monsen | Geological Model Building: A Hierarchical Segmentation Approach in Mathematical Methods and Modelling in Hydrocarbon Exploration and Production[END_REF], Randen et al. [START_REF] Randen | Three-dimensional texture attributes for seismic data analysis[END_REF]). More precisely, the fault extraction from 3D data is of key importance in reservoir characterization: detailed knowledge of the fault system may provide valuable information for production.

Let us note that in this work, we only focus on the segmentation process, and not in the reconstruction of extracted surfaces (as done in Apprato et al. [START_REF] Apprato | C k reconstruction of surfaces from partial data[END_REF], Gout [START_REF]Gout: C k surface reconstruction from surface patches[END_REF], Gout and Komatitsch [START_REF] Gout | Surface fitting of rapidly varying data using rank coding : application to geophysical surfaces[END_REF]) or inverse problems ( [START_REF] Symes | The Rice Inversion Project[END_REF])

We now introduce some classical tools used for segmentation processes. Parametric deformable models, originating from the active contour model introduced by Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF], explicitly represent the interface as parameterized contours in a Lagrangian framework. Active contour models use an energy-minimizing spline that is guided by internal and external energies in such a way that the spline is deformed by geometric shape forces and influenced by image forces. By optimizing the weights used in the internal energy and choosing the proper image forces (edges), one can use active contour models to evolve the curve toward the boundary of objects being segmented. Geodesic active contours basically consist in deforming an initial contour towards the boundary of the object to be detected. The deformation is obtained trying to minimize a functional defined so that its local minimum is obtained at the boundary of the object. This energy like functional minimization problem has led to many contributions ever since. Caselles et al. [START_REF] C A S E L L E S | G e o d e sic Active Contours[END_REF] (see also Kichenassamy et al. [23] ) have shown for instance, that by setting one of the regularization parameters to zero in the classical active contour model, one gets a problem equivalent to finding a geodesic curve in a Riemann space. The issue was then no longer seen as an energy-like minimization problem but as a curve evolution one.

With the introduction of the level set method developed by Osher and Sethian [START_REF] Osher | Fronts propagation with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], geometric deformable models allow automatic topological changes without using any other procedures. The technique of level set methods is to adopt an Eulerian approach to implicitly model the propagating interface using a level set function Φ, whose zero-level set always corresponds to the position of the interface Γ. The evolution of this propagating interface is governed by a partial differential equation in a higher dimensional space. The level set function can be constructed with high accuracy in space and time. The position of the zero-level set is evolved using a speed function that consists of a constant term and a curvature deformation in its normal direction (see [START_REF] Osher | Fronts propagation with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving interfaces in Computational Geometry[END_REF], [START_REF] Sethian | Evolution, Implementation and Application of Level Set and Fast Marching Methods for Advancing Fronts[END_REF] and [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] for more details). Image stopping terms, such as regional and gradient forces, are incorporated into this speed function for segmentation purposes. This representation is intrinsic, parameter free, topology free and allows for splits and merges.

The idea that governs the following model consists in seeking a geodesic curve in a Riemann space whose metric depends both on the image contents and on the geometrical constraints.

Constraints will be mathematically defined by a set of distinct points belonging to the edge and the goal will be defined as follows: "Localize the contour Γ at the points of maximum gradient while fitting the model on the given data points". This work will be done in the context of the Level Set method

In a first step (part 2), we define the shape optimization problem and give its level set formulation. The Euler-Lagrange theorem enables us to determine the PDE that the 3D function Φ must satisfy if it minimizes the above functional. The evolution equation is then given. We establish the existence and uniqueness of the viscosity solution of the associated parabolic problem.The discretization is given in Section 3 and numerical examples illustrate the matter in Section 4. Conclusions are given in Section 5.

Mathematical formulation 2.1 Model

Let I : Ω -→ < be a given bounded image function (for instance, I can be the attribut "seismic wave propagation velocity", or just the grey level pixel values on other applications...), with Ω an open bounded subset of < n . Let us consider n =2f o rt h e purpose of illustration. As mentioned in the introduction, we plan to introduce a geometrical approach in this new method by adding interpolation constraints. Thus let S =( x i ,y i ) i ∈ Ω be a finite set of given points which correspond the the interpolation condition (well data). We would like to find a contour Γ ⊂ Ω such that Γ is the boundary of the object under consideration, interpolating the points from S w h i c hb e l o n gt ot h i s boundary. Let g :[ 0 , +∞[-→ [0, +∞[ be an edge-function as in ( [START_REF] Caselles | A geometric model for active contours in image processing[END_REF], [START_REF] Chan | An Efficient Variational Multiphase Motion for the Mumford-Shah Segmentation Mode[END_REF]), such that g(0) = 1, g is positive, strictly decreasing and lim

By definition, d(a)=0ifandonlyifa ∈ S.I no r d e rt ofind a contour Γ such that g ' 0 or d ' 0o nΓ, we propose to minimize the following energy:

E(Γ)= Z Γ d.g(|∇I|)ds. (3)
We will start with an initial guess Γ 0 and we will apply gradient descent to the energy, in a Level Set approach. We will construct a family of curves Γ(t) decreasing the energy as t increases.

Remark 1 In the previous definition of the distance function d,i t se ffect is global, so, this model could be considered not appropiated because the evolution of points which are very far from the set of points S are going to be (in some cases) influenced by S because of the distance function d given in [START_REF] Apprato | C k reconstruction of surfaces from partial data[END_REF]. So, to overcome this fact and improve the model, it is possible to use the following distance function: at the point a, we consider the following distance function d a (x) defin e do na ni n t e r v a l[A, B] ,

d a (x)=1-D a (x)e -x 2 σ 2
where

D a (x)= ( 0 sur ]-∞, α[ ∪ ]α, +∞[ P a 5 sur [-α, α]
,

with P a 5 (x) ∈ P 5 (polynomial space, degree ≤ 5) and satisfying

     P a 5 (a)=1, (P a 5 ) 0 =0, P a 5 (a + α)=(P a 5 ) 0 (a + α)=0 P a 5 (a -α)=(P a 5 ) 0 (a -α)=0 so P a 5 (a)=(cx + b)(x -α) 2 (x + α) 2 ,
the parameter α allowing to control the support of the function

D a (x) (b, c ∈ <). We note that d a (x) ∈ C 1 .
In 2D , it is easy to get the following distance function:

d(x)= n Y k=1 Ã 1 -P a k,1 5
(x 1 )P a k,2 5
(x 2 )e -[ (x 1 -a 1 ) 2 +(x 2 -a 2 ) 2 ] σ 2 ! .

The effect of the function d is therefore local: in a neighborhood of the points from S, d ≈ 0, w h i c hm e a n st h a ti tw i lls t o pt h ee v o l u t i o no ft h ec u r v e ,w h i l ei tw i lln o ta ffect the motion

of the curve in remote points from S. ¥

The Level Set Approach

The Level Set approach ( [START_REF] Osher | Fronts propagation with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving interfaces in Computational Geometry[END_REF], [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF], [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF], [START_REF] Sethian | Evolution, Implementation and Application of Level Set and Fast Marching Methods for Advancing Fronts[END_REF]) consists in considering the evolving active contour Γ = Γ(t) as the zero level set of a function Φ, which is a Lipschitz continuous function defined by: (

Φ : Ω × [0, +∞[-→ < (x, y, t) 7 -→ Φ(x, y, t) (4)
such that

Γ(t)={(x, y) ∈ Ω : Φ(x, y, t)=0},
and Φ(•, •,t) takes opposite signs on each side of Γ(t). It enables us to re-write the energy in terms of Φ as follows

F (Φ)= R Ω d(x, y)g(|∇I(x, y)|)|∇H(Φ(x, y))|dxdy,
where H is the one-dimensional Heaviside function. By approximating H by a C 1 or C 2 regularization H ² ,a s² -→ 0 and letting δ ² = H 0 ² , the energy can be written as (see Chan and Vese [START_REF] Chan | Active Contours without edges[END_REF])

F ² (Φ)= R Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ)|∇Φ(x, y)|dxdy, (5)
where

Z Ω δ ² (Φ)|∇Φ(x, y)|dxdy
is an approximation of the length term.

Minimization of the energy

In this section, we minimize the energy F ² and we determine the associated partial differential equation satisfied by Φ. To this end, it is possible to use variational calculus with the extension of the classical Euler-Lagrange theorem to a function that depends on two variables (see le Guyader [START_REF] Guyader | Imagerie Mathématique: Segmentation sous contraintes géométriques[END_REF]). But we prefer proposing a proof based on the Gâteaux derivative which permits to get boundary conditions. We recall that F is differentiable in the Gâteaux sense at Φ ∈ X if the application

x 7 -→ F 0 Φ (Ψ)= lim h→0 F (Φ + hΨ) -F (Φ) h
is defined for any Ψ ∈ X and if it is linear and continuous. In this case, Riesz theorem provides the existence of F 0 (Φ) ∈ X such that

F 0 Φ (Ψ)=<F 0 (Φ), Ψ >,
F 0 (Φ) being the gradient of F at Φ. Coming back to the problem, let us determine the Gâteaux derivative of the energy F ² .T h eG âteaux derivative of F ² with respect to Φ in the direction of Ψ is

F 0 ² Φ (Ψ)= lim h→0 F ² (Φ + hΨ) -F ² (Φ) h
where

F ² (Φ + hΨ) -F ² (Φ)= Z Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ + hΨ)|∇Φ + h∇Ψ|dxdy - Z Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ)|∇Φ|dxdy.
Then

F ² (Φ + hΨ) -F ² (Φ) = Z Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ + hΨ)|∇Φ| q 1+h 2 |∇Ψ| 2 |∇Φ| 2 +2h <∇Φ,∇Ψ> |∇Φ| 2 dxdy - Z Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ)|∇Φ|dxdy
and we use a Taylor development (on h) to linearize the square root. Taking the limit when h → 0, we get

F 0 ² Φ (Ψ)= Z Ω d(x, y)g(|∇I(x, y)|)δ 0 ² (Φ)Ψ|∇Φ|dxdy + Z Ω d(x, y)g(|∇I(x, y)|)δ ² (Φ) < ∇Φ, ∇Ψ > |∇Φ| dxdy.
Hence

F 0 ² Φ (Ψ)= - Z Ω δ ² (Φ)div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| )Ψdxdy + Z ∂Ω d(x, y)g(|∇I(x, y)|) δ ² (Φ) |∇Φ| Ψ( ∂Φ ∂x ν x dσ + ∂Φ ∂y ν y dσ)dσ and
F 0 ² Φ (Ψ)= - Z Ω δ ² (Φ)div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| )Ψdxdy + Z ∂Ω d(x, y)g(|∇I(x, y)|) δ ² (Φ) |∇Φ| ∂Φ ∂ν Ψdσ.
This expression must vanish for all Ψ in order to satisfy the Euler-Lagrange equation. Therefore, we obtain the following problem:

       δ ² (Φ)div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| )=0.
with the boundary conditions

δ ² (Φ) |∇Φ| ∂Φ ∂ν =0. (6)

Evolution equation

From the previous section, we get the evolution equation

∂Φ ∂t = δ ² (Φ)div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| ). ( 7
)
As stressed by Zhao et al. [START_REF] Zhao | Implicit and non parametric shape reconstruction from unorganized data using a variational level set method[END_REF], there is a balance between the potential force and the surface tension. A parallel can be drawn with the classical deformable models, a model which shows off an equilibrium between the regularization energy and the energy linked to the image. The closer we are to the finite set of points or on edges, the more important is the flexibility in the model, since in this case the expression d(x, y)g(|∇I(x, y)|)v a n i s h e s .

Proposition 1 The energy F ² (Φ) is decreasing with time t.

Proof: We follow the same arguments as Zhao et al. [START_REF] Zhao | A variational level set approach to multiphase motion[END_REF].¥ When a local minimum is reached, then the quantity ∂Φ ∂t tends to 0, which means that the steady state is reached. A rescaling can be made so that the motion is applied to all level sets by replacing δ ² by |∇Φ|. As stressed by Zhao et al. [START_REF] Zhao | Implicit and non parametric shape reconstruction from unorganized data using a variational level set method[END_REF] and Alvarez et al. [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], it makes the flow independent of the scaling of Φ. Thus the proposed parabolic problem with the associated boundary conditions ∂Φ ∂ν = 0 can be written:

         Φ(x, 0) = Φ 0 (x), ∂Φ ∂t = |∇Φ| • div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| ) ¸, ∂Φ ∂ν =0 on ∂Ω. (8)
The evolution equation satisfied by Φ is defined by

∂Φ ∂t = |∇Φ| • div(d(x, y)g(|∇I(x, y)|) ∇Φ |∇Φ| ) ¸.
A l lt h el e v e ls e tm o v ea c c o r d i n gt o

Γ t = d(x, y)g(|∇I(x, y)|)κñ-< ∇(d(x, y)g(|∇I(x, y)|)),ñ>ñ, (9)
with κ = div( ∇Φ |∇Φ| )a n dñ = -∇Φ |∇Φ| (interior normal). Indeed, we have for all level sets Φ(Γ(t),t)=constant.

Calculating the derivatives with respect to the variable t, we have

d dt [Φ(Γ(t),t)] = 0 ⇔ Φ t + < ∇Φ, Γ t >=0,
where Φ t denotes ∂Φ ∂t and Γ t , ∂Γ ∂t . Hence, using the definition of Γ t and ñ,w eg e t

Φ t = d(x, y)g(|∇I(x, y)|)κ|∇Φ|+ < ∇(d(x, y)g(|∇I(x, y)|)), ∇Φ >.
A parallel can be drawn with the Caselles et al. model ( [START_REF] Caselles | Minimal surfaces: a geometric threedimensional segmentation approach[END_REF], [START_REF] C A S E L L E S | G e o d e sic Active Contours[END_REF]) in which the evolution equation of any level set is given by:

Γ t = g(|∇I(x, y)|)κñ-< ∇(g(|∇I(x, y)|)),ñ>ñ.
The convergence of the model can be increased by adding the component kdg(|∇I|),ñ in the evolution equation of a level set [START_REF] Caselles | Minimal surfaces: a geometric threedimensional segmentation approach[END_REF], k being a constant. This component can be seen as an area constraint. An analogy with the Balloon model developed by Cohen [START_REF] Cohen | On active contours models and balloons[END_REF] can be made: this constant motion term prevents the curve to stop on a non significative local minimum and is also of importance when starting from curves inside the object to be detected.

The associated parabolic problem becomes:

               Φ(x, 0) = Φ 0 (x), ∂Φ ∂t = |∇Φ|d(x)g(|∇I(x)|)div( ∇Φ |∇Φ| )+ < ∇(d(x)g(|∇I(x)|)), ∇Φ > +kd(x, y)g(|∇I(x, y)|)|∇Φ|, ∂Φ ∂ν =0 o n ∂Ω. (10)
In the next section, we aim at proving the existence and uniqueness of the solution of this parabolic problem which are the main results of this paper.

Remark 2 In the evolution equation satisfied by Φ, ∇d is well-defined, except at the points that are equidistant from at least two points of the given finite set S and at the points from S. Indeed, the function d(x)=d(x; S) is continuous as the inf of a finite number of continuous functions but is not differentiable at the points equidistant from at least two points of the given finite set S and at the points from S. The distance function d satisfies the Eikonal equation |∇d| =1. In the theoretical part devoted to the existence and uniqueness of the solution of our problem, we need a certain smoothness on the distance function d. Using the curvature as a regularizing or smoothing term enables us to get the desired properties on d. What follows is taken from Sethian's book [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving interfaces in Computational Geometry[END_REF]: the main conclusion that we use here is that "a front propagating at the speed 1 -²κ for ²>0 does not form corners and remains smooth for all time". Furthermore, as the dependence on curvature vanishes, the limit of this motion is the entropy-satisfying solution obtained for the constant speed case.

Existence, uniqueness of the solution

In this section, we use the notion of viscosity solutions (for more details, see the 'User's guide to viscosity solutions' by Crandall, Ishii and Lions [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF]). This theory applies to some partial differential equations that can formally be written in the form F (x, u, Du, D 2 u) where Du denotes the gradient and D 2 u the Hessian matrix (which is symmetric). In general, F is defined as

F : < n ×<×< n × S(n) -→ <
where S(n) denotes the set of symmetric (n × n) matrices. In our case, F will be defined by F : Ω ×<×< 2 × S(2) -→ <. A complete study of viscosity solutions for this problem is studied in le Guyader [START_REF] Guyader | Imagerie Mathématique: Segmentation sous contraintes géométriques[END_REF], where different approaches of viscosity solutions are given for our considered problem (based on works of Alvarez et al. [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], Barles [5] and [START_REF] Barles | Solutions de viscositéd e séquations de Hamilton-Jacobi[END_REF], Caselles et al. [START_REF] Caselles | A geometric model for active contours in image processing[END_REF], Chen et al. [START_REF] Chen | Uniqueness and existence of ciscosity solutions of generalized mean curvature flow equations[END_REF], Crandall and Lions [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF]).

Here is a first result that will be used in the sequel. Preliminary: Let p, q ∈ < n \{0} , then

| p |p| - q |q| | ≤ |p -q| min(|p|, |q|)
.

Proof: evident.

We use here the existence theorem for viscosity solutions introduced by Ishii and Sato in [START_REF] Ishii | Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain[END_REF]. This article treats the difficult case of singular parabolic equations with nonlinear oblique derivative boundary conditions while we wish to apply it to a problem with homogeneous Neumann boundary conditions. As formerly done in their article, we denote by ρ(p, q)=min µ |p -q| min(|p|, |q|) , 1 ¶ . We assume that Ω is a bounded domain in < n with a C 1 boundary. Let us consider the following conditions:

1. F ∈ C([0,T] × Ω ×<×(< n -{0}) × S n )
, where S n denotes the space of n × n symmetric matrices equipped with the usual ordering.

2. There exists a constant γ ∈ < such that for each

(t, x, p, X) ∈ [0,T] × Ω × (< n -{0}) × S n ,
the function u 7 → F (t, x, u, p, X)γu is non decreasing on <.

3. For each R>0, there exists a continuous function w

R :[ 0 , ∞[-→ [0, ∞[s a t i s f y i n g w R (0) = 0 such that if X, Y ∈ S n and µ 1 ,µ 2 ∈ [0, ∞[s a t i s f y : Ã X 0 0 Y ! ≤ µ 1 Ã I -I -II ! + µ 2 Ã I 0 0 I ! (11)
then

F (t, x, u, p, X) -F (t, y, u, q, -Y ) ≥-w R (µ 1 (|x -y| 2 + ρ(p, q) 2 )+µ 2 + |p -q| + |x -y|(1 + max(|p|, |q|))),
for all t ∈ [0,T] ,x,y ∈ Ω,u ∈ < with |u| ≤ R and p, q ∈ < n \{0} .

4. B ∈ C(< n ×< n ) T C 1,1 (< n × (< n \{0}))
5. For each x ∈ < n , the function p 7 → B(x, p) is positively homogeneous of degree one in p, i.e., B(x, λp)=λB(x, p), ∀λ ≥ 0,p ∈ < n \{0}.

6. There exists a positive constant θ such that < ν(z),D p B(z, p) >≥ θ for all z ∈ ∂Ω and p ∈ < n -{0}.H e r eν(z) denotes the unit outer normal vector of Ω at z ∈ ∂Ω.

We recall the following theorem taken from Ishii and Sato [START_REF] Ishii | Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain[END_REF].

Theorem 2 Consider the following problem :

( u t + F (t, x, u, Du, D 2 u)=0in ]0,T[×Ω B(x, Du)=0in ]0,T[×∂Ω ( 12
)
satisfying u(0,x)=g(x) for x ∈ Ω. Assume that conditions 1, 2, 3, 4, 5, 6 hold. Then for each g ∈ C( Ω) there is a unique viscosity solution u ∈ C([0,T[× Ω) of ( 10) satisfying u(0,x)=g(x) for x ∈ Ω.¥

We apply this theorem to the considered problem. The function F is defined by:

F (t, x, u, p, X)= -trace µ d(x)g(|∇I(x)|)(I - p N p |p| 2 )X ¶ -< ∇(d(x)g(|∇I(x)|)),p > .
Denoting by A(x, p) the symmetric positive matrix defined by

A(x, p)=d(x)g(|∇I(x)|)(I - p N p |p| 2 ), we get F (t, x, u, p, X)=-trace(A(x, p)X)-< ∇(d(x)g(|∇I(x)|)),p > .
F presents a singularity for p = 0 but is continuous otherwise. The first point is satisfied.

F does not depend explicitly on u so any negative constant γ satisfies the second condition.

For the third point, the inequality [START_REF] Chan | Active Contours without edges[END_REF] gives us that for all r, s ∈ < 2 , (Xr,r)+(Ys,s)

≤ µ 1 |r -s| 2 + µ 2 (|r| 2 + |s| 2 ).
Taking successively r = σ(x, p)e i and s = σ(y, q)e i with (e i ) i an orthonormal basis of < 2 (as done in [START_REF] Gout | Image Segmentation under Interpolation conditions[END_REF], A(x, p)=σ(x, p)σ T (x, p)), we get:

trace(A(x, p)X)+trace(A(y, q)Y ) ≤ µ 1 trace((σ(x, p) -σ(y, q))(σ(x, p) -σ(y, q)) T ) +µ 2 (d(x)g(|∇I(x)|)+d(y)g(|∇I(y)|)).
Hence

trace(A(x, p)X)+trace(A(y, q)Y ) ≤ µ 1 | p d(x)g(|∇I(x)|) p |p| - p d(y)g(|∇I(y)|) q |q| | 2 +2θµ 2 .
The function

x 7 → d(x)g(|∇I(x)|) is bounded by θ.
Moreover, one has:

| p d(x)g(|∇I(x)|) p |p| - p d(y)g(|∇I(y)|) q |q| | 2 = |( p d(x)g(|∇I(x)|) - p d(y)g(|∇I(y)|)) p |p| + p d(y)g(|∇I(y)|)( p |p| -q |q| )| 2 . Thus | p d(x)g(|∇I(x)|) p |p| - p d(y)g(|∇I(y)|) q |q| | 2 ≤ 2( p d(x)g(|∇I(x)|) - p d(y)g(|∇I(y)|)) 2 +2d(y)g(|∇I(y)|)| p |p| -q |q| | 2 .
Using the preliminary and properties of the functions

x 7 → d(x)g(|∇I(x)|)a n dx 7 → p d(x)g(|∇I(x)|)
as in [START_REF] Gout | Image Segmentation under Interpolation conditions[END_REF], we can conclude that

trace(A(x, p)X)+trace(A(y, q)Y ) ≤ µ 1 (2ζ|x -y| 2 +8θρ(p, q) 2 )+2θµ 2 .
N e x t ,w eh a v et oe v a l u a t et h ee x p r e s s i o n

F (t, x, u, p, X) -F (t, y, u, q, -Y )= -(trace(A(x, p)X)+trace(A(y, q)Y )) -(< ∇(d(x)g(|∇I(x)|)),p > -< ∇(d(y)g(|∇I(y)|)),q >).
Using the same arguments as in [START_REF] Gout | Image Segmentation under Interpolation conditions[END_REF], we have

| < ∇(d(x)g(|∇I(x)|)),p > -< ∇(d(y)g(|∇I(y)|)),q > | ≤ κ|x -y|max(|p|, |q|)+C 2 |p -q|,
from which we deduce

-(F (t, x, u, p, X) -F (t, y, u, q, -Y )) ≤ µ 1 [2ζ|x -y| 2 +8θρ(p, q) 2 ] +2θµ 2 + κ|x -y|max(|p|, |q|)+C 2 |p -q| and -(F (t, x, u, p, X) -F (t, y, u, q, -Y )) ≤ max(2ζ, 8θ,C 2 , κ)(µ 1 (ρ(p, q) 2 + |x -y| 2 ) +µ 2 + |p -q| + |x -y|(1 + max(|p|, |q|))).
We just have to take w R (l)=max(2ζ, 8θ,C 2 , κ)l. w R (0) = 0 and w R is non-decreasing on [0, ∞[. The fourth point is fulfilled with assumptions on ν (vector field of class C 1,1 ). Then, it is easy to check that B is positively homogeneous of degree one. For the last point, one can easily see that:

B(z, p)=< ν(z),p > (13)
and

< ν(z),D p B(z, p) >= |ν(z)| 2 =1.
We take θ = 1 and the last assumption is fulfilled.

Numerical algorithm

We recall the evolution equation

               ∂Φ ∂t = |∇Φ|d(x)g(|∇I(x)|)div( ∇Φ |∇Φ| ) + < ∇(d(x)g(|∇I(x)|)), ∇Φ > +kd(x)g(|∇I(x)|)|∇Φ| ,o n[ 0 , +∞[×Ω., Φ(0,x)=Φ 0 (x), ∂Φ ∂ν =0on∂Ω., (14)
Its discretization has the following characteristics:

• The distance function d is computed using the Fast Marching method (see Sethian [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving interfaces in Computational Geometry[END_REF], Hvistendahl Karlsen et al. [START_REF] Hvistendahl Karlsen | A fast marching method for reservoir simulation[END_REF]).

• We have chosen the Additive Operator Splitting (AOS, see Weickert and Kühne [START_REF] Weickert | Fast methods for implicit active contours models[END_REF]) scheme.

The classical spatial discretization of the term div(c|∇Φ|

)i s        div(c∇Φ) ' ∂ x (c ij Φ i+ 1 2 ,j -Φ i-1 2 ,j h )+∂ y (c ij Φ i,j+ 1 2 -Φ i,j-1 2 h
),

' c i+ 1 2 ,j Φ i+1,j -Φ i,j h 2 -c i-1 2 ,j Φ i,j -Φ i-1,j h 2 + c i,j+ 1 2 Φ i,j+1 -Φ i,j h 2 -c i,j-1 2 Φ i,j -Φ i,j-1 h 2 .
The terms

c i + - 1 2 ,j and c i,j + -1 2
can be determined by linear interpolation. To simplify the notation in the following, we use a vectorial representation of the function Φ via a concatenation of the rows of the image. So, Φ ∈ R N×M where N is the number of lines and M the one of columns. The center of gravity of a pixel i is associated w i t ht h en o d eo ft h em e s h i n go fc o o r d i n a t e sx i . Thus, the term Φ n i corresponds to an approximation of Φ(x i ,t n ). We first suppose that k =0 . The discretization of the evolution equation 14 is given by the semi-implicit scheme

Φ n+1 i = Φ n i + τ |∇Φ| n i X j∈Λ(i) ( dg |∇Φ| ) n i +( dg |∇Φ| ) n j 2 . Φ n+1 j -Φ n+1 i h 2 ,
where Λ(i) denotes the neighbourhood of the pixel i. However, to assure that |∇Φ| does not vanish, one can replace the arithmetic mean

( dg |∇Φ| ) n i +( dg |∇Φ| ) n j 2
by its harmonic counterpart. If |∇Φ| n i =0or(dg) i = 0, we set:

Φ n+1 i = Φ n i .
Then, we introduce the linear system

Φ n+1 = Φ n + τ X l∈{x,y} A l (Φ n )Φ n+1 .
where A x (Φ n )i sd e fined by

a ijx (Φ n )=          |∇Φ| n i 2 ( |∇Φ| dg ) n i +( |∇Φ| dg ) n j ,j∈ Λ x (i), -|∇Φ| n i P m∈Λx(i) 2 ( |∇Φ| dg ) n i +( |∇Φ| dg ) n m ,i= j, 0e l s e ,
where Λ x (i) represents the neighbouring pixels of i with respect to direction x. However, the solution Φ n+1 cannot be directly determined from this scheme, it requires to solve the following linear system: 

 Id -τ X l∈{x,y} A l (Φ n )   Φ n+1 = Φ n ,
where Id istheunitmatrix. Thematrix

à Id -τ P l∈{x,y} A l (Φ n )
! is strictly diagonally dominant, it follows from Gershgorin's theorem that it is invertible. To solve this linear system, iterative algorithms should be used. Classical methods like Gauss-Seidel does not need additional storage, and convergence can be guaranteed for the special structure of the system matrix. This convergence, however, may be rather slow (when parameter τ increases) since the condition number of the system matrix increases with the image resolution. An alternative discretization is proposed by Weickert and K ünhe [START_REF] Weickert | Fast methods for implicit active contours models[END_REF], they focus on a splitting-based alternative (AOS scheme).

It is easy to implement and does not require to specify any additional parameters. This may make it attractive in a number of practical applications. Instead of using the semi-implicit scheme

Φ n+1 =   Id -τ X l∈{x,y} A l (Φ n )   -1 Φ n ,
Weickert an K ühne [START_REF] Weickert | Fast methods for implicit active contours models[END_REF] consider the additive operator splitting (AOS) variant:

Φ n+1 = 1 2 X l∈{x,y} (Id -2τ A l (Φ n )) -1 Φ n . (15)
By means of a Taylor expansion, it is easy to see that the semi-implicit scheme and its AOS version differ by an O(τ 2 )t e r m . T h eA O Ss c h e m el e a d st o

2Φ n+1 =(Id -2τ A x (Φ n )) -1 Φ n +(Id -2τ A y (Φ n )) -1 Φ n . (16)
We set for l ∈ {x, y},

B l (Φ n )=Id -2τ A l (Φ n ).
The matrix B is strictly diagonally dominant tridiagonal linear system which can be solved very efficiently with a Gaussian algorithm. As shown in [START_REF] Weickert | Fast methods for implicit active contours models[END_REF], the fact that AOS schemes are based on an additive splitting guarantees that both axes are treated in exactly the same manner. AOS schemes are not only efficient, they are also unconditionally stable. ¥

Numerical results

• The user defines an initial condition (paraboloids, ellipses...).

• The user gives interpolation conditions (if there are ones).

• The stopping criteria is either the number of iteration, either a verification that the solution is stationnary while computing the relative error between Φ n+1 and Φ n .

Example 1 (velocity seismic model ). We give a numerical example on a velocity seismic model. There is a vertical fault and several layers. The time step is equal to 2, k =0 .02, and g : s 7 → 1 1+s 2 . CPU time is equal to 96 seconds ( p cI n t e lP e n t i u m42 . 8 G h z 256Mo). Example 2 (2D fault extraction). The time step is equal to 2, k =0 .02, and g : s 7 → 1 1+s 2 . CPU time is equal to 128 seconds (pc Intel Pentium 4 2.8Ghz 256Mo).T h i sa very difficult (extreme!) example because we want to locate a vertical fault which is very difficult to segment. In this example, the interpolation condition are not well data, but points given by the user in order to help the segmentation process. Without these (well chosen) interpolation condition, the segmentation is impossible to realize.

5C o n c l u d i n g r e m a r k s

We have developed a segmentation method under interpolation condition which gives good results.

But some works still remain to be done:

• The CPU time is suspiciously high for such a result, the main problem is that the fast marching method to approximate the distance function needs 80% of the CPU time. We are currently trying to improve this part of the algorithm.

• Of course, the final goal will be to segment more complex geophysical data. For now, to segment a 3D block (see Figure 7), it is needed to start the process with an initial guess, and the main problem is to be able to give a "correct" initial condition. And so to have a "good" visualization of the data set. Another approach would be to segment 2D slices (like example 2 of section 4) of the 3D block, and then apply a 3D reconstruction method (work in progress). We conclude with some remarks:

Remark 3 This method has been successfully tested on classical image segmentation without interpolation condition in Le Guyader [START_REF] Guyader | Imagerie Mathématique: Segmentation sous contraintes géométriques[END_REF].

Remark 4 Some comparisons with other methods (deformable models method under geometrical constraints developed by Le Guyader et al. [START_REF] Guyader | Using a Level Set Approach for Image segmentation under interpolation conditions[END_REF]) are given in [START_REF] Guyader | Imagerie Mathématique: Segmentation sous contraintes géométriques[END_REF].

Remark 5 This method has been validated on medical applications (see Apprato et al. [START_REF] Apprato | Segmentation of Medical Images Sequence for Non Invasive Assessment of Pulmonary Arterial Hypertension[END_REF]): in this case, the interpolation conditions are useful in order to help the process when some image data are missing. The goal was to compute the blood flow variations, throughout a cardiac cycle. With our colleagues from "C.H.U du Haut-Lévèque-Hôpital de Bordeaux", this method has been successfully applied on an patient image sequence to outline the cross-sectional area of a great thoracic vessel, namely the main pulmonary artery, in order to non invasively assess pulmonary arterial hypertension.

Remark 6 Let us note that in case of noisy images, it is possible to use usual filter or specific algorithm (see Sethian [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving interfaces in Computational Geometry[END_REF]).

  s→∞ g(s) = 0. The function g is applied to the gradient of the image |∇I(x, y)|. An example of such a function is given by gthe set of points S, we associate the distance function d(a) from every point a ∈ Ω to S d(a)= distance(a, S)=min c∈S |a -c| .

Figure 3 .

 3 Figure 3. 2D view of the initial condition. There are 6 interpolation conditions.

Figure 4 .

 4 Figure 4. The final contour after 140 iterations.

Figure 5 .

 5 Figure 5. Fault extraction. Many interpolation conditions (9 points) are needed on such an image, in order to get a vertical fault. User intervention is crucial in this process. (source data: confidential). The algorithm has been modified in order to define such an initial condition.

Figure 6 .

 6 Figure 6. Final result.

Figure 7 .

 7 Figure 7. Example of 3D geological data (Location: Pyrénées region near Pau, France) where several layers (A and B) are clearly located, but to segment the complex structure, one needs to give initial guesses and a scientific visualization tool must be used to look inside the data.

Acknowledgements: The authors are grateful to Luminita Vese (UC Los Angeles) and to Dominique Apprato (Université de Pau) for valuable comments. This work is related to the MAGIQUE 3D project (INRIA Futurs, Pau, France). A large part of this work has been done at the "Laboratoire de Mathématiques Appliquées (EA 3226)" at INSA de Rouen, France.