Nicolas Le Moine
email: nicolas.le_moine@upmc.fr

Pierre-Stéphane Gsell

A graph-based approach to glacier flowline extraction: an application to glaciers in Switzerland

Keywords: glacier flowline, skeleton, discrete curve evolution, optimum branching, directed minimum spanning tree

In this paper we propose a new, graph-based approach to glacier segmentation and flowline extraction. The method, which requires a set of glacier contours and a Digital Elevation Model (DEM), consists in finding an optimum branching that connects a set of vertices belonging to the topological skeleton of each glacier. First, the challenges associated with glacier flowline extraction are presented. Then, the three main steps of the method are described: the skeleton extraction and pruning algorithm, the definition and computation of a travel cost between all pairs of skeleton vertices, and the identification of the directed minimum spanning tree in the resulting directed graph. The method, which is mainly designed for valley glaciers, is applied to glaciers in Switzerland.

Introduction

Glacier morphology

Glaciers are moving ice bodies which flow under their own weight, due to the accumulation of solid precipitations on the higher slopes of a mountain range. As the strain rate increases, ice viscosity decreases and the accumulated ice literally 'flows' downslope. [START_REF] Bahr | Observations and analysis of selfsimilar branching topology in glacier networks[END_REF] first explicitly drew a parallel between rivers and valley glaciers (i.e. glaciers that are confined by topography, as opposed to ice caps). They showed that glaciers also exhibit branching topologies, and computed classical river network indices such as bifurcation and area ratios for glacier networks. This analogy stems from the fact that in most recent orogens where valley glaciers are found (the Alps, the Andes, the Himalayas, the Rocky Mountains, etc.), glacier inception took place in a topography previously shaped by fluvial erosion [START_REF] Gsell | Identifying the probabilistic structure of drained areas as a function of hypsometry in river networks[END_REF]. Bahr and Peckham also pointed out that self-similarity properties could provide a 'lever arm' for tackling glacier flow dynamics for complex geometries, just as these properties are used for treating subgrid, hydraulic propagation in complex river networks with concepts such as the Geomorphological Instantaneous Unit Hydrograph (Rodríguez-Iturbe and Valdés, 1979; [START_REF] Gupta | A representation of an instantaneous unit hydrograph from geomorphology[END_REF]. One of the reasons why this approach has not been given much attention is maybe the difficulty lying in the first step of identifying networks of glacier flowlines.

Limits of classical drainage network extraction methods

Figure 1 shows the downstream region of the Rhone glacier in Switzerland. In fluvial morphology, we typically find cross-sections such as A-A' with a concave topography in the talweg. This translates into V-or Ushaped (for former glacial valleys) elevation contours, with the lowest point roughly in the medial axis of the talweg. Hence, river network extraction from a DEM is relatively straightforward, except for problems such as flat areas or closed depressions (see e.g. [START_REF] Garbrecht | The assignment of drainage direction over flat surfaces in raster digital elevation models[END_REF][START_REF] Martz | The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models[END_REF].

Things are more complicated for ice-covered areas. In the accumulation (higher) area of the glacier, where hillslopes as well as valley floors are icecovered, the topography is still concave (B-B'): the surface of the ice is more or less homothetic to the bedrock surface (with lower roughness though). In contrary, in the ablation area the glacier is limited to a narrow ice tongue confined between lateral, ice-free hillslopes. Since ice thickness is maximum in the medial axis of the ice tongue, we have a convex cross-section (C-C') with seemingly two talwegs on each side of the glacier. Elevation contours in this area have the shape of a W with its two wedges pointing upstream, as opposed to the single wedge in concave topography. The central flowline of the glacier (i.e. the line of maximum ice thickness), which is also typically the line of the bedrock talweg, is a local maximum and not a local minimum of the ice surface (it looks like is a local water divide). Therefore, it cannot be extracted in a stable way from a DEM with classical algorithms.

Automatic methods for glacier flowline extraction

The problem of glacier flowline extraction has received some attention recently, due to the need of feeding glacier databases with attributes such as glacier length. A flowline or a set of flowlines has to essentially meet two requirements: (i) to stay as far as possible from the glacier boundary, and (ii) to cross elevation contours orthogonally. Le [START_REF] Bris | An automatic method to create flow lines for determination of glacier length: A pilot study with Alaskan glaciers[END_REF] propose to construct a set of waypoints located at the center of 'traverses' drawn perpendicular to a single, rectilinear 'main axis', and then connect them. However, the method can only extract one centerline per glacier.

Objectives of the study

In this paper, a new method is presented that aims at extracting glacier flowlines with an emphasis on preserving their tree-like structure, i.e. the structure of tributaries within the glacier.

As in Le Bris and Paul (2013), our method first identifies a set of waypoints (i.e., vertices of a graph) that are subsequently connected. However these waypoints are identified with a more general operation called skeletonization. Once these waypoints are identified (including special 'snout' vertices), we compute a travel cost between every pair of them: the cost function is designed so as to penalize displacements that deviate from the steepest slope direction. 2a), covering a total area of 1200 km 2 and mainly located in the headwaters of the Rhone, Rhine, and Danube rivers. The steps of the method are illustrated with a focus on a particular glacier complex in the Bernese Alps (Figure 2b), straddling the water divide between the Rhone and Rhine rivers.

Data

Digital Elevation Model

In this study we use the 25-meter Digital Elevation Model from the Swiss Federal Office of Topography (SwissTopo, 2004).

Randolph Glacier Inventory (RGI) glacier contours

Glacier outlines are taken from the Randolph Glacier Inventory (RGI, [START_REF] Arendt | Randolph Glacier Inventory -A Dataset of Global Glacier Outlines: Version 3.2. Global Land Ice Measurements from Space[END_REF]. The RGI provides a segmentation of glacier complexes (continuous ice bodies) into individual glaciers; we chose to dissolve (reaggregate) these elements and to work with the complexes in order test the ability of our approach to identify multiple snouts in such complexes. The segmentation is not a prerequisite and is even a by-product of our method.

Skeleton extraction and pruning

Glacier flowlines as a topological skeleton

The method proposed by Le Bris et al. (2013), which consists in picking the midpoints of 'traverses' drawn orthogonally to the glacier's 'main axis', is actually related to a topological operation called skeletonization, or medialaxis transform. A first, intuitive definition is the analogy with a grassfire [START_REF] Blum | A Transformation for Extracting New Descriptors of Shape[END_REF]: if one 'sets fire' simultaneously at all points on the border of a grass field enclosed in the object's boundary, then the skeleton is the set of points where two or more firefronts meet (see Figure 3a). The result (Figure 3b) is a 'thinned' version of the object (i.e., a set of edges, namely a graph) that preserves its essential geometrical and topological features (in the example of Figure 3, the skeleton edges look like the veins of the leaf).

We will see that this is a first interesting step to glacier flowline extraction.

Skeletonization using Voronoi tesselation

A second definition of the skeleton uses the concept of maximal disk (or maximal ball in dimensions higher than 2). A disk B is said to be maximal in set D if it is completely included in D, and such that if it is contained in any other disc B then B is not completely included in D. Mathematically,

B is a maximal disk in set D iif B ⊆ D B ⊂ B ⇒ B ⊆ D (1)
A maximal disk B(s) centered at point s is entirely contained in D and is interiorly tangent to the boundary ∂D in at least two different points, called generating points: these are two locations where firefronts originate from in the grassfire analogy, and they meet at the center of a maximal disk. The skeleton S(D) can be defined as the set of the centers of all maximal disks in D: We used the Qhull library [START_REF] Barber | The Quickhull algorithm for convex hulls[END_REF] to generate the tesselation. It allows to retrieve not only the Voronoi edges but also the pair of generating input sites (boundary points) for each edge, a highly valuable information for the subsequent steps. It is worth noting that the object may have inner boundaries ('holes' in the shape, such as inner rocks for a glacier), which will translate into cycles in the skeletal graph.

S(D) = s B(s) is a maximal disk in D (2)
RGI glacier contours were first densified so as to have at least one point every 50 m along the boundary, before running the tesselation with Qvoronoi.

All subsequent steps were then performed under Scilab (Scilab Entreprises, 2012), using the CLaRiNet library (Le Moine, 2013).

Skeleton pruning through Discrete Curve Evolution

As mentioned previously, the set of edges obtained with the Voronoi tesselation is only an approximation of the skeleton. Moreover, due to the noise in the boundaries, unsignificant boundary features can generate skeleton edges which do not reflect essential topological properties. Hence, we need to simplify, or 'prune', the skeleton i.e. remove many spurious edges.

In this study we use the pruning algorithm of [START_REF] Bai | Skeleton pruning by contour partitioning with discrete curve evolution[END_REF], who point out that not all n boundary points significantly contribute to the shape of the object: only a relatively small subset of these boundary (input) sites may be sufficient to describe the overall shape. Hence, if we select k ≤ n points on the boundary, we define a partition of the contour into k contour segments (subarcs). The pruning rule is to remove all skeleton edges whose generating points lie on a same contour segment, as illustrated in Figure 6.

It is important to note that this pruning process is not equivalent to reducing the number of input sites in the Voronoi tesselation (i.e. moving from right to left in Figure 5): the pruning removes some edges but the geometry of the remaining ones is not altered.

(a) (b) (c) The problem is to select the k most relevant points; [START_REF] Bai | Skeleton pruning by contour partitioning with discrete curve evolution[END_REF] propose to reduce the number of boundary points from n to k < n through Discrete Curve Evolution (DCE).

Let c 1 , . . . , c n be the n input sites sampled on the boundary, and let

u i-1,i = c i -c i-1
and u i,i+1 = c i+1c i be the vectors of the boundary edges incident to point c i (Figure 7). The contribution of this site to the overall shape is [START_REF] Latecki | Application of Planar Shape Comparison to Object Retrieval in Image Databases[END_REF]:

K(c i ) = θ i u i-1,i u i,i+1 u i-1,i + u i,i+1 (3)
where

θ i = (u i-1,i , u i,i+1 ) = arccos u i-1,i • u i,i+1 u i-1,i u i,i+1 sgn det(u i-1,i , u i,i+1 )
is the oriented turn angle at point c i (trigonometric, i.e. measured counterclockwise). This contribution is usually defined in 2D (planar computation), but there is no difficulty in extending it to a 3D curve. However, even if glaciers are located on steep topography, they remain objects in the geographical space i.e. with a rather flat aspect ratio, so that the addition of the vectors' z-coordinate would not change the contributions dramatically. • K is negative for a convex point (clockwise turn), and positive for a concave point (counterclockwise turn),

• the higher the absolute value |K|, the higher the contribution to the overall shape (through great segment lengths and/or large turn angle).

A point lying on a straight contour portion, i.e. such that θ = 0, will have a zero contribution (i.e. it can be removed without any loss of shape information).

We start the DCE with the n initial input sites and we removes the point If the shape has one or several inner boundaries, we simply add a loop on the boundaries in order to find the least-contributing site at each iteration.

c

Identification of snout vertices

In the previous section we illustrated the skeleton extraction and pruning with a simple shape (a leaf); we now apply these methods to RGI contours of Swiss glaciers. The main issues will be to decide the level k at which we should extract the skeleton for each glacier, and to use elevation data in order to identify snout and head vertices.

Choice of DCE level for each glacier

Consider a glacier with area A. Clearly, the larger the glacier, the more points will be needed on its boundary to correctly describe its shape or skeleton, i.e. the higher we will need to stay in the DCE hierarchy. Conversely, a small glacier with a single main axis without tributary could be correctly described at DCE level k = 2. A selection rule of thumb k = f (A) was devised empirically upon visual appreciation:

(k -2) = k 0 A γ
where is the floor function. We used k 0 = 13.5 and γ = 0.8 with A in km 2 . For example, the skeleton of Rhone-Trift glacier complex, with an area A = 35.8 km 2 , is still correctly described at DCE level k = 2+13.5(35.8) 0.8 = 238, i.e. with 238 sites on its boundary (Figure 9, left panel).

Hierarchical snout identification

The next step is the orientation of the skeletal edges in order to identify glacier snouts and heads. Indeed, skeletonization is a planar operation and we do not know if a 'leaf' vertex in the skeleton corresponds to the beginning of a flowline ('head' vertex), or to its end ('snout' vertex). Again, we will use the results of the DCE algorithm to create a hierarchy of potential snouts.

A skeletal vertex s is said to be a level-k snout if it meets the three following requirements:

(i) s is a leaf vertex of the skeletal graph of DCE level k, i.e. having a degree (number of incident edges) d = 1, (ii) s is at an elevation lower than its closest vertex s of degree d = 2 in the level k skeleton graph. s is such that all vertices on the path between s and s have degree 2: in the following we will call such a sequence of vertices of degree 2 (except the two extremal vertices) a stretch.

(iii) the generating points Gen(E s ) of its unique incident edge E s (pendant edge) lie on either side of an input site which is a local topographic minimum in the DCE level k. These seemingly complicated requirements translate the more intuitive notion that an edge is at a snout if it points downslope in a locally convex portion of the boundary that also points downslope, as shown in Figure 9.

On Figure 9, skeleton edges at level k are drawn in solid black lines, while edges that have been pruned at earlier levels are drawn with dotted lines.

There are several 'candidate' snouts in this area: skeleton vertices S 1 , S 2 , S 4 , S 6 and S 9 , which are all leaf vertices at the current pruning level, i.e. meeting requirement (i). All these 5 vertices also satisfy requirement (ii), since we have:

Leaf vertex

Nearest vertex of degree = 2 z(S 1 ) < z(S 3 ) z(S 2 ) < z(S 3 ) z(S 4 ) < z(S 5 ) z(S 6 ) < z(S 7 ) z(S 9 ) < z(S 8 ) However, Table 1 Hence, the only snout at level k in this part of the glacier is vertex S 2 . As input sites are iteratively removed from the glacier's outline in the DCE process, the combination of requirements allows a robust identification of convex portions of the boundary that point downslope, i.e. glacier tongues/snouts.

Finally, skeleton stretches which satisfy requirements (i) and (ii) but not requirement (iii) are further removed.

Hence, we obtain a snout hierarchy tied to the edge hierarchy: with a single threshold (the level k = f (A) of the DCE), we can extract both a set of skeletal edges and a set of snouts for a given continuous ice body of area A.

Definition of skeleton waypoints

Let us have a closer look at the final pruned skeleton in Figure 10. The solid black lines are the remaining edges while the black dots are just a subset of the vertices in the pruned skeleton that we downsampled for simplicity (these vertices are waypoints distant from at least 1000 m). Clearly the skeleton edges are far from a set of flowlines, except maybe in the ablation area. The location of glacier heads and snouts looks satisfying but major flaws appear:

• many edges in the skeleton significantly deviate from the local steepest slope direction, i.e. do not at all cross elevation contours at a right angle,

• the skeleton has cycles (around inner rocks). We will see that we can build another graph which links those same skeleton vertices/waypoints, but which do not exhibit these flaws, provided we can quantify what a 'significant deviation from the local steepest slope direction' is.

Construction of optimum branchings

In this section we define a cost function that will allow us to rate how much an edge between two vertices deviates from a slope line, and to select a suitable set of edges. The next section presents a cost function designed according to this rationale.

Definition of the cost function

Consider an elementary displacement δx = (dx, dy) from x = (x, y) to x = (x , y ) = x + δx = (x + dx, y + dy). The cost of this elementary displacement will be defined as:

C ↑ (x , x + δx) = z max -z(x) -z(x + δx) -z(x) + λ δx = ∇z δx - (∇z) • δx + λ δx (1) (2) (3)
where δx = dx 2 + dy 2 is the length of the displacement, ∇z is the gradient vector of the topographic surface at point x, and λ is a friction factor (scalar) in m • m -1 . This cost is expressed in meters of potential energy and is the sum of three terms:

(1) is the maximum elevation (relative to z(x)) we could reach with a displacement of length δx starting at x. By definition of the gradient, this maximum elevation is z max = ∇z δx , which is always positive.

(2) is the elevation we actually reached with the displacement δx, which is z = (∇z) • δx (note that z is algebraic and may be negative).

Consequently, the difference ( 1)-( 2), which is always positive whatever the displacement, is a measure of how much higher we could have gotten with a displacement of the same length. If the displacement is in the direction of the gradient (i.e. upslope) then ∇z and δx are colinear and (∇z) • δx = ∇z δx . Hence ( 1) and ( 2) cancel out in the case of a steepest ascent: the cost will be low. Conversely, if we move downslope, (∇z) • δx is negative and ( 1) -( 2) is largely positive: the displacement is very costly.

(3) is a friction term stabilizing the cost function: whatever the direction (upslope, downslope or along an elevation contour), there is always a cost λ for moving 1 meter across the surface. The effect is to smooth trajectories and to prevent paths from following too closely the smallscale irregularities on the surface.

We call C ↑ the upslope cost function, which is obviously not symmetric:

C ↑ (x, x ) = C ↑ (x , x)
The definition of C ↑ looks unusual, as classical cost functions (see e.g. Hence, every displacement of length δx around x has the same cost whatever the direction (C is necessarily isotropic). This kind of formulation is easily solved in classical GIS softwares (it only requires to specify the impedance raster, and the start and end points). However, we believe that our formulation is preferable for steepest ascent/descent problems, which are anisotropic in nature (see also [START_REF] Collischonn | A direction dependent least-costpath algorithm for roads and canals[END_REF].

Similarly, we can define a downslope cost function:

C ↓ (x , x + δx) = +(∇z) • δx - -∇z δx + λ δx (1) (2) (3)
This time, ( 2) is the minimum elevation we could reach with a displacement of length δx . The difference ( 1) -( 2) is always positive and is a measure of how much lower we could have gotten with a displacement of the same length: (1) and ( 2) cancel out in the case of a steepest descent.

Conversely, if we move upslope, ( 1) -( 2) is large and the cost is high.

Finally, these cost functions satisfy the requirement that the steepest ascent path from x to x is also the steepest descent path from x to x (assuming the same friction factor λ):

C ↑ (x, x ) = C ↓ (x , x)
In the following, we only use the upslope cost function C ↑ ; we will see why it is more convenient to treat the problem from downslope up. Moreover, in order to improve flow convergence near the glacier boundary, we slightly modify the cost function by increasing the friction factor λ near the boundary (this will penalize flowlines wandering along the boundary):

C ↑ (x , x + δx) = ∇z δx -(∇z) • δx+ λ ∞ +(λ 0 -λ ∞ )e -D(x) D λ δx λ(x)
where D(x) is the distance to the boundary at point x, λ ∞ is the friction cost far from the boundary, λ 0 > λ ∞ is the friction cost at the boundary and D λ is the scale of decrease with distance. Note that this modified cost is still anisotropic: λ(x) is an impedance, but it is not the dominant term of the

cost function. The values λ ∞ = 0.035 m • m -1 , λ 0 = 5λ ∞ = 0.175 m • m -1 ,
and D λ = 150 m were found to give very good results on all glaciers.

Finally, since we use a raster DEM (i.e. a square lattice with only 8

possible elementary displacements from a given grid point x i,j ), we define a discrete version of C ↑ :

C ↑ (x i,j , x i+δi,j+δj ) = z i,j + s max δi 2 + δj 2 -z i+δi,j+δj + λ(x i,j ) δi 2 + δj 2
where

s max = max 0<(δi 2 +δj 2 )≤2 z i+δi,j+δj -z i,j δi 2 + δj 2
is the estimated upward slope (norm of the gradient) at pixel x i,j . If the above expression for s max is negative, it is of course set to zero: there is a local maximum at x i,j .

Cost assignment for paths between waypoints

We use Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] to compute the least-cost path between each pair of waypoints. It uses the elementary cost function to produce an accumulated cost surface spreading around a start point S, and a backlink grid which gives the direction opposite to the gradient of the accumulated cost surface at each pixel. The least-cost path from S to any point is constructed in the reverse direction, starting from the destination and following the backlinks to S. In our case, the cost is defined from downslope up, so a backlink is defined from upslope down (see Figure 11). We use this algorithm to compute the costs between any pair of skeletal vertices/waypoints. In practice, we only compute the costs to the n nearest neighbors of each vertex with respect to the cost function, using a local run of Dijkstra's algorithm that is aborted as soon as n neighbor vertices have been visited (straightforward with an implementation based on a priority queue). The value n = 30 turned out to be largely sufficient for the following steps. We obtain the kind of graph shown in Figure 12. Each edge can be traversed in both directions, but at a different cost in each direction: the graph is directed.

Identification of minimum spanning branchings

The creation of a dense, redundant graph between the skeletal waypoints looks like a step backwards, compared to the simplicity of the skeleton. However, graph theory provides powerful tools which will help simplify this graph again, in the way we want.

Each waypoint (vertex) has several incoming and outgoing edges in the new directed graph. The set of flowlines we are looking for is a branching or, synonymously, a directed spanning tree i.e. a graph T such that:

• T contains no cycle,

• each vertex has one and only one incoming edge: each vertex is visited (the tree is spanning), but no two edges of T enter the same vertex (i.e., each vertex has only one downstream pixel, since we build the flowlines from downstream towards upstream).

Of course, such a tree T has to be rooted at a snout, i.e. a special vertex without an incoming edge. Moreover, we want the edges of this tree to be -as much as possible-steepest climb routes, i. Since there are several snouts in a glacier complex, the subgraph may actually be a forest i.e. a set of trees, each one rooted at a different snout.

The single-root version of the Edmonds algorithm is easily extended to the multiple-root case: a virtual 'master' root is created, and zero-cost edges are added from this root to each snout (see Figure 12).

Figure 13 shows the output of the algorithm for the Rhone-Trift glacier complex. It is worth noting that the segmentation of the complex into individual glaciers is a by-product of the method: each snout 'drains' a set of upstream vertices which form its catchment.

Conclusion and perspectives

In The method currently lacks a quality assessment, though this could be done on a small subset of glaciers by comparison with manually-extracted flowlines. The large-scale visual assessment is however very satisfying, and the resulting network can be used to compute indices such as Strahler orders, bifurcation ratios, etc. Many scaling properties of glaciers, such as volumeare scaling [START_REF] Bahr | The physical basis of glacier volume-area scaling[END_REF] or power-law behavior of accumulation basin areas [START_REF] Gsell | Identifying the probabilistic structure of drained areas as a function of hypsometry in river networks[END_REF], originate in glacier branching topology: such properties could act as a 'lever arm' for tackling the problem of catchment-scale glacier flow dynamics, much like other fundamental symmetries (plane-strain or radial), as advocated by [START_REF] Bahr | Observations and analysis of selfsimilar branching topology in glacier networks[END_REF]. We hope that this study will foster ideas in this field.

Acknowledgments

This work has been supported by an EC2CO grant from the French CNRS (INSU). The authors acknowledge the Swiss Federal Office of Topography (SwissTopo) for providing the 25-meter DEM.

Figure 1 :

 1 Figure1: Illustration of the spurious 'double talweg' in glacial landscape. This feature mainly appears in the glacier's ablation area where a narrow ice tongue is confined in a valley (C-C'). On a map, elevation contours in this area have the shape of a W with its two wedges pointing upstream, whereas contours in classical (ice-free) valleys are V-or U-shaped with a single wedge pointing upstream (A-A').

  [START_REF] Kienholz | A new method for deriving glacier centerlines applied to glaciers in Alaska and northwest Canada[END_REF] use a more complex approach based on a cost function which quantifies the trade-off between the two requirements; a set of flowlines is then extracted between glacier heads and a single snout (terminus) per glacier. Other methods apply alternate procedures in the accumulation and ablation zones[START_REF] Machguth | The length of the world's glaciers -a new approach for the global calculations of center lines[END_REF], also resulting in a large number of parameters.

Figure 2 :

 2 Figure 2: Area covered by this study. (a) Map of Switzerland with glacier contours in dark blue. (b) Zoom on the Rhone-Trift glacier complex illustrating the steps throughout the paper. This complex covers 35.8 km 2 and includes two large glaciers: the Rhone glacier (15.9 km 2 ) and the Trift glacier (16.6 km 2 ) as well as a number of smaller ones (e.g. Tiefen, Sidelen and Alpli glaciers).

Figure 3 :

 3 Figure 3: Grassfire analogy of the skeletonization: (a) propagation of the 'firefronts'; (b) resulting skeleton i.e. the set of points where two or more fronts meet.

Figure 4

 4 Figure 4 illustrates this definition. Gen(s) ⊂ ∂D denotes the generating points of the skeletal point s ∈ S(D).

Figure 4 :Figure 5 :

 45 Figure 4: Definition of skeleton S(D) as the set of centers of all maximal disks in D (gray shape). Gen(s) denotes the generating points of the skeletal point s.

Figure 6 :

 6 Figure 6: Pruning of the skeleton based on contour partitioning. Skeleton edges which have their generating points lying on a same contour segment are removed (in white). (a) Partition with k = 12 points; (b) k = 6 points; (c) k = 6 points but at different locations.

Figure 7 :

 7 Figure 7: Definition of the turn angle at a boundary point. Since contours are closed, we set c n+1 = c 1 and c 0 = c n . If boundary points are sorted clockwise (i.e. with the 'inside' on the right and the 'outside' on the left when moving along the boundary, as in the ESRI R shapefile format), then:

  i min having the lowest absolute value |K(c i min )|, thus leading to a simplified contour called DCE level n-1. The metric K is again computed for the n-1 remaining points, and the procedure is repeated p times to obtain a hierarchical set of contour partitions, called DCE level n -2, n -3, . . . , k = np.At each iteration, we remove the skeleton edges whose generating points lie on a same contour segment of DCE level k = np, and these edges are assigned the level k. Since we need at least two points on the contour to define a partition, the highest possible level is k = 2. Figure6cactually represents DCE level 6 of the leaf, and in Figure8we show the final hierarchy of skeleton edges.

Figure 8 :

 8 Figure 8: Hierarchy of skeleton edges obtained at the end of the pruning algorithm: we can stop at any level in this hierarchy.

Figure 9 :

 9 Figure 9: Definition of a snout vertex in the skeleton pruned at DCE level k (k = 238 here).

Figure 10 :

 10 Figure 10: Example of flaws in the pruned skeletal graph, notably the deviation of skeleton edges from steepest slope directions and the presence of cycles. Black dots are a downsampled subset of skeleton vertices (waypoints), and red dots are snout vertices.

5. 1 .) 2 .

 12 Rationale[START_REF] Douglas | Least-cost Path in GIS Using an Accumulated Cost Surface and Slopelines[END_REF] recalls that computing the least-cost path (with respect to various factors such as slope, soil cover, etc.) between a point A and a point B involves three steps:1. the definition of a cost of movement for an elementary movement from A to A' (close to A); 2. the computation of the accumulated cost surface spreading from the starting point, that sums the costs of all elementary movements; 3. the construction of the least-cost path as the slope line of the accumulated cost surface starting at B and ending at A[START_REF] Warntz | Transportation, social physics, and the law of refraction[END_REF]. If we want a least-cost path to look like a steepest ascent or descent line, we can figure out what the corresponding cost function has to look like by reversing the three steps: 3. The least-cost path has to look like a topographical slope line (it will never be a true one though, unless B is strictly upslope or downslope from AThis implies that the accumulated cost surface has to 'look like' the topographical surface z around A; 1. It means that the local cost of movement has to 'look like' a differential of the topographical surface, i.e. depend on its gradient ∇z.

Kienholz

  et al., 2014) are simply the product of a scalar impedance I(x) by the length of the elementary displacement: C(x , x + δx) = I(x) δx .

Figure 11 :

 11 Figure 11: Example of accumulated cost surface, spreading from a start point S. Arrows indicate the backlinks: the least-cost path from S to any destination (K, L, M, N, etc.) is constructed from the destination to the start S, following the backlink grid. White pixels are ice-free zones (inner rocks) across which no displacement is possible.

Figure 11 2 • 3 •

 1123 Figure 11 illustrates the effects of the definition of C ↑ : • destination K is located right upslope from start S, at an elevation of about 3130 m. Since edge (S,K) is almost a steepest ascent path, its cost is low: C ↑ (S, K) = 11.2 • destination L is located at about the same elevation, but one has to cross elevation contours at some angle to get there from S. The cost is higher: C ↑ (S, L) = 25.3 • destination M is located at an elevation lower than S, about 3030 m. Consequently, the edge from S to M is far from a steepest ascent (it is almost a steepest descent), and has a high cost: C ↑ (S, M) = 80.5 (note that this path would in contrary have a very low cost with respect to the downslope cost function C ↓ )• finally, destination N is located on the other side of a ridge. In order to reach N from S, one must first climb to a pass at about 3100 m (low cost), but then go down (high cost). We have C ↑ (S, N) = 169.3

Figure 12 :

 12 Figure 12: Edges and associated costs between skeleton waypoints in the Rhone-Trift glacier complex. Left: example of initial graph with numerous edges linking pairs of waypoints. For clarity, the plot displays only the edges between natural (planar) neighbors but a much denser graph can be created. Right: zoom on a region. Figures in red are the costs of each edge, in both directions.

Figure 13 :

 13 Figure 13: Output of the Chu-Liu/Edmonds algorithm with the set of edges and costs of Figure 12. The minimum branching (Directed Mininum Spanning Tree/Forest) allows to visit all waypoints, starting from a snout and always moving as upslope as possible.

Figure 14

 14 Figure 14 is a zoom on a pass across the topographical divide between Rhone and Trift glaciers (Undri Triftlimi). It explains why no edge in the DMST/DMSF should cross any major topographical ridge. Let us suppose that the ridge-crossing edge (K,L) is part of the DMSF, denoted F and of total cost C ↑ (F ). In this case, point L ultimately drains to Trift glacier's snout. However, the other incoming edge to L, (M,L), has a lower cost than (K,L): hence if we remove (K,L) from F and add (M,L), we obtain a forest F that is still spanning (L is visited i.e. has an incoming edge) and has a lower total cost C ↑ (F ) = C ↑ (F ) + C ↑ (M, L) -C ↑ (K, L) < C ↑ (F ). This means that F was not a DMSF in the first place (it is spanning but not of minimum cost): point L has to flow to Rhone glacier's snout, and no edge should go through the pass. Hence, our procedure is efficient not only for flowline extraction, but also for glacier segmentation.

Figure 14 :

 14 Figure 14: Zoom on the minimum branching in the vicinity of Undri Triftlimi, a pass on the divide between Rhone and Trift glaciers. Vertex L, which is located south of the ridge, has to be visited from M and not from K in the Directed Minimum Spanning Tree/Forest: no edge should go across any major divide.

  this paper we proposed a new method to extract glacier flowlines, based on Voronoi skeletonization of glacier boundaries, skeleton pruning using Discrete Curve Evolution (DCE), and the construction of a Directed Minimum Spanning Tree between skeletal vertices with respect to an anisotropic, upslope cost function C ↑ . The application of the method requires limited parameter tweaking: it only requires a selection rule k = f (A) for the level k of the DCE as a function of glacier area A (2 parameters), and 3 parameters for the friction law (i.e.the isotropic term of the cost function). Table2sums up the values used on the whole domain of Figure2a(1200 km 2 of glaciers, the largest contiguous icefield having an area of 261 km 2 ).DCE level selectionk = 2 + k 0 A γ k 0 13.5 γ 0.8 Cost function C ↑ (x , x + δx) = ∇z δx -(∇z) • δx + λ(x) δx with λ(x) = λ ∞ +(λ 0 -λ ∞ )e -D(x) D λ λ ∞ 0.035 m • m -1 λ 0 0.175 m • m -1 D λ 150 m

Table 1 :

 1 shows that only one candidate satisfies the third re-Check for requirement (iii) in the definition of snout vertices.

	quirement:
		Input site in DCE level k
	Leaf vertex that separates the vertex's Is requirement (iii) met ?
		generating pair
	S 1	I 3
	S 6	I 5	NO
	S 9	I 1	NO

NO: z(I 2 ) > z(I 3 ) > z(I 4 ) S 2 I 4 YES: z(I 3 ) > z(I 4 ) < z(I 5 ) S 4 I 2 NO: z(I 1 ) > z(I 2 ) > z(I 3 )

Table 2 :

 2 Parameters of the method.


Mohamed Hadded
email: mohamed.hadded@it-sudparis.eu

Paul Muhlethaler
email: paul.muhlethaler@inria.fr

Rachid Zagrouba
email: rachid.zagrouba@cristal.rnu.tn

Anis Laouiti
email: anis.laouiti@it-sudparis.eu

Leila Azouz Saidane
email: leila.saidane@ensi.rnu.tn

Using Road IDs to Enhance Clustering in Vehicular Ad hoc Networks

Keywords: VANET, MANET, Clustering, GPS, MAP, Road ID

Vehicular ad hoc networks (VANETs) where vehicles act as mobile nodes is an instance of Mobile Ad hoc NETworks (MANETs), which are essentially developed for intelligent transportation systems. A challenging problem when designing communication protocols in VANETs is coping with high vehicle mobility, which causes frequent changes in the network topology and leads to frequent breaks in communication. The clustering technique is being developed to reduce the impact of mobility between neighboring vehicles. In this paper, we propose an Adaptive Weighted Cluster Protocol for VANETs, which is a road map dependent and uses road IDs and movement direction in order to make the clusters structure as stable as possible. The experimental results reveal that AWCP outperforms four other most commonly used clustering protocols in terms of control packet overhead, the packet delivery ratio, and the average cluster lifetime, which are the most usual metrics used for comparing performance.

I. INTRODUCTION

Research on Vehicular Ad hoc Networks (VANETs) has attracted increasing interest over recent years due to its capability to improve road safety by using Vehicle To Vehicle (V2V) and/or Vehicle To Infrastructure (V2I) communications [START_REF] Tse | Improving message reception in vanets[END_REF]. It can also be used to improve traffic management conditions and to provide on-board infotainment such as Internet access, video streaming, etc. VANETs are characterized by the self-organization of the nodes, where their nodes can be vehicles, roadside units or sensors. Due to the varying vehicular densities caused by high vehicle mobility, supporting network connection requires a high communication overhead for exchanging and updating topology information. For instance, in a fully distributed VANET, each vehicle is required to maintain its own connectivity to its one hop neighboring vehicles. Without using expensive components such as central points, establishing a hierarchical clustering structure within the network can reduce the relative mobility between neighboring vehicles, thereby reducing communication overhead [START_REF] Rawashdeh | A novel algorithm to form stable clusters in vehicular ad hoc networks on highways[END_REF]. Clustering allows the formation of organized groups and is used to coordinate channel access, simplify routing, and security. In VANETs, each vehicle is equipped with a digital road map and a positioning system, e.g. Garmin Nuvi 50 GPS that allows it to obtain the time, its speed and position [START_REF] Carswell | A modular neural network approach to improve map-matched gps positioning[END_REF] and the ID of the road on which it is traveling. Therefore a clustering protocol for VANET network can be designed using a MAP digital and GPS based approach. However, the main challenge for clustering protocols in VANETs is ensuring topology stability by taking into account many mobility parameters. This should be performed with a minimum overhead. In this paper, we propose a multi-metrics based Adaptive Weighted Clustering Protocol (AWCP) that takes advantage of the geographic information of vehicles. The main contributions of the paper are listed below:

• We identify certain essential features that the clustering protocols must satisfy to build stable clusters in VANETs.

• We propose a Map and GPS based clustering protocol using the WCA algorithm [START_REF] Chatterjee | A weighted clustering algorithm for mobile ad hoc networks[END_REF] for VANET in which a vehicle only considers neighbors moving on the same road and in the same direction. To the best of our knowledge, this is the first attempt to design a clustering protocol by taking advantage ID of the road on which the vehicles are traveling.

• We evaluate the proposed protocol on realistic VANET mobility scenarios taken from the metropolitan area of San Jose (California) and we validate its performance by comparing it with other clustering protocols in the literature.

The rest of the paper is organized as follows. In Section 2, we present some relevant guidelines (in form of rules) to establish stable clusters in VANETs and a summary of related work. Section 3 describes our clustering protocol, called AWCP. Section 4 presents the simulation results and the performance evaluation. Finally, conclusions and future work are reported in Section 5.

II. CONTEXT AND RELATED WORK

A. Clustering in VANET

When a vehicle node wishes to participate in a cluster head election, it firstly collects all the necessary one-hop neighbors information. In this section, we identify the rules that the clustering protocol must satisfy in VANET networks with the aim to form stable clusters, where re-clustering is reduced, and cluster members lifetimes are prolonged. Rule 1. The cluster head vehicles should have sufficiently powerful radios to be able to communicate with the members of their clusters. This implies that the cluster heads should be close to the center of the cluster. Thus, the vehicle that has the minimum average Euclidean distance to their direct neighbors can be elected to act as the cluster head. The average Euclidean distance is:

δ(i, t) = ∑ j √ (xi-xj ) 2 +(yi-yj ) 2 n(i,t) (1)
n(i, t) is the number of vehicles connected directly to i at instant t, where j is any vehicle that is connected to i. Rule 2. To form a stable cluster, the cluster head should have similar mobility characteristics as the vehicles within its cluster. Indeed, if cluster heads are elected without taking speed into account, the number of vehicles that will quickly move out of communication range from their cluster head will increase. Thus, the vehicle whose current speed is the closet to the mean value is elected as the cluster head. The average value of speed is:

ρ(i, t) = ∑ j ν(j,t) n(i,t) (2)
ν(t, j) is the speed of vehicle j at instant t, where j is any vehicle that is in communication range of the cluster head candidate i. Rule 3. The formation of a high number of clusters increases the overhead and the inter-cluster interference and degrades the network performance. To overcome this, the cluster mechanism should group all vehicles in the network with a small number of cluster heads. Thus, we restrict the vehicle that has the maximum number of neighboring vehicles should be elected to act as the cluster head. However, this number i should be bounded by a capacity value β i which represents the maximum number of neighbors that vehicle i can optimally handle as a cluster head. Rule 4. To provide a stable cluster structure, the cluster members should always remain within the transmission range of their cluster heads. Thus, the clustering protocol should take into consideration the mobility features of VANETs (e.g. multiple roads including road junctions, opposite-direction of the traffic flow). The Mobility Direction (MD) is necessary information that can be used to form stable clusters. Indeed, the vehicles that are moving in different directions cannot remain within transmission range for a long period of time. Let us consider the VANET scenario shown in Figure 1. Vehicle v i can choose to join one of the two clusters y or n which are managed by the two vehicles v y and v n , respectively. If vehicle v i joins cluster y, it will quickly leave the cluster and it will need to choose another cluster head. However, if vehicle v i joins cluster n, they will both be moving in the left hand direction and as a result it will remain a member of this cluster for a longer period of time.

Rule 4.2. The mobility direction is not always sufficient to insure clustering stability in VANETs. As shown in Figure 1, based on the mobility direction metric the vehicle v j can join cluster k managed by the vehicle v k because they are moving on the same direction. However, they are not moving on the same road and vehicle v k may need to change its status and choose a new cluster head if the distance between the two roads becomes greater than the communication range of the cluster head. Moreover, based on the definitions of the mobility directions given in [START_REF] Omar | Vemac: A tdma-based mac protocol for reliable broadcast in vanets[END_REF], the two vehicles v i and v j are considered to be moving in the same "left" direction and thus these vehicles can be grouped together to form a cluster. Since the two vehicles are not moving on the same road, vehicle v i will be out of the cluster j after a short period of time and it will need to join a new cluster. Thus, the Road ID (RID) is critical information to provide a more stable cluster structure and to reduce the average number of times a vehicle must change clusters. In this paper, we impose that each vehicle only considers neighboring vehicles that are moving on the same road and in the same direction, and ignores control messages from vehicles on a different road and moving in the opposite direction.

B. Related work

Clustering is the best known method in VANETs to enable efficient resource allocation with low overhead and to reduce the relative mobility between neighboring vehicles. Several VANET research studies in the literature have focused on developing clustering protocols, most of which are based on MANET clustering techniques. However, none of the protocols proposed takes RID into account when forming clusters in VANETs. As a result, these protocols do not create stable clusters. Some of the most significant protocols are described below.

In [START_REF] Rawashdeh | A novel algorithm to form stable clusters in vehicular ad hoc networks on highways[END_REF], the authors proposed a multi-metric algorithm for cluster head elections, called Threshold-based Technique (TB), suitable for highway area. In addition to the position and the direction, this algorithm uses a speed difference metric as a new parameter to increase the cluster lifetime. The vehicles that are moving at high speed are regrouped into one cluster, while the vehicles moving at low speed are grouped into another cluster. An Adaptable Mobility-Aware Clustering Algorithm based on Destination positions, called AMACAD, is proposed and evaluated by Morales et al. [START_REF] Morales | An adaptable mobilityaware clusteringalgorithm in vehicular networks[END_REF]. The goal of this work has been to develop a clustering protocol with an efficient message exchange mechanism, which improves the clustering stability in VANETs. AMACAD performs clustering based upon information such as current location, vehicle velocity, relative destination and final destination of vehicles. A Multi-Head Clustering Algorithm, called Center-Position and Mobility (CPM), was proposed in [START_REF] Lo | A multi-head clustering algorithm in vehicular ad hoc networks[END_REF]. This technique aims to create stable clusters and reduce re-clustering overhead by supporting single and multiple cluster heads. In the cluster head election phase, vehicles within communication are organized into clusters and one vehicle for each cluster is elected to act as a Master Cluster Head (MCH). Then, some cluster members in the cluster are selected to be Slave Cluster Heads (SCHs). In order to form stable clusters, the authors imposed that all the vehicles in a cluster are moving in the same direction.

Several other clustering algorithms designed for MANETs also work in VANETs and are frequently used for comparison with other VANET clustering protocols. For instance, the Lowest-ID clustering algorithm (LID) [START_REF] Gerla | Multicluster, mobile, multimedia radio network[END_REF] is based on electing a node with the smallest ID as a cluster head. The Highest-Degree algorithm (HD) [START_REF] Ramalingam | Associativity based cluster formation and cluster management in ad hoc networks[END_REF] selects a node as a cluster head based on the nodes' connectivity. The node with the maximum number of neighbors becomes the cluster head. The Weighted Clustering Algorithm (WCA) [START_REF] Chatterjee | A weighted clustering algorithm for mobile ad hoc networks[END_REF] elects a node to act as a cluster head based on a combined weight which includes the number of its neighbors, their average distance and the node's average speed, and battery-life. MOBIC [START_REF] Basu | A mobility based metric for clustering in mobile ad hoc networks[END_REF] is a Mobilitybased algorithm designed for MANETs which is also used in VANETs. MOBIC is a mobility based version of the Lowest-ID algorithm and uses a signal power level metric to elect cluster heads. Table I compares the features of these protocols. It is clear from this table that while all of these protocols satisfy certain of these rules, none of them satisfies all of them. As result, these protocols may not work efficiently in VANET network. Thus, in this paper, we propose a clustering protocol, called AWCP, that uses road ID information in addition to the known parameters i.e., direction, connectivity, average speed. To the best of our knowledge, this is the first contribution to design a clustering protocol that uses road ID information and which satisfies all the above mentioned rules.

III. AWCP DESCRIPTION

A. System Model

Our algorithm is based on the assumption that each vehicle in a VANET can know its road ID via a digital road map and a positioning system, e.g. GPS (Global Positioning System) or a GALLILEO receiver that also allows it to obtain an accurate real-time three-dimensional geographic position (latitude, longitude and altitude), direction, speed and exact time.

B. Cluster Head Election

In this section, we present the cluster setup and maintenance mechanisms of AWCP in detail. Initially, all vehicles are in the Undecided State (US). To divide the network into clusters, each active vehicle changes its state to Cluster Head Candidate (CHC) and it starts to broadcast a HELLO message periodically containing all the necessary information ⟨V ID, RID, M D, position, speed⟩ to its One-Hop neighbors (OH). In order to form stable clusters, each vehicle uses RID and MD to filter out any vehicle that is moving on another road or in the opposite direction. Upon reception of a HELLO message from all its one-hop neighbors, each vehicle i calculates its current weight W (i, t) using the following function [START_REF] Carswell | A modular neural network approach to improve map-matched gps positioning[END_REF]. The weight function consists of three parts, i.e., average-distance weight factor (1), average-speed weight factor (2), number of neighbors weight factor.

W (i, t) = w 1 * δ(i, t)+w 2 * |ν(i, t)-ρ(i, t)|-w 3 * n(i, t) (3)
Where w 1 , w 2 and w 3 are the balancing factors such that ∑ 3 k=1 w k = 1. We note that the three weight factors are in conflict. For simplicity, we assume that all the factors should be minimized. In fact, the multiplication of the thirst weight factor by (-1) allows us to transform a maximization to minimization. Then, each node i broadcasts a beacon message containing all the necessary information for the CH election algorithm ⟨V ID, RID, M D, W, CH -ID⟩. Vehicle i announces itself as a CH by assigning its own ID to the CH-ID field of the election beacon. When a vehicle i receives beacons from its one-hop neighbors, it sorts its neighbor list OH i according to the weights received in the beacons, and then it executes the cluster head election algorithm to change its status from CH to Cluster Member (CM), Cluster Gateway (CG) or remain CH.

W (i, t) = {min W (j, t) ∀ j ∈ OH i } (4), n(i, t) ≤ β i (5)
The vehicle i that satisfies the two properties (4) and ( 5) at instant t is elected as the CH. Then, all vehicles that are within transmission range of the CH become CMs or CGs and are not allowed to participate in another cluster head election procedure. The CH election algorithm terminates once all the vehicles either become a CH, CM or a CG. Algorithm 1 outlines the details of the CH nodes' election. It is executed by each vehicle i having at least one neighboring vehicle. In Algorithm 1, i, j, and x represent three vehicles which are moving in the same road and on the same direction and are participating in the CH election process. In addition, ITJ Interval is the time interval for a CH vehicle to broadcast the Invite-To-Join (ITJ) message, PRE Interval is the time interval for a CM to signal its presence to its CH, while CH Timeout Interval is the time interval for a vehicle to elect itself as a CH, if it did not receive any ITJ messages during this period.

C. Cluster maintenance

In VANETs, a vehicle can join or leave a cluster at any time. These two operations will have only local effects on the topology of the cluster if the vehicle is a CM. However, if the vehicle is the CH, it must hand over the responsibility to one of the very close cluster members before leaving the cluster. The first reason for that is to maintain the cluster structure even if the current CH leaves. The second reason is to avoid using the re-clustering algorithm and thus no re-clustering overhead is generated when the CH leaves the cluster. Then, the current CH will order the CM to switch to CH and switch its own state to CM. Each cluster head periodically broadcasts an ITJ messages to its one-hop neighbors. Once a US or CHC vehicle receives an ITJ message, and if it wishes to join the cluster, it will check the received signal strength. The US or CHC vehicle will consider the ITJ message to be valid if its signal strength is greater than the predefined threshold denoted by Pr Threshold. When receiving a valid ITJ message, the vehicle sends a Request-To-Join (RTJ) message including the vehicle's ID, road ID and direction. When the CH receives the RTJ message, it checks the road ID on which the requesting vehicle moving and, if it is moving in the same direction, the CH sends an acknowledgment (ACK) including its ID number. After the reception of the ACK, the corresponding vehicle becomes a CM of this cluster. Once a US vehicle becomes a CM, it is not allowed to participate in another cluster head election procedure. Moreover, if a CM receives an ITJ message from another neighboring CH moving on the same road and in the same direction, the vehicle will switch from the CM state to the CG state. Figure 2 shows the vehicle state transitions diagram. A vehicle remains in the CM state as long as it receives an ITJ message from its CH every ITJ Interval. When the CM vehicle does not receive an ITJ message from its CH during CH Timeout Interval, it considers that it has lost contact with the CH and thus switches its state to CHC. Each CH updates a time stamp field for each CM based on the presence messages (PRE-MSG) received. The CH removes a CM from its cluster members list if the difference between the current time and the last time stamp of the PRE-MSG message received from it is greater than CM Timeout Interval. The CH will change its state to CHC, if its list of cluster members is empty. Moreover, when two or three CHs moving on the same road and in the same direction receive an ITJ message from each other, only one of them will keep its CH responsibility based on the the weight function (3), while the others will switch to a CM.

IV. SIMULATION RESULTS AND PERFORMANCE EVALUATION

This section presents NS2 simulation results to evaluate the performance of AWCP. The objectives of the evaluation are to: (1) Evaluate the performance of AWCP under realistic and different mobility scenarios, (2) Evaluate and compare the efficiency of AWCP with the most popular clustering protocols in the literature, (3) Test the efficiency of AWCP in reducing the number of states changed per vehicle.

A. Simulation Scenarios

We generated a realistic VANET environment by selecting a real highway area from a digital map which took into account road directions, road intersections, and traffic rules. Figure 3 shows a metropolitan area from the Map of San Jose (California) of size 4000m×4000m exported from Open-StreetMap (OSM) and edited using Java OpenStreetMap Editor (JOSM). Then MOVE and SUMO [START_REF] Karnadi | Rapid generation of realistic mobility models for vanet[END_REF] were used respectively to generate vehicle traffic scenarios and to simulate the area with vehicular traffic. To generate vehicular traffic by MOVE and SUMO, we defined a vehicle flow which described a swarm of vehicles in each direction. The parameters of each vehicle flow consisted of the maximum number of vehicles, the starting road and destination of the flow, the time to start three of the most widely used Quality of Service (QoS) metrics in this domain [START_REF] Rawashdeh | A novel algorithm to form stable clusters in vehicular ad hoc networks on highways[END_REF]: The Average Cluster Lifetime (ACL), which is the average time period from the moment when a vehicle becomes a CH, CM or CG to the time when it changes its state. The Control Packet Overhead (CPO), which is the rate of AWCP control packets used to form and maintain the cluster structures. Finally, the Packet Delivery Ratio (PDR), which is the ratio of the number of data packets that are correctly delivered to their destinations.

B. AWCP parameters

We have used the min-EUDT AWCP configuration found by NSGA-II in [START_REF] Hadded | A multi-objective genetic algorithm-based adaptive weighted clustering protocol in vanet[END_REF] which is the most balanced setting of parameters, that optimizes the AWCP QoS metrics. Table I shows the parameters of AWCP and their values. These parameters are four timers, four counters and three weighting factors.

C. AWCP performance evaluations

In this section, we evaluate and compare the performance of AWCP with other well known clustering protocols proposed in the literature, namely WCA, LID, HD and CPM. Figure 4 shows the ACL of all the algorithms for different VANET scenarios. This figure shows that the ACL is increased by respectively 63.6%, 62.2%, 59% and 45.5% on average when AWCP is used compared to the WCA, HD, LID, and CPM protocols. Therefore, we can conclude that the protocol proposed provides stable clusters which have a long lifetime. As AWCP takes into account road IDs and movement directions to form clusters, the CMs will be associated with their CHs for a longer period of time. We can also note that CPM performs better than WCA, LID and HD, because the CPM protocol forms clusters based on the mobility direction.

The Clustering Protocol Overhead (CPO) of AWCP and the other protocols for various of VANETs scenarios is shown in Figure 5. It is clear from this figure that our protocol has a lower overhead than the other protocols. In fact, AWCP reduces the CPO by respectively 38.8%, 37.2%, 37.1% and 47.2% on average compared to WCA, HD, LID and CPM. There are two reasons why AWCP decreases the overhead. Firstly because the maximization of the cluster heads' duration and the cluster members' duration decreases the number of control messages required to elect new cluster heads and to join a new cluster, respectively. Secondly, the minimization of the number of the clusters reduces the amount of ITJ messages broadcasted by the cluster heads. It clearly shows that for varying traffic densities the AWCP protocol gives the best performance in terms of PDR, except for the High scenario, where WCA, HD and LID deliver a higher data rate. Although the network performance in terms of throughput significantly decreases when the vehicle density increases, on average, AWCP guarantees a better PDR than the other protocols. This is due to the fact that AWCP does not generate an excessive clustering overhead and thus the data packets are transmitted to their destination vehicles with a lower collision rate. Fig. 6. PDR results for different VANET scenarios Figure 7 shows the average US duration (the average duration in which the vehicles are in the US state) with respect to road traffic density. We note that the AWCP protocol provides a smaller US average duration than the other protocols, except for Low and High scenarios where CPM and ACWP behave similarly. Moreover, it can be seen from this figure that the average US duration increases significantly for HD, WCA and LID in the high scenario, while it still remains reasonable for both the AWCP and CPM protocols. Figures 8 and9 show Fig. 7. The average duration in which the vehicles are in the US state for different VANET scenarios the number of changes of states for each vehicle during the simulation time for the S3 and S4 scenarios. In both figures 8 and 9, we note that AWCP causes the lowest number of transitions. For instance, vehicle 14 in Figure 8 kept its sate throughout the simulation time when AWCP was used, while it changes its state 3, 4, 5 and 8 times when CPM, WCA, LID and HD were used, respectively. These results can be explained by the fact that AWCP avoids the problem of merging multiple clusters into a single cluster at road junctions. V. CONCLUSION Designing of efficient clustering protocols is an important issue in VANETs due to the rapid changes in network topology and the lack of infrastructure. In this paper, we identify the essential properties that clustering protocols must satisfy to form stable clusters in VANETs, and we propose an optimized clustering protocol, called AWCP, whose objective is to maximize the lifetime of the cluster heads and cluster members. AWCP is a MAP and GPS based approach which takes advantage of knowing the road ID and the direction in which the vehicles are traveling. Our method of selecting cluster heads based on mobility features and a weight function is the key to achieving a more stable cluster. The experimental analysis shows that AWCP clearly improves the clustering performance in VANETs in terms of cluster lifetime duration and communication overhead compared to the well known clustering protocols WCA, LID, HD. This is also the case when AWCP is compared with another protocol CPM, which considers both the vehicles position and the direction of movement.

The ns2 simulation results presented in this paper are very encouraging and justify continuing in this research direction. We plan to design a cross-layer architecture (MAC/AWCP) to improve channel access efficiency in VANETs in which each cluster head is responsible for assigning bandwidth to all the members of its cluster.

Fig. 1 .

 1 Fig. 1. Road-based clustering Rule 4.1. The Mobility Direction (MD) is necessary information that can be used to form stable clusters. Indeed, the

Fig. 2 .

 2 Fig. 2. Vehicle state transition diagram

Fig. 3 .

 3 Fig. 3. San Jose (California) urban area captured from Google Maps (left) and exported to a VANET network topology by using MOVE/SUMO (right)

Fig. 4 .

 4 Fig. 4. ACL results for different VANET scenarios

Figure 6

 6 Figure 6 shows the Packet Delivery Ratio (PDR) achieved by each clustering protocol for various VANET scenarios.It clearly shows that for varying traffic densities the AWCP protocol gives the best performance in terms of PDR, except for the High scenario, where WCA, HD and LID deliver a higher data rate. Although the network performance in terms of throughput significantly decreases when the vehicle density increases, on average, AWCP guarantees a better PDR than the other protocols. This is due to the fact that AWCP does not generate an excessive clustering overhead and thus the data packets are transmitted to their destination vehicles with a lower collision rate.

Fig. 8 .

 8 Fig. 8. The number of vehicles' state transitions for the S3 VANET scenario

TABLE I .

 I Comparison of clustering protocols

	Protocols	WCA	HD	LID	CPM	AMACAD	TB	MOBIC
	Rule 1
	Rule 2
	Rule 3
	Rule 4.1
	Rule 4.2

TABLE II .

 II VANETs scenarios

	Scenario	Number of vehicles	Number of CBR sources
	Low (S1)	25	5
	Medium (S2)	50	15
	High (S3)	100	25
	Very High (S4)	150	35

TABLE III .

 III Simulation parameters in ns-2

	Parameter	Value/Protocol
	Simulation time	100 s
	Vehicle speed	120 -150 km/h
	Propagation model	T wo Ray Ground
	Medium Capacity	6 M bps
	PHY/MAC Layer	IEEE 802.11p
	Transmission range	1000 m
	Transport Layer	U DP
	CBR Packet Size	512 bytes
	CBR Data Rate	0.1M bps

TABLE IV .

 IV AWCP PARAMETERS

	Parameter	Type	Value
	Hello Interval	R	0.78s
	Election Interval	R	0.16s
	ITJ Interval	R	7.23s
	PRE Interval	R	9.16s
	CH Timeout Interval	R	12.75s
	CM Timeout Interval	R	12.7s
	Cluster Size	Z	50
	Pr Threshold	R	7.23E-16w
	Distance Weight factor (w1)	R	0.716
	Speed Weight factor (w2)	R	0.204
	Neig Weight factor (w3)	R	0.07