Cédric Demonceaux
email: cedric.demonceaux@u-picardie.fr

Pascal Vasseur
email: pascal.vasseur@u-picardie.fr

Claude Pégard
email: claude.pegard@u-picardie.fr

Robust Attitude Estimation with Catadioptric Vision

Attitude (roll and pitch) is an essential data for the navigation of a UAV. Rather than using inertial sensors, we propose a catadioptric vision system allowing a fast, robust and accurate estimation of these angles. We show that the optimization of a sky/ground partitioning criterion associated with the specific geometric characteristics of the catadioptric sensor provides very interesting results. Experimental results obtained on real sequences are presented and compared with inertial sensors measures.

I. INTRODUCTION

Autonomous navigation of a UAV requires the knowledge of its attitude (roll and pitch). Using computer vision simultaneously for attitude estimation and other applications can be very interesting in order to bound the number of embedded sensors. In this paper, we propose a catadioptric vision system for the attitude estimation. We demonstrate how the specific geometric characteristics of the sensor permit a formulation as an optimization problem which provides fast, robust and accurate results. Attitude computation methods based on vision generally consist in detecting the horizon and in estimating the angles in comparison with it [START_REF] Ettinger | Visionguided flight stability and control for micro air vehicles[END_REF] [8] [START_REF] Woo | Robust horizon and peak extraction for vision-based navigation[END_REF]. In [START_REF] Ettinger | Visionguided flight stability and control for micro air vehicles[END_REF] and [START_REF] Todorovic | Sky/ground modeling for autonomous mav flight[END_REF], the authors use a classical projective camera in order to detect the horizon. They have first presented a sky-ground modelization based only on colors [START_REF] Ettinger | Visionguided flight stability and control for micro air vehicles[END_REF]. The horizon is considered as a straight line and the method consists in searching in the image the straight line which maximizes the sky-ground separation criterion. This model has been improved in [START_REF] Todorovic | Sky/ground modeling for autonomous mav flight[END_REF] by adding texture information. In these approaches, the roll angle is given by the inverse tangent of the slope of the horizon line while the pitch angle is approximated by the percentage of sky in the image. In order to obtain a more accurate detection of the horizon line and of the peaks in the image, in [START_REF] Woo | Robust horizon and peak extraction for vision-based navigation[END_REF] the authors propose to use an infrared camera. Indeed, the difference in temperature between sky and ground is so large that the intensities of both regions are very different. In spite of the interesting results obtained with the previous approaches, the use of a single perspective camera generates several drawbacks. First, only a partial view of the environment is available and important occlusions in the horizon can have a serious influence on the final result. Second, the horizon is visible only in a particular interval of roll and pitch values. If the UAV gets out of this interval, the final image is exclusively made of sky or earth and the horizon can not be detected. Third, it is only possible to compute the roll angle while the pitch is only approximated thanks to an hypothesis on the altitude of the UAV [START_REF] Ettinger | Visionguided flight stability and control for micro air vehicles[END_REF]. In [START_REF] Wang | Complete pose determination for low altitude unmanned aerial vehicle using stereo vision[END_REF], the authors use a stereovision system in order to improve the computation of the attitude by determining the complete pose of the UAV. However, this system relies on the capture of ground targets/landmarks in both images which limits the environment in which the UAV can move. Rather than using perspective cameras, we propose in this paper a system based on a central catadioptric camera. Catadioptric vision consists in associating a convex mirror with a projective camera whose optical axis is aligned with the axis of the mirror [START_REF]Panoramic Vision: Sensors, Theory, Applications[END_REF]. The main advantage of these sensors is the acquisition of an omnidirectional image with a single shot. This kind of sensor has already been used for the navigation of UAVs in corridor [START_REF] Hrabar | Omnidirectional vision for an autonomous helicopter[END_REF]. They present the following advantages for the attitude computation : first, the surrounding of the UAV is completely captured and occlusions will then have a lower impact on the estimation of the final results. Second, whatever the attitude of the UAV, the horizon is always present in the image, even partially, and the angles can always be computed. Third, we are able to compute the roll angle but also the pitch angle without any prior hypothesis, contrary to the applications including a perspective camera. In [START_REF] Demonceaux | Omnidirectional vision on uav for attitude computation[END_REF], we proposed a first approach of attitude computation based on catadioptric vision. In this method, the algorithm consists first in segmenting the image in two classes respectively the earth and the sky thanks to Markov Random Fields segmentation based on color cues. The second step extracts the horizon pixels in the image which are projected onto the equivalent sphere. Then a M-robust estimation method is used in order to compute the best plane which passes through the horizon pixels and the attitude angles can be estimated. This method demonstrates very interesting results but do not use sufficiently the geometric characteristics of the catadioptric vision. Moreover, the MRFs segmentation and the M-robust estimator are particularly time consumers and do not permit a real time implementation. In this paper, we present a new method based on catadioptric vision which presents better accuracy, robustness and computation time. We use the particular geometric characteristics of the sensor for a formulation of the process as an optimization problem which is solved on the sphere of equivalence in order to compute directly the attitude angles.

In Section 2, we present the model of image formation and the catadioptric projection of the horizon. Section 3 is devoted to the formulation of the optimization problem for the horizon detection. In Section 4, we propose our algorithm for the solution of the optimization problem and the roll and pitch estimation. Experimental results are presented in Section 5. Finally, we conclude in Section 6 with a summary and perspectives.

II. CENTRAL CATADIOPTRIC PROJECTION OF THE HORIZON

Baker and Nayar classified catadioptric sensors into two categories depending on the number of viewpoints [START_REF] Baker | A theory of catadioptric image formation[END_REF]. Sensors with a single viewpoint, named central catadioptric sensors, permit a geometrically correct reconstruction of the perspective image from the original catadioptric image. Geyer and Daniilidis [START_REF] Geyer | Catadioptric projective geometry[END_REF] have demonstrated the equivalence for the single viewpoint category with a two-step projection via a unitary sphere centered on the focus of the mirror (the single viewpoint) (Fig. 1). This two-step projection consists first in projecting real 3D point P w to point P s from the center of the sphere O c . The second projection projects point P s to point P i in the image plane from point O p . In order to apply the equivalence, it is necessary to know the intrinsic parameters of the camera and two additional parameters namely ξ and ϕ which are respectively equal to distances |O c O p | and |O c O i | (Fig. 1). Parameters ξ and ϕ define the shape of the mirror (see [START_REF] Barreto | Issues on the geometry of central catadioptric image formation[END_REF] for further details on their signification) and can be estimated by calibration [START_REF] Ying | Catadioptric camera calibration using geometric invariants[END_REF]. This equivalence is particularly interesting because it permits to simplify the horizon detection and the attitude computation. As demonstrated in [START_REF] Ying | Catadioptric camera calibration using geometric invariants[END_REF], a 3D sphere projects on the equivalent sphere in a circle, and then on the catadioptric image plane in an ellipse. Consequently, the attitude computation consists in looking for an ellipse in the omnidirectional image or a circle on the equivalent sphere which corresponds to the horizon (Fig. 2). The geometrical properties of the equivalent sphere allow to deduce the roll and pitch angles. Indeed, if we consider Figure 3, we can note that the normal of the projected horizon on the sphere, which is also confounded with the line passing through the center of the sphere of equivalence and through the center of the earth represents in fact the attitude of the UAV depending on the position of the optical axis. Then, the computation of this normal permits to deduce the roll (ρ) and pitch (ψ) angles.

III. HORIZON DETECTION AS AN OPTIMIZATION PROBLEM IN EQUIVALENCE SPHERE

From the previous considerations, the horizon line in a central catadioptric image is equivalent to a circle on the equivalent sphere and then to a plane in 3D space which partitions the sphere into two regions (sky/ground). Let C a circle on the sphere S and P the plane which contains C (C = P∩S). We can show that C can be determined by a single point

(x c , y c , z c ) = (0, 0, 0) ∈ B(0, 1) = {(x, y, z)|x 2 + y 2 + z 2 ≤
1} and that the equation of P is:

x c (x -x c ) + y c (y -y c ) + z c (z -z c ) = 0. Moreover, the radius of C is R = 1 -(x 2 c + y 2 c + z 2 c
). Reciprocally, any point (x c , y c , z c ) = (0, 0, 0) ∈ B(0, 1) defines a single circle C. This property demonstrates that searching a circle (except the great circles) on the sphere is equivalent to searching a point = (0, 0, 0) inside the sphere.

Point (x c , y c , z c ) ∈ B(0, 1) partitions sphere S into two regions R s = {(x, y, z) ∈ S|x c (x -x c ) + y c (y -y c ) + z c (z - z c ) >= 0} and R g = {(x, y, z) ∈ S|x c (x -x c ) + y c (y - y c ) + z c (z -z c ) < 0}.
Searching this point which corresponds to the horizon in the image consists in classifying the points on the sphere {(x, y, z) ∈ S} into two classes. These classes represent respectively the sky and the ground and the color of the points defined in RGB space is used as measure of appearance.

Let R s and R g , these two classes characterized by their means

m s = (m R s , m G s , m B s ) and m g = (m R g , m G g , m B g
) and by their covariance matrices Γ s and Γ g . If we consider the Mahalanobis distance given by

d(R s , R g ) = (m s -m g ) T (Γ s + Γ g ) -1 (m s -m g )
the best partitioning into R s and R g is then defined by arg max

Rs,Rg d(R s , R g ).
In our case, the problem can be formulated as follows: we are looking for point (x c , y c , z c ) which gives the maximum of, B(0, 1)\(0, 0, 0) →

(x c , y c , z c ) → d(R s , R g ). (1)
Then, we can compute pitch (ψ) and roll (ρ) angles as follows :

ψ = x c | x c | × arccos |z c | z 2 c + x 2 c (2) ρ = -y c | y c | × arccos |z c | z 2 c + y 2 c ( 3
)
Searching this optimum of (1) is not trivial. Indeed, the optimization criterion is not convex and is particularly hard to control. Moreover, since an exhaustive research of the maximum by a very fine sampling of the sphere can not be performed because of the computation time, we propose an algorithm based on a multiscale sampling of the sphere.

IV. OPTIMIZATION ALGORITHM FOR HORIZON DETECTION IN EQUIVALENCE SPHERE

In order to solve (1), the search space B(0, 1) is discretized in 2 3N points :

{X ijk c = (-1 + i 2 2 N , -1 + j 2 2 N , -1 + k 2 2 N )|i, j, k = 0 . . . 2 N }.
The optimization algorithm consists first in evaluating the functional (1) for these 2 3N points. Next, we perform a research in the neighborhood of point X ijk c which maximizes the functional by refining successively the discretization until the desired precision has been reached. It is worth noting that the discretization of B(0, 1) is only necessary for the first image of the sequence in order to obtain the initialization of the point (x c , y c , z c ). Indeed, it appears reasonable to hypothesize that the searched point presents a very low variation between two consecutive images. Consequently, the circle search at time t+1 is performed in the neighborhood of the point computed at time t. The algorithm is summarized in Table I.

[ρ(t), ψ(t)] = Opt(I(t), , N )

% I(t) =image at time t % =precision term % N = initial cutting n ← N ρ(t) ← 0 ψ(t) ← 0 if t = 0 Xc(t) = (xc, yc, zc) ← (0, 0, 0) for i = -1 : 2 -n : 1 for j = -1 : 2 -n : 1 for k = -1 : 2 -n : 1 E(i, j, k) ← d(Rs, Rg) endfor endfor endfor Xc(t) = ( xc, yc, zc) ← arg max i,j,k E(i, j, k) estimation of ρ(t) and ψ(t) elseif Xc(t) ← Xc(t -1) endif While max(||ψ(t) -ψc(t)||, ||ρ(t) -ρ(t)||) ≥ Xc(t) ← Xc(t) ρ(t) ← ρ(t) ψ(t) ← ψ(t) Xc ← Xc for i = xc -2 -n : 2 -(n+1) : xc + 2 -n for j = yc -2 -n : 2 -(n+1) : yc + 2 -n for k = zc -2 -n : 2 -(n+1) : zc + 2 -n E(i, j, k) ← d(Rs, Rg) endfor endfor endfor Xc(t) ← arg max i,j,k E(i, j, k) estimation of ρ(t) and ψ(t) n ← n + 1 Loop TABLE I OPTIMIZATION ALGORITHM
controls the accuracy of the result, it is experimentally fixed to 0.5 • . Parameter N also has an important impact since it defines the size of the first discretization of the sphere. More this value is high and more the first discretization is fine, and consequently the risk of convergence in a local minimum is high. In contrary, more this parameter is low and more the algorithm is fast. It is then necessary to find a good compromise between computation time and convergence of the algorithm. Experimental results show that a value equal to 5 for N corresponds to this compromise.

V. EXPERIMENTAL RESULTS

In order to demonstrate the validity of our approach, we propose different tests which permit to evaluate the accuracy, the sensitivity and the robustness. We also present experimental results on a real video acquired from an airplane. In all experiments, we set the parameter N to 5 and to 0.5 • .

A. Sensor Calibration and Projection on the Sphere

In considering the model of the equivalent sphere, for a 3D point P w with coordinates (x w , y w , z w ), we obtain in the camera frame :

       x = (ξ+ϕ)xw ξ √ x 2 w +y 2 w +z 2 w -zw y = (ξ+ϕ)yw ξ √ x 2 w +y 2 w +z 2 w -zw z = -ϕ (4)
and in the image frame :

u v = α u 0 u 0 0 α v v 0   x y 1   , (5)
In order to perform the calibration with this model, we place the catadioptric system in a cube which forms a grid of points on each side. In this way, the points of the pattern are distributed regularly over the whole catadioptric image. The pattern contains 144 points and we estimate six extrinsic parameters, mirror parameters and intrinsic parameters of the camera. The estimation is performed by the minimization of the quadratic error between the selected points and those computed by the model. Thanks to this calibration, we are then able to project the catadioptric image on the equivalent sphere as shown in figure 4.

B. Accuracy Evaluation

In order to evaluate the accuracy of our system, we propose to compare our results with the measures given by an inertial sensor. In this way, 49 images were captured for which roll and pitch angles are included in the interval -30 • to 30 • with a step equal to 10 • (fig 5). Figure 6 shows the roll and pitch errors over the 49 images. The mean error for the roll is equal to 1, 3 • and to 2, 1 • for the pitch. It is worth noting that the images were taken from a very low altitude (on the roof of the university) in a urban environment. These conditions are typically the main reasons of the pitch error equal to 4 degrees when the catadioptric sensor is very inclined. Indeed, buildings in the neighborhood occult strongly the horizon line et the algorithm tends to find the best compromise between the sky in R g and the buildings in R s during the image partitioning as shown in figure 7. However, this problem is negligible in real sequences when the UAV has a higher altitude (fig 10).

C. Robustness Evaluation

In order to verify the robustness of the method, we apply our algorithm to an image degraded with a white noise. This image was taken with angles equal to -20 • for roll and -30 • for pitch. Results for attitude estimation are shown in figures 8 and 9 according to the standard deviation sigma (σ = 0 to 50) of the added white noise. This results show a very good stability of the angle estimation even when the white noise is very important in the image. For a white noise with a standard deviation included between 0 and 30, the roll angle has a variation less than 0. noting that the results are coherent in the sense that the angles curves are continuous (fig 11). The method has been implemented with Matlab and the computation time is about five seconds per image which can be easily improved in order to obtain a real time application. The result of this sequence is available at the following website address http://www.crea.upicardie.fr/ vasseur.

VI. CONCLUSION

In this paper, we were interested in the problem of attitude angle estimation for a UAV with a catadioptric visual sensor. For ground truth images, we obtain mean errors equal to 1.5 • for the pitch and roll angles. These relatively important errors can be explained by the very difficult experimental conditions (very low altitude and urban environment).

Fig. 1 .

 1 Fig. 1. Equivalence between the catadioptric projection and the two-step mapping via the sphere.

Fig. 2 .

 2 Fig. 2. Projection of the horizon in the sphere of equivalence and in the image plane.

Fig. 3 .

 3 Fig. 3. Relation between the projection of the horizon on the sphere and the roll (ρ) and pitch (ψ) angles.

Fig. 4 .

 4 Fig. 4. Example of a real image projection onto the sphere of equivalence.

Fig. 5 .

 5 Fig. 5. First row: Original image (a) ρ = 20 • and ψ = 10 • , (b) ρ = 30 • and ψ = -20 • . Second row: Horizon detection and estimation of the angles (c)ρ = 19.2 • and ψ = 8.9 • (d)ρ = 29.1 • and ψ = -20.6 •

Fig. 6 .

 6 Fig. 6. Pitch (a) and roll (b) errors in degrees over 49 images.

Fig. 7 .

 7 Fig. 7. Horizon detection with angles ρ = 20 • and ψ = 30 • . Horizon occultation by the buildings due to the low altitude generates a bad estimation ρ = 19.6 • et ψ = -34.2 • .

Fig. 8 .Fig. 9 .

 89 Fig. 8. Angle estimation with a noised image, real values ρ = -20 • ψ = -30 • . (a) Without noise (σ = 0), estimated angles ρ = -19.4 • ψ = -30.1 • . (b) With white noise σ = 30, estimated angles ρ = -20.15 • ψ = -28.8 • .

Fig. 10 . 7 •

 107 Fig. 10. Aerial sequence results (a) ρ = 26.2• ψ = 1.8 • , (b) ρ = 17.0 • ψ = 0.0 • , (c) ρ = 0.1 • ψ = 0.9 • , (d) ρ = -16.3 • ψ = 1.4 • , (e) ρ = -24.8 • ψ = 4.3 • , (f) ρ = 3.5 • ψ = -0.9 • , (g) ρ = -24.9 • ψ = 1.3 • , (h) ρ = 1.6 • ψ = 6.7 •

Fig. 11 .

 11 Fig. 11. Roll (red) and pitch (blue) estimation aver 4000 images of the sequence.

ACKNOWLEDGMENT

This work is supported by ANR -Project CAVIAR Ref: -France. We also would like to thank Pierre Détaille and Lionel Guyot for their help.