Hector Garcia De Marina
email: hgdemarina@ieee.org.

Zhiyong Sun
email: zhiyong.sun@anu.edu.au.

Murat Bronz

Gautier Hattenberger

H Garcia De Marina

Circular formation control of fixed-wing UAVs with constant speeds

In this paper we propose an algorithm for stabilizing circular formations of fixed-wing UAVs with constant speeds. The algorithm is based on the idea of tracking circles with different radii in order to control the inter-vehicle phases with respect to a target circumference. We prove that the desired equilibrium is exponentially stable and thanks to the guidance vector field that guides the vehicles, the algorithm can be extended to other closed trajectories. One of the main advantages of this approach is that the algorithm guarantees the confinement of the team in a specific area, even when communications or sensing among vehicles are lost. We show the effectiveness of the algorithm with an actual formation flight of three aircraft. The algorithm is ready to use for the general public in the open-source Paparazzi autopilot.

I. INTRODUCTION

There is a growing body of literature that recognises the importance of Unmanned Aerial Vehicles (UAVs) for aerial surveillance, search and rescue, patrol, inspection of facilities, precision agriculture, and atmosphere study [START_REF] Beard | Decentralized cooperative aerial surveillance using fixed-wing miniature UAVs[END_REF], [START_REF] Goodrich | Supporting wilderness search and rescue using a camera-equipped mini UAV[END_REF], [START_REF] Tokekar | Sensor planning for a symbiotic UAV and UGV system for precision agriculture[END_REF]. The current trend in multi-agent systems [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF] has led to a proliferation of works where individual and independent UAVs have gradually been replaced by a team of cooperative and coordinated ones [START_REF] Ryan | An overview of emerging results in cooperative UAV control[END_REF], [START_REF] Bertuccelli | Robust planning for coupled cooperative UAV missions[END_REF], [START_REF] Wang | Cooperative UAV formation flying with obstacle/collision avoidance[END_REF], [START_REF] Stipanović | Decentralized overlapping control of a formation of unmanned aerial vehicles[END_REF]. Formation control is one of the most well-known tools for assessing the cooperation and coordination of multi-agent systems [START_REF] Oh | A survey of multi-agent formation control[END_REF], [START_REF] Anderson | UAV formation control: Theory and application[END_REF], where the usage of such a tool aims at forming and maintaining a prescribed geometrical shape for a group of vehicles.

An important challenge in the formation control of UAVs is to consider demanding constraints for the vehicles. For example, fixed-wing UAVs are preferred to rotorcrafts in missions where the endurance or the capacity for traveling long distances are essential requirements. However, rather than point mass models [START_REF] Montenbruck | Fekete points, formation control, and the balancing problem[END_REF], [START_REF] Lee | Distributed formation control via global orientation estimation[END_REF], guidance algorithms for fixed-wing UAVs have to consider nonholonomic constraints, e.g., by modeling the aircraft as unicycles. Theoretical contributions on the coordination and formation control of unicycles include the consensus and rendezvous [START_REF] Lin | Necessary and sufficient graphical conditions for formation control of unicycles[END_REF], [START_REF] Dimarogonas | On the rendezvous problem for multiple nonholonomic agents[END_REF], and the circular formations [START_REF] Marshall | Pursuit formations of unicycles[END_REF], [START_REF] El-Hawwary | Distributed circular formation stabilization for dynamic unicycles[END_REF]. The circular or pursuit formation design is a practical method for steering the vehicles to coordinated and periodic trajectories, which provide means to sample data with a desired spatial and temporal separation [START_REF] Leonard | Coordinated control of an underwater glider fleet in an adaptive ocean sampling field experiment in monterey bay[END_REF]. However, all the above mentioned works consider that both the speed and orientation of the vehicles can be actuated by the formation control algorithm in order to accomplish the mission's goal.

Fixed-wing aircraft are usually designed for flying most efficiently at a given fixed airspeed [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]. Furthermore, an aircraft needs to fly with an airspeed over a certain lower bound or otherwise the aircraft stalls and falls down. Consequently, the control problem to be discussed in this paper considers unicycle-like vehicles with constant speed, i.e., we only actuate on the steering of the vehicle via coordinated turns by actuating on the bank angle of the aircraft. Note that having a constant airspeed does not imply to have a constant ground-speed because of the effect of the wind. Therefore, the wind causes the aircraft to travel with different ground speeds depending on its yaw and heading angles with respect to some frame of coordinates fixed on the ground, e.g., at center of the circular formation. Nevertheless, in this paper it will be assumed that the speed of the wind is much smaller than the desired airspeed, so the ground-speed can be considered almost constant during the vehicle's mission.

Circular formations for unicycles with the constraint of having constant speeds make the formation control problem more challenging as it has been shown in the early work in [START_REF] Justh | Equilibria and steering laws for planar formations[END_REF] for just two vehicles. In particular, the analysis of constantly moving vehicles is problematic if one wants to control geometrical relations between them. However, a clever strategy consisting in controlling the instantaneous center of rotation of the vehicle, instead of its position, for a given fixed angular velocity has been employed in different works [START_REF] Sepulchre | Stabilization of planar collective motion: All-to-all communication[END_REF], [START_REF]Stabilization of planar collective motion with limited communication[END_REF], [START_REF] Seyboth | Collective circular motion of unicycle type vehicles with non-identical constant velocities[END_REF], [START_REF] Briñon-Arranz | Cooperative control design for time-varying formations of multi-agent systems[END_REF], [START_REF] Sun | Collective control of multiple unicycle agents with non-identical constant speeds: Tracking control and performance limitation[END_REF]. The benefit of this approach is that such instantaneous center of rotation for a unit speed vehicle can be fixed with respect to some global frame of coordinates while the vehicle is still moving, i.e., it just circulates around a constant point. The cited works applying this approach include the rendezvous of such points or the control of geometric relations between them such as position-or distance-based formation control approaches [START_REF] Oh | A survey of multi-agent formation control[END_REF]. However, there are certain drawbacks associated with the control of such centers of rotation. For example, it is common that during the mission the steering control action of a vehicle could be close to or exactly zero. For such a situation, if communications or lines of sight for the sensors are lost, then the control action of the vehicle is not updated since the steering depends only on the relative states with respect to its neighbors, and consequently, the aircraft will continue flying straight. This is clearly a problem since nowadays the restrictions for drones in the airspace are tight in many countries. Therefore, the flight plan has to guarantee that the fixed-wing UAVs will not abandon their designated flying area.

In this paper we propose a distributed algorithm for controlling circular formations of fixed-wing UAVs with constant speeds, where each vehicle has a (feasible) prescribed inter-position with respect to its neighbors on the circle. In particular, the formation control algorithm does not directly actuate on steering the vehicle but setting the radius of the circle to be tracked by the vehicle, i.e., we actuate on the angular velocity of the vehicle around the center point of the circle. This approach can be related to algorithms where the agents are exclusively confined on the target circumference and they can change their phases [START_REF] Wang | Forming circle formations of anonymous mobile agents with order preservation[END_REF], [START_REF] Bullo | Distributed control of robotic networks: a mathematical approach to motion coordination algorithms[END_REF]. While several works controlling the inter-vehicle phases have considered constraints such as anonymity, restrictions in the communications or order preservation, to the best of the authors' knowledge, only the work by Wang et al [START_REF] Wang | Controlling anonymous mobile agents with unidirectional locomotion to form formations on a circle[END_REF] addresses the constraint of having vehicles that cannot move backwards. However, such condition implies that the vehicles can be stopped on the circumference, something impossible for a fixed-wing aircraft.

The proposed algorithm in this work has a series of advantages and features. First, it is distributed, i.e., the vehicles depend only on relative measurements, such as relative positions, with respect to their neighbors, and a complete graph is not necessary. Second, unlike many of the above cited works, the desired inter-vehicle angles on the circle can be prescribed. Third, the desired formation is exponentially stable. Fourth, an arbitrary bound of the maximum distance of the vehicles with respect to the center of the circle can be chosen by design, therefore it is guaranteed the confinement of the formation regardless of broken communications or sensing. Fifth, it is possible to extend the algorithm to not only circles, but at least to any smooth closed orbit homeomorphic to a circle.

The remaining parts of the paper are organized as follows. In Section II, the notation, the considered model of the fixedwing aircraft and the employed trajectory tracking algorithm for following a circle are introduced. In Section III we state the circular formation problem and propose the design of a controller as a solution together with a stability analysis. Experimental results with actual aircraft are presented in Section IV. The algorithm has been implemented in the popular open-source autopilot Paparazzi [START_REF] Hattenberger | Using the paparazzi UAV system for scientific research[END_REF] and it is ready to be used by the general public.

II. PRELIMINARIES

A. Notation

Consider a formation of n ≥ 2 fixed-wing UAVs whose positions are denoted by p i ∈ R 2 with i ∈ {1, . . . , n}. The vehicles are able to sense the relative positions with respect to their neighbors. The neighbors' relationships are described by an undirected graph G = (V, E) with the vertex set V = {1, . . . , n} and the ordered edge set E ⊆ V × V. The set N i of the neighbors of vehicle i is defined by N i ∆ = {j ∈ V : (i, j) ∈ E}. Two vertices are adjacent if (i, j) ∈ E. A path from a vertex i to a vertex j is a sequence starting at i and ending at j such that two consecutive vertices are adjacent, and if i = j then the path is called a cycle. We assume that the graph G is connected, i.e., there is a path between any two vertices i and j. We define the elements of the incidence matrix B ∈ R |V|×|E| , where |X | denotes the cardinality of the set X , for G by

b ik ∆ =      +1 if i = E tail k -1 if i = E head k 0 otherwise , (1)
where E tail k and E head k denote the tail and head nodes, respectively, of the edge E k , i.e., E k = (E tail k , E head k ). A circular trajectory with radius r ∈ R + can be described by the following implicit equation

C r ∆ = {p : ϕ(p) = 0}, (2)
where ϕ(p) = p 2 x +p 2 y -r 2 and p = p x p y T is a Cartesian position with respect to a frame of coordinates whose origin is at the center of C r . The plane R 2 can be covered by the following disjoint sets ϕ c (p)

∆ = ϕ(p) = c ∈ R,
where each level set is defined for a particular value of c such that the resulting circle's radius is non-negative, and in particular, the zero level set corresponds uniquely to C r . We define by n(p) ∆ = ∇ϕ(p) the normal vector to the circle corresponding to the level set ϕ(p) and the tangent vector at the same point p is given by the rotation

τ (p) = En(p), E = 0 1 -1 0 .
Note that C r belongs to the C 2 space and it is regular everywhere excepting at its center, i.e.,

∇ϕ(p) = 0 ⇐⇒ p = 0, (3)
and all the level sets ϕ c (p) can be parametrized. In particular, the vehicle i can calculate such parametrization associated to its position with the following expression

θ i (p) = atan2(p y , p x ) ∈ (-π, π]. (4)
Note that θ i (p) belongs to the circle group S1 .

B. Fixed-wing aircraft's model

Consider for the unit speed i'th fixed-wing aircraft the following nonholonomic model in 2D

ṗi = m(ψ i ) ψi = u ψi , (5)
where m = cos(ψ i ) sin(ψ i ) T with ψ i ∈ (-π, π] being the attitude yaw angle 1 and u ψi is the control action that will make the aircraft to turn. In particular, for coordinated turns where the altitude of the vehicle is kept constant and the pitch angle is close to zero, the control action u ψi corresponds to the following bank angle φ i to be tracked by the autopilot of the vehicle

φ i = arctan u ψi g , ( 6
)
where g is the gravity acceleration.

C. Trajectory tracking

One of the key points of the proposed formation control algorithm in this paper is to make sure that the aircraft is tracking C r . There exist many guidance algorithms in the literature [START_REF] Kendoul | Survey of advances in guidance, navigation, and control of unmanned rotorcraft systems[END_REF], [START_REF] Sujit | Unmanned aerial vehicle path following: A survey and analysis of algorithms for fixed-wing unmanned aerial vehicless[END_REF]. We have chosen the algorithm proposed in [START_REF] Kapitanyuk | A guiding vector field algorithm for path following control of nonholonomic mobile robots[END_REF] that has been successfully tested in real flights [START_REF] De Marina | Guidance algorithm for smooth trajectory tracking of a fixed wing uav flying in wind flows[END_REF] for two reasons. Firstly, the local exponential converge to the desired path is guaranteed. This property will help us later to support the convergence of the formation control algorithm under the argument of slow-fast dynamical systems in cascade. Secondly, the algorithm can be straightforwardly extended to other C 2 curves that are homeomorphic to C r , such as ellipses or the (possibly concave curve) Cassini ovals.

The trajectory tracking algorithm employs the level sets e(p) ∆ = ϕ(p) for the notion of error distance between the aircraft and C r . Note that for circular paths, the error e has a clear relation with the Euclidean distance, but for more general trajectories, such as ellipses, this is not always true. The vehicle has to follow the vector field defined by

ṗd (p) ∆ = τ (p) -k e e(p)n(p), (7)
where k e ∈ R + is a gain that defines how aggressive the vector field is, in order to converge to traveling on C r . Let us define x as the unit vector constructed from the nonzero vector x. Theorem 2.1: [START_REF] Kapitanyuk | A guiding vector field algorithm for path following control of nonholonomic mobile robots[END_REF], [START_REF] De Marina | Guidance algorithm for smooth trajectory tracking of a fixed wing uav flying in wind flows[END_REF] Consider the system (5), then the control action

u ψ = -E pd pT d E (E -k e e)H(ϕ) ṗ -k e n T ṗn T E ṗd || ṗd || 2 + k d pT E pd , (8)
where H(•) is the Hessian operator and k d ∈ R + , makes the aircraft to converge (locally) exponentially fast to travel over C r , i.e., for |e(0)| ≤ c * we have that |e(t)| ≤ a exp(-bt) with t → ∞ for some constants a, b, c * ∈ R + . The first term in [START_REF] Stipanović | Decentralized overlapping control of a formation of unmanned aerial vehicles[END_REF] makes the aircraft to stay on the guidance vector field [START_REF] Wang | Cooperative UAV formation flying with obstacle/collision avoidance[END_REF] while the second term makes the vehicle to converge to the guidance vector field in case that the vehicle is not aligned with it.

III. CIRCULAR FORMATION CONTROL

A. Problem definition

Given the neighbors' relationship described by the graph G, the stacked vector of inter-vehicle angles can be calculated as

z = B T θ, (9)
where θ ∈ T |V| = S 1 × • • • × S 1 (the n-torus) is the stacked vector of parameters for each vehicle as in (4), and z ∈ T |E| , therefore no necessarily all the inter-vehicle angles are calculated. Note that z k ∈ S 1 can be calculated from the relative measurement p ip j by just following trigonometric arguments in Figure 1. Consider a collection of desired z * k ∈ S 1 , k ∈ {1, . . . , |E|} on the circle. We define the formation error e θ ∈ T |E| as the stacked vector of signals

e θ k (t) = z k (t) -z * k , (10)
ϕ(p * ) = e > 0 ϕ(p) < 0 C r := ϕ(p) = 0 θ(p * ) O N p * τ n -k e en pd ṗψ
Fig. 1: The direction to be followed by the UAV at the point p * for converging to C r is given by pd . The tangent and normal vectors τ and n are calculated from ∇ϕ(p * ). The error distance e is calculated as ϕ(p * ). All the circles ϕ c can be parametrized by an angle θ with respect to the horizontal axis of a frame of coordinates at the center of C r .

where e θ k ∈ S 1 . The objective of the proposed algorithm for the team of fixed-wing UAVs in the next subsection is to achieve simultaneously e θ (t) → 0 and p i (t) → C r , ∀i ∈ N as t → ∞.

B. Controller design and stability

Consider that the unit speed aircraft i is tracking correctly C r , therefore its angular velocity around the center of C r is θi

= 1 r . ( 11
)
The idea is to control the inter-vehicle angles in z by changing in the vehicles the desired trajectory to be tracked, i.e., instead of (2), the vehicle i has to track

i C (r, i c) ∆ = {p : ϕ(p) = i c}, (12)
where i c ∈ R is the formation control signal to be designed and the superindex i denotes for the vehicle i ∈ V. Note that the smaller the i c (possibly negative), the bigger the radius of i C (r, i c) , and thus the smaller the angular velocity θi . For the sake of simplicity in the following analysis we define

i c ∆ = i u 2 r + 2r i u r , (13)
where i u r ∈ R is a control action with a more straightforward physical meaning than i c, i.e., we set the radius (around the desired r) of the circumference i C (r, i c) (or simply i C) because

x 2 + y 2 -r 2 = i u 2 r + 2r i u r ⇐⇒ x 2 + y 2 -(r + i u r ) 2 = 0. ( 14
)
Remark 3.1: Note that for a generic closed trajectory, if a vehicle tracks it on a negative level set, then it travels less distance than one that tracks the same trajectory on a positive level set after one loop.

We propose the following control action for achieving the desired circular formation

i u r = k r B i e, (15)
where B i stands for the i'th row of the incidence matrix B as in [START_REF] Beard | Decentralized cooperative aerial surveillance using fixed-wing miniature UAVs[END_REF], and k r ∈ R + . Since e θ k ∈ (-π, π], we impose to k r the following condition

r -πk r max i∈V (|N i |) > 0, (16)
i.e., we avoid the possibility of setting a negative radius2 in i C. Note that the control action ( 15) is based on the popular consensus algorithm in formation control [START_REF] Oh | A survey of multi-agent formation control[END_REF]. Before presenting the main result, we need the following technical lemma that will define the neighbors' relationships.

Lemma 3.2: If G does not contain any cycles, then the matrix A ∆ = -B T B is Hurwitz. Proof: If G does not contain any cycles, it has been show in [START_REF] Guattery | Graph embeddings and laplacian eigenvalues[END_REF] that Bx = 0 for any nonzero vector x ∈ R |E| . Note that ||Bx|| 2 = x T B T Bx > 0, implying that B T B is positive definite. Hence, A is Hurwitz if G does not contain any cycles.

We also make use of the following assumption. Assumption 3.3: A vehicle i is always tracking and traveling over i C as in [START_REF] Lee | Distributed formation control via global orientation estimation[END_REF].

The Assumption 3.3 considers that if there is a change in the radius of i C, then the vehicle instantaneously jumps to the required level set. As we will show, the circular formation controller [START_REF] Marshall | Pursuit formations of unicycles[END_REF] guarantees the exponential stability of the origin of the error signal e θ under the Assumption 3.3. Since the trajectory error tracking e in Theorem 2.1 is also locally exponentially stable, one may consider a slow-fast dynamics in cascade by tuning appropriately the gains k e and k r in ( 7) and ( 15) respectively [START_REF] Sepulchre | Constructive nonlinear control[END_REF]Chapter 4]. Informally, the controller (8) provides a fast transient process of the vehicle to i C if the aircraft is sufficiently close to it, while the whole formation slowly follows the formation controller [START_REF] Marshall | Pursuit formations of unicycles[END_REF]. In fact, if k r is sufficiently small, the set of possible trajectories i C will be very close to C r , therefore making reasonable the Assumption 3.3. Furthermore, since for circles the convergence of the algorithm in Theorem 2.1 is almost globally stable (with the exception of starting at the center of C r ), even if the vehicles start far away from the trajectories to be tracked, they will approach a situation where Assumption 3.3 can eventually be considered.

Theorem 3.4: Consider a team of n unit speeds aircraft modeled as in [START_REF] Ryan | An overview of emerging results in cooperative UAV control[END_REF], and the graph G defining the neighbors' relationships does not contain any cycles. All the vehicles are tracking [START_REF] Lee | Distributed formation control via global orientation estimation[END_REF] by employing [START_REF] Stipanović | Decentralized overlapping control of a formation of unmanned aerial vehicles[END_REF]. If Assumption 3.3 holds and the level sets i c in [START_REF] Lee | Distributed formation control via global orientation estimation[END_REF] are controlled by ( 15) via [START_REF] Lin | Necessary and sufficient graphical conditions for formation control of unicycles[END_REF], then the origin of the error e θ as in [START_REF] Anderson | UAV formation control: Theory and application[END_REF] is locally exponentially stable for the desired z * k , k ∈ {1, . . . , |E|}.

Proof: The proof is based on checking the stability of the linearization of the error dynamics e θ around the origin. First note that ėθ = ż = B T θ. According to Assumption 3.3, we also have that for each edge E k = (i, j) the agents i and j are tracking a circle of radius r + k r B i e θ and r + k r B j e θ respectively, so from [START_REF] Montenbruck | Fekete points, formation control, and the balancing problem[END_REF] it holds that

ėθ k = żk = 1 r + k r B i e θ - 1 r + k r B j e θ = k r (B j -B i )e θ (r + k r B i e θ )(r + k r B j e θ ) = k r A k e θ (r + k r B i e θ )(r + k r B j e θ ) , k ∈ {1, . . . , |E|}, (17)
where A k is the k'th row of the matrix A as in Lemma 3.2. We linearize ( 17) around e θ = 0, therefore the dynamics of small variations θ k of the error are given by

˙ θ k = ∂ ėθ k ∂e θ e θ =0 θ = k r r 2 A k θ , (18)
which leads to the following compact form

˙ θ = k r r 2 A θ , (19)
and because A is Hurwitz according to Lemma 3.2 since G has not any cycles, we can conclude that the equilibrium e θ = 0 is locally exponentially stable. Remark 3.5: Note that since the convergence to the trajectories i C is asymptotic, one can guarantee that all the vehicles will be confined in a disc D of radius (r + πk r max i∈V

(|N i |),
which corresponds to the worst case radius to be tracked, for all time t, even if the control u r is not updated, e.g., the vehicles are not exchanging or sensing their positions.

It is interesting to highlight that if G does not contain any cycles, then in such a disc D the only equilibrium point for the system is at e θ = 0, which has been proven stable. Since the vehicles are eventually confined in D according to Theorem 2.1, it seems reasonable to conjecture that an estimation of the region of attraction for the exponentially stable e θ = 0 is indeed D. Furthermore, for a proof of the convergence of the overall system without Assumption 3.3, one can use the stability theory of cascade systems [START_REF] Sepulchre | Constructive nonlinear control[END_REF]Chapter 4], while the exponential stability of the partial system (17) could guarantee the (locally) asymptotic stability of the overall system. We will present a rigorous proof in the extended journal version.

IV. IMPLEMENTATION AND FLIGHT PERFORMANCE

A. Experimental platform

The validity of Theorem 3.4 has been tested with the three fixed-wing UAVs shown in Figure 2. The aircraft have about 600 grams of weight, 1.2 m of wingspan, and they are actuated by two elevons and one motor. The electronics include a battery that allows about 45 minutes of autonomy at the nominal flight, which corresponds to an airspeed of 13 m/s. The chosen board for running the Paparazzi autopilot stack is the Apogee [START_REF] Paparazzi | UAV open-source project[END_REF], which includes the usual sensors of

B. Circular formation flight experiment

The formation flight 3 has been taken place at the aero model club of Eole at Muret, close to the city of Toulose in France. The flights were performed on the 17th of February, 2017 between the 10:00 and the 12:00 hours local time. The wind coming from the south was about 2 m/s according to MeteoFrance, therefore we can consider that ground-speed and airspeed are approximately equal. The airplanes 1, 2 and 3 are tagged with the colors blue, pink and red respectively at the ground station captions in Figure 3. The chosen incidence matrix for the communication between vehicles is

B =   1 0 -1 1 0 -1   , (20)
which clearly does not define any cycles. The desired formation is defined by z * 1 = z * 2 = 0, i.e., all the aircraft achieve consensus for their corresponding θ i . Potential collisions are avoided by making the airplanes to fly at different altitudes, which are 30, 32 and 35 meters above the ground for planes 1, 2 and 3 respectively. The target circle C r is set at the waypoint CIRCLE in Figure 3 with a desired radius r equal to 80 meters. The gains k d , k e and k r in Theorems 2.1 and 3.4 have been set to 1, 1 and 8 respectively. The airplanes exchange their positions with a frequency of 2 Hz, although lost communications are expected. Note that each airplane has a different understanding of e θ {1,2} , i.e., each airplane calculates on board the error signal and it might be different among neighbors due to lost communications or GPS delays. Interesting work studying the effects of this fact in formation control can be found at [START_REF] Mou | Undirected rigid formations are problematic[END_REF], [START_REF] De Marina | Design and implementation of formation control algorithms for fully distributed multi-robot systems[END_REF].

Before the formation control algorithm begins the three aircraft are orbiting at different standby points. The experiment starts at time 131 seconds in Figure 3. At this moment 3 The video of the experiment with HD quality can be watched at https://www.youtube.com/c/HectorGarciadeMarina . the algorithm commands the airplane 2 (red) to follow a circumference with a much smaller radius than airplanes 1 (blue) and 3 (pink) in order to catch them up. In fact, airplanes 1 and 3 are tracking a circumference with a bigger radius than C r in order to wait for airplane 2. In about 15 seconds, the errors e θ {1,2} have been reduced half. Some lost communications have been experienced between times 150 and 170 seconds. However, the algorithm seems robust against such an issue and the formation achieves consensus within a band of ±10 degrees of error, and continues stable until the end of the experiment after seven laps.

V. CONCLUSIONS

This paper has presented an algorithm for achieving circular formations with fixed-wing UAVs traveling with constant speeds. The strategy consists in controlling the angular velocities around the center of the desired circle. For that, we design a control action that is applied to the level set to be tracked around the desired circle. These level sets are tracked by a guidance vector field algorithm that is locally exponentially stable. Since the presented algorithm for circular formations is also exponentially stable if the vehicles are perfectly tracking the level sets, we employ the argument of slow-fast systems in cascade in order to show the compatibility of both algorithms, which has been demonstrated in practice with three aircraft. The algorithm can be potentially extended to other closed non-circle trajectories. A more rigorous analysis will be presented in an extension of this work. degrees is between the dashed lines. Each dot represents the calculated error by one of the airplanes once it receives the position of its neighbor. The density of points is different due to lost communications. However, the algorithm seems robust enough against such an issue as it drives both errors close to zero within a band ±10 degrees.

Fig. 2 :

 2 Fig. 2: The three fixed-wing UAVs employed for the circular formation at the aero club of Eole at Muret (Toulouse).

  (a) t = 131secs (b) t = 145secs (c) t = 162secs (d) t = 410secs

Fig. 3 :eθ 2 Fig. 4 :

 324 Fig. 3: Screenshots from the Paparazzi ground control station showing the evolution of the circular formation.

For our setup, the yaw angle and heading angle can be considered equal due to the absence of wind.

Note that while a level set can be negative, the radius of a circumference cannot.


Vincent Vidal
email: vincent.vidal@liris.cnrs.fr

Christian Wolf
email: christian.wolf@liris.cnrs.fr

Florent Dupont
email: florent.dupont@liris.cnrs.fr

Combinatorial mesh optimization

Keywords: Triangular meshes, mesh optimization, discrete optimization, graph cuts

. The desired mesh properties are controlled through a global energy function including data attached terms measuring the fidelity to the original mesh, shape potentials favoring high quality triangles and connectivity as well as budget terms controlling the sampling density. The optimization algorithm modifies mesh connectivity as well as the vertex positions. Solutions for the vertex repositioning step are obtained by a discrete graph cut algorithm examining global combinations of local candidates.

Results on various 3D meshes compare favorably to recent state-of-the-art algorithms. Applications consist in optimizing triangular meshes and in simplifying meshes, while maintaining high mesh quality. Targeted areas are the improvement of the accuracy of numerical simulations, the convergence of numerical schemes, improvements of mesh rendering (normal field smoothness) or improvements of the geometric prediction in mesh compression techniques.

Introduction

Nowadays, 3D triangular meshes are commonly used in many fields, as the hundreds of thousands of existing 3D triangular models can attest. They are of interest in visual effects, video games, scientific visualization, 3D animation and medical surgery simulation based on finite element methods, to name a few. Most of the existing triangular meshes are of unsatisfying quality because of their inappropriate vertex sampling, which is responsible for inequilateral triangles and irregular connectivity. The origin of this irregular sampling may be due to the scanning device, to 3D interactive solid modeling software or to simplification algorithms.

This poor quality causes instability and divergence of various mesh processing applications [START_REF] Shewchuk | What is a good linear element? Interpolation, conditioning, and quality measures[END_REF]. In that context, several remeshing techniques have been introduced in the literature. They consist in improving some quality requirements under soft and/or hard constraint conditions. Targeted goals vary according to the application [START_REF] Alliez | Recent Advances in Remeshing of Surfaces[END_REF]. Simplification techniques tend to preserve the overall shape of the mesh while removing as many triangles as possible. Mesh smoothing methods consist in removing high frequency noise so as to fair the mesh. Finally, Mesh optimization aims at improving the mesh quality, i.e. the regularity of the sampling and of the connectivity.

In this paper, the mesh optimization problem is stated as an energy minimization problem. The proposed energy function (see equation 4) captures the fidelity to the initial mesh geometry, the mesh quality (expressed in terms of triangle shape and vertex valence) and the number of vertices (simplification versus remeshing). Hence, according to the initial mesh and desired goals of the method, it is possible to improve the compactness of the representation, the triangle shape, the vertex valence, while controlling the geometric error introduced by such competing desires.

We solve the proposed energy minimization problem with an iterative two step minimization method. The first step consists of a greedy mesh connectivity energy minimization process. Edge flips and edge collapses/splits are done only if they decrease the energy function. The second step performs global vertex repositioning, which is formulated as a discrete max-flow/min-cut problem in a directed graph. More specifically, this problem is solved using a variant of QPBO (Quadratic Pseudo-Boolean Optimization), a recent graph cut technique which guarantees a nonincrease of the energy function after the vertex relocation step. As a result, the energy minimization problem is solved with guarantees on its convergence. Figures 4 and5 show some results of our regularization algorithm.

Related work

Two classes of approaches have been identified in the vast mesh optimization literature. The first group consists of methods which do not offer a high control over the geometric error and/or over the mesh connectivity when new points are sampled on the surface. Variational partitioning frameworks [START_REF] Cohen-Steiner | Variational shape approximation[END_REF][START_REF] Valette | Approximated centroidal voronoi diagrams for uniform polygonal mesh coarsening[END_REF][START_REF] Valette | Generic remeshing of 3d triangular meshes with metric-dependent discrete voronoi diagrams[END_REF] are based on vertex/triangle clustering and are often used to produce a coarser mesh with a high approximation quality. This coarse mesh is then usually retiled or refined using local re-triangulation to improve the point sampling. Parameterization methods [START_REF] Alliez | Isotropic surface remeshing[END_REF] optimize 2D patches instead of a 3D surface, but they generally suffer either from distortion produced by the global parameterization, or need patch boundary post-processing resulting from local parameterization. Semi-regular remeshing [START_REF] Guskov | Hybrid meshes: multiresolution using regular and irregular refinement[END_REF] uses an initial coarse mesh partition and treats each patch separately using subdivision rules. This produces a few irregular vertices (those of the coarse mesh), well-shaped triangles, and a small geometric error. However, it is sensitive to the patch structure and the resulting vertex sampling is difficult to control. Geodesic front propagation techniques [START_REF] Sifri | Geodesic-based surface remeshing[END_REF][START_REF] Peyré | Geodesic remeshing using front propagation[END_REF] consider geodesic equidistant curves over the surface and allow to get well-shaped triangles (vertices can be distributed according to the local curvature). However, some post-processing steps are needed to avoid artifacts when the curve topology is complex and the geometric approximation error is located near sharp features. All these methods give an overall good triangle shape, but they suffer from a lack of geometric error control when new points are sampled on the surface, since they cannot locally control the geometric approximation of some high frequency features.

The second group of methods works directly on the initial mesh simplices (vertices, edges and resulting triangles), which allows a better control of the geometric fidelity to the original 3D surface. Local approaches to mesh optimization consist in using a set of local legal moves (e.g. towards barycenters or angle-bisectors) and connectivity modifications (topological operators) to decrease an energy function. Such an approach may lead to a local minima configuration, especially in combination with greedy optimization (e.g. gradient descent). Surazhsky and Gotsman [START_REF] Surazhsky | Explicit surface remeshing[END_REF] use local operations such as edge-collapse, edge-split and edge-flip to regularize the mesh connectivity. Surazhsky et al. [START_REF] Surazhsky | Explicit surface remeshing[END_REF][START_REF] Surazhsky | Isotropic remeshing of surfaces: a local parameterization approach[END_REF] apply area-based smoothing to control both triangle quality and vertex sampling over the mesh. To achieve a precise isotropic vertex placement, Surazhsky et al. [START_REF] Surazhsky | Isotropic remeshing of surfaces: a local parameterization approach[END_REF] use a centroidal Voronoi tessellation (CVT). Global approaches to mesh optimization attempt to resolve the vertex repositioning problem in a global way, most of the time solving a sparse linear system [START_REF] Winkler | Mesh massage: A versatile mesh optimization framework[END_REF][START_REF] Liu | Non-iterative approach for global mesh optimization[END_REF] or using least squares approximation [START_REF] Nealen | Laplacian mesh optimization[END_REF]. The main idea in the global approaches using a so-called Laplacian global operator [START_REF] Winkler | Mesh massage: A versatile mesh optimization framework[END_REF][START_REF] Liu | Non-iterative approach for global mesh optimization[END_REF][START_REF] Nealen | Laplacian mesh optimization[END_REF], is to infinitely apply a Laplacian operator such that applying it one more time will not change the current vertex positions. That allows a direct formulation as a linear system. Then other constraints are added to take into account invariant vertex positions or to avoid the shrinking effect due to Laplacian smoothing. These global vertex repositioning techniques depend on the initial sampling and connectivity, and the triangle shape may therefore be difficult to improve in case of irregular configurations (e.g. irregular vertex degrees and/or low sampling). In most local approaches, connectivity and the vertex positions are optimized separately as in [START_REF] Hoppe | Mesh optimization[END_REF]. That is essentially due to the fact that the combinatorial complexity of mesh connectivity optimization does not allow to solve it neither globally nor jointly with the vertex repositioning problem.

Other works related to the surface sampling improvement (e.g. direct sampling) or to the approximation error control (e.g. volume computation between two meshes) may be of interest to the reader, but are beyond the scope of this paper.

The proposed method is characterized by the following advantages:

-Many objects present a noticeable amount of sharp features (sharp edges and corners), which are robustly detected and preserved during our optimization process (see section 2.3). The coherency of geometric features is exploited in the feature edge detection process to be more robust to the presence of noise.

-Geometric error is kept low, while triangle quality and compactness of the representation are improved. No high frequency noise is introduced in the smooth parts of the processed objects during all optimization stages. -Vertex regularity is improved only when it does not penalize the triangle quality and the geometric error.

-The locality of our approach allows a high control on the proposed vertex new positions, which can be exploited to favor vertex positions alignment along a curve.

Our main contributions are two-fold:

-A discrete-continuous mesh position optimization algorithm for triangular meshes. New vertex positions are proposed locally and the decisions which keep either the new vertex position or the old one are globally taken based on the calculation of the minimum st-cut/maximum flow in a graph. This ensures a local control over the candidate positions while minimizing the energy globally and avoiding oscillation problems. -A method for feature edge detection based on an improved Potts model, which favors contiguous lines of sharp edges with approximately the same dihedral angles and directions for adjacent sharp edges.

The paper is organized as follows: section 2 introduces the notations used throughout the paper, presents an outline of our regularization algorithm and explains how the method copes with sharp features. In section 3, we present our global energy minimization problem and in particular we define the objective function. Section 4 deals with the minimization of the objective function and more precisely the optimization of mesh vertices and connectivity. In section 5, we give experimental results and compare them with other recent algorithms from the state of the art. In section 6, we conclude and discuss some future work.

2 Outline of the proposed mesh optimization algorithm

Notations

The following notations are used throughout the paper:

-M (resp. M') denotes the initial mesh (resp. optimized/remeshed model -T represents the remeshed model triangles. (s, r, q) ∈ T iff s ∈ N r and s ∈ N q and r ∈ N q .

Algorithm description

The global scheme of our method for optimizing 2manifold triangular meshes is illustrated in figure 1 and in algorithm 1. In the initialization phase, we compute principle curvatures and detect geometric features, which remain constant during the subsequent steps. Curvature extraction is based on normal cycles calculus [START_REF] Cohen-Steiner | Restricted delaunay triangulations and normal cycle[END_REF]. The feature detection procedure is explained in section 2. Fig. 1 The remeshing pipeline: at every iteration i, new vertex positions are chosen using a global energy minimum approximation; less frequently, the mesh connectivity is improved to favor better global vertex configurations.

(e.g. [START_REF] Jiao | Feature detection for surface meshes[END_REF]). Geometric features are generally defined as discontinuities in the normal direction (i.e. first derivatives) or boundaries of the object.

In our method, feature detection allows to keep corner positions unchanged and to restrict vertices on feature lines to stick to these lines in the vertex optimization step. In addition, during the connectivity optimization stage, feature edges are prohibited from modifying their direction. For instance, feature edge flipping is forbidden. By doing so, the intrinsic geometric properties of the initial triangular mesh are preserved.

Intuitively, discontinuity edges or feature edges could be detected by thresholding the edge dihedral angles θ i to classify them as either normal or sharp using a threshold θ T h . However, due to the presence of noise, feature edges may be difficult to characterize with this technique. For instance, some geometric models do not have any dihedral angle threshold capable of distinguishing feature edges from others without connecting true features with geometric artifacts or vice versa.

In the same lines of our vertex repositioning optimization procedure described in the rest of the paper, the classification performance of the feature edge detection process can be significantly boosted by taking the decisions for all edges of the 3D model jointly, i.e. globally. This can be achieved through a model which complements the information residing in the dihedral angles of all edges with additional terms favoring consistent feature lines, i.e. coherent neighboring edge labels. We propose an improved Ising/Potts model, known from image restoration [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], which in our case minimizes a global energy function over all edge labels:

Ŵ = arg min W E(W ; θ, φ) = arg min W i E d (W i ; θ i ) -µ i,j:i∈N j S(φ ij )δ(W i , W j ) ( 1
)
where W is the set of all indexed edge variables W i , θ is the set of all indexed dihedral angles θ i , and φ is the set of all indexed turning angles φ ij measured on pairs of edges. Each binary edge variable W i can take a value/label in L={0, 1}, which correspond to normal edge and feature edge. N i is a set of edge indices which correspond to edge variables adjacent to the edge indexed by i. The energy function is characterized by data attached terms E d which push the solution into the direction determined by the dihedral angles, as well as regularizing pairwise terms which favor geometrically coherent feature lines. The weight µ sets the relative strength of the pairwise terms compared to the data attached terms.

E d (W i , θ i ) is a data attached term taking into account the dihedral angle θ i of the edge and which has been designed as a softened thresholding with threshold θ T h and linear dependency of the energy on the angle:

E d (W i ; θ i ) = cos(θ i ) if W i = 1 ⇔ feature 2cos(θ T h )-cos(θ i ) else W i = 0 ⇔ normal (2)
δ is the Kronecker symbol equal to one whenever the condition W i = W j is satisfied and zero otherwise. δ can be interpreted as a term favoring the same labeling for neighboring edges. S(φ ij ) is a function which controls the strength of the regularizing term through a measure of geometric coherency, i.e. homogeneous labels are favored more if the neighboring edges are similarly oriented:

S(φ ij ) = exp {-50.(1 -cos(φ ij ))} . (3)
All pairwise terms in (1) are submodular (see also section 4.1.2 on submodularity), therefore the global minimum can be very efficiently computed using graph cuts and Kolmogorov et al.'s st-graph construction [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF].

Once the feature edges are detected, vertices can be labeled given the decisions on the edges: a vertex is considered as a corner if it has at least three adjacent feature edges, or if it has two adjacent and non-aligned boundary edges.

The global energy minimization problem

In this work the initial mesh M is copied and kept as reference geometry during the remeshing process. The optimized mesh M' starts with the same vertices and connectivity as M. All vertex positions Y of M remain constant during the optimization process, but the set X may evolve (addition/removal of vertices) since the initial sampling may not be the one which minimizes the objective function (4) given below, which depends both on the vertex positions and on the mesh connectivity. To each vertex of the remeshed model we keep a link to the closest position of the original mesh, which allows us to keep track of the geometric distance between a vertex of the optimized mesh and the surface represented by the original mesh (cf. figure 2).

Definition of the objective function

The specific form of our objective energy function U (X; Y ) is defined as a scalar combination of sums of positive and unitless energy potentials. Each potential, which we also call feature function, locally measures a criteria (geometric error or a quality) and decreases (resp. increases) when the criteria is locally more (resp.

Algorithm 1: The whole method, including the continuous-discrete solution to vertex relocation. K (0) , C and i max are, respectively, start temperature, cooling speed and number of iterations. Remaining notations are defined in section 2.1. The graph cut algorithm takes the global decision minimizing the energy presented in section 3 for the whole set of vertices, considering for each vertex of the remeshed model the current position X s and a new candidate position X new s .

Input: M(Y ,F ), K (0) , C, imax Output: M'(X,E) Compute curvatures and feature edges X ← Y ; E ← F; K ← K (0) ; for i ← 0 to imax -1 do if i mod 5 = 0 then Regularize connectivity end for Xs ∈ X do σ f reedom ← 0.5 1 + e -K .||Xs, X Ns ||g.min(1, ρmax) X new s ←
first with ∆sU < 0 and in sphere, out of the following:

       angle based smooth. Laplacian smooth. guided uniform random end X ← globally optimize QPBO-P (X, X new , Y ) K ← K • C end less) respected.
Unitless energy potentials have the advantage of being scale invariant.

U (X; Y ) = λ s (s,r,q)∈T ψ s (X s , X r , X q ) Shape + λ d (s,r,q)∈T ψ d (X s , X r , X q , Y ) Fidelity + λ v s ψ v (X s ) Valence + λ p s 1 Penalty (4)
Here λ s , λ d , λ v and λ p are positive scalar weights.

The respective label subscripts s, d, v, and p denote shape quality, data fidelity, valence quality, and vertex penalty, respectively. The other symbols have been defined in section 2.1. U (X; Y ) evaluates the global configuration of the optimized mesh, and assigns a scalar energy value to each possible solution, i.e. each possible result mesh. To make the equation easier to read, we abusively simplified the dependency on the mesh connectivity in the notation; in particular, there are no variables associated to the connectivity information. It should nevertheless be clear that the energy function depends on the vertex positions as well as the mesh connectivity.

As a consequence, at each vertex relocation or at each mesh connectivity modification, U (X; Y ) may vary. By setting the parameters λ, a user can efficiently and easily create meshes with the desired properties.

The feature functions

Our model contains four different kinds of positive feature functions: a shape function measuring the quality of the mesh triangles, a data attached function encoding the approximation quality of the new surface, a valence potential function encoding the quality of the mesh connectivity, as well as a term penalizing a high vertex budget, which allows the user to control the number of vertices of the resulting mesh.

The role of the first type is to favor equilateral mesh triangles, therefore it is calculated on triangles: ψs(Xs, Xr, Xq) = R(Xs, Xr, Xq) min(||Xs -Xr||, ||Xs -Xq||, ||Xr -Xq||) [START_REF] Valette | Approximated centroidal voronoi diagrams for uniform polygonal mesh coarsening[END_REF] where R(X s , X r , X q ) denotes the circumradius associated with the triangle (X s , X r , X q ) and ||.|| is the usual Euclidean norm. Note that this feature function does not depend on the initial vertices Y and that it is scale invariant. ψ s (X s , X r , X q ) can be extended to take the value +∞ when the triangle (X s ,X r ,X q ) is degenerated (singularity of the function). The minimum of this feature function for one triangle is reached when the triangle is equilateral.

The data attached feature function ψ d (X s , X r , X q , Y ) measures the geometric error between a triangle (X s , X r , X q ) of the optimized model and the initial mesh. Ideally, ψ d should be equal to the absolute volume error produced by a local operation (either vertex repositioning or topological operation on edges). However, exact volume error computations may add extra cost to the whole iterative energy minimization process, especially for the vertex repositioning step, where we calculate two different candidate positions for each vertex (see section 4), which gives eight possibilities to check for each triangle. We therefore calculate the exact direct volume error only during the optimization of mesh connectivity. In particular, in the latter case the tetrahedron volume error is calculated over the subsequent local edge operations (and not between the optimized and initial meshes for computation time purposes). More precisely, the volume error associated with an edge flip, an edge collapse or an edge split are, respectively, the tetrahedron volume induced by the 4 points of the 2 adjacent triangles, the sum of the tetrahedron volumes induced by each one-ring neighboring edge and its adjacent old and new center vertex position, and zero.

During the optimization of the vertex positions, we approximate the fidelity term ψ d by a point-to-surface distance, whose main advantage is that it can be calculated vertexwise:

ψ d (X s , X r , X q , Y ) = F (X s , Y ) + F (X r , Y ) + F (X q , Y ) (6)
where F (X s , Y ) is the square of the shortest distance between X s and the initial mesh, i.e. the geometric distance. The distance F could be approximated up to second order with the Pottmann distance [START_REF] Pottmann | Geometry of the squared distance function to curves and surfaces[END_REF][START_REF] Pottmann | Fitting B-spline curves to point clouds by curvature-based squared distance minimization[END_REF], which allows to directly compute the gradient of the energy. Since this is not necessary in our discrete framework, a direct computation of F as the orthogonal projection distance on the initial mesh was possible, which is slightly faster and more robust in the presence of sharp features.

The proposed ψ d (X s , X r , X q , Y ) approximation works very well provided that the weighting factor λ d is high, which is the case for the targeted applications, namely regularization where geometric fidelity must be high.

The valence potential function is calculated on optimized mesh vertices:

ψ v (X s ) = (|N opt | -|N s |) 2 (7)
where |N s | is the number of edges adjacent to vertex s, i.e. the vertex degree, and |N opt | is the optimal (desired) vertex degree: 4 on borders and 6 otherwise. The vertex budget energy cost is linear in the number of remeshed model vertices, and as can be seen from equation (4), each vertex costs λ p . In general, λ p controls the desired amount of vertices.

Minimization of the objective function

The next two subsections describe how the unique global energy function is minimized through two dif-ferent steps, which update the vertex positions and the mesh connectivity.

Global vertex repositioning

Updating the vertex positions of the current mesh requires minimizing equation ( 4) over all variables X s ∈R 3 , a continuous optimization problem. The energy potentials associated with the vertex valence and with the vertex budget are constant at this stage and therefore omitted.

Unfortunately, the function U (X; Y ) is not convex and standard gradient descent methods will most likely return a sub-optimal solution. In addition, a leastsquare or linear system approach will return an oversmooth mesh, even if some constraints are added. This suggests a discrete approach for obtaining high fidelity meshes. Our work benefits from recent advances in optimization theory for discrete Markov Random Fields (MRFs) [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF] by transforming the continuous problem into a discrete problem, similar to the technique proposed for optical flow by Lempitsky et al. [START_REF] Lempitsky | Fusionflow: Discretecontinuous optimization for optical flow estimation[END_REF]. However, instead of employing a global discrete optimizer to merge several solutions obtained by existing techniques applied with different parameters, in our case the global discrete optimizer takes decisions on candidates calculated at each step in an iterative process. At each iteration and for each vertex of the remeshed model, a new candidate position is proposed and the optimal decision for the whole set of vertices is calculated, i.e. the decision minimizing (4). In the following two subsections 4.1.1 and 4.1.2, we explain how we generate new candidates and how we globally decide to keep or not a new position.

Local candidate proposals

We force each valid candidate position to stay within a small freedom sphere to prevent too much geometric error and to avoid creating a geometric fold-over. A new vertex candidate position will be rejected if it is outside of the corresponding freedom sphere, leaving the current vertex position unchanged. The more the region around a vertex X s is curved, the smaller its freedom radius σ f reedom is forced to be, which avoids large moves around a point of high curvature and thus limits the introduced geometric error. Moreover, like the stepwidth in gradient descent, this radius decreases with time to avoid big moves at the end of the optimization process. Similarly to simulated annealing techniques [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], a temperature parameter K decreases at each iteration (introduced in the algorithm 1 and in section 5.1). The radius σ f reedom is related to this temperature K through a sigmoidal function, as well as to the local geodesic radius of the X s one-ring, and to the local maximal curvature radius measured on the initial mesh:

σ f reedom (Xs, X Ns , Y ) = 0.5 1 + e -K .||Xs, X Ns ||g.min(1, ρmax) (8)
where ρ max is the maximum absolute curvature radius at the closest initial mesh vertex Y s from X s (σ f reedom thus depends on Y ). ||X s , X Ns || g is the maximum radius such that the sphere centered on X s with that radius does not intersect the X s one-ring neighborhood of edges. In other words, it is the local geodesic radius, computed as the minimum Euclidean distance from X s to its one-ring neighborhood of edges:

||X s , X Ns || g = sup ρ∈R +    ρ : ρ < ||X s , u|| ∀u on (r, q) ∈ E, ∀r ∈ N s , ∀q ∈ N s    (9)
It is easy to see that a freedom radius strictly less than 0.5||X s , X Ns || g for each optimized mesh vertex X s prevents the creation of new fold-overs (see figure 3). For the same reason, the factor min(1, ρ max ) in equation ( 8) limits the radius to the value described above. Since the initial mesh is normalized (its coordinates are divided by its bounding box diagonal) before processing and unnormalized at the end, the diagonal of the mesh bounding box does not appear in the equation involving the curvatures radius.

A good candidate position X new s for replacing X s must decrease the energy. The local energy variation ∆ s U due to the new vertex position X new s can be computed quickly from the vertex X s and its one-ring neighborhood N s . Because of the form of the energy function (4),

∆ s U = U X \ {X s } ∪ {X new s }; Y -U (X; Y ) (10)
can also be computed much quicker from only a small number of terms as

∆ s U = U {X new s , X Ns }; Y -U {X s , X Ns }; Y (11)
where the so-called local evidence U {X s , X Ns }; Y = U (X; Y ) |{Xs,Ns} only contains the terms in U (X; Y ) which involve the modified mesh vertex X s or the one ring X Ns .

The global algorithm which jointly takes decisions on candidate pairs (i.e. keeping either X s or X new s ) for the whole mesh, is given below in section (4.1.2). The convergence of the algorithm depends on the quality of the new candidate positions, which leads us to choose one candidate for each vertex among several candidates calculated by different methods. They have been ranked empirically in experiments, and the chosen candidate is the first one satisfying the freedom sphere constraints given above, i.e. X new s must strictly be in the freedom sphere and ∆ s U < 0, cf. ( 8), ( 10) and [START_REF] Surazhsky | Explicit surface remeshing[END_REF]. If no such candidate is found, the old position is kept. The techniques are, in decreasing rate of empirically measured convergence, angle-based smoothing, Laplacian smoothing, uniformly taken and guided random candidates. They consist in computing a displacement vector -→ v with different rules according to the method, from which the new candidate position X new s is set:

X new s = X s + - → v .
For angle-based (resp. Laplacian) smoothing, -

→ v = γ - → v (γ ∈ R + )
, where -→ v is computed as the mean direction towards surrounding angle bisectors (resp. the umbrella-operator U (X s ) = 1 r∈Ns zr r∈Ns z r X r -X s where the weights z r are set to the local area dispersion). For uniformly taken and random candidate positions, -→ v has the global energy gradient direction for gradient-guided iterations and is a direction on the local tangent plane otherwise. Special care is taken for vertex positions on feature edges, for which new candidates are forced to lie on the feature edge and which are kept only if the surrounding feature edge directions do not change with this new candidate. Vertices on corners remain unchanged.

Global candidate decisions

The global decision (i.e. keeping the current position or choosing the new candidate) on the whole set of vertices is taken by a graph cut technique. This involves constructing an st-graph representing the energy function (4) such that the minimum cut/maximum flow on this graph will give the solution which globally minimizes the energy. Although there are known graph cuts techniques able to optimize functions of some restricted classes where each variable may take values from sets larger than 2, the eventual gain in quality is far outweighed by the high computational complexity. We therefore concentrated on a single new candidate proposal for each vertex, which leads to binary valued functions, i.e. each variable associated to vertex may take values from a set of two values. For that purpose, it is necessary to associate a binary labeling to {X s , X new s }, for instance 0 could denote the current vertex position X s and 1 the proposed position X new s . Unfortunately, the energy function ( 4) is not submodular, which makes its exact global minimization with graph cuts difficult [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. Submodularity is an important concept in discrete optimization theory. Its meaning is the equivalent of the "convex" concept, restricted to discrete sets. This criterion requires that, taking one of the vertices of a given triangle (X r , X s , X t ) and fixing the decision to either 0 or 1, the projection onto the two other vertices satisfies the following constraint (without loss of generality we suppose that X r has been fixed):

ψ s (X r , 0, 0) + ψ s (X r , 1, 1) ≤ ψ s (X r , 0, 1) + ψ s (X r , 1, 0) (12)
In the case of ( 4), this criterion is generally not satisfied for the shape term ψ s (X s , X r , X q ) of every triangle, although a subset of triangles may satisfy the criterion.

In the proposed method, ( 4) is approximatively minimized with a variant of QPBO (Quadratic Pseudo-Boolean Optimization), a graph cuts algorithm which can deal with non-submodular terms (see section 5.1). However, the output of QPBO is a partial labeling, i.e. for some binary variables, the algorithm does not know if it is better to choose zero or one. On the other hand, convergence is guaranteed since, starting from an initial configuration, QPBO guarantees that the new labeling does not increase the energy (robustness). For a detailed explanation of QPBO and its properties, see [START_REF] Kolmogorov | Minimizing nonsubmodular functions with graph cuts-a review[END_REF].

Iterative mesh connectivity optimization

Our goal is to optimize the remeshed model connectivity by minimizing the objective function (4). Unfortunately a global optimization of the connectivity is intractable -the optimal scheduling of local topological operations cannot be computed in a reasonable amount of time. Here, a greedy scheme has been adapted to allow control through the global energy function (4), i.e. a local topological operator will be applied only if it decreases the energy.

Three priority queues handle local topological operations, respectively for edge-flips, edge-splits and edgecollapses. Each priority queue only contains valid local operations, which do neither create a vertex fold-over, nor change the direction of a feature edge and which will decrease the energy if applied. The priority of each queue is directly related to the energy decrease.

Again, calculating the optimal order of operation types -and therefore of priority queues -is intractable. Experimental tests (and intuition) led to the choice of giving higher priorities to flips, then to splits, and last to collapses.

Each priority queue is mutable: after the application of a local operation, each neighboring edge that already was in the priority queue is updated, i.e. removed if the local operation is no more applicable or its priority is updated. To make sure that the connectivity optimization step will terminate in a finite number of local operations, we limit the maximum number of consecutive operations of the same topological operator.

Topological changes and their energy -edge collapses and splits change the number of triangles of the mesh and therefore also the number of terms of the energy function. As a consequence, a change of the energy function is not necessarily an image of a quality/fidelity change of the mesh, but related to differences in the amount of mesh simplices. The budget term controlled through weight λ p has been designed to compensate for this. Its role is actually two-fold:

-To compensate for the lost terms in topological changes, leading to roughly unchanged vertex counts before and after optimization, if λ p has been judiciously chosen. -To control the vertex budget, if desired. A mesh simplification algorithm producing high quality meshes is obtained simply by increasing λ p .

Experimental results

In order to demonstrate the efficiency of our method, we applied it to several mesh models with unadapted sampling and very irregular connectivity, and which contain both, smooth parts and sharp features. Some visual results are given in figures 4, 5, 6, 7 and 8. The presented results have been obtained on an Intel Core 2 Duo P8400 (2.26 GHz) with 4 GB RAM with running time between 31 s and 3 min 56s. About 50% (resp. 50%) of the total running time is used for the connectivity optimization (resp. vertex repositioning). Clearly, the complexity of the connectivity optimization is linear in the number of edges of the input mesh, while the complexity of vertex repositioning is directly related to the complexity of the chosen graph-cut technique. The selected graph-cut technique is based on a push-relabel version of the max-flow algorithm with polynomial complexity in the worst case. Fortunately, the empirical (or average) complexity is nearly-linear with respect to the size of the input [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]. The average complexity of our method is nearly-linear with respect to the size of the input mesh (e.g. number of vertices, edges and triangles), which has been confirmed by our experiments (cf. figure 11).

Implementation details and settings

The temperature and cooling parameters mentioned in section 4.1.1 and detailed in algorithm 1 allow to decrease the radius of the freedom sphere associated with a given vertex when the number of iterations/the elapsed time increases. Their names are borrowed from physics, more precisely from the solidification of a melting metal. For the cooling schedule we used the suggestions in [START_REF] Duda | Pattern Classification, 2nd Edition[END_REF] (page 356), setting the temperature K to K (i) = K (0) • C i where i denotes the current iteration (0 ≤ i ≤ i max -1). The initial temperature K (0) has been set to 100 and the constant controlling the speed of the cooling process C has been set to 0.95. The chosen total number of iterations of our algorithm is i max = 170: a higher number of iterations does not improve significantly the results while adding extra computation time on the presented models (cf. figure 12). If the energy gradient can be computed for the vertex optimization step, the last 20 iterations are gradientguided which tends to fine-tune the final result.

Concerning the detection of robust features (cf. figure 13), we experimentally set θ T h = 35 o and µ = 0.1. We also experimentally fixed the quality/fidelity tradeoff of our method by setting the scalar weighting factors of the objective function with the following values: λ s = 1, λ d = 10 5 , λ v = 0.1 and λ p = 1. For candidate proposals in the vertex repositioning stage, γ is set to 0.1 for angle-based or Laplacian candidates.

In the connectivity optimization step, the maximum number of local topological change iterations is set to 10.

For the vertex repositioning step and its related energy minimization, we used the implementation of QPBO-P given in Rother et al. [START_REF] Rother | Optimizing binary MFRs via extended roof duality[END_REF], which is available online1 . We also performed experiments with other approximative graph cut techniques designed for nonsubmodular functions, namely classical QPBO and QPBO-I (both introduced in [START_REF] Rother | Optimizing binary MFRs via extended roof duality[END_REF]) and energy truncation to create a submodular function [START_REF] Rother | Digital tapestry[END_REF]. The results were slightly lower with these alternative techniques. Belief Propagation (BP) [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems[END_REF] is also an alternative to graph cuts techniques like QPBO-P. While both types of methods give approximate results only on this type of energy functions, their strengths and weaknesses are different. BP can be applied to non-submodular functions on multiple labels, it's complexity being O(|L| |C| ), where L is the set of labels and |C| the maximum clique size (3 in our case). The exactness of BP depends on the absence of cycles in the graph, whereas the exactness of QPBO-P depends on the submodularity of the energy function. Our choice towards QPBO-P, rather than towards BP, is motivated by the fact that QPBO-P gives an exact solution if the function is submodular, and high quality solutions if the number of non-submodular terms is low, which is the case in our problem, whereas the number of cycles is extremely high. Furthermore, nodes for which QPBO-P is able to find a solution are optimal (part of the exact solution of the problem). In our problem, at each step more than 99% of the vertex repositioning decisions are globally optimal, which gives in practice excellent results, in addition to the complexity advantages of the method. In addition, QPBO-P guarantees that the global energy does not increase after each vertex repositioning step. This last point is no more guaranteed with the BP in the presence of cycles, even the convergence of the energy minimization is not guaranteed in that case [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF].

Figure 12 shows the evolution of the global energy during 170 steps for a single CAD model.

Discussion

To illustrate the mesh quality, average triangle minimum and maximum angles are presented in table 1. Large angles cause discretization errors and large errors in interpolated derivatives while small angles are responsible for poor conditioning [START_REF] Shewchuk | What is a good linear element? Interpolation, conditioning, and quality measures[END_REF]. For high-quality meshes, the average minimum (resp. maximum) angle should be greater than 30 o (resp. less than 90 o ). The closer to 60 o the mean minimum and maximum angles are, the better are the results. According to table 1 and figure 8 and regarding the mean minimum and maximum triangle angles, regularized meshes obtained by our method are high-quality meshes. Moreover, figures 4, 5, 6 and 8 confirm this. Note that the regularized mesh vertex sampling is similar to the original mesh vertex sampling. That is convenient for regularizing meshes for which the initial vertex sampling must be preserved. Table 1 also shows that the vertex valence is slightly improved by our method, i.e. when it does not penalize neither triangle shape nor geometric fidelity. The number of irregular vertices can be significantly reduced by increasing λ v (e.g. to 0.6). However, the improvement in valence can degrade the fidelity to the original surface and may need additional global iterations to keep the mean triangle quality high.

To evaluate the surface fidelity of the remeshed models, the Hausdorff distance and the maximum of the two RMS (Root Mean Square) distances normalized to the bounding box diagonal are presented in table 1. These distances have been obtained using the Metro tool [START_REF] Cignoni | Metro: Measuring error on simplifed surfaces[END_REF]. According to these distances, to table 1 and to figure 8, the geometric error introduced by our method is small.

Our method preserves high frequency mesh features (cf. figure 7), while considerably improving triangle quality. Let's note that the number of vertices does not need to be chosen. It adapts itself to the geometry while maintaining the same order of magnitude, given the proposed setting of the vertex budget weight λ p .

We compared our method to those of Valette et al. [START_REF] Valette | Generic remeshing of 3d triangular meshes with metric-dependent discrete voronoi diagrams[END_REF], Surazhsky et al. [START_REF] Surazhsky | Explicit surface remeshing[END_REF][START_REF] Surazhsky | Isotropic remeshing of surfaces: a local parameterization approach[END_REF], and Liu et al. [START_REF] Liu | Non-iterative approach for global mesh optimization[END_REF]. As can be seen in table 1, our method gives better results in terms of triangle shape and surface fidelity when compared to Valette et al. and Liu et al.; Surazhsky et al.'s methods generate more regular triangles (better mean min and max angles), but our method better approximates the original surface. For instance, in reference [START_REF] Surazhsky | Explicit surface remeshing[END_REF] their proposed method smoothes the triceratops eye (cf. figure 7) resulting in a significant loss of details.

Our method can deal with high genus (i.e. > 1) models (cf. figure 10), and thus avoids the stitching problem that occurs in parameterization-based approaches (e.g. [START_REF] Alliez | Isotropic surface remeshing[END_REF]).

Mesh simplification

The proposed energy minimization framework allows to simplify meshes by simply changing the setting of the weight parameters, such that edge collapses are favored. This is achieved by setting the vertex penalty weight λ p to higher values (λ p = 25 in our simplification examples). Some simplified meshes are presented in figure 9. These results show a high fidelity to the original triangular surface mesh, while the number of vertices is significantly decreased. However, by favoring the removal of triangles, the average triangle quality is degraded.

In the same manner, it is possible to refine a coarse mesh by changing the vertex penalty weight λ p to negative values. To avoid varying vertex density over the refined mesh, refinement steps should alternate with regularization steps.

Conclusion and future work

In this paper we have presented a method for mesh optimization which includes robust feature detection with an improved Potts model and an original way of computing vertex positions using global combination of lo- [START_REF] Liu | Non-iterative approach for global mesh optimization[END_REF], [START_REF] Valette | Generic remeshing of 3d triangular meshes with metric-dependent discrete voronoi diagrams[END_REF], [START_REF] Surazhsky | Isotropic remeshing of surfaces: a local parameterization approach[END_REF], and [START_REF] Surazhsky | Explicit surface remeshing[END_REF]. Displayed times for Sur2 have been computed on a Pentium 4 PC (2.4 GHz) with 512 RAM [START_REF] Surazhsky | Explicit surface remeshing[END_REF], while others are from an Intel Core 2 Duo P8400 (2.26 GHz) with 4 GB RAM.

cally proposed candidates. Its main advantages are its feature sensitiveness and its ability to improve triangle shapes while preserving the original surface fidelity. The obtained results are better than other methods in terms of surface fidelity and surface fidelity/triangle quality trade-off. Our method is quite general, because by setting well-chosen weights in our objective function, a user can improve the vertex valence, improve the compactness of the representation or improve the quality of the triangles.

As future work, we will investigate quadrangular and anisotropic remeshing. We will also tackle combinatorial optimization (when computationally tractable) in the connectivity processing to improve mesh connectivity configuration. The robustness against local shape variations in the edges detection process will be improved by setting θ T h automatically (using statistics on geometrical measures) in equation 2, rather than by asking the end-user.

Fig. 2

 2 Fig.2Representation of the graph used for constructing the remeshed model: each vertex of the optimized mesh Xs is linked to its current closest initial mesh vertex Y s and has access to the whole initial mesh if needed.

Fig. 3

 3 Fig. 3 One vertex relocation step: take the best global decision for each vertex between the current (blue) and the new candidate position (orange) located in the vertex freedom sphere.

Fig. 5 Fig. 6

 56 Fig. 5 Results obtained with our method. From top to bottom: bimba, egea, triceratops. Left: input model, right: remeshed model.

Fig. 7 Fig. 8 Fig. 9 Fig. 10 Fig. 11 Fig. 12 Fig. 13

 78910111213 Fig. 7 Comparisons for the triceratops model (top) and the egea model: the original (left), our remeshed version (middle), the Surazhsky and Gotsman [11] remeshed version (right). The features surrounding the triceratops eye are well-preserved by our method.

Table 1

 1 Statistics on the remeshed models: number of vertices, percentage of irregular vertices, mean minimal angle, mean maximal angle, Hausdorff distance, maximum between the 2 RMS distances measured by Metro normalized to the bounding box diagonal, and running time. Liu, Val, Sur1, and Sur2 correspond respectively to

	Model	#v	Irreg	Amin	Amax	Er Haus	Er RM S	Time
			(%)	(deg)	(deg)	(10 -3 )	(10 -3 )	(sec)
	Fandisk (init)	6495	20	43.4	86.1	-	-	-
	Fandisk (Liu)	6495	20	44.7	82.0	3.3	0.8	n/a
	Fandisk (our)	5905	12	49.0	75.9	1.6	0.03	232
	Cow (init)	2904	53	30.2	93.7	-	-	-
	Cow (Liu)	2904	53	35.1	88.2	5.3	0.9	n/a
	Cow (our)	2695	39	41.0	81.0	5.5	0.5	59
	Shark (init)	2560	32	20.8	97.4	-	-	-
	Shark (Liu)	2560	32	26.2	107.5	3.0	0.3	n/a
	Shark (Sur1)	2560	31	50.6	71.1	6.8	0.8	n/a
	Shark (our)	1719	47	36.2	84.8	4.0	0.6	42
	Hand (init)	7950	58	32.4	94.1	-	-	-
	Hand (Liu)	7950	58	34.3	92.2	8.8	0.4	n/a
	Hand (Val)	6802	45	46.1	77.5	2.6	0.2	9
	Hand (our)	5847	33	50.2	72.3	1.7	0.2	193
	Bimba (init)	8857	62	34.2	92.8	-	-	-
	Bimba (Liu)	8857	62	38.1	87.0	4.9	0.5	n/a
	Bimba (Sur1)	8857	20	53.6	67.6	6.0	0.5	n/a
	Bimba (Val)	8143	48	45.2	78.1	6.0	0.4	10
	Bimba (our)	7986	41	47.6	75.3	3.0	0.2	232
	Egea (init)	8268	75	34.7	93.5	-	-	-
	Egea (Liu)	8268	75	38.2	88.3	2.6	0.2	n/a
	Egea (Sur2)	8705	7	52.4	69.1	2.7	0.2	15
	Egea (our)	7783	43	48.8	74.1	2.6	0.2	236
	Triceratops (init)	2832	59	29.6	95.5	-	-	-
	Triceratops (Sur2)	2758	13	42.2	82.5	8.4	1.1	12
	Triceratops (our)	2412	44	41.5	81.0	3.6	0.5	55

http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software.html

Acknowledgments We wish to thank Vitaly Surazhsky and Ligang Liu for their kind help in applying their respective remeshing methods on our 3D models and for sharing their results with us.

This work has been supported by the French National Research Agency (ANR) through the MDCO program (project MADRAS No. ANR-07-MDCO-015)