Mojdeh Rastgoo

Cedric Demonceaux

Ralph Seulin

Olivier Morel

Attitude estimation from polarimetric cameras

In the robotic field, navigation and path planning applications benefit from a wide range of visual systems (e.g. perspective cameras, depth cameras, catadioptric cameras, etc.). In outdoor conditions, these systems capture information in which sky regions cover a major segment of the images acquired. However, sky regions are discarded and are not considered as visual cue in vision applications. In this paper, we propose to estimate attitude of Unmanned Aerial Vehicle (UAV) from sky information using a polarimetric camera. Theoretically, we provide a framework estimating the attitude from the skylight polarized patterns. We showcase this formulation on both simulated and real-word data sets which proved the benefit of using polarimetric sensors along with other visual sensors in robotic applications.

I. INTRODUCTION

Large-field cameras and lenses (e.g. omnidirectional and fisheye cameras) are popular in robotic applications due to their ability to provide large field of view (up to 360 • ), extending the amount of visual information. It is the main reason for which they have been adopted for a broad range of tasks such as visual odometry [START_REF] Scaramuzza | Appearance-guided monocular omnidirectional visual odometry for outdoor ground vehicles[END_REF], navigation [START_REF] Winters | Omnidirectional vision for robot navigation[END_REF], simultaneous localization and mapping (SLAM) [START_REF] Kim | Slam with omni-directional stereo vision sensor[END_REF], and tracking [START_REF] Kobilarov | People tracking and following with mobile robot using an omnidirectional camera and a laser[END_REF]. With those systems, sky regions in the images acquired represent a large segment of information which are usually discarded. Here, we show that polarimetric information can be extracted from those regions and used in robotic applications.

Sun position, stars and sky patterns are hold as navigational cues for the past centuries. Indeed, before the discovery of magnetic compass, these natural cues have been the solitary source of navigation used by our ancestors [START_REF] Barta | Psychophysical study of the visual sun location in pictures of cloudy and twilight skies inspired by viking navigation[END_REF][START_REF] Horváth | On the trail of vikings with polarized skylight: experimental study of the atmospheric optical prerequisites allowing polarimetric navigation by viking seafarers[END_REF]. Similarly, some insects used the skylight polarized pattern created by the scattered sunlight to navigate in their environment [START_REF] Wehner | Desert ant navigation: how miniature brains solve complex tasks[END_REF][START_REF] Labhart | Neural mechanisms in insect navigation: polarization compass and odometer[END_REF]. For instance, desert ants (cataglyphis), butterflies and dragonflies among others, are able to navigate through their paths, efficiently and robustly by using the polarized pattern of sky, despite their small brains [START_REF] Labhart | Neural mechanisms in insect navigation: polarization compass and odometer[END_REF][START_REF] Wehner | Desert ant navigation: how miniature brains solve complex tasks[END_REF][START_REF] Hamaoui | Polarized skylight navigation[END_REF].

Acknowledging the nature, numerous studies have been conducted on polarized skylight pattern [START_REF] Lambrinos | A mobile robot employing insect strategies for navigation[END_REF][START_REF] Chu | Application of a novel polarization sensor to mobile robot navigation[END_REF][START_REF] Zhao | A novel angle algorithm of polarization sensor for navigation[END_REF][START_REF] Wang | A novel autonomous real-time position method based on polarized light and geomagnetic field[END_REF][START_REF] Chahl | Integration and flight test of a biomimetic heading sensor[END_REF][START_REF] Ashkanazy | Bio-inspired absolute heading sensing based on atmospheric scattering[END_REF][START_REF] Stürzl | A fisheye camera system for polarisation detection on uavs[END_REF][START_REF] Ma | An evaluation of skylight polarization patterns for navigation[END_REF][START_REF] Miyazaki | Estimating sunlight polarization using fish-eye lens[END_REF][START_REF] Wang | A bionic camera-based polarization navigation sensor[END_REF][START_REF] Lu | Angle algorithm based on hough transform for imaging polarization navigation sensor[END_REF][START_REF] Hamaoui | Polarized skylight navigation[END_REF]. These studies are generally reported in the optic field. They focus on estimating the solar azimuth angle by creating a sort of compass. Estimating polarized patterns have been, however, a difficult and complex task. The primary studies report the use of several photodiodes [START_REF] Lambrinos | A mobile robot employing insect strategies for navigation[END_REF][START_REF] Chu | Application of a novel polarization sensor to mobile robot navigation[END_REF][START_REF] Zhao | A novel angle algorithm of polarization sensor for navigation[END_REF][START_REF] Wang | A novel autonomous real-time position method based on polarized light and geomagnetic field[END_REF][START_REF] Chahl | Integration and flight test of a biomimetic heading sensor[END_REF], or of multiple cameras [START_REF] Ashkanazy | Bio-inspired absolute heading sensing based on atmospheric scattering[END_REF][START_REF] Stürzl | A fisheye camera system for polarisation detection on uavs[END_REF][START_REF] Wang | A bionic camera-based polarization navigation sensor[END_REF] or manually rotating filters [START_REF] Ma | An evaluation of skylight polarization patterns for navigation[END_REF][START_REF] Miyazaki | Estimating sunlight polarization using fish-eye lens[END_REF][START_REF] Lu | Angle algorithm based on hough transform for imaging polarization navigation sensor[END_REF][START_REF] Hamaoui | Polarized skylight navigation[END_REF]. As a consequence of those troublesome setups, robotic applications are not benefiting from the advantages of polarized patterns, as 1 Université de Bourgogne Franche-Comté, Le2i VIBOT ERL-CNRS 6000, France, mojdeh.rastgoo@gmail.com attested by the lack of polarized sensors used in Unmanned Aerial Vehicle (UAV). However, the recent introduction of division-of-focal-plane (DoFP) micropolarizer cameras has offered an alternative solution [START_REF]Micropolarizer array for infrared imaging polarimetry[END_REF][START_REF] Nordin | Diffractive optical element for stokes vector measurement with a focal plane array[END_REF][START_REF] Millerd | Pixelated phase-mask dynamic interferometers[END_REF]. In such cameras a micropolarizer filter array, composed of a pixelated polarized filters oriented at different angles, is aligned with a detector array. Thus, linear polarization information are simultaneously acquired taking a single image. Here, we use a DoFP coupled with a fisheye lens to exploit the polarized information of sky region to estimate vehicle attitude.

In this paper, Sect. II presents the specificity of the camera used and the adaptation required for our robotic application. The remainder of the paper is organized as follows: Sect. III introduces the concepts of polarization by scattering, Rayleigh model and its relation with attitude estimation. Our formulation to estimate attitude is presented in Sect. IV. Experiments and implementation details are given in Sect. V, and finally discussions and conclusions are drawn in Sect. VI.

II. SETTING THE POLARIMETRIC CAMERA FOR ROBOTICS

In this work, visual information is captured using the IMPREX Bobcat GEV polarimetric camera which is a DoFP polarimetric camera. In a single shot, the camera captures four different linearly polarized measures by using a micropolarizer with pixelated polarized filter array as illustrated in Fig. 1. Hence, each acquired image is subdivided into four linearly polarized images I 0 , I 45 , I 135 , and I 90 . Subsequently, the polarized state of the incident light is computed from these images by means of the Stokes' parameters [START_REF] Goldstein | Polarized light[END_REF], referred as s 0 , s 1 , and s 2 in Eq. (1). In addition, the polarized parameters angle of polarization (AoP) and degree of linear polarization (DoPl), respectively referred α and ρ l in Eq. (1) are computed. ) ,

ρ l = s 2 2 + s 2 1 s 0 . (1)
The camera captures images with a resolution of 640 px × 460 px and each polarized image is reconstructed from the interlaced pixels. In addition, we used a 180 • fisheye lens for the experiment reported in this paper to benefit from a large field-of-view. An example of an acquired raw image, the extracted polarized images and the computed polarized information is shown in Fig. 2 The IMPREX Bobcat GEV camera is controlled using the eBus SDK provided by Pleora Technologies Inc [START_REF]eBUS SDK | Pleora Technologies Inc[END_REF].

To enable the interaction with other robotic devices, we implemented a ROS package publicly available 1 enabling the usage of roslaunch and rosrun. In addition, our package allows to store and stream the raw data as well as computing the Stokes' and polarized parameters.

III. POLARIZED CUES USED FOR ATTITUDE ESTIMATION

Polarized cues used for attitude estimation are based on three main concepts which are presented in this section: (i) the Rayleigh scattering model and its implications on the polarization by scattering, (ii) the polarization parameters in pixel frame and its relation to camera, and (iii) the connection between the polarized parameters in the pixel frame and the parameters used to estimate the vehicle attitude.

A. Rayleigh scattering model

The unpolarized sunlight passing through our atmosphere gets scattered by different particles within the atmosphere. Besides deviating the direction of a propagated wave, this transition also changes the polarization state of the incident 1 https://github.com/I2Cvb/pleora_polarcam. This package is derived from earlier work [START_REF] Iralab | ROS device driver for PhotonFocus cameras based on Pleoras eBUS Software Development Kit (SDK)[END_REF]. light which can be explained using the Rayleigh scattering model. Rayleigh scattering describes the scattering of light (or any electromagnetic waves) by particles much smaller than their transmission wavelength. Accordingly, it assumes that scattering particles of the atmosphere are homogeneous and smaller than the wavelength of the sunlight. Despite these assumptions, this model proved to be sufficient for describing skylight scattering and polarization patterns [START_REF] Pomozi | How the clear-sky angle of polarization pattern continues underneath clouds: full-sky measurements and implications for animal orientation[END_REF][START_REF] Horvath | Groundbased full-sky imaging polarimetry of rapidly changing skies and its use for polarimetric cloud detection[END_REF].

The Rayleigh model predicts that the unpolarized sunlight becomes linearly polarized after being scattered by the atmosphere. Based on this model, two main outcomes are drawn. On the one hand, DoPl is directly linked to the scattering angle γ according to:

ρ l = ρ lmax 1 -cos 2 (γ) 1 + cos 2 (γ) , (2)
where ρ lmax is a constant equal to 1 in theory but slightly less than 1 in practice due to some atmospheric disturbances [START_REF] Pomozi | How the clear-sky angle of polarization pattern continues underneath clouds: full-sky measurements and implications for animal orientation[END_REF]. The scattering angle γ is defined by the angle between the observed celestial vector c and the sun vector s as presented in Fig. 3. Note that DoPl is 0 in the sun direction and maximum when the scattering angle is π 2 [START_REF] Smith | The polarization of skylight: An example from nature[END_REF][START_REF] Miyazaki | Estimating sunlight polarization using fish-eye lens[END_REF]. On the other hand, the scattered light is considered to be polarized and orthogonal to the scattering plane. Consequently, the Angle of Polarization is directly related to the orientation of the scattering plane.

B. Polarization by scattering model in pixel frame

As presented Fig. 4, an image is considered as a collection of pixels and each pixel measures the polarization parameters of the light traveling along a ray associated with that pixel. The pixel frame P is defined accordingly with the ray which coincides with c. The camera calibration determines the relationship between pixels and these 3D rays.

Let consider one pixel of the image with its associated pixel frame P ( obc). Based on Rayleigh scattering, the electric field of incident light after scattering is perpendicular to the scattering plane that is defined by the observer, celestial point, and the sun. Accordingly, the normalized electric field vector E in the world frame is presented as the normalized cross product of s and c as shown in Eq. ( 3).

E = s ∧ c s ∧ c , (3)
The same measurement in the pixel frame P is represented as:

E obc =   E o E b 0   =   cos α sin α 0   , (4)
where α is the measured AoP associated to the corresponding pixel. Combining Eq. ( 3) & (4) and using the scattering angle γ, between s and c lead to:

(s ∧ c) • o = sin γ cos α (s ∧ c) • b = sin γ sin α . ( 5
)
Applying the vector triplet cross product rule on Eq. ( 5) results in:

s • b = sin γ cos α s • o = -sin γ sin α . (6)
Using Eq. ( 2), the scattering angle γ is expressed as:

cos γ = s • c = ± 1 -ρ l 1 + ρ l , (7)
with ρ l = ρ l ρ lmax .

Equations ( 6) & ( 7) finally lead to a representation of the sun vector in pixel frame P which express a direct relation between the AoP, the scattering angle, and the sun position:

s p =   -sin γ sin α sin γ cos α cos γ   . (8)
In other words, the sun vector is expressed in the pixel frame as a vector depending only on the polarization parameters AoP and DoPl, which is directly linked to the scattering angle γ.

C. UAV attitude and polarized sky pattern

Figure 5 illustrates all transformations and frame conventions to estimate the attitude of a UAV. In addition, an inertial measurement unit (IMU) frame is added for later comparisons. In Fig. 5, W, W , I, V, C, and P refer to the world frame, global frame of IMU, IMU frame, vehicle frame, camera frame, and pixel frame, respectively. The rotation from one frame to another is presented with lowercase alphabet. In this scenario, a vector v p in pixel frame is expressed in the world frame v w as:

W {EN U } V {F LU } I {RLD} C {RDF } R ww R wv R vi R vc W R w i
v w = R wv • R vc • R cp • v p , (9)
where the rotation from the camera to the pixel frame R cp is obtained by camera calibration. The rotation is defined as the yaw and pitch rotations by the zenith and azimuth angles of the celestial point (θ c , φ c ) as shown in Eq. [START_REF] Horvath | Groundbased full-sky imaging polarimetry of rapidly changing skies and its use for polarimetric cloud detection[END_REF].

R cp =   cos θ c cos φ c -sin φ c sin θ c cos φ c cos θ c sin φ c cos φ c sin θ c sin φ c -sin θ c 0 cos θ c   , = R zc (φ c ) • R yc (θ c ) . (10)
In Eq. ( 8), we expressed the sun position in the pixel frame. Indeed, this representation can be applied to any point from the world frame, ergo:

s w = R wv • R vc • R cp •   -sin γ sin α sin γ cos α cos γ   , = R wv • R vc • v , (11)
thus

R T • s w = v . ( 12
)
The above equation shows a direct relation between the rotation matrix of the vehicle R, the AoP α -measured by the polarimetric camera at one pixel -and the angle of scattering γ for the corresponding pixel. In addition, if ρ lmax is known, the angle of scattering is directly obtained by inverting Eq. ( 7) providing a direct relation between polarization parameters and the rotation of the vehicle.

IV. ATTITUDE ESTIMATION

In this section, we present 2 approaches to estimate the attitude: (i) in the former method named absolute rotation, the absolute rotation and attitude of the vehicle is estimated under the assumption that the sun position is known or deduced from time and the GPS location of the vehicle and (ii) in the latter method named relative rotation, the relative rotation of the vehicle from its initial position is estimated without additional assumption regarding the sun position. Both approaches required an estimate of the scattering angle γ which will be presented beforehand.

A. γ estimation

By only measuring the AoP α in scattering effects, γ needs to be estimated to get the vector v defined in Eq. ( 11). This equation is valid for all points in sky region. However, only 2 celestial points are required to estimate γ such as:

                   R t • s = R cp1 •    -sin γ 1 sin α 1 sin γ 1 cos α 1 cos γ 1    R t • s = R cp2 •    -sin γ 2 sin α 2 sin γ 2 cos α 2 cos γ 2    . ( 13
)
Using the product of R cp and R z (α), Eq. ( 13) is rewritten as:

M 1 •   0 sin γ 1 cos 1   = M 2 •   0 sin γ 2 cos γ 2   . ( 14
)
By defining the matrix M such that M = M t 2 • M 1 , γ 1 and γ 2 are found as:

γ 1 = -arctan M02 M01 γ 2 = -arctan M20 M10 . ( 15
)
The AoP is 2π modulus, while the γ found in Eq. ( 15) is π modulus leading to two possible solutions for the vector v: (α 1 , γ 1 ) and (α 1 + π, -γ 1 ).

B. Absolute rotation

In order to estimate the absolute rotation and attitude of the UAV, it is assumed that: (i) the sun position is known (ii) the vector v is estimated using the AoP measures of the sky (2 points) and (iii) the vertical in the pixel frame is known or a second w is estimated using the AoP from horizontal reflected areas. In this study, the vertical in the pixel frame is assumed to be known.

The aforementioned assumptions lead to the following expression:

[s, z, s ∧ z] = R(t) • [v(t), w(t), v(t) ∧ w(t)] = R wv (t) • R vc • [v(t), w(t), v(t) ∧ w(t)] . ( 16
)
where z is the vertical in world frame ([0, 0, 1]) and t is the time instance.

Solving Eq. ( 16) enables to get R wv (t). However, due to γ ambiguities, v and therefore R wv have 2 solutions. At each iteration, the rotation R wv selected is the one the closest from previous rotation, assuming that the motion between two frames is smoothed.

C. Relative rotation

The relative rotation is estimated between two time stamps (t 1 , t 2 ). Let v(t 1 ) and v(t 2 ) referring to v 1 and v 2 to simplify the expression. Therefore, Eq. ( 16) becomes:

[s, z, s ∧ z] = R wv1 • R vc • [v 1 , w 1 , v 1 ∧ w 1 ] [s, z, s ∧ z] = R wv2 • R vc • [v 2 , w 2 , v 2 ∧ w 2 ] . (17)
Leading to:

R wv2 = R wv1 • R vc • [v 1 , w 1 , v 1 ∧ w 1 ] • [v 2 , w 2 , v 2 ∧ w 2 ] -1 • R T vc , (18)
Using the above equation, the relative rotation R v1v2 is equal to:

R v1v2 = R vc • [v 1 , w 1 , v 1 ∧ w 1 ] • [v 2 , w 2 , v 2 ∧ w 2 ] -1 • R t vc . (19)
As previously explained, only 2 points are required to compute the scattering angle γ. In practice (see Sect. V), more points can be used which leads to a more robust estimate of the vehicle rotation. Fig. 8. Absolute and relative rotation obtained from synthetic data in ideal case, noisy conditions without RANSAC optimization and noisy case with RANSAC optimization. The mean and standard deviation of the difference between each predicted and GT is illustrated in the legend as well.

V. EXPERIMENTS AND RESULTS

This section presents our experimental setup, the designed experiments and the obtained results. The setup used in our experiment is illustrated in Fig. 5. Instead of using a UAV, a camera was manually moved (see Fig. 6). The IMU was calibrated with the camera using kalibr toolbox [START_REF] Furgale | Unified temporal and spatial calibration for multi-sensor systems[END_REF][START_REF] Furgale | Continuous-time batch estimation using temporal basis functions[END_REF], to use the recordings as GT. The polarimetric camera with fisheye lens was also calibrated according to [START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF]. Using the above setup two data sets of synthetic and real images were created and the results obtained are presented in Exp. V-A and Exp. V-B, respectively.

A. Experiment 1

The synthetic data containing AoP and DoPl images of sky regions were created using the IMU recordings obtained during real acquisition. Figures 7(a Applying our framework on ideal synthetic data, perfect results were obtained for absolute and relative rotations, while γ was estimated using only 2 random points from the sky region (blue dotted curve in Fig. 8).

Although using our proposed framework, we were able to achieve perfect results on ideal synthetic data, in reality it is rare to obtain the perfect skylight polarization pattern. Variety of causes clutter the desired skylight pattern, the main one being pollution. To account for such cases, a second test was performed while significant level of noise was added to the created synthetic data. Figures 7(c Performing the same 2 random-points algorithm as before on noisy dataset leads to the results illustrated in Fig. 8, the orange curve. As expected the performance decline, simply due to the noise.

To solve this problem, Random sample consensu (RANSAC) was used to estimate the attitude in the absolute and relative rotation methods. For the absolute rotation, since the sun position is assumed to be known the model optimizes the full rotation of each frame in comparison with the origin considering the difference between the predicted and real sun positions.

However, in the relative rotation method there is no information about the original position, or sun position, and the algorithm only depends on the polarized vector, v = R cp .v p between two different frames. Therefore, using the model infers the optimal vector representing each frame is obtained.

Estimating the attitude on the noisy dataset with RANSAC significantly improved the results as shown in Fig. 8 (greendashed line). The parameters used were: an error threshold of 0.07, 10 random points (2 points for defining the model and the rest as test), and 2000 iterations.

The quantitative results in terms of mean difference (µ) and standard deviation (σ) between the predicted rotations and GT, for all the conditions are also shown in the Fig. 8. Note that all the angles are reported in radians.

As illustrated in the obtained results, using the RANSAC model the outliers are ignored and satisfactory results are achieved.

B. Experiment 2

This section presents the results obtained using real data. The same experimental setup was used using the IMU results to create the GT for the vehicle in the world frame. The original data set contains 593 recordings which was undersampled to a sampling rate of 20 frames. Both absolute and relative methods were ran with RANSAC with the same parameters than in the previous experiment. The results are presented in Fig. 9.

Even though the difference between the predicted rotation and GT using the real data is higher than for the synthetic measurements the results are promising and the pose trajectory of the vehicle is respected.

VI. DISCUSSION AND CONCLUSION

This paper presented a new method to estimate attitude of a vehicle using the polarization pattern of the sky. Contrary to conventional cameras, polarimetric cameras exploit part of images that represent the sky and we demonstrated how to take advantage of this property to estimate attitude. We first derived all equations that describe the relationship between the rotation matrix of the vehicle and the polarization parameters. Herein, we proposed a model based on AoP measurements of the light beam scattered by the sky, subsequently two approaches of the absolute attitude and relative attitude estimation were proposed. The former estimated the rotation the vehicle in comparison to the origin taking into account the sun position while the latter did not consider this assumption and estimated the position of the vehicle in the world frame considering two consecutive frames. Finally, in order to cope with the undesired artifacts and outliers that can occur during the measurements, a RANSAC model was integrated within our framework. Promising results were achieved after using RANSAC optimization, illustrating the potential and capacity of a polarimetric camera to be integrated in the robotic field. As future work, we will focus our attention to improve these preliminary results including a minimization process of the accumulated error of prediction using filtering. Aware that the polarimetric camera cannot be a standalone system for robust attitude estimation, we also plan to combine this modality with geometric information to improve the quality of the estimation.

VII. ACKNOWLEDGMENT

This work is part from a project entitled VIPeR (Polarimetric Vision Applied to Robotics Nav-igation) funded by the French National Research Agency ANR-15-CE22-0009-VIPeR.

Fig. 1 .

 1 Fig. 1. Structure of DoFP sensors: in a single shot, four polarized images are acquired, each of them with a different polarized angles.

Fig. 2 .

 2 Fig. 2. Polarimetric images: left: a raw images in which the different polarimetric images are interlaced; center: the four extracted linearly polarized images (I 0 , I 45 , I 135 , I 90 ); right: the AoP and DoPl images. For visualization purpose, AoP is represented in HSV colorspace.

  s 0 = I 0 + I 45 + I 135 + I 90 4 , s 1 = I 0 -I 90 , s 2 = I 45 -I 135 , α = 0.5 arctan( s 2 s 1

Fig. 3 .Fig. 4 .

 34 Fig. 3. Skylight polarization by scattering. Scattering plane is highlighted by light shade of red. (θs, φs) and (θc, φc) define the zenith and azimuth angle of sun and celestial point respectively. obc defines the pixel frame, P, and E is the electrical field orthogonal to the scattering plane.

Fig. 5 .

 5 Fig. 5. Frame conventions and rotations for attitude estimation of an UAV.

Fig. 6 .

 6 Fig. 6. Experimental setup.

Fig. 7 .

 7 Fig. 7. AoP and DoPl images synthetically created. All images have been generated with sky region for yaw, pitch and roll angle of 1.8 rad, -0.2 rad and 0.1 rad, respectively. (a)-(b) No noise added, (c)-(d) Noise level of 0.1 added.

  ) & 7(b) show an example of this dataset at optimal conditions. This dataset based on IMU recordings contains rotations along roll, pitch, and yaw. This dataset has originally 856 samples but has been downsampled by sampling rate of 30 samples.

  ) & 7(d) show an example of synthetic data with an additional Gaussian noise with an

Fig. 9 .

 9 Fig. 9. Absolute and relative rotation obtained from real data. The black line represents the GT, the orange line and the blue line represent the absolute and relative predicted rotations, respectively.


Imeen Ben Salah

Sébastien Kramm

Cédric Demonceaux

Pascal Vasseur

Summarizing Large Scale 3D Mesh

Recent progress in 3D sensor devices and in semantic mapping allows to build very rich HD 3D maps very useful for autonomous navigation and localization. However, these maps are particularly huge and require important memory capabilities as well computational resources. In this paper, we propose a new method for summarizing a 3D map (Mesh) as a set of compact spheres in order to facilitate its use by systems with limited resources (smartphones, robots, UAVs, ...). This vision-based summarizing process is applied in a fully automatic way using jointly photometric, geometric and semantic information of the studied environment. The main contribution of this research is to provide a very compact map that maximizes the significance of its content while maintaining the full visibility of the environment. Experimental results in summarizing large-scale 3D map demonstrate the feasibility of our approach and evaluate the performance of the algorithm.

I. INTRODUCTION

Environment mapping is a challenging and important question that has been widely discussed in mobile robotics research. Most of the early vision-based navigation systems need a map that contains a sufficient level of detail for an accurate navigation. Over the last few years, the introduction of High-Definition (HD) [START_REF] Heiko | Autonomous driving in the icityhd maps as a key challenge of the automotive industry[END_REF] and semantic maps have greatly participated in the large commercial success of navigation and mapping products and also in the enhancement of data fusion based localization algorithms. Several digital map suppliers like TomTom and HERE are now providing HD maps with very high navigation accuracy, especially in challenging urban environments. These HD maps provide a detailed representation of the environment thus they require a high processing capacity with severe time constraints as well as a large storage requirements. In order to use this type of maps by navigation systems with limited resources (computation / memory), the size of these maps must be reduced. Hence, we need to summarize such maps while preserving the essential information for navigation.

II. PREVIOUS WORKS

High quality maps are required for a wide range of applications, particularly for applications where constant localization or feature tracking is crucial. Hence understanding a maps structure and evaluating map quality is very useful for accurate navigation [START_REF] Merzić | Map quality evaluation for visual localization[END_REF], [START_REF] Zhen | Robust localization and localizability estimation with a rotating laser scanner[END_REF]. Furthermore, in some navigation tasks, setting a full-size map on a mobile device (car, robot, etc.) poses several difficulties in time and memory consumption. To solve this problem there are two possible approaches: map compression or map summarizing using only significant information. Over the past decade, several compression approaches have been proposed in literature. Some of them employ special data structures such as octree [START_REF] Devillers | Geometric compression for interactive transmission[END_REF] for progressive encoding of point clouds. Schnabel et al. [START_REF] Schnabel | Octree-based point-cloud compression[END_REF] proposed a prediction scheme to achieve the compression using the octree approach. A novel compression method has been proposed in [START_REF] Kammerl | Real-time compression of point cloud streams[END_REF] to code only the spatial and temporal differences within an octree structure. Jae-Kyun et al. [START_REF] Ahn | Large-scale 3d point cloud compression using adaptive radial distance prediction in hybrid coordinate domains[END_REF] proposed a geometry compression algorithm for large scale point cloud to encode radial distances in a range image. Compression algorithms are essential for efficient storage and transmission but they are not sufficient for navigation indeed a step of decompression is necessary to process the data. Several approaches based on feature selection to summarize a map for localization purposes were presented in [START_REF] Dymczyk | The gist of maps-summarizing experience for lifelong localization[END_REF], [START_REF] Mühlfellner | Summary maps for lifelong visual localization[END_REF], [START_REF] Soo | 3d point cloud reduction using mixed-integer quadratic programming[END_REF]. They propose some scoring functions to order the map landmarks according to their significance while guaranteeing a maximum coverage of the scene. An approach to map reduction was proposed in [START_REF] Steiner | Location utilitybased map reduction[END_REF]. It aims to select only the places that are particularly suitable for localization using the location utility metric. To simplify the process of appearance-based navigation, a selection process is applied to choose the key/reference features in the environment. For instance, in visual memory based approaches, a set of relevant and distinctive areas (images) are acquired and used during navigation for comparison with the current position. In the work of Cobzas [START_REF] Cobzas | Image-based localization with depth-enhanced image map[END_REF], a panoramic memory of images is created by combining acquired images with depth information extracted from a laser scanner. In this image database, only the salient information will be retained [START_REF] Benhimane | Linear and quadratic subsets for template-based tracking[END_REF] without degrading the performance during navigation. In order to build this image database, some techniques have been developed to guarantee the maximum efficiency in the choice of useful information. A spherical representation has been proposed by M. Meilland et al. [START_REF] Meilland | Dense omnidirectional rgb-d mapping of large-scale outdoor environments for real-time localization and autonomous navigation[END_REF]. This spherical representation is built by merging different images acquired by a set of cameras with the depth information extracted from a laser scanner. Methods based on Bag of Words (BoW) are widely used for localization. BoW methods can efficiently represent a huge amount of data using the occurrences of several visual vocabulary. By applying hierarchical dictionary to the visual navigation problem [START_REF] Fraundorfer | Topological mapping, localization and navigation using image collections[END_REF], BoW methods proved a high scalability and accuracy in vision-based localization and mapping processes. All existing methods for summarizing maps are based mainly on geometric or photometric characteristics to select the most salient information. However, these characteristics are insufficient for good perception and understanding of the environment. A combination of geometric, photometric and semantic characteristics when selecting information allows us to have a compact, precise and useful summary. Our work aims to perform several navigation tasks using only a map summary of the environment. This map should be not only compact but also coherent with the perception of the agent. To provide this map summary, we propose a new method dealing with largescale 3D meshs. The output of our summarizing method is a set of spherical images. the spherical images are constructed using the technique of ray tracing. For each pixel of the image, a ray is launched and the intersections of this ray with the faces of the mesh are calculated. Ray tracing is a rendering technique for generating an image by tracing the path of light as pixels in an image plane and simulating the effects of its encounters with virtual objects. In our case, for each pixel, we launch a ray, and we look for the triangles of the mesh intersected by the ray. We keep the nearest triangle. We have proposed to use the Mller-Trumbore [START_REF] Möller | Fast, minimum storage ray/triangle intersection[END_REF] algorithm to determine if a ray intersects a given triangle or not while keeping track of the nearest intersected object. Our main contributions are:

• We propose an efficient algorithm for summarizing maps based on geometric, photometric and semantic characteristics.

• We formulate map summarizing process as a multiobjective optimization problem. • We propose a new scalable representation of summarized map using augmented and labelled RGB-D-L spherical images. Where RGB is the color information, D is the depth information and L is the semantic information The remainder of this paper is organized as follows. Section III is dedicated to the problem formulation and presents an overview of our approach. The complete map summarization process is described in Sections IV and V. Before conclusion and perspectives, experimental results and discussions are presented in Section VI.

III. PROBLEM MODELING

In this article, we aim at summarizing a 3D mesh to a compact set of spherical images. Our objective is to reduce the number of spheres while maintening a high level of significant information for navigation and localization. To modelize this problem, we introduce some notations. Let us suppose that the 3D map is a 3D mesh with color and semantic information. In this 3D mesh we define the navigable areas Γ as the 2D polygon whose edges represent the borders of the navigable areas. We define the visible areas Σ as the 2D polygon whose edges represent the borders of the visible areas. To reduce the data, we propose to represent the different information into some augmented spherical images. The goal of this paper is to minimize the number 𝑛 of augmented spheres 𝑋 = {𝑋 𝑖 , 𝑖 = 1..𝑛} needed to describe the 3D maps. The final set of spheres must guarantee a compromise between two objectives:

• Maximize the Visibility • Maximize the Entropy .

The first goal is to maximize the visibility 𝑉 of the studied environment. We define the visibility of a 3D point 𝑐 in the navigable areas Γ as the number of all vertices of Σ visible from this point. Geometrically, the point 𝑝 is said to be visible from 𝑐 if the segment 𝑝𝑐 does not intersect any obstacle or any edge of Σ. We define the visibility 𝑉 as:

𝑉 (𝑐) = #{𝑝 ∈ vertices of Σ/𝑝 is visible from 𝑐} (1)
where #𝐴 is the cardinality of the set 𝐴. The second objective is the entropy 𝐸 which expresses the quantity of visible salient information from a viewpoint. We define the entropy of a sphere as the amount of significant information projected onto it. Using this entropy we are able to select the optimal viewpoint which represents as best as possible a 3D map. The idea is to compute the number of salient points in the spherical image. To decide the saliency of 3D point on the sphere, we propose to define two levels of saliency. The low-level saliency that is based on lowlevel characteristics (photometric and geometric) of a 3D point. The high-level saliency is based on the semantic information of each 3D point. In our summarizing process, points labelled "buildings", are considered among the most salient points for localization. Using these two types of saliency, we calculate the number of points of interest on the sphere according to their relevance. Therefore, we have four possible combinations as following.

• 𝑛 00 : number of non relevant points semantically, photometrically and geometrically. • 𝑛 10 : number of points relevant only photometrically and geometrically.

• 𝑛 01 : number of points relevant only semantically.

• 𝑛 11 : number of points relevant semantically, photometrically and geometrically. The entropy of a sphere will be characterized by the entropy of its center 𝐶. The entropy is given by the equation:

𝐸(𝐶) = -𝑛00 ℎ log 𝑛00 ℎ -𝑛01 ℎ log 𝑛01 ℎ -𝑛10 ℎ log 𝑛10 ℎ -𝑛11 ℎ log 𝑛11 ℎ
where ℎ is the number of 3D points on the sphere. We then propose to formulate this summary process of 3D mesh as a multi-objective optimization problem using a simultaneous optimization strategy. Thus we propose to optimise the visibility and the entropy in the same time. Algorithm 1 finds the optimal set of spherical images 𝑋 = {𝑋 𝑖 , 𝑖 = 1..𝑛} which allows to describe efficiently a 3D mesh. Here we represent each sphere by its center 𝑋 𝑖 .

A. Barycenters graph

The first goal of our summarizing process is to reduce the size of a 3D map while maintaining a global view of the entire scene. Thus, the optimal set of spherical images must guarantee maximum visibility. All spheres must be placed on the navigable areas (street). To simplify the search space we will look for these spheres on a starting graph called Barycenters graph. On this graph 𝑇 = (𝐶, 𝛾) we will place our spherical images. The nodes of this graph will be a set of points of view allowing to visualize the totality of the scene. All these points of view will be placed on the navigable area. In the literature several methods have been proposed to treat visibility in a polygon. The art gallery problem was presented as a question by the German painter Paul Klee: how much guards do we need to guard a room? This problem consists in finding in a polygon the optimal number of guards required to view the totality of the polygon. A variety of this problem is known as the "Watchman route". This problem consists in finding a path inside a polygon Γ so that each point in Γ will be visible from a certain point of this path. Several algorithms have been proposed to calculate the watchman route by taking as input a group of target points to visit [START_REF] Carlsson | Finding the shortest watchman route in a simple polygon[END_REF], [START_REF] Dror | Touring a sequence of polygons[END_REF], [START_REF] Packer | Computing multiple watchman routes[END_REF]. Another method known as Zookeeper's problem was proposed [START_REF] Jonsson | An approximative solution to the zookeeper's problem[END_REF]. Given a simple polygon (the zoo) with a set of 𝑘 convex disjointed polygons (the cages) inside, this solution consists in finding the shortest path inside the zoo that touches every cage without getting inside it. The majority of all these algorithms looks for the optimal path that minimizes the total distance traveled while guarding a polygon. In our case, we need to find a path in the navigability polygon that ensures the guard of the visibility polygon. With this intention, we propose to find a graph inside the navigability polygon which guarantees the total oversight of the visibility polygon. To visualize the totality of a convex area, only one point of view (the center) is sufficient. From this statement, the graph will be built from the points of view 𝐶 𝑘 selected in each convex area. Thus Fig. 1: Polygons Intersection we propose to decompose the visibility polygon into several convex zones. Since the navigability polygon is included in the visibility polygon, we will calculate the inter-polygons between each convex polygon of Σ and the navigability polygon. Fig 1 shows the decomposition of the visibility polygon into several convex polygons and the intersections with the navigability polygon. Several algorithms have been proposed to decompose a polygon in multiple areas. Keil's algorithm [START_REF] Mark | Polygon decomposition[END_REF], [START_REF] Keil | On the time bound for convex decomposition of simple polygons[END_REF]] consists in splitting a polygon into many convex sub-polygons by eliminating all reflex vertices. We suggest the use of Mark Bayasit's algorithm [START_REF] Bayazit | Decomposing polygons into convex regions[END_REF] which is an optimized version of the Keil's algorithm. Therefore we connect all the nodes of 𝑇 by a set of arcs 𝛾 so that these arcs do not intersect the navigability polygon.

B. Optimization Process

We intend to find out on the previous graph 𝑇 = (𝐶, 𝛾) a set of viewpoints which maximizes the entropy 𝐸 while ensuring that the majority of the environment areas are still visible. As shown previously in Algorithm 1 the final spheres can belong to the vertices of the starting graph 𝑇 or on its edges. While browsing the initial graph, we will select the best viewpoints satisfying the criteria of entropy and visibility. The optimization algorithm consists in selecting a start point and an end point on the initial graph 𝑇 . Starting from the start point, we will select points of view with entropy greater than a threshold 𝛼 and covering all vertices in visibility polygon. To select the next point on the graph we propose a bijection function Π between each two linked nodes. This function allows to split each edge into 𝑁 points. For each edge 𝛾 𝑘 connecting two nodes 𝐶 𝑘 and 𝐶 𝑘+1 . Π the bijection between 𝛾 𝑘 and the interval [0 1] Π :

⃒ ⃒ ⃒ ⃒ [01] -→ 𝛾 𝑘 𝑟 ↦ -→ Π(𝑟)
𝑟 𝑖 = 𝑖 𝑁 -1 , 𝑖 = 0 to 𝑁 -1 (2)
Fig 2 shows an example of edge discretization. The optimization process will be repeated until cover the entire visibility area. At the end of this process we will obtain an optimal graph of spheres. In the next section we will detail the process of evaluating the obtained results.

IV. EXPERIMENTAL RESULTS

A. Evaluation

To evaluate the quality of a map and estimate its performances during the localization, several solutions have been proposed. In [START_REF] Merzić | Map quality evaluation for visual localization[END_REF], the authors propose to predict performances of localization in context of the visual cartography thanks to a score that is estimated for each position 𝑣 according to the map structure. This score is calculated by counting the number of points of interest visible from 𝑣 and its closest neighbors. The authors showed that the score provided by their method is strongly correlated with the true performance of localization. The solution suggested by Zhen et al. [START_REF] Zhen | Robust localization and localizability estimation with a rotating laser scanner[END_REF] consists in calculating a criterion called localizability. This criterion is calculated from the geometric characteristics in each direction for each point in the map. To estimate the localizability of a given position first, they determine the set of visible points from this position. Then they estimate the surface normals for each visible point. After that they accumulate the normals in one matrix 𝑁 which describes the set of observable constraints of the given position.

𝑁 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑛 1𝑥 𝑛 1𝑦 𝑛 1𝑧 𝑛 2𝑥 𝑛 2𝑦 𝑛 2𝑧 . . . . . . . . . 𝑛 𝑘𝑥 𝑛 𝑘𝑦 𝑛 𝑘𝑧 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
In [START_REF] Zhen | Robust localization and localizability estimation with a rotating laser scanner[END_REF] they showed that to precisely locate a robot we needs to be able to adequately constrain its pose in all three dimensions. They proposed to examine the singular values 𝑁 with SVD as 𝑈 Σ𝑉 . Theoretically, we should be Fig. 3: Evolution of localizability while moving along the x-axis able to locate as long as all three of the singular values are non-zero. In practice, they have calculated the localizability 𝐿 as the minimum singular value of 𝑁 :

𝐿 = 𝑚𝑖𝑛(𝑑𝑖𝑎𝑔(Σ)) (3)
To test this criterion, we have selected several positions on a trajectory and then calculated localizability for each sphere in each position. All spheres are in the same reference (𝑥, 𝑦, 𝑧) and have the same height 𝑧 = 1.5𝑚 corresponding to the average height of the agent. Fig 3 shows the evolution of localizability by moving along the x-axis. We show that the localizability is better when the spherical image contains fewer empty areas. These empty areas are the invisible areas from this point of view. Therefore, this criterion reflects the quantity of geometric constraints in a given position. If this point of view allows sufficient visualization of the geometric constraints, it will be easily localizable. In our case, we need that each sphere in the optimal graph 𝑋 must have a good localizability value.

Fig. 4: Large-Scale 3D mesh

B. Results

We have applied our method on two databases which describes two parts of Rouen city. They cover about 33000 𝑚 2 of a developed urban area that contains urban structures such as houses, buildings and roads. The first database is a 3D textured, triangulated and labeled mesh ( To evaluate the gain in memory we have proposed to compute the compression ratio which is defined as the ratio between the size of the final set of spheres and the size of the initial mesh and it is notated as a percentage. This dataset permitted the evaluation of our solution's performance by using all semantic, photometric and geometric characteristics together. Table I shows the final results on both databases. and entropy of the studied environment, thus an optimal image registration. Reducing the size of the initial map up to 1% allows to gain in term of memory and time of computation. Indeed we can integrate easily this minimal graph on navigation system with low resources or send it over a network from a server to an agent. We have optimized the computation time of our algorithm. Indeed, the construction of a sphere of resolution 360𝑥180 took 20𝑚𝑖𝑛 because of the volume of the initial data. To do this, we have used parallel computing and organized the original data in order to facilitate research. As a result, we reduced the calculation time of a sphere to 0.9 𝑚𝑖𝑛. This computation time can be further improved to facilitate the use of our algorithm in real-time applications.

V. CONCLUSIONS

The developped method throughout this paper allows to summarize efficiently a large-scale 3D Mesh. The summarizing process is based on the extraction of several spherical views representing a part of the initial map. This spherical representation contains semantic, photometric and geometric information. This new method of summarizing 3D maps allows us to facilitate several navigation tasks when applied in intelligent transportation systems (localization, route planning, obstacle avoidance, ...) by reducing significantly the calculation time and the required memory size. We also believe that using semantic information permits the development of a precise summary map by rejecting unnecessary Fig. 7: Optimal RGB-D-L Spherical Image Graph Fig. 8: Second Data Set localization data such as points belonging to dynamic objects (cars, pedestrians, ... ). In our future work, we aim to extend our algorithm to treat the case where the navigable areas that would contain holes. The next important step will also consist in using the set of spherical views and the semantic information of the navigable area in order to produce a navigation graph. This graph will propose an adapted path to follow by the agent according to its mobility capabilities (car, pedestrian, bike, mobile robot, ...). In addition, we aim to add other criteria in the optimization process. For example, for the convergence of localization algorithms, we can impose a maximum distance 𝐷 between spherical images. In addition, we aim to develop our localization system in order to evaluate the localizability of the summarized map.

Algorithm 1

 1 Compute the sphere graph 𝑋 = {𝑋 𝑖 , 𝑖 = 1..𝑛} Require: 𝑇 = (𝐶, 𝛾) Barycenters graph Σ = {𝑉 𝑗 , 𝑗 = 1..

Fig. 2 :

 2 Fig. 2: Edge Discretization

  Fig 4) and contains over 5 millions faces. In this dataset, we have 15 classes of labels, namely {0: Misc, 1: Road,2:RoadMarks, 3: Terrain sidewalk, 4: Building, 5: vegetation, 6: Tree, 7: Pedestrian, 8: cars,9: Bicycle, 10: Motorbike,11: Traf-ficLight, 12: TrafficSign, 13: Pole, 14: sky}.

Fig. 5 :

 5 Fig. 5: Polygons Extraction

Fig. 6 :

 6 Fig. 6: RGB-D-L Spherical Image

  𝑘}𝑉 𝑖 vertices of Σ Starting point 𝐶 𝑠 ∈ 𝐶 End point 𝐶 𝑒 ∈ 𝐶 Ensure: Mark all vertices in Σ as uncovered Mark all vertices in Σ visible from 𝐶 𝑠 and 𝐶 𝑒 as covered 𝑋 ← {𝐶 𝑠 , 𝐶 𝑒 } 𝑞 ← 𝐶 𝑠 𝛾 𝑘 ← the edge in 𝛾 whose starting point is 𝑞 while ∃𝑥 an uncovered point in Σ do Π(𝛾 𝑘 ) =Discretization of the edge 𝛾 𝑘 𝑝 ← first point in Π(𝛾 𝑘 ) to the right of 𝑞 such that 𝑝 cover some uncovered 𝑥 ∈ Σ and 𝐸(𝑝) ≥ 𝛼 𝑋 ← 𝑋 ∪ {𝑞} Mark all points ∈ 𝑥 as covered 𝑞 ← 𝑝 if(𝑞 ∈ 𝐶) 𝛾 𝑘 ← the edge in 𝛾 whose starting point is 𝑞 end while return 𝑋

TABLE I

 I

		: Final Results
	Data Sets	Data Set 1 Data Set 2
	Trajectory Length	100m	500m
	Number of Triangles 5 millions 30 millions
	Number of Spheres	10	52
	Compression Ratio	99%	99%
	The final set of spherical images allows to summarize the
	initial 3D mesh. Most of the points visible from thes spheres
	belong to relevant points such as buildings, vegetation and
	TrafficSign. As a result, we have significantly decreased
	the size of the map. Nevertheless, we have succeeded to
	keep a maximum number of salient points. To evaluate

our solution, we have computed the Localizability of each spherical image, as defined in the previous section. All the spheres have a high localizability values which are greater than 0.5. Based on what has been shown in

[START_REF] Zhen | Robust localization and localizability estimation with a rotating laser scanner[END_REF]

, all spheres are easily localizable. This ensures a good localization in the sphere graph. Our algorithm guarantees maximum visibility

VI. ACKNOWLEDGMENTS

This work takes part in the ANR-15-CE23-0010-01 pLaT-INUM project. This project has been funded with the support from the French National Research Agency (NRA).


Nicolas Guillou
email: nicolas.guillou@cerema.fr

Georges Chapalain
email: georges.chapalain@cerema.fr

Wave Energy Potential in the Sea of Iroise

Keywords: renewable energy, wave power, SWAN, unstructured grid, Brittany, western Europe

The wave energy resource in coastal areas of the Sea of Iroise (western Europe) is evaluated with an unstructured version of the phase-averaged wave model SWAN (Simulating WAves Nearshore) for a eight-year period between 2004 and 2011. Numerical predictions are calibrated and evaluated against available measurements of significant wave height and peak period at nine wave buoys in offshore and nearshore waters. In spite of strong energy dissipation in shallow water, up to 60 % of wave power, the medium-term evaluation of the resource reveals major coastal energetic patterns, between 20 and 35 kW m -1 , off the isles of Ushant and Sein and in the nearshore areas of the Crozon Peninsula, the bay of Audierne and the northern coastline. The variability of wave power production is estimated revealing, in accordance with previous numerical estimations over European shelf seas, significant inter-annual and inter-seasonal evolutions at the scale of the Sea of Iroise. Changes are particularly noticeable during the winter period with opposite situations in the distribution of monthly averaged wave energy flux. In the perspective of implementation of wave energy converters, the local distributions of energy flux against periods and directions are finally investigated in areas of maximum mean wave power.

I. INTRODUCTION

Marine renewable energy (MRE) is recognised by many countries as a promising alternative to mitigate the effects of climate change induced by human activities and achieve future energy security [START_REF] Shields | Marine renewable energy: The ecological implications of altering the hydrodynamics of the marine environment[END_REF]. Among the different MREs, wave energy constitutes an abundant resource with (1) a global worldwide power estimated around 2 TW [START_REF] Gun | Quantifying the global wave power resource[END_REF] and (2) high power density in the nearshore areas [START_REF] Clément | Wave energy in europe: current status and perspectives[END_REF]. Numerous technologies are currently in development to improve extraction performance for electricity generation [START_REF] Lõpez | Review of wave energy technologies and the necessary power-equipment[END_REF]. Accurate energy assessments are thus requested by potential developers to optimise location and design of wave energy converters (WEC).

The available resource is commonly evaluated with the wave energy flux (also denominated the wave power or potential) characterising the transport of energy per unit length of wave front and expressed in kW m -1 . Among the different methods implemented to assess wave power (e.g., buoys records, advanced combination of measurements and large-scale numerical simulations) [START_REF] Soares | Numerical evaluation of the wave energy resource along the atlantic european coast[END_REF], third-generation spectral wave models are traditionally implemented to extrapolate the resource at extended time scales and locations approaching, in particular, the associated variability in coastal areas [START_REF] Angelis-Dimakis | Methods and tools to evaluate the avaibility of renewable energy sources[END_REF]. Numerous numerical studies have thus been conducted to refine wavepower assessments in the most energetic regions of the world.

With a total resource amount estimated over 230 GW ( [START_REF] Gun | Quantifying the global wave power resource[END_REF], [START_REF] Pontes | Assessing the european wave energy resource[END_REF]), particular attention has been given to European shelf seas with numerical applications at spatial scales ranging from regional domains (e.g., [START_REF] Mattarolo | Wave energy resource of the french coasts: the anemoc databse applied to the energy yield evaluation of wave energy converters[END_REF]- [START_REF] Neill | Wave power variability over the northwest european shelf seas[END_REF]) to shallow-water coastal areas with increased computational resolutions (e.g., [START_REF] Iglesias | Wave energy potential along the death coast (spain)[END_REF]- [START_REF] Gonc ¸alves | Wave energy conditions in the western french coast[END_REF]). A detailed review of these modelling has recently been established by Guedes Soares et al. [START_REF] Soares | Numerical evaluation of the wave energy resource along the atlantic european coast[END_REF].

Whereas coastal applications, mostly based on structured regular or curvilinear computational meshes, identify generally well the distribution of nearshore energetic patterns, unstructured grid computation offers attractive options [START_REF] Shields | Marine renewable energy: The ecological implications of altering the hydrodynamics of the marine environment[END_REF] increasing spatial resolution at the coast, (2) solving numerics and physics mismatches boundaries problems of embedded domains and (3) optimising CPU performance with a reduced number of grid nodes.

The present study extends these coastal numerical evaluations relying on an unstructured version of the phaseaveraged wave model SWAN (Simulating WAves Nearshore) [START_REF] Booij | A third-generation wave model for coastal regions, 1. model description and validation[END_REF] implemented in the Sea of Iroise at the western extend of Brittany (Fig. 1). The site of application is considered as one of the most energetic region along the French coasts with a mean offshore wave energy flux estimated around 50 kW m -1 [START_REF] Mattarolo | Wave energy resource of the french coasts: the anemoc databse applied to the energy yield evaluation of wave energy converters[END_REF].

The approach retained here relies on a comparison between numerical predictions and available observations of the significant wave height and the peak period. Field measurements are collected at nine offshore and nearshore locations including archived data of long-term observation systems and field data acquired during short-term campaigns in coastal areas (section II-A). SWAN (version 40.91), modified to integrate an enhanced dissipation term for current-induced whitecapping [START_REF] Van Der Westhuysen | Spectral modeling of wave dissipation on negative current gradients[END_REF] (section II-B), includes heterogeneous parameterisations of bottom roughness and variations of tidal free-surface elevation and depth-averaged currents (section II-D). A medium-term evaluation of the wave energy resource for a eight-year period, between 2004 and 2011 (section III-A), reveals major coastal energetic patterns exhibiting significant inter-annual and interseasonal variabilities (section III-B). In the perspective of WEC implementation, the local distributions of energy flux against periods and directions is finally investigated in areas of maximum mean wave power (section III-C).

II. MATERIALS AND METHODS

A. Measurements Description

Wave buoy data here used consist of the archived measurements of the French CANDHIS database ("Centre d'Archivage National de Données de Houle In Situ", Cerema, France) (points 02902, 02911 and 02914) complemented by observations acquired during campaigns conducted off and in the bay of Douarnenez (points E1 to E3 and T1 to T3) (Fig. 1). Long-term measurements are only available at the offshore wave buoys 02902 and 02911. Whereas the associated data may lack during periods over three months in relation to measuring-system malfunction, it covers globally well the period of interest between 2004 and 2011. Medium term observations are available at the wave buoy 02914 off Penmarc'h headland over the period November 2009-February 2010. Complementary measurements were acquired at six offshore and nearshore wave buoys (E1-E3, T1-T3) during three short-term campaigns: in 15-22 April 2005, in 13-23 September 2005 and in 10 April-10 May 2006 [START_REF] Guillou | Rôles de l'hétérogénéité des sédiments de fond et des interactions houle-courant sur l'hydrodynamique et la dynamique sédimentaire en zone subtidale -applications en manche orientale et à la pointe de la bretagne[END_REF]. The associated instrumentation network is deployed in water depths ranging from 10-15 m off beaches of the bay of Douarnenez (points E3, T2 and T3) to 110 m at the wave buoy 02902 off the isle of Ushant.

B. Model Theoretical Formulation

SWAN computes the evolution of the wave action density N (= E/σ with E the wave energy density distributed over intrinsic frequencies σ and propagation directions θ) using the time-dependent spectral action balance equation:

dN dt = S tot σ (1)
where t denotes times. The right-hand side of this equation contains the source and sink terms of physical processes which generate, dissipate or redistribute wave energy:

S tot = S nl4 + S nl3 + S in + S wc + S bot + S brk + S wc,cur . (2)
Parameterisations here adopted for each terms are briefly detailed hereafter. The redistribution of energy by nonlinear quadruplet wave-wave interactions S nl4 is computed with the Discrete Interaction Approximation of Hasselmann et al. [START_REF] Hasselmann | Computations and parameterizations of the nonlinear energy transfer in a gravity wave spectrum. part ii: Parameterizations of the nonlinear transfer for application in wave models[END_REF].

The non-linear triad re-distribution of wave energy S nl3 is approached with the Lumped Triad Approximation derived by Eldeberky [START_REF] Eldeberky | Nonlinear transformation of wave spectra in the nearshore zone[END_REF]. The transfer of energy from the wind to the waves S in and the dissipation of wave energy due to whitecapping S wc are approached with the saturationbased model of van der Whesthuysen [START_REF] Van Der Westhuysen | Advances in the spectral modelling of wind waves in the nearshore[END_REF] combined with the wind input formulation proposed by Yan [START_REF] Yan | An improved wind input source term for third generation ocean wave modelling[END_REF]. The sink term of energy dissipation by bottom friction S bot is computed according to the formulation proposed by Madsen et al. [START_REF] Madsen | Spectral wave attenuation by bottom friction: theory[END_REF].

Energy dissipation in random waves due to depth-induced breaking S brk is quantified according to Battjes and Janssen [START_REF] Battjes | Energy loss and set-up due to breaking of random waves[END_REF]. An additional dissipation term S wc,cur recently proposed by van der Westhuysen [START_REF] Van Der Westhuysen | Spectral modeling of wave dissipation on negative current gradients[END_REF] is finally included to limit the overprediction of wave height on negative current gradients (accelerating opposing currents or decelerating following currents). Its implementation is conducted following previous recent calibrations of this dissipation term in SWAN (e.g., [START_REF] Van Der Westhuysen | Spectral modeling of wave dissipation on negative current gradients[END_REF], [START_REF] Fairley | The effects of a severn barrage on wave conditions in the bristol channel[END_REF]). In SWAN, wave power is approximated with default quantities outputs of energy transport components along (x) and (y) directions:

P SW AN = P 2 SW AN,x + P 2 SW AN,y 1/2 . ( 3
)
with

P SW AN,x = ρg 2π 0 ∞ 0 c x Edσdθ (4)
and

P SW AN,y = ρg 2π 0 ∞ 0 c y Edσdθ ( 5
)
where ρ is the water density, g is the acceleration of gravity, θ is the wave direction and c x and c y are the propagation velocities of wave energy in spatial space [START_REF] Team | Swan cycle iii. version 40.91 abc -scientific and technical documentation[END_REF]. This computation of wave power based on the summation of squared (x) and (y) components neglects additional terms included when directly estimating wave potential from the amplitude of wave energy:

P = ρg 2π 0 ∞ 0 |c g + ū|Edσdθ (6)
where c g is the group velocity and ū is the depth-averaged current. Nevertheless, as demonstrated by wave-power assessments based on SWAN (e.g., [START_REF] Iglesias | Wave energy potential along the death coast (spain)[END_REF]- [START_REF] Gonc ¸alves | Wave energy conditions in the western french coast[END_REF], [START_REF] Rusu | Wave energy pattern around the madeira islands[END_REF]), this computational method provides accurate estimation of the wave energy flux in nearshore areas.

The wave action balance equation is expressed on a cartesian coordinate system and solved on an unstructured grid, a constant directional resolution and an exponential frequency distribution. Further details about the mathematical expressions of sources and sinks are available in SWAN technical documentation [START_REF] Team | Swan cycle iii. version 40.91 abc -scientific and technical documentation[END_REF] and associated scientific literature [START_REF] Zijlema | Computation of wind-wave spectra in coastal waters with swan on unstructured grids[END_REF].

C. Numerical Resolution

The wave action balance equation is solved on an unstructured computational grid adopting a constant directional resolution and an exponential frequency distribution. The time integration of the action equation is performed with an implicit first order Euler scheme improving the stability of the resolution with the time step retained. The geographic propagation terms are approximated with an upwind difference scheme while discretisation in spectral space is performed with a hybrid central/upwind scheme. Source-sink terms are integrated semi-implicitly following Patankar's [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] rules and linearised with the Newton-Raphson iteration if strongly nonlinear. The solution is finally found by means of an active solver integrating a sweeping algorithm to update the solution at each vertex. Further details about the numerical resolution are available in [START_REF] Zijlema | On convergence behaviour and numerical accuracy in stationary swan simulations of nearshore wind wave spectra[END_REF], [START_REF] Zijlema | Computation of wind-wave spectra in coastal waters with swan on unstructured grids[END_REF] and [START_REF] Team | Swan cycle iii. version 40.91 abc -scientific and technical documentation[END_REF].

D. Model Setup

SWAN is set up on an unstructured computational grid covering the Sea of Iroise and comprising 9971 nodes and 18443 elements with a size of 10 km offshore to less than 300 m nearshore (Fig. 2). The model runs with 30 exponentially spaced frequencies ranging from 0.05 to 1 Hz, 30 evenly spaced directions and a time step of 15 min. In nearshore areas, wave energy dissipation by bottom friction is parameterised with a heterogeneous roughness length scale determined combining bottom types with associated observations of roughness parameter compiled by Soulsby [START_REF] Soulsby | The bottom boundary layer of shelf seas[END_REF]. The offshore bottom roughness is set to an uniform value of k n = 10.5 mm (e.g., [START_REF] Muller | High-resolution atmospheric forcing for regional oceanic model: the iroise sea[END_REF], [START_REF] Guillou | Modelling impact of bottom roughness on sea surface temperature in the sea of iroise[END_REF]). Wind velocity components at 10 m above the free surface are provided at a time step of three hours and a spatial resolution of 10 km by the meteorological model ALADIN ("Aire Limitée, Adaptation dynamique, Développement Inter-National", Météo-France). SWAN integrates variations of tidal free-surface elevations and depth-averaged currents predicted by the bidimensional horizontal circulation model TELEMAC 2D [START_REF] Hervouet | Hydrodynamics of free surface flows, modelling with the finite element method[END_REF] set up at an extended computational grid covering the initial SWAN unstructured mesh. The wave model is finally driven by wave components (significant wave height, peak period, direction and spreading) predicted by a regional run of Wave Watch III at the scale the north-eastern Atlantic ocean with a spatial resolution of 18 km in the context of the IOWAGA (Integrated Ocean WAves for Geophysical and other Applications, Ifremer) project.

The wave model is run during eight years between 2004 and 2011 which corresponds to a period when most measurements and forcings were available. Model performances are assessed with the standard statistical parameters of the mean absolute bias the index of agreement introduced by Willmott [START_REF] Willmott | On the validation of models[END_REF] as

BIAS = 1 N i=N i=1 |x i -y i | (7)
RE = 1 - i=N i=1 (x i -y i ) 2 i=N i=1 (|x i -x| + |y i -x|) 2 (8)
and the Pearson's correlation coefficient

R = i=N i=1 (x i -x) (y i -ȳ) i=N i=1 (x i -x) 2 i=N i=1 (y i -ȳ) 2 1/2 (9)
where N is the number of data in the discretised time series considered, (x i ) and (y i ) represent the two sets of measured and simulated values and x and ȳ are the mean values of observed and modeled data, respectively.

III. RESULTS AND DISCUSSION

A. Comparison with Point Measurements

The quality of numerical results issued from the depthaveraged circulation model TELEMAC 2D has first been estimated. A preliminary confirmation of model predictions has been performed against water depth's observations in harbors of Le Conquet and Brest of the tide gauge network RONIM ("Réseau d'Observation du NIveau de la Mer"). This evaluation has been extended to mean near-surface spring tidal currents' measurements compiled by the SHOM ("Service Hydrographique et Océanographique de la Marine") at eight points titled S1 to S8 evenly spaced over the computational domain (Fig. 1). Predictions reproduce generally well the temporal variations of the amplitude and direction of the currents at the eight sites considered (Fig. 3). Whereas the model tends to overestimate currents amplitude at point S5, differences are globally restricted to less than 15 %. The currents direction is furthermore fairly well approached even in areas influenced by the formation of headland-associated eddies (point S7).

Wave model predictions are then evaluated on the basis of statistics computed for the significant wave height H s and the peak period T p at the nine locations considered over the different periods of measurements (Fig. 1). Predictions reproduce fairly well the temporal evolutions of long-term observations at offshore locations 02902 and 02911 (Fig. 4). The associated statistics (Table I) fall in the range of estimations performed by Gonc ¸alves et al. [START_REF] Gonc ¸alves | Wave energy conditions in the western french coast[END_REF] at point 02902 and Boudiere et al. [START_REF] Boudiere | A suitable metocean hindcast database for the design of marine energy converters[END_REF] at point 02911. At both measurement sites, good agreement is obtained for H s with indexes RE over 0.93. The quality of model predictions at these offshore locations is also exhibited with computed correlation coefficients over 94 %. Whereas increased differences are obtained in peak period estimations, the associated indexes of agreement remain over 0.81. On medium and short-term measurements, comparisons between numerical results and observations reveal slightly better estimations of H s in deep waters (points 02914 and T1) than in coastal waters (points T2 and T3). This comparison can not be accurately established between points E1, E2 and E3 as measurements cover different periods in offshore and nearshore waters. Differences obtained at coastal locations appear primarily in the bay of Douarnenez where the model tends to overestimate H s measurements at point T2 while underestimating it at points E3 and T3 (Fig. 5). Nevertheless, predictions of the significant wave height remain satisfactory with a minimum index of agreement RE equal to 0.85 at point T3. The model approaches also the observed semi-diurnal modulation of the significant wave height induced by the tidal current and particularly noticeable at point E2 in September 2005. Whereas the correlation coefficient R reaches values under 60 % at points E1 and T3 for the estimation of the peak period, T p predictions are in general good agreement with measurements reproducing the observed increase during storm events. No particular bias is thus denoted on predictions of H s and T p at the nine available measurements sites. 6). As pointed out by Rusu and Guedes Soares [START_REF] Rusu | Numerical modelling to estimate the spatial distribution of the wave energy in the portuguese nearshore[END_REF], whereas a close correlation exists between significant wave height H s and wave power P , the resulting field may present local differences revealing especially the influence of the group velocity on computation of wave energy flux (Eqs. 4 and5).

In the present investigation, slight differences are obtained in the nearshore areas, locations with the most important H s matching sites with maximum wave power. Strongest shallow-water mean wave power P = 35 kW m -1 is reached off the isle of Ushant characterised by the largest exposure to North-Atlantic incoming waves. Waves experience however more significant energy dissipation, mainly by bottom friction and wave breaking, in the eastern nearshore areas of the Sea of Iroise. Mean wave power decreases thus from offshore values around 27 kW m -1 off the Crozon Peninsula to less than 12 kW m -1 at the entrance of the bay of Douarnenez. This dissipation results in a strong spatial variability of coastal wave energy flux with values (1) reaching 20 kW m -1 in exposed shallow-water areas of the northern coastline, the Crozon peninsula and the bay of Audierne and (2) decreasing below 10 kW m -1 in sheltered areas behind the isles of Ushant and Sein or in the bays of Brest and Douarnenez. 2) Inter-annual and Inter-seasonal Variabilities: In the present investigation, the mean offshore wave energy flux estimated around P = 40 kW m -1 is lower than quantification of P = 50 kW m -1 established by Mattarolo et al. [START_REF] Mattarolo | Wave energy resource of the french coasts: the anemoc databse applied to the energy yield evaluation of wave energy converters[END_REF] on the basis of large-scale modelling at the scale of the European continental shelf over a period of 23 full years . It lies furthermore over the evaluation of P = 28 kW m -1 recently performed by Gonc ¸alves et al. [START_REF] Gonc ¸alves | Wave energy conditions in the western french coast[END_REF] in the western French coast for the three-years period between 1998 and 2000. Although these differences may be attributed to the numerical methods retained including various spatial and temporal resolutions, the variability of wave climate between simulation periods appears also to be taken into account. Indeed, as exhibited by Neill and Hashemi [START_REF] Neill | Wave power variability over the northwest european shelf seas[END_REF], wave power at the scale of the North-western European shelf seas is characterised by a strong inter-annual variability particularly noticeable during the winter period in close correlation with the evolution of the North Atlantic Oscillation (NAO). Nevertheless, assessment of wave-energy resource provides generally potential WEC developers with only averaged quantities neglecting the relative uncertainty associated with the variability of wave climate (e.g., [START_REF] Mackay | Uncertainty in wave energy resource assessment. part 1: historic data[END_REF], [START_REF]Uncertainty in wave energy resource assessment. part 2: variability and predictability[END_REF]). Further investigation is thus conducted about the spatial and temporal variabilities of wave power in the Sea of Iroise.

The seasonal evolution of the wave energy resource in the Sea of Iroise (Fig. 7) is in accordance with previous largescale estimations (e.g., [START_REF] Neill | Wave power variability over the northwest european shelf seas[END_REF], [START_REF] Gonc ¸alves | Wave energy conditions in the western french coast[END_REF]) exhibiting a clear contrast between winter energetic and summer low-energetic months. Between November and March, monthly offshore wave power remains thus globally over 45 kW m -1 reaching peak values over 75 kW m -1 in December and January. For the rest of the year, deep-waters predictions are globally restricted to values below 20 kW m -1 with an exception for the month of October when the offshore mean wave energy flux reaches 35 kW m -1 . This seasonal evolution is however characterised by a strong temporal variability as exhibited with monthly average wave power in early 2008 and 2011 (Fig. 8). An opposite situation is thus obtained between these two years with energetic conditions appearing (1) in 2008 during the months of January and March and (2) in 2011 during the month of February. Time-series of monthly wave power averaged over computational domain (Fig. 9) identify more precisely the inter-annual variability of energy flux in the Sea of Iroise. Following largescale estimations performed by Neill and Hashemi [START_REF] Neill | Wave power variability over the northwest european shelf seas[END_REF], a clear contrast is exhibited between (1) Spring and Summer when P remains nearly stable and (2) Autumn and Winter when P shows significant annual differences. The wave energy flux experiences thus prominent inter-annual variations during the month of November. Indeed, whereas November is classified at rank four among the most energetic months (Fig. 7), the most energetic period is obtained in November 2009 with mean wave power reaching 38 kW m -1 .

As pointed out by Neill and Hashemi [START_REF] Neill | Wave power variability over the northwest european shelf seas[END_REF], when considering a typical year, differences may also appear between the most energetic months in the geographic locations of peak wave power. Previsions obtained in 2007 at the scale of the European shelf seas were thus exhibiting peak wave power on the northwest of Scotland and Ireland in January and March while displaying focused wave energy flux on the Celtic Sea and Bay of Biscay in February. At the local scale of the Sea of Iroise, smoothed differences are obtained for the period of interest between 2004 and 2011. Synoptic investigations of monthly predictions confirm that nearshore spatial distribution of wave energy flux remains nearly the same with peak wave power matching sites identified on the overall average field (Fig. 6).

C. Local Analysis of Wave Power

1) General Selection: Several different technologies are currently in development for transforming wave power into electricity [START_REF] Falcão | Wave energy utilisation: a review of the technologies[END_REF] with operational range reaching maximum efficiency in restricted intervals of periods and directions. In the perspective of WEC selection, design and optimisation, the present assessment of wave energy flux is investigated further at nearshore locations with the largest energy along the coast of western Brittany. Taking into account the population density in the proximity of these areas, sites located off the isles of Ushant and Sein are not considered here. Four locations titled #1 to #4 are retained in water depths between 15 and 32 m at a distance to the coastline around 2.1 km (Fig. 6, Table II). The associated areas correspond to major expositions to North Atlantic incoming waves: along the northern coastline (point #1), off Crozon Peninsula (point #2), in the bay of Audierne (point #3) and off Penmarc'h headland (point #4). While the average nearshore wave power lies below 12 kW m -1 , these hots spots concentrate wave energy flux with mean values varying from 15.1 kW m -1 off Crozon Peninsula to 23.3 kW m -1 off Penmarc'h headland.

2) Distributions against periods and directions: These four locations present nearly the same distributions of wave power against peak periods (Fig. 10). Waves with longest periods over 12 s contribute during nearly 35 % of the time to the maximum values of the wave energy flux. The local mean wave power exceeds, in this range of periods, 31 kW m -1 reaching 47 kW m -1 off Penmarc'h headland. Waves with periods between 8 and 12 s are the most frequent with a percentage of occurrence estimated around 54 %. The associated average wave energy flux is however reduced to values between 8 and 13 kW m -1 at the four locations considered.

In comparison, short-period waves appearing during 12 % of the time have a negligible contribution reduced to mean values below 4.5 kW m -1 . Whereas North Atlantic waves come predominantly from west and northwest [START_REF] Soares | Numerical evaluation of the wave energy resource along the atlantic european coast[END_REF], the distribution of wave power against the incoming direction presents more variabilities (Fig. 11) in relation to nearshore processes like depth and/or currentinduced refraction and dissipation by bottom friction and wave breaking. Maximum energy density is obtained off Penmarc'h headland from the western direction with wave power exceeding 40 kW m -1 during more than 16 % of the time. Minimum wave power occurs off Crozon Peninsula in relation to increased dissipation of wave energy by bottom friction. At this site, more than 50 % of wave power lies below 20 kW m -1 . Points #1 and #3 are characterised by intermediate densities of wave energy with incoming waves contributing to values over 20 kW m -1 during more than 16 % of the time. At these two locations, the wave energy flux is nearly concentrated along one direction, the north-western direction at point #1 and the south-western direction at point #3. WEC selection should thus aim for maximum efficiency in these ranges of periods and directions.

IV. CONCLUSIONS

The wave propagation model SWAN has been set up on an unstructured computational grid to investigate and evaluate the nearshore wave energy resource in the Sea of Iroise. Numerical results have been compared with available in-situ measurements of significant wave height and peak period at nine locations between 2004 and 2011. The main outcomes of the present study are the following:

1) Whereas the Sea of Iroise is a high energetic area with mean offshore wave power estimated around 40 kW m -1 over the period of interest, a strong energy dissipation is exhibited in shallow water (from 50 m depth) with mean values decreasing below 15 kW m -1 in coastal areas.

2) The variability of the wave energy resource has also been identified. Predictions exhibit strong inter-annual and inter-seasonal variabilities of wave power particularly noticeable over the winter months. At the scale of the Sea of Iroise, monthly variations of wave energy flux may thus present opposite situations during the most energetic periods.

3) The energy resource was investigated further at four locations with the largest average wave power exhibiting significant variabilities against the incoming wave directions. Whereas the site located off Penmarc'h headland presents the maximum energy density with wave power exceeding 40 kW m -1 during more than 16 % of the time, the locations identified along the northern coastline and off Audierne appear also very interesting for WEC implementation as the wave energy flux is nearly concentrated along one primary direction. The implementation of an unstructured version of SWAN gives promising results for the quantification of the wave energy flux in the coastal areas. Whereas the present investigation will benefit from extending comparisons of numerical predictions with nearshore measurements, modelling exhibits the remarkable energy resource in the area of the Sea of Iroise. This refined assessment of wave power provides finally potential developers with relevant results for (1) selecting the implementation and (2) optimising the design of WEC projects.

Fig. 1 .

 1 Fig.1. Bathymetry of the Sea of Iroise with locations of available measurements points (red circles for wave buoys and black squares for currents stations).

Fig. 2 .

 2 Fig. 2. Computational unstructured grid for the (blue line) TELEMAC2D and (red line) SWAN models.

Fig. 3 .

 3 Fig.3. Measured (black line) and computed (blue line) time series of (top) the near-surface tidal current amplitude and (bottom) its direction (anticlockwise convention from the East) at points S1 to S8 for spring tidal conditions.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Measured (black line) and computed (red line) time series of significant wave height and peak period at wave buoys 02902 and 02911 in 2009.

Fig. 6 .

 6 Fig. 6. Predicted annual average (top) significant wave height and (bottom) wave power in the Sea of Iroise over the period 2004-2011 with locations of points #1 to #4.

Fig. 7 .

 7 Fig. 7. Average monthly evolution of wave power in the Sea of Iroise over the period 2004-2011.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Monthly average wave power in the Sea of Iroise in January, February and March 2008 and 2011.

Fig. 10 .

 10 Fig. 10. (Top) Predicted average wave power and (Bottom) percentage of occurrence for peak periods between 0 and 8 s, 8 and 12 s and over 12 s at the four locations #1 to #4.

Fig. 11 .

 11 Fig. 11. Predicted directional distribution of wave power at points #1 to #4 over the period 2004-2011.

TABLE I OVERALL

 I STATISTICS FOR THE SIGNIFICANT WAVE HEIGHT h m0 AND THE PEAK PERIOD Tp AT THE NINE MEASUREMENTS POINTS CONSIDERED

							02902
	m0 (m)
	h
	T p (s)
	Jan	Feb	Mar	Apr	May	Jun	Jul 2009	Aug	Sep	Oct	Nov	Dec
							02911
	h m0 (m)
	T p (s)
	Jan	Feb	Mar	Apr	May	Jun	Jul 2009	Aug	Sep	Oct	Nov	Dec
												Wave	h m0	Tp
												buoys	BIAS	RE	R	BIAS	RE	R
													(m)	(s)
												02902	0.47	0.93 0.94	1.21	0.81 0.68
												02911	0.35	0.96 0.96	1.07	0.84 0.73
												02914	0.39	0.98 0.96	1.29	0.83 0.71
												E1	0.32	0.93 0.89	1.21	0.75 0.59
												E2	0.16	0.96 0.95	0.89	0.87 0.78
												E3	0.09	0.94 0.94	1.08	0.80 0.64
												T1	0.22	0.95 0.92	0.96	0.86 0.74
												T2	0.08	0.91 0.89	1.17	0.76 0.59
												T3	0.09	0.85 0.88	1.53	0.69 0.50
												B. Global Wave Power Assessment
												1) Spatial Distribution of Energetic Patterns: Coastal wave
												energetic patterns in the Sea of Iroise are exhibited displaying
												the mean annual wave power over the period 2004-2011 (Fig.

TABLE II LOCATION

 II , WATER DEPTH, DISTANCE TO THE COASTLINE AND MEAN WAVE ENERGY FLUX COMPUTED IN 2004-2011 AT SITES #1 TO #4.

	Sites	X	Y	Depth	Distance to the	Pmean
		(m)	(m)	(m)	coastline (km)	kW m -1
	#1	80159 2422126	15.5	2.2	16.6
	#2	86038 2377260	27.5	2.3	15.1
	#3	84396 2354266	31.0	2.1	17.8
	#4	96045 2333391	24.3	2.1	23.3

ACKNOWLEDGMENT

The authors are particularly grateful to Météo-France for supplying wind fields data obtained with the meteorological model ALADIN over the computational domain. Wave data used at the open boundary of the wave propagation model are issued from regional runs of WaveWatch III in the curse of the IOWAGA project. Wave buoy data used for model validation were provided from the French CANDHIS database (Cerema). Simulations were performed on computer facilities CAPARMOR (CAlcul PARallèle Mutualisé pour l'Océanographie et la Recherche). The present paper is a contribution to the research programs DIADEMS and FLUSED of the Laboratory of Coastal Engineering and Environment (Cerema, http://www.cerema.fr, http://memphys-lgce.fr.ht).