Damien Eynard
email: damien.eynard@u-picardie.fr

Pascal Vasseur
email: pascal.vasseur@univ-rouen.fr

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

Vincent Frémont
email: vincent.fremont@hds.utc.fr

Real time UAV Altitude, Attitude and Motion Estimation from Hybrid Stereovision

Keywords: UAV, hybrid stereovision, motion, attitude, altitude

Knowledge of altitude, attitude and motion is essential for an Unmanned Aerial Vehicle during critical maneuvers such as landing and take-off. In this paper we present a hybrid stereoscopic rig composed of a fisheye and a perspective camera for vision-based navigation. In contrast to classical stereoscopic systems based on feature matching, we propose methods which avoid matching between hybrid views. A plane-sweeping approach is proposed for estimating altitude and detecting the ground plane. Rotation and translation are then estimated by decoupling: the fisheye camera contributes to evaluating attitude, while the perspective camera contributes to estimating the scale of the translation. The motion can be estimated robustly at the scale, thanks to the knowledge of the altitude.

We propose a robust, real-time, accurate, exclusively vision-based approach with an embedded C++ imple-

Introduction

Unmanned Aerial Vehicles (UAV) have received a lot of attention over the last decade, in relation to command and on-board computer vision. Interest in these topics has focused largely on increasing UAV autonomy, which includes maneuvers such as landing and take-off. In this context a fast, robust, accurate estimation of critical parameters such as altitude, attitude, motion and velocity is clearly crucial for the control loop.

A number of sensors have been used to estimate these parameters. Altitude can be determined by GPS, altimeter, radar, or laser. Attitude can be provided by Inertial Measurement Unit (IMU), and motion by GPS or radar. Raw Global Positioning Systems (GPS) have a poor accuracy both vertically (accuracy between 25 and 50 meters) and horizontally (less than 15 meters). Moreover, GPS is sensitive to interruptions in transmission in urban environments. Radar can estimate both altitude and velocity with greater accuracy, but it requires active sensors that consume power. Altimeters are very widely used, but they depend on pressure variations, which implies an accuracy error of between 6% and 7%. Laser altimeters, on the other hand, are highly accurate, but they have very specific requirements concerning reflecting surfaces. Finally, IMU can provide indications of velocity, acceleration, attitude and orientation, but it is subject to drift and error accumulation.

As an alternative to these sensors, vision-based systems, thanks to the growth of computational capacity, can now rapidly estimate all these parameters, as well as performing other visual tasks. Furthermore, cameras remain compact, passive systems with correspondingly low energy consumption, and can provide information at high rates (up to 200Hz).

Most of the works using computer vision for UAV are based on the optical flow analysis obtained with a single camera [START_REF] Barrows | Optic flow sensors for mav navigation[END_REF]; [START_REF] Beyeler | Visionbased altitude and pitch estimation for ultra-light indoor aircraft[END_REF]; [START_REF] Chahl | Landing strategies in honeybees and applications to uninhabited airborne vehicles[END_REF]; [START_REF] Green | Autonomous landing for indoor flying robots using optic flow[END_REF]. Indeed, optical flow can be efficiently computed and consequently embedded for navigation purpose. These systems were bio-inspired by information processing mechanisms that have evolved in bees, and consist in deducing the altitude according to the optical flow, given knowledge of the camera's motion. [START_REF] Beyeler | Visionbased altitude and pitch estimation for ultra-light indoor aircraft[END_REF] is alone in proposing an estimation of pitch to cor-rect the optical flow, without which an unstable system might result. An original approach using a single perspective camera is to be found in [START_REF] Cheriany | Autonomous altitude estimation of a uav using a single onboard camera[END_REF]. The authors use a technique based on the learning of a mapping between texture information contained in a top-down aerial image and altitude values. This learning stage is performed over different kinds of terrain, and a spatio-temporal Markov Random Field is used. However, optical flow information is clearly not sufficient for a complete control of a UAV.

Related Work

As explained before, orientation and position parameters are essential for the control of the UAV and vision based methods can be very interesting because of the richness of the provided information. A hierarchy between these parameters can be established according to their importance in the control loop. The most important one deals with the attitude knowledge for the stabilization followed by the altitude for hovering and landing maneuvers and finally position and yaw angle for trajectory. Consequently, this ranking imposes the refreshing rate requirement for the control. In this work, even if we do not manage the control of the UAV, our aim is to propose a complete onboard vision system able to provide orientation and position parameters at the required frequency rate for a future loop control. Many works based on vision solely or vision fused with other sensors have been already proposed. A complete recent review is presented in [START_REF] Weiss | Monocularslam-based navigation for autonomous micro helicopters in gps-denied environments[END_REF]. These works can be divided in two categories which are respectively based on cameras placed in the environment or based on embedded cameras on the UAV. The former consists generally in a fixed camera network that supervises completely a 3D scene. The most well known system is Vicon and has been successfully employed in UAV control [START_REF] Mellinger | Recent advances in quadrotor capabilities[END_REF]. Such systems are particularly efficient with a high measure frequency, very accurate with no drift and also robust. However, the fixed configuration does not allow any exploration of a new and unknown environment. The second category deals with onboard camera and can also be divided according to two sections. The first one is based on the use of a known pattern placed in the environment. In Garciapardo et al. (2000); [START_REF] Saripalli | Vision-based autonomous landing of an unmanned aerial vehicle[END_REF]; [START_REF] Sharp | A vision system for landing an unmanned aerial vehicle[END_REF], the authors make use of a downward-pointing perspective camera to estimate altitude according to a predefined pattern laid out on the ground. Recently, a circular pattern has been used in order to measure the pose of a UAV [START_REF] Eberli | Vision based position control for mavs using one single circular landmark[END_REF]. This kind of approach is interesting, since it requires only one camera, provides a complete pose and can be used in real time.

This solution is the complementary of Vicon approach since the camera looks an a priori known environment while with Vicon, cameras of the environment tracks a known pattern placed on the UAV. The advantages and drawbacks are then similar. The other solution consists of onboard cameras that capture and treat data from an unknown and unstructured environment (no artificial landmark). In this way, it is possible to estimate orientation and position according to features detected in the environment. If some absolute features such as horizon, vertical direction or ground are used, it is then possible to perform the estimation in an absolute reference frame and to obtain absolute measures of attitude and altitude. For example to estimate the attitude from horizon line, [START_REF] Demonceaux | Robust attitude estimation with catadioptric vision[END_REF] propose a method based on an omnidirectional view while [START_REF] Dusha | Fixedwing attitude estimation using computer vision based horizon detection[END_REF] use a perspective camera. [START_REF] Thurrowgood | A vision based system for attitude estimation of uavs[END_REF] estimates the attitude using a fast histogram based method which can be adapted for any type of camera. For a complete absolute six degrees of freedom estimation, an initial geo-localized reference is then necessary. If used features have no particular identity, the orientation and position estimation is performed relatively to the structure formed by these features. Such methods (SLAM, SFM, Visual Odometry, ego-motion) are now well-known and many different algorithms have been already proposed successfully in [START_REF] Davison | Monoslam: Real-time single camera slam[END_REF]; [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]; [START_REF] Nistér | An efficient solution to the five-point relative pose problem[END_REF]. Then, works of [START_REF] Li | A linear approach to motion estimation using generalized camera models[END_REF]; [START_REF] Pless | Using many cameras as one[END_REF] can determine the motion from non-overlapping cameras. Nevertheless they need a least of 14-points or 17-points to retrieve the motion that increases the computation time and noise sensitivity. However, if these methods are generally well adapted for mobile robots or vehicles, a direct transfer to micro-aerial vehicles is not possible and requires specific adaptations and developments [START_REF] Lee | A benchmarking tool for mav visual pose estimation[END_REF]. Since most of the efficient methods are based on feature points, a possible adaptation consists in including some a priori knowledges or hypotheses in order to reduce the number of necessary points. In [START_REF] Kalantari | A new solution to the relative orientation problem using only 3 points and the vertical direction[END_REF]), Fraundorfer et al. (2010), [START_REF] Naroditsky | Two efficient solutions for visual odometry using directional correspondence[END_REF], the authors propose a relative pose estimation between two views by a three-plus-one algorithm, based on the correspondence of three points, given a common direction in the two views. This common direction can then be obtained by IMU for example. A recent work proposed also to use only one point for structure from motion [START_REF] Scaramuzza | 1-point-ransac structure from motion for vehicle-mounted cameras by exploiting non-holonomic constraints[END_REF]. However, the used constraint for reducing the number of points is only valid for mobile robots.

The most recent and successful autonomous microaerial vehicle based on vision with adapted methods is described in [START_REF] Blösch | Vision based mav navigation in unknown and unstructured environments[END_REF]; [START_REF] Weiss | Monocularslam-based navigation for autonomous micro helicopters in gps-denied environments[END_REF].

The authors propose a monocular Simultaneous Localization and Mapping (SLAM) framework from which they extract pose parameters for a complete control of the six degrees of freedom. The use of a SLAM approach allows to correct the possible drift by loop closing. Since their system is based on a single wide field of view camera, the depth information can not be directly recovered. This is why they propose the initialization of the map scale by hand and maintain the consistency of this scale by fusing IMU data and camera poses through an EKF [START_REF] Nützi | Fusion of imu and vision for absolute scale estimation in monocular slam[END_REF]. The complete pose from the SLAM is computed at 15-30 Hz and the comparison of results with a VICON shows the effectiveness of the approach. In [START_REF] Artieda | Visual 3-d slam from uavs[END_REF], a monocular vision sensor is also proposed in order to perform the visual SLAM in partially structured environments. Other sensors are just used in order to propose comparisons. In their work, the authors deal with a deep study of the different steps of the visual SLAM such as the calibration, the influence of depth and the image processing techniques for feature detection and matching for example. Their results are also shown to be globally satisfactory but the treatment is made off-line at an average of 12 FPS. A discussion of the scale estimation is proposed but no real results on this parameter are shown. Recently, a visual monocular SLAM algorithm has been proposed especially for MAVs by integrating a plane constraint [START_REF] Lee | Mav visual slam with plane constraint[END_REF]. This constraint allows to obtain an approximately constant time algorithm that has been tested in indoor environments. A work quite more closer to ours consists in using horizon and points of the ground in order to estimate the ego-motion of a UAV [START_REF] Oreifej | Horizon constraint for unambiguous uav navigation in planar scenes[END_REF]. The use of the horizon line allows to linearize and disambiguate the planar flow, and consequently to obtain a unique solution for the UAV motion estimation. However, this work has only been tested on a balloon-based UAV with motion particularly different from a quad-rotor.

All these works propose very interesting results but they are based on a single camera and either are up to scale or require a manual initialization and the use of other sensors for maintaining this scale [START_REF] Weiss | Monocularslam-based navigation for autonomous micro helicopters in gps-denied environments[END_REF]. In order to obtain a complete vision approach, we propose in this paper to use a stereovision sensor for the estimation of the orientation and position parameters at scale. Due to the limited field of view of the perspective cameras, it has been demonstrated that the use of omnidirectional cameras is more suitable for motion estimation [START_REF] Lee | Large motion estimation for omnidirectional vision[END_REF]; [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF]. Indeed, omnidirectional cameras can gather more information from the environment but are less accurate owing to their limited resolution. In order to have both the wide field of view and the accuracy, we then propose an embedded hybrid stereovision sensor. However, classical stereovision based methods employ feature point matching between views, which generally constitutes a bottleneck in the whole process. In our case, this difficulty is compounded because of the heterogeneous nature of our images that prevents a direct comparison between points without any adaptation. Thus, in order to reduce this difficulty, a first solution consists in using the epipolar geometry between views and to limit the search area for each point to match. Recently, Puig et al. proposed an hybrid fundamental matrix and its use for hybrid matching [START_REF] Puig | Matching of omindirectional and perspective images using the hybrid fundamental matrix[END_REF]. However, their results show that very few points are positively matched and a complete and accurate pose estimation seems to be complicated. A second way for improving the results consists in adapting the point features to the geometry of the sensor. Many different works have been already proposed in order to adapt Harris corner detector [START_REF] Demonceaux | Central catadioptric image processing with geodesic metric[END_REF] or SIFT descriptors [START_REF] Bastanlar | Effective structure-from-motion for hybrid camera systems[END_REF], [START_REF] Lourenço | Feature detection and matching in images with radial distortion[END_REF]. However, the principal drawback is the expensive computation time that totally prevents a real time implementation on an UAV. We then propose in this paper an hybrid stereovision approach that avoids point feature matching between hybrid views and which is based on a direct comparison of images by assuming that the planar ground surface is visible in both images. The principal contributions of this paper are fourfold. The first is the hybrid sensor comprising fisheye and perspective views. Secondly, we have the use of correspondence-less methods. This leads to the third contribution, namely the real-time implementation of altitude, attitude and motion, adapted for future embedded applications. Finally, we have an estimation of the motion performed independently in each of the views, at a metric scale. The estimated motion is then merged and filtered by a Kalman filter.

The organization of the article is as follows. Part III deals with the hybrid sensor and its modeling. Part IV describes the estimation of attitude, altitude and motion. Part V is devoted to a quantitative evaluation of data obtained using a small UAV.

3 Hybrid sensor

Motivation

In this paper, we present a new approach for state (orientation and position) estimation of a MAV based on an hybrid stereovision sensor. Contrary to SLAM frameworks that provide simultaneously all the parameters, our method proposes a cascading approach that allows to obtain attitude, altitude and motion at different rates (see Fig. 1(c)) . Attitude is then first obtained from the omnidirectional camera by use of the horizon or vertical direction but can also be obtained by IMU if available. This parameter being the most important for flight stability, our system is able to provide it at 30Hz for the vanishing points method and 100Hz for the horizon line based method. The use of a stereovision sensor allows also to directly estimate the altitude at scale and consequently to avoid any manual initialization or any knowledge about the dimension of the scene or the use of an additional sensor for maintaining this information during the flight. We assume a consistent planar ground similarly than in [START_REF] Lee | Mav visual slam with plane constraint[END_REF] which is reasonable in most of the cases. Therefore, we will show in the results than imperfect grounds have low impact on the altitude estimation and that little obstacles can be present in the image. In the same way, we show experiments in order to show the accuracy of the system according to the ratio between altitude and baseline.

Thus, our hybrid sensor thus combines nested methods.:

-First, we calibrate the hybrid system using [START_REF] Caron | Multiple camera types simultaneous stereo calibration[END_REF]. -Second, attitude is estimated either by vanishing point detection in urban and indoor environments [START_REF] Demonceaux | Uav attitude computation by omnidirectional vision in urban environment[END_REF], or by the horizon fitting in other cases [START_REF] Demonceaux | Robust attitude estimation with catadioptric vision[END_REF]; -Third, as a result of the previous steps, knowing that a homography exists between the two views, altitude estimation and ground plane segmentation are performed by plane-sweeping [START_REF] Eynard | Uav altitude estimation by mixed stereoscopic vision[END_REF]. -Finally, motion is decoupled in rotation and translation. Rotation is known by attitude or IMU, and translation is estimated by a two-point algorithm in hybrid views [START_REF] Eynard | Uav motion estimation using hybrid stereoscopic vision[END_REF]. By knowing altitude previously estimated, the translation is provided at the metric scale.

By implementing these methods in C++ we show that our methods are real time and ready to be embedded in a UAV.

Modeling

Perspective Model

The perspective projection (Fig. 2(a)) models pinhole cameras. Let X = X Y Z T be a 3D point expressed in the camera frame. It is projected onto the image plane as x = x y 1 T :

x = pr(X) with x = X Z y = Y Z .
(1)

x is the point on the normalized image plane and u = u v 1

T , the pixel point, is obtained by the relation u = Kx. K is the intrinsic matrix, knowing parameters γ p = {p x , p y , u 0 , v 0 }:

K =   p x 0 u 0 0 p y v 0 0 0 1   (2)
The full perspective projection of a 3D point to the pixelic image plane is therefore pr γp (X) = Kpr(X).

Spherical Model

Fisheye lenses cannot be classified as single viewpoint sensors [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. Nevertheless, we can show that this modeling is a good approximation of [START_REF] Ying | Can we consider central catadioptric cameras and fisheye cameras within a unified imaging model[END_REF]. [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF] have proposed a calibration method based on a unitary sphere. This model is particularly accurate, and allows the radial and tangential distortions caused by the fisheye lens to be modeled. Using the spherical model (Fig. 2(b)), a 3D point X is first projected onto a unitary sphere, centered at 0 0 0 T . The obtained point, x s , is then perspectively projected onto the image plane as x, from the second center of projection 0 0 -ξ T (2(b)):

x s = X ρ , x = pr ξ (X) with x = X Z+ξρ y = Y Z+ξρ , (3)
and

ρ = √ X 2 + Y 2 + Z 2 .
The image point is obtained from a 3D point, knowing the intrinsic parameters γ s = {p x , p y , u 0 , v 0 , ξ} and using pr γs (X) = Kpr ξ (X) (eq. 1).

Calibration

The hybrid calibration is based on Virtual Visual Servoing work, which is in turn based on the calibration of perspective cameras Marchand and Chaumette. [START_REF] Caron | Multiple camera types simultaneous stereo calibration[END_REF] simultaneously estimate the projection parameters and relative poses between N cameras of N models composing the stereo rig. They have demonstrated that calibration results for standard cameras are similar to Bouguet's toolbox or Mei's toolbox using same points in the image. This method is implemented in the Hyscas software [START_REF] Hyscas | Hybrid stereoscopic calibration software[END_REF]. We use this software to calibrate our hybrid model. Then, the extrinsic parameters obtained from hybrid calibration are the rotation R c and the translation t c between the two cameras.

In the case of attitude estimation using IMU, the IMU and fisheye calibration is estimated by the InerVis software Ine. The matrices provided by IMU are directly corrected in order to be expressed in the fisheye view.

Attitude

We can distinguish three main approaches for estimating attitude using two sensors.

Vision in urban environments: attitude and orientation can be estimated in catadioptric views [START_REF] Bazin | Uav attitude estimation by vanishing points in catadioptric image[END_REF]; [START_REF] Demonceaux | Uav attitude computation by omnidirectional vision in urban environment[END_REF]. These works are based on vanishing points resulting from the projection of vertical and horizontal lines in 3D space onto the unitary sphere. The algorithm consists of first detecting edges (using Canny for instance) in the omnidirectional image. Then a step of split and merge keeps long edges by deleting short lines. Once lines have been extracted, their normals are computed from the unitary sphere, which is known by calibration. Finally, using a voting approach, the three main directions are computed from both vertical and horizontal lines. Views are modeled by the proposition of [START_REF] Barreto | Issues on the geometry of central catadioptric image formation[END_REF].

Vision in natural environments: Demonceaux et al. ( 2006) also proposes estimating attitude via a projection of the horizon line onto a catadioptric view. This method is able to estimate pitch and roll. By using one of [START_REF] Kalantari | A new solution to the relative orientation problem using only 3 points and the vertical direction[END_REF]; [START_REF] Montiel | A visual compass based on slam[END_REF]; [START_REF] Scaramuzza | Correcting vehicle heading in visual odometry by using image appearance[END_REF] methods, it is also possible to retrieve the heading.

IMU: attitude can be obtained from Euler angles.

In our experiments we estimate attitude using either a fisheye view or IMU.

Hardware

The hybrid stereo rig is composed of two IDS uEye cameras with M12 lenses, the first is a perspective lens while the second is a fisheye. The stereo baseline is 32 cm. Each camera provides 752 × 480 RGB images, and the two cameras are located at opposite extremities of the rig, pointing downwards. Frame grabbers are triggered by hardware, and ring-buffer acquisition allows continuous data transfer and non-stop frame acquisition.

An XSens IMU is attached to the rig to determine the attitude ground truth. An embedded board with an Atom CPU receives and records images as well as IMU data. Images are sent and the UAV acquisition controlled via a wifi interface. The Atom board, the cameras and IMU are embedded on an MD4-1000 UAV (see Fig. 1(a)). For processing, we use a Macbook Pro with Core 2 Duo at 2.4 GHz.

Altitude

Plane-sweeping was introduced by [START_REF] Collins | A space-sweep approach to true multi-image matching[END_REF], proposed in real time by [START_REF] Gallup | Real-time plane-sweeping stereo with multiple sweeping directions[END_REF] and then adapted to hybrid views in [START_REF] Eynard | Uav altitude estimation by mixed stereoscopic vision[END_REF]. First a reference view has to be defined. Then, for each normal n and for each distance to a 3D plane, the different warped images are compared (Eq. 4) to the reference image by homography. Let I and I ref be respectively the current and the reference image, I(x) the intensity of the normalized pixel x. Then, I * is the homography of I by the homography matrix H. For a pixel x at the altitude d we get:

I * (x, d) = I(H x,d ) = I((R c -t c n T d )x) (4)
The best estimation of the homography H corresponds to the minimum global error of the difference between the warped image and the reference view |I * -I ref |. In our application we extend this aspect. We take the perspective view as the reference I ref .

Neighborhood manipulations are easier to perform and faster than on a sphere. Images I obtained via the fisheye camera are projected onto the sphere and then onto the reference plane by homography, and are denoted I * . Since the cameras are pointing towards the ground, the centered working surface on the fisheye is subject to fewer distortions and gives a better resolution than at the edge of the image. Finally, plane-sweeping can be performed (see Fig. 3(a)). We propose estimating the altitude d and segmenting the ground plane by mixed plane-sweeping, with R c , t c , n known by calibration and attitude estimation. We first present the homography used in our hybrid case, before describing the mixed plane-sweeping algorithm.

4.1 Sphere to plane homography Given a mixed stereo rig modeled by a plane and a unitary sphere, we shall now define the homography of the 3D ground plane that exists between two views (see Fig. 3(b)). In [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], a homography links two projections of a 3D plane onto two planes. In [START_REF] Mei | Homography-based tracking for central catadioptric cameras[END_REF], homography relates two projections of a 3D plane onto two spheres. The homography is defined by the following equation:

H = R c -t c n T d (5)
We obtain x * ∼ H -1 x for a homography between two planar projections of a 3D plane. These two planar projections are then replaced by a planar and a spherical projection, giving Eq.6.

x ∼ x * s = ρ Hx s (6)

Homography H depends on R c , t c , the attitude n and the altitude d. R c , t c and n being known, the next step is to obtain the estimation of altitude that minimizes the error defined by Eq.7. It will be remarked that in this step we can extract the ground plane using the mask G. All the points in the perspective view satisfying H can be considered as belonging to the ground plane (Eq.8). In the following subsection we present the algorithm for segmenting the ground plane and estimating the altitude d.

d = argmin d ( x |I ref (x) -I * (x, d)|) (7) G = {x, x |I ref (x) -I * (x, d)| x (I ref (x)) < thres} (8) 4.2 Algorithm
Altitude is determined by an exhaustive search. Let [d min , d max ] be the altitude range in the estimation. For each iteration, the best altitude d is estimated from this range, then the mask G is updated (see Algorithm 1). Pixels corresponding to the ground plane are in white in figures 8(c) and 10(e). In order to obtain a real-time method we estimate the altitude at time t using the altitude at (t -1).

The estimation is performed in two steps:

Algorithm 1 Altitude and ground plane segmentation algorithm -initialization

Estimation(d min , d max , s, ∆d) {Initialization} a 0 = d -1 = d min b 0 = d 0 = d max G = {p P ixels} while | d k -d k-1 | > ∆ d do {Estimation of the best altitude} d k+1 = argmin d∈{ b k -a k s-1 t+a k ;t [0,s-1]} x G |I ref (x)- I * (x, d k )| {Estimation of the mask of inliers/outliers} G = {x, x 1 Wx |I ref (x 1 )-I * (x 1 , d)| x 1 Wx I ref (x 1 )
< thres} {Estimation of the new range depending on sampling}

a k+1 = d k+1 -b k -a k s-1 b k+1 = d k+1 + b k -a k s-1 k = k + 1 end while Return d
-Initialization: we estimate the best altitude from within a wide altitude range. -During flight we use the altitude estimated in the initialization phase to obtain a narrower range in Eq.10) with r d in Eq.9. This range depends on the vertical velocity v v of UAV (about ±5m/s) and the hardware's computational power in frames per second denoted f ps.

r d = v v f ps 10 3 (mm) (9) d t [d t-1 -r d , d t-1 + r d ] (10)
We now show how calibration, attitude, estimated altitude and the ground plane segmentation, once obtained, can be used to estimate autonomously the motion using the hybrid views, without any feature matching between hybrid views.

Motion estimation

To estimate the motion we need the calibration (R c , t c ), the altitude d at the meter scale, the attitude and the ground plane segmentation G. The attitude can be determined either by IMU or by vision using vanishing points in urban environments, or the horizon line in natural environments. In the case of natural environments the UAV's heading has to be computed a posteriori, using such methods as [START_REF] Kalantari | A new solution to the relative orientation problem using only 3 points and the vertical direction[END_REF]; [START_REF] Montiel | A visual compass based on slam[END_REF]; [START_REF] Scaramuzza | Correcting vehicle heading in visual odometry by using image appearance[END_REF]. Since altitude, attitude, ground plane and calibration are determined by vision or IMU, we decouple the rotation and the translation of the motion. Below we present an algorithm for estimating the translation from two sets of points (Fig. 4(b)), that is to say points located on the ground plane G and points located randomly in the environment with depth either known or unknown. Our algorithm is based on the technique of point tracking in hybrid views for estimating the translation from two views. The estimated translation t is then filtered and smoothed by a Kalman filter.

Motion of the stereo rig

In each image, tracking is performed by the method given in [START_REF] Bouguet | Pyramidal implementation of the lucas kanade feature tracker description of the algorithm[END_REF]. We define x t a tracked point in the image acquired at time t. Each point is related by a rotation R t relative to the world reference X w estimated by the IMU (see Eq.11).

x t = t R w X w ⇐⇒ X w = t R -1 w x t (11)
Then, for a couple of points in the image (x t ; x t+1 ) that indicate motion between times t and t + 1, we can express the point x t from the frame t to the frame t + 1 by Eq.12. We thus obtain a rotation corresponding to motion (eq.13).

x t+1 = t+1 R w t R -1 w x t (12) R t+1 = t+1 R w t R -1 w (13)
The points are related by a motion composed of the rotation R t+1 previously defined and the translation t t+1 (Fig. 4(a), Eq.14). The distance of the image point to the 3D point is defined by the altitude d (see Fig. 4(a)) known by plane-sweeping in the case of a point belonging to the ground plane.

x t+1 = d R t+1 x t + t t+1 (14)

Motion from ground points

As mentioned previously, plane-sweeping is used to estimate both the altitude and the segmentation of the ground plane G. Knowledge of these parameters entails knowledge of tracked points in the two views belonging to the ground plane, together with their depth (see Fig. 4(a)). Then, the motion is estimated at the metric scale and defined as follows, with × the cross product of two vectors:

(x t+1 × t t+1 ) = -d(x t+1 × R t+1 x t ) ( 15
)
We now show how motion can be estimated in each view and then in the mixed view. First, in the spherical an referential case, we obtain the matrices A s Eq.16 and B s Eq.17. Secondly, in the perspective case, images are considered to be rectified by the rotation R c , and we obtain the matrices A p Eq.16 and B p Eq.18. In the hybrid case A s and A p are concatenated in A h , and B s and B p are concatenated in B h . Finally, we obtain Eq.22.

A s/p =   0 -x t+1 x t+1 z t+1 0 -x t+1 -y t+1 x t+1 0   (16) B s = -d s (x t+1 × R t+1 x t ) T (17) B p = -d p ((x t+1 × R t+1 x t ) + x t+1 × ((R t+1 -I d )t c )) T (18)

Motion from environment points

Concerning points located randomly in the environment, without possessing any knowledge of their depth we can extend to the mixed view the motion estimation proposed by [START_REF] Bazin | Motion estimation by decoupling rotation and translation in catadioptric vision[END_REF] for the spherical view (see Fig. 4(a)). The translation t t+1 from two points is defined as in Eq.19 . As mentioned previously, perspective points and spherical points are concatenated in order to estimate the translation t t+1 for obtaining the relation in Eq.23.

(R

t+1 x t × x t+1 ) T • t t+1 = 0 (19)

Combining 3D points from mixed views

Note: In this part, we adapt the generic formulation to the case of perspective case. x pi,t denotes the i th point tracked in the perspective view at time t, and x si,t denotes the i th point tracked in the sphere at the time t.

C h =           R t+1 x s0,t × x s0,t+1 . . . R t+1 x sm,t × x sm,t+1 R t+1 x p0,t × x p0,t+1 . . . R t+1 x pn,t × x pn,t+1           (20)
Considering perspective images rectified by R c known by calibration, we get:

D h =           0 . . . 0 (R t+1 x p0,t × x p0,t+1 ).(R t+1 t c -t c ) . . . (R t+1 x pn,t × x pn,t+1 ).(R t+1 t c -t c )           (21) A h t t+1 = B h (22) C h t t+1 = D h (23)
We have already presented a generic formulation for translation estimation, either from ground points or from environment points. When translating from planar points, the translation is estimated at the metric scale. However, the main drawback in the case of 3D motion is the pixel projection noise sensitivity. When translating from environment points, on the other hand, motion estimation has the advantage of being more robust to noise than on the plane, but the estimation is performed up to scale. One of the contributions of this paper is thus the combination of the two methods to increase both accuracy and robustness. The first method is defined by Eq. 22, and the second by Eq. 23. C h is defined by Eq. 21. Concatenating the matrices A h and C h yields A f us , and concatenating B h and D h gives us B f us . Finally we obtain the A f us t t+1 = B f us solved by least squares. points and ground points. The introduced error is a random angle, applied to view-lines whichever the projection model used. The angle error is computed between the real vector of translation and the estimated vector of translation.

Kalman filtering

Once the motion has been estimated at the metric scale we observe certain discontinuities and brutal variations. In order to reduce bad estimations and to refine the trajectory, we have chosen to use a linear Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. The considered state is simply the translation vector of the ego-motion i.e. X k = t t+1 , and is modeled as a linear Gaussian system given by Eq. 24:

X k+1 = EX k + V k Y k = OX k + W k ( 24
)
where E = Id 3 is the linear state transition model and O = Id 3 the observation model, assuming constant speed during a sample interval. The vectors V k and W k correspond respectively to the error model and the observation noise. They are assumed to be additive and white zero-mean Gaussian, with user-defined covariance matrices. In order to retain the dynamic nature of the measurement, more uncertainty is given to the measures, i.e. 10 6 and 10 2 for the model. Given the above considerations, the Kalman filter consists in predicting the translation vector t t+1 and then obtaining a refined value using an update step when a new observation is available. We therefore obtain the translation vector of the ego-motion and its estimated accuracy, from all past observations up to the current time.

Results

In this section we present the results of the different methods. First, a ground truth comparison is used to validate the effectiveness of each method and also its accuracy, then each method is tested on UAV image sequences. Finally, we analyze the performance of the algorithm.

Attitude

First, attitude can be provided either by vision or by IMU. For the case of vision, two types of methods are analyzed: in urban environments and natural environments. Results are compared to an IMU.

Attitude estimation by vanishing points

We employ the method of [START_REF] Demonceaux | Uav attitude computation by omnidirectional vision in urban environment[END_REF] and [START_REF] Bazin | Uav attitude estimation by vanishing points in catadioptric image[END_REF]. The error introduced by this method does not exceed 3 • . Our algorithm is insensitive to low attitude errors. For example, with real images, an error of 5 • for the roll and the pitch will introduce respectively an error of 0.87% and 0.39% (fig. 6(b)) for altitude estimation. This method, which has been tested on our fisheye lens, is as good as a catadioptric lens at estimating attitude. In the works presented here, attitude can be estimated in real time (around 30 Hz).

Indoor experiments were carried out using a Pelican UAV, where vanishing points can be determined. Fig. 6(a) and Table 1 present the results the from visionbased method (blue line) and IMU (red line). Notice that the red curve is smoothed by the IMU Kalman filter. Finally, for the yaw, the absolute offset between the vision estimation and the IMU is removed by : y i,vp = y i,vp -ȳvp + ȳimu , with y i,vp the i th value estimated by vanishing points, ȳimu the mean of yaw angles estimated by IMU.

Altitude

Once the calibration and attitude have been estimated we proceed to estimate the altitude. We compare our algorithm to a laser telemeter, and UAV altitude estimation validates the effectiveness of this approach.

Comparison with the ground truth

We present two cases of experimental results where real altitude is estimated by a laser telemeter. The error between actual and estimated altitude is expressed as = ( d-dgt) dgt , with d gt the ground truth of the altitude.

-Case 1: two cameras with a 447mm baseline are fixed on a pneumatic telescopic mast. Altitude and ground plane estimation are performed offline on a GPU. -Case 2: two micro cameras with a 314mm baseline are embedded on a compact UAV. Altitude estimation is performed online by CPU processing.

In the first experiment we observe an accurate estimation of altitude on a free ground plane (tab. 2), with an error between 0.18% and 3.14%. When there are obstacles on the ground plane, we observe a higher error, between 7.52% and 8.82%. The higher the altitude, the less accurate the estimation because of the decrease in resolution as altitude increases. Moreover, the accuracy depends on the size of the baseline. This fact does not impact our application negatively, given that accuracy is most critical during the two phases of landing and take-off, i.e. where the UAV is close to the ground plane. The second aspect of our algorithm is the segmentation of the ground plane, which is well estimated for contrasted areas. In case of a plane without obstacles, the pneumatic telescopic mast, where cameras are fixed, is well represented by outliers (dark on the image) (Fig. 8(c)). For an image composed of a dominant ground plane and walls, the ground plane is segmented as inliers, while walls are segmented as outliers. Our algorithm allows inliers/outliers to be segmented globally when estimating the dominant ground plane for our application.

The aspect of our algorithm in greatest need of improvement is where there are poorly textured planes and reflections which depend on camera orientation. When the ground plane or outliers (walls, objects) are homogeneous or poorly textured, the segmentation of outliers/inliers becomes difficult. For the second experiment we implemented our system on a small quadrotor UAV. Micro cameras embed-ded on the UAV are linked to an external laptop to perform online altitude estimation. We tested the accuracy by comparing altitudes estimated by plane-sweeping to altitudes estimated by laser telemeter. Fig. 9 shows that altitude is well estimated for the range of altitudes corresponding to the landing and take-off phases of a UAV. The mean error is 2.41%.

UAV Altitude Estimation

UAV Ground Plane Segmentation

By estimating the altitude from the two views (Fig. 10(a),10(b) and the homography (Fig. 10(c)), the ground plane (Fig. 10(e)) is segmented using the error (Fig. 10(d)) and becomes the mask G. Figure 10(e) shows the good segmentation of planar (white pattern) and non-planar objects (cables in black) after the threshold.

Performances

First we developed the algorithm on GPU with brook+ for ATI that gives a real-time (30 Hz) frame rate in order to estimate altitude and to segment the ground plane. The algorithm was tested on ATI 4850 with an E8400 3Ghz CPU.

We then implemented the algorithm on a CPU with a good deal of optimizing and without segmenting the ground plane. With this implementation we obtained min : 80 Hz, mean : 180 Hz, max : 250 Hz, that is to say values higher than the video frame rate, which meant that we could use our algorithm for online processing. The platform for our tests is a Macbook Pro with a C2D P8400 2.26Ghz CPU. A demonstration can be seen at Eynard (2010). We use a stereo rig with uEye cameras and get the normal of the ground plane with an IMU. During this demonstration, the system is able to estimate altitude in real time with robustness and accuracy.

An embedded version of our algorithm was exported to the ARM of a Gumstix Overo Fire with OMAP3530 ARM @600Mhz based processor. With this implementation we get a frame rate of around 5 Hz, which is not enough for real-time applications, but relatively interesting given the power/size ratio.

By developing these algorithms on GPU, CPU and an embedded board we obtained interesting results. For simultaneous ground plane segmentation and altitude estimation, results are real time and can be implemented on a UAV with GPU. For altitude estimation only, computation is faster and can be implemented on smaller quadrotor UAV.

Motion

We propose estimating motion using a linear method, namely least squares (LS). This method is robust to Gaussian noise, but sensitive to outliers [START_REF] Malis | Méthodes robustes d'estimation pour la vision robotique[END_REF]. For this reason, outliers are rejected by the RANSAC method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. This method was first validated on robot arm performing a loop, as we now describe.

Comparison with the ground truth

The system is presented in Fig. 11(b), and images are processed offline. In each view the number of tracked points is between 50 and 200. In 11(c), we compare our algorithm with the ground truth, up to the micrometer, obtained by the robot arm. The trajectory is a square of side 330 mm, with a total motion of 1320 mm. All trajectories are expressed within the same coordinate system of the robot arm. First we estimate the altitude by plane-sweeping (Fig. 10(f)). Then the motion is estimated using the derotation obtained from IMU (blue line) or from the ground truth (green line). In the blue case, the error is 19.084 mm in XYZ, or 1.45%. In the green case, the error is 15.865 mm in XYZ, or 1.20%.

UAV motion estimation

Once it had been validated, we tested our algorithm on a quadrotor with two uEye cameras with images processed offline. This validation was performed both indoors and outdoors. A Xsens IMU provides the attitude and the rotation of the motion while altitude is estimated via plane-sweeping. In each view, the number of tracked points is between 50 and 200. We assume that the perspective camera is pointing at the ground and looking at planar points, while the fisheye observes environment points. In an indoor environment, Fig. 12 In the outdoor environment the ground texture, mainly grass, means that the plane-sweeping algorithm cannot give accurate results, because of ambiguities in the search range. The poorly textured ground, the high altitude (15 m), and the short baseline make altitude estimation difficult. In this case we hand altitude estimation over to a pressure altimeter, embedded on the UAV. Fig.14(a) shows the motion estimated by GPS filtered by the altitude provided by the pressure altimeter. We then estimate the motion using embedded IMU to estimate the attitude and the rotation of the motion.

Performances

The algorithm was implemented in C++ and tested on a Macbook Pro. Estimating rotation by IMU and with tracking points on two images of 752 × 480, we could estimate the translation up to 24 Hz. The Lucas Kanade tracker has been adapted [START_REF] Sanahuja | Commande et localisation embarquée d'un drone aérien en utilisant la vision[END_REF] for the OMAP 3530 DSP, that allows real-time tracking for one image per Gumstix.

Conclusions and future works

This paper presents a hybrid stereo system composed of fisheye and perspective views, modeled respectively by a sphere and a plane, which is able to estimate navigation parameters for a UAV without the use of any other sensor. Navigation parameters are attitude, altitude, ground plane segmentation and motion, assuming calibrated cameras.

Projecting the horizon line can determine the attitude, as well as vanishing points obtained from lines in an urban environment. The projections of the ground plane onto the sphere and the plane are related by a homography which allows both altitude and ground plane to estimated using plane-sweeping. Once the altitude has been estimated and the ground plane segmented, tracked points on the ground provide the information for metric translation, while environment points give an accurate estimation of translation up to scale. By combining points from mixed views, and by merging two sets of points, translation is estimated accurately using the metric information, and smoothed by Kalman filter.

The methods we propose are real-time methods that do not require matching between hybrid views. We show in this article that our sensor is able to solve the main drawbacks of conventional sensors by using a compact, low power-consumption, passive, accurate, robust vision sensor. Furthermore, once navigation parameters have been estimated, this system is available to perform other visual tasks such as surveillance and environment discovery.

Perspectives of this work will be the development of a new hardware architecture able to grab stereo images and provide all the parameters in parallel, as well as integrating the estimated parameters into the control loop.

Fig. 1

 1 Fig. 1 (a) UAV with hybrid system outdoors. (b) Pelican in an indoor environment. (c) General method proposed.

Fig. 2

 2 Fig. 2 (a) Perspective model used for pinhole cameras. (b) Spherical model used for omnidirectional cameras.

  Fig. 3 (a) Mixed plane-sweeping. (b) Sphere/plane homography.

Fig. 4

 4 Fig. 4 fig. 4(a) 3D points located on the ground plane and in the environment. Fig. 4(b) Combining points located in the environment and the ground plane from hybrid views.

Fig. 5

 5 Fig. 5 Comparison of Ground Points, 3D Points and 3D + Ground Points. Dashed lines represent estimation from 3D of fused points while plain lines represent estimation from planar points. A smaller slope corresponds to a low sensitivity of the corresponding algorithm to the noise.

Figure 5

 5 Figure5illustrates the advantages of combining 3D points and ground points. The introduced error is a random angle, applied to view-lines whichever the projection model used. The angle error is computed between the real vector of translation and the estimated vector of translation.

Fig. 6

 6 Fig. 6 (a) Comparison between attitude estimated by IMU (red) and by vision (blue). (b) Influence of attitude error (pitch and roll) on the estimation of altitude.

  Fig. 7 (a) View from the UAV's fisheye. (b) Segmented ground/sky.

Fig. 9

 9 Fig. 9 Altitudes (cross) are estimated by plane-sweeping and the line represents measurements from the laser range finder based on himself, thus ideal measures.

Fig. 8 Fig

 8 Fig. 8 Altitude and ground plane segmentation -4.8% of inliers -Fisheye view (a), perspective view (b), ground plane segmentation (c), sphere to plane homography (d), reference and homography comparison (e)

  (a) presents the final 3D trajectory of the motion estimation. Raw data sensitive to noise and bad estimations are shown in red, while data filtered and smoothed by Kalman filter are shown in blue. In other experiment, Fig.13(a) presents the final trajectory in green, superimposed to an external view.

Fig. 11

 11 Fig. 11 Loop of a robot 11(c) with hybrid and merged views, Fig.11(a) fisheye and perspective views of the mobile robot.

  Fig. 12 3D trajectory of the UAV.

  Fig. 14 Taking-off in natural environment with MD4-1000. The red line represents GPS ground truth. The purple line is the estimated motion and green line the estimated motion filtered by Kalman filter (for more details, please refer to corresponding video).

Table 1

 1 Attitude Comparison for IMU vs. the Vision Method

	Type	Roll	Pitch	Yaw
	Error max	3.97 •	2.43 •	5 •
	Mean of absolute error 1.32 •	1.96 •	0 •
	Std of error	1.65 •	2.07 •	0.78 •

Table 2

 2 Altitude ground plane estimation with and without obstacles -Algorithm parameters for this test: s = 6, thres = 25

	Type	Grnd truth Estim. altitude	Error
	Free Ground	2187mm	2200mm	0.59%
	Free Ground	3244mm	3250mm	0.18%
	Gnd + obstacles
	(low contrast)	3244mm	3488mm	7.52%
	Free Ground	4072mm	4200mm	3.14%
	Free Ground	5076mm	5202mm	2.48%
	Gnd + obstacles	4080mm	4440mm	8.82%

Acknowledgments

This work is supported by the European FEDER (Fonds Européen de Développement Régional) and Région Picardie Project ALTO (Automatic Landing and Take-Off).


Jean-Luc Boudenne
email: jean-luc.boudenne@univ-amu.fr

Julien Parinet

Carine Demelas

Tarek Manasfi
email: tarek.manasfi@univ-amu.fr

Bruno Coulomb
email: bruno.coulomb@univ-amu.fr

Jean - Luc

MONITORING AND FACTORS AFFECTING LEVELS OF AIRBORNE AND WATER BROMOFORM IN CHLORINATED SEAWATER SWIMMING POOLS

Keywords: Trihalomethanes, bromoform, air, seawater, swimming pool, brominated compounds

Br-DPBs (brominated disinfection by-products); CHBr 3 (bromoform); CHCl 3 (chloroform), CHCl 2 Br (dichlorobromomethane), CHClBr 2 (dibromochloromethane); Clf (free residual chlorine); NPOC (Non Purgeable Organic Carbon); PCA (Principal Component Analysis); T water (temperature of water); T air (temperature of air); TTHMs (total trihalomethanes = sum of THM-4);

Introduction

Thalassotherapy (term coined by the French doctor de la Bonnardière in 1869) has begun establishing health therapy centers for early 1920's, at first in France and then in Great Britain and Germany [START_REF] Schwartz | Encyclopedia of coastal science[END_REF]. This practice has been defined by the French Medical Academy in 1960 as being the use of seawater, seaweed, sea mud or other sea resources and/or the marine climate for the purpose of medical treatment with a medicinal effect [START_REF] Crecente | A multicriteria approach to support the location of thalassotherapy (seawater therapy) resorts: Application to Galicia region[END_REF]. Till 1998, thalassotherapy was covered by French medical insurance as a standard treatment for sore throats, digestive problems, rheumatic arthritis, musculoskeletal injury, development disorders in children, respiratory problems (asthma and bronchitis), endocrine imbalances and skin diseases (e.g. psoriasis) [START_REF] Roullier | Thalassothérapie et crénothérapie pédiatrique[END_REF][START_REF] Shani | Sustained decrease of blood pressure in psoriatic patients during treatment at the dead sea[END_REF][START_REF] Charlier | The healing sea: a sustainable costal ocean resource: thalassotherapy[END_REF]. The decision to stop reimbursing these cures was not linked to the questioning of the healing properties of this practice, but rather because it has turned into tourism and leisure activities at the end of 1990's. Thalassotherapy is indeed now included in health or spa tourism, which already represents a $156 billion global industry but is still an emerging and rapidly growing sector [START_REF] Johnston | Global Spa Summit, wellness tourism and medical tourism: Where do spas fit? Global Spa Summit LLC[END_REF][START_REF] Cockerell | Spas and health resorts in Europe[END_REF]. Furthermore, health or spa tourism is currently one of the tourism products with the highest rates of growth [START_REF] Farnos | Condiciones óptimas para un centro de talasoterapia (Optimal conditions for a thalassotherapy resort)[END_REF].

Whereas some researchers question the healing properties of seawater cures [START_REF] Robiner | Psychological and Physical Reactions to Whirlpool Baths[END_REF][START_REF] Lotti | Mineral waters: Instead of soap or better than soap?[END_REF]Tsuchiya et al., 2003a;Tsuchiya et al., 2003b;[START_REF] Boulet | Algae-induced occupational asthma in a thalassotherapist[END_REF], another area of concern is about the impact of chemical treatment on the physical-chemical content of seawater. Indeed, supposed healthy benefits of seawater are partially connected to their mineral salt contents that are believed to be absorbed through the skin during immersion of users [START_REF] Kazandjieva | Climatotherapy of psoriasis[END_REF][START_REF] Moss | Water and health: a forgotten connection?[END_REF]. However, as there is no specific legislation concerning pools and other medical therapy equipments relating to hygiene and security rules, in France and in most of European countries, the thalassotherapy pools use the public swimming pool legislation [START_REF] Chapuis | Management of infectious risk associated with therapeutic pools[END_REF]ANSES, 2013). This lack of specific legislation on seawater swimming pools implies thus the use of chlorine-based compounds (in France, only these compounds are authorized as chemical disinfectants since November 2011) to prevent outbreak diseases caused by feces-associated microbes, such as viruses and bacteria (WHO, 2006;[START_REF] Papadopoulou | Microbiological quality of indoor and outdoor swimming pools in Greece: Investigation of the antibiotic resistance of the bacterial isolates[END_REF][START_REF] Schets | Exposure assessment for swimmers in bathing waters and swimming pools[END_REF]. Users represents not only a major source of micro-organisms but also inputs of organic compounds (human body fluids coming from urine, sweat, hair, mucus, skin and/or personal-care products) [START_REF] Seredyńska-Sobecka | Monitoring organic loading to swimming pools by fluorescence excitation-emission matrix with parallel factor analysis (PARAFAC)[END_REF][START_REF] Weng | Disinfection by-product dynamics in a chlorinated, indoor swimming pool under conditions of heavy use: National swimming competition[END_REF][START_REF] Keuten | Definition and quantification of initial anthropogenic pollutant release in swimming pools[END_REF][START_REF] Duirk | Reaction of Benzophenone UV Filters in the Presence of Aqueous Chlorine: Kinetics and Chloroform Formation[END_REF][START_REF] Santos | The degradation products of UV filters in aqueous and chlorinated aqueous solutions[END_REF]. The combination of a residual disinfectant need (residual of free chlorine higher or equal to 2 mg/L Cl 2 in France) with continuous inputs of organic compounds from users leads to the formation of disinfection by-products [START_REF] Judd | The fate of chlorine and organic materials in swimming pools[END_REF][START_REF] Weaver | Volatile disinfection by-product analysis from chlorinated indoor swimming pools[END_REF][START_REF] Kanan | Formation of disinfection byproducts in indoor swimming pool water: the contribution from filling water natural organic matter and swimmer body fluids[END_REF]Manasfi et al., 2017a).

Adding of chlorine in rich-bromide waters leads to the rapid oxidation of bromide ions (Br -) into hypobromous acid (HOBr), that is readily more reactive with organic matter than hypochlorous acid (HClO) [START_REF] Kampioti | The impact of bromide on the formation of neutral and acidic disinfection by-products (DBPs) in Mediterranean chlorinated drinking water[END_REF][START_REF] Fabbricino | Modelling disinfection by-products formation in bromide-containing waters[END_REF][START_REF] Tian | Chlorination and chloramination of high-bromide natural water: DBPs species transformation[END_REF]. The occurrence of brominated compounds in seawater swimming pools is poorly documented whereas the toxicity of brominated compounds is known to be greater than their chlorinated analogs [START_REF] Erdinger | Pathways of trihalomethane uptake in swimming pools[END_REF][START_REF] Richardson | Occurrence, genotoxicity and carcinogenicity of regulated and emerging disinfection byproducts in drinking water: a review and roadmap for research[END_REF][START_REF] Kogevinas | Genotoxic effects in swimmers exposed to disinfection by-products in indoor swimming pools[END_REF][START_REF] Plewa | Mammalian cell cytotoxicity and genotoxicity of the haloacetic acids, a major class of drinking water disinfection by-products[END_REF]. In a previous paper [START_REF] Parinet | Exposure levels to brominated compounds in seawater swimming pools treated with chlorine[END_REF], occurrence of brominated compounds were investigated in water of seawater swimming pools : bromoform (CHBr 3 ) was the most abundant trihalomethane (THM) compound found with levels ranging from 28.7 mg/L to 930.7 mg/L (272.1 mg/L median) followed by dibromochloromethane (CHCl 2 Br) with a median level of 14.85 mg/L (63.58 mg/L maximum). Bromoform has been classified as a group B2 (probable human carcinogen) by US.EPA (ATSDR, 1990). Because bromoform can be absorbed into the body through the skin but also by inhalation [START_REF] Beech | Nitrates, chlorates and trihalomethanes in swimming pool water[END_REF][START_REF] Xu | Percutaneous absorption of trihalomethanes, haloacetic acids, and haloketones[END_REF], it seems thus important to determine the exposure levels of seawater pools' users to bromoform (and to the other THM) not only in water but also in air.

This study has thus been carried out to determine the THM-4 (Total THM, TTHM) levels of exposure in air and water of four thalassotherapy centres located in Southeast of France during four sampling campaigns, spread on two days in summer and two days in winter.

Materials and methods

Thalassotherapy centres.

Four seawater establishments (called after E1, E2, E3 and E4) from Southeast of France, using chlorination agents (bleach or dichloroisocynauric acid) for water disinfection were selected for this study. E1 and E2 possess 2 pools (E1.A, E1.B and E2.A and E2.B); E3 possesses only one pool (E3.A) and E4 possesses 3 pools (E3.A, E3.B and E3.C). Pools E1.A, E2A, E4.A and E4.B are used only for relaxing or swimming activities (notice that E4.A is not connected to other basins -E4.B and E4.C-in this establishment and is not heated), whereas E1.B, E2.B and E4.C are used only for sport activities. E3.A is used equally for relaxing or sport activities. Globally, water temperature of pools devoted to swimming or relaxing are lower (29-32 °C) than those used for sport activities (33-34 °C). These temperature are preset and kept constant, whatever the season, in the range measured. All the swimming pools are disinfected continuously by bleach except E4.A and E4.C, which are disinfected once a week by dichloroisocyanuric acid. The hydraulic characteristics and mean daily frequentation of the whole pools are presented in Table S1.

These indoor-swimming pools received seawater from the Mediterranean Sea. Physicalchemical and chemical quality of raw seawaters feeding these pools are reported in Table S2.

Physico-chemical parameters (conductivity, pH, turbidity, water and air temperature, hygrometry, total residual chlorine) were recorded on-site at the time of sampling. Total residual chlorine was measured [START_REF] Urbansky | Ascorbic acid reduction of residual active chlorine in potable water prior to halocarboxylate determination[END_REF] and the other physico-chemical parameters were determined with specific electrodes. Free and combined bromine were determined by the N-N-diethyl-p-phenylenediamine (DPD) colorimetric method [START_REF] Lourencetti | Trihalomethanes in chlorine and bromine disinfected swimming pools: Air-water distributions and human exposure[END_REF]. Bromide and chloride were measured with an ICS-3000 Dionex ion chromatography system using a 30 mM NaOH eluent at 1.5 ml.min -1 and 30 °C.

Water and air were sampled in triplicate at the opening and, at the closing of the basins, during two consecutive days in winter (February 2016) and two consecutive days in summer (August 2016), to study the potential influences of basin frequentation (higher frequentation leading to higher organic load in the pools), and of potential supplementary organic inputs from bathers in summer (e.g., sweat, sunscreens). Indeed, as observed during previous campaigns, thalassotherapy user behaviour pattern in summer, alternatively consists in selftanning outside the pools and in cooling himself down inside the pools, bringing thus in the pools a cocktail of organic compounds ready to react with HOBr [START_REF] Manasfi | Identification of disinfection by-products in Freshwater and Seawater Swimming Pools and Evaluation of Genotoxicity[END_REF]Manasfi et al., 2017b). Moreover, people who attend thalassotherapy centres to be older in winter ("curist" period) than in summer ("non-curist" period, including mainly families accompanied with children). Overall, the number of bathers is the same in winter than in summer (Table S1).

Sampling and analysis of water and air

3 air and 3 water samples were collected per day and per pool (1×2 in the morning, 1×2 at noon and 1×2 at the closing of the pools. The total number of water (or air) samples is thus 6 (×2 for duplicates) in winter and 6 (×2 for duplicates) in summer.

Samples for TTHM analysis in the air were collected during 30 min on a 150 mg NIOSH charcoal tube (MSA AUER Gmbh, Berlin, Germany) connected to a pump Escort ELF (MSA AUER Gmbh, Berlin, Germany) at a flow rate of 1L.min -1 . The air sampling was spiked at 20 cm above the surface of the swimming pool. The charcoal was then eluted in 5 mL of toluene (Analytical grade, Sigma-Aldrich) in an ultrasonic bath for 5 min.

TTHM in air were analyzed by Gas Chromatography coupled to an Electron Capture Detection (GC-ECD). Separation of chloroform (CHCl 3 ), dichlorobromomethane (CHCl 2 Br), dibromochloromethane (CHClBr 2 ) and bromoform (CHBr 3 ) in air was carried out using a Perkin Elmer (Norwalk, CT, USA) GC Model Clarus 580 equipped with 63 Ni ECD system and split-splitless injector was used. Helium 5.0 was used as carrier gas at a programmed flow of 1 mL/min. For separation, a capillary column DB5-ms (30 m × 0.25 mm × 1 µm) Perkin Elmer (Norwalk, CT, USA) was employed. The calibration ranges, the limits of detection and the RSD were evaluated based on the addition of 50 µL of a THM's solution at 1 mg/ml on 7 charcoal tubes and where then extracted by 2 mL of Toluene. The detection limits were determined for an air sample of 30 L. The extraction efficiency was calculated based on the addition of 0.1 mL of a THM's solution at 1 mg.mL -1 on 4 charcoal tubes and then eluted by 2 mL of toluene. All these analytical features are given in Table S3.

Sample for TTHM analysis in water were collected in 60 mL glass bottles, filled until overflowing with pool water and tightly sealed with a PTFE-lined screw cap. The water sampling was carried out at 20 cm under the surface of the swimming pool. Prior to sampling, 0.1 mL of 100 g.L -1 ascorbic acid was added to each 60 mL bottles to quench any residual chlorine reaction [START_REF] Lee | Production of various disinfection byproducts in indoor swimming pool waters treated with different disinfection methods[END_REF]. Samples were also collected in 1-L glass bottles for analysis of Total Organic Carbon (TOC), Total Nitrogen (TN), UV absorbance at 254 nm (UV 254 ), bromide and chloride. All samples were refrigerated at 4 °C until analysis.

Analysis of TTHM in seawater swimming pools was carried out by headspace-gas chromatography/mass spectrometry (HS-GC-MS) in Electron Impact (EI) mode. Separation of chloroform (CHCl 3 ), dichlorobromomethane (CHCl 2 Br), dibromochloromethane (CHClBr 2 ) and bromoform (CHBr 3 ) in water was carried out using a Varian Saturn 2000 gas chromatograph (Varian Inc., Walnut Creek, CA, USA) equipped with electronic flow control (EFC) and a 1079 universal capillary injector that allows programmed temperature injection (a PTV injection port). The gas chromatograph was also equipped with an autosampler (CombiPAL autosampler, CTC Analytics) for 32 headspace vials composed of an oven for sample heating/headspace generation and a robotic arm where the headspace syringe was placed. Helium was employed as carrier gas. The operating conditions of the headspace are presented in Table S3. Then 1 mL of the vial headspace was injected in the GC inlet heated at 140 °C, operated in splitless mode and equipped with a split open deactivated insert liner of 5 mm OD × 54 mm × 3.4 mm ID (Varian Inc., Walnut Creek, CA, USA). The column temperature was set at 50 °C and remains constant during 5 min. After this time, the temperature raised up to 220 °C at 10 °C/min. A constant column flow of 1.5 mL.min -1 of helium was used. An Agilent FactorFourVF-624 ms fused silica column of 60 m × 0.32 mm i.d. and 1.4 µm of film size (Agilent, Santa Clara, CA, USA) was used for chromatographic separation. Full scan data were acquired under the following: mass range 45-260 m/z, scan time 1 s, solvent delay 10 min, manifold temperature 120 °C, and emission current 20 µA.

The limits of detection and the RSD, presented in Tables S3 andS4, have been evaluated based on 7 samples at 0.1µg/L of each THM. The recoveries have been evaluated based on 4 seawater samples at 1 and 20 µg/L of each THM.

Details about oven program, sources of chromatographic standards, quantifications, recoveries and methods of limits of detection of the analytes are reported in the Supporting Information (Table S4).

Statistical Analysis

Statistical analyses were performed using SPSS statistics 17.0 (49), downloadable software for analyzing data, distributed by Softonic®. The taking into account of a large number of variables and the comparative study of different waters from seawater swimming pools has imposed us to use multidimensional statistical data processing. Principal Component Analysis, which is the tool of choice for this type of study, was therefore selected.

The objective of Principal Component Analysis (PCA) was to help differentiate between the different physical-chemical parameters involved in the appearance of disinfection by-products in water and air. Principal Component Analysis aims to turn the original variables (physicalchemical parameters) into new variables, principal components (PC), which are orthogonal and non-inter-correlated. The PCs are ordered in such a way that the first component represents the largest fraction of the total data variance. In order to evaluate the influence of each variable in the PCs, a varimax rotation is generally used to obtain the rotator factor loadings, which represent the contribution of each variable to each of the principal components obtained. The contribution of each variable to the construction of a component can be evaluated through the coefficient of correlation (r) between the variable and component. Together these coefficients make up the matrix of components. This matrix helps us to understand the various processes involved.

All statistical processing of the data collected was carried out using the statistics software v. 17 from SPSS®.

Results and discussion

Occurrence of TTHM

TTHM in water

As previously observed in such seawater swimming pools [START_REF] Parinet | Exposure levels to brominated compounds in seawater swimming pools treated with chlorine[END_REF], CHBr 3 represent the THM with the highest amount found (300.0 µg.L -1 mean, 1029.0 µg.L -1 maximum), about 95% of the sum of TTHM, followed by CHClBr 2 , 4 to 5% of the sum of THM-4 (18.9 µg.L -1 mean, 81.0 maximum) and CHCl 3 (0.1 presence of brominated compounds as compared to chlorinated compounds results from the rapid formation of hypobromous acid (HBrO) and hypobromite based disinfectant are added into bromide halogenating agent, bromine is usually more reactive than chlorine especially with compounds, brought by swimmers or bathers, (Manasfi et al., 2017b). -E3.A (Figure 1B) with TTHM (955.8 µg.m -3 mean, 1780.0 µg.m -3 maximum), CHCl 3 (7.5 µg.m -3 mean, 29.0 µg.m -3 maximum), CHBrCl 2 (1.6 µg.m -3 mean, 2.6 µg.m -3 maximum), CHClBr 2 (48.9 µg.m -3 mean, 150.0 µg.m -3 maximum) and CHBr 3 (897.7 µg.m -3 mean, 1600.0 µg.m -3 maximum). CHBr 3 levels is 3 to 10-fold greater than the others pools.

-E4.C (Figure 1B) with TTHM (69.5 µg.m -3 mean, 207.0 µg.m -3 maximum), CHCl 3

(5.0 µg.m -3 mean, 15.0 µg.m -3 maximum), CHBrCl 2 (3.6 µg.m -3 mean, 18 µg.m -3 maximum), CHClBr 2 (2.8 µg.m -3 mean, 12.0 µg.m -3 maximum) and CHBr 3 (58.0 µg.m -3 mean, 200.0 µg.m -3 maximum).

We interpreted the discrepancy of the establishment 3 by a low level of daily renewal of air and the low levels of TTHM in the air of the pool E4.C is due to the chlorination once a week.

Whereas in swimming pools supplied with tap water the most commonly detected THM in the air of swimming pools is chloroform -with averaged levels typically ranging from 3.1 to 320 µg.m -3 (Table 1), the predominant molecule found during the present study was bromoform (Figure 1B) with very high levels ranging from 80.0 µg.m -3 (pool E2.A) to 1380.0 µg.m -3 (pool E3.A) (964.7 µg.m -3 mean). Bromoform represented between about 95% of TTHM with a mean of 266.1 µg.m -3 , followed by dibromochloromethane (CHBr 2 Cl) representing about 4-5 % of the sum of THM-4.

Finally, dichlorobromomethane (CHCl 2 Br) and chloroform (CHCl 3 ) were present at low levels, representing less of 1 % of TTHM. As observed in waters, levels of CHBr 3 are higher at the closing of pools, and at higher concentrations in summer than in winter. However, levels of CHBr 3 in air are poorly correlated with levels found in water. These trends have also been observed in other studies (Aggazzotti et al., 1995;[START_REF] Bessonneau | Determinants of chlorination by-products in indoor swimming pools[END_REF].

Volatilization of THMs could be affected by several factors including water temperature, air temperature, concentration in water, intensity of air circulation (ventilation), and the turbulence in water caused by vigorous activity of pool users (Hsu et al., 2009). In this study, the high THM level in the air in basin E3.A seems related to poor ventilation in this facility (partially defective ventilation system as described by the pool operator).

Seasonal variation

Two temporal variability factors also affect the levels of TTHM found in waters (and air). The levels of CHBr 3 found in water (and air) are clearly correlated to time of sampling, with high increase of these levels between opening and closing of the pools (Figure 2), except for pool E4.A, most likely because this basin is not heated and consequently not very frequented. Moreover, levels of CHBr 3 are also strongly related to the season, their levels being found to be 1 to 40-fold higher in summer than in winter. This seasonal influence has already been demonstrated in freshwater chlorinated swimming pools [START_REF] Bessonneau | Determinants of chlorination by-products in indoor swimming pools[END_REF], where levels of chloroform and chloramines were systematically higher in summer than in winter.

This trends is often linked to the highest frequentation of swimming pools in hot season, but this cannot be the only explanation when considering thalassotherapy. Indeed, in such establishments, frequentation is the same in winter than in summer. Input of other compounds than body fluids must thus be raised in summer and, use of sunscreens combined to the absence of pre-swim showering may explained highest levels of DBPs (bromoform) found in water. Indeed, several constituents of sunscreens have been proven to be reactive with chlorine [START_REF] Duirk | Reaction of Benzophenone UV Filters in the Presence of Aqueous Chlorine: Kinetics and Chloroform Formation[END_REF][START_REF] Manasfi | Degradation products of benzophenone-3 in chlorinated seawater swimming pools[END_REF]. This relationships is also found when comparing activities carried out in these pools, whatever the season. CHBr 3 levels are indeed higher in pools where aquatic activities (aquagym) are offered daily to bathers (E1B vs E1.A; E2B. vs E2.A). Several past studies have already shown that heavy use of pools lead to a highest release of body fluids such as sweat that act as precursors of DBPs [START_REF] Judd | The fate of chlorine and organic materials in swimming pools[END_REF][START_REF] Keuten | Quantification of continual anthropogenic pollutants released in swimming pools[END_REF][START_REF] Afifi | Seasonal dynamics of water and air chemistry in an indoor chlorinated swimming pool[END_REF].

However, pools E4.B and E4.C do not show the same pattern of evolution : levels of CHBr 3 in the "relaxing" basin are higher than in the "sporting" basin. A possible explanation is the use of botanicals and essential oils when curists (spa patients) are massaged in thalassotherapy centres. The reactivity of these compounds with chlorine has not been yet reported in the literature, but because of their aromatic structure they are supposed to be highly reactive with chlorine (or bromine) and to be very efficient precursors of DBPs.

General study on TTHM in air and water using PCA

The objective of the PCA approach is to identify, in real environments, the water characteristics and the operating conditions that may be involved in the appearance of DBPs and to help better understand the origins of these molecules.

In general, PCAs only produce graphic interpretation resulting from the extraction of the first two components. However, for most aquatic environments, the first four or five components provide other information on the physical-chemical, biological and even climatic processes that govern the environment, for which the first two components are not always enough to explain [START_REF] Legendre | Numerical Ecology[END_REF][START_REF] Shirodkar | Factors controlling physico-chemical characteristics in the coastal waters off Mangalore-A multivariate approach[END_REF].

In this study, to allow each component to be represented by robust variables, a varimax rotation of different varifactors with factor loading was calculated using eigenvalues greater than 0.3. PCA was applied to normalized data, and so the covariance matrix coincides with the correlation matrix. To examine the suitability of these data for factor analysis, Kaiser-Meyer-Olkin (KMO) and Bartlett's tests were performed. KMO is a measure of sampling adequacy that indicates the proportion of variance which is common variance, i.e. which might be caused by underlying factors. High KMO values (close to 1) generally indicate that factor analysis may be useful, which is the case in this study: KMO = 0.567; if the KMO test value is less than 0.5, factor analysis will not be useful. Bartlett's test of sphericity indicates whether the correlation matrix is an identity matrix, which would indicate that variables are unrelated. The significance level is 0 in this study; less than 0.05 indicates that there are significant relationships among variables.

The PCA of the 19 parameters and rotated component matrix obtained for the varimax factor analysis indicate that altogether five components (Kaiser normalization) explain 86% of the total variance. The Kaiser normalization is a common strategy to select the "right" number of components by selecting components with an eigenvalue > 1. This value is significant given number of variables taken into account at the outset for the PCA. The interpretation of the processes described for each component can be carried out based on the variables that are used in their construction, we have retained uniquely the parameters in bold in the rotated matrix component (Table 2): -Component 1, which is built upon the variables: THM water , CHBr 3water , CHBr 2 Cl wate r, NPOC, Cl f , UV254 (strong positive loadings). This component may explain the presence of the major Br-DBPs in water in relation with precursors like: TOC, UV254 and Cl f . The presence of Br-DBPs is positively correlated to the aromaticity (UV254), carbonaceous organic matter concentration (NPOC) and chlorine (Cl f ). In this study, bromine seems to be not correlated to the component 1. This is probably explained by the fact that bromine is in large excess (65 mg.L -1 as compared to the µg.L -1 levels of halogenated compounds). -Component 2 is based on the following variables: CHBrCl 2water , CHCl 3water (strong negative loadings) and conductivity and Br -(strong positive loading). It is a component that may explain the minor chlorinated DBPs concentrations found in the seawater swimming pools are inversely correlated to bromide (Br -) and consequently to conductivity.

-Component 3, which is comprised of the variables THM air , CHBr 2 Cl ai r, CHBr 3air (strong positive loadings). This component may explain the major Br-DBPs concentrations found in the air of seawaters swimming pools.

-Component 4, which is comprised of the variables CHBrCl 2air , CHCl 3air and hygrometry (strong positive loadings) and pH (moderate negative loading). This component may explain the minor chlorinated DBPs concentrations found in the air of seawaters swimming pools.

These DBPs are correlated to hygrometry and inversely correlated to pH in water.

-Component 5, where we noted the presence of T water and T air is the component of temperature. Temperature of the air seems to be independent to DBPs formation, that is not the case for the temperature of the water which is positively correlated to component 1 (loading: 0.425) and component 2 (loading: 0.403). Regarding the establishment E3 and its swimming pool E3.A, we observe that the plots of this swimming pool are scattered along the component 3, so there is a wide variation in concentrations of DBPs in the air of this swimming pool. The same observations can be made for seawater swimming pools E4.A and E4.B, which are scattered along the component 3 too.

Based on this statistical study, we can conclude that PCA allows to distinguish pools compared to DBPs concentrations found in air and water and help to establish where the pollutions are highest, and how they vary. Moreover, PCA allows to observe the impact of ventilation, and the influence of the raw water quality, size of the swimming pools and attendance.

NPOC, Cl f and UV254 are found to be the precursors of THM found in water, particularly to the major ones (CHBr and CHBr 2 Cl water ). It is not the case for the minor ones (CHCl 3water and CHBrCl 2water ), that are only correlated to NPOC. Bromide concentrations are not linked to the formation of Br-DBPs probably due to their excess in the seawater swimming pools (NPOC -organic loads brought by bathers-act as the limiting reagent). TTHM concentrations found in the hall air of seawater swimming pools are presumably linked by a non-linear correlation to their concentrations in waters. This conclusion has been done by the fact that they are on four different components (component 1 and 2 for THM in water and component 3 and 4 for THM in air). Perhaps this is due to the operating conditions, particularly the air renewal.

Temperature of water seems to be positively correlated to the formation of Br-DBPs in water (the majors) and negatively to the chlorinated DBPs (the minors).

Conclusions

This study shows that in chlorinated seawater pools the speciation of DBPs in the water and air is predominantly brominated. Bromoform is the most abundant species (among THM) found in air and water of seawater swimming pools. These high concentrations related to the suspected carcinogenic effects of THM-Br are high of concern in seawater swimming pools treated with chlorine-based disinfectants. The THM-Br are genotoxic, mutagenic and are activated to mutagens by glutathion-transferase-theta-1 (GSTT1-1) [START_REF] Demarini | Glutathione S-transferase-mediated induction of GC-->AT transitions by halomethanes in Salmonella[END_REF][START_REF] Kogevinas | Genotoxic effects in swimmers exposed to disinfection by-products in indoor swimming pools[END_REF]) and a recent study found higher bladder cancer risks for people expressing GSTT1-1 [START_REF] Cantor | Polymorphisms in GSTT1, GSTZ1, AND CYP2E1, disinfection by-products, and risk of bladder cancer in Spain[END_REF]. Bromoform induced sister chromatid exchanges and micronuclei and is considered as a probable human carcinogen by US-EPA (IRIS, 1999;ATDSR, 1990). Dibromochloromethane induced chromosomal aberrations and sister chromatid exchanges [START_REF] Kogevinas | Genotoxic effects in swimmers exposed to disinfection by-products in indoor swimming pools[END_REF] and is considered to be a probable human carcinogen by US-EPA (ATDSR, 1990).

This study also emphasizes the need to not only improve water treatment systems or to decrease hydraulic retention time of water in pools (rate of renewal of water, rate of filtration rate), but also to improve air treatment systems, as it has already been demonstrated that there is a strong correlation between levels of Br-DBP found in ambient air and in exhaled breath of pool attendants [START_REF] Fantuzzi | Occupational exposure to trihalomethanes in indoor swimming pools[END_REF]. Moreover, the levels of inhaled Br-DBPs are also affected by physical activities of bathers in these pools [START_REF] Font-Ribera | Environmental and personal determinants of the uptake of disinfection by-products during swimming[END_REF]. Whatever these improvements in air or water treatments, nothing can replace bathers' hygienic behaviour, consisting in a 60-s pre-swim shower with soapy water, and by wearing clean bathing suits and swimming caps [START_REF] Keuten | Definition and quantification of initial anthropogenic pollutant release in swimming pools[END_REF].

Figure 1A :

 1A Figure 1A : THMs concentrations (logarithmic scale) in the waters of the different swimming pools studied

Figure 1B : 0 :M

 1B0 Figure 1B : THMs concentrations (logarithmic scale) in the airs of the different swimming pools studied

Figure 2 :

 2 Figure 2 : Temporal variations of CHBr 3 in water and air seawater swimming pools under study

Figure 3

 3 Figure 3 represents the score plots of the seawater swimming pools for different campaigns. The component 1 is the component of the major DBPs found in waters and component 3 is the component of the major DBPs found in airs. As we can see on this figure, the swimming

Figure

  Figure

Table 1 :

 1 Classical THM levels found in indoor air swimming pools

	References

Acknowledgements

This research was supported by a grant from the French Agency for Food, Environmental and Occupational Health and Safety (Grant 2009-CRD-21). Tarek Manasfi acknowledges the Doctoral School of "Environmental Sciences" (ED251) at Aix-Marseille University and the French Ministry of Higher Education and Research for the doctoral scholarship.