Olga Rodionova

Mohamed Sbihi

Daniel Delahaye

Marcel Mongeau

Optimization of aicraft trajectories in North Atlantic oceanic airspace

Keywords: Aircraft trajectory optimization, North Atlantic oceanic airspace (NAT), Automatic Dependent Surveillance-Broadcast (ADS-B), Simulated Annealing I

Air traffic control in non-radar oceanic airspace has always been relatively limited. Current American and European projects aim at transiting from ground-based systems to satellitebased systems. New technologies will allow to reduce significantly the present aircraft separation standards. As a consequence, the aircraft crossing the North Atlantic will be able to follow better routes, which will improve the air traffic situation by decreasing of flight durations and congestion in pre-oceanic continental airspace. In this paper an optimization model is introduced and a simple methodology to solve this problem is proposed that yields encouraging preliminary results.

INTRODUCTION

A. North Atlantic airspace

The North Atlantic (NAT) is the busiest oceanic airspace in the world. For the most part in the North Atlantic, Air Traffic Services (ATS) radar surveillance is unavailable, and Air Traffic Control (ATC) must rely significantly upon the aircraft HF voice position reports. Aircraft separation assurance and safety are nevertheless highly demanded.

Due to the constraints of large horizontal separation criteria and a limited economical flight level (FL) band (FL310-400), the airspace is congested at peak hours. In order to provide the best traffic control service, a system of tracks referred to as the Organized Track System (OTS) is constructed to accommodate as many flights as possible on, or close to, their minimum time tracks and altitude profiles [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF].

As a result of passenger demand and time zone differences, much of the NAT air traffic contributes to two major alternating flows: westbound and eastbound, whose effect is to concentrate most of the traffic unidirectionally. Due to the meteorological conditions, including the presence of jet streams, consecutive eastbound and westbound minimum-time tracks are seldom identical. The creation of a different OTS is therefore necessary for each of the major flows. In general the east-west tracks are situated more northerly than the west-east tracks (Figure 1). Currently, about half of NAT flights are planned on the OTS [START_REF]Guidance concerning air navigation in and above the North Atlantic MNPS airspace[END_REF]. OTS flights operate on great circle tracks joining successive significant waypoints and are commonly planned so that the specified ten degrees of longitudes (10°W, 20°W, etc.) are crossed at integer degrees of latitude. The aircraft may be planned at any of the flight levels published for that track. Unless otherwise requested by ATC, the flights on tracks should report their position at the designated reporting points listed in the flight plan. To mitigate against communication equipment failures and poor propagation conditions, ATC often employs strategic traffic planning, and issues Oceanic Clearances (authorization to follow the specified track at a defined speed and flight level). Flights that continue to follow such pre-coordinated strategic oceanic clearances are thereby guaranteed conflictfree progress to oceanic exit. Every effort is therefore made to clear aircraft as per their flight plans. However, this is not always possible, especially during peak traffic flow periods.

B. Separation standards

The current situation shows that the traffic often follows a routing which is not optimal in terms of flow management from an overall ATC point of view [START_REF]ASEP-ITM simulations from traffic data[END_REF], and therefore some rerouting is necessary when entering continental radar control areas. This leads to increasing length and duration of the flights as well as to increasing the congestion in the continental airspace. It seams reasonable to enable the re-routing of aircraft directly within the OTS-covered area from the initial predefined track to another oceanic exit point that is more suitable in view of its final destination. This would reduce ATC complexity in continental airspace near the exit of the track system, and improve the aircraft's remaining route towards destination.

A major factor in preventing traffic re-routing with current oceanic separation standards is the additional longitudinal separation required for flights which have not reported over a common point. Two consecutive aircraft following the same track should be separated 10 minutes apart (Figure 2), as their longitudinal relationship is established by their position reports and any errors in forward position estimates can be assumed to cancel out since they both experience the same weather. But this assumption can not be applied in the case where a particular aircraft desires to re-route onto the adjacent track (Figure 2), as it has experienced weather different from that of any aircraft following this track. As a consequence, current regulations impose an increased longitudinal separation of 15 minutes in this case. Due to the high density of traffic on tracks, this re-routing maneuver can therefore rarely be applied. These and other data are then broadcast to other aircraft and to ATC (Figure 3).

The new approach will allow to decrease significantly the oceanic separation standards [START_REF] Massoglia | The Use of Satellite Technology for Oceanic Air Traffic Control[END_REF]. Indeed, with ADS-B systems the consecutive aircraft following the same track would be separated only 2 minutes apart, and the aircraft performing a re-routing to the adjacent track would be separated only 3 minutes from the aircraft on its new track. Obviously, the new separation standards will raise the limit on the total number of flights allowed for crossing the NAT airspace, and generally improve the aircraft routes by decreasing their length.

Several papers [START_REF] Mohleji | Performance Analysis of North Pacific Operations Using an Automated Air Traffic Control System Simulation[END_REF][START_REF] Grabbe | Central East Pacific Flight Routing[END_REF][START_REF] Grabbe | Central East Pacific Flight Scheduling. AIAA Guidance, Navigation, and Control Conference and Exhibition[END_REF][START_REF] Grabbe | Scheduling Wind-Optimal Central East Pacific Flights[END_REF] are devoted to treating the related problem in the Pacific oceanic airspace. However, to our knowledge, no research works have been published concerning the NAT airspace. This paper describes a preliminary study aiming at modeling the improvement of oceanic air traffic situation when using the ADS-B technology.

II. PROBLEM FORMULATION

A. OTS modeling

Our aim is to provide a mathematical formulation of the problem consisting in searching for optimal flight routes for a set of flights within OTS. As the subject of the study is the NAT airspace, only the parts of the routes belonging to the OTS will be taken into account. Moreover, as the eastbound and westbound traffics are well separated in time and space, they can be considered independently. Thus, only the eastbound traffic (departing from North America and arriving in Europe) will be taken into consideration here. Obviously, the westbound traffic can be treated in the same way.

To define the flight route, first it is necessary to define the OTS structure. The OTS can be represented by a Nx×Ny×Nz grid of waypoints, where Ny is the number of OTS tracks, Nxthe number of waypoints on each track, and Nzthe number of flight levels for each track. Thus, the 3D position of an eastbound aircraft located on track j at waypoint i at flight level k is completely specified by the vector (i, j, k).

A flight entering a predefined track at a predefined flight level is supposed to follow this track at the same level until a change of trajectory is made. In this preliminary study such changes are only allowed at the waypoints. Thus, from its current 3D position (i, j, k), the flight has several possibilities to pursue its route:

 to follow the same track at the same flight level,  to re-route to an adjacent track at the same flight level,  to change the flight level.

In the first two cases, the flight should follow the straight line connecting the current waypoint to the next one, where the next waypoint can be one of the following (see Figure 4):

 (i+1, j, k) (same track);  (i+1, j+1, k) (northern reroute);  (i+1, j-1, k) (southern reroute).
This grid of waypoints represents a network for which nodes are waypoints and links are route segments connecting such nodes (Figure 4). Some extra nodes have been added for modeling route segment crossings. These extra nodes are considered to be outside the network.

When an aircraft changes its flight level, it is only authorized to climb (not to descend), in order to satisfy as much as possible the optimal kerosene consumption flight profile provided by the company. Moreover, according to the current ATC rules, aircraft may change flight level only when crossing a waypoint (Figure 5). In our model the distance between the real horizontal position of the aircraft at the new flight level and the previous waypoint (see Figure 5) is neglected, as well as the time required to reach the new flight level. This instantaneous-climbing hypothesis is reasonable as the distance between the waypoints is much greater than the distance needed for flight level changing, so the aircraft is supposed to climb almost instantly. Thus, the next aircraft position in this case will be assumed to be (i, j, k+1). The number of options given to a flight depends on its entry and exit track numbers. For instance, if the entry and exit tracks are the same, the aircraft has no opportunity to change its trajectory. The flight-level changes are fixed according to the optimal performance of the aircraft and happen at some computed waypoint. They are not decision variables.

In this paper a simplified model has been considered. For each flight f the following input data are given:

     Ny f Tr in ,..., 2 , 1 
-the desired entry track;

     Ny f Tr out ,..., 2 , 1 
-the desired exit track;

   f t in
-the desired entry time;

     f const f v  -the aircraft speed;      Nz f FL in ,..., 2 , 1 
-the entry flight level;

     Nz f FL out ,..., 2 , 1  -the exit flight level     f FL f FL in out  ;    1 ,..., 2 , 1 ,   Nx i f z i
-binary parameters defining the flight altitude profile, that depends on the aircraft weight, and satisfying the following conditions:

      otherwise; , 0 , t at waypoin FL its changes flight if , 1 i f f z i       f FL f FL f z in out Nx i i      1 1 .
For each flight f, we define decision variables as follows:

   1 ,..., 2 , 1 ,   Nx i f x i
-binary variables defining the flight re-routing maneuvers:

      otherwise. , 0 , t at waypoin track changes flight if , 1 i f f x i
For each flight f, these variables form a binary vector of size 1  Nx , where the number of ones is equal to the distance between actual entry and exit tracks. When

  1  f x i
, the aircraft leaves the current occupied track at waypoint i, and reroutes to the adjacent track (the next track towards the exit track).

Based on our first experiments, we have remarked that, for many situations, such a search space definition does not guarantee existence of conflict-free trajectories. In order to avoid this situation, entry and exit track numbers and entry time have been relaxed by allowing aircraft to enter (or to exit) an adjacent track with some delay. The delay (denoted as   Decision variables must satisfy the following constraints:

           1 , , 1 _    f Tr f Tr f Tr f Tr in in in act in
, (tolerance with regards to required entry track);

           1 , , 1 _    f Tr f Tr f Tr f Tr out out out act out
, (tolerance with regards to required entry track);

     slot f f d    , where     d N f ,..., 1 , 0   , (acceptable delays);        f Tr f Tr f x act in act out Nx i i _ _ 1 1      ,
(total number of re-routing maneuvers).

C. Objective function

For a single flight, different route optimality criteria can be chosen, e.g. the total trajectory length, the total flight duration, the fuel consumption, etc. While constructing the optimal routes for a set of flights, we must ensure the additional collision avoidance requirement based on the separation standards.

In this preliminary study, only a collision-avoidance criterion is taken into account in order to simplify modeling, as the aim of this first step of the study is to verify whether a conflict-free solution exists when applying the new separation standards. In future work alternative measure characteristics (such as fuel consumption, flight duration, etc.) will be considered.

Based on the route network structure, conflicts may happen only at nodes and links. The following auxiliary variables are introduced:

 n
Cf -the number of conflicts on nodes;  l Cf -the number of conflicts on links. Thus, in the current model, the objective function is simply given by  

l n Cf  .
Given an instantiation of decision variables, the aircraft trajectories can be computed on the network based on a flight simulator. Each time an aircraft passes over a node, the passing time is recorded. The same recording process is applied for link entry and exit times.

A conflict is detected on a node if the longitudinal separation constraint is violated for two aircraft passing this node. To calculate the number of conflicts at a given node n, all the flights f passing this node are selected and sorted according to the time ) ( f t n at which they pass this node. To satisfy the longitudinal separation requirements, the condition

      1 , 1      f f t f t f t s n n must be fulfilled for 1 ,..., 2 , 1   n N f
, where n N is the total number of flights passing the node, and

  1 ,   f f t s
is the imposed longitudinal separation value depending on the flights' maneuvers. If for any pair of flights f and f +1 this condition is violated, then a conflict is detected and n Cf is increased by one.

Based on the sequence order of aircraft at the link entry and exit, it is possible to detect conflict happening on this link. It happens if an aircraft is slower than the one following it on the same track. To compute the number of conflicts at a particular link l, all flights passing this link are selected and sorted into two lists according to the time at which they entry and exit this link. Then, the entry and exit lists are compared. If two aircraft are found swapped, then a conflict is detected and l Cf is increased by the number of flights having different ranks in both lists.

III. OPTIMIZATION METHOD

The problem we have to solve is fully discrete with many decision variables. If we consider N flights with a number of discrete delays

  N c d T Nx N S                     1 1 . For instance, if 5  d N , 4  c T , 10  Nx and 500  N
, then the number of different potential solutions to be considered is:

500 500 756 4 9 6                    S .
Due to this high combinatorics, we have developed a methodology based on stochastic optimization and implemented the Simulated Annealing (SA) method. This optimization method is a metaheuristic inspired from the thermodynamics theory [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF][START_REF] Aarts | Simulated Annealing and Boltzmann Machine[END_REF]. It imitates the annealing of the metal, involving heating and controlled cooling. Each point s of the search space is analogous to a state of some physical system, and a function E(s) to be minimized is analogous to the internal energy of the system in that state. The goal is to bring the system, from an arbitrary initial state, to a state with minimal energy (the optimal solution of the problem).

At each step, the SA heuristic considers some neighboring state (solution) s' of the current state (solution) s, and probabilistically decides between moving the system to the new state s' or remaining in state s. This probability depends both on the difference between the corresponding function values E(s) and E(s') and also on a global parameter T (called the temperature) that is gradually decreased during the process. The dependency is such that the choice between the previous and current solution is almost random when T is large (during the first iterations), but increasingly selects the better solution as T goes to zero (as the number of iterations grows). Typically, this step is repeated until the system reaches a state that is considered good enough, or until a given computation budget has been exhausted.  the number of transitions performed at a same temperature level;

 the stopping criterion for the algorithm.

In our context, the state of the system will represent the trajectories of the set of N flights using N vectors, each of which corresponds to a particular flight f and contains the decision variables defined above:

  f  ,   f Tr act in _ ,   f Tr act out _ ,   1 ,..., 2 , 1 ,   Nx i f x i
(see Figure 6). The flight trajectory is completely defined by these variables together with the flight predefined data such as the desired track entry time, the flight altitude profile and the aircraft velocity. The energy function is represented by our objective function: the total number of conflicts,  

l n Cf Cf s E   .
To generate a neighboring state ' s from the current state s, the mutation operator has to be defined. In the current model, a mutation consists in choosing randomly one flight f, and changing randomly some of its parameters. The mutation of rerouting variables involves choosing randomly two variables

  f x i and   f x j
having different values (0 and 1), and permuting their values (Figure 7). The goal is to explore widely the search space at high temperatures, and to reduce the search area as the temperature decreases in order to concentrate the search near the optimum. The probability of acceptance of the state transition is defined by the following formula:

                        s E s E e s E s E T s s pr T s E s E accept ' if , ; ' if , 1 , ' ' .
The initial temperature 0 T is to be adjusted from each particular problem in order to provide acceptance of most of the transitions. The process of the initial temperature setting is analogue to the process of the system heating. It starts with some small temperature value 0 0 T that can be obtained from examining a number M of system states. Then, at each step, the current temperature is increased by multiplying it with some value 1   :

1 0 0    i i T T



. These steps are repeated until the acceptance probability at the current temperature becomes sufficient (

     T s s pr accept , '
, for some user-defined value of  ). The last obtained temperature is then chosen as the initial temperature 0 T . Figure 6. System state representation The process of temperature decreasing is analogue to the system cooling. The cooling should be made rather slow in order to allow the system to reach the thermodynamic equilibrium. The exponential cooling schema has been chosen:

1    i i T T 
, where  is close to but smaller than 1. The number of transitions t N tested at each particular temperature level was supposed to be constant.

The stopping criterion is to stop either when an optimal (i.e., feasible) solution is found (   0

 s E
), or when the temperature goes below the predefined critical value f T .

IV. COMPUTATIONAL RESULTS

This section represents the experimental settings and obtained results in order to validate the proposed optimization formulation and the methodology introduced.

A. User-defined parameter values for SA

After the number of tests having been held, the following parameter value were adjusted empirically for SA algorithm:

 1000  M ;        M i i s E M T

B. Computational environment

The SA algorithm was implemented in Java and was run under Windows-32 operational system, on Intel® Core™ 2 CPU with 1.73 GHz.

C. Test problem description

To test our algorithm five random aircraft flight sets were generated. An aircraft flight set consists of 500 flights defined by the values The tests were performed with two longitudinal separation values: 10 minutes (for the set of aircraft not equipped with ADS-B) and 2 minutes (for the set of aircraft equipped with ADS-B). In the case where an aircraft re-routes from another track, the value of longitudinal separation at the targeted waypoint

  f Tr in ,   f Tr out ,   f t in ,   f v ,   f FL in ,   f FL out ,   f z i , for 500
  1 ,   f f t s
is multiplied by 1.5 as required by the standard (15 minutes and 3 minutes correspondingly). In the case where one of the consecutive aircraft changes its flight level, the value of longitudinal separation between these aircraft

  1 ,   f f t s
is multiplied by 1.1, in order to compensate the error caused by neglecting the climbing time.

D. Results

Two types of tests were performed on the same generated flight sets. In the first case, it was supposed that the aircraft do not possess the ADS-B system. In the second case, all the aircraft were considered to be equipped with ADS-B.

For the flight sets not equipped with ADS-B it was shown that no collision-free solution exists. This confirms that aircraft have to fly the routes that are not optimal. Different values of the temperature decrease parameter  were tested in this case in order to obtain better solution. The algorithm is able to decrease the total number of conflicts by a factor 3 when compared with the initial number. The results of simulations for flight set 1 are represented in Table I and in Figure 8. The initial number of conflicts in this case was equal to 2775. For the flight sets where all aircraft are equipped with ADS-B, the SA algorithm finds an optimal (collision-free) solution for every test set under study. The average number of iterations in heating for the 5 tests is 11. The maximum number of decreasing steps during the cooling process that the algorithm can perform with the temperature decrease parameter 9 . 0   is 88. For the tests with ADS-B, this number is sufficient, the value 9 . 0   is therefore preferred to be used. In this case, the average number of iterations in cooling for the Table II represents the average computational time of algorithm execution for both types of performed test (with and without ADS-B). V. CONCLUSION In this paper, a developed mathematical model of the oceanic air traffic is introduced. Based on the future performances of the ADS-B system, an optimization algorithm has been implemented in order to improve aircraft routes and to reduce the induced congestion at the exit continental airspace. Due to the high complexity of the combinatorial optimization problem, a Simulated Annealing algorithm has been applied.

The performed simulations prove that implementing the ADS-B technology has a strong positive effect on the current traffic situation in North Atlantic oceanic airspace. The reduction of the current separation standards makes it possible for the aircraft to perform re-routing within Organized Track System and therefore, to follow more optimal trajectories towards their destination. As a consequence, the total flight duration as well as the congestion level in the pre-oceanic airspace will significantly decrease.

In the future work it is planned to test the developed approach on real NAT traffic data and to compare the obtained optimized trajectories with the routes actually followed by aircraft, using appropriate aviation-focused metrics. Also it is planned to modify the objective function used in the model in order to take into account other performance measures (such as fuel consumption, flight duration, etc.). Future work could also consider other optimization algorithms to solve the described problem.

Figure 1 .

 1 Figure 1. North Atlantic oceanic airspace with Organized Track System

Figure 2 .

 2 Figure 2. Aircraft longitudinal separation criterion on the tracks The transition from present ATC tools to airborne-based systems and procedures proposed by the American and European projects NextGen (Next Generation Air Transport System) and SESAR (Single European Sky ATM Research) can help to overcome these drawbacks. The key component of these projects is the Automatic Dependent Surveillance-Broadcast (ADS-B) technology [3] that combines aircraft's positioning source, aircraft avionics, and a ground infrastructure to create an accurate surveillance interface between aircraft and ATC. Aircraft transponders receive GPS signals and use them to determine the aircraft's actual position.These and other data are then broadcast to other aircraft and to ATC (Figure3).

Figure 3 .

 3 Figure 3. Automatic Dependent Surveilllance-Broadcast operating principle

Figure 4 .Figure 5 .

 45 Figure 4. OTS grid model with nodes and links

  chosen among a number d N of discrete-valued multiples of a fixed slot duration (denoted as slot). The new associated decision variables are: time delay at track entry.



  Several notions and parameters are to be defined and adjusted for the SA algorithm:  the representation of state of the system (solution) s from the search space;  the energy function (the objective function of the minimization problem)   s E ;  the method to generate the neighboring state (solution) from the current one ( ' s s  );  the probability of acceptance of the state transition as a function of the temperature the law of temperature decreasing;

Figure 7 .

 7 Figure 7. System state mutation to perfrom transition to a neighbor state

  randomly generated on a 4-hour time period with realistic altitude profiles and velocities. The aircraft velocity   f v was chosen in the range from 450 knots to 500 knots. Aircraft were allowed to select an entry time delay in the range {0, 5, 10, 15, 20, 25, 30} minutes.The OTS was modeled with a regular grid consisting of 7 Nz flight levels. The distance between the tracks is 1º, which corresponds to 60 NM (Nautical Miles). The distance between waypoints is 10º, i.e. 300 NM at the latitude of 60º (taken as an example). The distance between the flight levels is 1000 feet (304.8 m).

Figure 9 .Figure 10 .

 910 Figure 9. Evolution of the number of conflicts with transitions performed during the process of cooling for flight set 1 with ADS-B

TABLE I .

 I COMPARING THE RESULTS OF SA FOR A TEST FLIGHT SET WITHOUT ADS-B 

		Number of iterations	Best solution, final
		of cooling	number of conflicts
	0.9	88	1135
	0.95	180	1068
	0.99	917	817

Figure 8. Evolution of the number of conflicts with the steps of cooling process performed for temperature decreasing for flight set 1 without ADS-B

tests is 37. The following graphs (Figures9 and 10) represent the algorithm progress.

ACKNOWLEDGMENT This work has been supported by French National Research Agency (ANR) through COSINUS program (project ID4CS n°ANR-09-COSI-005).


Cyril Allignol
email: allignol@recherche.enac.fr

Nicolas Barnier
email: barnier@recherche.enac.fr

Alexandre Gondran
email: gondran@recherche.enac.fr

Optimized Flight Level Allocation at the Continental Scale

Keywords: Deconfliction, Flight Level Allocation, Graph Coloring, Tabu Search

As acknowledged by the SESAR program, current ATC systems must be drastically improved to accommodate the predicted traffic growth in Europe. This study aims at assessing the performance of 4Dtrajectory planning and strategic deconfliction, two of the key concepts identified to meet SESAR capacity objectives.

Among the possible degrees of freedom, we focus here on the flight level (FL) optimization to avoid en-route conflicts between intersecting flights. The resulting problem can be reduced to Graph Coloring with a specific cost function minimizing the discrepancies to requested FLs.

This FL allocation leads to very large combinatorial optimization problems when applied at the continental scale, especially when considering temporal uncertainties. The instances were solved with a Tabu Search algorithm in a few seconds to a few minutes, depending on size and conflict density. Our results shows that the global conflict resolution workload is alleviated by at least 20 %, while bounding the individual FL discrepancies to three levels for a small proportion of the traffic.

I. INTRODUCTION

In an already saturated European sky, the predicted growth of air traffic volume (+50 % in 2025 compared to 2011 according to [START_REF] Statfor | Eurocontrol long-term forecast: Ifr flight movements 2010-2030[END_REF]) urges to improve Air Traffic Management (ATM) efficiency, as attested by the European Single Sky program SESAR. Current ATM optimization strategies, like reducing the size of control sectors or the distance of separation (RVSM, P-RNAV), seem to have reached the structural limits of the system, while the automation of Air Traffic Control (ATC) has known few significant improvements over the last decades [START_REF] Garot | Failures in the automation of air traffic control[END_REF].

Among the key concepts identified to meet SESAR performance objectives, the planning of 4D-trajectories would allow to increase en-route capacity, while preserving the current level of safety. Several degrees of freedom on the trajectories can be exploited to regulate the traffic in order to make the number of potential conflicts decrease and be able to improve the airspace capacity. Previous studies have considered re-routing [START_REF] Sorensen | Limited dynamic re-routing (LDRR)[END_REF], [START_REF] Gianazza | Assessment of the 3dseparation of air traffic flows[END_REF], ground holding to satisfy sector capacities [START_REF] Dalichampt | Innovative slot allocation (ISA)[END_REF], [START_REF] Barnier | Slot allocation with constraint programming: Models and results[END_REF] or to deconflict the traffic [START_REF] Barnier | 4D-Trajectory deconfliction through departure time adjustment[END_REF], but few are dedicated to the optimization of the cruise flight level (FL) to vertically separate intersecting flows [START_REF] Letrouit | Optimisation du réseau de routes aériennes en Europe[END_REF], [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF], avoiding any potential conflict between two aircraft of distinct flows. Moreover, recent studies [START_REF] Allignol | Planification de trajectoires pour l'optimisation du trafic aérien[END_REF] have assessed that more than 30 % of potential conflicts occur between pair of flights that are both in their cruise phase (for traffic data in the French airspace), so that an optimized FL allocation could significantly lower the conflict resolution workload of controllers.

In the current tactical regulation system managed by the Central Flow Management Unit (CFMU) of Eurocontrol, hardly any optimization is attempted 1 on the cruise Flight Levels. So, with each flight plan filed by airlines, a FL depending on the type of aircraft and the length of the trajectory considered is requested and then granted by Eurocontrol. As the performances of aircraft for the same range are likely to be similar, the same FLs are requested over and over, which generates high traffic loads on specific layers of airspace. However, the Requested Flight Level (RFL) of a flight plan usually corresponds to an optimal cruise altitude w.r.t. fuel consumption or CO 2 emission. Thus the distance between the allocated flight level and the RFL should be kept as low as possible, should such a scheme be accepted by airlines.

Of course, flow-based FL allocation schemes, which approximate trajectories in the horizontal plane without taking the climb or descent phases into account and gathers flights sharing the same origin, destination (and possibly RFL) only prevent conflicts occurring between distinct leveled flows in their cruise phase. So conflicts involving at least one vertically evolving trajectory, or catch up flights belonging to the same flow, will still have to be resolved by other means (ground holding, ATC actions...). Furthermore, [START_REF] Allignol | Planification de trajectoires pour l'optimisation du trafic aérien[END_REF] has shown that such a flow-based FL allocation corresponds to very dense conflict graphs and tends to generate very costly solution requiring large discrepancies from the RFLs.

To partially overcome these issues, we propose a FL allocation scheme considering each flight individually and taking into account the temporal dimension of conflicts, whereas previous works considered any given route corresponding to a flow as occupied all day long or at least during the whole flight time of each aircraft of the flow. Such a scheme allows to take catch up conflicts into account and generates much more sparse conflict graphs, depending on the uncertainties assumed on the temporal following of the trajectory.

Several optimization paradigms have been evaluated to solve the corresponding allocation problem: evolutionary algorithm, tabu search and constraint programming. We will focus here on the tabu approach as it appeared to consistently outperform the two others in our experiments, producing solutions with maximal discrepancies of only two or three FL for a limited proportion of the traffic (typically 30 %) in a few minutes of computation time.

In the following sections, we first briefly introduce the current ATC and Air Traffic Flow Management (ATFM) system in Europe and review the related research projects on FL allocation. Then we describe our conflict model and the resulting FL allocation combinatorial optimization problem taking each aircraft individually into account. Next, our results on instances of expanded sizes over one day of traffic in the ECAC 2 airspace, with varying amounts of temporal uncertainties, are presented. We end with planned further works to enhance the approach before concluding.

II. CONTEXT AND RELATED WORKS

The next section briefly describes current ATFM policies applied by Eurocontrol in the ECAC airspace, through its operational unit, the CFMU, which centralizes the filing of flight plans and issues ATFM regulations to prevent the overload of control sectors. However, no optimization is performed on the RFLs during this process, leading to saturated airspace layers. 2 European Civil Aviation Conference

In the following section, we review the few research projects that have attempted to address this issue with flow-based models in order to improve airspace capacity, and compare them with our approach.

A. ATC and ATFM

Air Traffic Control (ATC) is a ground-based service provided to ensure the safety and efficiency of the flow of aircraft. The first goal of ATC is to maintain aircraft separated: outside Terminal Areas (TMA) around airports (where specific procedures are applied), two aircraft should remain distant from each other at least by 5 NM horizontally and 1000 ft vertically, as illustrated by the safety volume of figure 1. Warning System). Our FL optimization approach would be integrated during the pre-tactical phase, with possible updates a few hours before takeoff, as the computation times of our resolution algorithm would be compatible with the current practice of the CFMU.

B. Research Projects on Flight Level Allocation

Even if CFMU experts can balance traffic load by separating flights of the same flows over alternative routes, no real optimization of the FL allocation is currently performed at the strategic or pre-tactical level in Europe, except for a few static rules attributing the parity of FL according to the heading of flights to prevent face-to-face conflicts. At the real time level, Air Traffic controllers can issue temporary FL changes to separate aircraft, and this kind of vertical manoeuvre has been integrated in the CATS traffic simulator and conflict solver [START_REF] Granger | Optimal resolution of en route conflicts[END_REF] developed at the French Civil Aviation R&D centre (DSNA/DTI).

Our work is more similar to the flow-based approaches presented in [START_REF] Letrouit | Optimisation du réseau de routes aériennes en Europe[END_REF], [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF] where graph coloring techniques are applied to optimize the FL allocation of direct routes for the French or European traffic, or in [START_REF] Gianazza | Assessment of the 3dseparation of air traffic flows[END_REF] where a Genetic Algorithm and an A ⋆ algorithm are used to optimize traffic flows with possible changes of FL along the route (as well as lateral deviations for direct routes). However, as our optimization scheme allocate a FL for each flight individually and take the temporal dimension accurately into account, our instances are much larger and the conflict graph much sparser.

Furthermore, contrary to these studies, our model does not fit well within standard graph coloring problems as the choice of FL is very much restricted for operational reasons (detailed in section III-B), so each flight will only have a couple of possible FL available above or below its RFL. Moreover, the FL allocation phase should handle over-constrained instances as well to obtain the best possible solution, even if some conflicts remain (which will have to be solved by other ATC or ATFM actions). Thus, classic coloring techniques providing lower bounds on the cost of solutions like the use of cliques or maximal flow algorithms (cf. the "alldifferent" constraints used in [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF]) cannot be used to speed up the search. Furthermore, the cost of solutions is measured w.r.t. the number of remaining conflicts and the sum of discrepancies from RFLs, the latter not being a conventional criteria for graph coloring problems and algorithms. Consequently, our tabu search algorithm has been tailored (as described in section IV) to handle such variations on the standard graph coloring problem, whose main optimization criterion is the number of colors (i.e. FLs).

III. DECONFLICTION WITH FLIGHT LEVEL ALLOCATION

The flight level allocation is aimed at vertically separating intersecting aircraft, in order to reduce the complexity of the traffic. This FL allocation is performed according to the following steps:

• computation of intersections between aircraft in the horizontal plane; • allocation of a FL to each flight in such a way that two intersecting flights are given a different FL (if possible). This allocation does not take the climbing or descending parts of trajectories into account, so that only conflicts occurring in the cruise phase will be solved, which accounts for approximately 25 % of the total number of conflicts in our experiments. The remaining conflicts can be avoided with other strategic methods such as the ground holding scheme presented in [START_REF] Barnier | 4D-Trajectory deconfliction through departure time adjustment[END_REF], or with standard ATC procedures used by controllers.

A. Conflict Detection

Our data are provided by the CATS3 simulator [START_REF] Alliot | CATS: A complete air traffic simulator[END_REF], which takes as input all filed flight plans related to the ECAC airspace for a given day of traffic and the relevant airspace configuration. It outputs 4D-trajectories sampled every 15 s, which is the largest time interval that guarantees the detection of the shortest possible conflict (i.e. facing aircraft at maximum speed).

To compute the constraints of the model, trajectories are pairwise probed for conflicting points. Given a flight i, we note:

• {p k i } the chronologically ordered sequence of the 3D-points of its trajectory; • t k i the time at which flight i will be at point p k i . For any geometrically conflicting points p k i and p l j such that the separation norm is violated (d h being the distance in the horizontal plane and d v in the vertical plane):

d h (p k i , p l j ) < 5 NM and d v (p k i , p l j )
< 1000 ft we must temporally ensure that t k i = t l j , which can be rewritten: 0 = t k it l j . If we consider the whole conflicting zone started at points p k i and p l j , we obtain the constraint:

0 ∈ [lb ij , ub ij ]
where lb ij and ub ij are respectively the lower and upper bounds of the t k it l j for the considered conflict.

A given pair of flights i and j may conflict several disjoint times over their trajectories, thus leading to several such intervals. The conflict constraint for flights i and j, conflicting σ times, is therefore:

0 ∈ [lb 1 ij , ub 1 ij ] ∪ • • • ∪ [lb σ ij , ub σ ij ] (1)
We note

C ij = [lb 1 ij , ub 1 ij ] ∪ • • • ∪ [lb σ ij , ub σ ij ].
The resulting conflict model is exact: a violated constraint necessarily entails a conflict during the simulation and, conversely, any conflict detected between two trajectories in the pre-processing phase will yield a C ij constraint. More details about this model can be found in [START_REF] Barnier | Combining flight level allocation with ground holding to optimize 4D-deconfliction[END_REF].

However, aircraft might not follow their 4D-trajectory precisely, so that the above conflict detection shall be carried out with the real transit time at point p k i :

θ k i = t k i + ε k i
where ε k i is a temporal uncertainty. The constraint for flights i and j becomes:

[-ε max , ε max ] ∩ C ij = ∅ (2)
where ε max is the maximum amplitude of temporal uncertainties: ∀i, k,

|ε k i | ≤ εmax 2 .
For example, ε max = 6 states that all aircraft are able to follow their 4Dtrajectory within a ±3 min uncertainty on each point. Such an accuracy is out of reach of the current flight management systems, but is within the performance objectives of the SESAR program.

B. Model

Standard models for flight level allocation, such as in [START_REF] Letrouit | Optimisation du réseau de routes aériennes en Europe[END_REF] or [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF], aggregate flights with the same route and possibly the same RFL into flows. FL allocation is then performed considering intersections among flows. However, we propose in this study to allocate FL to flights individually. This new model generates instances of higher dimensions but with a lower density conflict graph, leading to better solutions, i.e. with lower deviation from RFL, as depicted on figure 2 which presents the distribution of discrepancies from RFLs for several days of traffic in the French airspace (optimized by a Constraint Programming algorithm, see [START_REF] Allignol | Planification de trajectoires pour l'optimisation du trafic aérien[END_REF]).

The variables of the model are the flight levels FL i for each flight i. In order to keep the operational cost of the allocation as low as possible, flight levels are bounded by a tight interval around RFLs -typically 2 or 3 FL above or below -which induces acceptable over-consumption of a few percent for airlines (see [START_REF] Airbus | Getting to grips with fuel economy[END_REF]). If we note RFL i the RFL for flight i, we ensure that:

FL i ∈ [RFL i -∆ max , RFL i + ∆ max ] (3)
where ∆ max is an instance parameter defining the maximum deviation from RFL. Additionally, due to aircraft performances, each FL i is bound by a value FL max i corresponding to the aircraft's ceiling, so that:

FL i ∈ [RFL i -∆ max , min(RFL i + ∆ max , FL max i )]
Constraints for this model are directly derived from the conflicts presented in section III-A, i.e. whenever two flights i and j are in conflict (w.r.t. equation 1), they must be assigned to different flight levels:

0 ∈ C ij ⇒ FL i = FL j (4)
As aforementioned, the antecedent of the implication 4 can be replaced by equation 2 to take a given amount of temporal uncertainty into account, thus tightening the density of the conflict graph.

To solve this problem, the primary aim of our search algorithm is to satisfy the greatest possible number of constraints. Then, it tries to minimize a secondary cost that takes airlines preferences into account as described in the next section.

C. Cost of the FL Allocation

When registering a flight plan, airlines must indicate at which FL they wish to fly. This requested flight level usually corresponds to the optimal cruise altitude w.r.t. to the operational cost of the airlines. Flying at a much higher or lower FL might lead to a significant increase in fuel consumption or even be above the aircraft ceiling, thus allocated flight levels must be chosen as close as possible to the RFL.

In order to evaluate the performance of the FL allocation, we define a cost function that measures the deviations from RFL:

cost FL = i |FL i -RFL i |
It would also be possible, given a realistic performance model for aircraft, to define a cost function that takes fuel consumption (or CO 2 emission) into account:

cost consumption = i c i (FL i ) -c i (RFL i )
where c i gives the operational cost for flight i (depending on aircraft type, travel distance, load factor...) as a function of the FL.

IV. RESOLUTION WITH TABU SEARCH

The flight allocation problem described in the previous section can be considered as a Graph Coloring Problem (GCP): each flight corresponds to a vertex, each conflict detected between two flights (cf. section III-A) to an edge and two connected vertices must be allocated different FLs, which correspond to colors. Usually, the only added constraint is a maximal number of colors k and the problem is then named k-GCP. In our case, the number of FLs is limited for each vertex by operational constraints (cf. section III-B): each flight must be colored with a very tight FL range (compared to the global domain of FLs) around the RFL (relation 3). This problem is known as list-coloring problem.

We define in the following sections our slightly modified version of list-coloring, then we detail our adaptation of the Tabucol algorithm to solve the FL allocation problem.

A. List-Coloring

Given a graph G(V, E) with vertex set V and edge set E, and, for each vertex u, a color domain D u , a listcoloring of G is a function that allocates to each vertex u ∈ V a color c ∈ D u . If two adjacent vertices have the same color, the two vertices and the corresponding edge are said conflicting. A solution to our problem is a list-coloring without conflicting edges, named a legal solution, that minimizes the discrepancies from RFLs: minimize

S FLi∈S |FL i -RFL i | subject to FL i ∈ D FLi , ∀FL i ∈ S F L i = F L j , ∀(FL i , FL j ) ∈ E
where S = {FL i , ∀i}, E is the set of conflicting pairs of flights according to equation 2 and

D FLi = [RFL i - ∆ max , RFL i + ∆ max ]
as explained in relation 3.

B. Adaptation to the FL Allocation Problem

To solve our problem, we have adapted the wellknown Tabucol algorithm [START_REF] Hertz | Using tabu search techniques for graph coloring[END_REF], an instance of tabu search [START_REF] Glover | Tabu search -part i[END_REF] tailored for k-GCP. [START_REF] Galinier | A survey of local search methods for graph coloring[END_REF] improves this algorithm and presents a very good survey of local search methods for graph coloring. As the search space in Tabucol contains all k-colorings, legal or not, the domain constraint of our optimization problem must be relaxed in the objective function.

To adapt Tabucol to our non-conventional cost function, we modify the objective as follows:

minimize S   (FLi,FLj)∈E δ FLiFLj , FLi∈S |FL i -RFL i |   subject to FL i ∈ D FLi ∀FL i ∈ S
δ FLiFLj being the Kronecker function which returns 1 if FL i = FL j and 0 otherwise. This tailored strategy features two objectives which are considered in lexicographic order. The first one counts the number of conflicting edges in a list-coloring and the second is the sum of discrepancies from RFL.

This modification of the objective allows to find solutions that minimize the number of conflicts -which is the most important goal in practice -even if there is no legal list-coloring (i.e. the first objective equals to zero) or if Tabucol cannot reach one before a reasonable time limit.

C. Principle of the Tabucol Algorithm

As an instance of tabu search, Tabucol is a local search algorithm that attempts to iteratively improve the current solution by selecting (generally) the best candidate in its neighborhood, i.e. the set obtained by slightly modifying the current solution.

For Tabucol, the neighborhood structure is based on critical 1-moves: two list-colorings are neighbors if they have exactly the same FL allocation except for one flight. A critical 1-move then consists in changing the FL of only one flight involved in a conflicting vertex.

The algorithm is initialized by taking the coloring corresponding to the RFL, i.e. FL i = RFL i , as the current solution. Then, at each iteration of Tabucol, the best critical 1-move which is not tabu and minimizes the objective function is selected. Next, the reverse move of the chosen critical 1-move is added to the tabu list for a given number of iterations, called the tabu tenure. The tabu tenure is managed dynamically to prevent all moves from being forbidden as the number of critical neighbors may vary at each iteration.

We take the same parameters as [START_REF] Galinier | A survey of local search methods for graph coloring[END_REF] for the tabu tenure and use also the same aspiration criteria, i.e. a forbidden move in the tabu list is nevertheless selected if it improves on the currently-known best solution. The algorithm stops when the first objective is equal to zero (no remaining conflict) or after a given number of iterations without improvement of the best solution encountered so far.

V. RESULTS

This model has been implemented and experimented on three data sets of ECAC-wide traffic. The initial data set consists of one day of 2002 traffic, containing around 22 500 flights. The other data sets are build upon this day of traffic by making it uniformly 25 % and 50 % denser respectively, leading to instances with up to 32 000 flights and 700 000 constraints. Those three instances will be referred as ECAC-0, ECAC-25 and ECAC-50 in the following. Table I sums up a few dimensions of the instances, as the number of flights and conflicts. The number of constraints corresponds to the number of edges in the conflict graph ("potential" conflicts) and the number of conflicts is obtained by allocating each flight to its RFL. The resulting conflict graphs have a density on the order of 1 ‰. As Tabucol is a stochastic algorithm, it is run 10 times for each instance, with a fixed number of one million iterations without improvement, in order to show the stability of solutions. The results displayed in the following are, for each instance, the best solution encountered (in terms of solved conflicts) and the average over the 10 different runs.

A. FL Allocation

All instances where no uncertainty is considered (ε max = 0) are solved entirely for ∆ max = 30, with 73 % to 80 % flights obtaining their RFL and only 7 % to 12 % flights being deviated of more than one FL.

Table II gives trade-offs between the number of remaining conflicts and the value of the ∆ max parameter. As expected, the number of remaining conflicts increases when flexibility is reduced (i.e. instances cannot be solved with very low discrepancy from RFL). However, reducing ∆ max generally leads to lower values of cost F L as unsolved conflicts do not contribute to this criterion. The distributions of discrepancies to RFL for different values of the maximum deviation are shown on figure 3. We observe that the variation of ∆ max has little influence on the number of flights obtaining their RFL.

B. Robustness to Temporal Uncertainty

The results presented in the previous section consider no uncertainty on the trajectories. We tested our solutions by introducing random uncertainties (uniformly chosen in the interval [-ε c , ε c ]) on takeoff times in the CATS simulator and measuring the number of remaining conflicts.

Figure 4 shows the percentage of remaining en-route conflicts observed during the CATS simulations for different values of ε max and ε c .

It could be expected to observe no remaining conflicts for instances where ε c ≤ ε max (at least for fully solved instances, table III). However, changing flight levels leads to a slight modification of trajectories (cruise phase is shortened when FL is increased and the ground speed modified accordingly), which impacts the initial conflict detection. For higher values of the ε c (ε c ≥ 2), handling uncertainties as presented in section III-B gives good results as the number of remaining conflicts is reduced when ε max increases.

Nevertheless, increasing ε max leads to harder instances, which has a direct influence on the ability to solve all conflicts and on the cost of the solutions. constraint graph for handling a ±3 min uncertainty is two to three times denser than with no uncertainty, which explains the observed results: cost F L is twice higher and up to 9 % of the constraints cannot be satisfied. The distributions of discrepancies to RFL associated to these solutions are depicted on figure 5. For the biggest instances, there is still more than 50 % of the flights obtaining their RFL, but almost 30 % are deviated of two or three flight levels.

Computing times are highly correlated to the instance size and density, and vary from a few seconds to 50 min on a standard workstation (Opteron 6128 at 800 MHz).

VI. CONCLUSION AND FURTHER WORK

We have presented a novel flight level allocation scheme to avoid the conflicts occurring during the cruise phase of intersecting flights, which can amount to 20-30 % of all conflicts in a day of traffic in the ECAC airspace, and thus enhance the airspace capacity to accommodate the foreseen traffic growth in Europe. Our model is based on an exact temporal characterization of each potential conflict computed during a pre-processing phase and, unlike all other approaches known to the authors, allocates a level to each flight individually. This scheme yields much better results than flow-based aggregated ones as reported in [START_REF] Letrouit | Optimisation du réseau de routes aériennes en Europe[END_REF], [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF], and is therefore able to take into account, at a reasonable cost, a given amount of time uncertainty on the following of the trajectories. Furthermore, the discrepancy between the allocated flight level and the one requested in the flight plan is tightly bounded in our model to keep the operational cost of such regulation as low as possible for airlines.

The application of our scheme at the continental scale in Europe, with a traffic comprising more than 30 000 flights daily, leads to huge combinatorial optimization problems, with more than 700 000 conflict constraints for the largest instance. Nevertheless, we were able to solve all instances with at most three levels of discrepancy, from one day of traffic of 2002 up to an expanded instance with 50 % more flights, roughly corresponding to current forecast for 2025. Among the various optimization techniques tried, the best results in terms of number of shifted aircraft, overall cost and computation times were consistently obtained by a tailored version of the Tabucol graph coloring algorithm.

However, our approach does not yet take real-time updates into account, like the cancellation or postponement of a flight, but the pre-processing of conflicts and the low computation time of our resolution algorithm would allow our allocation scheme to be integrated within a dynamic framework like the sliding window algorithm presented in [START_REF] Allignol | Planification de trajectoires pour l'optimisation du trafic aérien[END_REF]. Another shortcoming of our model is the linear contribution of FL discrepancy in the objective function, which does not exactly correspond to the airlines operational cost in terms of fuel consumption or CO 2 emissions. Our algorithm could though very easily take into account more sensible cost functions, provided we are given enough configuration and performance data by airlines and aircraft manufacturers. Eventually, our approach could be combined with other regulation schemes, as proposed in [START_REF] Barnier | Combining flight level allocation with ground holding to optimize 4D-deconfliction[END_REF], or real-time conflict solvers, to further reduce the conflict resolution workload of ATC.

Fig. 1 . 3 )

 13 Fig. 1. Vertical and horizontal separation. Another aircraft cannot be inside the cylinder at the same time.

Fig. 2 .

 2 Fig. 2. Discrepancy to RFL for allocation of FL to flows (left) and individual flights FL allocation (right). Flights that obtained their RFL are depicted in light green, flights with a 10, 20, 30 discrepancy are depicted respectively in green, orange and red.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Deviation from RFL vs ∆max.

Fig. 5 .

 5 Fig. 5. Deviation from RFL vs εmax.

  Table III sums up the number of non-satisfied constraints and the operational cost for all instances processed. The

TABLE III IMPACT

 III OF HANDLING A εmax UNCERTAINTY ON REMAINING CONFLICTS AND COSTF L (∆max = 30).

	Instance	εmax	Remaining conflicts	costF L
			best (aver.)		best	(aver.)
		0	0	(0)	0 %	6637	(6656)
	ECAC-0	1 2	0 0	(0) (1.6)	0 % 10976 (11089) 0 % 17344 (17479)
		3	96	(110)	0.5 % 22915 (22620)
		0	0	(0)	0 %	9674	(9724)
	ECAC-25	1 2	4 80	(10) 0.02 % 17720 (17610) (96) 0.3 % 26907 (26758)
		3	790	(822)	2.4 % 31034 (30931)
		0	0	(0)	0 % 13553 (13634)
	ECAC-50	1 2	16 527	(28) 0.06 % 26481 (26359) (573) 1.5 % 35783 (35071)
		3	2184 (2313)	4.8 % 40321 (39189)

Only the static semi-circular rule that avoids face-to-face conflicts impacts flight level allocation.

The Complete Air Traffic Simulator developed at DSNA/DTI.


Michele A Quinto

Dominique Houzet
email: dominique.houzet@gipsa-lab.grenoble-inp.fr

Fanny Buyens
email: fanny.buyens@cea.fr

Tetrahedral Volume Reconstruction in X-Ray Tomography using GPU architecture

In this paper, we propose the use of the graphics processor unit (GPU) to accelerate a ray-tracing method in the framework of X-ray tomographic image reconstruction. We first describe an innovative iterative reconstruction method we have developed based on a tetrahedral volume with conjugate gradient. We do not use voxels here but instead tetrahedrons to increase the quality of reconstruction and the reduction of data as thus we need less resolution of the volume to fit the object reconstructed. This is an important point to use the GPU. We present here the algorithms adapted to the GPU and the results obtained compared to CPU.

I. INTRODUCTION

Tomography reconstruction from projections data is an inverse problem widely used in the medical imaging field. With sufficiently large number of projections over the required angle, the FBP (Filtered Back-Projection) algorithms allow fast and accurate reconstructions. However in the cases of limited views (low dose imaging) and/or limited angle (specific constrains of the setup), the data available for inversion are not complete, the problem becomes more ill-conditioned and the results show significant artifacts. In these situations, an alternative approach of reconstruction, based on a discrete model of the problem, consists in using an iterative algorithm or a statistical modelisation of the problem to compute an estimate of the unknown object. These methods are classicaly based on a volume discretization into a set of voxels and provide 3D maps of densities. High computation time and memory storage are their main disadvantages. Moreover, whatever the application is, the volumes are segmented for a quantitative analysis. Numerous methods of segmentation with different interpretations of the contours and various minimized energy functions are offered with results that can depend on their use by the application. This work presents a novel approach of tomographic reconstruction simultaneously to segmentation of the different materials of the object to reconstruct. The process of reconstruction is no more based on a regular grid of voxels but on a mesh composed of non regular tetraedra. After an initialization step, the method runs into three main steps: reconstruction, segmentation and adaptation of the mesh, that iteratively alternate until convergence. For that purpose, we have adapted and optimized two iterative algorithms of reconstruction (OSEM and Conjgate Gradient) usually used in a conventionnal way to be performed on irregular grids of triangular or tetraedric elements. For the segmentation step, a geometrical approach (level set) was implemented. The adaptation of the mesh to the content of the estimated image is then constrained by the contours segmented the step before that makes it progressively coarse from the edges of the object to the limits of the domain of reconstruction. At the end of the process, the result is a classical tomographic image in gray levels, whose representation by an adaptive mesh to its content provides a corresponding segmentation. The results show that the method provides reliable reconstructions and leads to drastically decrease the memory storage. In this context, both the 3D operators of projection and backprojection were implemented on parallel architectures namely the graphic processor units (GPU). These operators are based on the raytracing method used in image rendering.

II. X-RAY TOMOGRAPHY

X-rays tomography is a non invasive technique which provides 3D reconstruction of an object. Its principle is based on the interaction of the X-ray beams with the materials composing the object, and the counting of the transmitted photons by detectors placed in front of the X-ray tube. The acquired data are collected according to multiple orientations (angles) of the source-detector system over the object. The number of projections and the rotation step vary according to the application and the type of devices. By means of these data, an image in grey levels is mathematically reconstructed representing the internal density of the object.

Commonly, the algebraic algorithms used for reconstruction are based on a representation of the volume by a regular lattice of voxels. However, the spatial resolution of the detectors increases and consequently the number of voxels too to preserve the high spatial resolution. The number of unknowns to estimate can thus be very important, and the algebraic methods of reconstruction very time consuming.

The algebraic methods model the inverse problem by a discrete approach: the projections are interpreted as samples of a function with discrete variables, and the object is represented as an unknown function f defined as a linear combination of basic functions (commonly characterizing voxels). So the inverse problem is written as a linear system:

p = H f (1)
where the vector p contains the measures of the projections. The size of the vector p is equal to N D × N p , where N D is the number of the pixel detector and N p the number of projections. H is the matrix system ; it depends on the geometry of acquisition. The value of its elements h ij represents the contribution of element j to X-ray beam projection i. Commonly, h ij is equal to 1 if pixel j is intersected by X-ray beam i, otherwise 0, but other expressions of contribution can be used as, for example, the length or the surface of the X-ray beam intercepted in the pixel. Matrix H is very sparse : most of its elements are zero. Nevertheless, the computation of matrix H is very time consuming. In the case of small data sets, H can be pre-calculated and stored in memory. Although, as the size of the data sets enlarge, matrix H cannot be stored in memory and several approaches have been proposed to make it more compact. For example, [START_REF] Herraiz | FIRST: Fast iterative reconstruction software for (PET) tomography[END_REF] proposed a list-mode approach to only store the non-zero elements of the matrix. A more efficient solution is to calculate the coefficients of matrix H on the fly using ray-tracing methods. The optimization of matrix H computation improve the performances of both the projection and the back projection which are the main operators of algebraic algorithms of reconstruction. The approaches used to compute projection and backprojection are classified as [START_REF] Zhuang | Numerical evaluation of methods for computing tomographic projections[END_REF]:

• ray-driven: the X-ray beams are centered on the pixels of the detector

• pixel-driven: the X-ray beams are centered on the pixels of the volume of reconstruction.

In this work, the two approaches were implemented and optimized on GPU (Graphic Processor Unit, multi-threaded parallel architectures), to study the parallelisation and optimisation characteristics of both the projector and the back projector (the latest being the transposed operator of the first one). In the case of the ray-driven implementation, the projector performs gather operations by reading the data and writing the result to a specific memory address (corresponding to the thread IDs), while the backprojector performs scatter operations by updating the values of the estimated image shared by several threads. On the contrary, in the case of the pixel-driven implementation, the projector performs scatter operations by ... while the backprojector performs gather operations by ... . For more efficiency, gather operations are prefered to ensure no memory writing conflicts. In this way, we implemented a ray-driven approach for the two operators [START_REF] Siddon | Fast calculation of the exact radiological path for a threedimensional CT array[END_REF], the projection requiring more operations than the backprojection. The small number of detectors induced a limited amount of exploitable parallelism but a sufficient one for a single GPU.

III. PROPOSED METHOD

In this work, we proposed a method of tomography reconstruction combined with a stage of segmentation, that results to a classical image of reconstruction in grey levels paired with an image of segmentation of the studied object. In this method, we adopted a mesh representation based on tetrahedral elements to discretize the volume of reconstruction in order to obtain a

Image & Mesh initialization

Reconstruction Segmentation Mesh generation better representation of convex objets and to drastically reduce the number of elements during the process of reconstruction (the mesh being adapted to the content of the estimated image). Thus, as we do not consider anymore a three-dimensional volume described by voxels but a set of irregular tetrahedra, common used algebraic algorithms of reconstruction as EM and Conjugate Gradient were implemented and optimized to be used on irregular grids of tetraedra.

J (f ) <
The proposed method, we refer as ATM (Adaptive Triangular Mesh reconstruction), is composed of three main stages after initialization : reconstruction, segmentation and adaptation of the mesh which successively alternate until convergence (see Figure 1). At the end of the process, the estimated volume is represented by an optimized number of elements in the mesh, what drastically reduces the size of the reconstructed data and consequently turns down the problems of memory storage.

Among the three stages, only the reconstruction one is here concerned with a GPU implementation (through both the projection and the backprojection operators) because of its high time consuming and its large memory load. By the way, classical algorithms with a regular parallelism (independent nested loops) can be used for the two other stages.

In this work, the h i,j elements of matrix H represent the intercepted length of the X-ray beam i with the mesh element j (see Figure 2). Contrary to the conventional case of a fixed grid (common voxel discretization), the h i,j elements can not be stored throughout the process of reconstruction because each time the mesh is modified the elements of the matrix H have to be computed.

IV. IMPLEMENTATION OF THE PROJECTION AND THE BACKPROJECTION OPERATORS ON GPU

Both the projection (ALGORITHM 1) and the backprojection (ALGORITHM 2) operators were implemented on GPU using a ray-tracing method. Three methods designed for 3D rendering in computer graphics and adapted for our application were studied: one classification method (REFERENCE !!) and two parametric ones (REFERENCES !!). The first results have shown that one of the two parametric methods do not converge using single-precision floating point values because of the precision required to compute the intersection ray-tetrahedron. Also the classification method is slowest one. So, in the case of single-precision floating point, the second parametric method shows the best performances. The choice of single-precision processing is important as the ratio of performances can be very different according the used GPU. For example, there is a factor equal to 3 for the last Kepler GPU generation. The second parametric method thus used (REFERENCE!!) is based on the method of comparison of the parameters with the storage of the plane equations of the faces of the tetrahedra (see Figure 3): The outgoing face is identified by t = min{t i } and the intercepted length is computed as = t 1 -t 0 .

The used 3D data structures are the following ones:

• i th vertex: vertices[3 • i, • • • , 3 • i + 2]
• vertex indexes of i th tetrahedron: tetrahedron corners[4

• i, • • • , 4 • i + 3]
• neighbor opposite to corner of overall index j: tetrahedron neighbors[j]

• uniform density of i th tetrahedron: tetrahedron densities[i]

Algorithm 1 Projection operator end while 12: end for This method requires a initialization stage of the rays to identify the input tetrahedron into the tetrahedral volume. This stage is very expensive in terms of time consuming when compared to a CPU execution time of one projection. In order to speed up this stage of initialization, a quad-tree hierarchical structure was implemented on each face of the volume. A regular tetrahedral volume is then considered, and a leaf of the quadtree only contains the index of two tetrahedrons (see Figure 4).

To optimize the performance of the quad-tree hierarchical structure, the following subdivisions of the faces of the tetrahedral volume were used:

• 16 x 16 (512 peripheral tetrahedrons)

• 32 x 32 (2048 peripheral tetrahedrons) • 64 x 64 (8192 peripheral tetrahedrons)

We observe that the execution time only depends on the number of rays to initialize (see Figure 5). The stage of initialization of the rays with a hierarchical structure is 700× faster than without (see Figure 6). This optimized initialization step is thus now less time consuming than the projection itself.

V. RESULTS

Both the projection and backprojection algorithms were executed using 2D thread blocks. The size of the thread blocks was chosen in order to maximize the number of threads with respect to the material resources (registers) of the graphic card.

Each thread reads the first tetrahedron vector by its ID: if the value is equal to -1 the execution stops, otherwise the operator is performed. The latter is achieved when the next index of the tetrahedron is equal to -1.

A. Validation of the 3D projection operator

The validation of the projector was performed considering the object described in Table II (see Figure 7). A regular object was chosen to be perfectly described by the mesh so that it has no influence on the results. The results are shown in Figure 8. The 3D tetrahedral mesh is generated by using the library Tetgen [START_REF] Si | Tetgen: A quality tetrahedral mesh generator and a 3D Delaunay triangulator[END_REF]. The configuration of acquisition is described in Table III.

We studied the error introduced by a mesh which not perfectly describes an object. The test object is here an ellipsoid (see Figure 9) whose surface is discretized by a increasing number of points from 867 to 1622456. The error, expressed in percentage, is defined as:

Err = p T -p A p A • 100% (2)
where, p T is the vector of the tetrahedral projections and p A is the vector of the analytical projections. Figure 10 shows the error, Err, according to the number of points describing the surface of the ellipsoid. We observe that the error decreases when the number of points describing ellipsoid increases, up to 357612 points. After what, the error remains constant at 1%.

B. Performance CPU/GPU

The performance were estimated by calculating the projection of the numerical phantom previously described (see Figure 7). The hardware used for these experiments are: • CPU: Intel Xenon E5440 @ 2.83GHz

• GPU: Nvidia-Fermi Tesla C2070

To characterize the performance obtained on GPU, comparisons were leaden by considering the execution time according to the detector size and the number of tetrahedra (see Figure 11). The memory transfert was taken into account to compute the speed-up factor. As result, the speed-up factor is in the range of 40× to 80×. It depends on the size of the detector and the number of tetrahedra discretizing the volume. Large detectors are needed to take benefit from a high degree of parallelism. However, the sequential process of the tetrahedra limits the benefit. We studied the influence of the memory transfers on the execution time of the projector implemented on GPU. Figure 12 shows that bigger is the detector, more time consuming are the memory transfers. The memory transfer time reaches the same value as the time of execution of the projection for a detector size of 8192×8182 pixels. This is a main drawback that will impact speedup for larger detectors.

VI. CONCLUSION

In this paper we explored the feasibility of using Graphic Processor Units (GPUs) for tomographic reconstruction using tetrahedral discretization of the volume. We demonstrated fast projection/backprojection operators using a ray-tracing algorithm with low required memory size. We also demonstrated a speedup factor between 40× and 80× in the case of an implementation of the operators on GPU over the CPU implementation with a low impact on the precision (error less than 1%).

Compared with a classical voxel-based method, the ATM method implemented on GPU suffers from a loss of memory locality and element regularity. Indeed, the voxel-based approaches discretize the volume into a set of regular elements, and GPU implementations can take benefit of that by using for example the texture memory to get data locality that results in divergence free and bank conflict free memory operations. Besides, all the ray beams across the same number of voxels that balances the work load among all the threads in the warp [START_REF] Chou | A fast forward projection using multithreads for multirays on GPUs in medical image reconstruction[END_REF].

In our sutdy, an irregular discretization of the volume was chosen in order to fit the objet shape with a tetrahedral mesh. Because of the irregularity of the elements, two main disadavantages limit a highyield informatic optimization : the locality of the memory is lost (or almost), and the number of the tetrahedron processed by each thread belonging to the same warp can be different. That induces an unbalanced work load over a warp which means that threads still run while others remain idle. The minimization of the unbalanced work load is an important perspective to this work.

Figure 1 :

 1 Figure 1: Diagram of the proposed method.

ray projection i th linear detector Figure 2 : 1 nFigure 3 :

 detector213 Figure 2: Example of the 2D projection operator.

•

  compute the normal vector n of the face, for example AB × BC • compute s 1 = BD • n and s 2 = d • n if s 1 and s 2 have opposed signs, the ray r goes out of this face, and thus parameter t = s 0 B • n d • n

Figure 4 :Figure 5 :

 45 Figure 4: Decomposition of the surface of the tetrahedral volume.

Figure 6 :

 6 Figure 6: Gain between the stage of initialization without and with the hierarchical data structure.

Figure 7 :

 7 Figure 7: Object described to validate the 3D projector

Figure 8 :

 8 Figure 8: Projections of the objet at θ = 0 using (a) an analytical projector (intersection equations of the rays and the cubic elements) and a tetraedral projector (b) implemented on CPU and (c) implemented on GPU. The maximun error is about 10 -15 (d).

Figure 9 :

 9 Figure 9: Ellipsoid surface described by 14011 points.

  Error between analytical projections and tetrahedral projections Number of points of the ellipoid

Figure 10 :

 10 Figure 10: Error between the analytical projections and the tetrahedral ones with respect to the density of the tetrahedral mesh.

Figure 11 :Figure 12 :

 1112 Figure 11: Speed-up of the GPU projector with respect to the number of the tetrahedral elements and considering several sizes of detector.

Table I :

 I Memory write operations.

		ray-driven	voxel-driven
	Projector	Gather	Scatter
	Backprojector	Scatter	Gather

Table II :

 II Description of the test object used for the validation of the projector

	Object	xmin, xmax	ymin, ymax	zmin, zmax
	Externe cube	-0.25, 0.25	-0.25, 0.25	-0.25, 0.25
	1st hole	-0.1875, -0.03125	-0.1875, -0.03125	-0.078125, 0.078125
	2nd hole	0.046875, 0.085938	0.046875, 0.085938	-0.019531, 0.019531
	3th hole	0.082031, 0.09375	0.144531, 0.15625	-0.117188, 0.117188
	4th hole	0.160153, 0.167969	0.128906, 0.136719	-0.003906, 0.003906

Table III :

 III Configuration used for the validation of the projector on GPU.

	Parameter	Value
	source-object distance	98 mm
	object-detector distance	132 mm
	size of the detector	14.2mm × 14.2 mm
	resolution of the detector (pixels)	256 × 256
	number of the projections	64

ACKNOWLEDGMENT

The authors would like to thank Tomacz Toczek for design the ray-tracing algorithm on the CPU/GPU in the framework of the 3D rendering in computer graphic.


François Courvoisier

Amaury Mathis

Jinggui Zhang

Luc Froehly

Luca Furfaro

Maxime Jacquot

Remo Giust

Pierre-Ambroise Lacourt

John Michael Dudley

Femtosecond laser material processing with nondiffracting light

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Although versatile and widely used, femtosecond laser micro-nanomachining faces a challenge for the fabrication of high aspect ratio or deep structures. The control of the profile along the longitudinal dimension is extremely difficult. In this context, our approach is based on controlling the direction of light rather than shaping the intensity pattern in one plane. Here we review our recent work in this field using Bessel beams and accelerating beams.

In Bessel beams, light propagates along the generatrix of a cone, with a fixed radial wavevector. Contrary to Gaussian beams, the femtosecond filamentation regime of Bessel beams can be stationary, and can generate extended plasma tracks in dielectrics. This allowed us to generate nanochannels in glass with aspect ratio up to 100:1 [START_REF] Bhuyan | High aspect ratio nanochannel machining using single shot femtosecond Bessel beams[END_REF]. Numerical simulations of the nonlinear propagation in this regime show the importance of the conical structure in generating extremely dense plasmas [START_REF] Courvoisier | Applications of femtosecond Bessel beams to laser ablation[END_REF].

Airy beams and more generally accelerating beams constitute another family of beams that possess quasi-nondiffracting behaviour and their nonlinear propagation can also be stationary. Their primary intensity lobe propagates along a curved trajectory that can be arbitrarily shaped even in the nonparaxial regime. Since this lobe is adjacent to a region where no light propagates, accelerating beams can be used for direct curved edge profiling and trench processing in transparent and opaque materials [START_REF] Mathis | Micromachining along a curve: Femtosecond laser micromachining of curved profiles in diamond and silicon using accelerating beams[END_REF][START_REF] Mathis | Direct machining of curved trenches in silicon with f emtosecond accelerating beams[END_REF].

Our results show that controlling the linear and nonlinear propagation of ultrashort femtosecond laser pulses by nondiffracting beams is a key new technological approach for laser processing.

  (a) Intensity distribution of a circularly accelerating beam. (b) Edge profiling on silicon and diamond.

  Femtosecond Laser Material Processing with nondiffracting light F. Courvoisier*, A. Mathis, J. Zhang, L. Froehly, L. Furfaro, M. Jacquot, R. Giust, P.A. Lacourt, and J.M. Dudley Département d'Optique P.M. Duffieux, Institut FEMTO-ST, UMR 6174 CNRS Université de Franche-Comté, 25030 Besancon Cedex, France francois.courvoisier@femto-st.fr