Ying Huo

Daniel Delahaye

Ji Ma

Mohammed Sbihi

Integrated Traffic Flow Based Optimization of Airport and Terminal Area

Keywords: Integrated Optimization, Terminal Maneuvering Area, Airport, Simulated Annealing

As air traffic demand exceeds available capacity, severe congestion during peak periods occurring at airports and the surrounding terminal airspace leads to flight delays as well as potential safety issues. This paper addresses the optimization problem of integrated traffic management of airport and terminal airspace on macroscopic level. Instead of evaluating the flight conflicts as we did in our previous research, a congestion evaluation model is proposed to quantify the flights for traffic flows in the terminal airspace and the runways, meanwhile, the occupancy evaluations of taxiway network and terminals are conducted. An adapted heuristic optimization method simulated annealing is implemented to solve the problem. Optimization process is put forward based on the case studies of Paris Charles De-Gaulle airport with 2.5-hour data. The comparisons between initial traffic data and optimized solution are provided and the advantage of the proposed model is analyzed. In the aspect of time uncertainty, simulations based on the proposed model and the model of our previous work are conducted, the final results indicate that our model shows an advantage in uncertainty absorption over the previous work.

I. INTRODUCTION

As air traffic congestion increases, air traffic delays become more and more severe. The frequent communication all over the world urges the development of air transportation. The imbalance between demand and available resources of air transport causes a significant growth in delay, which increases air traffic service dissatisfaction and economic loss. Airport and its surrounding airspace are considered as a very complex region, thus the study of air traffic flow management (ATFM) concentrated in airport and Terminal Maneuvering Area (TMA) needs technical operation strategies and fully analytical evaluations.

Due to the distinct and complex characteristics of terminal airspace and airport, segregated researches have been studied widely. Murc ¸a et al. [START_REF] Murc | Control-based optimization approach for aircraft scheduling in a terminal area with alternative arrival routes[END_REF] studied the flight scheduling problem in the terminal area by considering the dynamic characteristic of flights in the airspace and the benefits of alternative routes. In [START_REF] Deau | Airport surface management and runways scheduling[END_REF][START_REF] Kjenstad | Integrated surface and departure management at airports by optimization[END_REF][START_REF] Balakrishnan | A framework for coordinated surface operations planning at dallas-fort worth international airport[END_REF], runway scheduling and taxiway scheduling were studied to enhance the coordination of the integrated airport management.

Although airport and TMA are usually considered as two individual areas due to their different configurations, the fact that airport and its surrounding airspace have a direct and interacted connection attracts considerable interests towards the integrated research of terminal area and airport including integration of arrival management, departure management and surface management. The integrated research related subproblems of routing, sequencing and conflict resolutions are investigated as well [START_REF] Kjenstad | Optimizing aman-sman-dman at hamburg and arlanda airport[END_REF][START_REF] Pavese | Dman-sman-aman optimisation at milano linate airport[END_REF].

Though there are lots of researches focusing on the microscopic level, the real-time management of airport and flight operations still rely on the decisions of air traffic controllers who provide feasible flight instructions according to operating rules and personal experience. Considering this situation, our study addresses the problem in a macroscopic perspective and tries to provide a solution for the integrated airport and terminal management without considering the specific instructions from the controllers.

In our previous work [START_REF] Ma | Integrated optimization of terminal manoeuvring area and airport[END_REF], the integrated problem of terminal airspace management and airport management on a macroscopic level was studied, however, the presented evaluation model took microscopic metrics of flight by flight conflicts into consideration. On macroscopic level, air traffic should be considered in a larger scale, while the investigation of detailed flight interactions are easily effected by minor uncertainty that generated during the flight operation. In order to emphasize the congestion elimination and robustness at the same time, traffic flow-based approach is proposed. The objective of our study is to mitigate the congestion and conflict resolution workloads in the terminal area and reduce the overload of each element in the airport. The proposed model focuses on the density of traffic flow which is less sensitive in terms of time variation comparing to the model in our previous work.

The paper is organized as follows: section II provides an overview of the proposed network abstraction. The mathematical problem is clarified, the evaluation method is described as well. Section III displays the case study based on real data of Paris-CDG airport. In section IV, optimization results are analyzed, and the results of the simulations regarding time uncertainty are displayed. Section V concludes the whole paper.

II. MATHEMATICAL MODELING

Our model is developed to contribute to both arrival and departure flights optimization in the airport and its surrounding airspace. A network abstraction is established to centralize the optimization process on network resources. The resources in the system are TMA routes, runways, taxiway network and terminals. With the proposed method, flight movements are grouped and quantified by time intervals.

A. Network model of TMA and airport surface

In the TMA, the route network is represented by links and nodes. The node-link structure can be presented as graph(N , L), where N is the node set and L is the link set. We have N = N e ∪ N ω ∪ {r} that constitutes the route network, where N e ={MOPAR, LORNI, OKIPA, BANOX} represents the set of entering points of the TMA, r represents the node at the runway threshold, and N ω = {n 2 e , ..., n r-1 e } is the set of nodes between the entry point e and the runway threshold. Similarly, we have L = L e ∪L ω ∪{l r }, where L e = {1, ..., l e } is the set of links that connects the entry nodes, l r is the last link of the route connecting runway threshold, L ω is the set of remaining links. Each entry point e ∈ N e corresponds to one route that is defined by r e = {e, l 1 e , n 2 e , l 2 e , ..., n r-1 e , l r , r} in the airspace with a slight change at the last node due to different runway assignment. Each route is defined by a succession of nodes and links, flights follow one of these routes according to their entry points and landing runways.

Fig. 1 shows the arrival route network abstraction of Paris Charles-De-Gaulle (CDG). One of the arrival routes from OKIPA to runway 26L is marked and relevant notations are displayed. We consider the links that converge to one node as the parent links of this node. In our network, we have |N | = 15 nodes and |L| = 15 links.

Three main resources in the airport are considered: runways, taxiway network and terminals. Each of the components will be considered as a resource with a certain capacity or relevant congestion constraints, the network resources are generally represented as S = N ∪ L ∪ R ∪ T e ∪ t a , where R and T e denote runway sets and terminals respectively. t a represents the taxi network.

B. Notations

Assuming that we are given a set of flights F = {1, ..., F }, among them, three types of flight operations are included, F = A∪D∪AD, where A stands for arrival, D stands for departure and AD represents connected flights. • C f : Wake turbulence category (f ∈ F).

• M f : Assigned terminal (f ∈ F).

• E f : Entry node of the TMA (f ∈ A ∪ AD). f is accessible. The exact time at which a flight arrives at a resource can be calculated with some assumptions:

• In TMA, we assume that the speed of flights will constantly decelerate until the Final Approach Fix (FAF) and then keep constant till the threshold of the runways. Final speeds for arriving at runways are set as 110 kt 130 kt and 150 kt for small, medium and large aircraft respectively. • The runway occupation time is obtained from the flight operation standards regarding the flight category. • The taxi duration is measured empirically. According to the airport operation information, the average taxi time for each flight can be obtained with regard to the combination of terminal and runway. Tab. I and Tab. II show the taxi-in and taxi-out time that we use in this paper. • The in-block time for each arrival flight can be obtained by adding the runway occupancy time and average taxi-in duration to the real landing time.

• For the departure flights, the evaluation is conducted until the take-off segment, the take-off time is estimated by adding up the taxi-out duration and the off-block time. • The take-off time of a connected flight can be computed by using the in-block time plus the turn around duration and average taxi-out time. By using the flight operation standards, we can obtain the runway occupancy time. Along with the given information, the operation time frame of all flights can be gathered.

C. Decision variables and constraints

The optimization model that we propose contains four types of decision variables: TMA entry time, TMA entry speed, landing runway and push back delay. The associated parameters are selected and interpreted as follows:

1) Entry time: in the TMA, flights (f ∈ A∪AD) can adjust the time to enter the TMA by changing the en-route speed. Discrete time is used and ∆T which denotes the time slot is considered as one unit to measure the total time range. The maximum tardiness and earliness are represented as ∆T max and ∆T min . The possible time range for each decision variable is composed of multiple ∆T . Therefore, for each flight (f ∈ A ∪ AD), a decision variable t f ∈ T f has a flexible range of: T f = {T o f + j∆T |∆T min /∆T j ∆T max /∆T, j ∈ Z}, where j denotes the interest of time slot deviation from the initial arrival time. Considering the flight operations in reality, the entry time range is set as ∆T max = 30 min and ∆T min = -5 min.

2) Entry speed: We discretize the entry speed in an similar way as the entry time. For flights f ∈ A ∪ AD, an entry speed decision variable v f ∈ V f has a constraint of:

V f = {V min f + j∆ v f |0 j (V max f -V min f )/∆ v f , j ∈ N}. In this study, we set V min f = 0.9V o f , V max f = 1.1V o f and ∆ v f = 0.01V o f .
We also guarantee a prior requirement that V max f do not exceed the maximum speed limited by the aircraft type.

3) Runway assignment: For f ∈ A ∪ AD, landing runway assignment decisions are represented by r a f ∈ R a f . Two landing runways are dedicated to arrival flights to balance the capacity of runways during the peak hours of operation. Departure runway assignment is not included in our research.

4) Pushback delay: For f ∈ D ∪ AD, it has been verified that the delay on the ground has more economic benefits than in the air. Therefore, pushback delay rather than air delay is considered as decision variable p f ∈ P f :

P f = {P o f + j∆T |0 j ∆T p max /∆T, j ∈ N},
where ∆T p max denotes the maximum allowed delay time. In our case, the maximum delay is chosen to be 20 min which is a reasonable time range in practice.

To summarize, the decision variables associated to flight f denoted as x f .

x f = (t f , v f , r a f , p f ) Let's denote x = (x 1 , ..., x F ) as a set of decision variables of all flights.

D. Objective

Our models are designed to eliminate the conflict resolution workload in the air and alleviate congestion in the whole system. The resources in this study are evaluated separately and the objective function is thus a weighted sum of evaluation metrics of five resources:

min s∈S γ s • G s (x), (1)
where γ s is a user-defined weighting coefficient for the congestion or exceeded overload of resource s in the model system. G s (x) is the corresponding objective for each resource s.

The objective is evaluated by two metrics, the first one can be interpreted as the average congestion during time horizon T which is the mean capacity exceeding value for all evaluation in time horizon. Another important metric measures the maximum overload, which is the difference between the maximum overload and the corresponding capacity. Mean congestion provides a general perspective of network situation in terms of flight numbers, and the maximum overload indicates the period at which severe congestion occurs. A general objective function G without evaluation details is formulated as follows:

G(x) = s∈S γ s 1 T t∈T [max(E s (t) -C s , 0) • ∆] + s∈S γ s max t∈T [max(E s (t) -C s , 0)] , (2)
where E s (t) represents the evaluation result of resource s at the evaluation time indicated by t, and C s indicate the resource capacity. We define E s to be the flight number distribution for resource s along T . ∆ is a time unit for the discrete time horizon.

Once the decision variables are set, we can calculate the time when aircraft arrive at each resource. Then we use the desired time values to compute E s . Considering the different characteristics of each resource, the evaluation of E s need to be specified. Evaluation process is introduced in the following section.

E. Traffic flow-based evaluation

Traffic flow-based approach is implemented on the operated traffic flow which is a series of scheduled flights specified with operation time. We investigate the traffic flow by focusing on specific time intervals. For each time interval, the number of flight that entry or exit the network resources are counted. In order to detail the evaluation and enhance the accuracy of our model without a huge loss of the flight information, a series of evaluations are conducted based on a fixed time period, for each evaluation, the starting time is shifted by one time step.

In the process of traffic flow-based evaluation, ∆ (e.g. 1 minute) is set as a time unit that constitutes the evaluation time interval and as a time step to proceed to the evaluation process. According to different attributes on each resource, the evaluation time interval L = k • ∆ should be set differently where k is a coefficient to define the length of the evaluation time.

Fig. 3 illustrates the evaluation time scheme based on the traffic flow-based approach. The blue strips indicate a time interval of L for each evaluation in which the number of flights that entry or exit a specific resource are counted. A continuous moving of the evaluation time is proceeded by shifting one time step afterward.

1) Link: The preliminary evaluation for link is the number of flights that enters a link. With the utilization of traffic flow-based approach, the evaluation time interval should be decided. Considering that a short evaluation time interval is not able to reflect the flight density and can not achieve a proper comparison with link capacity, thus k = 10 is used to define the link evaluation time.

For the link evaluation, the number of flights that enters link l from t to t + L is evaluated as follows:

E l (t) = t+L i=t f ∈F ω i l,f l ∈ L, (3)
where ω i l,f is a binary parameter that equals to 1 if flight f is supposed to enter the link l during time interval i to i + L, otherwise equals to 0.

Based on the aforementioned assumption, flights are conducted with a constant deceleration procedure which begin from the entry point of TMA until the FAF. Regarding the different spatial locations of links, flight speeds are distributed differently on each link, so the average speed v l of all passed flights on link l is estimated. Along with the flight separation requirements, the acceptable number of flights that enters link l during the evaluation period can be obtained and set as a capacity of link l for the following congestion evaluation.

The number of flights on each link in the TMA is constrained by M = 1 flight per 5NM which is derived from the flight operation rules. Based on the average speed of each link, the associated link capacity C l can be computed as:

C l = M • v l K ( 4
)
where K is a time scale. As in link, k equals to 10, ∆ is fixed to be one minute and the unit of link average speed is NM/h, K = 1 hour/L is a scale to convert the flight entry number of a specific link of 1 hour into a time of 10 minutes.

2) Node: In the TMA, nodes evaluation is focused on the estimation of potential conflicts resolution workload. Using the evaluation model proposed in [START_REF] Delahaye | Modeling and optimization of air traffic[END_REF], the quantification of conflicts resolution workload is determined by two converging traffic flows which perfectly match the scenario of nodes and their parent links. Suppose that two links l i and l j are parent links of node n, flight speed on the links are given as v li and v lj respectively. E out li and E out lj are the distributions of number of flights that exit the parent links during the time horizon, thus the conflicts resolution workload E n (t) on a specific time interval indicated by t can be computed as follows:

E n (t) = 2N s v 2 li -2v li v lj cos θ ij + v 2 lj v li v lj sin θ ij •E out li (t)E out lj (t)n ∈ N , (5)
where N s is a parameter for the standard separation norm, θ ij denotes the angle formed by two convergent flows associated with link l i and l j .

For nodes, the conflict resolution workload is set to be 6 per hour regarding the workload of air traffic controller.

3) Landing runway: Runway flight evaluation counts the number of flights that pass through the runway threshold from time t to time t + L, where k = 10 (L = k • ∆) ensures an adequate time interval for avoiding the potential risk of losing flight information. Along with the time going, flight entry distribution E r of the whole time horizon regarding arrival runway is formulated as Eq. 6 and the runway operation limitation can be imposed to calculate the sub-objective generated by runway.

E r (t) = t+L i=t f ∈F ω i r,f r ∈ R, (6)
same as in Eq. 3, ω i r,f is a binary parameter to indicate that flight f uses runway r during time interval i to i + L and equals to 0 otherwise. 4) Taxiway network and terminal: Taxiway network and terminal are resources with specific capacities. We can verify whether the number of flights exceeds the capacity or not at any time, therefore the congestion evaluation is conducted on each time step with a time interval of L = ∆. In this case, there is no overlap of evaluation time, evaluation intervals are separated and continuous. Flights that enter and exit the resource are recorded at each time interval to evaluate the number of flights that stay on terminal or taxiway network during time t to t + L.

E ta (t) = ( f ∈F ω t ta,f - f ∈F δ t ta,f ) + Q ta (7)
Eq. 7 takes taxiway as example. ω t ta,f is an mentioned binary parameter for taxiway. δ t ta,f is also a binary parameter that equals to 1 if flight f is supposed to proceed an exit operation of taxiway during time t to t + L and otherwise equals to 0. Q ta denotes the initial number of flights that stay in the taxi network at the time that we start the current evaluation. The evaluation process of terminals is the same as in the taxiway network. Based on flight number distribution, the two objective metrics regarding the specific resources can be computed.

III. SOLUTION APPROACH

It is known the computation for this problem is very expensive and NP-hard [START_REF] Beasley | Scheduling aircraft landingsthe static case[END_REF]

(|T f |•|V f |•|R f |•|P f |) |AD| +(|T f |•|V f |•|R f |) |A| + |P f |) |D| .

A. Simulated annealing

Simulated annealing is a meta-heuristic algorithm to approximate global optimization in a large search space for an optimization problem. It is used to find an approximate global optimum other than a precise local one in a fixed amount of time. This approach is interpreted as a slow decrease in the probability of accepting worse solutions as the solution space is explored. At each time step, the algorithm selects a solution close to the current one noted as its neighborhood solution, the quality of this solution is measured, and the decision is made to move to it or to stay with the current solution based on the fact that the new solution is better or worse than the current one. With this procedure repeating, the approximate global optimum solution can be achieved [START_REF] Aarts | Simulated annealing and boltzmann machines[END_REF].

In our case, in order to mitigate computational burden, a performance indicator is introduced to increase the probability of finding a better neighborhood solution. During the evaluation, flights who have a chance to violate the capacity constraint in any resource are marked with the overload. Total performance indicator is represented as m f for flight f which is the sum of exceeded capacity that flight f has experienced in the network. That is to say, if on a specific resource s, E s (t) exceeds the capacity at evaluation time interval t, and flight f is one of the flights that contributes to this overload, then this overload will considered as a part of performance indicator and added to m f .

In our algorithm, indicators are categorized into airspace performance a f and ground performance g f according to Algorithm 1 Neighborhood solution selecting method Require: For each flight, compute a f and g f and m f = a f + g f , generate a random number n=random(0,1);

p f a =a f /m f ; p f g = g f /m f ; 1: Choose one flight f based on its total performance; 2: if f ∈ A then 3: if n < p f a then 4:
one of t f , v f is chosen with equal probability and change its value. one of t f , v f and r a f is chosen with equal probability and change its value. 7:

end if 8: else if f ∈ D then 9: change P f 10: else if f ∈ A ∪ D then 11: if n < p f a then 12:
one of t f , v f is chosen with equal probability and change its value. Choosing and changing the decision variable of flight f as neighborhood according to neighborhood selecting rules; 6:

Calculate the new objective s i (x);

7: if s b (x) > s i (x) then 8: s b (x) ← s i (x); 9: else if β < exp( si(x)-s b (x)

Tc

) then able so as to generate the neighborhood solution. Algorithm 1 describes the procedure to find the neighborhood solution. The change of decision variable is the change of time deviation or speed deviation corresponding to the initial value with a range stated in the explanation for decision variables. Algorithm 2 presents the implementation of simulated annealing with the consideration of Algorithm 1 in which T c denotes the current temperature of the optimization process and s b (x) is the obtained best solution considering all the finished iteration.

IV. RESULTS

A. Case study background

The instances that we used in this paper are the real scheduling data of Paris CDG airport. In this airport, there are two runway configurations: West-flow (26L, 27R-26R, 27L) and East-flow ( 09L, 08R-09R, 08L). This paper focuses on the more frequently used West-flow configuration. The airport consists of three terminals, in which terminal 2 handles the most amount of flights.

The congestion evaluation standard of the related capacity of each component will be chosen according to the flight regulation and empirical estimation. The capacities of resources in the airport are estimated by using historical data which is shown in Tab. III. We would like to mention that in the taxiway, the acceptable number of flights is fluctuating due to different flight types, speeds or taxiway changing in practice. Thus the capacity of taxiway is confirmed during the process of problem solving by giving other resources priority. Then we choose the lowest flight accumulation in the taxi network as the taxiway capacity.

Actual flight data on 18th February 2016 from 7:00 AM to nearly 9:40 AM which was the peak hour of the day are studied in this work. Flights are mostly composed of heavy and medium types with a mix ratio of Heavy : Medium = 23% : 77%. The model was coded in Java and run on a 2.50GHz core i7 CPU, under Linux operating system. The required time for computation is around 5 minutes.

The resources in the network system are evaluated separately with different capacities, thus the respective analysis is made for each resource set. 2) Node: Since the evaluation time interval is set as 10 minutes, for each interval, the obtained data is finally transferred to an hour. As the conflicts calculation is highly related to the number of flights coming from the parent links, thus for some nodes, no overload is appeared during the evaluation time, but still, there are three nodes who have overloads from the initial scheduling data. Fig. 4 displays nodes that have their capacity exceeded with initial flight data are marked with the green circles.

Fig. 6 displays the comparison of conflicts resolution workload distributions before and after optimization, the indicated quantities of the red lines show that hourly conflicts resolution workloads are reduced to the required standard.

3) Runway: In each runway, similar evaluation process is conducted as in link. The flight number distributions are compared directly. From Fig. 7, we can see that the red

C. Robustness simulation

As a model focus on traffic flow evolution and evaluate a series of flight status at the same time, the advantage of uncertainty resistance can be found. Thus a simulation is conducted with an reference case of our previous study. The optimal solutions of traffic flow-based model and flight by flight conflict evaluation model which is proposed in [START_REF] Ma | Integrated optimization of terminal manoeuvring area and airport[END_REF] are obtained with the same input data. since the two models have different targets and distinctive perspectives, the robustness tests are conducted respectively with their own optimal solutions with respect to their own objective. In each simulation, random perturbations are drawn from uniform distribution of [-ε, ε] for three different values of ε. Each simulation has 5000 experiments and the average values of objective and subobjectives of each model are obtained and displayed in Tab. IV. For the two models, the sub-objective of airside evaluations are the same. In TMA, flight by flight conflict evaluation model computes the conflicts generated on each resource, while the proposed model take the average congestion and maximum overload into account. Two parts in the table are highlighted: the simulation results for the traffic flow-based model and for the flight by flight model. For each model, as the range of perturbation increased, the derived congestion or conflicts augment. One exception is revealed due to the inherent attributes that terminal is a resource with specific capacity and the small perturbations have minor effects on the number of flights in terminals. As we compare the result of two models with the same perturbation, regardless the values of the results or the growth rates, In this paper, network abstraction is proposed based on the distinctive characteristics of different resources so as to conduct separated investigation on each kind of resource. In order to mitigate the congestion in the network as well as increase the stability against time variation, traffic flow-based approach is put forward. In TMA as well as runway, flights are grouped and counted in a specific time interval, corresponding capacity is set so as to ensure the general safety and separation. In airside, the taxiway and terminal congestion are considered. An heuristic method of simulated annealing is proposed to solve the optimization problem. Practical instances of Paris CDG airport are used. Promising results are demonstrated in detail in terms of each resource. For verifying the model, a robustness test towards the proposed model and the model in our previous work is conducted, the results indicated that the traffic flow-based model is less sensitive than the other one.

As for future works, multiple airports system contains more uncertainties which yields a complex and imprecise study on microscopic level, thus the further research could be dedicated to a larger range of research of including other airports in the system.

Figure 1 .

 1 Figure 1. Arrival route structure

Figure 2 .

 2 Figure 2. Network model of TMA and airport surface

  • T o f : Initial RTA (Required Time of Arrival) of entering the TMA through the entry node (f ∈ A ∪ AD). • V o f : Initial speed of entering the TMA through the entry node (f ∈ A ∪ AD). • R a f : Initial assigned landing runway (f ∈ A ∪ AD). • t ad f : Turn around duration (f ∈ AD). • T d f : Scheduled off block time (f ∈ D). • R d f : Initial assigned departure runway(f ∈ AD ∪ D). The common information of all flights are C f and M f . For arrival flights, the information of E f , T o f , V o f , R a f are given. For departure flights, T d f and R d f are provided and for connected flights, except the information of arrival and departure flights, an additional data of t ad

Figure 3 .

 3 Figure 3. Traffic flow-based evaluation

  Considering an example of 100 flights containing 25 arrival flights, 50 connected flights and 25 departure flights. In this example, we assume that |T f | = 400, |V f | = 20, |R f | = 2 and |P f | = 200, then there would be a huge possible solutions to be considered. Due to this high combinatorics, simulated annealing algorithm is put forward to address the problem.

Algorithm 2

 2 of t f , v f , r a f and p f is chosen with equal probability and change its value.15:end if 16: end if the location of the constraints violation. The values of subperformance decide the selection of decision variable due to a fact that for a flight which encounters an airspace overload is useless to change the value of pushback delay. The designed algorithm considers all the marked flights. One of them will be chosen and assigned a new value for the picked decision vari-Simulated annealing application Require: Initialize( initial temperature T 0 , iteration I = 200, random number β ∈ [0, 1]); 1: Calculate the initial sequence objective s(x); 2: s b (x) ← s(x) 3: while T c > 0.0001 • T 0 , s b (x) is not optimum solution do 4:for i = 1 to I do 5:

T

  c = T c • 0.99; 14: end while 15: return s b (x)

B. Results analysis 1 )

 1 Link: Regarding the network abstraction, the link number is counted up to 15, by analyzing the initial flight data, only six of them experience congestion. Two of them are displayed as examples. In Fig.5, blue lines represent the original scheduling data derived evaluation distributions and the red lines are the optimized results obtained by heuristic

Figure 4 .

 4 Figure 4. Network congestion indication

Figure 9 .

 9 Figure 9. Initial and optimized gate occupancy comparison of terminal 2

  traffic flow-based model has a better performance than flight by flight conflict evaluation model. The results indicate that flight by flight conflict evaluation model is very sensitive to the time variation yet traffic flow-based model shows a strong resilience in the face of uncertainty. V. CONCLUSION This work contributes to the terminal and airside integrated flight scheduling problem in a macroscopic perspective. Based on the past deterministic work, new considerations are included in the research for each network resource. The goal of this work is to reduce the congestion and controller intervention workload in the terminal area and the overload in the airport.

  . If we consider |A| + |AD| + |D| flights with the entry time changes |T f |, entry speed changes |V f |, landing runways |R f | and pushback delay changes |P f |, the total number of possible combinations of decision variables is equal to

TABLE IV ROBUSTNESS

 IV TEST RESULTS OF TWO MODELS WITH THREE KINDS OF PERTURBATION.

			Traffic flow-based model		Flight by flight conflict evaluation model
			Perturbations according to uniform distribution on [-ε, ε]
	Average congestion/conflicts	Initial	ε= 0.5min	ε=1min	ε=2mins	Initial	ε=0.5min	ε=1min	ε= 2min
		congestion				conflict
	Node	0	6.585	10.1	16.56	0	21.07	44.882	78.37
	Link	0	0.474	0.963	1.62	0	4.94	43.975	144.55
	Runway	0	0	0.064	0.47	0	14.998	31.05	52.86
	Taxiway	0	0	0.00134	0.033	0	0.519	0.596	0.835
	Terminal	0	0	0	0	0	0	0	0
	Sum	0	7.06	11.136	18.69	0	41.54	120.5	276.6

ACKNOWLEDGEMENT

This work has been supported by China Scholarship Council (CSC). The author is gratefully acknowledge the support of Mr. Serge Roux for his assistance with technical improvements.

4) Taxiway: In Fig. 8, the blue line is the initial flight distribution on taxiway with respect of time which shows large fluctuation. The capacity setting is done in the process of optimization considering the performances of other resources. As shown in the red line, the occupancy of the taxiway can be effectively reduced with a smoother distribution which will provide a much safer environment during the operation of flights.

5) Terminal: Due to the airport configuration, although terminal 2 has most of the flight gates, it is the busiest terminal and easy to generate a congestion. Fig. 9 shows the comparison between initial gate occupancy distribution and the one after optimization. The capacity is set to be 130, while The flight distribution perform a tendency to manipulate the flights from the peak hour to an idle time period so as to ensure safety and make better use of the resources.


n , it can be geometrically represented by an ellipsoid E A = {x ∈ R n , x T Ax ≤ 1}, the lengths of E A 's semi-axis being 1/ λ i (A), if λ 1 , . . . , λ n denote the eigenvalues of A (the ellipsoid can be degenerated, with infinitely long semi-axis when λ i = 0). Then the Loewner ordering corresponds to a reversed inclusion ordering on ellipsoids: for A, B ∈ S + n ,

A ≥ B ⇐⇒ E A ⊆ E B .
To get a more intuitive geometrical correspondance where the ellipsoid associated to the smaller element is included in the one associated to the bigger element, we can consider the ellipsoid whose semi-axis's lenghts are √ λ i . When A ∈ S ++ n this corresponds to E A -1 , thus for A, B ∈ S ++ n ,

A ≥ B ⇐⇒ A -1 ≤ B -1 ⇐⇒ E B -1 ⊆ E A -1 ,
and this can also be extended for non invertible matrices 1 . S + n with the Loewner ordering is an anti-lattice An anti-lattice [START_REF] Kadison | Order properties of bounded self-adjoint operators[END_REF] is a poset in which, for any two elements A and B, a greatest lower bound exists if and only if A and B are comparable. In [START_REF] Kadison | Order properties of bounded self-adjoint operators[END_REF] Kadison shows that S n with the Loewner ordering is an anti-lattice. Moreland and Gudder [START_REF] Moreland | Infima of hilbert space effects[END_REF] proved that this is also true for S + n (that is, a greatest lower bound exists in S + n if and only if A ≥ 0 and B ≥ 0 are comparable), and Ando [START_REF] Ando | Problem of Infimum in the Positive Cone[END_REF] generalizes the result to infinite dimension.

The misleading geometrical interpretation Despite the aforementioned results, one can argue that it is possible to find a unique smallest (largest) ellipsoid containing (contained in) a set of centered ellipsoids. The problem in this definition is that it refers to two different orderings: whereas the set of upper/lower bounds is defined by the Loewner ordering (or equivalently the inclusion ordering on ellipsoids), the smallest/biggest element of this set is defined by a volume total ordering. Unfortunately, the volume ordering does not induce the inclusion ordering, and one can find upper/lower bounds in the Loewner order that are not comparable with the "smallest"/"biggest" (in the volume sense) upper/lower boundsee Figure 1.

The (questionable) proposition of [2]

The authors of [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] acknowledge that the Loewner ordering is not a lattice ordering. However, the paper builds another ordering that is not proved to induce a lattice and that, in fact, seems to suffer from the same problem (the new order seems to be an anti-lattice as well).

In [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF], computing the sup of a set of matrices in S + n boils down to finding the smallest sphere covering a set of non centered spheres. As in the case of the ellipsoids, here two orderings are mixed: the inclusion ordering and the volume ordering. Again, we can find examples, as in Figure 2, where the resulting "smallest upper bound" is not comparable to other upper bounds in the ordering defined by the authors. These counter-examples seem to prove that this kind of ordering does not define a lattice. Indeed, if it did, the two blue spheres of Figure 2 would admit a unique sup. Since inclusion implies the ordering on volume, this sup would be the sphere of minimal volume enclosing the two blue spheres, that is the green sphere. The red sphere shows the existence of upper bounds not comparabale to the green one, and contradicts the existence of a unique minimal upper bound. In the following I give a recap and reformulation of the construction presented in [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF].

Bases of S +

n and their extreme points Figure 3 provides an illustration of the structure of the cone S + n , with the main features used for the construction of a new ordering in [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF].

For any λ ∈ R, we note T λ = {M ∈ S n ,tr(M) = λ }. Then [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] recalls that B 1 := S + n ∩ T 1 , the set of positive matrices with trace 1, is a base of the cone S + n , that is to say: for any M ∈ S + n , M = 0, there exists a unique M ∈ B 1 and a unique µ > 0 such that M = µ M. Indeed, here µ = tr(M) and M = M tr(M) . Moreover, B 1 is convex and its set of extreme points is known: it is the set of matrices

E 1 := {vv T , v ∈ R n , ||v|| = 1}.
For a convex set B, extreme points ext(B) are defined as the points such that for any x ∈ ext(B), B \ {x} is still convex, and they have the property to belong to the boundary of the convex set.

In this precise case, E 1 has another interesting property: its elements lie on the Euclidean sphere of centre 1 n I n and radius 1 -

1 n : for any M = vv T ∈ E 1 , ||M -1 n I n || 2 = tr((M -1 n I n ) 2 ) = tr(M 2 ) -2 n tr(M) + 1 n 2 tr(I n ) = tr(M) -1 n = 1 -1 n .
Since extreme points are always included in the boundary of the convex set, and since B 1 is a base of S + n , knowing E 1 = ext(B 1 ) gives a good idea of the general shape of the cone. Indeed, we can view B 1 as some kind of convex polygone inscribed in the latter sphere, and then any "slice" B λ = S + n ∩ T λ , λ > 0 can be deduced from B 1 by homothety. We get that B λ is also a base of S + n , its extremal points are exactly λ • E 1 , and they lie on the Euclidean sphere S λ of center λ n I n and radius λ 1 -1 n . We note C λ = S λ ∩ T λ the intersection of that sphere with T λ .

The new cone [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] Just as

S + n = ∪ λ ≥0 B λ , a new cone C n can be defined as C n = ∪ λ ≥0 C λ .
Although this is not explicit in the paper, the new ordering defined in [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] is the one induced by C n : n . The gray polygones represent two bases of the cone, B 1 and an arbitrary B λ with λ > 1, intersections of the cone with the hyperplanes T 1 and T λ . The extreme points of the bases are marked by black dots (there are infinitely many of them in reality), and the central dashed line indicates the span of the identity matrix I n , orthogonal to each T µ . On each base B µ , the extreme points lie on a sphere S µ centered on µ n I n (marked by a black cross) and with radius r µ = µ 1 -1 n . Whereas S + n is the union of all the gray polygones B µ , the cone that defines the new ordering in [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] is the union of the C µ = S µ ∩ T µ , µ ≥ 0. Readily, we have S + n ⊆ C n , which means that the Loewner ordering implies the new ordering. Furthermore, we get an easy characterization of C n :

A ≥ B ⇐⇒ A -B ∈ C n . B 1 B λ
M ∈ C n ⇐⇒ M ∈ C tr(M) ⇐⇒ ||M -tr(M) n I n || ≤ tr(M) 1 -1 n ⇐⇒ tr(M 2 ) -1 n tr(M) 2 ≤ tr(M) 2 (1 -1 n ) and tr(M) ≥ 0 M ∈ C n ⇐⇒ ||M|| ≤ tr(M).
Recalling that a matrix M = a c c b is in S + 2 if and only if a ≥ 0, b ≥ 0 and det(M) ≥ 0, we can see that the Loewner ordering and the new ordering are equivalent for n = 2:

S + 2 = C 2 .
In practice: inclusion ordering on spheres The ordering in [START_REF] Burgeth | Mathematical morphology for matrix fields induced by the Loewner ordering in higher dimensions[END_REF] is directly defined in terms of inclusion of spheres. Using the notations introduced above, it can be written as

A ≥ B ⇐⇒ (B -C n ) ∩ T 0 ⊆ (A -C n ) ∩ T 0 ⇐⇒ B -C tr(B) ⊆ A -C tr(A) .
(B -C n ) ∩ T 0 is the intersection of the "ground" plane T 0 , and the reversed cone -C n whose vertex has been placed in B (see Figure 4). It is therefore the sphere of T 0 centered in m B = B -tr(B) n I n and of radius r B = tr(B) 1 -1 n . Similarly, (A -C n ) ∩ T 0 is the sphere of T 0 centered in m A = A -tr(A) n I n and of radius r A = tr(A) 1 -1 n . By applying a criterion on the distance between centres, we can check that this definition of the new ordering, based on the inclusion of spheres, is equivalent to the one we gave earlier:

A ≥ B ⇐⇒ ||m A -m B || ≤ r A -r B ⇐⇒ ||A -B -tr(A-B) n I n || ≤ tr(A -B) 1 -1 n ⇐⇒ ||A -B|| ≤ tr(A -B) ⇐⇒ A -B ∈ C n .

Conclusion on [2]

The paper proposes a new ordering but does not check whether it induces a lattice or not. From simple geometrical considerations it seems clear that it does not, and [START_REF] Joris Van De Gronde | Beyond Scalar Morphology: A Theoretical Framework for Mathematical Morphology on Colour and Tensor-valued Images[END_REF] asserts this ordering does not produce a lattice. I wonder if it is worth trying to show it is actually an anti-lattice.

This ordering and the definition of smallest upper bound may still be used as an approximation for morphological methods on S + n . The thesis [START_REF] Joris Van De Gronde | Beyond Scalar Morphology: A Theoretical Framework for Mathematical Morphology on Colour and Tensor-valued Images[END_REF] gives a quality measure of this approximation, I still need to have a look at it. However, it is not clear how it is better than the original Loewner ordering and the inclusion of centered ellipsoids. Therefore, if we decide to work with approximations of dilations and erosion, I would stay with the Loewner ordering.

2 Analysis of images of vessels based on structure tensors

Example of images

See Figure 5. In the following, we note f the image and suppose it is of size N × N.

Structure tensor

The definition and computation of the structure tensor follows G. Peyré's numerical tour: http://www. numerical-tours.com/matlab/pde_3_diffusion_tensor/.

Gradient The following centered finite difference approximation of ∇ f is used:

∇ f (x, y) = ( f X (x, y), f Y (x, y)) T , where f X (x, y) = 1 2 ( f (x + 1, y) -f (x -1, y)) for any y and 2 ≤ x ≤ N -1, f X (1, y) = 1 2 ( f (2, y) -f (N, y)) and f X (N, y) = 1 2 ( f (1, y) -f (N -1, y)) (and similarly for f Y (x, y)).
Tensor at scale zero The scale zero tensor T 0 maps each (x, y) to the symmetric matrix with rank ≤ 1

T 0 (x, y) = ∇ f (x, y) • ∇ f (x, y) T = f X (x, y) 2 f X (x, y) f Y (x, y) f X (x, y) f Y (x, y) f Y (x, y) 2 .
It is straightforward that ∇ f (x, y) is eigenvector of T 0 (x, y) with ||∇ f (x, y)|| 2 as corresponding eigenvalue, and 0 is the other eigenvalue. Hence for any (x, y), T 0 (x, y) is an extreme point of the cone of positive semi-definite matrices S + 2 .

Tensor T σ at scale σ > 0 It is a smoothed version of T 0 , obtained by convolving the latter with the σ -scale Gaussian kernel

T σ = G σ T 0 = G σ ( f 2 X ) G σ ( f X f Y ) G σ ( f X f Y ) G σ ( f 2 Y )
where G σ is the smoothing operator. Note that, for any u = (x, y), T σ (u) is a weighted sum of tensors T 0 (v) with positive weights, and is therefore a positive semi-definite matrix.

Eigen-decomposition of T σ For each u = (x, y), T σ (u) can be decomposed as

T σ (u) = λ 1 (u)e 1 (u) • e 1 (u) T + λ 2 (u)e 2 (u) • e 2 (u) T
where 0 ≤ λ 2 (u) ≤ λ 1 (u) are T σ (u)'s eigenvalues and (e 1 (u), e 2 (u)) its basis of orthogonal eigenvectors. If we represent T σ (u) by its corresponding ellipse

E u = {X ∈ R 2 , X T T σ (u)X ≤ 1}
then the main direction of E u is given by the second eigenvector e 2 (u). For σ = 0, this vector is orthogonal to the gradient ∇ f (u); more generally, for small σ , Peyré points out the Taylor expansion

T σ (u) = T 0 (u) + σ 2 H f (u) 2 + O(σ 3 ),
where H f 2 is the Hessian of f 2 . It is not clear what e2(u) represents for larger σ , but it seems that inside a vessel the ellipse's main direction roughly coincides with the vessel's main direction.

Trace and anisotropy images (Figures 7 and8) From the tensor field T σ we can build two scalar images: the anisotropy image I an and the trace image I tr , respectively defined as

I an (u) = 1 - 2λ 2 (u) tr(T σ (u))
and

I tr (u) = tr(T σ (u)) = λ 1 (u) + λ 2 (u).
We have 0 ≤ I an ≤ 1, and the closer I an (u) to 1, the more T σ (u) is anisotropic. I tr can be seen as a smoothed version of the square norm of the ∇ f . ) where p 0 is the origin and v = (1, 0) T ; Right: the elements form the left hand image that comply with the additional constraint that the (smaller) angle between the directions defined by --→ p 0 p and -→ v is smaller than π 6 .

Structuring elements

Here we define for each pixel u the neighboorhood to be taken into account in the computation of flat dilation and erosion. The idea is to choose in a spatial window around u, those pixels v with for which the main orientation of T σ (v) is consistent, in some way, with T σ (u). So far I have used co-circularity as consistency criterion.

Co-circularity Given two points p 1 , p 2 ∈ R 2 and two vectors -

→ v 1 , - → v 2 , we say that (p 1 , - → v 1 ) and (p 2 , - → v 2 )
are co-circular if there is a circle tangent in p 1 and p 2 to -→ v 1 and -→ v 2 respectively. An infinite radius for the circle corresponds to the case when --→ p 1 p 2 , -→ v 1 and -→ v 2 are colinear, which we consider a particular case of cocircularity.

Additional constraint As shown on the left of Figure 9, co-circularity allows "ladder" configurations, which we may want to discard in the following. If so, one can impose an additional constraint, namely that the angle between --→ p 1 p 2 and -→ v 2 is below a certain threshold.

Neighbourhood graph From the above we define an (undirected) graph (G, E) as follows: the set of nodes G is the set of pixels in the image {u 1 , u 2 , . . . , u N 2 }; (u i , u j ) ∈ E iff u j is in a square window of fixed size 2p + 1, centered on u i (that is, ||u iu j || ∞ ≤ p), and (u i , e 2 (u i )) and (u j , e 2 (u j )) are co-circular up to a certain angular tolerance, with the additional constraint described earlier to avoid ladder configurations -see Figure 9.

Morphological filters

Max-plus convolution Building on the graph (G, E) of the previous section, let W the N 2 × N 2 adjacency matrix defined by W i j = 1 ⇐⇒ (u i , u j ) ∈ E or i = j, and W i j = 0 otherwise. It encodes the structuring elements for any pixel, and we can now use it to compute max-plus convolutions (dilations

δ (I)(u i ) = j I(u j ) + log(W i j ) and ε(I)(u i ) = j I(u j ) -log(W ji ) . (1)
In particular, I can be the original image f , the anisotropy image or the trace image (Figures 5, 7 and8).

Geodesic reconstruction by dilation

We want to reconstruct a structure in an image I, e. g. a vessel, and possibly only this structure. If the structure is bright on a darker background, a possible strategy is to start from a marker image R representing a small part of the structure, and dilate it recursively under the constraint that the produced image remains smaller than I. By doing so, we hope to recover the bright structure as it is in I, while the background should remain relatively flat and smaller than it is in I.

More formally, given a mask image I and a marker image R, both N × N, the geodesic reconstruction by dilation consists in building recursively a sequence of images (I n ) n≥0 ,

I 0 = R ∧ I and I n+1 = δ (I n ) ∧ I,
and finally take the sup Ĩ = n≥0 I n . Since in our case W ii = 1, the dilation is extensive: δ (I)(u i ) ≥ I(u i ), and the sequence (I n ) n≥0 is increasing. Furthermore, I n ≤ I for any n, hence the sequence converges to its maximal element I n max after a finite number n max of iterations, producing the final reconstruction Ĩ = I n max .

Figure 10 shows an example of such a reconstruction, in which the mask image I is an original image of vessels, and the markers have been defined manually to match bright regions in the vessels. The structuring elements (or equivalently the adjacency matrix W ) was calculated on the tensor field T σ associated with the original image, with σ = 3.

Ideas

• In the geodesic reconstruction, the mask image could be the anisotropy or trace image as well.

Interestingly, they are independent on the contrast (it should work for bright vessels on dark background as well as for dark vessels on bright background). I have tried to work with the anisotropy but it is quite a noisy image and the reconstruction is not very accurate with respect to the original shape of the vessel -more work needs to be done on that. Experiments with the trace are on going as well.

• In the definition of the structuring elements, more information can be included than the spatial and angular consistency.I have tried to include the trace information, to avoid the existence of a path between pixels inside a vessel and pixels outside a vessel. So far this has not really worked out.

• The good part of working with the tensor field T σ is that it provides positive semi-definite matrices, on which we can test methods to be adapted later on diffusion tensors. However, it raises some issues regarding the analysis of vessel images. First, it contains a scale parameter σ , which means that either a multi-scale approach or an automatic scale definition is required. Second, as said earlier, a structure tensor T σ (u) is merely a weighted sum of (rank 1) positive semi-definite matrices.

It is not clear how this averaging behaves with respect to orientation information, especially in bifurcations.

• It seems that Frangi's vessel enhancement [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], based on the analysis of the Hessian, gives accurate segmentation of vessels. It may be interesting to figure out how our approach can add to Frangi's.

• As suggested in previous discussions, we should compare the results obtained with our definition of the structuring elements, to more classic ones (non-adaptive structuring elements, sets of differently oriented segments as structuring elements...).

Tropical and morphological operators for signal on graphs

This section aims at setting or recalling ( [START_REF] Velasco | Nonlinear Operators on Graphs via Stacks[END_REF]) some results about morphological operators for signals on graphs, defined as max/min-plus convolutions as in Equation 1. These results may help to get an insight on all the operators we will define based on this setting.

Setting and definitions

The first two definitions are taken from [START_REF] Velasco | Nonlinear Operators on Graphs via Stacks[END_REF], and deal with a class of real matrices W .

Definition 1 (Morphological weight matrix). A N × N real matrix W = (w i j ) 1≤i, j≤N is called a morphological weight matrix if -∞ ≤ w i j ≤ 0 for any (i, j), 1 ≤ i, j ≤ N.

Definition 2 (Conservative morphological weight matrix). A morphological weight matrix W is said conservative if w ii = 0 for any i,

1 ≤ i ≤ N.
We consider a weighted and directed graph G = (V, E) containing N nodes, whose N × N adjacency matrix, noted W , is a morphological weight matrix:

-∞ < w i j ≤ 0 if (i, j) ∈ E w i j = -∞ if (i, j) / ∈ E.
Note this can be interpreted as taking the log of weights in [0, 1], where the absence of edge from node i to node j would be represented by w i j = 0.

Definition 3 (Path). We call path from node i to node j in G, a tuple of nodes (k 1 , . . . , k n ) such that k 1 = i, k n = j, and (k p , k p+1 ) ∈ E (or equivalently w k p k p+1 > -∞) for 1 ≤ p ≤ n -1.

We will note Γ i j the set of paths from i to j in G, and Γ (p) i j the set of paths from i to j in G containing at most p nodes, for p ≥ 1.

Definition 4 (Weight of a path). Given the weight matrix W , the weight of a path γ = (k 1 , . . . , k n ) in G, noted ω(γ), is the sum ω(γ) = ∑ n-1 p=1 w k p k p+1 . It is easy to see that the set {ω(γ), γ ∈ Γ i j } has a maximum value whenever Γ i j is non empty. Indeed every w k p k p+1 is non positive, and so is the weight of a path. Hence, for a non-empty set of paths Γ i j , the set {ω(γ), γ ∈ Γ i j } has a least upper bound, noted ω * (Γ i j ). Clearly, Γ i j must contain at least a path without a cycle (that is, in which a node is present at most once). Noting Γ i j the set of paths without a cylce in Γ i j , then Γ i j is non empty and we notice that ω * (Γ i j ) = ω * (Γ i j ). This is because one can associate to any path with a cycle a path without a cycle, and the latter has a larger weight. Since Γ i j is finite (it contains no more than (N -2)! paths), ω * (Γ i j ) is a maximum, achieved by a at least a path γ * ∈ Γ i j . Since Γ i j ⊂ Γ i j , ω * (Γ i j ) is also a maximum. It follows that there is a maximal weight for any non empty set of paths, which allows Definition 5.

Definition 5 (Maximal weight, maximal path). We call maximal weight for a non-empty set of paths Γ the number ω * (Γ) := max γ∈Γ ω(γ).

Then a maximal path in Γ is a path γ * ∈ Γ such that ω(γ * ) = ω * (Γ), i.e. achieving the maximal weight among the paths in Γ.

For completeness, we will use the convention ω * ( / 0) = -∞. Now, let x = (x 1 , . . . , x N ) a signal on G, with values in the lattice L = ([0, 1], ≤). We define its dilation and erosion induced by W as ∀i ∈ {1, . . . , N}, δ W (x) i = N j=1 x j + w i j and ε W (x) i = N j=1

x jw ji .

(

) 2
As max-plus products, the dilation δ W and erosion ε W may be noted:

δ W (x) = W x and ε W (x) = 1 -W T (1 -x) .
Similarly, we will note W 2 the max-plus product of W by itself, that is

(W 2 ) i j = N k=1 w ik + w k j ,
and more generally W p the p-th power of W in the max-plus sense, for any integer p ≥ 0 (for p = 0, W p is the max-plus identity matrix, for which w ii = 0 and w i j = -∞ for i = j). Finally, since -∞ is the neutral element for the max operator and is absorbant for the sum, only the terms for which w i j > -∞ count in the previous definition of dilation (Eq. 2), that is to say only the nodes j for which (i, j) ∈ E -the neighbours of i in G. The same remark holds for the erosion, replacing max by min and switching i and j. We will note N i = j ∈ {1, . . . , N}, (i, j) ∈ E the set of neighbours of i in G. Then as we just noticed

δ W (x) i = j∈N i x j + w i j and ε W (x) i = j,i∈N j x j -w ji .
In all the remaining, we will always assume i ∈ N i , with w ii = 0. In other words, we assume W is a conservative morphological weight matrix (See Definitions 1 and 2).

Properties of W p

Proposition 1. Let W a conservative morphological weight matrix and G = (V, E) the associated graph. Then for any p ∈ N, 1 ≤ i, j ≤ N, 1. (W p ) i j > -∞ if and only if there is at least a path from node i to node j in G containing at most p + 1 node(s) (i.e. iff Γ (p+1) i j is non empty)

2. (W p ) i j is the maximal weight of the set of paths from node i to node j containing at most p + 1 node(s):

(W p ) i j = ω * (Γ (p+1) i j
).

At this point one should remember that the w i j are non positive, w i j = 0 modeling the strongest possible effect of node i on node j and w i j = -∞ meaning no effect at all. Therefore, a large weight ω(γ) for path γ ∈ Γ i j represents a strong link from node i to node j.

Proof. Note that once point 1 is proven2 , we only need to prove point 2 in the case where (W p ) i j > -∞.

Indeed if (W p ) i j = -∞ then, according to point 1, Γ (p+1) i j = / 0 and by convention ω * (Γ (p+1) i j ) = -∞ = (W p ) i j .
The case p = 0 is clear although not very informative. There is a path from i to j containing at most p + 1 = 1 node if and only if i = j, and if so there is only one path of weight w ii = 0. Recall that W 0 is the max-plus identity matrix, for which indeed

(W 0 ) i j = -∞ if i = j and (W 0 ) ii = 0.
We assume the proposition true for a given p ≥ 0. By definition

(W p+1 ) i j = (W p W ) i j = N k=1 (W p ) ik + w k j .
Hence (W p+1 ) i j > -∞ if and only if, for at least one k, (W p ) ik > -∞ and w k j > -∞, which means there is a path from i to k containing at most p + 1 nodes, and an edge from k to j, and therefore at least one path from i to j containing p + 2 nodes. This shows point 1.

Now we suppose (W p+1 ) i j > -∞. By assumption (W p ) ik is optimal for any k reachable from i along a path containing at most p + 1 nodes. Then the optimality of N k=1 (W p ) ik + w k j among the paths of

Γ (p+2) i j
is a consequence of a classical dynamic programming argument.

Remark Interestingly, W p is also a conservative morphological weight matrix. Its corresponding graph, that we note G p , can be seen as the original graph G to which a direct edge has been added between nodes i and j whenever there is a path containing at most p + 1 nodes from i to j in G. The weight associated with this new edge is

(W p ) i j = ω * (Γ (p+1) i j
), the maximal weight for the paths from i to j containing at most p + 1 nodes.

Corollary 1. Let W a N × N conservative morphological weight matrix. Then W N = W N-1 and therefore there exists an integer p max = min{p ∈ N,W p = W p+1 }.

Proof. This result is a direct consequence of the interpretation of (W N ) i j as the length of the largest path from i to j contanining at most N + 1 nodes. Since the graph contains N nodes, the latter path actually contains at most N nodes and therefore (W N ) i j = (W N-1 ) i j . This shows that the set of integers p such that W p = W p+1 is not empty and as such has a smallest element p max .

We will note W ∞ = W p max = W N-1 . It is clear from Proposition 1 that (W ∞ ) i j > -∞ if and only if there is a path from i to j, and (W ∞ ) i j is the maximal weight for the paths from i to j: (W ∞ ) i j = ω * (Γ i j ).

Corollary 2. Let W a conservative morphological weight matrix, and p ∈ N. Then (W T ) p = (W p ) T .

Proof. Again, this is clear from the graph interpretation. If G = (V, E) is the graph associated with W , we note G = (V, E ) the graph associated with W T . G has the same nodes as G but its edges are reversed with respect to G: (i, j) ∈ E ⇐⇒ ( j, i) ∈ E , and the associated weight is the same. Threfore, since (W T ) p i j is the maximal weight for the paths in G from i to j containing at most p + 1 nodes, it is also the maximal weight for the paths in G from j to i containing at most p + 1 nodes. Hence (W T ) p i j = (W p ) ji , and this holds for any (i, j).

Note that the result can also be shown by induction exactly like in the linear matrix power case. Indeed, here + law verify the same sufficient conditions as the × law in the linear case: it is associative and commutative.

Operators and filters

From now on, W is a N × N conservative morphological weight matrix, and the coefficient on the i-th row and j-th column of the matrix W p (noted (W p ) i j in the previous section) will be noted w (p) i j . The graph associated with W is noted G = (V, E) and the one associated with W p is noted G p = (V, E p ).

Iterated erosions and dilations

As it is well known in mathematical morphology, if a dilation δ and an erosion ε form an adjunction, then

δ p = δ • • • • • δ and ε p = ε • • • • • ε are
also adjoint dilation and erosion. The associativity of the max-plus product of matrices yields the following result for the dilation and erosion defined by Eq. 2.

Proposition 2. Let W a conservative morphological weight matrix, and p ∈ N. Then δ p W = δ (W p ) and

ε p W = ε (W p ) .
Proof. For the dilation, it is a direct consequence of the max-plus product associativity. For the erosion, we can use the fact that

(ε W ) p (x) = 1 -(δ W T ) p (1 -x) and that (W T ) p = (W p ) T .
A consequence of this is that we can simply express the iterated erosions and dilations as ∀i ∈ {1, . . . , N},

δ p W (x) i = N j=1 x j + w (p) i j and ε p W (x) i = N j=1 x j -w (p) ji , (3)
where, as shown earlier, w

i j represents the strongest link from i to j among the paths containing at most p + 1 nodes. Another consequence is the semi-group property:

δ (W p ) δ (W q ) = δ (W p+q ) .
Extensivity, anti-extensivity Since W and any of its max-plus powers W p are conservative morphological weight matrices, it follows that δ p W is extensive and ε p W is anti-extensive for any p ≥ 0. This holds for ε

∞ W = ε p max W and δ ∞ W = δ p max W . A consequence is that δ ∞ W ε ∞ W = ε ∞ W and ε ∞ W δ ∞ W = δ ∞ W .
Invariant nodes For a signal x on G and any p ≥ 0 or p = ∞, we note

F p (x) ={i ∈ V, ε p W (x) i = x i } and F p (x) ={i ∈ V, δ p W (x) i = x i
} the sets of nodes invariant under ε p W and δ p W respectively.

Proposition 3. For a signal x and any 0 ≤ p ≤ ∞, the sets F p (x) and F p (x) are non-empty.

Proof. We prove the case of F p (x), the other one is then obtained by duality. The anti-extensivity of ε p W yields ε p W (x) i 0 < x i 0 for any node i 0 that is not in F p (x). From equation 3, we deduce the existence of a node i 1 such that ε p W (x)

i 0 = x i 1 -w (p) i 1 i 0 ≥ x i 1 . Therefore, i 1 = i 0 and x i 0 > x i 1 . If F p (x)
were empty, we could proceed likewise to order the N values of x into a decreasing sequence

x i 0 > x i 1 > • • • > x i N .
But then x i N must be in F p (x) otherwise there would be a node j such that x i N > x j , which would contradict the previous ordering.

Idempotence of ε ∞ W and δ ∞ W This is straightforward from Corollary 1: ε p max W = ε p max +1 W = ε ∞ W .
The erosion and dilation ε ∞ W and δ ∞ W are therefore morphological filters, and this is one of the reasons for studying them.

We will now focus on ε

ε ∞ W (x) i = x j -w (∞) ji
and for any node k on a path from i to j achieving w

(∞) ji , ε ∞ W (x) k = x j -w (∞) jk ε ∞ W (x) i = ε ∞ W (x) k -w (∞) ki .
This means that the tranformed signal y = ε ∞ W (x) is composed of some minima y j 0 = x j 0 , j 0 ∈ F ∞ (x), and increasing sequences y j 0 ≤ y j 1 ≤ • • • ≤ y j k along "strong" paths of G, each step between two consecutive values being exactly y j p+1y j p = w (∞) j p j p+1 = w j p j p+1 .

Proof. Let us start with the first statement. Let i ∈ V be a node of G. Case 1: i ∈ F ∞ (x). Then we already have the result:

ε ∞ W (x) i = x i = x i -w (∞) ii . Case 2: i / ∈ F ∞ (x). Then there is a node j = i such that ε ∞ W (x) i = x j -w (∞)
ji , and we need to show that j ∈ F ∞ (x). Assume the contrary, that j / ∈ F ∞ (x). Then there is a node j = j such that ε ∞ W (x) j = x jw (∞) j j < x j , the last inequality coming from the anti-extensitvity of ε ∞ W and the fact that ε ∞ W (x) j = x j by assumption. Then we get

ε ∞ W (x) i = x j -w (∞) ji > x j -w (∞) j j -w (∞) ji . But w (∞) j j + w (∞)
ji is the weight of a path from i to j in G, and as such it is not larger than w (∞) ji , the largest weight of a path from i to j. So -w

(∞) j j -w (∞) ji ≥ -w (∞)
ji and therefore

ε ∞ W (x) i > x j -w (∞) j i which contradicts the expression of ε ∞ W (x) i as an infimum: ε ∞ W (x) i = N k=1 x k -w (∞) ki .
For the second statement, let k a node on a path from i to j achieving w (∞) ji . By the same argument as just above, if we assume that ε ∞ W (x) k = x jw (∞)

j k = x j -w (∞)
jk for a certain j = j, then we find again that ε ∞ W (x) i > x jw (∞) j i . This proves the first line of the second statement, and the second line is then the consequence of the equality w

(∞) ji = w (∞) jk + w (∞)
ki when k is a node on a path from i to j achieving w (∞)

ji .

An implication of Proposition 4 is that ε ∞ W (x) is fully determined by its values on F ∞ (x), in the following sense. Given a signal y ≥ ε ∞ W (x) such that y coincides with ε

∞ W (x) on F ∞ (x) then ε ∞ W (x) = ε ∞ W (y) and F ∞ (y) = {i, y i = ε ∞ W (x) i } (therefore F ∞ (x) ⊂ F ∞ (y)).
Hence knowing the set of nodes F ∞ (x) and the values of ε ∞ W (x) on that set, it is possible to reconstruct ε ∞ W (x) everywhere, for example by computing ε ∞ W (y) where y = ε ∞ W (x) on F ∞ (x) and y = 1 everywhere else.

Particular cases

1. G is an undirected graph with binary weights: w i j = w ji ∈ {0, -∞}. In this case, ε ∞ W (x) is a signal that is constant on each connected component of G, equal to the minimum of x on that connected component.

2.

G is a graph with binary weights:w i j ∈ {0, -∞}. In this case, ε ∞ W (x) is a signal that is constant on each equivalence class defined by the equivalence relation i j ⇐⇒ w ∞ i j > 0 and w ∞ ji > 0 (that is to say, iff there is a path from i to j and from j to i in G), and its value is equal to the minimum of x on that equivalence class. What is more, if there is a path from a node i in a class C i to a node j in a different class C j , then the constant value on C i is greater than the value on C j .

Openings

W ), p ≥ 0 forms a granulometry: γ (p+1) = δ p γε p ≤ δ p ε p = γ (p) (since γ is anti-extensive). What is more, as we saw earlier,

γ (∞) W = ε ∞ W .
Openings with changing weights Another kind of sequence of openings can be defined from a sequence (W n ) n≥1 of morphological conservative weight matrices. Then starting from a signal x we compute a sequence of signals defined by

x n = γ W n • γ W n-1 • • • • γ W 1 (x).

Application to structure tensors

We work on the same images and use the same definitions of structure tensor and anisotropy as in Section 2.

In what folows T is a n × n image of ellipsoids (or tensors). We will index its entries by a single index. The tensor T (i) is characterized by its coordinates in the image (x i , y i ) its eigenvalues 0 ≤ λ 2 (i) ≤ λ 1 (i) and orhtonormal eigenvectors e 1 (i) and e 2 (i). The anisotropy at pixel i is a

i = λ 1 (i)-λ 2 (i) λ 1 (i)+λ 2 (i) . Finally, we note θ i the angle in (-π 2 , π 2 ] such that e 2 (i) is colinear to [cos(θ i ), sin(θ i )] T .

Binary matrix based on co-circularity

Here we present operators based on a binary conservative morphological weight matrix, computed as explained in Section 2.3 and 2.4. We recall the idea of this computation and provide more details.

Definition of W

The graph G = (V, E) we consider is undericted and contains N = n 2 nodes. The signal x i on node i is the anisotropy of the ellipse at pixel i. There is an edge between i and j iff the point u j = (x j , y j ) is in a square window of fixed size 2k + 1, centered on u i (that is, d i j =||u iu j || ∞ ≤ k), and (u i , e 2 (i)) and (u j , e 2 ( j)) are co-circular up to a certain angular tolerance, with the additional constraint described earlier to avoid ladder configurations -see Figure 9. More precisely, given a point u ∈ R 2 and an angle θ ∈ (-π 2 , π 2 ] representing the orientation of a line passing through u, let f u,θ : R 2 → (-π 2 , π 2 ] be the function associating to each point v ∈ R 2 the orientation f u,θ (v) that makes (u, θ ) and (v, f u,θ (v)) cocircular. Then we define the degree of cocircularity between nodes i and j as

c i j = | cos(θ j -f u i ,θ i (u j ))| = | cos(θ i -f u j ,θ j (u i ))| = c ji .
It is a number in [0, 1], and the closer it is to 1, the "more cocircular" i and j are. We say that i and j are cocircular up to the angular precision α ∈ [0, π 2 ] if c i j ≥ cos(α). The additional constraint to avoid ladder configurations is also controlled by an angular parameter β . For i and j to be neighbours in G, we require that max(b i j , b ji ) ≥ cos(β ), with the following definition of b i j b

i j = | --→ u i u j || --→ u i u j || • e 1 ( j)|.
As shown in Figure 9, this restricts the space of possible neighbours to a kind of double cone. We will therefore refer to this constraint as the "conic constraint". Finally, the adjacency matrix W is defined by w i j = 0 if the two constraints are fullfilled or if i = j, and

w i j = -∞ otherwise. Noting bi j = 1 [cos(β ),1] (max(b i j , b ji )), ci j =1 [cos(α),1] (c i j ) and di j = 1 [0,2k+1] (d i j ),
we can also write for i = j w i j = log( bi j • ci j • di j ).

Hence, here W depends on the variables (x i , y i , θ i ) 1≤i≤N , but not on λ 1 (i) and λ 2 (i). W (as defined in Section 3.3.2) on an anisotropy image. The bright parts (high anisotropy) vanish first in the background, that is to say outside the vessels, but even the vessels eventually disappear. This is due to the fact that the graph G is connected (it has only one connected component). This kind of results motivates the search for an operator that would converge towards an anisotropy image where the brightness is preserved in the vessels but not in the background.

Non binary matrix including anisotropy information

The original motivation for including anisotropy in the definition of W was to apply a new opening on a tensor field resulting from a previous opening.

Definition of W

We use the same notations as in Section 4.1.1. In addition, a i and a j denote the anisotropy at pixels (or nodes) i and j. Then as before we set w ii = 0, and for i = j      wi j = di j • bi j • max c i j • min(a i , a j ), 1a i 1/4 w i j = log( wi j ).

Not only W is now non-binary (a i , a j , c i j ∈ [0, 1]), but it is also asymmetric. Indeed, the idea is that a node with low anisotropy value should have a strong effect on its spatial neighbourhood during the erosion (that is the meaning of max(., 1a i )); and a node which had originally a large anisotropy and is cocircular with other bright nodes, should recover a large value during dilation if it has been darkened during erosion (that is the meaning of max( C i j • min(a i , a j ), .)). The power 1/4 was set empirically to see significant effects of erosion and dilation.

Here W depends on the whole information containes in T : it depends on the variables (x i , y i , θ i ) 1≤i≤N and also on λ 1 (i) and λ 2 (i) through the anisotropy a i . We note g the function that associates W to these variables, that is: W = g(u, θ , a) where u = (x i , y i ) 1≤i≤N , θ = (θ i ) 1≤i≤N and a = (a i ) 1≤i≤N .

Granulometry γ

Openings with changing weights

Starting from T = (u, θ , a, λ 1 + λ 2 ) we define a sequence by a 0 = a, W 0 = g(u, θ , a 0 ), and then

a n+1 = γ (p)
W n (a n ), W n+1 = g(u, θ , a n+1 ).

An example of results is shown on Figure 13.

1

  Lattice on the cone of semi positive definite matrices 1.1 Notations • S n : the set of n × n symmetric real valued matrices • S + n : the set of semi-positive definite (SPD) matrices • S ++ n : the set of positive definite matrices • I n : the identity matrix • tr(M): the trace of matrix M • M T : transposition of matrix M • ||M||: the Euclidean norm of matrix M, for the canonical scalar product A, B = tr(A T B) • ∀λ ∈ R, T λ = {M ∈ S n ,tr(M) = λ }: the set of symmetric matrices with trace λ 1.2 The Loewner ordering is not a lattice ordering Definition and geometric interpretation The Loewner ordering is defined on S n as follows: for any A, B ∈ S n , A ≥ B ⇐⇒ A -B ∈ S + n . If A ∈ S +

1 IFigure 1 :

 11 Figure 1: An example of minimal-volume enclosing ellipsoid (yellow), not comparable (in the Loewner sense) to another enclosing ellipsoid (in purple).

Figure 2 :

 2 Figure 2: An example of minimal-volume enclosing ball (green), not comparable in the inclusion sense to another enclosing ball (red).

Figure 3 :

 3 Figure 3: Illustration of the cone S +n . The gray polygones represent two bases of the cone, B 1 and an arbitrary B λ with λ > 1, intersections of the cone with the hyperplanes T 1 and T λ . The extreme points of the bases are marked by black dots (there are infinitely many of them in reality), and the central dashed line indicates the span of the identity matrix I n , orthogonal to each T µ . On each base B µ , the extreme

Figure 4 :

 4 Figure 4: Examples of translated and reversed cones M -C n for several M. The light gray cone corresponds to the smallest ball (M -C n ) ∩ T 0 containing all the others. Figure extracted from [2].

Figure 5 :

 5 Figure 5: Example of input images.

Figure 6 :

 6 Figure 6: Structure tensors T σ for the bottom left hand image of Figure, with σ = 0.1, σ = 1 and σ = 3.The shape of the represented ellipses show their anisotropy and main direction, whereas their colors must encode their trace (this is not clear yet).

Figure 7 :

 7 Figure 7: Trace images corresponding to the input of Figure 5, and the tensor field T σ with σ = 3. They can be seen as smoothed version of the square norm of the ∇ f .

Figure 8 :Figure 9 :

 89 Figure 8: Anisotropy images corresponding to the input of Figure 5, and the tensor field T σ with σ = 3. The closer I an (u) to 1, the more T σ (u) is anisotropic.

Figure 10 :

 10 Figure 10: Illustration of the geodesic reconstruction by dilation. From left to right, top to bottom: mask image I (in this case, the original image); marker image R (defined manually); I 5 , I 25 , I 35 , I 50 , final reconstruction Ĩ = I 178 ; residual image I -Ĩ.

W

  As a family of openings decreasing with p, (γ

4. 1 . 2

 12 Figure11shows the result of openings γ (p)

Figure 11 :

 11 Figure 11: Example of openings γ (p) W on an anisotropy image, in the case of a binary and symmetric matrix W , computed as explained in Section 4.1.1. From left to right, top to bottom: p = 0 (original anisotropy image), p = 1, 5, 15, 20, 30, 40, 60. The parameters used here are k = 7 (spatial window), α = π 20 (cocircularity precision) and β = π 6 (conic constraint).

Figure 12 :

 12 Figure 12: Example of openings γ (p) W on an anisotropy image, in the case of a non-binary and non symmetric matrix W , computed as explained in Section 4.2. From left to right, top to bottom: p = 0 (original anisotropy image), p = 1, 5, 15, 20, 25, 30, 37. The value p = 37 is the p max for W . The parameters used here are k = 7 (spatial window), α = π 20 (cocircularity precision) and β = π 6 (conic constraint).

Figure 13 :

 13 Figure 13: Sequence of anisotropy images defined as in Section 4.2.3. In this example p = 4, and from left to right, top to bottom: n = 0 (original anisotropy image), n = 1, 2, 3, 5, 7, 10, 13. After n = 13, no significant change occurs.

  ∞ W but what follows holds for δ ∞ W by duality. The following result gives an idea of what ε ∞ W should look like. Proposition 4. Let x a signal on G. For any node i, there is a node j ∈ F ∞ (x) such that

Point 1 can be inferred from a classical graph theory result on the powers (in the usual algebra this time) of the adjaceny matrix, and taking the log. We will nevertheless prove it here too.