Raúl Sáez

Homeyra Khaledian
email: homeyra.khaledian@upc.edu

Xavier Prats
email: xavier.prats@upc.edu

Andréas Guitart
email: andreas.guitart@enac.fr

Daniel Delahaye
email: daniel.delahaye@recherche.enac.fr

Eric Feron
email: eric.feron@aerospace.gatech.edu

A Fast and Flexible Emergency Trajectory Generator Enhancing Emergency Geometric Planning with Aircraft Dynamics

Keywords: Aircraft emergency planning, Safety, Trajectory prediction, Dubins paths, Path generation, Geometric route planning

We present an automated emergency trajectory generator to compute the best emergency trajectory for a given landing site. A combination of the optimized version of the rapidly exploring random tree algorithm and Dubins paths is used to compute a path connecting the aircraft position with the landing site, which avoids the obstacles in the way. Then, this path is used as input for a trajectory prediction (TP) algorithm, which computes a four-dimensional trajectory by taking into account the current aircraft performance and weather. The set of vertical profiles considered in the TP has been designed in order to cover the widest possible range of emergency situations. Moreover, the aircraft intents considered in these profiles are chosen by taking into account the operational requirements of the air traffic operation system and the input of the flight crew. Among these profiles, two have been tested during the study, to verify the result of the proposed algorithm and its computing time, which is one of the main success criteria. This concept is expected to be part of an advanced flight management system on-board function to help the pilot take efficient and effective decisions in emergency situations and adverse conditions.

I. INTRODUCTION

Nowadays, emergency trajectories for airliners in a degraded flyability mode do not exist except for engine loss situations in SIDs (standard instrumental departures), for which airliners have to specifically design the corresponding flight procedure. The current process of defining emergency trajectories and landing sites remains completely manual and fully relies on the capabilities of the pilot for situation analysis. In emergency situations, an automated support to the pilot could suppose a clear advantage by providing a trajectory to safely bring the aircraft from the location where the emergency takes place until a safe and appropriate landing site.

We can observe in history several cases of successful landings with very degraded aircraft capabilities-e.g., the U.S. flight 1549, ditching in the Hudson river after a double engine failure caused by bird strike-in which, thanks to a good situation analysis of the aircraft crew, the safety of the passengers and integrity of the aircraft structure were preserved. However, in other situations the consequences were fatal and, even if the outcome is positive in some cases, having an automated emergency trajectory generator could help to ensure avoiding compromising the safety of the operation in abnormal situations.

In this work, we present an emergency-trajectory-generation function for commercial aircraft allowing a safer return to ground when normal operation is interrupted, in the context of future automated operations with reduced crew. Emergency trajectories are generated by using a combination of the optimized version of the rapidly exploring random tree algorithm (RRT*) and Dubins paths, which are used to generate the route (i.e., lateral path), avoiding the obstacles present from the aircraft position to the landing site. Then, a fourdimensional (4D) trajectory is generated with a trajectory predictor (TP), where both the aircraft performance and the most updated weather forecast available are taken into account. A comprehensive set of vertical profiles-used when generating the 4D trajectory-have been identified in this work, designed to be as generic as possible in order to cover the widest possible range of emergency situations. Ultimately, the outcome of this work is to bring a support to the pilot-in both flight management and decision making-in emergency situations with degraded aircraft capabilities, where emergency trajectories would be injected within the flight management system (FMS). In this work, we focus only on the generation of the trajectory, assuming a proper landing site has been previously identified. In future work, an insight on the landing site selection process will be given, and a wider range of situations will be investigated in order to assess the feasibility of our concept.

Similar works have dealt with the generation of emergency trajectories for civil aircraft. For instance, in [START_REF] Atkins | Emergency Flight Planning Applied to Total Loss of Thrust[END_REF], the authors presented an automatic post-failure flight planning to enable safe emergency landings. More recently, in [START_REF] Atkins | Emergency landing automation aids: an evaluation inspired by US airways flight 1549[END_REF], the authors also investigated the problem of the generation of a safe landing trajectory for an airplane with structural damage to its wing and flying very close to local terrain. In [START_REF] Di Donato | An Off-Runway Emergency Landing Aid for a Small Aircraft Experiencing Loss of Thrust[END_REF], an emergency planning technique focusing on small aircraft was presented. However, in all these works, only one type of emergency was considered (i.e., loss of thrust) and no obstacles were taken into account when generating the trajectory. Similarly, in [START_REF] Paul | Flight Trajectory Planning for Fixed-Wing Aircraft in Loss of Thrust Emergencies[END_REF], trajectories were generated for a loss of thrust emergency, but only considering purely geometric criteria for the generation of the trajectory, thus, neglecting aircraft dynamics. Finally, other remarkable works considering path planning under an emergency can be found in [START_REF] Meuleau | An emergency landing planner for damaged aircraft[END_REF], where road-maps are used to generate the path, or in [START_REF] Miwa | Real-Time Trajectory Generation Applicable to Emergency Landing Approach[END_REF], where a real-time flight trajectory generator was presented. Furthermore, Garmin Ltd. developed the Autoland system [START_REF]Autonomí™Autonomous Safety-enhancing Technologies[END_REF], exclusively dedicated to general aviation. It consists in an autopilot-based functionality which, in case of an emergency where the pilot is unable to fly, determines the most suitable airport and runway and leads the aircraft to that runway, by taking into account weather, terrain, obstacles and aircraft performance capabilities.

To the best of the authors' knowledge, although there are many works dealing with the generation of emergency trajectories, they are mostly based on the geometric planning of these trajectories. Furthermore, most of the existing works do not take into account the operational requirements of the air traffic operations system (or the flight crew point of view), or they just compute trajectories for one specific emergency. In some cases, the existing solutions are only applied to general aviation. Finally, a lot of research is based on very simple scenarios, where no obstacles are found in the vicinity. The methodology proposed in this paper tries to bridge all these gaps, producing operationally sound trajectories tailored to specific aircraft (degraded) performance, considering weather and obstacle environments and covering a wide range of emergency or urgency situations.

II. BACKGROUND

In this Section, we describe the different techniques we used in this work to generate the emergency trajectory. In Subsections II-A and II-B, we describe the methods used to generate the path: the RRT algorithm and Dubins path respectively. Then, in Subsection II-C, we focus on the trajectory prediction problem.

A. Path Planning Algorithm: Rapidly Exploring Random Tree Algorithm

Many methods propose to generate a graph in order to then find the optimal path between two points. The most used methods are sampling-based path planning algorithms, which consist in generating a graph from a sampling [START_REF] Karaman | Sampling-based Algorithms for Optimal Motion Planning[END_REF]. In this paper, we use one of these algorithms to generate the path, the RRT. Before discussing the algorithm, it is essential to define the problem and the main functions of this algorithm.

Let χ = (0, 1) d be the configuration space, where d ∈ N is the space dimension, (d ≥ 2). Let χ obs be an open set, which denotes the obstacle-free space as χ f ree = cl(χ\χ obs ), where cl(χ) denotes the closure of a set χ. The initial condition is denoted by x init ∈ χ f ree and the goal region, χ goal , is an open set of χ f ree .

This algorithm is composed of four main functions. The first one is the sampling function, which generates a sequence of points in χ. The second one is the nearest neighbor function, which determines the vertex that is closest to a point x ∈ χ. The third function is the near vertices function, which returns a set of vertices that are contained in a ball of radius r centered at a point x ∈ χ. Finally, the last function is the collision test function, used to know whether the straight line between two points lies in χ f ree or not.

At the beginning of this algorithm, the graph is composed of only the initial vertex and no edges. Then, at each iteration, a new point x rand ∈ χ f ree is generated. After that, this point is steered towards its nearest point. Given two points x, y ∈ χ, the function Steer : (x, y) → z returns a point z ∈ χ such that z is "closer" to y than x. The point returned by this function is connected to the graph if there is no obstacle between this point and the nearest point (Figure 1(a)). An improvement of such algorithm is the RRT* algorithm. It connects vertex in the same way as RRT. Moreover, at each step it tries to improve the graph around this new vertex. If the addition of the new vertex can reduce the cost of vertices that are within a distance r, a new connection replaces the previous connection (Figure 1(b)). The distance r is computed as a function of V, which is an open set composed of nodes, and a coefficient γ RRT * , whose value depends on the search space. Figure 2 shows, for the same number of iterations, the graph generated by the RRT algorithm (left) and the RRT* algorithm (right). This figure shows that with the phase of tree improvement, the RRT* tree is more ordered than the RRT tree and, therefore, the path computed by the RRT* algorithm is closer to the shortest path. The RRT* is a very efficient algorithm because its time complexity is O(n log n). However, the generated paths are not flyable. Therefore, it is necessary to modify this algorithm in order to take into account the minimum turn radius and the flight path angles computed by the path bounds generator. The resulting trajectory will not be completely flyable, but it will ensure that all obstacles are avoided. Then, the motion planner (Section III-D) will take into account the aircraft performance to generate the flyable trajectory that the aircraft will follow.

B. Dubins Path

Dubins path is the shortest curve that connects two points with a constraint on the curvature of the path and with initial and final orientation angles. There are two different types of Dubins curve [START_REF] Shkel | Classification of the Dubins set[END_REF] [START_REF] Satyanarayana | Shortest Dubins Path to a circle[END_REF]. The first one is the Circle-Segment-Circle (CSC), which is composed of a first turn, a segment, and a second turn. The second type of Dubins curve is Circle-Circle-Circle, which is composed of three different turns. In Sections II-B1 and II-B1, we show two examples of Dubins curves, a Left-Segment-Left Dubins curve and a Left-Right-Left Dubins curve respectively.

Let P 1 = (x 1 , x 2 ) and P 2 = (x 2 , y 2 ) be the initial position and the final position respectively. Let χ 1 and χ 2 be their track direction. An associated value of track is the orientation angle. Let α and β be the initial and final orientation angle defined as follows:

       α = mod( π 2 -χ 1 , 2π) β = mod( π 2 -χ 2 , 2π) (1)
1) Left-Segment-Left: Figure 3 depicts a Left-Segment-Left Dubins curve. The curve connects a starting point with an orientation angle α and an end point with an orientation angle β. It is composed of a turn to the left around the first green circle, a segment of length L s and a second turn to the left around the second green circle.

The length of the segment L S and the two angles φ 1 and φ 2 are defined as follows :

L s C 1 C 2 × × × L s π/2 -φ2 α φ1 φ1 φ2 β O D Fig. 3. Left-Segment-Left Dubins Curve F irstT urnAngle = φ 1 -α if α ≤ φ 1 φ 1 + (2π -α) otherwise (2)
L S = (x 2 -r sin β -x 1 + r sin α) 2 + (y 2 + r cos β -y 1 -r cos α) 2 (3)
φ 1 = mod(arctan y 2 + r cos β -y 1 -r cos α x 2 -r sin β -x 1 + r sin α (4)
φ 2 = mod(β -φ 1 , 2π) (5)
2) Left-Right-Left: The three angles, φ 1 , φ 2 and φ 3 , are defined as follows:

φ 1 = mod(-α + arctan y 2 -y 1 -r cos α + r cos β x 2 -x 1 + r sin α -r sin β + φ 2 2 , 2π) (6)
φ 2 = arccos 1 8r 2 (6r 2 -(y 2 -y 1 ) 2 -(x 2 -x 1 ) 2 + 2r 2 cos(α -β) +2(x 2 -x 1 )r(sin α -sin β) + 2(y 2 -y 1 )r(cos β -cos α)) (7)
φ 3 = mod(mod(β, 2π) -α + 2φ 2 , 2π) (8)

C. The Trajectory Prediction Problem

Let us divide the aircraft trajectory into N flight phases. For each phase i, defined over the time period [t

(i) 0 , t (i)
f ], a state vector x (i) (t) and a control vector u (i) (t) are defined. In this paper, the state vector x = [v, h, s, m] is composed, respectively, by the True Airspeed (TAS), the geometric altitude, the distance to go and the mass of the aircraft; while the control vector u = [T, γ, β] is composed, respectively, by the thrust, the aerodynamic flight path angle and the speed-brakes deployment.

The dynamics of x are expressed by the following set of ordinary differential equations (ODE), which consider a point-mass representation of the aircraft reduced to the well known "gamma-command" model (i.e. vertical equilibrium is assumed):

dv dt = v = T (v, h, π) -D(v, h, m, β) m -g sin γ (9a) dh dt = ḣ = v sin γ (9b) ds dt = ṡ = v 2 cos 2 γ -W x (s, h) 2 + W s (s, h) (9c) dm dt = ṁ = -f (v, h, π) (9d)
where π is the throttle; D is the aerodynamic drag; W x and W s are, respectively, the cross and along path wind components; g is the local gravity acceleration; and f is the total fuel flow.

From an initial aircraft state (x 0 ), a trajectory prediction algorithm [START_REF] Dalmau | A Fast and Flexible Aircraft Trajectory Predictor and Optimiser for ATM Research Applications[END_REF] aims at computing future states, based on the aircraft intent inputs, weather forecast data and an aircraft performance model.

Mathematically, the aircraft intent of each flight phase i could be given as a certain control vector closing the two degrees of freedom of Eq. ( 9), plus an end condition [START_REF] González-Arribas | Automated optimal flight planning based on the aircraft intent description language[END_REF]. In practise, however, aircraft are not operated following specific T and γ profiles, and these controls are not known beforehand. Instead, climbs and descents are typically composed by constant Mach (M ), calibrated airspeed (CAS) or energy share Factor2 (k) segments performed at maximum climb or idle thrust, respectively; while in the cruise phase aircraft typically fly at constant pressure altitude (h p ) and M . Thus, in a generic formulation, two path constraints close the mathematical problem for each phase defining the aircraft trajectory:

h i (x(t), u(t), p) = 0 ( 10
)
For instance, the first path constraint of an hypothetical phase could enforce to fly at constant pressure altitude (i.e. ḣp = 0), while the second one could restrict the CAS to be constant (i.e. vCAS = 0). In such case, the two parameters defining the pressure altitude and CAS values of the concerned phase must be specified to perform the numerical integration.

An excellent description of the different path constraints that could be defined in an ATM context is given in [START_REF] Gallo | Trajectory computation infrastructure based on BADA aircraft performance model[END_REF].

The dynamics of x (9) together with the two path constraints (10) form a system of differential algebraic equations (DAE). Provided that the path constraints are meaningful (i.e. they close the mathematical problem), u can be explicitly determined, reducing the DAE system to an ODE system suitable for numerical integration using standard numerical procedures. The integration is performed until reaching the end condition, which triggers the switch to the following phase.

III. METHODOLOGY

In this Section, we describe the methodology followed in this work. In Subsection III-A, we give a general description of the methodology. Then, in Subsections III-B, III-C and III-D, we give more details about the different modules in charge of generating the emergency trajectory.

A. High-Level Methodology

The methodology followed to generate a 4D emergency trajectory is depicted in the diagram of Figure 5. It consists of three main stages, developed by the following modules: and a set of minimum radius of turn (R min,1 ...R min,n ) as a function of the altitude and the track (h, χ g ). In order to compute these parameters, both aircraft performance and weather-i.e., temperature, pressure and winds (τ , p and w respectively)-are considered. At this point, a maximum value for the bank angle (which depends on the emergency) is also taken into account (φ max ). More details regarding the path bounds generator module can be found in Section III-B.

• Path generator: the path is generated by taking into account the bounds computed in the previous stage. The method used is based on a combination of the RRT*used to generate a first version of the path from a set of nodes in a grid-and Dubin curves, used to replace the straight-lines path obtained with RRT* with curves connecting each pair of points. This path links the aircraft position with the landing site, ensuring that all obstacles (and/or convective weather) in the way are avoided. More details regarding the path generator module can be found in Section III-C. • Motion planner: by using the path and the distance to go to the landing site, a 4D trajectory is generated with a TP, which refers to the process of computing a trajectory given a known sequence of control inputs.

From an initial state-with associated initial altitude and speed (h 0 and v 0 respectively)-a TP algorithm aims at computing future states, based on the flight intent, weather data and an aircraft performance model (more details in Section III-D).

B. Path Bounds Generator

The path bounds generator, for each of the profiles used by the motion planner (Section III-D), derives a finite number of flight path angles and minimum turn radius as a function of the altitude and the track. Then, these parameters are sent to the path generator. These parameters are computed as follows:

• Flight path angles (γ 1 to γ n ): for any given trajectory, the flight path angle is continuously changing as a function of the altitude (h) and the track (χ)-except if the aircraft follows a fixed-flight-path-angle profile, which usually happens only during the last part of the flight, close to the airport. In order to speed up the computational time of the path generator, the corresponding emergency trajectory is divided in several sections, and a constant flight path angle is derived for each of these sections, as depicted in Figure 6. Thus, the path generator needs to comply with a lower number of flight path angles, which do not continuously change along the trajectory. Still, the path bounds generator computes the flight path angles as a function of h and χ.

C. Path Generator

In this section, we describe the path generator, which is the module in charge of generating the path linking the aircraft position and the landing site. It is important to highlight the fact that this path is a three-dimensional path (3D). However, the path generator does not take into account the aircraft performance, so it is the task of the next module, the motion planner (Section III-D), to generate a flyable trajectory by the aircraft. The path computed by the path generator, however, ensures that all the obstacles in the way are avoided. In order to accomplish it, the path generator needs the terrain data, which in this work is represented by a cube. In this cube, each point p generated by the RRT* algorithm is a state vector containing the latitude ψ, the longitude λ, the altitude h and the track χ. In addition, both the track and the altitude of the aircraft are defined within certain bounds:

h terrain (λ, ψ) < h < h max (11) χ ∈ [0, 360[ (12)
where h terrain is the terrain altitude and h max is the aircraft ceiling altitude. In this study, the path is represented by a set of n Dubins curves (d 1 , d 2 , ..., d n ) (Figure 7) such as ∀i, 1 ≤ i < n, the final point of the d i is the first point of d i+1 . The selected objective function is the sum of the distance of each Dubins curve. For a path (d 1 , d 2 , ..., d n ), the cost is defined as follows :

cost(d 1 , d 2 , ..., d n ) = n i=1 distance(d i ) (13)
where distance(d i ) is computed as follows :

distance(d i ) = r 1 * φ 1 + L S + r 2 * φ 2 if CSC Dubins curve r 1 * φ 1 + r 2 * φ 2 + r 3 * φ 3 otherwise (14
) One of the main functions of sampling-based path planning algorithms is the free space checking function. The performance of this function is essential to have a very efficient algorithm in terms of computing time. The proposed function only verifies the straight line between the start point and the final point of the Dubins curve. This implies enlarging the obstacles in order to be sure that Dubins curves will be in the free space. In some cases, the straight line is in the free space but the Dubins curve passes through an obstacle (Figure 8). In order to avoid this problem, the obstacles have to be enlarged horizontally by a distance corresponding to the turn radius of the aircraft trajectory. The turn radius, as explained in the previous Subsection (Subsection III-B), depends on the altitude and the track, and the chosen radius is equal to the maximum of all radii. Therefore, the cells located at a horizontal distance lower than this given radius of an obstacle cell become obstacles (Figure 9(a)). Moreover, all cells under these cells are considered also as obstacles (Figure 9(b)).

The path generator, apart from the terrain data, needs to take into account the following constrains:

• Descent constraints: the path generator has to take into account the flight path angles computed by the path bounds generator. This means that, in addition to the obstacle constraints, two points can be connected only if the flight path angle between these two points-for the current altitude and track-corresponds to the flight path angle computed with the path bounds generator. The associated constraints can be written as follows (where γ prof ile refers to the flight path angle computed by the path bounds generator):

∀i < n arctan d horizontal (d i ) |z end -z start | = γ prof ile (15)
• Track constraints: in this study, the track constraints taken into account are the following :

-Initial track: at the start position, the aircraft has a given track. The emergency trajectory has to start with this orientation. -Turn radius: throughout the flight, the pilot changes the track of the aircraft to avoid the obstacles, by following a given turn radius, whose minimum value is computed by the path bounds generator. This radius is used when modeling the turns with Dubins curves. -Landing site track: the selected landing site has also a given orientation, which means that the final track of the trajectory has to match this orientation. The proposed method to generate the path is based on the RRT* algorithm. Each node N generated during the sampling phase is defined by its coordinates in free space. The nodes contain three additional pieces of information to construct the structure of the tree. The first piece of information is the track, which is randomly generated. The second piece of information is the cost of the node N . Finally, the third piece of information is the parent node N parent . It represents the previous node through which the shortest path to reach N passes. At each step of the algorithm, a new random point is generated, located in the free space. Then, this point is steered to its nearest point to be within a distance lower than the neighborhood radius r (Figure 10). Next, the Dubins curve horizontal distance is computed and the altitude of the point x horizontalSteering is modified to have a descent angle which corresponds to the flight path angle generated by the path bounds generator (Figure 11). Finally, the point is connected to the tree if there is no obstacle in the way.

D. Motion Planner

The motion planner is the module in charge of generating a 4D trajectory by using the 3D path obtained from the path generator, together with the current aircraft performance and weather (i.e., temperature, pressure and wind). Different predefined vertical profiles are considered in the motion planner and each one is considered independently. Given a triggering event (i.e. emergency situation for a particular scenario), a profile from this predetermined catalogue is automatically chosen by the tool and the corresponding 4D trajectory is computed. Furthermore, the path bounds generator (Section III-B) also considers these predefined profiles in order to compute the set of flight path angles that define the profile and the associated minimum turn radius. It is worth noting that these pre-computed profiles are not given in terms of 4D trajectories, but as a sequence of aircraft intents, tailored to each case.

The set of vertical profiles (Table I) considered in this work have been designed to be as generic as possible in order to cover the widest possible range of emergency situations. However, most of emergencies are aircraft model dependent, so the procedure described in each aircraft's flight crew operating manual (FCOM) would ultimately affect the definition of the vertical profiles.

In this work, the TP logic implemented in both the motion planner and the path bounds generator considers the following aspects:

• Trajectory phases: the vertical profile of the aircraft trajectory is split in a finite number of phases. Each phase is specified by an end condition and two aircraft intents. • Aircraft intents: operational intents, similar to those one could find in state-of-the-art FMSs, such as flying at constant Mach, CAS, altitude, throttle setting, flight path angle, etc. • End condition: needed to define the transition from one phase to the other where intents might change and/or aerodynamic conditions might change (i.e. usage of flaps, landing gear, etc). Example: a phase is flown at constant Mach and Idle thrust (2 intents) until the moment the CAS achieves a given value (end condition of the phase given in terms of CAS); then, the next phase is flown at constant CAS and constant vertical speed (2 different intents). • Trajectory computation (trajectory prediction): a full 4D trajectory is unequivocally specified given a sequence of phases, with the 2 intents per phase and the end condition specified. This allows to close the 2 degrees of freedom of the mathematical model describing the movement of the aircraft in the vertical plane (i.e. numerically integrate the equations of movement). The list of profiles, including how the trajectory is computed in each case are fixed by the motion planner and the pilot will not be able to change them. However, in some cases, more than one outcome could be presented by default to the pilot, who could select the desired trajectory. Moreover, for some phases the pilot could request the system to change the values for (some of) the intents. For instance, in a Mach descent, the pilot could decide the value for Mach, but could not replace that phase for a CAS descent, for instance.

The output of the motion planner component is a 4D trajectory plan. Then, this plan will have to be executed in flight by the (auto) pilot. In nominal operations, the FMS of the aircraft also computes this kind of trajectory plan. A typical example would be the aircraft descent, in where the top of descent (TOD) is fixed. This plan becomes the guidance command for the autopilot, which will try to follow the plan (in managed mode) and adjust the aircraft controls to compensate for uncertainty (mainly due to inaccuracies in the weather forecast and to a lesser extent, inaccuracies in the aircraft performance).

Although the aircraft crew is responsible for deploying speed-brakes, hyper-lift devices (i.e. flaps/slats) and landing gear at the moment they consider the most appropriate, a descent plan computed by a FMS makes some hypotheses on the moments these devices will be deployed. In this work, the same philosophy is followed, taking into account that the descent trajectory is not computed until the runway threshold, but down to a predefined point in space that depends on the type of emergency, health of the aircraft and characteristics of the landing site. Thus, some hypotheses are made regarding the usage of speed-brakes, flap/slats and landing gear, which are specified in each profile definition. The vertical profiles considered in this work depend on the following factors:

• Engines available: all engines out or at least one engine is operative. • Type of emergency: Land ASAP (as soon as possible, e.g., fire on cabin) or ANSA (at the nearest suitable airport, e.g. depressurization in cabin). • Approach procedure: airport with or without a suitable IFR (instrumental flight rules) approach procedure. • Maneuverability: aircraft fully or not fully maneuverable. In addition to these factors, in case no engines are available, the availability of fuel on board is also considered. The combination of these factors, taking into account that some combinations for the "all engines out" cases lead to the same profile, gives as a result the list of profiles of Table I. 4Yes ANSA --Yes

IV. RESULTS

In this Section, we show the results obtained in this work. In Subsection IV-A, we describe the scenario and we explain the inputs needed to apply our methodology. Then, in Subsection IV-B, we describe the case-studies. Finally, in Subsection IV-C, we show the computed trajectories.

A. Scenario Description and Inputs

The scenario chosen for this work is set in a challenging mountainous area around Grenoble airport, in France, where 15 available landing sites were identified. In order to successfully generate 4D emergency trajectories, the different modules need a series of inputs:

• Aircraft performance model (APM): in this work, we use the base of aircraft data (BADA) APM by Eurocontrol [START_REF] Eurocontrol | User Manual for the Base of Aircraft Data (BADA) Family[END_REF], which is a widely and recognised APM in the air traffic management community, providing performance models for a large number of aircraft with very accurate results. Fort all case-studies an Airbus A320-232 has been used. • Weather data: in this work, the longitudinal component of the wind is modelled by a smoothing spline [START_REF] De | On calculating with B-splines[END_REF], w : R → R, such that:

w(h) = nc i=1 c i B i (h), (16)
where B i , i = 1, . . . , n c are the B-spline basis functions and c = [c 1 , . . . , c nc ] are the control points of the smoothing spline.

It should be noted that the longitudinal wind has been modelled as a function of the altitude only, as done in similar works [START_REF] Sáez | Traffic Synchronization in Terminal Airspace to Enable Continuous Descent Operations in Trombone Sequencing and Merging Procedures: An Implementation Study for Frankfurt Airport[END_REF]. The control points of the spline approximating the longitudinal wind profile are obtained by fitting historical weather data. This data is obtained from gridded binary (GRIB) formatted files provided by the global forecast system (GFS) of the National Oceanic and Atmospheric Administration (NOAA) [START_REF]Global Forecast System (GFS)[END_REF]. Apart from the wind, this data is used as well to obtain the values for temperature and pressure as a function of altitude.

• Terrain data: this data is given as a 3D matrix, where altitude data is given every 30 meters. In addition, the location of the landing sites is included in this matrix, each one represented as a 3D point. Furthermore, if the landing site corresponds to an airport, the orientation and the length of the runway is detailed.

B. Case-Studies

Two case-studies have been analyzed in this work. For each case, the corresponding profile from the list of profiles (Table I) has been chosen. In Tables II andIII, the different aircraft intents that appear are the following: ALT, maintain a constant altitude; MACH, maintain a constant Mach number; CAS, maintain a constant CAS; ACC: accelerate at a given acceleration or with a given load factor; DEC, decelerate at a given deceleration or with a given load factor; THR: keep a given throttle setting. For the ACC and DEC cases, we have used the energy share factor value proposed by BADA [START_REF] Eurocontrol | User Manual for the Base of Aircraft Data (BADA) Family[END_REF], The two case-studies considered in this work are the following:

• Dual-engine failure: this case-study assumes a dualengine failure, which leads to a total loss of thrust (as the aircraft considered is a twin-engine aircraft, an A320). The set of phases and intents considered in this case are described in Table II. The aircraft flies at the recommended speed for maximum range until being sure to "make" the best landing site available. Once this is ensured, speed may be increased (to a certain CAS value, CAS n ) to increase rate of descent if necessary.

Then, a final deceleration should deliver the aircraft in landing configuration at V app -approach speed, which is the final approach speed when the flaps/slats are in landing configuration and the landing gear is extendedon a 650ft/Nm glide path, and at least 5NM before touchdown area. • ANSA case: this emergency could correspond, for instance, to a fuel leak emergency when the aircraft is in cruise. The set of phases and intents considered in this case are described in Table III. We propose that the aircraft keeps its current (cruise) speed and altitude. Then, a typical Mach/CAS descent (cruise Mach and a certain CAS, CAS n ) is performed until reaching the IAF (initial approach fix), with the aircraft aligned with the initial approach segment and at V F 1 , which is the target speed for flaps deployed in configuration 1. For both profiles, the final point is set at a distance of 5NM from the runway, at an altitude of 3250ft with respect to the runway altitude 6 and assuming a final aircraft mass of 60,000 kg. In addition, we enforce the aircraft to be aligned with the runway at this point. Finally, the initial altitude is set to 37,000 ft. By taking into account all these constraints, it was impossible to find a feasible trajectory linking the aircraft position with all landing sites. There were cases in which there was an obstacle between the final point and the landing site, making it impossible for the aircraft to reach the landing site by following this type of profiles. In the end, only 3 landing sites in the Grenoble area were reachable. The trajectories depicted in Subsection IV-C were generated for one of these landing sites, located at 9514 ft above mean sea level.

C. Final Trajectories

Figures 12 and 13 show the resulting trajectories for the profiles specified in Section IV-B. In both cases, the first part of both profiles-the cruise phase and the green dot descent in profiles 2 and 7 respectively-is not flown due to the proximity of the landing site. The resulting trajectory follows in both cases an MMO/VMO (maximum operating Mach/maximum operating CAS) profile until the final deceleration to V app and V F 1 for profiles 2 and 7 respectively.

The computational time for the generation of the trajectory is composed of the computational time spent by the different modules considered in this work. Regarding the path bounds generator, the flight path angles and the minimum turn radius are generated in 7 seconds. A higher computational time is usually spent on the path generator. Depending on the number of iterations, the algorithm computes one, two or three trajectories. Then, in order to ensure that a solution is always found for each airport, the number of iterations needs to be very high (at least 30000). This is due to the fact that the constraints are very restrictive. However, if one solution is sufficient (i.e. one trajectory for one airport), about 10000 iterations are enough. The computational time for the RRT* algorithm increases very fast when the number of iterations increases, leading to a computational time of 3s, 20s, 100s and 500s if the number of iterations is 10000, 30000, 50000, 100000 respectively.

Finally, the motion planner takes between 2 to 5 seconds to generate the 4D trajectory, depending on the profile considered. In general, the emergency trajectories are generated in less than 15 seconds. Still, we are still using a very simple prototype to test our concept; in the future, we expect to optimize our code, as well as to run the tests in a more powerful computer or server.

V. CONCLUSIONS

In this work, we showed a novel methodology to generate 4D trajectory plans for emergency situations, designed to safely lead the aircraft to a designated landing site. For that purpose, we used a combination of different methods: the RRT algorithm, Dubins paths and a TP algorithm.

We showed that we are capable of generating emergency trajectories for a wide range of emergency situations in a short amount of time (less than 15 seconds). Ultimately, the concept presented in this paper aims to be part of an in-flight support tool for the FMS. The 4D trajectory plans, after being injected in the FMS, would be executed in-flight by the (auto) pilot.

This work focuses only on the generation of the trajectory; in the future, an insight on the landing site selection process will be given, and a wider range of situations will be investigated in order to assess the feasibility of our concept.

Fig. 1 .

 1 Fig. 1. New connection building (a) and Connection replacement (b)

Fig. 2 .

 2 Fig. 2. Comparison between the graph generated by the RRT algorithm (left) and the one generated by the RRT* algorithm (right)

Figure 4 Fig. 4 .

 44 Fig. 4. Left-Right-Left Dubins curve

Fig. 5 .

 5 Fig. 5. Generation process of the emergency trajectory

Fig. 6 .

 6 Fig. 6. Example of descent profile with a finite number of flight path angles • Minimum turn radius (R min ): it is achieved when the aircraft flies with the maximum bank angle and it is also affected by the aircraft current ground speed. The path bounds generator computes the minimum turn radius as a function of altitude, track, and ground flight path angle (γ g ). A bank angle of 30 degrees is considered if the

Fig. 7 .

 7 Fig. 7. Path example

Fig. 8 .

 8 Fig.8. Example of collision with an obstacle after the addition of Dubins curves : The straight line path is drawn in green and is in the free space. The Dubins path is in blue and passes through the obstacle in red.

Fig. 9 .

 9 Fig. 9. Enlargement of obstacles

  Fig. 10. Horizontal Steering

  Fig. 12. Profile 2

TABLE II PROFILE 2 :

 II2 ALL ENGINES OUT + ASAP + NO FUEL ON BOARD

	Phase		a/c intent #1	a/c intent #2
	Deceleration to	ALT	Current altitude	THR	IDLE
	green dot 5
	CAS descent	CAS	Green dot	THR	IDLE
	Acceleration to increase	ACC	-	THR	IDLE
	rate of descent
	Descent at higher speed to	CAS	CAS n	THR	IDLE
	increase rate of descent
	Deceleration to V app	DEC	-	THR	IDLE
	which varies for acceleration and deceleration and depends on
	the flight phase (descent or climb).

The work presented in this paper has received funding from the Clean Sky

Joint Undertaking (JU) under grant agreement No 864771, corresponding to the SafeNcy project (https://cordis.europa.eu/project/id/864771). The JU receives support from the European Union's Horizon 2020 research and innovation programme and the Clean Sky 2 JU members other than the Union. The opinions expressed herein reflect the authors view only. Under no circumstances shall the Clean Sky 2 JU be responsible for any use that may be made of the information contained herein.

For an aircraft that is changing both altitude and speed, the energy share factor specifies how much the available thrust is allocated to the vertical and speed evolution.

Regardless of the availability of an IFR approach procedure and regardless of the maneuverability status of the aircraft. Note that big aircraft relying on hydraulic power are not fully maneuverable without engine power.

Depressurization or smoke in cabin. For the last part of this emergency trajectory, the aircraft intents from Profiles 7 to 10 are used.

For the Airbus A320, the green dot speed is the minimum operating speed in managed mode and clean configuration, being approximately the best liftto-drag ratio speed.

For profile

[START_REF]Autonomí™Autonomous Safety-enhancing Technologies[END_REF], as the engines are still working, it might not be needed to finish at such a steep glide path. In this case, the final point could be set at 5NM from the runway and at 1500ft above the runway altitude. Still, in order to simplify the problem, we considered the same final point for both cases.

ACKNOWLEDGMENTS

We would like to thank Capt. Claude Godel for his valuable advice when defining the vertical profiles used by the motion planner.


Zihao Wang
email: zihao.wang@inria.fr

Clair Vandersteen

Thomas Demarcy

Dan Gnansia

Charles Raffaelli

Nicolas Guevara

Hervé Delingette

Inner-ear Augmented Metal Artifact Reduction with Simulation-based 3D Generative Adversarial Networks ⋆

Keywords: Artifact Reduction, Deep Learning, GAN

Metal Artifacts creates often difficulties for a high quality visual assessment of post-operative imaging in computed tomography (CT). A vast body of methods have been proposed to tackle this issue, but these methods were designed for regular CT scans and their performance is usually insufficient when imaging tiny implants. In the context of post-operative high-resolution CT imaging, we propose a 3D metal artifact reduction algorithm based on a generative adversarial neural network. It is based on the simulation of physically realistic CT metal artifacts created by cochlea implant electrodes on preoperative images. The generated images serve to train a 3D generative adversarial networks for artifacts reduction. The proposed approach was assessed qualitatively and quantitatively on clinical conventional and cone beam CT of cochlear implant postoperative images. These experiments show that the proposed method outperforms other general metal artifact reduction approaches.

Introduction

Computed Tomography (CT) is one of the most widely used imaging techniques in clinical practice. The physical principles of CT lead to the unavoidable creation of artifacts in the reconstructed images in the presence of dense materials, i.e., those composed of atoms with high atomic numbers. Several physical phenomena contribute to the creation of such artifacts, including X-ray beam hardening, X-ray scatter, electronic noise, edge effects and also the geometrical characteristics of metal parts. The artifacts are commonly found in routine clinical postoperative imaging, for instance due to fixation plates in orthopaedics, cochlear electrode implants in otology, contrast agents, etc. These spurious signals in CT images may impair postoperative analysis. For instance, during cochlear implant (CI) surgery, an electrode array inserted along the cochlear scala tympani is usually comprised of a metal alloy, for its high electrical conductivity. The existence of metal artifacts in postoperative CT makes the evaluation of the position of the electrodes along the scala difficult. The knowledge of the relative position of the cochlear implant is one of the main determinants for assessing the success of the surgery and leads to appropriate and more personalized patient care.

Metal artifact reduction (MAR) methods aim to decrease the extent of such artifacts (see: 1). Classical non-learning-based MAR algorithms are divided into two groups: corrupted projection recovery and iterative image reconstructionbased methods [START_REF] Mehranian | X-ray ct metal artifact reduction using wavelet domain 0 sparse regularization[END_REF]. In the former case, projections corrupted by the presence of metal absorbing the X-rays are detected and then replaced by predicted or interpolated values, based on prior knowledge (Kalender et al., 1987a). The efficiency of the approach is related to the ability to recover the projected signals in the absence of metal parts [START_REF] Mehranian | X-ray ct metal artifact reduction using wavelet domain 0 sparse regularization[END_REF]. In the case of iterative methods, the missing data in image or projection space is estimated on the basis of statistical principles, possibly including prior knowledge. Aside from Filtered Back Projection (FBP) based methods, [START_REF] Naranjo | Metal artifact reduction in dental ct images using polar mathematical morphology[END_REF] introduced mathematical morphology algorithms for MAR by converting the image to polar coordinates centered on the metal artifact. [START_REF] Zhang | Convolutional neural network based metal artifact reduction in x-ray computed tomography[END_REF] Sinogram NO 2D Simulation (BH) YES RL-ARCNN [START_REF] Huang | Metal artifact reduction on cervical ct images by deep residual learning[END_REF] Image NO 2D Simulation (BH) YES DestreakNet [START_REF] Gjesteby | A dual-stream deep convolutional network for reducing metal streak artifacts in CT images[END_REF] Image NO 2D Simulation (BH) YES DudoNet++ [START_REF] Lyu | Dudonet++: Encoding mask projection to reduce CT metal artifacts[END_REF] Sinogram+Image NO 2D Simulation (BH;SC;EN) NO CycleGAN [START_REF] Nakao | Regularized three-dimensional generative adversarial nets for unsupervised metal artifact reduction in head and neck ct images[END_REF] Image NO 3D Unsupervised NO cGAN (Wang et al., 2019a[START_REF] Wang | Metal artifact reduction and intra cochlear anatomy segmentation inct images of the ear with a multi-resolution multi-task 3d network[END_REF] Image YES 2D Paired Data YES MARGAN (Proposed) Image YES 3D Simulation (BH;SC;EN) YES methods have been proposed that are summarized in Table 1. [START_REF] Zhang | Convolutional neural network based metal artifact reduction in x-ray computed tomography[END_REF] introduced CNNs as prior information in the sinogram (projection) space for the inpainting or sinogram completion task using a simulated dataset in the training stage. However, this method needs either the CT sinograms raw data (usually unavailable to the typical user) or to project back the input image in order to fill-in the missing traces. This limits its application to datasets that include CT raw data, and the sinogram-based MAR algorithms tend to generate over-smoothed images due to their filtering effect. [START_REF] Huang | Metal artifact reduction on cervical ct images by deep residual learning[END_REF] developed a deep learning network, RL-ARCNN, in image space to predict residual images (the difference between the images with and without artifacts) to remove metal artifacts in cervical CT images. The 2D network DestreakNet was proposed in [START_REF] Gjesteby | A dual-stream deep convolutional network for reducing metal streak artifacts in CT images[END_REF] for streak artifact reduction as a post-processing step in order to recover the details lost after the application of the interpolation-based normalized MAR Meyer (2010) algorithm. [START_REF] Lyu | Dudonet++: Encoding mask projection to reduce CT metal artifacts[END_REF] proposed Dudonet++ for 2D CT metal artifact reduction. Their approach relies on processing the image with artifacts (henceforth referred to as artifact image) in both sinogram and image spaces in order to restore fine details in the image. Their quantitative evaluation shows that the Dudonet++ is effective for artifact reduction on simulated CT images but it lacks a quantitative evaluation on a clinical dataset. Furthermore, the method of [START_REF] Verburg | Ct metal artifact reduction method correcting for beam hardening and missing projections[END_REF] uses a sinogram replacement method for beam hardening correction, which is constrained to a linear model for interpolation thus not always optimum for images containing non-linear effects, for instance in the case of CI MAR (see Tab. 4).

MARGAN approach

Recently, generative adversarial networks (GAN) were devised for solving MAR problems instead of CNN classification or regression networks, owing to their ability to generate high quality images. Wang et al. (2019a[START_REF] Wang | Metal artifact reduction and intra cochlear anatomy segmentation inct images of the ear with a multi-resolution multi-task 3d network[END_REF] proposed a conditional GAN (cGAN)approach for CT images with cochlear implants (CI), using a collection of paired and registered post-and preoperative cochlear implant volumes to train 2D/3D cGANs for inner ear MAR. A difficulty in this approach is to collect and, most importantly, to register (preoperative) artifact-free and (postoperative) artifact im-ages. This registration problem must be able to cope with the presence of outliers due to the presence of artifacts.

Nakao and et al. also proposed a MAR method based on CycleGANs for artifact reduction in dental filling and neck CT images. The approach is unsupervised and aims to achieve a cross-domain (artifact and artifact-free dataset) style transformation through feature swapping. This approach does not require training on paired datasets, i.e., with and without artifacts, but CycleGAN performance significantly worsens when unpaired data is used [START_REF] Zhu | Unpaired imageto-image translation using cycle-consistent adversarial networks[END_REF] for training instead of paired data. This approach was qualitatively compared with the manual corrections available in commercial CT scans and quantitatively assessed on synthetic datasets. While the output of the CycleGANs seems effective, this method may not be useful for the reduction of tiny artifacts like cochlear implants, due to the difficult separability of artifacts in feature space.

Simulating the presence of metal parts

In this paper, we propose a GAN-based MAR method that relies on simulated training data and is suitable for preand postoperative images. To the best of our knowledge, our approach is the first MAR algorithm that combines the physical simulation of metal artifacts with 3D GAN networks. While classical GAN-based methods such as (Wang et al., 2019a) rely on the existence of paired images with and without artifacts for training, our approach has several advantages. First, only preoperative images (without artifacts) are required for the training stage, because the generation of the corresponding artifact image is based on physical simulation. This allows a large set of training images (800 images) to be used, without the need for registering the pre-and postoperative images. Second, the nature of artifacts can be easily modulated by controlling the complexity of the artifact simulation model complexity. Third, we introduce the concept of augmented metal artifact reduction by optionally adding landmarks in the corrected image that indicate the central location of metal parts. More precisely, in this paper, we show that for the postoperative cochlear implant CT images, the location of each electrode center can be identified in the corrected image such that ENT (ear, nose and throat) surgeons can assess the quality of the implantation surgery. Compared to CycleGANs [START_REF] Nakao | Regularized three-dimensional generative adversarial nets for unsupervised metal artifact reduction in head and neck ct images[END_REF], the MARGAN approach allows artifacts to be easily disentangled from the background. This is why we believe this approach is probably more appropriate to attenuate artifacts created by tiny implants. Fourth, MARGAN was evaluated on postoperative, cone beam CT images. Finally, MARGAN was developed as a 3D GAN since metal artifacts usually vary continuously between slices. The contributions of the MARGAN framework is summarized in Tab. 1.

The MARGAN method is based on two main stages (see Fig. 3). In the first stage (Fig. 5), given a preoperative image from the training set, one or several images with metal artifacts are generated. This requires a rough segmentation of the structures of interest, the position of metal parts (e.g., electrode arrays) and the simulation of artifacts based on a CT image formation model. Furthermore, the location of the electrode arrays is added to the generated images. In the second stage (Fig. 3), a 3D GAN is trained using preoperative and corresponding simulated artifact images. The GAN loss is improved by adding a term based on Retinex theory to decrease the image blur in generated images. After training, the GAN is applied on a postoperative image without any segmentation or other preprocessing. It results in images with attenuated metal artifacts but also with landmarks corresponding to electrode centers.

The MARGAN method was applied to a set of inner ear CT images to reduce the artifacts created by cochlear implants. Qualitative and quantitative results are provided for 33 paired pre-and postoperative CT images, including a comparison with two classical open source MAR algorithms. Qualitative evaluation of cone beam CT (CBCT) postoperative images is also provided.

This paper extends the initial work published in [START_REF] Wang | Deep learning based metal artifacts reduction in post-operative cochlear implant ct imaging[END_REF] in several ways. The artifact simulation model is more sophisticated, including scattering effects and electronic noise of the CT system detectors. The algorithm evaluation is more comprehensive, with the addition of paired CT images, CBCT images and a study of the impact of the Retinex loss. The postoperative electrode position is assessed in a few cases with postmortem photographic views of the cochlea.

The paper is organized as follows: In section 2, we introduce the CI and CI metal artifact simulation procedures (the gray box in Fig. 3). In section 3, the network implementation is described (the green box in Fig. 3). Results of the MAR-GAN algorithm are presented in section 4. Sections 5 and 6 discuss the contributions and limitations of the proposed approach.

Simulation of metal artifacts in CT images

The simulation of metal artifacts in CT images from artifactfree images entails i) simulating the presence of metal parts in the images as shown in Fig: 2 and ii) simulating the creation of artifacts caused by those metal parts as shown in Fig: 4. The former algorithm is completely dependent on the organ or implant considered, while the latter is far more generic, based on the physics of image formation.

The processing pipeline to generate the training set for the MARGANs is displayed in Fig. 2. In this section, we consider the case of preoperative CT images of the inner ear prior to cochlear implant surgery. The objective is therefore to simulate, in these preoperative images, the addition of metal electrode arrays associated with the implant.

The 3D CT volumes of the inner ear, written as ( ), are first rigidly registered on a template image by a block matching algorithm [START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF]. The template is a sample CT image that has been manually cropped around the temporal bone. The registration is necessary to cope with the variations of field of view and pose in the input image dataset. A region of interest (ROI) is then cropped to get a cochlear volume suitable for further processing. We then fit a parametric shape model [START_REF] Demarcy | Automated analysis of human cochlea shape variability from segmented CT images[END_REF] to automatically reconstruct the shape of the cochlea (step (II) of Fig: 2). The accuracy required for the registration and segmentation steps is limited.

The signed distance map [START_REF] Wang | Metal artifact reduction and intra cochlear anatomy segmentation inct images of the ear with a multi-resolution multi-task 3d network[END_REF] from the fitted triangular mesh of the parametric shape model is generated as shown in step (III). It is then thresholded (step (IV)) to create a 3D tubular binary mask near the center-line of the cochlea. This mask corresponds to the probable location of the electrodes after a CI intervention. Finally, in step (V), the voxel values in Hounsfield units (HU) of the mask region are then set to 3071 which is the maximum detectable HU of the CI metal artifacts. This creates the image train ( ) used for training the GAN network.

Simulation of beam hardening, scattering and electronic noise due to metal parts

The metal parts have large absorption rates of X-ray energy which is the cause of the visible artifacts in CT images. It impacts the whole image formation process through several physical effects. Our previous work [START_REF] Wang | Deep learning based metal artifacts reduction in post-operative cochlear implant ct imaging[END_REF] only considered the simulation of the beam hardening effect, inspired by the work of [START_REF] Zhang | Convolutional neural network based metal artifact reduction in x-ray computed tomography[END_REF]. In this paper, we improve the realism of the simulated artifacts by also including the X-ray scatter effect through Monte Carlo simulation and the detector electronic noise. The three main physical effects governing the generation of metal artifacts are described below, along with the processing pipeline.

Beam hardening effect For a monoenergetic X-ray entering a homogeneous material of thickness along direction at position , , the number of photons ( , , ) is given by the Beer-Lambert law :

( , , ) = 0 -( , )
where 0 is the initial photon number and ( , ) is the linear attenuation coefficient of the material.

The attenuation coefficient depends on the energy of the input photon ( ), and therefore for a polychromatic X-ray beam having the energy distribution (or spectrum), ( ), the number of photons received by the entire detector surface is then:

= ∫ 0 ( ( ) -∭ ( , , , ) + ( )) (1)
where 0 and are the minimum and maximum energies for a fixed tube peak voltage, and ( ) is an additive offset that captures X-ray scattering.

Scattering effect

The Compton effect applies to incoming X-ray photons that interact with the free electrons in the traversed materials. This effect results in random changes (scatter) in the directions of the photons, which may still reach the detector plate despite collimator devices. The Compton scatter is enhanced in the presence of metal parts, thus resulting in an offset in the number of photons ( ) and leading to a reduction in the image contrast. Computing this additional scatter is very complex as it depends on the projected plane and the material and geometry of the tissue surrounding the metal parts. To this end, we use Monte Carlo simulation to estimate the offset value ( ) for different detector positions and orientations. The governing equation for the simulation provides the emission energy ( ) of a polychromatic ray deviating by an angle from its initial trajectory :

( ) = ∫ 0 (1 + ∕ 2 )(1 - ( )) (2)
where is the electron mass and the speed of light. To estimate the scatter effect inside the cochlea on X-ray detectors, we use the [START_REF] Zubal | Computerized three-dimensional segmented human anatomy[END_REF] head phantom where metal parts are roughly positioned inside the temporal bone. Based on the MCGPU software [START_REF] Arjovsky | Accelerating monte carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit[END_REF] performing GPU Monte Carlo simulations of photon transport in voxelized geometry, we simulate thousands of X-ray photon trajectories at different energies, positions and orientations through the head and produce both the scatter-free sinogram ( ) and the scatter sinogram offset ̃ ( ). The scatter sinogram offset is corrected by a scale factor such that the resulting scatter to primary ratio = mean( ̃ ( )) mean( ( )) , falls within the range of 0.1% to 2%, which was experimentally found by Glover et al. [START_REF] Gh | Compton scatter effects in ct reconstructions[END_REF]. This is simply done by randomly picking a ratio within 0.1% to 2% and computing

( ) = mean( ( )) mean( ̃ ( )) ̃ ( ) (3)
such that mean( ( )) mean( ( )) = . The same ratio is used for simulating all sinograms of the same image to obtain spatially consistent artifacts.

The computation of the scatter offset is dependent on the X-ray energy, position, and orientation but is independent of the input image as it relies on the digital head phantom augmented with metal parts next to the temporal bone. Only the scatter to primary ratio varies between different volumes. This implies that the scatter sinograms can be precomputed, thus alleviating the computational load when generating images with metal artifacts.

Detector Noise Once photons hit the x-ray detector, the scintillator transforms the deposited energy into visible light, while a photomultiplier translates this light into an electric signal. In this process, some electronic noise is introduced which can be modeled by a zero mean Gaussian distribution [START_REF] Benson | Synthetic ct noise emulation in the raw data domain[END_REF] with standard deviation :

(0, 2 ). The signal measured in each sinogram f inal can then be written as: f inal = +  (0, 2 ) where is the energy deposited as described Eq. 1 and 2 = 0.04 Žabić et al. (2013); [START_REF] Duan | Electronic noise in ct detectors: Impact on image noise and artifacts[END_REF]; [START_REF] Fuchs | Direct comparison of a xenon and a solid-state ct detector system: measurements under working conditions[END_REF]; [START_REF] Rui | Analysis of optimal ct spectrum for pet attenuation correction[END_REF].

Simulation pipeline

The overall metal artifact simulation pipeline is described in Fig. 5. In the first step, we use an X-ray energy spectrum ( ) extracted from a CT manufacturer dedicated site1 for a tungsten anode tube at 140 . The spectrum is sampled at five energies from which attenuation maps ( , , ,

) are generated. This computation is based on the Hounsfield unit formula and the water absorption coefficients as a function of energy. We then perform fan-beam projection (Step III) of the five attenuation maps to produce sinogram-like images representing absorbed energy on the CT detectors. The scattering and attenuation sinograms are precomputed on a head phantom for various orientations and positions of the source. The projection of the ROI of the head where metal parts have been inserted creates a sine trace on the scattering and attenuation sinograms. This trace is randomly sampled, then normalized as in Eq. 3 to obtain a plausible scatter to primary ratio. It is then added to the electronic noise and to the weighted sum of the five sinograms (Step IV) and a discretization of Eq. 1. Finally, inverse fan-beam projection produces the output image with metallic artifacts (Step V).

The difference between simulated images with and without scattering noise is shown in Fig. 6 with a subtraction map. We see that scattering and electronic noise can introduce significant new artifacts.

GAN-based Metal Artifact Reduction

Given pairs of preoperative and simulated postoperative images, we aim to train a network that generates the former given the latter as a way to reduce metal artifacts. The use of a GAN to tackle the MAR issue is motivated by the successful use of 2D and 3D GANs such as SRGAN [START_REF] Ledig | Photorealistic single image super-resolution using a generative adversarial network[END_REF][START_REF] Sanchez | Brain mri super-resolution using 3d generative adversarial networks[END_REF] to solve imaging Super-Resolution (SR) problems.

Network Overview

In MARGAN, two neural networks are used: the generator network produces the MAR images and the discriminator network indicates whether the input image contains metal artifacts or not. The generator network, , with weights, , aims at modeling the mapping between the image with artifacts, m , and the simulated artifact-free image, train . We denote by the 3D image created by the generator network, which should be as close as possible to train . The discriminator neural network, , tries to detect the presence of artifacts in the generated MAR images, . To train and networks, the sum of discriminator and generator losses is optimized as detailed below.

Network Architecture

The generator network architecture is similar to U-Net with convolution and deconvolution layers, skip connections and batch normalization layers to improve the training efficiency (see Fig 3). Moreover, unlike [START_REF] Sanchez | Brain mri super-resolution using 3d generative adversarial networks[END_REF] which is patch based, the input to the network consists of full 3D images as it is compatible with GPU memory. The number of filters increases gradually from 1 to 512, a number of feature maps that can fit on an 11 Gb video-memory GPU card. The discriminator network follows that of [START_REF] Sanchez | Brain mri super-resolution using 3d generative adversarial networks[END_REF] with eight groups of convolution layers and batch normalization layers combined sequentially.

Loss Functions

Discriminator Loss The discriminator network, , is trained using output images from the generator network, = ( m ), and images without any metal artifacts, nm . Following [START_REF] Sanchez | Brain mri super-resolution using 3d generative adversarial networks[END_REF], the discriminator loss enforces the ability of the discriminator network to distinguish the artifact-free images, nm , from the generated ones, :

arg max = ∼ nm log ( ) + ∼ m log 1 - ( ( ))
Generator Loss The objective of the generator network is to produce an image, = ( m ), as close as possible to the target image, train . This is why the first loss term is the mean square error (MSE), = ∼ m (| train -( )| 2 ), to encourage a similarity between generated and target voxels. But using only the MSE loss leads to blurred MAR images with a lack image detail at high frequencies.

To avoid this excessive smoothing, we propose a new loss term based on Retinex theory [START_REF] Land | Lightness and retinex theory[END_REF]). This theory is mostly used to improve images seriously affected by environmental illumination. The Retinex theory assumes that a given image can be considered as the product of environmental brightness (or illumination), ( , ), and the object reflectance, ( , ). This reflectance map contains high frequency details and is unaffected by the illumination condition, a property referred to as the color constancy phenomenon. The objective of Retinex-based algorithms is to recover the reflectance image from the original one. In single-scale Retinex approaches [START_REF] Zhang | A ultrasound liver image enhancement algorithm based on multi-scale retinex theory[END_REF], the environmental brightness is simply a Gaussian blur version of the input image and therefore log( ( , )) = log( ( , )) -log( ( , ) *  (0, )) where  (0, ) is a Gaussian function of standard deviation , and * is the convolution operator. This leads us to introduce the following Retinex loss to make its illumination part as close to 1 as possible :

= ∼ m | ( ) - log ( )-log ( ) *  (0, ) | | | (4)
where the expectation is taken over the image domain. This loss definition ensures numerically stable evaluations and enforces salient features in the image that would otherwise be attenuated. Combining it with the adversarial term = 1 2 | ( ) -1| 2 as in [START_REF] Sanchez | Brain mri super-resolution using 3d generative adversarial networks[END_REF], the full optimization target of the generator is:

arg min = ⋅ + + (5)
where is a parameter controlling the influence of the Retinex loss.

Results

Dataset

Training data

The cochlea dataset was collected from the Radiology Department of the Nice University Hospital with a GE Light-Speed CT scanner without any metal artifact reduction filters. The preoperative dataset includes 1000 temporal bone images (493 left and 597 right) from 597 patients. The original CT volumes are registered to a sample image by a pyramidal block-matching algorithm in order to spatially normalize all images, then they are resampled with 0.2 × 0.2 × 0.2 3 voxel size. They were then cropped to volumes of 60 × 50 × 50 voxels around the cochlea region. We then simulated on all volumes, the insertion of CI electrodes and the generation of metal artifacts as described in section 2. This created a set of 1000 pairs of images, with and without metal artifacts.

Evaluation Data

The evaluation dataset #1 includes 33 cadaver temporal bone CT images collected from the same site from different bodies. The imaging protocol was the same as for the training dataset but was performed before and after the implantation of CI, thus leading to 33 pre-and postoperative image pairs. The temporal bones were ground by an ENT (ear, nose and throat) surgeon, approximately along a plane perpendicular to the cochlear modiolar axis at the bottom of the scala tympani as shown in Fig. 12. Pictures of the ground bones were acquired in order to visualize the electrode array.

Finally, the second evaluation dataset includes 8 postoperative images that were acquired on a Carestream 9600 cone beam CT (CBCT) following the CI surgery. These images were resampled, registered and cropped following the same processing pipeline as the training set.

Implementation details

Artifact simulation A polychromatic X-ray source was simulated with MC-GPU v1.3, a GPU-based Monte Carlo simulator of photon transport in voxelized geometry [START_REF] Arjovsky | Accelerating monte carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit[END_REF]. To simulate the scatter effects, we simplified the contents of the human head by assuming it consists of air, water, soft tissue, bone, muscle and unalloyed titanium. Cochlear CT voxel values were converted to MC-GPU v1.3 units based on the material mapping in Table 2. The simulation of scatter was performed offline on a GPU parallel computing cluster. The beam hardening maps and the final simulation volumes were computed with Matlab 2017a on a Dell Mobile Workstation with Intel(R) Core(TM) i7-7820HQ @ 2.90GHz CPU.

Neural Networks

The networks were trained with a RM-Sprop optimizer [START_REF] Arjovsky | Accelerating monte carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit[END_REF] with learning rate = 1e-4 for the generator and = 1e-3 for the dis- criminator. The MARGAN was implemented with Tensorflow and the weight of Retinex loss was set to = 5e-5.

The batch size and number of epochs were set as 1 and 20 respectively.

Computational Time To assess the computational time, we ran the trained MARGAN for a single 3D volume on a Dell 2x Xeon Silver 4110 with a Tesla T4 GPU and on a Dell Precision 7520 Mobile Workstation with an Intel(R) Core(TM) i7-7820HQ @ 2.90GHz CPU. It took about 18.3 secs on Tesla T4 GPU node and 25.1 secs on i7-7820HQ CPU, respectively. Yet, the computational time cost can be reduced for processing a batch of volumes. This is because the I/O process occupies the major GPU/CPU computational time for processing a single volume.

Clinical Evaluation

Qualitative Study Fig. 9 shows the output of the MAR-GAN network for four patients on two selected slices together with pre-and postoperative CT images. The streak artifact patterns were largely suppressed by the MARGAN algorithm. As shown inside the yellow boxes, the artifact patterns were significantly reduced compared to postoperative images. The cochlear structures that were slightly distorted by the artifacts (indicated by yellow arrows) were mostly recovered in comparison to the preoperative image slices. Finally, the MARGAN-generated images include by design, high intensity pixels at the potential locations of electrode centers. The yellow circles are clearly positioned in the centers of the electrodes and can help otologists visualize the relative positions of electrodes with respect to the scala tympani.

Quantitative Comparison with other

Similar to and Yu, 2018), we compared our approach with three open source MAR algorithms: MAR with projection linear interpolated replacement (marLI) [START_REF] Kalender | Reduction of ct artifacts caused by metallic implants[END_REF], beam hardening correction (marBHC) (Ver- burg and Seco, 2012) and NMAR [START_REF] Meyer | Normalized metal artifact reduction (nmar) in computed tomography[END_REF]. The visual assessment of the different MAR algorithms is shown in Fig. 7. The MARGAN approach clearly outperforms the other three MAR methods in its ability to decrease the texture changes of artifacts and to generate an image similar to the preoperative image. All 33 postoperative images were processed by marLI, marBHC, and NMAR. Three global similarity indices, the root mean square error (RMSE), structural similarity index (SSIM) and peak signal to noise ratio (PSNR), were computed between the preoperative images and the MAR images generated by the three comparison methods and our proposed approach. These three indices are reported in Table 4 and capture the preservation of visible structures, the errors and the quality of the reconstructed images. Our method outperforms the other MAR methods for all three metrics (lowest RMSE and largest SSIM and PSNR). In Fig. 8, the same indices were computed for all patient #1 image slices to evaluate the spatial consistency of the reconstruction. Clearly the MARGAN approach exhibits the best performance, with a lower mean value and much lower variance. This can be explained by the fact that it is the only MAR algorithm working directly on 3D images.

Impacts of methodological contributions

We assess the importance of our methodological contributions by evaluating their impact on the generated MAR-GAN images when they are removed from the computational pipeline. More precisely, we consider the following two contributions:

• Retinex Loss When zeroing the Retinex scale factor = 0 (instead of setting = 5e-5) during the MAR-GAN training, only the loss term is used, which is equivalent to minimizing the 2 norm between the generated and ground truth images. We also include in the ablation study the replacement of with the • Simulation of scatter and electronic noise in artifact simulation We simulated the image training set with only the beam hardening effect (as in [START_REF] Wang | Deep learning based metal artifacts reduction in post-operative cochlear implant ct imaging[END_REF]) or with the full pipeline as described in section 2.1.

1 norm involving | - ( )| terms.
In Table 5, we used the three similarity measures PSNR, RMSE and SSIM with respect to the preoperative images as a way to quantify the impact of those contributions.

Table 5 shows that both the addition of scatter and electronic noise in the simulation and the addition of the Retinex loss can improve the performance of MAR for all three different metrics. We also see that using a single 1 loss function performs worse than the proposed loss combination approach. A visualization of the image difference output obtained using different training loss functions is shown in Fig. 10 with subtraction maps between different output images and the ground truth image. We see from the yellow and red marks in those subtraction maps the effectiveness of the proposed Retinex loss function in comparison with using pure L1 and L2 losses.

Out-of-sample Test

To assess the generalization ability of this MARGAN approach, we explore its performance on 8 postoperative CBCT images, noting that the network was trained on CT images.

In Fig. 11, we see that metal artifacts in CBCT are more extensive and complex than in CT images. Yet, the MAR-GAN can cope well with those CBCT images and is able to recover most of the cochlear structures.

CI Electrode Position Prediction

The positioning of CI electrodes in postoperative imaging provides important information for establishing a hearing prognosis [START_REF] Kós | Measurements of electrode position inside the cochlea for different cochlear implant systems[END_REF][START_REF] Todt | Evaluation of cochlear implant electrode position after a modified round window insertion by means of a 64-multislice ct[END_REF] and can be used to improve the cochlear implant programming strategy [START_REF] Noble | Clinical evaluation of an image-guided cochlear implant programming strategy[END_REF][START_REF] Aebischer | Intraoperative impedance-based estimation of cochlear implant electrode array insertion depth[END_REF]. The proposed MAR-GAN algorithm output images where the electrode centers are outlined by voxels in hypersignal as shown in Fig. 9 and 7. To qualitatively evaluate the positional accuracy of those electrode centers in generated MARGAN images, we use pictures of the cochlea acquired after the dissection and grind-

DISCUSSION

Our MARGAN approach combines an artifact simulation pipeline with a 3D GAN network that generates augmented preoperative images from postoperative images. The artifact generation algorithm relies on three physical phenomena: beam hardening, scatter and electronic noise. The scatter and noise effects clearly have less impact on the output image compared to beam hardening. Yet, these effects were shown in Table 5 to improve the realism of the output of MARGAN when compared to preoperative images. The simulation pipeline could easily be refined in many ways, for example, using a more hardware-specific energy spectrum, increasing the number of sample energies in the approximation, or including more application-dependent scatter to primary ratios. This approach could also be extended to other imaging modalities, such as cone beam CT, dual energy CT or trimodal low-dose X-ray tomography [START_REF] Zanette | Trimodal low-dose x-ray tomography[END_REF] by introducing suitable artifacts simulation pipelines. The use of 3D GANs allowed us to generate MAR images with spatial coherence across neighboring slices, which is not guaranteed when using 2D slice-by-slice MAR methods. Furthermore, Retinex loss was introduced to improve the sharpness of the MAR images. We show in Tab. 5 that the Retinex loss can improve the performance of the MAR-GAN with a scale coefficient = 5e-5. However, an inappropriate value can introduce distortions in the MARGAN output. Furthermore, the influence of other hyperparameters in the simulation pipeline on the artifact reduction needs to be further investigated.

The MARGAN approach is both data driven (for the generation of MAR images) and model driven (for the generation of training image pairs). This is in contrast to purely data-driven MAR methods that either rely on pairs of preand postoperative images [START_REF] Jia | Automatically segmenting the left atrium from cardiac images using successive 3d u-nets and a contour loss[END_REF] or on nonpaired data [START_REF] Nakao | Regularized three-dimensional generative adversarial nets for unsupervised metal artifact reduction in head and neck ct images[END_REF]. The collection of image pairs, with and without artifacts, is mostly restricted to images acquired before and after an intervention such as CI insertion. The use of such pairs makes the 3D GAN fairly effective at removing artifacts in postoperative images. However, the collection of those images may be difficult and an intra-patient image registration is required. The MARGAN pipeline aims at reaching the same efficiency but by replacing the postoperative image with a simulated one. This makes the MARGAN algorithm applicable to a larger set of clinical cases where such image pairs cannot be gathered, for instance in the case of hip, shoulder or knee prostheses. The use of CycleGAN on non-paired images as in [START_REF] Nakao | Regularized three-dimensional generative adversarial nets for unsupervised metal artifact reduction in head and neck ct images[END_REF] is very appealing, because it avoids both artifact simulation and collection of paired images. However, it has only been tested to remove large artifacts, such as those caused by dental fillings, and with limited quantitative assessment.

Another advantage of the artifact simulation approach in MARGAN is its ability to augment the generated MAR image with voxels indicating the location of the metal part. In the case of CI postoperative images, it enables visualization in the same image of both the cochlea and the implant electrode centers. Note that the augmentation of the MARGAN image is only optional in this framework, because the metalfree image can replace the augmented image as train in the loss function of the 3D GAN.

Specifically, in the cochlear metal artifact reduction problem, we see from Fig. 7 and Tab. 4 that almost all the traditional MAR approaches have degradation problems in terms of reconstruction image quality. It was reported in [START_REF] Meyer | Normalized metal artifact reduction (nmar) in computed tomography[END_REF] and [START_REF] Diehn | Ct dental artifact: Comparison of an iterative metal artifact reduction technique with weighted filtered back-projection[END_REF] that sinogram inpaintingbased methods can introduce new artifacts. These artifacts can have a severe impact on image quality if the metallic parts and artifacts occupy a large area of space in the image, which is typically the case for the CI electrodes discussed here. However, the risk of quality degradation is not applicable for MARGAN as the image domain-based methods do not need to access the sinogram and the Radon transform.

A limitation of MARGAN lies in the relative complexity of implementing the simulation of metal artifacts in CT images. This is especially true for the scatter effect, which only adds a marginal gain in realism to the generated images. A more thorough study should be performed to evaluate the level of realism required in the simulation pipeline to improve the MARGAN output. For example, the real artifacts and their simulation images can be used as a testing set for evaluating the similarity between the realistic and the simulation. The metal artifacts that appear on the images are typically non-linear phenomena, which lead to the difficulty in designing a proper similarity metric for quantification. It is unclear yet how realism the simulated artifacts need to be in order to train the MARGAN to learn the metal artifacts reduction efficiently, although the experiments show that the MARGAN have the ability to generalize the metal artifacts reduction. Simulating the insertion of metal parts can also be complex as it requires a segmentation algorithm to locate the region of insertion. But this complexity is rewarded by the ability to generate a vast training set accounting for variations in patient anatomy or implant design.

Learning-based MAR methods were shown to outperform traditional sinogram-based MAR approaches in several previous works (Wang et al., 2019a;[START_REF] Zhang | Convolutional neural network based metal artifact reduction in x-ray computed tomography[END_REF][START_REF] Wang | Deep learning based metal artifacts reduction in post-operative cochlear implant ct imaging[END_REF]. But by design, the performance of those supervised methods depends on the chosen training set and they are application-specific algorithms. Their integration into a clinical workflow remains to be demonstrated, in particular due to their potential lack of robustness. For instance, we rely on a specific type of electrode array for the CI electrode positioning simulation. It would be of interest to evaluate the impact of CI characteristics (e.g. diameter, material, electrode spacing etc.) on the generation of MAR images. Yet, the successful application of MARGAN on CBCT images unseen during training is an encouraging sign of the generalization ability of MARGAN, though further studies are required.

Finally a limitation common to all MAR methods is the difficulty of evaluating performances quantitatively, due to the lack of ground truth data. The use of paired pre-and postoperative image data enables quantitative comparison through global similarity indices (such as PNSR, RMSE) but is also dependent on the registration quality of the two images. Images with synthetic artifacts created by image processing were also considered in [START_REF] Nakao | Regularized three-dimensional generative adversarial nets for unsupervised metal artifact reduction in head and neck ct images[END_REF], for instance, but they are computationally intensive to reach sufficient realism. Physical anthropomorphic phantoms are a useful alternative for MAR assessment [START_REF] Bolstad | Metal artifact reduction in ct, a phantom study: subjective and objective evaluation of four commercial metal artifact reduction algorithms when used on three different orthopedic metal implants[END_REF] but are limited by the number of phantoms considered.

CONCLUSION

In this paper, we have introduced a simulation-based 3D GAN to attenuate metal artifacts in CT images. The network is trained on a thousand regular CT images without any artifacts and their corresponding images where metal artifacts have been simulated. We have demonstrated the introduction of scatter and electronic noise effects in addition to beam hardening in an efficient computational pipeline. The complexity of scatter simulation has been alleviated by precomputing the impact of scatter on a generic head phantom where metal parts have been introduced. A Retinex loss was introduced to enhance visible edges in the generated images. The MARGAN approach was evaluated on CT and CBCT images of the inner ear with cochlear implants inserted. The proposed approach provided images close to preoperative images and outperformed open source MAR methods. Furthermore, images generated by MARGAN included the location of the electrode centers, which is useful for assess-ing the quality of implant surgery.

The trade-off between the complexity of artifact simulation and MARGAN output requires additional study, and we will also investigate the impact of MARGAN images on the automatic registration of pre-and postoperative images.

Figure 1 :

 1 Figure 1: Sketch of the MARGAN algorithm applied on postoperative images

Figure 2 :

 2 Figure 2: Cochlear implant electrode positioning simulation; (I) Registration of CT image on a template image; (II) Cochlear shape fitting; (III) Signed distance map generation; (IV) Electrode positioning; (V) Image fusion with electrodes.

Figure 3 :Figure 4 :

 34 Figure 3: The framework of MARGAN for metal artifact reduction. (a) The cochlear implant positioning simulation; (b) CI metal artifact physical simulation. (c) A 3D GAN is trained with simulated and preoperative datasets. The generator outputs an image with reduced metal artifacts. The discriminator network aims to identify whether or not the input image includes metal artifacts.

Figure 5 :Figure 6 :

 56 Figure 5: Pipeline of metal artifact simulation. Given a preoperative image with simulated implants, the simulation starts from the computation of attenuation maps (steps I -II) for the cochlea ROI volume based on the energy spectrum of the X-ray tube. Step (III) performs fan-beam projection to simulate the sinograms of the attenuation map. (IV) Monte Carlo simulation of scattering effects is performed offline on a head phantom for the generation of the scattering sinograms whose traces in the ROI are randomly chosen, then normalized and added to the combined attenuation map sinograms. (V) Gaussian electronic noise is added and then inverse fan-beam projection is performed to get the final simulated artifact images.

Figure 7 :

 7 Figure 7: Metal artifact reduction visualization of MARGAN in comparison with other approaches for patients #1 -5.

Figure 8 :

 8 Figure 8: The 3D consistency between slices from patient #1 for three different metrics. We see the MARGAN algorithm achieves the best slice consistency in comparison to other approaches.

Figure 9 :

 9 Figure 9: Results from patients #2 (top left), #3 (top right), #4 (bottom left) and #6 (bottom right) for two middle slices (first and second rows). The four columns correspond to: original postoperative images, output of MARGANs, registered preoperative images with manually positioned electrodes in red and postoperative images with electrodes appearing in yellow.

Figure 10 :

 10 Figure 10: Qualitative ablation study of Retinex loss effectiveness. The first column is a middle slice of patient #1, the second column is the corresponding outputs of MARGAN with different loss functions and the last column shows subtraction maps between the first two columns.

Figure 11 :

 11 Figure 11: Performance of MARGAN on 8 CBCT postoperative images. The yellow box shows three views of postoperative images and MARGAN-processed images for patient #1.

Figure 12 :

 12 Figure 12: Evaluation of the electrode position after the application of the MARGAN algorithm on 2 subjects (top and bottom); (a) Reformat of a 3D MARGAN image along a plane orthogonal to the modiolar axis. Red circles were manually added at the location of high intensity voxels; (b) Image of the cochlea with electrodes inserted after dissection and grinding of the temporal bone; (c) fusion of images (a) and (b) after an affine transform based on the manual correspondence of the centers of the two circles outlined by black squares. A good overlap of green and red squares is observed. Note that some inner circle electrodes shown in (c) are not visible because they are displayed from a two-dimensional cross-sectional plane whereas the bone grinding photography is in 3D top view.

Table 1

 1 Summary of major MAR approaches. (In the dataset collection column: BH, SC and EN indicate Beam Hardening, Scattering and Electronic Noise, respectively)

		Processing Domain Inner Ear MAR 2D/3D	Dataset Collection	Quantitative evaluation on clinical data
	marBHC	Sinogram	NO	2D	Non-Learning	YES
	marLI	Sinogram	NO	2D	Non-Learning	YES
	NMAR	Sinogram	NO	2D	Non-Learning	YES
	CNN Prior

Table 2

 2 Material Mapping Table for Voxel Conversion to MCGPU File

		air	water bones muscle titanium soft tissue fat
	MC-GPU MATERIAL 1	15	4	2	16	3	6
	DENSITY [g/cm 3 ]	0.001205 1.000 1.990 1.041	4.506	1.038	0.916

Table 3

 3 Dataset Summary: Preoperative and postoperative refer to images collected before and after Cochlear Implant, respectively.

	Dataset	Pre-Op Post-Op Photography
	Training	800	0	0
	Validation	200	0	0
	Evaluation CT	33	33	33
	Evaluation CBCT 0	8	0

Table 4

 4 Quantitative evaluation of the MARGAN approach compared to marBHC, marLI and Nmar. It shows the performance gain of MARGAN compared to other methods.

	Metric Preoperative marBHC marLI Nmar MARGAN
	PSNR	16.33	11.59	16.53 13.58	18.31
	RMSE	0.15	0.28	0.15	0.52	0.12
	SSIM	0.58	0.56	0.55	0.52	0.64

Table 5

 5 Ablation Study of Retinex and Physical Simulation

	Dataset	PSNR RMSE	SSIM
	MARGAN L1 Scatter	16.67	0.1490 0.56
	MARGAN L2 Scatter	18.17	0.1257 0.64
	MARGAN L2+Retinex No-Scatter 18.02	0.1277 0.63
	MARGAN L2+Retinex Scatter	18.31	0.1242 0.64

https://www.oem-xray-components.siemens.com/x-ray-spectrasimulation

⋆ This work was partially funded by the regional council of Provence Alpes Côte d'Azur, by the French government through the UCA JEDI and 3IA Côte d'Azur "Investments in the Future" project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01 and ANR-19-P3IA-0002, and was supported by the grant AAP Santé 06 2017-260 DGA-DSH.