Thibaud Briand
email: thibaud.briand@enpc@ens-cachan.fr

Low Memory Image Reconstruction Algorithm from RAW Images

Keywords: Image reconstruction, kernel regression, denoising, demosaicking, registration

In this paper, we present a fast and low memory image reconstruction algorithm from a burst of RAW images. Existing algorithms produce high-quality images but the number of input images is limited by severe computational and memory costs. Our algorithm processes the images sequentially so that the memory cost only depends on the size of the output image. Data are combined using classical kernel regression and blur is removed by applying the inverse of the corresponding asymptotic equivalent filter that we introduce. In addition, we propose an accurate and efficient registration method for mosaicked images. We verify the performance of our algorithm on synthetic and real images. For a large amount of data, the results are similar to slower memory greedy methods.

I. INTRODUCTION

Many computational imaging methods produce high-quality images from a set of lower quality images. They compensate for the shortcomings of the imaging system by an adequate digital processing of the accumulated data. This idea has been sucessfully exploited for increasing resolution [START_REF] Tsai | Multiframe image restoration and registration[END_REF], [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF], denoising [START_REF] Buades | A note on multi-image denoising[END_REF], [START_REF] Liu | Fast burst images denoising[END_REF], reducing blur [START_REF] Delbracio | Removing camera shake via weighted Fourier burst accumulation[END_REF] or extending depth-of-field [START_REF] Mertens | Exposure fusion: A simple and practical alternative to high dynamic range photography[END_REF]. A particularly interesting problem is image reconstruction from a burst of RAW images since it solves simultaneously demosaicking and denoising, and eventually super-resolution, problems [START_REF] Vandewalle | Joint demosaicing and super-resolution imaging from a set of unregistered aliased images[END_REF], [START_REF] Farsiu | Multiframe demosaicing and superresolution of color images[END_REF]. The reconstructed image does not suffer from artefacts introduced by uncontrolled processing such as quantization, JPEG compression or demosaicking. However in practice the use of accumulated data is limited by huge computational and memory costs.

In this work, we focus on the low memory image reconstruction problem from a burst of RAW images. By burst we mean a set of digital images taken from the same camera, in the same state, and quasi instantaneously. The images differ because of small motions of the camera, noise and illumination variations. The scene is assumed to be static. The color intensity values are available through a color filter array (CFA), commonly the Bayer filter array [START_REF] Bayer | Color imaging array[END_REF].

Building an image from a burst is generally decomposed into two main steps: the registration, where the images are expressed in a common system of coordinates, and the combination, where the data are combined to form an image. An accurate subpixel motion estimation for CFA images is necessary to reconstruct a high-quality image. According to [START_REF] Vandewalle | Joint demosaicing and super-resolution imaging from a set of unregistered aliased images[END_REF] it can be done by applying a registration method designed for grayscale images on the luminance part of the luminance and chrominance decomposition of [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF]. However, the proposed Fourier-based method [START_REF] Vandewalle | A frequency domain approach to registration of aliased images with application to superresolution[END_REF] for estimating a rigid transformation is not accurate enough. In addition, a homographic transformation model may be preferable.

In burst denoising methods [START_REF] Buades | A note on multi-image denoising[END_REF], [START_REF] Liu | Fast burst images denoising[END_REF] aligned RGB images are computed from the estimated transformations by interpolation and a denoised image is obtained by averaging. In a demosaicking or super-resolution context [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF], the aliased images cannot be interpolated. A color image, expressed on a regular grid, is generally computed from the irregularly sampled data obtained by gathering all the pixel values in the reference system of coordinates. A band-limited color image can be computed by applying the adaptive weights conjugate gradient Toeplitz (ACT) method [START_REF] Feichtinger | Efficient numerical methods in non-uniform sampling theory[END_REF] independently on the three channels. The method proposed by [START_REF] Farsiu | Multiframe demosaicing and superresolution of color images[END_REF] is based on a maximum a posteriori estimation technique by minimizing a multiterm cost function, which handles the channel correlation. In [START_REF] Vandewalle | Joint demosaicing and super-resolution imaging from a set of unregistered aliased images[END_REF] the luminance and chrominance parts of the reconstructed image are computed using an adaptive normalized convolution [START_REF] Pham | Robust fusion of irregularly sampled data using adaptive normalized convolution[END_REF]. This is a robust variant of normalized convolution [START_REF] Knutsson | Normalized and differential convolution[END_REF] using polynomial bases and taking into account local linear structures. It can also be described as an iterative steering kernel regression method [START_REF] Takeda | Kernel Regression for Image Processing and Reconstruction[END_REF], [START_REF] Milanfar | A Tour of Modern Image Filtering: New Insights and Methods, Both Practical and Theoretical[END_REF].

While registration is performed sequentially on the images, all the above mentioned combination methods use an iterative scheme and require the availability of all the data in memory at once. A fast and low memory method requires a simpler combination. For instance a blurry image can be computed pixel-by-pixel without any iterative scheme using classical kernel regression, which takes into account the spatial distance but not the photometric distance. As shown in [START_REF] Friedman | The elements of statistical learning[END_REF], each pixel value is a local weighted averaging of the data, i.e., a local linear filtering. The equivalent filter depends on the data spatial repartition.

This paper proposes a fast and low memory image reconstruction algorithm from a burst of RAW images. The combination, using classical kernel regression of second order with Gaussian kernel, is split into an accumulation part and an image computation part. Using the particular structure of the linear systems involved for computing the pixel values, the input images are processed sequentially during the accumulation. An asympotic equivalent filter is introduced and its inverse is applied during a sharpening step. In this paper, the reconstructed image has the same size as the input CFA images so that the super-resolution potential of the algorithm is not discussed. The main contributions of this work are:

(1) an accurate registration method for CFA images, (2) the introduction of the asymptotic equivalent filter used to remove blur and (3) a fast and low memory combination method whose performances are similar to more complex and costly methods.

The paper is organized as follows: Section II presents the computation of pixel values using classical kernel regression and introduces the corresponding asymptotic equivalent filter. Our low memory image reconstruction algorithm from a burst of RAW images is presented in Section III and is experimentally validated in Section IV. The work is concluded in Section V.

II. CLASSICAL KERNEL REGRESSION

Let S be a set of irregularly sampled data and x 0 = (x 0 , y 0 ) be a pixel position. Denote by w the two-dimensional Gaussian function of scale σ s > 0, i.e.,

w(x) = 1 2πσ 2 s exp - x 2 2σ 2 s . (1)
Pixel value computation. Within a local neighborhood B(x 0 , r) of radius r = 4σ s around x 0 , the intensity value at position x = (x + x 0 , y + y 0 ) is approximated by a second order polynomial expansion

P (x, v) = v 1 + v 2 x + v 3 y + v 4 x 2 + v 5 xy + v 6 y 2 . (2)
Denote by (I i , x i ) 1≤i≤N the samples of S located in B(x 0 , r) at position x i = (x i + x 0 , y i + y 0 ) and with intensity value I i . Write to simplify w i = w(x i -x 0 ). The parameter v = (v 1 , . . . , v 6 ) T is chosen as the solution to a weighted linear regression on the intensities I i , which amounts to minimizing the energy

E(x 0 , v) = N i=1 w i (I i -P (x i , v)) 2 . ( 3
)
Set X i = (1, x i , y i , x 2 i , x i y i , y 2 i ) T and A i = X i X T i .
The minimization problem can be interpreted as a weighted least squares problem whose solution satisfies the 6-by-6 symmetric linear system Av = b where

A = N i=1 w i A i and b = N i=1 w i I i X i . (4)
In a non-degenerate case A is positive semi-definite and the system admits a unique solution. Finally, the pixel value at position x 0 is given by P

(x 0 , v) = v 1 .
Asymptotic equivalent filter. Consider the ideal case where the data in S correspond to the sampling, with a uniform spatial distribution, of a bounded function I : R 2 → R . As N and r go to infinity, the energy in (3) is approximated, after renormalization by 1 N , by

E(x 0 , v) = R 2 w(x -x 0 ) (I(x) -P (x, v)) 2 dx. ( 5
)
One can easily show that the above asymptotic energy has a unique minimum v. The corresponding polynomial expansion verifies

P (x 0 , v) = (k * I)(x 0 ) (6)
where k is called the asymptotic equivalent filter and is defined by

k(x) = w(x) 2 - x 2 2σ 2 s . (7)
This kernel does not depend on x 0 and its Fourier transform is given by

F(k)(x) = 1 + σ 2 s 2 x 2 exp -σ 2 s x 2 2 . ( 8
)
To summarize, an image computed by classical kernel regression can be asymptotically interpreted as the convolution between the ideal unknown target image and k. The image may seem blurry since k attenuates the high-frequencies. This justifies the application of an enhancement filter reverting the blur.

III. ALGORITHM

From a burst of N im RAW images of size W × H, our image reconstruction algorithm builds an RGB image of size W × H. An additional reference RAW image is used during the preprocessing and registration steps. The main steps are summarized in Fig. 1

(b).

Preprocessing. A variance stabilizing transformation (VST) is applied to the RAW images to approximate a homoscedastic noise. For the channel c, using the classical assumption that the variance V c of the noise is an affine function of the intensity V c (u) = a c u + b c , the VST is given by f c (u) = √ a c u + b c . The noise curve of the reference RAW image is estimated using the Ponomarenko et al. method [START_REF] Colom | Analysis and Extension of the Ponomarenko et al. Method, Estimating a Noise Curve from a Single Image[END_REF], [START_REF] Ponomarenko | An automatic approach to lossy compression of AVIRIS images[END_REF] and the parameters a c and b c are computed by linear least squares.

Then, a multiplicative mean equalization of the reference image's channels is done and its maximal value is arbitrarily set to 255. Finally, multiplicative mean equalizations between all the images and the reference are performed.

Registration. First, the burst and reference CFA images are converted to grayscale images by applying the 11 × 11 filter of [10, Fig. 4(b)]. Then, homographic transformations linking the images to the reference are estimated by applying the inverse compositional algorithm [START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework: Part 2[END_REF] to the corresponding filtered images. More precisely, the registration is done using a robust error function and a multiscale coarse-to-fine approach as described in [START_REF] Sánchez | The Inverse Compositional Algorithm for Parametric Registration[END_REF].

Accumulation. Computing pixel values by classical kernel regression is equivalent to solving a symmetric linear system Av = b. Thanks to the redundancies in A the system is determined by only 21 coefficients (instead of 42). The accumulation step consists in computing and storing the 21 system coefficients for each output pixel. It takes as input the images and their corresponding homographic transformations. Thanks to the summative expressions in (4) the contribution of a sample to the system coefficients can be computed by accumulation. The system coefficients are initialized to 0. Then, for each of the N im CFA images the contributions of its W × H samples are added to the system coefficients. Note that a sample only contributes to pixels of the same channel and at a distance lower than r.

In the following the value σ s = 1/ √ 2 is used. It insures a non-negligible contribution to a sample at distance 2σ s = √ 2, which is approximately the distance separating red (or blue) samples of a CFA image. Each sample contributes to at most 20 surrounding pixels, which represents 420 updates.

Image computation and sharpening. For each pixel of the color regular grid of size W ×H, the symmetric linear system, corresponding to the computed system coefficients, is solved using a Cholesky decomposition. The blur of the resulting RGB image is corrected in the sharpening step. It consists in applying the inverse of the asymptotic equivalent filter k to the RGB image. This is done using the inverse of ( 8) and the discrete cosine transform (DCT) convolution method.

Low memory requirement. Theoretically the accumulation and image computation steps are equivalent to a single combination step taking as input the three sets, one per channel, of irregularly sampled data. The structure of the resulting algorithm is shown in Fig. 1(a). The storage of the data makes this algorithm memory greedy. It requires to store three floating point numbers per sample with a total amount of 3N im W H.

On the contrary the memory cost of the proposed algorithm is constant with the number of images N im since the images and their samples are processed sequentially. The main cost comes from the storage of the system coefficients, which only represents 63W H floating point numbers.

IV. EXPERIMENTS

Our image reconstruction algorithm was compared to an algorithm using the memory greedy structure of Fig. 1(a) where the combination is done using the ACT method [START_REF] Feichtinger | Efficient numerical methods in non-uniform sampling theory[END_REF]. It uses the same preprocessing and registration but requires no sharpening.

When the reference target image is known we define E ref as the root mean square error (RMSE) of the reconstructed image.

We define E set as the root mean square difference (RMSD) between reconstructed images from two separated subsets of the burst. In the case of an ideal denoising from n images we would have

E ref = 10 3 σ √ n and E set = √ 2E ref .
The factor 10 3 is due to the different number of samples available in the channels (two times more in the green channel). Note that before any image comparison a band of 20 boundary pixels is discarded and a mean equalization is performed.

The first goal of the experiments is to show that we are able to provide similar results as the ACT with less computation time. Secondly, we verify that the decay of E ref and E set empirically follows the ideal denoising error decay.

In the figures, we display details on a small relevant zone. For the difference of color images we display the root mean square (RMS) over the channels. The experiments were made using C-language implementations and an Intel(R) Core(TM) i7-7820HQ CPU (2.90GHz).

On synthetic data. We generated a burst of N im = 200 CFA images from the Rubberwhale image [START_REF] Baker | A Database and Evaluation Methodology for Optical Flow[END_REF], which is a color image of size 584 × 388. The images were linked by homographies obtained by randomly moving the four corners of the input image. Each displacement was uniformly distributed in [-3, 3] 2 . After resampling with a high-order interpolation method, the color filter was applied and Gaussian white noise of standard deviation σ = 5 was added. A noisy CFA reference image was also generated.

The registration step was performed in 30s, i.e, 150ms per image. The end-point error was on average 0.00909 pixel. With such a subpixellic precision, the results with or without estimating the transformations are similar. Therefore, we present the results using the estimated transformations.

The reconstructed images using our algorithm and ACT are displayed in Fig. 2(a)-(b). It is hardly possible to detect differences to the naked eye. As shown in Fig. 2(c)-(d), the residuals using both methods are mainly composed of noise. Our result has higher errors at the discontinuities. This effect is less important with ACT since it is Fourier based. As shown in Fig. 2 denoising error. The error E ref using our algorithm is slightly higher than using ACT but the RGB image was computed 8.5 times faster using about Nim 21 10 times less memory. The evolutions of E ref and E set with the number of input images are shown in Fig. 3. In both cases the curve shape is similar to the ideal denoising one.

On real data. We built images from a burst of N im = 200 RAW images and an additional reference RAW image. The images were acquired using an Olympus E-M5 Mark II camera set in sequential mode and then cropped to have a size of 1024 × 1024.

The registration was performed in 167s, i.e., 835ms per image. The combination and sharpening were done in 367s using our algorithm and in 2705s using ACT. The constructed images are shown in Fig. 4(a)-(b) and compared in Fig. 4(c). The result using ACT is slightly better at the discontinuities but our result is computed 7 times faster and required about V. CONCLUSION This work proposes a fast and low memory image reconstruction algorithm from a burst of RAW images. Thanks to a simple combination method and a sequential processing of the input images, the memory cost only depends on the size of the output image. The registration method used for CFA images is proven experimentally to be efficient and accurate. Classical kernel regression has a particular structure that allows for the accumulation of the data and admits an asymptotic equivalent filter, which is successfully inverted to remove the blur. For a large amount of data our algorithm provides similar results as slower and memory greedy methods. The residual noise decreases as expected.

The super-resolution potential will be discussed in a coming paper. This efficient algorithm opens the way to real time image formation.

Fig. 1 .

 1 Fig. 1. Main steps of our image formation method. (a) and (b) are theoretically equivalent.

  (e), both reconstruction errors are close to the ideal (E ref and computation time

Fig. 2 .

 2 Fig. 2. Image reconstruction on synthetic data. The computation times in (e) do not take into account the registration.

Fig. 3 .

 3 Fig.3. Error evolution with the number of input images using our algorithm on synthetic data.

Nim 21 10

 21 times less memory.

Fig. 4 .

 4 Fig. 4. Image reconstruction on real data. (c) is the difference between (a) and (b), and (d) is the difference between reconstructed images built using our algorithm on two separated sets of 100 images.


Richard Alligier

Cyril Allignol

Nicolas Barnier

Nicolas Durand

Ruixin Wang

ICRAT 2018

Keywords: UAS, self-separation, conflict resolution, geometrical algorithm

In this article, we extend the 2D-framework introduced in 2016 to implement a horizontal detect and avoid algorithm for UASs flying in Terminal Control Areas. First, we introduce a model able to detect conflicting trajectories in 3D and select combined horizontal and vertical maneuvers, taking the urgency of the conflict into account.

We use a large data set of recorded real commercial traffic trajectories to evaluate the ability of our improved algorithm to avoid any loss of separation with commercial airliners. We test two different types of UASs, flying at 80 kt or 160 kt, with six different missions: constant heading or turning and leveled, climbing or descending. We consider both heading change and vertical maneuver so that UASs have more freedom to avoid conflicts, but keep the speed constant as they mostly have poor speed-up performance.

The article investigates the influence of the various parameters on the separation achieved and the amount of maneuvers required, especially the strategy used to select the best maneuver among the allowed headings. The analysis of our results shows that, amid two "basic" and "extreme" strategies that respectively favor minimal heading and vertical changes or the expectation to escape the conflict, the combination of both, switching from the first one to the second whenever the time before the conflict falls under a given threshold, gives the best results with very few remaining close encounters, while keeping low the amount and amplitude of maneuvers.

I. INTRODUCTION

Civilian UASs are used for diverse missions (ranging from fire detection and river bed surveillance to small items delivery, agriculture, etc). The demand for opening spaces to them is very high [START_REF] Abiresearch | Small unmanned aerial systems market exceeds US$8.4 billion by 2019, dominated by the commercial sector and driven by commercial applications[END_REF]. Most of Civilian UASs will fly in lower airspace (under FL180) and sometimes close to airports where commercial aviation traffic is dense. Developing new algorithms to help separating UASs from the rest of the traffic is critical for safety reasons. In this context, different approaches are investigated.

In a centralized approach, separation could be entirely managed by Air Navigation Service Providers. Experimental results [START_REF] Dcf | ODREA demonstration report[END_REF] have shown that mixing UASs with conventional traffic creates difficult situations for Air Traffic Controllers because their performances and time responses are different.

In an autonomous approach, the separation task could be shared by both commercial aircraft and UASs. However, in such a context, air traffic controllers would have to validate the maneuvers of commercial aircraft before they are executed.

It seems more reasonable to let UASs deal with their separation with the commercial traffic. This requires that the positions and speeds of surrounding aircraft are available (through ADS-B for example) and that the performances of UASs enable them to react early enough to handle the whole separation task.

In a previous article [START_REF] Allignol | Detect & Avoid, UAV Integration in the Lower Airspace Traffic[END_REF], we showed that a self-separation algorithm used in robotics in the horizontal plane could be adapted to our problem and we checked different scenarios on thousands of real commercial traffic samples to validate our model. In a small percentage of cases, our approach failed, despite a mixed strategy that focused either on keeping the maximal distance with the traffic (Safest) or tried to minimize the separation (Closest).

We then made our approach more realistic by adding noise to the detection signal recorded by the UASs [START_REF] Allignol | Assessing the Robustness of a UAS Detect & Avoid Algorithm[END_REF]. We were able to show that our approach is generally robust to uncertainties.

In this article, we generalize our approach by adding vertical maneuvers to UASs. Our approach is close to Snape et al.'s [START_REF] Snape | Navigating multiple simple-airplanes in 3D workspace[END_REF] who first introduced a geometrical approach for self-separation for aircraft conflict resolution. We model the forbidden zones with cylinders instead of spheres or ellipsoids in order to comply with various separation standards. We also adapt the hybrid strategy presented in our horizontal separation model that switches between the Safest and Closest modes. We also use a fallback strategy to minimize the constraint violation in case the separation standards are not satisfied.

A first detect and avoid algorithm was designed by van den Berg et al. [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] and tested with different speed constraint hypotheses by Durand et al. [START_REF] Durand | Does ATM need centralized coordination? Autonomous conflict resolution analysis in a constrained speed environment[END_REF] in the context of autonomous air conflict resolution. In [START_REF] Allignol | Detect & Avoid, UAV Integration in the Lower Airspace Traffic[END_REF] we tailored van den Berg's geometrical approach to model the performances of UASs and consider specific fallback strategies to handle cases for which the first approach failed to maintain separation.

As in our previous approach:

• The UAS assumes the whole avoidance maneuver.

• UASs considered fly with low speeds compared to surrounding commercial aircraft. The ratio used in this article can go from 1.5 to 5. We focus on the lower airspace where the aircraft speed is theoretically limited at 250 kt, but recorded data show that in practice some aircraft fly much faster (up to 400 kt). More specifically, we consider as in [START_REF] Allignol | Detect & Avoid, UAV Integration in the Lower Airspace Traffic[END_REF] two types of UASs: Fast UASs flying at 160 kt and Slow UASs flying at 80 kt. However, we now also allow UASs to climb at a rate up to 2200 ft min -1 and descend at a rate up to 4200 ft min -1 . • Because most civilian UASs have very poor speed-up performances compared to conventional aircraft, we only consider maneuvers at constant speed for UASs. Our approach could easily be generalized with variable speed, but relaxing this constraint would have little effect on the resolution process with realistic traffic. • Commercial aircraft flying in the lower airspace are generally climbing or descending and their speeds are constantly changing, either increasing when climbing or decreasing when descending, and changing direction as well. This factor has a great influence on the detect and avoid strategy in order to ensure that a reasonable distance to the encountered traffic can be maintained. Using real traffic data is therefore essential to validate a resolution algorithm for such evolving and intricate traffic. 1) Related Works: The Free-flight concept emerged in the early 90s. In Europe, the main idea was then to equip every aircraft with a detect and avoid algorithm able to ensure separation with the rest of the traffic.

Approaches inspired by physical laws were tested. Zeghal used sliding forces to coordinate maneuvers between aircraft [START_REF]A comparison of different approaches based on force fields for coordination among multiple mobiles[END_REF]. Potential or vortex fields [START_REF] Kosecka | 2-1/2 D Conflict Resolution Maneuvers for ATMS[END_REF], as well as a model based on an analogy with electrical particle repulsion [START_REF] Eby | Free flight separation assurance using distributed algorithms[END_REF], were also used. In 2001, we proposed a token allocation strategy combined with an A * algorithm to solve conflicts with realistic maneuvers [START_REF] Granger | Token allocation strategy for free-flight conflict solving[END_REF], [START_REF] Granger | Optimal resolution of en route conflicts[END_REF]. Even if some maneuvers could be simultaneously decided, a complete ranking of aircraft was necessary and finding an optimal ranking has been shown to be problem-dependent [START_REF] Archambault | Scheduling heuristics for on-board sequential air conflict solving[END_REF]. We also tried artificial Neural Networks on the two-aircraft problem [START_REF] Durand | Neural nets trained by genetic algorithms for collision avoidance[END_REF] but they could not be generalized to handle more aircraft. All these approaches have been tested on en-route traffic, mainly with leveled aircraft.

Geometrical algorithms have also been widely studied in robotics [START_REF] Snape | Navigating multiple simple-airplanes in 3D workspace[END_REF], [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF], [START_REF] Hwang | Protocol-based conflict resolution for air traffic control[END_REF], [START_REF] Ny | Geometric programming and mechanism design for air traffic conflict resolution[END_REF]. ORCA, the powerful technique developed by van den Berg et al. [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] to coordinate maneuvers in a distributed fashion, can handle thousands of agents in a small space. It was applied to aircraft by Snape et al. [START_REF] Snape | Navigating multiple simple-airplanes in 3D workspace[END_REF], but the hypotheses of the algorithm require simultaneous vertical and horizontal speed changes. We also tested them [START_REF] Durand | Does ATM need centralized coordination? Autonomous conflict resolution analysis in a constrained speed environment[END_REF] in the horizontal plane with speed constraints and showed that this algorithm is unable to deal with high densities of traffic when the speed norm cannot be changed.

More generally, conflict resolution has been proven to be a highly combinatorial optimization problem [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF]. Most centralized approaches that have been proposed to solve conflicts can be broadly divided into two main categories. The first ones [START_REF] Krella | Arc 2000 scenario (version 4.3)[END_REF]- [START_REF] Hu | Optimal coordinated maneuvers for three dimensional aircraft conflict resolution[END_REF] use greedy sequential algorithms to optimize trajectories one by one after ranking the aircraft (ordering aircraft is however very challenging [START_REF] Archambault | Scheduling heuristics for on-board sequential air conflict solving[END_REF]). The others try to find the global optimum without the need to prioritize aircraft. Among this second category, many models define aircraft trajectories through simple analytic expressions that introduce strong limits on the type of situations that can be dealt with, as the ones described in [START_REF] Oh | Design and analysis of conflict resolution algorithms via positive semidefinite programming [aircraft conflict resolution[END_REF]- [START_REF] Gariel | 3D conflict avoidance under uncertainties[END_REF]. In [START_REF] Durand | Optimal resolution of en route conflicts[END_REF], [START_REF] Allignol | A New Framework for Solving En-Route Conflicts[END_REF], we proposed a model to solve multiple aircraft conflicts based on Metaheuristics (Genetic Algorithm and Tabu Search) using trajectory simulation with uncertainties. However, these works mainly targeted en-route traffic control and used simulated traffic only with the BADA model on real flight plans.

2) Outline: In this article, we focus on a realistic traffic environment. We consider UASs flying in the lower airspace (under FL180) and design various conflict scenarios with real recorded commercial aircraft trajectories in TCAs (Terminal Control Area). The aim of the study is to assess the performances of a "3D -detect and avoid" strategy for UASs to maintain a reasonable horizontal and vertical separation with commercial traffic. We first consider conflicts involving one UAS and one commercial aircraft only, then we introduce multiple aircraft encountering the UAS.

In section II of this article, we detail our geometrical approach inspired by the concept of Velocity Obstacles [START_REF] Fiorini | Motion planning in dynamic environments using velocity obstacles[END_REF] in the case of a single UAS avoiding non cooperative aircraft. We show how our previous 2D algorithm can be adapted to handle arbitrary separation volume in 3D by discretizing the maneuvers space. In section III-A, we describe how the conflict scenarios were built from real traffic data and the hypotheses that were chosen for the UAS. In section III-B, we give some results obtained on the different scenarios and show the influence of the various parameters on the quality of the results in terms of separations achieved and amounts of maneuvers imposed on the UAS. The last section draws some conclusions about the results obtained with the simulations and highlights directions for future work.

II. DETECT AND AVOID MODEL IN 3D

Building upon the 2D horizontal model of conflict resolution between a UAS and intruding aircraft described in [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF], we present a somewhat simpler geometric view which can be naturally extended in 3D. This new model allows to efficiently discriminate conflicting maneuvers among a discrete set of alternatives and find the best one according to the current strategy.

We present as well a new maneuver selection strategy improving on the distance-based one described in [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF] to take into account the time before conflict to switch between the so-called Closest and Safest modes.

As the former one, our new model can be easily generalized to handle multiple aircraft.

A. Single Aircraft

In [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF], we presented an adaptation of the ORCA algorithm [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] tailored to the context of the integration of a UAS with limited speed change ability in a TCA with evolving airliners. Instead of equally dividing among the two aircraft an arbitrary maneuver that avoids the conflict, we constrained the scheme such that only the UAS is assigned a maneuver and the norm of its speed is kept constant. In the following, we resort to the more simple concept of Velocity Obstacles introduced by [START_REF] Fiorini | Motion planning in dynamic environments using velocity obstacles[END_REF] and at the core of the ORCA algorithm [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF].

As depicted in figure 1, let d h be the standard separation in the horizontal plane between UAS A and intruding aircraft B with respective velocities

# » v A and # » v B , #» v r = # » v A -# » v B their
relative speed and τ the look-ahead time. Whereas ORCA represents the future conflicting zone at time τ by a homothetic transformation of center A and ratio 1 τ of the B-centered protection disc, we simply scale the relative velocity #» v r by τ . In the former, a conflict would occur within the τ lookahead time-window whenever the endpoint of #» v r lies in the rounded cone formed by the protection disc and its transform. However, in our equivalent version, we only need to check if the scaled vector τ #» v r crosses the protection disc around B, avoiding the need to check intersections with a composite shape. Indeed, as can be seen on the top of figure 1, the position of a point p along the line defined by A and #» v r can be parameterized by the time t:

p = A + t #» v r
The time-window [t 1 h , t 2 h ] of the horizontal conflict, if any, can therefore be computed by the values of t at the intersection between this line and the protection circle. Thus, a conflict will occur in the horizontal plane within τ iff:

[0, τ ] ∩ [t 1 h , t 2 h ] = ∅ (1)
To extend our approach to the third dimension, the same reasoning is applied in the vertical plane where the protection volume is defined as a minimal vertical separation d v (usually 1000 ft). The 3D protection volume is therefore a cylinder of radius d h and height d v , but any shape for which an intersection function can be provided could be used instead. As shown on figure 1, the intersection time-window

[t 1 v , t 2 v ] in the vertical plane is computed, if any.
In 4D, a conflict will occur iff

[t 1 h , t 2 h ] and [t 1 v , t 2 v ] intersects within [0, τ ]: [0, τ ] ∩ [t 1 h , t 2 h ] ∩ [t 1 v , t 2 v ] = ∅ (2)
To avoid the conflict within τ , a maneuver must then be chosen such that this intersection is empty.

In the following, we note cf (A, B, # » v A , # » v B ) a boolean function corresponding to the negation of inequality (2), i.e. indicating that the UAS A with a speed # » v A is free of conflict until τ with the aircraft B at the speed # » v B .

B. Maneuvers Range

In our former 2D approach [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF], the whole set of constant norm maneuvers could be easily computed as a union of allowed angle ranges for the new velocity # » v A . With the complex geometry of dynamic cylinders, the non-convex surface of all conflict-free combinations of heading and path angle is much more difficult to compute. Therefore, we chose to discretize the set of maneuvers to check, trading complex geometric computations with brute force intersection checks (which could easily be done in parallel with several cores or a GPU).

To assess the performance of this new 3D model over the previous one, the norm of the UAS speed is maintained constant -even if this constraint could easily be relaxed at the price of multiplying the intersection checks described in the previous section. The maneuvers are therefore described by two parameters, the heading θ and the path angle φ. Horizontally, the heading deviation is symmetrically constrained by θ max on both side:

A B D h t 1 h t 2 h # » v A # » v B τ #» v r x y x z A B D h D v t 1 v t 2 v # » v A # » v B τ #» v r
θ ∈ [θ -θ max , θ + θ max ] = Θ
(with θ max = 30 • typically for a 10 s time step) whereas the vertical angle must be constrained both absolutely and relatively around the current path angle to roughly model the capacities of typical medium-size UAS:

φ ∈ [φ abs min , φ abs max ] ∩ [φ -φ rel min , φ + φ rel max ] = Φ In the end, the possible maneuvers is a set PM containing all the speeds obtained after a deviation has been applied to

# » v A : PM = # » v A | # » v A heading ∈ Θ, # » v A path ∈ Φ, # » v A = # » v A
These maneuver alternatives are then approximated by discretizing the aforementioned ranges with angular steps η h and η v , leading to a worst-case O( Θ η h × Φ ηv ) intersections computation (e.g. around 1000 with typical 60 • and 20 • of horizontal and vertical angular openings respectively, and η h = η v = 1 • ). Note that once an angular deviation is fixed on one of the parameters φ or θ , if the intersection with the lookahead time window [0, τ ] is empty (e.g. violating inequation 1 for the horizontal plane), then the intersection processing in the other dimension can be spared (as inequation 2 will be necessarily violated as well).

C. Heading Change Strategies

With the previous sections, we can easily compute the set of conflict-free maneuvers noted CF . More specifically, this set can be defined by using the boolean function cf defined in Section II-A. This function indicates that the aircraft B and the UAS A are free of conflict within τ minutes. With this function, the mathematical definition of CF is easy:

CF = { # » v A ∈ PM | cf (A, B, # » v A , # » v B )}
We need to chose one maneuver in CF . To do so, there are a number of strategies. As observed in our previous study ( [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF]), one strategy alone does not offer a good compromise between conflicts avoidance and maneuver quantity (cf. section III-B). But a hybrid strategy able to switch between minimal maneuvering and fastest escaping when the conflict becomes dire seems to be the most promising. To determine the current strategy, we present a new criterion that takes the dynamic of the conflict into account instead of the distance alone.

If the current speed does not generate a conflict within τ minutes (i.e. # » v A ∈ CF ) then the UAS can keep its current speed. If the current speed generates a conflict at a date t conflict τ (i.e. # » v A / ∈ CF ), then we have to choose a maneuver. If we can escape the conflict in one maneuver (i.e. CF = ∅) then we may have enough time to escape the conflict with a minimal deviation if t conflict is large (Closest). If t conflict is small then the actual conflict is imminent and we chose the fastest escaping maneuver among the conflict-free maneuvers (Safest). If the conflict cannot be escaped in one maneuver (i.e. CF = ∅) then the situation is dire and we opt for the fastest escaping maneuver (Fallback).

The difficulty in the conception of this hybrid strategy is to build a balanced criteria used to switch from strategy to one another. As a consequence this hybrid strategy is parameterized by a parameter γ. The value of γτ is the threshold that decides if the conflict is imminent (i.e. t conflict < γτ ). Therefore, this parameter allows to continuously slide from a pure Closest strategy for γ = 0 to a pure Safest strategy for γ = 1. The influence of parameter γ is discussed in section III-B.

To summarize, the hybrid strategy choosing # » v A is parameterized by γ ∈ R which varies from 0 to 1:

# » v A =        # » v A if # » v A ∈ CF Closest if γτ t conflict τ and CF = ∅ Safest if t conflict < γτ and CF = ∅ Fallback otherwise (i.e. CF = ∅)
This hybrid strategy uses different strategies that are described in the following sub-sections:

1) Closest: In order to minimize the maneuver burden assigned to the UAS when the threat is low, the conflict-free maneuver closest to the UAS current velocity is chosen. In 2D, it was sufficient to minimize the heading deviation, but in 3D, the euclidean distance between current and new speed must be minimized:

Closest: arg min # » v A ∈CF # » v A -# » v A (3)
where # » v A is defined by a feasible conflict-free maneuver in CF , i.e. for which inequation 2 is false. As the only global minimum of this distance trivially is the current speed # » v A (or any reference velocity), the search for the best conflict-free maneuver could consider them by increasing distance and select the first admissible one encountered to speed up the sampling process.

A B (A -# » v B ) # » v A # » v B #» v r Safest #» v r # » v A escape angle
Note that this criterion can be used as well if the speed norm is allowed to change. Other criterion, for example biased towards horizontal or vertical maneuvers to emphasize the least costly change, or taking into account the flight plan of the UAS with a reference velocity, could be used in the same fashion.

2) Safest: When the conflict is imminent or when unexpected changes of the intruding aircraft trajectory are worsening the threat, the Closest strategy is not always able to solve the conflict. More radical maneuvers should then be selected to maximize the expectation of escaping the conflict. In [START_REF] Allignol | Integration of UAS in Terminal Control Area[END_REF], we proposed to choose the most "robust" maneuver defined as the maneuver authorized by the turn rate of the UAS that is the closest to the middle angle of the largest conflict-free angle range (without taking maneuverability into account).

With discretized 3D maneuvers, we could define a similar criterion based on the conflict-free maneuvers graph regardless of the UAS maneuverability, where nodes are maneuvers defined by (θ , φ ) ∈ [-π, π]×[-π/2, π/2] and edges connect nodes that exactly differ of one discretization step on one component, i.e. the spherical grid graph that discretizes (θ , φ ) space. We could then select the center of its largest connected component as the analog of the most robust maneuver in the horizontal plane.

However, it may prove too costly to compute the set of all conflict-free maneuver over [-π, π] × [-π/2, π/2] and we resort to a much more simple criterion to avoid the worst threats: the escape angle between #» v r and # » AB, which should be chosen to avoid the conflict. Therefore, the Safest maneuver can be selected as the one that minimizes the cosine of the escape angle, which can be computed by:

Safest: arg min # » v A ∈CF #» v r • # » AB #» v r # » AB (4)
This is illustrated in figure 2 with an UAS at position A and an intruding aircraft at position B.

Similarly to the Closest strategy based on a distance from a reference vector, the search for the optimal maneuver w.r.t. the Safest one could be sped up by starting with potential minima, regardless of the conflict constraints. However, as the 3D dot product involves two cosine functions, there can be up to four local minima from which the search should start to find the global optimum.

Moreover, these schemes can easily be generalized for a UAS with variable speed and the graph-based one naturally extends to the multiple intruding aircraft case.

3) Fallback: If the conflict-free maneuvers set CF is empty, it means that no turning angle nor vertical move can guarantee the horizontal and vertical separation distance d h and d v for the next τ minutes.

In such a case, the hybrid strategy resorts to a fallback strategy consisting in choosing the maneuver among PM which gives the biggest escape angle between #» v r and # » AB. Thus, the maximal escape angle among the possible maneuvers PM is chosen:

Fallback: arg min # » v A ∈PM #» v r • # » AB #» v r # » AB (5)

D. Multiple Aircraft

In real traffic situations, conflicts involving more than two aircraft often occur and cannot always be solved pairwise in sequence, but should be globally handled. Fortunately, our geometrical model can be easily extended to take into account several aircraft. For n intruding aircraft B i , we just need to compute the intersection of all individual

CF i = { # » v A ∈ PM | cf (A, B i , # » v A , # » v Bi )
} to obtain the global conflict-free maneuvers set CF:

CF = n i=1 CF i
For the Closest strategy, the selection of the best maneuver is then identical to the single aircraft case with this new definition of CF . It is not so straightforward for the Safest strategy as the criterion depends on #» v r and # » AB, so will change for each intruding aircraft B i considered.

Various techniques could be used to solve this "multiobjective" optimization problem. For example the best conflict-free maneuver for the aircraft involved in the most dire conflict could be selected. The "direness" of a conflict could be estimated by the time before collision for closing trajectories:

# » AB i 2 #» v r • # » AB i
or by the size of the conflict-free maneuvers set CF for each intruding aircraft B i .

III. EXPERIMENTS

Our geometrical model has been implemented and thoroughly tested for the single aircraft case on recorded real TCA traffic enhanced by injecting various conflicting UASs scenarios. The following section describes the data and scenarios generation, and the next one reports the results and analysis of our experiments.

A. Description of simulations

Simulations have been carried out on recorded traffic data from a terminal maneuvering area in the south-west of the French airspace on 2013/09/14. The use of recorded tracks in TCA provides a most realistic picture of the trajectories that UASs might be confronted to in operations, as many UAS missions are flown in the lower airspace. Besides, such scenarios put the detect and avoid algorithm to stress because the trajectories of surrounding aircraft are convoluted and exhibit more complex 3D profiles than en-route traffic. The traffic sample contains 478 flights, evenly shared between departures and arrivals.

The simulated trajectories for UASs are derived from recorded tracks of aircraft with similar performances. Two types of UAS have been tested in our experiments, representing a small UAS with relatively poor flying performance and a more maneuverable UAS respectively. The specifications used in the simulations are described in table I. For each UAS model, six missions are implemented, corresponding to all possible combinations of two horizontal profiles (constant heading or circle around a fixed point) and three vertical profiles (constant altitude, climb or descent). Thus, for a given set of parameters and a given UAS, 2868 simulation scenarios are simulated. Each of these scenarios is built in such a way that, if no maneuver is issued at all, a collision would occur between aircraft and UAS.

A fast-time simulator enables us to play the trajectories (both recorded and built) and to modify them by sending maneuver messages consisting of a heading change and a turn rate for horizontal maneuvers, and a flight level change and a climb (or descent) rate for vertical maneuvers. Those messages are sent to UASs only, other aircraft are left unchanged.

To evaluate the outcomes of our experiments, we use several separation volumes around the aircraft corresponding to different objectives of the detect and avoid function, as illustrated in figure 3. The smallest one (in red) is the collision volume, in which the UAS should never enter. It is defined by ICAO [START_REF][END_REF] as a cylinder of radius 500 ft and height 200 ft centered on the aircraft position. The second volume (in green) around the collision volume corresponds to the remain well clear function as defined in [START_REF] Cook | Defining well clear for unmanned aircraft systems[END_REF]. This volume is also an aircraft-centered cylinder, with radius 4000 ft and height 900 ft. Finally, we use a third cylindrical volume (in blue) corresponding to the target separation of our algorithm. The dimensions of this target volume vary, in this study, from the dimensions of the well We also record the number of maneuvers given by the algorithm, as well as the cumulated quantity of deviation from original flight direction (both in horizontal and vertical planes).

For each set of parameters (i.e. γ and the separation target) and for each UAS and mission, we perform three different simulations, allowing only horizontal maneuvers, only vertical maneuvers or both horizontal and vertical maneuvers (these are referred to as H, V and 3D respectively in the following results).

B. Results and analysis

In this section, we analyze the outcomes of many simulations: two UASs, eight values for γ and four target separation volumes have been tested. For each set of parameters, all 2868 traffic scenarios are played, leading to a total of more than 180 000 single simulations.

1) Impact of γ on the maneuvers: In order to quantify the impact of γ on the maneuvers, we have summed all the deviations (both vertical and horizontal) given to the UAS. This sum is then divided by the number of scenarios giving us the average maneuver quantity. This value is depicted in figure 4, which also depicts the average number of maneuvers. γ is increased from 0 to 1, which corresponds to a shift from a Closest-only strategy to a Safest-only strategy.

We see that the average maneuver quantity increases slowly with γ, with a higher rate low γ values and for the 3D maneuvers. This is expected, as the Safest maneuvers tend to produce be more intense than Closest maneuvers.

Interestingly, the average number of maneuvers decreases with γ. This means that, even if the Safest strategy gives larger deviations than the Closest strategy, the maneuvers tend to be more efficient at maintaining separation, thus needing less further corrections in the following time steps. With the closest strategy, however, minimal deviations are applied that will more likely need to be adapted.

2) Separation and Deviation: The main objective of the detect and avoid algorithm is to maintain a given separation distance between the UAS and surrounding aircraft. Besides, the deviation of the UAS from its mission also matters significantly for the sake of efficiency. However, these two criterion are often antagonistic.

Our hybrid method is parameterized by γ and the separation volume. These parameters must be tuned to find the best combination. As two opposed criterion must be taken into account, it is not possible to find a parameter combination that is best for both criterion. Actually, we observe that several sets of parameters lead to solutions that are not dominated for both criteria.

These non dominated criteria constitute a Pareto frontier, which is depicted in figures 5 and 6. More specifically, figure 5 depicts the percentage of success of "remain well clear" with respect to the quantity of maneuver that is necessary to achieve this result.

On figure 5, we can see that we can remain well clear in 95 % of the scenarios if the parameters are well chosen. This can be achieved with a maneuver quantity inferior to 40 • except for UAS1 when only heading changes are authorized.

On figure 6, we can see that we can avoid the collision volume in 99.8 % of the scenarios if the parameters are well chosen. This can be achieved with a maneuver quantity inferior to 40 • except, again, for UAS1 when only heading changes are authorized. This corresponds to only a handful of situations (one or even zero depending on the parameters with 3D maneuvers, compared to the 2868 scenarios), that need further analyze to adapt our algorithm. For both figures, the best results were obtained with UAS2, the fastest UAS. This can be explained by the fact that a deviation on the heading or the path have a different impact on the new relative speed depending on the speed of the UAS. The higher the UAS speed, the higher is the impact on the relative speed. Also, the vertical maneuvers are more efficient to solve the conflicts than the horizontal maneuvers. This is especially true for UAS2. This can be explained by the fact that the vertical norm is smaller than the horizontal one. In terms of separation success percentage, the 3D maneuvers are the most efficient. In that respect, the use of vertical maneuvers seems mandatory for UAS1, the slowest UAS.

3D -80kt H -80kt V -80kt 3D -160kt H -160kt V -160kt
3D -80kt H -80kt V -80kt 3D -160kt H -160kt V -160kt
After analyzing our results with respect to both separation and efficiency objectives, we found that targeting a separation volume 50 % larger than the well clear volume and setting γ to 0.2 was a particularly interesting setting. This enables for an achievement of well clear separation of 95 % for UAS1 and almost 99 % for UAS2, while keeping the deviation to a very low level. Also, γ = 0.2 corresponds, in our settings, to triggering a maneuver only one minute before the potential loss of separation.

IV. CONCLUSION AND FUTURE WORK

In this study, we proposed a three-dimensional detect and avoid algorithm for UAS integration in air traffic, based on previous work on a 2D approach. This algorithm aims at providing separation between a UAS and commercial traffic at departure or arrival, but proved to be efficient for short-term collision avoidance as well.

The possible maneuvers for the UAS, either horizontal, vertical or both, are sampled and each one of them is geometrically confronted to the constraint induced by the intruder aircraft, in order to obtain a set of feasible trajectories. Then, two strategies have been defined and then combined to choose the best maneuver. The first strategy (Closest) follows the mission as close as possible, while the second one (Safest) tries to escape from the conflicting area in the fastest possible way. The algorithm is set up to switch between strategies depending on the seriousness of the situation.

This approach has been intensively tested in simulations involving real traffic recorded in a French TCA. With the proper set of parameters, the required separation was provided with a very high success rate, and collision avoidance was ensured at almost 100 %. The underlying trajectory deviations are kept at an operationally acceptable level, and the algorithm can be set up so that the maneuvers are triggered only one minute before the predicted loss of separation.

In further works, we plan to increase the possibilities of the UAS by authorizing speed changes, which might help in finding solutions with less deviations from the mission.

One of the pitfalls of our method is that it only takes into account the current state, so that any further change in the aircraft state could break the resolution. In order to improve the robustness of the maneuvers, we plan to try and anticipate better, both on the aircraft intentions and the UAS capabilities.

Knowing the past positions of the aircraft, it is possible to build a short-term predicted trajectory, based on the analysis of the derivatives of its speed and turn angle. For example, the beginning or the end of a turn, a climb or a descent could be inferred. Particular care would have to be taken during the calibration phase, especially when choosing the number of past states to consider: too much states would create some latency in the prediction, whereas too few states would yield unreliable ones.

If the aircraft trajectory could be predicted this way, then it becomes particularly interesting to anticipate several maneuvers for the UAS. This could be planned optimally with an A * or Dijkstra algorithm, using the geometrical algorithm at each step to prune the search tree or validate the existence of partial solutions, at the cost of a significantly longer computation. It could also be performed geometrically by an approximation of a few maneuvers aggregated into a single one.

Fig. 1 .

 1 Fig. 1. Horizontal (above) and vertical (below) projections of the geometric model in 4D: a conflict will occur within time τ if and only if the relative speed #» vr scaled by τ crosses the protection volume.

Fig. 2 .

 2 Fig. 2. Safest maneuver and underlying escape angle.

Fig. 4 .

 4 Fig. 4. Average number of maneuvers and average maneuvering quantity w.r.t. hybridization parameter γ.

Fig. 5 .

 5 Fig. 5. Percentage of success of "remain well clear" w.r.t. the necessary quantity of maneuver.

Fig. 6 .

 6 Fig. 6. Percentage of success of collision avoidance w.r.t. the necessary quantity of maneuver.

TABLE I UAS

 I SPECIFICATIONS FOR SIMULATION.

  ST) VOLUMES USED IN SIMULATION. ST1 CORRESPONDS TO THE WELL CLEAR VOLUME AND ST4 CORRESPONDS TO ATC SEPARATION IN TCA. ST2 AND ST3 ARE INTERMEDIATE VALUES.Fig. 3. Description of the aircraft-centered separation volumes: in red the collision volume, in green the well clear volume and in blue the target separation volume. ST h and STv are the horizontal and vertical separation targets respectively.clear volume up to 3 NM and 2000 ft (which corresponds to ATC separation in TCA in France). The various values used for the target separation are provided in table II. For each scenario where the algorithm could not avoid it, we keep track the fact that the UAS entered the well clear volume or the collision volume.

		ST1	ST2	ST3	ST4
	Radius 4000 ft	1 NM	2 NM	3 NM
	Height	900 ft 1000 ft 1400 ft 2000 ft
			TABLE II
	500 ft SEPARATION TARGETS (200 ft ST v ST h	900 ft
				4000 ft