Pierre-Jean Barjhoux
email: pierre-jean.barjhoux@irt-saintexupery.com

Youssef Diouane
email: youssef.diouane@isae.fr

Stéphane Grihon
email: stephane.grihon@airbus.com

Dimitri Bettebghor
email: dimitri.bettebghor@onera.fr

Joseph Morlier
email: joseph.morlier@isae.fr

A Bilevel Methodology for Solving a Structural Optimization Problem with both Continuous and Categorical Variables

In industrial structural optimization problems, two kinds of variables are involved : continuous variables, like geometrical parameters, and categorical variables like choices of material and cross-sectional shapes. A key difficulty of this kind of problem is the sensitivity of the computation cost with respect to the number of categorical design variables. In this article, a new methodology is proposed to address this mixed continuous-categorical structural optimization problem. For this purpose, a 10-bar truss test case representative for industrial design is introduced. Promising results have been obtained in terms of computation cost, compared to approaches based on combinatorial algorithms.

Nomenclature n

= number of structural elements or bars in the structure p = number of catalogs available for each element of the structure M n, p (R)

= set of matrices of size n × p with real coefficients B n, p = set of matrices of size n × p with binary coefficients that describe choices 0 n, p = set of matrices of size n × p with null coefficients Γ = a set of catalogs, i.e.; the list of available combination of choices for each element, Γ = {1, . . . , p} a = vector of areas of all elements or bars, a ∈ R n t = vector of thicknesses of all elements or bars, a ∈ R n c = categorical variables, ie a choice of catalog for each bar, c ∈ Γ n B = matrix of binary variables, B ∈ B n, p B i = i th line of the matrix of binary variables B, B i ∈ B 1, p Φ = vector of internal forces in the structure, Φ(a, t, c) ∈ R n m = number of stress constraints per element s = stress constraints in initial problem description, s(a, t, c) ∈ M n,m (R) g = stress constraints after problem reformulation, g(a, t, B) ∈ M n,m (R) w = weight of the whole structure, w(a, c) ∈ R

I. Introduction

In the field of structural optimization, two kinds of design variables are involved. Continuous or sizing variables describe the size of aircraft structural parts: in the case of thin-sheet stiffened sizing, they represent panel thicknesses and stiffening cross-sectional areas. On the other hand, technological choices (e.g. material choices, composite stacking sequences choices, cross-section shapes) or even computation assumptions (e.g. buckling margin policy) are depicted by categorical variables. These categorical variables are non-ordered. For example there is no ordering between combinations of material and cross-section shapes choices.

The challenge of this kind of optimization problem lies in the mixed continuous-categorical nature of the design space. More precisely, the large dimension of the categorical design space makes this kind of problem very hard to solve. Because of the curse of dimensionality, existing generic method solving Mixed Integer Non-Linear Programming (MINLP) are not efficient enough to face this problem in industrial context.

In this paper, details on related work are provided in a first place. Then the formulation of the optimization problem, as encountered in the industry, is described. It is followed by a description of the new approach proposed in this article. This approach is based on the decomposition of the mixed structure optimization problem in two levels: one where the sizing variables are optimized at given categorical choices, and a master problem that handles the categorical variables formulated as integers. Finally, an adapted 10-bar truss mixed structural optimization test case is presented and the results obtained are detailed. Promising results have been obtained compared to an hybrid branch and bound algorithm [START_REF] Barjhoux | Mixed Variable Structural Optimization: toward an Efficient Hybrid Algorithm[END_REF] in terms of computation cost. Indeed, the number of optimizations is proportional to the number of elements times the number of available choices per element, instead of a computation cost proportional to the number of elements at the power of the number of choices.

II. Related work

To handle this kind of problem, mixed-variable optimization algorithms have been proposed, such as Genetic Algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF], Differential Evolution [START_REF] Storn | Differential Evolution A Simple and Efficient Heuristic for global Optimization over Continuous Spaces[END_REF], Particle Swarm Optimization [START_REF] Kennedy | Swarm Intelligence[END_REF], Ant Colony Optimization [START_REF] Liao | Ant Colony Optimization for Mixed-Variable Optimization Problems[END_REF], and Pattern Search [START_REF] Torczon | On the convergence of pattern search algorithms[END_REF]. Some algorithms rely on continuous relaxation, or consider discrete variables as ordered [START_REF] Lampinen | Mixed Integer-Discrete-Continuous Optimization, by Differential Evolution, Part 1: the optimization method[END_REF][START_REF] Guo | Swarm intelligence for mixed-variable design optimization[END_REF][START_REF] Dimopoulos | Mixed-variable engineering optimization based on evolutionary and social metaphors[END_REF]. In other cases, categorical variables are handled through a native mixed-variable optimization approach for both categorical and continuous variables without relaxation [START_REF] Deb | A Flexible Optimization Procedure for Mechanical Component Design Based on Genetic Adaptive Search[END_REF][START_REF] Audet | Pattern Search Algorithms for Mixed Variable Programming[END_REF]. Other strategies consist in turning categorical variables into continuous ones, allowing the use of continuous optimization algorithms. For example, continuous indices containing shape and material information are built [START_REF] Rakshit | Simultaneous material selection and geometry design of statically determinate trusses using continuous optimization[END_REF], meta-models are constructed through Multiple Kernel Regression [START_REF] Herrera | Metamodel-assisted optimization based on multiple kernel regression for mixed variables[END_REF]. Other approaches rely on the structure of the problem, and suggest to handle discrete and continuous variables in two different problems, for example in [START_REF] Geoffrion | Generalized Benders decomposition[END_REF][START_REF] Duran | An outer-approximation algorithm for a class of mixed-integer nonlinear programs[END_REF]. Decomposition theory for multidisciplinary optimization has been adapted to mixed integer quasiseparable subsystems [START_REF] Haftka | Decomposition theory for multidisciplinary design optimization problems with mixed integer quasiseparable subsystems[END_REF] and then applied to structural optimization with discrete variables [START_REF] Schutte | Decomposition and Two-level Optimization of Structures with Discrete Sizing Variables[END_REF]. However, in most of the existing approaches that can handle MINLP problems, the algorithms are not able to solve a mixed structural optimization problem at an industrial scale. Furthermore, in many cases the algorithm relies on an ordering or relaxation of the discrete variables, that is not straightforward in the case presented in this article.

In industry, structural optimization software solutions handle continuous variables while categorical ones are often ignored. In fact, discrete variables always introduce a combinatorial behavior leading to an exponential number of configurations to be explored. This makes difficult to maintain the optimization performance while managing all variables types simultaneously. For instance, at Airbus, two structural optimization software were developed to split up the problem. The first one is inspired by classical continuous optimization approaches using gradients. The second solver, PRE STO [START_REF] Grihon | PRESTO: A rapid sizing tool for airframe conceptual design studies[END_REF], focuses on discrete and categorical variables in a preliminary phase of the design process where trade-offs are evaluated [START_REF] Grihon | Structure Sizing Optimization Capabilities at AIRBUS[END_REF]. In PRE STO, the overall optimization process alternates between local (element-wise) optimizations at given internal loads, and an internal loads update at given choices by Finite Element Model (FEM) analysis. The advantage of this method lies in the fact that the categorical variables are considered independent from each other during the element-wise optimization process. This makes this approach very efficient, and affordable in terms of computation costs even for industrial applications. However, one of the drawbacks is that since the optimizations are not performed on the whole structure, it is not possible to take into account global constraints like rigidity constraints. Furthermore, this optimization process is driven by local choices. In [START_REF] Barjhoux | Mixed Variable Structural Optimization: toward an Efficient Hybrid Algorithm[END_REF], the problem is decomposed in two levels. In a master level, the categorical variables are handled. In the slave problem, the continuous sizing variables are optimized, the choices being fixed. The categorical variables, at the upper level, are handled through a branch and bound algorithm. Usually, branch and bound lower bound computations rely on convex relaxation of the discrete variable. Since the non-ordered nature of the categorical variables does not allow to relax the categorical variable, an innovative formulation of a sub-problem for lower bound computations has been proposed. Thus, in this approach, the optimization is performed on the whole structure ; the FEM computations are nested into the optimizations. This approach can thus handle displacements constraints. However, the number of optimizations remains prohibitive for industrial applications.

III. Problem description

Let us name a the areas where a ∈ R n , and t the thicknesses where t ∈ R n . The set Γ is an enumerated set that contains all combinations of choices, for example the combination of choices of a material and the stiffener type of an element. To each combination of these choices is attributed an indice, that is the value of the catalog belonging to Γ. If p is the number of catalogs, we can note Γ = {1, . . . , p}. The categorical variables c can take values of the set Γ n (number of categorical choices per element at the power of the number of elements). Each value of this set is known as a catalog. In the industrial case, the structural model counts around 100 elements and from 10 up to 100 catalogs. This means that the continuous design space counts 200 continuous variables, and the categorical design space counts 100 variables. Furthermore, the categorical design space Γ n counts 100 10 up to 100 100 combinations describing the choices for the whole structure. This high combinatory dimension enforces the need for a methodology to solve efficiently such problems. The optimization problem is formulated as follows:

minimize c ∈Γ n ,a ∈R n ,t ∈R n w(a, t, c) subject to s(a, t, c) ≤ 0 d(a, t, c) ≤ 0 a ¯≤ a ≤ ā t ¯≤ t ≤ t (1)
with a ¯, t ¯, the lower bounds on areas and thicknesses, and ā, t their upper bounds. In this formulation, the constraints s ensure the structural strength element per element. It can be noted that s(a, t, c) is a matrix with n rows and m columns corresponding to the type of constraints per element. The topology of the structure to optimize is assumed to be fixed.

Furthermore, it is essential to describe the dependency of the constraints to state variables like internal loads Φ, and design parameters (a,t,c):

s(a, t, c) = Constraint type
Note that the continuous and categorical variables (i.e.; areas a, thicknesses t and categorical variables c hiding material or stiffening principle choices) have a significant role in this problem. The categorical variables affect the weight w, internal forces Φ, rigidity constraints d, stress constraints s. On the other hand, continuous variables affect the weight, internal loads, stress constraints and rigidity constraints. It is worth to note that a change in a categorical variable or area will modify the loads distribution Φ along the structure. Since the stresses s require the value of Φ, each component of s vector depends on the whole structure description. Internal forces Φ are computed using a finite element model (FEM). During this process the linear elasticity equation is solved:

Ku = F (3)
with K ∈ M q,q (R) the stiffness matrix of the considered structure, u ∈ R q the vector of displacements, F ∈ R q the vector of forces applied on the nodes. q is the dimension of the space times the number of nodes. Then, the strain is computed from the displacements u, and the constraints s from the Hooke's law.

The rigidity constraint d on displacements ensures that at one node of the truss, the displacement will not exceed a predefined value.

IV. A bi-level approach

To tackle the problem (1) at an industrial scale, a reformulation is proposed. Let B n, p describe a set of matrices whose coefficients have integer values. Each line of the matrices has only one binary coefficient equal to 1, the other coefficients are 0. Physically, each line of the matrix describes a choice of catalog applied to a given element, and only one choice can be made for each element. Let be g the function that outputs a matrix g(a, t, B) ∈ M n,m (R) of stress constraints. B i is the i th line of the matrix B. The expression of the coefficients of this matrix is thus:

B n, p = X = (x i j ) 1≤i ≤n 1≤ j ≤p x i j ∈ {0,
g i j (a, t, B) = p k=1 B ik s i j (a i , t i , B i γ, Φ(a, t, Bγ)) ∀(i, j) ∈ 1, n × 1, m (4)
The mixed variable optimization problem can be written as follows : minimize

B ∈B n, p a ∈R n t ∈R n w(a, t, B) subject to g(a, t, B) ≤ 0 d(a, t, B) ≤ 0 a ¯≤ a ≤ ā t ¯≤ t ≤ t (5)
The weight w is defined as:

w(a, t, B) = n i=1 p j=1 B i j ρ j a i L i (6)
with the vector of densities ρ ∈ R p , corresponding to the density associated to each catalog described in γ, in the same ordering.

A. Algorithm 1

This problem is based on a decomposition in two optimization problems. The slave problem consists in a continuous optimization with respect to continuous variables and relaxed components of the categorical variables. The optimal value of the objective of this problem is then linearized with respect to the categorical variables, and optimized in a master problem that computes the optimal categorical choices. This work is still in progress (see Appendix) and it is preferred to provide in this paper a detailed Bi-level methodology with new results.

B. Algorithm 2

In this approach, the post-optimal sensitivity analysis of the slave problem in the previous algorithm is replaced by a computation of the slope of the optimal Lagrangian between two discrete choices.

The problem ( 5) is decomposed in two problems, (P a,t ) and (P B ). The choices are optimized at the master level (P B ), relying on optimization with respect to a,t in the slave problem (P a,t ), where B is a parameter. The slave problem (P a,t ) consists thus in an optimization with respect to continuous design variables that are thicknesses and areas (resp. t, a), while the binary variable is fixed at a given B:

(P a,t ) B minimize a,t ∈R n w(a, t, B) subject to g(a, t, B) ≤ 0 d(a, t, B) ≤ 0 a ¯≤ a ≤ ā t ¯≤ t ≤ t (7)
This corresponds to a sizing step, where the choices of materials and shape are fixed. At the first step (0) of the algorithm, the problem (P a,t ) B (0) is solved at a given B (0) . Let be w * (B) the optimal weight of (P a,t ) B taken at optimal areas and thicknesses, respectively a * (B) and t * (B).

Let be L(B) the value of the Lagrangian of the problem (P a,t ) B , i.e.;

L(B) = w(a(B), t(B), B) + n i=1 m j=1 λ g i j (B) g i j (a(B), t(B), B) + λ d (B)d(a(B), t(B), B) + n i=1 λ ub a i (B)(a i -a ¯i) + λ lb a i (B)(a ¯i -a i ) + λ ub t i (B)(t i -t ¯i) + λ lb t i (B)(t ¯i -t i ) (8)
Let be g * (B) and d * (B) the constraints (on stress and displacement, resp.) at the optimum, and w * (B) the optimal weight taken at the optimum of (P (a,t) ), such that:

g * (B) = g(a * (B), t * (B), B), d * (B) = d(a * (B), t * (B), B), w * (B) = w(a * (B), t * (B), B). (9)
The associated optimal Lagrangian L * (B) is then given by:

L * (B) = w * (B) + n i=1 m j=1 λ g i j (B) g * i j (B) + λ d (B)d * (B) + n i=1 λ ub a i (B)(a * i (B) -a ¯i) + λ lb a i (B)(a ¯i -a * i (B)) + λ ub t i (B)(t * i (B) -t ¯i) + λ lb t i (B)(t ¯i -t * i (B)) . (10)
A linearization of L * around B (0) is given by:

L * (B) = L * (B (0) ) + n i=1 p j=1 ∆L * (B (0) ) ∆B i j i j (B i j -B (0) i j ) (11)
with :

• ∆B i j a matrix of perturbation of B (0) on the coefficient B (0) i j , • B (0) i j the binary design variable that drives the choice of the j th catalog on the i th element, at a given step (0) of the methodology,

• ∆L(B (0) ) ∆B i j i j
the coefficient i j of the matrix ∆L(B (0) ) ∆B i j .

At the master level, the linearization of the Lagrangian ( 11) of (P a,t ) B is minimized with respect to B :

(P B ) minimize B ∈B n L * (B) . (12)
Considering equation [START_REF] Audet | Pattern Search Algorithms for Mixed Variable Programming[END_REF], solving problem ( 12) is equivalent to find, for each element i, the binary coefficient B i j such that the discrete slope of the Lagrangian with respect to B i j is minimal. The problem ( 12) is thus equivalent to compute L * (B) so that :

B (1) i j = 1 if j = l i (∀i ∈ {1, . . . , n}) 0 else with l i = argmin j ∈ 1, p ∆L * (B (0) ) ∆B i j i j (∀i ∈ {1, . . . , n}) . (13)
The discrete slope of the Lagrangian ∆L * (B (0) ) ∆B i j kl is given by :

∆L * (B (0) ) ∆B i j kl =          L * (B (0) + ∆B i j ) -L * (B (0) ) ∆B i j kl if B (0) i j 1 0 else (14)
with ∆B i j ∈ M n, p (R), so that ∆B i j = 0 if B (0) i j

1, otherwise k ∈ {1, . . . , n} and l ∈ {1, . . . , p} :

∆B i j kl =          1 if (k, l) = (i, j) -1 if k = i and B (0) k,l = 1 . 0 else (15)
A numerical example of the computation of the perturbation ∆B i j , for a given B (0) is depicted in Fig. 2. According to KKT conditions, the Lagrangian at the optimum is strictly equal to the optimal value of the objective. Assuming that L * (B) is a good approximation of L(B) around B (0) :

B 0 = 1 0 0 0 1 0 ∆B 11 = 0 0 0 0 0 0 ∆B 12 = -1 1 0 0 0 0 ∆B 11 = -1 0 1 0 0 0
L * (B) = w * (B) . ( 16
)
According to equations ( 14) and ( 15), solving problem ( 13) is thus equivalent to find the binary variable such that, for each element, it minimizes the comparison between the optimal weight computed at a perturbed choice B (0) + ∆B i j , and the optimal weight at B (0) . The problem [START_REF] Herrera | Metamodel-assisted optimization based on multiple kernel regression for mixed variables[END_REF] is thus equivalent to compute L * (B) so that :

B (1) i j = 1 if j = l i (∀i ∈ {1, . . . , n}) 0 else with l i = argmin j ∈ 1, p w(B (0) + ∆B i j ) -w(B (0) ) (∀i ∈ {1, . . . , n}). (17)
Once B (1) has been computed, the next iteration of the algorithm starts again with the computation of the solution of the sizing problem (P a,t ) B (1) taken at B (1) . The full process of the algorithm is depicted in Algorithm 1.

Algorithm 1 Bi-level algorithm

1: B (0) ← initialize B initialize B 2: w hist , t hist , a hist , B hist ← list(), list(), list(), list()
initialize lists that store results history 3: w (0) , t (0) , a (0) ← solve (P a,t ) B (0) 4: add w (0) , t (0) , a (0) , B (0) to w hist , t hist , a hist , B hist 5: k ← 1 6: while |w (k)w (k-1) | < do check convergence on weight 7:

for each bar i from 1 to n do 8:

r ∈ R p 9:
for each catalog j from 1 to p do

10: if B (k-1) i j = 1 then avoids to repeat optimization at B (k-1)
11:

r j ← w (k-1)
12:

else 13:
r j ← solve (P a,t ) B (k-1) +∆B i j solve sizing taken at perturbed choice, r j is the optimal weight 14:

l ← argmin j ∈ 1, p r j find choice that minimizes weight 15:

B (k) i j = 1 if j = l i (∀i ∈ {1, . . . , n}), 0 otherwise build binary solution of iteration (k) 16:

w (k) , t (k) , a (k) ← solve (P a,t ) B (k)
sizing problem taken at current optimal choices (for all elements)

V. Test Case

In order to test the new methodology, the well-known 10-bar truss optimization test-case often introduced as a continuous optimization problem [START_REF] Merval | L'Optimisation Multiniveaux D'une Structure[END_REF][START_REF] Samuelides | Modèles réduits en optimisation multiniveau de structures aéronautiques[END_REF], has been adapted [START_REF] Barjhoux | Mixed Variable Structural Optimization: toward an Efficient Hybrid Algorithm[END_REF] to match the needs of a mixed optimization problem. The truss is illustrated in Fig. 3.

In this case, the sizing variables are the areas a of the bars cross-sections. Each s i j (a, Bγ) is computed as the difference between the stress in the current element i and allowable in tension, compression, and the computed values of local and Euler buckling. It is worth noting that the internal forces computation requires all the components of the variable B, and all areas a, i.e.; s(a, Bγ) =

Allowable Tension Allowable Compression Euler Buckling

Local Buckling

elt 1 s 11 (a 1 , B 1 γ, Φ(a, Bγ)) s 12 (a 1 , B 1 γ, Φ(a, Bγ)) s 13 (a 1 , B 1 γ, Φ(a, Bγ)) s 14 (a 1 , B 1 γ, Φ(a, Bγ)) elt 2 s 21 (a 2 , B 2 γ, Φ(a, Bγ)) s 22 (a 2 , B 2 γ, Φ(a, Bγ)) s 23 (a 2 , B 2 γ, Φ(a, Bγ)) s 24 (a 2 , B 2 γ, Φ(a,
The definition of the constraint g becomes, with B k the k th row of the matrix B:

g i j (a, B) = p k=1 B ik s i j (a i , B k γ, Φ(a, Bγ)) ∀(i, j) ∈ 1, n × 1, m (19)
Thus, in this test case the number of elements n is 10 and the number of stress constraints m is 4. The constraints g ∈ M 10, p are defined using equation (4). Here is the problem formulation :

B i j ρ j a i L i subject to g(a, B) ≤ 0 d(a, B) ≤ 0 100 ≤ a ≤ 1300 (mm) (20)
The sizing problem (P a ) is thus :

(P a ) minimize a ∈R n w(a, B) subject to g(a, B) ≤ 0 d(a, B) ≤ 0 100 ≤ a ≤ 1300 (mm) (21)
The implementation of this test-case relies on GEMS, a Python library for programing MDO processes [START_REF] Gallard | GEMS: A Python Library for Automation of Multidisciplinary Design Optimization Process Generation[END_REF]. Each elementary function like binary variable encoding, weight, internal forces and constraints computations has been catalog material shape implemented as a GEMS MDO discipline object. GEMS automatically chains all of them, depending on the design variables, objective and constraints declared to the GEMS MDO scenario object. The jacobian of each MDO discipline has been implemented, so that GEMS automatically builds the jacobian of objective and constraints. The gradients of objective and constraints have been checked with finite differences. Due to the quasi-separable structure of the problem, the implementation of the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | A class of globally convergent optimization methods based on conservative convex separable approximations[END_REF] from the Python package NLopt has been used to perform the continuous optimizations.

0 M 1 I 1 M 2 I 2 M 1 C 3 M 2 C

VI. Results

The catalogs definition, with respect to the case number, is detailed in Table 1. The characteristics of materials M 1 and M 2 are presented in Table 2.

The results obtained by the bi-level optimization can be compared to results obtained by 3 other approaches : PRE STO [START_REF] Grihon | Structure Sizing Optimization Capabilities at AIRBUS[END_REF], hybrid branch and bound [START_REF] Barjhoux | Mixed Variable Structural Optimization: toward an Efficient Hybrid Algorithm[END_REF] (h-B&B), and a complete enumeration. Each case is thus solved 4 times. All cases are different in terms of categorical design space Γ and bound on displacements. The complete enumeration is an enumeration of all continuous optimizations taken at all possible discrete choices. This is why for p = 2 and n = 10, it required 2 10 = 1024 optimizations on areas. This means that the complete enumeration that would solve problem [START_REF] Samuelides | Modèles réduits en optimisation multiniveau de structures aéronautiques[END_REF] requires 1024 different matrices B fixed. All the results are exposed in Table 3.

In cases 1 to 3 and 7 to 9, the bound constraint on displacements is active. These cases can not be solved by the PRE STO approach. On the other hand, cases 4 to 6 and 10 can be solved by a PRE STO approach since this constraint is not active (and can be omitted in the optimization problem). It is worth to note that an optimization during the PRE STO approach is less costly than an optimization with the other 3 methods. Indeed, with PRE STO the FEM is evaluated at each global iteration during the internal loads computation, i.e. outside element-wise optimizations. For example the case 4 required a total 200 optimizations but only 10 evaluations to the FEM, one per global iteration. First, it can be observed that the optimal choices and weight (and thus areas) computed by the 4 approaches in all cases are equal. Furthermore, the number of optimizations required by the Bi-level approach seems promizing. Indeed, it can be remarked that the total number of optimizations required to solve this problem with a complete enumeration is 1024 for p = 2. The number of optimizations required by the h-B&B approach varies from 72 to 252, i.e. around 7% to 25% of the computation cost of a complete enumeration. The Bi-level approach requires 33 optimizations (around 3% of the computation cost), and is converged after 3 global iterations. The efficiency of the Bi-level approach becomes more interesting when the number of catalogs is increased. Indeed, for cases 7 to 10, the Bi-level requires from 1.7% to 6.8% of the computations required by the h-B&B approach. The computation cost reaches 0.002% when it is compared to the complete enumeration. Thus, according to these experiments, the h-B&B computation cost is proportional to the number of catalogs p at the power of the number of n structural elements (p n ). The computation cost of the Bi-level is proportional to the number of catalogs times the number of structural elements (p × n). Furthermore, all the p × n optimizations performed at each global iteration of the Bi-level approach can be performed in parallel. With the h-B&B approach, p optimizations per branching step can be launched in parallel without modifying the overall methodology. The optimal weights obtained by the Bi-level approach are very close to those obtained with h-B&B and complete enumeration approaches. In particular for cases 4 to 6, the optimal catalog 0 is chosen at the 4 th bar instead of choice 1, leading to a higher weight. For cases 7 to 10, it can be remarked that even if optimal choices seem very different from the results obtained by the h-B&B approach, the optimal weights are very close. In the other hand, it is worth to note that when the number of catalogs p is increased, convergence difficulties can be observed with the Bi-level approach.

Cases

VII. Conclusion

In this work, a new methodology is proposed to solve a mixed categorical-continuous structural optimization problem. This methodology relies on the separation of the continuous variables from the categorical ones, that are treated in two different optimization problems. Compared to existing approaches like the hybrid branch and bound, the industrial approach PRE STO and a complete enumeration, the computation cost seems promising, and similar optimal weights are obtained. Furthermore, the whole process could benefit from the parallelization of all optimizations performed during the overall iterations. This observed efficiency is due to the linearization of the Lagrangian of the sizing optimization problem, making the choices of each structural element independent from each other. However, in a future work it will be interesting to assess the drawbacks of this linearization. In addition, another way to compute the discrete Lagrangian slope, based on post-optimal sensitivity definition, is under investigation. Furthermore, it would be interesting to scale up the model and increase the design space dimension in order to check the behavior of the convergence. Finally, it has been observed that in some cases, this current formulation of the Bi-level algorithm retrieves a higher optimal weight than a complete enumeration. This is why another approach, based on the framework AMIEGO [START_REF] Roy | Next generation aircraft design considering airline operations and economics[END_REF], is under investigation. The main idea of this framework is to leverage the efficiency of a gradient-based optimizer while globally exploring the integer design space with the help of the Efficient Global Optimization (EGO) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] algorithm.

Appendix

In this first version of the bi-level approach, it is proposed to apply a relaxation on binary variables in problem [START_REF] Liao | Ant Colony Optimization for Mixed-Variable Optimization Problems[END_REF]. The relaxed binary variable is afterwards ensured to be discrete thanks to a particular formulation of the master optimization problem. Let be C n, p such as : C n, p = X = (x i j ) 1≤i ≤n 1≤ j ≤p x i j ∈ [0, 1], and p j=1

x i j = 1, ∀i ∈ 1, n
The overall problem (5) becomes, with relaxed binary variables U belonging to C n, p : minimize

U ∈ C n, p (a,t)∈R 2n w(a, t, U) = n i=1 p j=1 U i j ρ j a i L i subject to g(a, t, U) ≤ 0 a ¯≤ a ≤ ā t ¯≤ t ≤ t ( 22
)
It can be remarked that solving problem ( 5) is equivalent to solve problem [START_REF] Gallard | GEMS: A Python Library for Automation of Multidisciplinary Design Optimization Process Generation[END_REF] on the edges of C n, p : Problem ( 5) is the restriction of problem ( 22) on B n, p .The constraint on displacements d is not considered in this first approach. Now that the problem has been defined, the bi-level algorithm can be described. First, the variable U is decomposed into two variables. A variable α drives the local choices through the limit stress computation and the choice of density, while β drives the global effects through the choices involved in the internal loads computations. Both α and β belong to C n, p . Let be E ∈ M p, f a matrix of material properties. On each column are listed the p values of one material property. There are f material properties, including for example Young modulus and Poisson coefficient. Since the algorithm will ensure that α never takes intermediate values, only β can take intermediate values between 0 and 1. This is why β can not directly multiply the vector of choices γ, but the matrix of continuous properties E. An example of definition of E is given in Fig. 4. The constraint formulation (4) is thus rewritten as:

g i j (a, t, α, β) = p k=1 α ik s i j (a i , t i , γ k , Φ(a, t, βE)) ∀(i, j) ∈ 1, n × 1, m (23)
and the weight definition is given by :

w(a, t, α) = n i=1 p j=1 α i j ρ j a i L i (24)
The slave problem (P a,t,β ), consists in an optimization with respect to the sizing variables a and t, and the variable β. The variable α is fixed.

(P a,t,β ) Given α : minimize (a,t)∈R 2n β ∈ C n, p w(a, t, α) subject to g(a, t, α, β) ≤ 0 d(a, t, α, β) ≤ 0 a ¯≤ a ≤ ā t ¯≤ t ≤ t ≤ α -β ≤ (25)
This corresponds to a sizing step, where the choices of materials and shape are fixed. Let be L(α) the value of the Lagrangian of the problem P (a,t,β) , i.e.;

L(α) = w(a(α), t(α), α) + n i=1 m j=1 λ g i j (α) g i j (a(α), t(α), α, β(α)) + n i=1 λ ub a i (α)(a i -a ¯i) + λ lb a i (α)(a ¯i -a i ) + λ ub t i (α)(t -t ¯i) + λ lb t i (α)(t ¯i -t i ) + n i=1 p j=1 λ ub α-β i j (α)(α i j -α ¯ij ) - n i=1 λ lb α-β i (α)(α ¯ij -α i j ) (26)
Let be L * the Lagrangian at the optimum of problem 25. A linearization of L * around α 0 is given by:

L * (α) = L * (α 0 ) + n i=1 p j=1 dL * (α 0 ) dα i j (α i j -α 0 i j ) (27) with :
• dL * (α 0 ) dα i j the coefficient i j of the jacobian of the Lagrangian taken at α 0 ,

• α 0 i j the relaxed binary design variable that drives the choice of the j th catalog on the i th element, at a given step 0 of the methodology. At the master level, the Lagrangian (27) of (P a,t,β ) is minimized with respect to α :

(P α ) Given a, t, β : minimize α∈ C n, p L * (α) (28)
Using equation ( 27), it can be demonstrated that solving problem (28) is equivalent to find, for each element i, the binary coefficient α i j = 1 such that the slope of the Lagrangian is minimal. The problem (28) becomes :

(∀i ∈ 1, n ), l i = argmin j ∈ 1, p dL * (α 0 ) dα i j with α * given by α * i j = 1 if j = l i (∀i ∈ 1, n ) 0 else (29)
The slope of the Lagrangian dL * (α 0 ) dα i j taken at α 0 can be computed using post-optimal sensitivity definition of problem [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF], i.e.;

dL * (α 0 )

dα i j = n k=1 ∂L * (α 0 ) ∂a k da k dα i j + n k=1 ∂L * (α 0 ) ∂t k dt k dα i j + n k=1 ∂L * (α 0 ) ∂ β k dβ k dα i j + n k=1 m l=1 ∂L * (α 0 ) ∂λ g kl dλ g kl dα i j + ∂L * (α 0 ) ∂λ d dλ d dα i j + ∂L * (α 0 ) ∂α i j (30)
Furthermore,

• the KKT conditions of problem [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] give

∀i ∈ 1, n ∂L * (α 0 ) ∂a i = 0, ∂L * (α 0 ) ∂t i = 0, ∂L * (α 0 ) ∂ β i = 0 (31)
• under hypothesis of Fiacco's theorem, a * (α), t * (α) and β * (α) are local solutions in the vicinity of α. Thus, a strong assumption is made: in a neighborhood of α the set of active constraints A(α) does not change, in particular no inactive constraint is activated (or the contrary) with a change of α :

∀(i, j) ∈ 1, n × 1, p , ∀(k, l) ∈ 1, n × 1, m dλ g k j dα i j = 0 (32)
Thus, under hypothesis of Fiacco's theorem, the expression 30 becomes:

dL * (α 0 ) dα i j = ∂L * (α 0 ) ∂α i j (33)

and the derivative of the Lagrangian, at the optimum, is:

dL * (α 0 ) dα i j = ∂w(a * (α 0 ), t * (α 0 ), α 0 ) ∂α i j +

The algorithm procedure is described in Algorithm 2:

Algorithm 2 Bi-level algorithm with post-optimal sensitivities 1: α (0) ← initialize α initialize α 2: w (0) , t (0) , a (0) , β (0) ← solve (P a,t,β ) α (0) 3: k ← 1 4: while |w (k)w (k-1) | < δ do check convergence on weight for each bar i from 1 to n do 7:

l i = argmin j ∈ 1, p ∂L(α (k) ) ∂α i j 8:
α (k) ← α (k-1) + ∆α il build binary solution 9:

w (k) , t (k) , a (k) , β (k) ← solve (P a,t,β ) α (k) 10:

k ← k + 1 11: k ← k-1 /2
Decrease the tolerancy between α and β 12: return w * , t * , a * , α * returns optimal weight, thicknesses, areas, and catalogs Thus in this algorithm, although α belongs to C n, p , the linearization of the Lagrangian leads to pure binary values. Theoretically, α belongs to C n, p , but the algorithm (in particular problem (29)) ensures that α coefficients are such that α belongs to B n, p . Later in the process, when problem (25) is solved, the relaxed variable β is forced to stay in a vicinity around α. This vicinity is decreased iteration per iteration, so that at the end α = β, with (α, β) ∈ B 2 n, p . However, numerical experiments show that the assumptions leading to equation (33) are not verified. Indeed, in practice, a change of the binary variable α induces a change in the set of active constraints. This is responsible of a non-linearity as illustrated in Fig. 5b. Fig. 5 Optimization with n = 1 and p = 2. This simplified case illustrates the importance of the behaviour of the Lagrangian between two solutions P 1 and P 2 . The post-optimal sensitivity can lead to wrong choices if, for example, there is a change in the set of active constraints (in Fig. 5b). In both situations, the optimal solution is α 1 = 1 (and thus α 0 = 0).

4 Fig. 1

 41 Fig. 1 Specific coding of the categorical design variable in the case of a 5 elements structure : n = 4. In this case Γ is a set of 4 catalogs : p = 4. Both B and c describe the same choices : choice of catalog 4 for bar 1, catalog 2 for bar 2, catalog 1 for bar 3, etc. Fig. a shows the categorical variables description used in the approach presented in Section III : c ∈ Γ 5 . Fig. b illustrates the binary coding B ∈ B 5,4 of choices c.

Fig. 2 A

 2 Fig. 2 A numerical example of the computation of ∆B 1j (∀ j ∈ 1, p ) with respect to B 0 , with n = 2, p = 3.

FFig. 3

 3 Fig. 3 Illustration of a 10-bar truss structure. Each member works in tension-compression only. A vertical load is applied on node 4. The length of horizontal and vertical bars is l = 1m.

minimizeB

  ∈B 10, p a ∈R 10 w(a, B) =

E = e 1 ν 1 e 2 ν 2 e 3 ν 3 e 4 ν 4 Fig. 4 2 .

 442 Fig. 4 Example of matrix of material properties, with p = 4 and f = 2. It provides the values of Young modulus (e) and Poisson coefficient (ν) in the same order than the list of choices γ.

g

  kl (a(α 0 ), t(α 0 ), α 0 , β(α 0 ))∂α i j -λ lb α-β i j (α 0 ) + λ ub α-β i j (α 0 ) = ρ j a * i (α 0 )L i + kl s k j (a * k (α 0 ), t * k (α 0 ), γ j , Φ(a * (α 0 ), t * (α 0 ), β * (α 0 ) E)) -λ

  (a) The Lagrangian is linear w.r.t. α and thus the post-optimal sensitivity returns a correct Lagrangian slope. (b) The misleading post-optimal sensitivity computed in P 2 makes the algorithm choose α 0 = 1, however the optimal weight in P 1 is greater than the one in P 2 . A change of constraints is responsible of a non-linearity in P 12 .

  n s n1 (a n , t n , c n , Φ(a, t, c)) s n2 (a n , t n , c n , Φ(a, t, c)) . . . s nm (a n , t n , c n , Φ(a, t, c))

		1	Constraint type 2	. . .	Constraint type m
	elt 1	s 11 (a 1 , t 1 , c 1 , Φ(a, t, c)) s 12 (a 1 , t 1 , c 1 , Φ(a, t, c)) . . . s 1m (a 1 , t 1 , c 1 , Φ(a, t, c))
	elt 2	s 21 (a 2 , t 2 , c 2 , Φ(a, t, c)) s 22 (a 2 , t 2 , c 2 , Φ(a, t, c)) . . . s 2m (a 2 , t 2 , c 2 , Φ(a, t, c))
	. . .	. . .	. . .	. . .	. . .

elt

Table 1

 1

			M 1	M 2
		density (kg/mm 3 ) 2.8 10 -6 2.7 10 -6
		young modulus (M Pa) 7.1 10 4	7.4 10 4
		poisson coefficient (-)	0.3	0.33
		tension allowable (M Pa) 1.6 10 2	1.5 10 2
	Catalogs definition with	compression allowable (M Pa) 2.1 10 2	2.0 10 2
	respect to case numbers, c ∈ Γ 4

Table 2 Numerical details on Material definitions

 2

Table 3 Number of optimizations, optimal weights, and optimal catalogs obtained for 10 different cases. Each case takes different values of bounds on displacements and/or number of catalogs.

 3

	Number of optimizations

Acknowledgments

This work is part of the MDA-MDO project of the French Institute of Technology IRT Saint Exupery. We wish to acknowledge the PIA framework (CGI, ANR) and the project industrial members for their support, financial backing and/or own knowledge: Airbus, Airbus Group Innovations, SOGETI High Tech, Altran Technologies, CERFACS.


Long Chen

Thomas Rodet

Nicolas Gac

A simple and efficient super-short-scan algorithm of Fan-beam reconstruction for multiple circular trajectories: solution towards the truncated data

The sufficiency condition of the accurate reconstruction of the region of interest (ROI) in fanbeam tomography with non-truncated data was introduced by F. Noo([1]). When the detector does not cover the whole objet, R. Clackdoyle extended this condition of the accurate reconstruction of ROI in presence of the truncated data using virtual fan-beam method (VFB) [2]. In this paper, we are interested in the image reconstruction of the whole object from truncated data, since the whole object reconstruction is always preferred even with small detector, instead of the accurate reconstruction of only a region. We propose a simple and efficient super-short-scan algorithm for multiple circular trajectories to reconstruct the whole object from truncated fan-beam projections. This algorithm is validated by simulation studies.

I. Introduction

In fan-beam tomography, the short-scan (180 • plus the fan angle) is a minimal complet data set to reconstruct the whole object [START_REF] Parker | Optimal short scan convolution reconstruction for fanbeam ct[END_REF]. This had been also hold even for the reconstruction of only a ROI [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF]. A new sufficient condition of the accurate reconstruction of ROI was introduced in [START_REF] Noo | Image reconstruction from fan-beam projections on less than a short scan[END_REF]. The ROI can be exactly reconstructed only and if only all the lines through the ROI are known. All the fan-beam projections should be non-truncated as the short scan.

When the detector is too small to cover all the object, the projection data is inevitably truncated. The accurate reconstruction of ROI is still able to be achieved in some cases with VFB in [START_REF] Clackdoyle | Quantitative reconstruction from truncated projections in classical tomography[END_REF]. An optimal virtual fan-beam circular trajectory is found by rebinning the truncated projections to non-truncated projections. Additionally, the whole object reconstruction from the truncated data has been interesting, as in dental imaging, the panoramic view of the whole jaw from a medium detector. In this paper, we propose a simple and efficient super-short-scan (SSscan) algorithm on twice and triple circular trajectories to achiev the reconstruction of the whole object. 'Simple' means the movement of the scanner is not complexe, and 'efficient' stands for as less fan-beam projections as possible to reconstruct the whole object. The super-shortscan projections here are not necessarily non-truncated, which are different from those in [START_REF] Kudo | New super-short-scan algorithms for fan-beam and cone-beam reconstruction[END_REF].

The rest of this paper is organized as follows: the proposed super-short-scans on twice and triple circular trajectories are described in section 2, we use a standard iterative least-squares method for image reconstruction in section 3. the proposed super-short-scan algorithm is evaluated by the simulation studis in section 4, and we draw the conclusion in the final section.

II. Super-short-scans on multiple circle trajectories

In CT tomography, how to reconstruct a whole object from the truncated data is a challenge, when the detector is too small to cover all the object. Obviously, a single circular scan is able to reconstruct only a region, not the whole object. More scans are needed to achieve the whole object reconstruction. We present the proposed supershort-scans on two and three circular trajectories below.

A. Twice super-short-scans

Usually, the object support is known a priori. Without loss of generality, let the object support Ω be an ellipse depicted by tow semi-axis parametres (a and b, a > b), which is centered at origine in figure 1.(a). Two iso-centers O 1 and O 2 of scans are located symetrically to the axis Y, with distance c to the origine. The system geometry is decribed in the figure 1.(a) with distance from the Xray source to the iso-center R, and from the source to the detector D. The length of detector is L. We have r = L * R/D, where r is the radius of the field of view (FOV). Obviously, r should be greater than c for the whole object reconstruction.

Intuitively, the first choice of the whole object reconstruction is to reconstruct a half ellipse from each scan. In presence of the truncated data, an optimum virtual trajectory of the X-ray is found in red solid circle using VFB ( [START_REF] Clackdoyle | Quantitative reconstruction from truncated projections in classical tomography[END_REF]), in the figure 2.(a). According to the sufficient condition in [START_REF] Noo | Image reconstruction from fan-beam projections on less than a short scan[END_REF], two reasonable trajectories are drawn in the figure 2.(a), the scanning angle of each scan Q 1 or Q 2 is π + 2 * γ 1 , where γ 1 = arcsin(c/R). As c < r, the angle is smaller than half of the fan angle γ m with γ m = arcsin(r/R). These tracjectories are refered to the reduced scan (R-scan). When we look at the two reduced scans, a random pixel M in the left half of object is mesured twice in one line in the figure 2.(a). The fanbeam projections of the whole object are double from the scanning range between the lines t 1 and t 2 and between t 3 and t 4 , which are the shared tangents of the object support Ω and two fileds of view (FOV) of the twice scans.

λ1 s = π + ζ 1 , λ1 e = κ 1 = ξ + τ 1 ( 1
)
λ2 s = -ζ 2 , λ2 e = π -κ 2 = π -(ξ + τ 1 ) (2)
where the tangent angle ξ is given by

tanξ = 1/ ( r * c + (a 2 -b 2 ) * (r 2 -b 2 ) + b 2 * c 2 r 2 -b 2 ) 2 -1,
which is derived from the solution of the tangent to the FOV of Q 1 or Q 2 and the object support function.

Compared to the reduced scan, the SS scan requires less fan-beam projections. Quantitively, the reduced projection angle is given below :

∆ Π = π + 2 * γ 1 -|λ1 e -λ1 s | = 2 * γ 1 + ξ + τ 1 -ζ 1 (3)
A numerical result will be given in the simulation studies, ∆ Π = 30 • .

B. Triple super-short-scans

When the object is much larger than the field of view of the scanner, triple scans are required for the whole object reconstruction from the truncated data. For simplification, let the object support Ω be an equilateral triangle, as in the figure 1.(b). The other parameters are defined as above.

As the object takes a great part of the FOV of each scan, no less than short-scan is able to reconstruct the part of the object inside the FOV, seperately, according to Noo's sufficient condition. A good choice is short-scan. The scanning angle of short-scan is π plus the fan angle 2 * γ m . The data redundancy here is larger than in case of twice scans. We propose a triple super-short-scans to decrease

λ0 s = π + κ 1 , λ0 e = 2 * π -κ 2 (4)
λ1 s = π 2 + ∠O2O1a(λ1s), λ1 e = 3π 2 -κ 3 (5) λ2 s = - π 2 + ∠O1O2a(λ2s), λ2 e = π 2 + ζ 3 (6)
where

h = r * sin (∠O0V0a(λ0s)) = 1/2 * r, ζ 1 = ζ 2 = ζ 3 = ∠O2O1a(λ1s) = arcsin(h/R) = arcsin( r 2 * R ), κ 1 = ∠O0V0a(λ0s)+ζ 1 , κ 2 = ∠V2V0O0-ζ 2 , κ 3 = π/3-arcsin( r 2 * R
), and ∠O1O2a(λ2s) = π/3 + arcsin( r 2 * R ). All the three scanning angles |λi e -λi s | of Π i are 2/3π, with i=0,1,2. For the whole object reconstruction, our triple super-short scans only needs 2π fan-beam projections, which are π + 6 * γ m less than those of the triple short-scans, since each shortscan needs π + 2 * γ m fan-beam projections.

We verify the sufficient condition where every lines through the object support should be known in the figure 3: Let M be an arbitrary pixel inside the object. All the lines passing M are intersected with the three colorful solid arcs (black, red, and blue). Therefore, M is able to be reconstructed exactly. This will be validated in the results of image reconstruction from the simulated data.

III. Image reconstruction

In this paper, we consider the image reconstruction as an optimisation problem. A standard iterative leastsquares (ILS) method is used for image reconstruction. The analytic methods suitable for the multiple scans are beyond this topic.

In fan-beam tomography, a forward model of data acquisition is given as follows:

g = Hf + ( 7
)
where g is the mesurement vector, H represents the system matrix, whose element h ij means the contribution of j th pixel of the object on the i th mesurement unit, the vector f describes the unknown object, and contains the detector noise and modelling errors. f is estimated by minimizing the criteria function J(f ) defined below:

J(f ) = 1/2 g -Hf 2 (8) f = arg min f J(f ) (9)
where . 2 is L 2 norm. In the case of the twice scans or triple scans, H is rewritten as

H= H 1 H 2 , or H=   H 0 H 1 H 2  
respectively, where H 0 , H 1 , and H 2 are the system matrices of the scans Q 0 , Q 1 , and Q 2 , separately. In case of the twice, the gradient ∇J(f ) and the optimal step α of ILS are deducted as follows:

∇J(f ) = H T 1 (H 1 f -g 1 ) + H T 2 (H 2 f -g 2 ) (10) α = ∇J(f ) 2 ∇J(f ) T (H T 1 H 1 + H T 2 H 2 )∇J(f ) (11)
where H T i are the backprojections, i = 1, 2, g 1 and g 2 are the mesurements of the scans Q 1 and Q 2 .

IV. Simulation studies

The proposed methods of image reconstruction of the whole object from the truncated data are evaluated by the noise-free and noisy simulated data. The geometry configuration of our CT scanner fits for dental imaging, R and D are 440 and 690 mms, there are 680 units of the detector with size of 120 µm, as a result, the achievable image size by a single circular scan is 52 × 52 mm 2 . However, in our simulation, the support size of the standard modified Shepp-Logan (MSL) phantom (Figure 5) is 72 × 24 mm 2 , whose attenuation coefficients are defined in matlab, and a larger simulated jaw (Figure 6) is included in an equilateral triangle support with side length 90 mm. Twice or triple circular scans are required to reconstruct the whole object. Our jaw contains 3 materials to simulate the tissu, teeth and high attenuated implants with attenuation coefficient of 0.02, 0.06 and 0.3 mm -1 . The fan-beam data was acquired using our ray-driven projector. To verify a stable recontruction with the proposed method, the noisy data was generated with addition of a gaussian noise in the noise-free fan-beam projections.

The images are reconstructed in 384 × 384 with pixel size 200 µm. The reconstructed images from the noise-free fan-beam data of MSL are presented in the left column of Figure 5. Each half of MSL is accurately reconstructed, even more the pixels inside the FOV (red dashed arc in 5.(a) and .(c)) from the reduced scan. The ILS method is able to reconstruct the MSL inside the FOV not only half of MSL from the reduced scan, as OSEM in [START_REF] Clackdoyle | Quantitative reconstruction from truncated projections in classical tomography[END_REF]. Both twice reduced scans and twice super-short scans allow to reconstruct the whole object exactly, in the figure 5.(e) and .(g). Moreover, the scanning angle of the proposed supershort-scan decreases by 16.3%, from 184 • of the reduced scan (π + 2 * γ 1 ) to 154 • (equations ( 1) and ( 2)). To assure the efficacy of our proposed super-short-scan, we reduce a bit of its scanning rangs, which is refered to the SS --scan, [λ1 s -ε λ1 e + ε] (clockwise) and [λ1 s + ε λ1 e -ε], with a small angle ε (one or two angular steps), in the figure 2 .(e)). However the proposed SS-scan requires only the fanbeam projections of 120 • , which is much less than 187 • (π+the fan-beam angle 2 * γ m ) of the short-scan. The scanning rang using SS-scan decreases by 35.8% compared to the short-scan. Moreover, the proposed SS-scan is the minimal scan for the whole jaw reconstruction. If we reduce a few of fan-beam projections (SS --scan) with the scanning ranges [λi s +ε λi e -ε], i = 0, 1, 2 (Figure 3), some artifacts appear in the figure 6.(c) and .(f).

a) Q 1 (b) Q 1 (c) Q 2 (d) Q 2 (e) R-scan (f) R-scan (

V. Conclusion

We propose a simple and efficient super-short-scan algorithm on twice and triple circular trajectories for the whole object image reconstruction from the truncated data. The scanning angle decreases by 16.3% in the case of twice super-short scans compared to the reduced scans, and by 35.8% in the case of the proposed triple supershort scans compared to the short-scans. Our proposed super-short scans algorithms could be applied in dental imaging for a panoramic view of the whole jaw with a small detector, as well as for a fast scan and reduction of dose in other applications of CT tomography. Furthermore, these super-short scans will be easily extended in cone-beam CT tomography.

Fig. 1 :Fig. 2 :

 12 Fig. 1: System geometries of twice scans (Q 1 and Q 2 ) with ellipse support of the object (a). and triple scans (Q 0 , Q 1 and Q 2 ) for the triangular object support (b). O i denotes the iso-center, i = 0, 1, 2, In (b), V 0 , V 1 , and V 2 are the three vertexes, their three iso-centers, O 0 , O 1 and O 2 , are located on the middle of the segments between the triangle gravity G and each vertex, respectively. The phantoms of Shepp-Logan (a) and simlulated jaw (b) inside Ω are given here only as exemple.

Fig. 3 :

 3 Fig. 3: Illustration of the triple super-short-scans. M is an arbitrary pixel in the object. The scanning ranges Π 0 (M ), Π 1 (M ) and Π 2 (M ) are sufficient for the accurate reconstruction of M.

Fig. 4 :

 4 Fig. 4: Profiles of the reconstructed images of the whole MSL from SS --scan (a) noise-free data and (b) noisy data. The profiles are drawn along the yellow solid line in 5.(a).

(

Fig. 5 :

 5 Fig. 5: Reconstructed images from the nois-free fan-beam data (left column)and the noisy fan-beam data of MSL (right column). The display window is [0.1 0.3].

  Fig. 6: Reconstructed images from the noise-free fan-beam data of the simulated jaw of the different acquisitions (a) short-scan, (b) SS-scan and (c) SS --scan. and from the noisy data (d) short-scan, (e) SS-scan and (f) SS --scan. The display window is [0.01 0.07].


Xi Wang
email: xi.wang@inria.frmarc.christie

Marc Christie

Eric Marchand
email: eric.marchand@irisa.fr

Optimized Contrast Enhancements to Improve Robustness of Visual Tracking in a SLAM Relocalisation Context

Robustness of indirect SLAM techniques to light changing conditions remains a central issue in the robotics community. With the change in the illumination of a scene, feature points are either not extracted properly due to low contrasts, or not matched due to large differences in descriptors. In this paper, we propose a multi-layered image representation (MLI) in which each layer holds a contrast enhanced version of the current image in the tracking process in order to improve detection and matching. We show how Mutual Information can be used to compute dynamic contrast enhancements on each layer. We demonstrate how this approach dramatically improves the robustness in dynamic light changing conditions on both synthetic and real environments compared to default ORB-SLAM. This work focalises on the specific case of SLAM relocalisation in which a first pass on a reference video constructs a map, and a second pass with a light changed condition relocalizes the camera in the map.

I. INTRODUCTION

Visual tracking systems such as SLAM and visual odometry are widely used in industrial and consumer devices. Except for direct methods, eg. [START_REF] Engel | LSD-SLAM: Large-scale direct monocular SLAM[END_REF] working on the analysis of changes in pixel gradients, most visual SLAMs rely on corner detection with extractors that extract keypoints (KP) and descriptors that identify and match the extracted KPs over different frames.

Unfortunately, the corner detection process and consequently the matching problem are strongly dependent on the illumination condition at the moment of capturing images. Although the matching process usually relies on gradient information that is more or less independent from intensity, SLAM methods still suffer from illumination changes at different degrees and may yield inaccurate maps and even tracking failures [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF], [START_REF] Seonwook | Illumination change robustness in direct visual slam[END_REF].

In this paper, we propose to improve the robustness of indirect SLAM techniques to light changing conditions by using a multi-layered image representation (MLI). We target the specific case of relocalisation, in which a first pass in a given lighting condition is used to construct a map, and a second pass in a different and dynamic lighting condition is performed to relocalise the camera in a way similar to [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF] or more recently [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF]. The idea of MLI is to dynamically generate k contrast enhancements of an input image into k layers and apply KP detection on each image layer to improve detection and hence KP matching and tracking (see Fig. 1). The challenge therefore holds in how to compute the optimal parameters of each contrast enhancement to maximize keypoint detection and matching. In this paper, we rely on the information theory approach of Mutual Information to compute the optimal contrast enhancements on each frame of a video sequence.

The contributions of this paper are:

• a multi-layered image representation using different contrast enhancements to improve detection and matching of keypoints • an efficient process to compute optimal parameters for these contrast enhancements using an information theoretic approach • a dramatic improvement of ORB-SLAM tracking in light changing conditions.

Results displayed in Fig. 1 shows that our approach outperforms the default ORB-SLAM technique in terms of percentage of tracked frames in a sequence, hence demonstrating a stronger robustness to light changing conditions.

II. RELATED WORK

Robustness to light changing conditions is a central issue that has received increased attention. The issue has often been tackled at the extractor level by searching an optimal contrast threshold in the KP extractor with respect to the current lighting condition. For example, in SuperFast [START_REF] Florentz | Superfast: Model-based adaptive corner detection scalable robotic vision[END_REF] the FAST contrast threshold -a threshold value that triggers a brighter, darker or similar decision on per-pixel comparisonis dynamically computed using a feedback-like optimization method that yields a new threshold value per region in the image. Lowering the threshold however tends to generate a large number of KPs that influence the computational capac-ity of other processes, and the proposed technique requires specific adaptations to be applied to other KP detectors.

Another possibility consists in applying image transformations (eg. contrast enhancers) on captured images before applying KP detectors. Interestingly, it has been demonstrated that KP extractors gain significant performance by using HDR images as input, converted to SDR images through tone-mapping operators [START_REF] Rana | Evaluation of feature detection in hdr based imaging under changes in illumination conditions[END_REF], [START_REF] Přibyl | Feature point detection under extreme lighting conditions[END_REF]. Among these techniques, a learning-based optimal tone-mapping operator has been proposed for SIFT-like detectors [START_REF] Rana | Learning-Based Tone Mapping Operator for Image Matching[END_REF]. But the high computational cost and specific HDR devices required, as well as HDRcustomized extractors hamper the wider applicability of such approaches. In comparison, for SDR images, research has mainly focused on contrast enhancement operators for aesthetic and perceptional goals through changes in the exposure times [START_REF] Yuan | Automatic exposure correction of consumer photographs[END_REF] which remain limited in addressing robustness of KP tracking.

For direct and semi-direct SLAM methods, i.e. methods that rely on analysis of pixel intensities rather than extracting intermediate features, robustness to illumination changes has been addressed by optimizing an affine brightness transfer transformation between consecutive frames [START_REF] Engel | Direct sparse odometry[END_REF], [START_REF] Engel | LSD-SLAM: Large-scale direct monocular SLAM[END_REF] or matching a dimension reduced deep-learning feature [START_REF] Zhang | Dynamic environments localization via dimensions reduction of deep learning features[END_REF]. Using mutual information instead of photometric error as the metric during the optimization process of pose estimation has also demonstrated its benefits [START_REF] Caron | Direct model based visual tracking and pose estimation using mutual information[END_REF], [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF]. While exhibiting a good robustness to illumination changes, these methods remain computationally expensive.

III. MULTI-LAYERED IMAGES

Our approach consists in computing, for every frame, a number of contrast enhancements of the original camera image into different images (called layers) before applying keypoint detection on each layer. The idea is inspired by the significant gain obtained in keypoint tracking on tone mapped images of High Dynamic Range (HDR) images [START_REF] Rana | Evaluation of feature detection in hdr based imaging under changes in illumination conditions[END_REF], [START_REF] Přibyl | Feature point detection under extreme lighting conditions[END_REF], but replacing an HDR image by artificially-exposed multiple images in a Multi-Layered Image representation called MLI. These enhancements are computed in a way to improve the detection and matching of keypoints. The parameters of the first contrast enhancement are searched by maximising the shared information between a reference image (in a good lit condition) and the contrast-enhanced image. We rely on a Mutual Information metric, an information theory measure of dependence between two random variables (images), and demonstrate the relevance of this metric in keypoint tracking (see Section IV). The parameters of the second layer are then searched by maximising the mutual information with the reference image, without the information already provided by the first layer. Other layers are computed in a similar incremental way. Given that landscapes of the Mutual Information metrics are difficult to optimize, a specific smoothing process is proposed that enables the use of straightforward gradient descent optimization techniques.

The contrast enhancement technique relies on a saturated affine brightness transfer per-pixel function (SAT). We use a SAT form that defines a contrast band u = (a, b) > which conveniently models the lower cut point (a) and higher cut point (b) of the saturation, with a linear interpolation between a and b on pixel intensity i (see Fig. 2). A given contrast u = (a, b) > is defined in a contrast space ✓ R 2 , where is the space of all contrast bands where b > a. This MLI representation can be seamlessly integrated in a state of the art SLAM tracking technique such as ORB-SLAM, replacing the keypoint detection of the camera image, by (i) a set of contrast enhancements of the camera image, and (ii) keypoint detection applied on each enhanced image (see Section V).

f SAT (i, (a, b) > ) = min(max(0, i/(b a)), 1) (1)

IV. OPTIMAL IMAGE ENHANCEMENT

The challenge consists in computing the best parameters for each contrast enhancement on each camera image in a way to improve detection and matching of keypoints. Our hypothesis is that we can compute a close to optimal value of u k by maximising Mutual Information between a well lit reference image I ⇤ and a transformed test image

f SAT (I, u k ).

A. Mutual Information

Mutual information (MI) was initially introduced in information theory [START_REF] Shannon | A mathematical theory of communication[END_REF], and then widely applied in the field of computer vision for image alignment, model registration as well as visual tracking and SLAM [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF], [START_REF] Dame | Second-order optimization of mutual information for real-time image registration[END_REF], [START_REF] Thevenaz | Optimization of mutual information for multiresolution image registration[END_REF]. The MI built from the image entropy of two different images provides a measure of their mutual dependence. In image alignment tasks, the higher the mutual information, the better the alignment since mutual information considers the distribution of the intensities as well as the intensities themselves.

Entropy h(I) is a variability measure of a random variable I. In image alignment or illumination evaluation scenarios, I is regarded as one image with r the possible values (graylevel intensities) of I. Equation p I (r)=P (I=r) therefore expresses the probability distribution function of r, in other words the normalized histogram of the image. The Shannon entropy h(I) of an image I is expressed as:

h(I) = X r p I (r)log(p I (r)) (2)
With the same principle, the joint entropy h(I, I ⇤ ) of two images I and I ⇤ can be defined in the following way:

h(I, I ⇤ ) = X t,r p II ⇤ (t, r)log(p II ⇤ (t, r)) (3)
where t and r are the possible grey-level intensities of the I and I ⇤ . The joint probability distribution function is defined as p

II ⇤ (t, r)=P (I=t \ I ⇤
=r), which can also be regarded as a normalized bi-dimensional histogram of images

I and I ⇤ .
With the above notations of entropy and joint entropy, the mutual information (MI) is expressed as the intersection of two random variables I and I ⇤ (see Fig. 3):

MI(I, I ⇤ ) = h(I) + h(I ⇤ ) h(I, I ⇤ ) (4)

B. Optimal enhancement for the first layer

We need to search for the optimal parameter u ⇤ that maximises the MI between a reference image I ⇤ (eg. an image lit in normal lighting conditions) and a contrast enhanced version of camera image f SAT (I, u).

u ⇤ = argmax u MI(f SAT (I, u), I ⇤ ) (5)
We can empirically show that the MI has similar behavior to the ground truth wrt illumination changes. Given a reference image I ⇤ under a given light condition, and a test image I in a different lighting condition, the computation of the ground truth (i.e. the absolute optimal enhancement) can be performed by an exhaustive sampling of the contrast band parameter u k , applying corresponding image transforms on I and evaluating the number of matched keypoints between I ⇤ and f SAT (I, u k ), as displayed in Fig. 7. We illustrate this on an example from the NewTsukuba data set [START_REF] Peris | Towards a simulation driven stereo vision system[END_REF]. We compare the landscapes generated by sampling u k on (1) the ground truth ORB detector and on (2) mutual information Eq. ( 4). Despite differences, we observe the optimums are positioned at similar contrast band values. In an obvious way, the more information is shared between a reference image and a contrast-enhanced image, the better are the detection and matching.

C. Optimal enhancements for other layers

The parameters of the second layer are searched by maximising the Mutual Information with the reference image, as well as the information already provided by the first layer (see Fig. 3).

This can be expressed as multivariate mutual information with the definition of higher dimensional joint probability distribution, to account for multiple image layers. From Eq. ( 2) and (3), a joint entropy of 3 random image variables is obtained with a definition of normalized tri-dimensional histogram p Given I ⇤ an image under reference light conditions, I the test image for which the contrast bands need to be computed, and I 0 = f SAT (I, u 0 ) the first contrast band computed by Eq. ( 9), we can express the tri-variable mutual information to represent the low-correlated information generated in second layer (see Fig. 3 right red line).

MI(I, I ⇤

)

MI(I, I ⇤ , I 0 ) = MI(I ⇤ , I|I 0 ) = h(I, I 0 ) + h(I ⇤ , I 0 ) h(I, I ⇤ , I 0 ) h(I 0 ) (8)
Given the contrast band from first layer u 0 , the optimization of second layer is carried out as follows:

u ⇤ = argmax u MI(I ⇤ , f SAT (I, u)|f SAT (I, u 0 )) (9)
The computation of further layers can be expressed in a similar way. The mutual information between four images

MI(I, I ⇤

, I 0 , I 1 ) or more is computationally expensive to achieve. However a reasonable approximation can be computed using mutual information of previous optimal image MI(I, I ⇤ , I 1 ) to replace MI(I, I ⇤ , I 0 , I 1 ), which balances the computational cost and preciseness.

D. Smoothing Mutual Information

Derivative optimization approaches favor smoother objective function landscape to make sure an efficient descent to the optimum. Unlike image alignment where reducing the number of bins is important to smooth the cost function (image alignment indeed concentrates more on the geometric information instead of illumination information in one image [START_REF] Josien | Mutual information matching and interpolation artifacts[END_REF], [START_REF] Dame | Second-order optimization of mutual information for real-time image registration[END_REF], [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF]). For illumination estimation, lowering the histogram bins during the estimation does not gain any benefits but loses information. This is illustrated on another example from the NewTsukuba data set in Fig. 4 by presenting the landscape of the cost function MI(I ⇤ , f SAT (I, u)) (see Fig. 5). In our work, we selected 256 bins for the purpose of illumination estimation combined with an sigmoid-smoothed SAT function Eq. ( 11). f

SgSAT (I, u) based on Eq. ( 1) supporting a more curved transition at cutting points which leads to a smoother and less aliased cost function landscape (see Fig. 6). Here we show our sigmoid function Sg(x) with k = 1/(b a):

Sg(x) = 1 1 + 8ke ( (x (b+a)/2))) with u = (a, b) > (10)
With the definition of d as the length of contrast band, d = b a, we have the sigmoid-smoothed SAT function f SgSAT (I, u) defined as: ) computed by Eq. ( 8), notations keep as aforecited for I ⇤ and I with I 0 generated by the optimal contrast band from MI(I ⇤ , f SAT (I, u)). We see clearly that standard SAT causes aliasing-like effect in the landscape, due to derivative discontinuity at cutting points a, b

f SgSAT (I, u) = max(1 d, 0) ⇥ Sg(I) +min(d, 1) ⇥ f SAT (I, u) (11)

E. Optimization Framework

Using the concepts introduced, we propose an optimization framework that computes a multi-layered image representation for each frame of a video sequence.

As illustrated in Algo. 1, the first step relies on the cost function of standard mutual information Eq. ( 5) to compute the optimal result of the first layer. The following layers are then computed by optimizing a cost function measuring conditional mutual information with the previous result. This aims at finding the best contrast band (i.e. with the lowest information correlation to the previously computed contrast bands). I ⇤ represents the image under reference light condition and I the current test image to optimize, u i , i = 0..N is referring to the computed optimal contrast band of each layer with a layer number N .

Algorithm 1 Optimal MLI Generated by MI

1: i 0 2: u 0 argmax u (MI(I ⇤ , f SgSAT (I, u))) 3: while i < N do 4: u i argmax u (MI(I ⇤ , f SgSAT (I, u)|f SgSAT (I, u i 1 ))) 5: i i + 1 6: end while 7: return {u k } k=1..N
A demonstration with the NewTsukuba data set (see Fig. 4) in Fig. 7 illustrates the idea of our multiple step optimization framework. First step is the computation of MI between two images, shown in Fig. 7 representing the first layer (layer 1). The second and third layer rely on the computation of the first layer and instead of optimizing a standard MI cost function, a conditional mutual information cost function is optimized using Eq. ( 9). Comparing with the ground truth generated by ORB [START_REF] Rublee | Orb: An efficient alternative to sift or surf[END_REF] detector, our proposed method presents a highly similar behavior as well as a characteristic of derivability for gradient optimization methods. Ground truth for layers 2 and 3 are generated by a subtraction between the common keypoints detected from the previous optimal contrast band u ⇤ and the keypoints from all others contrast bands i.e. the ground truth of layer 2 is computed by removing all keypoints common with keypoints detected in previous optimum u ⇤ = (0, 0.15) > in layer 1. Empirical results also show that even with relatively different reference images (c), the landscapes are similar.

V. EVALUATION AND EXPERIMENTS

To evaluate the benefits of our approach in visual SLAM relocalisation tasks, we first select a synthetic scene benchmark under different static and dynamic lighting conditions [START_REF] Peris | Towards a simulation driven stereo vision system[END_REF]. This dataset encompasses four videos rendered with identical virtual camera trajectories in a synthetic scene with different illumination conditions (Daylight, Fluorescent, Lamps, Flashlight).

We then designed a real scene benchmark in different static and dynamically changing lighting conditions by executing a same camera trajectory using a robotic arm (see companion videos). In both benchmarks, using a reference video in a given lighting condition, we tested the robustness of our approach compared to default ORB-SLAM in localising the camera from the second video sequence against the keyframes generated from the first video sequence, in a way similar to [START_REF] Royer | Monocular vision for mobile robot localization and autonomous navigation[END_REF] or more recently NID-SLAM [START_REF] Pascoe | Nidslam: Robust monocular slam using normalised information distance[END_REF]. In each benchmark, we report the success rate, i.e. the percentage of the frames from second video successfully relocated. Our implementation is integrated in ORB-SLAM [START_REF] Mur-Artal | Orb-slam: A versatile and accurate monocular slam system[END_REF].

Results of NewTsukuba data set are displayed in Table . I. The table reports an improved success rate against illumination changing environments in all but one condition, and provides a 96.5% success rate where both default ORB-SLAM and NID-SLAM fail (0%) in the Lamps to Flashlight condition. The reason related to the failure comparison needs to be investigated in the future work. Optimal contrast bands for the sequences are computed with a relative low sampling frequency wrt to image acquisition frequency (renew a contrast band every 5 to 10 frames). In the case of real scenes, we placed a monocular camera on a trajectory memory 7 DoF Franka robot arm to guarantee that the camera movement in each video is identical. In comparison with the synthetic scene, a strongly dynamic lighting condition is introduced, by having two operators randomly move spots lights in the experiment scene. Using the same success rate criterion, our MLI ORB implementation managed to track 100% of the dynamically lit scene against a keyframe map generated under normal lighting condition, while normal ORB-SLAM only retrieved 52.12% of the keyframes. Fig. 9 shows the inlier keypoints after graph optimization process (details see ORB-SLAM [START_REF] Mur-Artal | Orb-slam: A versatile and accurate monocular slam system[END_REF]), which can be regarded as trusty tracked points generated in the current frame. It demonstrates that MLI performs dramatically better than default ORB which frequently lost tracking during the video. Screenshots of the experiment video are displayed in Fig. 8. A better tracking quality especially around dark or over-exposed area of non-uniformly light images can be observed.

VI. CONCLUSIONS

We have introduced a novel multi-layered image representation based on mutual information optimization to tackle the illumination robustness problem in SLAM relocalization tasks. Each layer in MLI provides low-correlated information which helps to enhance the contrast and therefore increase the robustness during keypoints tracking process under varying illumination conditions. The optimal parameters are computed using a multiple steps mutual information optimization framework. The proposed method shows significant improvements on both synthetic and real videos.

Fig. 1 :

 1 Fig. 1: Keypoint tracking results in different lighting conditions. Only a few keypoints are matched between reference condition (image a) and standard ORB-SLAM (image b), compared to our novel MLI method (image c).

Fig. 2 :

 2 Fig. 2: Using SAT function with different contrast bands to generate a multi-layered image representation (MLI).

  Parameters a and b naturally represent the band region where the contrast is enhanced, which motivated the choice of this operator. To ensure enhancement or compression of contrasts, we define the range of values for u = (a, b) > as a 2 [ 1, 1] and b 2 [0, 1]. The computation of a layer k in our MLI representation is performed by applying the following operator MLI k on all pixel intensities of the image using a contrast band u k . A MLI is therefore represented as a set of k image layers where MLI k (I) = f SAT (I, u k ) for an image I, where f SAT (I, u k ) is the application of f SAT on all pixels of I.

7 )Fig. 3 :

 73 Fig. 3: Mutual information between two images (left), and three images (right) defined as the shared entropy between images.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Images I ⇤ and I from NewTsukuba data set [16]: a synthetic data set with identical camera trajectories and variant illumination conditions

Fig. 6 :

 6 Fig. 6: Comparison between standard SAT function and sigmoid-smoothed SAT function wrt u = (0.3, 0.5) > . Second row shows the conditional mutual information MI(I ⇤, I|I 0 ) computed by Eq. (8), notations keep as aforecited for I ⇤ and I with I 0 generated by the optimal contrast band from MI(I ⇤ , f SAT (I, u)). We see clearly that standard SAT causes aliasing-like effect in the landscape, due to derivative discontinuity at cutting points a, b

Fig. 7 :

 7 Fig. 7: In each layer, similar behaviors are shown with regard to optimums (see image a,b), compared with ORB detector.Ground truth for layers 2 and 3 are generated by a subtraction between the common keypoints detected from the previous optimal contrast band u ⇤ and the keypoints from all others contrast bands i.e. the ground truth of layer 2 is computed by removing all keypoints common with keypoints detected in previous optimum u ⇤ = (0, 0.15) > in layer 1. Empirical results also show that even with relatively different reference images (c), the landscapes are similar.

Fig. 8 :

 8 Fig. 8: Results of real scene against dynamic illumination variance. With a keyframe map generated under reference condition, MLI shows a better retrieve capacity especially when encountering non-uniform illumination variance. In contrast, standard ORB-SLAM only tracks the well lighted parts in the image.

Fig. 9 :

 9 Fig. 9: The number of trusty inlier keypoints after graph optimization of ORB-SLAM which is a critical indicator displaying the tracking quality. MLI generates significantly better results than standard ORB under the dynamic light changing environment.

  NID ORB MLI NID ORB MLI NID ORB MLI Daylight 99.3 100 100 96.7 96.2 100 73.9 97.6 99.7 74.6 79.8 90.7 Fluo 95.0 88.1 99.8 99.7 100 100 85.3 93.9 100 95.8 100 90.5 Lamps 88.3 55.7 99.0 93.6 79.8 94.1 93.1 100 100 84.3 37.9 92.4 Flash 23.8 30.7 92.8 92.2 90.6 94.6 0.00 0.00 96.5 92.0 100 99.3

	V 2	V 1	Daylight	Fluo	Lamps	Flash
			NID ORB MLI

TABLE I :

 I

SLAM keyframe retrieval success rate between our MLI implementation, default ORB-SLAM and NID-SLAM.