David Vicente

Anisotropic Minkowski Content and Application to minimizers of free discontinuity problems

come

Anisotropic Minkowski Content and Application to minimizers of free discontinuity problems

Introduction

For S ⊂ R n , the (n -1)-dimensional upper and lower Minkowski contents of S are defined by

M (S) = lim sup ρ→0 + L n ({x : dist(x, S) < ρ}) 2ρ , M (S) = lim inf ρ→0 + L n ({x : dist(x, S) < ρ}) 2ρ ,
where L n is the Lebesgue measure. If M (S) = M (S), their common value, denoted by M(S), is called the (n -1)-dimensional Minkowski content. In [START_REF] Federer | Geometric Measure Theory[END_REF] (Theorem 3.2.39), the following result is given.

Theorem 0.1 (Federer). If S is a closed and (n -1)-rectifiable subset of R n , then

M(S) = H n-1 (S),
where H n-1 is the (n -1)-dimensional Hausdorff measure.

In this article, we focus on a class of anisotropic perturbation of H n-1 defined on (n -1)rectifiable surfaces S ⊂ R n by S Mν, ν 1/2 dH n-1 , (0.1)

where M : R n → S + n (R) is a continuous field of symmetric definite positive matrices, ν is an unitary vector orthogonal to the tangent plane of S. As S is rectifiable, then ν is well defined H n-1 -almost everywhere in S. In order to get an approximation of (0.1) which involves the Lebesgue measure, as in Theorem 0.1, for any (x, v) ∈ R n × R n , we consider the following Riemannian metric:

φ(x, v) = M -1 (x)v, v 1/2
and for any (x, y) ∈ (R n ) 2 its associated integrated distance: L n ({x : dist φ (x, S) < ρ}) 2ρ .

dist φ (x, y) = inf
If M M (S) = M M (S), we call their common value the (n -1)-dimensional anisotropic Minkowski content M M (S). In this paper, we prove the following Theorem 0.2. If S is a closed and (n -1)-rectifiable subset of R n and M : R n → S + n (R) is continuous, then we have

M M (S) = S Mν, ν 1/2 dH n-1 ,
where ν is an unitary vector, normal to S.

A comparable result is given in [START_REF] Bellettini | Some results on surface measures in calculus of variations[END_REF] (Theorem 6.1) for S = ∂E and E is a set of finite perimeter. In [START_REF] Chambolle | A remark on the anisotropic Minkowski content[END_REF], the study focuses on anisotropic outer Minkowski content for the same class of sets. In our case, S is not necessary the boundary of a set. This generalization is of interest for the case where S is the jump set J u of a function u with bounded variation as we will see in last section.

Our result allows us to extend a regularity result for free boundary problems to a larger class of energies which involve a continuous perturbation of the Hausdorff measure. More precisely, as it has been introduced in [START_REF] Bucur | Monocity formula and regularity for general free discontinuity problems[END_REF], u ∈ SBV(Ω) is an almost quasi-minimizer of a free discontinuity problem if there exist σ ≥ 1, α > 0 and

c α > 0 such that: if v ∈ SBV(Ω), B(x, r) ⊂ Ω and [w = v] ⊂ B(x, r), then we have B(x,r) |∇w| 2 dx + H n-1 (J w ∩ B(x, r)) ≤ B(x,r) |∇v| 2 dx + σH n-1 (J v ∩ B(x, r)) + c α r n-1+α ,
where J u is the jump set of u. The difference of this definition with respect to the classical notion of quasi-minimizer consists in the coefficient σ ≥ 1 which allows bounded perturbation of H n-1 . We prove the following Theorem 0.3. Let M : R n → S + n (R) be continuous and u ∈ SBV(Ω) be an almost quasi-minimizer of a free discontinuity problem, then we have

M M (J u ) = Ju Mν u , ν u 1/2 dH n-1 .
This result is applied to the Mumford-Shah model in an anisotropic setting. Indeed, a minimizer this model is an almost quasi-minimizer of a free boundary problem and the approximation result of Theorem 3.1 is satisfied. In a joint paper [START_REF] Vicente | Anisotropic Mumford-Shah model and its approximation with Γ-convergence[END_REF], this is the key tool to perform a Γ-convergence approximation of our model.

In section 1, we recall some classical results of Geometric Measure Theory. Section 2 is devoted to the proof of Theorem 0.2. In section 3, we give the proof of Theorem 0.3 and we apply this result to give a rigorous setting for the approximation of an anisotropic Mumford-Shah model.

Geometric Measure Theory Framework

We adopt the notations:

• Ω for an open and bounded subset of

R n , • v 1 , v 2 ∈ R for the canonical scalar product of v 1 , v 2 ∈ R n , • n-1 i=1 v i ∈ R n for the canonical vectorial product of v 1 , . . . , v n-1 ∈ R n , • |v| for the euclidean norm of v ∈ R n ,
• B(x, r) for the open ball in R n with radius r and center x, The geometric measure theory framework is mainly extracted from [START_REF] Federer | Geometric Measure Theory[END_REF] and [START_REF] Bellettini | Some results on surface measures in calculus of variations[END_REF]. Definition 1.1 (Federer,3.2.1). Let f maps a subset of R n-1 onto R n , the (n -1)-dimensional Jacobian is defined by

• M n (R) for the space of n × n real matrices, • M for an induced norm of M ∈ M n (R), • S + n (R) ⊂ M n (R)
J n-1 (f )(a) = n-1 i=1 ∂f ∂x i (a)
whenever f is differentiable at a.

Theorem 1.1 (Area Formula, Federer, 3.2.3). Suppose f :

R n-1 → R n is Lipshitzian. If u is an L n-1 -integrable function, then R n-1 u(x)J n-1 (f )(x)dL n-1 (x) = R n x∈f -1 {y} u(x)dH n-1 (y).
We are interested in the class of subset of R n which are (n -1)-rectifiable.

Definition 1.2 (Federer, 3.2.14). The set E ⊂ R n is (n -1)-rectifiable if there exists a Lipshitzian function f : R n-1 → R n mapping some bounded subset of R n onto E.
The following result is useful to get an univalent parametrization of a rectifiable set.

Theorem 1.2 (Federer,3.2.4). For every L n-1 measurable set A, there exists a Borel set

B ⊂ A ∩ {x : J n-1 (f )(x) > 0} such that f B is univalent and H n-1 (f (A) \ f (B)) = 0.
We say that ϕ : Ω × R n → R + is a Finsler metric if ϕ is continuous and if there exists λ, Λ > 0 such that, for any (x, v, t)

∈ Ω × R n × R, it satisfies φ(x, tv) = |t|φ(x, v), λ|v| ≤ φ(x, v) ≤ Λ|v|.
We define the integrated distance associated to ϕ as

∀(x, y) ∈ (R n ) 2 , dist ϕ (x, y) = inf 1 0 ϕ γ, dγ dt dt : γ ∈ W 1,1 ([0; 1]; R n ), γ(0) = x, γ(1) = y .
For any E ⊂ R n , we define the associated diameter as diam ϕ (E) = sup{dist ϕ (x, y) : x, y ∈ E} and the associated k-dimensional Hausdorff measure by

H k ϕ (E) = lim δ→0 + inf ω k 2 k i∈I (diam ϕ (A i )) k : E ⊂ i∈I A i , ∀i ∈ I, diam ϕ (A i ) ≤ δ
where ω k is the volume of the unit ball in R k .

Theorem 1.3 (Bellettini). Let ϕ be a Finsler metric, then, for all Borel set E ⊂ R n , we have

H n ϕ (E) = E H n (B(0, 1)) H n (B ϕ (x, 1)) dx,
where B ϕ (x, 1) is the unit ball centered at x for the metric ϕ.

2 Proof of Theorem 0.2

The proof is divided in two steps. First, we assume in section 2.1 that the metric does not depend on x. Then, in section 2.2, we remove this assumption and we prove the Theorem 0.2 in the general setting.

Homogeneous Case

In this section, we assume that the metric is homogeneous. So, we denote by M 0 its common value and, for any v ∈ R n , we set

φ 0 (v) = M -1 0 v, v 1/2 , ( 2.1)
so that (R n , φ 0 ) is an Euclidean space. We prove the following Theorem 2.1. If S is a closed and (n -1)-rectifiable subset of R n and M 0 ∈ S + n (R). Then we have

M M0 (S) = S M 0 ν, ν 1/2 dH n-1 ,
where ν is an unitary and normal vector to S.

We decompose the proof in three Lemmas. As M -1 0 is symmetric positive definite, R n may be viewed as an euclidean space according to the scalar product M -1 0 •, • . So, Theorem 0.1 may be directly applied. The three following Lemmas consist in the change of variable from R n endowed with the scalar product M -1 0 •, • to R n endowed with the canonical scalar product. Lemma 2.1. For M 0 ∈ S + n (R), φ 0 defined by (2.1) and for any E ∈ B(R n ), we have

L n (E) = det(M 0 )H n φ0 (E).
Proof. As M 0 is a symmetric positive definite matrix, there exists P ∈ O n (R) and (λ i ) i=1...n such that λ i > 0 for any i ∈ {1, . . . , n} and

M 0 = P       λ 1 0 • • • 0 0 λ 2 . . . . . . . . . . . . . . . 0 0 • • • 0 λ n       P -1 . ( 2.2)
We denote by (c i ) i=1...n the canonical basis of R n , then ( √ λ i P c i ) i=1...n is an orthonormal basis for the scalar product M -1 0 •, • . So, the linear application L : R n → R n characterized by L(c i ) = √ λ i P c i for any i ∈ {1, . . . , n}, is an isomorphism which satisfies L(B(0, 1)) = B φ0 (0, 1) and det(L) = det(M 0 ). It gives

L n (B φ0 (0, 1)) = B(0,1) |det(L)|dx, = det(M 0 )L n (B(0, 1))
According to isodiametric inequality (see 2.10.35 in [START_REF] Federer | Geometric Measure Theory[END_REF]), we have H n = L n and then, for any

E ∈ B(R n ), Theorem 1.3 yields H n φ0 (E) = L n (E) det(M 0 ) .
which concludes the proof of Lemma 2.1.

Setting

e i = √ λ i P c i , then (e i ) i=1,...,n is an orthonormal basis of R n for the scalar product M -1 0 •, • . In the Euclidean space R n , φ 0 ), the vectorial product of (v i ) i=1,...,n-1 , denoted by n-1 φ0,i=1 v i , is characterized by the following equality M -1 0 v 1 , e 1 . . . M -1 0 v n-1 , e 1 M -1 0 w, e 1 . . . . . . . . . . . . . . . . . . M -1 0 v 1 , e n . . . M -1 0 v n-1 , e n M -1 0 w, e n = M -1 0 n-1 φ0,i=1 v i , w , (2.3)
which must be true for any w ∈ R n .

Lemma 2.2. If M 0 ∈ S + n (R) and φ 0 defined by (2.1), then we have

n-1 φ0,i=1 v i = M 0 (∧ n-1 i=1 v i ) det(M 0 )
.

Proof. As P ∈ O n (R), then (P c i ) n i=1 is an orthonormal basis for the usual scalar product and, for any w ∈ R n , we have

v 1 , P c 1 . . . v n-1 , P c 1 w, P c 1 . . . . . . . . . . . . . . . . . . v 1 , P c n . . . v n-1 , P c n w, P c n = ∧ n-1 i=1 v i , w .
According to (2.2), for any v ∈ R n , we have

v, P c i = λ i M -1 0 v, λ i P c i and then λ 1 . . . λ n M -1 0 v 1 , e 1 . . . M -1 0 v n-1 , e 1 M -1 0 w, e 1 . . . . . . . . . . . . . . . . . . M -1 0 v 1 , e n . . . M -1 0 v n-1 , e n M -1 0 w, e n = M -1 0 M 0 ∧ n-1 i=1 v i , w .
According to (2.3), it concludes the proof of Lemma 2.2.

As for the vectorial product, we may define an anisotropic (n -1)-dimensional Jacobian by

J φ0 n-1 (f )(a) = n-1 φ0,i=1 ∂f ∂x i (a)
whenever f is differentiable at a. The following result is a straightforward consequence of Lemma 2.2.

Lemma 2.3. Let M 0 ∈ S + n (R) be fixed and φ 0 defined by (2.1), then we have

J φ0 n-1 (f )(a) = 1 det(M 0 ) M 0 n-1 i=1 ∂f ∂x i (a) , n-1 i=1 ∂f ∂x i (a) 1/2
. Now we give the proof of Theorem 2.1.

Proof. Let S be a (n -1)-rectifiable closed set. As L, introduced in the proof of Lemma 2.1, is an automorphism, then S is also (n -1)-rectifiable and closed for R n endowed with the scalar product M -1 0 •, • . According to Lemma 2.1, we have

L n ({x : dist φ0 (x, S) < ρ}) 2ρ = det(M 0 ) H n φ0 ({x : dist φ0 (x, S) < ρ}) 2ρ . (2.4)
We may apply Theorem 0.1 in the euclidean space (R n ,

M -1 0 •, • ), it gives lim ρ→0 + H n φ0 ({x : dist φ0 (x, S) < ρ}) 2ρ = H n-1 φ0 (S). (2.5)
As S is rectifiable, there exists a Lipschitzian function f : R n-1 → R n and a bounded subset A ⊂ R n-1 such that S ⊂ f (A). According to Corollary 1.2, there exists a Borel set

B ⊂ A ∩ {x : J φ0 n-1 (f )(x) > 0} such that f B is univalent and H n-1 (f (A) \ f (B)) = 0. We denote by C = f -1 (S) ∩ B and then H n-1 φ0 (f (C)) = H n-1 φ0 (S)
. As f C is univalent, according to Area formula 1.1 with u = 1 C , we have

H n-1 φ0 (S) = C J φ0 n-1 (f )(x)dx. (2.6)
For any x ∈ C we have J φ0 n-1 (f )(x) = 0 which implies that ∂f ∂x1 (x), . . . , ∂f ∂xn-1 (x) is free and then J n-1 (f )(x) = 0. According to Lemma 2.3, we may write

J φ0 n-1 (f )(x) = J φ0 n-1 (f )(x) J n-1 (f )(x) J n-1 (f )(x), = 1 det(M 0 ) M 0   n-1 i=1 ∂f ∂xi (x) n-1 i=1 ∂f ∂xi (x)   , n-1 i=1 ∂f ∂xi (x) n-1 i=1 ∂f ∂xi (x) 1/2 J n-1 (f )(x).
We set

u(x) = M 0   n-1 i=1 ∂f ∂xi (x) n-1 i=1 ∂f ∂xi (x)   , n-1 i=1 ∂f ∂xi (x) n-1 i=1 ∂f ∂xi (x) 1/2 .
According to Lemma 3.2.25 in [START_REF] Federer | Geometric Measure Theory[END_REF], we have

u(x) = M 0 ν(f (x)), ν(f (x)) 1/2 where ν(f (x)
) is an unitary and orthogonal vector to f (C) at f (x). Applying Area formula 1.1, this time in R n endowed with its canonical euclidean structure, gives

C J φ0 n-1 (f )(x)dx = 1 det(M 0 ) C u(x)J n-1 (f )(x)dx, = 1 det(M 0 ) f (C) M 0 ν, ν 1/2 dH n-1 , = 1 det(M 0 ) S M 0 ν, ν 1/2 dH n-1 .
According to (2.4), (2.5) and (2.6) we may conclude

lim ρ→0 + L n ({x : dist φ0 (x, S) < ρ}) 2ρ = S M 0 ν, ν 1/2 dH n-1 .

Inhomogeneous setting

In this section, we remove the homogeneity assumption and, using a piecewise constant approximation of M, we prove Theorem 0.2. The proof is divided in three Propositions.

Proposition 2.1. Let M be as in Theorem 0.2. Then, for any compact neighborhood C of S, there exists λ, Λ > 0, such that the following ellipticity condition is satisfied for any

(x, v) ∈ C × R n λ|v| 2 ≤ M(x)v, v ≤ Λ|v| 2 .
Proof. As M(x) is symmetric positive definite for any x ∈ C, there exists λ(x), Λ(x) > 0 such that, for any v ∈ R n , we have

λ(x)|v| 2 ≤ M(x)v, v ≤ Λ(x)|v| 2 .
As M is continuous, the previous inequalities remain true in any compact C with (λ, Λ) only depending on C.

Proposition 2.2. Let S and M be as in Theorem 0.2. Then, we have

M M (S) ≤ S Mν, ν 1/2 dH n-1 .
Proof. We may assume that S Mν, ν 1/2 dH n-1 is finite, otherwise the result is ensured. Let C be a compact neighborhood of S. Let δ > 0 be fixed. As M is continuous, there exists

η > 0 such that |x -y| ≤ η ⇒ M(x) -M(y) ≤ δ (2.7)
for any (x, y) ∈ C 2 . According to Proposition 2.1, there exists λ > 0 such that

λ 1/2 H n-1 (S) ≤ S Mν, ν 1/2 dH n-1
and then H n-1 (S) is also finite. For t ∈ R and i ∈ {1, . . . , n}, we set

Π i t = {x ∈ C : x, e i = t}. Thus, for k ∈ N fixed, t ∈ R : H n-1 (S ∩ Π i t ) > 1 k
is finite and then, for any i ∈ {1, . . . , n}, the set t ∈ R : H n-1 (S ∩ Π i t ) > 0 is at most countable. So, if we consider the cubes with its faces orthogonal to the vectors of the orthogonal basis, there exists a partition K of C by cubes with diameter less than η, such that for any K ∈ K, we have

H n-1 (S ∩ ∂K) = 0.
(2.8)

For any K ∈ K, we fix a K ∈ K and, for any (x, v) ∈ K × R n , we set

M(x) = M(a K ), φ(x, v) = (M(a K )) -1 v, v .
(2.9)

For K ∈ K and r > 0, we set K r = {x : dist(x, ∂K) ≤ r}. We have the following decomposition

M M (S) ≤ K∈K M M (S ∩ K r ) + K∈K M M (S ∩ (K \ K r )).
(2.10)

In Claim 1 and Claim 2, we will determine an upper bound for the two previous sums.

Claim 1: We have K∈K M M (S ∩ K r ) ≤ λ -n 2 Λ 1/2 K∈K H n-1 (S ∩ K r ).
According to Proposition 2.1, there exists λ > 0 and Λ > 0 such that, for any

(x, v) ∈ C × R n , we have Λ -1 2 |v| ≤ M(x) -1 v, v 1/2 ≤ λ -1 2 |v|.
(2.11)

Let (x, y) ∈ C 2 and γ ∈ W 1,1 ([0; 1]; R n ) such that γ(0) = x and γ(1) = y, then we have

Λ -1 2 1 0 | γ|dt ≤ 1 0 φ(γ, γ)dt ≤ λ -1 2 1 0 | γ|dt.
According to (0.2) and (2.11), for any (x, y) ∈ C 2 , it gives

Λ -1 2 dist(x, y) ≤ dist φ (x, y) ≤ λ -1 2 dist(x, y) (2.12)
and then, for any x ∈ C, the following inclusions are satisfied

B(x, λ 1/2 ) ⊂ B φ (x, 1) ⊂ B(x, Λ 1/2 ).
For any x ∈ C, , as

H n (B(x, t)) = t n H n (B(0, 1)), we get Λ -n 2 ≤ H n (B(0, 1)) H n (B φ (x, 1)) ≤ λ -n 2
and with Theorem 1.3, for any E ∈ B(C), it yields

Λ -n 2 H n (E) ≤ H n φ (E) ≤ λ -n 2 H n (E).
Moreover, (2.12) implies {x : dist φ (x, S ∩ K r ) < ρ} ⊂ {x : dist(x, S ∩ K r ) < Λ 1/2 ρ} and then we have

M M (S ∩ K r ) = lim sup ρ→0 + H n φ ({x : dist φ (x, S ∩ K r ) < ρ}) 2ρ , ≤ λ -n 2 lim sup ρ→0 + H n ({x : dist(x, S ∩ K r ) < Λ 1/2 ρ}) 2ρ , ≤ λ -n 2 Λ 1/2 M (S ∩ K r ).
As S is closed and (n-1)-rectifiable, so is S ∩K r . We may apply Theorem 0.1 to get M (S ∩K r ) = H n-1 (S ∩ K r ) and Claim 1 is proved.

Claim 2: We have

K∈K M M (S ∩ (K \ K r )) ≤ (1 + θ δ) S Mν, ν 1/2 dH n-1 ,
where θ is a constant which depends only on the restriction of M in C.

Let C be a compact neighborhood of S. As M is continuous, then there exists a constant m > 0 which depends only on C such that, for any (x, y) ∈ C 2 , we have

M -1 (x) -M -1 (y) ≤ m M(x) -M(y) .
Let (x, y) ∈ C 2 and γ ∈ W 1,1 ([0; 1]; C) such that γ(0) = x and γ(1) = y, according to (2.7), we have the inequalities

1 0 M -1 (γ) γ, γ 1/2 dt - 1 0 M -1 (γ) γ, γ 1/2 dt ≤ 1 0 | (M -1 (γ) -M -1 (γ)) γ, γ | M -1 (γ) γ, γ 1/2 + M -1 (γ) γ, γ 1/2 dt, ≤ 1 0 M -1 (γ) -M -1 (γ)) γ, γ M -1 (γ) γ, γ 1/2 + M -1 (γ) γ, γ 1/2 dt, ≤ mδ 2Λ -1 1 0 γ, γ 1/2 dt, ≤ mδ 2Λ -3 2 1 0 M -1 (γ) γ, γ 1/2 dt.
We denote m 2Λ -3 2 by θ , it gives

1 0 M -1 (γ) γ, γ 1/2 dt ≤ (1 + θ δ) 1 0 M -1 (γ) γ, γ 1/2 dt, so that {x : dist φ (x, S ∩ (K \ K r )) < ρ} ⊂ x : dist φ (x, S ∩ (K \ K r )) < (1 + θ δ) ρ , (2.13)
and then

M M (S ∩ (K \ K r )) ≤ (1 + θ δ)M M (S ∩ (K \ K r )). ( 2

.14)

As φ is an homogeneous metric in a neighborhood of K \ K r as in Section 2.1, we may apply Theorem 2.1 to obtain

M M (S ∩ (K \ K r )) = S∩(K\Kr) Mν, ν 1/2 dH n-1 and M M (S ∩ (K \ K r )) ≤ S∩K Mν, ν 1/2 dH n-1 .
According to (2.14), we get

K∈K M M (S ∩ (K \ K r )) ≤ (1 + θ δ) S Mν, ν 1/2 dH n-1 ,
so it concludes the proof of Claim 2.

Applying Claim 1 and Claim 2 to the decomposition (2.10), we have

M M (S) ≤ λ -n 2 Λ 1/2 K∈K H n-1 (S ∩ K r ) + (1 + θ δ) S Mν, ν 1/2 dH n-1
Taking the limit r → 0 + gives

M M (S) ≤ λ -n 2 Λ 1/2 K∈K H n-1 (S ∩ ∂K) + (1 + θ δ) S Mν, ν 1/2 dH n-1 .
Applying (2.8), and then taking the limit as δ → 0 + concludes the proof of Proposition 2.2

M M (S) ≤ S Mν, ν 1/2 dH n-1 .
Proposition 2.3. Let S and M as in Theorem 0.2. Then, we have

M M (S) ≥ S Mν, ν 1/2 dH n-1 .
Proof. Let δ > 0, η > 0 and K as in (2.8) and (2.9). For r > 0 and K ∈ K, we still set K r = {x : dist(x, ∂K) ≤ r}. As K is finite and

(K, L) ∈ K 2 , K = L ⇒ dist(K r , L r ) > 0, then we have M M (S) ≥ K∈K M M (S ∩ (K \ K r )) (2.15)
With the same proof as for (2.13), we have

{x : dist φ (x, S ∩ (K \ K r )) < ρ} ⊂ {x : dist φ (x, S ∩ (K \ K r )) < (1 + θ δ)ρ}
and then

L n ({x : dist φ (x, S ∩ (K \ K r )) < ρ}) 2ρ ≤ L n ({x : dist φ (x, S ∩ (K \ K r )) < (1 + θ δ)ρ}) 2ρ .
Passing to the lim inf with ρ → 0 + gives

M M (S ∩ (K \ K r )) ≤ (1 + θ δ)M M (S ∩ (K \ K r )).
As M is constant in a neighborhood of K \ K r as in Section 2.1, we may apply Theorem 2.1 to obtain

M M (S ∩ (K \ K r )) = S∩(K\Kr) Mν, ν 1/2 dH n-1 .
According to (2.15), we get

M M (S) ≥ (1 + θ δ) -1 K∈K S∩(K\Kr) Mν, ν 1/2 dH n-1 , ≥ (1 + θ δ) -1 S∩Cr Mν, ν 1/2 dH n-1 ,
where we have set C r = K∈K (K \ K r ). Remark that

r 1 < r 2 ⇒ C r2 ⊂ C r1 , r>0 C r = C \ K∈K ∂K so, we have lim r→0 + S∩Cr Mν, ν 1/2 dH n-1 = S\∪∂K Mν, ν 1/2 dH n-1 .
As H n-1 (∂K) = 0 for any K ∈ K (2.8), then lim

r→0 + S∩Cr Mν, ν 1/2 dH n-1 = S Mν, ν 1/2 dH n-1 .
We deduce that

M M (S) ≥ (1 + θ δ) -1 S Mν, ν 1/2 dH n-1
and passing to the limit for η → 0 + and then for δ → 0 + concludes the proof.

so, H n-1 (J u ) is finite too.

According to [START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF], Section 5.9, the set J u is rectifiable up to a H n-1 -negligible set N . More precisely, there exists a countable family (K i ) i∈N of compact C 1 -hypersurfaces such that

J u = N ∪ i∈N K i ,
where H n-1 (N ) = 0. Let K ⊂ J u ∩ Ω be compact, we have the decomposition

K = K ∩ J u \ J u ∪ K ∩ q i=1 K i ∪   K ∩ ∞ i=q+1 K i   ∪ [K ∩ N ] . Theorem 3.2 gives H n-1 (K ∩ J u \ J u ) = 0. Let δ > 0 be fixed. As H n-1 (J u ) is finite, there exists q ∈ N such that H n-1   K ∩ ∞ i=q+1 K i   ≤ δ. (3.1)
In the sequel, we omit the dependance with δ for the sake of simplicity. We set S = K ∩ q i=1 K i and for A ⊂ R n we adopt the following notation

A ρ := {x : dist φ (x, A) < ρ}.
Let τ > 0 and ρ > 0 be fixed, we decompose

K = (K \ S τ ρ ) ∪ S τ ρ and set E = K \ S τ ρ . So, we have K ρ ⊂ E ρ ∪ S (1+τ )ρ ,
and then

L n (K ρ ) 2ρ ≤ L n (E ρ ) 2ρ + L n (S (1+τ )ρ ) 2ρ . (3.2)
The rest of the proof consists in computing an upper bound, when ρ → 0 + , for the two terms in the right hand side of (3.2). As S is a closed and rectifiable set, Theorem 0.2 gives

lim ρ→0 + L n (S (1+τ )ρ ) 2ρ = (1 + τ ) S Mν u , ν u 1/2 dH n-1
and then

lim ρ→0 + L n (S (1+τ )ρ ) 2ρ = (1 + τ ) Ju Mν u , ν u 1/2 dH n-1 . (3.3)
In the fourth following Claims, we prove that lim sup

L n (Eρ) 2ρ
converges to 0 when δ converges to 0. The main tool is the regularity result given by Theorem 3.3.

Claim 1: There exists p ∈ N and x 1 , . . . ,

x p ∈ E ∩ J u such that E ⊂ p i=1 B φ (x i , τ ρ) and (i, j) ∈ {1, . . . , p} 2 , i = j ⇒ dist φ (x i , x j ) ≥ τ ρ. ( 3.4)
We construct (x i ) i by an iterative way and we show that the number of iterations is finite. If E = ∅, then the result is obvious. Otherwise, there exists x1 ∈ E. As E ⊂ J u , there exists (y k ) k ⊂ J u converging to x1 . As dist φ ( x1 , S) > τ ρ, there exists k 0 ∈ N such that dist φ (y k0 , S) > τ ρ and we set x 1 = y k0 .

Let us assume that there exits x 1 , . . . , x p ∈ E ∩ J u which satisfy (3.4). If E ⊂ p i=1 B φ (x i , τ ρ), then the iterative process stops. Otherwise, there exists xp+1 ∈ E \ p i=1 B φ (x i , τ ρ). As E ⊂ J u , there exists (y k ) k ⊂ J u converging to xp+1 . As dist φ (x p+1 , S) > τ ρ and dist φ (x p+1 , xi ) > τ ρ for any i ∈ {1, . . . , p}, there exists k p ∈ N such that dist φ (y kp , S) > τ ρ and dist φ (y kp , xi ) > τ ρ. We set x p = y kp .

If the iterative process does not finish, then there exists a sequence (x i ) i∈N ⊂ Ω which satisfies (3.4). According to Proposition 2.1, we have

(i, j) ∈ {1, . . . , p} 2 , i = j ⇒ |x i -x j | ≥ τ ρ √ Λ .
As Ω is bounded there exists a converging subsequence which is a contradiction. So, the iterative process is finite.

Claim 2: There exists a constant c = c(n, λ, Λ) which only depends on the dimension n and the ellipticity coefficients λ, Λ such that, for any x ∈ Ω, we have

# {i ∈ {1, . . . , p} : x ∈ B φ (x i , τ ρ)} ≤ c(n, λ, Λ). We set c(n, λ, Λ) = sup q ∈ N : ∃y 1 , . . . , y q ∈ R n , ∀i, |y i | ≤ 1 √ λ , i = j ⇒ |y i -y j | ≥ 1 √ Λ .
We consider a finite partition of B 0, 1 √ λ of parallelepipeds whose diameter is less than 1 √ Λ . Then, c(n, λ, Λ) is finite and less than the cardinality of such partition. For x ∈ Ω, we set y i = xi-x τ ρ . According to Proposition 2.1, we have

x ∈ B φ (x i , τ ρ) ⇒ y i ∈ B 0, 1 √ λ , i = j ⇒ |y i -y j | ≥ 1 √ Λ .
and then # {i ∈ {1, . . . , p} :

x ∈ B φ (x i , τ ρ)} ≤ c(n, λ, Λ).
Claim 3: We still denote by c a generic constant depending on (n, λ, Λ). We have

pβρ n-1 ≤ c τ n-1 H n-1 (J u ∩ E τ ρ ).
According to Claim 2, we have

p i=1 1 Ju∩B φ (xi,τ ρ) ≤ c1 Ju∩Eτρ . Integrating with respect to H n-1 gives p i=1 H n-1 (J u ∩ B φ (x i , τ ρ)) ≤ cH n-1 (J u ∩ E τ ρ ). Proposition 2.1 condition gives B(x i , λ 1/2 τ ρ)) ⊂ B φ (x i , τ ρ)) and then p i=1 H n-1 (J u ∩ B(x i , λ 1/2 τ ρ)) ≤ cH n-1 (J u ∩ E τ ρ ). (3.5)
According to Theorem 3.3, there exists β > 0 such that, for any i ∈ {1, . . . , p}, we have

H n-1 (J u ∩ B(x i , λ 1/2 τ ρ)) ≥ βλ n-1 2 τ n-1 ρ n-1 . ( 3.6)
Inequalities (3.5) and (3.6) conclude the proof of Claim 3.

Claim 4: We have

lim sup ρ→0 + L n (E ρ ) 2ρ ≤ (1 + τ ) n c τ n-1 δ
where δ is given by (3.1), β by Theorem 3.3 and ω n is the volume of the unit ball of R n .

According to Claim 1, we have E ⊂ p i=1 B φ (x i , τ ρ) and then

E ρ ⊂ p i=1 B φ (x i , (1 + τ )ρ). Proposition 2.1 gives B φ (x i , (1 + τ )ρ) ⊂ B(x i , Λ 1/2 (1 + τ )ρ), it yields L n (E ρ ) ≤ pL n (B φ (x i , (1 + τ )ρ)), ≤ pL n (B(x i , Λ 1/2 (1 + τ )ρ)), ≤ pΛ n 2 (1 + τ ) n ω n ρ n .
As E = K \ S τ ρ , then we deduce E τ ρ ⊂ K τ ρ \ S and Claim 3 gives

L n (E ρ ) ≤ Λ n 2 (1 + τ ) n ω n c βτ n-1 ρH n-1 (J u ∩ E τ ρ ), ≤ Λ n 2 (1 + τ ) n ω n c βτ n-1 ρH n-1 (J u ∩ (K τ ρ \ S)).
According to (3.1), we deduce lim sup

ρ→0 + L n (E ρ ) 2ρ ≤ Λ n 2 (1 + τ ) n ω n c 2βτ n-1 H n-1 (J u ∩ (K \ S)), ≤ Λ n 2 (1 + τ ) n ω n c 2βτ n-1 δ.

Conclusion of the proof.

According to (3.2), (3.3) and Claim 4, we have lim sup

ρ→0 + L n (K ρ ) 2ρ ≤ (1 + τ ) n c τ n-1 δ + (1 + τ ) Ju Mν u , ν u 1/2 dH n-1 .
Letting δ, τ → 0 + successively gives lim sup

ρ→0 + L n (K ρ ) 2ρ ≤ Ju Mν u , ν u 1/2 dH n-1 ,
which concludes the proof of Lemma 3.1. Now, we give the proof of Theorem 3.1.

Proof. We divide the proof in two inequalities.

Claim 1: We have M M (J u ) ≤ Ju Mν u , ν u 1/2 dH n-1 .
We assume that Ju Mν u , ν u 1/2 dH n-1 is finite, otherwise the result is ensured. We set Ω r = {x ∈ Ω : dist(x, ∂Ω) > r}. So, the following set is at most countable Π = r > 0 :

Ju∩∂Ωr Mν u , ν u 1/2 dH n-1 > 0 .
In particular, for any δ > 0, there exists a sequence (r i ) i ⊂ Π strictly decreasing to 0 + and such that Ju∩(Ω\Ωr 0 )

Mν u , ν u 1/2 dH n-1 ≤ δ.

Application to an anisotropic Mumford-Shah model

This work was initially motivated by the problem of segmentation of thin structures in an image.

Let Ω ⊂ R n be the domain of an image g ∈ L ∞ (Ω). In order to favor thin details which are elongated in one or two directions (as tubes or thin plates), we introduce a field of symmetric definite positive matrices M whose principal directions correspond to the elongation of the sets we want to detect. So, we introduce an anisotropic version of the Mumford-Shah model (see [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variationnal problems[END_REF] for the seminal paper)

E(u, M) = Ω (u -g) 2 dx + Ω |∇u| 2 dx + Ju Mν u , ν u 1/2 dH n-1 + Ω DM n+α dx,
where u ∈ SBV(Ω), J u is its jump set and ν u is an unitary and normal vector to J u . The field of matrices M belongs to W 1,n+α (Ω) and α > 0 ensures the continuity of M. Moreover, M is constrained to take its values in a compact subset G of S + n (R). If M ≡ Id n , our model corresponds to the classical Mumford-Shah energy. The novelty of our approach consists in introducing this new variable which must capture the local anisotropy of the image. For example, if we want to detect thin tubes in a 3-D image, we fix µ > 1, the elongation parameter, and we define G as the subset of symmetric matrices whose spectra is {1; 1; µ}. By this way, in order to minimize E, ν u must be orthogonal to the principal directions of M and then this energy will favor set elongated in the principal direction of M.

In order to compute the Euler-Lagrange equation associated to E, we need to perform an approximation of Ju Mν u , ν u 1/2 dH n-1 which involves the Lebesgue measure. That is the reason why we are interested by the previous results on Minkowski contents. We prove the following

Theorem 3.4. Let M : Ω → S + n (R) be continuous, h ∈ L ∞ (Ω), α > 0 and u ∈ SBV(Ω) a minimizer of E α,h,M (v) = α Ω (v -h) 2 dx + Ω |∇v| 2 dx + Jv Mν v , ν v 1/2 dH n-1 : v ∈ SBV(Ω) .
Then, u is an almost quasi-minimizer and we have

M M (J u ) = Ju Mν u , ν u 1/2 dH n-1 .
Proof. We prove in Claim 1 that, up to a scaling, a minimizer is an almost quasi-minimizer. Then, in Claim 2, we compute the re-scaled problem.

Claim 1: For β > 0 and f function defined in Ω we denote by f β the function defined for any x ∈ βΩ by

f β (x) = f x β .
Then, there exists β > 0 such that v β ∈ SBV (βΩ) is an almost quasi-minimizer of a free discontinuity problem.

Let ṽ ∈ SBV (Ω) be a minimizer of E α,h M in SBV(Ω). We introduce

ψ(t) =    -h L ∞ (Ω) if t ≤ -h L ∞ (Ω) , t if |t| ≤ h L ∞ (Ω) , h L ∞ (Ω) if t ≥ h L ∞ (Ω) .
According to the chain rule in SBV (see [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF], Theorem 3.99), for u ∈ SBV(Ω), we have

D(ψ • u) = 1 {x : |u(x)|≤ h L ∞ } ∇uL n + ((ψ • u) + -(ψ • u) -)ν u H n-1 J u and then L n ({x ∈ Ω : |u(x)| > h L ∞ }) > 0 ⇒ E α,h M (ψ • u) < E α,h M (
u). We deduce that |ṽ(x)| ≤ h L ∞ (Ω) for any x ∈ Ω. By an homothetic change of variable, ṽβ is a minimizer of the following rescaled problem

αβ 2 βΩ (v -h β ) 2 dx + βΩ |∇v| 2 dx + β Jv M β ν v , ν v 1/2 dH n-1 : v ∈ SBV(βΩ) . Assuming that v ∈ SBV (βΩ), x ∈ βΩ, B(x, r) ⊂ βΩ, [ṽ β = v] ⊂ B(x, r),
then, the minimality assumption on ṽβ gives

B(x,r) |∇ṽ β | 2 dx + β Jṽ β ∩B(x,r) M β ν ṽβ , ν ṽβ 1/2 dH n-1 ≤ αβ 2 B(x,r) (φ • v -h β ) 2 dx + B(x,r) |∇(φ • v)| 2 dx + β J φ•v ∩B(x,r) M β ν φ•v , ν φ•v 1/2 dH n-1 .
We deduce

B(x,r) |∇ṽ β | 2 dx + β Jṽ β ∩B(x,r) M β ν ṽβ , ν ṽβ 1/2 dH n-1 ≤ 4αβ 2 h 2 L ∞ ω n r n + B(x,r) |∇v| 2 dx + β Jv∩B(x,r) M β ν v , ν v 1/2 dH n-1 ,
where ω n = L n (B(x, 1)). Now, we set β = λ -1 where λ is the ellipticity coefficient introduced in Proposition 2.1, it gives So, we may conclude that ṽβ satisfies the definition 3.1 and then it is an Λ λ , 1,

4α h 2 L ∞ ωn λ 2
-almost quasi-minimizer of a free boundary problem.

Claim 2: We have

M M (J ṽ ) = Jṽ Mν ṽ , ν ṽ 1/2 dH n-1 .

According to Claim 1, we may apply Theorem 3.1 to ṽβ and it gives φ β (x) < ρ}) = β n L n ({x ∈ Ω : J ṽ (x) < βρ}) and then M M β (J ṽβ ) = β n-1 M M (J ṽ ).

M M β (J ṽβ ) =
(3.8)

The following application f : J ṽβ → J ṽ x → x β for any x ∈ J ṽβ satisfies

J n-1 (f )(x) = 1 β n-1 ,
where J n-1 (f ) corresponds to the (n -1)-dimensional Jacobian of f (3.2.1 in [START_REF] Federer | Geometric Measure Theory[END_REF]). According to the change of variable formula between (n -1)-rectifiable subsets (3.2.20 in [START_REF] Federer | Geometric Measure Theory[END_REF]), we deduce that 1 β n-1 Jṽ β M β ν ṽβ , ν ṽβ 1/2 dH n-1 = Jṽ Mν ṽ , ν ṽ 1/2 dH n-1 .

(3.9)

From equations (3.7), (3.8) and (3.9) we conclude that M M (J ṽ ) = Jṽ Mν ṽ , ν ṽ 1/2 dH n-1 .

Conclusion

We have proved that the approximation of Hausdorff measure by Minkowski content may be generalized to continuous and anisotropic perturbation of the Hausdorff measure. It has a natural application in free boundary problems. Classical results of regularity of the jump set may also be generalized in this anisotropic setting. Those results are the key tools to perform Γ-convergence approximation. For example, in a forthcoming work, we will prove that the following energy

E

W 1 , 1 (

 11 [0; 1]; R n ), γ(0) = x, γ(1) = y , dist φ (x, S) = inf {dist φ (x, y) : y ∈ S} . (0.2)Then, the associated anisotropic Minkowski (n -1)-dimensional upper and lower contents are defined by the limitsM M (S) = lim sup ρ→0 + L n ({x : dist φ (x, S) < ρ}) 2ρ , M M (S) = lim inf ρ→0 +

  for the subset of symmetric definite positive matrices, • GL n (R) ⊂ M n (R) for the subset of invertible matrices, • O n (R) ⊂ GL n (R) for the subgroup of orthogonal matrices, • B(Ω) the class of Borelian subsets of Ω, • L n for the Lebesgue measure in R n , • H k for the k-dimensional Hausdorff measure.

H n- 1 (

 1 J ṽβ ∩ B(x, r)) = βλH n-1 (J ṽβ ∩ B(x, r)), ≤ β Jṽ β ∩B(x,r) M β ν ṽβ , ν ṽβ 1/2 dH n-1 , so B(x,r) |∇ṽ β | 2 dx + H n-1 (J ṽβ ∩ B(x, r)) ≤ B(x,r) |∇v| 2 dx + β Jv∩B(x,r) M β ν v , ν v 1/2 dH n-1 + 4αβ 2 h 2 L ∞ ω n r n .The right hand side of Proposition 2.1 givesB(x,r) |∇ṽ β | 2 dx + H n-1 (J ṽβ ∩ B(x, r)) ≤ B(x,r) |∇v| 2 dx + βΛH n-1 (J v ∩ B(x, r)) + 4αβ 2 h 2 L ∞ ω n r n .

Jṽ β M β ν ṽβ , ν ṽβ 1 / 2

 12 dH n-1 . (3.7)We denote by φ β the metric which is defined, for any (x, v) ∈ Ω × R n , byφ β (x, v) = φ x β , vand by d φ β the associated integrated distance. Then, we haveM M β (J) = lim ρ→0 + L n ({x ∈ Ω : d J φ β (x) < ρ}) 2ρ.

ε

  (u, M) = Ω (u -g) 2 dx + Ω |∇u| 2 dx + Ju Mν u , ν u 1/2 dH n-1 + Ω DM n+α dxmay be approached, in the sense of the Γ-convergence, byE ε (u, z, M) = Ω (u -g) 2 dx+ Ω |∇u| 2 (1 -z) 2 dx+ Ω M∇z, ∇z + z 2 4ε dx+ Ω DM n+α dx, when ε → 0 + .

Regularity result for almost quasi-minimizers of a free boundary problem

The standard functional framework is the theory of functions with bounded variation and in particular the special set of functions with bounded variations SBV that may be found in [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF].

Definition and main result

In [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via gamma-convergence[END_REF], it is proved that if u ∈ SBV(Ω) is a minimizer of the well known Mumford-Shah functional

where J u is the jump set of u, then we have H n-1 (J u ) = M(J u ). We want to get the same result when H n-1 is perturbed by an anisotropic and continuous function. So, for u ∈ SBV(Ω) and B(x, r) ⊂ Ω, we set

In [START_REF] Bucur | Monocity formula and regularity for general free discontinuity problems[END_REF], the following definition is introduced. Definition 3.1. For σ ≥ 1, α > 0 and c α > 0, we say that w ∈ SBV(Ω) is a (σ, α, c α )-almost quasi-minimizer of a free discontinuity problem, if there exists α > 0 and c α ≥ 0 such that for any v ∈ SBV(Ω), we have

The main result we introduce in this section is the following.

Proof of Theorem 3.1

To prove Thorem 3.1, we need the two following regularity results for the jump set of almost quasi-minimizers which are extracted from [START_REF] Bucur | Monocity formula and regularity for general free discontinuity problems[END_REF].

Theorem 3.2. Let u be an almost quasi-minimizer of a free discontinuity problem, then

Theorem 3.3. There exist constants β, ρ 0 such that for every (σ, α, c α ) almost quasi-minimizer u, for every x ∈ J u and for every 0 < ρ < ρ 0 such that B(x, ρ) ⊂ Ω, we have

First, we prove the following Lemma 3.1. Let u be an almost quasi-minimizer of a free discontinuity problem (3.1), then for any compact set K ⊂ J u ∩ Ω, we have

Proof. We may assume that Ju Mν u , ν u 1/2 dH n-1 is finite, otherwise the result is ensured. According to Proposition 2.1, we have

We consider the following partition

We may apply Lemma 3.1, with

As r i ∈ Π for any i ≥ 1, we have

As δ > 0 is arbitrary, it proves Claim 1.

Claim 2: We have

As for the proof of Lemma 3.1, there exists a countable family (K i ) i∈N of compact C 1 -hypersurfaces such that

where H n-1 (N ) = 0. As q i=0 K i is rectifiable and closed, Theorem 0.2 gives

So, we deduce that

and then, as H n-1 (N ) = 0, we conclude that


Renaud Barate
email: renaud.barate@ensta.fr

Antoine Manzanera
email: antoine.manzanera@ensta.fr

Evolving Vision Controllers with a Two-Phase Genetic Programming System Using Imitation

We present a system that automatically selects and parameterizes a vision based obstacle avoidance method adapted to a given visual context. This system uses genetic programming and a robotic simulation to evaluate the candidate algorithms. As the number of evaluations is restricted, we introduce a novel method using imitation to guide the evolution toward promising solutions. We show that for this problem, our two-phase evolution process performs better than other techniques.

Introduction

Our goal is to design vision based obstacle avoidance controllers for mobile robots. The most popular method to perform obstacle avoidance with a single camera uses optical flow and is directly inspired by the flight of insects [START_REF] Muratet | A contribution to vision-based autonomous helicopter flight in urban environments[END_REF][START_REF] Zufferey | Fly-Inspired Visual Steering of an Ultralight Indoor Aircraft[END_REF]. However, systems based on optical flow don't cope well with thin or lowly textured obstacles. There is also evidence that information based on appearance rather than movement, like texture information, is extracted in the first stages of the vision chain [START_REF] Malik | Preattentive texture discrimination with early vision mechanisms[END_REF]. This kind of information can also be used for obstacle avoidance. For instance, Michels implemented a system to estimate depth from texture information in outdoor scenes [START_REF] Michels | High speed obstacle avoidance using monocular vision and reinforcement learning[END_REF]. Other systems use this kind of information to discriminate the floor from the rest of the scene and calculate obstacle distances in several directions [START_REF] Ulrich | Appearance-based obstacle detection with monocular color vision[END_REF]. Nevertheless those methods suppose that the floor may be clearly discriminated and they neglect potentially useful information from the rest of the scene.

As there is no method that can deal with all contexts, we want our robot to automatically select and adapt an obstacle avoidance method for the current environment. For now, this adaptation is an offline process based on genetic programming which creates original controllers adapted to a given simulation environment. The next step will be either to install this system on a real robot for online evolution or to design a higher level controller able to select in real-time an algorithm adapted to the current context in a database of evolved algorithms. As we only use artificial evolution as an optimization technique, we won't make assumptions about implications on the development of vision processing in the brain of animals by natural selection and evolution. However there is an important issue addressed by this work which is of interest in a bio-inspired framework: how can Imitation increase the learning speed, and what are the consequences on performances and generalization?

Evolutionary techniques have already been widely used for robotic navigation and the design of obstacle avoidance controllers [START_REF] Walker | Evolving controllers for real robots: A survey of the literature[END_REF] but in general vision is either overly simplified or not used at all. For instance, Marocco used only a 5 × 5 pixels retina as visual input [START_REF] Marocco | Active vision and feature selection in evolutionary behavioral systems[END_REF]. On the other hand, genetic programming has been proved to achieve human-competitive results in image processing systems, e.g. for the detection of interest points [START_REF] Trujillo | Synthesis of interest point detectors through genetic programming[END_REF]. Parisian evolution has also been shown to produce very good results for obstacle detection and 3D reconstruction but those systems need two calibrated cameras [START_REF] Pauplin | Evolutionary Optimisation for Obstacle Detection and Avoidance in Mobile Robotics[END_REF].

To our knowledge, only Martin tried evolutionary techniques with monocular images for obstacle avoidance [START_REF] Martin | Evolving visual sonar: Depth from monocular images[END_REF]. The structure of his algorithm is based on the floor segmentation technique and the evaluation is done with a database of hand labeled real world image. The advantage of such an approach is that the evolved algorithms are more likely to work well with real images than those evolved with computer rendered images. Nevertheless, it introduces an important bias since the algorithms are only selected on their ability to label images in the database correctly and not on their ability to avoid obstacles.

Evolution of the Vision Algorithms

Structure of the Vision Algorithms

Generally speaking, a vision algorithm can be divided in three main parts: First, the algorithm will process the input image with a number of filters to highlight some features. In our case, this filter chain consists of spatial and temporal filters, optical flow calculation and projection that will produce an image highlighting the desired features. Then these features are extracted, i.e. represented by a small set of scalar values. We implemented this with a mean computation of the pixel values on several windows of the transformed image(s). Finally these values are used for a domain dependent task, here to generate motor commands to avoid obstacles. Fig. 1 shows an example program for a simple obstacle avoidance behavior based on optical flow. Here is the list of all the primitives (transformation steps) that can be used in the programs and the data types they manipulate:

-Spatial filters (input: image, output: image): Gaussian, Laplacian, threshold, Gabor, difference of Gaussians, Sobel and subsampling filter. -Temporal filters (input: image, output: image): pixel-to-pixel min, max, sum and difference of the last two frames, and recursive mean operator. -Optical flow (input: image, output: vector field ): Horn and Schunck global regularization method, Lucas and Kanade local least squares calculation and simple block matching method. The rotation movement is first eliminated by a transformation of the two images in order to facilitate further use of the optical flow. -Projection (input: vector field, output: image): Projection on the horizontal or vertical axis, Euclidean or Manhattan norm computation, and time to contact calculation using the flow divergence. -Windows integral computation (input: image, output: scalar ): For this transformation, we define a global coefficient α 0 and several windows on the left half of the image with different positions and sizes. With each window is paired a second window defined by symmetry along the vertical axis.

A coefficient α i and an operator (+ or -) are defined for each pair. The resulting scalar value R is a simple linear combination calculated with the following formula:

R = α 0 + n i=1 α i µ i µ i = µ Li + µ Ri or µ i = µ Li -µ Ri
where n is the number of windows and µ Li and µ Ri are the means of the pixel values over respectively the left and right window of pair i.

-Scalar operators (input: scalar(s), output: scalar ): Addition, subtraction, multiplication and division operators, temporal mean calculation and simple if-then-else test. -Command generation (input: two scalars, output: command ): The motor command is represented by two scalar values: the requested linear and angular speeds.

Most of those primitives use parameters along with the input data to do their calculations (for example, the standard deviation value for the Gaussian filter or the number and position of windows for the windows integral computation). Those parameters are specific to each algorithm; they are randomly generated when the corresponding primitive is created by the evolution process.

Evaluation in the Simulation Environment

For the evaluation of the different obstacle avoidance algorithms, we use a simulation environment in which the robot moves freely during each experiment. The simulation is based on the open-source robot simulator Gazebo. The simulated camera produces 8-bits gray-value images of size 320×160 (representing a field of view of approximately 100 • × 60 • ) at a rate of 10 / sec. The simulation environment is a closed room of 36 m 2 area (6 m × 6 m) containing three bookshelves (Fig. 1). All the obstacles are immovable to prevent the robot from just pushing them instead of avoiding them. In each experiment, the goal of the robot is to go from a given starting point to a goal location without hitting obstacles.

Due to the complexity of the vision algorithms, we're limited to about 40,000 evaluations to keep the evolution time acceptable (a few days at most). We therefore face a common problem with genetic programming systems, that is the state space is immense compared to the number of evaluations so large parts of it will never be explored. More, the fitness landscape is very chaotic so the evolution can easily get stuck in a local minima of the fitness function.

In previous work we showed that a classical evolution process often produces controllers with a seemingly random trajectory, even if the visual features they use are coherent with the environment [START_REF] Barate | Automatic Design of Vision-based Obstacle Avoidance Controllers using Genetic Programming[END_REF]. To overcome those problems, we propose a novel approach based on the imitation of a given behavior in a first phase to guide the evolution toward more efficient solutions in the second phase.

First Phase: Evolution of Algorithms that Imitate a Behavior

In this first phase, the population is initialized with random algorithms and the evolution lasts for 50 generations. The goal will be to match a recorded example behavior. More precisely, we first record the video sequence and parameters of an experiment where we manually guide the robot from the starting point to the target point while avoiding obstacles. We make the trajectory as short and smooth as possible to limit the difficulty of the matching task.

For the evaluation of the algorithms, we replay this sequence and compare the command issued by the evaluated algorithm with the command recorded during the manual control of the robot. The goal is to minimize the difference between these two commands along the recorded sequence. Formally, we try to minimize two variables F and Y defined by the formulas:

F = n i=1 (f Ri -f Ai ) 2 and Y = n i=1 (y Ri -y Ai ) 2
where f Ri and y Ri are the recorded forward and yaw speed commands for frame i, f Ai and y Ai are the forward and yaw speed commands from the tested algorithm for frame i and n is the number of frames in the video sequence.

Second Phase: Evolution of Efficient Solutions

In this second phase, the population is initialized with the final population of the first phase and the evolution lasts for 50 more generations. The algorithms will be evaluated on their ability to really avoid obstacles and reach the target location. For that, we place the robot at the fixed starting point and let it move in the environment during 30 s driven by the obstacle avoidance algorithm. Two scores are attributed to the algorithm depending on its performance: a goalreaching score G rewards algorithms reaching or approaching the goal location, whereas contact score C rewards the individuals that didn't hit obstacles on their way. Those scores are calculated with the following formulas:

G = t G if the goal is reached t max + d min /V otherwise C = t C
where t G is the time needed to reach the goal in seconds, t max is the maximum time in seconds (here 30 s), d min is the minimum distance to the goal achieved in meters, V is a constant of 0.1 m/s and t C is the time spent near an obstacle (i.e. less than 18 cm, which forces the robot to keep some distance away from obstacles). The goal is hence to minimize those two scores G and C.

For the two phases of evolution, the evaluations consist in fact in two runs with a different starting point and a different goal location. Final scores are the sum of the scores obtained for the two runs. Performing two different runs favors algorithms with a real obstacle avoidance strategy while not increasing evaluation time too much. The starting points are fixed because we want to evaluate all algorithms on the same problem.

The Genetic Programming system

We use grammar based genetic programming to evolve the vision algorithms [START_REF] Whigham | Grammatically-based genetic programming[END_REF]. As usual with artificial evolution, the population is initially filled with randomly generated individuals. In the same way that a grammar can be used to generate syntactically correct random sentences, a genetic programming grammar is used to generate valid algorithms. The grammar defines the primitives, data and the rules that describe how to combine them. The generation process consists in successively transforming each non-terminal node of the tree with one of the rules. This grammar is used for the initial generation of the algorithms and for the transformation operators. The crossover consists in swapping two subtrees issuing from identical non-terminal nodes in two different individuals. The mutation consists in replacing a subtree by a newly generated one. Table 1 presents the exhaustive grammar that we used in all our experiments.

The numbers in brackets are the probability of selection for each rule. A major advantage of this system is that we can bias the search toward the usage of more promising primitives by setting a high probability for the rules that generate them. We can also control the size of the tree by setting small probabilities Table 1. Grammar used in the genetic programming system for the generation of the algorithms.

for the rules that are likely to cause an exponential growth (rules like real → ifThenElse(real,real,real,real) for example).

As described previously, we wish to minimize two criteria (F and Y in the first phase, G and C in the second phase). There are different ways to use evolutionary algorithms to perform optimization on several and sometimes conflicting criteria. For the experiments described in this paper, we chose the widely used multiobjective evolutionary algorithm called NSGA-II. It is an elitist algorithm based on the non-dominance principle. A diversity metric called "crowding distance" is used to promote diversity among the evolved individuals. All the details of the implementation can be found in the paper by K. Deb [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF].

In order to prevent problems of premature convergence, we separate the population of algorithms in 4 islands, each containing 100 individuals. Those islands are connected with a ring topology; every 10 generations, 5 individuals selected with binary tournament will migrate to the neighbor island while 5 other individuals are received from the other neighbor island. For the parameters of the evolution, we use a crossover rate of 0.8 and a probability of mutation of 0.01 for each non-terminal node. We use a classical binary tournament selection in all our experiments. Those parameters were determined empirically with a few tests using different values. Because of the length of the experiments, we didn't proceed to a thorough statistical analysis of the influence of those parameters.

Experiments and Results

Comparison with Other Evolution Strategies

We compare here our own system with a classical one-phase evolution process and with two methods commonly used to guide the evolution toward promising solutions, namely incremental evolution and seeding. For these three experiments we only use the evaluation function described in 2.4, to select the controllers on their ability to reach the goal while avoiding obstacles.

Seeding an evolution consists in inserting in the initial population one or several individuals (generally hand written solutions) performing reasonably well on the target problem. This seeding method guides the evolution towards a specific part of the state space where we know that there is a good probability to find interesting solutions. Here we simply inserted our manually designed controller presented before in the initial population of one of the four islands. The evolution then proceeds normally for 100 generations.

Incremental evolution [14] is based on the decomposition of a problem in a set of problems of increasing complexity. Starting with a simple problem, the evolution quickly finds good solutions. The difficulty is then progressively increased and new individuals to adapt to those more and more complex conditions. Here, we proceed for the 20 first generations with a simplified environment with only one bookshelf (the one that is closest to the wall). This way, there is no real obstacle between the starting points and the target points. Then we proceed for 40 more generations with one more bookshelf, and finally for the 40 last generations with the complete environment containing three bookshelves.

In total, all the experiments last for 100 generations, that is 40,000 evaluations. Fig. 2 presents the Pareto fronts at different times of the evolution for the different methods. With the classical one-phase evolution, the first generation contains only individuals with poor performance. After 50 generations we obtain controllers with interesting behaviors, performing better than the hand designed controller. The end of the evolution doesn't improve these results much, we can say here that the evolution process has become stuck in a local minimum of the fitness function. With seeding, the first generations contain better controllers due to the seed individual but the evolution doesn't manage to further improve this behavior. In the incremental evolution experiment, the evolution finds very efficient controllers in the two first environments. But when we introduce the third bookshelf, the evolution doesn't find a way to adapt to those new conditions. Nevertheless the final controllers perform better than in the previous cases. With two-phase evolution, we obtain very good individuals even at the end of the first phase. The second phase will then optimize those controllers and the final ones definitely surpass the ones obtained with the other methods.

Analysis of the Controllers

We represent on Fig. 3 the trajectories of the robot in the environment with controllers issued from the different evolution processes. The hand written controller manages to avoid obstacles but the robot goes too slowly and never reaches the target points. With the one-phase evolution process, the evolved controllers sometimes manage to reach one of the target points without hitting obstacles but generally they move in an almost random way and with a lot of rotations. Turning quickly on themselves prevents the robots from getting stuck against obstacles and thus generally improves their goal reaching score in the first generations. This kind of controllers quickly overcomes the rest of the population and prevents the evolution from finding better controllers. As shown before, the seeded evolution didn't manage to improve the seed individual much. The trajectory obtained with evolved controllers is very close to the trajectory of the hand written controller used as seed.

The incremental evolution produces efficient solutions with one and two bookshelves but when we add the final one, the evolution fails to adapt the controllers. The trajectory clearly shows this problem: one of the bookshelf is quickly avoided but the robot just stops before hitting the second one without trying to avoid it. The first phase here selects very simple controllers going straight on to the target points. To achieve this, they just use the target heading value for the yaw speed command. Although it is a good strategy without obstacles, these controllers are difficult to adapt afterward for more complex problems.

With two-phase evolution, the first phase evolves controllers with smooth trajectories very similar to the recorded sequence. Some of those controllers even reach the two target points although most of them get stuck on an obstacle as shown on Fig. 3. With this initialization, the second phase quickly evolves very efficient controllers reaching the two target points without hitting obstacles and even faster than when the robot is manually controlled. The structure of the controllers evolved in this environment is generally quite simple (mainly based on contrast information obtained with simple threshold filters) and there's almost no structure difference between the end of the first phase and the end of the second phase. This means that the imitation phase is crucial for the evolution of the structure to produce efficient controllers. The second phase mainly optimizes the parameters of the algorithms to obtain a faster and more robust behavior.

The main limitation of all those controllers is that they don't generalize well to different conditions. If we move around the obstacles, the starting point and the target point, the performance of the controllers decreases quickly. The controller issued from two-phase evolution still shows some obstacle avoidance abilities but it doesn't move toward the target point. Some other controllers just seem to move randomly in the environment. There are two reasons for this mediocre generalization behavior: First, it seems that the choice of only two trajectories in the evaluation environment was too optimistic. The controllers overlearn these trajectories and this decreases their generalization performance. Second, the evolved controllers often use the target heading information in quite complicated ways or even not at all. In this case when we move the target point they are not able to find it unless by chance. To overcome these problems we have designed experiments with more trajectories during the evolution and with a different grammar which facilitates an efficient use of the target heading for the controllers. The generalization performance of the controllers evolved this way is presented in [START_REF] Barate | Generalization Performance of Vision Based Controllers for Mobile Robots Evolved with Genetic Programming[END_REF].

Conclusion

In this paper, we present a system that selects and adapts automatically an appropriate obstacle avoidance method to a given visual context. This system is based on genetic programming and we introduced a novel evolution process to guide the evolution toward promising solutions in the state space. We have shown that for our target application, it improves the performance of the final controllers greatly compared to classical evolution, seeding or incremental evolu-tion. This method could be adapted for many evolutionary robotic systems and it would be interesting to validate it on other applications. Our goal is now to test this method on a wide range of environments (including real-world environments) to check if some primitives are selected more often than others and how they are combined by the evolution. This will lead to a better understanding on how to combine different vision information in a robotic navigation system in order to obtain more adaptive and robust behaviors.

Fig. 1 .

 1 Fig. 1. Left: Algorithmic tree of a program example for obstacle avoidance. Rectangles represent primitives and ellipses represent data. Right: Snapshots of the simulation environment.

Fig. 2 .

 2 Fig.2. Pareto fronts at different times of evolution for the different evolution strategies. Those fronts are created by taking the best individuals in the generation in all the four islands and keeping only non-dominated ones. In some cases one individual dominates all others, reducing the front to one point.

Fig. 3 .

 3 Fig.3. a) Trajectory of the hand written controller. b) Trajectory of a controller evolved with one-phase evolution. c) Trajectory of a controller evolved with seeded evolution. d) Trajectory of a controller evolved with incremental evolution. e) Trajectory of the recorded sequence. f) Trajectory of a controller evolved with only example matching evaluation during 50 generations (end of the first phase). g) Trajectory of the controller evolved with two-phases evolution process (end of the second phase). The small dots on the trajectories are placed every 2 seconds and thus indicate the speed of the robot.


Renaud Barate
email: renaud.barate@ensta.fr

Antoine Manzanera
email: antoine.manzanera@ensta.fr

Generalization Performance of Vision Based Controllers for Mobile Robots Evolved with Genetic Programming

Keywords: I.2.2 [Artificial Intelligence]: Automatic Programmingprogram synthesis; I.4.8 [Image Processing and Computer Vision]: Scene Analysis-depth cues Vision, obstacle avoidance, robotic simulation, generalization

We present a genetic programming system to evolve vision based obstacle avoidance algorithms. In order to develop autonomous behavior in a mobile robot, our purpose is to design automatically an obstacle avoidance controller adapted to the current context. Using a simulation environment, several trajectories are used to evaluate the performance of an evolving population of controllers. Two different methods are proposed and compared. In the first one (structurefree controller), we discard any a priori in the structure of the algorithm. We show that the evolved controllers behave well in the evolution environment, but have very limited generalization capabilities. In the second method (structure restricted controller), the compromise between the two antagonist functions of the robot navigation, i.e. obstacle avoidance and target reaching, is explicitely integrated in the structure of the algorithms. We show that controllers evolved with this second method are more generic. We give hints in the discussion to detect overlearning earlier in the evolution process, so as to quickly evolve controllers with good generalization abilities.

INTRODUCTION

Our goal is to automatically design vision based controllers for robots. We chose obstacle avoidance as target application as it is a necessary task for any mobile robot. One commonly used method to perform obstacle avoidance using a single camera is based on the calculation of optical flow. Optical flow represents the perceived movement of the objects and can be used to compute a fairly accurate approximation of their distance. It has been used recently for obstacle avoidance with an autonomous helicopter for instance [START_REF] Muratet | A contribution to vision-based autonomous helicopter flight in urban environments[END_REF]. However, systems based on optical flow don't cope well with thin or lowly textured obstacles. On the other hand, more simple primitives can be used to extract useful informations in the image, generally based on texture. For instance, Michels implemented a system to estimate depth from texture information in outdoor scenes [START_REF] Michels | High speed obstacle avoidance using monocular vision and reinforcement learning[END_REF]. Other obstacle avoidance systems use this kind of information to discriminate the floor from the rest of the scene and calculate obstacle distances in several directions [START_REF] Horswill | Polly: A vision-based artificial agent[END_REF][START_REF] Lorigo | Visually-guided obstacle avoidance in unstructured environments[END_REF][START_REF] Ulrich | Appearance-based obstacle detection with monocular color vision[END_REF]. Nevertheless those methods suppose that the floor may be clearly discriminated and they neglect potentially useful information from the rest of the scene.

As none of those methods is robust enough to deal with all kind of environments, our system should be able to automatically select the kind of visual features that allow obstacle detection in a given environment. More, it should adapt and parameterize the whole obstacle avoidance controller to use those features efficiently. For that, we propose to use genetic programming as it allows us to evolve algorithms with little a priori on their structure and parameters. For now, the evolution is an offline process using a simulation environment for the evaluations where the robot must reach a target point while avoiding obstacles.

Evolutionary techniques have already been used for robotic navigation and the design of obstacle avoidance controllers but in general vision is either overly simplified or not used at all. For instance, Harvey used a 64×64 pixels camera but simplified the input to a few average values on the image [START_REF] Harvey | Seeing the Light: Artificial Evolution, Real Vision[END_REF]. Marocco used only a 5 × 5 pixels retina as visual input [START_REF] Marocco | Active vision and feature selection in evolutionary behavioral systems[END_REF]. A description of several other evolutionary robotics systems can be found in [START_REF] Mataric | Challenges in evolving controllers for physical robots[END_REF][START_REF] Walker | Evolving controllers for real robots: A survey of the literature[END_REF] but most of them rely only on range sensors. On the other hand, genetic programming has been proved to achieve human-competitive results in image processing systems, e.g. for the detection of interest points [START_REF] Ebner | Evolving a task specific image operator[END_REF]15]. Parisian evolution has also been shown to produce very good results for obstacle detection and 3D reconstruction but those systems need two calibrated cameras [14,[START_REF] Olague | Parisian evolution with honeybees for three-dimensional reconstruction[END_REF].

To our knowledge, only Martin tried evolutionary tech-niques with monocular images for obstacle avoidance [START_REF] Martin | Evolving visual sonar: Depth from monocular images[END_REF]. The structure of his algorithm is based on the floor segmentation technique and the evaluation is done with a database of hand labeled real world image. The advantage of such an approach is that the evolved algorithms are more likely to work well with real images than those evolved with computer rendered images. Nevertheless, it introduces an important bias since the algorithms are only selected on their ability to label images in the database correctly and not on their ability to avoid obstacles. We first tried to evolve controllers with very little a priori on their algorithmic structure. In this case, all the algorithm is designed by the evolutionary process. The advantage of this freedom is that every kind of algorithm can potentially be evolved this way, even those that we didn't think of while designing the system. Nevertheless the compromise between moving toward the target point and avoiding obstacles is very difficult to achieve with this kind of free structure. In general the evolution will either find controllers moving toward the target point but with limited obstacle avoidance abilities, or controllers that avoid obstacles but that won't move toward the target location. Therefore this kind of controllers won't generalize well to slightly different conditions.

To overcome this problem, we designed a more restricted algorithmic structure that facilitates the compromise between avoiding obstacles and moving toward the target location. We will show that those controllers are more generic as they can reach the target point while avoiding obstacles even when the environment changes.

STRUCTURE OF THE CONTROLLERS

The Vision Algorithms

Generally speaking, a vision algorithm can be divided in three main parts: First, the algorithm will process the input image with a number of filters to highlight some features. Then these features are extracted, i.e. represented by a small set of scalar values. Finally these values are used for a domain dependent task, here to generate motor commands to avoid obstacles. We designed the structure of our algorithms according to this general scheme. First, the filter chain consists of spatial and temporal filters, optical flow calculation and projection that will produce an image highlighting the desired features. Then we compute the mean of the pixel values on several windows of this transformed image (feature extraction step). Finally those means are used to compute a single scalar value by a linear combination. We will use this scalar value to determine the presence of an obstacle and to generate a motor command to avoid it.

An algorithm is represented as a tree, the leaves being input data, the root being output command, and the internal nodes being primitives (transformation steps). The program can use different types of data internally, i.e. scalar values, images, optical flow vector fields or motor commands. For each primitive, the input and output data types are fixed. Some primitives can internally store information from previous states, thus allowing temporal computations like the calculation of the optical flow. Here is the list of all the primitives that can be used in the programs and the data types they manipulate:

• Spatial filters (input: image, output: image): Gaussian, Laplacian, threshold, Gabor, difference of Gaus-sians, Sobel and subsampling filter.

• Temporal filters (input: image, output: image): pixelto-pixel min, max, sum and difference of the last two frames, and recursive mean operator.

• Optical flow (input: image, output: vector field ): Horn and Schunck global regularization method, Lucas and Kanade local least squares calculation and simple block matching method [START_REF] Barron | Performance of optical flow techniques[END_REF]. The rotation movement is first eliminated by a transformation of the two images in order to facilitate further use of the optical flow.

• Projection (input: vector field, output: image): Projection on the horizontal or vertical axis, Euclidean or Manhattan norm computation, and time to contact calculation using the flow divergence.

• Windows integral computation (input: image, output: scalar ): The method used for this transformation is:

1. A global coefficient α0 is defined for the primitive.

2. Several windows are defined on the left half of the image with different positions and sizes. With each window is paired a second window defined by symmetry along the vertical axis. A coefficient αi and an operator (+ or -) are defined for each pair.

3. The resulting scalar value R is a simple linear combination calculated with the following formula:

R = α0 + n i=1 αiµi µi = µLi + µRi or µi = µLi -µRi
where n is the number of windows and µLi and µRi are the means of the pixel values over respectively the left and right window of pair i.

The number of windows pairs, their positions, sizes, operator and coefficient along with the global coefficient are characteristic parts of the primitive and will be customized by the evolutionary process.

• Scalar operators (input: scalar(s), output: scalar ):

Addition, subtraction, multiplication and division operators and a temporal mean calculation.

• Command generation (input: two scalars, output: command ): The motor command is represented by two scalar values: the requested linear and angular speeds.

Most of those primitives use parameters along with the input data to do their calculations (for example, the standard deviation value for the Gaussian filter). Those parameters are specific to each algorithm; they are randomly generated when the corresponding primitive is created by the genetic programming system described in Sect. 3.4.

Structure-Free Controllers

For this first kind of controllers, all the algorithmic tree will be built by the evolutionary process. We add two scalar input variables: the target distance and the target direction, and one scalar operator: an if-then-else test. This allows the evolutionary process to create more complex and non-linear combinations between those scalar input variables and the scalars issuing from the vision part of the algorithm. Fig. 1 shows an example controller built with this structure. The algorithm is in this case a single tree, and all the primitives between the input data and the resulting motor command are created and parameterized by the evolution.

Structure-Restricted Controllers

Here, we restrict the global structure of the controllers to facilitate the compromise between avoiding obstacles and going toward the target point. First, a vision algorithm is used to detect nearby obstacles. If no obstacle is around, the robot will just go straight to the target point. If an obstacle is detected, another vision algorithm will be used to generate a command to avoid this obstacle. The parts that will be customized by the evolution are the two vision algorithms along with a few parameters like the speed of the robot when going straight to the goal.

The first vision algorithm (obstacle detection) is in fact always linear, starting with the video image and computing a Boolean as output. We obtain this Boolean value by comparing the scalar produced by the feature extraction step with a scalar threshold. This result will indicate the presence or absence of a nearby obstacle.

The second vision algorithm (command generation) can use as input either the original video image or the image filtered by the obstacle detection algorithm. This allows the reuse of interesting features between the two algorithms. Here, several scalar values can be produced with different filter chains and these values are combined using scalar operators. The final step will generate a motor command using the resulting scalar value(s). Fig. 2 illustrates this restricted structure with an example.

THE EVOLUTION PROCESS

Simulation Environment

For the evaluation of the different obstacle avoidance algorithms, we use a simulation environment in which the The simulator uses ODE physics engine for the movement of the robot and collisions detection and OpenGL for the rendering of the camera images. The physics engine update rate is 50 Hz while the camera update rate is 10 Hz. In all the experiments presented in this paper, the simulated camera produces 8-bits gray-value images of size 320 × 160 representing a field of view of approximately 100 • × 60 • . This large field of view reduces the dead angles and hence facilitates obstacle detection and avoidance. The simulation environment is a closed room of 36 m 2 area (6 m × 6 m) containing three bookshelves (Fig. 3). All the obstacles are immovable to prevent the robot from just pushing them instead of avoiding them. In each experiment, the goal of the robot is to go from a given starting point to a goal location without hitting obstacles.

First Phase: Evolution of Algorithms that Imitate a Behavior

Due to the simulation time and mainly to the complexity of the vision algorithms, we're limited to about 40,000 evaluations to keep the evolution time acceptable (a few days at most). We therefore face a common problem with genetic programming systems, that is the state space is immense compared to the number of evaluations. More, the fitness landscape is very chaotic so the evolution can easily get stuck in a local minima of the fitness function. To overcome those problems, we use a two-phase evolution process: the first phase is based on imitation to quickly guide the evolution toward promising solutions, while the second phase will evolve more robust controllers by using a functional evaluation in the simulation environment.

With this kind of problem, it is very difficult to manually design even a mediocre controller but it is very easy to manually guide the robot toward the target point while avoiding obstacles. So we can obtain a good example of an efficient behavior by recording the video sequence and command parameters while we guide the robot. In this first phase, we will try to evolve algorithms that imitate this efficient example behavior.

The population is initialized with random algorithms and the evolution lasts for 50 generations. For the evaluation of the algorithms, we replay the recorded sequence and compare the command issued by the evaluated algorithm with the command recorded during the manual control of the robot. The goal is to minimize the difference between these two commands along the recorded sequence. Formally, we try to minimize two variables F and Y defined by the formulas:

F = n i=1 (fRi -fAi) 2 and Y = n i=1 (yRi -yAi) 2
where fRi and yRi are the recorded forward and yaw speed commands for frame i, fAi and yAi are the forward and yaw speed commands from the tested algorithm for frame i and n is the number of frames in the video sequence.

Second Phase: Evolution of Efficient Solutions

The goal of this second phase is to produce more robust and generic controllers. For that, we will evaluate them on their ability to really avoid obstacles and reach the target location in the simulation environment. The population is initialized with the final population of the first phase and the evolution lasts for 50 more generations. For the evaluation, we place the robot at the fixed starting point and let it move in the environment during 30 s driven by the obstacle avoidance algorithm. Two scores are attributed to the algorithm depending on its performance: a goal-reaching score G rewards algorithms reaching or approaching the goal location, whereas contact score C rewards the individuals that didn't hit obstacles on their way. Those scores are calculated with the following formulas:

G = tG if the goal is reached tmax + dmin/V otherwise C = tC
where tG is the time needed to reach the goal in seconds, tmax is the maximum time in seconds (here 30 s), dmin is the minimum distance to the goal achieved in meters, V is a constant of 0.1 m/s and tC is the time spent near an obstacle (i.e. less than 18 cm, which forces the robot to keep some distance away from obstacles). The goal is hence to minimize those two scores G and C.

For the two phases of evolution, the evaluations consist in fact in four runs with different starting points and goal locations. Final scores are the means of the scores obtained for the four runs. The starting points are fixed because we want to evaluate all algorithms on the same problem. We chose four runs as a tradeoff between generalization behavior and evaluation time but more tests would be needed to determine the optimal number of runs and the precise consequence on generalization.

The Genetic Programming System

We use genetic programming to evolve vision algorithms with little a priori on their structure. As usual with evolutionary algorithms, the population is initially filled with randomly generated individuals. We use the grammar based genetic programming system introduced by Whigham [START_REF] Whigham | Grammatically-based genetic programming[END_REF] to overcome the data typing problem. It also allows us to bias the search toward more promising primitives and to control the growth of the algorithmic tree.

In the same way that a grammar can be used to generate syntactically correct random sentences, a genetic programming grammar is used to generate valid algorithms. The grammar defines the primitives and data (the bricks of the algorithm) and the rules that describe how to combine them.

The generation process consists in successively transforming

Table 1: Grammar used in the genetic programming system for the creation and transformation of the command generation algorithm for structurerestricted controllers. The grammars used for the obstacle detection algorithm and for the structurefree controllers are very similar to this one. each non-terminal node of the tree with one of the rules. This grammar is used for the initial generation of the algorithms and for the transformation operators. The crossover consists in swapping two subtrees issuing from identical nonterminal nodes in two different individuals. The mutation consists in replacing a subtree by a newly generated one. Table 1 presents the exhaustive grammar that we used in all our experiments.

The numbers in brackets are the probability of selection for each rule. A major advantage of this system is that we can bias the search toward the usage of more promising primitives by setting a high probability for the rules that generate them. We can also control the size of the tree by setting small probabilities for the rules that are likely to cause an exponential growth (rules like real → add(real,real) for example).

As described previously, we wish to minimize two criteria (F and Y in the first phase, G and C in the second phase). There are different ways to use evolutionary algorithms to perform optimization on several and sometimes conflicting criteria. For the experiments described in this paper, we chose the widely used multi-objective evolutionary algorithm called NSGA-II. This algorithm is based on the Pareto dominance principle. Individuals are sorted by non-dominance rank, so that non-dominated individuals get a higher probability of being selected for breeding. This algorithm is elitist and a "crowding distance" is used to promote diversity among the individuals. More details can be found in the paper by K. Deb [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF].

In order to prevent problems of premature convergence, we separate the population of algorithms in 4 islands, each containing 100 individuals. Those islands are connected with a ring topology; every 15 generations, 5 individuals selected with binary tournament will migrate to the neighbor island while 5 other individuals are received from the other neighbor island. For the parameters of the evolution, we use a crossover rate of 0.8 and a probability of mutation of 0.01 for each non-terminal node. We use a classical binary tournament selection in all our experiments. Those parameters were determined empirically with a few tests using different values. Due to the length of the experiments, we didn't proceed to a thorough statistical analysis of the influence of those parameters.

EXPERIMENTS AND RESULTS

Analysis of the Evolved Controllers

In total, each experiment lasts for 100 generations, that is 40,000 evaluations. Fig. 4 presents the final Pareto fronts for the evolution of structure-free and structure-restricted controllers. Those final results are comparable in performance for the two kinds of controllers. Structure-restricted controllers generally come closer to the target points but structure-free controllers are better for avoiding obstacles. The real difference appear when we compare the trajectories of the controllers. For that we take an example solution in the final population of each experiment (individuals A and B in Fig. 4). Those individuals are not on the Pareto front but they show good generalization abilities (see Sec. 4.2). Fig. 5 represents the trajectories of the robot in the simulation environment used during the evolution with the structure-restricted controller A. We can see that it applies a very simple back and forth strategy: it goes forward toward the goal as long as there is no obstacle in front of it, and it goes backward turning on itself when it faces an obstacle. Nevertheless this simple strategy is sufficient as the robot reaches the target point in each case and with only a small number of contacts.

Starting point

Target point Trajectory Contact points Fig. 6 represents the trajectories in the same environment with the structure-free controller B. Here the robot moves much faster than in the previous case but it doesn't use the target direction information at all, so it will never find the target point unless by chance.

Generalization Performance

The previous experiment showed the performance of the evolved controllers in the environment that was used for the evolution process but it does not prove that they really can avoid obstacles in a generic way. They could just as well be adapted only to reproduce the trajectory used during the evolution. In this case, they would have learned a trajectory and not a generic obstacle avoidance behavior.

In order to test the genericity of the evolved controllers, we designed a second environment with the same appearance but where the obstacles and the start and target points have been moved around. We also have four different runs with different start and target points. We evaluate all the individuals of the last generation from the previous evolution

Starting point

Target point Trajectory Contact points In this environment, the structure-restricted controllers generally outperform the structure-free ones. This is even clearer when we consider the trajectories of the robot with the two kinds of controllers (Fig. 8 and9). The structurerestricted controller A adapts quite well to those new conditions and reaches the target point in three out of four test cases, with only a few more contacts with the obstacles than before. The structure-free controller B still avoids the obstacles correctly but since it does not move toward the target point, it finds it by chance only in one of the test cases.

Starting point

Target point Trajectory Contact points Those experiments show that using a customized structure for the evolution of algorithms can bring significant benefits, especially when we must find a compromise between two different behaviors. Indeed we can evolve very efficient obstacle avoidance controllers or controllers that move toward the goal with a free structure, but only by using a kind of rulebased structure we can evolve controllers that move toward a target point while avoiding obstacles.

Another interesting result is that the best individuals in the environment used during the evolution process are generally not the best ones when we change the environment. This can be explained by the fact that those individuals overlearn the environment used to evaluate them and this decreases their generalization performance. Therefore the question is: how to evolve and select individuals with good generalization abilities? We partially answered this question by using four different test cases in the evaluation environment. This limits the overlearning problem as the evolved controllers must be adapted for the four different trajectories. By evaluating the whole final population on a new environment and selecting individuals performing well with those new conditions, we also choose the solutions that are more likely to generalize well.

But this does not tell us when we should stop the evolution process to avoid overlearning. Gagné [START_REF] Gagné | Genetic Programming, Validation Sets, and Parsimony Pressure[END_REF] proposed to use a three data sets methodology, already commonly used in the Machine Learning community. This consists in testing the individuals on the Pareto front after each generation on a validation set to detect when overlearning begins to occur and therefore to stop the evolution process. But as we have

Starting point

Target point Trajectory Contact points

Figure 9: Trajectories with the structure-free controller B in the generalization environment seen before, the more generic individuals are not always on the Pareto front. We could evaluate the whole population on a new environment after each generation but this would double the evolution time and we would select individuals that are "by chance" adapted to our validation set. Nevertheless we plan to test and adapt this kind of technique in future work in order to select more generic individuals and to shorten the evolution time.

CONCLUSION

In this paper, we presented the benefits of using a specific algorithmic structure for the evolution of vision based obstacle avoidance controllers with genetic programming. This structure allows us to facilitate the compromise between going toward a target point and avoiding obstacles. Therefore the controllers evolved with this structure are more generic than structure-free controllers.

Our goal is now to further improve this system to facilitate the detection of different kind of obstacles and to detect overlearning earlier in the evolution process. We will then adapt the system to create controllers for a real robot. Finally, we want to design a high-level controller able to select its behavior depending on the context in real-time, in order to have a fully autonomous and adaptive robot.

Figure 1 :

 1 Figure 1: Example of a structure-free controller for obstacle avoidance based on optical flow. Rectangles represent primitives and ellipses represent data.

Figure 2 :

 2 Figure 2: Example of a structure-restricted controller. Top: overview of the global fixed structure. Bottom left: an example of obstacle detection algorithm. Bottom right: an example of command generation algorithm. Rectangles represent primitives and ellipses represent data.

Figure 3 :

 3 Figure 3: Snapshot of the simulation environment, containing three bookshelves in a 6 m × 6 m closed room.

Final

  Pareto front, structure-restricted Final Pareto front, structure-free

Figure 4 :

 4 Figure 4: Final Pareto fronts of the solutions obtained with the two kinds of controllers. Those fronts are created by taking the best individuals in the final population in all the four islands and keeping only non-dominated ones. Individuals A and B are the examples used for the trajectories in the next figures.

Figure 5 :

 5 Figure 5:Trajectories of the robot with the structure-restricted controller A in the environment used during the evolution. The small dots are placed every 2 seconds and thus indicate the speed of the robot.

Figure 6 :

 6 Figure 6: Trajectories with the structure-free controller B in the evolution environment

Figure 7 :

 7 Figure 7: Pareto front of the evolved individuals in the generalization environment.

Figure 8 :

 8 Figure 8: Trajectories in the generalization environment with the structure-restricted controller A