Léo Agélas
email: leo.agelas@ifpen.fr

Global regularity of Quantum Navier-Stokes equations in R 2

come

Introduction

Quantum hydrodynamic models become important and necessary to model and simulate electron transport, affected by extremely high electric fields, in ultra-small sub-micron semiconductor devices, such as resonant tunnelling diodes, where quantum effects (like particle tunnelling through potential barriers and built-up in quantum wells [START_REF] Gardner | The quantum hydrodynamic model for semiconductor devices[END_REF][START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF] take place and dominate the process. They arise in semiclassical mechanics in the study of semiconductor devices, in which case being derived from the Wigner-Boltzmann equation [START_REF] Ferry | Form of the quantum potential for use in hydrodynamic equations for semiconductor device modeling[END_REF][START_REF] Markowich | Semiconductor Equations[END_REF][START_REF] Gardner | The quantum hydrodynamic model for semiconductor devices[END_REF][START_REF] Gasser | Quantum hydrodynamics, Wigner transforms and the classical limit[END_REF][START_REF] Jüngel | Quasi-hydrodynamic semiconductor equations[END_REF][START_REF] Gardner | Dispersive/Hyperbolic Hydrodynamic Models for Quantum Transport (in Semiconductor Devices)[END_REF]. Quantum hydrodynamics has engendered substantial activity in the field of theoretical chemical dynamics, where one may refer to Wyatt et. al.([36]) for a comprehensive introductory overview of the numerous recent results emerging from this blossoming field. In quantum chemistry, they arise as solutions to chemical kinetic systems, in which case they are derived from the Schrödinger equation by way of Madelung equations [START_REF] Madelung | Quantentheorie in hydrodynamischer Form[END_REF]. The basic idea emerging from quantum chemistry is to employ the time-dependent Schrödinger equation to solve dynamical properties (probability densities, "particle" velocities, etc.) of chemical systems. In the same spirit in which the de Broglie-Bohm interpretation (see [START_REF] Bohm | Quantum Theory[END_REF][START_REF] Bohm | A suggested interpretation of the quantum theory in terms of "hidden variables", i and ii[END_REF][START_REF] Brunton | Kinetic applications of electron paramagnetic resonance spectroscopy. 26. Quantum-mechnical tunneling in the isomerization of sterically hindered aryl radicals[END_REF]) of quantum mechanics may be used to recover "trajectories" of individual fluid elements along the characteristics of motion of the solution (see [START_REF] Wyatt | Quantum Dynamics with Trajectories[END_REF] and [START_REF] Jüngel | Quasi-hydrodynamic semiconductor equations[END_REF] for a comprehensive overview). They are also used to describe superfluids [START_REF] Loffredo | On the creation of quantum vortex lines in rotating He II[END_REF], weakly interacting Bose gases [START_REF] Grant | Pressure and stress tensor expressions in the fluid mechanical formulation of the Bose condensate equations[END_REF]. Some other topics of interest in quantum hydrodynamics are quantum turbulence, quantized vortices, second and third sound, and quantum solvents. Later, so-called quantum hydrodynamic equations have been derived by Ferry and Zhou [START_REF] Ferry | Form of the quantum potential for use in hydrodynamic equations for semiconductor device modeling[END_REF] from the Bloch equation for the density matrix and by Gardner [START_REF] Gardner | The quantum hydrodynamic model for semiconductor devices[END_REF] from the Wigner equation by a moment method. More recently, dissipative quantum fluid models have been proposed. For instance, the moment method applied to the Wigner-Fokker-Planck equation leads to viscous quantum Euler models [START_REF] Gualdani | Analysis of the viscous quantum hydrodynamic equations for semiconductors[END_REF], and a Chapman-Enskog expansion in the Wigner equation leads under certain assumptions to quantum Navier-Stokes equations [START_REF] Brull | Derivation of viscous correction terms for the isothermal quantum Euler model[END_REF]. In this paper, we consider the Cauchy problem for Quantum Navier Stokes equations as follows:

     ∂n ∂t + ∇ • (nu) = 0, ∂nu ∂t + ∇ • (nu ⊗ u) -2ε 2 n∇ ∆ √ n √ n -2ε∇ • (nDu) + ∇p(n) = 0 (1)
with initial conditions, n(x, 0) = n 0 , u(x, 0) = u 0 ,

on O = R 2 /L 0 Z 2 , L 0 > 0 a given real and with periodic boundary conditions and we require in addition that n 0 and u 0 are periodic functions of period L 0 where Du = ∇u + ∇u T

2

, p(n) = gn 2 with g > 0 and ε > 0 is the scaled Planck constant. From [START_REF] Gualdani | Analysis of the viscous quantum hydrodynamic equations for semiconductors[END_REF][START_REF] Jüngel | Derivation of New Quantum Hydrodynamic Equations Using Entropy Minimization[END_REF], we can assume that the scaled Planck constant ε is of order 10 -2 . For the initial data n 0 , we suppose that it is a small perturbation of some positive constant n 0 > 0 . There are only few mathematical results for these viscous quantum hydrodynamic model due to difficulties coming from the third-order derivatives. The existence of classical solutions to the one-dimensional stationary model with ε = 0 and with physical boundary conditions was shown in [START_REF] Jüngel | Physical and numerical viscosity for quantum hydrodynamics[END_REF]. The transient equations are considered in [START_REF] Chen | The viscous model of quantum hydrodynamics in several dimensions[END_REF][START_REF] Chen | Viscous quantum hydrodynamics and parameter-elliptic systems[END_REF][START_REF] Dreher | The transient equations of viscous quantum hydrodynamics[END_REF], and the local-in-time existence and exponential stability of solutions were proved. Global-in-time solutions in one space dimension are obtained if the initial energy is assumed to be sufficiently small. We also mention that in the inviscid case (ε = 0) there is a recent proof of nonglobal-in-time existence for a quantum hydrodynamic equation in bounded domains with prescribed data corresponding to high boundary and initial energy [START_REF] Gamba | On the blowing up of solutions to the quantum hydrodynamic equations in a bounded domain[END_REF]. Later, Existence of global-in-time weak solutions in one space dimension without smallness conditions is proved in [START_REF] Gamba | Global existence of solutions to one-dimensional viscous quantum hydrodynamic equations[END_REF]. Concerning the multidimensional case, local-in-time existence theorems have been obtained in [START_REF] Chen | The viscous model of quantum hydrodynamics in several dimensions[END_REF][START_REF] Dreher | The transient equations of viscous quantum hydrodynamics[END_REF]. Global-in-time existence of weak solutions in a three-dimensional torus for large data is achieved in [START_REF] Jüngel | Global weak solutions to compressible Navier-Stokes equations for quantum fluids[END_REF][START_REF] Dong | A note on barotropic compressible quantum navier-Stokes equations[END_REF][START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF]. This proof of existence given in [START_REF] Jüngel | Global weak solutions to compressible Navier-Stokes equations for quantum fluids[END_REF] relies on the reformulation of the Quantum Navier-Stokes model as a viscous quantum Euler system and vice versa by introducing a new velocity variable ,w, involving gradients of the particle density, w = u + ε∇ log n. Following the results obtained in [START_REF] Jüngel | Effective velocity in Navier-Stokes equations with third-order derivatives[END_REF], it is shown provided that the particle density n = 0, that the particle density n and the new velocity w solve a viscous Euler system and this new formulation is equivalent to the Quantum Navier-Stokes equations ( 1)- [START_REF] Ball | Remarks on blow-up and nonexistence theorems for nonlinear evolution equations[END_REF]. In this paper, by using this new viscous Euler system with variables (n, m), m ≡ nw = nu + ε∇n (see [START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF]), we prove existence and uniqueness of global strong solutions in C([0, +∞[; H s (O)), s > 1 of the Quantum Navier-Stokes equations for large initial data for both n 0 and u 0 . The proof comes from results obtained after introducing the variables ρ (new particle density) and v (new velocity) defined from n and m as follows ρ(x, t) = . In Proposition 4.1, under the assumption that ρ L ∞ < 1, we show that the blow-up of smooth solutions (ρ, v) to the viscous Euler equations [START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF] is controlled by the time integral of the maximum magnitude of the velocity v to the power four. More precisely, in Lemma 4.1, we show the following energy estimates in the homogeneous Sobolev space Ḣs (Ω), with

Ω def = √ gn 0 O ( def = R 2 /( √ gn 0 L 0 )Z 2 ) and s ≥ 0, (ρ(t), v(t)) Ḣs ≤ (ρ 0 , v 0 ) Ḣs e C R t 0 ( ρ(τ ) 2 L ∞ + v(τ ) 2 L ∞ + v(τ ) 4 L ∞ )dτ , (3)
where (ρ 0 , v 0 ) are the initial data corresponding to the variables (ρ, v). Furthermore, from a BD entropy estimates and under the assumption that ρ

L ∞ ≤ 3 4 , in Lemma 4.2, we show, (ρ(t), v(t)) 2 L 2 + ε t 0 ∇ v(τ ) 1 + ρ(τ ) 2 L 2 + ∇ρ(τ ) 2 L 2 dτ (ρ 0 , v 0 ) 2 L 2 . ( 4
)
Since these Inequalities are valid only for ρ L ∞ ≤ 3 4 , the challenge was then to find bounds for both ρ L ∞ and v L ∞ allowing to ensure that ρ L ∞ ≤ 3 4 and v L ∞ ≤ 1 as soon as n 0 is sufficiently large relative to (n 0 -n 0 , u 0 ) in some sense. It was thus important to obtain inequalities invariant in respect of any scaling. For this, thanks to Gagliardo-Nirenberg inequalities and under the assumption that ρ L ∞ ≤ 3 4 and v L ∞ ≤ 1, from (3) and (4), we show that for all s > 1,

(ρ(t), v(t)) L ∞ (ρ 0 , v 0 ) s-1 s L 2 (ρ 0 , v 0 ) 1 s Ḣs e C R t 0 ( ρ(τ ) 2 L ∞ + v(τ ) 2 L ∞ )dτ . (5)
It remained then the difficulty to obtain a bound of

t 0 ( ρ(τ ) 2 L ∞ + v(τ ) 2 L ∞
)dτ depending only on the initial data and not on the time. We overcome this difficulty thanks to Lemma 4.3 used with σ = r, σ = -r, for 0 < r < 1, a fractional Gagliardo-Nirenberg inequality [START_REF] Chen | Viscous quantum hydrodynamics and parameter-elliptic systems[END_REF] and inequality (4), to obtain under the assumption that ρ

L ∞ ≤ 3 4 and v L ∞ ≤ 1, t 0 ( ρ(τ ) 2 L ∞ + v(τ ) 2 L ∞ )dτ (ρ 0 , v 0 ) Ḣr (ρ 0 , v 0 ) Ḣ-r e C1 (ρ0,v0) 2 L 2 .
Therefore, from (5), we obtain,

(ρ(t), v(t)) L ∞ (ρ 0 , v 0 ) s-1 s L 2 (ρ 0 , v 0 ) 1 s Ḣs e C2 (ρ0,v0) Ḣr (ρ0,v0) Ḣ-r e C 1 (ρ 0 ,v 0 ) 2 L 2 ,
and then infer the condition on the initial data ensuring that ρ L ∞ ≤ 3 4 and v L ∞ ≤ 1 as soon as n 0 is sufficiently large relative to (n 0 -n 0 , u 0 ) in some sense, thus yielding to existence and uniqueness of global strong solutions of (1)-( 2). This paper is organized as follows: In section 2, we introduce some notations. In section 3, we establish some crucial estimates for the proof of Theorem 4.1. In section 4, we give the proof of the Lemmata mentioned previously to obtain the proof of Theorem 4.1. From the latter, we deduce our main Theorem 4.2 with the condition on the initial data (n 0 -n 0 , u 0 ) ensuring existence of global strong solutions of the Quantum Navier Stokes equations ( 1)-(2).

Some notations

We denote A B, the estimate A ≤ C B where C > 0 is a absolute constant. Given a function f which is periodic with period L, and thus representable as a function on the torus R 2 /LZ 2 , we define the discrete Fourier transform f : Z 2 -→ C by the formula,

f (k) = 1 L 2 R 2 /LZ 2 e -2π L ix•k f (x) dx, when f is absolutely integrable on R 2 /LZ 2 . If f ∈ L 2 (R 2 /LZ 2 )
, then from Parseval equality, we have,

k∈Z 2 | f (k)| 2 = 1 L L 0 |f (x)| 2 dx. Let s ∈ R, we define the Sobolev norm f H s (R 2 /LZ 2 ) of a tempered distribution f : R 2 /LZ 2 -→ R by, f H s (R 2 /LZ 2 ) = k∈Z 2 1 + 2π|k| L 2 s | f (k)| 2 1 2
, and then we denote by H s (R 2 /LZ 2 ) the space of tempered distributions with finite H s (R 2 /LZ 2 ) norm. On the torus R 2 /LZ 2 , for s > -1, we also define the homogeneous Sobolev norm,

f Ḣs (R 2 /LZ 2 ) = k∈Z 2 2π|k| L 2s | f (k)| 2 1 2
, and then we denote by Ḣs (R 2 /LZ 2 ) the space of tempered distributions with finite Ḣs (R 2 /LZ 2 ) norm.

We use the Fourier transform to define the fractional Laplacian operator (-∆) α , α > -1 on R 2 /LZ 2 and we define it as follows,

(-∆) α f (k) = 2π|k| L 2α f (k).
3 Some estimates

In this section, we give a serie of Lemmata which will be used in the next section. Let Ω = R 2 \LZ 2 with L > 0 a real. We begin with the following inequality proved in [START_REF] Chemin | Perfect Incompressible Fluids[END_REF] which states that for all f, g ∈ L ∞ (Ω) ∩ Ḣs (Ω),

f g Ḣs ( f L ∞ g Ḣs + f Ḣs g L ∞ ). ( 6
)
We have also the following inequality proved in [START_REF] Chemin | Perfect Incompressible Fluids[END_REF] for nonlinear composition which states that for s > 0, I an open subset of R with 0 ∈ I and for any real function f ∈ BC [s]+2 (I) such that f (0) = 0, we have for all u ∈ H s (Ω) such that u(x) ∈ I, f (u) ∈ H s (Ω). More precisely there exists a non decreasing continuous function C depending only on s and max

0≤k≤[s]+2 f (k) L ∞ (I) such that, f (u) H s ≤ C( u L ∞ ) u H s . (7)
Let us mention also the following Interpolation inequality: for all (α, β) such that 0

≤ α < 1 < β and for all u ∈ H β (Ω), u L ∞ (Ω) u γ Ḣα (Ω) u δ Ḣβ (Ω) , (8)
where γ > 0, δ > 0, γ + δ = 1 and γα + βδ = 1. Now, we give three crucial Lemmata.

Lemma 3.1 Let σ ∈ R, 0 < |σ| < 1, f, h ∈ Ḣσ (Ω) ∩ Ḣσ+1 (Ω) and g ∈ Ḣ1 (Ω). Then, we have, | f Dg, h Ḣσ | ( f Ḣσ ∇h Ḣσ + h Ḣσ ∇f Ḣσ ) ∇g L 2 ,
where D is a derivative of first order.

Proof. We set

A = -∆. Let δ ≥ 0 such that σ 2 + δ > 0, σ 2 + δ < 1 2 and σ 2 + 1 2 -δ > σ 2 that means -σ 2 < δ < min( 1 2 , 1 2 -σ 2 )
. Then, we have

| f Dg, h Ḣσ | = | A σ 2 (f Dg), A σ 2 h | = A σ 2 +δ-1 2 (f Dg), A σ 2 + 1 2 -δ h ≤ A σ 2 +δ-1 2 (f Dg) L 2 A σ 2 + 1 2 -δ h L 2 ,
where we have used Cauchy-Schwarz inequality. Since 0 < σ 2 + δ < 1 2 , we have (see [START_REF] Giga | Solutions in L r to the Navier-Stokes initial value problem[END_REF], see also (A5) in [START_REF] Kato | On the nonstationary Navier-Stokes system[END_REF], [START_REF] Fujita | On the Navier-Stokes initial value problem I[END_REF]),

A σ 2 +δ-1 2 (f Dg) L 2 A σ 2 +δ f L 2 A 1 2 g L 2 .
Then, we deduce,

| f Dg, h Ḣσ | A σ 2 +δ f L 2 A σ 2 + 1 2 -δ h L 2 A 1 2 g L 2 . ( 9
)
Thanks to the fractional Gagliardo-Nirenberg inequality given by Corollary 1.5 in [START_REF] Hajaiej | Necessary and sufficient conditions for the fractional Gargliardo-Nirenberg inequalities and applications to Navier-Stokes and generalized boson equations[END_REF], we have for any

u ∈ H 1 (Ω) and 0 ≤ θ ≤ 1, u Ḣθ u 1-θ L 2 u θ Ḣ1 . (10)
We use Inequality [START_REF] Dreher | The transient equations of viscous quantum hydrodynamics[END_REF] with u = A σ 2 ρ and θ = 2δ, then we get,

A σ 2 +δ f L 2 = A σ 2 f Ḣ2δ A σ 2 f 1-2δ L 2 A σ 2 f 2δ Ḣ1 = f 1-2δ Ḣσ ∇f 2δ Ḣσ . (11)
We use Inequality [START_REF] Dreher | The transient equations of viscous quantum hydrodynamics[END_REF] with u = A σ 2 h and θ = 1 -2δ to obtain,

A σ 2 + 1 2 -δ h L 2 = A σ 2 h Ḣ1-2δ h 2δ Ḣσ ∇h 1-2δ Ḣσ . (12)
Therefore, from (9), using [START_REF] Dong | A note on barotropic compressible quantum navier-Stokes equations[END_REF], [START_REF] Fujita | On the Navier-Stokes initial value problem I[END_REF] and Young inequality, we deduce,

| f Dg, h Ḣσ | ( f Ḣσ ∇h Ḣσ + h Ḣσ ∇f Ḣσ ) ∇g L 2 ,
which concludes the proof.

Lemma 3.2 Let σ ≥ 0, 0 < γ < 1, ∈ H σ (Ω) ∩ L ∞ (Ω) such that L ∞ (Ω)
≤ γ, then there exists a constant C > 0 depending only on σ and γ such that,

1 + Ḣσ ≤ C Ḣσ .
Proof. If σ = 0, the proof holds. Let us assume that σ > 0. Let λ > 0, λ (x) = (λx) and f the function defined by

f (x) = x 1 + x then f ∈ C ∞ ([-γ, γ]
) and λ L ∞ ≤ γ therefore thanks to (7), we deduce that there exists a constant C > 0 depending only on σ, γ such that,

λ 1 + λ H σ ≤ C λ H σ . (13)
From ( 13), we deduce that,

λ 1 + λ L 2 + λ 1 + λ Ḣσ λ L 2 + λ Ḣσ . (14)
Then, we get, 1

λ 1 + L 2 + λ σ λ 1 + Ḣσ 1 λ L 2 + λ σ λ Ḣσ . (15)
We multiply Inequality (15) by λ λ σ to obtain,

1 λ σ 1 + L 2 + 1 + Ḣσ 1 λ σ L 2 + Ḣσ . (16)
Then, taking the limit as λ → ∞ in ( 16), we conclude the proof.

Lemma 3.3 Let σ ∈ R, 0 < |σ| < 1, 0 < γ < 1, f, h ∈ Ḣσ (Ω) ∩ Ḣσ+1 (Ω), g ∈ Ḣ1 (Ω) and ∈ Ḣσ (Ω) ∩ H σ+1 (Ω) ∩ L ∞ (Ω) such that L ∞ ≤ γ. Then, we have, f 1 + Dg, h Ḣσ ( f Ḣσ ∇h Ḣσ + h Ḣσ ∇f Ḣσ ) ∇g L 2 + f L ∞ ( Ḣσ ∇h Ḣσ + h Ḣσ ∇ Ḣσ ) ∇g L 2 ,
where D is a derivative of first order.

Proof. We set A = -∆. Let δ ≥ 0 such that σ 2 + δ > 0, σ 2 + δ < 1 2 and σ 2 + 1 2 -δ > σ that means -σ 2 < δ < min( 1 2 , 1 2 -σ 2 )
. Similarly as (9), we have,

f 1 + Dg, h Ḣσ A σ 2 +δ f 1 + L 2 A σ 2 + 1 2 -δ h L 2 A 1 2 g L 2 .
From ( 12), we get,

A σ 2 + 1 2 -δ h L 2 h 2δ Ḣσ ∇h 1-2δ Ḣσ .
Then, we deduce,

f 1 + Dg, h Ḣσ A σ 2 +δ f 1 + L 2 h 2δ Ḣσ ∇h 1-2δ Ḣσ ∇g L 2 . ( 17
) Since f 1 + = f -f 1 + , then A σ 2 +δ f 1 + L 2 = f 1 + Ḣσ+2δ ≤ f Ḣσ+2δ + f 1 + Ḣσ+2δ .
Since σ + 2δ > 0, thanks to (6), we have,

f 1 + Ḣσ+2δ f Ḣσ+2δ 1 + L ∞ + f L ∞ 1 + Ḣσ+2δ ≤ f Ḣσ+2δ L ∞ 1 - L ∞ + f L ∞ 1 + Ḣσ+2δ ≤ γ 1 -γ f Ḣσ+2δ + f L ∞ 1 + Ḣσ+2δ , (18)
where we have used the fact that L ∞ ≤ γ. Then, we deduce,

A σ 2 +δ f 1 + L 2 f Ḣσ+2δ + f L ∞ 1 + Ḣσ+2δ . ( 19
)
Thanks to Lemma 3.2, we have,

From [START_REF] Giga | Solutions in L r to the Navier-Stokes initial value problem[END_REF], we infer,

A σ 2 +δ f 1 + L 2 f Ḣσ+2δ + f L ∞ Ḣσ+2δ . (20)
Using the same arguments as in [START_REF] Dong | A note on barotropic compressible quantum navier-Stokes equations[END_REF], we get,

f Ḣσ+2δ f 1-2δ Ḣσ ∇f 2δ Ḣσ Ḣσ+2δ 1-2δ Ḣσ ∇ 2δ Ḣσ (21)
From [START_REF] Gamba | Global existence of solutions to one-dimensional viscous quantum hydrodynamic equations[END_REF], using [START_REF] Gardner | Dispersive/Hyperbolic Hydrodynamic Models for Quantum Transport (in Semiconductor Devices)[END_REF], [START_REF] Grant | Pressure and stress tensor expressions in the fluid mechanical formulation of the Bose condensate equations[END_REF] and Young inequalities, we deduce,

f 1 + Dg, h Ḣσ ( f Ḣσ ∇h Ḣσ + h Ḣσ ∇f Ḣσ ) ∇g L 2 + f L ∞ ( Ḣσ ∇h Ḣσ + h Ḣσ ∇ Ḣσ ) ∇g L 2 , (22)
which concludes the proof.

Global regularity

This section is devoted to the proof of Theorem 4.2. We assume that n 0 > 0 and we introduce m 0 = n 0 u 0 + ε∇n 0 . We begin by noting that the system of equations ( 1) is equivalent to a system of equations of type [START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF]. Indeed, thanks to Theorem 2.1 in [START_REF] Jüngel | Effective velocity in Navier-Stokes equations with third-order derivatives[END_REF], under the assumption that n > 0, if (n, u) is solution of the shallow water equations (1) for the initial data (n 0 , u 0 ) then (n, m) is solution of the system of equations ( 23) for the initial (n 0 , m 0 ), with m = nu + ε∇n. Thanks again to Theorem 2.1 in [START_REF] Jüngel | Effective velocity in Navier-Stokes equations with third-order derivatives[END_REF], under the assumption that n > 0, if (n, m) is solution of the system of equations [START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF] for the initial data (n 0 , m 0 ) then (n, u) is solution of the shallow water equations [START_REF] Adams | Sobolev spaces[END_REF] for the initial data (n 0 , u 0 ) with u = m -ε∇n n . Then, we consider the following system of equations,

     ∂n ∂t + ∇ • m -ε∆n = 0, ∂m ∂t + ∇ • m n ⊗ m -ε∆m + gn∇n = 0. ( 23
)
We introduce ñ defined by ñ ≡ n -n 0 n 0 and m ≡ m n 0 , from [START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF], we obtain the following equivalent system of equations,

     ∂ ñ ∂t + ∇ • m -ε∆ñ = 0, ∂ m ∂t + ∇ • m 1 + ñ ⊗ m -ε∆ m + gn 0 ñ∇ñ + gn 0 ∇ñ = 0. ( 24
)
For any λ ∈ R * , by using the rescaled solutions ñλ (x, t) = ñ(λx, λ 2 t) and mλ (x, t) = λ m(λx, λ 2 t), we obtain the equivalent system of equations,

     ∂ ñλ ∂t + ∇ • mλ -ε∆ñ λ = 0, ∂ mλ ∂t + ∇ • mλ 1 + ñλ ⊗ mλ -ε∆ mλ + gn 0 λ 2 ñλ ∇ñ λ + gn 0 λ 2 ∇ñ λ = 0. ( 25
) Setting λ = 1 √ gn 0 in (25), then with ρ(x, t) = n x √ gn0 , t gn0 -n 0 n 0 and v(x, t) = 1 n 0 √ gn 0 m x √ gn 0 , t gn 0 ,
provided that ρ > -1, we deduce that the system of equations [START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF] on O × [0, T 0 ] is equivalent to the following system of equations on Ω × [0, T ] with

Ω def = √ gn 0 O and T = gn 0 T 0 ,      ∂ρ ∂t + ∇ • v -ε∆ρ = 0, ∂v ∂t + ∇ • v 1 + ρ ⊗ v -ε∆v + ρ∇ρ + ∇ρ = 0. ( 26
)
with initial conditions,

ρ(x, 0) = ρ 0 (x), v(x, 0) = v 0 (x), (27)
where

ρ 0 (x) = n 0 x √ gn0 -n 0 n 0 (28)
and

v 0 (x) = 1 n 0 √ gn 0 m 0 x √ gn 0 . (29)
Then, in what follows, we study the following system of equations,

     ∂ρ ∂t + ∇ • v -ε∆ρ = 0, ∂v ∂t + ∇ • v 1 + ρ ⊗ v -ε∆v + (1 + ρ)∇ρ = 0. ( 30
)
with initial conditions,

ρ(x, 0) = ρ 0 (x), v(x, 0) = v 0 (x). ( 31
)
We thus establish the following energy estimate in H s for any s > 1.

Lemma 4.1 Let 0 < δ < 1, ω 0 ≡ (ρ 0 , v 0 ) ∈ H s (Ω) 3 , s ≥ 0. If ω ≡ (ρ, v) ∈ C([0, T ]; H s (Ω) 3 ) ∩ C(]0, T ]; H s+1 (Ω) 3 ) with ρ L ∞ (Ω×[0,T ])
≤ δ is a solution of the system of Equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF], then there exists a constant c 1 > 0 depending only on δ, s such that for all τ ∈ [0, T ],

ω(τ ) 2 Ḣs + ε 2 τ 0 ω(σ) 2 Ḣs+1 dσ ≤ ω 0 2 Ḣs e c 1 ε R τ 0 a(σ)dσ ,
where

a(σ) = ω(σ) 2 L ∞ (1 + v(σ) 2 L ∞ ).
Proof. We take the inner product in Ḣs (Ω) of the first equation of ( 30) with ρ, use integrations by parts to obtain, 1 2

d dt ρ 2 Ḣs + ε ∇ρ 2 Ḣs = v, ∇ρ Ḣs . (32)
Now, we take the inner product in Ḣs (Ω) 2 of the second equation of ( 30) with v, use integrations by parts to obtain, 1 2

d dt v 2 Ḣs + ε ∇v 2 Ḣs = v 1 + ρ ⊗ v, ∇v Ḣs + ρ 2 2 , ∇ • v Ḣs -∇ρ, v Ḣs . ( 33
)
Thanks to Cauchy-Schwarz inequality and Young inequality, we have,

v 1 + ρ ⊗ v, ∇v Ḣs ≤ v 1 + ρ ⊗ v Ḣs ∇v Ḣs ≤ 1 2ε v 1 + ρ ⊗ v 2 Ḣs + ε 2 ∇v 2 Ḣs and ρ 2 2 , ∇ • v Ḣs ≤ ρ 2 Ḣs ∇ • v Ḣs 2 ≤ 1 2ε ρ 2 2 Ḣs + ε 8 ∇ • v 2 Ḣs ≤ 1 2ε ρ 2 2 Ḣs + ε 4 ∇v 2 Ḣs .
Then, we deduce, 1 2

d dt v 2 Ḣs + ε 4 ∇v 2 Ḣs ≤ 1 2ε v 1 + ρ ⊗ v 2 Ḣs + 1 2ε ρ 2 2 Ḣs -∇ρ, v Ḣs . ( 34
)
We sum Equations ( 32) and [START_REF] Markowich | Semiconductor Equations[END_REF], to obtain,

1 2 d dt ( ρ 2 Ḣs + v 2 Ḣs ) + ε ∇ρ 2 Ḣs + ε 4 ∇v 2 Ḣs ≤ 1 2ε v 1 + ρ ⊗ v 2 Ḣs + 1 2ε ρ 2 2 Ḣs . (35)
Thanks to (6), we get,

ρ 2 Ḣs ρ L ∞ ρ Ḣs (36)
and also

v 1 + ρ ⊗ v Ḣs v 1 + ρ Ḣs v L ∞ + v 1 + ρ L ∞ v Ḣs . (37)
Furthermore, we have,

v 1 + ρ Ḣs ≤ v Ḣs + vρ 1 + ρ Ḣs . ( 38
)
Thanks to (6), we have also,

vρ 1 + ρ Ḣs v L ∞ ρ 1 + ρ Ḣs + v Ḣs ρ 1 + ρ L ∞ . ( 39
)
Then with (38) and (39), we deduce,

v 1 + ρ Ḣs v Ḣs 1 + ρ 1 + ρ L ∞ + v L ∞ ρ 1 + ρ Ḣs . (40)
Therefore, using (40) and (37), we get,

v 1 + ρ ⊗ v Ḣs 1 + ρ 1 + ρ L ∞ v Ḣs v L ∞ + v 2 L ∞ ρ 1 + ρ Ḣs + v 1 + ρ L ∞ v Ḣs ≤ 2 1 -δ v Ḣs v L ∞ + ρ 1 + ρ Ḣs v 2 L ∞ , (41)
where we have used the fact that ρ L ∞ ≤ δ < 1. Then, using (41) and [START_REF] Wyatt | Quantum Dynamics with Trajectories[END_REF], from [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF], we deduce,

1 2 d dt ( ρ 2 Ḣs + v 2 Ḣs ) + ε ∇ρ 2 Ḣs + ε 4 ∇v 2 Ḣs 1 ε v 2 Ḣs v 2 L ∞ + 1 ε ρ 1 + ρ 2 Ḣs v 4 L ∞ + 1 2ε ρ 2 L ∞ ρ 2 Ḣs . (42)
Thanks to Lemma 3.2, from (42), we deduce that for all t ∈]0, T ],

1 2 d dt ( ρ(t) 2 Ḣs + v(t) 2 Ḣs ) + ε ∇ρ(t) 2 Ḣs + ε 4 ∇v(t) 2 Ḣs 1 ε v(t) 2 Ḣs v(t) 2 L ∞ + 1 ε ρ(t) 2 Ḣs v(t) 4 L ∞ + 1 2ε ρ(t) 2 L ∞ ρ(t) 2 Ḣs
(43) which implies,

1 2 d dτ ( ρ(τ ) 2 Ḣs + v(τ ) 2 Ḣs ) +ε ∇ρ(τ ) 2 Ḣs + ε 4 ∇v(τ ) 2 Ḣs 1 ε ( v(τ ) 2 L ∞ + v(τ ) 4 L ∞ + ρ(τ ) 2 L ∞ )( ρ(τ ) 2 Ḣs + v(τ ) 2 Ḣs ). (44)
Using Gronwall inequality, from (44), we deduce that there exists a constant C > 0 depending only on δ, s such that for all τ ∈ [0, T ],

ρ(τ ) 2 Ḣs + v(τ ) 2 Ḣs + τ 0 2ε ∇ρ(σ) 2 Ḣs + ε 2 ∇v(σ) 2 Ḣs dσ ≤ ( ρ 0 2 Ḣs + v 0 2 Ḣs )e C R τ 0 a(σ)dσ , where a(σ) = 1 ε ( v(σ) 2 L ∞ + v(σ) 4 L ∞ + ρ(σ) 2 L ∞ ).
Thanks to Plancherel Theorem, we have ∇ρ(τ ) Ḣs = ρ(τ ) Ḣs+1 and also ∇v(τ ) Ḣs = v(τ ) Ḣs+1 , then we conclude the proof.

In the following Proposition, we prove existence and uniqueness of local strong solutions of [START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF].

Proposition 4.1 Let ω 0 ≡ (ρ 0 , v 0 ) ∈ Ḣ-r (Ω) 3 ∩ H s (Ω) 3 with 0 ≤ r < 1 < s, ρ 0 L ∞ < 1.
Then there exists a maximal time of existence T * > 0 such that there exists a unique strong solution

ω ≡ (ρ, v) ∈ C([0, T * [; H -r (Ω) 3 ∩ H s (Ω) 3 ) ∩ C 1 (]0, T * [; H s-1 (Ω) 3 ) ∩ C(]0, T * [; H s+1 (Ω)
3 ) of the system of Equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF]. Moreover if T * < ∞, then either ρ L ∞ (Ω×[0,T * ]) ≥ 1 or either

T * 0 v(τ ) 4 L ∞ dτ = ∞. (45)
Proof.

Let a = 3 + ρ 0 L ∞ 2(1 + ρ 0 L ∞ ) > 1 and δ = 1 2 (1 + ρ 0 L ∞ ) < 1.
We introduce χ(λ) a smooth bump function with values in the interval [0, 1], identically equal to one for -1 ≤ λ ≤ 1 and identically equal to zero for |λ| ≥ a.

We notice for all λ δ ≥ a, 1 + λχ λ δ = 1, and for all λ δ < a,

1 + λχ λ δ ≥ 1 -|λ| ≥ 1 -δa = 1 4 (1 -ρ 0 L ∞ ) > 0.
Then, we get for all λ ∈ R,

1 + λχ λ δ ≥ 1 4 (1 -ρ 0 L ∞ ) > 0. ( 46
)
For the proof, we use some results which deal with existence, uniqueness, regularity of solutions ω = (ρ, v) for nonlinear evolution equations of the form

∂ t ω = Aω + f (ω), (47)
with initial conditions ω(0) = ω 0 .

More precisely, we use Proposition 2.1 in [START_REF] Ball | Remarks on blow-up and nonexistence theorems for nonlinear evolution equations[END_REF] with X = H s-1 (Ω) 3 , A = ε∆ for our generator of holomorphic semigroup T (t) = e -tA of bounded linear operators on X and f our locally Lipschitz continuous function on

X α = H s (Ω) 3 , α = 1 2 defined by f (ω) = (f 1 (ω), f 2 (ω)) with f 1 (ω) = -∇ • v and f 2 (ω) = -∇ • v 1 + ρχ( ρ δ ) ⊗ v -ρ∇ρ -∇ρ. (49)
Indeed, f is locally Lipschitz continuous on H s (Ω) 3 , since for all ω 1 = (ρ 1 , v 1 ) ∈ H s (Ω) 3 and ω 2 = (ρ 2 , v 2 ) ∈ H s (Ω) 3 , we have,

∇ • v 1 -∇ • v 2 H s-1 ≤ v 1 -v 2 H s ∇ρ 1 -∇ρ 2 H s-1 ≤ ρ 1 -ρ 2 H s ρ 1 ∇ρ 1 -ρ 2 ∇ρ 2 H s-1 = 1 2 ∇((ρ 1 -ρ 2 )(ρ 1 + ρ 2 )) H s-1 ≤ 1 2 (ρ 1 -ρ 2 )(ρ 1 + ρ 2 ) H s ρ 1 -ρ 2 H s ( ρ 1 H s + ρ 2 H s ), (50)
where we have used the fact that for all h 1 ∈ H s (Ω), h 2 ∈ H s (Ω) with s > 1 (see [START_REF] Chemin | Perfect Incompressible Fluids[END_REF]),

h 1 h 2 H s h 1 H s h 2 H s , (51)
We introduce the function g ∈ C ∞ defined by g(λ

) = 1 1 + λχ( λ δ )
, thanks to (46) we notice that for all k ∈ N, there exists a real C k > 0 depending only on δ such that

g (k) L ∞ ≤ C k . (52)
Then, we get,

∇ • (g(ρ 1 )v 1 ⊗ v 1 ) -∇ • (g(ρ 2 )v 2 ⊗ v 2 ) H s-1 ≤ g(ρ 1 )v 1 ⊗ v 1 -g(ρ 2 )v 2 ⊗ v 2 H s ≤ (g(ρ 1 )v 1 -g(ρ 2 )v 2 ) ⊗ v 1 H s + g(ρ 2 )v 2 ⊗ (v 1 -v 2 ) H s ≤ (v 1 -v 2 ) ⊗ v 1 g(ρ 1 ) H s + v 2 ⊗ v 1 (g(ρ 1 ) -g(ρ 2 )) H s + g(ρ 2 )v 2 ⊗ (v 1 -v 2 ) H s . (53)
Using (51), we estimate each term on the right hand side of Inequality (53),

(v 1 -v 2 ) ⊗ v 1 g(ρ 1 ) H s v 1 -v 2 H s v 1 H s (1 + g(ρ 1 ) -1 H s ) v 2 ⊗ v 1 (g(ρ 1 ) -g(ρ 2 )) H s v 1 H s v 2 H s g(ρ 1 ) -g(ρ 2 ) H s g(ρ 2 )v 2 ⊗ (v 1 -v 2 ) H s v 1 -v 2 H s v 2 H s (1 + g(ρ 2 ) -1 H s ). ( 54
)
Thanks to ( 7) and (52), we deduce that there exists a non-decreasing function C > 0 depending only on s and δ such that,

g(ρ 1 ) -1 H s C( ρ 1 L ∞ ) ρ 1 H s g(ρ 2 ) -1 H s C( ρ 2 L ∞ ) ρ 2 H s .
Furthermore, using Taylor formula at order 1, we get g(ρ 2 ) -g(ρ 1 ) = (ρ 2 -ρ 1 )

1 0 g ((1 -σ)ρ 1 + σρ 2 )dσ, then we deduce, g(ρ 1 ) -g(ρ 2 ) H s ρ 1 -ρ 2 H s 1 0 g ((1 -σ)ρ 1 + σρ 2 ) H s dσ ρ 1 -ρ 2 H s 1 0 C 1 ( (1 -σ)ρ 1 + σρ 2 L ∞ ) (1 -σ)ρ 1 + σρ 2 ) H s dσ,
where we have used ( 7) and C 1 > 0 is a non-decreasing function depending only on s and δ. Therefore, we get,

g(ρ 1 ) -g(ρ 2 ) H s ρ 1 -ρ 2 H s C 1 ( ρ 1 L ∞ + ρ 2 L ∞ )( ρ 1 H s + σρ 2 ) H s ).
Thanks to (54) and the Sobolev embedding H s (Ω) → L ∞ (Ω) since s > 1, from (53), we deduce therefore,

∇ • (g(ρ 1 )v 1 ⊗ v 1 ) -∇ • (g(ρ 2 )v 2 ⊗ v 2 ) H s C 2 ( ρ 1 H s , ρ 2 H s , v 1 H s , v 2 H s )( v 1 -v 2 H s + ρ 1 -ρ 2 H s ), ( 55
)
where C 2 > 0 is a continuous function on R 4 . Using (50) and (55), we get that for all ω 1 ∈ H s (Ω) 3 ,

ω 2 ∈ H s (Ω) 3 , f (ω 1 ) -f (ω 2 ) H s-1 C 3 ( ω 1 H s , ω 2 H s ) ω 1 -ω 2 H s , (56)
where C 3 > 0 is a continuous function on Ω, which proves that f is well locally Lipschitz continuous on H s (Ω) 3 . Then, we deduce thanks to Proposition 2.1 combined with Theorem 3.1 in [START_REF] Ball | Remarks on blow-up and nonexistence theorems for nonlinear evolution equations[END_REF], that there exists a maximal time T > 0 such that there exists an unique solution ω

= (ρ, v) ∈ C([0, T [; H s (Ω) 3 ) ∩ C 1 (]0, T [; H s-1 (Ω) 3 ) of the system of Equations (47)-(48). Moreover if T < ∞ then lim t→T * ω(t) H s = ∞. From (56), we get that for all w ∈ H s (Ω) 3 , f (w) H s-1 C 3 ( w H s , 0) w H s , then we deduce from (47) that ω ∈ C(]0, T [; H s+1 (Ω)).
Then, we write ω under its integral form,

ω(t) = e -εt∆ ω 0 + t 0 e -ε(t-σ)∆ f (ω(σ))dσ. ( 57
)
Since s > 1 -r > 0, then by Interpolation inequality, we have H s (Ω) → Ḣ1-r (Ω), then using the same arguments as for (56), we get that for all w ∈ Ḣ-r ∩ H s , f (w) Ḣ-r C 3 ( w H s , 0) w H s . Since ω 0 ∈ H -r (Ω), therefore, from (57), we deduce also that ω ∈ C([0, T [;

H -r (Ω)). Since ρ 0 L ∞ < 1 then ρ 0 L ∞ < δ, moreover ρ ∈ C([0, T [; L ∞ (Ω)) due to the Sobolev embedding H s (Ω) → L ∞ (Ω)
, hence we deduce that there exists a time 0 < T < T such that ρ L ∞ (Ω×[0,T ]) ≤ δ. Since ρ L ∞ (Ω×[0,T ]) ≤ δ then we get χ( ρ δ ) = 1 on [0, T ] and from (49), we deduce in fact that ω is the unique solution of ( 30)-( 31) on [0, T ]. Then we deduce that there exists a maximal time of existence 0 < T * < ∞ such that there exists an unique solution ω

= (ρ , v ) ∈ C([0, T * [; H -r ∩ H s (Ω)) ∩ C 1 (]0, T * [; H s-1 (Ω)) ∩ C(]0, T * [; H s+1 (Ω)
) of the system of equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF]. We get that if T * < ∞, then either ρ L ∞ (Ω×[0,T * ]) ≥ 1 or either lim t→T * ω (t) H s = ∞ (by using the same arguments as ω) and in this case thanks to Lemma 4.1 and the fact that for any 0 < T < T * ,

T 0 v(τ ) 2 L ∞ dτ ≤ √ T T 0 v(τ ) 4 L ∞ dτ 1 2
, we get

T * 0 v(τ ) 4 L ∞ dτ = ∞
, then, we conclude the proof.

We derive now a Bresch-Desjardin Entropy type, as in [START_REF] Jüngel | Effective velocity in Navier-Stokes equations with third-order derivatives[END_REF].

Lemma 4.2 Let ω 0 ≡ (ρ 0 , v 0 ) ∈ L 2 (Ω) 3 with ρ 0 ∈ L ∞ (Ω) such that ρ 0 L ∞ ≤ 3 4 . If ω ≡ (ρ, v) ∈ C([0, T ]; L 2 (Ω)) 3 ∩ C 1 (]0, T ]; L 2 (Ω) 3 ) with ρ ∈ L ∞ (Ω × [0, T ]) such that ρ L ∞ (Ω×[0,T ]) ≤ 3
4 is a solution of the system of Equations ( 30)-( 31), then we have for all t ∈ [0, T ],

1 4 ω(t) 2 L 2 + ε t 0 1 4 ∇ v(τ ) 1 + ρ(τ ) 2 L 2 + ∇ρ(τ ) 2 L 2 dτ ≤ 2 ω 0 2 L 2 .
Proof. We can write the system of equations ( 30) as follows,

     ∂η ∂t + ∇ • (ηw) -ε∆η = 0, ∂(ηw) ∂t + ∇ • (ηw ⊗ w) -ε∆(ηw) + ∇p(η) = 0, (58)
where

w = v 1 + ρ , η = 1 + ρ and p(η) = η 2 2
. By using the enthalpy h(η) = η -1, we notice ηh (η) = p (η), then from Equation ( 15) of [START_REF] Jüngel | Effective velocity in Navier-Stokes equations with third-order derivatives[END_REF], we obtain,

d dt Ω 1 2 η|w| 2 + H(η) + Ω (εη|∇w| 2 + |∇η| 2 ) = 0,
where H(η) = η 1 h(τ )dτ , which yields to,

d dt Ω |v| 2 2(1 + ρ) + ρ 2 2 + ε Ω (1 + ρ) ∇ v 1 + ρ 2 + |∇ρ| 2 = 0.
We integrate Equation just above on [0, t] to obtain, for all t ∈ [0, T ],

1 2 v(t) 1 + ρ(t) 2 L 2 + 1 2 ρ(t) 2 L 2 + ε t 0 1 + ρ(τ )∇ v(τ ) 1 + ρ(τ ) 2 L 2 + ∇ρ(τ ) 2 L 2 dτ = 1 2 v 0 √ 1 + ρ 0 2 L 2 + 1 2 ρ 0 2 L 2 . Since ρ L ∞ (Ω×[0,T ]) ≤ 3 4 and ρ 0 L ∞ ≤ 3 4 , then we deduce that for all t ∈ [0, T ], 1 4 v(t) 2 L 2 + 1 2 ρ(t) 2 L 2 + ε t 0 1 4 ∇ v(τ ) 1 + ρ(τ ) 2 L 2 + ∇ρ(τ ) 2 L 2 dτ ≤ 2 v 0 2 L 2 + 1 2 ρ 0 2 L 2 ,
which concludes the proof.

The following Lemma will help us to express

t 0 (ρ, v)(σ) 2 L ∞ dσ in terms of t 0 ( ∇ρ(τ ) 2 L 2 + ∇ v(τ ) 1+ρ(τ ) 2 L 2
)dτ . 3 is a solution of the system of Equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF] with ω L ∞ (Ω×[0,T ]) 3 ≤ 3 4 and such that T 0 a(τ )dτ < ∞, then there exists a constant C > 0 such that for all s ∈ [0, T ],

Lemma 4.3 Let ω 0 ≡ (ρ 0 , v 0 ) ∈ Ḣσ (Ω) 3 ∩ L 2 (Ω) 3 , σ ∈ R, 0 < |σ| < 1. If ω ≡ (ρ, v) ∈ C([0, T ]; Ḣσ (Ω) 3 ) ∩ C(]0, T ]; H 1+σ (Ω) 3 ) ∩ L ∞ (Ω × [0, T ])
ω(s) 2 Ḣσ + ε s 0 ω(τ ) 2 Ḣσ+1 dτ ≤ ω 0 2 Ḣσ e C R t 0 a(τ )dτ ,
where a(τ

) = ∇ρ(τ ) 2 L 2 + ∇ v(τ ) 1+ρ(τ ) 2 L 2 ε .
Proof.

We take the inner product in Ḣσ (Ω) of the first equation of ( 30) with ρ, use integrations by parts to obtain, 1 2

d dt ρ 2 Ḣσ + ε ∇ρ 2 Ḣσ = v, ∇ρ Ḣσ . (59)
Now, we take the inner product in Ḣσ (Ω) of the second equation of ( 30) with v, to obtain, 1 2

d dt v 2 Ḣσ + ε ∇v 2 Ḣσ = -∇ • v 1 + ρ ⊗ v , v Ḣσ -ρ∇ρ, v Ḣσ -∇ρ, v Ḣσ . (60)
We sum Equations ( 59) and (60), to obtain,

1 2 d dt ( ρ 2 Ḣσ + v 2 Ḣσ ) + ε ∇ρ 2 Ḣσ + ε ∇v 2 Ḣσ = -∇ • v 1 + ρ ⊗ v , v Ḣσ -ρ∇ρ, v Ḣσ . (61)
Thanks to Lemma 3.1, we get,

-ρ∇ρ, v Ḣσ ( ρ Ḣσ ∇v Ḣσ + v Ḣσ ∇ρ Ḣσ ) ∇ρ L 2 . ( 62
)
We estimate now the first term at the right hand side of Equation (61). We notice

∇ • v 1+ρ ⊗ v = v 1+ρ • ∇ v + ∇ • v 1+ρ v, then we deduce, -∇ • v 1 + ρ ⊗ v , v Ḣσ = - v 1 + ρ • ∇ v, v Ḣσ -∇ • v 1 + ρ v, v Ḣσ .
Thanks to Lemma 3.1, we get,

-∇ • v 1 + ρ v, v Ḣσ v Ḣσ ∇v Ḣσ ∇ v 1 + ρ L 2 . (63)
Thanks to Lemma 3.3, we have,

- v 1 + ρ • ∇ v, v Ḣσ v Ḣσ ∇v Ḣσ ∇v L 2 + v L ∞ ( ρ Ḣσ ∇v Ḣσ + v Ḣσ ∇ρ Ḣσ ) ∇v L 2 . ( 64
)
Using inequalities (62)-(64) and Young inequalities, from (61), we deduce that for all t ∈ [0, T ],

1 2 d dt ( ρ(t) 2 Ḣσ + v(t) 2 Ḣσ ) + ε 2 ∇ρ(t) 2 Ḣσ + ε 2 ∇v(t) 2 Ḣσ 1 ε ρ(t) 2 Ḣσ ∇ρ(t) 2 L 2 + 1 ε v(t) 2 Ḣσ ∇ρ(t) 2 L 2 + 1 ε v(t) 2 Ḣσ ∇v(t) 2 L 2 + ∇ v(t) 1 + ρ(t) 2 L 2 + 1 ε v(t) 2 L ∞ ρ(t) 2 Ḣσ ∇v(t) 2 L 2 + v(t) 2 L ∞ ε v(t) 2 Ḣσ ∇v(t) 2 L 2 . ( 65
)
Using the fact that v L ∞ (Ω×[0,T ]) 2 ≤ 3 4 , from (65), we infer that there exists a constant C > 0 depending only on σ such that,

1 2 d ds ( ρ(s) 2 Ḣσ + v(s) 2 Ḣσ ) + ε 2 ∇ρ(s) 2 Ḣσ + ε 2 ∇v(s) 2 Ḣσ ≤ C ∇ρ(s) 2 L 2 + ∇ v(s) 1+ρ(s) 2 L 2 + ∇v(s) 2 L 2 ε ( ρ(s) 2 Ḣσ + v(s) 2 Ḣσ ), ( 66
) which implies, 1 2 d ds ( ρ(s) 2 Ḣσ + v(s) 2 Ḣσ ) ≤ C ∇ρ(s) 2 L 2 + ∇ v(s) 1+ρ(s) 2 L 2 + ∇v(s) 2 L 2 ε ( ρ(s) 2 Ḣσ + v(s) 2 Ḣσ ).
(67) Then, using Gronwall Inequality, from (67), we deduce for all s ∈ [0, T ],

ρ(s) 2 Ḣσ + v(s) 2 Ḣσ ≤ ( ρ 0 2 Ḣσ + v 0 2 Ḣσ )e 2C R t 0 b(τ )dτ , (68)
where b(τ

) = ∇ρ(τ ) 2 L 2 + ∇ v(τ ) 1+ρ(τ ) 2 L 2 + ∇v(τ ) 2 L 2
ε . After integrating (66) and using (68), we deduce that for all s ∈ [0, T ],

ρ(s) 2 Ḣσ + v(s) 2 Ḣσ + s 0 ε ∇ρ(τ ) 2 Ḣσ + ε ∇v(τ ) 2 Ḣσ dτ ≤ ( ρ 0 2 Ḣσ + v 0 2 Ḣσ )e 2C R t 0 b(τ )dτ . (69)
We notice that ∇v(τ Then, gathering the previous results, we prove our Theorem 4.1 which deals with existence and uniqueness of global strong solutions of [START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF].

) = ∇ v(τ ) 1 + ρ(τ ) (1 + ρ(τ )) = ∇ v(τ ) 1 + ρ(τ ) (1 + ρ(τ )) + v(τ ) 1 + ρ(τ ) ⊗ ∇ρ(τ ), then we get, |∇v(τ )| ≤ ∇ v(τ ) 1 + ρ(τ ) (1 + |ρ(τ )|) + |v(τ )| 1 -|ρ(τ )| |∇ρ(τ )|.
Theorem 4.1 Let ω 0 ≡ (ρ 0 , v 0 ) ∈ Ḣ-r (Ω) 3 ∩ H s (Ω) 3 with 0 < r < 1 < s and ω 0 L ∞ ≤ 1 2 .
Then there exists a real α 1 > 0 depending only on r, s, v 0 , ρ 0 such if α 1 ≤ 1 2 , then there exists a unique global strong

solution ω ≡ (ρ, v) ∈ C([0, ∞[; Ḣ-r (Ω) 3 ∩ H s (Ω) 3 ) ∩ C 1 (]0, ∞[; H s-1 (Ω) 3 ) ∩ C(]0, ∞[; H s+1 (Ω)
3 ) of the system of Equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF] with 3 and the real α 1 is given by,

ω L ∞ (Ω×[0,+∞[) 3 ≤ 1 2 . Furthermore, ω ∈ C ∞ (Ω×]0, ∞[)
α 1 = C 1 ω 0 s-1 s L 2 ω 0 1 s Ḣs e c 1 ε α0 ,
where,

α 0 = ω 0 Ḣr ω 0 Ḣ-r e C ε 2 ω0 2 L 2 ,
with c 1 > 0, C > 0 and C 1 > 0 constants depending only on r, s. Moreoever, for all t ≥ 0,

t 0 ω(τ ) 2 L ∞ (Ω) dτ ≤ c 1 α 0 .
Proof.

Thanks to Proposition 4.1, there exists a maximal time of existence T * > 0 such that there exists a unique strong solution ω

≡ (ρ, v) ∈ C([0, T * [; Ḣ-r (Ω) 3 ∩ H s (Ω) 3 ) ∩ C 1 (]0, T * [; H s-1 (Ω) 3 ) ∩ C(]0, T * [; H s+1 (Ω) 3 )
of the system of Equations ( 30)- [START_REF] Kato | Commutator Estimates and the Euler and Navier-Stokes Equations[END_REF]. Moreover if T * < ∞, then either ρ L ∞ (Ω×[0,T * ]) ≥ 1 or either

T * 0 v(τ ) 4 L ∞ = ∞. (70)
Let us assume that T * < ∞.

Since s > 1, we get the Sobolev embedding H s (Ω) → L ∞ (Ω) and then ω ∈ C([0, T * [; L ∞ (Ω) 3 ). Furthermore, since ω 0 L ∞ < 3 4 , then we deduce that there exists a maximal time 0 < T ≤ T * such that for all t ∈ [0, T [, ω(t) L ∞ < 3 4 . Let us show that T = T * . Indeed if T < T * , then we get,

ω(T ) L ∞ = 3 4 . ( 71
)
Thanks to Lemma 4.1 and the fact that v L ∞ (Ω×[0,T ]) ≤ 1 , there exists a constant c 1 > 0 depending only on s such that for all t ∈ [0, T ],

ω(t) 2 Ḣs + ε 2 t 0 ω(τ ) 2 Ḣs+1 dτ ≤ ω 0 2 Ḣs e c 1 ε R t 0 ω(τ ) 2 L ∞ dτ , (72)
Thanks to (8), we have for a.e t ∈ [0, T ],

ω(t) L ∞ ω(t) 1 2
Ḣ1+r ω(t)

1 2
Ḣ1-r .

Thanks to Cauchy-Schwarz inequality, we infer,

t 0 ω(τ ) 2 L ∞ dτ t 0 ω(τ ) 2 Ḣ1+r dτ 1 2 t 0 ω(τ ) 2 Ḣ1-r (Ω) dτ 1 2
.

Thanks to Lemma 4.3 used first with σ = r and after with σ = -r, we deduce,

t 0 ω(τ ) 2 L ∞ dτ ω 0 Ḣr ω 0 Ḣ-r e C R t 0 b(τ )dτ , (73)
where C > 0 is a constant and b(τ

) = ∇ρ(τ ) 2 L 2 + ∇ v(τ ) 1+ρ(τ ) 2 L 2 ε
. Thanks to (73) and Lemma 4.2, we get that there exists a constant c 2 > 0 depending only on r, s such that for all t ∈ [0, T ],

t 0 ω(τ ) 2 L ∞ dτ ≤ c 2 α 0 ,
where,

α 0 = ω 0 Ḣr ω 0 Ḣ-r e 4C ε 2 ω0 2 L 2 .
Then, from (72), we deduce that for all t ∈ [0, T ],

ω(t) Ḣs ≤ ω 0 Ḣs e c 1 ε c2α0 . (74)
Thanks again to Lemma 4.2, we have for all t ∈ [0, T ],

ω(t) L 2 ω 0 L 2 . ( 75
)
Thanks to (8), we have for all t ∈ [0, T ],

ω(t) L ∞ ω(t) s-1 s L 2 ω(t) 1 s
Ḣs .

Then, using (75) and (74), from (76), we deduce that there exists two constant C 1 > 0 and C 2 > 0 depending only on r, s such that for all t ∈ [0, T ],

ω(t) L ∞ ≤ C 1 ω 0 s-1 s L 2 ω 0 1 s Ḣs e C 2 ε α0 . (77)
We assume in what follows that,

C 1 ω 0 s-1 s L 2 ω 0 1 s Ḣs e C 2 ε α0 ≤ 1 2 . ( 78
)
Then thanks to (77) and ( 78), we obtain a contradiction with (71), therefore T = T * . Owing to (78), this means for all t ∈ [0, T * [,

ω(t) L ∞ ≤ 1 2 , (79)
which leads to a contradiction with (70) and then T * = ∞. Therefore, under the assumption (78), we deduce that there exists a unique global strong solution ω

≡ (ρ, v) ∈ C([0, ∞[; Ḣ-r (Ω) 3 ∩ H s (Ω) 3 ) ∩ C 1 (]0, ∞[; H s-1 (Ω) 3 ) ∩ C(]0, ∞[; H s+1 (Ω)
3 ) of the system of Equations ( 30) for the initial data

ω 0 ≡ (ρ 0 , v 0 ) ∈ Ḣ-r (Ω) 3 ∩ H s (Ω) 3 , moreover ω L ∞ (Ω×[0,+∞[) ≤ 1 2 . It remains to prove that ω ∈ C ∞ (Ω×]0, ∞[) 3 . Let > 0, k ∈ N and k = (1 -2 -k ), notice k+1 > k .
By considering for each k ∈ N, the system of equations [START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF] for the initial data ω( k ), then using a recurrence argument combined with Proposition 4.1 and uniqueness of the solution ω, we infer

ω ∈ C([ k , ∞[; Ḣ-r (Ω) 3 ∩ H s+k (Ω) 3 ) ∩ C 1 (] k , ∞[; H s+k-1 (Ω) 3 ) ∩ C(] k , ∞[; H s+k+1 (Ω) 3 ) for any k ∈ N. Since for all k ∈ N, k < , then we get ω ∈ C([ , ∞[; Ḣ-r (Ω) 3 ∩ H s+k (Ω) 3 ) ∩ C 1 ([ , ∞[; H s+k-1 (Ω) 3 ) ∩ C([ , ∞[; H s+k+1 (Ω) 3
) for all k ∈ N. Then, thanks to Sobolev embedding and the system of equations [START_REF] Klusdahl | Self-consistent study of the resonanttunneling diode[END_REF], we deduce that ω ∈ C ∞ (Ω×]0, ∞[) 3 , which allows us to conclude the proof.

Let us compute the real α 1 given in Theorem 4.1 with ω 0 ≡ (ρ 0 , v 0 ) given by ( 28) and [START_REF] Kato | On the nonstationary Navier-Stokes system[END_REF].

Let ñ0 = n 0 -n 0 n 0 and m0 = m 0 n 0 . We recall that m 0 = n 0 u 0 + ε∇n 0 . Then, we get, m0 = (1 + ñ0 )u 0 + ε∇ñ 0 . We recall that Ω def = √ gn 0 O. Then, we notice, for all x ∈ Ω,

ρ 0 (x) = ñ0 x √ gn 0 and v 0 (x) = 1 √ gn 0 m0 x √ gn 0 . By introducing, ω0 def = ñ0 gn 0 , m0 , (80)
we deduce that for all x ∈ Ω,

ω 0 (x) = 1 √ gn 0 ω0 x √ gn 0 .
Then, we get α 1 = α1 √ gn 0 , with, Then, an immediate consequence of Theorem 4.1 with initial data ω 0 ∈ Ḣ-r (Ω) 3 ∩ H s (Ω) 3 with 0 < r < 1 < s is the result which follows. (O) 2 with 0 < r < 1 < s and we assume also that there exists a constant C > 0,

α1 def = C 1 ω0 s-1 s L 2
n 0 -n 0 H s+1 + n 0 -n 0 W 1, 2 1+r ≤ C min(n 0 , √ n 0 ). ( 82
)
We recall that α1 is given by,

α1 = C 1 ω0 s-1 s L 2
where C 1 > 0 and c 1 > 0 are constant.

Since n 0 u 0 n 0 = u 0 + n 0 -n 0 n 0 u 0 and s > 1, we get,

n 0 u 0 n 0 H s u 0 H s 1 + n 0 -n 0 n 0 H s u 0 H s .
Then, we have, .

ω0 H s ≤ n 0 -n 0 √ n 0 H s √ g + n 0 u 0 n 0 H s + ε n 0 ∇n 0 H s C( √ g + ε) + u 0 H s . ( 84
Then, we get, ω0 Ḣ-r ω0 L 2 1+r

≤ n 0 -n 0 √ n 0 L 2 1+r √ g + n 0 u 0 n 0 L 2 1+r + ε n 0 ∇n 0 L 2 1+r C( √ g + ε) + (1 + C) u 0 L 2 1+r , (85)
where we have used the fact that

n 0 n 0 L ∞ ≤ 1 + n 0 -n 0 n 0 L ∞ 1 + n 0 -n 0 n 0 H s ≤ 1 + C.
Under the assumption (82), thanks to (84) and (85), from (83) we deduce that there exists a real C 2 > 0 depending only on r, s, g, ε, u 0 H s , u 0 L

  Since ρ L ∞ (Ω×[0,T ]) ≤3 4 and v L ∞ (Ω×[0,T ]) 2 ≤ 3 4 , then we deduce for a.e τ ∈ [0, T ],|∇v(τ )| ≤ 7 4 ∇ v(τ ) 1 + ρ(τ )+ 3|∇ρ(τ )|. Therefore, we deduce that b(τ ) Theorem, we have ∇ρ(τ ) Ḣσ = ρ(τ ) Ḣσ+1 , ∇v(τ ) Ḣσ = v(τ ) Ḣσ+1 , then from (69), we conclude the proof.

4C ε 2 ω0 2 L 2 .

 22 Ḣr ω0 Ḣ-r e

Corollary 4 . 1 2 √ 2 √ 3 . 2 √

 412232 Let (n 0 -n 0 , m 0 ) ∈ Ḣ-r (O)3 ∩H s (O)3 with 0 < r < 1 < s. If ω0 L ∞ ≤ 1 gn 0 (ω 0 given by (80)) and if α1 ≤ 1 gn 0 (α 1 given by (81)), then there exists a unique global strong solution (n, m) of the system of Equations (23) for the initial data(n 0 , m 0 ) such that (n -n 0 , m) ∈ C([0, ∞[; Ḣ-r (O) 3 ∩ H s (O) 3 ) ∩ C 1 (]0, ∞[; H s-1 (O) 3 ) ∩ C(]0, ∞[; H s+1 (O) 3 ) and n -n 0 L ∞ (O×[0,+∞[) ≤ n 0 2 . Furthermore, (n -n 0 , m) ∈ C ∞ (O×]0, ∞[)Let us give a sufficient condition ensuring that α1 ≤ 1 gn 0 .For this, we assume n0 -n 0 ∈ H s+1 (O) ∩ W 1, 2 1+r (O), u 0 ∈ H s (O) 2 ∩ L 2 1+r

  n 0 u 0 + ε∇n 0 n 0 ,

) 2 1ϕ

 2 For any f ∈ Ḣ-r (O), we can write f Ḣ-r (O) as follows f Ḣ-r (O) = supϕ∈C ∞ 0 (O)∩ Ḣr (O), ϕ Ḣr (O) =1 | f, ϕ |.Since for all ϕ ∈ L -r (O) ∩ Ḣr (O), we have ϕ Ḣr (see[START_REF] Adams | Sobolev spaces[END_REF]), then we deduce that for allf ∈ L 2 1+r (O), f ∈ Ḣ-r (O) and f Ḣ-r (O) f

2 1+r such that α1 ≤ C 2 . 2 .

 222 We get also ω0L ∞ ω0 H s ≤ C( √ g + ε) + u 0 H s .By using Corollary 4.1 and the first equation of[START_REF] Jiang | A remark on weak solutions to the barotropic compressible quantum Navier-Stokes equations[END_REF] to obtain again more regularity on h, we deduce our main Theorem,Theorem 4.2 Let n 0 -n 0 ∈ W 1, 2 1+r (O) ∩ H s+1 (O), u 0 ∈ L 2 1+r (O) 2 ∩ H s (O) 2 with 0 < r < 1 < s. If n 0 -n 0 H s+1 + n 0 -n 0 W 1, 21+r ≤ C 0 min(n 0 , √ n 0 ) with C 0 > 0 a constant, then there exists a real C > 0 depending only on C 0 , r, s, g, ε, u 0 H s , u 0 L 2 1+r such that if n 0 ≥ C then there exists a unique global strong solution (n, u) of the system of Equations (1) for the initial data (n 0 , u 0 ) such that(n -n 0 , u) ∈ C([0, ∞[; Ḣ-r (O) 3 ∩ H s (O) 3 ) ∩ C 1 (]0, ∞[; H s-1 (O) 3 ) and n -n 0 L ∞ (O×[0,+∞[) ≤ n 0Moreover, (n -n 0 , u) ∈ C ∞ (O×]0, ∞[)3 .


Silvère Gauthier
email: silvere.gauthier@lirmm.fr

William Puech

Roseline Bénière

Gérard Subsol

R B Éni Ère

Digitized 3D mesh segmentation based on curvature analysis

Today, it is increasingly frequent and easy to digitize the surface of real 3D objects. However, the obtained meshes are often inaccurate and noisy. In this paper, we present a method to segment a digitized 3D surface from a real object by analyzing its curvature. Experimental results -applied on real digitized 3D meshes -show the efficiency of our proposed analysis, particularly in a reverse engineering process.

Introduction

The availability of 3D scanners has increased the fast development of applications in Computed-Aided Design (CAD), reverse engineering, medical research and inspection. Many 3D processes need to segment an object as a preprocess. In the production line of manufactured objects, steps can be distributed to many partners, and some data can be missing at the end of the line. In reverse engineering, we try to reconstruct the continuous geometry of an object from a 3D mesh, which is discrete. This can lead to quality control or object modification issues for example. In fact, the mesh can be a discretization of a CAD object or a digitized one. To reconstruct the initial geometry, we must take into account the shape of the objects as they are related. But an object shape can be very complex, and the measured data is often noisy. So, we need robust 3D descriptors to accurately define the objects shape.

In previous works, geometry descriptors like curvatures [START_REF] Chen | Intrinsic surface properties from surface triangulation[END_REF][START_REF] Dong | Curvatures estimation on triangular mesh[END_REF] allow us to deal with the 3D mesh shape. But curvatures are locally computed, while it is often necessary to characterize the shape globally. To do this, we can construct curvature distributions [START_REF] Chen | Automatic prismatic feature segmentation of scanning-derived meshes utilising mean curvature histograms[END_REF][START_REF] Masmoudi | A finite mixture model of geometric distributions for lossless image compression[END_REF] and analyze them.

In this paper, we propose a method based on the analysis of a digitized 3D mesh curvature histogram. We use the curvature approximation from Bénière et al. [START_REF] Bénière | A comprehensive process of reverse engineering from 3d meshes to cad models[END_REF] who incorporates two other methods [START_REF] Chen | Intrinsic surface properties from surface triangulation[END_REF][START_REF] Dong | Curvatures estimation on triangular mesh[END_REF]. Then, a distribution is constructed continuously by a kernel estimation from all of the curvature values. Finally, an accurate curvature distribution analysis is realized. In the distribution, we propose to search for peaks and valleys, and automatically compute segmentation thresholds. Indeed, curvature distribution approximately describes the object shape, so these thresholds allow us to extract points belonging to object edges. In the proposed approach, we use the curvature histogram to segment digitized 3D meshes. This paper is organized as follow. We first present previous works in this topic. Then, we expose in detail our proposed method of 3D mesh segmentation. After, we apply it on digitized 3D surfaces of real objects and show that our curvature analysis hugely improves the obtained results, particularly in reverse engineering.

Finally, we conclude and propose directions for future research.

Previous works

In this section, we first present previous research in 3D mesh segmentation. Next, we introduce curvature and previous works on 3D curvatures. Then, we develop previous research in distributions.

Segmentation

A segmentation is the partitioning of a digital image, a 3D mesh or a 3D point cloud in several regions, as illustrated in Fig. 1. Many different 3D mesh segmentation algorithms have been proposed in the literature [START_REF] Petitjean | A survey of methods for recovering quadrics in triangle meshes[END_REF][START_REF] Shamir | A survey on mesh segmentation techniques[END_REF][START_REF] Theologou | A comprehensive overview of methodologies and performance evaluation frameworks in 3d mesh segmentation[END_REF], but each segmentation gives more or less good results depending on the chosen application. Most of the time, a segmentation brings together points with similar criteria. For example, the segmentation can be based on a waterfall, hierarchical clustering, iterative merging or remeshing. In our case, the best results are reached using curvatures, because in reverse engineering, features are extracted from curvature analysis. Some methods use curvatures to segment by discontinuities [START_REF] Benkő | Algorithms for reverse engineering boundary representation models[END_REF][START_REF] Chen | Automatic prismatic feature segmentation of scanning-derived meshes utilising mean curvature histograms[END_REF] or clustering [START_REF] Angelo | Geometric segmentation of 3d scanned surfaces[END_REF], but are often not robust enough around object edges or are designed for CAD meshes [START_REF] Lavoué | A new cad mesh segmentation method, based on curvature tensor analysis[END_REF][START_REF] Bénière | Recovering Primitives in 3D CAD meshes[END_REF]. Indeed, curvatures are often wrong around object edges because adjacent points can run over many different features. Moreover, we are searching for a fully automatic method, so we cannot use parameters like a cluster number or ask for user help.

Curvatures

Intuitively, curvature quantifies the deviation between a curve and a straight line, or between a surface and a plane in 3D. The curvature of a 2D curve at a point P equals the inverse of the osculating circle radius r at P. The osculating circle is the circular arc which best approximates the curve around P (Fig. 2.a). On a 3D surface, an infinity of curvatures exists, each computed on a slice containing the normal vector of point P (Fig. 2.b). So, we need to distinguish particular curvatures. Principal curvatures are minimum and maximum curvatures. If principal curvatures are different, their matching planes are orthogonal. The intersections between these matching planes and the tangent plane at P define the principal directions. Mean curvature and Gaussian curvature equal respectively the mean and the product of principal curvatures. Dong and Wang [START_REF] Dong | Curvatures estimation on triangular mesh[END_REF] compute discrete curvatures with:

k n (t) = < P i -P, N i -N > ||P i -P|| 2 , ( 1
)
where P i is a neighbor of P with a normal N i , and t the projection of P i P on the tangent plane of P. A linear regression is also applied on all discrete curvatures, but with a coefficient fixed to the maximum computed value. Bénière et al. [START_REF] Bénière | A comprehensive process of reverse engineering from 3d meshes to cad models[END_REF] compute discrete curvatures with formula 1 for each neighbor of P, and apply the linear regression of [START_REF] Chen | Intrinsic surface properties from surface triangulation[END_REF]. In our work, we prefer to use this approximation because it is more accurate. Indeed, equation 1 uses each neighbor independently and avoids some curve distortion. Moreover, fixing a linear regression coefficient when we do not know if we have computed the real maximum curvature is dangerous. Curvatures are often used to caracterize surface shape. So, for example, it is possible to analyze a shape to detect saliency or apply segmentation.

Distributions

There are many different types of distributions. For example, we can find probability, frequency or cumulative distributions. But overall, we can distinguish discrete and continous distributions. A discrete distribution is often represented by a histogram, as illustrated in Fig. 3.a. A histogram is a set of intervals with the same width, usually named "bins", which group homogeneous values. Each bin represents an occurence frequency of values from its interval. On the other hand, continous distributions are numerous because they correspond to mathematical models. Common models are Gaussian or Normal distributions (Fig. 3.b), which are couples of a mean value µ and a standard deviation σ .

In computer science, most of the value sets are discrete. So, to generalize methods for continuous spaces, it is necessary to approximate a discrete distribution with a continuous model. Most of the time, models are optimized from the measured data. Model optimization is a minimisation of the difference between the model and the data, updating the model until convergence. A common optimization algorithm is the Expectation-Maximisation (EM) [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. For example, these distributions can be used in image processing to improve compression [START_REF] Masmoudi | A finite mixture model of geometric distributions for lossless image compression[END_REF]. Some methods in 3D mesh processing also begin to use distributions, for example to apply segmentation [START_REF] Chen | Automatic prismatic feature segmentation of scanning-derived meshes utilising mean curvature histograms[END_REF] or extract object edges [START_REF] Demarsin | Meshless extraction of closed feature lines using histogram thresholding[END_REF].

Proposed digitized 3D mesh segmentation

Our proposed method first computes a discrete curvature on each point of the initial 3D mesh. The method constructs a continuous estimated and normalized histogram for each curvature, then analyzes it. Next, a histogram analysis leads to an automatic curvature threshold computation, allowing us to extract object edges. Finally, our proposed method retrieves all disconnected regions using many processes. An overview of our proposed method is illustrated in Fig. 4. In fact, we provide region growing to retrieve disconnected regions, an artefact merging and bridge removal to handling inaccuracies leading to wrong regions, and finally a recursivity handling objects which are a composition of mechanical parts.

Probability curvature distribution

In probability and statistics, a probability distribution assigns a probability to each measurable subset of the possible outcomes of a random experiment, survey, or procedure of statistical inference. We can represent a probability distribution by a histogram. But to define the probability distributions for the simplest cases, we need to distinguish between discrete and continuous random variables. In the discrete case, we can easily assign a probability to each possible value. By contrast, when a random variable like curvature takes values from a continuum then, typically, probabilities can be nonzero only if they refer to intervals.

To approximate continuous curvatures on a discrete 3D mesh, we use the method from Bénière et al. [START_REF] Bénière | A comprehensive process of reverse engineering from 3d meshes to cad models[END_REF]. In our case, meshes can be defined with different scales and can give different curvature ranges. So, to construct curvature histograms which can be compared between many objects, we must normalize curvature values with each mesh. Moreover, curvature values are real, so it is more suitable to estimate a continuous histogram, with kernelestimation for example.

Normalized curvature

Histograms must have the same range to make a comparison. To homogenize histogram range, curvature values have to be normalized. For example, our method normalizes curvature values by multiplying them by the mean edge length of the mesh.

If the object is correctly meshed, we can deduce minimum and maximum possible curvature values, as illustrated Fig. 5. So we propose to limit the histogram range according to extreme possible values, since we normalize curvature by edge length.

Kernel estimation

Curvature values are real, so it is more suitable to compute a histogram with a continuous estimation. Our method computes histograms with a kernel-type estimation. We chose a gaussian kernel because digitized meshes with only one feature often give gaussian-type curvature distributions. We compute the histogram with:

fh (x) = 1 n n ∑ i=1 K h (x -x i ), (2)
where x is the central value of a bin, fh (x) is the quantity inside the bin, n is the number of points, x i is the i th point and h is the kernel standard deviation. The gaussian kernel is defined by:

K h (x -x i ) = 1 h √ 2π e -1 2 ( x-x i h ) 2 . ( 3
)
Each bin of the histogram is computed by centering the kernel on it. Then, kernel density estimation (KDE) is applied on each curvature value and added to the bin. So each bin is computed with a neighborhood defined by the kernel standard deviation. Fig. 6 shows an example of normalized mean curvature distribution with a kernel deviation h = 0.01. This continuous estimation makes histograms less sensible to noise and bin number. Indeed, high frequency fluctuations are naturally smoothed by the kernel.

Distribution analysis

Many caracteristics of a histogram can be useful to many applications. We can for example search for modal number and positions (named "peaks" here), pattern, sparsity and statistics. In our method, we essentialy provide robust peaks and valleys detection, and use those to compute segmentation thresholds.

To detect homogeneous curvature intervals, we need to detect peaks and valleys in the histogram, as illustrated in Fig. 7.

Here, a peak defines a dominant curvature value and a couple of two consecutive valleys defines an homogeneous interval of curvature.

We begin by detecting robust peaks (Fig. 7.b). A peak is a bin with a higher probability than the two adjacent bins. A robust peak is a bin with the higher probability in a window. Thus, we can detect valleys between peaks (Fig. 7.c). For each adjacent peaks couple, we compute the line passing through the two peaks. Then, a valley is a bin with a higher distance to the line than the two adjacent bins. A robust valley is a bin with the higher distance in a window. Our method uses sliding windows of 10 bins for both peaks and valleys detection.

Segmentation

In reverse engineering, a segmentation is often needed. In our case, we work on digitized meshes, with point coordinate inaccuracies and noise. To correctly segment those, we aim to extract object edges matching with intersections between geometric features. We can retrieve disconnected and homogeneous regions matching with these features. Finally, we also apply some postprocessing and recursivity to improve results. Our proposed seg-mentation, based on curvature histogram analysis, is completely automatic.

Edge extraction

Edges of an object are mesh areas caracterized by a high curvature. We propose to extract those edges with curvature histogram analysis. Indeed, high curvatures are on the extremities of curvature histograms. To detect and extract edges, we apply a thresholding on curvature values. Our method defines two thresholds, matching with the extreme left and right valleys of the histogram (Fig. 9). Points with curvature between the two thresholds are labelled "uniform", and others are labelled "edge". We can also realize multiple thresholding by taking each valley couple as thresholds (Fig. 10). This method directly isolates homogeneous intervals, but we need to be aware of potential curvature inaccuracy. To improve edge detection, we may have to remove some valleys, which are due to low sampling on some curvature intervals.

Region growing

After curvature thresholding, we can retrieve connected points by region growing. Our method retrieves triangles instead of points, but it is based on the same principle. A triangle type is defined by its dominant point type. Region growing is a common algorithm that groups similar adjacent elements: take a "seed" triangle and assign a unique ID (Fig. 11.a) and propagate seed ID on its neighbors (Fig. 11.b) until it reaches the borders (Fig. 11.c).

Post-processing

Because of digitization inaccuracy, some regions obtained by region growing can be wrong. To improve segmentation quality, we need to take those into account. Indeed, we often have many regions with an insufficient area, so we chose to merge them with adjacent regions. A small region in our method is defined by an area smaller than:

T hreshold area = 20 × 1 n n ∑ i=1 Area(T i ), (4)
where n is the number of triangles and T i the i th triangle. If a small region has only one adjacent region, it can be merged (Fig. 12). Otherwise, we prefer to remove it to avoid region connectivity modifications. Sometimes, because of digitization inaccuracies, region connectivity can be corrupted by small bridges (Fig. 13.a). We chose to remove these bridges to disconnect the regions into many subregions. To do this, we use an erosion-dilatation algorithm. First, the connected region is eroded many times and disconnected subregions are retrieved by region growing (Fig. 13.b and Fig. 13.c). Then, subregions are dilated as many times as the erosion (Fig. 13.d). Finally, we can merge bridges with only one adjacent region, and remove the others (Fig. 13.e).

Recursivity

A mechanical object can be composed of many parts with different scales. In this case, it is not possible to compute a unique threshold to obtain optimal results. So in our method, we choose to apply recursive segmentation. In fact, we segment the input mesh, then segment again each submesh with the same method, until we obtain only one submesh after segmentation. Since our method is based on curvature histogram analysis, each submesh has a different histogram and we can detect object edges with many scales.

Experimental results

This section presents three meshes from different scanners. We apply segmentation on these meshes and show the results.

To have an accuracy with an order of magnitude of -3 and center values around zero in the same bin, our method uses histograms with 1001 bins. Choosing a good kernel standard deviation is a difficult problem, so we prefer to choose it empiricaly. In fact, our method uses a kernel with a standard deviation h = 0.01. Moreover, we limit our histograms in the interval [-2; 2] since we normalize curvature by edge length. The whole parameter set of our method has been validated experimentally.

Presentation of the used meshes

For experimental results, we used three digitized meshes from two different structured light scanners. The first two come from a first scanner and are illustrated in Fig. 14 and Fig. 15 respectively. The third comes from a second scanner and is illustrated in Fig. 16.

Segmentation result

This section presents results of our segmentation using curvature histograms analysis. Each figure shows edge extraction and final set of submeshes. For a reverse engineering application, our results are very good, since most of the features are correctly disconnected. Indeed, edge extraction is accurate since curvature thresholds are computed from distribution and so are adaptative.

The sharp edges of Aerospace are properly detected (Fig. 17 The sharp edges of Moldy are properly detected (Fig. 18.a) and the features are correctly disconnected. We note that freeforms are not over-segmented. We obtain 48 submeshes whose 100% match with only one feature (Fig. 18.b). We have segmented about 30 meshes, with a processor Intel R Core TM i7-4710 CPU @ 2.50GHz. Some results are presented in Table 1.

We can see that our segmentation is fast: less than one minute, except for very large amounts of triangles. Morover, these times also include curvatures and mesh topology computation, which represent a large part of the time.

To validate our approach, we count the number of submeshes that contain only one feature. We can see in Table 1 that about 96% of submeshes match with only one feature (the remaining 4% can contain similar tangent features).

Since features are correctly disconnected, obtained results are suitable for reverse engineering applications. Indeed, it is more accurate and easy to extract only one feature, than many on the same mesh, since we do not encounter curvature neighborhood problems or feature intersections.

Conclusion

In this paper, we propose a new digitized 3D mesh segmentation based on curvature anaylsis. Our proposed method is fully automatic, which is an advantage for industrial applications like reverse engineering. Our proposed segmentation first constructs a continuous normalized curvature distribution and analyzes it. Then, thresholds are computed from the distribution to extract object edges. In this way, it is possible to retrieve disconnected regions corresponding to object features. We also provide some post-processing like artefact merging, bridge removal and recursivity to improve results.

The most important part of our proposed method is the edge extraction based on curvature distribution analysis. Indeed, final submeshes are more homogeneous, and provide important information for reverse engineering, like neighborhood of the feature.

In future research, we propose to analyze more precisely our curvature distribution construction parameters, like bin number or kernel standard derivation, to improve threshold computing accuracy. We can also improve post-processing in order to take more parameters into account. A further objective consists of extending our segmentation to natural objects.

Figure 1 .

 1 Figure 1. Example of a segmentation: a) A 3D point cloud, b) Segmentation of the point cloud [12].

Figure 2 .

 2 Figure 2. Curvature representation on: a) 2D curves and b) 3D surfaces.

Figure 3 .

 3 Figure 3. a) Discrete and b) Continuous distribution samples.

Figure 4 .

 4 Figure 4. Method overview.

Figure 5 .

 5 Figure 5. Minimum and maximum possible for absolute curvature values.

Figure 6 .

 6 Figure 6. Normalized mean curvature kernel-estimated histogram sample.

Figure 7 .

 7 Figure 7. Histogram peaks and valleys.

Figure 8 .

 8 Figure 8. Valley computing in the histogram extremities.

Figure 9 .

 9 Figure 9. Edge extraction.

Figure 10 .

 10 Figure 10. Homogeneous areas extraction.

Figure 11 .

 11 Figure 11. Region growing algorithm.

Figure 12 .

 12 Figure 12. Small regions (grey) are merged if they only have one adjacent region.

Figure 13 .

 13 Figure 13. Small bridges corrupting connectivy are removed.

Figure 14 .

 14 Figure 14. Initial mesh from Scanner 1: Aerospace.

Figure 15 .

 15 Figure 15. Initial mesh from Scanner 1: Moldy.

Figure 16 .

 16 Figure 16. Initial mesh from Scanner 2: Outlet.

  .a) and the features are correctly disconnected, except for a few tangent ones. We obtain 70 submeshes whose 94.3% match with only one feature (Fig. 17.b).

Figure 17 .

 17 Figure 17. a) Edges extraction and b) Segmentation of Aerospace.

Figure 18 .

 18 Figure 18. a) Edges extraction and b) Segmentation of Moldy.

Figure 19 .

 19 Figure 19. a) Edges extraction and b) Segmentation of Outlet.

Table 1 : Segmentation performances and associated features.

 1

	Mesh	Triangles	Time Regions One feature regions in %
	Vase	20 000	<1s	6	6	100
	Fandisk	23 964	<1s	21	20 95.2
	Lego	24 748	<1s	35	35	100
	Lego small	26 371	<1s	10	10	100
	Cup	55 552	1s	32	32	100
	Yoke	62 276	1s	8	7 87.5
	Manique	65 090	1s	40	33 82.5
	Nespresso	71 012	1s	5	5	100
	MediumBolt	89 000	2s	11	10 90.9
	StripedShoe	100 000	1s	16	16	100
	Connector	195 424	4s	36	33 91.7
	Outlet	195 853	5s	72	72	100
	Etui	210 963	5s	3	3	100
	Shoe	258 994	3s	4	4	100
	Czslowakei	400 026	8s	162	162	100
	Part2	414 823	12s	20	20	100
	Chair	500 000	7s	85	79 92.9
	Gear	500 000	6s	283	279 98.6
	Aerospace	799 296	16s	70	66 94.3
	MasterCylinder	820 793	29s	53	49 92.5
	Moldy	851 194	13s	48	48	100
	Watertight	921 216	16s	40	35 87.5
	OilPump	1 064 031	22s	175	169 96.6
	Carter	1 067 079	34s	108	106 98.1
	Pump	1 105 570	21s	518	502 96.9
	MachinedBlock 1 125 832	33s	113	111 98.2
	Te	1 297 428	40s	39	36 92.3
	Splint	2 095 079 1m09s	21	19 90.5
	Metrologic	2 159 724 1m27s	14	13 92.9
	ProductPart	3 427 245 2m16s	191	182 95.3


Alexandre Meyer

Hector M Briceño

Saïda Bouakaz

User-guided Shape from Shading to Reconstruct Fine Details from a Single Photograph

Many real objects, such as faces, sculptures, or low-reliefs are composed of many detailed parts that can not be easily modeled by an artist nor by 3D scanning. In this paper, we propose a new shape from shading (SfS) approach to rapidly model details of these objects such as wrinkles and reliefs of surfaces from one photograph. The method first determines the surface's flat areas in the photograph. Then, it constructs a graph of relative altitudes between each of these flat areas. We circumvent the ill-posed problem of shape from shading by having the user set if some of these flat areas are a local maximum or a local minimum; additional points can be added by the user (e.g. at discontinuous creases) -this is the only user input. We use an intuitive mass-spring based minimization to determine the final position of these flat areas and a fast-marching method to generate the surface. This process can be iterated until the user is satisfied with the resulting surface. We illustrate our approach on real faces and low-relief photographs.

Introduction

Despite recent advances in surface modeling and deformation, creating photorealistic 3D models remains a difficult and time consuming task. Many real objects, such as people, faces, sculptures, masks or low-reliefs are composed of many detailed parts that can not be easily modeled by an artist. Alternatively, 3D scanning technology is still an expensive process. While much work has been devoted to using several photographs to build 3D models, or to rendering new views from many photographs, little work has been done to address the problem of modeling objects from a single photograph.

The fine aspects of these surfaces appear on photographs as a variation of shading. The methods that recover these features from the shading are called Shape from Shading (SfS). Nevertheless, it has been shown that this is an ill-posed problem [START_REF] Durou | Non-Visible Deformations[END_REF][START_REF] Prados | Shape from shading: a well-posed problem ?[END_REF]: a solution does not necessarily exist and when it exists, it is not unique meaning that different surfaces may have produced a given image. Figure 1 illustrates this point: the image on the left may have been produced by both objects on the right.

In this paper, we propose a new practical Shape from Shading method which can be applied to a real photograph to help on the difficult problem of modeling fine aspect such as wrinkles and reliefs of surfaces. The ill-posed of shape from shading is solved by asking an user to decide whether some areas orthogonal to the viewing direction are a local extrema (local maximum or local minimum) or not. The user information is propagated during a minimization process based on a mass-spring simulation. This massspring minimization has the advantage of presenting a graphical visualization which allows the user to interact during the computation according to his knowledge. Indeed, we have noticed that many minimization approaches like simulated annealing [START_REF] Courteille | A global solution to the sfs problem using b-spline and simulated annealing[END_REF] which are fully automatic do not offer any convenient way of correcting reconstructed surface if it is incorrect. With our approach, the user can visually follow the intuitive massspring minimization and correct any errors in the subsequent reconstruction. Moreover, each time he wants to see the potential reconstructed surface, the fast marching method [START_REF] Sethian | A Fast Marching Level Set Method for Monotonically Advancing Fronts[END_REF][START_REF] Kimmel | Optimal Algorithm for Shape from Shading and Path Planning[END_REF] generates it in few seconds.

Related Work

The problem of surface reconstruction from single image can vary in the degree of user interaction: from fully automatic approaches to interactive modeling. Hoiem et al. in [START_REF] Hoiem | Automatic photo pop-up[END_REF] proposed a fully automatic system that creates a 3D model of the scene made up of several texture-mapped planar billboards. Their approach captures the global geometry of the scene as planar. On a single image without any underlying lines (or planes), for example a face, it seems difficult to obtain better results without using shading information. The Shape from Shading (SfS) problem has been widely studied in the computer vision area. See [START_REF] Kozera | An overview of the shape from shading problem[END_REF][START_REF] Zhang | Shape from Shading: A Survey[END_REF][START_REF] Durou | A Survey of Numerical Methods for Shape from Shading[END_REF] for a survey on SfS methods. The SfS problem is known as a difficult problem because of his ill-posedness [START_REF] Prados | Shape from shading: a well-posed problem ?[END_REF]. Consequently, few approaches have been tried on real photographs. Courteille et al. in [10] propose a method to set flat a photograph of a curved sheet of paper in order to facilitate the character recognition. Prados et al. in [START_REF] Prados | Shape from shading: a well-posed problem ?[END_REF][START_REF] Prados | A unifying and rigorous shape from shading method adapted to realistic data and applications[END_REF] have tried a method based on taking into account the light attenuation on face photographs with relatively good results. Zhu et al. in [START_REF] Zhu | Shape from shading: Recognizing the mountains through a global view[END_REF] tackle the ambiguities of shape from shading by a semi-definite programming relaxation process which flip patches and adjust heights until the result surface has no kinks.

To our knowledge, beside Zeng et al. [START_REF] Zeng | Interactive Shape from Shading[END_REF] SfS approaches are mostly fully automatic in spite of the ill-posed nature of the problem. Zeng's method asks the user to enter some normal in order to determine in which direction the slope is going up to a local maximum. Once all local maxima are computed, they compute the relative altitude between each of them and the fast-marching algorithm generates the surface. Similarly to Zeng et al. , we believe that the ill-posed aspect may only be tackled with user interaction. Comparing to Zeng's approach we differ in the user input aspect. In Zeng's approach, the user has to enter enough normals to capture the small variation of surface such as wrinkles on a forehead. Our approach automatically computes all flat areas and mass-spring minimization is more visual. Thus, the user may interact more directly on the data (the mass-spring graph for us) and may explore different solutions as illustrated on Fig. 1 with the two plausible configurations.

Formulation and Overview of our Approach

Our technique takes as input a color image and produces as output an heightmap which is an image where each pixel stores a distance to the underlying plane. Input images are obtained with a camera using a flash as then only light source.

The first step of our approach is to compute the shading image from the RGB colored photograph. For that, we convert each RGB pixel into YUV color space and keep the luminance Y as shading. A more elaborate solution based on assumption of shading continuity is one proposed by Funt et al. in [START_REF] Funt | Recovering shading from color images[END_REF]: reflectance changes are located and removed from the luminance image by thresholding the gradient at locations of abrupt chromaticity change. More recently, Tappen et al. in [START_REF] Tappen | Integrating Surface Normal Vectors Using Fast Marching Method[END_REF] proposed a solution based on both color information and a classifier trained to recognize gray-scale patterns.

Once the shading image is computed, we have the classical SfS problem. In this part, we assume that the shading image is photographed orthographically, and the scene is composed of Lambertian surfaces which exhibit single-bounce reflections and are illuminated from the camera direction by a point light source at infinity.

Our approach is decomposed as follows:

1. For each pixel of the RGB photograph, extract the shading value by converting it into YUV and taking the Y value (or with [START_REF] Funt | Recovering shading from color images[END_REF]'s or [START_REF] Tappen | Integrating Surface Normal Vectors Using Fast Marching Method[END_REF]'s methods). 2. Detect the flat regions (highlights) in the image which will become the vertices to our relative altitude graph which will guide the user-interaction and the reconstruction (Sect. 4). 3. Using fast-marching (for a recap see Sect. 3.1) we compute the relative altitude difference between the vertices in the relative-altitude graph (Sect. 4). 4. The user defines few vertices as peaks or saddles. 5. The position of the remaining vertices is computed by solving a spring-mass system over the vertices (Sect. 5), and a new surface is quickly computed 6. If features on the surface do not have a vertex associated with them, it is possible for the user to add new vertices to the relative-altitude graph. 7. These three last steps can me re-iterated until the user is satisfied with the reconstructed surface.

Shading Image Formation Model and Fast-Marching

We first review the shading image I(x, y) formation model for a 3D Lambertian object. In our approach the camera and the light-source have the same position which we consider at infinity from the object; we then define the light source direction as L = (0, 0, 1). The surface normal direction is given by N x,y = ( ∂zx,y ∂x , ∂zx,y ∂y , 1). Notice that N is not normalized. The shading image is the dot product of the light and the normalized surface normal, it is computed by:

I x,y = L. N x,y ||N x,y || = 1 ∂zx,y ∂x 2 + ∂zx,y ∂y 2 + 1 With ∇z x,y = ( ∂zx,y ∂x , ∂zx,y ∂y ) we get ||∇z x,y || = I -2 x,y -1
This equation is known as the Eikonal equation which can be solved by the numerical algorithm fast marching which we recap coarsely here. More detailed information can be found in [START_REF] Sethian | A Fast Marching Level Set Method for Monotonically Advancing Fronts[END_REF][START_REF] Kimmel | Optimal Algorithm for Shape from Shading and Path Planning[END_REF]16].

At initialization, all pixel's altitude z x,y are set to ∞ beside few pixels whose altitudes are known. All known pixels are put into a priority queue ordered by their altitude: smaller altitude first. The algorithm extracts pixels from the priority queue until it is empty. Starting from a known pixel, the altitude of its four-connected neighbors is updated and added to the queue. The altitude z x,y of pixel (x, y) is updated as follows:

• Let z 1 = min(z x-1,y , z x+1,y ) and z 2 = min(z x,y-1 , z x,y+1 ).

• If |z 1 -z 2 | < ||∇z x,y || then z x,y = z1+z2+ √ 2×∇z 2 x,y -(z1-z2) 2 2 else z x,y = min(z 1 , z 2 ) + ||∇z x,y ||;

Relative Altitude Graph

Our minimization process described in the next section is based on a relative-altitude graph mapped on the shading image. This section is dedicated to this graph computation.

Vertices Detection

We define a graph over the photograph where the vertices correspond to highlights (high intensity values) on the shading image. These points correspond to singular points on the surface where the gradiant is 0 (e.g. local minima, local maxima, saddle) On the shading image, several adjacent pixels may have the same intensity, we thus consider only one area, and thus one vertex in the graph. Since we assume an orthographic camera, a vertex represent a flat area of the surface orthogonal to the viewing direction. We named them OVD areas for Orthogonal to the Viewing Direction. All these OVD areas are parallel because of the orthographic assumption. Theoretically, these OVD areas have a maximal shading value. Since on real photograph, light attenuation may appear, we define them by a threshold T shading . To compute these OVD areas we consider all regions of 4-connected pixels with equivalent values.

The OVD areas are computed by a depth-first search algorithm on the shading image interpreted as a graph: two pixels are neighbor if their shading values are equal. Among all these areas, we keep as OVD areas only those with a shading value greater than T shading and where its neighboring area is less bright. On the vase image of Fig. 2, this algorithm finds four OVD areas which intuitively are the four highlight spots.

Edges and Relative Altitude Computation

We now determine edges and their weight which are the relative-altitude between vertices. For each vertex v i of the graph, we compute the relative altitude to all other vertices by fast marching: We set to zero the altitude of the pixel p i under the considered vertex v i . The altitude of all other pixels is set to ∞. We run the fast marching algorithm as described on Sect. 3.1. The altitude of all other pixels will be lower than v i because we only descend from this pixel. We then look at the altitude of pixels that correspond to the other vertices v j , j = i. We use this difference to set the weight of the edge e ij to be the relative altitude difference between vertex v i and v j . Notice the fact that relative altitude difference between two vertices might be wrong if there is a inflection point between them, this situation will be addressed in next section by a simplification of the graph. We iterate this process for each vertex until we get the weights of all the edges between all the vertices. After this process, we do not know if a vertex is, for example, a local maximum, local minimum, saddle; we only know its relative altitude to its neighbors.

Mass-spring Simulation and User Interaction

The relative altitude graph is mapped onto the shading image, then it is simplified and converted into a mass-spring network which will serve as a visual aid and a way for the user to correct the minimization process.

Our initial complete graph is composed of C n 2 = 1 2 n(n -1) edges, n being the number of vertices (OVD areas). However, the relative altitude between two far vertices (in the sense of Euclidean distance) is probably incorrect and should be discarded: the monotonic descent assumption used for the relative altitude calculation does not hold if the path between the two vertices crosses a valley, a saddle or a ridge. Thus, a relative altitude is only meaningful for adjacent vertices, the graph can be reduced to a subset of edges. For the same reasons that Zeng et al. in [START_REF] Zeng | Interactive Shape from Shading[END_REF], we simplify the graph by its minimum spanning tree using Prim's algorithm [START_REF] Prim | Shortest connection networks and some generalizations[END_REF]. The number of edges is thus reduced to n. Since all vertices are directly/indirectly connected to each other by the tree, the user can seed a minimization process to find their absolute altitude. Notice, that our approach of computing a complete graph which is then simplified is simple to setup whereas determining directly which vertices are neighbors would have been error-prone.

The weight of each edge is the (relative-)altitude difference between its two vertices but the sign of this relative altitude is unknown: we do not know which vertex is above the other. In others words, we do not which vertices are local maximum, which are local minimum, and which ones are saddle. Thus, we ask a user to select some vertices and to move them up or to the down according to his knowledge of the target surface, this will serve as the initial condition to the minimization process. For example on a face, user will move up the vertex corresponding to the nose. In order to respect the relative altitude constrains between vertices and to propagate user's information to the remaining vertices we build a mass-spring network. This saves the user from having to adjust the altitude of all vertices. Each vertex becomes a mass which will be able to move only in the z direction as its (x,y) position is fixed. Indeed, we consider only relative-altitude between vertices (and not Euclidean distance) this simulation is like a single column of mass linked by springs as illustrates on Fig. 2. All vertices have a mass of 1 meaning any vertex is more important than another. Each edge of the spanning tree becomes a spring having a rest-length equal to the relative altitude between its two vertices.

This mass-spring network is animated by an Euler-explicit integration [START_REF] Desbrun | Interactive animation of structured deformable objects[END_REF] until it stabilizes. This simulation has the two advantages: propagating user inputs and being visually intuitive for the user. Indeed, even before the stabilization, the user may want to interact by moving vertices according to his knowledge of the surface. Moreover, each time he wants to see the potential reconstructed surface, the fast marching method generates it in less than a second. This iterative process "mass-spring simulation/user interaction on the vertices" is running until the reconstructed surface satisfies the user.

Our interactive method allows to deal with surfaces with sharp edges: meaning local minimum or maximum without highlight. Sharp edges are points where the surface is C 0 continuous but piecewise C 1 continuous meaning where the gradient is not smooth. For instance, at the line of the junction of the lips, the surface changes its orientation without producing a flat area with a highlight. Thus, our method allows the user to add a vertex to the relative-altitude graph which will allows a change in the surface orientation. Fig. 3 illustrates this feature: on the left the graph without user intervention produce incoherent lips whereas on the right, after the addition of the blue vertex, our method produces correct lips. Fig. 3. A surface may have sharp edges corresponding to local minimum or maximum without producing highlights. For example, at the line of the junction of the lips, the surface changes its orientation without forming a flat area. To correct this, the user can add vertices to the graph which will allow a change in the orientation. On the left, the graph without user intervention produces incoherent lips whereas on the right, after adding a vertex (blue), our method produces correct lips.

Results

On Fig. 5 we show results from real photographs to demonstrate that our technique is well suited to image-based modeling. This surface is reconstructed from one input photograph in few minutes by an user: graph generation takes around 30 seconds on a Intel Centrino Laptop for approximately 70 vertices, surface generation by fast-marching takes less than a second for an image of 300 × 200. For the faces, the user has to interact with fewer vertices, between 2 and 10. We illustrate our concept on face photographs to show the capability of our technique to capture fine wrinkles of the skin which are difficult to obtain with multi-view approaches. Once the surface is reconstructed as an heightmap, we use the surface normal to extract a pseudo intrinsic color of each pixel by solving the diffuse equation (R, G, B) image = (R, G, B) intrinsic × N.L with L = (0, 0, 1). Thus, a textured image of the surface is obtained by combining the intrinsic color and the shading computed with surface normal. Notice, that if the light is similar to the original photograph (position and color) we should obtain similar results.

The heightmap produced by our system is easily triangulated to a mesh. Nevertheless, it can also be directly included as displacement map, for instance to produce realistic scene of low-relief walls as illustrate on Fig. 4.

Limitations.

Our method allows the reconstruction of fine details, nevertheless there are some limitations with our system. First, the reconstructed surface might have some kinks at the surface jonctions during the fast marching process. On Fig. 5, kinks near the eyebrows are present due to the difference of albedo between the skin and the eyebrows. At the end of [START_REF] Zeng | Interactive Shape from Shading[END_REF], Zeng et al. propose a method to fix this kind of surface incoherency by a minimization process. Second, if there are too many fine details, the amount of user input can become important. It is concevaible to add heuristics, to modify the mass-spring simulation or to use hierarchical approach to alleviate this limitation. Additionnally, the algorithm supposes that the surface is C1 continuous, thus discontinuities in the surface can be hard to capture. This problem is partially mitigated by allowing the user to add vertices to the relative altitude graph. Fig. 4. On the left, the original photograph and the computed shading image. On the right, the reconstructed surface representing a low-relief hand. After an empirical test on this image, we do not perform any particular process to manage the specular aspect, nor to take into account that the wall behind the hand has probably a different albedo than the hand.

Conclusion and future work

Starting from a single color image (a photograph), we have presented an intuitive method for user reconstruction of surfaces which may have produced the image. Our method is interactive, guided by the user it may reconstruct different surfaces for a same input image. It allows to explore different SfS solutions in case of doubt. For example, the image on the left of Fig. 1 which may have been produced by the two surfaces on the right. This exploration facility allows the user to interactively, quickly, and easily reconstruct the surface of a given object with only one photograph. The ambiguity around the global shape of the photographed object is hard to resolve automatically without any a priori knowledge, so we ask the user to specify few local extrema (maximum or minimum). Since reconstructed surface is computed in less than few seconds, it is easy for the user to converge to a surface. Manual intervention is only needed to reconstruct the global shape whereas fine part of the surface is automatically extracted.

Finally, we believe that a little user interaction can help to reconstruct many real objects. Thus, SfS approaches may be practically included in 3D mesh modelers1 by defining a shape by example paradigm. In the future, we also would like to combine SfS approaches to global surface reconstruction based on multiple-views. Using shading information, our technique allows to capture them from a single photograph. We show the original image, the computed shading image used to reconstruct the surface, a rendering of the reconstructed image with only the shading computed using the normal and some results with the color texture. Top left is a photo downloaded from the web. Others are extracted from a 640 × 480 video sequence. Notice that Zeng et al. in [START_REF] Zeng | Interactive Shape from Shading[END_REF] propose a minimization method to fix the kinks of the surface due to fast marching imprecision like the ones present near the eyebrows (bottom).

Fig. 1 .

 1 Fig. 1. Shape from Shading is an ill-posed problem: two different surfaces may produce the same shaded image. The left image may have been produced by several surfaces shown on the on the right. The highlights correspond to the flat areas. The upper right image was produced considering highlights A and C as peaks, and the lower right image considering only highlight B as a peak. Other combinations are possible.

Fig. 2 .

 2 Fig. 2. (a) Each highlight of the shading image gives a vertex in the relative-distance graph. Fast marching algorithm computes the relative distance between each pair of vertices. The graph is simplified into a minimum spanning tree. (b) An iterative process of mass-spring simulation/user inputs on the graph (c) runs until the user is satisfied by the reconstructed surface (d). Notice that, since our graph represents the relative-altitude (and not euclidean distance), each vertex can move only in a column (change altitude) during the mass-spring simulation.

Fig. 5 .

 5 Fig. 5. Fine detail of facial expression are hard to capture because of the wrinkles of the skin.Using shading information, our technique allows to capture them from a single photograph. We show the original image, the computed shading image used to reconstruct the surface, a rendering of the reconstructed image with only the shading computed using the normal and some results with the color texture. Top left is a photo downloaded from the web. Others are extracted from a 640 × 480 video sequence. Notice that Zeng et al. in[START_REF] Zeng | Interactive Shape from Shading[END_REF] propose a minimization method to fix the kinks of the surface due to fast marching imprecision like the ones present near the eyebrows (bottom).

such as Maya(Alias Wavefront), 3D Studio Max(Discreet) or Image Modeler(Realviz).