Amaury Dame

Eric Marchand
email: eric.marchand@irisa.fr

Using mutual information for appearance-based visual path following

Keywords: Navigation, mutual information, visual servoing

In this paper we propose a new way to achieve a navigation task (visual path following) for a non-holonomic vehicle. We consider an image-based navigation process. We show that it is possible to navigate along a visual path without relying on the extraction, matching and tracking of geometric visual features such as keypoint. The new proposed approach relies directly on the information (entropy) contained in the image signal. We show that it is possible to build a control law directly from the maximisation of the shared information between the current image and the next key image in the visual path. The shared information between those two images is obtained using mutual information that is known to be robust to illumination variations and occlusions. Moreover the generally complex task of features extraction and matching is avoided. Both simulations and experiments on a real vehicle are presented and show the possibilities and advantages offered by the proposed method.

Introduction

In recent years, robot localization and navigation have made considerable progress. Navigation can be seen as the ability for a robot to move autonomously from an initial position to a desired one (which may be far away from the initial one). Thanks to sensor based navigation, we have seen autonomous robots in various challenging areas (from highways to deserts and even on Mars). Nevertheless the design of these autonomous robots usually relies on more than one sensor (camera, stereo sensors, lidar, GPS,...). In this paper, we propose a new method that demonstrates the capability of a mobile robot to navigate autonomously using the information provided by a monocular camera. Furthermore we will show that the proposed approach does not require any tracking nor matching process which is usually a bottleneck for the development of such an approach.

Most navigation approaches consider a (partial) 3D reconstruction of the environment [START_REF] Ohno | Autonomous navigation for mobile robots referring pre-recorded image sequence[END_REF][START_REF] Sim | Mobile robot localization from learned landmarks[END_REF][START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF][START_REF] Clemente | Mapping large loops with a single hand-held camera[END_REF][START_REF] Se | Mobile robot localization and mapping with uncertainty using scale-invariant visual landmarks[END_REF][START_REF] Furgale | Visual path following on a manifold in unstructured three-dimensional terrain[END_REF], leading to SLAM-like techniques. Such solutions are attractive, since the navigation task will be achieved using a classical posebased control of the robot in the metric space. Within this context, during a learning step the environment is reconstructed using bundle adjustment approaches [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF], Kalman/particle filters based approaches [START_REF] Clemente | Mapping large loops with a single hand-held camera[END_REF], visual odometry [START_REF] Furgale | Stereo mapping and localization for long-range path following on rough terrain[END_REF][START_REF] Furgale | Visual path following on a manifold in unstructured three-dimensional terrain[END_REF]. Despite the complexity of the underlying problem, SLAM has proved to be a viable solution to create accurate maps of the environment [START_REF] Clemente | Mapping large loops with a single hand-held camera[END_REF] [START_REF] Se | Mobile robot localization and mapping with uncertainty using scale-invariant visual landmarks[END_REF] even in large ones [START_REF] Frese | Closing a million-landmarks loop[END_REF]. In this context, the control of the robot during the navigation task is a well known problem and the main difficulties here are i) the complexity of the initial reconstruction step and ii) the matching of visual features ob-served during the learning step with current observations. With a monocular camera as unique sensor, these are mainly computer vision issues.

Another class of techniques relies on the definition of a visual path: the appearance-based approaches [START_REF] Matsumoto | Visual navigation using viewsequenced route representation[END_REF][START_REF] Burschka | Vision-based control of mobile robots[END_REF][START_REF] Remazeilles | Image-based robot navigation from an image memory[END_REF][33] [START_REF] Courbon | Autonomous navigation of vehicles from a visual memory using a generic camera model[END_REF][7] [START_REF] Zhang | Robust appearance based visual route following for navigation in large-scale outdoor environments[END_REF]. The trajectory is no longer described in the metric space but as a set of reference images. A 2D visual servoing step allows the robot to navigate from its current position to the next key image. When the robot gets close to this image, a new key image is selected. In this context, the environment can be modeled by a graph whose nodes are the key images. A visual path in the environment is nothing but a path in the graph [START_REF] Remazeilles | Image-based robot navigation from an image memory[END_REF]. Working directly in the sensor space, such approaches do not require prior 3D reconstruction step. In some cases, partial reconstruction has to be considered. In [START_REF] Segvic | A mapping and localization framework for scalable appearance-based navigation[END_REF] [START_REF] Diosi | Outdoor visual path following experiments[END_REF][START_REF] Booij | Navigation using an appearance based topological map[END_REF] a part of the epipolar geometry that links the current and key images is considered in order to predict the location of currently not visible features and ensure a robust tracking. In [START_REF] Courbon | Autonomous navigation of vehicles from a visual memory using a generic camera model[END_REF] homography computation wrt. the reference images allows to precisely localize the robot. In any case, the learning step of these appearance-based approaches is far less complex since it does not require any prior 3D reconstruction.

Nevertheless at navigation level, for both pose-based or most of the image-based visual navigation approaches, features have to be extracted or tracked in the image stream and matched with either the 3D database or key images to design the control law. Robust extraction and real-time spatio-temporal tracking or matching of these visual cues are non trivial tasks and also one of the bottlenecks of the expansion of visual navigation. It could be then very interesting to consider directly a comparison with current image and next keyframe to control the robot. In [START_REF] Zhang | Robust appearance based visual route following for navigation in large-scale outdoor environments[END_REF], this relative orientation is computed in the Fourrier space and the control law is directly computed from the orientation difference. Similar approach is proposed by [START_REF] Matsumoto | Visual navigation using viewsequenced route representation[END_REF] but used cross-correlation. In [10][9], it has been shown that no other information than the image intensity (the pure image signal) need to be considered to control the robot motion and that these difficult tracking and matching processes can be totally removed. Although very efficient, this approach is sensitive to light variations and, thus, can hardly be considered in outdoor environment. In this paper, we propose a new approach that no longer relies on geometrical features nor on pixels intensity [START_REF] Collewet | Photometric visual servoing[END_REF] but uses directly the information (entropy) contained in the image signal as proposed in [START_REF] Dame | Entropy-based visual servoing[END_REF][START_REF] Dame | Mutual information-based visual servoing[END_REF]. More precisely we will consider mutual information [START_REF] Ce | A mathematical theory of communication[END_REF][39] as a similarity criterion.

MI has been introduced in the context of information theory by Shannon [START_REF] Ce | A mathematical theory of communication[END_REF]. It has been later considered as an image similarity measure back in the mid ninety's independently by Collignon [START_REF] Collignon | Automated multi-modality image registration based on information theory[END_REF] for tomographic image registration, Studholme [START_REF] Studholme | Automated 3D registration of truncated mr and ct images of the head[END_REF] for MR and CT image, an by Viola [START_REF] Viola | Alignment by maximization of mutual information[END_REF] for projection image. Since then MI has become a classical similarity measure especially for multi-modal registration techniques [START_REF] Pluim | Mutual-informationbased registration of medical images: a survey[END_REF] (eg, for medical or remote sensing applications). Being closer to the signal, we will show that this approach is robust to very important illumination variations and robust to large occlusions.

Our goal is then to propose a control law that allows the robot to maximize the mutual information between the current acquired image and the next image in its visual path. This is an optimization process. We show that it is possible to compute the interaction matrix that relates the variation of the mutual information to the vehicle velocity leading to the definition of the control law. Let us emphasize that since mutual information is computed from the whole images (current and key images) it is possible to directly control the motion of the vehicle along a given path without any feature extraction or matching. Furthermore no 3D reconstruction of the environment is necessary.

We will demonstrate the efficiency of this new approach on a navigation task carried out at 0.5 m/s over 400 meters. Images are acquired at 30Hz (nearly 25.000 images were acquired and processed in real-time during this navigation task).

In this paper we will first present a general overview of the method with the learning and the navigation steps. Then section 2 and 3 will focus on the two parts of the navigation steps which consists of the visual servoing task and the key images selection task. Finally, simulated and experimental results are presented in section 4.

Navigation process overview

In this work, we consider a non-holonomic robot with a camera mounted on the front. Our goal is not to localize the robot within its environment (visual odometry) but only to ensure that it is able to reproduce a visual path defined as a set of images previously acquired by the camera.

Learning step: definition of the visual path

With respect to previous approaches that rely on 3D reconstruction (eg, [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF]) or even on appearance-based approaches [START_REF] Segvic | A mapping and localization framework for scalable appearance-based navigation[END_REF], the learning step of the approach is simple. It does not require any feature extraction nor scene reconstruction: no image processing is done, only raw images are stored. The vehicle is driven manually along a desired path. While the vehicle is moving, the images acquired by the camera are stored chronologically thus defining a trajectory in the image space. Let us call I * 0 , . . . , I * N the key images that define this visual path.

Navigation step: following the visual path

The vehicle is initially positioned close to the initial position of the learned visual path (defined by the image I * 0 ). The navigation is performed using a visual servoing task. Figure 2 shows the general control scheme used for the navigation. In [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF][33] [START_REF] Chen | Qualitative vision-based path following[END_REF] the considered control scheme are either posebased control law or consider classical visual servoing process based on the use of visual features extracted from the current and key images (I and I * k ). In this work the definition of a new control law is proposed. One of the originality of this work is that, rather than relying on features extraction and tracking, we build the control law directly from the information shared by I and I * k measured using the mutual information [START_REF] Ce | A mathematical theory of communication[END_REF]. When the mutual information between two images is maximized, the two images are similar. We then control the robot in order to maximize the mutual information between I and I * k . As for any visual servoing scheme it is then necessary to exhibit the Jacobian that links the variation of the mutual information to the control input of the robot (that is the steering angle ψ or the camera rotational velocity ρ) needed to follow the path with a constant translational velocity v. This process is presented in the next section. In the same time, when the vehicle reaches the neighboring key image I * k , a new one I * k+1 is selected in the visual path. To achieve a seamless switching between key images, a specific process described in section 4 is proposed.

Mutual information based navigation

In [START_REF] Dame | Mutual information-based visual servoing[END_REF], it has been shown that it is possible to achieve 6 degrees of freedom (d.o.f) visual servoing task using only the information contained in the images acquired from the camera mounted on a robot and one reference image. The desired position of the robot is reached by maximizing the mutual information between the two images. Since mutual information is robust to illumination variations and occlusions, the use of mutual information-based visual servoing is well suited for outdoor navigation tasks.

Mutual information

In this section, a brief reminder of the definition of mutual information is given. Mutual information is the information shared by two signals (here, images). For the two signals X and Y, mutual information is given by the following equation [START_REF] Ce | A mathematical theory of communication[END_REF]:

MI(X, Y) = H(X) + H(Y) -H(X, Y) (1)
where H(X) denotes the entropy of the signal X, that means its variability. H(X, Y) denotes the joint entropy of the signals X and Y, that is the joint variability of the system defined by the two signals. By substracting the joint variability from the variabilities, as in equation 1, we obtain the shared information of the two signals that is mutual information.

In the present work we focus on mutual information between two images. The desired image (in the navigation context, desired image is the current key image of the visual path) is noted I * and the current image acquired by the vehicle camera is noted I. To address this definitions, let us consider that I is now a random variable and that the actual pixel intensities are samples of this random variable (I(x) being the intensity of the pixel x).

The entropy H(I) is a measure of variability of a random variable I. If i are the possible values of I(x) (i ∈ [0, N c I ] with N c I = 255) and p I (i) = Pr(I(x) = i) is the probability distribution function of i, then the Shannon entropy H(I) of a discrete variable I is given by the following expression:

H(I) = - N c I i=0 p I (i) log (p I (i)) . (2)
The formulation can be seen as follows: sincelog (p I (i)) is a measure of the uncertainty of the event i, then H(I) is a weighted mean of the uncertainties. H(I) is then the variability of I. Since a sample of I is in our case given by the pixel intensities I(x), the probability distribution function can be estimated using the normalized histogram of this image. The entropy can therefore be considered as a dispersion measure of the image histogram.

Following the same principle, joint entropy H(I, I * ) of two random variables I and I * can be defined as the variability of the couple of variables (I, I * ). The Shannon joint entropy expression is given by:

H(I, I * ) = - N c I i=0 N c I * j=0 p II * (i, j) log (p II * (i, j)) (3)
where i and j are respectively the possible values of the variables I and I * , and The original definition of mutual information given in [START_REF] Ce | A mathematical theory of communication[END_REF] [START_REF] Viola | Alignment by maximization of mutual information[END_REF] can thus be used as:

p II * (i, j) = Pr(I(x) = i ∩ I * (x) = j)
MI (I, I * ) = i, j p II * (i, j) log p II * (i, j) p I (i)p I * ( j) . (4)
The probabilities p II * , p I and p I * involved in the computation of MI are obtained by normalizing the histograms and joint his-togram of the images. Their analytical formulations are given by:

p II * (i, j) = 1 N x x φ (i -I(x)) φ ( j -I * (x)) (5)
p

I (i) = 1 N x x φ (i -I(x)) (6)
p I * ( j) = 1 N x x φ ( j -I * (x)) (7)
where N x is the number of points x in the region of interest (the complete image in our case). φ(ξ) is the function used to fill the histogram. Typically p I * ( j) is incremented each time I * (x) = j.

Then φ(ξ) = 1 if ξ = 0 and null otherwise. The equation of p I (i) is similar to the one of p I * ( j).

Since the two images are typical 8 bits images with 256 gray level values, the initial definition of mutual information is given for 256 entries for i and j. This definition gives a cost function that is subject to noise, artifacts and then local extrema that may induce issues in a non-linear optimization process [START_REF] Panin | Mutual information-based 3D object tracking[END_REF].

Several modifications on the computation allow to have a smooth cost function with a large convergence domain. The first is to consider a smaller number of entries for each histograms. The effect is to smooth the extremum and enlarge the convergence domain. To do so, image intensities used in equation ( 5), ( 6) and ( 7) are scaled to fit in the new number of bin N c . Let us note I and I * respectively the scaled images I and I * :

I(x) = I(x) N c -1 255 I * (x) = I * (x) N c -1 255 . ( 8
)
The intensities of this images are no more integer values. Thus, the original φ function has to be modified to update the histograms entries using real values (e.g. to compute φ i -I * (x) ). A solution of this problem is given by the Partial Volume Interpolation [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF] that defines φ as a first order B-spline (corresponding to a bilinear interpolation). An other operation is obtained by improving this interpolation. Instead of using a simple bilinear interpolation to compute the histograms (φ = B 1 ), B-splines of higher order are used [START_REF] Dame | Mutual information-based visual servoing[END_REF]. In the present work, we consider φ = B 3 (this will also be necessary to compute the second order derivatives of mutual information). B-splines functions are represented in Figure 3.

0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 B 1 B 2 B 3 (a) -1 -0.5 0 0.5 1 -3 -2 -1 0 1 2 3 B 1 B 2 B 3 (b)
Finally, the images I and I * can be smoothed (using for example a gaussian filter) to increase the domain of convergence by smoothing mutual information [START_REF] Pluim | Mutual-informationbased registration of medical images: a survey[END_REF]. Moreover an interesting effect of the image filtering is a convexification of the cost function. The filter variance has however to be controled to keep the accuracy of the maximum.

Combining the three operations, the obtained MI is smooth with a wide and accurate maximum and thus adapted for the optimization problem that will be used in the next section. Figure 4 shows the results obtained using the comparision between the original computation (φ = B 1 and 256 probability bins) and the proposed computation.

Navigation using Mutual information 3.2.1. Visual Servoing using Mutual information

In the general visual servoing formulation, the goal is to minimize a dissimilarity measure (generally the difference) between some desired and current features using a non linear minimization [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. Such approaches have been already used in navigation [START_REF] Segvic | A mapping and localization framework for scalable appearance-based navigation[END_REF] [START_REF] Chen | Qualitative vision-based path following[END_REF]. Here, our goal is to propose a more direct approach that uses the image as a whole and that does not rely on geometric features and, hence, avoid the tracking and matching steps. Rather than minimizing, as usual, the error between current and desired features, the goal is to maximize the amount of information shared by the two images.

Regarding the optimization of mutual information, various approaches have been proposed: first-order gradient descent [START_REF] Viola | Alignment by maximization of mutual information[END_REF], multi-resolution hill climbing algorithm [START_REF] Studholme | Automated 3D registration of truncated mr and ct images of the head[END_REF] or simulated annealing techniques [START_REF] Ritter | Registration of stereo and temporal images of the retina[END_REF]. Powell's [START_REF] Collignon | Automated multi-modality image registration based on information theory[END_REF][START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF] or Simplex [START_REF] Meyer | Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin-plate spline warped geometric deformations[END_REF][START_REF] Chen | On registration of regions of interest (roi) in video sequences[END_REF] methods (which do not require function derivatives to be analytically expressed) have been very popular in MI optimization but the former is sensitive to local optima or computationally inefficient and is not adapted to our navigation problem. Considering that MI computation is evaluated from the joint image intensity histogram, an analytic derivative of the mutual information is difficult to obtain. In order to compute MI derivatives, [START_REF] Maes | Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information[END_REF] introduces partial volume interpolation for the construction of the joint histogram leading to an analytic computation of MI gradients. In [START_REF] Thévenaz | Optimization of Mutual Information for Multiresolution Image Registration[END_REF], the authors formulate the mutual-information criterion as a continuous and differentiable function of the registration parameters using B-Spline Parzen windows. These derivatives can then be considered within a Newton like approach. Such approach has been considered here. Nevertheless we propose a dedicated control law specifically adapted to the MI cost function. We propose an inverse compositional optimization approach [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF] where an important part of the required derivatives can be precomputed, resulting in small computation times. A precise, complete and efficient computation of the Hessian matrix [START_REF] Dame | Second order optimization of mutual information for real-time image registration[END_REF] is also described.

Considering that our vehicle has a constant translational velocity v, we will control only the vehicle steering angle ρ. The navigation task toward the next key image I * can then be seen as a gradient descent optimization process where the cost function is defined as the mutual information between I and I * wrt. the angle ρ :

ρ * = arg max ρ f (ρ) with f (ρ) = MI(I * , I(ρ)). ( 9
)
This maximization is performed by updating the parameter ρ to find a null derivative of mutual information using a non linear optimization. Using the parameter update δρ, the expression to maximize can typically be rewritten as :

MI(I * , I(ρ t+1 )) = MI(I * , I(ρ t + δρ)). ( 10
)
where ρ t is the current steering angle and ρ t+1 is the steering angle at the next iteration. As it is demonstrated in [START_REF] Baker | Equivalence and efficiency of image alignment algorithms[END_REF], this formulation is equivalent to the inverse compositional formulation, where the expression to maximize is the following one:

f = MI(I * (-δρ), I(ρ t )). ( 11
)
As it will be explained later, this formulation (also classical in visual servoing [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]) allows to precompute some terms and then have a faster computation. Using a classical Newton's method, the parameter update can be computed using:

δρ = -αH -1 MI L MI ( 12
)
where L MI and H MI are respectively the interaction matrix of MI and its Hessian matrix. α ∈ [0, 1] is a scalar factor that allows to set the speed of the convergence (in our navigation task the goal is not to get to the maximum in one iteration).

Since I * is now depending on ρ, the expressions defined in equations ( 6) and ( 5) depends on ρ too. The derivatives of mutual information are then given by the following equations:

L MI = i, j L p II * 1 + log p II * p I * (13)
H MI = i, j L p II * L p II * 1 p II * - 1 p I * +H p II * 1 + log p II * p I * ( 14
)
Although it is classical to consider the second term of ( 14) as null [17][37], in this work the computation of the exact second derivative of mutual information is used. The two previous expressions depend on the joint probability derivatives. Using equation ( 5), the formulation of the the joint probability p II * between the key image I * and current image I and its variations are given by:

p II * (i, j, ρ) = 1 N x x φ i -I(x, ρ) φ j -I * (x) L p II * (i, j,ρ) = 1 N x x L φ i-I(x,ρ) φ j -I * (x) H p II * (i, j,ρ) = 1 N x x H φ i-I(x,ρ) φ j -I * (x)
where N x is the number of pixels considered in the images I and I * . φ is a B-spline function differentiable twice (see the MI definition in the previous section). The interaction matrix and Hessian of φ are given by:

L φ i-I(x,ρ) = - ∂φ ∂i L I (15)
H φ i-I(x,ρ) = ∂ 2 φ ∂i 2 L I L I - ∂φ ∂i H I (16)
and:

L I = ∇I L x (17)
H I = L x ∇ 2 I L x + ∇IH x (18)
where ∇I and ∇ 2 I are respectivelly the gradient and the second order gradient of the image. Since the only degree of freedom considered is the rotation around the vertical axis (the y axis of the camera), the interaction and Hessian matrices are then:

L x = -(1 + x 2 ) -xy (19)
and

H x = 2x(1 + x 2 ) y(1 + 2x 2 )] (20)
The final update δρ values have been computed on the example presented in Figure 4. Figure 5 shows the value of the derivatives of mutual information and the corresponding value δρ depending on the rotation between the desired and the current position. The relation between the real rotation and the computed update is quasi linear. The result of this proposed update will then cause a quasi exponential decreasing of the error, that is the ideal goal of typical visual servoing tasks [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. Controlling both steering angle and translational velocity.. Increasing the number of d.o.f to be controlled is quite straightforward. In [START_REF] Dame | Mutual information-based visual servoing[END_REF], experiment with 6 controlled d.o.f have been reported. In our case if two d.o.f have to be considered (translational velocity v and steering angle ρ the main difference w.r.t the presented approach is the computation of the control law that is then given by:

v δρ = -αH -1 MI L MI ( 21
)
Such a control law would allow forward and backward motions of the vehicle (although not point stabilization). Formulation of L MI and H MI is not modified and equation ( 13) to [START_REF] Frese | Closing a million-landmarks loop[END_REF] remain valid (although size are modified, L MI is now a 2 × 1 vector and H MI is a 2 × 2 matrix). Nevertheless it is necessary to modify equation ( 19) and [START_REF] Furgale | Visual path following on a manifold in unstructured three-dimensional terrain[END_REF]. L x given by

L x = x/Z -(1 + x 2 ) y/Z -xy (22)
and the complete derivation of the Hessian matrices can be found in [START_REF] Lapresté | A Hessian approach to visual servoing[END_REF]. Let us note that this formulation of the problem required a knowledge of depth information Z in [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF]. Furthermore since Z can be large in the case of outdoor scene this may lead to instabilities in the control law. This is why we have considered only the steering angle in this paper.

Navigation using visual servoing

For every acquisition of an image I, an update δρ is computed in order to move the camera and increase the mutual information between I and I * . Using the model of the vehicle it is possible to go back to the steering angle that will give the estimated update.

Firstly the update of the rotation of the camera is linked to the camera rotational velocity ρ (around the y axis) by the following equation:

ρ = δρ ∆t ( 23
)
where ∆t is the processing time (30Hz in our case).

The velocity is directly linked to the steering angle ψ of the wheels. Using the model of the non-holonomic vehicle used in our experiments (See the car-like model Figure 6) the general steering angle is computed as follows:

ψ = arctan L ρ v = arctan L δρ v ∆t ( 24
)
where v is the translational velocity of the vehicle (along the z axis) and L is the distance between the front and rear wheels.

Key images switch in the visual path

The previous section shows how to control a vehicle toward one key image using information shared by the current and key images. To be able to follow the learned trajectory a switching process between the key images of the visual path has to be defined.

Various switching solutions

Several solutions can be considered. One could propose to simply analyse the cost function evolution and check if the function (equation ( 9)) is increasing since mutual information is supposed to increase in the visual servoing task. If it is no longer the case it could mean that the key image is outdated. This solution is unfortunately limited to nominal or simple cases. In an outdoor environment such conditions are unpractical: if the illumination is changing then the cost function value will be affected and the desired image selection will fail.

Another solution could be to consider only the rotation required (given by the parameter update δρ) to reach the alignment position. If the computed rotation is smaller than a given threshold it could mean that the vehicle is next to the desired position, then the next key image can be loaded. But such a simple solution will obviously fail when the tracked trajectory is a straight line.

Proposed key image selection process

The proposed approach is based on this solution coupled with a translation estimation. The key image is updated each time the remaining rotational error is low and that the translation to reach the desired position is small (so that there is no more problem in straight lines).

The rotational error is directly given by the parameter update δρ. However no estimation is given on the translational error (that is the remaining distance between the current position and the position corresponding to the key image).

The approach is to consider the variation of the mutual information between the key image and the current image depending on the variation of translational position t z of the vehicle. If the variation of MI is null, it means that the maximum of MI is reached and that the robot is at the desired position. If the variation is positive (respectively negative) it means that mutual information is increasing (respectively decreasing) and that the robot is getting close to (respectively moving away from) the desired position.

This variation is simply computed as the derivation of mutual information wrt. the translational velocity v z of the camera. The formulation of the problem is similar to the one proposed in section 3.2.1 as the difference that the current image is now depending on v z . The derivative of the mutual information is now expressed with the interaction matrix corresponding to the translational degree of freedom that is L t z = [x/Z y/Z] t with Z the depth of each points. Since an accurate estimation of the translation is not needed, Z is approximated to be constant with Z = 20 meters.

To validate the proposed approach, some simulations have been performed using a strongly rough environment. Figure 7 illustrates the performed experiment. The value of the derivative of mutual information is shown depending on the translation between the current and the key positions along the z axis. We can see that the choice of the depth value is not critical (in fact using the previous equations, it can be seen that changing Z is only modifying the derivative by a scale factor). Considering a strongly non flat scene, mutual information derivative with respect to the translation remains accurate with a null value when the robot reaches the desired translation.

Using two given thresholds on both the parameter update and the translation estimation allows to update the key image each time the robot is close to the current desired position.

Experimental results

This section presents navigation experiments performed with the vehicle represented in Figure 6 using the mutual information-based navigation process.

Simulation

The first experiment is a simulation that describes the behavior of the proposed navigation task in nominal conditions. Since this is a simulation, the acquired trajectory and the resulting one are perfectly known. The simulation is performed in an urbanlike environment that is shown in Figure 8(a). The ground is flat and there are no occlusions nor illumination variations between the environment in the learning and in the navigation phases. The buildings of the environment have various type of textures with low, high and repetitive textures. The experiment has been Considering the steering angle sent to the vehicle (red plot on Figure 9), we can see that the control law is obviously not continuous. Each time the key image is changed, the computed rotation is abruptly changed and then decreases exponentially. The effect on a real vehicle may be hard to endure. To solve this issue, we propose to filter the previous result to have smoother changes of the steering angle. A simple Kalman filter with a constant velocity model has been applied to the computed angle. The result of the Kalman filter is shown on previous computed steering angles (see Figure 9). This result is adapted for the control of the vehicle that keeps on following properly the path with smoother changes in its direction.

Navigation in natural environment

The experiments have been carried out in real-time using a camera mounted on an electric car-like robot named Cycab (see Figure 10) . The presented experiments consider realistic physical paths for which no common landmarks are visible from the initial and the desired position. The navigation process has been achieved with the same vehicle/camera that were used during the learning step (construction of the visual path). Therefore potential camera distortions were not an issue in the computation of MI between two images. All computation are performed in the normalized space. Nevertheless, if two cameras are considered, then it will be necessary to undistort both images (using known/estimated calibration parameters).

The mutual information navigation scheme has been tested on a non-holonomic vehicle (see Figure 6) in an outdoor envi- ronment. The final approach presented in the previous paragraph has been used. Let us emphasize that the vehicle is equipped with a monocular camera and that no other sensor such as GPS, radar or odometry are considered in these experiments. Furthermore, the 3D structure of the scene remains fully unknown during the learning and navigation steps. Aerial views of the environment, where the navigation task takes place, are shown in Figure 11 along with the considered trajectory (about 400 meters). As seen on the pictures, the environment is semi-urban with both trees and buildings (with windows acting as repetitive textures). Let us note that the vehicle crosses a covered parking lot (green part of the trajectory in Figure 11) and that the ground is no longer perfectly flat (mainly in the first 100 meters of the trajectory). When learning the path the vehicle is manually driven at a roughly constant velocity. For this experiment we consider 1200 key images (that is around three key images per meter). The navigation task itself is carried out at 0.5 m/s. Images are acquired at 30Hz (nearly 25.000 images are acquired and processed in real-time during this navigation task). Some pictures of one navigation task are shown in Figure 12. By comparing the current and key images (and the image error on the third row), we can see that the robot is qualitatively (as defined in [START_REF] Remazeilles | Image-based robot navigation from an image memory[END_REF]) following the same path. The navigation task has been tested with both cloudy and sunny weather using the same learned visual path). Since time had passed between the acquisition of the visual path and the navigation task, there have been very large illumination changes between the current and the key images as it is highlighted in Figure 14(a). The task has even been tested with a ground recovered with snow still using the same initial visual path (See Figure 14(b)). Despite those illumination variations the navigation task was still converging. That shows the robustness of the proposed control law to illumination variations and the efficiency of considered mutual-information similarity criterion to perturbation. Let us also note that the road on which experiments have been conducted is not completely flat. As can be seen on the video, just entering under the cover parking, the road features a quite important slope (at least 1 meter down in 5 meters), the camera is then pitching at this point.

Since no obstacle avoidance process is considered, the navigation task has been performed in quiet conditions. Nevertheless several vehicles have overtaken our experimental vehicle and appeared in the camera view. Despite this perturbation, and thanks to the robustness of the similarity criterion, the navigation task has never failed showing the robustness of mutual information to occlusions. One of these moments is represented in Figure 14(c) (the van in the current image was not present in the key image). The complete video is provided with this paper. Let us note that MI between the current image and the next key image can be monitored in order to detect failures (deviation from the visual path). In nominal conditions (few scene modification), if this MI is smaller than the MI between two successive key frames it is possible that the robot deviates from its visual path and security actions has to be considered.

Conclusion

In this paper, we have presented a new way to achieve image based navigation task. We show that our vehicle is able to track a previously learned trajectory using only the information provided by a monocular camera. The navigation task can be achieved despite important variation in the lighting condition and possible perturbations. This can be achieved thanks to various elements related to the use of mutual information which is new in vision-based control:

• Our approach does not rely on features extracted from the image. Therefore we do not need to track or match features (eg, keypoints) which has proved to be a difficult and not always reliable process. Furthermore no 3D information related to the scene structure is required.

• To avoid this tracking and matching processes, the visionbased control law of the non-holonomic vehicle is directly linked to the optimization of a similarity criterion based on the information shared by two images.

• We propose a control law that directly links the vehicle motion to the variation of the mutual information. The proposed approach that considers a derivation of the MI, up to the second order, allows a large convergence domain along with fast (video rate) computation.

• Considering information contained in the images and not features extracted from the image or the image intensities induces a natural robustness to perturbation that is essential in our navigation context.

We also considered a key images selection process which is efficient regarding the considered navigation taks. Future work is planned on the post processing of the steering angle update. For the moment, it is simply filtered, while using the knowledge on the vehicle kinematics and smoothness of the trajectory could improve the navigation task. Moreover, future work will devoted to the definition of efficient navigation task that require more degrees of freedom and more complex control models such as aerial drones.

Figure 2 :

 2 Figure 2: Navigation based on multiple visual servoing tasks.

Figure 1 :

 1 Figure 1: Key images that define the visual path. This visual path is learned prior to the navigation step.

  is the joint probability distribution function. Here, I and I * being images, i and j are the pixel intensities of the two images and the joint probability distribution function is a normalized bidimensional histogram of the two images. As for entropy, joint entropy measures the dispersion of the joint histogram of I and I * .

Figure 3 :

 3 Figure 3: B-splines from order one to three (a) and their derivatives (b).

Figure 4 :

 4 Figure 4: Method used to compute mutual information with respect to the rotation around the vertical axis ρ. (a) External view of the camera at the desired position in red (used to acquire the desired image I * ) and at the various positions of rotation in green. (b) shows the computed mutual information between I(ρ) and I * with respect to the rotation ρ computed using various formulation.

Figure 5 :

 5 Figure 5: Derivatives of mutual information and corresponding rotational update. First row: derivative of mutual information with respect to the rotational error ρ ( o ), second row: second derivative and third row: update δρ.

Figure 6 :

 6 Figure 6: Model of the Cycab vehicle (Cycab can be seen on Figure 10).

Figure 7 :

 7 Figure 7: Translation estimation between the current and desired image. (a) Desired image, (b) acquired image with a 4 meter translation, (c) scene depth and (d) derivative of mutual information with respect to the translation along the z axis (in meters) with various fixed scene's depth Z.

Figure 8 :

 8 Figure 8: Simulation experiment. (a) Aerial view of the environment, (b) 2D representation of map with the learned trajectory in green and the resulting path in red (gray rectangles are buildings).

Figure 9 :

 9 Figure 9: Simulation experiment: steering angle (in degrees) in the first turn of the path with respect to the time (in second). The computed steering angle in red shows two exponential decreasement corresponding to two visual servoing tasks, the filtered steering angle is in green.

Figure 10 :

 10 Figure 10: The Cycab vehicle considered in this experiment.

Figure 11 :

 11 Figure 11: Outdoor environment with an approximation of the learned trajectory (the green part is the trajectory executed under a covered parking lot). (a) First aerial view, (b) second aerial view.

Figure 14 :

 14 Figure 14: Mutual information robustness. (a) robutness to illumination variations, (b) illumination variations and snow on the ground, (c): robustness to occlusions. First column: desired image, second column: current image and third column: difference of the current and desired images.

Figure 12 :

 12 Figure 12: Outdoor navigation experiment, begining of the experiment. First row: current image acquired by the vehicle, second row: desired image and third row: difference between the current and desired image. The first top column and the last bottom column show respectively the first images and the last images used in the navigation task.

Figure 13 :

 13 Figure 13: Outdoor navigation experiment, end of the experiment.

Preprint submitted to Robotics and Autonomous Systems December 7, 2012

Acknowledgment

This work was supported by DGA under contribution to student grant.


Y Grandati
email: grandati@yahoo.fr

C Quesne
email: cquesne@ulb.ac.be

Disconjugacy, regularity of multi-indexed rationally-extended potentials, and Laguerre exceptional polynomials

Keywords: Quantum mechanics, supersymmetry, orthogonal polynomials PACS Nos.: 03.65.Fd, 03.65.Ge

The power of the disconjugacy properties of second-order differential equations of Schrödinger type to check the regularity of rationally-extended quantum potentials connected with exceptional orthogonal polynomials is illustrated by re-examining the extensions of the isotonic oscillator (or radial oscillator) potential derived in kth-order supersymmetric quantum mechanics or multistep Darboux-Bäcklund transformation method. The function arising in the potential denominator is proved to be a polynomial with a nonvanishing constant term, whose value is calculated by induction over k. The sign of this term being the same as that of the already known highestdegree term, the potential denominator has the same sign at both extremities of the definition interval, a property that is shared by the seed eigenfunction used in the potential construction. By virtue of disconjugacy, such a property implies the nodeless character of both the eigenfunction and the resulting potential.

I INTRODUCTION

The introduction of the concept of exceptional orthogonal polynomials (EOP) [START_REF] Gómez-Ullate | The Darboux transformation and algebraic deformations of shape-invariant potentials[END_REF][START_REF] Gómez-Ullate | Supersymmetry and algebraic Darboux transformations[END_REF] and the construction of the first X 1 families of Laguerre-and Jacobi-type EOP in the context of Sturm-Liouville theory [START_REF] Gómez-Ullate | An extended class of orthogonal polynomials defined by a Sturm-Liouville problem[END_REF][START_REF] Gómez-Ullate | An extension of Bochner's problem: exceptional invariant subspaces[END_REF] have aroused a lot of interest both in mathematics and in physics. On the one hand, the EOP families are novel systems of complete orthogonal polynomials with respect to some positive-definite measure, generalizing the classical families of Hermite, Laguerre, and Jacobi, but differing from these among others by the presence of gaps in the degrees n of the polynomials appearing in the sequences (e.g., n = 0 for the X 1 families). On the other hand, they were shown [START_REF] Quesne | Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry[END_REF] to provide a keystone for rationally extending some well-known (translationally) shape-invariant quantum potentials [START_REF] Gendenshtein | Derivation of exact spectra of the Schrödinger equation by means of supersymmetry[END_REF][START_REF] Cooper | Supersymmetry and quantum mechanics[END_REF][START_REF] Cariñena | Shape-invariant potentials depending on n parameters transformed by translation[END_REF] in such a way that shape invariance is preserved. In this context, the usefulness of an approach based on supersymmetric quantum mechanics (SUSYQM) or, equivalently, the Darboux transformation was pointed out [START_REF] Bagchi | Isospectrality of conventional and new extended potentials, second-order supersymmetry and role of PT symmetry[END_REF][START_REF] Quesne | Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics[END_REF].

An important step was the description of infinitely many shape-invariant potentials associated with X m EOP families, where the number of gaps m, called codimension (and corresponding here to n = 0, 1, . . . , m -1), may be arbitrarily large [START_REF] Odake | Infinitely many shape invariant potentials and new orthogonal polynomials[END_REF][START_REF] Odake | Another set of infinitely many exceptional (X ℓ ) Laguerre polynomials[END_REF][START_REF] Odake | Infinitely many shape-invariant potentials and cubic identities of the Laguerre and Jacobi polynomials[END_REF][START_REF] Ho | Properties of the exceptional (X ℓ ) Laguerre and Jacobi polynomials[END_REF]. The existence of two families of X m -Laguerre and X m -Jacobi EOP was also reported (thereby extending an observation made for m = 2 in Ref. [START_REF] Quesne | Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics[END_REF]) and it was later on explained through Darboux-Crum transformation [START_REF] Sasaki | Exceptional Laguerre and Jacobi polynomials and the corresponding potentials through Darboux-Crum transformations[END_REF][START_REF] Gómez-Ullate | Exceptional orthogonal polynomials and the Darboux transformation[END_REF][START_REF] Gómez-Ullate | On orthogonal polynomials spanning a non-standard flag[END_REF]. In addition, these families were shown to be obtainable through two alternative approaches, namely the Darboux-Bäcklund transformation [START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF] and the prepotential method [START_REF] Ho | Prepotential approach to solvable rational potentials and exceptional orthogonal polynomials[END_REF].

Counterparts of the classical family of Hermite polynomials have also been known in physical applications since the early 90s [START_REF] Yu | Equidistant spectra of anharmonic oscillators[END_REF][START_REF] Yu | Equidistant spectra of anharmonic oscillators[END_REF][START_REF] Bagrov | Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics[END_REF][START_REF] Samsonov | New features in supersymmetry breakdown in quantum mechanics[END_REF][START_REF] Junker | Conditionally exactly solvable problems and non-linear algebras[END_REF][START_REF] Junker | Conditionally exactly solvable potentials: a supersymmetric construction method[END_REF][START_REF] Cariñena | A quantum exactly solvable nonlinear oscillator related to the isotonic oscillator[END_REF][START_REF] Fellows | Factorization solution of a family of quantum nonlinear oscillators[END_REF], but these X m EOP are of a slightly different nature in the sense that the m gaps in the degree sequence correspond to n = 1, 2, . . . , m, instead of n = 0, 1, . . . , m -1. Since then, similar X m -Laguerre and X m -Jacobi EOP have been constructed and form the so-called third families [START_REF] Quesne | Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF][START_REF] Ho | Prepotential approach to solvable rational potentials and exceptional orthogonal polynomials[END_REF].

A significant advance was then the introduction of multi-indexed X m 1 ,m 2 ,...,m k EOP families and associated rationally-extended potentials through multistep Darboux algebraic transformations [START_REF] Gómez-Ullate | Two-step Darboux transformations and exceptional Laguerre polynomials[END_REF]. The same concept was developed through several approaches, such as the Crum-Adler mechanism [START_REF] Odake | Exactly solvable quantum mechanics and infinite families of multi-indexed orthogonal polynomials[END_REF], higher-order SUSYQM [START_REF] Quesne | Higher-order SUSY, exactly solvable potentials, and exceptional orthogonal polynomials[END_REF][START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF], and multistep Darboux-Bäcklund transformations [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF].

In such a context, there remain two open essential questions. The first one is to know whether the multi-indexed EOP families exhaust all the possibilities of highercodimensional complete orthogonal polynomial systems or, in other words, whether all the higher-codimensional complete orthogonal polynomial systems are generated by the application of successive algebraic Darboux transformations (or one of the equivalent methods quoted above). This was recently proved to be true for codimension two [START_REF] Gómez-Ullate | A conjecture on exceptional orthogonal polynomials[END_REF].

The second problem has to do with the existence of a well-behaved measure defining an EOP system, which is directly connected with the regularity of the associated rationallyextended quantum potentials. In the one-indexed X m case, this problem is usually solved by showing that the classical polynomial occurring in the denominator has no zero in the domain of the variable, which is an easy task since the distribution of zeros of classical polynomials is well known [START_REF] Erdélyi | Higher Transcendental Functions[END_REF][START_REF] Szegö | Orthogonal Polynomials[END_REF]. In the multi-indexed X m 1 ,m 2 ,...,m k case, however, the denominator contains a Wronskian of functions written in terms of classical polynomials, so that the problem becomes much more difficult to solve. While still feasible for k = 2, the study of zeros remains outside scope for most higher k cases, hence the construction of multi-indexed EOP is still rather formal.

A notable progress towards the solution of the regularity question was recently made [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF] by making use of the disconjugacy properties of the Schrödinger equation for eigenvalues below the ground state [START_REF] Hartman | Ordinary Differential Equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF][START_REF] Böcher | Leçons sur les méthodes de Sturm[END_REF]. As a consequence of such properties, the nodeless character of the function present in the denominator can be inferred from the signs it takes on both ends of the domain. Nevertheless when applied to the rationally-extended isotonic oscillator (or radial oscillator) potentials [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF], except in the k = 2 case, the equality of signs was obtained without exhibiting explicit formulas for the denominator in the asymptotic regions.

The purpose of the present paper is to give such explicit expressions in the general multiindexed case. We then extend some results previously obtained in [START_REF] Quesne | Higher-order SUSY, exactly solvable potentials, and exceptional orthogonal polynomials[END_REF][START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF], confirming in particular a conjecture about the structure of the polynomial factor in the denominator. This allows a more direct proof of the regular character of the considered extensions.

The paper is organized as follows. After recalling in Sec. II the disconjugacy properties of the Schrödinger equation, we apply it to the construction of rationally-extended radial oscillator potentials and corresponding X m -Laguerre EOP in first-order SUSYQM. Section III summarizes the corresponding multistep construction in kth-order SUSYQM, where each step may be associated with any one of the first two families of X m -Laguerre EOP.

An improved proof of the regularity of these extensions is then given in Sec. IV. Finally, Section V contains the conclusion.

II RATIONAL EXTENSIONS OF THE ISOTONIC POTENTIAL AND X m -LAGUERRE EOP A Disconjugacy

Let us briefly recall some essential elements about the disconjugacy properties of secondorder differential equations of Schrödinger type,

ψ ′′ (x) + E -V (x) ψ(x) = 0. ( 2.1)
To any solution ψ(x) of Eq. (2.1) we can associate the corresponding Ricatti-Schrödinger function w(x) = -ψ ′ (x)/ψ(x), which is a solution of

-w ′ (x) + w 2 (x) = V (x) -E. (2.2)
Equation (2.1) is said to be disconjugated on I ⊂ R (V (x) being supposed continuous on I) if every solution of this equation has at most one (necessarily simple) zero on I [START_REF] Hartman | Ordinary Differential Equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF].

For a closed or open interval I, the disconjugacy of Eq. (2.1) is equivalent to the existence of solutions that are everywhere non zero on I [START_REF] Hartman | Ordinary Differential Equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF]. In the following we will consider

I =]0, +∞[.

Let us also quote the following result:

Theorem [START_REF] Hartman | Ordinary Differential Equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF] If there exists a continuously differentiable solution on I of the Riccati inequality

-w ′ (x) + w 2 (x) ≤ V (x) -E, (2.3)
then Eq. (2.1) is disconjugated on I.

In the case we consider, the Hamiltonian corresponding to Eq. (2.1) and associated to Dirichlet boundary conditions on I has an infinite bound-state spectrum {E ν , ψ ν } ν∈N .

Consequently, since ψ 0 is nodeless on I, w 0 (x) is a continuously differentiable solution of (2.2) on this domain and condition (2.3) is fulfilled for every value of E not greater than E 0 .

The domain of values ] -∞, E 0 ] of the spectral parameter E will be called the disconjugacy sector of Eq. (2.1). In this sector, except for E = E 0 , ψ(x) does not satisfy the required Dirichlet boundary conditions to be a bound-state wavefunction and diverges at one or both extremities of I. To ensure that this eigenfunction is nodeless, it is sufficient to verify that it has the same sign at both extremities of the definition interval.

B Isotonic oscillator and X m -Laguerre EOP

The isotonic oscillator potential is defined on the positive half-line x > 0 as

V l (x) = 1 4 ω 2 x 2 + l(l + 1) x 2 . (2.4)
The spectrum of the associated Schrödinger equation has an infinite number of bound states, whose (unnormalized) wavefunctions are given by

ψ (l) ν (x) ∝ η l (z)L (α) ν (z), ν ∈ N, (2.5)
where z = ωx 2 /2, α = l + 1/2, and

η l (z) = z (α+1/2)/2 e -z/2 , (2.6)
L (α)
ν (z) being a generalized Laguerre polynomial [START_REF] Erdélyi | Higher Transcendental Functions[END_REF]. The associated energies are

E (l) ν = ω(2ν + α + 1). (2.7)
In first-order SUSYQM, the superpartner Hamiltonians are defined as [START_REF] Cooper | Supersymmetry and quantum mechanics[END_REF]   

H (+) = A † A = -d 2 dx 2 + V (+) (x) -E, H (-) = AA † = -d 2 dx 2 + V (-) (x) -E, (2.8) with    A † = -d dx + W (x), A = d dx + W (x), (2.9)
where the superpotential W (x) satisfies the Riccati equations

∓W ′ (x) + W 2 (x) = V (±) (x) -E.
(2.10) W (x) is obtained from a solution φ(x) of the initial Schrödinger equation

H (+) φ(x) = 0, (2.11)
via

W (x) = -φ ′ (x)/φ(x).
(2.12)

The new potential V (-) (x) = V (+) (x) + 2W ′ (x) is regular as soon as the "seed" solution φ(x) is nodeless on ]0, +∞[. This can be achieved if we choose the factorization energy in the disconjugacy sector of H (+) , that is, when E is not greater than the ground-state energy

E (+) 0 of V (+) (x)
. This is true in particular if we choose as seed solution the ground-state wavefunction ψ

(l) 0 (x) = η l (z) itself. For E < E (+)
0 , as discussed above, the analysis of the sign of φ(x) at 0 + and +∞ is then sufficient to control the regularity of W (x). (2.13)

In the first two cases, 1/φ(x) diverges at one extremity of the positive half-line and, although we have

H (-) 1 φ(x) = 0, (2.14)
it is not a bound-state wavefunction of H (-) . Then H (+) and H (-) turn out to be isospectral.

If we are looking for seed functions which are polynomials (up to some gauge factor), we obtain them from the eigenstates by applying the following discrete symmetries acting on the parameters of the initial potential (2.4) [18]

     Γ I : (ω, α) → (-ω, α), Γ II : (ω, α) → (ω, -α), Γ III : (ω, α) → (-ω, -α).
(2.15) By applying the symmetry Γ i , i = I, II, III, to ψ (l) ν , we obtain a seed eigenfunction φ ν , which satisfies boundary conditions of type i (see (2.13)). In particular, for i = I or II, we obtain Γ I :

ψ (l) m (x) → φ I lm (x) = χ I l (z)L (α) m (-z), Γ II : ψ (l) m (x) → φ II lm (x) = χ II l (z)L (-α) m (z), (2.16)
where

χ I l (z) = z (α+1/2)/2 e z/2 , χ II l (z) = z -(α-1/2)/2 e -z/2 , (2.17)
for the corresponding energies

E I lm = -ω(α + 2m + 1), E II lm = -ω(α -2m -1).
(2.18)

If φ I lm is always in the disconjugacy sector, for φ II lm to reach it we need to assume the condition α > m. Starting from the isotonic potential V (+) (x) = V l ′ (x), the SUSYQM partnership based on seed eigenfunctions φ I l ′ m or φ II l ′ m with respectively l ′ = l -1 and l ′ = l + 1 gives then rational extensions of the form [31]

V (-) (x) = V l (x) + V l,rat (x) + C, (2.19)
where

V l,rat (x) = -2ω    ġ(α) m g (α) m + 2z   g(α) m g (α) m - ġα) m g (α) m 2      (2.20)
(the dot denoting a derivative with respect to z), with in the first case

g (α) m (z) = L (α-1) m (-z), C = -ω (2.21)
and in the second case

g (α) m (z) = L (-α-1) m (z), C = ω. (2.22)
Due to disconjugacy properties, the regularity of these extensions on the positive halfline is ensured for every m = 1, 2, . . . (and α large enough for type II).

The bound-state wavefunctions of the superpartner potential V (-) (x) are obtained by applying the operator A (with

W (x) = -φ i′ l ′ m (x)/φ i l ′ m (x), i = I, II) on those of V (+) (x) (ψ (+) ν (x) ∝ η l ′ (z)L (α ′ ) ν (z) with α ′ = l ′ + 1/2). We can write ψ (-) ν (x) ∝ η l (z) g (α) m (z) y (α) m+ν (z), ν = 0, 1, 2, . . . , (2.23)
where y

(α) m+ν (z) is a polynomial of degree n = m + ν, y (α) n (z) = L I α,m,n (z) if i = I, L II α,m,n (z) if i = II, (2.24)
satisfying the second-order differential equation

z d 2 dz 2 + α + 1 -z -2z ġ(α) m g (α) m d dz + (z -α) ġ(α) m g (α) m + z g(α) m g (α) m y (α) m+ν (z) = -νy (α) m+ν (z). (2.25)
The set of polynomials y (α) m+ν (z) (ν = 0, 1, 2, . . . ), defined on the positive half-line, is orthogonal and complete with respect to the positive-definite measure of density

z α e -z / g (α) m (z) 2 .
(2.26)

III LAGUERRE EOP AND REDUCIBLE k-TH OR-DER SUSYQM

In kth-order SUSYQM, the first-order differential operators A and A † are replaced by kthorder ones A and A † [START_REF] Bagrov | Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics[END_REF][START_REF] Samsonov | New features in supersymmetry breakdown in quantum mechanics[END_REF][START_REF] Andrianov | Higher-derivative supersymmetry and the Witten index[END_REF][START_REF] Andrianov | Second order derivative supersymmetry, q deformations and the scattering problem[END_REF][START_REF] Fernández | Higher-order supersymmetric quantum mechanics[END_REF]. The correspondence between the initial and final

Hamiltonians H (1) and H (2) is determined by the intertwining relations

AH (1) = H (2) A, A † H (2) = H (1) A † . (3.27)
In the reducible case, the operators A and A † can be factorized into products of firstorder differential operators as

A = A (k) A (k-1) • • • A (1) , (3.28)
where

A (i) = d dx + W (i) (x), W (i) (x) = -log φ (i) (x) ′ , (3.29)
the seed eigenfunctions being given by Crum's formula [START_REF] Crum | Associated Sturm-Liouville systems[END_REF]

φ (i) (x) = φ (i) 1,2,...,i (x) = W(φ 1 , φ 2 , . . . , φ i | x) W(φ 1 , φ 2 , . . . , φ i-1 | x) , (3.30)
with W(φ 1 , φ 2 , . . . , φ i |x) denoting the Wronskian of φ 1 (x), φ 2 (x), . . . , φ i (x). Here we have introduced lower indices 1, 2, . . . , i to specify which seed functions have been used. The partner potentials are then linked by the second Crum's formula

V (2) (x) = V (1) (x) -2 d 2 dx 2 log W(φ 1 , φ 2 , . . . , φ k |x). (3.31)
As in Sec. II, we start from the isotonic potential V (1) (x) = V l ′ (x) to get an extended potential of the form

V (2) (x) = V l (x) + V l,rat (x) + C. (3.32)
The final potential is of type I q II k-q if the chosen set of seed eigenfunctions includes q ≤ k seed functions of type I and k -q seed functions of type II. The order of the seed functions being irrelevant, we can assume that

φ i (x) = φ I l ′ m i (x) = χ I l ′ (z)L (α ′ ) m i (-z), i = 1, . . . , q, φ II l ′ m i (x) = χ II l ′ (z)L (-α ′ ) m i (z), i = q + 1, . . . , k, (3.33)
with

l ′ = l + k -2q (hence α ′ = α + k -2q), 0 < m 1 < • • • < m q , and 0 < m q+1 < • • • < m k ,
these functions being nodeless for α ′ > sup i=q+1,...,k (m i ). In [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF], it has been shown that the Wronskian appearing in Crum's formula (3.31) can be written as

W(φ 1 , φ 2 , . . . , φ k | x) = (ωx) k(k-1)/2 z -q(k-q) χ I l ′ q χ II l ′ k-q g (α) µ (z), (3.34)
where

g (α) µ (z) = z -(k-q)(k-q-1) det Γ(α) µ , (3.35)
with Γ(α)

µ ij =                    L (α ′ +i-1) m j -i+1 (-z), 1 ≤ i ≤ q + 1, 1 ≤ j ≤ q, (m j + 1) i-1 z k-i L (-α ′ -i+1) m j +i-1 (z), 1 ≤ i ≤ q + 1, q + 1 ≤ j ≤ k, L (α ′ +i-1) m j -q (-z), q + 2 ≤ i ≤ k, 1 ≤ j ≤ q, (m j + 1) q (m j -α ′ -i + q + 2) i-q-1 × z k-i L (-α ′ -i+1) m j +q (z), q + 2 ≤ i ≤ k, q + 1 ≤ j ≤ k, (3.36)
and (a) n = a(a + 1) • • • (a + n -1) denoting the usual Pochhammer symbol [START_REF] Erdélyi | Higher Transcendental Functions[END_REF].

Observe that on taking (3.31) and (3.34) into account, we can express V l,rat (x) and C of Eq. (3.32) as

V l,rat (x) = -2ω    ġ(α) µ g (α) µ + 2z   g(α) µ g (α) µ - ġα) µ g (α) µ 2      , C = (k -2q)ω. (3.37)
In [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF], it has been proved that

g (α) µ (z) = C (α) µ z µ + lower-order terms, (3.38) with        µ = k i=1 m i -q(q -1)/2 -(k -q)(k -q -1)/2 + q(k -q), C (α) µ = (-1) σ ∆(m 1 , . . . , m q )∆(m q+1 , . . . , m k )/(m 1 ! . . . m k !), σ = k i=q+1 m i + q(k -q), (3.39)
∆(n 1 , . . . , n l ) being a Vandermonde determinant of order l,

∆(n 1 , . . . , n l ) = 1≤i<j≤l (n j -n i ). (3.40)
As shown in [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF], the Wronskian (3.34) has no node on the positive half-line, which ensures the regularity of the rationally-extended potential V (2) (x) on this domain. The bound-state wavefunctions of the latter are then given by [31]

ψ (2) ν (x) = η l (z) g (α) µ (z) y (α) µ+ν (z), ν = 0, 1, 2, . . . , (3.41)
where y

(α)
µ+ν (z) is a polynomial of degree n = µ + ν satisfying the second-order differential equation

z d 2 dz 2 + α + 1 -z -2z ġ(α) µ g (α) µ d dz + (z -α) ġ(α) µ g (α) µ + z g(α) µ g (α) µ y (α) µ+ν (z) = -νy (α) µ+ν (z). (3.42)
The y (α)

µ+ν (ν = 0, 1, 2, . . . ) polynomials constitute a complete orthogonal family on the positive half-line with respect to the positive-definite measure of density z α e -z / g (α) µ (z) 2 .

(3.43)

IV POLYNOMIAL CHARACTER OF g (α) µ AND BE-HAVIOUR NEAR THE ORIGIN

In [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF] it has been conjectured that g

(α)
µ (z) is a polynomial, which can be directly checked in the k = 2 case, but is not obvious for higher k values whenever k -q ≥ 2 (see Eq. (3.35)).

In the following we will give a complete proof of this assertion. In [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF] the regularity of the extended potential, which is equivalent to the absence of node of g

(α)
µ , has been established without giving the exact highest and lowest degrees coefficients of this polynomial. We will also determine the value of the lowest-degree coefficient and then verify explicitly the nodeless character of the Crum Wronskian.

First note that from Eq. (3.34) we can write for k -q ≥ 2

W(φ 1 , φ 2 , . . . , φ k | x) = (ωx) k(k-1)/2 z -q(k-q) (χ I l ′ ) q (χ II l ′ ) k-q g (α) µ (z), W(φ 1 , φ 2 , . . . , φ k-1 | x) = (ωx) (k-1)(k-2)/2 z -q(k-q-1) (χ I l ′ ) q (χ II l ′ ) k-q-1 g (α+1) µ ′ (z), (4.44)
where

l ′ = l + k -2q (hence α ′ = α + k -2q) and        µ = k i=1 m i -q(q -1)/2 -(k -q)(k -q -1)/2 + q(k -q), µ ′ = k-1 i=1 m i -q(q -1)/2 -(k -q -1)(k -q -2)/2 + q(k -q -1) = µ -m k + k -2q -1. (4.45)
Inserting these expressions in Eq. (3.30), we arrive at

φ (k) 1,2,...,k (x) = (ωx) k-1 z -q χ II l ′ (z) g (α) µ (z) g (α+1) µ ′ (z) = (2ω) (k-1)/2 z -(2α ′ -2k+4q+1)/4 e -z/2 g (α) µ (z) g (α+1) µ ′ (z) . (4.46)
On the other hand, we can obtain a different expression for φ

(k)
1,2,...,k by using some properties of Wronskians. Indeed, from Sylvester's theorem [START_REF] Muir | A treatise on the theory of determinants[END_REF], we can write

W(φ 1 , . . . , φ k-1 , φ k | x) = W W(φ 1 , . . . , φ k-2 , φ k-1 | x), W(φ 1 , . . . , φ k-2 , φ k | x) | x W(φ 1 , . . . , φ k-2 | x) (4.47)
and using

W(f φ 1 , . . . , f φ k | x) = f k W(φ 1 , . . . , φ k | x), (4.48)
we have

W(φ 1 , . . . , φ k-1 , φ k | x) = W(φ 1 , . . . , φ k-2 | x) × W W(φ 1 , . . . , φ k-2 , φ k-1 | x) W(φ 1 , . . . , φ k-2 | x) , W(φ 1 , . . . , φ k-2 , φ k | x) W(φ 1 , . . . , φ k-2 | x) x , (4.49)
that is

φ (k) 1,2,...,k (x) = W φ (k-1) 1,...,k-2,k-1 , φ (k-1) 1,...,k-2,k x φ (k-1) 1,...,k-2,k-1 (x)
.

(4.50)

From the Wronskian theorem [START_REF] Messiah | Mécanique quantique[END_REF], it follows that the Wronskian in the numerator of this expression satisfies the property

W φ (k-1) 1,...,k-2,k-1 , φ (k-1) 1,...,k-2,k ξ +∞ x = (E k-1 -E k ) +∞ x φ (k-1) 1,...,k-2,k-1 (ξ)φ (k-1) 1,...,k-2,k (ξ)dξ, (4.51)
because φ

(k-1) 1,...,k-2,k-1 and φ (k-1) 1,...,k-2,k correspond to the energies k-1 i=1 E i and k-2 i=1 E i + E k , respectively.
As a consequence of the exponential factor in Eq. (4.46), both φ

W φ (k-1) 1,...,k-2,k-1 , φ (k-1) 1,...,k-2,k x = -(2ω) k-5/2 (E k-1 -E k ) × +∞ z ζ k-2q-α ′ -2 e -ζ g (α+1) µ ′ (ζ)g (α+1) μ′ (ζ) g (α+2) µ ′′ (ζ) 2 dζ, (4.53)
where

μ′ = µ -m k-1 + k -2q -1, µ ′′ = µ -m k-1 -m k + 2(k -2q) -3. ( 4
φ (k) 1,2,...,k (x) = -(2ω) (k-3)/2 (E k-1 -E k ) × +∞ z ζ k-2q-α ′ -2 e -ζ g (α+1) µ ′ (ζ)g (α+1) μ′ (ζ)/ g (α+2) µ ′′ (ζ) 2 dζ × z (2k-4q-2α ′ -3)/4 e -z/2 g (α+1) µ ′ (z)/g (α+2) µ ′′ (z) -1 , (4.55)
that is (see Eq. (3.40))

g (α) µ (0) = (α ′ + 1) m 1 • • • (α ′ + k) m k -k+1 m 1 ! • • • m k ! ∆(m 1 , . . . , m k ) = (α -k + 1) m 1 • • • (α) m k -k+1 m 1 ! • • • m k ! ∆(m 1 , . . . , m k ). (4.61)
Consider now the first mixed case I k-1 II (k = q + 1), where α ′ = α -k + 2 [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF]. We then get

g (α) µ (z) = L (α ′ ) m 1 (-z) . . . L (α ′ ) m k-1 (-z) z k-1 L (-α ′ ) m k (z) . . . . . . . . . . . . L (α ′ +k-2) m 1 -k+2 (-z) . . . L (α ′ +k-2) m k-1 -k+2 (-z) (m k + 1) k-2 zL (-α ′ -k+2) m k +k-2 (z) L (α ′ +k-1) m 1 -k+1 (-z) . . . L (α ′ +k-1) m k-1 -k+1 (-z) (m k + 1) k-1 L (-α ′ -k+1) m k +k-1 (z) , (4.62)
which leads to

g (α) µ (0) = (m k + 1) k-1 L (-α ′ -k+1) m k +k-1 (0) L (α ′ ) m 1 (0) . . . L (α ′ ) m k-1 (0) . . . . . . . . . L (α ′ +k-2) m 1 -k+2 (0) . . . L (α ′ +k-2) m k-1 -k+2 (0) . (4.63)
The right-hand member determinant is exactly the constant term of the g polynomial in the pure case I k-1 and we deduce from the previous result that

g (α) µ (0) = (α ′ + 1) m 1 • • • (α ′ + k -1) m k-1 -k+2 m 1 ! • • • m k-1 ! ∆(m 1 , . . . , m k-1 ) × (m k + 1) k-1 (-1) m k +k-1 (α ′ -m k ) m k +k-1 (m k + k -1)! , (4.64)
or

g (α) µ (0) = (-1) m k +k-1 (α ′ + 1) m 1 • • • (α ′ + k -1) m k-1 -k+2 (α ′ -m k ) m k +k-1 m 1 ! • • • m k ! × ∆(m 1 , . . . , m k-1 ) = (-1) m k +k-1 (α -k + 3) m 1 • • • (α + 1) m k-1 -k+2 (α -k + 2 -m k ) m k +k-1 m 1 ! • • • m k ! × ∆(m 1 , . . . , m k-1 ). ( 4

.65)

For the general mixed case I q II k-q , we now make the following hypothesis

g (α) µ (0) = (-1) σ ∆(m 1 , . . . , m q )∆(m q+1 , . . . , m k ) m 1 ! • • • m k ! × (α + k -2q + 1) m 1 • • • (α + k -q) mq-q+1 × (α + k -2q -m q+1 ) m q+1 +q • • • (α + k -2q -m k ) m k +2q-k+1 , (4.66)
which combined with the previous expressions yields g (α) µ

(0) = (-1) σ ′ +σ ′ -σ ′′ (m k -m k-1 ) (α + 1)m 1 ! • • • m k ! × ∆(m 1 , . . . , m q )∆(m q+1 , . . . , m k-1 )∆(m q+1 , . . . , m k-2 , m k ) ∆(m q+1 , . . . , m k-2 ) × (α + k -2q + 1) m 1 • • • (α + k -q) mq-q+1 × (α + k -2q -m q+1 ) m q+1 +q • • • (α + k -2q -m k-2 ) m k-2 +2q-k+3 × (α + k -2q -m k-1 ) m k-1 +2q-k+2 (α + k -2q -m k ) m k +2q-k+2 . (4.72)
Here we note that

σ ′ + σ′ -σ ′′ = k-2 l=q+1 m l + m k-1 + m k + q(k -q) = σ, (4.73)
(α + k -2q -m k ) m k +2q-k+2 α + 1 = (α + k -2q -m k ) m k +2q-k+1 , (4.74)
and

(m k -m k-1 ) ∆(m q+1 , . . . , m k-1 )∆(m q+1 , . . . , m k-2 , m k ) ∆(m q+1 , . . . , m k-2 ) = (m k -m k-1 ) k-2 i=q+1 k-1 j=i+1 (m j -m i ) k-2 i=q+1 (m k -m i ) = k-1 i=q+1 k j=i+1 (m j -m i ) = ∆(m q+1 , . . . , m k ). ( 4

.75)

The final result for g

µ (0) is therefore given by Eq. (4.66), showing that the conjecture is still verified for k, which ends its proof.

Let us observe that in terms of α ′ , Equation (4.66) can be rewritten as We conclude that the sign of g (α)

g (α) µ (0) = (-1) σ ∆(m 1 , . . . , m q )∆(m q+1 , . . . , m k ) m 1 ! • • • m k ! × (α ′ + 1) m 1 • • • (α ′ + q) mq-q+1 × (α ′ -m q+1 ) m q+1 +q • • • (α ′ -m k ) m k +2q-k+1 . ( 4
µ (0) is given by (-1) σ , hence is identical with the one of C

V CONCLUSION

In the present work, we have established some further properties of the function g

(α)
µ (z) appearing in the denominator of rationally-extended isotonic oscillator (or radial oscillator) potentials constructed in kth-order SUSYQM [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF] or using the multistep Darboux-Bäcklund transformation method [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF]. We have first confirmed that g (α) µ (z) cannot have any pole at the origin and is a polynomial with non-zero constant term, a property that was conjectured in Ref. [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF]. We have then determined the value of this constant term by induction over k and shown that it has the same sign as the coefficient of the previously determined highest-degree term z µ [START_REF] Quesne | Rationally-extended radial oscillators and Laguerre exceptional orthogonal polynomials in kth-order SUSYQM[END_REF].

These results have been combined with the disconjugacy theorem applied to Schrödinger equation [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF] to prove the absence of zeros of g (α) µ (z) on the positive half-line and, therefore, the regularity of the corresponding rationally-extended potential. This has provided a demonstration of the latter property, which is more explicit than that already given in Ref. [START_REF] Grandati | Multistep DBT and regular rational extensions of the isotonic oscillator[END_REF].

Our work has illustrated the power of the disconjugacy properties of second-order differential equations of Schrödinger type to check the regularity of rational extensions of well-known quantum potentials. Other examples of application might be found in cases where no direct information on the denominator zeros is available. This will be the object of further investigations.

  Suppose that we are in this case. Then we have three possibilities (up to a global change of sign), namely      I : φ(0 + ) = 0 + , φ(+∞) = +∞, II : φ(0 + ) = +∞, φ(+∞) = 0 + , III : φ(0 + ) = +∞, φ(+∞) = +∞.

  ..,k-2,k-1 (ξ) and φ (k-1) 1,...,k-2,k (ξ) vanish for ξ → ∞, which ensures the vanishing of their Wronskian in the same limit. Hence Equation (4.51) reduces toW φ (k-1) 1,...,k-2,k-1 , φ (k-1) 1,...,k-2,k x = -(E k-1 -E k ) ..,k-2,k (ξ)dξ, (4.52) or, on taking Eq. (4.46) into account and defining ζ = ωξ 2 /2,

  .54) Inserting Eqs. (4.46) and (4.53) in the right-hand side of Eq. (4.50), we get the searched for second expression of φ (k) 1,2,...,k (x),

  .76) Since the angular momentum l ′ in the starting potential V (1) (x) is necessarily nonnegative, the inequality α ′ > -1/2 is fulfilled. Furthermore, we have assumed 0< m 1 < • • • < m q and 0 < m q+1 < • • • < m k < α ′ .As a consequence, the two Vandermonde determinants in Eq. (4.76) are positive and(α ′ + 1) m 1 • • • (α ′ + q) mq-q+1 (α ′ -m q+1 ) m q+1 +q • • • (α ′ -m k ) m k +2q-k+1 > 0.(4.77)

  .(3.39)). The combination of this result with the disconjugacy theorem of Sec. II ensures the regularity of the extended potential V (2) (x).

which can also be rewritten as a recursion relation for the g (α) µ (z) functions, g (α) µ (z) =

On the right-hand side of this relation, g

µ (z) can be easily obtained by induction. This is already verified for k = 1 and k = 2. Assuming that for k ′ < k, g (α) µ is a polynomial with a non-zero constant term, it turns out that in the neighbourhood of the origin, the right-hand side of Eq. ( 4.56) has a well-defined value given by

where we have used

From Eq. ( 4.57) we immediately deduce that g

µ (z) cannot have a pole at the origin and is then a polynomial with non-zero constant term.

To completely describe the behaviour of g (α) µ (z) at the origin, we need to determine the exact expression of this constant term. For such a purpose, we will proceed by induction again. Consider first the pure case I k (k = q), where α ′ = α -k. From (3.35) and (3.36), we have

.58)

On using [START_REF] Erdélyi | Higher Transcendental Functions[END_REF]]

we easily get

with σ defined in Eq. (3.39), and prove it by induction over k.

This conjecture is verified in the k = 1 and k = 2 cases. Since it is also valid for the two previous examples (k -q = 0 or 1), it only remains to prove it whenever k -q ≥ 2. In such cases, the last two seed functions φ k-1 and φ k , used in kth-order SUSYQM, are necessarily of type II. Suppose that Eq. (4.66) is verified for k ′ < k. We may therefore write

σ ′ = k-1 l=q+1 m l + q(k -q -1), σ ′′ = k-2 l=q+1 m l + q(k -q -2), σ′ = k-2 l=q+1 m l + m k + q(k -q -1).