Keywords: Nonlinear subdivision scheme, convergence of subdivision schemes, box splines

We study in this paper nonlinear subdivision schemes in a multivariate setting allowing arbitrary dilation matrix. We investigate the convergence of such iterative process to some limit function. Our analysis is based on some conditions on the contractivity of the associated scheme for the differences.

In particular, we show the regularity of the limit function, in L p and Sobolev spaces.

Introduction

Subdivision schemes have been the subject of active research in recent years. In such algorithms, discrete data are recursively generated from coarse to fine by means of local rules. When the local rules are independent of the data, the underlying refinement process is linear. This case is extensively studied in literature. The convergence of this process and the existence of the limit function was studied in [START_REF] Cavaretta | Stationary Subdivision[END_REF] and [START_REF] Dyn | Subdivision Schemes in computer aided geometric design, Advances in Numerical AnalysisII[END_REF] when the scales are dyadic. When the scales are related to a dilation matrix M, the convergence to a limit function in L p was studied in [START_REF] Han | Multivariate Refinement Equation and Convergence of Subdivision Schemes[END_REF] and generalized to Sobolev spaces in [START_REF] Chen | Convergence of Vector Subdivision Schemes in Sobolev Spaces[END_REF] and [START_REF] Jia | Convergence of Cascade Algorithms in Sobolev Spaces and Integrals of Wavelets[END_REF]. In the linear case, the stability is a consequence of the smoothness of the limit function.

The nonlinearity arises naturally when one needs to adapt locally the refinement rules to the data such as in image or geometry processing. Nonlinear subdivision schemes based on dyadic scales were originally introduced by Harten [START_REF] Harten | Discrete Multiresolution Analysis and Generalized Wavelets[END_REF][10] through the so-called essentially non-oscillatory (ENO) methods. These methods have recently been adapted to image processing into essentially non-oscillatory edge adapted (ENO-EA) methods. Different versions of ENO methods exist either based on polynomial interpolation as in [START_REF] Cohen | Matei Quasi-linear Subdivision Schemes with Applications to ENO Interpolation[END_REF] [START_REF] Arandiga | Matei Approximation of Piecewise Smooth Functions and Images by edge-adapted (ENO-EA) multiresolution Techniques[END_REF] or in a wavelet framework [START_REF] Chan | ENO-Wavelet Transforms for Piecewise Smooth Functions[END_REF], corresponding to interpolatory or noninterpolatory subdivision schemes respectively.

In the present paper, we study nonlinear subdivision schemes associated to dilation matrix M. After recalling the definitions on nonlinear subdivision schemes in that context, we give sufficient conditions for convergence in Sobolev and L p spaces.

General Setting

Notations

Before we start, let us introduce some notations that will be used throughout the paper. We denote #Q the cardinal of the set Q. For a multi-index

µ = (µ 1 , µ 2 , • • • , µ d ) ∈ N d and a vector x = (x 1 , x 2 , • • • , x d ) ∈ R d we define |µ| = d i=1 µ i , µ! = d i=1 µ i ! and x µ = d i=1 x i µ i .
For two multi-index m, µ ∈ N d we also define

  µ m   =   µ 1 m 1   • • •   µ d m d   .
Let ℓ(Z d ) be the space of all sequences indexed by Z d . The subspace of bounded sequences is denoted by ℓ ∞ (Z d ) and u ℓ ∞ (Z d ) is the supremum of

{|u k | : k ∈ Z d }.
We denote ℓ 0 (Z d ) the subspace of all sequences with finite support (i.e. the number of non-zero components of a sequence is finite).

As usual, let ℓ p (Z d ) be the Banach space of sequences u on Z d such that u ℓ p (Z d ) < ∞, where

u ℓ p (Z d ) := k∈Z d |u k | p 1 p for 1 ≤ p < ∞.
As in the discrete case, we denote by L p (R d ) the space of all measurable functions f such that f L p (R d ) < ∞, where

f L p (R d ) := R d |f (x)| p dx 1 p for 1 ≤ p < ∞ and f L ∞ (R d ) is the essential supremum of |f | on R d .
Let µ ∈ N d be a multi-index, we define ∇ µ the difference operator

∇ µ 1 1 • • • ∇ µ d d ,
where ∇ µ j j is the µ j th difference operator with respect to the jth coordinate of the canonical basis. We define

D µ as D µ 1 1 • • • D µ d d
, where D j is the differential operator with respect to the jth coordinate of the canonical basis. Similarly, for a vector x ∈ R d the differential operator with respect to x is denoted by D x .

A matrix M is called a dilation matrix if it has integer entries and if lim n→∞ M -n = 0. In the following, the invertible dilation matrix is always denoted by M and m stands for |det(M)|.

For a dilation matrix M and any arbitrary function Φ we put Φ j,k (x) = Φ(M j x -k).

We also recall that a compactly supported function Φ is called L p -stable if there exist two constants C 1 , C 2 > 0 satisfying

C 1 c ℓ p (Z d ) ≤ k∈Z d c k Φ(x -k) L p (R d ) ≤ C 2 c ℓ p (Z d ) .
Finally, for two positive quantities A and B depending on a set of parameters, the relation A < ∼ B implies the existence of a positive constant C, independent of the parameters, such that A ≤ CB. Also A ∼ B means A < ∼ B and B < ∼ A.

Local, Bounded and Data Dependent Subdivision Operators, Uniform Convergence Definition

In the sequel, we will consider the general class of local, bounded and data dependent subdivision operators which are defined as follows:

Definition 1. For v ∈ ℓ ∞ (Z d
), a local, bounded and data dependent subdivision operator is defined by

S(v)w k = l∈Z d a k-M l (v)w l , (1)
for any w in ℓ ∞ (Z d ) and where the real coefficients

a k-M l (v) ∈ R are such that a k-M l (v) = 0, if k -Ml ℓ ∞ (Z d ) > K (2)
for a fixed constant K. The coefficients a k (v) are assumed to be uniformly bounded by a constant C, i.e. there is C > 0 independent of v such that:

|a k (v)| ≤ C.
Note that the definition of the coefficients depends on some sequence v, while S(v) acts on the sequence w.

Note also that, from

v j = Sv j-1 = S(v j-1 )v j-1 , j ≥ 1.
(3)

Polynomial Reproduction for Data Dependent Subdivision Operators

The study of the convergence of data dependent subdivision operators will involve the polynomial reproduction property. We recall the definition of the space P N of polynomials of total degree N:

P N := {P ; P (x) = |µ|≤N a µ x µ }.
With these notations, the polynomial reproduction properties read:

Definition 2. Let N ≥ 0 be a fixed integer.

1. The data dependent subdivision operator S has the property of reproduction of polynomials of total degree N if for all u ∈ ℓ ∞ (Z d ) and P ∈ P N there exists P ∈ P N with P -P ∈ P N -1 such that S(u)p = p where p and p are defined by p k = P (k) and pk = P (M -1 k).

The data dependent subdivision operator S has the property of exact reproduction of polynomials of total degree

N if for all u ∈ ℓ ∞ (Z d )
and P ∈ P N , S(u)p = p where p and p are defined by p k = P (k) and pk = P (M -1 k).

Remark:

The case N = 0 is the so-called "constant reproduction property". For a data dependent subdivision operator defined as in (1), the constant reproduction property reads

k∈Z d a k-M l (v) = 1, for all v ∈ ℓ ∞ (Z d ).

Definition of Schemes for the Differences

Another ingredient for our study is the schemes for the differences associated to the data dependent subdivision operator. The existence of schemes for the differences is obtained by using the polynomial reproduction property of the data dependent subdivision operator.

Let us denote ∆ l = (∇ µ , |µ| = l) and then state the following result on the existence of schemes for the differences:

Proposition 1.
Let S be a data dependent subdivision operator which reproduces polynomials up to total degree N. Then for 1 ≤ l ≤ N + 1 there exists a data dependent subdivision rule S l with the property that for all v,w

in ℓ ∞ (Z d ), ∆ l S(v)w := S l (v)∆ l w Proof:
Let l be an integer such that 1 ≤ l ≤ N + 1. By using the definition of ∇ µ with |µ| = l, we write:

∇ µ S(v)w k = ∇ µ 1 1 • • • ∇ µ d d S(v)w k .
From the definition of S(v)w we infer that

∇ µ S(v)w k = max(µ 1 ,••• ,µ d ) m 1 ,••• ,m d =0 (-1) l   µ m   p∈Z d a k-m•e-M p (v)w p ,
where we have used the notation m

• e = m 1 e 1 + • • • + m d e d . Straightforward computations give ∇ µ S(v)w k = p∈Z d w p max(µ 1 ,••• ,µ d ) m 1 ,••• ,m d =0 (-1) l   µ m   a k-m•e-M p (v) = p∈Z d w p f k,p (v, µ). (4)
Let us clarify the definition of f k,p (v, µ). Since the data dependent subdivision operator is local we have a k-M p (v) = 0 for any data v ∈ ℓ ∞ (Z d ) and any index

k such that k -Mp ℓ ∞ (Z d ) > K. Now by putting k = ε + Mn, we get that f k,p (v, µ) is defined for p in the set V µ (k) := p : n -p + M -1 (ε -m • e) ∞ ≤ K M -1 ∞ , 0 ≤ m i ≤ µ i ∀i
Then, we define V (k) := {p : k -Mp ∞ ≤ K}. Since the data dependent subdivision scheme reproduces polynomials up to total degree N, we have for any |ν| = r ≤ N:

p∈V (k) a k-M p (v)p ν = P ν (k) for all k ∈ Z d , (5)
where P ν is a polynomial of total degree r. By tacking the differences of

order |ν ′ | = r + 1 in (5) we get p∈V ν ′ (k) f k,p (v, ν ′ )p ν = 0.
Note that the above equality is true for any ν such that |ν| = r. We deduce

that (f k,p (v, ν ′ )) k ∈ Z d is orthogonal to (p q ) p∈V ν ′ (k) where |q| ≤ r. Note that (∇ ν δ n-β ) n∈V ν ′ (k) , |ν| = r + 1, β ∈ Z d spans (p q ) p∈V ν ′ (k)
and we may thus write for any p ∈ V ν ′ (k):

f k,p (v, ν ′ ) = |ν|=r+1 r∈Z d c µ k,r (v)∇ µ δ p-r .
Now, by using (4) we obtain for any |µ| ≤ N + 1:

∇ µ S(v)w k = p∈V µ (k) w p |ν|=l r∈Z d c ν k,r (v)∇ ν δ p-r = p∈V µ (k) |ν|=l c ν k,r (v)∇ ν w p
If we now make µ vary, we obtain the desired relation.

Now that we have proved the existence of schemes for the differences, we introduce the notion of joint spectral radius for these schemes, which is a generalization of the one dimensional case which can be found in [START_REF] Harizanov | Stability of Nonlinear Subdivision and Multiscale Transforms[END_REF].

Definition 3. Let S(v) : ℓ p (Z d ) → ℓ p (Z d ) be a data dependent subdivision operator such that the difference operators S l (v) : ℓ p (Z d ) q l → ℓ p (Z d ) q l , with q l = #{µ, |µ| = l} exists for l ≤ N + 1. Then, to each operator S l , l = 0, • • • , N + 1 (putting S 0 = S)
we can associate the joint spectral radius given by

ρ p,l (S) := inf j≥1 (S l ) j 1 j ℓ p (Z d ) q l .
In other words, ρ p (S) is the infimum of all ρ > 0 such that for all v ∈

ℓ p (Z d ), one has ∆ l S j v ℓ p (Z d ) q l < ∼ ρ j ∆ l v ℓ p (Z d ) q l , (6)
for all j ≥ 0.

Remark: Let us define a set of vectors

{x 1 , • • • , x n } such that [x 1 , • • • , x n ]Z n = Z d , n ≥ d (i.e
. a set such that the linear combinations of its elements with coefficients in Z spans Z d ). We use the bold notation in the definition of the set so as to avoid the confusion with the coordinates of vector x. Then, consider the differences in the directions x 1 , • • • , x n . One can show that there exists a scheme for that differences which we call Sl for l ≤ N + 1 provided the data dependent subdivision operator reproduces polynomials up to degree N (the proof is similar to that using the canonical directions). If we denote by ∆l the difference operator of order l in the directions

x 1 , • • • , x n , one can see that ∆l v ℓ p (Z d ) ql ∼ ∆ l v ℓ p (Z d ) q l for all v in ℓ p (Z d ) and where ql = #{µ, |µ| = l, µ = (µ i ) i=1,••• ,n }.Then, following (6)
, one can deduce that the joint spectral radius of Sl is the same as that of S l .

Convergence in L p spaces

In the following, we study the convergence of data dependent subdivision schemes in L p which corresponds to the following definition:

Definition 4. The subdivision scheme v j = Sv j-1 converges in L p (R d ), if
for every set of initial control points v 0 ∈ ℓ p (Z d ), there exists a non-trivial function v in L p (R d ), called the limit function, such that

lim j→∞ v j -v L p (R d ) = 0.
where

v j (x) = k∈Z d v j k φ j,k (x) with φ(x) = d i=1 max(0, 1 -|x i |).

Convergence in the Linear Case

When S is independent of v, the rule (1) defines a linear subdivision scheme:

Sv k = l∈Z d a k-M l v l .
If the linear subdivision scheme converges for any v ∈ ℓ p (Z d ) to some function in L p (R d ) and if there exists v 0 such that lim

j→+∞ v j = 0, then {a k , k ∈ Z d } determines a unique continuous compactly supported function Φ satisfying Φ(x) = k∈Z d a k Φ(Mx -k) and k∈Z d Φ(x -k) = 1. Moreover, v(x) = k∈Z d v 0 k Φ(x -k).

Convergence of Nonlinear Subdivision Schemes in L p Spaces

In the sequel, we give a sufficient condition for the convergence of nonlinear subdivision schemes in L p (R d ). This result will be a generalization of the existing result in the linear context established in [START_REF] Han | Multivariate Refinement Equation and Convergence of Subdivision Schemes[END_REF] and only uses the operator S 1 .

Theorem 1. Let S be a data dependent subdivision operator that reproduces the constants. If ρ p,1 (S) < m 1 p , then Sv j converges to a L p limit function.

Proof: Let us consider

v j (x) := k∈Z d v j k φ j,k (x), (7)
where

φ(x) = d i=1 max(0, 1 -|x i |) is the hat function.
With this choice, one can easily check that

k∈Z d φ(x -k) = 1. Let ρ p,1 (S) < ρ < m 1 p , it follows that ∆ 1 v j (ℓ p (Z d )) d < ∼ ρ j ∆ 1 v 0 (ℓ p (Z d )) d ,
since q 1 = d. We now show that the sequence v j is a Cauchy sequence in L p :

v j+1 (x) -v j (x) = k∈Z d v j+1 k φ j+1,k (x) - p∈Z d v j p φ j,p (x) = k∈Z d p∈Z d (v j+1 k -v j p )φ j+1,k (x)φ j,p (x)
where we have used

k∈Z d φ(• -k) = 1
. Now, since the subdivision operator reproduces the constants:

v j+1 (x) -v j (x) = p∈Z d k∈Z d l∈Z d a k-M l (v j )(v j l -v j p )φ j+1,k (x)φ j,p (x).
Note that

l∈Z d a k-M l (v j )(v j l -v j p ) = l∈Z d a k-M l (v j )v j l -v j p = l∈Z d (a k-M l (v j ) -δ p-l )v j l .
Since

l∈Z d a k-M l (v j )-δ p-l = 0, ∇ i δ l-β , l ∈ {F (k) ∪ {p}} , β ∈ Z d , i = 1, • • • , d spans (a k-M l -δ p-l ) l∈{F (k)∪{p}}
. This enables us to write:

v j+1 (x) -v j (x) = p∈Z d k∈Z d l∈V (k) {p} d i=1 d i k,p,l ∇ i v j l φ j+1,k (x)φ j,p (x).
Since | k∈Z d φ j+1,k (x)| = 1 following the same argument as in Theorem 3.2 of [START_REF] Han | Multivariate Refinement Equation and Convergence of Subdivision Schemes[END_REF], we may write:

v j+1 -v j L p (R d ) < ∼ m -j p max 1≤i≤d ∇ i v j ℓ p (Z d ) < ∼ m -j p ∆ 1 v j ℓ p (Z d ) ∼ ( ρ m 1 p ) j ∆ 1 v 0 ℓ p (Z d ) (8)
which proves that v j converges in L p , since ρ < m 1 p . Note that, for p = ∞, we obtain that the limit function is continuous.

Furthermore, the above proof is valid for any function Φ 0 satisfying the property of partition of unity when p = ∞. In general, we could show, following Theorem 3.4 of [START_REF] Han | Multivariate Refinement Equation and Convergence of Subdivision Schemes[END_REF], that the limit function in L p is independent of the choice of a continuous and compactly supported Φ 0 .

Uniform Convergence of the Subdivision Schemes to

C s functions (s < 1)
We are now ready to establish a sufficient condition for the C s smoothness of the limit function with s < 1.

Theorem 2. Let S(v) be a data dependent subdivision operator which reproduces the constants. If the scheme for the differences satisfies ρ p,1 (S) < m -s+ 1 p , for some 0 < s < 1 then Sv j is convergent in L p and the limit function is C s .

Proof: First, the convergence in L p is a consequence of ρ p,1 (S) < m

1 p .
In order to prove that the limit function v be in C s , it suffices to evaluate

|v(x) -v(y)| for x -y ∞ ≤ 1. Let j be such that m -j-1 ≤ x -y ∞ ≤ m -j .
We then write :

|v(x) -v(y)| ≤ |v(x) -v j (x)| + |v(y) -v j (y)| + |v j (x) -v j (y)| ≤ 2 v -v j L ∞ (R d ) + |v j (x) -v j (y)| Note that (8) implies that v -v j L ∞ (R d ) < ∼ ρ j ∆ 1 v 0 ℓ ∞ (Z d )
. Since v j is absolutely continuous, it is almost everywhere differentiable, so putting y =

x + M -j h, with h = (h i ) i=1,••• ,d satisfying h ∞ ≤ 1 we get: |v j (x + M -j h) -v j (x)| ≤ |v j (x + M -j h) -v j (x + M -j (h -h d e d ))| +|v j (x + M -j (h -h d e d )) -v j (x + M -j (h -h d e d -h d-1 e d-1 ))| + • • • + |v j (x + M -j (h 1 e 1 )) -v j (x)|
Then, using a Taylor expansion we remark that, there exists θ ∈] -h d , h d [ such that:

|v j (x + M -j h) -v j (x + M -j (h -h d e d ))| = k∈Z d v j k h d D d φ(M j x + h -k + θ d e d ) If we denote Ψ d (x) = Φ(x 1 ) • • • Φ(x d-1 )Ψ(x d )
, where Ψ is the characteristic function of [0, 1] and Φ(x i ) = max(0, 1 -|x i |), we may write:

|v j (x + M -j h) -v j (x + M -j (h -h d e d ))| ∼ k∈Z d ∇ d v j k h d Ψ d (y)
where

y i = (M -j x + h -k + θ d e d ) i if i < d and y d = 2(M -j x + h + θ d e d ) d -k d
(we have used the fact that the differential of the hat function Φ is the Haar wavelet). Iterating the procedure for other differences in the sum, we get:

|v j (x + M -j h) -v j (x)| < ∼ d i=1 ∇ i v j ℓ ∞ (Z d ) < ∼ ∆ 1 v j ℓ ∞ (Z d ) .
Combining these results we may finally write:

|v(x) -v(y)| ≤ |v(x) -v j (x)| + |v(y) -v j (y)| + |v j (x) -v j (y)| ≤ 2 v -v j L ∞ (R d ) + |v j (x) -v j (y)| < ∼ ρ j ∆ 1 v 0 ℓ ∞ (Z d ) + ∆ 1 v j ℓ ∞ (Z d ) < ∼ ρ j ∆ 1 v 0 ℓ ∞ (Z d ) < ∼ x -y s ∞ with s < -log(ρ ∞,1
)/ log m.

Examples of Bidimensional Subdivision Schemes

In the first part of this section, we construct an interpolatory subdivision scheme having the dilation matrix the hexagonal matrix which is:

M =   2 1 0 -2   ,
For the hexagonal dilation matrix, the coset vectors are ε 0 = (0, 0) T , ε 1 =

(1, 0) T , ε 2 = (1, -1) T , ε 3 = (2, -1) T . The coset vector ε i ,i = 0, • • • , 3 of M defines a partition of Z 2 as follows:

Z 2 = 3 i=0 Mk + ε i , k ∈ Z 2 .
The discrete data at the level j, v j is defined on the grid Γ j = M -j Z 2 , the value v j k is then associated to the location M -j k. We now define our bidimensional interpolatory subdivision scheme based on the data dependent subdivision operator which acts from the coarse grid Γ j-1 to the fine grid grid Γ j . To this end, we will compute v j at the different coset points on the fine grid Γ j using the existing values v j-1 of the coarse grid Γ j-1 , as follows: for the first coset vector ε 0 = (0, 0) T we simply put v j M k+ε 0 = v j-1 k , for the coset vectors ε i , i = 1, • • • , 3. the value v j M k+ε i , i = 1, • • • , 3 is defined by affine interpolation of the values on the coarse grid. To do so, we define four different stencils on Γ j-1 as follows:

V j,1 k = {M -j+1 k, M -j+1 (k + e 1 ), M -j+1 (k + e 2 )}, V j,2 k = {M -j+1 k, M -j+1 (k + e 2 ), M -j+1 (k + e 1 + e 2 )}, W j,1 k = {M -j+1 (k + e 1 ), M -j+1 (k + e 2 ), M -j+1 (k + e 1 + e 2 )}, W 2 k = {M -j+1 k, M -j+1 (k + e 1 ), M -j+1 (k + e 1 + e 2 )}.
We determine to which stencils each point of Γ j belongs to, and we then define the prediction as its barycentric coordinates. Since we use an affine interpolant we have:

v j M k = v j-1 k and v j M k+ε 1 = 1 2 v j-1 k + 1 2 v j-1 k+e 1 . (9)
To compute the rules for the coset point ε 2 , V 1 k or V 2 k can be used leading respectively to:

v j,1 M k+ε 2 = 1 4 v j-1 k+e 1 + 1 2 v j-1 k+e 2 + 1 4 v j-1 k v j,2 M k+ε 2 = 1 2 v j-1 k + 1 4 v j-1 k+e 2 + 1 4 v j-1 k+e 1 +e 2 . ( 10
)
When one considers the rules for the coset point ε 3 , W 1 k or W 2 k can be used leading respectively to:

v j,1 M k+ε 3 = 1 4 v j-1 k+e 2 + 1 4 v j-1 k+e 1 +e 2 + 1 2 v j-1 k+e 1 v j,2 M k+ε 3 = 1 4 v j-1 k + 1 4 v j-1 k+e 1 + 1 2 v j-1 k+e 1 +e 2 . ( 11
)
This nonlinear scheme is converging in L ∞ since we have the following result:

Proposition 2. The prediction defined by ( 9), ( 10), [START_REF] Goodman | Spectral Radius formulas for subdivision operators, Recent Advances in Wavelet Analysis[END_REF] satisfies:

∆ 1 v j M.+ε i ℓ ∞ (Z d ) ≤ 3 4 ∆ 1 v j-1 ℓ ∞ (Z d )
We do not detail the proof here but the result is obtained by computing the differences in the canonical directions at each coset points. If we then use Theorem 2 we can find the regularity of the corresponding limit function is C s with s < -log(3/4) log(4) ≈ 0.207. In the second pat of the section, we build an example of bidimensional subdivision scheme based on the same philosophy but this time using as dilation matrix the quincunx matrix defined by:

M =   -1 1 1 1   ,
whose coset vectors are ε 0 = (0, 0) T and ε 1 = (0, 1) T . Note that a 0,0 = 1 and since the nonlinear subdivision operator reproduces the constants we have i a M i+ε = 1 for all coset vectors ε. To build the subdivision operator, we consider the subdivision rules based on interpolation by of first degree polynomials on the grid Γ j-1 . v j M k+ǫ 1 corresponds to a point inside the cell delimited by M -j+1 {k, k + e 1 , k + e 2 , k + e 1 + e 2 }. There are four potential stencils, leading in this case only to two subdivision rules:

vj,1 M k+ǫ 1 = 1 2 (v j-1 k + v j-1 k+e 1 +e 2 ) (12) vj,2 M k+ǫ 1 = 1 2 (v j-1 k+e 1 + v j-1 k+e 2 ) ( 13
)
Note also, that since the scheme is interpolatory we have the relation:

v j M k = v j-1
k . Let us now prove a contraction property for the above scheme.

Proposition 3. The nonlinear subdivision scheme defined by ( 12) and [START_REF] Jia | Convergence of Cascade Algorithms in Sobolev Spaces and Integrals of Wavelets[END_REF] satisfies the following property:

1. when k = Mk ′ : v j,1 M.+ε 1 -v j M. ℓ ∞ (Z d ) ≤ 1 2 ∆ 1 v j-2 . ℓ ∞ (Z d ) v j,2 M.+ε 1 -v j M. ℓ ∞ (Z d ) ≤ ∆ 1 v j-2 . ℓ ∞ (Z d )
2. when k = Mk ′ + ε 1 , we can show that:

v j,2 M.+ε 1 -v j M. ℓ ∞ (Z d ) ≤ 1 2 ∆ 1 v j-2 . ℓ ∞ (Z d ) v j,1 M.+ε 1 -v j M. ℓ ∞ (Z d ) ≤ ∆ 1 v j-2 . ℓ ∞ (Z d )
The proof of this theorem is obtained computing all the potential differences. This theorem shows that the nonlinear subdivision scheme converges in L ∞ since ρ 1,∞ (S) < 1.

Convergence in Sobolev Spaces

In this section, we extend the result established in [START_REF] Jia | Convergence of Cascade Algorithms in Sobolev Spaces and Integrals of Wavelets[END_REF] on the convergence of linear subdivision scheme to our nonlinear setting. We will first recall the notion of convergence in Sobolev spaces in the linear case. Following [START_REF] Jia | Characterisation of Smoothness of Multivariate Refinable Functions in Sobolev Spaces[END_REF] Theorem 4.2, when Φ 0 (x) =

k∈Z d a k Φ 0 (Mx -k) is L p -stable, the so-called
"moment condition of order k + 1 for a" is equivalent to the polynomial reproduction property of polynomial of total degree k for the subdivision scheme associated to a. In what follows, we will say that Φ 0 reproduces polynomial of total degree k. When the subdivision associated to a exactly reproduces polynomials, we will say that Φ 0 exactly reproduces polynomials. We then have the following definition for the convergence of subdivision schemes in Sobolev spaces in the linear case [START_REF] Jia | Convergence of Cascade Algorithms in Sobolev Spaces and Integrals of Wavelets[END_REF]:

Definition 5. We say that v j = Sv j-1 converges in the Sobolev space

W k N (R d ) if there exists a function v in W k N (R d ) satisfying: lim j→+∞ v j -v W k N (R d ) = 0 where v is in W k N (R d ), and v j = k∈Z d v j k Φ 0 (M j x -k)
for any Φ 0 reproducing polynomials of total degree k.

We are going to see that in the nonlinear case, to ensure the convergence we are obliged to make a restriction on the choice of Φ 0 . We will first give some results when the matrix M is an isotropic dilation matrix, we will also emphasize a particular class of isotropic matrices, very useful in image processing.

Definitions and Preliminary Results

Definition 6. We say that a matrix M is isotropic if it is similar to the diagonal matrix diag(σ 1 , . . . , σ d ), i.e. there exists an invertible matrix Λ such that

M = Λ -1 diag(σ 1 , . . . , σ d )Λ, with |σ 1 | = . . . = |σ d | being the eigenvalues of matrix M.

Evidently, for an isotropic matrix holds |σ

1 | = . . . = |σ d | = σ = m 1 d .
Moreover, for any given norm in R d , any integer n and any v ∈ R d we have

σ n u < ∼ M n u < ∼ σ n u .
A particular class of isotropic matrices is when there exists a set ẽ1 , ẽ2 , • • • , ẽq such that:

M ẽi = λ i ẽγ(i) ( 14
)
where γ is a permutation of {1, • • • , q}. Such matrices are particular cases of isotropic matrices since M q = λI where I is the identity matrix and where

λ = d i=1 λ i .
For instance, when d = 2, the quincunx (resp. hexagonal) matrix satisfies M 2 = 2I (resp. M 2 = 4I).

We establish the following property on joint spectral radii that will be useful when dealing with the convergence in Sobolev spaces.

Proposition 4. Assume that S reproduces polynomials up to total degree N.

Then,

ρ p,n+1 (S) ≥ 1 M ∞ ρ p,n (S),
for all n = 0, . . . , N.

Remark: If M is an isotropic matrix and S reproduces polynomials up to total degree N, then

ρ p,n+1 (S) ≥ σ -1 ρ p,n (S),
for all n = 0, . . . , N.

Proof: It is enough to prove

ρ p,1 (S) ≥ 1 M ∞ ρ p (S).
According to the definition of spectral radius there exists ρ > ρ p,1 (S) such that for any u 0

S 1 (u j-1 ) . . . S 1 (u 0 )∇u ℓ p (Z d ) < ∼ ρ j ∇u (ℓ p (Z d ) .
Using the notation ω j := S(u j-1 ) • . . . • S(u 0 )u we obtain

∇ω j ℓ p (Z d ) < ∼ ρ j ∇u ℓ p (Z d ) .
Since

ω j l = n A j l,n u n ,
where

A j l,n = l 1 ,...,l j-1 a l-M l j-1 (u j-1 )a l j-1 -M l j-2 • . . . • a l 1 -M n (u 0 ).
We can write down the ℓ p -norm as follows:

ω j p ℓ p (Z d ) = k∈Z d m j i=1 |ω j M j k+ε j i | p ,
where {ε j i } m j i=1 are the representatives of cosets of M j . First note that:

k -n ∞ ≤ k -n + M -j ε j i ∞ + M -j ε j i ∞ .
Note that M -j ε j i belongs to the unit square so that

M -j ε j i ∞ ≤ K 1 . When A j M j k+ε j i ,n = 0, one can prove that there exists K 2 > 0 such that k -n + M -j ε j i ∞ ≤ K 2 ,
the proof being similar to that of Lemma 2 in [START_REF] Goodman | Spectral Radius formulas for subdivision operators, Recent Advances in Wavelet Analysis[END_REF]. From these inequalities it follows that if A j M j k+ε j i ,n = 0 there exists

K 3 > 0 such that k -n ∞ ≤ K 3 ,
that is, for a fixed k, the values of ω j l for l ∈ {M j k + ε j i } m j i=1 depend only on

u n with n : { k -n ∞ ≤ K 3 }.
Let us now fix k and define ũ such that ũl =

     u l , if k -l ∞ ≤ K 3 ; 0, otherwise.
Let ωj := S(u j-1 ) • . . . • S(u 0 )ũ, then

ωj l =      ω j l , if l ∈ {M j k + ε j i } m j i=1 ; 0, if k -M -j l ∞ ≥ K 4 , since if A j l,n = 0, then k -M -j l ∞ ≤ k -n ∞ + n -M -j l ∞ ≤ K 3 + K 2 := K 4 . Moreover, from k -M -j l ∞ ≤ K 4 , it follows that M j k -l ∞ ≤ K 4 M j ∞ .
Taking all this into account, we get

k∈Z d l∈{M j k+ε j i } |ω j l | p = k∈Z d l∈{M j k+ε j i } |ω j l | p ≤ k∈Z d M j k-l ∞≤K4 M j ∞ |ω j l | p < ∼ M j ∞ ∆ 1 ωj l ℓ p (Z d ) < ∼ ( M ∞ ρ) j ∆ 1 ũ ℓ p (Z d ) .
That is,

ω j ℓ p (Z d ) < ∼ ( M ∞ ρ) j u ℓ p (Z d ) , consequently ρ p (S) < ∼ M ∞ ρ. Now, if ρ → ρ p,1 (
S) we get ρ p (S) ≤ M ∞ ρ p,1 (S).

Convergence in Sobolev Spaces When M is Isotropic

First, Let us recall that the Sobolev norm on W p N (R d ) is defined by:

f W p N (R d ) = f L p (R d ) + |µ|≤N D µ f L p (R d ) . (15)
If one considers a set

x 1 , • • • , x n such that [x 1 , • • • , x n ]Z n = Z d , an equivalent
norm is given by:

f W p N (R d ) = f L p (R d ) + |µ|≤N Dµ f L p (R d ) . ( 16
)
where Dµ = D µ 1 x 1 • • • D µn xn . We then enounce a convergence theorem for general isotropic matrix M: Theorem 3. Let S be a data dependent nonlinear subdivision scheme which exactly reproduces polynomials up to total degree N -1, then the subdivision scheme Sv j converges in W p N (R d ), provided Φ 0 is compactly supported and exactly reproduces polynomials up to total degree N -1 and ρ p,N (S) < m 1 p -s d for some s > N.

(17)

Proof: Note that because of Proposition 4, the hypotheses of Theorem 3 imply that ρ p,k (S) < m

1 d ρ p,k+1 (S) < m 1 p -s-1
d , which means that (17) is also true for k < N. Let us now show that v j is a Cauchy sequence in L p . To do so, let us define

q j (x) = d l=1 λ j,l x l ,
where Λ = (λ j,l ) is defined in [START_REF] Cohen | Matei Quasi-linear Subdivision Schemes with Applications to ENO Interpolation[END_REF]. For a multi-index µ = (µ 1 , . . . , µ d ) ∈ Z d let q µ (x) = q µ 1 1 (x) . . . q µ d d (x).

Since Λ is invertible, the set {q µ : |µ| = N} forms a basis of the space of all polynomials of exact degree N, which proves that

D µ (v j+1 -v j ) L p (R d ) ∼ q µ (D)(v j+1 -v j ) L p (R d )
Now, we use the fact that, since M is isotropic, q µ (D)(f (M j x)) = σ jµ (q µ (D)f )(M j x)

where σµ = d i=1 σ µ i i ( [START_REF] Chen | Convergence of Vector Subdivision Schemes in Sobolev Spaces[END_REF]). We can thus write:

q µ (D)(v j+1 -v j ) = q µ (D) l∈Z d v j+1 l Φ 0 (M j+1 x -l) - l∈Z d v j l Φ 0 (M j x -l) .
that, we need to recall the definition of box splines and some properties that we will use. Let us define a set of n vectors, not necessarily distinct:

X n = {x 1 , • • • , x n } ⊂ Z d \ {0}.
We assume that d vectors of X n are linearly independent. Let us rearrange the family X n such that

X d = {x 1 , • • • , x d } are linearly independent. We denote by [x 1 , • • • , x d ][0, 1[ d the collection of linear combinations d i=1 λ i x i with λ i ∈ [0, 1[
. Then, we define multivariate box splines as follows [START_REF] Chui | Multivariate Splines[END_REF][16]:

β 0 (x, X d ) =    1 | det(x 1 ,••• ,x d )| if x ∈ [x 1 , • • • , x d ][0, 1[ d 0 otherwise β 0 (x, X k ) = 1 0 β 0 (x -tx k , X k-1 )dt, n ≥ k > d. (18)
One can check by induction that the support of β

0 (x, X n ) is [x 1 , x 2 , • • • , x n ][0, 1] n .
The regularity of box splines is then given by the following theorem [START_REF] Prautsch | Bézier and B-spline Techniques[END_REF]:

Proposition 5. β 0 (x, X n ) is r times continuously differentiable if all subsets of X n obtained by deleting r + 1 vectors spans R d .

We recall a property on the directional derivatives of box splines, which we use in the convergence theorem that follows: Proposition 6. Assume that X n \ x r spans R d , and consider the following box spline function s

(x) = k∈Z d c k β 0 (x -k, X n ) then the directional derivative
of s in the direction x r reads:

D xr s(x) = k∈Z d ∇ xr c k β 0 (x -k, X n \ x r ).
We will also need the property of polynomial reproduction which is [START_REF] Prautsch | Bézier and B-spline Techniques[END_REF]: Theorem 4. Let S be a data dependent nonlinear subdivision scheme which reproduces polynomials up to total degree N -1 and assume that M satisfies relation [START_REF] Meignen | Application of the Convergence of the Control Net of Box Splines to Scale-Space Filtering[END_REF], then the subdivision scheme Sv j converges in (20)

Proposition 7. If β 0 (x, X n ) is r times continuously differentiable then, for any polynomial c(x) of total degree d ≤ r + 1, p(x) = i∈Z d c(i)β 0 (x -i, X n ) ( 19
W p N (R d ), if when N ≥ 2, Φ 0 is a C N -2 box spline generated by x 1 , • • • , x n satisfying Φ 0 (x) = k g k Φ 0 (Mx-k) and if N = 1 Φ 0 (x) = k g k Φ 0 (Mx-k) and
Proof: We here prove the case N ≥ 2, the case N = 1 can be proved similarly. First note that since Φ 0 (x) is a C N -2 box spline, we can write for any polynomial p of total degree N -1 at most:

p(M -1 x) = i∈Z d p(i)Φ 0 (M -1 x -i, X n ) = q∈Z d i∈Z d g q-M i p(i)Φ 0 (x -q, X n )
Using Proposition 7 we get p and p have the same leading coefficients, and that i∈Z d g q-M i p(i) is a polynomial evaluated in M -1 i having the same leading coefficients as p. That is to say the subdivision scheme (Sv j ) q = i∈Z d g q-M i v j i reproduces polynomials up to degree N -1.

As already noticed, the joint spectral radius of difference operator is independent of the choice of the directions x 1 , • • • , x n that spans Z d . Furthermore, it is shown in [START_REF] Meignen | Application of the Convergence of the Control Net of Box Splines to Scale-Space Filtering[END_REF], that the existence of a scaling equation for Φ 0 implies that the vectors x i , i = 1, • • • , n satisfy a relation of type [START_REF] Meignen | Application of the Convergence of the Control Net of Box Splines to Scale-Space Filtering[END_REF]. We consider such a set {x i } i=1,••• ,n and then define Φ 0 (x) = β 0 (x, Y N ) the box spline associated to the set

Y N := N x 1 , • • • , x 1 , • • • , N x n , • • • , x n .
which is C N -2 by definition. We then define the differential operator Dµ M -j := Dµ

Dµ M -j-1 (v j+1 (x) -v j (x)) = k∈Z d i∈Z d (a k-M i (v j ) -g k-M i )v j i ( Dµ β 0 )(M j+1 x -k, Y N ) = k∈Z d ∇µ ( i∈Z d (a k-M i (v j ) -g k-M i )v j i )β 0 (M j+1 x -k, Y µ N )
where Y µ N is obtained by removing µ i vector x i , i = 1, • • • , d to Y N and ∇µ = ∇ From this, we conclude that:

Dµ M -j-1 (v j+1 (x) -v j (x)) L p (R d ) < ∼ ρ p,|µ| (S) j m -j+1 p ∆|µ| v j 0 (ℓ p (Z d )) q|µ| .
Now, consider a sufficiently differentiable function f and remark that D M -j-1 x 1 f (x) = (Df )(x).M -j-1 x 1 , where Df is the differential of the function f . We also note that M q = λI which implies that λ = σ q and we then put j + 1 = q × ⌊ j+1 q ⌋ + r with r < q and where ⌊.⌋ denotes the integer part. From this we may write: D M -j-1 x 1 f (x) = σ -q⌊ j+1 q ⌋ (Df )(x).M -r x 1

and then D M -j-1 x 1 f (x) ∼ σ -q⌊ j+1 q ⌋ (Df )(x).x r j where r j depends on j. Making the same reasoning for any order µ of differentiation and any direction x i , we get, in L p :

( Dµ

M -j-1 f )(x) L p (R d ) ∼ σ -q|µ|⌊ j+1 q ⌋ ( Dµ f )(x) L p (R d ) .
We may thus conclude that Dµ (v j+1 (x) -v j (x)) L p (R d ) ∼ Dµ M -j-1 (v j+1 (x) -v j (x)) L p (R d ) σ q|µ|(×⌊ j+1 q ⌋)| < ∼ ρ p,|µ| (S) j m -j+1 p σ q|µ|⌊ j+1 q ⌋ ∆|µ| v 0 (ℓ p (Z d )) q|µ| .

( 1 )

 1 and (2) the new defined value S(v)w k depends only on those values l satisfying k -Ml ℓ ∞ (Z d ) > K. The subdivision operator in this sense is local. To simplify, in what follows a data dependent subdivision operator is an operator in the sense of Definition 1. With this definition, the associated subdivision scheme is the recursive action of the data dependent rule Sv = S(v)v on an initial set of data v 0 , according to:

d

  for some s > N.

1 M

 1 -j x 1 • • • Dµn M -j xn . We will use the characterization[START_REF] Prautsch | Bézier and B-spline Techniques[END_REF] of Sobolev spaces therefore µ = (µ i ) i=1,••• ,n . For any |µ| ≤ N we may write:Dµ M -j-1 (v j+1 (x) -v j (x)) = k∈Z d v j+1 k ( Dµ β 0 )(M j+1 x -k, Y N ) -p∈Z d i∈Z d v j i g p-M i ( Dµ β 0 )(M j+1 x -p, Y N ),using the scaling property satisfied by β 0 . Then, we get:

  µ i x i i=1,••• ,n . As both a k-M. (v j ) and g k-M. reproduce polynomials up to total degree N -1, there exist a finite sequence c k,p such that:∇µ ( i∈Z d (a k-M i (v j ) -g k-M i )v j i ) = p∈V (k) Ṽ (k) |ν|=|µ| c k,p (ν) ∇ν v j p ,where Ṽ (k) = {i, k -Mi ≤ K}, where g k-M i = 0 if k -Mi > K. We finally deduce that:Dµ M -j-1 (v j+1 (x) -v j (x)) = k∈Z d p∈V (k) Ṽ (k) |ν|=|µ| c k,p (ν) ∇ν v j p β 0 (M j+1 x -k, Y µ N ).

We use now the scaling equation of Φ 0 to get q µ (D)(v j+1 -v j ) = q µ (D)

Since S and Φ 0 exactly reproduce polynomials up to total degree N -1, we have for |µ| ≤ N -1:

Remark that g l-M r = 0 for l -Mr > K since Φ 0 is compactly supported.

Since

, we deduce:

Consequently,

Since ρ p,N (S) < m 1/p-s/d , with s > N we obtain

From this we deduce that q µ (D)(v j+1 -v j ) L p (R d ) tends to 0 with j. Making µ vary, we deduce the convergence in

We now show that when the matrix M satisfies ( 14) and when Φ 0 is a box spline satisfying certain properties, the limit function is in

To state the above result, we have used the fact that the joint spectral radius is independent of the directions used for its computation. Since we have the

which tends to zero with j, and thus the limit function is in W p N (R d ) . A comparison between Theorem 3 and 4 shows that when the subdivision scheme reproduces exactly polynomials, which is the case of interpolatory subdivision schemes, the convergence is ensured provided Φ 0 also exactly reproduces polynomials. When the subdivision scheme only reproduces polynomial the convergence is ensured provided that Φ 0 is a box spline. Note also that the condition on the joint spectral radius is the same. We are currently investigating illustrative examples which involve the adaptation of the local averaging subdivision scheme proposed in [START_REF] Cohen | Matei Quasi-linear Subdivision Schemes with Applications to ENO Interpolation[END_REF] to our non-separable context.

Conclusion

We have addressed the issue of the definition of nonlinear subdivision schemes associated to isotropic dilation matrix M. After the definition of the convergence concept of such operators, we have studied the convergence of these subdivision schemes in L p and in Sobolev spaces. Based on the study of the joint spectral radius of these operators, we have exhibited sufficient conditions for the convergence of the proposed subdivision schemes. This study has also brought into light the importance of an appropriate choice of Φ 0 to define the limit function. In that context, box splines functions have shown to be a very interesting tool.


Alessandro Nordio
email: alessandro.nordio@polito.it

Carla-Fabiana Chiasserini
email: carla.chiasserini@polito.it

Emanuele Viterbo
email: viterbo@deis.unical.it

Linear Signal Reconstruction from Jittered Sampling

This paper presents an accurate and simple method to evaluate the performance of AD/DA converters affected by clock jitter, which is based on the analysis of the mean square error (MSE) between the reconstructed signal and the original one. Using an approximation of the linear minimum MSE (LMMSE) filter as reconstruction technique, we derive analytic expressions of the MSE. Through asymptotic analysis, we evaluate the performance of digital signal reconstruction as a function of the clock jitter, number of quantization bits, signal bandwidth and sampling rate.

Introduction

A significant problem in Analog Digital Conversion (ADC) of wide-band signals is clock jitter and its impact on the quality of signal reconstruction. Indeed, even small amounts of jitter can measurably degrade the performance of analog to digital and digital to analog converters. Clock jitter is typically detrimental because the analog to digital process relies upon a sample clock to indicate when a sample or snapshot of the analog signal is taken. The sample clock must be evenly spaced in time; any deviation will result in a distortion of the digitization process. If one had a perfect ADC and a perfect DAC and used the same clock to drive both units, then jitter would not have any impact on the reconstructed signal. In a real world system, however, a digitized signal travels through multiple processors, usually it is stored on a disk or piece of tape for a while, and then goes through more processing before being converted back to analog. Thus, during reconstruction, the clock pulses used to sample the signal are replaced with newer ones with their own subtle variations. Jitter may have different probability distributions which may have different effects on the quality of the reconstructed signal. While several results are available in the literature on jittered sampling [START_REF] Liu | Error bounds for jittered sampling[END_REF][START_REF] Tourabaly | A jittered-sampling correction technique for ADCs[END_REF] as well as on experimental measurements and instruments performance [START_REF] Viterbo | SMT4-CT98[END_REF][START_REF] Arpaia | Characterization of digitizer timebase jitter by means of the Allan variance[END_REF][START_REF] Rubiola | Noise induced time interval measurement biases[END_REF][START_REF] Verspecht | Accurate spectral estimation based on measurements with a distorted-timebase digitizer[END_REF], an analytical methodology for the performance study of the AD/DA conversion is still missing. In this paper we fill this gap and propose a method for evaluating the performance of AD/DC converters affected by This work was supported by Regione Piemonte through the VICSUM project.

jitter, which is based on the analysis of the mean square error (MSE) between the reconstructed signal and the original one [START_REF] Verspecht | Accurate spectral estimation based on measurements with a distorted-timebase digitizer[END_REF]. As reconstruction technique, we consider linear filtering methods, which typically have low complexity and are used in a wide variety of fields. If jitter were known exactly, the linear minimum MSE (LMMSE) reconstruction technique would be optimal, since it minimizes the MSE of the reconstructed signal. In practice this is not the case, hence we apply a reconstruction filter with the same structure of the LMMSE filter, where we let the jitter vanish. Then, we apply asymptotic analysis to derive analytical expressions of the MSE on the quality of the reconstructed signal. We then show that our asymptotic expressions provide an excellent approximation of the MSE even for small values of the system parameters, with the advantage of greatly reducing the computation complexity. We apply our method to study the performance of the AD/DA conversion system as a function of the clock jitter, number of quantization bits, signal bandwidth and sampling rate.

System model

Throughout the paper we use the following notations. Column vectors are denoted by bold lowercase letters and matrices are denoted by bold upper case letters. The (k, q)th entry of the generic matrix Z is denoted by (Z) k,q . The n × n identity matrix is denoted by I n , while I is the generic identity matrix. (•) T is the transpose operator, while (•) † the conjugate transpose operator. We denote by f x (z) the probability density function (pdf) of the generic random variable x, and by E[•] the average operator.

Signal sampling and reconstruction

We consider an analog signal s(t) sampled at constant rate f s = 1/T s over the finite interval [0, M T s ). T s is the sample spacing. When observed over a finite interval, s(t) admits an infinite Fourier series expansion. Let N ′ denote the largest index of the non-negligible Fourier coefficients, then N ′ /T s can be considered as the approximate onesided bandwidth of the signal. We therefore represent the signal by using a truncated Fourier series with

N = 2N ′ + s(t) = 1 √ N N ′ ℓ=-N ′ a ℓ exp j2πℓ t M T s , (1)
0 ≤ t < M T s . The vector a = [a -N ′ , .
. . , a 0 , . . . , a N ′ ] T represents the complex discrete spectrum of the signal.

Observe that the signal representation given in (1) includes sine waves of any fractional frequency f 0 = f s N ′ /M (when a ℓ = 0 for -N ′ < ℓ < N ′ and a -N ′ = a * N ′ ), which are frequently used as reference signal for calibration of ADC [START_REF] Viterbo | SMT4-CT98[END_REF][START_REF]IEEE Standard for Terminology and Test Methods for Analog-to-Digital Converters[END_REF]. We note that when the signal s(t) is observed in the frequency domain through its M samples, the spectral resolution is given by ∆f = 1/(M T s ). Therefore, considering the expression in (1), the signal bandwidth is given by B = N ∆f 2 = N 2M Ts . By defining the parameter

β = M N (2)
as the oversampling factor of the signal s(t) with respect to the Nyquist rate, we can also write:

B = f s /2 β (3)
In work, we consider that sampling locations suffer from jitter, i.e., the m-th sampling location is given by

t m = mT s + d m , (4) m
s = V † a
where V is an N × M random Vandermonde matrix defined as

(V) ℓ,m = 1 √ N exp -j2πℓ t m M T s (5)
ℓ = -N ′ , . . . , N ′ , and m = 0, . . . , N -1. Note that V accounts for the jitter in the AD/DA conversion process, and that the parameter β defined in (2) also represents the aspect ratio of matrix V. Furthermore, in addition to jittered sampling, we assume that signal samples are affected by some additive noise and are therefore given by

y = s + n
where n is a vector of M noise samples, modeled as zero mean i.i.d. random variables. In practice, the dominant additive noise error is due to the n-bit quantization process [START_REF] Ergen | Effects of A-D conversion nonidealities on distributed sampling in dense sensor networks[END_REF].

In order to reconstruct the signal we consider a reconstruction technique that provides an estimate â of the discrete spectrum a. The reconstruction ŝ(t) of s(t) obtained from â is given by

ŝ(t) = 1 √ N N ′ ℓ=-N ′ âℓ exp j2πℓ t M T s

Reconstruction error

We consider as performance metric of the AD/DA conversion process the mean square error (MSE) associated to the estimate. The MSE, evaluated in the observation interval [0, M T s ), can be equivalently computed in both time and frequency domains as:

MSE = E M Ts 0 |s(t) -ŝ(t)| 2 dt = E a -â 2 N
More specifically, we consider the MSE relative to the signal average power, i.e.,

J = MSE σ 2 a
which can be thought of as a noise to signal ratio and will be plotted using a dB scale in our results. Among the possible techniques that can be applied to reconstruct the original signal, we focus on linear filters that provide an estimate of a through the linear operation â = By where B is an N × M matrix.

Jittered AD/DA conversion with linear filtering

Let us assume a 2 = σ 2 a N and E[nn † ] = σ 2 n I, then we define the signal to noise ratio (SNR) in absence of jitter as

γ = σ 2 a σ 2 n
Under the assumption that E[aa † ] = σ 2 a I, the linear filter that provides the best performance in terms of MSE is the linear minimum mean square error (LMMSE) filter, which is given by

B opt = VV † + 1 γ I -1 V (6)
In [START_REF] Nordio | Performance of linear field reconstruction techniques with noise and uncertain sensor locations[END_REF], it has been shown that, by applying the LMMSE filter, we obtain:

J = 1 σ 2 a N E a -â 2 = E tr γVV † + I -1
where tr{•} is the normalized matrix trace operator and the average is over the randomness in V. Note, however, that the filter in ( 6) cannot be employed in practice, since the jitters d m (hence the matrix V) are unknown (see the definition of V in ( 5)). We therefore resort to an approximation of the optimum filter B opt , based on the assumption that jitter has a zero mean.

In particular, we approximate V with the matrix F defined as,

F = V| dm=0
with the generic element of F given by, (F) ℓ,m = exp -j2πℓ m M / √ N , ℓ = -N ′ , . . . , N ′ , and m = 0, . . . , N -1. We observe that F is such that: FF † = βI and it is related to the discrete Fourier transform matrix. Substituting the approximation of V in (6), we obtain:

B = β + 1 γ -1 F (7)
Notice that the filter in [START_REF] Verspecht | Accurate spectral estimation based on measurements with a distorted-timebase digitizer[END_REF] is the LMMSE filter adapted to the linear model y = F † a + n. By letting ω = (β + 1/γ) -1 , the noise to signal ratio J provided by the approximate filter ( 7) is given by

J = 1 σ 2 a N E a -ωFy 2 = tr ω 2 E d FV † VF † -2ωℜ{E d FV † } +1 + ω 2 β γ ( 8
)
where the operator E Assuming the jitters to be independent [START_REF] Viterbo | SMT4-CT98[END_REF] and with characteristic function

C d (w) = E d [exp(jwz)],
the first two terms in ( 8) are given by

tr E d FV † = β N N ′ ℓ=-N ′ C d 2πℓ M T s E d FV † VF † = β + β (β -1) N N ′ ℓ=-N ′ C d 2πℓ M T s 2
Hence, we can write:

J = + ω 2 β 1 + 1 γ -2ω β N N ′ ℓ=-N ′ C d 2πℓ M T s +ω 2 β (β -1) N N ′ ℓ=-N ′ C d 2πℓ M T s 2 (9)
In order to reduce the complexity of the computation of the reconstruction error and provide simple but accurate analytical tools, in the next section we let the parameters N and M go to infinity, while the ratio β = M/N is kept constant. We therefore derive an asymptotic expression of J, which we will show well approximates the expression in ( 9) even for small N and M .

Asymptotic analysis

When N and M grow to infinity while β is kept constant, we define the asymptotic noise to signal ratio J as:

J (β,γ) ∞ = lim N,M →+∞ β J
In [START_REF] Nordio | Performance of linear field reconstruction techniques with noise and uncertain sensor locations[END_REF], it has been shown that J

(β,γ) ∞
provides an excellent approximation of MSE/σ 2 a even for small values of N and M , with the advantage of greatly simplifying the computation. In the limit N, M → ∞ with constant β, we compute

µ 1 = lim N,M →+∞ β 1 N N ′ ℓ=-N ′ C d 2πℓ M T s = 1/2 -1/2 C d (4πBx) dx (10)
where, from (3), we used the fact that 1/βT s = f s /β = 2B. Similarly, we define

µ 2 = lim N,M →+∞ β 1 N N ′ ℓ=-N ′ C d 2πℓ M T s 2 = 1/2 -1/2 |C d (4πBx)| 2 dx (11)
By using (10) (11), and ( 9), the asymptotic expression of J is given by

J (β,γ) ∞ = 1+ω 2 β(1+1/γ)-2ωβµ 1 +ω 2 β(β-1)µ 2 (12)
It is worth mentioning that for large SNRs (i.e., in absence of measurement noise), J

(β,γ) ∞ reduces to J (β) ∞ = lim γ→∞ J (β,γ) ∞ = 1 + 1 β -2µ 1 + 1 - 1 β µ 2 (13)
Equation ( 13) provides us with a floor that represent the best quality of the reconstructed signal (minimum MSE) we can hope for.

Example: uniform jitter distribution

Let us now assume the jitter to be uniformly distributed with pdf given by

f d (z) = 1 2dmax -d max ≤ z ≤ d max 0 elsewhere
where d max is the maximum jitter, independent of the sampling frequency f s . In this case, the characteristic function of the jitter is given by C d (w) = sin(d max w)/(d max w). Then,

µ 1 = Si(2πη u ) 2πη u and µ 2 = cos 2 (2πη u ) + 2πη u Si(4πη u ) -1 4π 2 η 2 u
where Si(•) is the integral sine function and η u = d max B is a dimensionless parameter which relates maximum jitter and signal bandwidth.

Results

For the ease of representation, we assume that the dominant component of the additive noise is due to quantization, and we express the SNR in absence of jitter, γ, as a function of the number of quantization bits n of the ADC [START_REF] Gielen | Analog building blocks for signal processing[END_REF]: (γ) dB = 6.02n + 1.76

Then, in the following plots we show the value of J as a function of γ or, equivalently, of the number of quantization bits n. Figure 1 compares the value of J obtained through its asymptotic expression against the performance of a system with finite parameters values (i.e., the value of J computed using ( 9)). The results are derived for η u = 10 -1 , 10 -2 , 10 -3 , and β = 10. Solid lines refer to the asymptotic expression (12), while markers represent the values of J computed through [START_REF] Gielen | Analog building blocks for signal processing[END_REF], with N ′ = 100. We observe an excellent matching between our approximation of J (β,γ) ∞ and the results computed through [START_REF] Gielen | Analog building blocks for signal processing[END_REF], even for small values of N and M . We point out that this tight match can be observed for any β > 1 and η u ≪ 1.

We also notice that J shows a floor, whose expression is given by (13). This floor is due to the mismatch between the filter F employed in the reconstruction and the matrix V characterizing the sampling system. Furthermore, in the case of unknown jitter, and, thus, of a floor in the behavior of J, there exists a number of quantization bits n = n * beyond which a further increase in the ADC precision does not provide a noticeable decrease in the reconstruction error J. The relation between η u , β, and n * is shown in Figure 2. Note that n * is lightly affected by an increase of β, provided that β > 1, and a good compromise for choosing the oversampling rate is β = 5.

Conclusions

We studied the performance of AD/DA converters, in presence of clock jitter and quantization errors. We considered that a linear filter approximating the LMMSE filter is used for signal reconstruction, and evaluated the system performance in terms of MSE. Through asymptotic analysis, we derived analytical expressions of the MSE which provide an accurate and simple method to evaluate the behavior of AD/DA converters as clock jitter, number of quantization bits, signal bandwidth and sampling rate vary. We showed that our asymptotic approach provides an excellent approximation of the MSE even for small values of the system parameters. Furthermore, we derived the MSE floor, which represents the best reconstruction quality level we can hope for and gives useful insights for the design of AD/DA converters.

  = 0, . . . , M -1, where d m is the associated independent random jitter whose distribution is denoted by f d (z). Typically, we have |d m | ≪ T s . Let the signal samples be s = [s 0 , . . . , s M -1 ] T where s m = s(t m ), 0 ≤ m ≤ M -1. Using (1), the set of signal samples can be written as

d

  [•] averages over the random jitters d m , m = 0, . . . , M -1.

Figure 1 :

 1 Figure 1: Comparison between the reconstruction error J derived through (9), the approximation of J (β,γ) ∞ and the floor J (β) ∞ in (13).

Figure 2 :

 2 Figure 2: Minimum number of bits n * required to reach the floor of J (β,γ) ∞ as a function of β and η u .


Alireza Entezari
email: entezari@cise.ufl.edu

Uniform Sampling and Reconstruction of Trivariate Functions

The Body Centered Cubic (BCC) and Face Centered Cubic (FCC) lattices have been known to outperform the commonly-used Cartesian sampling lattice due to their improved spectral sphere packing properties. However, the Cartesian lattice has been widely used for sampling of trivariate functions with applications in areas such as biomedical imaging, scientific data visualization and computer graphics. The widespread use of Cartesian lattice is partly due to the availability of tensor-product approach that readily extend the univariate reconstruction methods to trivariate setting. In this paper we report on recent advances on non-separable reconstruction algorithms, based on box splines, for reconstruction of data sampled on the BCC and FCC lattices. It turns out that these box spline reconstructions are faster than the corresponding tensorproduct B-spline reconstructions on the Cartesian lattice. This suggests that not only the BCC and FCC lattices are more accurate sampling patterns, their respective reconstruction methods are also more computationally efficient than the tensor-product reconstructions -a fact which is contrary to the common assumption among practitioners.

Introduction

Sampling and reconstruction play a vital role in visualization and computer graphics. Various volume rendering algorithms rely on accurate reconstruction as a key step since the quality and fidelity of the rendered image heavily depends on reconstruction. In image processing reconstruction is used in resampling, resizing, conversion, and manipulation of sampled data. In the realm of sampling, the term regular is often used to refer to the case that the sampling grid is uniform. Although there has been significant research, recently, in non-uniform sampling (e.g., sparse sampling, compressed sensing), the regular sampling is the most commonly-used sampling scheme in practice [START_REF] Unser | Sampling-50 Years after Shannon[END_REF]. When it comes to sampling multivariate functions, the tensor-product of uniform sampling, which forms a Cartesian lattice, is almost always the choice. The simple structure of the Cartesian lattice and its separable nature allows one to readily apply a tensor-product paradigm to many problems in a multi-dimensional setting. The power of the dimensionality reduction will remain the major reason that the Cartesian lattice is the preferred tool in numerical algorithms. The other attraction of the Cartesian lattice is that it simply exists in any dimension and often tools and theory extend to problems in a higher dimensional setting in a trivial manner. However, the Cartesian lattice has been known to be an inefficient lattice from the sampling-theoretic point of view. Miyakawa [START_REF] Miyakawa | Sampling theorem of stationary stochastic variables in multidimensional space[END_REF] and then Petersen and Middleton [START_REF] Petersen | Sampling and Reconstruction of Wave-Number-Limited Functions in N -Dimensional Euclidean Spaces[END_REF] were among the first people to discover the superiority of sphere-packing and sphere-covering lattices for sampling multivariate functions. In particular they have demonstrated that Cartesian lattice is very inefficient for sampling multivariate functions.

Optimal Sampling Lattices

When sampling a multivariate function with a lattice, generated by (integer linear combinations of the columns of) a sampling matrix, M , the spectrum of the signal is contained in the Brillouin zone. Brillouin zone is the Voronoi cell of the reciprocal lattice. The reciprocal lattice to the lattice M is generated by the columns of the matrix 2πM -⊤ . The multivariate version of the Nyquist frequency is the boundary of the Brillouin zone. Without a priori knowledge when sampling multivariate functions, one often assumes that the underlying function has features possibly in all directions. Therefore, without knowledge about particular orientations of high-frequency features, we need to capture an isotropic spectrum during the sampling process. Therefore, the objective of optimal sampling is to maximize the isotropic content of the Brillouin zone. In other words, the sampling lattice whose Brillouin zone has the largest inscribing (hyper) sphere is the best sampling lattice. Therefore, the optimal sampling lattice in any dimension is the lattice whose reciprocal lattice allows for the densest packing of spheres. In the bivariate setting the hexagonal lattice is the best sampling lattice since its reciprocal lattice, which happens to be the dual hexagonal lattice, allows for the best packing of 2-D with disks. When compared to the commonlyused Cartesian lattice with the same sampling density, the hexagonal lattice allows for about 14% more information to be captured in the spectrum of the underlying signal. This is illustrated in Figure 1 as the area of inscribing disc to the Brillouin zone of the hexagonal lattice (i.e., hexagon) is larger than the area of inscribing disc to the Brillouin zone of the Cartesian lattice (i.e., square), even though the two Brillouin zones have the same area.

In the trivariate setting, the optimal sampling lattice is the BCC lattice whose reciprocal lattice (i.e., the FCC lattice) is the densest sphere packing lattice. The sampling efficiency of the BCC lattice, when compared to the commonly-used Cartesian lattice is about 30% higher. Appendix A in [START_REF] Entezari | Optimal Sampling Lattices and Trivariate Box Splines[END_REF] presents a thorough comparison of the Brillouin zone of the Cartesian, BCC and FCC lattices.

The FCC lattice, is also superior to the Cartesian lattice as its efficiency compared to the Cartesian lattice is about 27% higher. Although among the FCC and BCC lattices the BCC wins, by a small margin, for optimal sampling, the FCC lattice appears to have good resistance to aliasing. This can be justified since its reciprocal lattice (i.e., the BCC lattice) allows for the best sphere covering of the space. The best covering of the space translates to replication of isotropic spectrum with minimal overlap between them-minimizing the aliasing for that sampling resolution. These facts about comparison of the Cartesian, BCC and FCC lattices together with their higher-dimensional counter parts are discussed for sampling stationary isotropic random processes [START_REF] Hr Kunsch | Optimal lattices for sampling[END_REF]. The arguments of the optimal sampling (BCC) and resilience to aliasing (FCC) is generalized to the notion that the reciprocal lattice for optimal sphere-packing lattice is the best choice for sampling functions at relatively high resolutions, while the spherepacking lattice is the best option for sampling functions at relatively low resolutions [START_REF] Hr Kunsch | Optimal lattices for sampling[END_REF].

Reconstruction

There is abundant research on reconstruction (i.e., interpolation or approximation) of data based on univariate filtering methods [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]. Various 1-D filters have a low-pass behavior and approximate the ideal kernel (i.e., sinc) for reconstruction into the space of band-limited functions. Bsplines, offer a framework for representation of piecewise polynomial functions and thus are widely used in reconstruction of univariate functions [START_REF] De Boor | A practical guide to splines[END_REF].

There are two common methods for extending the univariate reconstruction 'kernels' to multivariate setting. The separable approach builds the multivariate kernel by a simple tensor-product of univariate kernels. The separable approach is obviously suitable for reconstruction of data on the Cartesian lattice since the lattice itself is also separable. The radial basis approaches construct the multivariate reconstruction kernel by spherical extension of univariate kernel. Due to the spherical extension, the radial basis approach ignores the underlying geometry of the sampling lattice and is often used for scattered data interpolation/approximation. Splines have been widely accepted for image processing [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF]. In the context of image processing, splines are often constructed as a tensor-product of two univariate splines. Mitchell and Netravali [START_REF] Mitchell | Reconstruction Filters in Computer Graphics[END_REF], demonstrated the advantages of using splines for image processing. Recently, Van De Ville [START_REF] Van De Ville | Hex-Splines: A Novel Spline Family for Hexagonal Lattices[END_REF], developed the so called Hexsplines that are used for reconstruction of hexagonal images. Hex-splines can not be constructed as a tensorproduct of univariate splines. Due to the non-separable structure of hexagonal lattice, the tensor-product splines can not be applied for processing of hexagonal data.

Reconstruction of trivariate functions

In the visualization community reconstruction filters have received a lot of attention since accurate reconstruction of trivariate functions and their gradients is crucial in fidelity of rendering algorithms [START_REF] Möller | Design of Accurate and Smooth Filters for Function and Derivative Reconstruction[END_REF][START_REF] Carlbom | Optimal Filter Design for Volume Reconstruction and Visualization[END_REF][START_REF] Dutta Roy | Handbook of Statistics, volume 10, chapter Digital Differentiators[END_REF][START_REF] Möller | A Comparison of Normal Estimation Schemes[END_REF]. Similar to image processing, in volume visualization algorithms, often the tensor-product approach is used for reconstruction of Cartesian sampled data.

Theußl [START_REF] Theußl | Optimal Regular Volume Sampling[END_REF] introduced the BCC sampling in volume rendering. However, since the BCC lattice is a nonseparable lattice, various ad-hoc tensor-product [START_REF] Theußl | Reconstruction schemes for high quality raycasting of the body-centered cubic grid[END_REF] and radial basis [START_REF] Theußl | Optimal Regular Volume Sampling[END_REF] algorithms fail to provide satisfactory reconstruction algorithms and they exhibit blurry artifacts.

Csébfalvi [START_REF] Csébfalvi | Prefiltered gaussian reconstruction for high-quality rendering of volumetric data sampled Cartesian, C 2 , fourth order BCC, C 2[END_REF] proposed a global pre-processing algorithm (based on generalized interpolation [START_REF] Thévenaz | Interpolation revisited[END_REF]) that reconstructs the BCC lattice based on its two Cartesian sub-lattices. This approach is computationally inefficient and does not guarantee approximation order. The author's recent work in this area establishes the relationship between box splines and the above-mentioned sampling lattices. The box splines have been developed as a generalization of B-splines to the multivariate setting. While box splines have been considered as non-separable basis functions for approximation based on their shifts on the Cartesian lattice [START_REF] De Boor | Box Splines[END_REF], here their shifts on BCC and FCC lattices are considered. The interesting fact about these box splines is that while their shifts on the Cartesian lattice do not form a linearly independent set of functions, their shifts on the FCC and BCC lattices are linearly independent -a rare and useful property for the spline space!

Four direction box splines on BCC

The relation of box splines with the BCC lattice was established based on the fact that the immediate neighborhood of a lattice point on the BCC pattern forms a rhombic dodecahedron (see Figure 2). This polyhedron has the special property that is a projection of a four-dimensional hypercube (tesseract). This makes it a perfect match to be the support of a box spline since the geometric definition of box splines precisely amounts to projecting hypercubes (i.e., box) down to lower dimensional spaces. Generally, the class of polytopes that are the shadow of higher dimensional hypercubes are referred to as zonotopes. This linear box spline is defined by the four direction and is a C 0 kernel. The shifts of this box spline on the BCC lattice generate a spline space whose approximation order is two. By convolving this box spline by itself, one obtains a smoother, C 2 , quintic box spline that is specified by a repetition of the four principal directions. The shifts of this box spline generate a spline space whose approximation order is four [START_REF] Entezari | Linear and Cubic Box Splines for the Body Centered Cubic Lattice[END_REF][START_REF] Entezari | Practical box splines for volume rendering on the body centered cubic lattice[END_REF]. This smoothness and approximation order match that of the tricubic B-spline on the Cartesian lattice and hence we compare the two on a Carp fish dataset in first row in Figure 4. The piecewise polynomial representation of these box splines along with efficient evaluation methods can be found in [START_REF] Entezari | Practical box splines for volume rendering on the body centered cubic lattice[END_REF].

The six direction box spline on FCC

Unlike the BCC lattice, the immediate neighborhood in the FCC lattice is not a zonohedron. However, by enlarging the neighborhood one finds the truncated octahedron which is a zonohedron Figure 5. This polyhedron is a projection of a six-dimensional hypercube and the corresponding box spline is a cubic six-direction box spline [START_REF] Entezari | Optimal Sampling Lattices and Trivariate Box Splines[END_REF].

The spline space that is generated by shifts of this cubic box spline on the FCC lattice is a C 1 space whose approximation order is three. These characteristics match the triquadratic B-spline on the Cartesian lattice which is the base for our comparisons in second row in Figure 4.

The piecewise polynomial representation of the cubic box spline along with efficient spline evaluation method on the FCC lattice is demonstrated in [START_REF] Kim | Box Spline Reconstruction on the Face Centered Cubic Lattice[END_REF].

Computational advantages

Once efficient evaluation algorithms are derived for the four-direction box splines [START_REF] Entezari | Practical box splines for volume rendering on the body centered cubic lattice[END_REF] and the six direction box spline [START_REF] Kim | Box Spline Reconstruction on the Face Centered Cubic Lattice[END_REF], one can compare these box spline reconstructions to the commonly-used tensor-product B-spline reconstructions on the Cartesian lattice.

For the C 2 , fourth-order method the tricubic B-spline uses a neighborhood of 4 × 4 × 4 = 64 points for reconstruction, while the quintic box spline only uses a total of 32 points for reconstruction. Therefore as documented in [START_REF] Entezari | Practical box splines for volume rendering on the body centered cubic lattice[END_REF] the BCC non-separable box spline approach outperforms the comparable tensor-product B-spline approach by a factor of two. Similarly the triquadratic B-spline uses a neighborhood of 3 × 3 × 3 = 27 Cartesian data points, while the cubic box spline only requires a total of 16 FCC data points for the reconstruction. Therefore, the non-separable box spline reconstruction outperforms the comparable tensor-product B-spline approach as documented in [START_REF] Kim | Box Spline Reconstruction on the Face Centered Cubic Lattice[END_REF].

Conclusions

The recent research on optimal sampling lattices suggests that not only the FCC and BCC lattices offer higherfidelity sampling schemes, but also their reconstruction algorithms outperform the corresponding tensor-product reconstructions on the traditionally-popular Cartesian lattice. These encouraging results are crucial for acceptance of these efficient lattices in practical applications.

Figure 1 :

 1 Figure1: A square and a hexagon with unit area corresponding to the Brillouin zone of Cartesian and hexagonal sampling. The area of inscribing disk to a square is about 14% less than the area of the inscribing disk to the hexagon.

Figure 2 :

 2 Figure 2: The neighborhood of a BCC lattice point forms a rhombic dodecahedron. This polyhedron is a zonohedron which is the support of a linear box spline.

Figure 3 :

 3 Figure 3: Benchmark example dataset. The CT dataset of a carp fish at a high resolution of 256 × 256 × 256.

Figure 5 :

 5 Figure 5: The neighborhood of a FCC lattice point forms a truncated octahedron. This polyhedron is another zonohedron which is the support of a six-direction box spline.

Acknowledgments

The author would like to thank Dimitri Van De Ville, Torsten Möller and Carl de Boor for valuable insight and advice at various stages of the work.


Marc Louis

Maurice François

Marc François

A new yield criterion for the concrete materials

Keywords: endommagement, critère, béton

The proposed yield criterion depends upon two material constants and is proven to be smooth and convex under a simple condition. These properties induce a mathematical robustness that allows a further use in a damage mechanics model. The analytical gradient and hessian are given. The obtained yield surface is relevant to Kupfer's biaxial testings on concrete. The identification procedure, with respect to the classical uniaxial tension and compression testings, is detailed.

Résumé

Un nouveau critère de limite d'élasticité pour les bétons. Le critère proposé utilise deux constantes matériaux. Son gradient et son hessien sont donnés, sa convexité est démontrée, cette dernière ne dépendant que de la positivité des constantes matériau. Ces propriétés induisent une robustesse mathématique qui autorise son emploi futur au sein d'un modèle d'endommagement. La surface de charge obtenue est en bonne adéquation avec les essais biaxiaux de Kupfer. La procédure d'identification, à partir des contraintes limites d'élasticité identifiées en traction et en compression, est donnée.

keyword damage ; yield function ; concrete

Introduction

Kupfer [START_REF] Kupfer | Behavior of concrete under biaxial stresses[END_REF] performed biaxial testings on concrete (σ I , σ II , σ III = 0 in principal stresses) that are recognized as a reference and relevant for the stress states met in classical civil engineering. Although his best known results are the peak stress curves, the figure 15 of this article exhibits the elastic domain of concrete and the present work is referred to it. Kupfer measured the elastic limit from the loss of linearity of the stress to strain curve: this identification depends upon the precision of the measure. Some authors obtain, in pure tension, a yield stress very close to the peak stress i.e. a fragile behavior; on the contrary, Kupfer identified the yield stress in tension approximately as the half of the peak stress. A similar result is also obtained by Terrien [START_REF] Terrien | Acoustic emission and post-critical mechanical behaviour of a concrete under tensile stress[END_REF] thanks to a very accurate method, confirming Kupfer's results. Most of old elasticity criteria, such as Mohr-Coulomb's one, or more recent approaches [START_REF] Hoek | Empirical strength criterion for rock masses[END_REF], are based on stress vector considerations. They may be accurate, have a strong mechanical sense and allow to define the crack orientation, but are non smooth. These corners are not experimental evidences and lead to numerical difficulties when the criterion is used in a damage model. Actually, the most used criteria are simple functions of the stress tensor σ and the stress deviator σ d classical invariants I 1 = tr(σ), J ′ 2 = tr(σ d .σ d )/2 and J ′ 3 = tr(σ d .σ d .σ d )/3 whose polynomial structure insure the smoothness of the yield surface. Among them, we shall compare to the present model the De Vree [START_REF] De Vree | Comparison of nonlocal approaches in continuum damage mechanics[END_REF], Willam [START_REF] Etse | Fracture energy formulations for inelastic behavior of plain concrete[END_REF] and Maïolino [START_REF] Maïolino | Proposition of a general yield function in geomechanics[END_REF] expressions (all using two constants). The identification will be made preferentially with respect to the uniaxial compression and tension yield stresses, respectively denoted by points C and T on figure [START_REF] Kupfer | Behavior of concrete under biaxial stresses[END_REF]. The De Vree criterion depends upon (I 1 , J ′

2 ); its identification, with respect to the points C and T, gives the constants (γ = 7.75, κ = 3.9 10 -5 ). The Maïolino's criterion has been recently proposed for rock materials; it takes into account the effect of the third invariant J ′ 3 and, in order to keep a twoconstants form, a linear dependancy (σ + = I 1m -I 1 ) has been retained here. The model, in this form, cannot fit both points C and T (the identified constant L s = 0, 5 reaches the bound allowed for the convexity of the criterion) then the identification has been done with respect to the points C and P, giving a good compromise and leading to the constants (I 1m = 2.40 MPa, L s = 0, 62). The Willam's criterion also makes use of the third invariant of the deviator J ′ 3 . The identification with respect to the points C and T leads to the constant e = 0.5 which is the limit for the convexity of the model (dashed curve in figure 1); again, a better identification has been found from points C and P (plain curve in figure 1), with the constants (e = 0.68, m o = 5.3 MPa). It can be seen that these models (and many others, see [START_REF] Ragueneau | Anisotropic damage modelling of biaxial behaviour and rupture of concrete structures[END_REF] for a more complete comparison) fail to describe precisely the particular shape of the concrete's yield surface.

The proposed criterion

In compression, i.e. when (every) σ I 0, the yielding of concrete involves a diffuse microcracking, and high stresses. The mechanism in tension, i.e. when (at least one) σ I 0, is very different, with low stress levels and a quick localization of the microcracks [START_REF] Fokwa | Matériaux hétérogènes : analyse expérimentale et modélisation numérique par une approche hiérarchique[END_REF]. A part of the present criterion is related to the von Mises norm of the stress deviator that appears to be relevant to the compression states, in the low confinement range. The other part uses a stress tensor exponential (recalled in equation 12) in order to describe the dramatic effect of the positiveness of a principal stress, in the sense of the Rankine criterion, via an euclidean norm ( A = (A ij A ij ) 1/2 ). The constant σ y possesses a limiting role; the second constant σ 0 rules the dissymmetry of the yield surface and limits the numerical value of the exponential term. The term √ 3 makes the second member equal to zero when σ = 0.

f (σ) = ||σ d || + σ 0 exp σ σ 0 - √ 3 -σ y (1)
Rewritten in term of principal stresses (σ I , σ II , σ III ), the equation ( 2) exhibits clearly the role of the second member that is positive if at least one principal stress is positive. One can remark that the numerical calculus of a tensor exponential does not require to diagonalize it (it is an integrated function in many softwares).

f (σ I , σ II , σ III ) = (σ II -σ III ) 2 + (σ III -σ I ) 2 + (σ I -σ II ) 2 3 + σ 0 exp 2σ I σ 0 + exp 2σ II σ 0 + exp 2σ III σ 0 - √ 3 -σ y (2)
The constants (σ 0 , σ y ) can be identified with respect to the classical (and simplest) uniaxial compression and tension tests. Calling σ t > 0 and σ c < 0 the experimental elastic limit stresses in pure tension and compression, corresponding respectively to points T and C on figure [START_REF] Kupfer | Behavior of concrete under biaxial stresses[END_REF], the difference f (σ t , 0, 0) -f (σ c , 0, 0) = 0 leads to:

f (σ t , 0, 0) -f (σ c , 0, 0) = 0 σ 0 exp 2σ t σ 0 + 2 + σ t 2 3 = σ 0 exp 2σ c σ 0 + 2 -σ c 2 3 ( 3
)
Each member is a monotonic function of σ 0 : their unique intersection, numerically calculated, gives the value of σ 0 . The equation f (σ c , 0, 0) = 0, or f (σ t , 0, 0) = 0, gives the second constant σ y . Here, the identified values are σ 0 = 0, 455 MPa and σ y = 8, 00 MPa. The figure [START_REF] Kupfer | Behavior of concrete under biaxial stresses[END_REF] shows the good agreement of this criterion with the biaxial testings. The strong curvature around the point C of simple compression is well depicted. The figure [START_REF] Terrien | Acoustic emission and post-critical mechanical behaviour of a concrete under tensile stress[END_REF], in the principal stresses space, shows that this criterion can be seen as a softened version of the combination of von Mises (for the cylindrical shape) and Rankine (for flattened faces) ones. The figure [START_REF] Hoek | Empirical strength criterion for rock masses[END_REF] shows the deviatoric sections of the criterion that evolve from the triangular shape of Rankine envelope to the circular shape of von Mises.

Gradient, hessian and convexity

The use of a yield surface in a damage (or plasticity) model is helped by the knowledge of its gradient. We can already calculate the gradient of f with respect to the stress tensor σ. The derivative of the first member of equation ( 1) is straightforward. For the second one, it is useful to use the property (13). The obtained gradient presents a simple (and intrinsic) expression:

∂f ∂σ = σ d σ d + exp(2σ/σ 0 ) exp(σ/σ 0 ) (4)
The hessian of f may also be useful in numerical calculus. For the first member of equation ( 1) we obtain the following expression, in which ⊗ represents the dyadic product and P d is the fourth rank projector onto the deviatoric subspace:

∂ 2 σ d ∂σ 2 = 1 σ d P d - σ d σ d ⊗ σ d σ d with P d ijkl = 1 2 (δ ik δ jl + δ il δ jk ) - 1 3 δ ij δ kl (5)
The fourth rank operator between parenthesis in the equation ( 5) is a projector onto the four dimensions subspace of the second rank symmetric and deviatoric tensors orthogonal to σ d ; its eigenvalues are then (1, 1, 1, 1, 0, 0) which proves the well known convexity of the von Mises criterion. For the second member of equation ( 1), it is convenient to have a change of variable σ/σ 0 = σ ′ :

σ 0 ∂ 2 ∂σ 2 exp σ σ 0 = 1 σ 0 ∂ 2 exp(σ ′ ) ∂σ ′2 (6)
Using the gradient (4) leads to:

∂ 2 exp(σ ′ ) ∂σ ′2 = 2 exp(σ ′ ) ∂ exp(2σ ′ ) ∂(2σ ′ ) - 1 exp(σ ′ ) 3 exp(2σ ′ ) ⊗ exp(2σ ′ ) (7)
Let us call -→ u i the eigenvectors and λ i the eigenvalues of σ. Then -→ u i are also the eigenvectors of σ ′ and exp(σ ′ ) (these tensors are coaxial with σ) and λ i /σ 0 are the eigenvalues of σ ′ . We consider now the orthonormal base of the second order symmetric tensors B I , with I ∈ (1..6), i.e. B I : B J = δ IJ (the symbol ":" represents the tensor contraction (inner product) defined, for any symmetric second order tensors A and A ′ , as A : A ′ = A ij A ′ ij and δ is the Kronecker delta).

B 1 = -→ u 1 ⊗ -→ u 1 , B 4 = ( -→ u 2 ⊗ -→ u 3 + -→ u 3 ⊗ -→ u 2 ) / √ 2 B 2 = -→ u 2 ⊗ -→ u 2 , B 5 = ( -→ u 3 ⊗ -→ u 1 + -→ u 1 ⊗ -→ u 3 ) / √ 2 B 3 = -→ u 3 ⊗ -→ u 3 , B 6 = ( -→ u 1 ⊗ -→ u 2 + -→ u 2 ⊗ -→ u 1 ) / √ 2 (8)
From formula (14), the second derivative [START_REF] Ragueneau | Anisotropic damage modelling of biaxial behaviour and rupture of concrete structures[END_REF] expresses in the base B I ⊗ B J (in the Voigt notation, the terms 44, 55 et 66 would be divided by 2), in the general case (λ 1 = λ 2 = λ 3 ), as:

∂ 2 exp(σ ′ ) ∂σ ′2 = s -3/2         2p 1 s -p 2 1 -p 1 p 2 -p 1 p 3 0 0 0 -p 1 p 2 2p 2 s -p 2 2 -p 2 p 3 0 0 0 -p 1 p 3 -p 2 p 3 2p 3 s -p 2 3 0 0 0 0 0 0 r 23 0 0 0 0 0 0 r 31 0 0 0 0 0 0 r 12         B I with      r ij = σ 0 s p i -p j λ i -λ j p i = e 2 λ i σ 0 s = 3 i=1 p i (9)
From the convexity of the exponential function, the r ij are positive if σ 0 is positive. The case λ i = λ j is treated similarly from equation ( 14) and it implies again r ij > 0. The bilinear form associated to this expression can be written as follows, where A is a general symmetric second order tensor:

s 3/2 A. ∂ 2 exp(σ ′ ) ∂σ ′2 .A = p 1 A 2 11 + p 2 A 2 22 + p 3 A 2 33 + 2r 23 A 2 23 + 2r 31 A 2 31 + 2r 12 A 2 12 +2p 2 p 3 (A 2 22 -A 22 A 33 + A 2 33 ) + 2p 3 p 1 (A 2 33 -A 33 A 11 + A 2 11 ) + 2p 1 p 2 (A 2 11 -A 11 A 22 + A 2 22 )(10)
This expression is obviously positive and, added with the convexity of the von Mises term of the criterion, it proves the (strict) convexity of the proposed criterion f . This result can also be obtained without the hessian: it is easy to prove that the second partial derivatives ∂ 2 f /∂σ 2 I are positive thus, from the Ball's theorem (number 5.1) [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], f is convex.

4 Conclusions tension curves. Its convexity is assured by a unique condition (σ 0 > 0), leading to mathematical robustness. Future work may consists in the development of a full damage model using this criterion as a yield surface. Another possible enhancement is to take the high confinement effects into account by replacing the von Mises part of the criterion by, for example, the Drucker-Prager or an elliptic closed expression.

A Mathematics

Let A = A ij -→ e i ⊗ -→ e j be a symmetric second order tensor expressed in an orthonormal basis { -→ e i }; indexes (i, j, k, l) refer to it. Let λ p be the eigenvalues and { -→ v p } the orthonormal basis of eigenvectors of A; indexes (p, q, r, s) refer to it. The derivative of the power low A n expresses as:

∂A n ∂A ijkl = n-1 p=0 A p ik A n-p-1 jl ⇒ ∂A n ∂A pqrs =         
nλ n-1 p if (p = q = r = s) nλ n-1 p if (p = r) = (q = s) et λ p = λ q λ n p -λ n q λp-λq if (p = r) = (q = s) et λ p = λ q 0 else (11)

This allows to express the tensor exponential derivative as: In the base of eigenvectors ( -→ v p ⊗ -→ v q ⊗ -→ v r ⊗ -→ v s ), the derivative (12) has the following components:

exp(A) = ∞ n=0 A n n! ⇒ ∂ exp(A) ∂A ijkl = ∞ n=1 1 n! n-1 p=0 A p ik A n-p-1 jl (12)
∂ exp(A) ∂A pqrs         
exp(λ p ) if (p = q = r = s) exp(λ p ) if (p = r) = (q = s) et λ p = λ q exp(λp)-exp(λq) λp-λq if (p = r) = (q = s) et λ p = λ q 0 else (14)

Figure 1 :

 1 Figure 1: Yield surfaces (MPa) of De Vree (a), Maïolino (b), Willam (c) and the present one (d). Circular dots represent Kupfer's data.

Figure 2 :Figure 3 :

 23 Figure 2: Yield surface in principal stresses space (MPa). The plain lines correspond to the intersections with planes σ I = 0.

  Two interesting properties are obtained from the (double) contraction of the previous expression of the exponential with A on the right side (i.e. A kl ∂(exp(A)) ij /∂A kl ) and with exp(A) on the left side (i.e. (exp(A)) ij ∂(exp(A)) ij /∂A kl ).∂ exp(A)∂A : A = A. exp(