Y Grandati

A short proof of the Gaillard-Matveev theorem based on shape invariance arguments

Keywords:

We propose a simple alternative proof of the Wronskian representation formula obtained by Gaillard and Matveev for the Darboux-Pöschl-Teller (DPT) potentials. It rests on the use of singular Darboux-Bäcklund transformations applied to the free particle system combined to the shape invariance properties of the DPT.

I. INTRODUCTION

In 2002, Gaillard and Matveev established an interesting connection between the Darboux-Crum dressing formula and the Darboux-Pöschl-Teller potentials. More precisely, they proved the following theorem:

Gaillard-Matveev's theorem: The trigonometric Darboux-Pöschl-Teller potential (TDPT) with integer parameters m ≥ n defined on x ∈ ]0, π/2[ by

V (x; m, n) = m(m + 1) sin 2 x + n(n + 1) cos 2 x , (1)
admits the following Wronskian representation

V (x; m, n) = -2 (ln W (u 1 , ..., u m | x)) ′′ , (2)
where W (u 1 , ..., u m | x) denotes the Wronskian of the functions u 1 (x) , ..., u m (x) which are given by u k (x) = sin(lx), if 1 ≤ l ≤ m -n u m-n+l (x) = sin((m -n + 2l) x), if 1 ≤ l ≤ n.

(

) 3
Gaillard and Matveev underlined that this result traduces the existence of a chain of Darboux transformations connecting the constant potential to the TDPT V (x; m, n). If such a construction was already envisaged in Darboux's seminal works [START_REF] Darboux | Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitésimal[END_REF] , the theorem above precises what is the correct choice of "seeds" eigenfunctions of the free Hamiltonian to reach this last potential.

The proof given in 2 and 3 rests on direct evaluations of the Wronskian W (u 1 , ..., u m | x) and is quite sophisticated and technically involved. In this letter, we show that this result can be obtained in a shorter way directly by following Darboux's original idea. In a modern formulation the proof makes use of shape invariance and SUSY QM arguments associated to singular Darboux-Bäcklund Transformations (DBT) built from excited eigenstates. If the use of singular DBT has already been envisaged in several papers [START_REF] Casahorran | A family of supersymmetric quantum mechanics models with singular superpotentials[END_REF][START_REF] Panigrahi | Singular superpotentials in supersymmetric quantum mechanics[END_REF][START_REF] Robnik | Supersymmetric quantum mechanics based on higher excited states[END_REF][START_REF] Marquez | Factorization method and singular hamiltonians[END_REF] , they are generally considered of less interest, although recently proven to be crucial in the study of some two-dimensional superintegrable systems [START_REF] Marquette | Singular isotonic oscillator, supersymmetry and superintegrability[END_REF] .

We also examine the case of the Bessel potentials, which is the first example envisaged by Darboux 1 . This is a confluent case which necessitate to employ Matveev's generalized Wronskian leading then to a Wronskian version of the Rayleigh formula for Bessel functions.

II. DARBOUX-B ÄCKLUND TRANSFORMATIONS

We begin to briefly recall ther essential features of DBT and shape invariance. If ψ λ (x; a) is an eigenstate of the hamiltonian

H(a) = -d 2 /dx 2 + V (x; a), x ∈ I ⊂ R, (a ∈ R m being a multiparameter) associated to the eigenvalue E λ (a) (E 0 (a) = 0) ψ ′′ λ (x; a) + (E λ (a) -V (x; a)) ψ λ (x; a) = 0, (4)
then the Riccati-Schrödinger (RS) function w λ (x; a) = -ψ ′ λ (x; a)/ψ λ (x; a) satisfies the corresponding Riccati-Schrödinger (RS) equation 9

-w ′ λ (x; a) + w 2 λ (x; a) = V (x; a) -E λ (a). (5)
The class of RS equations is preserved by a specific subset of the group G of smooth SL(2, R)-valued curves M ap(R, SL(2, R)) . These transformations, called Darboux-Bäcklund Transformations (DBT), are built from any solution w ν (x; a) of the initial RS equation Eq(5) as [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF][START_REF] Cariñena | Integrability of Riccati equation from a group theoretical viewpoint[END_REF][START_REF] Cariñena | Group theoretical approach to the intertwined hamiltonians[END_REF] w λ (x; a)

A(wν ) → w (ν) λ (x; a) = -w ν (x; a) + E λ (a) -E ν (a) w ν (x; a) -w λ (x; a) , (6)
where

E λ (a) = E ν (a). w (ν)
λ is then a solution of the RS equation:

-w

(ν)′ λ (x; a) + w (ν) λ (x; a) 2 = V (ν) (x; a) -E λ (a), (7)
with the same energy E λ (a) as in Eq( 5) but with a modified potential

V (ν) (x; a) = V (x; a) + 2w ′ ν (x; a) (8)
called an extension of V . If V (ν) is a rational function of an a priori given variable (x, tan x, sin x, cosh x...) we call it a rational extension of V . This can be schematically summarized as

   w λ A(wν ) w (ν) λ V A(wν ) V (ν) . (9)
The corresponding eigenstate of H (ν) (a) = -d 2 /dx 2 + V (ν) (x; a) can be written

ψ (ν) λ (x; a) = exp -dxw (ν) λ (x; a) ∼ 1 E λ (a) -E ν (a) A (w ν ) ψ λ (x; a), (10)
where A (w ν ) is a first order operator given by A (w ν ) = d/dx + w ν (x; a). Eq (10) can still be written as

ψ (ν) λ (x; a) ∼ W (ψ ν , ψ λ | x) ψ ν (x; a) , (11)
where W (y 1 , ..., y n | x) is the Wronskian of the functions y 1 , ..., y n . From V , the DBT generates a new potential V (ν) (quasi) isospectral to the original one and its eigenfunctions are directly obtained from those of V via Eq (10). Nevertheless, in general, w ν (x; a) and then the transformed potential V (ν) (x; a) is singular at the nodes of ψ ν (x; a). For instance, if H(a) admits a bound state spectrum (E n , ψ n ) n∈N (x; a), V (n) is regular only when n = 0, that is when ψ n=0 is the ground state of H, and we recover the usual SUSY partnership in quantum mechanics [START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF] . Note that A(w 0 ) is a "state-deleting" DBT, the spectrum of the superpartner hamiltonian H = H (0) having for fundamental level and ground state E 1 and ψ (0) 1 respectively. V is said to be a translationally shape invariant potential (TSIP) if its superpartner has the form

V (x; a) = V (0) (x; a) = V (x; a + α) + R(a), (12)
where α ∈ R m . The question of successive iterations of DBT is very natural and is at the center of the construction of the hierarchy of hamiltonians in the usual SUSY QM scheme [START_REF] Sukumar | Supersymmetry, factorization of the Schrödinger equation and an hamiltonian hierarchy[END_REF] . Staying at the formal level, it can be simply described by the following generalization of Eq(9) (N j denotes the j-uple (ν 1 , , ..., ν j ), with

N 1 = ν 1 and N m = (N m-1 , ν m ))      w λ A(wν 1 ) w (ν1) λ A(w (N 1 ) ν 2 ) w (N2) λ ... A(w (N m-1 ) νm ) w (Nm) λ V A(wν 1 ) V (ν1) A(w (N 1 ) ν 2 ) V (N2) ... A(w (N m-1 ) νm ) V (Nm) , (13)
where w

(Nm) λ is a RS function associated to the eigenvalue E λ of the potential

V (Nm) (x; a) = V (x; a) + 2 m j=1 w (Nj-1) νj (x; a) ′ . ( 14
)
The corresponding eigenfunction is given by (cf Eq (10) and Eq (11))

ψ (Nm) λ (x; a) = A w (Nm-1) νm ψ (Nm-1) λ (x; a) = W ψ (Nm-1) νm , ψ (Nm-1) λ | x ψ (Nm-1) νm (x; a) . (15)
or

ψ (Nm) λ (x; a) = A w (Nm-1) νm ... A (w ν1 ) ψ λ (x; a). (16)
By induction, we can prove the Crum formulas 13,14

ψ (Nm) λ (x; a) = W (ψ ν1 , ..., ψ νm , ψ λ | x) /W (ψ ν1 , ..., ψ νm | x) V (Nm) (x; a) = V (x; a) -2 (log W (ψ ν1 , ..., ψ νm | x)) ′′ . (17)
A succession of m SUSY QM partnerships gives rise to a hierarchy of regular potentials V (Nj ) , j = 1, ..., m, associated to the m-uple N m = (0, , ..., m -1) , that is, to a set of seeds functions (ψ 0 , ..., ψ m-1 ) constituted by the m first bound states of the initial hamiltonian. In the case where the initial potential is a TSIP, we obtain simply

V (Nm) (x; a) = V (x; a) -2 (log W (ψ 0 , ..., ψ m-1 | x)) ′′ = V (x; a + mα) + m-1 j=0 R(a j ), (18)
with

ψ (Nm) k+m (x; a) = W (ψ 0 , ..., ψ m-1 , ψ k+m | x) W (ψ 0 , ..., ψ m-1 | x) , k ≥ 0, (19)
which is the k-th excited state of V (Nm) with the associated energy E k+m .

III. FROM THE CONSTANT POTENTIAL TO THE TDPT POTENTIAL

Following Darboux 1 , our starting point is the constant potential on the real line V (x) = 0, x ∈ R. The physical energy spectrum of the corresponding hamiltonian H is constituted by an open band

E k = k 2 ∈ ]0, +∞[ , k > 0,
where the two-dimensional eigenspace associated to E k is spanned by the (unnormalized) scattering eigenstates

φ odd k (x) = sin(kx)/k, φ even k (x) = cos(kx)/k. ( 20
)
We then build a chain of extensions of V where the first DBT we choose is based on a scattering state satisfying the Dirichlet boundary condition at the origin, namely

ψ 0 (x) = φ odd 1 (x) = sin (x) , w 0 (x) = -cot (x) , (21)
where the choice of the value k = 1 is for pure convenience. The first extended potential is then

V (0) (x) = V (x) + 2w ′ 0 (x) = 2 sin 2 x ( 22
)
and the corresponding Schrödinger equation

-d 2 /dx 2 + V (0) (x) ψ(x) = Eψ(x). (23)
At this step, Darboux [START_REF] Darboux | Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitésimal[END_REF] writes "Pour des valeurs particulières de E, cette équation admet les solutions sin 2 x, sin 2 x cos x. En les employant et en poursuivant l'application de la méthode, on parviendra successivement à des équations de la forme

-d 2 /dx 2 + n(n + 1) cos 2 x + m(m + 1) sin 2 x ψ(x) = Eψ(x)."
As emphasized by Gaillard and Matveev, the argument needs to be completed and specified in particular to determine in a systematic way what is the correct choice of the successive seeds functions to reach this equation.

First note that Darboux's procedure does not consider the spectral properties of the built hamiltonians which can then become more singular at each step. This is the case of V (0) which is singular at every node ξ l = lπ of ψ 0 (x). The corresponding singularities are of the centrifugal barrier λ/ (x -ξ l ) 2 with λ > 3/4, that is, strong singularities for which the transmission probability is zero [START_REF] Panigrahi | Singular superpotentials in supersymmetric quantum mechanics[END_REF][START_REF] Frank | Singular potentials[END_REF][START_REF] Lathouwers | The hamiltonian H = -1 2 d 2 dx 2 + x 2 2 + λ x 2 reobserved[END_REF] . This confines the particle in a single interval ]ξ l , ξ l+1 [. If we take for instance l = 0, the existence of the singularities in 0 and π can then be considered as equivalent to the addition of Dirichlet boundary conditions at the origin and at x = π.

ψ 0 + = ψ π -= 0. ( 24
)
These lasts act as a filter which select among all the formal eigenfunctions ψ (0) k only those which satisfy Eq (24). The question is now to determine what are the appropriate values of k which allows to verify this constraint.

From a more general point of view, consider the problem to determine what is the condition to which the solution ψ λ (x) of Eq(4) is submitted in order that its image ψ (ν) λ via the DBT A (w ν ) satifies

ψ λ (x 0 ) = 0, (25)
where x 0 ∈ I is a zero of ψ ν . In the vinicity of x 0 we can write

ψ µ (x) = a µ,0 + a µ,1 (x -x 0 ) + a µ,2 (x -x 0 ) 2 + O (x -x 0 ) 3 . ( 26
)
x 0 being a simple zero of ψ ν we have a ν,0 = 0 and a ν,1 = 0, which gives

w ν (x) ≃ x→x0 -1 x -x 0 - a ν,2 a ν,1 + O((x -x 0 )). ( 27
)
Inserting these results in Eq (11), we obtain

ψ (ν) λ (x) ∼ x→x0 -a λ,0 x -x 0 - a ν,2 a ν,1 a λ,0 + O((x -x 0 )) (28)
and the condition ( 25) is satisfied iff a λ,0 = 0, that is, iff x 0 is also a zero of ψ λ .

In brief, a DBT A (w ν ) associated to a seed function ψ ν having nodes on the definition interval generates an extended hamiltonian which is singular on this interval. The restriction of this extension between two singularities has a spectrum which consists in a filtered version of the initial hamiltonian spectrum: we only keep the levels E λ for which the initial eigenstates ψ λ have nodes at the singular points. This point of view is equivalent to the one of Marquez et al [START_REF] Marquez | Factorization method and singular hamiltonians[END_REF] .

Applying this analysis to to the system considered above, we obtain that the potential V (0) restricted to the positive half line and subjected to the boundary conditions Eq (24) has a discrete non degenerate spectrum given by E l , l ∈ N * , with the corresponding eigenstates

ψ (0) l (x) = A (w 0 ) ψ l (x) ∼ (l + 1) cos ((l + 1) x) -cot x sin ((l + 1) x) , (29)
where we have noted ψ l (x) = φ odd l+1 (x) and E l = E l+1 . In particular, the fundamental and first excited eigenstates are ψ (0) 1 (x) = sin 2 (x) and ψ (0) 2 (x) = sin 2 x cos x. V (0) can also be considered as the image via A (w 0 ) of the infinite square well

V SW (x) = 0, if x ∈ ]0, π[ +∞, if x / ∈ ]0, π[ . (30)
V (0) being a particular TDPT potential (1) with n = 0, m = 1, we can consider that V SW correspond to the limit case

V SW (x) = V (x; 0, 0) = lim µ,ν→0 V (x; µ, ν).
As it is well known [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF][START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF][START_REF] Gendenshtein | Derivation of exact spectra of the Schrodinger equation by means of supersymmetry[END_REF] , the TDPT potentials V (x; µ, ν) are second category translationally shape invariant potentials (TSIP), the SUSY QM partner of V (x; µ, ν) being given by

V (x; µ, ν) = V (x; µ + 1, ν + 1), (31)
for µ, ν > 0 and

V (x; µ, 0) = V (x; µ + 1, 0), µ > 0 V (x; 0, ν) = V (x; 0, ν + 1), ν > 0. (32)
Starting from V (0) (x) = V (x; 1, 0), we can build a hierarchy of p -1 extensions of via successive SUSY QM partnerships, that is, via a chain of p -1 DBT based on the successive ground states. From Eq(32) and Eq (18), we deduce that the final potential is V (x; p, 0) and can be seen as the end of a chain of p extensions which starts from the zero potential V (x) and which is associated to the p-uple N p = (0, ..., p -1), ie to the set of seeds eigenfunctions (ψ 0 (x) , ..., ψ p-1 (x)) = (sin (x) , ..., sin (px)). Using a straightforward recurrence, we deduce that the ground state of V (Np) (x) = V (x; p, 0) = p(p + 1)/ sin 2 x at energy E p is

ψ (Np) p (x) = sin p+1 x. (33)
Since is even with respect to π/2, its eigenfunctions ψ (Np) n are repectively symmetric (if n-p is even) or antisymmetric (if n -p is odd) with respect to this point. An immediate consequence is that its first excited eigenstate (at energy

E p+1 ) ψ (Np)
p+1 has only one node on ]0, π[ which is always located at x = π/2 independently of p. From Eq(33) we deduce

ψ (Np) p+1 (x) ∼ A + w (Np+1) p+1 ψ (Np+1) p+1 (x) ∼ cos x sin p+1 x. ( 34
)
Suppose that we continue the chain of extension and build the following DBT from this first excited eigenstate. The obtained potential is then singular and presents a pole of second order at x = π/2. ψ

p+1 (x) = A w (Np-1) p ... A (w 0 ) ψ p+1 (x), ( 35
)
the new chain is associated to the (n + 1)-uple (N p , p + 1) and the final extension is given by

V (Np,p+1) (x) = 2 cos 2 x + (p + 1) (p + 2) sin 2 x . ( 36
)
The presence of the additional "strong" singularity in π/2 imposes a Dirichlet boundary condition at this point. In the spectrum of the restriction of V (Np,p+1) to ]0, π/2[, half of the levels E (Np) n of the preceding potential V (Np) are eliminated. These are all those associated to "even" (with respect to π/2) eigenstates, that is to the even values of the quantum number n -p.

The spectrum of V (Np,p+1) (x) contains then only the levels E p+1+2(j+1) , j ≥ 0, and the associated eigenstates

ψ (Np,p+1)
p+1+2(j+1) are the images of the initial eigenstates ψ p+1+2(j+1) via the chain of BDT corresponding to the (p + 1)uple (N p , p + 1) = (0, ..., p -1, p + 1). Starting from V (Np,p+1) we can continue the chain by using standard SUSY QM partnerships, all the extended potentials thus generated being perfectly regular on ]0, π/2[. But V (Np,p+1) is also a TDPT

V (Np,p+1) (x) = V (x; p + 1, 1). ( 37
)
Consequently the shape invariance property Eq(31) implies that after q steps, we obtain as final extension

V (Np,p+1,...,p+q) (x) = V (x; p + q, q). (38)
This shows in an explicite way what is exactly the chain of DBT which permit to build the general TDPT potential of integer parameters as a rational (in sin x) extension of the free particle system.

The Crum formula gives then

V (x; p + q, q) = -2 (ln W (ψ 0 , ..., ψ p-1 , ψ p+1, ..., ψ p-1+2j , ..., ψ p-1+2q | x)) ′′ , (39)
or, noting l = p + q, V (x; l, q) = q(q + 1) cos 2 x + l(l + 1) sin 2 x = -2 (ln W (sin x, ..., sin((l -q) x), sin ((l -q + 2)x) , ..., sin ((l -q + 2j)x) , ..., sin ((l + q) x) | x))

′′ ,(40)

which is precisely the Gaillard-Matveev result.

IV. FROM THE CONSTANT POTENTIAL TO THE BESSEL POTENTIAL AND A WRONSKIAN RAYLEIGH FORMULA.

To generate the TDPT potentials from the constant one, we start with a DBT built from any (physical) diffusion state of the constant potential. Nevertheless, we are not limited to this case and on a formal point of view we can use any eigenfunction, physical or not to build this DBT. In particular, we can choose as seed eigenfunction ψ 0 (x) = x = lim k→0 ψ k (x) whose associated eigenvalue E 0 = 0 is located at the lower boundary of the physical spectrum of V . The extended potential generated by the DBT A(w 0 ) is

V (0) (x) = 2/x 2 = V B (x, 1) , (41)
that is, the first Bessel potential of integer parameter, the general form of which being V B (x; a) = a(a + 1)/x 2 . The two partner hamiltonians H and H (0) are strictly isospectral but the levels of are no more degenerated since the strong singularity of V (0) at the origin imposes to retain in the spectrum of H only the eigenstates which satisfy the Dirichlet boundary condition at the origin

ψ 0 + = 0. ( 42
)
The physical spectrum of H (0) is then

E (0) k = k 2 ∈ ]0, +∞[ , k > 0 ψ (0) k (x) = 1 E k -E0 A(w 0 )ψ k (x) = (-sin(kx) + kx cos (kx)) /k 3 x, (43)
the unphysical "fundamental" eigenfunction for k = 0 being ψ

(0) 0 (x) = x 2 ∼ lim k→0 ψ (0) k (x).
We can start to build a chain of isospectral potentials by successive DBT based on the "fundamental" unphysical eigenfunctions of zero energy which are recessive at the origin (satisfying the Dirichlet boundary conditions at the origin). An immediate recurrence gives, for N m = (0, ..., 0)

m times ψ (Nm) 0 (x) = x m+1 V (Nm) (x) = V B (x; m) = m(m + 1)/x 2 . ( 44
)
All the extensions V (Nm) subject to the Dirichlet boundary condition at the origin are then strictly isospectral to V (0) with a unique energy band ]0, +∞[. These results are naturally well-known and this constitutes the first example of application that Darboux gave of his method in his first article on the subject [START_REF] Darboux | Sur une proposition relative aux équations linéaires[END_REF] . They are also one of the key ingredient given by Duistermaat and Grünbaum 21 in the proof that every trivial monodromy potential decaying at infinity can be considered as a rational extension of the free particle potential.

Nevertheless, the question of the Wronskian representation of the extended potentials Eq(44) and of their eigenfunctions is less trivial. Indeed, we are now in a case of confluency where all the successive DBT are built from eigenfunctions associated to the same energy. Then we cannot use the Crum formulas to express V (Nm) and its eigenstates ψ (Nm) k in terms of Wronskians of distinct "seeds functions" (note that A(w (Nm) 0

) is an annihilator for ψ )). In this case, we can however use Matveev's formulas [START_REF] Matveev | Generalized Wronskian formula for solutions of the KdV equations: first applications[END_REF][START_REF] Matveev | Positons: slowly decreasing analogues of solitons[END_REF] and express V (Nm) and ψ (Nm) k in terms of "generalized Wronskians": Matveev's formulas: In the confluent case where the repeated DBT are built on eigenfunctions associated to the same value k 0 of the spectral parameter k (N m = (k 0 , ..., k 0 ) ) m times , we can write

                 V (Nm) (x) = V (x) -2 ln W ψ κ , ∂ i 1 ψκ ∂κ i 1 κ=k0 , ..., ∂ i m-1 ψκ ∂κ m-1 κ=k0 | x ′′ ψ (Nm) k (x) = W ψ k 0 , ∂ i 1 ψκ ∂κ i 1 κ=k 0 ,..., ∂ i m-1 ψκ ∂κ m-1 κ=k 0 ,ψ k |x W ψ k 0 , ∂ i 1 ψκ ∂κ i 1 κ=k 0 ,..., ∂ i m-1 ψκ ∂κ m-1 κ=k 0 |x , (45)
where the ∂ i j ψκ ∂κ i j κ=k0

, j = 0, ..., m -1, are the m first non-zero derivative of ψ κ with respect to the spectral parameter κ at the value κ = k 0 . In our case, for k 0 = 0, we have

ψ k (x) = ∞ n=0 (-1) n x 2n+1 (2n + 1)! (k) 2n , (46)
The sequence of the i j is given by i j = 2j, j ≥ 0 with

∂ 2j ψ κ (x) ∂κ 2j κ=0 = (-1) j x 2j+1 2j + 1 (47)
and the Matveev formula for the potential gives

V (x; m) = 2 ln W x, x 3 , ..., x 2m-1 | x ′′ . (48)
Darboux's result is readily recovered by using standard properties of the Wronskians 26 which allow us to write

W x, x 3 , ..., x 2m-1 | x = x m dy dx m(m-1)/2 W 1, y, ..., y m-1 | y y=x 2 = x m(m+1)/2 2 m(m-1)/2 m-1 j=1 j!, (49)
that is,

V (Nm) (x) = V B (x; m) = m (m + 1) /x 2 .
More interesting is the case of the eigenfunctions. Since ψ

(Nm) k is a solution of ψ ′′ (x) + k 2 - m(m + 1) x 2 ψ (x) = 0, (50)
satisfying the Dirichlet boundary condition at the origin, it can then be written as

ψ (Nm) k (x) ∼ x 1/2 J m+1/2 (kx) = 2 π xj m (kx), (51)
where J m+1/2 (x) and j m (x) are the usual Bessel and spherical Bessel functions [START_REF] Abramowitz | Hanbook of mathematical functions[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF] . But from the second Matveev's formula we also have

This can be considered as a Wronskian version of the Rayleigh formula for Bessel functions [START_REF] Abramowitz | Hanbook of mathematical functions[END_REF] . Indeed, the same handling as before [START_REF] Muir | A treatise on the theory of determinants[END_REF] applied to the Wronskian at the numerator of Eq(52) gives W x, x 3 , ..., x 2m-1 , sin(kx)/k | x = 2 m(m+1)/2 x (m+1)(m+2)/2 W 1, y, ..., y m-1 , g(y) | y ,

where y = x 2 and g(y) = sinc(kx) = sin(kx)/kx. Then

V. CONCLUSION

In this letter we have shown how to recover in a simple way the content of the Gaillard-Matveev theorem which provides a Wronskian representation for the TDPT potentials. This is directly achieved by combining the Crum formula and shape invariance arguments applied to specific singular extensions of the constant potential. The same reasoning can be adapted to obtain Wronskian representations for the eigenfunctions of the Bessel potentials. This confluent case necessitates to employ Matveev's generalized Wronskians, obtaining then a Wronskian version of the Rayleigh formula.

In 3 , Gaillard and Matveev also consider the case of the discrete Darboux-Pöschl-Teller potentials (DDPT), for which they give Casoratian representation formulas. A treatment of these discrete systems by the preceding approach will be the object of further investigations.

  p+1 being the image of ψ p+1 (x) = sin((p + 2) x) by the chain associated to N p ψ (Np)

  W x, x 3 , ..., x 2m-1 , sin(kx)/k | x W (x, x 3 , ..., x 2m-1 | x).

W x, x 3

 3 , ..., x 2m-1 , sin(kx)/k | x = 2 m(m+1)/2 x (m+1)(m+2)/2

VI. ACKNOWLEDGEMENTS

I would like to thank R. Milson and D. Gomez-Ullate for stimulating discussions.


Céline Teulière

Eric Marchand

Direct 3D servoing using dense depth maps

This paper proposes a novel 3D servoing approach using dense depth maps to achieve robotic tasks. With respect to position-based approaches, our method does not require the estimation of the 3D pose (direct), nor the extraction and matching of 3D features (dense) and only requires dense depth maps provided by 3D sensors. Our approach has been validated in servoing experiments using the depth information from a low cost RGB-D sensor. Positioning tasks are properly achieved despite the noisy measurements, even when partial occlusions or scene modifications occur.

I. INTRODUCTION

Most of the robotic positioning tasks are achieved by estimating first the relative pose between the robot and the scene or the object of interest, and then using a position-based control scheme [START_REF] Wilson | Relative End-Effector Control Using Cartesian Position Based Visual Servoing[END_REF]. However, the pose estimation problem itself is complex in its general formulation. Also known as the 3D localization problem, this problem has been widely investigated by the computer vision community [START_REF] Dementhon | Model-based object pose in 25 lines of code[END_REF] [START_REF] Lowe | Three-dimensional object recognition from single twodimensional images[END_REF] but remains non-trivial for unknown environments. Using range data, a range flow formulation has been proposed [START_REF] Horn | Rigid body motion from range image sequences[END_REF] [START_REF] Gharavi | 3-D Motion Estimation Using Range Data[END_REF] to estimate the 3D pose of a mobile robot. Alternatively, the alignement of successive 3D point clouds using ICP [START_REF] Besl | A method for registration of 3-D shapes[END_REF] [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] has become a very popular method. Many variants have been proposed in the litterature [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF] and the development of the so-called RGB-D cameras attracted lot of attention on these methods [START_REF] Swadzba | A comprehensive system for 3D modeling from range images acquired from a 3D ToF sensor[END_REF] [16] [START_REF] Henry | RGB-D Mapping: Using depth cameras for dense 3D modeling of indoor environments[END_REF] [START_REF] Ra | KinectFusion: Real-time dense surface mapping and tracking[END_REF] in the recent years.

In this paper, we propose to perform robotic tasks without reconstructing the full 3D pose between the robot and its environment, but using a sensor-based servoing scheme, the considered data being directly the depth map obtained from a range sensor. Our approach is thus related to other sensor-based methods, such as image-based visual servoing (IBVS) [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], where a robotic task is expressed directly as the regulation of a visual error. In IBVS, the visual error is usually defined as the difference between a current and a desired set of geometric features (points, straigth lines, etc.) selected from the image, to control the desired degrees of freedom. Therefore, IBVS schemes usually require the extraction of visual features from image measurements, and their matching in successive frames. However, those steps, based on image processing techniques, are often considered as the bottleneck of visual servoing methods.

Recently, some work proposed to use all the image directly, without any extraction or matching step, by minimizing the difference between the current image and a reference image. This approach is referenced as photometric C. Teulière is with INRIA, Lagadic Project, Rennes, France E. Marchand is with Université de Rennes 1, IRISA, INRIA, Lagadic Project, Rennes, France firstname.name@irisa.fr visual servoing [START_REF] Collewet | Photometric visual servoing[END_REF]. However, since it is based on the luminance consistancy assumption, it is sensitive to illumination changes. In our work we propose to use the depth map obtained from a range sensor as a visual feature, without any feature extraction or matching step, and to control a robot with this feature directly. Our approach is thus both direct (without any 3D pose estimation) and dense (without feature extraction). To our knowledge this is the first time dense depth information is used in a direct servoing approach.

The remainder of this paper is organized as follows: our dense depth map servoing framework is described in section II. In section III we propose simple solutions to increase its robustness to some practical issues such as incomplete measurements or occlusions. Positioning experiments have been conducted to validate the approach. The results are given and discussed in section IV.

II. DIRECT DENSE DEPTH MAP SERVOING

This section presents the heart of our approach, i.e. how to control a robot using dense depth maps. We first introduce what we call a depth map and what it means to use it as a feature to regulate (section II-A). Then we derive the fondamental equations necessary to compute our control law (section II-B and II-C). In section II-D we underline the main differences between our approach and sparse 3D approaches.

A. Depth map sensing

There are multiple technologies of sensors capable of providing depth information, (or range). Most range sensors without contact are active, and based on the time of flight (ToF) principle: the idea is to send waves of known velocity and measure the time it takes them to go from the sensor and come back after reflection on the scene. This can be achieved by sending light pulses. Another approach would consist in using a modulated signal and measuring the phase shift. In each case, knowing the velocity of the sent signal, the depth information is derived (eg: Laser scans, sonars, radars, ToF or RGB-D cameras).

Another existing technology for active range sensing is based on structured light: known patterns (stripes, dots, ...) are projected onto the scene and the depth information is deduced from their deformation. This technology is used for instance in the recent Microsoft Kinect or Asus Xtion pro devices, based on PrimeSense technology [START_REF] Freedman | Depth mapping using projected patterns[END_REF].

Depth can also be measured with passive sensors such as cameras: by matching image features in two different views of a calibrated stereo rig, depth can be computed from geometry. The depth information is sparse when a finite set [22] or elevation rate constraint equation [START_REF] Horn | Rigid body motion from range image sequences[END_REF]. It is analogous to the brightness change constraint equation that is used in the computation of optical flow [START_REF] Horn | Determining optical flow[END_REF] and used in direct photometric visual servoing methods [START_REF] Collewet | Photometric visual servoing[END_REF].

From equation ( 5) we know that the temporal variation of the depth is:

∂Z ∂t = Ż - ∂Z ∂x ẋ - ∂Z ∂y ẏ, (6)
where Ż = dZ dt denotes the Z velocity expressed in the sensor frame.

Therefore, the interaction matrix L Z related to one depth value is expressed by:

L Z = L P Z - ∂Z ∂x L x - ∂Z ∂y L y . (7)
The matrices L x , L y defined such that ẋ = L x v and ẏ = L y v are the well-known interaction matrices of image point coordinates, given by:

L x = -1 Z 0 x Z xy -(1 + x 2 ) y (8)
L y = 0 -1 Z y Z -(1 + y 2 ) -xy -x , (9)
and L P Z is the interaction matrix related to the coordinate Z of a 3D point, such that Ż = L P Z v. It is given by:

L P Z = 0 0 -1 -yZ xZ 0 . ( 10
)
More details on the derivation of those interaction matrices can be found in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. The full interaction matrix L Z of size N ×6 corresponding to the entire depth map is thus the stack of the 1×6 matrices L Zi :

L Z =    L Z1 . . . L Z N    . (11)

D. Dense vs sparse 3D servoing

Depth information has already been used in position-based visual servoing. For example, [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] proposed to use the 3D coordinates (X, Y, Z) of a set of 3D points as features to be regulated in a proportional control law. In other words, the positioning task was expressed as the regulation of the feature

P = (X 1 , Y 1 , Z 1 , ..., X N , Y N , Z N ) to a reference feature P * = (X * 1 , Y * 1 , Z * 1 , ..., X * N , Y * N , Z * N ) corresponding
to the 3D coordinates of the set of points at the desired robot position. The interaction matrix related to a single 3D point is then given by [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]:

L P =   -1 0 0 0 -Z Y 0 -1 0 Z 0 -X 0 0 -1 -Y X 0   . ( 12
)
At first sight, the formulation of this kind of 3D feature

(X 1 , Y 1 , Z 1 , ..., X N , Y N , Z N
) can seem very close to the vector formulation Z = (Z 1 , ..., Z N ) that we defined in section II-B. However, a key difference with respect to our approach is that [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] uses a sparse set of 3D features. Consequently, in [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] a matching step is required to determine the feature values through the sequence, and the range flow equation ( 5), based on a smoothness assumption, does not hold in the sparse case. On the contrary, one of the key advantages of the method we propose, is that it does not require any feature extraction or matching step and uses directly the dense depth information from the range sensor thanks to the range flow equation.

III. PRACTICAL ISSUES AND ROBUSTNESS

IMPROVEMENTS

In the previous section, we presented our depth map based servoing method. When testing it, we found that this method was efficient in simulation sequences, with perfect data, but we had to face some practical issues in real conditions, in particular, in our case, using a Kinect sensor. This section presents the modifications we had to undertake in order to improve the robustness of the servoing task with respect to noisy and incomplete measurements (section III-A) and to scene perturbations and occlusions (section III-B).

A. Noisy and incomplete measurements

As illustrated in Figure 2-b the depth map acquired by a Kinect sensor is noisy and incomplete. In practice, we only considered the pixels for which a depth value was available both in the reference Z * and the current Z depth maps. This means that the number N of depth values in Z and [START_REF] Horn | Determining optical flow[END_REF], is inferior to the size of the depth map (320 × 240). In the experiments presented in this paper, about 80% of the total number of pixels could typically be used.

In addition, we reduced the noise by applying a simple 3 × 3 Gaussian filter on the depth maps, the convolution being computed only with the valid neighbors.

Similarly, the spatial gradient was computed using a simple 3 × 3 derivative kernel taking into account the valid neighbors only.

B. Occlusions and scene modifications

Another issue to take into account is the possibility of partial occlusions or scene modifications during the servoing process. To reduce the effect of such events on the task achievement, we use robust M-estimation [START_REF] Huber | Robust Statistics[END_REF] [START_REF] Comport | Statistically robust 2-D visual servoing[END_REF]. We thus introduce a modification of our task objective (1) allowing uncertain measures to be less likely considered or in some cases completely rejected. The new task error is given by:

e = D(Z -Z * ) ( 13
)
where D is diagonal weighting matrix given by: D = diag (w 1 , ..., w N )), the weights w i depending on their distance to the median of the error vector e according to a robust function [START_REF] Huber | Robust Statistics[END_REF]. Different functions are possible for the robust estimation. In practice, we used Tukey's estimator to completely reject the least likely values. Using ( 13), the new control law becomes [START_REF] Comport | Statistically robust 2-D visual servoing[END_REF]:

v = -λ(DL Z ) + D(Z -Z * ). ( 14
)
Experimental results using this control scheme are presented in the next section.

IV. EXPERIMENTAL RESULTS

In this section we present the results of our approach for positioning tasks. A Kinect sensor has been mounted on a ADEPT Viper robot (see Figure 3). In each experiment, the robot first acquires the reference depth map at the desired position. It is then moved to an initial position from which the control scheme is launched, aiming at going back to the desired one. A fixed gain λ = 2.5 is used in these experiments.

The first experiment illustrates the behavior of our system in a nominal case, namely with a static scene and no occlusions. The depth maps are acquired using the LibFreenect 1 driver through the ViSP library [START_REF] Marchand | ViSP for visual servoing: A generic software platform with a wide class of robot control skills[END_REF], with a resolution of 320 × 240. The scene is composed of various objects of different shapes and materials (Figure 5 1-a).

The initial and final states are illustrated in Figure 5. The first row shows the RGB views provided by the Kinect for the initial (1-a) and final (1-b) positions. Those images are never used in the control scheme but are useful for a better understanding of the setup. The depth maps are shown in the second row, and the last row gives the corresponding error, i.e. the difference between the desired and the current depth maps, unavailable data being discarded as explained in III-A. The difference images are scaled so that a plain grey frame (3-b) corresponds to a null error, and thus to the good achievement of the task. In Figure 5 (3-a) we see that the initial error was significant. The corresponding quantitative values for the task error, the 3D positioning errors and the velocities are given in Figure 6. Figures (b) and (c) show that with an initial error of 5 to 15cm in translation and 5 to 22deg in rotation, the positioning task is achieved with a remaining error of less than 3mm in translation and 0.4deg in rotation. Given the low depth resolution of the sensor and its noisy measurements, this corresponds to a good achievement of the task.

Note that in this scene the smoothness assumption was not verified since large depth discontinuities exist at the border of the objects, for example between the table and the floor. This experiment thus shows that the method is successful beyond its initial assumption.

In the second experiment, we evaluate the robustness of our approach with respect to partial occlusions or modifications of the observed scene. The initial scene is illustrated in Figure 8 (1-a). During the task achievement, someone entered the sensor field, removed an object and put it back several times. Some selected frames of this sequence are shown in Figure 8. The initial and final positions are illustrated in the first and last columns, while columns (b) and (c) show examples of occlusions. Note that at the end of the sequence the white bear has been completely removed from the scene, and the final depth map (Figure 8 (2-d)) is thus different from the desired one (Figure 7 (b)). This difference appears in the final difference image (Figure 8 (4-d)) and the task error function (Figure 7 (a)). However, despite the scene modifications and occlusions, the positioning task is successfully achieved, as shown by the convergence of the positioning errors in Figure 7 (b) and (c). The robustness of our control scheme to perturbations is the result of the use of M-estimation (see III-B). The effect of M-estimation is illustrated on the third row of Figure 8, where we represented the relative weights of each data in equation [START_REF] Lowe | Three-dimensional object recognition from single twodimensional images[END_REF]. Black pixels correspond to rejected values and brightest ones to inliers. Figure 8 (3-b), (3-c), and (3-d) show that the perturbations are correctly detected since the corresponding pixels are given a smaller weight. the positioning accuracy for this experiment is similar to the first one. The videos of these experiments are provided with this paper.

V. CONCLUSIONS

We have demonstrated that it is possible to use a dense depth map directly to achieve a robotic task. The main advantage of this approach is that it does not require any pose estimation, feature extraction or matching step. Moreover, when the depth map is obtained from an active sensor, the resulting approach is not sensitive to illumination changes as photometric approaches can be. Some limitations can appear with the use of active sensors such as Kinect RGB-D camera, in particular the noise and the absence of some measurements. We shaw however that those issues can be overcome thanks to the use of M-estimators and basic image pre-processing.

Fig. 4 .

 4 Fig. 4. The task error is the difference of depth maps Z -Z * .

Fig. 5 .

 5 Fig. 5. First experiment. First column corresponds to the initial position. The RGB view from the Kinect (1-a) is not used in the algorithm. (2-a) Initial depth map, where white parts correspond to unavailable data. (3a) Difference between the initial and desired depth maps. Second column corresponds to the end of the motion. The final depth map (2-b) corresponds to the desired one, since their difference (3-b) is a uniform grey.

Fig. 6 .

 6 Fig. 6. First experiment. (a) Task error, (b) translational part of positioning error, (c) rotational part of positioning error, (d) translational velocities, (e) rotational velocities.

Fig. 7 .

 7 Fig. 7. Second experiment. (a) Task error, (b) desired depth map, (c) translational part of positioning error, (d) rotational part of positioning error, (e) translational velocities, (f) rotational velocities.


Bertrand Delabarre
email: bertrand.delabarre@irisa.fr

Guillaume Caron
email: guillaume.caron@u-picardie.fr

Eric Marchand
email: eric.marchand@irisa.fr

Omnidirectional Visual Servoing using the Normalized Mutual Information

Keywords: Vision, visual servoing, mutual information, information theory, robotics

Geometrical features have been a key element of visual servoing for several years. Recently, different works have shown how using all the information within the images can lead to successful servoing tasks. In particular, works using mutual information have been proposed and tested on perspective cameras. In this paper we propose to adapt this approach to vision systems following the unified sphere model for central cameras using a normalized version of the mutual information. This will permit to apply the technique to large fields of view with a more reliable similarity function. Several experiments are performed on a cartesian robot with a fisheye camera to validate our approaches.

INTRODUCTION

Visual servoing uses the information provided by a vision sensor to control the movements of a dynamic system [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]; [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]; Chesi and Hashimoto (2010)). This approach requires to extract and track visual information (usually geometric features) from the image in order to design the control law. This tracking process is one of the bottlenecks in the development of visual servoing techniques.

Other works have tried to circumvent these problems by using directly the information provided by the entire image [START_REF] Nayar | Subspace methods for robot vision[END_REF]; [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]; [START_REF] Kallem | Kernelbased visual servoing[END_REF]; [START_REF] Collewet | Photometric visual servoing[END_REF]; [START_REF] Dame | Mutual information-based visual servoing[END_REF]). Features are no longer extracted from the image. Those works have began with [START_REF] Nayar | Subspace methods for robot vision[END_REF] and later on [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF], where the images were reduced to eigenspaces. Later works have used directly the whole images. In [START_REF] Collewet | Photometric visual servoing[END_REF], a control law was proposed that minimizes the error between the current image and the desired one. In that case the vector of visual feature is nothing but the image itself and the error to be regulated is the sum of squared differences (the SSD). This approach features many advantages: it does not require any matching or tracking process. Furthermore since the image measurements are nothing but the pixel intensity, there is no error in the feature extraction process leading to a very precise realization of the task. The method was later extended to omnidirectional camera in Caron et al. (2010). [START_REF] Kallem | Kernelbased visual servoing[END_REF] also considered the pixels intensity with a kernel-based method that leads to a highly decoupled control law. However, this approach can not control the 6 degrees of freedom of the robot and it is also very limited in the case of appearance variations.

As previously stated, considering image intensities [START_REF] Collewet | Photometric visual servoing[END_REF]) is quite sensitive to modification of the environment and more robust registration functions should be considered. To solve this problem, the considered approach does not use directly the luminance of the pixels but the information contained in the images. The visual feature is the mutual information defined by [START_REF] Shannon | A mathematical theory of communication[END_REF]. The mutual information (built from the image entropy) of two random variables (images) measures their mutual dependence. This function does not directly compare intensities of the two images but the distribution of the information in the images. Considering two images, the higher the mutual information (MI), the better the alignment between the two images is. Mutual information was considered for positioning and navigation task using visual servoing in [START_REF] Dame | Mutual information-based visual servoing[END_REF]. Since this is an entropy based measure, it is very robust towards illumination changes, occlusions or multimodality.

In this paper, we extend the approach proposed in [START_REF] Dame | Mutual information-based visual servoing[END_REF] in two directions. First, we consider a new expression for mutual information which is normalized. This measure, proposed by [START_REF] Studholme | An overlap invariant entropy mesure of 3d medical image alignment[END_REF] to perform image alignment, is bounded which makes it easier to interpret when used in an optimization process. Then, we extend the proposed method to the use of the unified model of central projection cameras (Barreto (2001)). Omnidirectional cameras allow to consider more information on the surrounding environment which improves the efficiency of the approach. This paper is organized as follows. First, direct visual servoing approaches are discussed, before describing our approach based on normalized mutual information. The adaptation to the use of the unified model of central projection cameras is then exposed. Finally, several experiments are exposed to validate our work.

DIRECT VISUAL SERVOING

Visual servoing tasks have used geometrical features for a long time [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]). Recently, works have proposed new direct approaches. They achieve a direct visual servoing using the information provided by the whole image, thus eliminating the need to extract and track features. In [START_REF] Collewet | Photometric visual servoing[END_REF] the authors used directly the pixel luminance whereas in [START_REF] Dame | Mutual information-based visual servoing[END_REF] a measure of mutual information was chosen.

Positioning Task

The aim of a positioning task is to reach a desired pose of the camera r * , starting from an arbitrary initial pose. To achieve that goal, one needs to define a cost function which will indicate whether or not the camera is going in the right direction. Most of the time this cost function f is a dissimilarity measure, function of the camera pose, which needs to be minimized by controlling the camera. Considering the actual pose of the camera r the visual servoing problem can therefore be written as an optimization process:

r = arg min r f (r, r * ) (1)
where r, the pose reached after the optimization (servoing process), is the closest possible to r * (optimally r = r * ). For example, considering a set of geometrical features s, the task will typically have to minimize the difference between s(r) and the desired configuration s* which leads to:

r = arg min r (s(r) -s * ). (2)
This visual servoing task is achieved by iteratively applying a velocity to the camera. This requires the knowledge of the interaction matrix (also known as image Jacobian) L s of s(r) that links the variation of ṡ to the camera velocity and which is defined as:

ṡ(r) = L s v
(3) where v is the camera velocity. This equation leads to the expression of the velocity that is applied to the robot as it is defined in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] by:

v = -λ L s + (s(r) -s * )
(4) where λ is a convergence factor tuned to improve convergence speed.

Using the mutual information

In [START_REF] Dame | Mutual information-based visual servoing[END_REF], it was proposed to use the mutual information defined by [START_REF] Shannon | A mathematical theory of communication[END_REF] as the cost function of the positioning task: r = arg max r MI (I(r), I * ) .

(5)

The mutual information can be defined as the quantity of information shared by two signals. It is an entropy-based measure and its main advantage is that it is very robust. Since it is based on the information itself rather than on its representation it is resistant to many variations as occlusions or illumination changes. That is the reason why it can be very useful, as the conditions when effecting a visual servoing task can change with time. The computation of mutual information between two images is given by the following equation: MI(I(r), I * ) = H(I(r)) + H(I * ) -H(I(r), I * ) (6) where H(I(r)) is the measure of entropy of the image I(r) and H(I(r), I * ) the joint entropy of the images I(r) and I * . They are defined as:

H(I(r)) = - Nc I i=0 p I (i) log (p I (i)) (7)
H(I(r), I * ) = -

Nc I i=0 Nc I * j=0 p II * (i, j) log (p II * (i, j)) (8)
where N c is the dynamic of the considered image (typically 255), p I (i)=Pr(I(x)=i) the probability distribution of i and p II * (i, j)=Pr(I(x)=i ∩ I * (x)=j) the joint probability distribution function. The problem of that cost function is that it has no fixed upper bound. If the values of entropy of I(r) and I * change, the maximum possible value of MI(I(r), I * ) changes thus making it impossible to compare values in different situations. Using a normalized mutual information, the task has an identifiable goal that it is trying to reach. This makes it easier to evaluate and makes comparison possible between different situations.

Using the Normalized Mutual Information

To get rid of the drawbacks of mutual information we propose in this paper to use a normalized measure of mutual information as defined in [START_REF] Studholme | An overlap invariant entropy mesure of 3d medical image alignment[END_REF]. It is expressed as:

NMI(I(r), I * ) = H(I(r)) + H(I * ) H(I(r), I * ) . ( 9
)
This measure is bounded. To prove it, let us first consider I(r) and I * not sharing any information. This results in

H(I(r), I * ) = H(I(r)) + H(I * ) leading to NMI(I(r), I * ) = 1.
For the upper bound, by definition:

H(I(r), I * ) ≥ max [H(I(r)), H(I * )] (10)
which gives:

1 H(I(r), I * ) ≤ 1 max [H(I(r)), H(I * )] . (11)
Knowing this, the upper bound of the NMI can be shown as:

NMI ≤ max [H(I(r)), H(I * )] (1 + β) max [H(I(r)), H(I * )] with β ≤ 1 (12) NMI ≤ 1 + β. ( 13
)
The best case scenario being I(r) and I * sharing all their information, therefore β = 1, equation ( 13) becomes:

NMI(I(r), I * ) = 2. ( 14
)

Shape of the cost function

The goal of this subsection is to see if the shape of the chosen cost function makes it suitable for an optimization process. In order to do that, a comparison between two images is made. On a reference image, a template of 100x100 pixels is extracted from a determined chosen starting position in the image. Then, a current image is chosen and a first patch is extracted beginning with a translation in the image of tx=-10 pixels and ty=-10 pixels with relation to the starting point of the template. The value of NMI between those two images of 100x100 pixels is then computed which gives the coordinate of the point (-10,-10) on the 3D plot. That process is iterated incrementing tx and ty until they both reach a value of 10 pixels. Figure 1 shows the different comparisons made when images are reduced to a dynamic of 8 grey levels. A first comparison was made to show the shape of the NMI in a nominal case. The function shows in that conditions a marked maximum and a smooth shape with no local minima, making the NMI fit for an optimization process in these conditions. Then, the target image underwent an illumination variation to see the effects on the NMI. Here again, the function is smooth and possesses a marked maximum when (tx,ty) equals to (0,0). This shows why the NMI can be used even if the light changes in the scene. Finally, the target image was occluded to show the effects of an occlusion on the NMI. The resulting shape is impacted by this large occlusion but even in those conditions no clear local maxima can be seen and the global maxima in (0,0) is still very marked which shows the suitability of the NMI for visual servoing in varying conditions.

Control law

e being the task function which needs to be regulated to zero, the control law giving the velocity to apply is defined as in [START_REF] Samson | Robot Control: The Task Function Approach[END_REF] by:

v = -λ L e + e ⊤ ( 15
)
where λ is a convergence factor and L e is the interaction matrix associated to the task. As the normalized mutual information needs to be maximized, the task is defined as the regulation to zero of the NMI derivative that is nothing but the interaction matrix L NMI related to NMI :

e = L NMI . (16)
The task and the velocity having the same dimension, the control law becomes, as in [START_REF] Dame | Mutual information-based visual servoing[END_REF]:

v = -λH -1 NMI L ⊤ NMI ( 17
)
where H NMI is interaction matrix of L NMI , hence the hessian of the normalized mutual information regarding the position of the camera. To compute those matrices let us express u(r) the numerator and v(r) the denominator of NMI as in equation ( 9):

u(r) = H(I) + H(I * ) (18) = - i,j p II * (i, j) log (p I (i)p I * (j)) v(r) = H(I, I * ) (19) = - i,j p II * (i, j) log (p II * (i, j))
That being set, L NMI and H NMI are obtained by chain rules of derivation:

∂u(r) ∂r = - i,j ∂p II * (i, j) ∂r log(p I (i)p I * (j)) (20) ∂v(r) ∂r = - i,j ∂p II * (i, j) ∂r (1 + log(p II * (i, j))) (21) ∂ 2 u(r) ∂r 2 = - i,j ∂ 2 p II * (i, j) ∂r 2 log(p I (i)p I * (j)) + 1 p II * (i, j) ∂p II * (i, j) ∂r 2 (22) ∂ 2 v(r) ∂r 2 = - i,j ∂ 2 p II * (i, j) ∂r 2
(1 + log(p II * (i, j)))

+ 1 p II * (i, j) ∂p II * (i, j) ∂r 2 (23)
With first order derivatives, L NMI can be expressed as:

L NMI (r) = ∂u(r) ∂r v(r) -∂v(r) ∂r u(r) v(r) 2 . ( 24
)
To simplify the expression of H NMI let us denote:

α(r) = ∂u(r) ∂r v(r) -u(r) ∂v(r) ∂r (25)
β(r) = v(r) 2 .
(26) This gives:

H NMI (r) = ∂α(r) ∂r β(r) -α(r) ∂β(r) ∂r β(r) 2 (27)
with:

∂α(r) ∂r = ∂ 2 u(r) ∂r 2 v(r) -u(r) ∂ 2 v(r) ∂r 2
(28)

∂β(r) ∂r = 2 ∂v(r) ∂r v(r). ( 29
)
To compute the derivatives of the probabilities, [START_REF] Dame | Mutual information-based visual servoing[END_REF] used B-spline functions φ to perform histogram binning. The joint probability used in equation ( 8) therefore becomes:

p II * (i, j, r) = 1 N x x φ i -I(r, x) φ j -I * (x) (30)
where N x is the wanted number of bins in the histogram and I(r, x) and I * are the images processed to contain only N x grey levels. Deriving this joint probability yields :

∂p II * (i, j, r) ∂r = 1 N x x ∂φ ∂r i -I(r, x) φ j -I * (x) (31) ∂ 2 p II * (i, j, r) ∂r 2 = 1 N x x ∂ 2 φ ∂r 2 i -I(r, x) φ j -I * (x) (32)
The Jacobian and Hessian of φ are finally given by:

∂φ ∂r = - ∂φ ∂i L I (33) ∂ 2 φ ∂ 2 r = ∂ 2 φ ∂ 2 i L T I L I - ∂φ ∂i H I ( 34
)
where L I and H I are respectively the interaction matrix and the hessian of I. They are given by:

L I = ∇I L x (35)
H I = L x T ∇ 2 I L x + ∇I H x ( 36
)
where ∇I is the gradient of I, L x is the interaction of a point [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]) and H x its hessian matrix.

CASE OF CENTRAL CAMERAS

Defining L I requires to compute the image gradient and the definition of the interaction matrix related to a point. To be able to consider a wide range of cameras, the unified model of central projection cameras (Barreto ( 2001)) is used. With such a projection model there exists various ways to compute these matrices.

Projection Model

Since the work of Barreto ( 2001), a unified projection model for central projection cameras was designed. This model describes a family of cameras from perspective to catadioptric ones with particular shape mirrors. Furthermore, [START_REF] Ying | Can We Consider Central Catadioptric Cameras and Fisheye Cameras within a Unified Imaging Model[END_REF] showed this model could be used for fish-eye lenses.

According to this model, a central projection camera can be modelled by a first projection on a sphere with coordinates (0, 0, ξ) in the camera frame followed by a perspective projection on the image plane. Such a model can be defined using parameter ξ which depends intrinsically on the catadioptric camera mirror parameters.

Knowing intrinsic parameters γ = {p x , p y , u 0 , v 0 , ξ}, a 3D point X = (X, Y, Z) is first projected on a unitary sphere and then in the image plane as x = (x, y, 1). The relationship between X and x can be expressed as:

x = pr γ (X) with    x = X Z + ξ √ X 2 + Y 2 + Z 2 y = Y Z + ξ √ X 2 + Y 2 + Z 2 (37)
x is the point on the virtual normalized plane and the image point in pixelic coordinates is obtained by:

u = u v 1 = p x 0 u 0 0 p y v 0 0 0 1 x y 1 = Kx (38)
As this paper also deals with methods using data lying on the equivalent sphere, let us define the spherical projection function of a 3D point:

X S = pr S (X) with          X S = X √ X 2 + Y 2 + Z 2 Y S = Y √ X 2 + Y 2 + Z 2 Z S = Z √ X 2 + Y 2 + Z 2 (39)
where X S = (X S Y S Z S ) T . On the contrary, the inverse projection function pr -1 γ allows to retrieve the point on the sphere corresponding to the spherical projection of the 3D point, from x:

X S = pr -1 γ (x) =        ξ + 1 + (1 -ξ 2 )(x 2 + y 2 ) x 2 + y 2 + 1 x ξ + 1 + (1 -ξ 2 )(x 2 + y 2 ) x 2 + y 2 + 1 y ξ + 1 + (1 -ξ 2 )(x 2 + y 2 ) x 2 + y 2 + 1 -ξ        (40)

Image Plane Visual Servoing (IPVS)

Considering image plane visual servoing (ie considering the point x on the image plane), gradients computation are the same as with perspective cameras. The interaction matrix however changes. It was shown in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] that the interaction matrix can be expressed as the product of two Jacobians:

L x = ∂x ∂X ∂X ∂r ( 41
)
where ∂x ∂X represents the movement of a point in the image with relation to the corresponding point in 3D and ∂X ∂r the movement of the 3D point with relation to the camera pose. The advantage of this formulation is that the second term of the product does not depend on the type of projection and therefore only the first part needs to be redefined. This is done by deriving (37). The interaction matrix is eventually expressed as in Barreto (2001).

Cartesian Spherical Visual Servoing (CSVS)

As described in section 3.1, the projection model first projects the points on a sphere before projecting them again, this time on the image plane. To compute gradients adapted to the omnidirectional image shape, one can therefore decide to compute them based on X S , the projection of X on the sphere. With this representation, L X S is defined as in [START_REF] Hamel | Visual servoing of an under-actuated dynamic rigid-body system: An image-based approach[END_REF]:

L X S = ∂X S ∂X ∂X ∂r = 1 ρ (X S X T S -I 3 ) [X S ] × (42)
where [X S ] × is the skew matrix of the vector X S .

Gradient computation

To evaluate the gradients which are necessary to compute ( 35) and ( 36) for CSVS, it was chosen to act as proposed in Caron et al. (2010). A sample step is first defined. Then each point is projected from the image to the sphere and the sampling is done to determinate its neighbors on the sphere. Those point are finally projected in the image to determinate the neighborhood to be used when computing the gradient. This step only has to be done once for the whole process and can therefore be done beforehand. As it provides real values, a bilinear interpolation was used to compute values at the neighbors' locations.

EXPERIMENTAL RESULTS

Several experiments were performed to validate the method. They were realized on a robot with six degrees of freedom with a fisheye camera mounted on its end effector. The camera provided 320x240 input images. The test program was implemented using the ViSP [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]) library. Three experiments are exposed in this section. The first one shows how the proposed method manages to complete a positioning task and compares the two modelizations seen in section 3. The second and third one respectively study the effects of illumination variations and of occlusions on the results of the positioning task. The initial and final position difference introduces a translation along X axes and a rotation along Y axes typically leading to projection ambiguities and a consequent translation along Y axes making it difficult for a servoing task to perform well. During the experiments, interaction matrix is computed at each iteration. Since it is an unknown parameter, computations are made assuming a constant Z for all the pixels. An adaptive gain was used to enhance speed when the camera is at the limits of its convergence cone. For readability purposes, graphics of NMI were created depicting the evolution of (NMI-1) instead of NMI, giving a measure between 0 and 1.

Experiment 1

The aim of this experiment is to compare the two representation methods exposed in section 3. To do so, both tasks were launched from the same sets of positions and the results of the servoing tasks were monitored (see figures 2 to 5). Let us first note, that in both cases, the positioning task is correctly achieved. Cost function and trajectory are more noisy with IPVS which shows an evolution of the NMI shaky when close from convergence. This is mainly due to the fact that more approximations are done in the interaction matrix computation (especially in the image gradient computation, see Caron et al. ( 2010)). In both cases repositioning precision was good, with final errors being around 0.1 millimeter.

Different experiments also showed CSVS to converge from further positions and being less affected by the modification of the scene or the depth approximation, therefore it is this method that will be used to demonstrate the robustness of the servoing task towards variations of the scene. The evolution of NMI on figure 4 is interesting as it can be related to the shape of the NMI function (see figure 1). Indeed, when far from the maximum, the function is rather monotonous and does not grow rapidly which makes the task advance slowly. On the other hand, when it comes close to the maximum the evolution becomes quicker and the NMI grows rapidly as it can be seen here around iteration 300. Also, the value of NMI does not reach its theoretical maximum of 2 (in the case of the evolution of (NMI-1) it should converge to 1). This is mainly due to the fact that to perform histogram binning Bspline functions were used, impacting the estimated probability distribution functions.

Experiment 2: light variations

The goal of this experiment is to assess the robustness of the method towards illumination changes. This was done using the CSVS approach and changing illumination conditions in several places of the room thanks to various spot lights. Even though the final positioning precision drops a little, averaging around 1 to 3 millimeters, the task is still successful and the robot stabilizes at the desired position even with distant starting points. The trajectory is moderately impacted, as the overall course of the camera stays the same even though small variations appear. Figure 6 also shows the evolution of the cost function and positioning errors which are very similar to the nominal case, which shows that the task is almost not impacted by the changes of conditions.

Experiment 3: impact of occlusions

This experiment was also designed to evaluate the robustness of the method, but this time with respect to occlusions (see fig 7).

In order to do that, several objects were moved and people walked in the sight of the camera during the servoing task.

The task succeeds and positioning precision averages around 0.5 to 2 millimeters. Here again the camera stabilizes with an acceptable final error and the trajectory is very similar to the nominal case, which shows that the task is almost not impacted by occlusions.

CONCLUSION

In this paper, a new way of achieving visual servoing using the information theory was presented. The method was detailed on an omnidirectional camera and different approaches taking into consideration the resulting projection model were detailed. Experiments were realized on a cartesian robot with six degrees of freedom and robustness towards occlusions and illumination variations were demonstrated. When compared to MI, NMI shows similar results but gives better readability of its results since the upper bound sets an identifiable goal to reach.

Fig. 1 .

 1 Fig. 1. Shape of the NMI cost function. Below the shapes are shown the compared patches when tx = 0 and ty = 0 pixel.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. On the upper line are the initial and desired images, respectively on the left and on the right. On the lower line are the initial and final error images, respectively on the left and on the right.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Trajectory of the camera in the 3D world. Blue curve represents the CSVS method and green curve the IPVS method.Reference image Initial image

  REFERENCESBarreto, J.P.Araujo, H. (2001). Issues on the geometry of central catadioptric images. In Int. Conf. on Computer Vision and Pattern Recognition. Hawai, USA. Caron, G., Marchand, E., and Mouaddib, E. (2010). Omnidirectional photometric visual servoing. In IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, IROS'10, 6202-6207. Taipei, Taiwan.

Fig. 7 .

 7 Fig. 7. Omnidirectional visual servoing using normalized mutual information: impact of occlusions.