Philippe Angot
email: angot@cmi.univ-mrs.fr

Rima Cheaytou
email: rima.cheaytou@gmail.com

VECTOR PENALTY-PROJECTION METHOD FOR INCOMPRESSIBLE FLUID FLOWS WITH OPEN BOUNDARY CONDITIONS

Keywords: Vector penalty-projection methods, Stokes equations, incompressible flows, open boundary conditions, second-order convergence AMS subject classifications. 35Q30, 35Q35, 49M30, 65M08, 65M85, 65N08, 65N85, 74F10, 76D05, 76D45, 76M12, 76R10, 76T10

 to compute the solution of unsteady incompressible fluid flows and to overcome most of the drawbacks of the usual incremental projection methods. In this work, we deal with the time-dependent incompressible Stokes equations with outflow boundary conditions using the present method. The spatial discretization is based on the finite volume scheme on a MAC staggered grid and the time discretization is based on the backward difference formula of second-order BDF2 (namely also Gear's scheme). We show that the (V P Pr,ε) method provides a second-order convergence for both velocity and pressure in space and time even in the presence of open boundary conditions with small values of the augmentation parameter r; typically 0 ≤ r ≤ 1 and a penalty parameter ε small enough; typically ε = 10 -10 . The resulting constraint on the discrete divergence of velocity is not exactly equal to zero but is satisfied approximately as O(ε δt) with a penalty parameter ε as small as desired. Finally, the efficiency and the second-order accuracy of our method are illustrated by several numerical cases.

1. Introduction. The numerical solution of incompressible fluid flows has always been one of important subjects in fluid dynamics. The major difficulty in numerically solving unsteady incompressible Navier-Stokes equations in primitive variable form arises from the fact that the velocity and the pressure are coupled by the incompressibility constraint at each time step. There are numerous ways to discretize these equations. The most widespread methods are projection methods like pressurecorrection methods. This family of methods has been introduced by Chorin and Temam [START_REF] Chorin | Numerical solution of the Navier-Stokes equationss[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II)[END_REF] in the late sixties. The interest of pressure-correction projection methods is that the velocity and the pressure are computed separately. Actually, in a first step, we solve the momentum balance equation to obtain an intermediate velocity and then, this intermediate velocity is projected on the space of solenoidal vector fields. Using this method, a numerical error named the splitting error appears and several papers have been written to estimate this error. In 1992, Shen [START_REF] Shen | On error estimates of some higher order projection and penalty-projection methods for Navier-Stokes equations[END_REF] has introduced a modified approach which consists of adding in the first step of the scheme a penalty term built from the divergence constraint, it is the so-called penalty-projection method. In the same way, Jobelin et al. [START_REF] Jobelin | A finite element penalty-projection method for incompressible flows[END_REF] proposed a numerical scheme which falls in the category of the penalty-projection method with an augmentation parameter r ≥ 0 independent of the time step δt. Recently, Angot et al. [START_REF] Ph | Vector penalty-projection methods for the solution of unsteady incompressible flows, Finite Volumes for Complex Applications V[END_REF] propose a new family of methods, the so-called two-parameters vector penalty-projection method (V P P r,ε ) r≥0,0<ε≤1 (r and ε represent the augmentation parameter and the penalty parameter respectively). Angot et al. have tested in [START_REF] Ph | Vector penalty-projection methods for the solution of unsteady incompressible flows, Finite Volumes for Complex Applications V[END_REF][START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] the (V P P r,ε ) method using a first-order scheme in time and in presence of Dirichlet boundary conditions. They found that the scheme is O(δt) in time for velocity and pressure. However, many applications such as free surface problems and channel flows have to deal with open boundary conditions on a part of the boundary which is still the challenging problem. Many recent progress were made in this direction including Guermond et al. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], Jobelin et al. [START_REF] Jobelin | A finite element penalty-projection method for incompressible flows[END_REF], Fevriére et al. [START_REF] Févrière | On the penalty-projection method for the Navier-Stokes equations with the MAC mesh[END_REF], Liu [START_REF] Liu | Open and traction boundary conditions for the incompressible Navier-Stokes equations[END_REF] and Poux et al. [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF].

The objective of the present work is to show the optimal convergence of the vector-penalty projection methods when solving the unsteady incompressible Stokes problems including open boundary conditions.

ρ ∂v ∂t -µ ∆v + ∇p = f in Ω×]0, T [ (2.1) ∇ • v = 0 in Ω×]0, T [ (2.2) v = v D on Γ D ×]0, T [ (2.3) -p n + µ ∇v • n = g on Γ N ×]0, T [ (2.4)
where v=(u, v) T stands for the fluid velocity of initial value v(t = 0) = v 0 , p for the pressure field, ρ for the fluid density (the density is taken to be one), µ for the dynamic viscosity (here, µ = 1 Re with Re a given Reynolds number). We impose Dirichlet boundary conditions (2.3) on Γ D , and open boundary condition (2.4) on Γ N . The external body forces f, the pseudo-stress vector g and the Dirichlet boundary condition v D are known.

Finally, the reader will keep in mind that bold letters like v, g, etc., indicate vector valued quantities.

Vector penalty-projection methods with open boundary condition.

Let 0 = t 0 < t 1 < ... < t N = T be a partition of the time interval of computation [0,T], which we suppose uniform for the sake of simplicity. We denote by δt the time step, i.e., the constant difference between two successive times t n and t n+1 .

3.1. Description of the (V P P r,ε ) method with OBC (3.6). We describe now, the two-step vector penalty-projection method (V P P r,ε ) using OBC (3.6) as an open boundary condition in the projection step with an augmentation parameter r ≥ 0 and a penalty parameter 0 < ε 1.

3.1.1. BDF2 time scheme. For all n ≥ 1 such that (n + 1)δt ≤ T and for v 0 , v 1 , v 0 , v 1 ∈ L 2 (Ω) and p 0 , p 1 ∈ L 2 0 (Ω) given, find (v n+1 , p n+1 ) such that: • Vector penalty-prediction step with an augmentation parameter r ≥ 0:

3 v n+1 -4 v n + v n-1 2δt -µ ∆ v n+1 -r ∇(∇ • v n+1 ) + ∇p ,n+1 = f n+1 in Ω (3.1) v n+1 = v n+1 D on Γ D (3.2) (-p ,n+1 + r ∇ • v n+1 ) n + µ ∇ v n+1 • n = g n+1 on Γ N (3.3)
where p ,n+1 = 2p n -p n-1 ( by a second-order Richardson's extrapolation)

• Vector penalty-projection step with a penalty parameter 0 < ε ≤ 1:

3 v n+1 -4 v n + v n-1 2δt -ε µ ∆ v n+1 - 1 ε ∇ ∇ • v n+1 = 1 ε ∇ ∇ • v n+1 in Ω(3.4) v n+1 = 0 on Γ D (3.5) µ ∇ v n+1 • n = 0 on Γ N (OBC) (3.6)
• Correction step for velocity and pressure:

v n+1 = v n+1 + v n+1 (3.7) p n+1 = 2p n -p n-1 - 1 ε (∇ • v n+1 ) -r ∇ • v n+1 (3.8)
Remark 1. Summing the prediction and projection steps gives the discrete problem which is effectively solved by the above splitting scheme:

3v n+1 -4v n + v n-1 2δt -µ (∆ v n+1 + ε ∆ v n+1 ) + ∇p n+1 = f n+1 in Ω×]0, T [ (εδt) p n+1 -p ,n+1 δt + ∇ • v n+1 + r ε ∇ • v n+1 = 0 in Ω×]0, T [ v n+1 = v n+1 D on Γ D ×]0, T [ (-p ,n+1 + r ∇ • v n+1 ) n + µ ∇v n+1 • n = g n+1 on Γ N ×]0, T [
As a consequence of the previous system, one would note that the (V P P r,ε ) method lost a little consistency. Nevertheless, the method can be ultra-fast and very cheap if η = ε/δt is sufficiently small due to the adapted right hand-side in the correction step which lies in the range of the left hand side as ε is taken small enough. This crucial property was already shown in [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF]Theorem 1.1 and Corollary 1.3] and in [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]Theorem 3.1] and numerically confirmed also in [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF][START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF].

3.1.2. Correction of pressure gradient for r = 0 and small values of ε. The augmentation parameter r in our method is kept constant and within small values (0 ≤ r ≤ 1) to avoid to excessively degrade the conditioning of the linear system associated to the prediction step. However, taking r = 0, we obtain a poor convergence convergence rate in time for velocity and pressure: see the numerical results in Section 4. To improve the convergence rate, we reconstruct the pressure field itself very fast from its gradient as already proposed by Angot et al. in [START_REF] Ph | A fast vector penalty-projection method for incompressible Darcy and multiphase Navier-Stokes problems[END_REF][START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] where they observed that it is numerically far better to update directly the pressure gradient to avoid the effect of round-off errors when ε is very small.

In this regard, the updating of the pressure in our method is as follows:

Taking the gradient of the equation (3.8), we get:

∇p n+1 = 2∇p n -∇p n-1 - 1 ε ∇(∇ • v n+1 ) -r ∇(∇ • v n+1 ) (3.9)
Moreover, we have the following equality in the penalty-projection step:

3 v n+1 -4 v n + v n-1 2δt -ε µ ∆ v n+1 = 1 ε ∇(∇ • v n+1 ) Now, replacing the term , 1 ε ∇(∇ • v n+1
) in (3.9) by the term deduced from the above equality. It yields the estimation of the gradient of the pressure:

∇p n+1 = 2∇p n -∇p n-1 - 3 v n+1 -4 v n + v n-1 2δt +ε µ ∆ v n+1 -r ∇(∇• v n+1 ) with r ≥ 0 (3.10)
3.2. Description of (V P P r,ε ) method with OBC1 (3.16). The (V P P r,ε ) method with OBC (3.6) yields good numerical results (see Section 4). However, we observe that the well-posedness of the penalty-projection step using OBC (3.6) is not straightforward. Thus, we propose to replace OBC (3.6) by OBC1 (3.16) which clearly yields a well-posed penalty-projection step (see Lemma 3.1).

BDF2 time scheme. For all

n ≥ 1 such that (n + 1)δt ≤ T and for v 0 , v 1 ,v 0 , v 1 ∈ L 2 (Ω) and p 0 , p 1 ∈ L 2 0 (Ω) given, find (v n+1 , p n+1
) such that: • Vector penalty-prediction step with an augmentation parameter r ≥ 0:

3 v n+1 -4 v n + v n-1 2δt -µ∆ v n+1 -r∇(∇ • v n+1 ) + ∇p ,n+1 = f n+1 in Ω (3.11) v n+1 = v D on Γ D (3.12) (-p ,n+1 + r∇ • v n+1 )n + µ∇ v n+1 • n = g n+1 on Γ N (3.13)
where p ,n+1 = 2p n -p n-1 (by a second-order Richardson's extrapolation)

• Vector penalty-projection step with a penalty parameter 0 < ε ≤ 1:

3 v n+1 -4 v n + v n-1 2δt -εµ∆ v n+1 - 1 ε ∇ ∇ • v n+1 = 1 ε ∇ ∇ • v n+1 in Ω (3.14) v n+1 = 0 on Γ D (3.15) v n+1 • n = 0 and (µ∇ v n+1 • n) ∧ n = 0 (OBC1) on Γ N (3.16)
• Correction step for velocity and reconstruction of the pressure from its gradient:

v n+1 = v n+1 + v n+1 ∇p n+1 = 2∇p n -∇p n-1 - 3 v n+1 -4 v n + v n-1 2δt + ε µ∆ v n+1 -r∇(∇ • v n+1 ) Remark 2. For any vector u ∈ R d defined on Γ, the tangential component u n| Γ is defined by: u n| Γ = u τ n| Γ with u τ = u -(u • n)n. Thus, in 2-D (d = 2
), we have simply: u n| Γ = u • τ , where τ denotes a unit tangential vector on Γ.

Lemma 3.1. (Well-posedness of the velocity correction step (3.14-3.16)) For all v n+1 given in H 1 (Ω) and for all ε > 0 and δt > 0, there exists at each time step a unique solution v n+1 ∈ W(Ω) to the velocity-correction step (3.14-3.16), where

W(Ω) = u ∈ H 1 (Ω) d ; u = 0 on Γ D and u • n = 0 on Γ N is an Hilbert space.
Sketch of Proof. For sake of simplicity, we take the discrete form of projection step and we use Euler scheme to discretize in time. Multiplying the projection step by a test function ϕ ∈ W(Ω) and applying the Green formula, we obtain the following bilinear form:

a( v, ϕ) = ε δt ( v, ϕ) 0 + ε 2 µ(∇ v, ∇ϕ) 0 + (∇ • v, ∇ • ϕ) 0 in W(Ω) × W(Ω). We note that a( v, ϕ) is a continuous and coercive form in W(Ω) × W(Ω) and L(ϕ) := (∇ • v, ∇ • ϕ) 0 is a linear continuous form in W(Ω).
Under these hypotheses, we can easily apply Lax-Milgram theorem which ensures that the penalty-projection step is well-posed and admits a unique solution v n+1 in the Hilbert space W(Ω) equipped with the usual norm of H 1 (Ω).

4. Numerical experiments. We consider the unit square as our computation domain Ω = ]0, 1[ 2 . The Dirichlet condition is prescribed on ∂Ω, except for the part included in the y-axis, where the open boundary conditions (2.4) are imposed. We present hereafter a convergence study for two manufactured test cases with an open boundary condition, homogeneous or not.

4.1. Homogeneous outflow boundary conditions ( g = 0 ). We choose a test case already used in the literature [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Jobelin | A finite element penalty-projection method for incompressible flows[END_REF]. It consists in unstationary Stokes problem, with a forcing term and initial and boundary conditions corresponding to the following expression for the velocity and pressure: u(x, y, t) = sin(x) sin(y + t) v(x, y, t) = cos(x) cos(y + t) p(x, y, t) = cos(x) sin(y + t) Convergence rate in space We are interested to study the space convergence rate for Stokes equations with open boundary conditions. In order to estimate the spatial error, we focus on the stationary solution of the above numerical experiment. We fix the time step at a small value: δt = 10 -2 , the penalty parameter ε at 10 -10 , the augmentation parameter r at 10 -4 and we run the algorithm for different values of the spatial mesh step h.

The convergence rate of the error on velocity and pressure is around O(h 2 ) (see Figs. [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]

.1 and 4.2).

As a conclusion on the spatial convergence rate, we can say that the results obtained here conform with those reported by Poux et al. [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. Besides, the (V P P r,ε ) method appears more accurate than the standard incremental scheme studied by Guermond et al. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. Particularly, our method improves the spatial convergence from O(h) to O(h 2 ) for the velocity and from O(h 1/2 ) to O(h 2 ) for pressure. Convergence rate in time We consider the unsteady homogeneous case. The convergence tests with respect to δt are reported in Figs 4.3 and 4.4. The convergence rate of the error on the velocity in the L 2 -norm is close to O(δt 2 ). The situation is somewhat similar for the pressure. The convergence results for this quantity in L 2norm seem to behave like O(δt 2 ). These rates seem similar to those obtained by Poux et al [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF] in the rotational form of their method. Also, the errors of velocity and pressure for the (V P P r,ε ) method seem to be smaller than that computed in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF] even if the mesh we have used is coarser. Moreover, the (V P P r,ε ) methods improves the poor convergence rates of the standard BDF2 pressure-correction [START_REF] Guermond | On the error estimates for the rotational pressure-correction projection methods[END_REF] scheme from O(δt) to O(δt 2 ) for the velocity and from O(δt 1/2 ) to O(δt 2 ) for the pressure. Indeed, in [START_REF] Guermond | On the error estimates for the rotational pressure-correction projection methods[END_REF], the second-order rotational pressure-correction yields O(δt 3/2 ) accuracy for the velocity in the L 2 -norm and O(δt) accuracy for the pressure.

It is useful to note that we checked numerically the value of v • n on the boundary Γ N in the projection step using (V P P r,ε ) method with OBC (3.6). We obtain that v • n is of order 10 -12 so, practically equal to zero.

Besides, Fig. 4.5 shows the L 2 -norm of the velocity divergence as a function of epsilon. We find that the constraint on the discrete velocity divergence is satisfied approximately as O(εδt) with 0 < ε ≤ 1 a small as desired. Moreover, in Fig. 4.6, we plot the L 2 -norm of the velocity divergence as a function of the time step δt. We obtain that the velocity divergence varies like O(ε δt) for ε quite small.

Finally, in order to check also the (V P P r,ε ) method using OBC1 (3.16) in the projection step, we simulate numerically the same test case already done. Figs. 4.7 and 4.8 display the errors of the computed velocity and gradient pressure in the L 2norm for ε = 10 -10 and r = 0. The numerical results show that a 2nd order accuracy in time is recovered for both velocity and gradient pressure. Convergence rate in time at r = 0 for the (V P P r,ε ) method with OBC (3.6) with the pressure gradient in (3.10) We treat now the temporal convergence of the velocity and pressure in the case when r = 0 and ε 1. For this purpose, we use BDF2 scheme for discretization in time, the open boundary condition OBC (3.6) in the projection step and we take the same above test. Fig. 4.9 displays the errors of the computed velocity in the L 2 -norm for r = 0 and ε = 10 -10 . The second-order convergence rate is not clearly observed contrary to the cases when 10 -4 ≤ r ≤ 1. Unfortunately, the slope of the velocity-error appears to be rather of first order. For the considered case also, Fig. 4.11 shows the distribution of the pressure error (using (13)). As expected, the behavior of the pressure error in L 2 -norm is strongly affected; we observe a sharp degradation of the pressure convergence (order 1/2 only) compared to the case when 10 -4 ≤ r ≤ 1. In fact, this deterioration of the convergence rate is due to the cumulation of round-off errors when ε is small enough.

Using the improvement of pressure explained in Section. 3.1.2, we repeat the same tests for r = 0. Fig. 4.10 presents the errors of the velocity in the L 2 -norm. The velocity exhibits now a second-order convergence rate in time as well as in the case when 10 -4 ≤ r ≤ 1. In Fig. 4.12, the pressure gradient field achieves successfully 1.8th-order rate of convergence. These good results confirm the interest of updating the pressure by its gradient and provide a second-order convergence in time for r ≥ 0 Homogeneous open boundary condition-OBC1-Velocity error L 2 -norm versus time step at mesh size: 1/h = 128, ε=10 -10 , r = 0. Homogeneous open boundary condition-OBC1-Gradient Pressure error L 2norm versus time step at mesh size: 1/h = 128, ε=10 -10 , r = 0. which is in agreement with [START_REF] Ph | A fast vector penalty-projection method for incompressible Darcy and multiphase Navier-Stokes problems[END_REF] nonhomogeneous case. For this aim, we take the same problem as in [START_REF] Liu | Open and traction boundary conditions for the incompressible Navier-Stokes equations[END_REF][START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]: ) is approximately observed for velocity and pressure. These results are in agreement with those reported in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. Indeed, it seems that the velocity divergence is almost of order (ε δt) for ε small enough. The results are not shown here for sake of shortness.

u(x, y, t) = cos 2 ( πx 2 ) sin(πy) cos(2πωt) v(x, y, t) = -cos 2 ( πy
We have also tested the convergence rate of velocity and gradient pressure using (V P P r,ε ) method with OBC1 (3.16). The results ensure also a second-order convergence for both unknowns (the results are not shown for brevity).

Conclusions.

In this work, we have extended the vector penalty-projection (V P P r,ε ) method to the case of incompressible viscous flows with open boundary conditions. The numerical experiments reported here show that the (V P P r,ε ) method yields a considerable gain in accuracy compared to incremental pressure-correction scheme as well as in the standard form as rotational form [START_REF] Guermond | On the error estimates for the rotational pressure-correction projection methods[END_REF] with small values of the augmentation parameter r; typically 0 ≤ r ≤ 1 and a penalty parameter ε small enough. We also show that using the BDF2 second-order scheme for time discretization, the (V P P r,ε ) method yields O(δt 2 ) accuracy for velocity and pressure as well as in the homogeneous and nonhomogeneous open boundary conditions. Besides, the loss of spatial convergence in case of open boundary conditions does not occur anymore with our method. In fact, we obtain O(h 2 ) convergence in the L 2 -norm of the velocity and pressure. The price to pay in the present method is that the constraint on the discrete divergence of velocity is not exactly equal to zero, but is satisfied approximately as O(ε δt) with a penalty parameter ε as small as desired up to machine precision. Finally, this method proves to be really efficient since it is fast, cheap, and provides very accurate results with an optimal spatial and temporal convergence despite the existence of outflow boundary conditions whereas other famous methods don't reach the same high convergence rate.

2 .

 2 Formulation of the continuous problem. Let Ω ⊂ R d (d = 2 or 3 in practice) be an open, bounded and connected domain with a Lipschitz continuous boundary Γ = ∂Ω. We suppose that Γ is partitioned in two subsets Γ D and Γ N , of outward unit normal vector n, such that Γ = Γ D ∪ Γ N , Γ D ∩ Γ N = Ø. The generic point of Ω is denoted x. For T > 0, we consider the time-dependent Stokes equations governing incompressible fluid flow on a finite time interval [0,T] which read:

Fig. 4 . 1 .Fig. 4 . 2 .

 4142 Fig. 4.1. Homogeneous open boundary condition -OBC -Velocity error L 2 -norm versus mesh step at T=2, δt=10 -2 , ε=10 -10 , r = 10 -4 .

Fig. 4 . 3 .Fig. 4 . 4 .Fig. 4 . 5 .Fig. 4 . 6 .

 43444546 Fig. 4.3. Homogeneous open boundary condition -OBC -Velocity error L 2 -norm versus time step at T=2, mesh size: 1/h = 128, ε=10 -10 .

  Fig. 4.7.Homogeneous open boundary condition-OBC1-Velocity error L 2 -norm versus time step at mesh size: 1/h = 128, ε=10 -10 , r = 0.

  Fig. 4.8.Homogeneous open boundary condition-OBC1-Gradient Pressure error L 2norm versus time step at mesh size: 1/h = 128, ε=10 -10 , r = 0.

Fig. 4 . 9 .Fig. 4 . 10 .

 49410 Fig. 4.9. Homogeneous open boundary condition -OBC -Velocity error L 2 -norm versus time step with standard pressure correction (3.8) at mesh size: 1/h = 128, r = 0, ε=10 -10 .

  .

4. 2 .Fig. 4 . 11 .Fig. 4 . 12 .

 2411412 Fig. 4.11. Homogeneous open boundary condition -OBC -Pressure error L 2 -norm versus time step with standard pressure correction (3.8) at mesh size: 1/h = 128, ε=10 -10 .

Fig. 4 . 13 .Fig. 4 . 14 .

 413414 Fig. 4.13. Nonhomogeneous open boundary condition -OBC -Velocity error L 2 -norm versus time step at T = 2, mesh size: 1/h = 128, ε=10 -10 .


R Martínez-Noriega

A Roumy

G Blanchard

EXEMPLAR-BASED IMAGE INPAINTING: FAST PRIORITY AND COHERENT NEAREST NEIGHBOR SEARCH

Keywords: exemplar, inpainting, coherence, priority

Greedy exemplar-based algorithms for inpainting face two main problems, decision of filling-in order and selection of good exemplars from which the missing region is synthesized. We propose an algorithm that tackle these problems with improvements in the preservation of linear edges, and reduction of error propagation compared to well-known algorithms from the literature. Our improvement in the filling-in order is based on a combination of priority terms, previously defined by Criminisi, that better encourages the early synthesis of linear structures. The second contribution helps reducing the error propagation thanks to a better detection of outliers from the candidate patches carried. This is obtained with a new metric that incorporates the whole information of the candidate patches. Moreover, our proposal has significant lower computational load than most of the algorithms used for comparison in this paper.

INTRODUCTION

The purpose of inpainting is to filling-in the missing part of an image. Diffusion-based inpainting preserves the structure (i.e. lines and object contours) by propagating the isophote (a line of equal luminance) in the unknown region (see the seminal work [START_REF] Bertalmio | Image inpainting[END_REF] and the overview in [START_REF] Bugeau | A comprehensive framework for image inpainting[END_REF]). By construction, this approach achieves excellent results when the missing region is small. In this paper, we focus on filling-in larger regions 1 . This problem occurs when blocks of data have been lost in transmission, or intentionally pruned for compression, or when an unwanted large object has to be removed. In this context, the missing part contains both image structure and texture. State-of-the-art methods are exemplar-based inpainting [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] and were inspired by texture synthesis algorithms [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Xu | Chaos mosaic: Fast and memory efficient texture synthesis[END_REF][START_REF] Ashikhmin | Synthesizing natural textures[END_REF].

Exemplar-based inpainting fills the holes in an image by searching similar information on the known region and simply copy it to the unknown region [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] as is done in texture synthesis [START_REF] Xu | Chaos mosaic: Fast and memory efficient texture synthesis[END_REF][START_REF] Ashikhmin | Synthesizing natural textures[END_REF]. The previous premise is based on the fact that natural images contain redundant or very similar information [START_REF] Zontak | Internal statistics of a single natural image[END_REF]. Greedy inpainting using exemplars, also known as patches, consists of two major steps: select the patch to be filled, and propagate the texture and structure. The former step selects patches with mostly linear structures by giving them higher priority. The latter is related to the selection of the 1 This problem is also known as image completion K-most similar patches from which the information is copied. These two steps are iterated until the holes in the image are fully restored.

The contribution presented in this paper is twofold. We first propose an improvement to the priority order proposed by [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF]. More precisely, [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] proposes a simple but efficient priority criterion based on the isophote strength and on the angle between the isophotes and the boundary (between the known and the unknown region). Different improvements have been proposed that can be classified with respect to their complexity. For instance, [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF]Chapter 6] proposes to add a term that counts the number of known pixels which belongs to an edge. Therefore an edge detection algorithm has to be added. [START_REF] Komodakis | Image completion using global optimization[END_REF] proposes to solve the global inpainting optimization problem with a message passing. To make the solution sequential, a scheduling is proposed that requires the computation of many distances at the initialization. More precisely, for each patch located at the boundary between the known and the unknown region, the distance between this patch and all known patches have to be computed. In [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF], only patches located in a neighborhood window are tested, in order to reduce the number of distances to compute. However, the priority term [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF] still increases significantly the complexity of the priority term computation. In this paper, we show that by better combining the two priority terms proposed by [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], we get significantly improvements in the structure propagation without any complexity increase.

Our second contribution concerns spatial coherence. Different attempts have been proposed to increase the spatial coherence of inpainted images. Some of them reduce the search of candidate patches i.e. the patches from which the information is copied (spatial neighborhood in [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF] or parents of the spatial neighborhood [START_REF] Ashikhmin | Synthesizing natural textures[END_REF]). Other methods, impose that candidates for overlapping patches are coherent either through an energy term [START_REF] Bugeau | A comprehensive framework for image inpainting[END_REF] or through a graph that represents the spatial connection between the underlying patches of an image [START_REF] Komodakis | Image completion using global optimization[END_REF]. Interestingly, all these methods, when they select a candidate patch, base their choice on partial data of the patch. More precisely, only the pixels that correspond to the known pixels of the patch to be filled in, are used. Instead, we propose a novel criterion that takes into account all the pixels in the candidate patch. The idea is to include the probability density function (pdf) of the patches to the metric used for measuring the similarity. With this improvement we are able to better detect outliers from the candidate patches and therefore reduce the error propagation which commonly affects algorithms based on exemplars.

The paper is organized as follows. In Sec. 2 the overall algorithm is introduced. Detail explanation of the two major steps: filling order 978-1-4673-1026-0/12/$31.00 c 2012 IEEE and propagation of textures-structures are described in Sec. 3 and Sec. 4 respectively. Sec. 5 presents experimental results, and Sec. 6 concludes the paper.

DEFINITIONS AND OVERALL ALGORITHM

Throughout the paper, the function I define an image [START_REF] Bertalmio | Image inpainting[END_REF]. Let pixel denotes the pixel position of the image, whereas pixel value refers to as the value taken by the image at this pixel. More precisely, a pixel is mapped to a scalar (d = 1) for a gray level image and a vector (d = 3) for a color image:

I : I -→ {0, ..., 255} d p -→ Ip (1)
The image domain I contains missing regions, i.e. holes on the image whose information is unknown. Define the unknown region as Ω, and the known region as Ω. The two regions define a partition of the whole image domain and are separated by a so called contour δΩ (also known as the "fill front"), see Fig. 1-(a). Our algorithm is patch based, as it was first introduced in [START_REF] Xu | Chaos mosaic: Fast and memory efficient texture synthesis[END_REF] for texture synthesis.

The patch Ψ i is a square subregion of I centered at pixel i. In practice, the size of Ψ i is decided to be slightly larger than the largest distinguishable texture element. In our case, and for the remaining of the paper, the size of the patches is 9 × 9 pixels. The algorithm aims to complete I by filling Ω with visually plausible information taken from Ω. This task is carried out by iterating two main steps: decision of the filling order and selection of suitable information to fill in the unknown part. More precisely, the algorithm first selects the patch Ψp, with p ∈ δΩ, whose filling priority is maximum among the patches of the fill front. Notice that Ψp contains a region already filled (known region) and a region where there is no information (unknown region). Then, the K-most similar patches Ψq i to Ψp are found, where i = 1, 2, . . . , K and Ψq ∈ Ω, Fig 1-(b). Finally, the Ψq i are linearly combined to predict the unknown region of Ψp. These steps are iterated until the unknown region Ω is filled. This two step approach (filling order and information propagation) has been originally proposed in [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF]. By construction, the algorithm is sequential. At each iteration n, the missing part of the current processed patch is inferred. Formally, let P be the matrix that selects the known pixels from the current patch Ψp. Let Ωn and Ωn denote the known and unknows regions at iteration n, respectively. The goal is to infer

f : {0, ..., 255} d×|Ψp∩Ωn| -→ {0, ..., 255} d×|Ψp| x = P Ψp -→ y (2)
In [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], the function f (2) is estimated with a 1-nearest neighbor (NN) approach. Further improvements have been proposed which are based on K-nearest neighbor (K-NN) search and Kernel smoothers [START_REF] Hastie | The Elements of Statistical Learning[END_REF]Chapter 6]. A Gaussian Kernel is used in [START_REF] Wong | A nonlocal-means approach to examplar based inpainting[END_REF], a Kernel with sparse constraint is used in [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF], whereas the weights results from a constrained least square optimization in [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF]Chapter 6]. We do not consider here the Kernel in [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF], since it is well suited to the sparsity filling order proposed in the same paper. Instead, we consider Kernels that use Criminisi's filling order and choose the one that give the best results that is [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF].

We now explain in detail the two main steps of the algorithm: selection of filling order, and propagation of texture-structure.

FILLING ORDER

Previous works, [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF][START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF][START_REF] Komodakis | Image completion using global optimization[END_REF], have shown that filling order is crucial for exemplar-based algorithms. The most popular idea hitherto is to propagate linear structures first and then textures, in such a way that the formation of prominent peninsulas in the fill front are avoided. In this Section, we introduce our definition of filling order that better encourage the propagation of structures before the synthesis of textures.

Given the patch Ψp for some p ∈ δΩ, Criminisi et al. [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] proposed a priority for the filling order P (p) based on the product of two terms

P (p) = C(p)D(p).
(3) The first term, C(p), is the confidence term and quantifies the belief we have in the known information contained in Ψp. The second term is the data term and it is a function of the strength of isophotes hitting the front δΩ at each iteration, and it is the responsible of encourage the propagation of linear structures. The terms are defined as follows:

C(p) = q∈Ψp∩Ω C(q) |Ψp| , D(p) = |∇I ⊥ p • np| α ( 4
)
where |Ψp| is the area of Ψp, α is a normalization factor (e.g., α = 255 for a typical gray-level image), and np is a unit vector orthogonal to the front δΩ in the point p. The function C(•) is initialized with C(p) = 1, ∀p ∈ Ω and once a patch has been filled, say Ψp, the confidence term of the filled pixels is set to the belief associated to the patch: C(q) = C(p), ∀q ∈ Ψp ∩ Ω. Therefore, the first term C(p) counts the ratio between the belief we have on the current patch and the belief associated with a perfect known patch (a known pixel has belief 1). The second term is defined on [0, 1], and it increases with the norm of the gradient |∇Ip| (strong isophote) and also if this isophote becomes more perpendicular to the contour. Therefore, when D(p) → 1, the patch is very likely to have a strong linear structure. Although (3) tends to propagate linear structures, it often fails to fully reconstruct the edges of natural images. Then, we propose to amplify the data term in a nonlinear fashion by modifying (3) in the following way:

P (p) = C(p) exp D(p) 2σ 2 . ( 5
)
The rationale for this exponential amplification is that we seek for a filling-order that promote the early propagation of strong structures. However we would like to keep the concentric filling order when regions are smooth. The exponential in [START_REF] Xu | Chaos mosaic: Fast and memory efficient texture synthesis[END_REF] reaches this compromise since it highly increases the influence of the data term when D(p) → 1 (strong linear structures), whereas no amplification is performed when D(p) → 0 (smooth region). Another rationale is that the confidence term C(p) (4) decreases exponentially with the distance between the current patch and the contour δΩ.

PROPAGATING TEXTURE AND STRUCTURE

Once the patch Ψp with maximum priority has been selected, the next step is to synthesize its missing region Ψp ∩Ω. This is achieved by searching K similar patches Ψq i , where Ψq i ∈ Ω and i = 1, 2, . . . , K, to the patch to be filled Ψp. Then, the unknown region of Ψp is filled with equivalent information extracted from the K similar patches Ψq i , e.g. using a linear combination or simply copy the information from the best match.

New metric for the selection of similar patches

The most popular way to chose K-similar patches is to select those patches that minimize the sum of square distance (SSD) [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], formally

Ψq i = arg min Ψq∈Ω dSSD(P Ψp, P Ψq), (6)
where dSSD is the SSD distance and the matrix P extracts the known region from the patch Ψp (pixels known or already filled). However as pointed out in [START_REF] Bugeau | A comprehensive framework for image inpainting[END_REF], this criterion can select patches that are visually different and the final result is not visually coherent. Moreover, dSSD(P Ψp, P Ψq) is a pixel-by-pixel metric constrained to the region defined by P and it does not use the whole information contained in Ψq.

We propose to replace dSSD in (6) by a more robust metric (7) that uses a combination of SSD and the Hellinger distance in the following way:

d(Ψp, Ψq) = dSSD(P Ψp, P Ψq) × dH (Ψp, Ψq), (7)
where dH stands for the Hellinger distance. The latter defines how similar the pdfs of two patches are. The motivation for such a metric is that with pdfs we can compare patches of different sizes and therefore include the whole information of the patch Ψq. Note that this is not the case with dSSD since the computation is performed on the known region defined by P Ψq. The values taken by dH are between [0, 1] and then, (7) can be seen as a weighted SSD. Hellinger distance is defined as:

dH (Ψp, Ψq) = 1 - N i=1 ρp(i)ρq(i), (8)
where ρp, ρq are the pdf of the patches Ψp, Ψq respectively. Those pdfs are computed with the normalized histograms of the patches. In total, N bins are used to compute the histograms (e.g. for color images in RGB format, N/3 bins are used for each channel). The normalization is given according to the number of elements of each patch and the number of color dimensions.

In [START_REF] Bugeau | A comprehensive framework for image inpainting[END_REF], a similar metric to (7) was proposed. However, the originality of our proposal is based on how the term ρp • ρq from ( 8) is computed, and used to achieve two goals: reduction of error propagation and increase the accuracy of the prediction. Moreover, we introduce a distinction between texture and structure patches that dictates the way in which the pdfs are computed. For textures patches we are interested to include a coherence constraint in order to reduce the error propagation, whereas with structures patches the pdfs are computed in such a way to increase the accuracy of the prediction.

Patch classification: texture vs. structure Let Ψp be the patch to be filled, and which has two well defined regions, the known part P Ψp and the unknown P Ψp. In order to decide if the patch is a texture or structure patch, we search for the three biggest gradient vectors inside P Ψp and average them. If the magnitude of the average gradient vector is higher than a threshold T the current patch Ψp is considered as structure patch, otherwise as texture one. We now define how to compute ρp • ρq depending on the nature of the patch, i.e. texture or structure.

Hellinger distance for Structure patches

We need to compare the pdfs of patches that are defined on different regions: known region for the current patch Ψp, whereas whole region for the candidate patches Ψq, taken from the dictionary. Therefore, in order to be comparable 3 , we need to infer the pdf of the current patch Ψp on the whole region. To do so, we need to review what a structure patch is.

A structure patch has gradient vectors with big magnitude that might indicate the presence of an edge, or in other words, two regions that contrast in colors. Our idea is to divide the current patch Ψp in two regions r1 and r2, such that pixels in the same region have similar colors, see Fig. 1-(c). For simplicity and also because the patches are rather small (9 × 9 is a current patch size), we assume that the boundary between the two regions is linear. In order to construct this linear edge across Ψp, we first compute the average vector of the three biggest gradient vectors4 of P Ψp. Then, the edge is the line defined by the average vector and xc, the center of the three pixels with biggest gradient vector.

The pdf of the current patch Ψp inferred on the whole domain can not be computed. Let R be the matrix that extracts the information contained in r1 and R for r2. Then, the pdf of Ψp is computed as

ρp = |RΨp| |Ψp| ρp 1 + |RΨp| |Ψp| ρp 2 , (9)
where:

ρp 1 = H(P Ψp ∩ RΨp) |P Ψp ∩ RΨp| , ρp 2 = H(P Ψp ∩ RΨp) |P Ψp ∩ RΨp| .
As for the candidate patch Ψq, since the patch is known on the whole region, the pdf is the usual one. Formally,

ρq = H(Ψq) |Ψq| . (10)
Finally, the combined metric [START_REF] Zontak | Internal statistics of a single natural image[END_REF] contains the Hellinger distance between the histograms ( 9) and [START_REF] Xu | Image inpainting by patch propagation using patch sparsity[END_REF]. Therefore, patches selected with ( 7) are more coherent with the current patch since it ensures that even on the unknown region, the candidate patch is similar to the current patch. The rationale to reinforce the coherence as shown above is to reduce the error propagation observed in many exemplarbased inpainting algorithms.

Hellinger distance for Texture patches Texture patches are rather smooth and there is no need to separate the patch into two regions. The pdf computation for a texture patch is as follows. Therefore, computations of the histograms of P Ψp (on known region) and Ψq (on the whole region) might sufficient. However, since the known regions of Ψp and Ψq are already compared in the SSD, a way to enforce coherence is to favor patches Ψq whose region P Ψq (in the unknown area of Ψp) is similar to P Ψp.

Formally, given the patch Ψp, the pdf of Ψp is computed as:

ρp = H(P Ψp) |P Ψp| , ( 11
)
where H is a function that computes the histogram, and | • | is the area of P Ψp. As for the candidate patch Ψq, where Ψq ∈ Ω, we define

ρq = H(P Ψq) |P Ψq| . (12)
Finally, the combined metric [START_REF] Zontak | Internal statistics of a single natural image[END_REF] can be seen as a weighted SSD distance that is biased toward patches Ψq which are coherent. Here again, the goal is to avoid the error propagation phenomenon observed in many exemplar-based inpainting algorithms.

Synthesis

Once the K nearest patches Ψq i are found either using the texture or structure definition, we use a linear combination of the patches to finally obtain the prediction of the missing pixels:

P Ψp = K i=1 ωiP Ψq i .
As in [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF], the coefficients ω are obtained by solving

min ω K i=1 ωiP Ψq i -P Ψp 2 s.t. K i=1 ωi = 1. (13)
The optimization problem given in ( 13) can be solved in the same fashion as applied in local linear embedding (LLE) for data dimensionality reduction. In our case, the Gram matrix is G = (P Ψp1 -X) T (P Ψp1 -X), where X is a matrix with columns {P Ψq 1 , P Ψq 2 , • • • , P Ψq K } and 1 is a column vector of ones. Then it has a close form solution

ω = G -1 1 1 T G -1 1 .

EXPERIMENTS AND COMPARISONS

In this Section, we present simulation results of our algorithm on natural images, and comparisons with other popular algorithms from the literature. We also show empirical evidence of the advantages conveyed by the ideas proposed in this work i.e. improved filling order and spatial coherence. Fig. 2 is a summary of inpainted images using different methods considering that the light green region shown in (a) is the missing information. We compare our proposal with three other patch based inpainting algorithms. (b) is the seminal proposal of [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] in which the filling order is decided with (3) and the synthesis is achieved by copying the information from the most similar patch only. Turkan [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF], (c), uses edge detection to augment the priority defined in [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] and the synthesis is produced by a linear combination of patches. Recently, Adobe Photoshop has introduced a tool called PatchMatch [START_REF] Barnes | Patchmatch: A randomized correspondence algorithm for structural image editing[END_REF] to inpaint missing regions on images. The inpainted results of this tool is reported in (d). Finally, the last three columns are our results. (e) is the result of the algorithm as exactly described in this paper. In (f), we explore a different metric, structural similarity index (SSIM) [START_REF] Zhou Wang | Image quality assessment: from error visibility to structural similarity[END_REF] that is more consistent with human eye perception than SSD. Therefore, SSIM replaces SSD for the selection of the most similar patches. As for the weight computation in the synthesis step, SSIM cannot be used directly by LLE, since LLE directly refers to euclidean distances. Therefore in (f), the K-most similar patches were selected with SSIM and then, their SSD were used to compute the weights in the synthesis. In order to involve only the SSIM in the weight computation, we decide to change LLE by nonlocal means (NLM) [START_REF] Wong | A nonlocal-means approach to examplar based inpainting[END_REF]. This result is shown in column (g). The first row corresponds to a block erasure scenario 5 . In this context, the ground truth is known, which allows to compute a PSNR value after reconstruction of the missing part. In all the other rows, an unwanted large object has to be removed and therefore no PSNR can be computed.

The strongest characteristic of our algorithm is to stop the error propagation that highly affect patch-based algorithms. This is more evident in the block erasure scenario, shown in the first row of Fig. 2, in which our algorithm outperforms other proposals by at least 12 dB.

In general our proposal provides a better completeness of linear edges and a more natural look of the missing region. For example in the Bungee jumping image where the man is complety erased, [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] and PatchMatch [START_REF] Barnes | Patchmatch: A randomized correspondence algorithm for structural image editing[END_REF] fail to complete the roof of the shelter. Instead, our result not only completes the edge but also provides more visually pleasant images than [START_REF] Turkan | Novel texture synthesis methods and their application to image prediction and image inpainting[END_REF]. Visual coherence is only achieved with our proposal for the Soldier image. In this case, the horizontal edge between the grass and the wall on the left part is not linearly propagated by our algorithm because there is a stronger edge on the right part (black plastic bag and grass) that is filled in first. The image Bench was the toughest for all the algorithms, but our proposal is the only one to complete more naturally the edges of the bench. In the Jaguar image, our result is the only one which fully completes the tree branch and eliminates the propagation of errors in the leftupper region of the branch.

SSIM metric is promising since it estimates perceived errors, a characteristic that matters in inpainting. However, SSIM is applied only on luminance value and therefore faces problems of color mismatch. This can be notice in (g) of Fig. 2-Soldier, where SSIM chooses patches with similar structure but different color for the wall. The same problem is shown in Soldier-(f) but reduced, because, in this case, SSIM is combined with SSD.

We introduced two novel ideas to improve the performance of exemplar-based inpainting. First, a modified priority for the filling order has been introduced in order to better reproduce the edges in the missing part. Second, a novel way to select the candidate patches has been proposed in order to reduce the propagation of errors in the inpainted image. We now focus on the former one. Fig. 3 shows the filling order used by [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] and our algorithm. We represent the filling order with a gray scale of colors. The darker patches correspond to places that are filled-in first. Unlike [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], our proposal better encourages the early propagation of linear structures without disturbing the concentric filling order in regions that are smooth. This trade-off between synthesis of structure-textures is important to obtain good performance. Since our proposal of filling order is based on increasing the influence of the data term D(p), one might think that the filling order which first propagates structures only, would produce better results. However allowing the structures or edges to fully drive the inpainting process is not efficient, because it causes the formation of peninsulas in the fill front, as observed in [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF].

Our second contribution aims to reduce the error propagation generated by the presence of outliers in the candidate patches when dSSD is used to select those patches. We propose a more robust metric [START_REF] Zontak | Internal statistics of a single natural image[END_REF] that account the whole information from the candidate patches. In Fig. 4 we show how our metric can better distinguish the outliers patches over the popular dSSD. We notice that patches chosen with dSSD contain artifacts that are not well suited for the neighborhood whereas our metric selects more visually pleasant patches.

CONCLUSION

We proposed an inpainting algorithm that produces better completeness of linear edges and reduces the error propagation problem. The empirical results on natural images show better performance than well-known greedy algorithms, and better or similar performance than highly complex algorithms that need to inpaint the full missing part several times until convergence. Instead, our algorithm only needs to inpaint the missing region once and therefore the complex-ity load is lower. Moreover, the improvement in the filling-in order does not increment the complexity. In this proposal, we used crossvalidation to compute the threshold T which distinguishes texture from structure patches, and the variance σ 2 that controls the amplification of the data term. An open problem is an easy way to tune the parameters or even avoid them. Further work includes improving the classification of the patches and the inpainting adapted to each class, in order to obtain a less parametric technique.

Fig. 1 .

 1 Exemplar based inpainting. (a) Original image I with the unknown region Ω, its contour δΩ, and the known region Ω. (b) K patches Ψ q i are found to synthesize the unknown area of the patch Ψ p , with p ∈ δΩ. (c) Example of the imaginary edge that helps the pdf computation for patches classified as structures.

  Therefore, in order to compensate for the high values of C(p) near the contour, an exponential amplification of data term D(p) (4) has to performed. Finally, in order to be an amplification, the function f of the data term has to satisfy f (D) ≥ D. The function expD(p)

2 .

 2 Missing part (b) Criminisi [3] (c) Turkan [8] (d) PatchMatch [13] (e) Ours SSD (f) Ours SSIM (g) Ours SSIM-NLM Fig. Comparison of inpainted images obtained with our proposals, (e)-(g), and popular algorithms from literature. In columns, (a) Light green color represents the missing region. (b) Algorithm proposed in [3]. (c) Inpainted images obtained with [8]. (d) Resulting images after inpainting with PatchMatch tool [13] from Adobe Photoshop.

Fig. 3 .Fig. 4 .

 34 Filling order for the "Bungee jumping". Darker patches are regions that are inpainted first. (a) Location of the patch (b) Patch to be filled. Our metric produces better selection of candidate patches than d SSD . (a) The red square show the location to be filled. (b) Light green color is the region to be filled. (c) From left-to-right the five most similar patches to (b) chosen with d SSD . (d) Five most similar patches chosen with our metric.

Note that other exponential functions exist, as D(p) n . But these functions do not satisfy the above mentioned amplification property.

Note that the usual way to make the patches comparable, is to restrict the computation to the known region. But this is exactly what we want to avoid.

For the sake of brevity, we refer to the gradient vector with biggest magnitude value as the biggest gradient vector.

This problem occurs when blocks of data have been lost in transmission, or intentionally pruned for compression.

ACKNOWLEDGMENTS

This work was supported by the ANR-09-VERS-019-02 grant (AR-SSO project). The authors would also like to thank the reviewers for their comments that help improve the manuscript.


Ariane Ducellier

Hiroshi Kawase

Shinichi Matsushima

VELOCITY STRUCTURE INVERSION FROM H/V SPECTRAL RATIOS OF EARTHQUAKE DATA: APPLICATION TO THE TOHOKU REGION, JAPAN

Keywords: H/V spectral ratio, site effects, inversion

come

Velocity structure inversion from H/V spectral ratios of earthquake data: Application to the Tohoku region, Japan

INTRODUCTION

Multiple scattering waves sample the medium along their multiple paths. In the acoustic case, a statistical approach for wave propagation led to the concept of diffuse fields. In the seismic case, a uniformly distributed set of random forces provides an equipartitioned illumination of the medium (e.g. Sánchez-Sesmaet al., 2008). The seismic field engendered by such equipartitioned illumination behaves as a diffusion-like regime. When the seismic motion is dominated by multiple scattering of waves, the energy densities have therefore a diffusionlike behavior (Margerinet al., 2009). Under these conditions, we can retrieve the Green's function from averaging cross-correlations of the recorded motion of such diffuse field [START_REF] Sánchez-Sesma | Retrieval of the Green's function from cross correlation: the canonical elastic problem[END_REF].

In this study, we assume that the seismic field engendered by the summation of a sufficiently large number of seismic events covering a large interval of azimuths and incidence angles can be considered as a diffuse-wave field. A new formulation for the average horizontal-tovertical (H/V) spectral ratio of earthquake motions was then constructed by Kawaseet al. (2011). This leads to a new method to retrieve the underneath velocity structure.

We consider strong ground motion in elastic layered media. Kawaseet al. (2011) explored the theoretical consequences of assuming both a "sufficiently" flat layered site and "sufficiently" deep earthquake sources, such that surface waves are negligible. Under these conditions, the illumination is produced by incident plane waves and a one dimensional description of the wave propagation can be applied. The ground motion is supposed to be spatially homogeneous in a statistical sense, that is to say the average of normalized ground motion spectral densities will depend only on depth. By summing a sufficiently large number of earthquake motions, the averaged autocorrelations can thus be related through the diffuse field formalism to the imaginary part of the one dimensional Green's function, when the source and the receiver are both at the same point. [START_REF] Claerbout | Synthesis of a layered medium from its acoustic transmission response[END_REF] established a relationship between reflection response and autocorrelation of surface motion for one dimensional layered media. Kawaseet al. (2011) interpreted these pioneering results under the light of diffuse-field concepts and proposed an extension of Claerbout's results. They assert that the imaginary part of the one dimensional Green's function at the surface for a surface source is equal to the square of the absolute value of the transfer function for an incoming unit displacement divided by four times the circular frequency and the half-space impedance.

By connecting this last result with the illumination due to earthquakes, Kawaseet al. (2011) established a relationship between the average horizontal-to-vertical (H/V) spectral ratio of earthquake motions and the ratio of transfer functions for P-and S-waves. In this study, we focus on three sites from the KiK-net network in the Tohoku area, Japan. We can compute the theoretical transfer function of a soil column from the soil layers' mechanical properties using the Thomson-Haskell propagator matrix method [START_REF] Thomson | Transmission of elastic waves through a stratified solid[END_REF][START_REF] Haskell | The dispersion of surface waves in multilayered media[END_REF]. We then carry out inversions of the velocity structure for the three sites considered and verify that there is a good match between the earthquake H/V spectral ratio and the transfer function corresponding to the new velocity structure obtained from the inversion.

INVERSION METHOD

From the relationship between the average autocorrelation of the recorded motions and the one dimensional Green's function, Kawaseet al. (2011) deduced that the average horizontalto-vertical (H/V) spectral ratio can be written as follows: (2.2)

The relationship established by Kawaseet al. (2011) between the imaginary part of the one dimensional Green's function and the transfer function for an incoming unit displacement can be written as:

𝐼𝑚 𝐺 𝑖𝑖 1𝐷 𝑥, 𝑥, 𝜔 = 4𝜔𝜌 𝐻 𝑐 𝐻 -1 𝑇𝐹 𝑖 𝑥, 𝜔 2 (2.3)
where 𝑇𝐹 𝑖 𝑥, 𝜔 is the i th component of the transfer function between the top of the half-space and the position x, ρ H is the density of the half-space and c H is eitherthe P-or S-wave velocity depending on whether we consider a vertical or an horizontal displacement.

By specializing Equation (2.3) for both horizontal and vertical motion, we can express the H/V spectral ratio at the soil surface (x = 0) and on top of the half-space (x = h) with the transfer function for the horizontal motion due to the S-wave TF 1 (ω), and for the vertical motion due to the P-wave TF 3 (ω). We get:

𝐻 0,𝜔 𝑉 0,𝜔 = 2𝛼 𝐻 𝛽 𝐻 𝑇𝐹 1 0,𝜔 𝑇𝐹 3 0,𝜔 and 𝐻 ℎ,𝜔 𝑉 ℎ,𝜔 = 2𝛼 𝐻 𝛽 𝐻 𝑇𝐹 1 ℎ,𝜔 𝑇𝐹 3 ℎ,𝜔 (2.4)
whereα H is the P-wave velocity and β H is the S-wave velocityof the half-space (bedrock).

The Thomson-Haskell propagator matrix method is a frequency domain method due to [START_REF] Thomson | Transmission of elastic waves through a stratified solid[END_REF] and corrected by [START_REF] Haskell | The dispersion of surface waves in multilayered media[END_REF]. A generalization of the method was introduced in seismology by [START_REF] Gilbert | Propagator matrices in elastic wave and vibration problems[END_REF]. The Thomson-Haskell method uses the equations of motion to establish the relationship between the motion-stress vector at depth z 1 and the one at depth z 2 . As long as we know the mechanical properties of the soil column (e.g. S-and P-wave velocities, density and thickness of each layer in the elastic case), we can compute easily both the horizontal and vertical transfer functions from the Thomson-Haskell method. By minimizing the difference between the ratio of transfer functions and the average H/V spectral ratio at the soil surface and on top of the half-space, it becomes possible to inverse the velocity structure below the station from the bedrock to the surface.

Genetic Algorithms use techniques inspired by evolutionary biology to find exact or approximate solutions to optimization problems. In this study, we use the Genetic Algorithm code developed in De Martin et al. (2010) and adapt it for the case where the H/V spectral ratio is the objective function.

INVERSION OF VELOCITY STRUCTURE FOR THREE STATIONS

For the three stations FKSH06, FKSH09 and IWTH19 of the KiK-net network, we invert the S-wave velocities, the P-wave velocities and the thicknesses of the soil layers. As the H/V spectral ratio depends also on the ratio of the P-and S-wave velocities of the half-space, we also invert the S-wave velocity of the half-space while the P-wave velocity is kept constant, for simplification. However, it should be noted that any couple of P-and S-wave velocities of the half-space for which the ratio between both velocities is the same, could be a valid solution.

The inversion code we use in this study is based on a Genetic Algorithm; such that the solution of the optimization problem is searched into a finite number of different possible soil columns. Each parameter of the soil column can take a finite number of discrete values.In most cases, the invertedparameter can take 32 different values, ranging from 0.25 to 1.80 timesthe initial value. However, forthe S-wave velocity of the last layer and the half-space, the rangeof variations is smaller, ranging from 0.50 to 1.25 the initial value, in order to avoid a negative Poisson ratio. These last parameters can take only 16 different discrete values. The initial values of the soil parameters for stations FKSH06, FKSH09 and IWTH19 are shown in Tables 1, 2 and 3 respectively (column 3). They correspond to the velocity log posted onthe KiK-net website. Tables 1, 2 and 3 show also the intervals of variations where we look for a solution for each parameter (column 4).

The thicknesses of the shallower layers are inverted but the total depth of each soil column remains constant. The thicknesses of layer 6 of FKSH09, layer 4 of FKSH09, and layer 4 of IWTH19 may then vary and are computed such that the total depth is constant. To evaluate if the ratio of transfer functions corresponding to a given set of soil column parameters fits well the objective H/V spectral ratio, we use the following residual:

𝑅𝑒𝑠 = 1 2 𝐻𝑉𝑇 𝜔 -𝐻𝑉𝑇 0 𝜔 2 𝑑𝜔 𝐻𝑉𝑇 0 𝜔 2 𝑑𝜔 + 𝐻𝑉𝐵 𝜔 -𝐻𝑉𝐵 0 𝜔 2 𝑑𝜔 𝐻𝑉𝐵 0 𝜔 2 𝑑𝜔 (3.1)
where HVT 0 (ω) and HVB 0 (ω) are the objective H/V spectral ratios obtained from the observation data at frequency ω on the top (soil surface) and the bottom (top of the halfspace) of the soil column, and HVT(ω) and HVB(ω) are the theoretical ratios of horizontal and vertical transfer functions on the top and the bottom of the soil column. The integration is done over the interval 0.1 to 25 Hz for all the inversions carried out in this study.

As the Genetic Algorithm generates random numbers, two successive inversions with the same input parameters will not necessarily give the same result. Therefore, for each station, we carry out ten independent inversions with the same input files to check whether we always obtain the same best soil column at the end of the inversion. The results of the inversion for stations FKSH06, FKSH09 and IWTH19 are shown respectively on Figure 1, 2 and 3. On the left panel, the abscissa is frequency and the ordinate is the H/V spectral ratio on top of the soil column. The dashed black line corresponds to the objective H/V spectral ratio computed from the observation data, while the plain black line is the best ratio of transfer functions obtained from the ten independent inversions. On the middle panel, the ordinate is the H/V spectral ratio on bottom of the soil column. On the right panel, the abscissa is velocity and the ordinate is the depth from the surface. The two plain black lines are the S-wave velocity and the Pwave velocity of the soil column corresponding to the best ratio of transfer functions on the leftand middlepanels. The two grey lines are the initial S-and P-wave velocities. The best values of the S-and P-wave velocities and the thicknesses of soil layers obtained from the inversion are also given in Tables 1, 2 and 3 (column 6). We have now obtained the soil column parameters corresponding to the transfer function fitting best the two H/V spectral ratios. However, all the parameters are not well constrained by the inversion process. It is possible to find soil columns whoseparameters take different values, yetcorresponding transfer function still stays very close to the best solution from the inversion. We thus take the best soil column obtained from the inversion and we make each parameter vary, one at a time, from 0.25 to 1.80 the best value, as we did during the inversion process. We then compute the corresponding ratio of transfer functions and compare it to the ratio of the best soil column.

We need to define a criterion to establish whether the ratio of a soil column is close enough to the best ratio. We use a bootstrap method to estimate the variance of the H/V spectral ratio. We select randomly NE accelerograms from the data set and compute the corresponding H/V spectral ratio. This step is repeated NB times. We then compute the mean and the variance of the NB H/V ratios thus obtained, for each frequency. In this study, we use NE = NT, NT being the total number of seismic events available in the data set. For all computations we take NB = 500. We then compute the residual between the mean H/V spectral ratio plus the square of the variance and the mean H/V spectral ratio, and take it as a target residual to compare soil columns. If the residual between the ratio of transfer functions of a given soil column and the ratio of the best soil column is lower than this target residual, we consider that both soil columns are as good a solution of the inversion problem. We check if each soil column is a good solution for both H/V spectral ratio on top and H/V spectral ratio on bottom of the soil column. The column 5 of Tables 1, 2 and 3 gives the range of values where each parameter can vary while keeping a good fit between the two objective H/V spectral ratios and the two ratios of the horizontal and vertical transfer functions. The grey lines on left and middle panels of Figures 1, 2 and 3 represent the corresponding H/V ratio on top and on bottom of the soil column. However, it should be noted that these ranges of values were established while the soil parameters varied one at a time. It does not necessarily imply that any set of parameters picked up in this table will give a ratio of transfer functions fitting well the H/V spectral ratio.

CONCLUSION

A series of inversions of underground 1D structure from the bedrock to the surface from three observations sites in the Tohoku region, Japan, has been carried out with a Genetic Algorithm code. A good agreement could be found between the H/V spectral ratio at the soil surface obtained with the inversion and the objective ratios. However, the agreement is not as good for the H/V spectral ratio on top of the half-space, particularly for the higher frequencies. It should be noted that the soil column parameters can take only a few different discrete values in our inversion scheme. Inverting the velocity structure with a larger number of possible values in a narrower interval of variations could lead to better results.

Figure 3 . 1 .Figure 3 . 2 .

 3132 Figure 3.1. Result of inversion for station FKSH06

Table 3

 3

			.1. Inversion parameters for station FKSH06
	Parameter	Layer	Initial	Interval of variations	Rangeof values	Best value
			value		obtained from the	obtained from
					inversion	the inversion
	V S (m/s)	1	200	0.25 -1.80 * Initial value	47.5 -52.5	50
		2	300	0.25 -1.80 * Initial value	283.5 -729	405
		3	650	0.25 -1.80 * Initial value	250.25 -643.5	357.5
		4	1200	0.25 -1.80 * Initial value	231 -630	420
		5	1400	0.25 -1.80 * Initial value	1330 -1470	1400

Table 3

 3

			.2. Inversion parameters for station FKSH09
	Parameter	Layer	Initial	Interval of variations	Range of values	Best value
			value		obtained from the	obtained from
					inversion	the inversion
	V S (m/s)	1	140	0.25 -1.80 * Initial value	166.25 -183.75	175
		2	300	0.25 -1.80 * Initial value	324 -396	360
		3	1930	0.25 -1.80 * Initial value	405.3 -926.4	579
		4	2540	0.25 -1.80 * Initial value	1447.8 -1600.2	1524
		5	1960	0.50 -1.25 * Initial value	1724.8 -2587.2	2156
		Half-space	1960	0.50 -1.25 * Initial value	1666 -2352	1960
	V P (m/s)	1	450	0.25 -1.80 * Initial value	567 -1458	810
		2	1500	0.25 -1.80 * Initial value	1721.25 -2632.50	2025
		3	3480	0.25 -1.80 * Initial value	4541.4 -5550.6	5046
		4	4950	0.25 -1.80 * Initial value	3056.625 -3378.375	3217.5
		5	4950	0.25 -1.80 * Initial value	7573.5 -11583	8910
	Thickness (m) 1	2	0.25 -1.80 * Initial value	2.945 -3.255	3.1
		2	8	0.25 -1.80 * Initial value	9 -11.5	10
		3	34	0.25 -1.80 * Initial value	47.6 -71.4	59.5
		4	126	Not inverted
		5	30	Not inverted

Table 3

 3

			.3. Inversion parameters for station IWTH19
	Parameter	Layer	Initial	Interval of variations	Range of values	Best value
			value		obtained from the	obtained from
					inversion	the inversion
	V S (m)	1	170	0.25 -1.80 * Initial value	110.5 -110.5	110.5
		2	550	0.25 -1.80 * Initial value	235.125 -940.5	522.5
		3	700	0.25 -1.80 * Initial value	308 -462	385

ACKNOWLEDGEMENT

Part of this work was carried out when one of the authors (A.D.) was staying at the Disaster Prevention Research Institute of Kyoto University, at Uji, Kyoto, Japan, supported by the Carnot funding program of the French National Research Agency (ANR) under grants Carnot 2009-ac17 and Carnot 2010-ac24.