Marina Sergeeva

Daniel Delahaye

Leïla Zerrouki

DYNAMIC AIRSPACE CONFIGURATIONS GENERATED BY EVOLUTIONARY ALGORITHMS

This paper focuses on the process of generating a sequence of sector configurations composed of two airspace component types Sharable Airspace Modules (SAMs) and Sectors Building Blocks (SBBs). An algorithm has been developed that manages the main features of the dynamic sectors configuration (including sector design criteria). In order to make it run efficiently a pre-processing step will be presented to create a graph modelling of the inputs. Based on this initial graph, a mathematical model is defined which can be summarized by a multi-periods geometric graph partitioning problem. State, space, objective function and constraints will be also presented. Due to the induced complexity, a stochastic optimization algorithm based on artificial evolution is then proposed. A two layer chromosome is used for such a genetic algorithm for which recombination operators are proposed. Evaluation of the algorithm will be presented with a comparison to existing tools and operational approach.

Introduction

This paper presents an ongoing research work performed in the scope of SESAR program to assess the concept of Dynamic Airspace Configuration (DAC) in highly complex free routing areas [START_REF]5.4 project, Dynamic Airspace Configuration[END_REF].

The first section of the paper presents an overview of past work related to dynamic sectorisation and dynamic airspace configurations. It helps to distinguish new elements of the approach proposed in this paper.

In order to set-up the operational context and identify the operational objectives and constraints, the second section provides a high level description of the DAC concept in the framework of SESAR program.

The technical solution presented in the third and fourth sections consists in a graph modelling of the inputs and a mathematical model which can be summarized by a multi-periods geometric graph partitioning problem. Due to the induced complexity, a stochastic optimization algorithm based on artificial evolution is then proposed.

The resulting airspace configurations are assessed in section V, by comparing the solution to currently used tools and operational methods. Section VI concludes the paper and gives some perspectives of future work within SESAR program and its potential applications.

Previous Related works

A number of researchers, mainly from Europe and USA have achieved certain results in the development of the solution of the DAC problem [START_REF] Flener | Automatic airspace sectorisation: A survey[END_REF]. Currently, airspace is partitioned into functional blocks or pre-defined sectors that may be combined into bigger sectors when the traffic is low. Existing support tools enable the computation for the flow management planning phase of the sector opening schemes, which minimize sector overloads based on pre-defined sector configurations. The ICO tool developed at EUROCONTOL is one of them; cf. [START_REF] Vehlac | Improved Configuration Optimizer[END_REF] and [START_REF]NEST[END_REF]. More in a research perspective, a methodology for dynamic sectorisation for highly complex Area has been developed [START_REF] Ehrmanntraut | Airspace design process for dynamic sectorisation, 26th DASC[END_REF]. In preparatory work of SESAR program, "Dynamically Shaped Sectors" early project proposed a heuristic for automated sector design which combines best practice of sectors design and clustering techniques [START_REF] Zerrouki | Dynamically Shaped Sectors[END_REF].

In some of the existing works the DAC problem was solved using approaches that are based on the initial division of the airspace into small 2D or 3D areas [START_REF] Drew | A method of optimally combining sectors[END_REF]- [START_REF] Yangzhou | Dynamic airspace sectorization via improved genetic algorithm[END_REF], such as elementary sectors for example. Other researches are concentrated on redesigning the airspace boundaries. Stochastic local search, Mathematical programming [START_REF] Drew | A method of optimally combining sectors[END_REF], [START_REF] Gianazza | Forecasting workload and airspace configuration with neural networks and tree search methods[END_REF], Computational geometry [START_REF]NEST[END_REF], Evolutionary algorithms [START_REF] Yangzhou | Dynamic airspace sectorization via improved genetic algorithm[END_REF], [START_REF] Delahaye | Genetic algorithms for partitioning airspace[END_REF]- [START_REF] Delahaye | Airspace sectoring by evolutionary computation[END_REF],

Global optimization and several other techniques have been used in prior research on airspace optimization. Several different complexity metrics have been developed [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management[END_REF], [START_REF] Kopardekar | Airspace complexity measurement: An air traffic control simulation analysis[END_REF], and applied [START_REF] Klein | Simplified dynamic density: A metric for dynamic airspace configuration and nextgen analysis[END_REF] in those works in order to evaluate resulting sector configurations.

In [START_REF] Yangzhou | Dynamic airspace sectorization via improved genetic algorithm[END_REF], the DAC problem is solved by using an improved genetic algorithm (iGA). On the first step of the algorithm a graph model is constructed that represents the airspace static structure. Then, the algorithm solves a graph-partitioning problem. It aims to balance the sector workload under the constraint of ensuring safety. The iGA includes multiple populations and hybrid coding. The algorithm aims to find balanced airspace configurations with an optimal number of sectors. New criteria were included in the algorithm: sector balancing, coordination workload and sector average flight time. The algorithm includes two steps. In the first step, the number of sectors is determined and airspace partition is accomplished. In the second step, sector boundary optimization is implemented. Both of the two steps are done using a hybrid encoding multi objective GA on the basis of the weighted graph model. The method has been proven and is being applied to the airspace of North China.

Authors of [START_REF] Gianazza | Forecasting workload and airspace configuration with neural networks and tree search methods[END_REF] have proposed to use relevant air traffic complexity metrics instead of flight counts and sector capacities. In this work possible combinations of elementary sectors are explored, instead of the small subset of pre-defined configurations currently being used, so as to offer the maximum capacity to the incoming traffic. Artificial neural networks are used to predict the sector status from the relevant complexity metrics. In this work the trained network [START_REF] Gianazza | Selection and evaluation of air traffic complexity metrics[END_REF] and the most relevant metrics are used to predict when the controlled sector will be normally loaded, overloaded, or under-loaded, and when an airspace configuration should be recombined. The decision to reconfigure sectors is driven by the prediction made by the neural network. A classical tree search algorithm is used to build all the valid sector configurations from an initial set of elementary sectors. The tree search algorithm explores all possible airspace configurations, among which only one is chosen using defined evaluation criterion. The computed configurations were compared to the actual configurations archived by ATC centers.

The next two methods are more freeform and have not used any elements of the current functional airspace block design.

In [START_REF] Martinez | A weighted graph approach for dynamic airspace configuration[END_REF] authors have described a method for partitioning airspace into smaller regions based on a peak traffic-counts metric. The method was used for partitioning the airspace over the United States into 466 sectors, once for each hour of the 24-hour day. The algorithm consists of two phases. In the setup phase, a weighted graph is created based on traffic flow patterns. Topological structure of the airspace (pathways of a flow pattern) can be described as a graph. Nodes of the graph are represented as an averaged position coordinates of flights that enter or leave a given sector through a boundary and links between them represent the flow. The flow is computed as a number of aircraft along a link of the graph, and is assigned as its weight. The weight of the node is computed as a number of aircraft located in the area close to it. The final phase consists of three steps that are repeated in a cycle. During the first step one spectral bisection method splits a given graph into two sub-graphs such that the node-link connectivity is preserved in each sub-graph. On the next step, graphs weight peaks are computed. The sub-graph with the highest weight is selected for further partitioning. The process is repeated until the termination criterion is met.

Authors of [START_REF] Li | A spectral clustering based algorithm for dynamic airspace configurations[END_REF] have proposed an airspace modeling technique based on the air route structure and a spectral clustering based dynamic airspace configuration algorithm to partition the graph model of the airspace into subgraphs. The performance of the algorithm was validated with real air traffic data. In this method, first, a weighted graph model using the airspace information and air traffic data is built. Then, a novel airspace sectorisation algorithm is applied. The algorithm is based on a graph search algorithm. Each designed sector contains one subgraph, and has a smooth polygonal boundary. Sector boundaries created by this algorithm should satisfy the minimum distance constraint, which means that the main flow paths and intersections should be at a minimum distance away from the sector boundaries. The authors have as well proposed various evaluation metrics of resulting sectorisation. The simulation results have shown that the new sectorisation computed by this algorithm outperforms the currently used sectorisation.

Most of the previously developed approaches do not take into account the stability of the generated configurations over time and some important airspace design aspects. As a matter of fact, one must be able to generate airspace configurations that have minimum changes from one time period to another with reduced airspace complexity, and more balanced workload among airspace sectors.

DAC operational concept in SESAR

In comparison with a sectorisation based on a fixed route network, the SESAR program introduces the User Preferred Routing/free routing concept to enable the airspace users to freely plan and fly the 4-D trajectory that best copes with their business or mission needs. Contrary to a fixed-route network a FRA environment will produce a much larger number of different trajectories for which the dynamic nature and flexibility of the DAC process will cater most efficiently.

Increasing adaptability of Sectors Configurations

The DAC concept in SESAR proposes increasing levels of adaptability of the sectors configuration to match with the capability of each Air Traffic Control Center regarding types of free routes implementation, ATC licensing, and airspace delegation agreements.

Four levels of configurations adaptability have been identified: DAC Level-1 is the lowest level of airspace adaptability, which provides an improvement to the current operational process at minimum cost. Instead of using predefined sectors configuration, Level-1 uses dynamically computed configurations by combining existing elementary sectors provided by an automated optimization process. In level-1 any sector combination which forms controllable airspace blocks is eligible and may be used during the day of operation. DAC level-2 increases the adaptability of the airspace to the traffic pattern by delineating from the nominal sectors, a limited number of airspace components that can be recombined laterally and/or vertically within an automated and optimized sector configuration process to form controllable airspace blocks. Two airspace structures are specified in level-2: the Sector Building Blocks (SBB), are permanently busy areas with a high traffic load delineated by recurring traffic patterns, while Sharable Airspace Modules (SAM) are less busy with a temporary high traffic load that can be used in the sector configuration computation process to equally balance the traffic load among the sectors. SBBs in level-2 can be sufficiently large to be workable and controlled separately. SAMs cannot be used separately in a configuration. The number of the elementary sectors resulting from the combination of the airspace components is known and limited. However, the number of possible configurations will be of a greater number than in level-1, which allows a better adaptation to changes of the traffic pattern.

Figure 2. DAC level-2 SBBs and SAMs components (lateral view)

DAC level-3 combines a larger number of SBBs and SAMs. The SBBs are defined by a less structured traffic pattern which makes them smaller and their workability as a separate entity not ensured.

However, the SBB component is preserved as a reference for controllers training and licensing purposes.

Figure 3. DAC level-3 SBBs and SAMs components (top view)

DAC Level-4 combines only SAMs which are not related to any particular traffic pattern. It is the most challenging level since there will be no more predefined airspace component known and full generic ATC validation is a prerequisite for the implementation of this level.

DAC processes overview

Each DAC level is associated with an airspace design approach which produces the airspace components that will be used in the airspace configuration process. Elementary sectors, SBBs and SAMs are the result of a sector design process dedicated to level 1, 2 and 3, while DAC level 4 is itself a sector design process performed dynamically.

Airspace design for level-1 is similar to the current process which produces elementary sectors; however supporting automated tools are proposed at the aim of increasing the efficiency of sector design process [START_REF] Sergeeva | 3D Airspace Sector Design by Genetic Algorithm[END_REF]. Sector design for level-2 and 3 which delineates SBBs and SAMs is a subject of future research.

Because of the huge number of possible sectors configurations enabled by more modularity of the sector design phase, the definition of sector configurations and opening schemes will require the use of automated support tools.

Given the organizational framework of the concerned airspace, the challenge is to organize, plan and manage airspace configurations to meet User Preferred Routing, in a Free Route environment with enough flexibility to respond to changes in traffic demand, to unexpected events including weather, and to any update in airspace reservation within the Advanced Flexible Use of Airspace concept (AFUA) in the optimum way, while maintaining the safety targets.

The Optimum Configuration process will be able to build Configurations of sectors based on the grouping/de-grouping of elementary sectors and reallocation of Sharable Airspace Modules. The process will at least provide either a configuration with the minimum number of ATC working positions without imbalance or, a configuration with the minimum of imbalance considering the number of ATC working positions available.

It is obvious that workload assessment is a key requirement for generating workable sector configurations and more so in a Free Route User preferred trajectory environment where there might be very little historical data available and due to its very nature no possibility of managing trajectories by using published routings.

The workload assessment in FRA must be achieved through evolution/development of tools. These tools must support load assessment in generic airspace defined and managed through SBBs and SAMs. Workload assessment in a free route environment is only possible based on a trajectory that is well described and shared by the NOP.

Figure 5 provides an overview of the sector design and sector configuration processes described above within the overall context of airspace organization and management in SESAR.

DAC models and algorithms

Airspace modelling

The first stage in the optimization process consists of modeling the real problem using a mathematical abstraction which should be as faithful as possible. The modeling stage then consists of characterizing the state space and the objective space.

The state space represents the set of parameters of the system upon which we may act in order to optimize one (or more) objective(s). Examination of the properties of the state space then assists us in the choice of a suitable optimization method. In most industrial optimization problems, variables of the state space must remain within a sub-domain defined by a set of constraints. One very important point that is characterizing the state space is its dimension. Generally, the higher the dimension n of the state space Χ is, the harder it is to find the optimum.

The objective space represents the set of criteria which we wish to optimize. Within the context of industrial optimization problems, we seek to determine global optima and to avoid being trapped on local optima.

First, we need to model the airspace. Let's consider a given airspace which consists of a known set of 3D airspace blocks (Figure 6).

Figure 6. Initial airspace blocks (2D projection). The green blocks are SBBs and the pink ones are SAMs

The airspace area is split into several altitude layers. The shape of those blocks can vary depending on the altitude level. A set of aircraft trajectories that are crossing this area is known as well. For the purpose of DAC, traffic data is taken just for one day. There are two type of blocks used in our model: sharable (SAMs) and non-sharable (SBBs). Each controlled sector is built of at least one SBB. Building of the controlled sector starts from choosing a central block, which can be chosen only among non-sharable blocks. Even if the number of controlled sectors is different in two successive configurations, centers of controlled sectors will be chosen among the same group of non-sharable blocks. This partly ensures continuity between successive airspace configurations.

State space modelling

Our model is based on a weighted graph model of the airspace (Figure 7). Let a graph G = (N, L) represent the airspace, where N is the set of nodes and L is the set of links.

Figure 7. Initial graph. Green nodes represent non-sharable blocks (SBBs)

In this graph each node represents a sharable or non-sharable block and each link represents the relation is neighbor with between two blocks, it means that when two blocks share a common vertical or horizontal border (Figure 7), a link is built between them. Weight of the node represents a workload and weight of the link represents the flow between two sectors.

As for the evaluation metric, several of them have been included in the model. In order to estimate the complexity of the traffic in each block a convergence metric [START_REF] Gianazza | Selection and evaluation of air traffic complexity metrics[END_REF] is used and for the evaluation of the traffic load -an occupancy count. Occupancy count is computed as a number of aircraft inside the block at each minute. Link's weight of the graph is computed as a number of aircraft crossing the border between two airspace blocks connected by this link. The workload of airspace blocks is computed for each chosen time period, as well as transfer traffic between each pair of neighboring blocks (Figure 8). In order to be able to compute the total number of overloads in configuration, for each minute of chosen day we are computing the number of aircraft in each block.

Figure 8. Graph time extension

In order to be able to compute re-entry events inside created controlled sectors, we register the list of airspace blocks crossed by each trajectory with the associated time horizon (Figure 9).

A graph partitioning problem

Based on the airspace model described above, the DAC problem consists in finding an optimal partitioning of the graph into several sub-graphs for each time period. The regular graph partition problem is defined on data represented in the form of a graph G = (N, L), where N is the set of nodes and L the set of links, such that it is possible to partition G into smaller components with specific properties. Having a connectivity constraint on elementary sectors belonging to the same sector, such graph partitioning is a geometric graph partitioning. In this case, nodes belonging to the same component have to be connected. The Figure 10 gives an example of partition of a graph into three connected components. As it can be seen on the figure, for each pair of nodes belonging to the same component, there is a path connected them for which all the elements of such path (links and nodes) belong to this component. If we note N k the set of nodes belonging to the component k, the connectivity constraint can be modelled by: pair of node (i, j) N k a path P = {l 1 , l 2 ,…l d } connecting node i and node j for which each link l i has its origin and destination nodes in N k .

Figure 10. Example of a partition of a graph into three connected components

For a given time period, the state space will be modelled the following way: X = {N 1 , N 2 ,…,N K } where N i represents the set of nodes belonging to the component i. Based on those sets of nodes, it is easy to determine the links shared by two components (black dash line of the Figure 10). Having a problem with several time periods, such set of nodes are computed for each time period:

where K i represents the number of component at time period t i .

Based on the state space definition, we must model the associated objective function.

Objective function

The objective function as well as the constraints included in our model is designed according to requirements developed in co-operation with operational experts from EUROCONTROL. Five criteria have been included in our objective function for evaluation of a solution.

The first criterion measures the level of the workload imbalance among all controlled sectors in configuration for each time period t i (i = 1..N). The workload of the controlled sector is computed as the sum of workload of airspace blocks that it is composed of. The workload imbalance of all components (controlled sectors) for the time period t i is then computed with the following formula:

where K i is the number of components in the configuration for period t i , Wc ik is a total weight of all nodes (blocks workload) included in the component k (a workload of the sector) for period i, C is the sector capacity (a value of the workload that we are trying to obtain for each controlled sector).

Then, the total imbalance associated to the whole planning is given by: where T is a number of periods.

The second criterion measures the total number of overloads in each controlled sector of the configuration. For each one minute of the period we are computing the number of aircraft in the controlled sector, and if this number exceeds a given value during several minutes consequently, we register an overload. Then we summarize the number of overloads of all controlled sectors for each time period:

where v i is the total number of overloads for period t i and vcs i is the number of overloads in one controlled sector.

Then, the total number of overloads associated to the whole planning is given by: The third criterion measures the transfer traffic between neighboring components. When two neighboring blocks belong to different controlled sectors, the traffic flow between them is cut by the sector border, this increases the coordination workload. The total flow cut for the time period t i is given by: where -the flow cut between nodes i and j, in both directions, for the time period t.

The total flow cut for a given planning is then computed with the following formula:

The fourth criterion counts the number of reentry events. Using the list of nodes (crossed airspace blocks) for each trajectory and the set of nodes of each component, it is possible to compute the number of re-entry events for each time period, which be noted N r .

The last criterion measures the number of short transits inside each sector (N s ). It is relatively easy to compute efficiently the number of aircraft that are staying too short time inside one sector using trajectories that are represented by its associated list of blocks with entering and exit times.

All those criteria are then aggregated into a single objective function:

Application of EAs to the sector design Problem Evaluation Algorithms

Typically, graph partitioning problem falls under the category of NP-hard problem and stochastic optimization is a good candidate to address it. For this type of problem there could not be only one solution, there may be several optimal (or near optimal) solutions and we should be able to find most of them, as long as they have to be later evaluated and refined by airspace experts. This last point makes us reject classical simulated re-annealing optimization which updates only one state variable. On the other hand, Evolutionary Algorithms (EAs) maintain and improve a population of numerous state variables according to their fitness and will be able to find several solutions. Thus, EAs appear to be relevant to solve the sector design problem.

Coding the chromosome

In the previous section we have proposed a way of modeling the airspace configuration as a set of connected components for each time period. The coding used for this problem consists in representing connected components by sub-sets of nodes.

EAs use techniques inspired by evolutionary biology such as inheritance, mutation, natural selection, and recombination (or crossover) to find an approximate solution of the optimization problem.

An individual (a solution of the problem), is represented by a list of parameters, called chromosome. The chromosome used in our algorithm is presented in Figure 11.

Figure 11. Chromosome structure

The chromosome consists of two parts. First part includes a table with permutation of root nodes (their indexes) and a table that contains temporal segments. For each time period this second table will contain the list of root nodes, chosen from nodes that represent non-sharable blocks, and that are used in construction of the controlled sectors. Root nodes are used as starting nodes for graph partitioning (see Figure 12). Second part contains connected components, coded as lists of nodes belonging to same component. The method ensures that nodes of the same sub-set are connected by at least one path.

For each time period several nodes are selected from the list of non-sharable nodes as shown in Figure 11. Those nodes will be considered as root nodes. The minimum number of root nodes that can be selected equals to 1 and the maximum equals to the number of all non-sharable blocks i.e., to the number of maximum controlled sectors in the configuration. Here, nodes 1, 4, 8 and 9 are nonsharable nodes and potential centers of controlled sectors. Three time periods are considered and three temporal time segments are generated first randomly. At the time period 1 there is one sub-set created with root node 8, at the second period -3, with root nodes 8, 9 and 1, and etc. (see Figure 11).

Figure 12. Example of the coding graph partitioning for one time period

The initial permutation of root nodes ensures a random mapping between temporal segments and root nodes (this avoids the same root node to be associated with the first temporal segment). The way of coding chromosomes with temporal segments ensures the stability in time of controlled sectors shapes. As a matter of fact, for successive time periods, similar list of root nodes will be used as a sector centers, ensuring those volume of airspace being controlled by the same controller.

Chromosome Initialization

In order to create an initial population of solutions, we should be able to generate chromosomes randomly. First, the root node order is selected by the mean of the initial permutation. Then, temporal segments are randomly built. For each time period several root nodes are chosen from the permutation table. The number of chosen root nodes per time period is first generated randomly and is equal to the number of controlled sectors in the configuration. The first root node in the permutation table participates in the partitioning process for each time period. Finally, based on those chosen root nodes a greedy heuristic is used to build the initial partitioning which is summarized in Figure 13.

Figure 13. Greedy heuristic is used to create initial partitions

In this example, there are two root nodes, marked with labels A and B (step 1). On a first step, all nodes connected with one of the root nodes are marked respectively, with the label A or B (step 2 and 3). Then, the algorithm makes several iterations, running each time through the permutation list of all other nodes, until every node receives its label. At each iteration, the algorithm checks all neighbors of a non-labeled node and if at least one of its neighbors is labeled, the node receives the same label (step 4). If several neighbors of the node have different labels, then the label that represents the lightest component (with the smallest workload) is chosen. If the node does not receive a label, it is left without it for the next iteration.

Recombination operators

After creating the first population, each solution is evaluated, and a value of fitness is returned by a fitness function. This initial population undergoes a selection process that identifies best solutions. Then, three following recombination operators are applied to individuals: nothing, crossover, or mutation. The associated probability of application are respectively (1-p c -p m )), p c and p m .

Crossover between two parent solutions results into two new child chromosomes, which are added to the next population. During the application of crossover operator, chromosomes of parents are mixed or exchange between each other. The purpose of mutation in EAs is to allow the algorithm to avoid local minimum by preventing the population of chromosomes from becoming too similar to each other, which could cause slowing or even stopping evolution. These processes ultimately result in the next population of chromosomes POP(k+1) that is different from the initial generation. This generational process is repeated in cycle N G times, where N G is a required number of generations.

Several recombination operators have been developed in our algorithm.

First type of operators has been designed to modify the first part of the chromosome. In order to change initial root nodes, a mutation operator has been developed. It consists in changing the order of root nodes by randomly exchanging two nodes in the permutation table. The mutation operator for temporal segments changes the number of segments by adding or removing one segment i.e. adding or removing one controlled sector into configuration. One slicing crossover operator has been developed for temporal segments as well. Two parents with time intervals are used to create two children by randomly selecting one gene position and by exchanging all genes from such position.

The second type of operators for the second part of the chromosome consists of two mutation operators. Before applying any mutation operator a time period is selected either randomly, either according to the associated graph partitioning performances, meaning that configuration with a worst objective function will have more chances to be selected. Then one of the graph partitioning mutation operators is applied.

First operator modifies permutation table of all nodes (sharable and non-sharable). In order to change a resulting graph partitioning without changing root nodes order we are simply changing the order of all other nodes. This ensures that the graph partitioning process will not start from the same nodes, and so resulting components will probably consist of different nodes.

Then, the second operator begins by selecting the "worst" component. After that, it starts to search the suitable neighbouring component. In case if the chosen component is not loaded enough, then it searches for the neighbour with the highest workload (the same for the flow cut). On the other hand, if the chosen component is overloaded, then it searches for the neighbour with the lowest workload. This step is carried out statistically (introducing a bias into a random selection). Then we obtain one link which connects these two components. One of the nodes belonging to this link is then moved from one component to the other one (depending if the component is overloaded or not), while verifying that the component which is losing a node remains connected.

Having described our algorithm, the next section presents several examples of results.

Simulation and results

This section presents the results of dynamic sector configurations algorithm applied to Maastricht Amsterdam FIR and free routes traffic.

Scenarios description

Three scenarios have been used to test the algorithm. Each scenario is based on free route trajectories simulated by Lido 1 and Sabre 2 systems, which provide a sample of 2332 full free route trajectories for the 14 th of July 2014 crossing Maastricht/Amsterdam Airspace (EDYYDUTA). The baseline scenario to which the algorithm will be compared characterizes current operational methods and tools. It uses a limited number of predefined sectors configurations, based on known 1 Lido: Lufthansa Systems flights planning tool 2 Sabre: Sabre Airlines Solutions, flight planning tool traffic pattern and current organizational framework, including ATC licensing constraints. The tool used to compute the optimum sectors configuration sequence over a traffic day (the opening scheme) is ICO, Improved Configuration Optimizer. ICO algorithm is implemented within NEST simulator (cf. [START_REF] Vehlac | Improved Configuration Optimizer[END_REF] and [START_REF]NEST[END_REF]), the Network Strategic Tool developed by EUROCONTROL.

Two solution scenarios have been tested. The first is a DAC level-1 case, which uses 6 initial elementary sectors as components inputs to the algorithm, cf. figure 16. The sectors traffic load is monitored through occupancy count metric. It is assumed that the capacity used for all the sectors for the baseline as well as solution scenarios is a constant value. It corresponds to a maximum number of flights per minute over a sustained period of time (e.g.: maximum of 8 flights per minute over a sustained period of 12 minutes). It is also assumed that the backbone (the number of available controllers) is a constant value of 8 controllers all over the day. Capacity and backbone values are used for technical testing purpose and do not reflect real operational values.

Figure 21 and 22 presents results about the quality of the computed configurations; i.e. criteria of workload balancing, number of re-entries and shortcrossing flights that are taken into account by the algorithm for the selection of the most acceptable configurations in the operational view point.

Figure 20 below shows that the proposed algorithm used with DAC level-2 scenario outperforms the baseline as for overload and number of controllers working hours. DAC level-2 decreases significantly the overload of about 15% compared to the baseline. While DAC level-1 scenario provides better results regarding the number of controllers hours. Nevertheless, the operational acceptance of the proposed configurations is not ensured. Further analysis of each new configuration and also the changes in the configuration sequence must be properly analyzed. An initial expert judgment is provided in the next section.

Initial analysis of the new configurations

The results presented previously are promising; they demonstrate that the algorithm performs in adequacy with the specified operational objectives and constraints. Nevertheless, further operational evaluation is required to assess the workability of the proposed configurations and their sequence (opening scheme). In the scope of this paper, an initial judgmental assessment is provided for three new configurations proposed in scenario DAC level-2.

Figure 23 presents sectors configuration 4AXE. This configuration is built by grouping all the sectors of the lower level to form an extended lower sector displayed in light-green. This configuration is not operationally used because the current system visualization imposes limit to the size of the sector to be controlled. The maximum length of a leg in a combined sector should be 250NM. This technical constraint is not taken into account by the algorithm. Figure 24 presents one often used configuration, "4AXC", in DAC-level 2 scenario. The lower sector displayed in light-blue color, is composed of two grouped sectors (SBBs) which currently require two different ATC licenses, which explains why it is not used in today's configurations. However, the algorithm considering a unique license for all the sectors proves that it can offer better performance regarding sectors load balancing and overload. In addition, due to the use of SAMs components, can be added to this combination of SBBs, different combinations of SAMs, which leads to new sectors configurations, with a slight modification, as it is presented in configuration '4AXB' displayed in figure 25. The difference between configurations 4AXC and 4AXB is one SAM component that is allocated to the upper sector displayed in pink color in figure 25.

Configurations 4AXC and 4AXB are used consecutively in the opening scheme proposed by the algorithm to better balance the workload, and there is no foreseen operational opposition to the use of either 4AXC and 4AXB configurations or sequence, providing that a unique ATC licensing is applied to their components.

Conclusion and next steps

The modelling approach and algorithm solution presented in this paper have been tested and compared to existing technical solution. The provided results demonstrate that the new solution fits with the requirements of the DAC concept and outperforms the baseline solution regarding the minimization of the overload, sector load balancing, as well as reentering and short-crossing flights, at least for DAC level-1 and level-2 cases. Initial tests have been performed for scenarios DAC level-3. The provided results show that the algorithm is not adapted to that case, since an increasing number of SBBs and SAMs, requires geometrical constraints to control the shape of the sectors, in terms of convexity, balconies, etc. Further investigations will be achieved in order to define a solution for DAC level-3 sectors configuration.

Figure 1 .

 1 Figure 1. DAC level-1 Elementary Sectors components (top/lateral view)

Figure 4 .

 4 Figure 4. DAC level-4 SAMs components (top/lateral view)

Figure 5 .

 5 Figure 5. DAC processes overview

Figure 9 .

 9 Figure 9. List of airspace blocks associated to a given trajectory

Figure 15

 15 Figure 15 representing free route trajectories, compared to the trajectories planned in today's Routes Network in Figure 14, shows a noticeable dispersion of the trajectories and a loss of traffic pattern.

Figure 14 .Figure 15 .

 1415 Figure 14. Planned trajectories in fixed-route network

Figure 16 .

 16 Figure 16. DAC level-1: 6 sectors scenario The second scenario is a DAC level-2 case. It uses 6 SAMs components obtained by cutting vertically the 6 original sectors at altitude levels 335 and 355 in order to increase flexibility of new sectors configurations built by the algorithm, cf. figures 17 and 18.

Figure 17 .Figure 18 .

 1718 Figure 17. DAC level-2: 6 sectors and 6 SBBs

Figure 20 .

 20 Figure 20. Overload and controllers hours Figure 21 shows a significant improvement of the quality of the configurations provided by the solution scenarios in terms of workload balancing compared to the baseline. The mean difference of the sectors load within a configuration of the solution scenario is about 20% while it reaches 33% for the baseline configurations. Those results demonstrate that the proposed configurations, which freely combines elementary sectors (DAC level-1) and 6 SBBs and 6 SAMs enable to propose sectors configuration with more balanced sectors load.

Figure 21 .

 21 Figure 21. Sector load balancing criteria Figure 22 presents the average number of reentering flights and short-crossing flights per configuration for the baseline and solutions scenarios.It clearly shows that the algorithm optimizes both criteria. Similarly to overload criteria, DAC level-2 is the scenario which also best optimizes re-entries and short-crossing criteria due to the increased number of airspace component combinations which allows the algorithm to find a greater number of solutions fitting the traffic structure.

Figure 22 .

 22 Figure 22. Re-entering and short-crossing flights

Figure 23 .

 23 Figure 23. DAC Level-2 Conf 4AXE

Figure 24 .

 24 Figure 24. DAC Level-2 Conf 4AXC

Figure 25 .

 25 Figure 25. DAC Level-2 Conf 4AXB

Some Quantitative Results

The algorithm has been executed for both solutions scenario and ICO for the baseline scenario. Figure 19 displays the Opening Schemes computed for each scenario and arranged in time tables.

The values in Y-axis are the configuration name, where the first digit informs about the number of sectors used in the configuration. It shows that the algorithm builds configurations composed of a number of sectors which increases and decreases with the traffic profile of the day. The maximum number of sectors in a configuration is 4, which is in line with the limit imposed by the number of available controllers [START_REF] Sabhnani | Flow conforming operational airspace sector design[END_REF]. The same is for ICO algorithm in the baseline scenario. It also shows that the number of configurations and configuration changes increase significantly for DAC level-2 scenario. For example, it is proposed 5 different configurations composed of 4 sectors and 3 configurations of 3 sectors and 19 configuration changes over the day.

The performance of the opening schemes proposed for the baseline and solutions scenarios are displayed in the next three figures. Figure 20 provides global results regarding the overload and the number of controllers hours generated by the opening scheme of each scenario.

Figure 19. Opening Scheme for baseline and solution scenarios

Baseline Opening Scheme: ICO using pre-defined configurations DAC level-1 scenario: solution algorithm using elementary sectors DAC level-2 scenario: solution algorithm using 6 SBBs and 6 SAMs


Arnaud Dessein

William A P Smith

Richard C Wilson

Edwin R Hancock

SEAMLESS TEXTURE STITCHING ON A 3D MESH BY POISSON BLENDING IN PATCHES

Keywords: 3D model, texture reconstruction, seamless stitching, Poisson blending, uniform overlapping patches

In this paper, we propose a novel approach to seamless texture stitching on a 3D mesh. The main idea is to blend the sampled images by least-angle selection of gradients in overlapping patches. This is in contrast to previous works which focus on vertex-or face-based strategies with additional heuristics for robustness. Patches are obtained by growing a uniform mesh segmentation via geodesic projections. Blending is achieved by formulating a screened Poisson equation using discrete differential operators. The processing pipeline further includes an optional step of color transformation for calibration and correction. This is applied to a zippered mesh based on range scans and a morphable model fitted to face photographs.

INTRODUCTION

The problem of stitching several textures to form a seamless texture arises from various applications in 2D image and 3D mesh processing. In particular, this problem is inherently tied to the process of image-based object modeling to capture real-world objects. Standard techniques in this context use an RGB-D scanner or camera to collect multiple view-dependent observations consisting of color and depth maps, with technologies including time-of-flight or structured light for range images and binocular stereo for intensity images. It is relatively easy to recover the 3D shape of the object from these depth maps, using classical methods such as the iterative closest point algorithm [START_REF] Besl | A method for registration of 3-D shapes[END_REF] to align observations, and zippering [START_REF] Turk | Zippered polygon meshes from range images[END_REF] or volumetric merging [START_REF] Curless | A volumetric method for building complex models from range images[END_REF] to integrate them on a base shape mesh. However, it is much more difficult to rebind color information to recover texture.

A related problem is shape-from-X, where the 3D shape is estimated from 2D images via cues such as shading or texture. Common methods include bundle adjustment [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] for image registration, and silhouette-based visual hull [START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF] or space carving [START_REF] Kutulakos | A theory of shape by space carving[END_REF] for shape recovery. Another example is when a statistical model for a specific object class, such as 3D morphable models for human faces [START_REF] Blanz | A morphable model for the synthesis of 3D faces[END_REF], is fitted to one or more 2D images. In these scenarios, the recovered shape must be textured from the available observations of the different scene points. Although established pipelines exist for reconstructing shape, texture stitching is, surprisingly, much less studied.

In this paper, we assume that geometry has already been processed and focus on texture. There is in general no exact correspondence between points on the images and on the recovered surface. Thus, one needs both to interpolate color and to integrate concurrent data from several images. For instance, a baseline method is to sample the images on the mesh by projection and color interpolation, and then merge the textures from overlapping images by averaging. However, such a back-projection of color is likely to contain seams. The main issues are that (i) the reconstructed geometry is only approximate, (ii) the views cannot be registered perfectly, (iii) the acquisition settings may vary across views, and (iv) the lighting may change between poses. Hence, texture reconstruction for 3D models is still a challenging problem. Several approaches have been developed to cope with these issues and obtain a seamless stitching.

Related work

Early works have focused on different weighting heuristics to average overlapping textures. Pulli et al. [START_REF] Pulli | Acquisition and visualization of colored 3D objects[END_REF] couple viewing angle and proximity to the projected outline. Neugebauer and Klein [START_REF] Neugebauer | Texturing 3D models of real world objects from multiple unregistered photographic views[END_REF] additionally consider texture deviation from a normal model. Matsushita and Kaneko [START_REF] Matsushita | Efficient and handy texture mapping on 3D surfaces[END_REF] or Johnson and Kang [START_REF] Johnson | Registration and integration of textured 3D data[END_REF] exploit projected face areas. Levoy et al. [START_REF] Levoy | The digital Michelangelo project: 3D scanning of large statues[END_REF] merge several cues: projected area, lighting angle, proximity to the projected outline, to the image border, and to the silhouette edge seen from the light source or the mirror direction, in a confidence index which is then smoothed conservatively among adjacent vertices. Bernardini et al. [START_REF] Bernardini | High-quality texture reconstruction from multiple scans[END_REF] also combine different factors: viewing angle, distance to an edge, deviation from a normal model and resolution. Wang et al. [START_REF] Wang | Optimal texture map reconstruction from multiple views[END_REF] rather formulate an optimal weighting based on sampling theory. Similarly, Baumberg [START_REF] Baumberg | Blending images for texturing 3D models[END_REF] proposes a multi-band approach to blend low and high frequencies over long and short ranges, respectively, where weights are smoothed ratios of projected area over actual area. More recently, Totz et al. [START_REF] Totz | Patient-specific texture blending on surfaces of arbitrary topology[END_REF] propose a multi-band weighting based on the viewing distance and angle, where textures are averaged after flattening based on diskhomeomorphic patches that overlap for increased robustness.

In the meantime, some alternatives to weighting schemes have also been devised. Rocchini et al. [START_REF] Rocchini | Multiple textures stitching and blending on 3D objects[END_REF] select textures by least angle, smooth out selection along adjacent vertices to obtain larger contiguous regions that map to the same image, and only average border faces using barycentric coordinates as weights after a local linear registration. Lensch et al. [START_REF] Lensch | A silhouettebased algorithm for texture registration and stitching[END_REF] refine this approach and introduce additional heuristics to discard points based on the viewing angle, depth variation and proximity to the projected outline.

A common drawback of all these approaches is that they do not correct for high color discrepancies due to different acquisition settings or light variations across views. Moreover, ghosting and blurring artifacts appear as soon as the textures are misaligned: if the transition width of the weighting function is too large, then a slight misalignment produces a double-exposure effect where some features appear twice, whereas if it is too small, then texture merging is made within a few pixels only which is prone to seams.

To tackle the former issue, Beauchesne and Roy [START_REF] Beauchesne | Automatic relighting of overlapping textures of a 3D model[END_REF] propose a preprocessing step of automatic relighting, assuming almost convex objects, where differences in illumination between two images are corrected based on a lighting ratio. Agathos and Fisher [START_REF] Agathos | Colour texture fusion of multiple range images[END_REF] account for lighting and camera calibration by applying both a global linear correction estimated from overlapping views, and a local correction to smooth out remaining seams based on iterative averaging along the mesh connectivity. Bannai et al. [START_REF] Bannai | Multiple color texture map fusion for 3D models[END_REF] extend this to remove propagation errors by accounting for a groupwise transform rather than pairwise transforms between images, where a single joint problem is solved instead of several binary problems. This approach is yet computationally demanding, requiring several hours to run. Moreover, these approaches do not deal with issues of misalignment.

Alternatively, Lempitsky and Ivanov [START_REF] Lempitsky | Seamless mosaicing of imagebased texture maps[END_REF] apply a global optimization based on a Markov random field to favor smoothness in texture assignment and penalize sharp seams, followed by a seamleveling scheme similar to gradient-domain stitching. Gal et al. [START_REF] Gal | Seamless montage for texturing models[END_REF] extend this by also optimizing with respect to local translations to overcome misalignment problems. Following the success of gradient-domain stitching methods, Chuang et al. [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF] formulate a robust stitching scheme via a screened Poisson equation, although the approach lacks a local alignment scheme. The interest of such methods is that the texture reconstruction error is optimized in the gradient domain, leading to minimized perceptual errors in respective frequency bands and accounting for large color offsets between views, without requiring any weighting heuristics. Moreover, they include texturing of unobserved points as part of the process, hence removing the need for postprocessing heuristics toward hole filling.

Contributions

In this work, we develop a novel patch-based approach for texture stitching on a 3D mesh. Our methods are fully intrinsic to the mesh surface, generic in that we make no restrictive assumptions on the mesh topology (e.g., watertight, closed) or geometry (e.g., convex, smooth), and only require a couple of minutes to run. After image sampling, we perform gradient-based stitching within a uniform mesh segmentation. Patch textures are first selected globally by least angle, and gradient textures are then retained locally in patch overlaps. We use an original variant on Poisson blending by coupling discrete differential operators with screening. We also consider an extension to handle color transformation via calibration and correction. Results are presented on both a zippered mesh based on range scans and a morphable model fitted to face photographs.

PROPOSED METHODS

In this section, we expose our methods for seamless texture stitching. The pipeline includes core steps of image sampling, mesh segmentation and gradient stitching, with an optional color transformation.

Image sampling

For each view, we determine the set of visible vertices on the 3D model. Occlusions can be tested by ray casting from the viewer, although we here use an approximate but much less demanding depth criterion with a z-buffer. The image is then sampled on the mesh by back-projection of textures for visible vertices after bilinear color interpolation within the pixel grid. Additionally, the viewing angles for each face and vertex are computed as part of the process and stored for later use. We notice finally that any of the heuristics discussed in the introduction could be employed here to discard some vertices that are potentially corrupted by considering them as occluded.

Mesh segmentation

We segment the mesh with a farthest-point strategy based on fast marching [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF], enhanced with an original patch growing step to form an overlapping structure. This produces a uniform segmentation compared to the region-growing scheme in [START_REF] Totz | Patient-specific texture blending on surfaces of arbitrary topology[END_REF]. Moreover, we do not require patches to be disk-homeomorphic since our methods are intrinsic to the mesh surface. Not only does it eliminate undesirable distortions inherent to flattening, but it also allows to consider various mesh topologies with arbitrary genus and number of boundary components. Lastly, instead of using extrinsic 3D ball radii as thresholds for patch dilation, overlaps are obtained by growing patches intrinsically within neighbors via geodesic projections.

We consider a triangular base shape mesh M and assume it describes a 2D manifold with geodesic distance map D. The connectivity is given by a simplicial complex K whose elements are vertices {i}, edges {i, j} or faces {i, j, k}, with indices i, j, k ∈

[1 .. N ],
where N is the number of vertices. We write a vertex {i} as i for simplicity. We first select vertices iteratively by adding a new sample one at a time. Denoting by D l (i) the geodesic distance map to the first l selected samples, we select sample i l+1 as the vertex that maximizes D l (i). The distance map D l+1 (i) can simply be updated as the minimum between D l (i) and D(i, i l+1 ). We continue this process until a desired number M of vertices have been sampled.

Patches P1, . . . , PM are then obtained via the geodesic Voronoi tessellation based on the samples. The segmentation thus defines a dual graph G = (V, E), where V = [1 .. M ], and (m, n) ∈ E if Pm and Pn are neighbors, i.e., are connected by an edge {i, j} ∈ K.

To grow a patch Pm, we consider separately each of its neighbor patches Pn with (m, n) ∈ E, and define thresholds dmn as follows:

dmn = σ × D(i m , i n ) , (1)
where σ ≥ 0 is set by the user and can be seen as an overlap ratio or factor, and the geodesic distance D is restricted to the union Pm∪Pn of the reference patch and considered neighbor. The overlap Omn of Pm onto Pn is then constructed by geodesic projections:

Omn = i ∈ Pn : min j∈Pm D(i, j) ≤ dmn . (2)
A given grown patch Qm is eventually constructed by concatenation of the reference patch Pm with the respective overlaps:

Qm = Pm ∪ n|(m,n)∈E Omn . (3)
The time complexity of the whole process is in O(N log N log M ).

Gradient stitching

We stitch textures in the gradient domain by modifying a principled approach to Poisson's equation on a mesh with Dirichlet boundary conditions [START_REF] Yu | Mesh editing with Poisson-based gradient field manipulation[END_REF] so as to incorporate a screening term instead. This allows both to consider meshes with no boundary without having to arbitrarily define a virtual one, and to relax hard constraints on texture values at boundary vertices via a soft regularization distributed over all vertices. Since we rely on discrete differential operators on the mesh surface, our approach completely preserves conservative vector fields compared to extrinsic 3D finite elements in [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF]. This makes our approach more natural from a theoretical perspective, even if non-conservative fields are rather formed in practice.

A discrete vector field V is a piecewise constant vector function defined for each triangle T l by a coplanar vector v l . A discrete potential field is a piecewise linear function φ(s) = i∈K φiBi(s) on the mesh surface, where Bi is the piecewise linear basis function valued 1 at vertex i and 0 at other vertices, and φi specifies the value of φ at vertex i. The discrete gradient of φ for triangle T l is ∇φ l = i∈K φi∇B il , where ∇B il is the gradient of Bi within T l . The divergence of V at vertex i is div

V(i) = T l ∈K i |T l | ∇B il v l ,
where Ki is the set of triangles sharing vertex i and |T l | is the area of triangle T l . Writing Poisson's equation div ∇φ = div V in this framework leads to a linear system of equations Ax = y for the unknown potential values xi = φi, where:

aij = T l ∈K i |T l | ∇B il ∇B jl , yi = T l ∈K i |T l | ∇B il v l . (4)
This system is sparse since the sum for coefficients aij is non-null iff {i, j} ∈ K (it is an edge). The sum is then simply over the triangles T l (two if not a boundary edge, one otherwise) sharing this edge.

This equation can be interpreted as seeking for a potential field φ whose gradient ∇φ matches the guide vector field V . If V is conservative, i.e., it is the gradient of an existing potential field φ, then φ is the exact solution. Otherwise, a more general minimizer can still be obtained by least squares but its gradient differs from V . In addition, we regularize the minimization via screening:

min x∈R N Ax -y 2 2 + λ x -x 2 2 , (5)
where λ > 0 and x defines a guide potential field φ as φ i = x i . We apply this to solve for texture by considering each color channel independently as a potential field φ. For each view v, we compute the mean viewing angle of vertices in the different patches. Unobserved vertices, due either to occlusion or missing information, are assumed to have a viewing angle of π/2. Hence, patches with unobserved data are penalized and no difference on the nature of non-observability is made. For each patch now, we select texture from the view where the patch has the smallest viewing angle. For unobserved vertices, we also select texture from subsequent sorted views. We end up with partial textures φ (v) that we stitch in overlaps by Poisson blending. To build up the guide vector field V , we select local texture gradients by least angle for each triangle T l :

v l = i∈K φ (v l ) i ∇B il , (6)
where v l is the view whose angle is minimal for triangle T l . We also fill in unobserved faces simply by setting their gradients to zero for smoothness. Screening is done via a rough estimate φ obtained by averaging textures φ (v) , unobserved textures being discarded from the regularization. We use a small penalty λ = 10 -6 to remove color offset indeterminacies since we did not observe dramatic color bleeding issues compared to [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF]. The time complexity for building up the linear system is in O(M + N ). A naive complexity for least squares optimization via Cholesky decomposition is in O(N 3 ), though efficient solvers that exploit sparsity can be used instead.

Color transformation

To reduce texture discrepancies between views, we introduce two optional steps of color transformation via calibration and correction.

Compared to [START_REF] Bannai | Multiple color texture map fusion for 3D models[END_REF], we consider pairwise interactions for computational efficiency, and assume that more subtle effects due to propagation errors are fixed by gradient-domain stitching. However, we allow for an affine instead of linear transform, which further accounts for differences in color offset and ambient light intensity, rather than just color gain and directional light intensity. On the one hand, calibration is done if color charts are available. We first average pixels for each checking color to produce color maps under the different views. We then choose the color map of an arbitrary view or a virtual one (e.g., by averaging color maps) as reference, and solve for multiple pairwise affine transforms from color maps to this target by least squares. We finally apply these respective transforms to the sampled textures. On the other hand, correction is estimated directly from the sampled textures. Here, a robust regression method (e.g., least absolute deviations, least trimmed squares) is necessary to cope with potential outliers. We choose a reference view and compute pairwise affine transforms where views overlap. This not only compensates for different lighting conditions, but also for acquisition settings if views have not been calibrated before.

EXPERIMENTAL RESULTS

In this section, we present experimental results to evaluate our methods on real-world datasets. We consider both a zippered mesh based on range scans and a morphable model fitted to face photographs.

Range scans

We compare our methods on range scans with the only available tool we have found: the state-of-the-art texture stitcher [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF]. We use their Rooster dataset, which consists of color information and a base shape mesh reconstructed from 8 depth maps by zippering (Fig. 1). The mesh features 68,612 vertices, a positive genus due to 1 visible handle, and 1 boundary component as a result of a small hole in the top. For a fair comparison, we do not apply color transformation since their tool does not include such a step. To assess the effect of the number of patches and benefits of using overlap, we try different values of M, σ. For illustration, we also compare to two baseline vertex-based strategies: averaging and least-angle selection.

On the one hand, these baseline strategies suffer from texture discrepancies across scans. Averaging textures produces seams at boundaries where some scans become or cease to be observable, and a general color offset resulting in a darker texture due to accounting for shadowed vertices (Fig. 2a). Selecting the less foreshortened vertices addresses these issues, but introduces more marked seams at transitions between selected scans (Fig. 2b). Moreover, inherent to these two approaches, many unobserved vertices are not textured.

On the other hand, our proposed approach tackles these issues at once and leads to much smoother color transitions. For a small number of patches (M = 50, σ = 0.0), however, we observe a blurring effect around the eye (Fig. 2c). Increasing the number of patches (M = 100, σ = 0.0), we are able to remove this artifact, yet we still observe slight seams (Fig. 2d). The best result is when further growing the patches so as to overlap (M= 100, σ= 0.5), where the obtained stitching is seamless (Fig. 2e). Looking at the frontal view, our result preserves sharp edges of specularities, details of the painted pattern, and color differences between the pattern and body (Fig. 2f). In contrast, the texture stitcher [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF] blurs the specularities with the painted pattern, misses in particular the central specularity, and features a noticeable color bleed from the base so that details in the bottom half of the painted pattern are lost (Fig. 2g). The improvement is confirmed by computing the root mean square errors between scans and rendered models, with an average of 0.099 across views for our methods against 0.109 for the texture stitcher [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF].

Face photographs

We now consider photographs from the CMU PIE face database [START_REF] Sim | The CMU Pose, Illumination, and Expression (PIE) database[END_REF]. We choose 3 views of a given subject under the same illumination (Fig. 3). We fit a state-of-the-art 3D morphable model [START_REF] Paysan | A 3D face model for pose and illumination invariant face recognition[END_REF] to estimate shape, texture, camera and lighting parameters for the 3 (f) Idem.

(g) Stitcher [START_REF] Chuang | Estimating the Laplace-Beltrami operator by restricting 3D functions[END_REF].

Fig. 2. Texture reconstruction of the Rooster dataset. We compare our patch-based approach to two baseline and one state-of-the-art vertexbased approaches. Merging vertex information is not sufficient on this difficult dataset, whereas patch integration improves quality by removing seams while preserving details, as soon as we blend the solutions with enough patches and overlaps for robustness and smoothness. The results show that rendering a model fit [START_REF] Paysan | A 3D face model for pose and illumination invariant face recognition[END_REF] provides an inaccurate estimate of the actual face appearance and misses local asperities because of statistical regularization and smoothing (Fig. 4a). In contrast, our core stitching pipeline succeeds in improving texture details, although subtle color discrepancies still remain in that the right side of the face is lighter than the left one (Fig. 4b). The optional step of color transformation allows to tackle this latter problem so that the colors of both sides are more consistent (Fig. 4c). The results are again corroborated by the root mean square errors, which average to 0.0710 for the morphable model [START_REF] Paysan | A 3D face model for pose and illumination invariant face recognition[END_REF], against 0.0577 for our core pipeline, and 0.0570 with the additional color calibration.

CONCLUSION

We presented a novel technique for seamless texture stitching on a 3D model by Poisson blending in overlapping patches on the mesh surface. The obtained results demonstrate the relevance of our approach. Several perspectives were however left out for future work.

A drawback of our affine color transform is to ignore nonlinearities due to specularities and shadows, so that an enhanced lighting model could improve this step. In addition, we would like to evaluate color correction in more difficult conditions where calibration charts are not available and lighting varies between views. It would also be worth integrating a shape correction step to improve robustness against dramatic misalignment issues. Lastly, a more ambitious line is to couple shape and texture reconstruction in a unified paradigm.

  This work was supported by grant DSTLX1000070369 from the Defence Science & Technology Laboratory.

Fig. 1 .

 1 Fig. 1. Scans of the Rooster dataset. Colors from the 8 depth images are back-projected to the base shape mesh. The different views reveal substantial amounts of unavailable information and variations in illumination, hence seamless stitching for texture reconstruction is nontrivial.

Fig. 3 .

 3 Fig. 3. Photographs from the CMU PIE face database. The 3 poses of a given subject are taken simultaneously with different cameras. Subtle texture discrepancies appear across views because of different acquisition settings, which might lead to artifacts such as seams or different eye colors if textures are not stitched properly.

Fig. 4 .

 4 Fig. 4. Texture reconstruction of face photographs on the morphable model. Rendering the texture from the fitted model does not provide a realistic appearance, whereas our result features plausible local asperities. However, care should be taken to compensate for texture discrepancies via color transformation and avoid related artifacts.