V Flammang
email: valerie.flammang@univ-lorraine.fr

Bounds for the trace of small Salem numbers

come

Introduction

A Salem number is an algebraic integer τ > 1 of degree n = 2d ≥ 4 whose other conjugates, except for τ -1 have modulus 1. Now, the Mahler measure of an algebraic integer α with conjugates

α 1 = α, • • • , α n is M(α) = n j=1
max(1, |α j |) and its trace is tr(α)= n j=1 α j . Note that, if τ is a Salem number then M(τ ) = τ .

If α is an algebraic integer and M(α) = 1, then a classical theorem of Kronecker [K] tells us that α is a root of unity. It suggests the question: inf

α not a root of unity M(α) > 1 ? It is known as the Lehmer's problem and it is still open. Another formulation can be given as follows. Does there exist an absolute constant c > 0 such that: if M(α) > 1 then M(α) > 1 + c ? The smallest known value for M(α) is due to Lehmer himself [L] and is M(λ) = 1.176280 . . . where λ is the positive root > 1 of the polynomial L(z) = z 10 + z 9 -z 7 -z 6 -z 5 -z 4 -z 3 + z + 1. Here is the link between small Salem numbers (with Mahler measure less than 1.3) and Lehmer's problem: λ is a Salem number and moreover is up to now the smallest known Salem number. Currently, there are 47 known small Salem numbers (see the table [M2]). D.W. Boyd compiled 43 of such numbers [B1], [B2] in the years [START_REF] Boyd | Small Salem numbers[END_REF][START_REF] Boyd | Pisot and Salem numbers in intervals of the real line[END_REF]. Later, in 1990, M. J. Mossinghoff discovered the last four ones [M1], including one of degree 46. In 1999, the author, M. Grandcolas and G. Rhin proved that this table was complete up to degree 40 [F2]. M.J. Mossinghoff, G. Rhin and Q. Wu extended the result up to the degree 44 [MRW]. For more details on Salem numbers, see [S1].

In this paper, we focus our attention on the trace of small Salem numbers and prove that the trace can only take a small number of values. More precisely, we get:

Theorem 1. Let τ be a small Salem number of degree n = 2d ≥ 4. Then, we have:

-0.207188d ≤ tr(τ ) ≤ 0.193466d + 0.765083.

As the list of small Salem numbers is exhaustive up to the degree 44, this result is interesting from the degree 46. We have computed the bounds from this degree to the degree 30. We obtain:

Corollary 1. Let τ be a small Salem number of degree n = 2d ≥ 46. Then, we have:

d = 23 -4 ≤ tr(τ ) ≤ 5 d = 27 -5 ≤ tr(τ ) ≤ 5 d = 24 -5 ≤ tr(τ ) ≤ 5 d = 28 -5 ≤ tr(τ ) ≤ 6 d = 25 -5 ≤ tr(τ ) ≤ 5 d = 29 -6 ≤ tr(τ ) ≤ 6 d = 26 -5 ≤ tr(τ ) ≤ 5 d = 30 -6 ≤ tr(τ ) ≤ 6
Obviously, we have tested our bounds on the 47 known small Salem numbers and the results are recorded in Table 2 at the end of the paper. The method used to obtain the bounds in Theorem 1 also allows us to obtain the following result, which is less efficient than the one on the trace but which nevertheless provides some information on small Salem numbers:

Theorem 2. Let τ be a small Salem number of degree n = 2d ≥ 4 with minimal polynomial P . Then, we have:

|P (1)| ≤ exp(0.377596d -0.153515).
Corollary 2. Let τ be a small Salem number of degree n = 2d ≥ 46 with minimal polynomial P . Then, we have:

d = 23 |P (1)| ≤ 5070 d = 27 |P (1)| ≤ 22961 d = 24 |P (1)| ≤ 7396 d = 28 |P (1)| ≤ 33495 d = 25 |P (1)| ≤ 10790 d = 29 |P (1)| ≤ 48862 d = 26 |P (1)| ≤ 15740 d = 30 |P (1)| ≤ 71279
Again, we have tested our bounds on the 47 known small Salem numbers and the results are recorded in Table 3 at the end of the paper.

Relations between Salem numbers and totally positive algebraic integers

Let τ be a Salem number with minimal polynomial P (z) of degree n = 2d. Then P (z) is a reciprocal polynomial i.e., P (z) = z n P (1/z) The change of variable x = z + 1 z + 2 transforms the polynomial P (z) into a monic integer polynomial Q(x) of degree d with the following properties:

• To the roots τ and τ -1 of the polynomial P is associated a root α > 4 of the polynomial Q. Moreover, we are interesting only in small Salem numbers. It means τ < 1.3 and thus α < 4.069230.

• To the roots of modulus 1 of the polynomial P are associated d -1 roots in (0, 4) of the polynomial Q. Thus, all the conjugates of α are positive real numbers so α is called a totally positive algebraic integer.

• tr(α)=tr(τ ) + 2d.

In [F3], we proved that if α is a totally positive algebraic integer of degree d whose minimal polynomial is different from

x -1, x 2 -3x + 1, x 3 -5x 2 + 6x -1, x 4 -7x 3 + 13x 2 -7x + 1 and x 4 -7x 3 + 14x 2 -8x + 1 then tr(α) ≥ 1.792812d
. Hence, we immediately deduce the lower bound in Theorem 1: tr(τ ) ≥ 1.792812d -2d.

• P(1)=Q(4).

From now, we will restrict our attention to positive algebraic integers α of degree d with 4 < α < 4.069230 and all of whose other conjugates lie in (0, 4).

Proof of the Theorems

We detail the proof for Theorem 1 and we only indicate in the last below subsection what changes to get Theorem 2.

The principle of auxiliary function

It was introduced into number theory by C. J. Smyth in [S2]. The auxiliary function considered here is :

for x ∈ (0, 4.069230), f (x) = -x - 1≤j≤J c j log |Q j (x)| (1)
where the c j s are positive real numbers and the polynomials Q j s are nonzero polynomials in

Z Z[x].
Let m 1 be the minimum of the function f on (0, 4) and m 2 be the minimum of the function f on (4, 4.069230). If Q does not divide any Q j , then we have

d i=1 f (α i ) ≥ (d -1)m 1 + m 2 i.e., -tr(α) ≥ (d -1)m 1 + m 2 + 1≤j≤J c j log | d i=1 Q j (α i )|. Since Q does not divide any Q j then d i=1 Q j (α i ) is a nonzero integer because it is the resultant of Q and Q j .
Hence, if α is not a root of Q j , we have

-(d -1)m 1 -m 2 ≥ tr(α).

Link between the integer transfinite diameter and the auxiliary functions

Let K be a compact subset of C and ϕ be a positive function defined on K. The ϕ integer transfinite diameter of K is defined by

t Z,ϕ (K) = lim inf inf sup |H(z)| 1 n ϕ(z) . n ≥ 1 H ∈ Z[X] z ∈ K n → ∞ deg(H) = n
This weighted version of the integer transfinite diameter was introduced by F. Amoroso [A].

In the auxiliary function (1), we replace the coefficients c j by rational numbers a j /q where q is a positive integer such that q.c j is an integer for all 1 ≤ j ≤ J. Then we can write:

for x ∈ (0, 4.069230), f (x) = -x - t r log |H(x)| ≥ (d -1)m 1 + m 2 (2)
where

H = J j=1 Q a j j ∈ Z Z[X] is of degree r = J j=1 a j deg Q j and t = J j=1 c j deg Q j (
|H(x)| t/r e x ≤ e -(d-1)m 1 -m 2 .
If we suppose that t is fixed, it is equivalent to find an effective upper bound for the weighted integer transfinite diameter over the interval (0, 4.069230) with the weight ϕ(x) = e x :

t Z Z,ϕ ((0, 4.069230)) = lim inf inf sup |H(x)| t r ϕ(x) r ≥ 1 H ∈ Z Z[X] x ∈ (0, 4.069230) r → ∞ deg(H) = r

Construction of an explicit auxiliary function

The main point is to find a set of "good "polynomials Q j , i.e., which gives the best possible value for m 1 and m 2 . Until 2003, the polynomials were found heuristically. In 2003, Q. Wu [Wu] developed an algorithm that allows a systematic search of "good "polynomials. In 2009 [F1], we made two improvements to this previous algorithm in the use of the LLL algorithm. The idea is to get good polynomials Q j by induction. Thus, we call this algorithm the recursive algorithm. For more details, see [F3]. Suppose that we have some polynomials Q 1 , Q 2 , . . . , Q J . Then we use the semi-infinite linear programming ( introduced into number theory by C. Smyth [S2]) to optimize f for this set of polynomials (i.e., to get the greatest possible m 1 and m 2 ). We obtain the numbers c 1 , c 2 , . . . , c J and f in the form (2) with t = x(r + k) t to be as small as possible. We apply the LLL algorithm to the linear forms

Q(x i )R(x i ) exp x i (r + k) t .
The x i are control points which are points uniformly distributed on the interval I = (0, 4.069230) to which we have added points where f has local minima. Thus we find a polynomial R whose irreducible factors R j are good candidates to enlarge the set {Q 1 , . . . , Q J }. We only keep the factors R j that have a nonzero coefficient in the newly optimized auxiliary function f .After optimization, some previous polynomials Q j may have a zero exponent and so are removed.

To obtain the upper bound in Theorem 1, k has varying from 3 to 7. The polynomials and their coefficient can be read off from Table 4.

Proof of Theorem 2

The only change is the type of the auxiliary function. In the proof of this theorem, the function involved is:

for x ∈ (0, 4.069230), f (x) = -log(x + 4) - 1≤j≤J c j log |Q j (x)|
where, as previously, the c j s are positive real numbers and the polynomials Q j s are nonzero polynomials in Z Z[x]. Then, the rest of the proof follows that of Theorem 1. The polynomials and their coefficients can be read off from Table 5.

Remarks on the computations

It is obvious that to be able to do the computations we had to set a value of d. As the table is complete up to degree 44, we started with d = 23. This gave us values for m 1 and m 2 and the corresponding upper bound for tr(τ ). With these values of m 1 and m 2 , we then calculated the different other bounds for d = 24 to d = 30 and this gave the results recorded in Corollary 1. Moreover, to be sure of having optimum results, we have redone each calculation entirely for each value of d, from 24 to 30. As the bounds obtained are not better than the previous ones, we preferred to present the results under the form of Theorem 1, rather than only in the form of Corollary 1.

Table 1: The 47 known small Salem numbers

N°Polynomial P 1 z 8 -z 5 -z 4 -z 3 + 1 2 z 10 -z 6 -z 5 -z 4 + 1 3 z 10 -z 7 -z 5 -z 3 + 1 4 z 10 -z 8 -z 5 -z 2 + 1 5 z 10 -z 8 -z 7 + z 5 -z 3 -z 2 + 1 6 z 10 + z 9 -z 7 -z 6 -z 5 -z 4 -z 3 + z + 1 7 z 12 -z 11 + z 10 -z 9 -z 6 -z 3 + z 2 -z + 1 8 z 14 -z 12 -z 7 -z 2 + 1 9 z 14 -z 13 -z 8 + z 7 -z 6 -z + 1 10 z 14 -z 11 -z 10 + z 7 -z 4 -z 3 + 1 11 z 14 -z 12 -z 11 + z 9 -z 7 + z 5 -z 3 -z 2 + 1 12 z 16 -z 15 -z 8 -z + 1 13 z 18 -z 17 -z 10 + z 9 -z 8 -z + 1 14 z 18 -z 17 -z 14 + z 13 -z 9 + z 5 -z 4 -z + 1 15 z 18 -z 12 -z 11 -z 10 -z 9 -z 8 -z 7 -z 6 + 1 16 z 18 -z 14 -z 12 -z 11 -z 9 -z 7 -z 6 -z 4 + 1 17 z 18 -z 17 -z 14 + z 13 -z 12 + z 10 -z 9 + z 8 -z 6 + z 5 -z 4 -z + 1 18 z 18 -z 17 + z 16 -z 15 -z 12 + z 11 -z 10 + z 9 -z 8 + z 7 -z 6 -z 3 + z 2 -z + 1 19
z 18 -2z 17 + 2z 16 -2z 15 + 2z 14 -2z 13 + 2z 12 -3z 11 + 3z 10 -3z 9 + 3z 8 -3z 7 + 2z 6 -2z 5 + 2z

4 -2z 3 + 2z 2 -2z + 1 20 z 20 -z 18 -z 15 -z 5 -z 2 + 1 21 z 20 -z 18 -z 15 -z 12 + z 10 -z 8 -z 5 -z 2 + 1 22 z 20 -z 19 -z 15 + z 14 -z 11 + z 10 -z 9 + z 6 -z 5 -z + 1 23 z 20 -2z 19 + 2z 18 -2z 17 + 2z 16 -2z 15 + z 14 -z 12 + z 11 -z 10 + z 9 -z 8 + z 6 -2z 5 + 2z 4 -2z 3 + 2z 2 -2z +1 24 z 22 -z 21 -z 17 + z 11 -z 5 -z + 1 25 z 22 -z 21 -z 20 + z 19 -z 13 + z 11 -z 9 + z 3 -z 2 -z + 1 26 z 22 -z 20 -z 19 + z 15 + z 14 -z 12 -z 11 -z 10 + z 8 + z 7 -z 3 -z 2 + 1 27 z 22 -z 21 -z 19 + z 18 -z 14 + z 13 -z 12 + z 11 -z 10 + z 9 -z 8 + z 4 -z 3 -z + 1 28 z 24 -z 23 -z 18 -z 6 -z + 1 29 z 24 -z 23 -z 20 + z 19 -z 17 + z 16 -z 15 + z 13 -z 12 + z 11 -z 9 + z 8 -z 7 + z 5 -z 4 -z + 1 30 z 26 -z 25 -z 20 + z 13 -z 6 -z + 1 31
z 26 -z 24 -z 21 -z 18 + z 16 + z 13 + z 10 -z 8 -z 5 -z 2 + 1 32 z 26 -z 24 -z 23 + z 19 -z 17 -z 16 + z 14 + z 13 + z 12 -z 10 -z 9 + z 7 -z 3 -z 2 + 1 33 z 26 -z 20 -z 19 -z 18 -z 17 -z 16 -z 15 -z 14 -z 13 -z 12 -z 11 -z 10 -z 9 -z 8 -z 7 -z 6 + 1 34 z 26 -z 25 -z 22 + z 21 -z 20 + z 18 -z 17 + z 16 -z 14 + z 13 -z 12 + z 10 -z 9 + z 8 -z 6 + z 5 -z 4 -z + 1 35 z 26 -2z 25 + z 24 + z 23 -2z 22 + z 21 -z 18 + z 17 -z 15 + z 14 -z 13 + z 12 -z 11 + z 9 -z 8 + z 5 -2z 4 + z 3 + z 2 -2z + 1 36 z 26 -z 25 -z 24 + 2z 22 -2z 20 -z 19 + 2z 18 + 2z 17 -2z 16 -2z 15 + 3z 13 -2z 11 -2z 10 + 2z 9 + 2z 8 -z 7 -2z 6 +2z 4 -z 2 -z + 1 37 z 28 -z 24 -z 23 -z 22 -z 21 -z 20 + z 16 + z 15 + z 14 + z 13 + z 12 -z 8 -z 7 -z 6 -z 5 -z 4 + 1 38 z 30 -z 29 -z 22 -z 18 -z 15 -z 12 -z 8 -z + 1 39 z 30 -z 28 -z 25 -z 24 + z 20 + z 17 -z 15 + z 13 + z 10 -z 6 -z 5 -z 2 + 1 40 z 30 -z 25 -z 24 -z 23 -z 22 -z 21 -z 20 + z 15 -z 10 -z 9 -z 8 -z 7 -z 6 -z 5 + 1 41 z 30 -2z 29 + 2z 28 -2z 27 + z 26 -z 24 + 2z 23 -2z 22 + z 21 -z 19 + z 18 -z 17 + z 16 -z 15 + z 14 -z 13 + z 12 -z 11 +z 9 -2z 8 + 2z 7 -z 6 + z 4 -2z 3 + 2z 2 -2z + 1 42 z 34 -z 33 -z 30 + z 29 -z 28 + z 26 -z 25 + z 24 -z 22 + z 21 -z 20 + z 18 -z 17 + z 16 -z 14 + z 13 -z 12 + z 10 -z 9 +z 8 - 0.0334017400 x 2 -5x + 5 0.1931207414 x 3 -6x 2 + 9x -1 0.0288203880 x 3 -6x 2 + 9x -3 0.0185605107 x 3 -7x 2 + 14x -7 0.1095768403 x 4 -9x 3 + 26x 2 -24x + 1 0.0514099247 x 5 -11x 4 + 43x 3 -69x 2 + 36x + 1 0.0551475719 x 5 -11x 4 + 44x 3 -77x 2 + 55x -11 0.0324052775 x 6 -14x 5 + 76x 4 -199x 3 + 252x 2 -132x + 17 0.0030224114 2x 7 -34x 6 + 236x 5 -856x 4 + 1719x 3 -1854x 2 + 945x -163 0.0058282451 0.0413670035 x 3 -5x 2 + 6x -1 0.1989461843 x 3 -6x 2 + 9x -1 0.1138658462 x 5 -9x 4 + 28x 3 -35x 2 + 15x -1 0.1023209293 x 6 -11x 5 + 45x 4 -84x 3 + 70x 2 -21x + 1 0.0544779065

z 6 + z 5 -z 4 -z + 1 43 z 34 -z 33 -z 31 + z 29 + z 27 -2z 26 + z 23 + z 22 -z 21 -z 20 -z 19 + z 18 + z 17 + z 16 -z 15 -z 14 -z 13 + z 12 +z 11 -2z 8 + z 7 + z 5 -z 3 -z + 1 44 z 36 + z 35 -z 33 -2z 32 -2z 31 -z 30 + z 28 + z 27 -z 25 -z 24 + z 22 + z 21 -z 19 -z 18 -z 17 + z 15 + z 14 -z 12 -z 11 + z 9 + z 8 -z 6 -2z 5 -2z 4 -z 3 + z + 1 45 z 40 -z 37 -z 35 -z 33 -z 31 -z 29 + z 26 + z 24 + z 22 + z 20 + z 18 + z 16 + z 14 -z 11 -z 9 -z 7 -z 5 -z 3 + 1 46 z 44 -z 43 -z 37 -z 33 + z 25 + z 22 + z 19 -z 11 -z 7 -z + 1 47 z 46 -z 42 -z 41 -z 40 -z 39 + z 25 + z 24 + z 23 + z 22 + z 21 -z 7 -z 6 -z 5 -z 4 + 1

c

  j deg(Q j ). We seek a polynomial R ∈ Z Z[x] of degree k (k = 3 for instance) such that sup x∈I |H(x)R(x)| t r+k e x ≤ e -(d-1)m 1 -m 2 where H = J j=1 Q j and I is an interval composed with points uniformly distributed on (0, 4.069230). We want the quantity sup x∈I |H(x)R(x)| exp

  this formulation was introduced by J. P. Serre). Thus we seek a polynomial H ∈ Z Z[X] such that

	sup
	x∈(0,4.069230)

Table 2 :

 2 The bounds of Theorem 1 for the 47 known small Salem numbers

	N°deg(P ) Lower bound tr(P ) Upper bound
	1	8	-0.82875200	0	1.5389470
	2	10	-1.0359400	0	1.7324130
	3	10	-1.0359400	0	1.7324130
	4	10	-1.0359400	0	1.7324130
	5	10	-1.0359400	0	1.7324130
	6	10	-1.0359400	-1	1.7324130
	7	12	-1.2431280	1	1.9258790
	8	14	-1.4503160	0	2.1193450
	9	14	-1.4503160	1	2.1193450
	10	14	-1.4503160	0	2.1193450
	11	14	-1.4503160	0	2.1193450
	12	16	-1.6575040	1	2.3128110
	13	18	-1.8646920	1	2.5062770
	14	18	-1.8646920	1	2.5062770
	15	18	-1.8646920	0	2.5062770
	16	18	-1.8646920	0	2.5062770
	17	18	-1.8646920	1	2.5062770
	18	18	-1.8646920	1	2.5062770
	19	18	-1.8646920	2	2.5062770
	20	20	-2.0718800	0	2.6997430
	21	20	-2.0718800	0	2.6997430
	22	20	-2.0718800	1	2.6997430
	23	20	-2.0718800	2	2.6997430
	24	22	-2.2790680	1	2.8932090
	25	22	-2.2790680	1	2.8932090
	26	22	-2.2790680	0	2.8932090
	27	22	-2.2790680	1	2.8932090
	28	24	-2.4862560	1	3.0866750
	29	24	-2.4862560	1	3.0866750
	30	26	-2.6934440	1	3.2801410
	31	26	-2.6934440	0	3.2801410
	32	26	-2.6934440	0	3.2801410
	33	26	-2.6934440	0	3.2801410
	34	26	-2.6934440	1	3.2801410
	35	26	-2.6934440	2	3.2801410
	36	26	-2.6934440	1	3.2801410
	37	28	-2.9006320	0	3.4736070
	38	30	-3.1078200	1	3.6670730
	39	30	-3.1078200	0	3.6670730
	40	30	-3.1078200	0	3.6670730
	41	30	-3.1078200	2	3.6670730
	42	34	-3.5221960	1	4.0540050
	43	34	-3.5221960	1	4.0540050
	44	36	-3.7293840	-1	4.2474710
	45	40	-4.1437600	0	4.6344030
	46	44	-4.5581360	1	5.0213350
	47	46	-4.7653240	0	5.2148010

Table 3 :

 3 The bound of Theorem 2 for the 47 known small Salem numbers

	N°deg(P )	|P (1)| Upper bound
	1	8	1	3.8840134
	2	10	1	5.6658957
	3	10	1	5.6658957
	4	10	1	5.6658957
	5	10	1	5.6658957
	6	10	1	5.6658957
	7	12	1	8.2652584
	8	14	1	12.057140
	9	14	1	12.057140
	10	14	1	12.057140
	11	14	1	12.057140
	12	16	1	17.588636
	13	18	1	25.657835

Table 4 :

 4 Polynomials and their coefficients used for the upper bound of Theorem 1

	Polynomials	Coefficients
	x	0.0019679665
	x -2	0.1898490761
	x -3	0.5308931151
	x -4	0.5815137153
	x 2 -4x + 1	0.0725582571
	x 2 -4x + 2

Table 5 :

 5 Polynomials and their coefficients used for the upper bound of Theorem 2

	Polynomials	Coefficients
	x	0.3620883142
	x -1	0.5053687655
	x -2	0.3014609143
	x -3	0.1287326889
	x 2 -3x + 1	0.2940495893
	x 2 -4x + 1	0.0027131692
	x 2 -4x + 2	0.1046760710
	x 2 -5x + 5


Antoine Vacavant
email: antoine.vacavant@liris.cnrs.fr

David Coeurjolly
email: david.coeurjolly@liris.cnrs.fr

Laure Tougne
email: laure.tougne@liris.cnrs.fr

Topological and Geometrical Reconstruction of Complex Objects on Irregular Isothetic Grids

HAL is

Introduction

The character and symbol representation, description and classification are necessary tasks in many current applications, and concern both research and industrial challenges. Those tasks are applied on images generally designed within a regular grid, i.e. all the pixels have the same size, and their position can be easily indexed. However, it is now common to successively divide an image into subimages, as in quadtree decomposition [START_REF] Samet | The quadtree and related hierarchical data structures[END_REF][START_REF] Samet | Hierarchical spatial data structures. In Design and Implementation of Large Spatial Databases[END_REF], to represent a part of an image in a more compact and adapted structure. These techniques describe interesting parts of an image, from different points of view, through a set of irregular pixels. In this paper, we introduce the concept of shape representation within an irregular isothetic grid (I-grid for short) [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF]. The pixels are defined by variable sizes and positions, and may be determined by subdivision rules. We propose to represent the topology of the elements contained in the irregular twodimensional (2-D) image by constructing their associated Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complément intégrable ou d'une fonction numérique[END_REF], then we represent them by a simple polygonal structure that respects the extended supercover digitization model defined in [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF]. This structure also preserves the topology that we reveal in the previous stage. We clearly address the problem of vectorization (or raster-to-vector) on irregular isothetic grids, and not only in the scope of document analysis. In our framework, we are interested in binary images containing irregular objects, i.e. k-objects in respect to the definition given in [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF], where k represents the considered relation of adjacency (see Section 2 for further details). Those complex objects may contain holes, and could represent characters, symbols, lines, etc. An application of such binary image processing is clearly document and line drawings analysis, but we can also consider a discrete subdivision of a part of R 2 representing the solutions of a given function f : R 2 → R. The algorithms designed in interval arithmetic are interesting approaches to address those problems [START_REF] De Figueiredo | Fast interval branch-and-bound methods for unconstrained global optimization with affine arithmetic[END_REF][START_REF] Kearfott | Interval computations: introduction, uses, and resources[END_REF][START_REF] Snyder | Interval analysis for computer graphics[END_REF].

The techniques of vectorization developed until now on the discrete regular domain can be divided into several classes, up to the final application of the method [START_REF] Cordella | Symbol recognition in documents: a collection of techniques ?[END_REF][START_REF] Hilaire | Robust and accurate vectorization of line drawings[END_REF][START_REF] Mertzios | On applying fast and efficient methods in pattern recognition[END_REF][START_REF] Wenyin | From raster to vectors: extracting visual information from line drawings[END_REF]. We will only focus on a few kinds of raster-to-vector methodologies, largely developped for document analysis applications. To our knowledge, there exists no generic extension of those approaches on irregular isothetic grids. The run length encoding (RLE) based methods first build a decomposition into elongated cells along an axis of the image where we can build a line adjacency graph (LAG) [START_REF] Burge | A minimal line property preserving representation of lines images[END_REF][START_REF] Elgammal | Graph-based segmentation and feature-extraction framework for arabic text recognition[END_REF]. Those methods aim to describe the topology of the encountered objects in the image, but the geometrical structure deduced from it has to be improved by many post-treatment processings. The skeletonization and thinning methods are surely the most widely employed methods in vectorization. We can notice that tools designed in mathematical morphology [START_REF] Soille | Morphological image analysis[END_REF] are a frequent choice to prepare the images before processing the skeletonization. A survey of vectorization methods based on skeleton can be found in [START_REF] Lam | Thinning methodologies -a comprehensive survey[END_REF], and another one about such techniques not using it in [START_REF] Wenyin | A Survey of non-thinning based vectorization methods[END_REF]. The aim is to compute a medial axis of the object that minimally represents its shape [START_REF] Klette | Digital geometry[END_REF]. However, those techniques modify the original geometry of the object to obtain a minimal representation of it. Besides, they need filtering or smoothing pre-treatment processings to reduce the noise that could pertubate the final medial axis. The k-object can contain holes, and so may be composed by thick arcs. In the work of Debled et al. [START_REF] Debled | A linear algorithm for segmentation of digital curves[END_REF][START_REF] Debled | Fast polygonal approximation of digital curves[END_REF][START_REF] Debled | Optimal blurred segments decomposition in linear time[END_REF], the definitions of discrete lines and blurred segments join the concept of thick regular arcs. But, beyond this geometrical representation of arcs, the global structure is not aborded, and thus there are no description of the topology of the recognized objects.

In this article, we first introduce the concepts of k-arcs and k-objects by recalling some definitions, then we present the extended supercover model on an I-grid. We also recall the invertible reconstruction of k-arcs described in [START_REF] Coeurjolly | Supercover model, digital straight line recognition and curve reconstruction on the irregular isothetic grids[END_REF]. In the third part, we give details about the two main phases of our system: the description of the topology of a complex object based on the Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complément intégrable ou d'une fonction numérique[END_REF], and its polygonal reconstruction. Then, we present some experiments and revealing results to illustrate the two phases of our algorithm. We also proove the robustness of the polygonal reconstruction by a test on a large image of technical drawing. We finally discuss the applications of our contribution, and the improvement on its global performance.

We first define an irregular isothetic grid, denoted I, as a tiling of the plane with isothetic rectangles. We shortly recall that each rectangle P (also called cell) of I is defined by its center (x P , y P ) ∈ R 2 and a size (l x P , l y P ) ∈ R 2 . The position and the size of P may be controlled by different level of constraints; e.g. in the case of quadtree decomposition [START_REF] Samet | The quadtree and related hierarchical data structures[END_REF][START_REF] Samet | Hierarchical spatial data structures. In Design and Implementation of Large Spatial Databases[END_REF], for a cell of level k, (x P , y P ) = ( m 2 k , n 2 k ) and l x P = l y P = 1 2 k-1 for some m, n ∈ Z [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF][START_REF] Coeurjolly | Supercover model, digital straight line recognition and curve reconstruction on the irregular isothetic grids[END_REF]. In our framework, adjacency relation is an important feature that we depict through the following definitions.

Definition 1 (ve-adjacency and e-adjacency). Let P and Q be two cells. P and Q are ve-adjacent (vertex and edge adjacent) if :

or    |x P -x Q | = l x P +l x Q 2 and |y P -y Q | ≤ l y P +l y Q 2 |y P -y Q | = l y P +l y Q 2 and |x P -x Q | ≤ l x P +l x Q 2
P and Q are e-adjacent (edge adjacent) if we consider an exclusive "or" and strict inequalities in the above ve-adjacency definition. k may be interpreted as e or ve in the following definitions.

Definition 2 (k-arc). Let E be a set of cells, E is a k-arc if and only if for each element of E = {P i , i ∈ {1, ..., n}}), P i has exactly two k-adjacent cells, except P 1 and P n which are called extremities of the k-arc.

Definition 3 (k-object). Let E be a set of cells, E is a k-object if and only if for each couple of cells (P, Q) belonging to E × E, there exists a k-path between P and Q in E.

We now consider the extension of the supercover model from [START_REF] Cohen-Or | Fundamentals of surface voxelization[END_REF] on irregular isothetic grids [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF] to digitize Euclidean objects on I.

Definition 4 (Supercover on irregular isothetic grids). Let F be an Euclidean object in R 2 . The supercover S(F ) is defined on an irregular isothetic grid I by:

S(F ) = P ∈ I | B ∞ (P ) ∩ F = ∅ = P ∈ I | ∃(x, y) ∈ F, |x P -x| ≤ l x P 2 and |y P -y| ≤ l y P 2
where B ∞ (P ) is the rectangle centered in (x P , y P ) of size (l x P , l y P ) (if

l x P = l y P , B ∞ (P ) is the ball centered in (x P , y P ) of size l x P for the L ∞ norm).
This model has several interesting properties, e.g. for F , G two Euclidean objects in R 2 , we have [START_REF] Coeurjolly | Supercover model and digital straight line recognition on irregular Isothetic Grids[END_REF] for more details).

S(F ∪G) = S(F )∪S(G) or S(F ∩G) ⊆ S(F )∩S(G) (see proposition 2 in
We now present the k-arc reconstruction algorithm we use in our complex object geometrical representation phase (Section 3.2). Moreover, this approach respects the supercover model we have just presented. The algorithm proposed in [START_REF] Coeurjolly | Supercover model, digital straight line recognition and curve reconstruction on the irregular isothetic grids[END_REF] to decompose a curve into segments is first based on the following definition of an irregular digital line.

Definition 5 (Irregular isothetic digital straight line). Let S be a set of cells in I, S is called a piece of irregular digital straight line (IDSL for short) iff there exists an Euclidean straight line l such that:

S ⊆ S(l)
In other words, S is a piece of IDSL iff there exists l such that for all P ∈ S,

B ∞ (P ) ∩ l = ∅.
The algorithm inspired from [START_REF] Sivignon | Discrete analytical curve reconstruction without patches[END_REF] principally uses the construction and update procedures of a visibility cone, and can be sketched as follows. We first fix the extremity p 0 of the first segment such that p 0 ∈ P 0 . We note e 0 the Euclidean segment shared by P 0 and P 1 , and we consider the first cone C 0 (p 0 , s, t) such that s and t coincide with the extremities of e 0 and {p 0 , s, t} is sorted counterclockwise. Then, for each cell P i , we consider the shared segment e i between P i-1 and P i , and the current cone C j (p j , s, t) is updated. When the update procedure fails, a new cone C j+1 (p j+1 , s, t) is set up, and we add the point p j+1 to the reconstruction: to compute the new cone, authors of [START_REF] Coeurjolly | Supercover model, digital straight line recognition and curve reconstruction on the irregular isothetic grids[END_REF] consider the bisector of the cone and define p j+1 as the midpoint of the intersection between the bisector and the pixel P i-1 . The Figure 1 illustrates the progressive construction of cones in a k-arc, and the resulting segmentation into lines.

Complex Objects Definition and Representation on Irregular Isothetic Grids

In this section, we present the two main phases of our system for object representation on irregular isothetic grids.

Representation of a Complex Topology

To represent the shape of a k-object E, we have chosen an incremental directional approach to build its associated Reeb graph G, as in continuous space (see Figure 2). It is an interesting structure introduced by G. Reeb [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complément intégrable ou d'une fonction numérique[END_REF] based on the Morse theory [START_REF] Gramain | Topologie des surfaces[END_REF][START_REF] Hart | Computational topology for shape modeling[END_REF]. This graph is also used in many applications for surface and curve description [START_REF] Hétroy | Méthodes de partitionnement de surfaces[END_REF][START_REF] Tung | Indexation 3D de bases de données d'objets 3D par graphes de Reeb améliorés[END_REF][START_REF] Xiao | A discrete Reeb graph approach for the segmentation of human body scans[END_REF]. The Reeb graph G is associated to a height function f defined on E, and nodes of G represent the critical points of f . Moreover, to have a minimal representation of the topological information of E, each edge of the Reeb graph corresponds to a k-arc. Those k-arcs will be segmented in the stage of polygonal description of E (Section 3.2). We denote the left, right, top and bottom borders of a cell P respectively P L , P R , P T and P B . We have, for example, the abscissa of P L equal to x P -(l X P /2) (that we denote P L = x P -(l X P /2)). We also abusively say that a k-arc A and a cell P are k-adjacent if there exists a cell Q in A such that P and Q are k-adjacent. Let E = {P i } i=1,...,n be a given 2-D set of cells. We first choose a direction to treat the cells of E. Without loss of generality, we can suppose that we choose the left-to-right orientation above X axis, i.e. the height function f is defined along X axis. At time t = 0, we merge together all the k-adjacent cells P of E with the smallest left border x t=0 = x 0 , e.g. P L = x 0 = 0. This merging task is processed by the update procedure described below. Those m collections of cells define the begin cells of the initial recognized k-arcs A 1 , A 2 , ..., A m . Update Procedure. Let A be a k-arc, and P 1 and P 2 two adjacent cells of E such that P 1 ∈ A, P L 1 < P L 2 , and P 2 should be added to A. If P L 2 = P R 1 , we just add P 2 to A, else the procedure updates the k-arc A with P 2 , and may recode A. For that, we first build the greatest common rectangle F 2 of P 1 and P 2 . Definition 6 (Greatest common rectangle). Let P 1 and P 2 be two adjacent rectangles. F 2 is the greatest common rectangle (or GCR) of P 1 and P 2 iff i)

F 2 ⊆ P 1 ∪ P 2 , ii) F 2 ∩ P 1 = ∅ , iii) F 2 ∩ P 2 = ∅ , iv
) there is no rectangle greater than F 2 by inclusion respecting i), ii) and iii). Then, we consider the rectangles P 1 -F 2 and P 2 -F 2 . If P R 1 < P R 2 , we denote P 1 -F 2 = F 1 and P 2 -F 2 = F 3 , else we prefer P 1 -F 2 = {F 1 , F 3 }. We can notice that those rectangles may be empty, e.g. 3 presents the five general configurations of update procedure (there are also five other configurations, obtained by symetry when P T 2 > P T 1 ), and the k-arc recoding that we have to consider. Besides, we propose to reduce the number of cells in A by joining the two rectangles F 1 and

F 3 = ∅ if P R 1 = P R 2 , since in that case F L 3 = F R 3 . Figure
F 3 if F T 1 = F T 3 , F B 1 = F B 3 and F 2 = ∅.
This junction is processed by replacing F 1 and F 3 by the rectangle F 1 ∪ F 3 . Finally, the procedure ends by removing P 1 from A, and by adding the cells corresponding to the rectangles F 1 and F 2 to A. F 3 is also pushed in E, and will be treated later; more exactely at time t such that x t = F L 3 . At time t + 1, our algorithm consists first in merging the adjacent cells with the same left border x t+1 in k cells C 1 , C 2 , ..., C k (see update procedure for details). Those candidate cells may be added to one or more k-arcs among A i , i ∈ {1, ..., m} if they are adjacent to A i . It is clear that only a cell Q built at time t and having its right border Q R equal to x t+1 may be adjacent with a cell C j , j ∈ {1, ..., k}. A cell C j can be treated by several manners: -C j is not adjacent with any k-arc A i . We initialize a new k-arc A m+1 with the cell C j . C j represents the begin cell of A m+1 . -If C j is adjacent with one k-arc A i , then we just update A i with C j .

-When C j is k-adjacent with p k-arcs A i , A i+1 , ..., Ai + p, it is a merge phase. First, we update each k-arc with C j . The cell C j is marked as a merge cell and indicates that each k-arc A i , ..., A i+p has a k-arc A m+1 = {C j } linked as a next arc. -The case where p cells C j , C j+1 , ..., C j+p are k-adjacent with an k-arc A i is called a split phase. We first update A i with C j by the update procedure.

Then we denote Q the cell in A i such that Q R = x t+1 . We also define p new next k-arcs A m+1 , ..., A m+p of A i such that A m+1 = {Q, C j }, ..., A m+p = {Q, C j+p }. In those p k-arcs and in A i , Q is marked as a split cell.

When the algorithm ends, at time t such that x t is the greatest left border in E, we define the last added cell in every k-arc A i as an end cell. In this stage of our algorithm, there may also appear a split phase and a merge phase for a cell C j . We do not detail this specific case but it can be easily handled. We depict in Figure 4 the progressive construction of the graph and the recoding of the k-object presented in Figure 2 (b) in five stages of the algorithm.

(deg(n)) - n∈V,(n=s)∨(n=m) (deg(n) -2)
where deg(n) is the degree of the node n in G, so deg(n) = 1 if n is a begin or end node. The Euler number permits to describe the topology of an object by an unique value. For example, for a torus, χ = 0, for a disc, χ = 2, and the object described in Figure 2 (b) has a Euler number χ = -4; we can also say that this shape is homeomorpheous to a torus with 3 holes where χ = 2 -2 × #(holes) = -4. Beside the topological invariants obtained by critical points, the structure of the graph clearly depends on the direction we choose for the height function f . A part of the nodes and the edges may change, but the information on the topology of E, i.e. internal nodes of G, is not modified. The Euler number is an example of the use of the Reeb graph for shape description. Let us consider now E ′ as the object drawn in the fourth image of Figure 4. The three cells added during the last iteration could be noise modifying the contour of E ′ . The Reeb graph is modified by a split phase, three nodes are created, whereas these cells are maybe noise. Actually, the problem of the perturbation of the contour of E ′ could be certainly reduced if the object was first filtered or smoothed. This kind of pre-treatment processings is often adopted, whatever the approach we may choose for shape representation, e.g. skeletonization. Finally, with the update procedure, we recode the cells in E so that a k-arc is always represented between two nodes of G. This geometrical rearrangement clearly depends on the direction of f , but does not change neither the topology nor the contour of the recognized k-arcs. The topological structure so described is simple, and prepares the next phase of our complex objects reconstruction system.

Polygonal Reconstruction of Thick Objects

Since the reconstruction into polylines always affects the first point p 0 as the center of the first treated cell, we propose to start the reconstruction of every k-arcs computed in the previous stage by the merge and split nodes detected in the Reeb graph G. This insures that each of those particular nodes of G will be represented by an unique point in the final polygonalization. The segments are recognized from intersections between several parts of the object E to its extremities, i.e. we consider the edges me, se, mb and sb of G. Moreover, since the recognition algorithm is greedy, the possible error induced by the visibility cone approach is propagated to the extremities of E, instead of those intersections that represent the shape of the object. For the ms, mm, ss and sm configurations of edges in G we propose to process a bidirectional reconstruction that begins from each node of the edge, and ends in its center. Thus, the error may be concentrated in the midpoint of those edges. This approach confirms that the nodes m and s of G represent the places of an object where the description of its geometry must be precise. Finally, we choose to treat edges b-e by the same bidirectional reconstruction, that seems to be the more efficient way to insure a robust reconstruction. We do not deal about the problem of linking the two reconstructions on the k-arc (reconstruction with patch), because an efficient and general joint technique between two discrete lines implies that our algorithm would not be linear anymore [START_REF] Breton | Reconstruction inversible d'objets discrets 2D[END_REF]. Hence we just add a segment between the two polylines. This phase of our system can not be handled without patch, since we use the internal points of the shape of E to guide the geometrical reconstruction.

In Figure 5 b, we illustrate the behaviour of our algorithm in the case of the object E presented in the previous section. We also show the interest of our approach for a symetrical complex object. This structure respects the supercover model, and the symetrical shape of this object. We also show the result of our algorithm on the k-object presented in the previous section (b)

Contrary to conventional vectorization methods, we propose a technique that respects the supercover model on an I-grid. We do not address the quality of the global polygonal structure deduced from this second phase of our system. To introduce the concept of quality in the framework of document analysis, we may refer to [START_REF] Wenyin | A protocol for performance evaluation of line detection algorithms[END_REF].

Experimentation and Results

In Figure 6, we present the polygonal structure obtained on an image first rearranged by a quadtree-based approach. The reconstruction of k-arcs stands inside the object, and the split and merge nodes are represented with one point in the reconstruction. The polygonal representation also permits to measure geometrical features (e.g. length) of a complex function f : R 2 → R (Figure 7). f is first discretized by an interval computing algorithm through a set of cells E, then we use our system to minimally describe the curves of E. Finally, to show the robustness of our system, we present in Figure 8 the polygonal and topological reconstructions of a large image of technical drawing. The representation by lines of an object described on a binary image is a classical problem often considered in the framework of document analysis. We have proposed to enlarge the scope of vectorization methodologies to irregular isothetic representation of binary data. Depending on the final application of our system, we can treat the initial image with pre-treatment processings, reorganize the Reeb graph (edge contraction, etc.), or rearrange the segments finally processed in the second phase. The geometrical reconstruction stands inside the object, i.e. it respects the irregular digitization supercover model. Moreover, this reconstruction preserves the topology described by the Reeb graph. Thus, our system is robust, and topologically and geometrically correct. The Reeb graph can be extended to three-dimensional (3-D) object description, with a similar incremental approach. However, visibity cone reconstruction is hardly adaptable to such irregular objects. Our system should be modified to provide a 3-D polygonalization based on the Reeb graph. Such technique would be convenient especially for medical imaging, e.g. organ representation in an irregular 3-D CT-scan image.

Fig. 1 .

 1 Fig. 1. An example of the progressive construction of cones in a k-arc (left), and the reconstruction into segments we obtain (right)

Fig. 2 .

 2 Fig. 2. (a): an example of the Reeb graph G of a continuous object E. The nodes of G represent the critical points of f (maxima, minima, inflection points), and an edge is a connected component of E between two critical points. (b): an example of an irregular object E (left), the final recoded structure with k-arcs (right) and the Reeb graph associated to the height function f defined on E (bottom). The notations b, e, m and s are given at the end of this section

Fig. 3 .

 3 Fig. 3. Description of rectangles F1, F2 and F3 in the update procedure (top), and the associated cells as result (bottom). When P R 1 < P R 2 (a and b), P1 -F2 = F1 and P2 -F2 = F3, else P1 -F2 = {F1, F3} (d and e). If P R 1 = P L 2 , F2 = ∅, when P R 1 = P R 2 , F3 = ∅ and finally F1 = ∅ in the case P L 1 = P L 2

Fig. 4 .

 4 Fig.4. The recognized k-arcs and the associated Reeb graph for some iterations of our algorithm on the object presented in Figure2 (b). First, we initialize a k-arc with the cell with the smallest left border. Then, we progressively update and recode k-arcs. The third and fourth images present merge and split phases. We can notice that in one hand the recoding stage is not detailed in this figure, and in the other hand the edges ms represent a k-arc with one cell in this example

Fig. 5 .

 5 Fig.5. If we consider the original orientation of the k-arcs, the shape of the k-object presented in the next section (a left) is not well defined since the symetry is not preserved. So, we propose to start the reconstruction by the nodes s and m (a right). This structure respects the supercover model, and the symetrical shape of this object. We also show the result of our algorithm on the k-object presented in the previous section (b)

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig.6. An image of a chinese character (left), compressed by a quadtree-based approach (center). We show the final k-arcs recoding and the polygonalization (right)


Dinh-Van Nguyen
email: dinh-van.nguyen@inria.fr

Raoul De Charette

Fawzi Nashashibi

Trung-Kien Dao
email: trung-kien.dao@mica.edu.vn

Eric Castelli

WiFi Fingerprinting Localization for Intelligent Vehicles in Car Park

Keywords: component, Intelligent Vehicles, Indoor localization, WiFi Fingerprinting, Ensemble Neural Network, Particle Filter, Sensor Fusion

In this paper, a novel method of WiFi fingerprinting for localizing intelligent vehicles in GPS-denied area, such as car parks, is proposed. Although the method itself is a popular approach for indoor localization application, adapting it to the speed of vehicles requires different treatment. By deploying an ensemble neural network for fingerprinting classification, the method shows a reasonable localization precision at car park speed. Furthermore, a Gaussian Mixture Model (GMM) Particle Filter is applied to increase localization frequency as well as accuracy. Experiments show promising results with average localization error of 0.6m.

I. INTRODUCTION

As intelligent vehicles become closer to reality, a vast majority of researching resources are focusing on outdoor localization where GNSS is available. Works on indoor scenarios as well as the transition between indoor and outdoor are not yet fully addressed while a report shows that vehicles spend 95% of the time in car parks [START_REF] Barter | Cars are parked 95% of the time[END_REF]. Moreover, a study in France suggests a high average searching time for a free slot in car parks from 10 -40minutes in major cities which results in 70 million hours or 700 million euro cost a year [START_REF] Gantelet | The time looking for a parking space: Strategies, associated nuisances and stakes of parking management in France[END_REF]. This indicates a big gap to fill for intelligent vehicles localization.

Due to lack of space in urban cities, more and more multistory car parks are built to meet drivers' need. Inside these car parks, there are certain constraints about vehicles localization such as:

 No GPS signal.
 Average speed in a car park is around 10-12km/h [START_REF] Belloche | On-street Parking Search Time Modelling and Validation with Survey-based Data[END_REF].

 Areas with identical physical features.

Normally, a localization system for an intelligent vehicle is a fusion of multiple sensors such as: GPS antenna, LiDAR, camera, Inertial Measurement Unit (IMU), etc. Among them, the GPS antenna is critical for two primary roles: a global coordinate reference for fusion of multiple sensors (each with its own local coordinate) and a semantic indicator for a vehicle to perceive the surrounding environment. Thus, the lack of GPS signal for the indoor environment such car parks will most likely interrupt the vehicle localization system. Consequently, an alternative of GPS signal for indoor environment is required.

Speed is one of the keys constraints to the indoor localization of vehicles also. Although sensors like camera, LiDAR, IMU, etc. are doing well at much higher speed, those sensors do not provide a replacement for GPS signal or require a complex and expensive map of the environment to function. Other solutions for indoor movement tracking including beacons, WiFi signal, Ultra-Wideband, etc. often operate at low sampling frequency hence not suitable for the speed of a vehicle.

Another constraint mentioned above is the symmetrically identical structure of a car park. This characteristic of car parks increases the uncertainty of localization methods which use physical features such as laser/visual SLAM, image processing, etc. causing them to converge slower to true position. Also, symmetricity will add more confusion to the kidnapped robot problem [START_REF] Engelson | Error correction in mobile robot map learning[END_REF]. This further emphasizes the need for a global positioning system that can give absolute positioning feedback to the vehicle. This paper will address the issue of car park localization and the transition between the car park and outdoor environment. By introducing a bagging ensemble neural network and a hybrid offline phase learning to WiFi fingerprinting localization method, the method is capable of replacing GPS signal for the indoor environment localization. Also, a fusion with IMU using Gaussian Mixture Model Particle Filter will be presented to increase sampling frequency as well as accuracy of the whole system.

The paper is organized as follows: Section I gives a brief introduction of the problem and the idea for a solution. Section II summarizes state of the art of car park localization for vehicles. Section III explains WiFi Fingerprinting method as well as adaptation to vehicles localization. Section IV presents a Gaussian Mixture Model Particle Filter to fuse information from IMU and WiFi Fingerprinting method to increase sampling frequency and accuracy of localization. Section V gives details about experiment setup and results.

II. STATE OF THE ART

In the last few years, research community in Intelligent Vehicles has been developing several dedicated systems for car park localization and navigation. These systems can be categorized into two major classes: an absolute positioning approach with a pre-built map of the environment and a relative positioning approach with local feature detection.

In absolute positioning approach, two conditions are assumed: a pre-built map of the environment is available and little to no change in the environment static obstacles. This method often makes use of vision sensors such as LiDAR or camera. Work from [START_REF] Groh | Advanced Real-time Indoor Parking Localization based on Semi-Static Objects[END_REF] divides a map into known static objects like walls, parking slots, etc. and semi-static objects like parked cars. Instead of performing a full SLAM problem, a particle filter for localization on a pre-calculated probability grid map is deployed. A full SLAM approach can be found in [START_REF] Wahl | Localization inside a populated parking garage by using particle filters with a map of the static environment[END_REF] with Rao-Blackwellized particle filter SLAM together with a static map of the environment. The idea is to reduce the uncertainty of the SLAM with a known static map thus improving the accuracy of localization. Instead of mounting LiDAR sensor on vehicles, the approach in [START_REF] Ibisch | Towards autonomous driving in a parking garage: Vehicle localization and tracking using environmentembedded LIDAR sensors[END_REF] proposes to embed LiDAR sensors in the environment. Using RANSAC algorithm, the authors have successfully localized vehicles within 8cm of mean errors. A camera approach in this category is found in [START_REF] Bojja | Indoor Localization Methods Using Dead Reckoning and 3D Map Matching[END_REF]. The method builds a detailed 3D map of the environment together with dead reckoning algorithm and particle filter to localize vehicles inside a car park. This method also requires a GPS position initialization for the 3D matching process. In general, methods fall into this category often perform well for a carefully mapped environment. However, these methods are expensive in computation as well as in map building process and sometimes also costly sensors. They are dedicated to one environment without any or little consideration of transition phase for the system from indoor to outdoor or vice-versa.

For relative positioning approach, vehicles are located relatively to surrounding local objects such as other vehicles, lane marking, parking lot positions etc.. An example of this method can be found in [START_REF] Levinson | Map-Based Precision Vehicle Localization in Urban Environments[END_REF] with lane extraction and relative positioning. A vehicle to vehicle communication method to improve localization accuracy is proposed in [START_REF] Amini | Improving GPS-Based Vehicle Positioning for Intelligent Transportation Systems[END_REF]. A method for detecting parking slot position in [START_REF] Houben | On-Vehicle Video-Based Parking Lot Recognition with Fisheye Optics[END_REF] using multiples cameras on vehicles yields 0.2 meters of error in localization. One notable characteristic of this approach is the uncertainty of relative positioning increase significantly as the uncertainty of each agent adds up. Also, this method suffers heavily from accumulated error as there is little chance to do correction once error is made.

III. ENSEMBLE BAGGING WIFI FINGERPRINTING

WiFi Fingerprinting localization is a popular method for indoor localization typically for pedestrian walking. The method makes use of existing WiFi infrastructure (WiFi Access Points -APs, WiFi receiver) to determine location based on an offline mapping phase. The main argument of this method is that the combination of WiFi signal strengths from multiple static APs in the environment for one position is unique. By learning these unique features for several key positions in the environment, one can estimate its location just by scanning WiFi signal strength. Generally, WiFi Fingerprinting has two phases: offline phase and online phase.

Offline phase of WiFi Fingerprinting method is where signal strengths from multiple locations are gathered and stored with where 𝑥 𝑖,𝑗 is WiFi signal strength from jth WiFi APs that is recorded from the absolute coordinate 𝜌 𝑙 of 𝐹𝑃 𝑙 in ith scan. A collection of multiple scans from each FP in this offline phase will act as a reference for estimating position in online phase. One advantage of this method is the coordinate of fingerprint location can be measured in global frame (i.e. GPS frame) so that the transition from the indoor to the outdoor environment can be done seamlessly. At the same time, this method does not suffer from symmetricity of car park and highly feasible and scalable due to availability of APs in the environment.

In the online phase where localization estimation is carried out, the vehicle will move inside the environment while scanning for WiFi signal strength from surrounding APs. A likelihood function based on data from offline phase is defined as:

ℎ(𝑋) = 𝛾 𝜌 𝑙 (2)
Here, 𝑋 is input signal scan and 𝛾 𝜌 𝑙 is the likelihood score of 𝑋 to be scanned at 𝐹𝑃 𝑙 in the environment with regard to offline collected mapping coordinate 𝜌 𝑙 . In most of the cases, the position with the highest likelihood score will be the output of localization process.

WiFi Fingerprinting appears to perform well for a human user in multiple indoor scenarios with modification in either hardware or software. A 20-40cm of error are reported in [START_REF] Kotaru | SpotFi : Decimeter Level Localization Using WiFi[END_REF]- [START_REF] Kumar | Accurate indoor localization with zero start-up cost[END_REF]. However, there are several issues with this method when it comes to vehicles movement:  Low sampling frequency: the time to complete a scan of WiFi signals in a particular environment depends on multiple factors but generally falls into 1 second range (1Hz). This sampling frequency is neglected with human walking speed around 1.4 m/salthough in practice, most systems require user to walk at a lower speed for a better precision. In contrast, a vehicles moving in a car park is expected to move at 3 -3.3 m/s [START_REF] Belloche | On-street Parking Search Time Modelling and Validation with Survey-based Data[END_REF]. This results in a single measurement every 3 meters. In addition, localization system for a vehicle is much more demanding in term of precision and sampling frequency. Consequently, a classic WiFi fingerprinting system alone is inadequate for vehicle localization.

 Due to the high speed movement of vehicles, the signal scan of WiFi is expected to have higher variance and noise. In fact, with one measurement every 3 meters, the signal scan vector does not reflect any single position characteristic but rather a range.

To solve these issues, changes in both offline phase and online phases are proposed. A hybrid learning database in offline phase and an ensemble neural network for online phase is implemented as follows.

A. Hybrid Offline Phase

In classical approach of WiFi Fingerprinting localization method, offline collection of mapping between position (Fingerprint -FP) and Received Signal Strength Indicator (RSSI) from Access Points (APs) in the environment is done statically. This means at each FP, a number of scans for RSSI are performed and stored to model the signal strength feature of that particular position. However, as suggested above, this method likely produces biased observation of signal strength as vehicles moving at much faster speed thus results in a range of signal scans instead of static position signal scan. We then define a scan range as the distance traveled by a vehicle from initiation to complete of a WiFi signal scan.

Consider a typical scan of RSSI in online phases as in Fig. 2. For 1 scan of signal strength at 1Hz sampling frequency, the vehicle is likely to pass by a Fingerprint defined in offline phase. In the best case, the vehicle is likely to travel 3.3m (at 3.3m/s), passing by the fingerprint at exact centre of the path as shown in the figure. This is a significant range that cannot be represented with only static scans at fingerprint position in offline phase.

One way to address this issue is to include dynamic signal scans in offline phase to emulate the character of a real-time online phase observation. Instead of only collecting scan statically at fingerprint position in offline phase, the vehicle will also move at different speeds passing by fingerprint positions. Multiple scans within a representative area around a fingerprint will also be recorded.

{𝑋 𝑖 , 𝜌 𝑙 } = {𝑥 𝑖,1 , 𝑥 𝑖,2 , 𝑥 𝑖,3 , … , 𝑥 𝑖,𝑛 , 𝜌 𝑙 ′ } (3)
Where

𝜌 𝑙 ′ ∈ (𝜌 -∆𝜌, 𝜌 + ∆𝜌) (4)
The representative area for a fingerprint position is defined as in [START_REF] Engelson | Error correction in mobile robot map learning[END_REF]. Note that, as in (3), although the real position 𝜌 ′ that is mapped with the scan is different from fingerprint position 𝜌, the final mapping will still be the static fingerprint position to preserve fingerprinting characteristic. Thus all scans collected within the representative area of a fingerprint will be mapped to the fingerprint position only.

∆𝜌 = 𝑣̅ 𝑣𝑒ℎ𝑖𝑐𝑙𝑒 2 ⁄ (5)
The fingerprint representative area parameter ∆𝜌 will be determined as in [START_REF] Groh | Advanced Real-time Indoor Parking Localization based on Semi-Static Objects[END_REF] with average speed 𝑣̅ 𝑣𝑒ℎ𝑖𝑐𝑙𝑒 since in the best case scenario, the vehicle will pass by the fingerprint location at exactly the middle of scanning. One reason for that is if worst case scenario where the vehicle passed by the fingerprint location at either end of the scanning cycle were considered, then the representative area would be much wider and therefore would introduce more noise to the collected data. The worst case scenario also did not reflect the desired feature for targeted fingerprint location. In case of the car park with an average speed of a vehicle is 3.3m/s, ∆𝜌 is approximately 1.6𝑚. This closely brings back the problem to human walking speed if we consider the representative area as a fingerprint.

Finally, a normalization of collected RSSI in each scan is performed to reduce the numerical impact of signal strength feature. A detected AP RSSI would be normalized in the range [0, 1) with 0 as the weakest possible signal strength and 1 as the strongest possible signal strength.

𝑥 𝑖 = { -1, 𝐴𝑃 𝑖 𝑢𝑛𝑑𝑒𝑡𝑒𝑐𝑡𝑒𝑑 1 - (-1)×𝑅𝑆𝑆𝐼 100 , 𝐴𝑃 𝑖 𝑑𝑒𝑡𝑒𝑐𝑡𝑒𝑑 (6)

B. Ensemble Neural Network with Bagging

In online phase, a likelihood function is designed to estimate which fingerprinting is closest to input scanning signals. There are several algorithms to estimate likelihood score of input scanning signal. In this paper, an Ensemble Neural Network with Bagging is proposed to be the classifier. This method is known to perform well with high variance and noisy data. In this case, the WiFi signal strength collected during online phase is noisy with high variance as suggested in the previous section.

Ensemble neural network with bagging (Bootstrap Aggregating) is a method of combining multiple learning models to derive better prediction results [START_REF] Breiman | Bagging Predictors[END_REF], [START_REF] Dietterich | Ensem ble Methods in Mac hine Learning[END_REF]. Consider a mapping {𝑋, 𝑌} where 𝑋 is a vector of features and 𝐶 is corresponding class. In this paper, 𝑋 is a vector of signal strengths scanned at a particular position in the environment and 𝑌 is the label of FP position defined in offline phase. A function estimator which results from a set of training samples and a classification model is formed.

Human speed scan range Scan Initiation Scan complete

Vehicle speed scan range Fingerprint 𝑔(•) = ℎ((𝑋 1 , 𝑌 1 ), (𝑋 2 , 𝑌 2 ), … , (𝑋 𝑛 , 𝑌 𝑚 ) ))

Bagging algorithm is executed as follows:

 Construct bootstrap sample by randomly sampling with replacement n times from original training data. Let 𝑋 ^i and 𝑌 ^i is a bootstrap sample data from collected database {𝑋, 𝑌}:

(𝑋 ^1, 𝑌 ^1), (𝑋 ^2, 𝑌 ^2), … , (𝑋 ^n, 𝑌 ^m)

 Compute bootstrap estimator 𝑔 ^(•) by applying same classification model to newly formed bootstrap sample:

𝑔 ^(•) = ℎ((𝑋 ^1, 𝑌 ^1), (𝑋 ^2, 𝑌 ^2), … , (𝑋 ^n, 𝑌 ^m))) (9)
 Repeat two steps above for K times with K is large the bagging estimator is:

𝑔 ^𝑏𝑎𝑔𝑔 (•) = 1 𝐾 (∑ 𝑔 ^𝑖(•) 𝐾 𝑖=1 ) (10)
Theoretically, as K goes to infinity, the estimator converges to the expected value. In practice, a finite large K is expected to improve Monte Carlo accuracy.

Since vehicles are running at different speed, scanned RSSI is expected to be noisy and high variance. Ensemble Neural Network appears to be functioning well with high variance and noisy data. By picking randomly with replacement training data from whole dataset, the combination of final learning models is able to overcome bias and eliminate noise.

IV. GUASSIAN MIXTURE MODEL PARTICLE FILTER

As mentioned in previous section, WiFi fingerprinting localization suffers from low sampling frequency. To enhance the performance, a fusion with other sensors information is recommended. In this paper, a Gaussian mixture model particle filter is deployed to fuse WiFi fingerprinting localization with other sensor such as IMU information. The fusion strategy is illustrated in Fig. 3.

Given a signal scan, a classification result from the ensemble neural network is defined as a list of fingerprints and their corresponding confidence {𝐹𝑃 𝑖 , 𝐶 𝑖 }. In classic fingerprint solution, it is likely that the highest confidence fingerprint will be the final prediction of the method. However, due to environment noise and multipath effect of signal transmission, the highest confidence fingerprint does not always represent the best estimation. One way to improve the prediction accuracy is to consider few fingerprints with highest confidence.

A particle filter with GMM model is established as following:

 Initialization: The filter is initialized after the first frame of WiFi Fingerprinting. Particles are generated around the highest confidence fingerprint that results from the ensemble neural network with a standard deviation 𝜎. (𝜎 is estimated from the statistic of WiFi fingerprinting localization experiments).

 Updating transition: A constant speed model that receives velocity and heading angle of the vehicle as inputs. This model moves all particles with an assumption that at 10Hz frequency of IMU, movement of a vehicle can be considered as constant speed.

 Predicting transition: Top 3 highest confidence FPs returned from network will be considered as possible observation. Score of each particle will be updated as in [START_REF] Houben | On-Vehicle Video-Based Parking Lot Recognition with Fisheye Optics[END_REF]. Here, 𝑥 𝑡 is a particle at time t, 𝑐 𝑖 is the confidence of fingerprint ith returned by ensemble neural network, 𝜎 𝑖 is standard deviation of observation at fingerprint ith and 𝜇 𝑖 is coordinate of Fingerprint ith.

𝐺𝑀𝑀(𝑥 𝑡 ) = ∑ 𝑐 𝑖 𝑃(𝑥 𝑡 | 𝜎 𝑖 , 𝜇 𝑖 ) 𝑖 (11)
A Gaussian distribution is assumed with multiple observations are considered. Each observation has a posteriori confidence score. The higher the confidence, the smaller variance from mean is expected. Thus, given an estimated standard deviation 𝜎 for WiFi fingerprinting localization observation, the standard deviation for each fingerprinting observation in this case will be calculated as:

𝜎 𝑖 = 1 -𝑐 𝑖 𝜎 (12)
Since WiFi scan frequency goes at 1Hz and IMU data is filtered at 10Hz, there are cases where there is no observation is fed into particle filter except IMU information. At these cases, all particles will be updated with constant speed model with equal score. Final result of particle filter estimation is weighted sum of all particles with its score.

𝑥 = 1 𝑛 ∑ 𝑥 𝑡 ∑ 𝑐 𝑖 1 √2𝜋(1-𝑐 𝑖 𝜎) 2 𝑒 - (𝑥 𝑡 -𝜇 𝑖 ) 2 2(1-𝑐 𝑖 𝜎) 2 ⁄ 𝑖 𝑡 (13)
Note that by talking weighted sum of all particles instead of just the highest score particle as the final prediction, this allows room for recovery from biased observation. After receiving predicted result, a random proportion of particles is sampled with replacement using multinomial resampling method.

V. EXPERIMENTS AND RESULTS

During nearly one year period, experiments have been carried out in Rocquencourt campus -INRIA Paris. The setup of experiments is explained as follows.

Fig. 3. Fusion of IMU and WiFi using GMM Particle Filter

There are two testing platforms: a Cybercars (blue) and an electric car Citroën C1 (red) (Fig. 4). The Cybercar is used in both online and offline phase of WiFi fingerprinting localization. The C1 is used in testing phase only. Sensors on both vehicles are:

 RTK GPS antenna

VI. CONCLUSION

In this paper, a novel method of localizing intelligent vehicles in a GPS-denied environment such as car park is proposed. A number of constraints are identified such as lack of GPS, vehicle speed and symmetricity of the environment. The proposed method has two main components: (1) an ensemble neural network for WiFi fingerprinting localization and (2) a Gaussian Mixture Model Particle Filter to fuse WiFi fingerprinting with IMU information.

In the first component, by applying a hybrid offline phase for data collection and introducing an ensemble neural network, localization results with 2.25m mean error were obtained. By mapping signal scan at fingerprinting location to GPS frame coordinate, it also allows seamless transition between environments.

The second component provides improvements for localization system in both sampling frequency and accuracy. By deploying a Gaussian mixture model into a bootstrap Particle Filter, it allows the system to overcome WiFi fingerprinting limitations.

Experiments in the span of nearly a year from March 2017 to January 2018 shows promising results of localization for intelligent vehicles. A 99.3 % of error estimation is under one sigma (with σ=1.5). Mean error of localization is 0.62m with low standard variation of 0.175. These results prove the correctness and stability of the proposed system in a long run.

Still, there are several key factors which need to be addressed in this work such as: a real multistory carpark situation for experiments or changes in the environment condition (i.e. weather, physical changes, etc.). The experiments in this work did not address z-axis of the localization problem (for multistory carpark case) and there was no experiment conducted under extreme weather condition due to the nature of RTK GPS ground truth. Also, no changes in physical condition such as new structure, changes in APs location was recorded. These are factors which could potentially hurt the performance of the proposed system. In the future, these factors will be tackled.

Fig. 1 .

 1 Fig. 1. Illustration of offline phases in WiFi Fingerprinting corresponding to coordinate of each location. As shown in Fig. 1, each blue dot will be called a fingerprint (FP). For each fingerprint, a mapping of scanned WiFi signal strength and its coordinate is stored. {𝑋 𝑖 , 𝜌 𝑙 } = {𝑥 𝑖,1 , 𝑥 𝑖,2 , 𝑥 𝑖,3 , … , 𝑥 𝑖,𝑛 , 𝜌 𝑙 } (1)

Fig. 2 .

 2 Fig. 2. Scan range of vehicles in online phase.

Fig. 7 .

 7 Fig. 7. (a) Single run with initial position out of fingerprint range. (b) cumulative distribution of error.

Fig. 8 .

 8 Fig. 8. (a) Distribution of error estimation in all test runs; (b) Cumulative distribution of error for all test runs.

All experiments are performed in the outdoor environment, near parking lot and surrounded by buildings (instead of an indoor car park). The reason for adopting the outdoor scenario is to take advantage of centimetre precision RTK GPS as the ground truth for experiments' results evaluation later. The area of experiments is shown in Fig. 5.

In the offline phase of WiFi Fingerprinting method, the Cybercar collected static RSSI scan at fingerprints position which was chosen to be at roughly every 5m apart on the travelling route. Dynamic scans were collected in next 3 days, after static scans, and at some random times of the day.

In the online phase, both Cybercar and C1 are running tests on two localization modes: only WiFi Fingerprinting and fusion of WiFi and IMU using GMM Particle Filter. Both car are moving within 2.5 -3.5m/s of average speed in all experiments to comply with the average speed of vehicles moving inside a car park.

For WiFi only localization method using fingerprinting, results obtained were shown in [START_REF] Dinh-Van And | Indoor Intelligent Vehicle localization using WiFi Received Signal Strength Indicator[END_REF]. Within 1012.9m of the travel path, the maximum error recorded is 6.84m, and the average error is 2.25m. However, the sampling frequency of this method is restricted to the sampling frequency of WiFi scan which is 1Hz.

By fusing this result with IMU information using GMM Particle Filter, the localization performance is drastically improved. Fig. 5 shows a smoother and much more accurate localization outcome of a single test run is obtained in comparision to the result of WiFi fingerprinting method alone. With GMM fusion, the localization system is capable of running at 10Hz frequency instead of just 1Hz. Not only that, the error of localization is significantly reduced. Having statistic from WiFi fingerprinting localization results in [START_REF] Dinh-Van And | Indoor Intelligent Vehicle localization using WiFi Received Signal Strength Indicator[END_REF], a standard deviation 𝜎 = 1.5 is used for the fusion.

From March 2017 to January 2018, 64 test runs were recorded with Cybercar and C1. Fig. 6a shows a single run statistic with maximum error is 1.223m and mean error is 0.473m. The mean error estimation and variance of each experiment error are evaluated. A portion of those data is presented in Fig. 6b. For all experiments, the maximum mean error is 0.62m and the maximum variance is 0.175. This statistic shows the stability of the method over the course of one-year experiments. Note that, during this time, there was no significant change in the WiFi infrastructure and the environment structure. We also targeted only infrastructure mode WIFI networks in the area to reduce noise in the signal collection process.

In addition, with a random initial position, there are cases where the first position of the vehicle is not within the range of any fingerprint. These cases are examined and results are shown in Fig. 7a. Here, due to the initial position, the maximum error is large at 3.28m. However, after 23 frames, the particle filter quickly converges to reduce the error by 47% to reach 1.722m of maximum error. This happens due to the observation from WiFi fingerprinting localization helps converging the particle filter. Once a correct observation is made, the bias from starting position is quickly reduced. The results also prove that WiFi GMM fusion plays a critical role in correcting bias and accumulated error from the IMU motion model. The mean error for the entire path is recorded at 0.77m. The cumulative distribution of error for this case is shown in Fig. 7b. At 𝜎 = 1.5, 95.3% of error estimations are below one sigma. Furthermore, a distribution of all error estimations for every conducted experiments is shown in Fig. 8a. The distribution confirms the observation made from Fig. 6b with most of the error estimations are around 0.5m. The cumulative distribution of these estimation demonstrates similar result in Fig. 8b. With 𝜎 is estimated around 1.5, 99.3% of the error estimation from all experiments is below one sigma. This suggests a reliable prediction in localization process of the system in general.