+ S Abba Rapaya

P Notry

D Delahaye

[EN-A-37] Coordinated Sequencing of Traffic on Multiple En-route Constraint Points

Keywords: Air Traffic Management, Traffic sequencing, Arrival Management, Conflict resolution, Multi-objective optimization, Simulated annealing, Sliding window

Air transportation traffic is progressively increasing over the years and dealing with it is an essential task to guarantee fluid flights in the future. Several works already indexed multiple aspects of aviation, among them, the E-MAN system. It introduced the sequencing of arriving traffic, starting from early stages of the En-route phase. This change facilitated the work for the approach controllers but also increased the workload of the En-route controllers. To handle that workload, controllers are now assisted by tools that consider the new constraints introduced by the arrival management system and propose advisories. From that same perspective, our project focuses on an algorithm for a helper tool that will combine both aspects of traffic sequencing in the En-route phase and conflict resolution. With this novel approach, we automatically generate near-to-optimal flight decisions, given that we can modify the speed and the flight level to respect the sequencing constraints and cut down potential conflicts. We categorize the problem as a mathematical optimization case. Thus, we describe a detailed mathematical model which covers all the aspects of the problem. This model gives a basis for the implementation of the flight optimizer. Later, we propose a solution based on a sliding window simulated annealing algorithm which reduces the complexity and takes into account uncertainties. Finally, we successfully test an implementation of the solution with real-life traffic data. It corresponds to flights within France going towards Paris CDG airport over a period of 24 hours. The results demonstrate a total conflicts resolution with satisfying compliance with sequencing constraints.

Introduction

The air traffic has been increasing lately. In its Global Market Forecast (GMF) [1], Airbus has predicted that the global air traffic will encounter an annual growth of 4.4% for the next coming 20 years. With that increase of traffic, some major airports will now rely on arrival management systems such as the XMAN system to help them manage their flow of traffic and generate efficient flights sequencing. However, most of the workload is now transferred to the En-route controllers since they will have to deal with occurring conflicts and also try to respect new constraints imposed by the XMAN to flights.In En-route, the airspace is characterized by several route flows crossing into potential conflict points. To avoid those conflicts, some separation constraints are defined for the flights and controllers need to make sure they are respected by using classic methods like vectoring, speed controlling or flight level changes. Each decision applied to an aircraft can have an impact on surrounding flights.This research aims to develop an algorithm that will assist the controllers. Based on the traffic flows analysis, this algorithm will provide near optimum flight decisions in order to remove conflicts at the crossing points.From the controller's side, the main objective is to reduce the workload and minimize the frequency occupation. This will be done by automatically eliminating conflicts due to lack of separation at crossing points and route links. It also implies reducing the number of flight level changes and simplifying the instructions given to pilots.Another objective is to integrate our tool with an XMAN system by generating decisions compatible with XMAN constraints.Finally, from the airlines point of view, we will need to take into consideration the preferred flight profile by generating decisions not far from that profile.The first part of the paper introduces some previous related works linked to arrival manager optimization in a wide sense. The second part presents the associated mathematical model. The meta-heuristic (simulated annealing) used for solving the underlying optimization problem is presented in the third part. Finally, the fourth part presents the results given by the algorithm on a real case at Paris CDG. This section presents the concepts and previous works related to our problem. It focuses on points around air traffic sequencing, conflict resolution and mathematical optimization methods. It is concluded by a synthesis which serves as baselines for our modeling and solution approach.

Extended AMAN

The Cross-border SESAR Trials for Enhanced Arrival Management (X-STREAM) project main objective is to extend the current horizon of the Arrival Management System beyond 200 NM towards the upstream Area Control Center (ACC). Extended-AMAN (E-AMAN) [START_REF] Sesar | Extended arrival management (aman) horizon[END_REF] allows for the sequencing of arrival traffic much earlier than is currently the case, by extending the AMAN horizon from the airspace close to the airport to further upstream and so allowing more smooth traffic management. Controllers in the upstream sectors, which may be in a different control center or even a different Functional Airspace Block (FAB), obtain system advisories to support an earlier pre-sequencing of aircraft. Controllers implement those advisories by, for example, instructing pilots to adjust the aircraft speed along the descent or even before top-of-descent, thus reducing the need for holding and decreasing fuel consumption.

Characteristics

The traffic begins to be processed at the Eligibility Horizon (EH) and controllers are provided with AMAN advisories from the Active Advisory Horizon (AAH) onward [START_REF] Sesar | Integrated sequence building/optimization of queues[END_REF]. With the extension of the horizon of the AMAN to the En-route phase, as illustrated on Figure 1, the EH goes further to around 180-200NM which supports the controllers in applying more efficient arrival management techniques. The typical optimum top of descent is approximately 100-120NM from touchdown which implies that up to 100NM of flight within the extended AAH will be in En-route airspace.

Decisions on flights

AMAN advice may be in the form of a Target Time at the Initial Metering Point (IMP) or Time To Lose (TTL)/Time To Gain (TTG) advisories to controllers calculated by AMAN working back from the runway time [START_REF] Sesar | Integrated sequence building/optimization of queues[END_REF]. Upstream ACC controllers can then provide instructions to pilots to make adjustment of speed to meet TTL/TTG need, by executing advisories for speed coming from tools like Speed And Route Advisory (SARA) or to delegate responsibility for adherence to a Controlled Time of Arrival (CTA) to flight crew.

Previous works

Optimization problems involve the minimization (or maximization) of a function by choosing input values from a given set of data and computing the value of the function. Nowadays, many studies have been conducted on the resolution of air traffic conflicts and the sequencing using various mathematical optimization techniques. This section will present a few of those works. [START_REF] Vela | A mixed integer program for flight-level assignment and speed control for conflict resolution[END_REF] presented a new approach based on concepts of speed control and flight-level assignment for conflict resolution over predefined routes. It was based on a Mixed Integer Linear Program (MILP). The way it was formulated helped to reduce fuel burn over time horizons between 15-45 minutes. [START_REF] Li | A flowbased flight scheduler for en-route air traffic management[END_REF] worked on a flight scheduler for En-route air traffic. It was based on the application of delays or route changes to flights to produce an ordered merged sequence of flights at the exit nodes. They didn't tackle the conflict resolution aspect in that paper. During the 2017 SESAR Innovation Days, [START_REF] Courchelle | Simulated annealing for strategic traffic deconfliction by subliminal speed control under wind uncertainties[END_REF] presented a method to perform small changes on aircraft speed in order to resolve conflicts. In their model, they considered that the True Air Speed (TAS) was constant for each aircraft in the airspace and the uncertainties appeared due to wind components, with wind data collected from MétéoFrance PEARP (Prévision d'Ensemble ARPège). A protection area, as illustrated in Figure 2, is defined around each aircraft to represent its potential positions. Thus, two overlapping protection areas implies a potential conflict between the two aircraft that must be dealt with. The solution was implemented in Python using a simulated annealing algorithm. As a result, the algorithm was able to solve at least 70% of virtual conflicts with a computing time of 30 minutes to 2 hours for a 4 hours traffic data depending on the refining of the parameters. On the other hand, [START_REF] Barragan | Optimised algorithms for extended arrival managers[END_REF] presented their work at the 29th Congress of the International Council of the Aeronautical Sciences (ICAS). It was done on a study for a collaborative tool that aimed at helping controllers with the integration of the Extended Arrival Manager. They studied the aircraft sequencing problem as an optimization problem under different resolution methods (Linear programming, non-linear programming, heuristic methods). The objective function took into consideration the runway capacity, the delays and the fuel consumption. They also didn't consider the conflict resolution aspect.

Synthesis

As seen in this section, the Extended Arrival Management tool gives advisory directives to controllers to help create an efficient sequence of flights starting from a distant En-route point. However, it is still up to the controller to decide how he will implement those directives taking into account the potential occurring conflicts. A tool like SARA, in combination with XMAN, can helps controllers by giving speed directives but no global conflict resolution is provided. Thus, it will be an added value for our solution to provide both sequencing and conflict resolution capabilities. Moreover, most of the previous works done on the optimization of the traffic focused on changes made on aircraft speed or route to delay or advance a flight. Therefore, not many alternatives were given for conflict resolution when there is a dense traffic. In that case, we can introduce changes on flight levels to broaden the solutions possibilities as suggested by [START_REF] Courchelle | Simulated annealing for strategic traffic deconfliction by subliminal speed control under wind uncertainties[END_REF]. As for the optimization method to apply, the nature of our problem makes it more appropriate for a meta-heuristic method, given the dimension of the problem. Also, the simulated annealing appears as a good choice regarding its stability, its memory consumption profile and its speed. [START_REF] Ji | Aircraft merging and sequencing problems in tma[END_REF] introduced an implementation of the simulated annealing algorithm coupled with the sliding window concept, the base of our implementation will be that algorithm. Another added value we can integrate in our solution is the collaboration with multiple XMAN tools at the same time. The state of the art review enabled us to define basis for our study case. In the next section, we will define a detailed mathematical model for the problem. This will highlight abstract concepts useful for the design of our solution.

Mathematical model

This section focuses on an abstract mathematical modeling of our problem. Each aspect of the problem is represented using mathematical notations that can be interpreted into any optimization algorithm later. It presents the global network representation, how uncertainties are modeled, equations related to conflict detection, the optimization objective function definition and an evaluation of the complexity of the problem.

Constraints

We have identified some constraints for our model. They are elements that will influence the way the tool will handle the optimization process.

Sequencing constraints : those constraints are coming from the sequencing system. The XMAN generates a TTL or TTG for flights. It must be considered during the optimization process.

Aircraft performances : those correspond to the physical performances of the aircraft. It consists of the minimum (maximum) speed and the minimum (maximum) flight level the aircraft can fly. Basically, our tool should not advise a flight to perform outside of its physical limits.

Flight preferred profile : It is a constraint that can be proposed by the airline. In our case, we will consider the preferred cruising flight level. It is a flight level which if cleared to the flight, will help the airline to better achieve their business objective.

Separation constraints : we have three types of separation : the distance separation which is a minimum defined distance should be maintained between two aircraft when crossing a specific point, the time separation for which two consecutive aircraft should reach a given location with a minimum defined time interval and the wake turbulence separation which is the separation distance between two consecutive aircraft flying on the same link at the same flight level and which depends on their wake turbulence categorization as defined by ICAO.

Network modeling

The network of routes is modeled as an oriented graph with nodes and oriented links.

Nodes A node represents any point of interest on a route. Most of the time, it will be where two routes are crossing. We can define a separation constraint on them (distance or time separation). A node n is characterized by its cartesian coordinates (x, y), x, y ∈ R, a minimum separation distance d sep , in case of distance separation constraint and a minimum separation time t sep , in case of time separation constraint.The set of nodes will be noted N : N = {n i / i ∈ N, i ≤ nb nodes } where nb nodes is the total number of nodes in the network Links A link is a portion of route between two nodes. It is oriented from an origin node to a destination node. Each link is characterized by an entry node n ori , an exit node n dest and a length d l = dist(n ori , n dest ). The set of links will be noted

L : L = {l i / i ∈ N, i ≤ nb links }. ∀l i ∈ L, l i = (n ori , n dest ) / n ori
,n dest ∈ N, n ori n dest where nb links is the total number of links in the network.

Routes A route is an ordered list of adjacent links and is characterized by a parity p which can be odd or even: p ∈ {ODD, EV EN} and a list of available flight levels f ls = { f l / f l ∈ N}. We will use the semi-circular rule to determine those flight levels depending on the parity of the route. The set of routes will be noted R : R = {r i / i ∈ N, i ≤ nb routes }, ∀r ∈ R, r = {l i / l i ∈ L} where nb routes is the total number of routes in the network.

Flights modeling

A flight f is characterized by a speed V f , a flight level FL f , the time of arrival at the first node (entry in the airspace) t f init , a wake turbulence category wtCat f and a defined route r f ∈ R.Apart from those characteristics, a flight is also subject to some constraints: a maximum speed V max f , a maximum flight level FL max f , a preferred flight level FL pre f f given by the airline and a time constraint ttl f given by the sequencing system (a positive value means a time to gain and a negative value refers to a time to lose both before reaching the final exit node). The set of flights will be noted F .

F = { f i / i ∈ N, i ≤ nb f lights }
where nb f lights is the total number of flights in the network.

Decision variables

During the optimization process, two features of the flights are modified: the speed and the flight level. For that, we are using two decision variables which are the speed delta and the flight level delta.

Speed delta As expressed by [START_REF] Courchelle | Simulated annealing for strategic traffic deconfliction by subliminal speed control under wind uncertainties[END_REF] in their article, we will also use small speed changes to control flights speed. Each speed delta ∆V is an integer in the range of -6 to +3, expressed in percentage.

∀ f i ∈ F , ∆V i ∈ Z ∩ [-6, +3].The new TAS V i of flight f i can be ex- pressed as follow: V i = ∆V i 100 × V 0 i where V 0 i is the initial TAS.
Flight level delta The second decision variable we are using is the flight level delta. It represents the number of levels an aircraft should climb or descend. It is expressed as an integer varying from -2 to 2. A positive value represents a climb whereas a negative value represents a descend.

Uncertainties

In our assumptions, we supposed that the TAS of flights is constant all along its route. However, in the real world, each flight is subject to wind influence in the air which can affect the ground speed. This section presents how we model the wind and how it generates uncertainties in our model.

Wind modeling : we use a simplified model for the wind. The wind vector is characterized by a constant direction all over the airspace (it is represented as an angle α w relative to the geographical north) and an intensity W s randomly varying within a range of values, with W s ∈ [W s min , W s max ].

Flight uncertainty : the ground speed Gs of the aircraft is the resultant from both the TAS component and the Wind Speed component (

-→ Gs = ---→ T AS + --→ W s
). Let's assume a flight arriving at the entry node n ori of a link l at t ori , its time interval of arrival at the exit node n dest will be:

[t dest min , t dest max ],= [t ori + d l V fmax , t ori + d l V f min
]. This interval is built with the earliest time of arrival at the exit node t dest min which depends on V f max which is the maximum ground speed of the aircraft and the latest time of arrival at the exit node t dest max which depends on V f min which is the minimum ground speed of the aircraft.

Conflicts identification

We will consider two types of conflicts: node conflicts and link conflicts.

Node conflicts They are detected using the separation constraints. Depending on the type of separation constraint at the node (time or distance separation), the equations to detect node conflicts will differ. To detect a distance separation conflict, a circular protection area with a diameter of d sep is defined around the node. Each flight arriving at the node will be monitored using four different times : an early time of entry in the protection area t in min , a late time of entry in the protection area t in max , an early time of exit from the protection area t out min and a late time of exit from the protection area t out max . Two consecutive flights are considered to be in conflict at the node when they are both present within the protection area at the same time. In other words, the trailing flight enters the protection area before the leading flight exits it. Let's consider two flights f and g with f leading.

∀ f, g ∈ F , c f lt( f, g) =        1 if t f,out max -t g,in min < 0 0 otherwise
To detect a time separation conflict, a separation time t sep is defined on the node. Each flight arriving at the node will be monitored using two different times: an early time of arrival at the node t min and a late time of arrival at the node t max .Two consecutive flights are considered to be in conflict at the node when the time interval between the two of them crossing the node is less than the separation time. Let's consider two flights f and g with f leading.

∀ f, g ∈ F , c f lt( f, g) =        1 if |t f max -t g min | < t sep 0 otherwise
Link conflicts Link conflicts are essentially detected using wake turbulence category separation. Three situations are analyzed to detect those conflicts: at the entry node of the link, at the exit node of the link and catch-ups along the link itself. To detect entry conflicts, two different times are monitored at the entry of the link: an early time of arrival at the entry node t in min and a late time of arrival at the entry node t in max . Two consecutive flights are considered to be in conflict at the entry of a link when the distance between them is less than the required wake turbulence separation distance. Let's consider two flights f and g with f leading.

∀ f, g ∈ F , c f lt( f, g) =        1 if dist ori ( f, g) < wTCat sep ( f, g) 0 otherwise dist ori ( f, g) = V f min × (t f,in max -t g,in min
) With wTCat sep ( f, g) the wake turbulence separation distance between f and g. To detect exit conflicts, two different times are monitored at the exit of the link: an early time of arrival at the exit node t out min and a late time of arrival at the exit node t out max . Conflicts at the exit of the link are detected the same way as at the entry of the link. The distance between two consecutive flights is computed by using the exit times : dist dest ( f, g) = V g max × (t g,out mint f,out max ) Conflict detection at entry and exit only checks if the aircraft are properly separated when entering or exiting the link. However, within the link, faster aircraft can catch-up slower ones and cross them, thus implicating a loss of separation during the crossing. To detect catch-up conflicts, a comparison is made between the entry and the exit sequence of flights. Let's assume S l entry and S l exit respectively the sequences of flights at the entry and at the exit of the link l, pos : (F , S) → N the function that gives the position of a flight in a sequence. There is a catch-up conflict if :

∃ f ∈ F / pos( f, S l entry ) pos( f, S l exit )

Objective function

The mathematical modeling of the problem leads to a multi-objective optimization problem which aims at minimizing parameters linked to conflict occurrences, flight level modification, speed modification and also flights timing. The main objective is to reduce the global conflict count both on nodes and links. Let's consider c ob j the total conflict objective value.

c ob j = n∈N c n + l∈L c l
Two values related to flight level modifications are evaluated in the objective function:

f l ob j = f l changes + f l gap
where f l changes represents the total number of aircraft which decision implies FL change and f l gap is the summation of the difference between the final flight level and the preferred flight level of each flight. As for the flight level gap objective, we also want to minimize the gap between the initial speed and the final speed for each flight. For that, the speed objective s ob j is expressed as follows.

s ob j = f ∈F |∆V f |
Finally, we also want to make sure that the flights will arrive at their last node at the requested time with a reasonable (minimal) error margin. Thus, the time objective is expressed as follows.

t ob j = f ∈F |t f,last node max -t f,last node planned,max |
The global objective function f is the sum of all the objectives listed above.

f = c ob j + f l ob j + s ob j + t ob j
The complexity of our model mainly depends on the dimension of the problem (number of flights involved) and the range of the decision variables. If we consider n s options for the speed changes and n f l options for the FL changes, then for N f flight the size of the solution space is given by (n s × n f l ) N f . In our model, ∆V ∈ [-6, +3] and ∆FL ∈ [-2, +2], the size of the solution space can be expressed by (10 × 5) N f . To address such combinatorial complexity, we have developed a meta-heuristic based on a sliding-window simulated annealing.

Simulated annealing

This method, as described by [START_REF] -H Soliman | Modern Optimization Techniques with Applications in Electric Power Systems[END_REF], is based on the annealing process which is a physical process consisting of heating up a solid until it melts, then cooling it down in order to obtain a perfect crystallized form. During this process, the resulting structural properties depend on the residual energy in the material which is influenced by the rate of cooling. Thus, the cooling phase must be controlled in order not to get trapped in a locally optimal structure with high energy crystals, resulting in imperfections. In combinatorial optimization, the equivalent process is the simulated annealing which aims at finding a solution with minimal cost.

Principle

Let's consider the following elements [START_REF] Henderson | The theory and practice of simulated annealing[END_REF]. Let Ω be the solution space (the set of all possible solutions). Let f : Ω → R be an objective function defined on the solution space. The goal is to find a global minimum, ω * (ω * ∈ Ω such that f (ω) ≥ f (ω * ) for all ω ∈ Ω). The objective function must be bounded to ensure that P{Accept

ω as next solution} =        exp[-( f (ω -f (ω))/t k )] if f (ω -f (ω) > 0 1 if f (ω -f (ω) ≤ 0.
As illustrated in figure 3, each solution vector X which gather together the decision variables of aircraft is evaluated in a simulator process. This simulation is very powerful in the sense it can take into account many realistic aspects of the operational system. The result of the simulation generates the objective function which will be used by the optimizer (simulated annealing). The initial temperature is essential to define the behavior of the algorithm when it comes to the acceptance of solutions. We have set the initial temperature in order to get 80% of solutions to be initially accepted. We have used a geometric cooling law (t k+1 = α.t k with α in the range 0.8 -0.99.Efficient stopping criterion can prevent the simulated annealing algorithm from performing unnecessary computations [START_REF] Henderson | The theory and practice of simulated annealing[END_REF], thus reducing its global execution time.

In our case, three criteria are used to stop the cooling procedure: when the final temperature T f inal reaches a defined fraction of the initial temperature T init , when the evaluated criteria reach the value zero, or when the evaluated criteria remain constant over a defined number of iterations.

Neighborhood function

Small changes on a local solution are performed using a neighborhood function. The efficiency of simulated annealing is highly influenced by the neighborhood function used. In our case, each solution is a set of flight decisions and determining a neighbor solution consists only on randomly modifying a flight decision in our set. The method used to choose the random flight is similar to the roulette wheel selection [START_REF] Ji | Aircraft merging and sequencing problems in tma[END_REF]. On each flight decision, a performance criterion is evaluated, it corresponds to the sum of the total number of conflicts on the flight and the resulting delay. All the decision performances are added up together then a target value is randomly chosen between zero (0) and the total performance sum. Then, the cumulative sum of decision performances is calculated starting from the first decision until the cumulative sum reaches the previously chosen target value. The index at which it stops marks the flight decision on which a change will be done to generate a neighbor solution. Choosing the decision to change with this method guarantees that the selected flight is the one which is most likely to generate more conflicts and a large absolute delay. Given a flight decision, generating a neighbor solution consists in modifying the flight level delta ∆FL and the speed delta ∆V. Both modifications are done within their respective discrete range of values [-2, +2] and [-6, +3], with the same probability of occurrence p = 0.5.

Sliding window

The sliding window approach as proposed by [START_REF] Ji | Aircraft merging and sequencing problems in tma[END_REF] is a technique based on the receding horizon control. It consists into dividing the entire time horizon into smaller equal intervals and thus evaluating the state of the network within the small time intervals. The evaluation interval starts from the earliest time and progressively moves forward in time with a defined step until reaching the latest times. With this approach, it is not necessary to evaluate all the flights at once during the optimization since not all flights are involved in each sliding window interval. This method is more convenient in a dynamic environment with uncertainties and improves the computational time of the optimization process. This section presents the results of our algorithm on real-life traffic data at Paris CDG. Paris CDG TMA is accessible through four entry points: OKIPA, MOPAR, LORNI and BANOX. Our experimentation airspace is focused on the En-route portion of traffic arriving at CDG airport which corresponds to four areas each having a TMA entry point as last node as we can see in figure 5. All nodes (waypoints) except the four entry points are constrained with a 5NM separation minimum. As for the TMA entry points, the minimum separation between aircraft is 8NM. Concerning the wind data, we considered wind components with an angle of 30 degrees relative to the geographical north. Its intensity has been set within the range 7-10 kts all over the test airspace. A simple analysis of the network structure shows 742 different nodes with 817 links. This will contribute to increase the complexity and so the resolution time. In our scenario, we will consider a 24-hours traffic data of 690 flights corresponding to the filed flight plans arriving at LFPG on July 28th, 2016. The arrivals are distributed among the four entry points of CDG Terminal Area (TMA) with the majority of flights arriving at LORNI and OKIPA nodes (Table 1). There is also a distribution of medium and heavy wake turbulence category aircraft with a majority of medium aircraft.

Statistics

Our algorithm written in Java 8 has been executed on an Ubuntu 18.04 operating system PC equipped with an Intel Core i5-3230M processor (4 × 2.60GHz) and 8GB of memory. The simulated annealing algorithm has been tested with the parameters in Table 2. Given the length of the sliding window (2 hours) and the distribution of flights (Figure ??), each resolution step on a sliding window will handle around 100 flights on average. A first analysis of the flight data detected 407 conflicts. The optimization algorithm analyzed the 690 flights and has solved the sequencing and the conflicts within a computation time of 11 minutes and 23 seconds. It solved all the 407 initially identified potential conflicts, which gives an efficiency rate of 100%. It changed flight levels on a total of 127 flights: 67 climbs and 66 descents. As for the speed, 04 flights have been accelerated and 04 slowed down which make a total of 08 speed changes. At the end, a total absolute delays value of approximately 611 seconds has been generated for the flights with speed changes which gives an average of approximately 01 minute per modified flight (the 08 flights with speed modification). Very few flights have been affected by speed variations, also, the generated absolute delays, an average of 76 seconds per modified flight is still acceptable. Table 3 summarizes all those changes.

Conclusion

This paper introduced the work done on the sequencing of the air traffic in En-route segments influenced by constraint points. We saw that En-route controllers are getting more involved in the pre-sequencing of arriving flows when they are still in the cruising phase which causes an increase of their workload. On the other hand, airlines wish to have efficient flights with few flight level changes around a certain preferred vertical profile and also less maneuverings due to conflict avoidance. To solve those issues, we have developed a decision support tool which can assist controllers in their tasks for sequencing traffic and solving conflicts in En-route airspace. After reviewing the concepts and previous works related to our subject, we based our study on a mathematical modeling of the problem followed by an optimization algorithm in order to extract traffic sequences. Using the simulated annealing algorithm for optimizing flights decisions appeared to be a good choice given its efficiency and adaptability properties. A first trial of our solution on real traffic data over Paris airspace displayed a resolution ratio of 100% for conflict solving and an acceptable level of speed and flight level changes. Also, the generated delays due to the compliance with sequencing constraints and conflict resolution appeared within an acceptable range. Moreover, a preliminary version of the algorithm was able to generate flights instructions that can be directly applicable by controllers. Apart from this first test scenario, the solution as it has been designed is able to simultaneously provide En-route sequencing for several airports arrival flows. Even if it has b=not been used in the Paris CDG case, the algorithm can also manage time constraints for some points in the airspace (points connecting countries with letter of agreement). On the other hand, it can also be helpful for airports not yet equipped with an arrival management system as long as the constraints at TMA entry points are well defined. As another advantage, our solution will facilitate novel arrival techniques and procedures such as Continuous Descent Operation (CDO) and Point Merge System (PMS).

The authors confirm that they, and/or their company or institution, hold copyright of all original material included in their paper. They also confirm they have obtained permission, from the copyright holder of any third party material included in their paper, to publish it as part of their paper. The authors grant full permission for the publication and distribution of their paper as part of the EIWAC2019 proceedings or as individual off-prints from the proceedings.

Figure 1

 1 Figure 1 Simple Horizon Extension

Figure 2

 2 Figure 2 Protected area around the aircraft. If two areas overlap, there is a potential conflict

Figure 4 Figure 3

 43 Figure 4 Sliding window principle [8],W: the length of the sliding window, S : time shift of the sliding window

Figure 4

 4 illustrate how the sliding windows are generated all along the global optimization time interval.

Table 2

 2 Simulated annealing / Sliding window test parameters

	Criteria	Value
	Initial conflict count	407
	Solved conflicts	407
	Solved conflicts ratio	100%
	Climbed flights	61
	Descended flights	66
	Total FL changes	127
	FL changes ratio	18.41%
	Accelerated flights	4
	Slowed down flights	4
	Speed changes	8
	Speed change ratio	1.16%
	Total speed variation	35.029 kts
	Average speed variation	4.38 kts
	Total absolute delays	611 seconds
	Average absolute delays	76 seconds
	Execution time	11min 23sec

Table 3

 3 Optimization results on LFPG data

	Figure 5 Routes towards Paris CDG airport


Jim Sharples
email: jim.sharples@enac.fr

Gautier Hattenberger
email: gautier.hattenberger@enac.fr

Keywords: Unmanned Air Systems (UAS), Education, University

The French Civil Aviation University's (ENAC, École Nationale de l'Aviation Civile) main objectives are to train personnel specialised in air transport, and to help these specialists adapt to civil aviation's evolution, such as the advent of Unmanned Air Systems (UAS) in our skies. This paper shows that in ENAC, education in UAS is carried out in the same way than education in manned aviation topics. The training is divided in main subjects, which are then allocated to each student population according to their future roles: engineers, air traffic controllers, technicians, etc. Teaching is led by skilled personnel and backed up by modern facilities and equipment. This allows ENAC to provide a competent, specialised work force, ready to work on all aspects of the emerging and disruptive UAS-based operations.

I. INTRODUCTION

Many universities around the world do provide education on Unmanned Air Systems (UAS) or Unmanned Systems in general. This is for instance the case at Cranfield (UK), TU Delft (The Netherlands), Embry-Riddle Aeronautical University (USA) or RMIT (Australia), among others. Most of these education programs are focusing on Engineering aspects. Thus, UAS can be used to teach e.g. flight control [START_REF] Zhang | Using cutting-edge unmanned aerial vehicles ( uavs ) technology for flight controls courses teaching[END_REF] or data fusion [START_REF] Guerrero | Instrumentation of an array of ultrasonic sensors and data processing for unmanned aerial vehicle (uav) for teaching the application of the kalman filter[END_REF]. Drones can also be seen as a tool, therefore some of the teaching focuses on the end-user application, such as cartography [START_REF] Al-Tahir | Integrating uav into geomatics curriculum[END_REF], or meteorological instrumentation integration as done by the French Meteorological Office [START_REF] Cayez | Active teaching and learning of the use of RPA in meteorology[END_REF].

This paper shows how skills and knowledge in UAS are instilled into the French Civil Aviation University's (ENAC, École Nationale de l'Aviation Civile) 1 diverse student populations. The term "UAS" is here used to encompass the airborne vehicle (all sizes are considered), the command and control links with the ground station and/or remote pilot, and the ground station and/or remote pilot. UAS can be used for an ever-wider range of applications, some of which will be highlighted in this paper.

ENAC was created post-World War II as a training school for the French Civil Aviation Authority's (CAA, or DGAC in French) civil servants. It thus trained engineers, managers, air Work for this paper was supported by Groupe ADP and Sopra Steria 1 http://www.enac.fr/en traffic controllers and technicians. Since that time, ENAC has kept its original function but also grown to now train a wide variety of professionals in the field of civil aviation. It has also developed a strong research program, whose results are then integrated within the French CAA, but also worldwide.

One of ENAC's overarching aims is to help the different categories of professionals it trains to adapt to civil aviation's evolution.

ENAC divides its training in civil aviation in three main pillars: management (rules and regulations, Air Traffic Control); system engineering; and operations (of airlines, airports, ground systems). These competencies are then allocated to the student populations according to their future roles and diplomas, as described in the subsections below, and in figure 1.

For each of the student groups, we will detail how we craft UAS education tailored to their specific needs, thanks to the competencies of the teaching staff and to appropriate methods and tools.

II. ENAC GRADUATE ENGINEER (IENAC)

A. Grade and Roles

The French degree Diplôme d'Ingénieur is equivalent to a Master of Science (MSc) in Europe. These engineers will either serve as civil servants (10 % of a given year group) or work in the industry. Those that are civil servants can be involved in the definition, writing and/or application of the rules and regulations that govern civil aviation. Those that work for the industry will mostly be employed as Systems Engineers or experts in Aviation Operations.

B. UAS education for IENAC

1) Rules and Regulations: UAS rules and regulations are a relatively new topic (the first UAS laws were introduced in France in 2012) and an ongoing effort. In France, the CAA is the regulatory body in charge. It therefore stands to reason that the IENAC graduate engineers called to become civil servants within the CAA are taught this subject. But rules and regulations are just one of the pillars supporting UAS operations and market growth. A UAS traffic management system must also be built. Of course, the development of such a system intertwines deeply with the deployment of 2) UAS Design for IENAC: The teaching of system design from the UAS point of view has two purposes: firstly, provide knowledge and skills to future engineers targeting employment in the drone market; secondly, bolster education by giving students challenging and appealing practical projects.

The teaching methods used span a wide range of formats: talks (1 or 2 hours) on a specific topic, classic lectures (4 to 10 hours per topic) to provide more in-depth knowledge, practical and lab work (up to 20 hours) to gain practical skills and finally group projects (up to 80 hours) to fuse and apply to concrete problems the knowledge gained. IENAC students can thus benefit from a dedicated program of 80 hours to build foundation knowledge and, for those who wish and are selected, a research-oriented project of 80 hours. The selection process takes into account student motivation and academic results.

The knowledge imparted relies on several technical domains that are already covered in depth by the school's IENAC program. These include:

• system engineering, 2 http://drone-chair.enac.fr/

• control of dynamic systems,

• communication and control links,

• embedded systems and programming,

• conceptual design,

• payload and instrumentation.

To provide this broad spectrum of teaching, ENAC can rely on the professors and researchers of its UAS Research Program, who have been developing open-source hardware and software for UAS called Paparazzi3 since 2005 [START_REF] Hattenberger | Using the Paparazzi UAV System for Scientific Research[END_REF]. This system is of course deeply embedded in the lectures and practical projects and is one of their main pillars. In addition to that, an indoor flight arena is available on the campus. It has an effective volume of 11x9x9 meters and is equipped with an Optitrack4 camera-based positioning system, that enables 3D positioning with millimetric precision. This facility allows for safe flights and demonstrations. It is well suited to student projects, as a basic training of 2 hours is enough for them to become autonomous and fully use the system unsupervised. Figure 2 shows students testing a dedicated human-machine interface design for a system engineering project.

An example of one of the research-oriented projects, involving four students for 80 hours, is the study of autonomous landing of a fixed-wing aircraft in a net using ground-based ranging modules. Their tasks included writing up a state of the art report on the topic, leading a theoretical study of the expected positioning accuracy, and integrating the different components on a test aircraft as well as analysing the flight test results. The experimental plane, partially assembled by the students, can be seen in figure 3. This work served as the preliminary study for a paper published in a UAS-dedicated conference [START_REF] Hattenberger | Precision Landing for Fixed-wing UAV using Ultra-Wide-Band Ranging[END_REF].

III. AIR TRAFFIC SAFETY ELECTRONICS PERSONNEL (ATSEP)

A. Grade and Roles

Air navigation and air traffic management rely on safetycritical ground infrastructure: radars, radio navigation beacons, Control Tower IT networks, etc. ATSEP are responsible for this infrastructure: from assessing its needs and following industrial developments, to deploying, maintaining and upgrading its components, whilst also contributing to the regulations governing it. In France, ATSEP are civil servants, and they are trained exclusively at ENAC. Their diploma is equivalent to a Master's degree.

B. UAS education for ATSEP

1) UAS Design for ATSEP: ATSEP can be tasked with conducting inspection missions. For some of these operations, UAS are coming to the fore as potentially cost and timecutting tools, generating exhaustive data sets that complement traditional methods. For these inspection missions, ATSEP need to be competent in UAS design, particularly in the payload and instrumentation topics. Practical projects of 80 hours are proposed to groups of two students each, where the aim is to design and implement a payload for a designated mission. The UAS is provided by the UAS Research Program, along with the constraints for the payload (in terms of size, weight and power). For these topics, since the instrumentation can be quite specific, the ATSEP's teachers also take part in the projects alongside staff from the UAS Research Program.

2) UAS-based operations: For an inspection mission, the drone must of course be equipped for the task at hand (see Section III-B1), but the setup of the operation must also be carefully planned. This is where training in UAS-based operations fits in. As of now, it is limited to ATSEP who carry out specific projects on drone instrumentation (as mentioned in Section III-B1), and it encompasses most of the topics detailed in Section II-B1 and IV-B. These topics are tailored to fit a specific mission's requirements and are delivered by an engineer from the ENAC Research Chair on Drone Systems, who is also a certified drone pilot.

IV. AIR TRAFFIC CONTROLLERS (ATCOS) A. Grade and Roles

ATCos manage the airspace and the air traffic. Their aims are: to prevent collisions in the air and on the ground; to accelerate and organise air traffic; to provide information vital for maintaining flight safety; to take part in search-and-rescue efforts. In France, ATCos are civil servants, and they are trained exclusively at ENAC. Their diploma is equivalent to a Master's degree. ENAC also partly trains ATCos working for EUROCONTROL's Maastricht Upper Area Control Centre (MUAC).

B. UAS education for ATCos

In France, leisure drones are constrained by law to fly in segregated airspaces, where no other manned flights take place (under 500 ft above ground level and out of cities and airports). However, registered professional drones can be authorised, by the relevant bodies, to fly close to airports. Among the relevant bodies, Air Traffic Control can be called for. This means that air traffic controllers benefit from being educated in these very new systems: how they are operated, what are typical size weight and power and trajectory characteristics, how they are regulated, tracked, managed. Currently at ENAC, only ATCos for EUROCONTROL's Maastricht Upper Area Control Centre are educated in UAS by way of a two-hour lecture that brushes on all the topics detailed in the previous sentence. This lecture is delivered by an expert from the ENAC Research Chair on Drone Systems, who is also a certified drone pilot. The lecture encompasses a flight demonstration in the school's indoor flight arena based on ENAC's own open-source UAS. Additionally, members of ENAC's UAS Research Program can be called upon as experts for the assessment of ATCos' (either for the French CAA or for Maastricht) school projects.

V. ADVANCED TECHNICIAN DEGREE IN AVIATION (TSA) A. Grade and Roles

These technicians (TSAs in French) can either serve as civil servants or work in the industry. Civil servants are involved in a wide array of tasks for the French CAA: air traffic control on small aerodromes, defining aircraft arrival and departure procedures, airport inspection and surveillance, supporting the department in charge of incident/accident investigations, etc. In the private sector, TSAs can work for airport management companies, airlines, aircraft manufacturers. They lead tasks such as runway use optimisation or airworthiness checks.

B. UAS education for TSAs

Some TSAs might control air traffic in small aerodromes once they graduate; so they benefit from being educated in the same UAS topics as the ATCOs. Additionally, TSAs might be tasked with conducting inspection missions; so, like the ATSEPs, they benefit from being educated in UAS design and operations. Currently, this UAS education for TSAs is carried out by members of ENAC's UAS Research Program, who assist students with ongoing school projects, and can be called upon as experts for the assessment of said school projects.

VI. UAS SERVICES AND MANAGEMENT ADVANCED MASTER A. Grade and Roles

The French Advanced Master is a one-year course available to students that already hold a MSc. It is a means to specialise in a given topic. It is not open to ATSEP nor ATCos. ENAC developed this Advanced Master to meet UAS operators, manufacturers and air transport regulatory bodies' demand in high-level experts. Graduates are not civil servants, and can work as auditors, regulatory officers, research managers, business developers, designers, test managers, operations managers, data processing supervisors; in sectors such as air transport, audio-visual media, agriculture, energy production and distribution, etc.

B. UAS education for Advanced Master

ENAC is the first aeronautical university in Europe to open an Advanced Master on UAS management, applications and services. The course is divided in three blocks:

• Strategy (regulatory framework and UAS market analysis), • Technical Aspects (all the basic knowledge to understand UAS design from an engineering point of view), • Field Operations (pilot training, vehicle maintenance, mission setup and roll-out). During the first 6 months, the knowledge is imparted thanks to talks, lectures and case studies, delivered by a wide range of experts: members of ENAC's UAS Research Program, professors in air transport economics (one of ENAC's specialities), and professionals of the field (drone operators and/or manufacturers). Also, during these 6 months, a practical project runs parallel to the classroom lectures. The second block of 6 months is an internship in a company of the field.

VII. CONCLUSION

In this paper, we have explained how ENAC provides the personnel needed by manned air transport: engineers, air traffic controllers, technicians, etc. But manned aviation is currently facing a newcomer in its airspaces: Unmanned Air Systems. In order to increase the competencies of the professionals it trains, and help them adapt to civil aviation's evolution, ENAC also provides education in UAS. This training is modelled on the training given in manned aviation topics. The UAS subject is divided in three main categories: management, system design, and operations, which are then allocated to the different student populations according to their future roles and responsibilities. At the core of this teaching stands the UAS Research Program, its staff and its equipment, and its open source hardware and software platforms. These resources are enhanced by personnel with complementary skills and knowledge: in instrumentation, air transport economics, UAS operations, etc. This allows ENAC to provide a competent, specialised work force, ready to work on all aspects of the emerging and disruptive UAS-based operations.

Fig. 2 .

 2 Fig. 2. Human-Machine Interface student project inside the flight arena.

Fig. 3 .

 3 Fig. 3. Experimental setup for a student project.

http://paparazziuav.org

http://www.optitrack.com

ACKNOWLEDGMENT

The authors would like to thank Yannick Jestin (Head of UAS Research Program), and all the ENAC teachers involved in the training courses in UAS. The authors would also like to thank researchers Torbjørn Cunis and Murat Bronz for reviewing this paper.


Xavier Chermain

Frédéric Claux

Stéphane Mérillou

Glint Rendering based on a Multiple-Scattering Patch BRDF

Multiple-scattering glint integrator Chermain et al. [CCM18] Yan et al. [YHMR16

Introduction

Glittery materials, metallic paints, rough, brushed and scratched plastics and metals are surface types that are challenging to render in computer graphics. These materials have tiny mirror-likes details visible only under sharp and powerful lighting. A major feature of these materials is the spatial coherence of their hundreds of high frequency specular lobes, where small light or camera translations reveal the micro-surface. Using smooth normal distribution functions (NDF) to model material appearance is common for rendering specular surfaces, because these functions integrate well with common, field-proven shading models and frameworks. They are however ineffective at representing real life surfaces having complex micro-structures. Traditional surface filtering in computer graphics uses low-pass filters to reduce calculation costs, restricting NDFs to simple definitions, poorly replicating material appearance.

To capture the glints featured by specular surfaces, the reflectance model cannot use traditional Bidirectional Reflectance Distribution Functions (BRDF) based on smooth NDFs. More sophisticated models have to be used, such as glint integrators. This family of methods never simplify the micro-surface, and guarantees the rendering of all the micro-details falling into a screen pixel. The complexity of the material is fully taken into account and its related appearance is meant to be faithfully reproduced.

A Recent Interest. Rendering of specular micro-structures has received a growing interest in the last five years. The contributions of Jakob et al. [JHY + 14] and Yan et al. [YHJ + 14] have inspired many recent works exploring procedural materials [AK16, ZK16, WWH18], primitive-based BRDF [START_REF] Raymond | Multi-scale rendering of scratched materials using a structured sv-brdf model[END_REF] or normal map-based glint integrators [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF][START_REF] Xavier Chermain | A microfacet-based brdf for the accurate and efficient rendering of high-definition specular normal maps[END_REF][START_REF] Gamboa | Scalable appearance filtering for complex lighting effects[END_REF]. Each kind of method has associated strengths and weaknesses. Some methods specialize into a particular type of micro-structure (scratches, mirror-like flakes), sometimes achieving multiple-scattering effects and proposing good sampling routines. Others are more versatile, handling a large range of materials (glittery surfaces, brushed and scratched metals), often at the expense of memory consumption and rendering quality.

Remaining Issues. Most glint integrators based on normal maps or procedural surfaces simulate one bounce of light on the micro-surface. The rendered images produced with these integrators come out unnaturally dark, especially when using very rough surfaces (Figure 1, middle and right). A lot of surface normals are back-facing when the viewing or lighting direction is grazing. These unfortunate cases prevent light from being reflected. The second issue is an inefficient sampling procedure, yielding high sample weights, introducing isolated but nonetheless disturbing firefly artifacts (Figure 12). Current state-of-the art glint integrators randomly pick a normal of the micro-surface without taking into account the masking and interreflections that are inherent to the geometry of the micro-surface. At very grazing angles, about half of the samples are wasted and high weights compensate for these losses.

Position of our Work. The aim of this paper is to address all these issues. We choose to represent the microsurface by a normal map, for its versatility. Our representation is based on a dual normal model, inspired by the work of Schüssler et al. [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF]. A key difference is that we use a symmetric dual, not a tangent one (vertical), so as to not alter material roughness and related overall appearance. We take into account multiple-scattering by re-introducing the energy lost by a sole one-bounce BRDF formulation, a quantity that we analytically derive. This lost energy quantity is approximated and compensated for by using an energy-compensation BRDF, which is added to the single-scattering expression, eventually forming our final local BRDF. For importance sampling, we sample the single-scattering BRDF and keep the sample only if it is not shadowed or transmitted (Figure 5). Otherwise, we sample the energy-compensation BRDF, to account for multi-scattering.

Our local BRDF is then used as part of a footprint or patch based glint integrator. Our patch-BRDF is a discrete weighted sum of local BRDF instances, and is energy conserving.

Our contributions in the field of glint rendering are the following:

• We propose the first multiple-scattering glint integrator based on a specular normal map.

• The proposed patch-BRDF does not have the rendering-quality related issues of the previous methods: black holes on the surface, dark overall appearance and fireflies.

• The sampling procedure is optimal: no samples are wasted, the probability density function (PDF) is analytic and sample weights are one when the surface is non-absorptive.

Our representation is Jacobian free, i.e. it does not take into account surface curvature. It is asymmetric, therefore it limits the use of our method to forward path tracing. Our BRDF compensation scheme is restricted to specular surfaces. The dual normal approach may change the look of some materials in unwanted ways, even though this has not been significant for the examples in the paper. After a detailed review of the related works (Section 2), we present our local multiple-scattering BRDF in Section 3. We use this local model to build the patch-BRDF in Section 4, which filters the surface and guarantees the spatial coherence of the glints. We evaluate our method and discuss the results in Section 5.

Related Work

Our work is related to three different areas of research in computer graphics: glint integrators, microfacet-based normal mapping and multiple-scattering BRDFs. We review here the most relevant works in these areas.

Glint Integrators

Glint integrators can be seen as special types of BRDFs where specular micro-structures are rendered using either geometrical primitives randomly placed on the surfaces or high definition normal maps. Instead of relying on stochastic Monte Carlo sampling to estimate shading, glint integrators calculate all the radiance leaving a ray footprint (pixel-bound for primary rays). The BRDF associated with the footprint is usually referred to as a patch BRDF.

Stochastic surfaces. The first method making use of stochastic specular micro-structures is the method of Jakob et al. [JHY + 14] and focuses on glittery materials. The micro-surface is composed of mirror flakes, randomly placed on the surface. These are used to compute a discrete, patch-related NDF used as part of a microfacet-based BRDF. An acceleration hierarchy speeds up contributing reflector queries. Their method only models single-scattering and micro-flakes can be back-facing and consequently dismissed. The BRDF normalization is global, not footprint (patch) based. Importance sampling is also not optimal, as non-contributing reflectors may inadvertently be sampled. The extension proposed by Atanasov and Koylazov [START_REF] Atanasov | A practical stochastic algorithm for rendering mirror-like flakes[END_REF] alleviates this issue, at the price of going through all the flakes falling into a footprint, which is inefficient. Other works concentrate on improving performance [START_REF] Wang | Fast global illumination with discrete stochastic microfacets using a filterable model[END_REF]. Selected methods offer real time performance [START_REF] Zirr | Real-time rendering of procedural multiscale materials[END_REF][START_REF] Wang | A robust and flexible real-time sparkle effect[END_REF] but cannot really be used for photorealistic rendering.

The work of Raymond et al. [START_REF] Raymond | Multi-scale rendering of scratched materials using a structured sv-brdf model[END_REF] models specular scratches randomly distributed over a surface. Users have control over the profile, the orientation and the density of the scratches. They simulate multiple-scattering by pre-computing the reflectance distributions of scratch profiles with a 2D ray tracer. The full 3D-BRDF is reconstructed by interpolating these 2D-BRDFs. Mipmapped filtering is used to calculate the relative area of the scratches in a footprint. Their method gives realistic results, its importance sampling procedure is efficient, multiple-scattering is taken care of, but they only handle scratched surfaces.

Glint Integrators based on a Normal Map. As opposed to methods using procedural surfaces, glint integrators based on normal maps provide more flexibility over material appearance. Scratched, brushed, rough metals and glittery objects can be modeled by changing the normal map contents. The work of Yan et al. 2014 [YHJ + 14] was the first to introduce this micro-surface representation, obtaining realistic renderings of a wide range of materials. Several papers improve the performance of their method by either using 4D NDF mixtures [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF] or spherical harmonics [START_REF] Gamboa | Scalable appearance filtering for complex lighting effects[END_REF]. The work of Chermain et al. [START_REF] Xavier Chermain | A microfacet-based brdf for the accurate and efficient rendering of high-definition specular normal maps[END_REF] focuses on rendering quality and proposes footprint-relative normalization factors. The importance sampling procedure is efficient since it uses footprint-based NDFs, however visibility information is not used, causing fireflies and zero-radiance paths. The method of Gamboa al. [START_REF] Gamboa | Scalable appearance filtering for complex lighting effects[END_REF] does not use importance sampling. Instead, they efficiently integrate all the incident radiance, at the cost of a large bunch of pre-computed data. The common shortcoming of these methods is their single-scattering formulation, causing black overall appearance for rough surfaces. They also rely on classic normal mapping, where back-facing normals lead to black holes on the surface. Our work tackles these issues.

Microfacet-Based Normal Mapping

Classic normal mapping causes black fringes on surfaces and requires complex re-normalization factors, as demonstrated in the work of Schüssler et al. [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF]. To solve these issues, they propose a microfacet-based BRDF faking normal perturbations, using tangent facets to seal the micro-surface. The tangent facets make really continuous an otherwise discontinuous normal map composed of unrelated, individual discrete facets. They rely on random walks to model infinite multiple-scattering on the micro-surface, and derive an interesting closed-form expression for the two-scattering order. They obtain much more realistic results than classic normal mapping. Our model is inspired by their work but uses a symmetric dual normal, more appropriate for glints rendering that would otherwise be negatively impacted by vertical facets (Figure 2). Another major difference is that Schüssler et al. use expensive random walks to account for multiple-scattering, whereas we estimate the energy going into indirect bounces and compensate for it by using an alternate (or energy-compensation) BRDF to increase performance. Figure 3 gives a quick hint at the differences between classic normal mapping, the model of Schüssler et al. and our model.

Multiple-Scattering BRDFs

Recently, the multiple-scattering of light over micro-surfaces has gained a lot of interest in the graphics community. This is thanks to the work of Heitz et al. [START_REF] Heitz | Multiple-scattering microfacet bsdfs with the smith model[END_REF] on microfacet BSDFs using the Smith model. They provide ground truth results by matching stochastic reflectance functions with simulated data. The bottleneck of their method is the random walk approach that requires several evaluations to provide converged results. Recent works propose more practical multiple-scattering reflectance models, using distributions of specular Vgrooves [LJJ + 18, XH18, SPSH19]. These V-cavity-based methods are too expensive for our purpose. Finally, to efficiently model multiple-scattering, severals works directly re-integrate the energy lost by a single-scattering formulation, focusing on either reciprocity [CK17, Hil18, FA19] or efficient sampling [START_REF] Turquin | Practical multiple scattering compensation for microfacet models[END_REF]. We also use energy compensation, but do it in the context of glint rendering using microfacet-based normal mapping. Our method offers robust importance sampling, and a customizable energy-compensation BRDF.

Local Multiple-Scattering BRDF

Our model is built around the definition of a local and a patch-wide BRDFs. We start with the formulation of the local BRDF.

Micro-Surface Model Overview

Our micro-surface model is based on an interpolated, continuous normal map typically built from a discrete texture map. Each position u over the map has a perturbed normal ω p as well as an associated dual normal ω s . This dual normal is symmetric to ω p with respect to the geometric normal ω g , forming a V-shaped cavity [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF]. Right: our model. The symmetric facets ω s do not change overall material roughness the way tangent facets ω t do. For glint rendering, where the normal map is highly repeated, tangent facets would alter too much the appearance.

u u 1 Classic u u 1 [SHHD17] u Ours ω p (u 1 ) ω p (u 1 ) u 1 ω p (u 1 ) ω s (u 1 ) ω t (u 1 ) ω g
(Figure 3). Locally, at any position u over the surface, the orientation of the micro-normals ω p and ω s are used to derive the probability that an incident or outgoing ray direction hits either normals, taking into account the shading induced by its symmetric counterpart. All V-cavities have the same micro-BRDF. We detect situations incurring a loss of energy caused by local single-scattering within the V-cavity and account for this loss by also evaluating an additionnal BRDF. Stochastic Monte-Carlo rendering leverages this model to importance sample the surface locally.

Normal Distribution Function and Projected Areas

Normal distribution functions (NDF) are widely used to model micro-surfaces. They give the density of microfacets ω m at a surface position u. In our case, this function is composed of two Dirac distributions centered around the perturbed normal ω p and the symmetric normal ω s (we drop the dependance on u for the rest of the section):

D(ω m ) = 1 2 δ ωp (ω m ) ω p • ω g + 1 2 δ ωs (ω m ) ω s • ω g . ( 1
)
This NDF is a Cook-Torrance V-cavity NDF. Its projected area onto the geometric normal is exactly one [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]. The projected areas of a facet and its dual towards the observation direction ω o are

a p (ω o ) = ω o , ω p ω p , ω g , a s (ω o ) = ω o , ω s ω s , ω g (2)
where ω 1 , ω 2 is the clamped dot product. From these definitions, the probabilities λ p (ω o ) and λ s (ω o ) that an observation ray intersects ω p or ω s can be derived (Figure 4):

λ p (ω o ) = a p (ω o ) a p (ω o ) + a s (ω o ) , λ s (ω o ) = a s (ω o ) a p (ω o ) + a s (ω o ) . (3)
The terms λ p and λ s will be used to evaluate (Section 3.4) and sample (Section 3.6) the BRDF.

λ p (ω o ) λ s (ω o ) ω i ω o G 1 (ω i , ω p )
Figure 4: Illustration of the intersection probabilities and the shadowing function.

Masking Function

The masking function gives the proportion of microfacets ω p or ω s that are visible from the observation direction ω o . It is called the shadowing function when using the incident direction ω i as parameter (Figure 4). The V-Cavity masking function is [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]:

G 1 (ω o , ω p )= H( ω o , ω p )min 1, ω o , ω g 1 2 a p (ω o ) + 1 2 a s (ω o ) = H( ω o , ω p )min 1, 2 ω o , ω g ω p , ω g ω o , ω p + ω o , ω s (4)
for the facet with normal ω p . For the micro-normal ω s , it is

G 1 (ω o , ω s ) = H( ω o , ω s )min 1, 2 ω o , ω g ω s , ω g ω o , ω p + ω o , ω s (5)
where H(-) is the Heaviside function.

As in Schüssler et al. [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF], we assume no correlation of masking and shadowing, giving a simple separable masking-shadowing function (here for the micro-normal ω p ):

G 2 (ω o , ω i , ω p ) = G 1 (ω o , ω p )G 1 (ω i , ω p ) (6)
where ω i is the incident direction. This choice has an important implication for the importance sampling procedure (Section 3.6). We have to determine when a ray is shadowed and with this assumption, we can directly use the shadowing function G 1 (ω i , ω m ) instead of using explicit ray-facet intersections [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF].

Single-Scattering BRDF

In the context of glint rendering, the material of the two facets is very specular. Assuming that the two facets of our micro-surface have the same micro-BRDF f m , the cosine-weighted single-scattering BRDF f 1 for our V-cavity is

f 1 (ω o , ω i , ω p ) ω i , ω g = λ p (ω o )f m (ω o , ω i ) ω i , ω p G 1 (ω i , ω p ) + λ s (ω o )f m (ω o , ω i ) ω i , ω s G 1 (ω i , ω s ). (7)
The BRDF expression f 1 is the result of the integral of the contribution of each micro-normal ω m over the hemisphere Ω [START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF]. The shadowing function and intersection probability of each facet (Figure 4) weight the corresponding facet's cosine-weighted micro-BRDF f m , rotated on either ω p or ω s . Note that f 1 does not create any energy and is symmetric as long as the same goes for f m .

Multiple-Scattering BRDF

In this section, we demonstrate how we estimate f ∞ , the multiple-scattering BRDF, for a single cavity. Its definition is based on f 1 and an energy-compensation BRDF f 2+ meant to represent second-or-more order scattering. We estimate the energy 1 -E 1 going through this latter BRDF, and weight its evaluation using this value.

Energy Term Assuming that the micro-BRDF f m is perfectly specular and non-absorptive, the directional albedo E 1 of f 1 , also called the energy term, is

E 1 (ω o , ω p ) = Ω f 1 (ω o , ω i , ω p ) ω i , ω g dω i = λ p (ω o )G 1 (ω r (ω o , ω p ), ω p ) + λ s (ω o )G 1 (ω r (ω o , ω s ), ω s ) (8)
where ω r (ω 1 , ω 2 ) is the reflected direction of ω 1 around the normal ω 2 . The term E 1 gives the radiance that escapes the micro-surface after the first bounce of light in the outgoing direction ω o . For the single-scattering BRDF f 1 , this value ranges from 0 to 1 for any different couple of (ω o , ω p ).

In order to achieve 100% energy conservation, we define our local multiple scattering BRDF f ∞ as

f ∞ (ω o ,ω i , ω p ) = f 1 (ω o , ω i , ω p ) + (1 -E 1 (ω o , ω p ))f 2+ (ω o , ω i ). (9)
Weighted by the cosine, the BRDF integrates to one over incident directions

Ω f ∞ (ω o , ω i , ω p ) ω i , ω g dω i = E 1 (ω o , ω p ) + (1 -E 1 (ω o , ω p )) Ω f 2+ (ω o , ω i ) ω i , ω g dω i = 1 + E 1 (ω o , ω p ) -E 1 (ω o , ω p ) = 1 (10)
only if the energy-compensation BRDF f 2+ (ω o , ω i ) ω i , ω g does. The choice of f 2+ is arbitrary and is discussed in Section 5.2. Our formulation for f ∞ has a simple closed-form expression and can be evaluated efficiently. It is non-symmetric as can be seen in Equation 9, swapping ω o and ω i yields a different

E 1 term (E 1 (ω o , ω p ) = E 1 (ω i , ω p ))
. This is only a deal breaker when using bidirectional path tracing, and does not negatively impact visual appearance otherwise.

In practice, we mollify the Dirac distribution of f m to a NDF with tiny roughness α m [YHJ + 14, YHMR16, CCM18] in order to evaluate f 1 (Equation 7). Note that the energy term of Equation 8 still guarantees energy conservation with a mollified Dirac. Section 5.1 shows that f ∞ passes the white furnace test when the f m NDF has a tiny roughness, validating this approximation.

Fresnel Term. The above model does not incorporate Fresnel effects, which give the object its color and is, in the case of multiple scattering, responsible for the saturation effect that is visible on rough surfaces. We therefore need to choose the absorption coefficient for the energy-compensation BRDF f 2+ , but it is not trivial. We choose to approximate it by the Fresnel term F 1 of the micro-BRDF f m . The local multiple-scattering BRDF f ∞ for an absorptive normal map is now

f ∞ (ω o , ω i , ω p ) = f 1 (ω o , ω i , ω p ) + (1 -E 1 (ω o , ω p ))F 1 (ω o , ω g )f 2+ (ω o , ω i ), (11)
where F 1 (ω o , ω g ) colors the non-absorptive f 2+ BRDF the same way it does for f 1 . Note that Turquin [START_REF] Turquin | Practical multiple scattering compensation for microfacet models[END_REF] performs a similar approximation for its compensated microfacet BRDF. His model exhibits similar results as the ground truth [START_REF] Heitz | Multiple-scattering microfacet bsdfs with the smith model[END_REF].

Importance Sampling the Local BRDF

To importance sample f ∞ , we first sample the cosine-weighted single-scattering BRDF f 1 by choosing a facet using the intersection probabilities λ p (ω o ) and λ s (ω o ). A direction ω i is then picked by importance sampling the micro-BRDF f m of the sampled micro-normal. If the sampled direction ω i is shadowed (Figure 5, middle) or goes inside the surface (Figure 5, left), multiple-scattering occurs and f 2+ is sampled. The probability that the sampling of the single-scattering BRDF f 1 fails, i.e. ω i is shadowed or transmitted, is exactly equal to 1 -E 1 (ω o , ω p ). Our importance sampling procedure is summarized in Algorithm 1 (U is a uniform random number in the [0, 1[ range) and is illustrated in Figure 5. It is guaranteed to give weights w between 0 and 1 as long as f m and f 2+ also do. The weight will always be 1 if there is no Fresnel effect, considering f m has a near-zero roughness. Otherwise, w is multiplied by the Fresnel term. In all cases, arbitrary large weights cannot happen and firefly artifacts are avoided. 4 Multiple-Scattering Patch BRDF

Algorithm 1: Sample f ∞ (ω o , ω i , ω p ) ω i , ω g ω m ← [ U < λ p (ω o ) ? ω p : ω s ] // Choose a facet (ω i , w) ← sample f m (ω o , ω i , ω m ) ω i , ω m // Sample f m if ω i • ω g ≤ 0 or U ≥ G 1 (ω i , ω m ) then /* Ray goes inside the surface or is shadowed */ (ω i , w) ← sample f 2+ (ω o , ω i ) ω i , ω g // Sample f 2+ end ω o ω o ω o Transmitted Unsuccessful sampling of f 1 → sampling of f 2+ Successful sampling of f 1 Shadowed Single-scattering
In this section, we introduce the theoretical background to define, evaluate and sample patch BRDFs (P-BRDFs) based on a normal map. We then explain how to efficiently compute our P-BRDF based on the local multiplescattering BRDF f ∞ introduced in the previous section. We close this section with our importance sampling procedure.

Patch BRDF Definition

Our P-BRDF computes the density of outgoing radiance for a flat geometric surface P instead of a surface point u. The patch P matches the ray footprint. It is represented by a normalized low-pass filter k P , i.e.

R 2 k P (u)du = 1 and antialiases the radiance of all the glints in the footprint. Several techniques for surface filtering [Ige99, SW01, BYRN17] give methods to construct P during path tracing.

Our P-BRDF f P is defined as an integral of local BRDFs weighted by the ray footprint k P (u) at each position u:

f P (ω o , ω i ) = R 2 f (ω o , ω i , ω p (u))k P (u)du (12)
where f (ω o , ω i , ω p (u)) is the surface's BRDF perturbed by ω p (u). Note that this sharply contrasts previous works [JHY + 14, YHJ + 14, YHMR16, CCM18, GGN18] for which the P-BRDF formulations use a P-NDF.

Patch-BRDF Evaluation

By replacing f with our multiple-scattering BRDF f ∞ perturbed by ω p (u), we get

f P (ω o , ω i ) = R 2 f ∞ (ω o , ω i , ω p (u))k P (u)du = f 1 (ω o , ω i , P) + E 2+ (ω o , P)F 1 (ω o , ω g )f 2+ (ω o , ω i ). (13)
The Fresnel term F 1 (ω o , ω g ) and the energy-compensation BRDF f 2+ (ω o , ω i ) can be taken out of the integral because these terms do not depend on u. The term f 1 (ω o , ω i , P) is the single scattering patch-BRDF, defined as

f 1 (ω o , ω i , P) = R 2 f 1 (ω o , ω i , ω p (u))k P (u)du. (14)
The term E 2+ (ω o , P) is the energy lost by the single-scattering BRDF f 1 in the footprint, and its expression is

E 2+ (ω o , P)= R 2 1 -E 1 (ω o , ω p (u))k P (u)du = 1 - R 2 E 1 (ω o , ω p (u))k P (u)du. (15)
We now develop how we calculate f 1 (ω o , ω i , P) and E 2+ (ω o , P).

Single-Scattering Patch BRDF Evaluation

Until now our model, based on Equation 14, implies an infinite set of local BRDFs. In practice, we use discrete local BRDFs and assume a finite set of m perturbed normals ω p :

f 1 (ω o , ω i , u) = m j f 1 (ω o , ω i , ω p (u j ))k j (u). (16)
Seed points u j separated by a distance h are uniformly distributed across the normal map. Each BRDF at position u j is weighted by its distance to u to ensure smooth transitions between two adjacent local BRDFs. As our footprint k P is Gaussian, we use a Gaussian weight k j for the BRDF in order to have an analytic weighting of f 1 . The standard deviation of k j is set to σ h = h/ √ 8 log 2 as in the methods of Yan et al. and Chermain et al. [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF][START_REF] Xavier Chermain | A microfacet-based brdf for the accurate and efficient rendering of high-definition specular normal maps[END_REF] to achieve an equal contribution of two adjacent BRDFs when the evaluation point u is halfway between the two (Figure 6). Substituting the sum expression of f 1 (Equation 16) into the single-scattering patch-BRDF (Equation 14) gives

k P (u) u 1 u k 1 (u) k 6 (u) u 6 . . . . . .
f 1 (ω o , ω i , P)= R 2 m j f 1 (ω o , ω i , ω p (u j ))k j (u)k P (u)du = m j f 1 (ω o , ω i , ω p (u j )) R 2 k j (u)k P (u)du = m j f 1 (ω o , ω i , ω p (u j ))W j ( 17
)
where the weight W i is the analytic solution of the integral of the product of two 2D Gaussians, here the local weight k j and the ray footprint k P .

Acceleration Hierarchy. As explained in Section 3, the local BRDF f 1 has two symmetric facets with a near-perfectly specular material f m . This micro-BRDF is a delta Dirac function, but is modeled for practical reasons by a microfacet-based BRDF with very low roughness α m . This implies that f 1 is different than zero only when the half vector ω h = ωo+ωi ||ωo+ωi|| approximately matches the perturbed normal ω p or its symmetric ω s , and is used to accelerate the sum evaluation of Equation 17 by pruning out non-contributing discrete BRDFs.

A 4D Bounding Volume Hierarchy (BVH) is built using the discrete set of m perturbed normals ω p . Each normal has a position u j , a slope pj and a 4D bounding box defining the (u, p) domain where the perturbed facet of f 1 reflects a non-negligible amount of light, for a given couple of position u and half vector slope h. The 4D bounding box is (u j ±3σ h , pj ±3σ m ), where σ m = α m / √ 2, assuming that the contribution is zero beyond 3σ. This is mostly similar to the hierarchies used in Yan et al. [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF] and Chermain et al. [START_REF] Xavier Chermain | A microfacet-based brdf for the accurate and efficient rendering of high-definition specular normal maps[END_REF]. Because of the use of a symmetric facet, we have to intersect the 4D BVH with (u, h) as in the previous methods, but also with (u, -h) to account for the contribution of the ω s micro-normals.

Energy Lost in the Patch

Evaluating the energy lost in the patch requires us to process all the perturbed normals ω p falling into the ray footprint. For performance reasons, we assume that normals ω p follow a Gaussian distribution D of normals ω p having an anisotropic footprint roughness (α x (P), α y (P)). In other words, we change the integral domain from spatial positions u to perturbed normals ω p , similarly to Heitz et al. [START_REF] Heitz | Filtering color mapped textures and surfaces[END_REF]

E 2+ (ω o , P) = 1 - R 2 E 1 (ω o , ω p (u))k P (u)du ≈ 1 - Ω E 1 (ω o , ω p )D(ω p , α x (P), α y (P))dω p ≈ 1 -cos 2 (φ o ) Ω E 1 (θ o , ω p )D(ω p , α x (P))dω p + sin 2 (φ o ) Ω E 1 (θ o , ω p )D(ω p , α y (P))dω p . ( 18
)
The terms θ o and φ o are respectively the polar and azimuthal angles of observation direction ω o . We use φ o to weight the spherical integrals implying an isotropic distribution D(ω p , α) of normals ω p . By doing so, we can precompute integrals of the form Ω E 1 (θ o , ω p )D(ω p , α)dω, and store the values into a small 2D array (128 × 128 floats, Figure 7) for different couples of roughness α and polar angles θ o . To calculate α x (P) and α y (P) at render time, we use the method described in LEADR mapping [DHI + 13], based on the first and second moments of a slope map.

Importance Sampling the Patch-BRDF

Importance sampling is used to efficiently solve the patch scattering equation:

L(ω o ) = Ω L(ω i )f P (ω o , ω i ) ω i , ω g dω i = Ω L(ω i ) R 2 f (ω o , ω i , ω p (u)) ω i , ω g k P (u)dudω i (19)
where L(ω i ) is the incident radiance from direction ω i . The sampling of a patch-BRDF f P (ω o , ω i ) is done by sampling (u, ω i ) with the following probability density function (PDF)

PDF(ω i , u) = f (ω o , ω i , ω p (u)) ω i , ω g k P (u). (20)
The marginalized random variable u is first chosen, followed by the conditional random variable ω i |u. For the first step, the marginalized PDF is

PDF(-, u)= Ω f (ω o , ω i , ω p (u)) ω i , ω g dω i k P (u) = E(ω o , ω p )k P (u) (21)
where E is the energy term. This simplifies to k P (u)

only if the local BRDF f is 100% energy conserving, that is, if Ω f (ω o , ω i , ω p (u)) ω i , ω g dω i = 1 for all u ∈ R 2 .
Because our local BRDF honors this condition, the sampling first consists in sampling the ray footprint k P (u). During the second step, the conditional PDF is f (ω o , ω i , ω p (u)) ω i , ω g at the sampled position u.

To our knowledge, no other method has such a simple and efficient sampling procedure. Other methods either disregard the observation direction when sampling the patch [YHJ + 14, JHY + 14, YHMR16, CCM18], or need to walk through all normals if they want to do it [AK16] (Figure 8), which is inefficient.

Results and Discussion

In this section, we evaluate our method and compare it with-state-of-the-art methods. All results have been rendered with the PBRT-v3 [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF] renderer on a dual Intel Xeon Gold 6138 CPU computer. Figures 2 and9 do not use pixel footprints, while all other figures in the paper do. After the first light bounce, we never compute any footprint and rely on a smooth BRDF using the overall roughness and Fresnel term of the normal map, as in Jakob et al. [JHY + 14]. We only use forward path tracing since our BRDF is asymmetric. All figures use h = 1 (Section 4.2.1) except Figures 14 and15 (scratches part only).

Validation of the Local BRDF

In this section, we validate our multiple-scattering local BRDF f ∞ . To do so, we compare renderings and white furnace tests of classic and microfacet-based normal mapped surfaces. In Figure 9, the normal map is a Beckmann surface (α = 0.5). We compare f 1 and f ∞ to classic normal mapping. Black fringes are visible on the surface with this latter method. The single-scattering BRDF f 1 partially eliminates these artifacts while the multiple-scattering BRDF f ∞ makes them completely disappear. Furthermore, the BRDF f ∞ passes the white furnace test. We use the following parameters: f m is a classic microfacet-based BRDF with a Beckmann NDF (α m = 0.02). Its Fresnel term F 1 models a gold metallic material. The energy-compensation BRDF f 2+ is lambertian (diffuse), with an albedo equal to F 1 . Notice how the diffuse compensation and coloring incurred by F 1 is natural and plausibly simulates multiple-scattering on the micro-surface. For all other results, we use α m = 0.02 and a diffuse f 2+ , unless otherwise stated. In Figure 9, bottom rows, the normal map is repeated 20 times more than in the top rows, with other parameters unchanged. The back-facing normals associated with classic normal mapping aggregate, leading to an overall dark appearance (Figure 9, bottom left). Our local multiple-scattering BRDF f ∞ does not have this issue (Figure 9, bottom right).

Patch-BRDF Results

Figure 10 shows different materials that our method can model. An animation is provided in the supplemental materials to evaluate the temporal coherence offered by our method. Brushed metals with dents can be obtained using a Beckmann surface instance with very high anisotropy (α x = 0.5, α y = 0.01, see the bigger spoon). The knife uses a texture with scratches in only one direction, while the dessert spoon uses random directions, reproducing a completely different appearance. The forks have a specular f 2+ matching a plastic appearance. The green fork uses a Beckmann normal map with high isotropic roughness α = 1, and the purple one uses a Voronoï texture with α = 1.41. Finally, the silver spoon has a Beckmann surface with α = 0.1. Memory footprint. For the normal map representation and the roughness estimation in the footprint, we need two mipmapped textures, one for the normal map slopes and another one for the second moments. All our results use 2048 × 2048 mipmapped textures, each requiring 96 MiB. The lost energy approximation of Section 4.2.2 takes up 128 × 128 floats (64 KiB). The acceleration hierarchy of Section 4.2 uses a similar size than the one used by the methods it is based on [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF][START_REF] Xavier Chermain | A microfacet-based brdf for the accurate and efficient rendering of high-definition specular normal maps[END_REF] for the same value of h (512 MiB for h = 1, 2048 MiB for h = 0.5).

Impact of the Patch-Energy Approximation Our P-BRDF relies on an approximation of the energy lost in the patch (Section 4.2.2) to speed up evaluations. Figure 11 compares this approximation with the accurate quantity, directly calculated with a discrete weighted sum of f ∞ across the whole ray footprint and fully passing the white furnace test. Our approximation method provides a 2× speed up, except for the scratched surface which includes a lot of planar normals. While the white furnace tests don't come out perfect, they are significantly better than those of previous methods, and we offer a plausible appearance, without any black holes on the surface, unsaturated colors or overly dark general appearance.

Importance Sampling. The weights of our importance sampling procedure are never greater than one, guaranteeing renderings without unnaturally bright pixels (Figure 12, right). On the contrary, previous methods In the bottom rows, the normal map repetition is twentyfold, leading to a coarse Beckman specular lobe. We use 131,072 spp to render the majority of glints and achieve noise-free renderings.

can give arbitrary large sampling weights to compensate for zero-radiance samples (Figure 12, left). When the surface is partially shadowed as in Figure 12, no tiny reflections should be visible. This is the case for our model and not for the previous method.

Choosing the Energy-Compensation BRDF. Most figures in the paper use a diffuse energy-compensation BRDF f 2+ to model multiple-scattering. This choice gives results comparable to the previous methods, where their overall, dark appearance can be seen as a diffuse material with a low albedo. Figure 13 shows the difference between the use of a diffuse f 2+ and a specular f 2+ (α = 0.02). Metallic paint can be modeled with a nearperfectly specular f 2+ , giving the impression that metallic flakes are coated with a transparent varnish. Using a lambertian f 2+ suits materials with rough appearance.

Comparison with Methods using Gaussian Mixtures. We compare the renderings of glittery surfaces using two different methods in Figure 1 (and supplemental video). The methods of Chermain et al. and Yan et al. suffer from backfacing normals, related black fringes and dark overall appearance. The outcome of the white furnace test is significantly better with our method. Figure 1 uses the flat variant of the Gaussian mixtures [START_REF] Yan | Position-normal distributions for efficient rendering of specular microstructure[END_REF], that is best suited to represent glittery surfaces (constant micro-BRDF roughness α m = 0.02).

Figure 14 (and supplemental video) shows an object with a structured micro-surface rendered using different methods. Gaussian mixtures-based methods exhibit unnatural dark or bright scratches. They can be observed in the right area of the orb, due to back-facing normals or incorrect P-BRDF normalization. Structured microsurface renderings use h = 0.5 (Figures 14 and15) for our method for a reason that we now explain.

The models of Chermain et al. and Yan et al. reconstruct an interpolated normal map with NDF mixtures, where each NDF is a Gaussian with a variance that best fit locally the curvature of the input surface. Per-Gaussian reconstruction kernels k j further interpolate BRDF values between one another at render time, resulting in a surface model fundamentally continuous, even when the surface has relatively high frequency contents locally. Our model resembles the Jacobian-free, flat element variant (with fixed Gaussian roughness) of Yan et al., where normals are not interpolated continuously. It defines two normals for a given surface position instead of one. Equation 8 of our model forbids the use of high local roughness for the micro-BRDF. Our method therefore solely relies on the k j reconstruction filter to interpolate neighboring BRDF values at render time. Consequently, local high frequency contents are approximated less accurately than when using mixtures of Gaussians with varying roughness. Scratches appear sharper with our model, less smoothed down. This issue can be alleviated by reducing the h sampling distance (increasing the density of V-cavities). On the upside, the intrinsic shadowing of V-cavities naturally weights the intensities of the glints, favoring glints continuity when animating. Note that high roughness severely impacts performance of previous methods, causing BVH node extents to be extremely large. For example, the orb of Figure 14 The quality of the lost energy approximation (Section 4.2.2) depends on the size of the ray footprint. Large ray footprints lead to favorable Gaussian approximations of normals distributions (Equation 18). For close-up views, smaller foorprints cause lost energies to be underestimated, leading to poor white furnace tests (Figure 15). At low scales, roughness is not a reliable indicator for the lost energy estimation. Other methods exhibit even poorer results. Yan et al. normalize their footprint-dependant BRDF in a footprint-independent way. Chermain et al. suffer from the lack of multiple-scattering.

Conclusion and Future Works

In this work, we have presented the first practical multiple-scattering glint integrator based on normal maps. Our method leverages a local energy preserving microfacet-based BRDF that fakes normal perturbations. Our patch-BRDF can be used in Monte Carlo forward path tracing using multiple importance sampling. A large collection of material appearances can be reproduced through the use of different texture maps. Compared to the previous methods based on classic normal mapping, we obtain renderings that are free of artifacts such as black holes or a too dark overall appearance.

Currently available glint integrators only model opaque surfaces. To our knowledge, no method addresses the problem of rendering rough or scratched glass. Transparent materials are more sensitive to energy loss. We hope that our new approach will make solving this problem more accessible.

Figure 2 :

 2 Figure 2: Comparison between a smooth Beckmann surface (left), a microfacet-based normal mapped surface using tangent normals (middle, [SHHD17]) and using symmetric normals (right, our micro-surface). The normal map being very dense in this Figure, tangent facets introduce unwanted roughness compared to the smooth model while our micro-surface does not. Even with 262,144 spp in this scene, there still are fireflies as random walks are used to model multiple-scattering.

Figure 3 :

 3 Figure3: Left: classic normal mapping. Middle: normal mapping of Schüssler et al.[START_REF] Schüssler | Microfacet-based normal mapping for robust monte carlo path tracing[END_REF]. Right: our model. The symmetric facets ω s do not change overall material roughness the way tangent facets ω t do. For glint rendering, where the normal map is highly repeated, tangent facets would alter too much the appearance.

Figure 5 :

 5 Figure 5: Different cases of local multiple-scattering BRDF sampling.

Figure 6 :

 6 Figure 6: Discrete set of Gaussian weighted BRDFs, and associated discrete set of normals ω p (u j ).

Figure 7 :

 7 Figure 7: Let D(ω p , α) be a Gaussian normal distribution with isotropic roughness α over the normal map space ω p . D approximates the accurate normal distribution within a ray footprint. We use it to precompute footprintdependent first-bounce energies Ω E 1 (θ o , ω p )D(ω p , α)dω. The values are stored in a 2D array, mapping energy quantities to (α, θ o ) value pairs.

Figure 8 :

 8 Figure8: Classic normal-mapped surface. Importance sampling visible normals can only be done if all normals in the footprint are examinated one by one. Our method does not have to go through this.

Figure 9 :

 9 Figure 9: Comparison between classic normal mapping (left) and either our single (middle) or multiple (right) scattering microfacet-based normal mapped surface. Our multiple-scattering BRDF f ∞ is 100% energy conserving and removes the black holes (top rows) and the overall dark appearance (bottom rows) on the surface.In the bottom rows, the normal map repetition is twentyfold, leading to a coarse Beckman specular lobe. We use 131,072 spp to render the majority of glints and achieve noise-free renderings.

  rendered with a fixed micro-BRDF roughness α = 0.02 (flat elements variant) requires about 90 s with the method of Yan et al. and Chermain et al. Decreasing the value of h does of course affect performance, but to a lesser extent.

Figure 10 :Figure 11 :Figure 15 :

 101115 Figure 10: Cutlery scene rendered at 1024 spp and exhibiting a wide range of materials that our method can handle. From left to the right: brushed and scratched metals, green coated rough plastic, machined metal, rough metal and finally purple coated glittery plastic.