Lucie Baudouin
email: lucie.baudouin@laas.fr

Eduardo Cerpa
email: eduardo.cerpa@usm.cl

Emmanuelle Crépeau
email: emmanuelle.crepeau@math.uvsq.fr

Alberto Mercado
email: alberto.mercado@usm.cl

On the determination of the principal coefficient from boundary measurements in a KdV equation

Keywords: Inverse problem, Korteweg-de Vries equation, Carleman estimate AMS subject classifications: 35R30, 35Q53

This paper concerns the inverse problem of retrieving the principal coefficient in a Korteweg-de Vries (KdV) equation from boundary measurements of a single solution. The Lipschitz stability of this inverse problem is obtained using a new global Carleman estimate for the linearized KdV equation. The proof is based on the Bukhgeȋm-Klibanov method.

Introduction

The Korteweg-de Vries (KdV) equation y t (t, x) + y xxx (t, x) + y x (t, x) + y(t, x)y x (t, x) = 0, is a nonlinear dispersive equation that serves as a mathematical model to study the propagation of long water waves in channels of relatively shallow depth and flat bottom [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF]. In this model, the function y = y(t, x) represents the surface elevation of the water wave at time t and at position x. Since the works by Johnson [START_REF] Johnson | On the development of a solitary wave moving over an uneven bottom[END_REF] and Grimshaw [START_REF] Grimshaw | Evolution equations for long, nonlinear internal waves in stratified shear flows[END_REF] (see also the recent survey [START_REF] Grimshaw | Solitary waves propagating over variable topography[END_REF]), the study of water waves moving over variable topography has been considered. If we denote h = h(x) the function describing the variations in depth of the channel, then the proposed model becomes (after scaling) y t (t, x) + h 2 (x)y xxx (t, x) + ( h(x)y(t, x)) x + 1 h(x) y(t, x)y x (t, x) = 0.

In a recent paper, Israwi [START_REF] Israwi | Variable depth KdV equations and generalizations to more nonlinear regimes[END_REF] proposes a correction to this model by saying that the main coefficient is not h 2 (x)y xxx (t, x) but h 5 2 (x)y xxx (t, x). Thus, we are leaded to consider variable coefficients KdV equations to model the water wave propagation in non-flat channels.

The original inverse problem which motivated this paper can be stated in a informal way as "to recover information on the topography of a channel (denoted by h) by means of boundary measurements of a water wave (given by y)". This kind of coefficient inverse problem is by now classical and referred to as the determination of coefficients in evolution partial differential equations (pde's) by means of the partial knowledge of a single solution.

Many other inverse problems on KdV equations could be studied with possible applications. For example the pulsatile blood flow in an arterial compartment can be modelized by KdV equation, see [START_REF] Crépeau | A reduced model of pulsatile flow in an arterial compartment[END_REF] for details. It leads to several open inverse problems, as the determination of the flow at the input of the arteria with only one measure at the output, or the determination of physical coefficients like stiffness of arterias. Those open inverse problems are very relevant for cardiovascular medicine.

Since the pioneer works by Bukhgeȋm, Klibanov and Malinsky [START_REF] Bukhgeȋm | Uniqueness in the large of a class of multidimensional inverse problems[END_REF][START_REF] Klibanov | Newton-Kantorovich method for threedimensional potential inverse scattering problem and stability of the hyperbolic Cauchy problem with time-dependent data[END_REF], there have been a number of papers devoted to the proof of uniqueness and stability for inverse problems concerning the determination of source or coefficients in pde's by means of the localized knowledge of a single solution.

The problem of recovering sources has been first addressed for hyperbolic equations in [START_REF] Puel | On a global estimate in a linear inverse hyperbolic problem[END_REF][START_REF] Puel | Generic well-posedness in a multidimensional hyperbolic inverse problem[END_REF] and for parabolic equations in [START_REF] Imanuvilov | Lipschitz stability in inverse parabolic problems by the Carleman estimate[END_REF]. The determination of coefficients has also drawn the attention of many authors. Among the most classical results for the wave equation, we can mention the determination of a potential from boundary [START_REF] Yamamoto | Uniqueness and stability in multidimensional hyperbolic inverse problems[END_REF] or internal [START_REF] Imanuvilov | Global Lipschitz stability in an inverse hyperbolic problem by interior observations[END_REF] measurements and the determination of the main coefficient from boundary [START_REF] Bellassoued | Uniqueness and stability in determining the speed of propagation of second-order hyperbolic equation with variable coefficients[END_REF][START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF] or internal [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF] observation. Regarding the problem of recovering coefficients in parabolic equations, we send the interested reader to the survey [START_REF] Yamamoto | Carleman estimates for parabolic equations and applications[END_REF] by Yamamoto and the references therein.

Concerning dispersive equations, a stability result was obtained in [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF] for the problem of recovering the potential in a Schrödinger equation. That paper motivated several further works as [START_REF] Cardoulis | Inverse problem for the Schrödinger operator in an unbounded strip[END_REF][START_REF] Mercado | Lionel Inverse problems for the Schrödinger equation via Carleman inequalities with degenerate weights[END_REF], where appropriate Carleman estimates were proven in different contexts. To the best of our knowledge, there are no previous works concerning coefficient inverse problems for KdV equations.

Most of the references given above are inspired by what is often called the Bukhgeȋm-Klibanov method (see the recent book [START_REF] Beilina | Approximate global convergence and adaptivity for coefficient inverse problems[END_REF] or the survey [START_REF] Klibanov | Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems[END_REF]), which was meant to prove global uniqueness results for single-measurement coefficients inverse problems. Roughly speaking, it consists in two steps. The first step comprises the time differentiation of the governing pde, getting an auxiliary system, in which the unknown coefficient appears now in the initial condition. The second step relies on the application of a Carleman estimate whose parameters and weight are crucial to drive to the conclusion. An idea of Imanuvilov and Yamamoto, presented for the first time in [START_REF] Imanuvilov | Lipschitz stability in inverse parabolic problems by the Carleman estimate[END_REF] (for parabolic equation) and then often used (for other pde's), allows to use the appropriate global Carleman inequality to estimate the unknown coefficient in terms of internal or boundary measurements of the solution. The main drawback of this method, that is not lifted yet is the hypothesis that at least one initial condition, (or some of its derivatives) never vanishes in the domain.

This article is a first step in order to address the problem of recovering the topography of the channel by means of boundary measurements of the water wave. Indeed, we will deal with a KdV equation posed on a bounded interval [0, L] with an unknown third-order coefficient:

   y t + a(x)y xxx + y x + yy x = 0, ∀(x, t) ∈ (0, L) × (0, T ), y(0, t) = y(L, t) = y x (L, t) = 0, ∀t ∈ (0, T ), y(x, 0) = y 0 (x), ∀x ∈ (0, L), (2)
where the initial data y 0 is known. The unknown coefficient a = a(x) is assumed to be time independent. The case where we have several unknown coefficients (as in the equation ( 1)) is very hard to deal with. There are few results in this direction and they ask many conditions on the coefficients. See for instance [START_REF] Cardoulis | Simultaneous identification of the diffusion coefficient and the potential for the Schrödinger operator with only one observation[END_REF] where a Schrödigner equation is considered on an unbounded strip. Up to our knowledge, this paper solves the first coefficient inverse problem for the KdV equation. As mentioned before, the Bukhgeȋm-Klibanov method requires a Carleman estimate for the equation. This type of inequalities has already been proven for the KdV equation on a bounded interval in the framework of control problems. The control of this equation from the left Dirichlet boundary condition have been studied in [START_REF] Rosier | Control of the surface of a fluid by a wavemaker[END_REF] and [START_REF] Glass | Some exact controllability results for the linear KdV equation and uniform controllability in the zero-dispersion limit[END_REF], where the authors get the null-controllability of this equation by proving one-parameter Carleman estimates. Unlike those papers, here we will prove a twoparameter estimate in order to be able to deal with a space dependent main coefficient.

The inverse problem we are interested in can be stated as follows.

Inverse problem: Retrieve the principal coefficient a = a(x) of equation ( 2) from the measurement of y x (0, t), y xx (0, t) and y xx (L, t) on (0, T ), where y is the solution of equation [START_REF] Baudouin | Global carleman estimate on a network for the wave equation and application to an inverse problem[END_REF].

A partial and local answer for this nonlinear inverse problem will be proved in this paper. To be more specific, let us denote by y and ỹ the solutions of Equation (2) corresponding to coefficients a and ã respectively (we shall also use the notation y[a], y[ã]).

Stability: Estimate ã -a L 2 (0,L) by suitable norms of the measurements ỹx (0, t)y x (0, t) , ỹxx (0, t) -y xx (0, t) and ỹxx (L, t) -y xx (L, t) .

The next theorem is the main result of this article and gives the stability result for the inverse problem under consideration.

Theorem 1 Let r 0 , a 0 , α and K be given positive constants. Assume that the initial data y 0 ∈ w ∈ H 7 (0, L); w(0) = w(L) = w (L) = 0 satisfies for every x ∈ [0, L]

y 0 (x) = y 0 (L -x) (3)
and

|y 0 (x)| ≥ r 0 > 0. ( 4
)
Define the set

Σ(a 0 , α) = a ∈ W 6,∞ (0, L) ∀x ∈ [0, L], a(x) ≥ a 0 > 0, a(x) = a(L -x),
and a W 6,∞ (0,L) ≤ α . (5)

Then, there exists a constant C = C(L, T, r 0 , a 0 , α, K) > 0 such that

C a -ã L 2 (0,L) ≤ y x (0, t) -ỹx (0, t) H 1 (0,T ) + y xx (0, t) -ỹxx (0, t) H 1 (0,T ) + y xx (L, t) -ỹxx (L, t) H 1 (0,T ) (6)
for every a, ã in Σ(a 0 , α) such that the corresponding solutions of Equation (2) satisfy max{ y W 1,∞ (0,T ;W 3,∞ (0,L)) , ỹ W 1,∞ (0,T ;W 3,∞ (0,L)) } ≤ K.

Remark 1 In the Appendix, we will sketch the proof of the fact that the required hypotheses on a ∈ W 6,∞ (0, L) and y 0 ∈ w ∈ H 7 (0, L); w(0) = w(L) = w (L) = 0 are enough to guarantee that the solution of (2) lies in W 1,∞ (0, T ; W 3,∞ (0, L)) (see Proposition 3). Thus, the set of trajectories where Theorem 1 is valid is not empty. We mention that hypothesis y 0 ∈ H 7 (0, L) is not sharp, see Remark 6 in the Appendix.

Remark 2

The symmetry hypothesis on the initial data and the main coefficient are technical points coming from our proof. More precisely, they allow us to extend the solution of KdV equation to negative times, in order to apply the Bukhgeȋm-Klibanov method.

Remark 3 Concerning the question of the minimal number of boundary measurement required in Theorem 1, one could expect to get estimate (6) with either one measurement at x = L or two measurements at x = 0, but not with just one measurement at x = 0. This is based on the fact that Carleman estimates imply some observability inequalities, which do not hold for the linear KdV equation with only one observation at x = 0. This lack of observability holds for some critical lengths L, as one can read in [START_REF] Rosier | Exact boundary controllability for the Korteweg-de Vries equation on a bounded domain[END_REF][START_REF] Coron | Exact boundary controllability of a nonlinear KdV equation with critical lengths[END_REF][START_REF] Cerpa | Exact controllability of a nonlinear Korteweg-de Vries equation on a critical spatial domain[END_REF][START_REF] Cerpa | Boundary controllability for the nonlinear Korteweg-de Vries equation on any critical domain[END_REF].

On the other hand, we need in this paper three measurements because of the symmetry assumptions and the extension to negative times used.

Other related papers on this type of inverse problem for pde's can be listed: [1, 26] (one dimensional fourth order parabolic equation), [START_REF] Benabdallah | Stability of discontinuous diffusion coefficients and initial conditions in an inverse problem for the heat equation[END_REF][START_REF] Baudouin | A global Carleman estimate in a transmission wave equation and application to a one-measurement inverse problem[END_REF][START_REF] Baudouin | An inverse problem for Schrödinger equations with discontinuous main coefficient[END_REF] (discontinuous coefficients), [START_REF] Baudouin | Global carleman estimate on a network for the wave equation and application to an inverse problem[END_REF] (network of one dimensional waves) [START_REF] Bellassoued | Logarithmic stability in determination of a coefficient in an acoustic equation by arbitrary boundary observation[END_REF][START_REF] Bellassoued | Logarithmic stability in the dynamical inverse problem for the Schrödinger equation by arbitrary boundary observation[END_REF] (logarithmic stability estimates), [START_REF] Cristofol | Identification of two coefficients with data of one component for a non linear parabolic system[END_REF] (parabolic system), [START_REF] Roques | On the determination of the nonlinearity from localized measurements in a reaction-diffusion equation[END_REF] (unknown coefficient in the nonlinearity), [START_REF] Yamamoto | Global uniqueness and stability for a class of multidimensional inverse hyperbolic problems with two unknowns[END_REF][START_REF] Cardoulis | Simultaneous identification of the diffusion coefficient and the potential for the Schrödinger operator with only one observation[END_REF][START_REF] Cristofol | Identification of two coefficients with data of one component for a non linear parabolic system[END_REF] (two unknown coefficients). We also mention some important books which can be the starting point to study inverse problems [START_REF] Isakov | Inverse problems for partial differential equations[END_REF], Carleman estimates [START_REF] Fursikov | Controllability of evolution equations[END_REF], and control theory for partial differential equations [START_REF] Coron | Control and nonlinearity[END_REF].

Outline: In section 2, a two-parameters global Carleman estimate for the linearized KdV equation with non-constant main coefficient is obtained. It is then used in section 3, following the Bukhgeȋm-Klibanov method, to prove the desired stability result. As mentioned before, the Appendix A is concerned with the Cauchy problem.

Carleman estimate

In this section, a global Carleman inequality will be proven for the linearized KdV equation on the domain denoted here by Q := (0, L) × (-T, T ). As mentioned in the introduction, there are in the literature some Carleman estimates for the KdV equation with constant main coefficient posed on a bounded interval [START_REF] Rosier | Control of the surface of a fluid by a wavemaker[END_REF][START_REF] Glass | Some exact controllability results for the linear KdV equation and uniform controllability in the zero-dispersion limit[END_REF]. In this work we deal with non-constant coefficients a = a(x). We prove a two-parameter (s, λ) Carleman estimate in order to avoid extra hypothesis on the coefficients when solving the inverse problem.

Let us consider

d ∈ L ∞ ((-T, T ) × (0, L)), b ∈ L ∞ (-T, T ; W 1,∞ (0, L)
) and a ∈ W 3,∞ (0, L), with a(x) ≥ a 0 for every x ∈ [0, L], and a W 3,∞ (0,L) ≤ α, where a 0 , α > 0. We define the operator

P = ∂ t + a(x)∂ xxx + b(x, t)∂ x + d(x, t) (7)
and the space

V = v ∈ L 2 (-T, T ; H 3 ∩ H 1 0 (0, L)) P v ∈ L 2 ((0, L) × (-T, T )) . (8)
Consider β ∈ C 3 ([0, L]) such that for some r > 0 we have

0 < r ≤ β(x) and 0 < r ≤ β (x), ∀x ∈ (0, L). (9)
We define, for each λ > 0, the functions

φ(x, t) = e 2λ β ∞ -e λβ(x) (T + t)(T -t) , θ(x, t) = e λβ(x) (T + t)(T -t) (10) for (x, t) ∈ [0, L] × (-T, T ). Let us notice that there exists a constant C = C(T ) such that 0 < 1 C ≤ θ(x, t) and |θ t (x, t)| ≤ Cθ 2 (x, t), ∀(x, t) ∈ [0, L] × (-T, T ). (11)
We will prove the following Carleman estimate.

Theorem 2 Let φ and θ be defined by [START_REF] Benabdallah | Stability of discontinuous diffusion coefficients and initial conditions in an inverse problem for the heat equation[END_REF], and let P be the linearized KdV operator defined in [START_REF] Bellassoued | Uniqueness and stability in determining the speed of propagation of second-order hyperbolic equation with variable coefficients[END_REF]. There exist s 0 > 0, λ 0 > 0 and a constant C(L, T, s 0 , λ 0 , α, a 0 , r) > 0 such that for every s ≥ s 0 and λ ≥ λ 0 ,

T -T L 0 e -2sφ s 5 λ 6 θ 5 |u| 2 + s 3 λ 4 θ 3 |u x | 2 + sλ 2 θ|u xx | 2 dxdt ≤ C T -T L 0 e -2sφ |P u| 2 dxdt + Cs 3 λ 3 T -T e -2sφ(L,t) θ 3 (L, t)|u x (L, t)| 2 dt + Csλ T -T e -2sφ(L,t) θ(L, t)|u xx (L, t)| 2 dt (12)
for all u ∈ V defined by [START_REF] Bellassoued | Logarithmic stability in the dynamical inverse problem for the Schrödinger equation by arbitrary boundary observation[END_REF].

Remark 4 When looking for a one-parameter Carleman estimate (as in [START_REF] Glass | Some exact controllability results for the linear KdV equation and uniform controllability in the zero-dispersion limit[END_REF][START_REF] Rosier | Control of the surface of a fluid by a wavemaker[END_REF] where

φ = β(x) (t+T )(T -t)
), some hypothesis on the second derivative of the function β have to be imposed. The use of the second parameter λ allows to omit them. Proof: Let s > 0 and define W s = {e -sφ u / u ∈ V}. For u ∈ V, we set w = e -sφ u and we define the operator P φ from W s to L 2 ((0, L) × (-T, T )) by P φ w = e -sφ P (e sφ w).

After straightforward computations, we write

P φ w = P 1 w + P 2 w + Rw
where, for a constant m > 0 to be chosen later,

P 1 w = w t + 3as 2 φ 2 x w x + aw xxx + 3ams 2 φ x φ xx w, (13)
P 2 w = as 3 φ 3 x w + 3asφ x w xx + 3sw x (aφ x ) x , (14)
and

Rw = bsφ x w + bw x + asφ xxx w + 3as 2 φ x φ xx w + dw +sφ t w -3sa x φ x w x -3ams 2 φ x φ xx w. (15)
Therefore,

P φ w -Rw 2 L 2 (Q) = P 1 w 2 L 2 (Q) + 2 P 1 w, P 2 w + P 2 w 2 L 2 (Q) (16)
where

•, • is the L 2 (Q) scalar product.
Step 1. Explicit calculations Notice that w(0, t) = w(L, t) = 0 for all t ∈ (-T, T ) and w(x, ±T ) = 0 for all x ∈ (0, L). This is due to the fact that w = e -sφ u, with u ∈ L 2 (-T, T ; H 1 0 (0, L)) and lim t→±T φ(x, t) = +∞.

By developing P 1 w, P 2 w , we get P 1 w, P 2 w = i=1,...,4; j=1,...,3

I i,j
where I i,j is the L 2 (Q) scalar product of the i-th term in [START_REF] Bukhgeȋm | Uniqueness in the large of a class of multidimensional inverse problems[END_REF] with the j-th term in [START_REF] Cardoulis | Inverse problem for the Schrödinger operator in an unbounded strip[END_REF]. By integrations by parts, we obtain:

• I 1,1 = - 3s 3 2 Q aφ 2 x φ xt |w| 2 dxdt = 3 2 s 3 λ 3 Q aβ 3 x θ 2 θ t |w| 2 dxdt ; • I 1,2 = -3s Q (aφ x ) x w x w t dxdt + 3s 2 Q aφ xt |w x | 2 dxdt = -3s Q (aφ x ) x w x w t dxdt - 3 2 sλ Q aβ x θ t |w x | 2 dxdt ; • I 1,3 = 3s Q (aφ x ) x w x w t dxdt ; • I 2,1 = - 3s 5 2 Q (a 2 φ 5 x ) x |w| 2 dxdt = 3 2 s 5 λ 5 Q (a 2 β 5 x ) x θ 5 |w| 2 dxdt + 15 2 s 5 λ 6 Q a 2 β 6 x θ 5 |w| 2 dxdt ; • I 2,2 = - 9s 3 2 Q (a 2 φ 3 x ) x |w x | 2 dxdt + 9s 3 2 T 0 a 2 φ 3 x |w x | 2 x=L x=0 dt = 9 2 s 3 λ 3 Q (a 2 β 3 x ) x θ 3 |w x | 2 dxdt + 27 2 s 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt - 9 2 s 3 λ 3 T 0 a 2 β 3 x θ 3 |w x | x=L x=0 dt ; • I 2,3 = 9s 3 Q aφ 2 x (aφ x ) x |w x | 2 dxdt = -9s 3 λ 3 Q aβ 2 x (aβ x ) x θ 3 |w x | 2 dxdt -9s 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt ; • I 3,1 = 3s 3 2 Q (a 2 φ 3 x ) x |w x | 2 dxdt - s 3 2 Q (a 2 φ 3 x ) xxx |w| 2 dxdt - s 3 2 T 0 a 2 φ 3 x |w x | 2 x=L x=0 dt = - 3 2 s 3 λ 3 Q (a 2 β 3 x ) x θ 3 |w x | 2 dxdt - 9 2 s 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt + s 3 λ 3 2 Q (a 2 β 3 x θ 3 ) xxx |w| 2 dxdt + s 3 λ 3 2 T 0 a 2 β 3 x θ 3 |w x | x=L x=0 dt ; • I 3,2 = - 3s 2 Q (a 2 φ x ) x |w xx | 2 dxdt + 3s 2 T 0 a 2 φ x |w xx | 2 x=L x=0 dt = 3 2 sλ Q (a 2 β x ) x θ|w xx | 2 dxdt + 3 2 sλ 2 Q a 2 β 2 x θ|w xx | 2 dxdt - 3 2 sλ T 0 a 2 β x θ|w xx | x=L x=0 dt ; • I 3,3 = -3s Q a(aφ x ) x |w xx | 2 dxdt + 3s 2 Q [a(aφ x ) x ] xx |w x | 2 dxdt - 3s 2 T 0 (a(aφ x ) x ) x |w x | 2 x=L x=0 dt + 3s T 0 a(aφ x ) x w x w xx x=L x=0 dt = 3sλ Q a(aβ x ) x θ|w xx | 2 dxdt + 3sλ 2 Q a 2 β 2 x θ|w xx | 2 dxdt + 3 2 sλ T 0 (a(aβ x θ) x ) x |w x | 2 x=L x=0 dt - 3 2 sλ Q (a(aβ x θ) x ) xx |w x | 2 dxdt -3sλ T 0 a(aβ x θ) x w x w xx x=L x=0 dt ; • I 4,1 = 3ms 5 Q a 2 φ 4 x φ xx |w| 2 dxdt = -3ms 5 λ 5 Q a 2 β 4 x β xx θ 5 |w| 2 dxdt -3ms 5 λ 6 Q a 2 β 6 x θ 5 |w| 2 dxdt ; • I 4,2 = -9ms 3 Q a 2 φ 2 x φ xx |w x | 2 dxdt + 9ms 3 2 Q (a 2 φ 2 x φ xx ) xx |w| 2 dxdt = 9ms 3 λ 3 Q a 2 β 2 x β xx θ 3 |w x | 2 dxdt + 9ms 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt - 9m 2 s 3 λ 3 Q a 2 β 2 x β xx θ 3 xx |w| 2 dxdt - 9m 2 s 3 λ 4 Q a 2 β 4 x θ 3 xx |w| 2 dxdt ; • I 4,3 = - 9ms 3 2 Q (aφ x φ xx (aφ x ) x ) x |w| 2 dxdt = 9m 2 s 3 λ 3 Q aβ x β xx θ 2 (aβ x θ) x x |w| 2 dxdt + 9m 2 s 3 λ 4 Q aβ 3 x θ 2 (aβ x θ) x x |w| 2 dxdt.
Therefore,

P 1 w, P 2 w = 15 2 -3m s 5 λ 6 Q a 2 β 6 x θ 5 |w| 2 dxdt ( 17
)
+ 9ms 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt + 9 2 sλ 2 Q a 2 β 2 x θ|w xx | 2 dxdt + B + X 1
where X 1 and B satisfy the following. The term X 1 gathers the non-dominating terms and satisfies that there exists a positive constant C = C(L, T, α, r) independent of s and λ, (which may change from line to line) such that, using [START_REF] Bona | Solutions of the Korteweg-de Vries equation in fractional order Sobolev spaces[END_REF], we get

|X 1 | ≤ C(s 3 λ 6 + s 5 λ 5 ) Q θ 4 |w| 2 dxdt + C(sλ 4 + s 3 λ 3 ) Q θ 2 |w x | 2 dxdt. ( 18
)
The term B gathers the boundary terms. We split them as follows

B := B + + B -+ B * (19)
with the notations

B + = 4s 3 λ 3 T -T a 2 β 3 x θ 3 x=0 |w x (0, t)| 2 dt + 3 2 sλ T -T a 2 β x θ x=0 |w xx (0, t)| 2 dt, B -= -4s 3 λ 3 T -T a 2 β 3 x θ 3 x=L |w x (L, t)| 2 dt - 3 2 sλ T -T a 2 β x θ x=L |w xx (L, t)| 2 dt, B * = + 3 2 sλ T -T (a(aβ x θ) x ) x x=L |w x (L, t)| 2 dt - 3 2 sλ T -T (a(aβ x θ) x ) x x=0 |w x (0, t)| 2 dt -3sλ T -T a(aβ x θ) x x=L w x (L, t)w xx (L, t)dt + 3sλ T -T a(aβ x θ) x x=0 w x (0, t)w xx (0, t)dt.
Step 2. Dominating terms By choosing, and fixing, m such that

0 < m < 5 2 , (20)
we obtain that there exists a constant C = C(L, T, a 0 , r) > 0 such that:

15 2 -3m s 5 λ 6 Q a 2 β 6 x θ 5 |w| 2 dxdt ≥ Cs 5 λ 6 Q θ 5 |w| 2 dxdt; 9ms 3 λ 4 Q a 2 β 4 x θ 3 |w x | 2 dxdt ≥ Cs 3 λ 4 Q θ 3 |w x | 2 dxdt; 9 2 sλ 2 Q a 2 β 2 x θ|w xx | 2 dxdt ≥ Csλ 2 Q θ|w xx | 2 dxdt.
We introduce the notation given by

w 2 s,λ,θ = s 5 λ 6 Q θ 5 |w| 2 dxdt + s 3 λ 4 Q θ 3 |w x | 2 dxdt + sλ 2 Q θ|w xx | 2 dxdt.
From [START_REF] Coron | Control and nonlinearity[END_REF] we obtain that

|X 1 | ≤ C 1 s 2 + 1 λ w 2 s,λ,θ .
Therefore, choosing s and λ large enough, from [START_REF] Cerpa | Boundary controllability for the nonlinear Korteweg-de Vries equation on any critical domain[END_REF] we get some C > 0 depending only on L, T, α, a 0 , r, s 0 and λ 0 and is such that, for all s > s 0 , λ > λ 0 ,

C w 2 s,λ,θ ≤ P 1 w, P 2 w -B. (21)
Step 3. Boundary terms From ( 19) and ( 21), considering that B + and -B -are non-negative, we write

C w 2 s,λ,θ + B + ≤ P 1 w, P 2 w -B --B * . (22)
Now, we shall deal with the boundary terms in B * , that has no sign. We will bound |B * | by above. For the four integrals one by one, using [START_REF] Bona | Solutions of the Korteweg-de Vries equation in fractional order Sobolev spaces[END_REF] and Cauchy-Schwarz inequality, this gives:

3 2 sλ T -T (a(aβ x θ) x ) x x=L |w x (L)| 2 dt ≤ Csλ 3 T -T a 2 (L)β 3 x (L)θ(L)|w x (L)| 2 dt, 3 2 sλ T -T (a(aβ x θ) x ) x x=0 |w x (0)| 2 dt ≤ Csλ 3 T -T a 2 (0)β 3 x (0)θ(0)|w x (0)| 2 dt, 3sλ T -T a(aβ x θ) x x=L w x (L)w xx (L)dt ≤ Cs 2 λ 3 T -T a 2 (L)β 3 x (L)θ(L)|w x (L)| 2 dt + Cλ T -T a 2 (L)β x (L)θ(L)|w xx (L)| 2 dt, and
3sλ T -T a(aβ x θ) x x=0 w x (0)w xx (0)dt ≤ Cs 2 λ 3 T -T a 2 (0)β 3 x (0)θ(0)|w x (0)| 2 dt + Cλ T -T a 2 (0)β x (0)θ(0)|w xx (0)| 2 dt.
From these last four inequalities we get that

|B * | ≤ C s B + + C s |B -|.
Therefore, taking s 0 large enough, inequality [START_REF] Fursikov | Controllability of evolution equations[END_REF] becomes, for all s > s 0 ,

w 2 s,λ,θ + B + ≤ C P 1 w, P 2 w + C|B -|. (23)
Step 4. Carleman estimate for w. Since B + ≥ 0, inequality (23) yields

w 2 s,λ,θ ≤ C P 1 w, P 2 w + Cs 3 λ 3 T -T a 2 (L)β 3 x (L)θ 3 (L, t)|w x (L, t)| 2 dt + Csλ T -T a 2 (L)β x (L)θ(L, t)|w xx (L, t)| 2 dt. ( 24
)
Moreover, from [START_REF] Cardoulis | Simultaneous identification of the diffusion coefficient and the potential for the Schrödinger operator with only one observation[END_REF] we have

P φ w -Rw 2 L 2 (Q) ≤ 2 P φ w 2 L 2 (Q) + 2 Rw 2 L 2 (Q) ≤ 2 P φ w 2 L 2 (Q) + Cs 4 λ 6 Q θ 2 |w| 2 dxdt + Cs 2 λ 2 Q |w x | 2 dxdt.
Therefore, choosing again s 0 and λ 0 large enough, we have proved the following.

Proposition 1 There exist s 0 , λ 0 > 0 and a constant C = C(L, T, s 0 , λ 0 , α, a 0 , r) > 0 such that for all s ≥ s 0 , for all λ ≥ λ 0 ,

s 5 λ 6 T -T L 0 θ 5 |w| 2 dxdt + s 3 λ 4 T -T L 0 θ 3 |w x | 2 dxdt + sλ 2 T -T L 0 θ|w xx | 2 dxdt + T -T L 0 |P 1 w| 2 dxdt + T -T L 0 |P 2 w| 2 dxdt ≤ C T -T L 0 |P φ w| 2 dxdt + Cs 3 λ 3 T -T θ 3 (L, t)|w x (L, t)| 2 dt + Csλ T -T θ(L, t)|w xx (L, t)| 2 dt (25)
for all w ∈ W s := {e -sφ v : v ∈ V}.

Step 5. Back to the variable u The Carleman estimate stated in Theorem 2 will now be deduced from Proposition 1.

It is only a matter of going back to the variable u ∈ V. Recall that u = e sφ w and u(0, t) = u(L, t) = 0 for all t ∈ (-T, T ). One easily checks that there exists a positive constant C = C(L, T, s 0 , λ 0 , α, a 0 , r) such that for all x ∈ (0, L) and t ∈ (-T, T ),

e -2sφ |u x | 2 ≤ 2|w x | 2 + Cs 2 λ 2 θ 2 |w| 2
, and

e -2sφ |u xx | 2 = 2|w xx | 2 + Cs 2 λ 2 θ 2 |w x | 2 + Cs 4 λ 4 θ 2 |w| 2 .
Hence the left hand side of ( 12) is estimated by the left hand side of [START_REF] Grimshaw | Solitary waves propagating over variable topography[END_REF]. Moreover, we have P φ w = e -sφ P u, and concerning the boundary terms, using u(0, t) = u(L, t) = 0 for all t ∈ (-T, T ) we obtain that in {0, L} × (-T, T ),

|w x | 2 ≤ Ce -2sφ |u x | 2 ,
and
|w xx | 2 ≤ Ce -2sφ |u xx | 2 + s 2 λ 2 θ 2 |u x | 2 .
Therefore, for s 0 and λ 0 large enough, Theorem 2 directly follows from Proposition 1.

Inverse Problem

This section is devoted to the proof of the Lipschitz stability result stated in Theorem 1 about the inverse problem of retreiving the main coefficient a in Equation ( 2) from boundary measurements of the solution. We recall here that we will use the Bukhgeȋm-Klibanov method. In the sake of clarity, we divide the proof in several steps.

Step 1. Local solution of the inverse problem

We consider two coefficients a, ã and the corresponding solutions of (2), y and ỹ, and we define u(x, t) := y(x, t) -ỹ(x, t), and σ(x) := ã(x) -a(x).

The function u satisfies

   u t + a(x)u xxx + (1 + ỹ)u x + y x u = σ ỹxxx , ∀(x, t) ∈ (0, L) × (0, T ), u(0, t) = 0, u(L, t) = 0, u x (L, t) = 0 ∀t ∈ (0, T ), u(x, 0) = 0, ∀x ∈ (0, L). (26)
We consider now z(x, t) := u t (x, t), which satisfies the system

   z t + a(x)z xxx + (1 + ỹ)z x + y x z = f, ∀(x, t) ∈ (0, L) × (0, T ), z(0, t) = 0, z(L, t) = 0, z x (L, t) = 0 ∀t ∈ (0, T ), z(x, 0) = σ(x)y 0 (x), ∀x ∈ (0, L), (27)
where f = σ(x)ỹ xxxt -y xt u -ỹt u x . One should notice that we have now σ, that we seek to estimate, appearing not only in the source term of the equation, but also into the initial condition.

Step 2. Extension to negative time In order to estimate σ through z(x, 0) using the Carleman estimate of the previous section, we shall extend the partial differential equations ( 26) in u and ( 27) in z to negative times. In order to extend these linearized KdV equations to the interval (-T, T ), we define the symmetric extension of any function g defined on [0, L] × [0, T ] by

ĝ(x, t) = g(x, t) if x ∈ [0, L], t ∈ [0, T ], g(L -x, -t) if x ∈ [0, L], t ∈ [-T, 0). ( 28
)
One should notice that this extension satisfies f ĝ = f g. We also define the following anti-symmetric extension to [0, L] × [-T, T ] of any function g defined on [0, L] × [0, T ] by

ǧ(x, t) = g(x, t) if x ∈ [0, L], t ∈ [0, T ], -g(L -x, -t) if x ∈ [0, L], t ∈ [-T, 0). ( 29
)
Therefore, we first assume that the initial data z 0 (x) = z(x, 0) of equation ( 27)

satisfies z 0 (x) = z 0 (L -x) ∀x ∈ [0, L], (30)
meaning that we have to assume this for a, ã and y 0 , as we did in assumptions ( 3) and ( 5). We define the extension of y (resp. ỹ) solution of equation ( 2) by ŷ (resp. ŷ). Now we can define v := ẑ on [0, L] × [-T, T ], which satisfies the equation given by

           v t + a(x)v xxx + (1 + ŷ)v x + ŷx v = f , ∀x ∈ (0, L), t ∈ (-T, T ), v(0, t) = 0, v(L, t) = 0, ∀t ∈ (-T, T ), v x (L, t) = 0, ∀t ∈ (0, T ), v x (L, t) = -z x (0, -t), ∀t ∈ (-T, 0), v(x, 0) = σ(x)y 0 (x), ∀x ∈ (0, L). (31)
The partial differential operator of this equation is eligible for the Carleman estimate of Theorem 2, according to [START_REF] Bellassoued | Uniqueness and stability in determining the speed of propagation of second-order hyperbolic equation with variable coefficients[END_REF] with b = 1 + ŷ and d = ŷx , provided that b ∈ L ∞ (-T, T ; W 1,∞ (0, L)) and d ∈ L ∞ ((-T, T ) × (0, L)). This regularity is fulfilled thanks to the hypothesis on the regularity of y 0 , a and ã, which ensures the corresponding solutions y and ỹ to be regular enough. See Proposition 3 in the Appendix.

Step 3. First use of the Carleman estimate Following the definitions stated in the proof of Theorem 1, we set w = e -sφ v and since we have w(0, t) = w(L, t) = 0 and w(x, ±T ) = 0 (because of the weight function φ), we can either compute or estimate the following integral

I := 0 -T L 0 w P 1 w dxdt = 0 -T L 0 w(w t + 3as 2 φ 2 x w x + aw xxx + 3as 2 φ x φ xx w) dxdt (32)
where P 1 has been defined in [START_REF] Bukhgeȋm | Uniqueness in the large of a class of multidimensional inverse problems[END_REF] and m = 1 chosen according to [START_REF] Crépeau | A reduced model of pulsatile flow in an arterial compartment[END_REF]. First of all, making integrations by parts, it is not difficult to see that

I = 1 2 L 0 |w(x, 0)| 2 dx + R, where R = 0 -T L 0 3as 2 φ 2 x w x w + aw xxx w + 3as 2 φ x φ xx |w| 2 dxdt = - 3 2 s 2 λ 2 0 -T L 0 (aβ 2 x θ 2 ) x |w| 2 dxdt - 0 -T L 0 aw xx w x dxdt + 0 -T L 0 a x |w x | 2 dxdt - 1 2 0 -T L 0 a xxx |w| 2 dxdt + 3ms 2 λ 2 0 -T L 0 aβ x (β xx + λβ 2 x )θ 2 |w| 2 dxdt
Using a ∈ Σ(a 0 , α), β ∈ C 3 ([0, L]), the property (11) of θ and Cauchy-Schwartz inequality, one obtains

|R| ≤ C T -T L 0 (s 2 λ 3 θ 3 |w| 2 + |w x | 2 + |w xx | 2 ) dxdt
for the same generic constant C > 0 as in the previous section. Hence, estimating the quantity s 5 2 I by Carleman inequality (25), and choosing s 0 large enough to absorb the terms of R by the dominant ones of the left hand side of (25), we prove

s 5 2 L 0 |w(x, 0)| 2 dx ≤ 2s 5 T -T L 0 |w| 2 dxdt + 2 T -T L 0 |P 1 w| 2 dxdt + 2s 5 2 |R| ≤ C T -T L 0 e -2sφ | f (x, t)| 2 dxdt + Cs 3 λ 3 T -T θ 3 (L, t)|w x (L, t)| 2 dt + Csλ T -T θ(L, t)|w xx (L, t)| 2 dt.
On the one hand, since we assume that |y 0 (x)| ≥ r 0 > 0, we have Therefore, by these two last estimates, fixing λ > λ 0 , we get

s 5 2 L 0 e -2sφ(x,0) |σ(x)| 2 dx ≤ C T -T L 0 e -2sφ | f (x, t)| 2 dxdt + Cs 3 T -T θ 3 (L, t)|w x (L, t)| 2 dt + Cs T -T θ(L, t)|w xx (L, t)| 2 dt
and using the definition of f from [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF], and the fact that θ and φ are even in time, w = e -sφ v = e -sφ ẑ with definition [START_REF] Imanuvilov | Global Lipschitz stability in an inverse hyperbolic problem by interior observations[END_REF], and the boundary properties of z in ( 27), we obtain

s 5 2 L 0 e -2sφ(x,0) |σ(x)| 2 dx ≤ C T 0 L 0 e -2sφ(x,t) + e -2sφ(L-x,t) |f (x, t)| 2 dxdt + Cs 3 T 0 e -2sφ(L,t) θ 3 (L, t)|z x (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|z xx (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|z xx (L, t)| 2 dt. (33)
On the other hand, since φ(x, t) ≥ φ(x, 0) for all (x, t) ∈ (0, L)×(-T, T ), and ỹt , y xt and ỹxxxt belong to L ∞ ((0, L)×(0, T )) (recall that y 0 ∈ {w ∈ H 7 (0, L); w(0) = w(L) = w (L) = 0} and a, ã ∈ Σ(a 0 , α), which implies that y, ỹ ∈ W 1,∞ (0, T ; W 3,∞ (0, L))), we can write

T 0 L 0 e -2sφ(x,t) |f (x, t)| 2 dxdt = T 0 L 0 e -2sφ(x,t) |σ(x)ỹ xxxt -y xt u -ỹt u x | 2 dxdt ≤ C L 0 e -2sφ(x,0) |σ(x)| 2 dx + C T 0 L 0 e -2sφ(x,t) (|u| 2 + |u x | 2 ) dxdt.
and also,

T 0 L 0 e -2sφ(L-x,t) |f (x, t)| 2 dxdt = T 0 L 0 e -2sφ(L-x,t) |σ(x)ỹ xxxt -y xt u -ỹt u x | 2 dxdt ≤ C L 0 e -2sφ(L-x,0) |σ(x)| 2 dx + C T 0 L 0 e -2sφ(L-x,t) (|u| 2 + |u x | 2 ) dxdt ≤ C L 0 e -2sφ(x,0) |σ(x)| 2 dx + C 0 -T L 0 e -2sφ(x,t) (|û| 2 + |û x | 2 ) dxdt Thus, we get, T -T L 0 e -2sφ | f (x, t)| 2 dxdt ≤ C L 0 e -2sφ(x,0) |σ(x)| 2 dx + C T -T L 0 e -2sφ(x,t) (|û| 2 + |û x | 2 ) dxdt (34)
Here, the second term in the right hand side has to be estimated, and this is done using again a Carleman inequality.

Step 4. Second use of the Carleman estimate Applying now Carleman estimate [START_REF] Bona | A Nonhomogeneous Boundary-Value Problem for the Korteweg-de Vries Equation Posed on a Finite Domain[END_REF] to the equation satisfied by û that is the extension of (26) to negative times, and where the first and zero-th order potentials are b(x, t) = 1 + ŷ(x, t), and d(x, t) = ŷx (x, t), we obtain

T -T L 0 e -2sφ (|û| 2 + |û x | 2 ) dxdt ≤ C T -T L 0 e -2sφ |σ ỹxxx | 2 dxdt + Cs 3 T -T e -2sφ(L,t) θ 3 (L, t)|û x (L, t)| 2 dt + Cs T -T e -2sφ(L,t) θ(L, t)|û xx (L, t)| 2 dt ≤ C L 0 e -2sφ(x,0) |σ(x)| 2 dxdt + Cs 3 T 0 e -2sφ(L,t) θ(L, t)|u x (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|u xx (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|u xx (L, t)| 2 dt.
Thus, from [START_REF] Klibanov | Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems[END_REF], we can write, [START_REF] Klibanov | Newton-Kantorovich method for threedimensional potential inverse scattering problem and stability of the hyperbolic Cauchy problem with time-dependent data[END_REF] and from [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF] and [START_REF] Klibanov | Newton-Kantorovich method for threedimensional potential inverse scattering problem and stability of the hyperbolic Cauchy problem with time-dependent data[END_REF], choosing s 0 large enough, we deduce that

T 0 L 0 e -2sφ | f (x, t)| 2 dxdt ≤ C L 0 e -2sφ(x,0) |σ(x)| 2 dxdt + Cs 3 T 0 e -2sφ(L,t) θ(L, t)|u x (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|u xx (0, t)| 2 dt + Cs T 0 e -2sφ(L,t) θ(L, t)|u xx (L, t)| 2 dt
s 5 2 L 0 e -2sφ(0) |σ(x)| 2 dx ≤ Cs 3 T 0 e -2sφ(L,t) θ(L, t)(|u x (0, t)| 2 + |z x (0, t)| 2 ) dt + Cs T 0 e -2sφ(L,t) θ(L, t)(|u xx (0, t)| 2 + |z xx (0, t)| 2 ) dt + Cs T 0 e -2sφ(L,t) (|u xx (L, t)| 2 + |z xx (L, t)| 2 ) dt (36)
Taking into account that z = u t = ∂ t (y -ỹ) the result of Theorem 1 directly follows from [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF].

A Cauchy problem

In this appendix, we state the well-posedness results for the KdV equation considered in the paper. We only give the main ideas in the proofs because the tools are pretty similar to those applied for the constant coefficient case, which is already standard in the literature (see [START_REF] Bona | Solutions of the Korteweg-de Vries equation in fractional order Sobolev spaces[END_REF], [START_REF] Bona | A Nonhomogeneous Boundary-Value Problem for the Korteweg-de Vries Equation Posed on a Finite Domain[END_REF], [START_REF] Rosier | Exact boundary controllability for the Korteweg-de Vries equation on a bounded domain[END_REF], [START_REF] Coron | Exact boundary controllability of a nonlinear KdV equation with critical lengths[END_REF] and [START_REF] Coron | Control and nonlinearity[END_REF]). It is worthwhile to emphasis that these results allow to get, under smoothness hypothesis on the data, the existence of solutions for the KdV equation with the regularity required for our approach in order to prove the stability of the inverse problem.

We consider a main coefficient satisfying

∃a 0 ∈ R, ∀x ∈ [0, L] a(x) ≥ a 0 > 0. ( 37
)
Let us introduce, for any s ≥ 0 the space

B s = C([0, T ], H s (0, L)) ∩ L 2 (0, T ; H s+1 (0, L))
and consider only the principal part of the equation, i.e. we study the problem

       y t + a(x)y xxx = f, ∀(x, t) ∈ (0, L) × (0, T ), y(0, t) = 0, y(L, t) = 0, ∀t ∈ (0, T ), y x (L, t) = 0, ∀t ∈ (0, T ), y(0, x) = y 0 (x), ∀x ∈ (0, L). (38)
In order to take into account the coefficient a = a(x) we work in the space X = L 2 (0, L) endowed with the inner product

w 1 , w 2 X := L 0 1 a(x) w 1 (x)w 2 (x)dx.
Notice that since a ∈ L ∞ (0, L) satisfies ( 37), the norm defined by the inner product •, • X is equivalent to the L 2 -norm.

In the domain D(A) = {w ∈ H 3 (0, L) w(0) = w(L) = w (L) = 0}, we define the operator A : D(A) ⊂ X → X as A(w) = -aw . It is easy to see that both A and its adjoint operator A * are dissipative, and therefore, by standard semigroup theory (for instance a corollary of Lumer-Phillips Theorem, see Chapter 1 in [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]), we get that A generates a strongly continuous semigroup in L 2 (0, L).

Thus, if y 0 ∈ L 2 (0, L), a ∈ L ∞ (0, L) satisfies (37) and f ∈ L 1 (0, T, L 2 (0, L)) then the linear KdV equation ( 38) has a unique solution (called mild solution) in the space C([0, T ], L 2 (0, L)). Applying the multipliers technique, we get a Kato smoothing effect, which implies y ∈ L 2 (0, T ; H 1 (0, L)) and thus y ∈ B 0 . This regularity is indeed obtained by multiplying formally equation ( 38) by x a(x) y(t, x) and integrating in (0, L) × (0, T ). This proves the following result.

Proposition 2 Let y 0 ∈ L 2 (0, L), a ∈ L ∞ (0, L) verifying (37) and f ∈ L 1 (0, T, L 2 (0, L)).
Then the linear KdV equation (38) has a unique solution in the space B 0 . Moreover, there exist C > 0 such that

y B 0 ≤ C y 0 L 2 (0,T ) + f L 1 (0,T,L 2 (0,L)) .
Depending on the regularity of the data, we can prove the existence of more regular solutions. Let us consider y solution of [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]. The function u := y t satisfies

       u t + a(x)u xxx = f t , ∀(x, t) ∈ (0, L) × (0, T ), u(0, t) = 0, u(L, t) = 0, ∀t ∈ (0, T ), u x (L, t) = 0, ∀t ∈ (0, T ), u(0, x) = f (0, x) -a(x)y 0xxx , ∀x ∈ (0, L). (39)
Assuming that f t ∈ L 1 (0, T ; L 2 (0, L)) and y 0 ∈ D(A), we have that (f (0, x)a(x)y 0xxx ) ∈ L 2 (0, L), and then Proposition 2 implies that u = y t ∈ B 0 .

If additionally we ask f ∈ B 0 , then by equation [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF] we have a(x)y xxx = f -y t ∈ B 0 , and if a ∈ W 1,∞ (0, L) satisfies [START_REF] Mercado | Lionel Inverse problems for the Schrödinger equation via Carleman inequalities with degenerate weights[END_REF] we get that the original solution satisfies y ∈ B 3 .

Analyzing in the same way each one of the equations fulfilled by y tt and y ttt , we get the following result.

Proposition 3 Let y 0 ∈ H 9 (0, L) ∩ D(A), a ∈ W 6,∞ (0, L) verifying (37), and f ∈ B 6 such that f t ∈ B 3 , f tt ∈ B 0 , and f ttt ∈ L 1 (0, T, L 2 (0, L)). Then the linear KdV equation (38) has a unique solution in the space B 9 . This solution also satisfies y t ∈ B 6 , y tt ∈ B 3 and y ttt ∈ B 0 .

Remark 5 The hypothesis a ∈ W 6,∞ (0, L) is used to ensure that the initial condition of the problem solved by y ttt belongs to the space L 2 (0, L).

In order to prove a similar result for the nonlinear system        y t + a(x)y xxx + y x + yy x = 0, ∀(x, t) ∈ (0, L) × (0, T ), y(0, t) = 0, y(L, t) = 0, ∀t ∈ (0, T ), y x (L, t) = 0, ∀t ∈ (0, T ), y(0, x) = y 0 (x), ∀x ∈ (0, L),

we define the map Π : ỹ ∈ B 0 → Π(ỹ) = y ∈ B 0

where Π(ỹ) = y is the solution of [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF] with f replaced by (-ỹ x -ỹ ỹx ). For ỹ ∈ B 0 , we have ỹ ∈ L 2 (0, T ; H 1 (0, L)), and therefore ỹ ỹx ∈ L 1 (0, T ; L 2 (0, L)). Thus we are able to use Proposition 2 with the right-hand side f = (-ỹ x -ỹ ỹx ) to define Π(ỹ) = y ∈ B 0 . Applying a fixed point argument in small time and global a priori estimates to deal with the long-time case, we obtain the following Proposition 4 Let y 0 ∈ L 2 (0, L) and a ∈ L ∞ (0, L) verifying [START_REF] Mercado | Lionel Inverse problems for the Schrödinger equation via Carleman inequalities with degenerate weights[END_REF]. Then, the nonlinear KdV equation (40) has a unique solution in the space B 0 .

We can also prove the nonlinear result in the more regular framework.

Proposition 5 Let y 0 ∈ H 9 (0, L) ∩ D(A) and a ∈ W 6,∞ (0, L) verifying [START_REF] Mercado | Lionel Inverse problems for the Schrödinger equation via Carleman inequalities with degenerate weights[END_REF]. Then, the nonlinear KdV equation (40) has a unique solution in the space B 9 .

The proof of Proposition 5 can be done as before (with a fixed point argument and global a priori estimates) with the same map Π but now defined in the space Y := {y ∈ B 9 y t ∈ B 6 , y tt ∈ B 3 and y ttt ∈ B 0 } instead of B 0 . Let us notice that the nonlinearity ỹ ỹx satisfies all the hypothesis required in Proposition 3 for the right hand side f provided that ỹ ∈ Y.

Using some interpolation arguments, we can prove the following result, which gives us the existence of solutions as we require in order to solve our Inverse Problem in Theorem 1.

Theorem 3 Let y 0 ∈ H 7 (0, L) ∩ D(A) and a ∈ W 6,∞ (0, L) satisfying [START_REF] Mercado | Lionel Inverse problems for the Schrödinger equation via Carleman inequalities with degenerate weights[END_REF]. Then the nonlinear KdV equation (40) has a unique solution in the space B 7 . Furthermore y t ∈ B 4 and therefore y ∈ W 1,∞ (0, T ; W 3,∞ (0, L)).

Indeed, the proof of this result can be deduced from the Interpolation Theorem by Bona and Scott ( [START_REF] Bona | Solutions of the Korteweg-de Vries equation in fractional order Sobolev spaces[END_REF], see Theorem 4.3 in [12]). This theorem is applied in the following way. Denoting by S the mapping which sends the initial condition y 0 to the solution y of the equation [START_REF] Puel | Generic well-posedness in a multidimensional hyperbolic inverse problem[END_REF], by Propositions 4 and 5 we have that S : L 2 (0, L) -→ B 0 and S : H 9 (0, L) ∩ D(A) -→ B 9 are well defined. Additionally, two inequalities are needed (see (i) and (ii) from Theorem 4.3 in [START_REF] Bona | A Nonhomogeneous Boundary-Value Problem for the Korteweg-de Vries Equation Posed on a Finite Domain[END_REF]):

• The first one is the estimation of S(y 01 )-S(y 02 ) B 0 in terms of y 01 -y 02 L 2 (0,L) , which can be proven exactly in the same way as in the proof of Proposition 4.2 of [START_REF] Bona | A Nonhomogeneous Boundary-Value Problem for the Korteweg-de Vries Equation Posed on a Finite Domain[END_REF], but now using our Proposition 2.

• The second one consists in estimating S(y 0 ) B 9 by y 0 H 9 (0,L) , which can be proven from Proposition 5 and following Proposition 15 in [START_REF] Coron | Exact boundary controllability of a nonlinear KdV equation with critical lengths[END_REF].

From the Interpolation Theorem we get that S : H 7 (0, L) ∩ D(A) → B 7 is well defined, with gives us the proof of Theorem 3.

Remark 6 A sharper result can be obtained with y 0 ∈ H s (0, L)∩D(A) where s > 6+ 1 2 . In that case, previous arguments give the existence of a solution y ∈ B s with y t ∈ B s-3 , and then y txxx ∈ C([0, T ], H q (0, L)) with q > 1/2. Thus, we obtain for system (40) a solution y ∈ W 1,∞ (0, T ; W 3,∞ (0, L)).

ee

  -2sφ(0) |σ(x)y 0 (x)| 2 dx ≥ s -2sφ(0) |σ(x)| 2 dx.

Acknowledgements. This work has been partially supported by MathAmsud project CIP-PDE, by the ANR Project CISIFS number NT09-437023 (L. Baudouin, E. Crépeau), by Fondecyt grant 1120610 (A. Mercado), by Basal CMM U. de Chile, and by Conicyt grant ACT 1106.


Jan Bartovsky
email: bartovsj@esiee.fr

Petr Dokládal
email: dokladae@esiee.fr

Eva Dokladalova

Vjaceslav Georgiev

Jan Bártovský

Eva Dokládalová

Parallel Implementation of Sequential Morphological Filters

Keywords: Mathematical Morphology, Serial Filters, Real-Time Implementation, Dedicated Hardware 1 Introduction Mathematical Morphology is very popular, self-contained, image processing framework providing

Many useful morphological filters are built as more or less long concatenations of erosions and dilations: openings, closings, size distributions, sequential filters, etc.

An efficient implementation of such concatenation would allow all the sequentially concatenated operators run simultaneously, on the time-delayed data. A recent algorithm (see below) for the morphological dilation/erosion allows such inter-operator parallelism.

This paper introduces an additional, intra-operator level of parallelism in this dilation/erosion algorithm. Realized in a dedicated hardware, for rectangular structuring elements with programmable size, such an implementation allows to obtain previously unachievable, real-time performances for these traditionally costly operators. Low latency and memory requirements are the main benefits when the performance is not deteriorated even for long concatenations or high-resolution images.

tools from filtering, multi-scale image analysis to pattern recognition. It has been used in a number of applications, including the biomedical and medical imaging, video surveillance, industrial control, video compression, stereology, or remote sensing ( [8, 16, 18, 19]).

In image-interpretation applications requiring a high correct-decision liability, one often use robust multi-criteria and/or multi-scale analysis. It generally consists of a serial concatenation of alternating atomic operators dilation and erosion with a progressively increasing computing window, the so-called structuring element (SE). Its examples include:

-Alternate Sequential Filters (ASF) -that are concatenations of openings and closings with a progressively increasing structuring element, useful for multi-scale analysis [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Serra | An overview of morphological filtering[END_REF]. -Size distributions -(aka granulometries) are concatenations of openings allowing measuring the size distribution in a population of objects [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF][START_REF] Sabourin | Off-line signature verification by local granulometric size distributions[END_REF][START_REF] Vincent | Granulometries and opening trees[END_REF]. -Statistical learning -a selected set of morphological operators ζ i can be separately applied to an image f . Then for every pixel f (x, y), the vector of values ζ i ( f )(x, y) can serve as a vector of descriptors for the pixel-wise learning and classification [START_REF] Cord | Segmentation of random textures by morphological and linear operators[END_REF].

Although built from basic blocks (the dilation and erosion), these operators are costly due to the number of iterations. The real-time capabilities (i.e., low latency) are even more difficult to achieve due to the sequential data dependence and high memory requirements.

The recently introduced algorithm for the dilation and erosion [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] shows how to handle efficiently the implementation of such concatenations. It enables an inter-operator level of parallelism where all the sequentially concatenated operators can run simultaneously, on time-delayed data. Obviously, it is fully exploited only if the algorithm is implemented in an adequate hardware (HW).

In this paper, we propose a HW implementation of the original algorithm and we introduce an additional, intra-operator level of parallelism. Such an implementation allows obtaining previously unachievable, real-time performances for these traditionally costly operators. Section 2 discusses the state of the art of existing dilation/erosion algorithms and concludes by the novelties presented in this paper. Section 3 and 4 recall the definitions and algorithmic principles and Section 5 illustrates the functional scheme of a sequential HW implementation. Then, Section 6 introduces the parallel implementation allowing obtaining an additional performance increase. Finally, Section 7 presents results obtained on FPGA.

Fast Implementations of Morphological Filters

During the last decades, propositions of optimized implementations concentrated on the efficiency of computing the dilation and erosion. The majority of authors measures this efficiency as a number of comparisons per pixel. Nevertheless, the minimization of comparisons can result in high memory requirements. It can even penalize the execution time since the overall latency issues are neglected.

For the following, we define as operator latency the latency introduced by the dependence of the result on future data samples. For example the max filter y i = max(x i-2 , x i-1 , . . . , x i+2 ) has operator latency 2. We define as algorithm latency any additional latency introduced by the algorithm, e.g., the necessity to perform a reverse scan on data. The latency is a time-less measure expressed in a number of data samples.

Note that there is also an additional delay, called computing latency, induced by the time needed to compute the result after all data are available. It is a platform dependent measure, independent of two previous latency definitions. In the example above, the polyadic max can either be executed sequentially on sequential machines, or in parallel on the dedicated hardware.

Then the overall latency of the system is the sum of these three terms.

Algorithmic Advances

The most efficient dilation/erosion algorithms are based on the SE decomposition to a set of basic, more easily optimized shapes, see [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF][START_REF] Xu | Decomposition of convex polygonal morphological structuring elements into neighborhood subsets[END_REF][START_REF] Zhuang | Morphological structuring element decomposition[END_REF]. A special attention is paid to 1-D algorithms obtaining a significant gain in the overall performance.

The most popular 1-D algorithm is called HGW (published by van Herk [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels[END_REF], and Gill and Werman [START_REF] Gil | Computing 2-d min, median, and max filters[END_REF]). The computation complexity per pixel is O(1), i.e., is independent of the SE size. Nonetheless, the algorithm requires two scans: forward and reverse. Lemonnier [START_REF] Lemonnier | Fast dilation by large 1D structuring elements[END_REF] proposes to identify local extrema and propagate their values as long as it is covered by the SE. Again, forward and reverse scans are needed. Notice that in 2-D the reverse scan of the vertical component multiplies the algorithm latency by a factor of the image width.

Lemire [START_REF] Lemire | Streaming maximum-minimum filter using no more than three comparisons per element[END_REF] proposes a fast, stream-processing algorithm for causal linear SE. It runs also on floating-point data, has low memory requirements and zero algorithm latency. However, an intermediate storage of local maxima results in a random access to the input data. This problem is solved in Dokladal and Dokladalova [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] using the strictly sequential access to the data. It allows the real on-the-fly computing and has zero algorithm latency.

A different approach represents the algorithm proposed by Buckley and Van Droogenbroeck [START_REF] Van Droogenbroeck | Morphological erosions and openings: Fast algorithms based on anchors[END_REF]. It detects the anchors -the portions of the signal unaffected by the operator -and updates only the parts to be modified by the operation. It has zero algorithm latency. However, the algorithm uses a histogram which makes memory requirements dependent on the number of gray levels.

Recently, Urbach and Wilkinson [START_REF] Urbach | Efficient 2-D grayscale morphological transformations with arbitrary flat structuring elements[END_REF] propose an algorithm for arbitrary shaped 2-D flat SEs based on the computation of multiple horizontal linear SEs for every pixel and storing them in a look-up table. The result is then computed by taking the maximum from the intermediate values (stored in the look-up table) corresponding to the shape of the SE. The horizontal linear SE can be computed with one of the above mentioned 1-D algorithms.

Implementations

In the beginning of the 70's, Klein and Serra [START_REF] Klein | The texture analyser[END_REF] propose a texture analyzer for linear and rectangular SE by decomposition based on the delay-line concept. More recently, Velten and Kummert [START_REF] Velten | Implementation of a high-performance hardware architecture for binary morphological image processing operations[END_REF] propose also a delay-line based architecture supporting arbitrary-shaped SEs. However, the complexity being quadratic O(W 2 ) (W denotes the length of the SE), it becomes penalizing for large SEs. In Chien et al. [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF], the authors show how to reduce the number of redundant comparisons within large SEs by merging adjacent smaller SEs. The complexity becomes O(⌈log 2 (W )⌉) with identical memory requirements.

Clienti et al. [START_REF] Ch | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF] proposes a highly parallel morphological System-on-Chip. It is based on a set of neighbourhood processors optimized for 3×3 SE, interconnected in a partially configurable pipeline. Larger SE is obtained by homothecy (see Basic Notions below) requiring to instantiate a deep pipe of these processors or multiple image scans.

A similar approach has been published by Deforges et al. [START_REF] Déforges | Fast recursive grayscale morphology operators: from the algorithm to the pipeline architecture[END_REF]. Based on Xu's [START_REF] Xu | Decomposition of convex polygonal morphological structuring elements into neighborhood subsets[END_REF] decomposition combined with a stream implementation, the authors propose a pipeline architecture composed of the elementary parametrizable blocks. It handles an arbitrary convex shape of structuring element in only one scan of the input image. However, using large SE will require the proportional increase of the atomic HW resources, concatenated in a deep pipe. The principal limitation comes from a limited programmability of such a pipe.

To complete this brief survey, we can also cite the systolic architectures proposed by Diamantaras and Kung [START_REF] Diamantaras | A linear systolic array for realtime morphological image processing[END_REF], Malamas et al. [START_REF] Malamas | Fast implementation of binary morphological operations on hardwareefficient systolic architectures[END_REF] or Shih et al. [START_REF] Shih | Pipeline architectures for recursive morphological operations[END_REF] for gray-scale or binary morphology. Their common inconvenience is the need of an intermediate storage for 2-D structuring element and a long response time of the system.

Novelty of this Paper

All previous algorithms optimize the dilation/erosion algorithm, rather than the entire operator chain. The performance will inevitably decrease for more complex applications with long loops (iterations, idempotence) or concatenations.

Consider some serial morphological filter ζ = δ ε . . . δ ε (with δ and ε standing for dilation and erosion, see Section 3 below for details). If the atomic operators δ and ε use sequential access to data then the entire ζ can run on pipelined, time-delayed data. If the atomic algorithms δ and ε -in addition -have zero algorithm latency, then the entire chain ζ inherits the same properties: sequential data access and zero algorithm latency. This is an interesting property, since computing ζ suddenly becomes very efficient: in stream, with only the (further irreducible) operator latency of ζ .

In comparison with the preceding state of the art, the Dokladal [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] algorithm extends the possibility to implement erosion/dilation filters with the arbitrarily large, 2-D SE in only one scan over the image, with the minimal algorithm latency and memory requirements. If implemented in a dedicated hardware, we can obtain the same features even for the implementation of long concatenations ζ .

This paper starts from the sequential HW implementation of the Dokladal algorithm (published in Bartovsky et al. [START_REF] Bartovský | Pipeline architecture for compound morphological operators[END_REF]). It describes more deeply the implementation features and optimization techniques. It shows how to exploit the inter-operator parallelism in ζ . Additionally, it introduces another intra-operator parallelism in the computation of the 2-D erosion/dilation. The 2-D erosion/dilation is implemented as a run-time programmable block. The operation (erosion or dilation), the size and the origin of the SE can be modified on-the-fly between two frames.

Several such blocks concatenated in a pipeline allow obtaining previously unachievable, real-time performances for operators in the form of ζ . We can reach almost 100Hz HDTV 1080p performance, independent of the length of ζ .

Basic principles

Let δ B , ε B : Z 2 → R be a dilation and an erosion on grey-scale images, parametrized by a structuring element B, assumed rectangular, flat (i.e., B ⊂ Z 2 ) and translation-invariant, defined as

δ B ( f ) = b∈B f b (1) ε B ( f ) = b∈ B f b (2)
The hat denotes the transposition of the structuring element, equal to the set reflection B = {x | -x ∈ B}, and f b denotes the translation of the function f by some vector b.

The SE B is equipped with an origin x ∈ B. Below, B(x) denotes B placed with its origin at x. The implementation of (1) and (2) consists of searching the extremum of f within the scope of B

[δ B ( f )](x) = max b∈B [ f (x -b)] (3) [ε B ( f )](x) = min b∈B [ f (x + b)] (4)
The dilation and erosion by convex structuring elements verify the homothecy. Let B be some convex structuring element, and rB the change of scale of B, with r > 1, r ∈ Z. Then for the dilation we have

δ rB f = δ B . . . δ B r-times f . (5)
The homothecy allows obtaining large-size dilations by repeating several times the dilation by a small SE.

Combinations of dilations and erosions form other operators. The basic concatenation products are opening γ B = δ B ε B and closing ϕ B = ε B δ B . From here we can form the Alternating Filters obtained as γϕ, ϕγ, γϕγ and ϕγϕ. The number of combinations obtained from two filters is rather limited. Other filters can be obtained by combining two families of filters. This leads to morphological Alternate Sequential Filters (ASF), originally proposed by Sternberg [START_REF] Sternberg | Grayscale morphology[END_REF], and studied in Serra [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], Chapter 10. In general, it is a family of operators parametrized by some λ ∈ Z + , obtained by the alternating concatenation of two families of increasing and decreasing filters {ξ i } and {ψ i }, respectively, such that

ψ n ≤ . . . ≤ ψ 1 ≤ ξ 1 ≤ . . . ≤ ξ n .
The most known ASF are those based on openings and closings, obtained by taking ψ = γ and ξ = ϕ :

ASF λ = γ λ ϕ λ . . . γ 1 ϕ 1 (6)
starting with a closing, and

ASF λ = ϕ λ γ λ . . . ϕ 1 γ 1 ( 7
)
starting with an opening.

The second application example is the size distribution of a population of objects [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF][START_REF] Sabourin | Off-line signature verification by local granulometric size distributions[END_REF][START_REF] Vincent | Granulometries and opening trees[END_REF]. One way to compute them is the residue from a sequence of openings

sd(λ ) = ||γ λ -γ λ -1 || (8)
The following section briefly recalls the principles of the used algorithm, [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF]. It starts by the 1-D dilation algorithm, followed by the principle of separation of the n-D dilation into perpendicular 1-D computations, preserving the stream aspects at all levels.

Algorithm Description

1-D Dilation Algorithm

The algorithm principles and properties have been originally published in [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF]. We briefly recall the main important principles for HW implementation.

For some 1-D input signal f : 1 . . . N → R, the algorithm 1 computes the value δ B f (wp) = f (rp). The SE B, B ⊂ Z, is a line segment, containing its origin, and not necessarily symmetric. Consequently, the size of B is given by the span from the centre to the left and to the right, SE1 and SE2. The length of B is SE1 + SE2 + 1. The coordinates wp and rp stand for the current writing and reading positions.

The algorithm uses a FIFO queue Q. The queue supports operations push, pop and dequeue (modifying the FIFO's content) and queries front and back. The input signal f is read sequentially. A newly read value f := f (rp) is inserted in the FIFO queue as a pair { f , rp}, the sample f and reading position rp (code line 3). In this pair, one can independently access either the value or the position by indexing. For example the last stored element's value can be accessed (without dequeing it) by a query Q.back() [START_REF] Bartovský | Pipeline architecture for compound morphological operators[END_REF].

The algorithm does not store non decreasing intervals (see [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] for details and proof). The values that appear to belong to increasing or constant intervals are dequeued (code lines 1-2). Consequently, the values stored in the queue are always ordered in a decreasing order.

The old values, uncovered by the SE, are retrieved from the queue (code lines 4-5). The result of the dilation δ B f (x) is read at the front of the queue (code line 7). The result becomes available as soon as enough input data have been read, otherwise the output is empty (code line 9).

2-D Dilation Algorithm

The separability of n-D morphological dilation into lower dimensions is a well known property. For example, a rectangular SE R decomposes as R = H ⊕V where H and V are horizontal and vertical segments and ⊕ is the Minkowski addition. Then the dilation by a rectangle R can be computed by concatenation of two perpendicular 1-D dilations

δ R = δ V δ H . (9)
The sequential access to the data in 1-D makes that two perpendicular 1-D computations can be assembled into 2-D with sequential access at both levels, 2-D and 1-D, for both input and output data. There is no additional latency and no intermediate storage (the data are pipelined).

See the example of dilation by a rectangle R=H ⊕ V of an N × M image f , Fig. 1. The image is sequentially read in the raster-scan mode, line by line from left to right. The various indices rp and wp denote reading and writing positions, respectively, for the segments H and V , and the rectangle R. The computation is illustrated for column i and line j, i.e., the result δ R f (i, j) is to be written at wp R .

The computation of δ R = δ V δ H decomposes as follows: The current reading position of δ R coincides with rp H , that is rp R = rp H . The result of the horizontal dilation, at wp H , is immediately read by the vertical dilation in the respective column, that is wp H = rp V . The result of the vertical dilation δ V is written at the writing position wp R , i.e., δ V = wp R .

Notes: -The rp r and wp r run over the image in the raster scan mode. The distance between rp R and wp R is the (further irreducible) operator latency.

-There is one instance of the horizontal dilation running at the current line j, and N instances of vertical dilation, i.e., one per each column.

Sequential Hardware Implementation

In this section, we firstly describe in details the implementation of the 1-D dilation in a basic block that (thanks to the separability) can be used as a building brick in any dimensional system. We illustrate this below on a 2-D dilation or erosion. Secondly, we show how the intra-operator parallelism can be introduced to increase the performance.

1-D Algorithm Implementation

The 1-D algorithm presented in Section 4 is a system with sequential behavior. It contains a while loop that can not be unrolled (uncertain number of iterations). The common way to implement such a system is the Mealy Finite-State Machine (FSM, see [START_REF] Mealy | A method for synthesizing sequential circuits[END_REF]). The FSM issues all the necessary operations over the memory as well as it controls the input and output data-flow. It consists of 2 states {S1, S2}. the while statement (code line 1). Consequently, its computation time varies from 1 to the SE size clock cycles in the worst case when all the previously stored pixels are unnecessary.

S1 S2 End Start Q.push({F, rp}); Q.pop(); return (Q.front()[1] ); return (Q.front()[1] ); Q.dequeue(); Q.back()[1] > F End of
The Enqueue current sample signal (code line 3) is issued upon the transition from S1 to S2.

The S2 state handles the code lines 4 and 5, Delete too old values, and the lines 6 to 9 Return valid value or Return empty. These instructions are independent and executed in parallel. Consequently, the execution of S2 takes only one clock cycle.

1-D Block Architecture

The HW implementation can be separated into 2 areas (Fig. 3), the FSM part and the memory part. The FSM manages entire computing procedure and temporarily stores values in the memory part. The memory instantiates one FIFO queue in the case of horizontal direction (horizontal scan) and N FIFO queues in the vertical case (N is the image width). The queues are addressed by a modulo N Page counter (active in the case of vertical direction).

The Control unit is a sequential circuit that manages the state transitions. It increments the rp, wp and manages the Page and position Stamp counter appropriately. The Control unit also performs the queue memory operations and handles the backward full flags used for data-flow control.

Principle

In the beginning of S1, the last queued pixel is invoked by Back() operation from the queue and fetched to the Comparator 1. The pixel value is compared with the value of the current sample. Notice that the comparator evaluates all three possible relations (>,<,=) at the time, for both dilation and erosion. The Control unit decides on the basis of comparison results and selected morphological function (dilation or erosion) whether the enqueued pixel is to be dequeued (lines 1-2). Otherwise, the current pixel is extended with the reading position stamp and enqueued (line 3).

The S2 invokes the oldest queued pair {pixel, stamp} by Front() operation. The read pixel is a correct result if rp has already reached or exceeded the SE2 parameter. This output allowing condition (line 6) is checked by Comparator 3. The deletion of outdated values is performed by comparing the current value of the reading-position stamp with the rp value of the oldest pair. Notice that the deletion has no impact on the output dilation value because Pop() operation (lines 4 and 5) issued by the Control unit has an effect only with the next clock edge.

The switch request logic is used only in the parallelized version of the architecture, see Section 6. It is a simple block containing several comparators which generate a signal with the last output value of each parallel segment. Its purpose is to inform the switch connected to the output that the end of the segment has been reached and the following segment is to be processed.

The entire set of parameters, i.e., SE dimensions and selection of the morphological function, is run-time programmable at the beginning of the line for 1-D, and of the frame for the 2-D implementation, respectively. In addition, no further controller is needed; the internal behavior is driven only by the regular scan order data-flow.

Reducing the impact of data-dependency

Hereafter, we briefly describe two techniques brought to the system for higher throughput and lesser area occupation.

Number of dequeue steps

The data-dependent number of dequeue steps (below denoted by Steps) has an unpleasant consequence on the HW design: longer balancing FIFOs (see Fig. 4), lower data throughput. For HW design it is important to minimize the worst case upper bound Steps orig =SE -1.

The number of stored pixels is within [1, SE[. Suppose that we are to dequeue D pixels. We know that the pixels are queued in a strictly decreasing order. Thus, if the DL-th pixel (DL< D) can be dequeued then also all previous pixels can be dequeued. This can be done at the same time. Hence, the worst-case number of dequeue steps is

Steps = max D<SE (D div DL + D mod DL) ( 10
)
where D denotes the number of pixels to be dequeued and div and mod the integer division and the remainder operations. D can be regarded as a uniformly distributed, random variable D ∈ [1, SE[. Then we need to find the optimal DL that minimizes Steps (Eq. 10) for all D such as

DL optim = arg min DL<D max D<SE (D div DL + D mod DL) (11)
The optimal DL optim brings us the minimal number of dequeue steps Steps optim

Steps optim = min

DL<D max D<SE (D div DL + D mod DL) ( 12
)
Table 1 exemplifies, for some SE widths, the original and reduced number of dequeue Steps, obtained with optimal DL. Notice that more than one optimal DL can exist.

The SE is user programmable. DL optim also is programmable, though it is useless to make it accessible to the user; it can instead be read from a LUT for every given user-specified SE.

Pixel addressing

The absolute pixel addressing in the queues can be advantageously replaced in the HW by using the modulo addressing. Instead of the absolute reading position rp, we use the relative modulo position stamp = rp mod SE. The pixels are enqueued by Q.push( f , stamp) (code line 3). The delete condition of line 4 changes accordingly. Using the modulo addressing, a stored pixel becomes outdated whenever its modulo address equals the current pixels' one (stamp = Q. f ront() [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF]).

The advantage of the modulo addressing is a smaller data width. It fits into ⌈log 2 (SE -1)⌉ bits, whereas the absolute addressing requires ⌈log 2 (N -1)⌉ bits. This is mainly advantageous for vertical orientation using N queues for a unit.

2-D Dilation Implementation

Recall that dilation is separable into lower dimensions, Eq. 9. The dilation by a rectangle can be implemented using two 1-D dilation blocks, Fig. 4. The computing latency of the dilation varies per each pixel. In order to preserve the input/output stream flow, one needs to compensate the different latencies by insertion of balancing FIFOs. The FIFO fills when the preceding block outputs data faster than the subsequent block can read. The depth of this FIFO directly defines the upper bound of the system latency of the 2-D block.

Obviously, the FIFOs should be as small as possible. The necessary depth infers from the dequeuing worst case

F input = Steps hor + 2 StreamRate -1 (13)
F balance = N Steps ver + 2 StreamRate -1 (14)
where Steps hor and Steps ver are numbers of the dequeue steps in horizontal and vertical directions [START_REF] Lemonnier | Fast dilation by large 1D structuring elements[END_REF]. The output FIFO ensures a permanent stream delay in all circumstances. Its maximal size is a sum of both FIFOs (input and balancing). The instantaneous filling of output FIFO is complementary to the filling of both FIFOs combined. The overall delay does not change. If more 2-D blocks are pipelined to form compound operators (e.g., opening, closing, ASF), only one output FIFO at the end is necessary. The output and balancing FIFOs can be merged (see Fig. 5) into one memory thanks to the following properties: 1) the vertical unit reads exactly one pixel from the balancing FIFO for each pixel written to the output FIFO. Consequently, filling of these two FIFOs is complementary; the occupied memory spaces can not collide with each other, 2) the read/write activity is at most 1 access per 2 clock cycles. Hence, reading ports of both FIFOs can use one memory port and the writing ports can use the other memory port (without overloading it). Merging both FIFOs reduces the memory to approximately one half. The result memory (see Fig. 5) has two pairs of standard FIFO ports, but it contains only one dual-port RAM.

Clock rate

The overall average clock rate stays in the interval from 2 clock cycles per pixel in the best case, up to 3 clock cycles per pixel in the worst case. The current rate between 2 and 3 clock cycles per pixel is data dependent.

A temporarily worst case arrives whenever a monotonously decreasing signal is followed by a high value. This makes a number of samples to be dequeued at the time (code lines 1 -2, and the S1 state of the FSM), and the computing latency temporarily increases. However, the average computing latency remains unchanged, compensated by the fact that during the entire monotonous decrease of the signal no values have been dequeued. Therefore, the average clock cycles per pixel rate remains constant.

Memory requirements

The memory requirements of the 2-D architecture consist of horizontal and vertical computation-involved memories and two balancing FIFOs, defined by ( 13) and ( 14).

In the vertical case, the algorithm uses a several queues. Instantiating N separated memories would be resource inefficient because the FPGA RAM blocks could not be exploited. Instead, these queues are gathered in a single dualport memory (see Fig. 6) since only one queue is accessed at the time (the others are idle). A single memory block also allows using an off-chip memory. Every queue has a related pair of front and back pointers which must be retained throughout the entire computation process. The appropriate pair is always read before the particular queue is used and the modified pointers are stored back after the computation left the queue. These pointers are stored in a separated pointer memory. The queues are efficiently packed into RAM blocks resulting in a small memory extension.

Let W×H denote the width×height of the rectangular SE, and bpp bits per pixel. The memory contribution per 2-D unit is given by: M hor =W (bpp + ⌈log 2 (W -1)⌉)

[bits] (15)

M ver =N(H(bpp + ⌈log 2 (H -1)⌉)+ + 2⌈log 2 (H -1)⌉) [bits] (16)
The following example illustrates the very low memory consumption achieved thanks to the stream processing. Neither the input, nor the output or any working image are buffered.

Example: Consider a dilation of 8bpp, SVGA image (i.e., 800×600=N×M) by a square, 31x31 SE.

The computation (the queues) requires ( 15) and ( 16) The input and the balancing FIFOs require ( 13) and ( 14)

M
F input + F balance = (N + 1)( Steps + 2 StreamRate -1)8bpp = = (800 + 1)( 9 + 2 3 -1)8bpp ≈ 17 kbits
The total memory needs to implement the 2-D dilation are 331+17=349 kbits. This is far below the mere size of the image itself 800 × 600 × 8bpp= 3.84 Mbits which, at any moment, does not need to be stored.

Parallel Hardware Implementation

This section develops and implements the concept of the previously mentioned intra-operator parallelism in the dilation/erosion operator. Its main objective is to increase the throughput while maintaining the beneficial properties of the proposed algorithm, namely the sequential data access and minimal algorithm latency as much as possible.

The principle is based on utilization of concurrently working units that process different parts of the image simultaneously. The number of units used in parallel for horizontal and vertical directions defines the parallelism degree (PD). Considering that the input data are fetched line by line, we propose a solution minimizing the waiting-for-data periods of all units.

The image partition for 2-D dilation conforms to the intersection of two horizontal and vertical partitions (Fig. 7). Its granularity is determined by the PD. The horizontal partition (partition of image among horizontal units) is interleaved, whereas the vertical units use the partition into compact blocks. During the parallel processing the computation runs simultaneously at multiple segments of the image, see Fig. 8. These segments must belong to different columns and lines,i.e., must be placed on a diagonal. The input data rate can be theoretically PD-times faster than the computational throughput of one unit. Therefore, each image line needs to be buffered in a line buffer. The line buffers are filled at the external (fast) pixel rate and read by the internal PD-times slower rate. Figure 8 gives an example for PD = 3. We have three horizontal (H1 .. H3), and tree vertical (V1 .. V3) processing units. As soon as the line buffer receives the first pixel, the first horizontal unit H1 starts the processing and feeds results to the first vertical unit V1. Its output is fed to the first output line, see Fig. 8(a). After N received pixels, the output of H1 is connected to V2 which belongs to output line 1. Since the H1 left V1 and line 2 is read, the H2 can start processing second line feeding V1 connected to output line 2, see Fig. 8(b). When the 2N input pixel is received, the H1 connects to V3, H2 connects to V2 and H3 connects to V1, see Fig. 8(c), and so on.

Architecture

The parallel architecture depicted in Fig. 9 contains four separable generic parts scalable by n ≡ PD: input buffer, horizontal and vertical parts and output buffer. The input buffer is mainly composed of the 1-to-n multiplexer and n line buffers (we omitted the control logic). It divides the fast input stream into n (n-times slower) streams processed by computational units as described above. The output buffer composes n slow streams of the processed data into a single, fast, output stream respecting the image horizontal scan order. The operator blocks can be concatenated into more complex functions (opening, closing, ASF, etc.). The buffers are used only at the beginning and at the end of the chain. Both horizontal and vertical parts instantiate n balancing FIFOs, n horizontal or vertical units, and one switch that manages the interconnection. Each horizontal unit along with the front-end FIFO conforms to Section 5.2.

The width of the processing area proportionally affects both vertical memories, see [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF] and [START_REF]Mathematical Morphology: From Theory to Applications[END_REF]. The area of every horizontal unit remains unchanged, since every unit processes the entire line. The overall memory of the horizontal part is a factor of n. Contrarily, the memory requirements of every vertical part is divided by n because it processes only a fraction of the original image width. The area of the FSM of vertical units increases linearly with n.

Switching

The routing of the computation units is handled by a switch block. Every switch contains n input ports from previous units and the same number of output ports linked to the subsequent units. The purpose of the switch is to manage up to n interconnection channels. Notice that they are bidirectional: forward data and backward FIFO full flag. As described in Fig. 8, the output switching of all input ports is circular, i.e., V1, V2 ... Vn, V1, V2, ... and so forth. This property makes the switching easier because the only condition to evaluate is when to switch and whether the requested output unit is available.

The moment when to switch a given port is provided by the preceding unit's Switch Request logic. It generates a request every time it crosses the border of adjacent segments. If the desired unit is free, the switch reconnects the channel. If not, the switch sets high the FIFO full flag of requesting unit to stall it until the desired destination unit is freed and the channel can be established. All the channels are switched independently so stalling one unit does not affect the others. Figure 10 depicts the basic unit of the switch for one pair of input/output ports referred to as A. For n pairs of ports this circuitry is instantiated n-times. Each input port possesses a related control unit block that manages all channel transitions considering the availability of the requested partition. If this is still occupied, the requesting computation unit is stalled by holding its FIFO full flag active.

Experimental results

The proposed 2-D stream processing architectures have been implemented in VHDL, and targeted to the Xilinx Virtex5 FPGA (XC5VSX95T-2) using the XST synthesis tool. The processing clock frequency is 100 MHz. Notice that the queues are gathered in a block RAM memory, and thus its access time augments the critical path delay.

The measured performance for non-parallel architectures (PD=1) in terms of overall latency, clock cycles per pixel and FPGA area are given by Tables 2 and3. One can observe that the overall latency is factor of the SE size, the image width (both caused by operator latency) and the pixel rate (computing latency). The average pixel rate (AR) remains constant (Table 2). The average pixel rate can be expressed by [START_REF] Sabourin | Off-line signature verification by local granulometric size distributions[END_REF] and the stream frame-per-second (FPS) ratio by [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. T proc is overall time consumed by processing and f clk = 100 MHz is clock frequency of computation units.

AR = T proc -2SE 2 (N + M + SE 2 )/PD N M [clk/px] (17) FPS = f clk PD AR N M + 2SE 2 (N + M + SE 2 ) [fr/s] (18)
M, N denote the width and height of the image, SE 2 denotes the width of the structuring element from the origin rightwards.

Concerning the area occupation (see the Xilinx documentation [START_REF]Xilinx. Virtex-5 family documentation[END_REF]), the number of registers is quasi-constant; the number of LUTs and BRAM blocks increases linearly with the SE and image sizes (Table 3). Although the vertical memory (size is given by ( 16)) is packed into the RAM block, the amount of the used memory always exceeds the theoretical value. It is caused by a different memory organizations; e.g., the required word is 13 bits whereas available memories are of width 36 bits and its fractions.

The experimental frames-per-second (FPS) rate is obtained on a natural test image (see Fig. 11). The worst-case FPS is a theoretical worst-case performance of the system expected on the synthetic saw-shaped data.

Table 4 presents the relative speed-up of the parallel architecture vs. the intra-operator parallelism PD. In terms of overall latency and average processing rate, the processing domain clock cycle is considered as a reference unit. Note The FPGA area results, Table 5, are separated into 2 groups: the area of computing parts and buffers. The area of input and output buffers is linear w.r.t. both N and PD since their essential components are PD line buffers (FIFO memories with independent ports of N elements). The area of the operator units in terms of Slice registers and LUTs is proportional to PD as well because n independent circuits are instantiated in a parallel manner. Although the overall vertical memory requirements remain unaffected by PD, practically the number of occupied RAM blocks slightly increases. It is caused by a different memory organization. The ultimate timing results (PD=6) versus the image size are listed in Table 6. It illustrates the real performance of the architecture. It allows to achieve at least 96 fps with 1080p image size (full HD TV image size).

The worst case occurs on artificial saw-shaped image with no constant plateaus. Such an image infers the maximal number of algorithm's while-loop iterations. The best case fps (not mentioned in the table) is obtained with a constant image. A real, unfiltered image containing textures or random noise achieves performance somewhere between best and worst cases. For instance at 1080p, the worst case is 96 fps, best case 140 fps, achieved experimental performance is 113 fps.

This frame rate remains constant for any morphological serial filter (such as ASF). Obviously, the FPGA area increases accordingly to the size of the ASF. The implementation is eased by the fact that one can use an off-chip memory.

Comparison with existing HW implementations

Table 7 presents a comparison other recent architectures. The table is divided into three sections. The processing unit section presents the features of a single 2D computational unit. The second part the HW specifications, and the third part the performance on a given application, an ASF filter.

One can see that Clienti [START_REF] Ch | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF] offers a high throughput for small 3×3 rectangular SEs. Similarly, the Chien ASIC chip [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF] provides very reasonable performance on small SEs. On the other hand, Déforges [START_REF] Déforges | Fast recursive grayscale morphology operators: from the algorithm to the pipeline architecture[END_REF] directly offers large, nonrectangular, convex SE, but with a lower processing rate. The programmability is not mentioned, namely, the possibility to control the SE shape after the synthesis is not clear.

Although all these solutions are efficient for small SE sizes or short concatenations, they become more or less penalized for longer filters. This issue is illustrated in an Example Application, Table 7. It estimates the performances on a five-stage ASF 5 = ϕ 11×11 γ 11×11 . . . ϕ 3×3 γ 3×3 . Decomposed into a sequence of dilations and erosions, it can be realized as ASF 5 = ε 11×11 δ 21×21 . . . ε 5×5 δ 3×3 . Notice that it makes use of a progressively increasing SE. On neighborhood processors, large SE can be obtained using the homothecy Eq. 5. The Clienti SPOC instantiates 16 of 3×3 processing units. Hence, the ASF 5 will require 5 image scans with the entire image necessarily buffered in the memory. Chien also uses the homothecy. This deteriorates the throughput.

One could immediately figure out to instantiate a longer pipe in order to reduce the number of image scans. Alas, a long, fixed-length pipe lacks the flexibility. Consider another application for the illustration of the problem: the size distributions, exemplified by Fig. 12. Contrarily to ASF, the size distributions are often sampled sparsely, the SE increments by more than one and, at the same time, one often goes to much larger SE sizes. Every opening {γ B i } in (8) needs to be output and stored in the memory to compute the subtraction. For small sizes, a long pipeline is underused and the workload of the processing units unbalanced, whereas for large λ one may still need several image scans.

For example, for sizes λ = 5, 10, 15, 20, 25, as in Fig. 12, the Clienti SPOC will require 7 image scans. The 16 processor pipe is underused for λ = 5, 10, 15, whereas it will require 2 scans for λ = 20, 25.

Our processing unit with programmable SE size avoids using the homothecy. This allows optimal workload distribution over the entire pipe, so important for processing large images in real-time systems.

Conclusions

This paper describes an efficient implementation of serial morphological filters with flat, rectangular structuring elements of arbitrary size. The efficiency is obtained through the following properties: -The computational complexity is linear w.r.t. the image size and independent of the SE size.

-The overall latency is mostly equal to the latency of the operator, inferred by the size of the used structuring element. -It uses strictly sequential access to the data at all algorithm levels.

-Low memory consumptions (far below the size of the image) allow embedding on a single chip complex operators able to process large images.

-Two levels of parallelism: i) the inter-operator parallelism in serial concatenations ζ = δ ε . . . δ ε, allow running all these atomic δ and ε operators simultaneously, and ii) the intra-operator parallelism in every atomic dilation/erosion. The intra-operator parallelism is scalable (tested up to six) and allows the decomposition of fast streams into several slower streams processed in parallel without altering the streaming property of the system.

The architecture serves as a basic building block to be used for construction of more complex operators such as ASF, granulometries, etc., with the same properties and performance. The performances obtained on an FPGA are approaching the 100Hz HDTV 1080p standard. These performances are far above what has been reported in the literature. These performances allied to the programmability are extremely interesting. They open the accessibility of advanced morphological operators in industrial systems running under severe time constraints. The number of examples includes the on-line production control, aging material defectoscopy, etc., wherever one requires processing of high resolution images and low latency.

Fig. 1

 1 Fig. 1 Decomposition of dilation by a rectangle R into two 1-D dilations by segments H and V , see (9).

Fig. 2 Fig. 3

 23 Fig. 2 State diagram of the 1-D Algorithm FSM. State transition conditions are typeset in bold; the output signals are given in a shadow bounding box.

Fig. 4 2

 4 Fig. 4 2-D implementation is composed of 1-D blocks for respective directions.

Fig. 5

 5 Fig. 5 Merged FIFO replaces the balancing and output FIFOs to reduce memory requirements.

hor = 31 ( 8 + 5 )

 3185 = 403 bits and M ver = 800(31(8 + 5) + 2 × 5) = 330.4 kbits resulting in a total of 331 kbits for the 2-D dilation.

Fig. 7

 7 Fig. 7 Example of image partition for PD=3: image is divided horizontally line by line and into PD equal stripes in a vertical direction. The final image partition is obtained by intersection.

Fig. 8

 8 Fig. 8 Image partitioning and switch routing in parallel processing for PD=3. Decomposed in time -(a) beginning of processing, (b .. d) after kN pixels, k=1..3. The shading denotes the state : Dark Gray -being computed, Light Gray -already computed, White -waiting.

Fig. 9

 9 Fig.[START_REF] Gil | Computing 2-d min, median, and max filters[END_REF] Overview parallel 2-D architecture. The horizontal and vertical stages can be instantiated several times between input/output buffers to create compound operators.

Fig. 10

 10 Fig.[START_REF] Klein | The texture analyser[END_REF] Basic unit of the switch. Every switch contains n basic units for a correct routing between n input/output ports.

Fig. 11

 11 Fig. 11 (a) Experimental 800×600 lotus image. (b) Zoom on the fine veinous texture disadvantageous for the algorithm. (c) Result of ASF 4 filter, see Eq. 6 or 7. (d) Zoom on the ASF 4 filtered image.

→Fig. 12

 12 Fig.[START_REF] Lemonnier | Fast dilation by large 1D structuring elements[END_REF] The size distribution of the texture grain.

Appendix: The 1 -D Dilation Pseudocode Algorithm 1 : 2 Q 4 Q

 1124 df←1D DILATION (rp, wp, f, SE1, SE2, N) Input: rp, wp -reading/writing position; f -input signal value f (rp); SE1, SE2 -SE size towards left and right; N -length of the signal; Q -a FIFO-like queue Result: output signal value δ B f (wp) while Q.back()[1] ≤ f do 1 Q.dequeue() ; // Dequeue useless valuess .push({f, rp}) ; // Enqueue the current sample 3 if wp -SE1 > Q.front()[2] then

Table 1

 1 Optimal dequeue length DL, original and reduced number of dequeue steps for selected SE widths

	SE width	3	11	21	31	41
	Steps orig	2	10	20	30	40
	DL optim	2 3, 4 4, 5, 6 5, 6, 7 6, 7
	Steps optim 2	4	7	9	10

  Vertical Queues are mapped into linear memory space side by side. The front and back pointers are stored at separated memory.

	Address:	0	H-1	2*H-1	((N-1)*H)-1	N*H-1
		Queue 1	Queue 2	. . .	Queue N
	Q1.back	Q2.back		QN.back
	Q1.front	Q2.front	QN.front
	Fig. 6

Table 2

 2 Timing and area vs. SE, SVGA image size, PD=1.

	Size of SE (sq.)	3x3	11x11 21x21 31x31 41x41
	Latency [clk]	1908	9474	18888 28351 37969
	Av. rate [clk/px] 2.344	2.356	2.360	2.361	2.361
	Registers	212	232	242	242	252
	LUTs	584	761	859	859	953
	Block RAMs	2	6	13	13	28

Table 3

 3 Timing, frame rate and area w.r.t. image, SE = 31x31 square, PD=1.

	Size of Image	CIF	VGA	SVGA	XGA	1080p
	Latency [clk]	12826 23465 28351 37472 69548
	Av. rate [clk/px]	2.371	2.376	2.361	2.383	2.368
	Experimental FPS	384	130	85	51	20.5
	Worst-case FPS	319	106	68	41	16
	Registers	231	237	242	242	253
	LUTs	761	853	859	859	1057
	Block RAMs	7	13	13	13	26

Table 4

 4 Timing vs. degree of intra-operator parallelism PD. SVGA image, SE = 31x31 square.

	PD	2	3	4	5	6
	Latency [clk]	14243	9561	7244	5818	4893
	Av. rate [clk/px]	1.220	0.824 0.625 0.505 0.426
	Exp. speed up	1.938	2.869 3.785 4.682 5.554

Table 5

 5 Area vs. degree of intra-operator parallelism PD. SVGA image, SE = 31x31 square.

	PD	2	3	4	5	6
	Registers	650	978	1280 1605 1938
	LUTs	2138 3227 3862 4875 6054
	Block RAMs	13	14	14	18	21
	Reg. buf	661	969	1279 1587 1896
	LUTs buf	1408 2086 2776 3459 4135
	that the latencies of parallel versions are merely fractions
	(divided by PD) of non-parallel values.

Table 6

 6 Timing and frame rate vs. image size, PD = 6, SE = 31x31 square

	Size of Image	CIF VGA SVGA XGA SXGA 1080p
	Latency [clk]	2208 3996 4893 6390 7391 11641
	Av. rate [clk/px] 0.443 0.431 0.426 0.426 0.427 0.418
	Experimental FPS 2075 724	472	290	174	113
	Worst-Case FPS	1915 640	411	246	151	96

Table 7

 7 Comparison of several FPGA and ASIC architectures concerning morphological dilation and erosion Number of available stages varies with size of used FPGA

			Processing unit			HW System	Example Application ASF 5
		Parallel Supported Throughput	f max Clock rate Number of Supported Image	FPS
		degree	SE	[Mpx/s]	[MHz] [clk/px]	units	image	scans
	Clienti [3]	4	arb. 3x3	403	100	0.25	16 *	1024x1024	5	80
	Chien [2]	1	disk 5x5	190	200	1.052	1	720x480	27	21.5
	Déforges [5]	1	arb. convex	50	50	1	1 *	512x512	11	17.2
	This paper	6	rectangles	234	100	0.426	11 *	1920x1080	1	113

*

See Appendix for the 1-D dilation pseudocode