Florent Benaych-Georges
email: florent.benaych@upmc.fr

FINITE DIMENSIONAL PROJECTIONS OF THE BROWNIAN MOTION ON LARGE UNITARY GROUPS

Keywords: . 2000 Mathematics Subject Classification. 15A52, 60B15, 60F05, 46L54. Key words and phrases. Unitary Brownian Motion, Heat Kernel, Random Matrices, Central Limit Theorem, Haar measure

In this paper, we are concerned with the large n limit of the distributions of linear combinations of the entries of a Brownian motion on the group of n × n unitary matrices. We prove that the process of such a linear combination converges to a Gaussian one. Various scales of time and various initial distributions are concerned, giving rise to various limit processes, related to the geometric construction of the unitary Brownian motion. As an application, we propose a quite short proof of the asymptotic Gaussian feature of the linear combinations of the entries of Haar distributed random unitary matrices, a result already proved by Diaconis et al.

Introduction

There is a natural definition of Brownian motion on any compact Lie group, whose distribution is sometimes called the heat kernel measure. Mainly due to its relations with the object from free probability theory called the free unitary Brownian motion and with the two-dimentional Yang-Mills theory, the Brownian motion on large unitary groups has appeared in several papers during the last decade. Rains, in [R97], Xu, in [X97], Biane, in [B97a, B97b] and Lévy and Maïda, in [START_REF]Lévy Schur-Weyl duality and the heat kernel measure on the unitary group[END_REF][START_REF] Lévy | Maïda Central limit theorem for the heat kernel measure on the unitary group[END_REF], are all concerned with the asymptotics of the spectral distribution of large random matrices distributed according to the heat kernel measure. Also, in [D08], Demni makes use of the unitary Brownian motion in the study of Jacobi processes, and, in [START_REF] Benaych-Georges | Lévy A continuous semigroup of notions of independence between the classical and the free one[END_REF], Lévy and the author of the present paper construct a continuum of convolutions between the classical and free ones based on the conjugation of random matrices with a unitary Brownian motion. In this paper, we are concerned with the asymptotic distributions of linear combinations of the entries of an n × n unitary Brownian motion as n tends to infinity. We first give the joint limit distribution, as n tends to infinity, of (possibly rescaled) random processes of the type (Tr[A(V t -I)]) t≥0 for (e -t/2 V t ) t≥0 a Brownian motion staring at I on the group of unitary n × n matrices and A an n × n matrix (Theorem 1.2). This theorem is the key result of the paper, since specifying the choice of the matrices A and randomizing them will then allow us to prove all other results. As a first example, it allows us to find out, for any sequence (α n ) n of positive numbers with a limit α ∈ [0, +∞], the limit distribution of any upper-left corner of n/α n (V log(αnt+1) -I) t≥0 (Corollary 1.4): for small scales of time (i.e. when α = 0), the limit process is purely skew-Hermitian, whereas for large scales of time (α = +∞), the limit process is a standard complex matricial Brownian motion. For intermediate scales of time (0 < α < +∞), the limit process is an interpolation between these extreme cases. The existence of these three asymptotic regimes can be explained by the fact that the unitary Brownian motion is the "wrapping", on the unitary group, of a Brownian motion on the tangent space of this group at I (which is the space of skew-Hermitian matrices), and that as the time goes to infinity, its distribution tends to the Haar measure (for which the upper-left corners are asymptotically distributed as standard complex Gaussian random matrices).

Secondly, we consider a unitary Brownian motion (e -t/2 V t ) t≥0 whose initial distribution is the uniform measure on the group of permutation matrices: its rows are exchangeable, as its columns. In this case, for any positive sequence (α n ) and any positive integer p, the p ×p upper left corner of ( n/α n V log(αnt+1) ) t≥0 converges to a standard complex matricial Brownian motion (Corollary 1.9).

Since the unitary Brownian motion distributed according to the Haar measure at time zero has a stationary distribution, our results allow us to give very short proofs of some well-known results of Diaconis et al, first proved in [START_REF] D'aristotile | Newman Brownian motion and the classical groups. With Probability, Statisitica and their applications[END_REF], about the asymptotic normality of linear combinations of the entries of uniform random unitary matrices (Theorem 1.11 and Corollary 1.12).

It is clear that the same analysis would give similar results for the Brownian motion on the orthogonal group. For notational brevity, we chose to focus on the unitary group.

Let us now present briefly what problems underlie the asymptotics of linear combinations of the entries of a unitary Brownian motion.

Asymptotic normality of random unit vectors and unitary matrices:

The historical first result in this direction is due to Émile Borel, who proved a century ago, in [B1906], that, for a uniformly distributed point (X 1 , . . . , X n ) on the unit euclidian sphere S n-1 , the scaled first coordinate √ nX 1 converges weakly to the standard Gaussian distribution as the dimension n tends to infinity. As explained in the introduction of the paper [ADN03] of Diaconis et al., this says that the features of the "microcanonical" ensemble in a certain model for statistical mecanics (uniform measure on the sphere) are captured by the "canonical" ensemble (Gaussian measure). Since then, a long list of further-reaching results about the entries of uniformly distributed random orthogonal or unitary matrices have been obtained. The most recent ones are the previously cited paper of Diaconis et al., the papers of Meckes and Chatterjee [START_REF] Meckes | Linear functions on the classical matrix groups[END_REF][START_REF] Chatterjee | Meckes Multivariate normal approximation using exchangeable pairs[END_REF], the paper of Collins and Stolz [START_REF] Collins | Borel theorems for random matrices from the classical compact symmetric spaces[END_REF] and the paper of Jiang [J06], where the point of view is slightly different. In the present paper, we give a new, quite short, proof of the asymptotic normality of the linear combinations of the entries of uniformly distributed random unitary matrices, but we also extend these investigations to the case where the distribution of the matrices is not the Haar measure but the heat kernel measure, with any initial distribution and any rescaling of the time.

Second order freeness: A theory has been developed these last five years about Gaussian fluctuations (called second order limits) of traces of large random matrices around their limits, the most emblematic articles in this theory being [START_REF] Mingo | Nica Annular noncrossing permutations and partitions, and second-order asymptotics for random matrices[END_REF][START_REF] Mingo | Speicher Second order freeness and fluctuations of random matrices. I. Gaussian and Wishart matrices and cyclic Fock spaces[END_REF][START_REF] Mingo | Speicher Second order freeness and fluctuations of random matrices. II. Unitary random matrices[END_REF][START_REF] Collins | Speicher Second order freeness and fluctuations of random matrices. III. Higher order freeness and free cumulants[END_REF]. The results of this paper can be related to this theory, even though, technically speaking, we do not consider the powers of the matrices here 1 .

Brownian motion on the Lie algebra and Itô map: The unitary Brownian motion is a continuous random process taking values on the unitary group, which has independent and stationary multiplicative increments. The most constructive way to define it is to consider a standard Brownian motion (B t ) t≥0 on the tangent space of the unitary group at the identity matrix and to consider its image by the Itô map (whose inverse is sometimes called the Cartan map), i.e. to wrap 2 it around the unitary group: the process (U t ) t≥0 obtained is a unitary Brownian motion starting at I. Our results give us an idea of the way the Itô map alterates the Brownian feature of the entries of (B t ) t≥0 at different scales of time. Moreover, the question of the choice of a rescaling of the time (depending on the dimension) raises interesting questions (see Remark 1.1).

Notation. For each n ≥ 1, U n shall denote the group of n × n unitary matrices. The identity matrix will always be denoted by I. For each complex matrix M, M * will denote the adjoint of M. We shall call a standard complex Brownian motion a complex-valued process whose real and imaginary parts are independent standard real Brownian motions divided by √ 2. For all k ≥ 1, the space of continuous functions from [0, +∞) to C k will be denoted by C([0, +∞), C k ) and will be endowed with the topology of the uniform convergence on every compact interval.

Statement of the results

1.1. Brief presentation of the Brownian motion on the unitary group. There are several ways to construct the Brownian motion on the unitary group 3 . For the one we choose here, all facts can easily be recovered by the use of the matricial Itô calculus, as exposed in Section 2.1.

Let n be a positive integer and ν 0 a probability measure on the group of unitary n × n matrices. We shall call a unitary Brownian motion with initial law ν 0 any random process (U t ) t≥0 with values on the space of n × n complex matrices such that U 0 is ν 0 -distributed 1 The reason is that the constant matrices we consider here, like √ n×(an elementary n × n matrix), have no bounded moments of order higher than two: our results are the best ones that one could obtain with such matrices.

2 For G a matricial Lie group with tangent space g at I, the "wrapping" w γ , on G, of a continuous and piecewise smooth path γ : [0, +∞) → g such that γ(0) = 0 is defined by w γ (0) = I and w ′ γ (t) = w γ (t) × γ ′ (t). It B is a Brownian motion on g and (B n ) n≥1 is a sequence of continuous, piecewise affine interpolations of B with a step tending to zero as n tends to infinity, then the sequence w Bn converges in probability to a process which doesn't depend on the choice of the interpolations and which is a Brownian motion on G [IW81, Sect. VI.7], [RW87, Eq. (35.6)], [START_REF] Friz | Rough path limits of Wong-Zakai type with modified drift term[END_REF].

3 See [H56, SV73, IW81, RW87]. A very concise and elementary definition is also given in [R97].

and (U t ) t≥0 is a strong solution of the stochastic differential equation

(1) dU t = idH t U t - 1 2 U t dt,
where (H t ) t≥0 is a Brownian motion4 on the space of n × n Hermitian matrices endowed with the scalar product A, B = n Tr(AB).

It can be proved that for such a process (U t ) t≥0 , for any t 0 ≥ 0, U t 0 is almost surely unitary and both processes (U t 0 +t U * t 0 ) t≥0 and (U * t 0 U t 0 +t ) t≥0 are unitary Brownian motions starting at I n and independent of the σ-algebra generated by (U s ) 0≤s≤t 0 .

Remark 1.1 (Communicated by Thierry Lévy). There are other ways to scale the time for the Brownian motion on the unitary group. Our scaling of the time is the one for which the three limit regimes correspond respectively to small values of t, finite values of t and large values of t and for which the limit non commutative distribution of (U t ) t≥0 is the one of a free unitary Brownian motion. It also has a heuristic geometrical meaning: with this scaling, for any fixed t, the distance5 between U 0 and U t has the same order as the diameter6 of the group. It means that for any fixed t > 0, large values of n, U t is probably no longer too close to its departure point, while it also probably hasn't "orbited" the unitary group too many times.

1.2. The three asymptotic regimes for the unitary Brownian motion starting at I n . Let (α n ) n≥1 be a sequence of positive numbers with a limit α ∈ [0, +∞]. Let us fix a positive integer k and let, for each n ≥ 1, A n 1 , . . . , A n k be a family of non-random n × n matrices and (e -t/2 V n t ) t≥0 be a Brownian motion on U n starting at I. Suppose that there exists complex matrices [a l ] k l=1 , [p l,l ′ ] k l,l ′ =1 and [q l,l ′ ] k l,l ′ =1 such that for all l, l ′ ∈ {1, . . . , k}, we have

1 n Tr(A n l ) -→ n→∞ a l , (2) 1 n Tr(A n l A n l ′ ) -→ n→∞ p l,l ′ , (3) 1 n Tr(A n l A n * l ′ ) -→ n→∞ q l,l ′ . (4)
For each n, let us define, for t ≥ 0,

X n t = α -1/2 n (Tr[A n 1 (V n log(αnt+1) -I)], . . . , Tr[A n k (V n log(αnt+1) -I)]
). Let µ be the probability measure on C([0, +∞), C k ) defined by the fact that any process X t = (X t,1 , . . . , X t,k ) distributed according to µ is a Gaussian centered process with independent increments such that for all t > 0, for all l, l ′ = 1, . . . , k,

(5) E(X t,l X t,l ′ ) = q l,l ′ t , E(X t,l X t,l ′ ) =      -p l,l ′ t if α = 0, -p l,l ′ log(αt+1) α + a l a l ′ log 2 (αt+1) 2α if 0 < α < +∞, 0 if α = +∞.
Theorem 1.2. As n tends to infinity, the distribution of the process (X n t ) t≥0 converges weakly to µ.

Remark 1.3. Theorem 1.2 can easily be extended, using standard topological arguments, to the case where for each n, the matrices A n 1 , . . . , A n k are random, independent of (V n t ) t≥0 and the convergences of (2), (3) and (4) hold in probability (with non-random limits). It will be useful in the proofs of Theorem 1.6 and Theorem 1.11.

Recall that a principal submatrix of a matrix is a matrix obtained by removing some columns, and the rows with the same indices.

Corollary 1.4. Let us fix p ≥ 1 and let (H t ), (S t ) be two independent standard Brownian motions on the euclidian spaces of p × p respectively Hermitian and skew-Hermitian matrices endowed with the respective scalar products X, Y = Tr(XY )/2, X, Y = -Tr(XY )/2. Then, as n tends to infinity, the distribution of the C p×p -valued process of the entries of any p × p principal submatrix of n/α n (V n log(αnt+1) -I) t≥0 converges to the one of the random process (H t-fα(t) + S t+fα(t) ) t≥0 , where

f α (t) =      t if α = 0, log(αt+1) α if 0 < α < +∞, 0 if α = +∞.
Remark 1.5. Note that when α = 0, the limit process is simply a standard Brownian motion on the space of p × p skew-Hermitian matrices, whereas, as α grows from zero to +∞, the Hermitian part of the limit process keeps growing, and at last, when α = +∞, the limit process is a standard Brownian motion on the space of p × p complex matrices. As said in the introduction, the existence of these three asymptotic regimes can be explained by the fact that the unitary Brownian motion is the "wrapping", on the unitary group, of a Brownian motion on the tangent space at I (which is the space of skew-Hermitian matrices), and that as the time goes to infinity, its distribution tends to the Haar measure (for which, as stated by Corollary 1.12, the upper-left corners are asymptotically distributed as standard complex Gaussian random matrices).

1.3. The particular case of unitary Brownian motions with exchangeable rows and columns. Let (α n ) n≥1 be a sequence of positive numbers (no hypothesis is made on its convergence). Let us fix a positive integer k and let, for each n ≥ 1, A n 1 , . . . , A n k be a family of non-random n × n matrices and (e -t/2 V n t ) t≥0 be a Brownian motion on U n such that V n 0 is uniformly distributed on the group of matrices of permutations of {1, . . . , n}.

We suppose that there exists a complex matrix [q l,l ′ ] k l,l ′ =1 for all l, l ′ ∈ {1, . . . , k},

1 n Tr[A n l A n * l ′ ] -→ n→∞ q l,l ′ , (6) 1 n 2 ♯{(i, j) ; ((i, j)-th entry of A n l ) = 0} -→ n→∞ 0 (7)
and that there exists a probability measure µ 0 on C k such that

α -1/2 n (Tr[A n 1 V n 0 ], . . . , Tr[A n k V n 0 ]) -→ n→∞ µ 0 (conv. in distribution). ( 8
)
For each n, let us define, for t ≥ 0,

X n t = α -1/2 n (Tr(A n 1 V n log(αnt+1) ), . . . , Tr(A n k V n log(αnt+1) )
). Let P be an Hermitian matrix such that P 2 = [q l,l ′ ] k l,l ′ =1 , (Z 1 , . . . , Z k ) be an independent family of standard complex Brownian motions and C be a µ 0 -distributed random variable, independent of the Z l 's. Let us define µ to be the distribution, on C([0, +∞), C k ), of the process (C + (Z 1,t , . . . , Z k,t )P ) t≥0 .

Theorem 1.6. As n tends to infinity, the distribution of X n converges weakly to µ.

Remark 1.7. As Theorem 1.2, Theorem 1.6 can be extended to the case where for each n, the matrices A n 1 , . . . , A n k are random, independent of (V n t ) t≥0 and the convergences of ( 6) and (7) hold in probability (with non-random limits). In several examples given below, the matrices A n l 's are actually random. Examples 1.8. Let us give a few examples of sequences A n which satisfy the hypotheses (6), ( 7) and (8) (or their probabilistic versions mentioned in Remark 1.7). a) Firstly, if (α n ) is bounded from below by a positive constant, if (6) holds and if for all l, n -1 ♯{(i, j) ; ((i, j)-th entry of A n l ) = 0} -→ 0, then (8) holds for µ 0 the Dirac mass at zero (this can easily be deduced from Lemma 2.3 bellow). b) Secondly, if, α n tends to one and if, for all n, the matrices A n l , 1 ≤ l ≤ k, are random real n × n matrices which satisfy (6) and (7) for the convergence in probability, then (8) holds for µ 0 the Gaussian measure with covariance matrix [q l,l ′ ] k l,l ′ =1 . This follows easily from [S88]. As an example, if k = 1 and if, for each n, the entries of A n 1 are i.i.d. with distribution τ n such that as n tends to infinity,

τ n ({0}) -→ 1 and n(1 -τ n ({0})) t 2 dτ n (t) -→ 1, then (6)
, ( 7) and (8) hold for q 1,1 = 1 and µ 0 the standard Gaussian law. c) Other examples can be found using [C78, Th. 5.1], where the laws µ 0 are other infinitely divisible laws.

Both following corollaries are direct applications of the previous theorem, the first one using implicitly the fact that any entry of V n 0 is null with probability 1 -n -1 , and the second one using implicitly the fact that the distribution of the number of fixed points of a uniform random permutation of {1, . . . , n} converges weakly, as n tends to infinity, to the Poisson distribution with mean one [START_REF] Diaconis | Shahshahani On the eigenvalues of random matrices[END_REF].

Corollary 1.9. Let (α n ) be a sequence of positive numbers. For any p, q ≥ 1, as n tends to infinity, the distribution of any p × q submatrix of n/α n V n log(αnt+1) t≥0 converges weakly to the one of an independent family of pq standard Brownian motions on the complex plane, i.e. a Brownian motion on the space of p × q complex matrices.

Corollary 1.10. As n tends to infinity, the distribution of (Tr(V n log(t+1) )) t≥0 converges weakly to the one of (C + Z t ) t≥0 , where C is a Poisson random variable with mean one and (Z t ) t≥0 is a standard complex Brownian motion, independent of C.

1.4. Application to the asymptotics of the uniform measure on the unitary group. Since the Brownian motion on the unitary group distributed according to the Haar measure at time zero has a stationary distribution, our results allow us to recover certain results of asymptotic normality of linear combinations of the entries of uniform random unitary matrices.

The following theorem is not new [START_REF] Collins | Borel theorems for random matrices from the classical compact symmetric spaces[END_REF][START_REF] Chatterjee | Meckes Multivariate normal approximation using exchangeable pairs[END_REF]. However, our method allows to give a very direct proof, even under these very general hypotheses.

Theorem 1.11. Let us fix k ≥ 1 and let, for each n ≥ 1, A n 1 , . . . , A n k be n × n nonrandom matrices and U n be a random matrix with uniform distribution on the group of n × n unitary matrices. Suppose that for all l, l ′ , there is q l,l ′ ∈ C such that

1 n Tr(A n l A n * l ′ ) -→ n→∞ q l,l ′ .
Then as n tends to infinity, the distribution of the random vector (Tr[A n 1 U n ], . . . , Tr[A n k U n ]) converges weakly to the one of a Gaussian centered family (Z 1 , . . . , Z k ) of complex random variables such that for all l, l ′ = 1, . . . , k, E(Z l Z l ′ ) = 0 and E(Z l Z l ′ ) = q l,l ′ .

The following corollary is immediate.

Corollary 1.12. Let us fix k ≥ 1 and let, for each n, Z n 1 , . . . , Z n k be k different entries of a random n × n matrix with uniform distribution on the unitary group. Then the joint distribution of √ n(Z n 1 , . . . , Z n k ) converges weakly, as n tends to infinity, to the one of a family of independent standard complex Gaussian random variables.

Proofs

2.1. Preliminaries on matricial Itô calculus. a) Let n be a positive integer. Let (F t ) t≥0 be a filtration and (H n t ) t≥0 be an (F t ) t≥0 -standard Brownian motion on the space of n × n Hermitian matrices endowed with the scalar product A, B = n Tr(AB), i.e. a process with values in the space of n × n Hermitian matrices such that the diagonal and upper diagonal entries of ( √ nH n t ) t≥0 are independent random processes, the ones on the diagonal being standard real Brownian motions and the ones above the diagonal being standard complex Brownian motions. We shall use these formulas many times in the paper, without citing them every time.

b) With the same notation, let us now consider a deterministic C 1 function f with positive derivative such that f (0) = 0. Then for X the process introduced above, the process Xt := X f (t) satisfies

d Xt = f ′ (t)A f (t) (d Hn t )B f (t) + f ′ (t)C f (t)
dt, where Hn is the (F f (t) ) t≥0 -Brownian motion defined by the formula Hn

t = t 0 1 √ f ′ (s) d(B f (•) ) s .
2.2. Proof of Theorem 1.2. Let us first state some matricial inequalities we shall often refer to in the following. Let X, Y be two complex matrices and G, H be two Hermitian nonnegative matrices. Then we have 11) Inequality (11) follows from ( 9) and (10), which are obvious.

| Tr(XY )| ≤ Tr(XX * ) Tr(Y Y * ), (9) Tr(G 2 ) ≤ (Tr G) 2 , (10) | Tr(GH)| ≤ Tr(G) Tr(H). (
Lemma 2.1. Let us fix n ≥ 3, an n × n matrix A and a Brownian motion (e -t/2 V t ) t≥0 on U n starting at I. Then there exists some real numbers C 1 , C 2 , C 3 , C 4 , independent of t, whose absolute values are bounded by 100(Tr AA * ) 2 /n 2 and such that for all t ≥ 0,

E[Tr(AV t AV t )] = Tr(A 2 ) cosh(t/n) -(Tr A) 2 sinh(t/n), (12) E[| Tr(AV t AV t )| 2 ] = | Tr(A 2 )| 2 + | Tr A| 4 2n 4 + C 1 n 3/2 cosh 2t n -1 (13) -ℜ[(Tr A) 2 Tr(AA)] + n 3/2 C 2 sinh 2t n +nC 3 (e t -1) + C 4 (e 2t -1).
Proof. Since the formulas we have to state are invariant under multiplication of A by a scalar, we can suppose that Tr(AA * ) = n.

Note that for H n t as in Section 2.1, (V t ) is a strong solution of dV t = i(dH n t )V t . Hence by the matricial Itô calculus,

d Tr(AV t AV t ) = 2i Tr(V t AV t AdH n t ) -n -1 Tr(AV t ) Tr(AV t )dt, d Tr(AV t ) Tr(AV t ) = 2i Tr(AV t ) Tr(V t AdH n t ) -n -1 Tr(AV t AV t )dt. It follows that for x(t) = E[Tr(AV t AV t )] and y(t) = E[Tr(AV t ) Tr(AV t )], we have x ′ = -n -1 y and y ′ = -n -1 x.
Equation (12) follows. Let us now prove (13). We introduce the functions

f (t) = E(| Tr(AV t AV t )| 2 ), g(t) = ℜ{E[Tr(AV t ) Tr(AV t ) Tr(A * V * t A * V * t )]}, h(t) = E[| Tr(AV t )| 4 ].
By the matricial Itô calculus again (using the hypothesis Tr(AA

* ) = n), n × f ′ (t) = -2g(t) + 4e t u A * A,AA * (t), (15) n × g ′ (t) = -f (t) -h(t) + 4e t ℜ[v A,A * AA * (t)], n × h ′ (t) = -2g(t) + 4ne t v A,A * (t).
It follows, by ( 14), that g

′′ (t) -4 n 2 g(t) = 8e 2t n κ + 4e t n θ, for κ = n -1 Tr(AA * AA * ) -1, θ = 2 -n -1 Tr(AA * AA * ) -n -1 Tr(AAA * A * ) -| Tr A| 2 + ℜ{Tr A Tr(A * AA * )}, hence g(t) = µ cosh 2t n + ν sinh 2t n + 2nκ n 2 -1 e 2t + 4nθ n 2 -4 e t , with µ = ℜ[Tr A Tr A Tr(A * A * )] - 2nκ n 2 -1 - 4nθ n 2 -4 , ν = - 1 2 | Tr A 2 | 2 - 1 2 | Tr A| 4 + 2ℜ[Tr A Tr(A * AA * )] - 2n 2 κ n 2 -1 - 2n 2 θ n 2 -4 .
From (15), it follows that

f (t) = -| Tr A 2 | 2 -µ sinh 2t n -ν cosh 2t n -1 + w(t),
where w(t) = -2κ n 2 -1 (e 2t -1) -8θ n 2 -4 (e t -1) + 2(e 2t -1) + 4(n -1 Tr(AA * AA * ) -1)(e t -1). Now, the conclusion follows from the fact that since Tr(AA * ) = n, the inequalities (9), (10) and (11) allow to prove that

| Tr A|, | Tr A 2 | ≤ n, Tr(AA * AA * ), | Tr(AAA * A * )| ≤ n 2 and | Tr(A * AA * )| ≤ n 3/2 . Lemma 2.2. Let (a n ), (b n ), (c n ) be sequences of real numbers such that (a n ) tends to +∞ and (b n ) and (c n ) are both bounded. Then we have (16) u n := n 2 b n a n cosh log a n n -1 + nc n a n sinh log a n n -→ n→∞ 0.
Proof. Let us define

K = max sup n≥1 |b n | , sup n≥1 |c n | , sup 0<x≤1 cosh(x) -1 x 2 , sup 0<x≤1 sinh(x) x .
Then (16) follows from the fact that, since u n can also be written

u n = n 2 b n 2a 1-1 n n + n 2 b n 2a 1+ 1 n n - n 2 b n a n + nc n 2a 1-1 n n - nc n 2a 1+ 1 n n ,
we have the upper-bound:

|u n | ≤ 5K n 2
e n-1 ½ an>e n + 2K 2 log 2 an+log an an ½ an≤e n .

Proof of Theorem 1.2. For all n, X n is a C k -valued continuous centered martingale. To prove that its distribution tends to µ, by Rebolledo's Theorem (see [START_REF] Papanicolaou | Martingale approach to some limit theorems[END_REF] or [AGZ09, Th. H.14]), it suffices to prove that the bracket of X n converges pointwise, in L 1 , to the one of a µ-distributed process. Hence it suffices to fix λ 1 , . . . , λ k ∈ C, to define the process

Y n t = λ 1 X n t,1 + • • • + λ n X n t
,k and to prove that as n tends to infinity,

(17) Y n , Y n t L 2 -→ qt and Y n , Y n t L 2 -→      -pt if α = 0, -p log(αt+1) α + a 2 log 2 (αt+1) 2α if 0 < α < +∞, 0 if α = +∞.
for q = k l,l ′ =1 λ l λ l ′ q l,l ′ , p = k l,l ′ =1 λ l λ l ′ p l,l ′ and a = k l=1 λ l a l . Let us define, for each n, A n = k l=1 λ l A n l . We have

Y n t = α -1/2 n Tr[A n (V n log(αnt+1) -I)].
Hence by Section 2.1 b), Y n satisfies

dY n t = i √ α n t + 1 Tr[V n log(αnt+1) A n dH n t ],
where H n t is an Hermitian Brownian motion as introduced in Section 2.1. Thus, since V n log(αnt+1) V n * log(αnt+1) = (α n t + 1)I, by the matricial Itô calculus, we have

d Y n , Y n t = 1 n Tr(A n A n * )dt,
so that the first part of (17) follows directly from (4). Let us now prove the second part. By the matricial Itô calculus again, we have

d Y n , Y n t = -1 n(α n t + 1) Tr[A n V n log(αnt+1) A n V n log(αnt+1) ]dt.
Hence it suffice to prove that as n tends to infinity, we have the convergence

(18) -1 n(α n t + 1) Tr[A n V n log(αnt+1) A n V n log(αnt+1) ] L 2 -→ -p αt+1 + a 2 log(αt+1) αt+1 if 0 ≤ α < +∞, 0 if α = +∞,
uniformly as t varies in any compact subset of [0, +∞). This follows easily from Lemma 2.1 (with Lemma 2.2 in the case where α = +∞).

2.3. Proof of corollary 1.4. By Theorem 1.2 and formula (5), applied with matrices A n l of the type √ n×(an elementary n × n matrix), the distribution of any p × p principal submatrix of n/α n (V n log(αnt+1) -I) t≥0 converges weakly to the one of the random process (M t ) t≥0 with independent increments such that for all p × p complex matrices X, Y and all t ≥ 0,

E[Tr(M t X)Tr(M t Y )] = t
≤ n -1 + √ C A C B n(n -1)
Tr(AA * ) Tr(BB * ),

where for each matrix X = [x i,j ] n i,j=1 , C X denotes ♯{(i, j) ; x i,j = 0}.

Proof. Let us denote by respectively [a i,j ], [b i,j ], [s i,j ] the entries of A, B, S. By Hölder's inequality, we have

E{| Tr(AS)|} ≤ i,j |a i,j |E(s i,j ) = 1 n i,j |a i,j | ≤ 1 n C A Tr(AA * ).
Moreover, we have

E{| Tr(ASBS)|} ≤ j,k |a k,j b k,j |E(s 2 j,k ) + i,j,k,l j =l or k =i |a i,j b k,l |E(s j,k s l,i ) ≤ 1 n j,k |a k,j b k,j | + 1 n(n -1) i,j,k,l |a i,j b k,l |,
and the conclusion follows from Hölder's inequality again.

Proof of Theorem 1.6. We consider C, Z 1 , . . . , Z k , P as introduced above the statement of the theorem. For each t ≥ 0, let us define K t = (Z 1,t , . . . , Z k,t )P . It suffices to prove that as n tends to infinity, the joint distribution of (X n 0 , (X n t -X n 0 ) t≥0 ) converges weakly to the one of (C, (K t ) t≥0 ).

First, Lemma 2.3 and the hypothesis (6) and (7) allow us to claim that as n tends to infinity, 1 n Tr(A n,l V n,0 ) and 1 n Tr(A n,l V n,0 A n,l ′ V n,0 ) both converge to zero in probability. Moreover, by a standard topological argument, one can suppose that (α n ) admits a limit α ∈ [0, +∞]. Now, note that for all n, (e -t/2 V n * 0 V n t ) t≥0 is a unitary Brownian motion starting at I, independent of V n 0 . By Theorem 1.2 and Remark 1.3, it implies that the joint distribution of (X n 0 , (X n t -X n 0 ) t≥0 ) converges weakly to the one of (C, (K t ) t≥0 ), which closes the proof of the theorem. Set U = [u i,j ] n i,j=1 . One can write A = V A ′ W , with V, W unitary matrices and A ′ a diagonal matrix whose diagonal entries a 1 , . . . , a n are the eigenvalues of √ * . Since the law of U is invariant under the left and right actions of the unitary group, one can suppose that A = A ′ . These invariances of the law of U also imply that for all i, j, E(u i,i u j,j ) = δ j i /n. Equation ( 19) follows. Equation (20) follows from (13) and the fact that the Haar measure on the unitary group is the limit of the distribution of e -t/2 V t as t tends to infinity.

Proof of Theorem 1.11. Step I. Firstly, by Lemma 2.4, for all l, l ′ , n -1 Tr(A n l U n ) and n -1 Tr(A n l U n A n l ′ U n ) both tend in probability to zero as n tends to infinity.

Step II. Following [DZ98, Th. D8], we fix a bounded real function f on C k which is 1-Lipschitz for the canonical hermitian norm || • || on C k , and we shall prove that Let us fix ε > 0.

a) Let, for each t ≥ 0, (Z 1,t , . . . , Z k,t ) be a Gaussian family of centered complex random variables such that for all l, l ′ , E[Z l,t Z l ′ ,t ] = 0 and E[Z l,t Z l ′ ,t ] = q l,l ′ (1 -e -t ) (such a family exists because the matrix [q l,l ′ ] l,l ′ is nonnegative). The distribution of (Z 1,t , . . . , Z k,t ) tends to the one of (Z 1 , . . . , Z k ) as t tends to infinity. Hence there is t 0 > 0 such that |E{f (Z 1 , . . . , Z k )} -E{f (Z 1,t 0 , . . . , Z k,t 0 )}| ≤ ε, (22)

e -t 0 2 sup n≥1 {n -1 Tr[A n 1 A n * 1 + • • • + A n k A n * k ]} 1/2 ≤ ε. (23)
b) For each n, up to an extension of the probability space where U n is defined, one can suppose that there exists a unitary Brownian motion (U n t ) t≥0 , starting at I, independent of U n . Let us define, for each n, 2 U n t 0 -I)]) By Step I and the randomized version of Theorem 1.2 stated in Remark 1.3, as n tends to infinity, the distribution of D n converges to the one of (Z 1,t 0 , . . . , Z k,t 0 ). It follows that for n large enough, (24) |E{f (Z 1,t 0 , . . . , Z k,t 0 )} -E{f (D n )}| ≤ ε. c) At last, since f is 1-Lipschitz for || • ||, for all n, we have

|E{f (D n )} -E{f (X n )}| ≤ E{||D n -X n ||} ≤ (e -t 0 E{| Tr[A n 1 U n ]| 2 + • • • + | Tr[A n k U n ]| 2 }) 1/2 ≤ e -t 0 /2 (n -1 Tr[A n 1 A n * 1 ] + • • • + n -1 Tr[A n k A n * k ])

  Now, let us define, for C, D some n× n matrices, u C,D (t) = E[Tr(V t CV * t D)] and v C,D (t) = E[Tr(V t C) Tr(V * t D)].By the matricial Itô calculus again, one has (14) u C,D (t) = 1 n (e t -1) Tr C Tr D+Tr(CD), v C,D (t) = 1 n (e t -1) Tr(CD)+Tr(C) Tr(D).

2. 5 .

 5 Proof of Theorem 1.11. Lemma 2.4. Let U be Haar-distributed on U n and A be an n × n matrix, with n ≥ 3. Then E{| Tr(AU)| 2 } = Tr(AA * )/n, (19) E{| Tr(AUAU)| 2 } ≤ 100(Tr(AA * )) 2 /n 2 . (20) Proof.

  E{f (Tr[A n 1 U n ], . . . , Tr[A n k U n ])} -→ n→∞ E{f (Z 1 , . . . , Z k )}.

X

  n = (Tr[A n 1 U n U n t 0 ], . . . , Tr[A n k U n U n t 0 ]) Y n = e -t 0 2 (Tr[A n 1 U n ], . . . , Tr[A n k U n ]) D n = X n -Y n = e -t 0 2 (Tr[A n 1 U n (e t 0 2 U n t 0 -I)], . . . , Tr[A n k U n (e t 0

  ) t≥0 -adapted matrix-valued processes A, B, C, D, E, F , then, by Itô's formula, d(XY ) t = (dX t )Y t + X t dY t + 1 n Tr(B t D t )A t E t dt,

	X, Y such that dX d Tr(X), Tr(Y ) t =	If one considers two matrix-valued semimartingales n 1 Tr(B

t = A t (dH n t )B t + C t dt, dY t = D t (dH n t )E t + F t dt,

for some (F t t A t E t D t )dt.

  Tr(XY * ), E[Tr(M t X)Tr(M t Y )] = -Tr(XY ) or 0. Since a standard Brownian motion (B t ) on an euclidian space (E, •, • ) satisfies, for all u, v ∈ E, E( B t , u B t , v ) = u, v t, the result can easily be verified. Let us fix n ≥ 2, let S be the matrix of a uniform random permutation of {1, . . . , n} and let A, B be n × n non-random matrices. Then we have

			log(αt + 1) α	Tr(XY ),
	where if α = 0 or +∞, the second right-hand term has to be replaced by respectively
	-t 2.4. Proof of Theorem 1.6.
	Lemma 2.3. E{| Tr(AS)|} ≤	1 n	C A Tr(AA

* ), E{| Tr(ASBS)|}

  1/2 , the last inequality following from(19). By (23), it allows us to claim that for all n,(25) |E{f (D n )} -E{f (X n )}| ≤ ε. d)To conclude, note that by the right invariance of the law of U n , X n and the random vector (Tr[A n 1 U n ], . . . ,Tr[A n k U n ]) have the same distribution. By (22), (24) and (25), it follows that for n large enough,|E{f (Z 1 , . . . , Z k )} -E{f (Tr[A n 1 U n ], . . . , Tr[A n k U n ])}| ≤ 3ε.It closes the proof of the theorem.

(iH t ) t≥0 is in fact the skew-Hermitian Brownian motion that the process (U t U * 0 ) t≥0 wraps around the unitary group, as explained in Footnote 2.

As explained in Footnotes 2 and 4, (U t ) t≥0 is the wrapping of (iH t ) t≥0 on the unitary group, hence the distance between U 0 and U t has the same order as the one between 0 and iH t , which, by the Law of Large Numbers, has order √ tn.

By definition, the diameter of the group is the supremum, over pairs U, V of unitary matrices, of the length of the shortest geodesic between U and V . Here, it is equal the length of the geodesic t ∈ [0, nπ] → exp(itI n /n) between I n and -I n , i.e. to nπ.

This work was partially supported by the Agence Nationale de la Recherche grant ANR-08-BLAN-0311-03.


Marek Bucki
email: marek.bucki@imag.fr

Claudio Lobos
email: clobos@inf.utfsm.cl

Yohan Payan
email: yohan.payan@imag.fr

Nancy Hitschfeld

Jacobian-based Repair Method for Finite Element Meshes after Registration

Keywords: Meshing, Mesh Repair, Registration, Finite Element Method, Quality Improvement

come

Introduction

The Finite Element Method (FEM) [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Hughes | The finite element method: Linear static and dynamic finite element analysis[END_REF] is a numerical technique computing an approximate solution to a set of Partial Differential Equations (PDEs) and therefore, it is widely used to simulate mechanical behaviors.

The FEM relies on a subdivision of the continuum domain Ω to be simulated. This subdivision, called mesh, is the partition of Ω into simpler geometrical bodies known as elements. The FEM provides an approximation of the solution to the system of PDEs at each node of the mesh. This approximate solution can be further computed at any point p ∈ Ω by interpolation of the values found at the nodes of the element containing p.

The context of our work is biomechanical modeling of human organs for biomedical applications. In this field, the most common way to produce a mesh is to perform the following steps:

-To obtain volumetric medical images of the patient's organ using Magnetic Resonance Imaging (MRI), Ultrasound (US), Computed Tomography scanner (CT) or other imaging techniques. -To achieve image segmentation in order to produce a cloud of points or an initial surface mesh defining the geometry of the modelled domain. -To add internal nodes to produce the 3D elements that will conform the final volumetric mesh.

Some examples of such mesh generators are [START_REF] Molino | A crystalline, red green strategy for meshing highly deformable objects with tetrahedra[END_REF][START_REF] George | An efficient algorithm for 3d adaptive meshing[END_REF][START_REF] Alliez | Variational tetrahedral meshing[END_REF] for tetrahedral meshes. Others like [START_REF] Zhang | Adaptive and quality quadrilateral/hexahedral meshing from volumetric data[END_REF][START_REF] Nesme | Milieu mecanique deformable multiresolution pour la simulation interactive[END_REF] are interested in hexahedral meshes and finally other works such as [START_REF] Lobos | Mixed-element mesh for an intra-operative modeling of the brain tumor extraction[END_REF][START_REF] Kwok | A simple and effective mesh quality metric for hexahedral and wedge elements[END_REF] focus on the production of mixed-element meshes.

Even though those techniques are suited for most of the cases, in the medical field it is sometimes preferred to work with "registration" techniques. Consider for instance the concept of image registration which is the process of transforming different sets of data into one coordinate system. In the case of mesh registration, the different sets of data come from a generic model, named the Atlas, and patient data. The Atlas is normally a predefined volume mesh conserving all the relevant data for the simulation. The patient data is normally a set of points or a surface mesh. The first step is to produce rigid registration, which is, as in the case of imaging, putting both in the same coordinate system. The second step is to produce elastic registration, which is to deform the Atlas in order to represent patient data. The overall process is then called "mesh registration".

There exist several mesh registration techniques, but to continue with the explanation, the so called Mesh-Matching (M-M) algorithm introduced by Couteau et al. in [START_REF] Couteau | The mesh-matching algorithm: an automatic 3d mesh generator for finite element structures[END_REF] will be used. The statement is that building a mesh of an organ is a complex task. To achieve accurate surface representation and good quality while limiting the number of elements in the mesh is somehow, a difficult compromise. Therefore in many cases, meshes are built by hand. As this task is time consuming and usually requires expert user interaction, it cannot be easily put in practice for a great number of patients. Therefore an interesting alternative is to use a pre-defined generic FE mesh (the Atlas mesh) of the target organ in order to achieve accurate representation of its structure.

FE mesh construction relies on patient data obtained from segmented medical images. The resulting data can either be a cloud of points or a surface mesh of the organ. Using this representation a registration process fits the surface nodes of the Atlas onto the target organ by applying a transformation function T . In other words, an elastic registration function is an application T : R 3 → R 3 that minimizes a given registration error, or better said 'registration energy', between a given source points set S and a destination points set D. The registration energy is a 'similarity' measure between the transformed source points set T (S ) and the target destination points set D. If the registration energy and function search space are not carefully defined, finding T is an ill-posed problem, i.e. the existence and uniqueness of the solution are not guaranteed.

The registration energy can be freely defined in order to reflect the nature of the registration problem. The optimal registration is achieved when this energy reaches zero i.e. when each transformed source point has found a counterpart within the destination set. Given the sparse and uneven nature of the data, the search for the elastic transformation must be constrained in order to avoid erratic registrations. Various approaches to produce a registration have been proposed in the literature [START_REF] Couteau | The mesh-matching algorithm: an automatic 3d mesh generator for finite element structures[END_REF][START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Schnabel | A generic framework for non-rigid registration based on non-uniform multi-level free-form deformations[END_REF][START_REF] Montagnat | Globally constrained deformable models for 3d object reconstruction[END_REF], the discussion of which falls out of the scope of this article. Note that T is computed using surface nodes S (of Atlas volume mesh) and the points set D (patient data). Moreover, T is a 3D transformation function so it can be applied to Atlas internal nodes and so a final 3D mesh of the target organ is produced (patient mesh). Unfortunately it has been reported [START_REF] Luboz | Orbital and maxillofacial computer aided surgery: patient-specific finite element models to predict surgical outcomes[END_REF] that, in some cases, the application of T to internal nodes can produce invalid or poor quality elements (to be defined in §2). Fig. 1 shows an example of M-M (here illustrated in 2D for clarity) and suggests how invalid and poor quality elements can be created.

Note that invalid and poor quality elements cannot only be produced by registration methods [START_REF] Castellano-Smith | Constructing patient specific models for correcting intraoperative brain deformation[END_REF][START_REF] Berar | 3d meshes registration : Application to statistical skull model[END_REF][START_REF] Montagnat | 4d deformable models with temporal constraints: application to 4d cardiac image segmentation[END_REF], but also by meshing techniques such as the whisker-weaving from CUBIT1 [START_REF] Knupp | Hexahedral mesh untangling & algebraic mesh quality metrics[END_REF].

The goal of the present work is to repair the invalid and poor quality elements that might result from the process of mesh registration. The proposed method is based on the computation of the determinant of the Jacobian matrix and to this end section 2 defines what are invalid and a poor quality elements. Section 3 gives the State of the Art regarding mesh repair and quality improvement in our scope. Section 4 explains in further detail our proposed solution, in particular the Jacobian formulation, the penalty functions that are used, our algorithm and a practical example. Finally sections 5 and 6 provide our conclusions and perspectives.

Definition of invalid and poor quality elements

This section introduces the concepts of invalid and poor quality elements. It is important to mention that the two definitions proposed here rest upon the Jacobian matrix of the mapping function F relating the parent or reference element frame to its actual coordinates within the mesh. In the remainder of the article this matrix is denoted as J. Note that other geometrical functions can be used to define the quality of an element [START_REF] Shewchuk | What is a good linear element? interpolation, conditioning, and quality measures[END_REF][START_REF] Field | Qualitative measures for initial meshes[END_REF][START_REF] Kwok | A simple and effective mesh quality metric for hexahedral and wedge elements[END_REF][START_REF] Lobos | Mixed-element mesh for an intra-operative modeling of the brain tumor extraction[END_REF].

Invalid element

As mentioned before, the FEM computes for each node of the underlying mesh an approximated solution to the system of PDEs describing the physical phenomenon under study. For this purpose, a perfect element is defined in a reference parent system (ξ 1 , ξ 2 , ξ 3 ) where each reference element point is related to its actual counterpart within the element in the modeled domain referential (x 1 , x 2 , x 3 ) through the mapping F , as shown in Fig. 2. The Jacobian matrix of mapping F considered at parent frame point ξ is defined as:

J(ξ) := ∂F ∂ξ (ξ)
The FE analysis can be carried out on a mesh as long as the mapping F , considered for each mesh element, remains a one-to-one application. This property can be assessed by considering the determinant of Jacobian matrix of F , |J(ξ)|, at a number of element "control points", such as the element nodes or Gaussian quadrature (integration) points, for all the elements forming the mesh. If a Jacobian is found to be negative,

P N 0 N 1 N 2 N 3 N 4 N 5 N 6 N 7 x x x 3 2 1 F ξ 1 ξ 2 ξ 3
Fig. 2 The reference-to-actual frame mapping from which the Jacobian matrix is derived. Left: the reference/parent space and right: the normal/actual mesh space.

then the mapping F for the considered element is not a bijection, and the element cannot be used for FE analysis. In this case, the mesh is not not suitable for Finite Element analysis as the singularity of the mapping F leads to modeling inconsistency.

Let's consider this situation from a geometrical point of view as illustrated in Fig. 3. If nodes b, c and d are connected to node a then an imaginary plane can be formed by nodes b, c and d. The normal n of this plane is computed as (cb) × (db). After normalizing the vector n the signed distance of a to the plane can be computed as: (ab) • n. If this distance is positive, then the position of a is valid within this element. If not, then this element is considered invalid for the FEM at node a. 3 If the signed distance of node a to the plane formed by "connected nodes" b, c and d is positive, the element is valid (left), in any other case the entire element is said to be invalid (right).

Note that the value of |J(ξ)| is computed for each element at each of its nodes. If node a is shared by two elements e 1 and e 2 , it is possible to devise a nodal configuration where the position of a is valid within e 1 but invalid within e 2 .

Poor quality element

The Jacobian value measures the distortion of the actual mesh element w.r.t. its reference configuration, but does not carry overall distortion information. To overcome this issue the quality of an element can be measured at each node n as a ratio between the nodal Jacobian value at n, and a global distortion information. The latter is given by the maximal Jacobian value among those measured at all element nodes. This element node quality criterion is called Jacobian Ratio or JR [START_REF] Kwok | A simple and effective mesh quality metric for hexahedral and wedge elements[END_REF].

Mathematically, the JR of node ξ 0 taken within element e, JR ξ 0 e , is computed as:

JR ξ 0 e := |J(ξ 0 )| / max ξ∈e {|J(ξ)|}
Considering the case presented in Fig. 3, the value of |J(a)|, computed at node a within the element, would be very small if a was close to plane bcd and therefore, so would be the value of the JR a . Actually, when JR a value is positive, but close to 0 then the element presents a poor quality at node a. When JR increases, so does the quality of the element until it reaches the ideal value of 1.

As it is impossible to converge to 1 all the JR values within the mesh, a quality acceptance threshold must be defined. To this end, we have chosen to follow the recommendations adopted in the commercial FE analysis software ANSYS (ANSYS Inc., USA), where an element e is said to present poor quality when any node ξ within e has a JR value lower than 0.03 [START_REF] Kelly | Element shape testing[END_REF]. A good quality element, thus, has all its JR values within the interval [0.03, 1].

Overview of mesh repair techniques proposed in the literature

In this work, our goal is to achieve validity and improve the quality of the elements in a mesh after applying a non-rigid registration in the context of mesh generation procedures such as the Mesh-Matching introduced in section 1. To do so, we focus specifically on the Jacobian matrix as it is at the root of many popular mesh rejection criteria implemented in commercial FE analysis softwares.

In [START_REF] Freitag | On combining laplacian and optimization-based mesh smoothing techniques[END_REF], Freitag proposed the use of Laplacian-based algorithms to improve the overall quality of meshes. In particular, the Smart-Laplacian algorithm proposed in [START_REF] Freitag | On combining laplacian and optimization-based mesh smoothing techniques[END_REF] changes the position of a given node only if it results in overall quality improvement. This filter is implemented in such a way that no invalid element can be generated within the mesh. Following the same idea, Freitag and Plassmann proposed in [START_REF] Freitag | Local optimization-based simplicial mesh untangling and improvement[END_REF] the combination of Smart-Laplacian with optimization algorithms, having a much higher computational cost than the Laplacian filter alone, in order to achieve mesh untangling. In this approach, invalid mesh configurations are eliminated by mesh nodes reallocation and the applied optimization maximizes the minimum volume of tetrahedra in a local sub-mesh.

Both techniques presented above are designed to repair tetrahedral meshes, and the idea of expanding this concept to hexahedra and wedges was introduced by Knupp in [START_REF] Knupp | Hexahedral mesh untangling & algebraic mesh quality metrics[END_REF][START_REF] Knupp | Achieving finite element mesh quality via optimization of the jacobian matrix norm and associated quantities. part ii -a framework for volume mesh optimization and the condition number of the jacobian matrix[END_REF]. The main principle here is that the reallocation of a given node will only change the Jacobian values of the nodes connected to it (as showed in subsection 2.1). Moreover, for the tetrahedron, wedge and hexahedron each node is connected to 3 nodes within each element it is part of, therefore for those elements, a Jacobian-based repair method can always be geometrically reduced to an optimization based on tetrahedra.

Paoletti in [START_REF] Paoletti | Polyhedral mesh optimization using the interpolation tensor[END_REF] expanded the work to polyhedral meshes by including the "interpolation tensor" to the optimization.

Another important work was carried out by Kwok and Chen in [START_REF] Kwok | A simple and effective mesh quality metric for hexahedral and wedge elements[END_REF] where the focus was on the quality of hexahedra and wedges. The Jacobian ratio was used to improve the quality in conjunction with other quality functions such as the aspect-ratio, the warping factor and the control number [START_REF] Knupp | Matrix norms & the condition number, a general framework to improve mesh quality via node-movement[END_REF].

The main issue with the above mesh untangling techniques is that they are not focused on the strategy of "mesh registration". Lets recall that this type of strategy intend to match a predefined mesh model to a target domain. If the "registration" is accurate, the goal is to displace as less as possible the nodes in order to achieve mesh validity and acceptable quality levels. The above techniques are focused on overall element quality and therefore they have no constraint on reallocating every node, which in our case may cause an accuracy loss in the representation of target domain. Moreover, they do not separate validity recovery and quality improvement procedures, intending to achieve it all at the same time. The problem with this approach is that if invalid node a and poor quality node b are connected, finding a new position to improve the quality of b might decrease the chances of finding a new valid position for a.

A completely different approach was proposed by Luboz et al. in [START_REF] Luboz | Automatic 3d finite element mesh generation: an atlas fitted onto patient data[END_REF]. This work was considered as a continuation of the Mesh Matching (M-M) technique, initially presented in [START_REF] Couteau | The mesh-matching algorithm: an automatic 3d mesh generator for finite element structures[END_REF], as it repairs invalid configurations after registration. First of all, it relies on the determinant of the Jacobian matrix of the mapping between parent and actual element configuration F , to detect the nodes that make the elements invalid. Let n be one of these nodes and {e i } i=1,...,k the set of elements for which n is invalid. The direction of displacement of n is deduced from the gradient of F computed from each e i . In other words, the direction of displacement of n is the one that maximizes the sum of the invalid Jacobians associated to node n. A "displacement step" is also computed for each node and finally the overall set of invalid nodes is iteratively shifted until the mesh is repaired.

Although the approach of Luboz was successfully used to achieve a valid state of several registered meshes, it didn't consider the JR to improve the quality of the elements. Another important remark is that Luboz computed the direction of displacement of invalid nodes n by considering only elements at which n was invalid, and it should be noted here that if two elements e 1 and e 2 share node n and e 1 is valid but e 2 is not, the displacement of n to improve e 2 might cause invalidity at e 1 .

In this article, we propose a new mesh repair technique that only considers the displacement of nodes that make the mesh invalid or of poor quality. To compute the new position of a given node, all the elements linked to it are taken into account. Moreover, cases where invalid configurations arise from the computation of new nodal positions are also managed. Once mesh validity is achieved, the overall quality is assessed and if necessary, our algorithm reallocates mesh nodes exhibiting poor quality. Here, the chosen quality indicator is the JR [START_REF] Kwok | A simple and effective mesh quality metric for hexahedral and wedge elements[END_REF], for which we have chosen a threshold value JR min = 0.03, as mentioned earlier in subsection 2.2. The details of our mesh repair technique are given in the next section.

The proposed mesh repair procedure

The mesh repair procedure proposed in this article is a twofold sequential process:

first all the elements in the mesh are inspected and, if necessary, their nodes positions are adapted so as to recover validity, -then a second procedure is carried out on the mesh if the achieved quality levels are unacceptable.

The corrections applied to the nodal positions are based on the validity and quality criteria, both derived from the Jacobian of the mapping between the parent and actual element frame, as described above. Section 4.1 gives details on the evaluation of the Jacobian for a given element.

Usually distorted elements are sparsely distributed within the mesh, therefore both validity and quality optimizations can be carried out on "repair regions" defined as connected subsets of the mesh encapsulating distorted elements where displacements may only be applied to "improper nodes", i.e. nodes being either invalid or presenting poor quality, leaving other regions unaffected. This local repair strategy significantly decreases computational complexity as the number of degrees of freedom to be considered in a region is usually small. The local relaxations strategies for both validity recovery and quality improvement are described in section 4.2.

The above mentioned relaxations are carried out as maximizations of validity and quality "energy" terms. The formulation of each energy term ensures that the mesh nodal configuration reached at the end of the maximization procedure meets the desired validity or quality standards. The formulations of both energy terms are discussed in section 4.3.

The complete implementation of the iterative relaxation algorithm, common to both validity and quality optimization phases, is given in section 4.4.

Finally, section 4.5 gives a complete registration and repair example from which a particular repaired region is explained in detail.

Jacobian computation

Parent-to-actual element mapping

Mesh relaxation relies on the iterative computation of the Jacobian matrix of the mapping between the element parent (or reference) coordinates system (ξ 1 , ξ 2 , ξ 3 ) and the actual element coordinates (x 1 , x 2 , x 3 ): F : R 3 → R 3 , ξ → x (see Fig. 2). The Jacobian matrix is computed at a number of "control points" which can be either mesh nodes or element integration points (Gaussian quadrature points).

For a given element type and mesh configuration the Jacobian matrix ∂F/∂ξ can be assembled and the exact value of its determinant quickly evaluated using pre-computed element shape function derivatives w.r.t. parent coordinates, as described below.

Let ξ = (ξ 1 , ξ 2 , ξ 3 ) be an element point in parent coordinates and x = (x 1 , x 2 , x 3 ) the same point in actual element coordinates. For an isoparametric element the mapping F can be expressed as a combination of the shape functions φ i : R 3 → R used for value interpolation of the PDE solution within the element:

x(ξ) := F (ξ) = N X n=1 φ n (ξ) x n ( 1
)
where N is the number of nodes of the considered element and x n = (x n 1 , x n 2 , x n 3 ) are the coordinates of the n th element node.

The Jacobian matrix J is defined at parent frame point ξ as:

J(ξ) := ∂x ∂ξ (ξ) = 0 @ J 11 J 12 J 13 J 21 J 22 J 23 J 31 J 32 J 33 1 A
where J ij = (∂x i /∂ξ j )(ξ).

The Jacobian matrix elements can be further rewritten using Eq. 1, as:

J ij = ∂ ∂ξ j N X n=1 φ n (ξ) x n i = N X n=1 ∂φ n ∂ξ j (ξ) x n i = N X n=1 P n j (ξ) x n i
In the above expression all P n j (ξ) := (∂φ n /∂ξ j )(ξ) are independent of the current element configuration {x n } n=1,...,N and can be evaluated beforehand for all points of interest ξ.

Finally, the Jacobian matrix of mapping F computed at point ξ is assembled using the 3 • N pre-computed shape function derivatives P n j as follows:

J(ξ) = N X n=1 0 @ x n 1 x n 2 x n 3 1 A `P n 1 (ξ) P n 2 (ξ) P n 3 (ξ) ´(2)
From the computational point of view, the assembly of the matrix J(ξ) requires 9 • N multiplications and 9 • N additions, and the computation of the 3 × 3 determinant |J(ξ)|, 9 multiplications and 5 additions.

If fast computation is mandatory, the technique proposed by Knupp in [START_REF] Knupp | Achieving finite element mesh quality via optimization of the jacobian matrix norm and associated quantities. part ii -a framework for volume mesh optimization and the condition number of the jacobian matrix[END_REF] and based on a geometrical decomposition of hexahedra and wedges into tetrahedra, can be considered.

Normalized Jacobian value

The Jacobian of the mapping F at parent point ξ represents the ratio between the measure of an infinitesimal volume located around ξ and its image by F in the element actual configuration, as suggested by the notation ∂x/∂ξ.

At element definition stage, parent coordinates magnitudes are usually not related to the physical world of the modeled domain. If the modeled domain has dimensions the order of the millimeter and the mesh nodes coordinates are expressed in meters, then these coordinates are the order of 0.001. The parent coordinates, on the other hand, are usually defined in a standard frame and their value range is in the order of 1.

As a consequence, if we consider a perfectly regular element, the Jacobian computed at one of the nodes already carries the information of the scale change in all 3 dimensions. Its value is thus 0.001 • 0.001 • 0.001 = 10 -9 . This small value, although mathematically correct, does not numerically reflect the evenness of the element shape. In this case, the expected Jacobian value, ignoring the scale change, would be 1.

In order to measure element distortion while eliminating the scale factor between parent and actual element referential, we use a "normalized" value of the Jacobian defined as the product of the original Jacobian value defined above and a correction term K computed as follows.

Let's consider the trihedral defined by the Jacobian under study. Let x be the mean edge length of the 3 edges forming the trihedral in the actual element frame. Now let ξ be the mean edge length of the corresponding trihedral in parent configuration. We define the correction factor in 3D as:

K = " ξ x « 3
In 2D the above fraction should be raised to the power of 2. The normalized value of the Jacobian computed at parent point ξ is now defined as:

J(ξ) := K J(ξ)
If we consider the initial example where a scale factor of 0.001 was applied to all 3 coordinates between parent and actual element frame, the normalized value of the Jacobian for all element points is now J = 1. This fact can be formally apprehended by considering the following rewriting of the normalized Jacobian and noting that the 3D determinant raises the contents of the Jacobian matrix to the power of 3:

J(ξ) = " ξ x « 3 ˛∂x ∂ξ ˛= 1
To simplify subsequent computations, we will assume that both mean edge lengths ξ and x computed beforehand using the input mesh configuration, do not significantly vary during validity or quality relaxations. The K scale factor computed for a given Jacobian is thus considered constant.

In the remainder of the article, for the sake of clarity, we will omit the bar placed over the J symbol, as in the above notation, and refer to the normalized Jacobian simply as "Jacobian".

Jacobian differentiation

As the relaxation procedure relies on variational computation of Jacobian or JR values w.r.t. actual element nodes coordinates, it is interesting to derive from the formulation given in §4.1.1 the required differential expression. To this end, let's consider a Jacobian value estimated at parent point ξ within an element built on the set of nodes {x n } n=1,...,N .

The Jacobian matrix depends now on the element nodes positions, hence J = J(x 1 , . . . , x N , ξ), and we shall seek the expression of the partial derivative of the determinant ˛J(x 1 , . . . , x N , ξ) ˛w.r.t. the i th coordinate of the n th node, i.e.:

∂ ˛J(x 1 , . . . , x N , ξ) ˛/∂x n i
For clarity, we denote this quantity as ∂J n i (ξ). Thus:

∂J n i (ξ) := ∂ ∂x n i ˛J(x 1 , . . . , x N , ξ) ˛= ∂ ∂x n i ˛N X n=1 0 @ x n 1 x n 2 x n 3 1 A `P n 1 (ξ) P n 2 (ξ) P n 3 (ξ)
´T he analytical expression of the above derivative can be easily obtained for i = 1, 2, 3. Let's start by deriving it for i = 1. Expression 2 can be rewritten as:

J(x 1 , . . . , x N , ξ) = 0 B B B B @ P n x n 1 P n 1 (ξ) P n x n 1 P n 2 (ξ) P n x n 1 P n 3 (ξ) P n x n 2 P n 1 (ξ) P n x n 2 P n 2 (ξ) P n x n 2 P n 3 (ξ) P n x n 3 P n 1 (ξ) P n x n 3 P n 2 (ξ) P n x n 3 P n 3 (ξ) 1 C C C C A (3)
Now let's denote M ij the minor of matrix J after removal of row i and column j. For example:

M 23 := ˛J11 J 12 J 31 J 32 ˛= J 11 J 32 -J 31 J 12
Using the minor notation, the determinant of the matrix in Eq. 3 can be developed w.r.t. the first row as:

˛J(x 1 , . . . , x N , ξ) ˛= M 11 X n x n 1 P n 1 (ξ) -M 12 X n x n 1 P n 2 (ξ) + M 13 X n x n 1 P n 3 (ξ)
After deriving the above expression w.r.t.

x n 1 , it follows that:

∂J n 1 (ξ) = M 11 P n 1 (ξ) -M 12 P n 2 (ξ) + M 13 P n 3 (ξ)
Finally, the general expression of ∂J n i (ξ) can be obtained by developing the determinant of the Jacobian matrix w.r.t. the i th row:

∂J n i (ξ) = (-1) i+1 ˆMi1 P n 1 (ξ) -M i2 P n 2 (ξ) + M i3 P n 3 (ξ) 4.2

Local relaxations

The goal of the relaxation procedures is to successively optimize the Jacobian and JR values, during validity recovery and quality optimization phases respectively, by moving the "improper nodes" (INs), i.e. the mesh nodes which Jacobian or JR values need improvement. As the INs are being moved by the relaxation procedure, the Jacobian or JR values at the neighboring nodes can also be altered and potentially degraded.

Given the large number of nodes often found in FE meshes, repairing procedures can quickly become computationally prohibitive. This complexity can be greatly reduced by grouping all INs and their neighborhood into "repair regions" defined in such a way that all nodal corrections applied inside a region leave the outside mesh configuration unchanged. This local repair approach makes it possible to perform all relaxation procedures independently on each repair region identified within the mesh. The construction of such regions is described below.

The topological definition of repair region is different between validity and quality optimization phases. This difference stems from the criteria being optimized: Jacobian or JR. During validity recovery, the displacement of a mesh node n affects the Jacobians of all nodes connected by an edge to n. While during quality optimization phase, the Jacobian at n can become the maximal Jacobian of any element n is part of, and thus the JR of any node sharing an element with n can be affected by its position shift.

Prior to mesh validity recovery, the repair regions are extracted as follows. All mesh nodes are examined. If a node n has a negative or null Jacobian within an element and is not already part of a region then a new region is grown starting from n.

All Jacobians influenced by the newly found improper node n are listed. For each of these Jacobians, the 4 nodes forming the Jacobian trihedral are evaluated. If one of these nodes appears to be invalid, it is assigned to the current region and all its related Jacobian trihedrals within the mesh are in turn examined. The repair region growth stops where all nodes under scrutiny have strictly positive Jacobian values within all elements they are part of.

The repair regions have the following properties. An IN is part of one and only one repair region. The same being true of Jacobian trihedrals, the regions are computationally independent. A repair region can comprise many INs along with their associated Jacobian trihedrals.

The repair region extraction performed before quality optimization only slightly differs from the above algorithm. All mesh nodes are examined. If a node n has a JR below the acceptance threshold JR min , and is not already part of a region then a new region is grown starting from the newly found improper node n.

All JR influenced by n are listed and for each one, the nodes forming the associated element are evaluated. If one of these nodes appears to have a poor quality level, it is assigned to the current region and all its related JR are in turn examined. The region growth stops where all nodes under scrutiny have satisfactory JR values within all elements they are part of. As before, the resulting repaired regions are computationally independent.

During mesh relaxation phases, it is not uncommon to find a repair region having a nodal configuration for which a solution cannot be found by only moving its INs. In such case, the number of degrees of freedom of the region needs to be extended by allowing other nodes to be displaced by the relaxation procedure.

In our mesh repairing algorithm, we have implemented a straightforward region extension method. A locking situation is detected in a repair region when the attractor of the region relaxation is a nodal configuration having persistent INs. Should this situation arise, the region degrees of freedom are extended by allowing the displacement of a mesh node connected to a region IN by an edge, during mesh validity recovery phase, or by an element, during quality optimization phase. This new node is tagged as "improper" although its associated Jacobian or JR values are acceptable, as it was not identified as IN in the first place. Validity or quality criteria are used to choose the new IN among all possible neighbors to the region INs. During mesh validity recovery or quality optimization, the node having the lowest Jacobian or JR respectively, is selected.

As the extended region has a wider influence on the surrounding mesh configuration and might overlap with another region, a "fusion" check is performed after a region has been extended. To this end the sets of Jacobians influenced by the INs of all regions are checked against each other. If a common Jacobian is found between two regions, these regions are merged into a single one. This way, the resulting repaired regions layout retains the initial numerical independence that allows all regions to be processed separately and efficiently.

-1.5 1.5 -2 2 ϕ(t) t k = 0.5 (a) -1.5 1.5 -2 2 ϕ(t) t k = 1 (b) -1.5 1.5 -2 2 ϕ(t) t k = 2 (c)
Fig. 4 Penalty function ϕ k plot, at the vicinity of 0, for k=0.5 (a), k=1 (b) and k=2 (c).

Validity and quality energy formulation

The goal of validity recovery for a mesh region R is to find a configuration where all the Jacobians2 in R, {J j } j∈R , are positive. This nodal relaxation can be formulated as a maximization procedure driven by a "validity energy" E V , expressed as the sum of all Jacobians within R.

As the maximization of a sum of terms does not guarantee that all terms end up being positive, all the Jacobian values within the sum are affected by penalty function ϕ k , which strength is controlled by an index k, giving the expression of validity energy:

E V = X j∈R ϕ k (J j ) , where ϕ k (t) = 1 -e -kt
As the parameter k is increased during the maximization process, the influence of negative values overbalances the positive ones thus favoring a solution where all Jacobians in the sum are positive.

Fig. 4 shows the plot of ϕ k for three different values of k. As can be seen in the figure, ϕ k (t) ≤ 1 for all values of t, the slope at 0, which is the acceptance threshold for the region Jacobians, is k and the function decreases rapidly as t reaches lower negative values.

The aim of the quality optimization procedure, on the other hand, is to find a nodal configuration where all the JR in a region R, {JR j } j∈R , are above a predefined value JR min which, in our implementation, was set to 0.03 in accordance with the quality standard requested by the commercial FE analysis software ANSYS Workbench (ANSYS Inc., USA) [START_REF] Kelly | Element shape testing[END_REF].

We define the "quality energy" E Q associated to this relaxation problem as the sum of the JR within R weighted by a different penalty function ψ l , thus:

E Q = X j∈R ψ l (JR j ) , where ψ l (t) = 1 -e l(JRmin-t)
Similarly to mesh validity recovery, a penalty parameter l is used to find a solution where all the JR contributing to the sum E Q are above the quality threshold JR min . The penalty functions ψ l behave similarly to ϕ k : all values are lower than 1, the slope at the JR acceptance threshold JR min is l and all values below this limit are increasingly penalized as the value of l is raised during relaxation.

Prior to region relaxation at validity recovery or quality optimization stage, the initial values of the parameter for the corresponding penalty function ϕ k or ψ l must be carefully chosen. Indeed, given their formulation, an excessive penalization of unwanted Jacobian or JR values might induce strong numerical instabilities in the optimization process. To avoid this issue, the starting value of k or l is determined by computing the slope of the considered penalty function at the most penalized energy term (minimal Jacobian during validity recovery; minimal JR during quality optimization) and ensuring that it does not exceed a predefined slope threshold.

As the optimization advances, the penalty level k or l is steadily risen to force all energy term values to converge towards the desired interval: ]0, 1] during validity recovery and ]JR min , 1] during quality optimization phase.

Both optimizations are carried out by gradient ascent, with a limited number of ascent steps, each having a limited amplitude. These restrictions ensure that the nodal configuration of the relaxed mesh retains a certain similarity with the input mesh. The detailed relaxation algorithm is discussed in §4.4.

Validity and quality energy maximization

Both relaxation procedures, i.e. the maximizations of the energies E V and E Q defined above, are performed for each region by gradient ascent. The optimal nodal correction directions are given by the gradient of the energy with respect to the coordinates of the INs in the considered region. These vectors are computed at each step using the derivatives of the penalty functions dϕ k /dt or dψ l /dt along with the expression of the Jacobian differential given in §4.1.3.

The procedure successfully stops once an acceptable nodal configuration has been found. Otherwise, if locking configurations are detected, successive region extensions are performed and the relaxation is repeated with a higher number of degrees of freedom in the region under study. The relaxation procedure fails if no suitable nodal configuration can be found.

The mesh repair technique presented here was developed as part of the Mesh-Match-and-Repair [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: Application to 60 clinical cases[END_REF] algorithm which relies on a non-rigid deformation of a generic mesh towards a specific geometry. In this context, the quality of the domain surface representation yielded by the non-rigid deformation procedure must be preserved as much as possible by the mesh relaxation.

In the biomedical field, the dimensions of the modeled organs are the order of the centimeter, and usually a sub-millimetric mean surface representation accuracy is required to ensure the relevance of the Finite Elements analysis carried upon. To control the alteration of surface representation produced by non-rigid registration, our implementation of mesh repair limits the number of relaxation steps to 50, and the maximal amplitude of nodal corrections at each step to 0.1 mm. The maximal allowed nodal displacement is thus 5 mm, a situation which seldom occurs as it requires the same node to be moved in a constant direction throughout the relaxation process.

The common structure for both relaxation procedures is given by Alg. 1 and the main steps (identified by a number) are commented below.

Critical steps within the above relaxation algorithm are commented below.

- [START_REF] Alliez | Variational tetrahedral meshing[END_REF] : To ensure numerical stability of the relaxation, find the maximal value of k or l, depending on the current relaxation phase, which complies with the penalty function slope limit, as discussed in §4.3.

the iterations stop when all mesh nodes connected to the initial improper nodes in the considered region have been included.

Repairing example

Fig. 5 shows an example of mesh registration and repair for a face model. At the end of the repair process, all the elements were valid and presented acceptable quality (JR > 0.03). In particular, Fig. 6 shows a sample validity recovery for a repaired region. The invalid Jacobian trihedrals are drawn in dark and the valid Jacobians in light gray-scale. The region repair depicted in Fig. 6 took approximately 600 ms. The region had initially 5 invalid nodes and comprised 24 elements. In order to find a suitable nodal configuration the region had to be extended by 2 nodes, therefore including 8 more elements in the region. The mean magnitude of the corrections applied to the 7 displaced mesh nodes was 0.5 mm (max = 0.8 mm, σ = 0.1 mm).

More examples like this can be found in [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: Application to 60 clinical cases[END_REF] where 10 femur and 50 face FE meshes were successfully registered and repaired in the context of orthopaedic and orthognathic hard and soft tissues modelling studies.

Conclusions

A novel mesh repair algorithm, based on a mesh relaxation driven by specific validity and quality energy formulations, has been presented. It first recovers mesh validity and further improves mesh elements quality in order to enhance numerical results reliability. Both procedures rely on the computation of Jacobian values within the mesh elements, which is a mesh acceptance criterion widely used in popular FE analysis softwares. Within both validity and quality energy formulations, the Jacobian measures are weighted by carefully designed penalty functions that lead the nodal corrections towards a configuration where both mesh validity and quality criteria are satisfied. The procedure proposed here was successfully used to produce FE meshes suitable for FE analysis within commercial softwares [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: Application to 60 clinical cases[END_REF][START_REF] Bucki | Finite element speaker-specific face model generation for the study of speech production[END_REF], such as ANSYS (Ansys Inc., USA).

The minimized energy function presented in this work is not smooth. To approximate the descent direction despite the slope discontinuities a Centered Finite Differences (CFD) scheme was used. From a numerical point of view, the CFDs act as a smoothing operator over the non-smooth energy landscape yielding a good estimate of the gradient in the smooth domain and managing continuous transitions through the slope discontinuity boundaries.

Our mesh repair procedure only affects invalid and poor quality regions, leaving other mesh areas unchanged. If an acceptable configuration cannot be found by applying corrections to invalid or poor quality nodes, other nodes are progressively included into the repair process, according to their validity and quality levels, following the region extension procedure described above, in §4.2. This strategy preserves valid mesh regions as well as those with higher quality level. In addition to the initially identified improper nodes, only nearly invalid or low quality nodes can be reallocated by our procedure which strongly limits the difference between the repaired and initial input mesh.

The region-based mesh untangling process proposed in the article is able to deal with Jacobian inter-locking situations where a solution can only be found by moving simultaneously a number of nodes involved in shared elements. These situations require the whole set of involved elements to be considered simultaneously and, from our experience, can hardly be solved using a node-based approach as antagonist nodal displacements can quickly lead to an infinite loop in the mesh untangling routine. Moreover, several mesh untangling techniques proposed in the literature do not separate validity recovery and quality improvement procedures, intending to achieve it all at the same time. The problem with this approach is that if invalid node a and poor quality node b are connected, finding a new position to improve the quality of b may decreases the chances of finding a new valid position for a. We consider that to focus on mesh validity first, increases the chances of finding a valid position for a. Once a and b are of poor quality, both will have the same chances to increase their quality, but the difference is that both start from a valid mesh.

To our knowledge, the only other strategy to perform such a task was introduced by Luboz et al. [START_REF] Luboz | Automatic 3d finite element mesh generation: an atlas fitted onto patient data[END_REF]. The differences with this work and the main contributions of our approach are summarized bellow.

-In order to compute the new position of an invalid node n, only the elements that were invalid at n were considered in the work of Luboz et al. Our technique, on the other hand, takes into consideration all the elements attached to n as its displacement can cause previously valid elements to become invalid. -In our approach, the nodes to be repaired are grouped in independent regions to avoid prohibitive computational costs. -As opposed to Luboz et al., our algorithm considers mesh quality as well as mesh validity indicators. It is thus capable of producing a mesh that guarantees the stability of the numerical solution yielded by the FE analysis carried out upon.

Perspectives

In the context of Mesh-Matching, repair strategies should take into account the surface representation constraint, in situations where surface nodes are being moved during relaxation. This can be accomplished by considering the gradient of the energy undergoing maximization at a given improper node not along all space directions, as described above, but tangentially to the represented organ surface. The surface representation constraint should only be implemented for surface mesh nodes, i.e. the ones that need to rest as close as possible to the modeled domain surface during the relaxation, to ensure the accuracy of geometry representation. Domain surface representation could be further maintained by correcting small registration degradations likely to occur during repairing. To this end, an Atlas-tosurface registration step might be performed after each relaxation step. Nevertheless, as registration and relaxation are independent and not necessarily collaborative procedures, the convergence of this alternating method should be carefully studied so as to avoid infinite loops and eventually reach a repaired mesh configuration.

The local repair technique presented in §4.2 has proved to be effective in all cases considered so far [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: Application to 60 clinical cases[END_REF]. Yet, nodal configurations featuring intricate element inter-locking can be devised so that the approach proposed here fails. To overcome this issue, alternate relaxation paradigms can be considered. Indeed, we have often observed that during relaxation complex repair regions transit through a configuration with a strongly reduced number of improper nodes. In cases where the relaxation fails, instead of extending the original repaired region and repeating the process from its initial configuration, a more efficient approach would be to start from the region configuration featuring the minimal number of INs encountered during the unsuccessful relaxation.

Should a repair region be restricted to this smaller subset, its influence within the mesh must be recomputed and possible region overlaps checked for. Furthermore, in order to control the nodal correction amplitudes, track must be kept of anterior nodal displacements in order to guarantee that the original surface representation does not get excessively perturbed. Although promising, this approach has not yet been implemented in our software.

Finally, future work should be focused on adding other quality criteria to repair the mesh. The JR is an excellent quality criterion that helps to detect invalid and poor quality elements, but it cannot detect, for instance, elements with small dihedral angles or elements with warped faces. For these reasons, other quality criteria like dihedral angles, aspect-ratio, volume, warping factor, etc. should be also considered to ensure a quality output mesh.

Fig. 1

 1 Fig.1The mesh-matching algorithm: (a) The Atlas mesh, (b) the target domain or patient data, (c) the Atlas and the target surface domain represented by a dashed line and (d) the final mesh of the target organ obtained by applying the transformation function T onto the entire Atlas mesh. The circled region shows the location of possible distorted element leading to mesh invalidity or poor quality.

  Fig.3If the signed distance of node a to the plane formed by "connected nodes" b, c and d is positive, the element is valid (left), in any other case the entire element is said to be invalid (right).

Fig. 5 Fig. 6

 56 Fig. 5 Left: input patient data (skull and skin surfaces segmented from medical images). Middle: patient data and Atlas mesh. Right: the final mesh registered onto the patient's skin and skull and repaired.

See http://cubit.sandia.gov for details.

Note that in the remainder of the article, all considered Jacobian values are assumed to be normalized as described in 4.1.2

Acknowledgements This project has been financially supported by FONDECYT 1061227, ALFA IPECA project, STIC-AmSud PLOMO, and ECOS-Sud C06E04.

Algorithm 1 Relaxation algorithm for mesh validity recovery and quality optimization

Identify repair regions within mesh : Rr , r = 1, . . . , R for r = 1 to R do {Regions loop} repeat {Region extension loop} Find initial penalty parameter for region Rr : k or l [START_REF] Alliez | Variational tetrahedral meshing[END_REF] for penalty = 1 to 5 do {Penalty loop} for step = 1 to 10 do {Energy maximization loop} Compute energy value and gradient at current region configuration Find optimal nodal correction in gradient direction within 0. - [START_REF] Berar | 3d meshes registration : Application to statistical skull model[END_REF] : The considered relaxation energy gradient vector being formed by the derivatives of the energy w.r.t. each 3D coordinate of the region improper nodes, the coordinates of each nodal correction appear in the gradient vector on the corresponding row. In our specific implementation, the magnitudes of the nodal corrections are clamped, at each step, to 0.1 mm by applying a scale factor to the overall gradient vector. The value of the maximal nodal correction amplitude used here is not universal and must be determined for each application field according to the dimensions of the modeled domain and maximal tolerance on the representation of its geometry. The local maximum line search in the gradient direction is performed using the golden section technique [START_REF] Press | Numerical Recipes in C[END_REF] between the current nodal configuration and the configuration obtained after applying the maximal allowed amplitude correction. - [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: Application to 60 clinical cases[END_REF] : Raise the value of k or l so as the slope at the most penalized term reaches its maximal allowed value. - [START_REF] Bucki | Finite element speaker-specific face model generation for the study of speech production[END_REF] : Increase the region number of degrees of freedom by including new nodes, edge or element-connected to the region improper nodes, as described in §4.2. - [START_REF] Castellano-Smith | Constructing patient specific models for correcting intraoperative brain deformation[END_REF] : The region extension can be carried out to the point where the region overlaps the entire mesh and hence all other regions. Yet it appears unnecessary to put into practice this exhaustive and computationally prohibitive solution as in most cases, acceptable region configurations can be found by including in the relaxation process the immediate neighborhood of the originally improper region nodes. Thus,