Yu Liu

Andrew Thomas Cameron

Marcial Gonzalez

Carl Wassgren
email: wassgren@purdue.edu

Modeling granular

Keywords: Blending, Granular material, Finite element method, Multi-scale model

Introduction

Granular material blending plays an important role in many industries ranging from those that manufacture pharmaceuticals to those producing agrochemicals. The ability to create homogeneous powder blends can be critical to the final product quality. For example, insufficient blending of a pharmaceutical formulation may have serious consequences on product efficacy and safety.

A wide variety of blenders are available in the marketplace [START_REF] Fan | Solids Mixing[END_REF][START_REF] Poux | Powder mixing: Some practical rules applied to agitated systems[END_REF]. However, predictive engineering design of industrial powder blenders remains underdeveloped due to the lack of quantitative modeling tools. As a result, design and scale-up of blending equipment often relies on empirical studies. Discrete element method (DEM) computational models are considered state-of-the-art for predicting blending of particulate materials; however, this modeling approach is not well-suited for modeling industrial-scale blenders due to computational limitations [START_REF] Moakher | Experimentally validated computations of flow, mixing and segregation of non-cohesive grains in 3D tumbling blenders[END_REF]. Moreover, determining particle-level input parameters to use in such models is challenging and not widely agreed upon [START_REF] Asaf | Determination of discrete element model parameters required for soil tillage[END_REF][START_REF] Barrios | Contact parameter estimation for DEM simulation of iron ore pellet handling[END_REF][START_REF] Coetzee | Discrete element parameter calibration and the modelling of dragline bucket filling[END_REF].

In recent work by Liu et al. [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF], a multi-scale model was developed by combining particle diffusion coefficient correlations with advective flow field information from finite element method (FEM) simulations. The model was able to quantitatively predict the magnitude and rate of powder blending in a two-dimensional rotating drum blender and was computationally faster than DEM simulations. The current work extends this multi-scale modeling approach to a more complex, three-dimensional blender geometry with transient powder dynamics.

Background

Most models that track the spatial and temporal evolution of particle blending can be categorized as being either of the discrete or continuum type. Discrete algorithms model the dynamics of each particle individually while continuum algorithms ignore individual particles and instead assume that the particle assembly is a continuum.

The computational discrete element method (DEM) has been used extensively for predicting particle blending dynamics and is considered state-of-the-art [START_REF] Sarkar | Effect of particle size on flow and mixing in a bladed granular mixer[END_REF][START_REF] Yamane | Discrete-element method application to mixing and segregation model in industrial blending system[END_REF][START_REF] Liu | DEM study of the transverse mixing of wet particles in rotating drums[END_REF][START_REF] Chandratilleke | Flow and mixing of cohesive particles in a vertical bladed mixer[END_REF]. Although these studies have provided valuable qualitative insight, they are limited in their quantitative predictions by computational requirements. Even a small, lab-scale blender containing 100 µm particles contains more than 10 9 particles, far exceeding what is possible to model with standard computational tools. Typically, modelers will increase the particle size or decrease the workspace size in order to reduce the number of particles that must be modeled; however, previous work [START_REF] Sarkar | Effect of particle size on flow and mixing in a bladed granular mixer[END_REF][START_REF] Bridgwater | Geometric and dynamic similarity in particle mixing[END_REF] has shown that particle size can have an influence on the rate and extent of blending. Hence, DEM may not produce quantitatively accurate results for systems of industrial interest.

Analytical continuum models have also been proposed for modeling mixing and segregation [START_REF] Bridgwater | Particle mixing and segregation in failure zonestheory and experiment[END_REF][START_REF] Savage | Particle size segregation in inclined chute flow of dry cohesionless granular solids[END_REF][START_REF] Gray | Particle-size segregation and diffusive remixing in shallow granular avalanches[END_REF]. To be tractable analytically, these early models oversimplified the effects of key physical properties and the predictions from these purely theoretical studies were not very accurate. To improve the accuracy of these models, recent works [START_REF] Fan | Modelling size segregation of granular materials: The roles of segregation, advection and diffusion[END_REF][START_REF] Schlick | Granular segregation in circular tumblers: Theoretical model and scaling laws[END_REF] combined DEMcalculated correlations for particle diffusion and segregation at a local scale with analyticallyderived flow fields. These models have been shown to provide quantitatively accurate predictions; however, since they depend on analytical solutions for the macroscopic flow field, their use is restricted to simple two-dimensional geometries, such as a simple heap flow.

Recently, researchers have utilized multi-scale modeling approaches to make quantitative blending predictions. Bertuola et al. [START_REF] Bertuola | Prediction of segregation in funnel and mass flow discharge[END_REF] combined velocity field data predicted from finite element method (FEM) simulations, which treat the granular material as a generalized Newtonian fluid, with local-scale mixing and segregation relations derived from previous works [START_REF] Fan | Modelling size segregation of granular materials: The roles of segregation, advection and diffusion[END_REF][START_REF] Hajra | Granular mixing and segregation in zigzag chute flow[END_REF]. Key model parameters needed to be fit to the experimental data for quantitative agreement. Bai et al. [START_REF] Bai | FEM simulation of particle flow and convective mixing in a cylindrical bladed mixer[END_REF] used an FEM model assuming Mohr-Coulomb constitutive behavior to predict the velocity field in a cylindrical, bladed mixer, which was shown to be qualitatively accurate compared to DEM simulations. They also predicted the degree of blending by assuming convective mixing only (neglecting diffusive mixing), but observed a dependence of the mixing rate on FEM mesh size. Liu et al. [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF] also investigated particle blending using a multi-scale modeling approach. An FEM model, which assumed Mohr-Coulomb material behavior, was used to predict the velocity field in a two-dimensional rotating drum. This information was then combined with particle diffusion correlations. The model was shown to quantitatively predict the magnitude and rate of blending when compared to DEM simulations without the need for backfitted parameters.

The current work extends the work by Liu et al. [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF] to investigate blending in a more industrially-relevant Tote blender. Several key implementation details are also different. First, the model utilizes transient velocity field information from the FEM simulations instead of a steady velocity field. Second, the governing equations are extended to three dimensions instead of two dimensions. Together, these two modifications greatly increase the flexibility of the model. In addition, rotating drum experiments are used to calibrate the diffusion coefficient used in the multi-scale model and the model predictions are compared against published experimental results. Section 3 of this paper introduces the FEM modeling approach and implementation for the current work. Section 4 describes the advection-diffusion equation used in the multi-scale model and the numerical method used to solve it. Section 5 describes the material calibration methods and experiments. And in Section 6, comparisons are made between published experimental results and the multi-scale model predictions.

Finite element method model

A three-dimensional, coupled Eulerian-Lagrangian, FEM model is used in the present work to provide predictions of the advective flow field in a Tote blender. Previous works [START_REF] Zheng | Modelling the granular flow in a rotating drum by the Eulerian finite element method[END_REF][START_REF] Zheng | Finite element investigation of the flow and stress patterns in conical hopper during discharge[END_REF][START_REF] Swaminathan | Modeling the Formation of Debossed Features on a Pharmaceutical Tablet[END_REF][START_REF] Swaminathan | Characterizing the powder punchface adhesive interaction during the unloading phase of powder compaction[END_REF][START_REF] Liu | Modifications to Johanson's roll compaction model for improved relative density predictions[END_REF] have shown that FEM models can accurately simulate granular material behavior well, including advective flow fields [START_REF] Zheng | Modelling the granular flow in a rotating drum by the Eulerian finite element method[END_REF][START_REF] Zheng | Finite element investigation of the flow and stress patterns in conical hopper during discharge[END_REF]. Details of the model implementation can be found in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. The following sub-sections describe the model geometry, boundary conditions, and initial conditions.

Model geometry and boundary conditions

The commercial FEM package Abaqus/Explicit v6.14 is used to perform the bulk flow simulations. The geometry of the simulated Tote blender is based on the experiments by Sudah et al. [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF]. Those experiments were carried out in a 14-L GEI Gallay Tote blender. Details of the geometry and dimensions of the blender are presented in their work [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF]. For convenience, the dimensions are also shown in Figure 1. A Mohr-Coulomb elastoplastic model is used in the current work to describe the stress-strain behavior of the particulate material. This model has been shown previously to accurately describe granular flow fields [START_REF] Zheng | Modelling the granular flow in a rotating drum by the Eulerian finite element method[END_REF][START_REF] Zheng | Finite element investigation of the flow and stress patterns in conical hopper during discharge[END_REF]. The material properties needed in the model are bulk density, Young's modulus, Poisson's ratio, and material internal friction angle. The methodology for obtaining those material properties was described in detail in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF].

Boundary conditions applied in the model include the material-wall friction angle, the rotational speed along the axis of rotation, and the gravitational acceleration (g = 9.8 m/s 2 directed in the negative y direction in Figure 1).

Initial conditions

The coupled Eulerian-Lagrangian approach implemented in Abaqus is adopted in the current model. Details of this approach were described in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. As shown in Figure 2, the Eulerian mesh covers the entire material domain and ensures that no material leaks outside of the mesh. The Eulerian Volume Fraction (EVF) value is used to determine the volume of material within each element. A value of EVF = 0 indicates that no material is present in the element while EVF = 1 indicates that the element is completely filled with material. The initial bed state is generated by filling a fraction of the elements with material. Details of the process can be found in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. For monodisperse mixing, material loadings were achieved by assigning user-defined field variables at each material point to represent the initial material concentration. Details of this field variable approach can be found in the Abaqus documentation [START_REF]Dassault Systèmes, Abaqus analysis user's manual[END_REF]. An example of left-right initial loading is shown in Figure 3. The color represents a field variable value from 0 (blue) to 1 (red). For monodisperse mixing with materials A and B, the field variable represents the material concentration of A. A field variable value of one indicates that the element is completely filled with material A, while a value of zero indicates that the element is completely filled with material B. The simulation process consists of two steps. First, the material is allowed to settle as the gravitational acceleration is slowly increased from zero to its final value. The blender remains stationary during this step. This procedure is used to eliminate the transient oscillations of material as it settles and results in a stable bed compressed under its own weight. Next, the blender is allowed to rotate immediately at the rotation speed and the simulation is considered started.

The multi-scale blending model

The current multi-scale model is extended from a previously published two-dimensional model [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. Details of the model description and development can be found in this previous work. The following sub-sections summarize the main aspects of the model and highlight its extension to transient, three-dimensional problems.

Advection-diffusion equation

The advection-diffusion equation is used in the current work to model the spatiotemporal evolution of the concentration of a particular material species, c. The governing equation is derived as,

( ) ( ) c c c t ∂ = ∇ ⋅ ∇ -∇ ⋅ ∂ D v , (4.1)
where c is the local concentration of a particular species of material, D is the diffusion coefficient tensor for that species, and v is the local advective velocity vector. Using the local mass conservation equation and assuming an incompressible material, i.e.,

0 ∇ ⋅ = v , (4.
2) the governing equation can be written in index notation form as,

ij i i j i v c c c D t x x x -   ∂ ∂ ∂ ∂ =     ∂ ∂ ∂ ∂   . ( 4.3)
The self-diffusion coefficient D is a tensor quantity with components Dij, and with off-diagonal components Dij ( ≠ ) an order of magnitude smaller than the diagonal components Dij ( = ) [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. By neglecting off-diagonal components of the self-diffusion coefficient tensor, the index notation form of Eq. (4.3) in three-dimensional form is

2 2 2 2 2 2 yy xx zz x xx y yy z zz D D D c c c c c c c v D v D v D t x x x y y y z z z ∂   ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂     = - + + - + + - +       ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂       . ( 4.4)
The self-diffusion coefficient D is also proportional to the local shear rate γɺ and the local mean particle diameter d [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. The particle diffusivity is approximately 1.9 times larger along the mean flow direction than it is in the perpendicular direction, according to work by Utter et al. [START_REF] Utter | Self-diffusion in dense granular shear flows[END_REF]. Thus, the shear rate-dependent diffusion coefficient D can be written as,

( ) ( ) ( ) 2 2 1 2 2 2 yy 1 2 2 2 1 2 xx x y z y x z zz z x y D k d k d D k d k d D k d k d γ γ γ γ γ γ γ γ γ = + + = + + = + + ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ , ( 4.5)
where ( )

x y z v x v x γ = ∂ ∂ + ∂ ∂ ɺ ,
( )

y x z v y v y γ = ∂ ∂ + ∂ ∂ ɺ , and
( ) z x y v z v z γ = ∂ ∂ + ∂ ∂ ɺ .
The constant k1 can be found from experiments or small-scale DEM simulations, with k2 = 1.9k1.

4.2

Numerical method A finite difference method using a central explicit scheme is used to solve Eq. (4.4) due to its simplicity and computational efficiency [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. Other numerical methods are also available to solve the system of partial differential equations, such as the finite element method, finite volume method, domain decomposition method, and the matrix mapping method [START_REF] Fan | Modelling size segregation of granular materials: The roles of segregation, advection and diffusion[END_REF][START_REF] Schlick | Granular segregation in circular tumblers: Theoretical model and scaling laws[END_REF].

A second-order Tylor Lax-Wendroff scheme is used to generate the finite difference expression of the governing equation [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF],

1 2 2 2 2 0 0 2 2 0 0 0 0 0 0 0 1 1 2 2 1 2 n n n n n n ijk ijk x x ijk x x x ijk y y ijk y y y ijk n n n n n z z ijk z z z ijk x y x y ijk y z y z ijk z x z x ijk c c c c c c c c c c c ν ν µ δ ν ν µ δ ν ν µ δ ν ν ν ν ν ν +         = -∆ - + -∆ - +                     -∆ - + + ∆ ∆ + ∆ ∆ + ∆ ∆         , ( 4.6)
where, ( )

0 x x xx ijk x v D t x ν   + ∆ ∆   = ∆ , (4.7)
( )

0 y y yy ijk y v D t y ν   + ∆ ∆   = ∆ , (4.8)
( )

0 z z zz ijk z v D t z ν   + ∆ ∆   = ∆ , (4.9)
( )

2 x xx ijk t D x µ ∆ = ∆ , (4.10) ( ) 2 y yy ijk t D y µ ∆ = ∆ , (4.11)
( )

2 z zz ijk t D z µ ∆ = ∆ , (4.12)
, 1, , 1, 0

2 n n i j k i j k n x ijk c c c + - - ∆ = , (4.13) 1, , 1, , 0 2 n n i j k i j k n y ijk c c c + - - ∆ = , (4.14)
, , 1 , , 1 0

2 n n i j k i j k n z ijk c c c + - - ∆ = , (4.15) 2 , 1, , 1, 2 n n n n x ijk i j k ijk i j k c c c c δ + - = - + , (4.16) 2 1, , 1, , 2 n n n n y ijk i j k ijk i j k c c c c δ + - = - + , (4.17) 2 , , 1 , , 1 2 n n n n z ijk i j k ijk i j k c c c c δ + - = - + , ( 4.18) ( ) 1, 1, 1, 1, 1, 1, 1, 1, 0 0 4
n n n n i j k i j k i j k i j k n x y ijk c c c c c + + -+ + - -- - - + ∆ ∆ = , (4.19) ( ) 1, , 1 1, , 1 1, , 1 1, , 1 0 0 4
n n n n i j k i j k i j k i j k n y z ijk c c c c c + + - + + - - - - - + ∆ ∆ = , (4.20)
( )

, 1, 1 , 1, 1 , 1, 1 , 1, 1 0 0 4 n n n n i j k i j k i j k i j k n z x ijk c c c c c + + -+ + - -- - - + ∆ ∆ = , (4.21) ( ) ( ) ( ) , 1, , 1, 0 2 xx xx i j k i j k x xx ijk D D D + - - ∆ = , (4.22) ( ) ( ) ( ) 1,
( ) ( ) ( ) , , 1 , , 1 0 2 zz zz i j k i j k z zz ijk D D D + - - ∆ = . (4.24)
This finite difference formula is illustrated using the computational molecule shown in Figure 4. Details of the Taylor series expansion and stability of the numerical computations can be found in the previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. It should be emphasized that transient velocity fields are included in the solution to the advection-diffusion equation, Eq. (4.4), by extending the work by Liu et al. [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. Therefore, a C++ post-processing script was developed to process the FEM output files (.obd files [START_REF]Dassault Systèmes, Abaqus analysis user's manual[END_REF]) and handle the large amount of data from these simulations.

As mentioned in Section 3, the initial particle concentrations used in the model are determined by the loading conditions, and the user-defined field variable used in the FEM simulation is set accordingly, as in Figure 3. Since transient velocity fields are used in the current model, the boundaries of the material domain change throughout the entire computation. These boundaries are computed within the FEM simulation from the EVF values and, similarly, are determined by the advection-diffusion solver, after reading EVF values from the FEM output file at each time step. Specifically, the advection-diffusion boundary conditions at the free surface are enforced by setting the material concentration of the boundary node equal to the value of the node that is one grid point inward in the direction normal to the surface. This approach is similar to the polynomial fitting approach often used in computational fluid dynamics (CFD) [START_REF] Chung | Computational Fluid Dynamics[END_REF].

Next, a MATLAB program is used to iterate the finite difference form of the advection-diffusion equation given in Eq. (4.6). After generating initial particle concentrations using the extracted field variable values, the material concentration evolution is iterated using the transient velocity field obtained from the FEM simulation for each time step. A threshold is set to ensure the material concentration value remains between zero and one, and a small time step is carefully chosen to ensure the stability of the explicit scheme. Details of the iteration algorithm were presented in the work by Liu et al. [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF].

To achieve a converged result, the number of elements used in the advection-diffusion MATLAB program should be much larger than the number of elements in the FEM simulation. Hence, a linear interpolation algorithm is introduced in the MATLAB program to generate enough elements and ensure convergence. The MATLAB numerical algorithms are parallelized to divide the computational domain onto different cores. By using the full processing power of an eight-core desktop, the iteration runs approximately 5-8 times faster than on a single core.

Calibration of material properties

Material properties for the FEM simulation

The published experimental work for monodisperse particle mixing by Sudah et al. [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF] was used to validate the multi-scale model predictions. Those experiments utilized a 14-L GEI Gallay Tote blender containing 12-mm diameter glass beads of two different colors.

The material density ρ, elastic modulus E, and Poisson's ratio ν are known to have little influence on the material flow behavior [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. Moreover, the mixing process was shown to be unaffected by the wall friction angle φw within a rotating blender since the free surface angle remains constant. Hence, the material density, elastic modulus, Poisson's ratio, and wall friction angle used in the current FEM simulation are all based on values from previous work for hard spheres [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. In addition, the Mohr-Coulomb dilation angle of the material was set to 0.1°, since the dilation of cohesionless granular materials is usually small. The internal friction angle has been shown to influence the mixing rate [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF] and, hence, must be calibrated. Previous work has shown that the internal friction angle is not sensitive to particle size or consolidation stress for the same material [START_REF] Liu | An investigation of the effect of particle size on the flow behavior of pulverized coal[END_REF][START_REF] Kara | Contribution of particles size ranges to sand friction[END_REF][START_REF] Podczeck | The influence of particle size and shape on the angle of internal friction and the flow factor of unlubricated and lubricated powders[END_REF]. Hence, 1-mm diameter glass beads, as opposed to the 12 mm beads used in the Sudah et al. experiments [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF], were used in the current work to calibrate the internal friction angle. A Schulze Ring Shear Tester (Model RST-XS) was used to make these measurements, obtaining values similar to those reported in the literature [START_REF] Castellanos | Cohesion and internal friction of fine glass beads as affected by small intensity vertical vibration[END_REF][START_REF] Mitarai | Wet granular materials[END_REF]. A summary of all of the material parameters used in the FEM simulation is given in Table 1. It is worth noting that these parameters are found from independent, standard material tests, rather than being back-fit to match experimental data with the multi-scale blending simulation results.

Experimental calibration of the diffusion constant

One significant parameter needed in the multi-scale model is the spanwise diffusion constant k 1 . To calibrate k 1 for glass beads, a lab-scale rotating drum experiment was performed. A photograph of the experiment setup is shown in Figure 5. An acrylic circular drum of diameter D = 150 mm and width W = 50 mm was used to contain the material, and two shafts were used to stabilize and rotate the drum. The driving torque was provided by a gear motor that rotated one of the shafts. Rubber bands covered the shafts to prevent slipping. Values for the experiment parameters are listed in Table 2. As shown in Eq. (4.5), the influence of particle diameter on the diffusion coefficient has been included explicitly and thus diffusion constants k 1 and k 2 should be independent of particle diameters. Hence, the drum was filled side by side with 1 mm red and blue glass spheres. To facilitate filling of the drum, the front side of the drum was made removable and a separate barrier was used to help fill each side of the drum with equal volumes of red and blue glass beads. Friction tape was used to seal the drum and to allow the container to rotate smoothly. To analyze the degree of mixing, a high-speed camera was used to film the front of the drum. Once the drum was filled, the drive shift started to rotate and the mixing process was recorded. Several snapshots of the system at different times are shown in Figure 6. Qualitatively, the mixing dynamics in the experiment followed the same trend as the simulations reported in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. The experimental results were compared with predictions from the 2-D multi-scale model developed in the previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. To justify a comparison, the system geometry and operating conditions, such as the drum diameter, drum width, particle diameter, filling level, and rotation speed, were consistent between the experiment and the multi-scale model. The material parameters used in the FEM simulation are shown in Table 1 for glass beads.

The segregation intensity I was computed to compare the experiment results with the 2-D multiscale model predictions and calibrate the spanwise diffusion constant k 1 . The segregation intensity I is defined as,

2 i 2 0 i I σ σ = , ( 5.1)
( )

2 2 i , 1 1 1 M i i m m c c M σ = = - -∑ , (5.2) ( ) 2 0 1 i i c c σ = -, (5.3) , 1 1 M i i m m c c M = = ∑ .
(5.4) In these relations, σi 2 is the measured variance of component i's concentration (here component i is the red particles), σ0 2 is the variance of component i's concentration for a fully segregated system, c i is the measured mean concentration of component i, and M is the total number of samples used to calculate the mean and variance. The scale of scrutiny does affect the calculation, with larger scales of scrutiny having smaller segregation intensities. Hence, the same grid size used to analyze the experiment results, as described in Appendix A, was used in the 2-D multi-scale model. Computationally, the concentration values of all the nodes within one grid cell were averaged to compute the material concentration for the cell. A cell size of five particle diameters was used.

Due to the fact that a continuum is assumed in the FEM simulations, the multi-scale model would predict an asymptotic segregation intensity of nearly zero, corresponding to a perfectly mixed state. However, a perfectly mixed state is generally not achievably in practice and instead a randomly mixed state is the expected asymptotic state. For a randomly mixed system, the segregation intensity is derived as, where N is the number of particles in the cell used to calculate the concentration. The segregation intensity for a randomly mixed system in the current work is I = 0.04, which is shown in Figure 7 as the dashed line. To save computation time, the simulation was stopped once the segregation intensity reached the randomly mixed state. Segregation intensity plotted as a function of the number of drum revolutions is plotted in Figure 7 for different assumed values of the spanwise diffusion constant k 1 . Figure 8 shows the sum of the absolute differences between the segregation intensities measured from experiments and computed from the 2-D multi-scale model for different k 1 values. It can be seen that the segregation intensities computed from the 2-D multi-scale model match best with the values measured from experiments when k 1 = 0.01. Hence, a calibrated spanwise diffusion constant k 1 for glass beads equal to 0.01 is adopted. Note that the asymptotic value of the experimental results is slightly larger than the value for the randomly mixed state. This difference is because the image correction algorithm described in Appendix A introduces some error into the system and the perfectly random mixed state cannot be reached.

Comparison of the experiments and multi-scale model results

An FEM simulation was performed to predict the advective flow field information for the Sudah et al. Tote blender [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF] using a rotation speed of 10 rpm and the material properties shown in Table 1. Furthermore, since transient velocity field information was used in the current work, the velocity components in all three directions needed to be outputted at every time step, which generates considerable amounts of data. Although the velocity field is not constant within a Tote blender single revolution, the material reaches a periodic steady state condition after a few initial revolutions. Table 3 shows the averaged velocity differences between subsequent revolutions, defined as,

ave

j j N i n i n j n i n i i v v v v v N + = - ∆ = ∑ , (6.1)
where i is the direction of the velocity component, j is the revolution number, and N is the total number of nodes within the material domain. It is evident that the periodic steady state condition is reached after the first revolution. Therefore, to save computational effort without losing accuracy, the FEM simulation was only performed for two revolutions and the velocity fields within the second revolution were used in the multi-scale model subsequently. In the experiments reported by Sudah et al. [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF], materials were initially loaded into the Tote blender in two different ways: top-bottom loading and left-right loading. Two filling levels were also studied -40% fill and 60% fill. The same loading conditions and filling levels were modeled in the current work, as shown in Figure 9. Figure 10 shows the evolution of the material domain within the first revolution for left-right loading and 40% fill. In the experiments by Sudah et al. [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF], mixing curves were constructed by plotting the relative standard deviation (RSD) as a function of mixing time. The relative standard deviation is calculated as,

i RSD c σ = , (6.2) ( ) 2 2
, 1 [START_REF] Fan | Solids Mixing[END_REF] 1

M i i m m c c M σ = = - -∑ . ( 6.3)
In the above equations, σ 2 is the measured variance of component i's concentration, c i is the mean concentration of component i, and M is the total number of samples used to calculate the mean and variance. The relative standard deviation varies from zero, corresponding to perfect mixing, to one, which is a fully segregated state. The RSD value was also computed using the multi-scale model to compare with the published experimental results. Note that in the multiscale model, every node at which a concentration is calculated is used in the evaluation of the segregation intensity, while in the experiments nine core samples were used [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF]. Although the scale of scrutiny plays a role in the calculated results, as shown in previous work [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF], the predicted values using different grid sizes are close to each other.

Figure 12 plots the RSD values with respect to the number of revolutions for top-bottom and leftright initial loading patterns for the multiscale-model and experiments. Results for 40% and 60% fill levels were compared. In the current work, one million cells were used to maintain accuracy and computational efficiency (i.e., 100 cells in each direction). A mesh dependency study was performed to ensure solution convergence and the error was found to be within 5% between one million cells (i.e., 100 cells in each direction) and 3.375 million cells (i.e., 150 cells in each direction). As shown in Figure 12, there is good quantitative agreement between the multi-scale model predictions and experimental measurements, although it does appear that there is some offset for the case with 40% fill and top-bottom loading. There is no error information given for the experimental results so it is difficult to determine how significant this difference is. The total wall-clock time required to run a single case on 10 cores, including the FEM simulation and MATLAB processing, was approximately 2.5 days.

Conclusions

In this work, a three-dimensional transient multi-scale modeling approach is described for predicting blending in particulate systems. This model is extended from a previously published two-dimensional steady model [START_REF] Liu | Modeling granular material blending in a rotating drum using a finite element method and advection-diffusion equation multi-scale model Yu Liu[END_REF]. This three-dimensional model combines, within the advection-diffusion equation, finite element method generated transient macroscopic velocity fields with experimentally-obtained particle diffusion correlations at a local scale. The model is applied here to a three-dimensional Tote blender. Predictions of the mixing rate from the multiscale model compare well quantitatively to published experimental data [START_REF] Sudah | Simulation and experiments of mixing and segregation in a tote blender[END_REF].

A significant advantage of the multi-scale blending model over a DEM-only model is that the multi-scale model is more computationally efficient for industrially-relevant system sizes. This is because the number of DEM particles increases with the cube of the ratio of the system size to particle size, while the FEM nodes do not necessarily need to increase since the mesh size increases linearly with the system size. Furthermore, if particle size is reduced in the DEM simulations, then the integration time step also decreases, further increasing the time required to complete DEM simulations. In addition, all of the parameters used in the multi-scale model can be measured from independent, standard tests or calibrated from simple two-dimensional experiments (such as the material internal friction angle, the wall friction angle, and the diffusion constant). No back-fitting of the Tote blender results was used in the current work.

Future work should focus on expanding the multi-scale model to consider segregation of different materials, which is important in industrial practice. In addition, a simple, standard method should be developed to calibrate the diffusion constants. The current work used a twodimensional blender experiment, but other systems, such as a 2-D heap flow, may be considered.

Appendix A Image Processing algorithm

A flowchart of the image processing algorithm used in the current work is shown in Figure A1.

For each image, the material domains (red and blue glass beads) were extracted using the freeware program ImageJ [START_REF] Schneider | NIH Image to ImageJ : 25 years of image analysis HISTORICAL commentary NIH Image to ImageJ : 25 years of image analysis[END_REF], which was used to identify only the glass beads among a white background. This RGB image was further analyzed to extract only the red component of each image's pixels, with a value ranging between 0 (no red) and 255 (all red). A threshold value of 80 was selected to differentiate between the red and blue pixels. A value larger than 80 indicated that a pixel corresponded to a red particle while a value smaller than 80 indicated that it was a blue particle. This threshold value was chosen to ensure that at any point in time the fraction of red pixels in the entire system was 50±5% since the system consisted of 50% red beads.

Figure A1. The image analysis algorithm for a single image. The steps proceed from 1 to 5.

After distinguishing red and blue pixels, a black and white binary image was generated, with black corresponding to the blue pixels and white corresponding to the red ones. Note that since spherical glass beads were used in the experiment, reflections and shadows were introduced due to the light source. An image correction algorithm was developed in the MATLAB program in an attempt to account for these effects.

The corrected binary image was used to compute the segregation intensity. A non-overlapping grid of square cells was used in the current analysis. The cell size varied from 3 to 10 particle diameters on a side. To avoid using samples containing data outside the cylinder boundary, boundary flags were used to make sure each sample contained at least 95% material. All of the qualified samples were then used to calculate the segregation intensity for the current time step.

Figure 1 .

 1 Figure 1. Dimensions of the GEA Gallay Tote blender used in the FEM simulation.

Figure 2 .

 2 Figure 2. A schematic of the FEM model domain. The Eulerian mesh is shown in grey, the Tote blender mesh is shown in blue, and the outlines of the initial material domain are shown in red.

Figure 3 .

 3 Figure 3. Initial material concentration for a simulation with left-right loading. The color represents field variable value from zero (blue) to one (red). The materials in this simulation have homogeneous properties, except for color.

Figure 4 .

 4 Figure 4. Computational molecule for the second-order, three-dimensional Taylor Lax-Wendroff scheme.

Figure 5 .

 5 Figure 5. The rotating drum experiment setup.

Figure 6 .

 6 Figure 6. Snapshots showing the state of particle mixing for different numbers of drum revolutions in the rotating drum experiments.

Figure 7 .

 7 Figure 7. Segregation intensity with respect to the number of drum revolution for the experiment and the 2-D multi-scale model using different k 1 values.

Figure 8 .

 8 Figure 8. Sum of the absolute differences between the segregation intensities measured from experiments and computed from the 2-D multi-scale model for different k 1 values.

Figure 9 .

 9 Figure 9. The left-right (a) and top-bottom (b) initial loadings for 40% fill.

Figure 10 .

 10 Figure 10. Snapshots showing the change of the material domain in the FEM simulation (a) and multi-scale model (b) for left-right loading and 40% fill. The vertical color scale in (b) is the red particle concentration.

Figure 11 .

 11 Figure 11. Snapshots showing the state of mixing at the end of each revolution for (a) left-right loading and (b) top-bottom loading and 40% fill. The vertical color scale is the red particle concentration for (a) and (b).

Figure 12 .

 12 Figure 12. Mixing state comparison between published experimental results and the multi-scale model predictions. The random mixing state is represented by horizontal dashed lines and the spanwise diffusion constant is k 1 = 0.01.

Table 1 .

 1 Parameters used in the FEM simulation.

	Parameter	Value
	Material density (kg/m 3 )	1500
	Young's modulus (MPa)	3.65
	Poisson's ratio (-)	0.065
	Internal friction angle (degree)	27.1
	Cohesion (Pa)	87
	Dilation angle (degree)	0.1
	Wall friction coefficient (-)	0.324

Table 2 .

 2 The rotating drum experiment parameters.

	Parameter	Value
	Inner drum diameter (mm)	150
	Inner drum width (mm)	50
	Glass sphere diameter (mm)	1
	Filling level (% of max level depth) 32
	Drum rotation speed (rpm)	3.26

Table 3 .

 3 The averaged velocity differences between subsequent revolutions.

	Revolution numbers (j, j+1)	(1,2)	(2,3)	(3,4)
	∆	/ v v x	x	ave	3.95%	0.56%	1.85%
	∆	/ v v y		y	ave	2.30%	0.38%	3.03%
	∆	/ v v z	z	ave	2.91%	1.18%	2.53%


Jean-Baptiste Bellet
email: jean-baptiste.bellet@univ-lorraine.fr

Gérard

Gérard Berginc
email: gerard.berginc@fr.thalesgroup.com

GEOMETRIC PROBLEMS FROM 3D REFLECTIVE TOMOGRAPHY

Keywords: Mathematics Subject Classication. 78A97, 94A12, 44A12 3D imaging, computational optics, reconstruction, Radon transform, geometric tomography

1. Introduction 1.1. Three-dimensional optical imaging. There is a considerable interest in the development of new optical imaging systems that are able to give three-dimensional images. Potential applications range across the eld of defense and security for the recognition of targets, the medical eld for the detection of subcutaneous and cutaneous tumors or the archaeological eld for the discovery of remains in forests. The framework is the following: given a set of bi-dimensional images in the visible to near-infrared band, construct a threedimensional model of the original scene. This eld is related to several topics. 1.1.1. Radiometry. The rst one deals with image formation [START_REF] Horn | Robot vision[END_REF][START_REF] Mohammad-Djafari | Problèmes inverses en imagerie et en vision[END_REF]. This subject, called radiometry, is about the modeling of an optical image. In this eld, a camera is often modeled as an ideal pinhole, which realizes an ideal perspective projection of the scene [START_REF] Ma | An invitation to 3-D vision: from images to geometric models[END_REF]; this provides a geometric model of an image by the means of rays of projection. In radiometry the irradiance represents the amount of light incident on a surface: it is the power per unit area (W.m -2 ) incident on the surface. The radiance represents the amount of light radiated from a surface: it is the power per unit area per unit solid angle (W.m -2 .sr -1 ). The sensors of a camera measure the (incident) irradiance. Since light is an electromagnetic wave, the Maxwell equations govern its propagation and its behavior at interfaces between media. For a scene of several meters and visible to near-infrared wavelengths (0.4µm-3µm), light is often scattered from surfaces, due to subsurface volumic diusion and due to the rugosity of the surface. The models of radiometry take often this phenomenon into account, by the means of the Bidirectional Reectance Distribution Function (BRDF): the BRDF is a surfacic function dened as the ratio of radiance to irradiance; the incident direction and the radiation direction are arguments of the BRDF. It tells how the surface reects light toward a radiation direction when it is illuminated from an incident direction. In particular the BRDF can encode the specular reections on a perfect mirror (Snell-Descartes law); while an ideal matte surface, called a Lambertian surface, is encoded by a BRDF which does not depend on the radiation angle. In order to take into account the eects due to the polarization of light, some models go further by using a vectorized version: the interaction of light with a surface is expressed by the means of a Mueller matrix, or a polarimetric BRDF, which acts on a Stokes vector [START_REF] Chipman | Polarimetry, chapter 22 in: Handbook of optics[END_REF][START_REF] Flynn | Polarized surface scattering expressed in terms of a bidirectional reectance distribution function matrix[END_REF][START_REF] Shell | Polarimetric remote sensing in the visible to near infrared[END_REF]. For a polarized light, incident on a rough surface, some other scattering models are based on the Small Slope Approximation [START_REF] Berginc | Scattering models for range proling and 2D-3D laser imagery[END_REF][START_REF] Voronovich | Wave Scattering from Rough Surfaces[END_REF].

1.1.2. Reconstruction. The second main topic deals with inverse problems: reconstruct the 3D scene by combining the 2D projections. In the eld of computer vision, this problem is the so-called multiple-view reconstruction; standard strategies establish a correspondence of features, such as corners or edges, between the frames; the reconstruction algorithm takes benet of these correspondence, thanks to linear algebra considerations [START_REF] Ma | An invitation to 3-D vision: from images to geometric models[END_REF]. Other methods are based on a voxelized 3D reconstruction. In the shape-from-silhouette problem, the input images are binary; volumic intersection of the corresponding viewing cones, which is a way of backprojecting the images, reconstructs the visual hull of the scene [START_REF] Baumgart | Geometric modeling for computer vision[END_REF][START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF]. Some other methods carve iteratively a 3D volume by checking consistency of the input irradiances [START_REF] Kutulakos | A theory of shape by space carving[END_REF]. In this paper we will discuss reective tomography, which enters also in the class of methods which computes a voxelized reconstruction.

And we mention a method which reconstructs a surface and its irradiance by an iterative tting algorithm [START_REF] Yezzi | Stereoscopic segmentation[END_REF].

1.1.3. Tomography. In transmission tomography, the records are projections containing the X-ray transform of an attenuation; they are processed to provide a reconstruction. The reconstruction method looks like a Radon inversion formula, which is often implemented under the form of a ltered backprojection algorithm [START_REF] Natterer | Mathematical methods in image reconstruction[END_REF]. This is the basis of tomography, while many variants exist. Exact or heuristic inversion formula dedicated to special geometry of acquisition are derived. For instance the famous FDK algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] and its renements such as [START_REF] Zheng | A New Heuristic Weighting Function for FDK-based reconstruction of Cone Beam Tomography[END_REF] are heuristic methods dedicated to cone beam (perspective) projections. Also some other imaging modalities are modeled by related transforms; for instance, in emission tomography, the considered transform is an X-ray transform with attenuation [START_REF] Natterer | Mathematical methods in image reconstruction[END_REF].

The so-called reective tomography draws a parallel with transmission tomography: the records are projections with similar geometries. Then the idea is to compute a reconstruction by the means of a ltered backprojection from transmission tomography [START_REF] Knight | Tomographic Techniques Applied to Laser Radar Reective Measurements[END_REF]. This heuristics is successful for several kinds of reective data [START_REF] Knight | Tomographic Techniques Applied to Laser Radar Reective Measurements[END_REF]; the principle has also been proposed in object modeling from photographies [START_REF] Gering | Object modeling using tomography and photography[END_REF]. A recent patented technology is based on 3D reective tomography [START_REF] Berechet | Method for 3D reconstruction of an object in a scene[END_REF][START_REF] Berechet | Procédé de discrimination et d'identication par imagerie 3D d'objets d'une scène[END_REF][START_REF] Berginc | Optronic system and method dedicated to identication for formulating three-dimensional images[END_REF]; it can reconstruct surfaces from 2D optical images [911], and it recovers partially occluded objects [START_REF] Berginc | Optical 3D imaging and visualization of concealed objects[END_REF].

1.1.4. Example. The following example illustrates the eciency of reective tomography. We consider a sequence of 360 images of size 342 × 181, measured by turning around the scene, one degree step (courtesy of Thales Optronique SA). These images have been recorded by a real laser system. See Figure 1 for samples of the sequence. Using a home-made software, the FDK reconstruction from these 360 images takes 2.6 seconds on a Graphics Processing Unit (Nvidia Tesla C2075). Then the reconstruction is displayed in real-time, using a volume rendering method (Maximum Intensity Projection [START_REF] Wallis | Three-Dimensional Display in Nuclear Medicine and Radiology[END_REF]). In the Figure 2, we represent snapshots of such re-projections. The representations of the reconstruction contains clearly features and details that are useful for the recognition. 1.2.1. X-ray inversion in the visible to near-infrared band. The usual imaging methods, such as multipleview reconstruction, shape-from-silhouette, or transmission tomography, tries to invert some forward model, to nd respectively the location of features, the visual hull, or the attenuation. This is not the case for reective tomography which is an empirical approach. This last heuristics inherits characteristics of the previous methods, without being an exact inversion procedure: it applies the inversion formula of transmission tomography, it contains a backprojection step as shape-from-silhouette, and it recovers features of the original This has a strong consequence. In general, a point of the scene is not visible on each plane of projection, due to occlusions: thus the data are incomplete [START_REF] Rigaud | Reective Imaging Solved by the Radon Transform[END_REF]. Not only applies reective tomography an X-ray inversion on other things than X-ray projections, but also the data are incomplete. By the way, it is known in the eld of transmission tomography that serious diculties appear for incomplete data, such as instabilities and formation of artifacts [START_REF] Natterer | Mathematical methods in image reconstruction[END_REF]. In particular, despite transmission tomography has been extensively studied over the past decades, there are still ongoing works about the Radon inversion of incomplete data: it can characterize the artifacts, and then strategies of artifact reduction can be derived [START_REF] Borg | Full Characterization of Reconstruction Artifacts from Arbitrary Incomplete X-ray CT Data[END_REF][START_REF] Frikel | Characterization and reduction of artifacts in limited angle tomography[END_REF]. For that point, the problem of reective tomography meets the trendy problem of limited Radon data. A study of reective tomography in a well-suited framework should also be a way to describe its artifacts.

1.3. Contributions of this work.

1.3.1. From Optics to Mathematics. In this work, we transfer results from optical engineering to applied mathematics. We dene geometric problems and a heuristic solver inspired from reective tomography. We formulate in a proper way a geometrical version of problems that arise in this eld.

We dene an original class of geometric problems. For that purpose we dene a new transform that we call the reective projection; it encodes a geometric version of the incident irradiance on sensors in optics.

The main geometric problem is to invert this transform: reconstruct the initial surfaces of a scene, from the (partial) knowledge of a reective projection of the surfaces. In this way, we have transferred a problem from optical imaging to a problem of geometric tomography: Geometric tomography is the area of mathematics dealing with the retrieval of information about a geometric object from data about its sections, or projections, or both, Gardner [START_REF] Gardner | Geometric tomography[END_REF].

Furthermore, concerning the inversion of the reective projection, we transfer the heuristic principle of reective tomography: apply an algorithm from transmission tomography, for the right geometry of acquisition. In this way, we propose a scheme, in general linear (ltered backprojection), which solves eciently an original non-linear geometric problem. Then the main question is to understand the result of this heuristics.

1.3.2. Discussion about the heuristics. We combine analytic calculations for a canonical case and numerical results in order to discuss the eect of the heuristics. Roughly speaking, this discussion reveals that the heuristics computes a volumic reconstruction, in which some features of the original surfaces are enlightened; these features come from the backprojection of the discontinuities of the data set.

The canonical case concerns the reconstruction of a sphere ∂K from its silhouettes, for a horizontal cone beam scan. We determine the exact structure of such a reective projection: we nd a surface S of singularity, based on the rays of projection that are tangent to the sphere. Concerning the reconstruction, we guess, from this surface S, three regimes representing respectively the noise, the useful reconstruction, and the artifacts:

• the reconstructed value should be in general low;

• the reconstructed value should be high at the points of the initial surface ∂K that are projected on a ray of S; • the reconstructed value can be high along some lines of projection which are tangent to the surface ∂K, with a condition (that we determine) on the tangent plane.

The numerical results consider three kinds of data:

• projections with discontinuities due to the geometry or due to the projected texture;

• projections which are smoothed and thus without discontinuities;

• projections of a uctuating texture which keep the discontinuities of the silhouettes.

The main conclusion of these tests is again that the heuristics looks essentially like a backprojection of the discontinuities of the projections, at their true location in space.

1.4. Impact of this work.

1.4.1. In Mathematics. This work does not aim at lling denitively the mathematical gap described in the previous section. But it introduces several ingredients for further studies. Geometric problems of interest have been dened in a proper way, and we have extended the observations of the discussion under the form of conjectures. We suggest approaches to explore this new kind of geometric problems. Proofs based on asymptotic expansions are currently under study [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF], while microlocal analysis [START_REF] Borg | Full Characterization of Reconstruction Artifacts from Arbitrary Incomplete X-ray CT Data[END_REF][START_REF] Ramm | Reconstructing singularities of a function from its Radon transform[END_REF] should be also a relevant framework: we have seen that the singularities play a central role. This work lays the foundation for more advanced studies of the subject.

By the way, geometric tomography establishes a correspondence between transmission tomography and some geometric problems [START_REF] Gardner | Geometric tomography[END_REF]. In this work, this correspondence goes further: the proposed heuristics directly uses algorithms from the rst class of problems to solve problems of the second class: algorithms of transmission tomography solve eciently geometric problems! More generally we could imagine extensions of the principle to solve eciently other geometric problems.

1.4.2. In Optics. The proposed geometric problems are the core of the solver in reective tomography. Hence understanding them is important for optical applications: it is a way to strengthen the subject, and it may yield improvements of the solver, such as denoising, artifact reduction, or acceleration.

Here we have seen that the structure of the projections, especially the singularities, is more important than the values themselves. This has a strong impact: reective tomography is very robust. It enables active/passive imagery, reconstruction of active surfaces, and so on.

1.5. Organization. The paper is organized as follows. In the rst part, we transfer the ideas of reective tomography from optics to geometry; in particular we propose a heuristics for problems of geometric tomography. In the second part, we examine in details the heuristics on the silhouettes of a ball. In the third part, we present several numerical tests which emphasize the impact of the discontinuities. To nish with, we formulate conjectures that must be investigated in further studies.

Mathematical framework

2.1. Geometric problems.

2.1.1. Geometric inversion. We represent a scene of one or several objects by a compact set K ⊂ R 3 . We will project such a scene along rays. An opacity constraint is taken into account using a notion of visibility: Denition 1. A ray is a half-line L = {x 0 + tu, t 0} with x 0 ∈ R 3 and u ∈ S 2 ; L and (x 0 , u) ∈ R 3 × S 2 are identied.

Denition 2 (Visible point). Let K ⊂ R 3 be a compact set, and let L(x 0 , u) ∈ R 3 × S 2 be a ray which

intersects K: L ∩ K = ∅. The visible point v K (L) of K along the ray L is dened by: v K (L) = arg min{(v - x 0 ) • u, v ∈ L ∩ K}; in particular, v K (L) ∈ ∂K ∩ L. We adopt the convention that v K (L) = ∞ if K ∩ L = ∅.
Roughly speaking, for every ray L(x 0 , u), the visible point of K is the rst point of K when we travel from x 0 in the direction u: see Figure 3.

L v K (L) K K x 0 u x 0 u x 0 + uu T (z -x 0 ) K r z v K (L) Figure 3. Left: visible point v K (L) ∈ ∂K ∩ L(x 0 , u). Right: v K (L) when K is the ball of radius r centered at z.
Denition 3 (Reective projection). Let K ⊂ R 3 be a compact set. Let f : ∂K × S 2 → R be a function; we adopt the convention that f (∞, •) = 0. The reective projection of f is

P[f ] : R 3 × S 2 -→ R L(x 0 , u) -→ f (v K (L), u).
Remark 1. The values P[f ](L), or more generally the restrictions of P[f ] will be also called reective projections.

Roughly speaking, along the ray L(x 0 , u), the reective projection P[f ](L) = f (v K (L), u) represents an information which travels from the visible point v K (L) toward the origin x 0 , in the direction -u. This information depends on the visible point v K (L) but also on the direction u. If L does not meet K, P[f ](L) = 0 due to the convention. We propose now a name dedicated to graphical representations of reective projections: We can now formulate a problem of reconstruction from a set of reective projections: Problem 1 (Geometric inversion of reective projections). Let K ⊂ R 3 be a compact set, f : ∂K × S 2 → R, and let L ⊂ R 3 × S 2 be a set of rays. For every L ∈ L, we record the reective projection P[f ](L). We assume that K and f are unknown, while L and P[f ]| L are known. Question: reconstruct the set ∂K from P

[f ]| L .
This problem is in general dicult. The transform of this problem:

P : [f : ∂K × S 2 → R] → P[f ]| L ,
is in general neither invertible, nor linear. Nonetheless we will see a heuristics which can eciently solve such a problem. Remark 2. Similar concepts can be analogously introduced for the two-dimensional case in R 2 . We use the same terminology and the same notations: Problem (Geometric inversion of reective projections in 2D). Let K ⊂ R 2 be a compact set, f : ∂K × S 1 → R, and let L ⊂ R 2 × S 1 be a set of rays. For every L ∈ L, we record the reective projection P[f ](L). We assume that K and f are unknown, while L and P[f ]| L are known. Question: reconstruct the set

∂K from P[f ]| L .
2.1.2. Scanning geometries. We now formulate two special choices for the set of rays L in Problem 1. The rst choice denes the cone beam scan. This case occurs when we consider a camera that turns around a scene on a circle: assuming a pinhole model for the camera, each image is a perspective projection, while the optical center scans a circle. Denition 5 (Cone beam scan). Let r, a, b > 0. For every β ∈ [0, 2π], we set: θ(β) = (cos β, sin β, 0), θ ⊥ (β) = (sin β, -cos β, 0), e 3 = (0, 0, 1), x 0 (β) = rθ(β). For every (y

⊥ , y 3 ) ∈ [-a, a] × [-b, b], we set y = y ⊥ θ ⊥ (β) + y 3 e 3 ; the unit vector that points from x 0 to y ∈ θ(β) ⊥ is u(β, y ⊥ , y 3 ) = y-x0 |y-x0| . The set of rays L CB (r, a, b) = {(x 0 (β), u(β, y ⊥ , y 3 )), β ∈ [0, 2π], (y ⊥ , y 3 ) ∈ [-a, a] × [-b, b]} denes the rays of a cone beam scan.
Here, for every β, the rays (x 0 (β), u(β, y ⊥ , y 3 )), (y ⊥ , y 3 ) ∈ [-a, a] × [-b, b], are rays of a perspective projection through x 0 (β): see Figure 4. Furthermore the optical center x 0 (β) scans the horizontal circle When r becomes large compared to the size of K, the perspective projections of K look like orthographic projections along parallel rays. Then we mention a scan with orthographic projections; this is often used in order to approximate perspective projections, when the distance scene-camera is large compared to the distances in the scene. Here we can see that the problem of 3D reconstruction from an orthographic scan L (r, a, b) can be reduced to problems of 2D reconstruction from orthographic scans in 2D. Indeed, if L = L (r, a, b), the knowledge of P[f ]| L is equivalent to the knowledge of the horizontal cross-sections:

{x 2 1 + x 2 2 = r 2 , x 3 = 0} when β scans [0, 2π]. θ θ ⊥ 0 e 3 u(β, y 2 , y 3 ) θ ⊥ x 0 (β) = rθ y = y ⊥ θ ⊥ + y 3 e 3 θ θ ⊥ 0 x 0 (β) = rθ(β) e 3 = θ ⊥ ∧ θ
P[f | Hy 3 ]| L∩Hy 3 , H y3 = {((x, y 3 ), (u, 0)), x ∈ R 2 , u ∈ S 1 }, y 3 ∈ [-b, b].
This motivates the problem of reconstruction from orthographic reective projections in 2D [START_REF] Bellet | Reective ltered backprojection[END_REF].

2.2. Tomography.

2.2.1. General formulation. The standard model of transmission tomography is based on the X-ray transform: Denition 7 (X-ray transform). Let f : R 3 → R (integrable over lines). The X-ray transform of f is:

R[f ] : R 3 × S 2 -→ R (x, u) -→ R f (x + tu)dt.
Here a point v gives a contribution for every ray of the set {L : v ∈ L}; a graphical representation of such tomographic data is often called a tomogram, or a sinogram.

Problem 2 (Inversion of X-rays). Let f : R 3 → R, and let L ⊂ R 3 × S 2 be a set of rays. For every L ∈ L,

we record R[f ](L). Question: reconstruct f from R[f ]| L .
This problem is now linear. We can nd in the literature ecient reconstruction procedures, at least for some sets of rays L, such as the rays of a cone beam scan, or the rays of an orthographic scan. The analogous bi-dimensional problem is the Radon inversion, which is often solved by the well-known ltered backprojection algorithm [START_REF] Natterer | Mathematical methods in image reconstruction[END_REF]. The FDK algorithm presented page 8 is a very popular heuristics which solves this problem. Roughly speaking, it provides a ltered backprojection operator A = BΦ such that

BΦR[f ]| LCB(r,a,b) ≈ f, f : R 3 → R.
It is also known that if r is large, then such an operator behaves like the ltered backprojection for orthographic projections in 2D, per horizontal cross-section. F w (β, y ⊥ , y 3 ) = w(y ⊥ , y 3 )F (x 0 (β), u(β, y ⊥ , y 3 )), w(y ⊥ , y 3 ) = r (r 2 +y ⊥ 2 +y3 2 ) 0.5 .

Let us consider now

Step 1.b) Compute the horizontal ltering ΦF : ΦF (β, y ⊥ , y 3 ) := F -1 (|σ| ĥ(σ)F(F w (β, y ⊥ , y 3 ))(σ))

y 3 ∈ [-b, b], β ∈ [0, 2π],
where F(g)(σ) = g(y ⊥ )e -iσy ⊥ dy ⊥ is the Fourier transform, and ĥ is an even windowing function with compact support.

Step 2) Compute the backprojection BΦF , where B is a weighted summation over lines through x:

BG(x) = 2π 0 r 2 (r -x • θ) 2 G(β, y ⊥ , y 3 )dβ, with y ⊥ = rx•θ ⊥ r-x•θ , y 3 = rx3 r-x•θ , θ = (cos β, sin β, 0).
Output. FDK reconstruction BΦF . 3.1.1. Reective projection. We consider Problem 1, where: K = {|x -z| r} is the ball of radius r > 0

centered at z ∈ R 3 , ∂K = {|x -z| = r} is the corresponding sphere, f : (v, u) ∈ ∂K × S 2 → 1 is constant, L = L CB (R, M, M ) with R, M > 0 (large enough).
In that case, the reective projection for a xed angle β 0 , P[f ]| L∩{β=β0} , represents the silhouette, or the shadow, of the object K: see Figure 5. It is clear that P[f ] and R[χ(K)] have the same support, and that P

[f ] = χ(supp R[χ(K)]
), where χ is the characteristic function and supp denotes the support [START_REF] Gardner | Geometric tomography[END_REF]. We use the parametrization of L CB (see Denition 5): a ray of the scan takes the form L(x 0 (β), u(β, y ⊥ , y

3 )) ∈ L, with (β, y ⊥ , y 3 ) ∈ [0, 2π] × [-M, M ] 2 . The projection is P[f ](β, y ⊥ , y 3 ) := P[f ](x 0 (β), u(β, y ⊥ , y 3 )).
As in Figure 3 ψ(β, y ⊥ , y 3 ) = |x 0 -z| 2 -|(x 0 -z) • u| 2 -r 2 ; we get:

(right), L∩K = ∅ ⇔ |x 0 -z| 2 -|(x 0 -z) • u| 2 r 2 ; in this case, the visible point is v K (x 0 , u) = x 0 + [(z -x 0 ) • u -C(x 0 , u)]u where C(x 0 , u) 0 is the unique number 0 such that v K (x 0 , u) ∈ ∂K. Let
P[f ](β, y ⊥ , y 3 ) = χ(ψ -1 ({t 0})(β, y ⊥ , y 3 ) = 1 if ψ(x 0 , u) 0, 0 otherwise.
In other words { P[f ] = 1} = {ψ 0}, and { P[f ] = 0} = {ψ > 0}.

3.1.2. Discontinuities. The interface between these two pieces is S := ∂{ψ 0} ∩ ∂{ψ > 0}. We show that S = {ψ = 0}. The inclusion S ⊂ {ψ = 0} is due to the continuity of ψ. By the way, ∂ψ ∂(y ⊥ ,y3) does not vanish on {ψ = 0}: by the chain rule,

∂ψ ∂(y ⊥ , y 3 ) = -2[(x 0 -z) • u](x 0 -z) T ∂u ∂(y ⊥ , y 3 ) , with ∂u ∂(y ⊥ , y 3 ) = 1 |y -x 0 | (I -uu T )[θ ⊥ , e 3 ],
where y = y ⊥ θ ⊥ + y 3 e 3 , and I -uu T is the orthogonal projection on u ⊥ . If R is large enough, then (x 0 -z) • u = 0; so ∂ψ ∂(y ⊥ ,y3) = 0 ⇔ (I -uu T )(x 0 -z) ⊥ span(θ ⊥ , e 3 ) = θ ⊥ . But due to the condition ψ = 0, x 0 -z and u are linearly independent and so (I

-uu T )(x 0 -z) ∈ u ⊥ \ {0}. Then ∂ψ ∂(y ⊥ , y 3 ) = 0 ⇔ u ∈ θ ⊥ .
But u / ∈ θ ⊥ and so ∂ψ ∂(y ⊥ ,y3) = 0 . Hence the sign of ψ changes across {ψ = 0}, and then {ψ = 0} ⊂ ∂{ψ 0} ⊂ S. This also proves that the rank of the jacobian of ψ is one on S. And thus S is a surface.

P [f ] is a piecewise constant function and its discontinuities occur at the surface S. In Figure 6, we get a reectogram of P[f ] by the means of three cross-sections of P[f ]. Here we have emphasized the interface S by drawing {|ψ| < τ } where τ > 0 is a small threshold. This interface S encodes the same information than P[f ], so it is an object of interest. One can check that

ψ = 0 ⇔ x 0 + uu T (z -x 0 ) ∈ ∂K. Hence S = {(β, y ⊥ , y 3 ) : (x 0 (β), u(β, y ⊥ , y 3 )) is tangent to ∂K}. If (β, y ⊥ , y 3 ) ∈ S, the visible point is v K = x 0 + uu T (z -x 0 ) (and u • (v K -z) = 0).
In Figure 3 (right), it means that v K and x 0 + uu T (z -x 0 ) are equal.

3.2. Heuristic reconstruction of the sphere. We now investigate the eect of the FDK algorithm on before ltering are on Figure 6.

P[f ].
S : ψ = 0 ψ > 0 ψ < 0 - - - - + + + + ΦP[f ] ≈ 0 ΦP[f ] :
(β) = (β, x⊥ (β), x3 (β)), β ∈ [0, 2π]; here, x⊥ (β) = Rx•θ ⊥ R-x•θ , x3 (β) = Rx3 R-x•
ψ • γ(β) = |x -z| 2 -|(x -z) • u| 2 -r 2 , with u = x-Rθ |x-Rθ| . Case a. Firstly, ψ • γ(β) = 0 ⇔ x + uu T (z -x) ∈ ∂K ⇔ the line (Rθ, x) is tangent to K. In particular, if the line (Rθ(β), x) is far from the tangents of K, then ΦP[f ](γ(β)
) is small and of arbitrary sign; it does not provide a signicant value at x.

Case b. Secondly, we study the case where the line (Rθ

(β 0 ), x) is tangent to K, i.e. (x, β 0 ) is such that ψ • γ(β 0 ) = 0. The derivative of ψ • γ is: d dβ ψ • γ(β) = -2[(x -z) • u] (x -z) • du dβ , du dβ = 1 |x -Rθ| (I -uu T )Rθ ⊥ . Then d dβ ψ(γ(β)) = 0 ⇔ (x -z) • u = 0 or θ ⊥ • (I -uu T )(x -z) = 0. Therefore d dβ ψ • γ(β 0 ) = 0 if,
and only if, x ∈ ∂K or the plane Rθ(β 0 ) + span(θ ⊥ (β 0 ), u(β 0 )) is tangent to K. See Figure 11. x 3

x 2

x 3

x 1 x 1

x 2 3): the ray (Rθ, u) and the plane Rθ + span(θ ⊥ , u) are tangent to K at x + uu T (z -u).

(1) If ψ(γ(β 0 )) = 0 but the other conditions are violated, then d dβ ψ • γ(β 0 ) = 0, the sign of ψ • γ(β) changes when β crosses β 0 , and the curve γ crosses the surface S at γ(β 0 ). In that case the (weighted) values of ΦP[f ](γ(β)) on either side of β 0 are expected to be high but with dierent signs; they should oset themselves by summation, and they should not contribute signicantly at x. See Figure 12 (left).

(2) On the contrary, if (Rθ(β 0 ), x) is tangent to K at x ∈ ∂K (left of Figure 11) and Rθ(β 0 ) + span(θ ⊥ (β 0 ), u(β 0 )) is not tangent to K, then the second derivative is:

d 2 ψ • γ dβ 2 | β=β0 = -2 (x -z) • du dβ | β=β0 2 < 0.
So the curve γ is tangent at S at the point γ(β 0 ); furthermore ψ • γ > 0 in a deleted neighborhood of β 0 . In that case the (weighted) values of ΦP[f ](γ(β)) on either side of β 0 are expected to be high but with the same sign; they should be accumulated by summation and give a signicant contribution (3) The last case occurs when the line (Rθ(β 0 ), x) and the plane Rθ(β 0 )+span(θ ⊥ (β 0 ), u(β 0 )) are tangent to K (right of Figure 11). Then ψ•γ(β 0 ) = d dβ ψ•γ(β 0 ) = 0 and we can expect signicant contributions as before, but now along the line (Rθ(β 0 ), x). This case describes eventual artifacts. Such lines are slightly remarkable in the vertical views of Figures 9 and10.

S : ψ = 0 ψ < 0 ψ > 0 ∇ψ + + + - - - γ γ - ΦP[f ] : S : ψ = 0 ψ < 0 ψ > 0 ∇ψ γ γ ΦP[f ] : + + + + Figure 12. Left -case b.(1): if ψ(γ(β 0 )) = 0 but dψ•γ dβ (β 0 ) = 0,

Numerical results

4.1. Contribution of the discontinuities. To observe the contribution of the discontinuities, we create synthetic data sets where the projected pattern f is piecewise constant with discontinuities. We increase the number of jumps from a data set to the next one.

Here, the shape ∂K is a sphere with a dent. To create this object, we deform the sphere |x| = 1, in spherical coordinates (ψ, ϕ, ρ), where ψ ∈ [-π, π] is the azimuth, ϕ ∈ [-π 2 , π 2 ] is the elevation, and ρ > 0 is the radius. For all points of the sphere (ψ, ϕ, ρ = 1), the point of the considered surface ∂K is (ψ, ϕ, ρ := 1 + 0.75(r -1)1 1 r<1 ), with 0.08r := ( ψ π + 1/4) 2 + ( 2ϕ π + 1/6) 2 . This object is computed from a discrete version of the sphere, discretized with 640 2 patches. For all integer m, we dene on this surface the following piecewise constant pattern, in spherical coordinates:

(ψ, ϕ) → p m (ψ)p m (ϕ), with p m (s) = 0.5 + 0.251 1 (ms-ms )<0.5 .

We project this pattern: for each ray L(x 0 , u), the projection is

P[f ](L) = f (v K , u) := p m (ψ)p m (ϕ),
where (ψ, ϕ, ρ) are the spherical coordinates of the visible point v K ∈ ∂K. We simulate an orthographic scan of this reective projection, using plots of surfaces with Matlab. We consider a uniform discretization (constant steps) of size 801 × 201 × 201, 801 being the number of angles of projection, and 201 × 201 being the size of each image.

Increasing m increases the number of jumps; we simulate data sets for several values of m: 0, 1, 2, 4, 8 and 16. On the rst line of Figure 13, we represent one image of the sequence, for the successive values of m.

Of course we distinguish here two kinds of jumps in the images: jumps due to discontinuities of the pattern f , and geometrical jumps due to the shape (interface object/background).

We apply the heuristics on these data sets. For the visualization, the reconstructions are re-projected, using a Maximum Intensity Projection (MIP). In order to improve the visual perception of the dent, we restrict the domain to a half-space before visualization: our full volume being a set of voxels (i, j, k) ∈ [1, 201] 3 , we keep only the j 91. On the second line of Figure 13, we represent a vertical view of the reconstructions, associated with the view of the data sets. More jumps in the input can improve the visual perception of the object, for both the input and the output. The heuristics computes much more than just a convex hull. A remarkable property here is that the dent is even more perceptible in the reconstruction than in the reectogram. This is because the boundary ( ψ π + 1/4) 2 + ( 2ϕ π + 1/6) 2 = 0.08 of the dent introduces geometrical jumps in many input images; the heuristics combines them, and so the boundary is emphasized in the reconstruction. We consider the sphere ∂K : |x| = 1, which is discretized using 640 2 patches as before. We project a smooth pattern, given in spherical coordinates by:

(ψ, ϕ) → 1 + 0.5 cos m(ψ + ϕ), with m = 2 i , 1 i 6, playing the role of a xed frequency. As before, by rotating over 360 degrees, we simulate an orthographic scan of size 801 × 201 × 201. The simulated images are smooth, except at the interface object/background where discontinuities occur. For the purpose of this subsection, we weight the data in order to obtain smooth images. For a pixel (i, j) ∈ [1, 201]× [START_REF] Baumgart | Geometric modeling for computer vision[END_REF]201], the weight is (r +1) 2 (r -1) 2 1 1 r<1 , with 99r(i, j) = ((i -100) 2 + (j -100) 2 ) 0.5 . On the rst line of Figure 14, we represent one smooth image of the sequence that we get, for the successive values of m.

We apply the heuristics on these data sets. On the second line of Figure 14, we represent a vertical view (MIP) of the reconstructions, associated with the view of the data sets. From the visual point of view, it is easier to identify the original sphere for large m. For smooth data, the result of the heuristic may be poor for low frequency components, but is meaningful when the frequency is large enough. In any case, smooth variations in the input produce contributions which add contrasts in the renderings; they can be useful for the visual perception. The results also emphasize the importance of singularities: here, we have seen that without the discontinuities at the interface object/background in the projections, it is more dicult to reconstruct the shape. 4.3. Test of robustness. We show how the method deals with changes in the forward problem during the acquisition, by considering a randomized uctuating pattern. This pattern could represent active surfaces whose reectance properties are random.

For all σ = 0, 2 j , -2 j 2, we consider an orthographic scan F σ (β, y ⊥ , y 3 ) of the Stanford Bunny ∂K [START_REF] Turk | Zippered polygon meshes from range images[END_REF]; the size of this scan is 801 × 200 × 157. For each ray of projection (x 0 , u) such that the angle of projection is β (and thus u = (-cos β, -sin β, 0), see Denition 6), the projected surfacic pattern is:

f σ (v, u) = 1 + (0.2 + ση 1 (β)) sin(πση 2 (β) + 20π |v|), v ∈ ∂K,
where the η i (β) are independent realizations of the Gaussian N (0, 1). The ση i (β) correspond to some modications of the amplitude and the phase of the pattern. Increasing σ increases the dependency in β for the projected pattern. If v ∈ ∂K is a xed point, then along the set of rays that sees v, i.e. {(β, y ⊥ , y 3 ) : v K (x 0 , u) = v}, the intensity level of v randomly varies with a standard deviation which increases when σ increases. To observe this, see the rst line of Figure 15, where we represent reectograms: horizontal slice y 3 = 0 in the reective projections.

We apply the heuristics on these data sets. On the second line of Figure 15, we represent a vertical view of the reconstructions (MIP). Despite changes in the forward problem, the silhouette of the object still appears under the form of discontinuities in the reective projection. As a result the heuristics successfully reconstructs the shape. By the way, the jumps at the discontinuities do not have the same values. Thus the reconstruction is better than the shape alone: surfacic contrasts are perceived. This is even true for large disturbances: see the last case, where the level of disturbance is about 4 times the intensity level of the original signal.

Conjectures

We would like to extend the previous discussions to more general congurations. The rst thing to do is to nd a framework which describes the structure of a reective projection, and more particularly its singularities. The following conjecture deals with an extension for several smooth convex objects:

Conjecture 1. Let (K i ) 1 i N be a family of compact sets such that Ki = K i and K i ∩ K j = ∅ if i = j. Let (f i ) 1 i N be a family of distinct numbers. We assume that the K i are convex and smooth. Let K = ∪ N i=1 K i , and let f : ∂K × S 2 → R, such that f | ∂Ki = f i , 1 i N . Let M, R > 0 be large enough. We conjecture that the reective projection P[f ]| LCB(R,M,M ) is piecewise constant. We conjecture that its singularities are included in the set of the rays of L CB (R, M, M ) that are tangent to one of the K i 's (at least), and that they describe pieces of surfaces. We assume that F is piecewise smooth. Let S ⊂ L CB (R, M, M ) denotes the set of rays L ∈ L CB (R, M, M ) such that F is singular at L. We conjecture that the FDK reconstruction BΦF emphasizes several points, including the points v K (L) ∈ ∂K, L ∈ S, which are the visible points at the singularities.

The meaning of emphasize must be precised. Essentially we have seen that the method consists in accumulating coherent information, and in cancelling incoherent one. We know that methods of stationary phase provide a mathematical framework to describe such phenomena. This option is explored in ongoing works: [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF] focuses on asymptotic expansions with respect to the cut-o frequency of the ltering (width of the spectral windowing). Another idea is to study the problem in a framework of microlocal analysis: describe the singularities of the reective projection and the singularities of the empirical reconstruction. Such studies will strengthen the heuristics by proving what it exactly does. In particular a full description of the emphasized points will also describe the artifacts: the emphasized points that do not belong to ∂K.

We see here that the subject meets again the problems of Radon inversion, where such ideas are currently developed to describe and reduce incompleteness artifacts [START_REF] Borg | Full Characterization of Reconstruction Artifacts from Arbitrary Incomplete X-ray CT Data[END_REF].

Figure 1 .

 1 Figure 1. A few images of a sequence of real laser images.

Figure 2 .

 2 Figure 2. Volume rendering of a 3D tomographic reconstruction from 2D laser images.

Denition 4 .

 4 Let F = P[f ] be a reective projection. A reectogram of F (or of f ) is a graphical representation of F . Such representations have peculiarities. If v ∈ ∂K, the set of rays whose visible point is v is {L : v K (L) = v}. If we follow this set through a reectogram of P[f ], we observe the intensities of projection of v: f (v, •). In general the set {L : v K (L) = v} does not cross the whole reectogram due to the denition of the visible point: occlusions appear.

Figure 4 .

 4 Figure 4. Cone beam scan: perspective rays (x 0 , u) through the optical center x 0 (left); the optical center x 0 (β) scans {x 2 1 + x 2 2 = r 2 , x 3 = 0} (right).

Denition 6 (

 6 Orthographic scan). Let r, a, b > 0. For every β ∈ [0, 2π], we set: θ(β) = (cos β, sin β, 0), θ ⊥ (β) = (sin β, -cos β, 0), e 3 = (0, 0, 1), u(β) = -θ(β). For every (y ⊥ , y 3 ) ∈ [-a, a] × [-b, b], we set x 0 (β, y ⊥ , y 3 ) = rθ(β)+y ⊥ θ ⊥ (β)+y 3 e 3 . The set of rays L (r, a, b) = {(x 0 (β, y ⊥ , y 3 ), u(β)), β ∈ [0, 2π], (y ⊥ , y 3 ) ∈ [-a, a] × [-b, b]} denes the rays of an orthographic scan.

Problem 1 and

 1 Problem 2 are not the same, but the geometry of the reective projection P[•]| L and the geometry of the X-ray transform R[•]| L have similarities: the scene is projected along the rays of L. The main idea of reective tomography is to use methods which solve Problem 2 in order to solve eciently Problem 1: Heuristics (Principle of reective tomography). We consider Problem 1. We assume that A is an operator which solves Problem 2 for the corresponding rays: AR[•]| L ≈ I (identity). Then the volumic function AP[f ]| L : R 3 → R is an empirical representation of the original set ∂K. The main question that arises now concerns the meaning and the description of AP[f ]| L . How does it represent the shape ∂K? 2.2.2. Case of a cone beam scan. We present the specic case of a cone beam scan. Problem (Cone Beam Computed Tomography). We measure partially the X-ray transform of f : R 3 → R, by cone beam scanning: R[f ]| LCB(r,a,b) . Question: reconstruct f from R[f ]| LCB(r,a,b) .

  Problem 1 with L = L CB (r, a, b) being the rays of a cone beam scan. We measure the reective projection P[f ]| LCB(r,a,b) where f : ∂K × S 2 → R and ∂K are unknown; we would like to recover ∂K. Following the principle of reective tomography, we use the reconstruction operator BΦ of Problem 2: the reconstructed ∂K is represented by the volumic function BΦP[f ]| LCB(r,a,b) : R 3 → R. This solves eciently a nonlinear geometric problem with a linear operator: this can be fastly computed on Graphics Processing Units. But proving mathematically in what extent BΦP[f ]| LCB(r,a,b) is a true representation of the original ∂K is still an open question. We now discuss this heuristic representation.

  FDK algorithm. Input. Projection F : L CB (r, a, b) → R, measured by a cone beam scan.Step 1.a) Compute the weighted data set F w :

Figure 5 .

 5 Figure 5. Silhouette of a ball: if K is a ball, with f = 1, the perspective projection P[f ](x 0 (β), u(β, •, •)) represents a silhouette and a shadow.

3. Reconstruction from silhouettes 3 . 1 .

 31 Cone beam scan of a ball.

Figure 6 .

 6 Figure 6. Reectogram from silhouettes of a ball: three cross-sections of P[f ](β, y ⊥ , y 3 ). In light gray, P[f ] = 0 and ψ > 0, in dark gray, P[f ] = 1 and ψ < 0; in black, the interface ψ = 0 (or more precisely: |ψ| < τ ).

Figure 7 .

 7 Figure 7. Cross-sections of the ltered projection ΦP[f ]; the same cross-sections for P [f ]

Figure 8 .

 8 Figure 8. The ltered projection ΦP[f ] looks like a zero-crossing detection of S : ψ = 0.

3. 2 . 1 .

 21 Filtering. The rst step computes the ltered projection ΦP[f ]: see Figure7for cross-sections. We observe that this step especially detects and emphasizes the singularity S. The lter F -1 (|σ| ĥ(σ))(y ⊥ ) is indeed the derivative ∂ y ⊥ , followed by a regularized Hilbert transform in y ⊥[START_REF] Natterer | Mathematical methods in image reconstruction[END_REF]. As a result, for all β and y 3 , this step emphasizes the contours of y ⊥ → P[f ](β, y ⊥ , y 3 ); it is a zero-crossing detection. The signicant values of ΦP[f ] are located near the surface S : ψ = 0 and the sign of ΦP[f ] is expected to change across S (zero-crossing). See Figure8. 3.2.2. Backprojection. The second step computes the backprojection BΦP[f ]; see Figure 9 for orthogonal cross-sections in this reconstruction, and Figure 10 for Maximum Intensity Projections. For every point x, the reconstruction BΦP[f ](x) is especially a summation through ΦP[f ], along the curve γ

  θ denotes the coordinates of the projection of x: Rθ, x and x⊥ θ ⊥ + x3 e 3 are aligned. In general this is a summation of arbitrary signed values. Theses values are small far from S, but are large in a neighborhood of S. The summation cancels the values in general. But for special circumstances, high values of the same sign are accumulated without being compensated and thus the sum is high. We want to identify such high contributions. So we study the cancellation of the function ψ along γ, by considering

2 Figure 9 .

 29 Figure 9. Cross-sections of the reconstruction BΦP[f ] of the sphere.

Figure 10 .Figure 11 .

 1011 Figure 10. Maximum Intensity Projections of the reconstruction BΦP[f ] of the sphere.

  Figure 12. Left -case b.(1): if ψ(γ(β 0 )) = 0 but dψ•γ dβ (β 0 ) = 0, the values of ΦP[f ](γ(β)) should oset themselves by summation. Right -case b.(2): if ψ(γ(β 0 )) = 0 and dψ•γ dβ (β 0 ) = 0, the values of ΦP[f ](γ(β)) should be accumulated by summation.

4. 2 .

 2 Contribution of the smooth variations. Realistic images of realistic scenes have often jumps and smooth variations. Here we would like to observe the contribution of the smooth variations alone, and to observe what happens without discontinuities. So we create synthetic data sets whose images are smooth, and whose frequency increases from a data set to the next one.

Figure 13 .

 13 Figure 13. Eect of jumps for a sphere with a dent. From the left to the right: the projected pattern has more and more discontinuities, m = 0, 2 i , 0 i 4. On the rst line: one image of the sequence; on the second line: associated heuristic reconstruction.

Figure 14 . 6 .

 146 Figure 14. Eect of smooth variations for a sphere. From the left to the right: the frequency parameter is larger and larger: m = 2 i , 1 i 6. On the rst line: one image of the sequence; on the second line: associated heuristic reconstruction.

Figure 15 . 2 .

 152 Figure 15. Test of robustness: from the left to the right, the level of disturbance is σ = 0, 2 j , -2 j 2. On the rst line, reectograms (β, y ⊥ ) → F σ (β, y ⊥ , 0). On the second line: a vertical view computed by the heuristics.

  The next step is to understand the eect of the heuristics. The following conjecture states that the method is based on the singularities of the data set: Conjecture 2. Let F = P[f ]| LCB(R,M,M ) be a cone-beam scan of the reective projection of f : ∂K ×S 2 → R.


Daniele Sportillo

Alexis Paljic

Luciano Ojeda

Philippe Fuchs

Vincent Roussarie

Psa Group

¶ Psa

Light Virtual Reality systems for the training of conditionally automated vehicle drivers

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Conditionally automated vehicles do not require drivers to constantly monitor their driving environment; drivers can, therefore, engage in secondary activities such as reading, writing emails and watching videos. However, when the automated system encounters unexpected situations, it will assume that drivers will adequately respond to a Take-Over Request (TOR).

In this context, Virtual Reality (VR) constitutes a potentially valuable learning and skill assessment tool because it would allow drivers to familiarize themselves with the automated vehicle and interact with the novel equipment involved in a free-risk environment.

The objective of this research is to explore the potential of light Virtual Reality systems, in particular, their role for the acquisition of skills for the Transfer of Control (ToC) in highly automated cars. By using the adjective light, we want to mark the difference between VR systems that are portable and/or easy to set up and systems that are cumbersome and require dedicated space to operate. The idea is that thanks to the portability and the cost-effectiveness, light VR systems could be easily deployed in car dealerships to train a large amount of people in a safe and reliable environment. This paper aims to compare the effectiveness of a training program based on a headmounted display with a user manual and with a fixed-base driving simulator. To validate the light VR system, user performances are evaluated during a test drive in a high-end driving simulator and self-reported measures are collected via questionnaires.

Related work

Although most of the studies in Driving Simulation (DS) uses static screens as the display system, recent studies prove that HMD-based * e-mail: daniele.sportillo@{mpsa.com,mines-paristech.fr} † e-mail: alexis.paljic@mines-paristech.fr ‡ e-mail: luciano.ojeda@mpsa.com § e-mail: philippe.fuchs@mines-paristech.fr ¶ e-mail: vincent.roussarie@mpsa.com DS leads to similar physiological response and lane change performance when compared to stereoscopic 3D or 2D screens [START_REF] Weidner | Comparing VR and non-VR driving simulations: An experimental user study[END_REF] [START_REF] Taheri | Development of a Driving Simulator with Analyzing Driver's Characteristics Based on a Virtual Reality Head Mounted Display[END_REF]. Even if the steering wheel is the most used driving interface, novel HMD systems usually come with wireless 6-DoF controllers which can be used to control a virtual car [START_REF] Sportillo | An immersive virtual reality system for semi-autonomous driving simulation: a comparison between realistic and 6-dof controller-based interaction[END_REF]. DS, in fact, provides the opportunity to implement critical scenarios and hazardous situations which are ethically not possible to evaluate on real roads [START_REF] Ihemedu-Steinke | Virtual reality driving simulator based on head-mounted displays[END_REF]. For this reason and to overcome the limited availability of physical prototypes for research purposes, DS is extensively used for studies on automated vehicles to design future automotive Human-Machine Interface (HMI) [START_REF] Melcher | Takeover requests for automated driving[END_REF] for TORs and to investigate the behavioral responses during the transition from automated to manual control [START_REF] Merat | Transition to manual: Driver behaviour when resuming control from a highly automated vehicle[END_REF]. Assessing the quality of the take-over performance remains an open problem. Usually reaction times (such as gaze reaction time, hands on wheel time, and intervention time), time to collision, lateral accelerations and minimum clearance are objective metrics analyzed in obstacle avoidance scenarios [START_REF] Happee | Take-over performance in evasive manoeuvres[END_REF].

The intense scientific production in recent years suggests that the transition of control in automated cars is a valuable research topic worth investigating from the design stage to the final implementation of the new systems. Moreover, the compelling need and interest of the car industry to train a large amount of people in a reliable and cost-effective way, without compromising security, make light virtual reality system tools a promising solution for this purpose.

METHODS

This study contained two parts: training and test drive. The aim of the training was to introduce the principles of the Level 3 Automated Driving System (ADS)-equipped vehicle, present the novel humanmachine interfaces, and describe the actions to perform in order to appropriately respond to unplanned requests to intervene. A between-subject study with 60 participants was designed in order to compare a light Virtual Reality systems to a user manual and a fixedbase driving simulator in terms of training effectiveness evaluated through a test drive which required the application of knowledge and skills acquired during the training.

The training

The aim of the training was to teach drivers how to interact with automated cars in three situations: the manual mode, automated mode and the take-over request. To do so, the training introduced the participants to the HMI for each situation, the actions they were free to perform during the automated driving and the best practice to respond to a take-over request. For all the participants, the training program started with an introductory video that briefly presented the main functionalities of a Level 3 ADS-equipped car.

In the study three different training systems were compared. The first one was a user manual (UM) consisting of a slide presentation displayed on a 13.3" screen of a laptop computer. The participants were asked to carefully read each of the 8 slides and to proceed when they felt ready, without a time limit. The slides used text and images to present the actions to be performed during the manual driving, the automated driving and the take-over requests.

The second system was a fixed-base simulator (FB) consisting of an actual car cockpit including a driving seat, a dashboard, a force-feedback steering wheel and a set of pedals. A 9.7" tablet used by the driver to perform the secondary activity was placed in the center console. To display the virtual environment a 65" plasma screen was positioned behind the cockpit at 1.5m from the driver.

The third system was a Light Virtual Reality (LVR) system including an HMD as a display system, and a racing wheel as driving system. Spatial sound was presented via headphones. To have a spatial correspondence between the real steering wheel and the virtual one, the steering wheel inside the virtual car was a 3D model of the real racing wheel with which the participants were interacting. Moreover, the position and the movements of the virtual model corresponded to the real one, allowing for co-located manipulation.

The Virtual Learning Environment

For the training using the LVR system and the fixed-base driving simulator, a step-by-step tutorial was developed in the form of a Virtual Learning Environment (VLE) (Fig. 1). The task of the participants consisted of interactions with the car following the instruction of a virtual vocal assistant. The messages announced by the assistant were also displayed on a yellow panel in front of the trainee. The driving scenario was a straight 2-lane road delimited by guardrails with no traffic. Before the driving scenario, an acclimatization virtual environment was proposed to the participants to help them locate and identify the controls of the car. This training also included a secondary activity (a video) that required the use of a tablet to distract the human driver from the driving task during the automated driving. The participants were asked, but not forced, to look at the tablet. The video was automatically played when the automated system was enabled and paused during the manual driving and the TORs.

Test drive

After the training, the participant performed a test drive designed to evaluate their performance in a more realistic driving scenario. The system used for this purpose was a high-end driving simulator consisting of the front part of a real car surrounded by a panoramic display. Data including position, speed and acceleration of the car, and current driving mode were recorded. Inside the car, a 10.8 inch tablet which provided 9 different secondary activities was placed in the center console. During the test drive 3 TORs were issued during the test drive: (A) a 10-second TOR caused by a road narrowing provoked by a stationary car on the right lane; (B) a 10-second TOR caused by a loss of ground marking; (C) a 5-second TOR caused by a sensor failure. During the autonomous driving, participants were asked to engage in one of the secondary activities proposed by the tablet.

RESULTS

To evaluate the training systems and the learning environment, objective and self-reported measures were collected anonymously and treated confidentially. A first outcome of the study is that the training allowed all the participants to respond to the TORs.

The quality of the take-over was evaluated in terms of reaction time (rt), the elapsed time from the TOR until the driver takes back control. A significant difference between the user manual group (rt = 5.15s) and the other two systems was observed. However, no differences were observed comparing the fixed-base simulator (rt = 3.17s) and the HMD (rt = 3.16s). Self-reported measures were collected through set of questions at the beginning of the test, after the training and after the test drive. To evaluate the favorability of the training, the participants filled out a 10 questions survey containing questions about perceived usefulness, easiness, pleasantness and realism. In summary, according to the objective metrics measured during the test drive, the group of participants trained with the Virtual Learning Environment (fixed-base driving simulator and light VR system) were able to respond to the take-over request in less time than the group of participants trained with the user manual. Furthermore, self-reported measures showed responses in favor of the light VR training system.

CONCLUSION

The results of this research persuade us that Light Virtual Reality systems represent a valuable tool for the acquisition of driving skills in conditionally automated vehicles. The proposed training system, composed of an HMD and a game racing wheel, is a portable and cost-effective device that provides an adequate level of immersion for teaching drivers how to respond to a take-over request in a safe environment. Therefore, this system could be employed for the training of future customers of automated cars before their first ride. A direct outcome of these results is the acknowledgment of VR as key player in the definition of the set of metrics for profiling driver interaction in automated vehicles. In future work, the training will be implemented in the form of a serious game in which the level of instruction adapts to the users' needs in order to assess the acquisition of skill during the training itself. Furthermore, test drives with the real vehicles are considered of primary importance to validate current results.

Figure 1 :

 1 Figure 1: Illustration of a user with HMD immersed in the Virtual Leaning Environment

ACKNOWLEDGMENTS

This research was supported by the French Foundation of Technological Research under grant CIFRE 2015/1392 for the doctoral work of D. Sportillo at PSA Group.