Viviana Contreras Moreno

Alain Ledoux
email: alain.ledoux@insa-rouen.fr

Lionel Estel

Rabah Azzouz

Vincent Levacher

Laurent Bischoff

Salim Derrouiche

Kinetic modeling of carbon dioxide valorization with epoxides

Keywords: Carbon dioxide valorization, cyclic carbonates, non-metallic catalyst, modeling, genetic algorithm

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Influence of atmospheric concentration of carbon dioxide in global warming is now well-established and is caused by an increase of greenhouse effect according to an acceleration of fossil energy consumption. Several ways are being implemented in the mitigation of carbon dioxide emissions like pre-combustion capture, oxyfuel combustion and post-combustion capture. These technologies are not at the same level of maturity but they can only be acceptable for a transition period. At present, only the long-term storage of carbon dioxide captured has been considered in the large scale. However the capture and storage of carbon dioxide have been questioned on the last years due to economics and environmental problems, which are associated with the energy consumed in the process and the stability of the stored gases. About 20-30% of energy in a power station is required to store carbon dioxide [START_REF] Metcalfe | An Integrated Approach to Energy and Chemicals Production[END_REF].

During the transition period valorization of CO2 should be developed. The objective is to consider carbon dioxide as a starting material needed by the chemical industry. In this way, carbon dioxide transformation would avoid the costs of carbon capture and storage associated [START_REF] Metcalfe | Influence of Reactor Design on Cyclic Carbonate Synthesis Catalysed by a Bimetallic Aluminium(salen) Complex[END_REF]. Several chemical process for carbon dioxide valorization have been studied, the most common processes include the synthesis of urea, urethane derivatives, carboxylic acid, dimethyl carbonate, polycarbonates, cyclic carbonates and the catalytic hydrogenation of carbon dioxide [START_REF] Olajire | Valorization of Greenhouse Carbon Dioxide Emissions into Value-Added Products by Catalytic Processes[END_REF].

Chemical valorization of carbon dioxide

Carbon dioxide transformation is not simple due to its high thermodynamic stability and low reactivity [START_REF] Sakakura | Transformation of Carbon Dioxide[END_REF]. In most cases high input energy is required to convert carbon dioxide or the use of complex catalytic systems is necessary to reduce the activation energy of the reaction. There are some exothermic systems which can be used for the chemical valorization of carbon dioxide such as epoxides, hydrogen, acetylene, amines, etc. Epoxides seem potential candidates for converting carbon dioxide into high value products.

The reaction between carbon dioxide and epoxides in presence of a catalyst (Figure 1) can produce two types of products: cyclic carbonates and polycarbonates. Typically, the formation of polycarbonates is less thermodynamically favorable than the synthesis of cyclic carbonates, but it requires lower activation energy. Therefore, the formation of any of these products will depend on the selectivity of catalytic systems, substrate used (type of epoxide) and conditions of the reaction. Cyclic carbonates has attracted much attention in the chemical industry due to the large number of applications [START_REF] Pescarmona | Challenges in the Catalytic Synthesis of Cyclic and Polymeric Carbonates from Epoxides and CO2[END_REF].

Figure 1 -Synthesis of cyclic carbonates from epoxides and CO2

Cyclic carbonates are mainly used as intermediates in the synthesis of fine or bulk chemicals (polymers), polar aprotic green solvents and as electrolyte solvents for lithium-ion batteries. The commercial production of organic carbonates mostly relies on the use of phosgene, which is toxic and dangerous [START_REF] Abbas | Organic Carbonates[END_REF]. Only minor amounts are produced by other pathways as oxidative carbonylation or by using epoxides and carbon dioxide (Aresta 2003). The last option requires high temperatures, above 150°C and high carbon dioxide pressures between 50 -80 bar [START_REF] Metcalfe | An Integrated Approach to Energy and Chemicals Production[END_REF].

Different types of homogeneous and heterogeneous catalytic systems have been suggested for the synthesis of cyclic carbonates. The recent works are focusing on new catalytic systems to avoid the use of large amounts of energy in the reaction between epoxides and carbon dioxide [START_REF] Pescarmona | Challenges in the Catalytic Synthesis of Cyclic and Polymeric Carbonates from Epoxides and CO2[END_REF].

In this work the kinetic of a non-metallic catalytic system was investigated using epichlorohydrin and carbon dioxide under mild conditions of pressure and temperature. Since the reaction is biphasic, the mass transfer and the reaction mechanism should be studied in order to model the system. The kinetic model was developed considering the reaction mechanism and the mass transfer in a stirred tank reactor. The modeling of system using the approach of quasi steady state was considered and compared with the complex model.

Experimental Set-up

The experiments between carbon dioxide and epoxides were carried out in a stirred tank reactor at several temperatures (50°C -90°C) and pressures (P < 5 bar). The reactor is composed of a jacket which keeps constant the temperature inside of the vessel by using a heat carrier fluid. The vessel (Tr) and jacket (Tj) temperatures are regulated to an accuracy of ±0.1°C. A calorimetric approach is possible by using a heater installed in the system that quantifies the amount of energy released or absorbed by the mixture. In order to obtain a heat balance in the reactor, a calibration is conducted to determine global heat coefficient (U W/m²/K) and heat capacity for the mixture (C J/kg/K). A stirrer in the reactor allows keeping the homogeneity in the liquid phase and it improves the gas transfer to the liquid. A gas inlet is available in the reactor which is linked to a gas chamber where the pressure and temperature of gas are measured as a function of time. The gas feeding is controlled from the regulation of pressure inside of the reactor. The scheme of the reactor is shown in Figure 2.

In a typical reaction sequence, the vacuum is established in the reactor, followed by the catalyst (99%) and epichlorohydrin (99%) introduction and then stirring of the mixture. The catalyst concentration used is lower than 7.0% mol. The reactor is heated to the desired temperature which is kept constant throughout the reaction. Carbon dioxide is introduced continuously to assure a constant pressure inside of the reactor. At specific time intervals, samples are withdrawn for further analysis. The samples are analyzed by gas chromatography (GC) and infrared spectroscopy (FTIR).

O R 1 R 2 CO 2 + O O O R 1 R 2 Catalyst Results
The substrate used was a mixture of epichlorohydrin with non-metallic catalyst. Effects of stirring, catalyst concentration and temperature were studied for modeling of the system.

Effect of stirring

The mass transfer limitation should be studied in order to correlate the kinetic parameters. The interaction between carbon dioxide and epichlorohydrin was studied at several stirring speeds and temperatures. In the absence of catalyst, carbon dioxide and epichlorohydrin present a low reactivity, therefore the effect of the mass transfer phenomena can be observed. We assumed that at low concentrations of catalyst, the physical properties of epichlorohydrin do not have a meaningful influence, thus the effect of stirring in the physical system can be extrapolated to the chemical system. The results indicate that the mass transfer is increased by stirring enhancement which is reflected on the mass transfer coefficient of the liquid phase (Figure 3). The mass transfer coefficient increases with temperature, first by a decrease of liquid viscosity and second by an increase of gas diffusion into the liquid film. Since the mass transfer phenomena in our case was faster than the rate of reaction, the reaction regime is directly controlled by the chemical kinetic but the hydrodynamic of the system must be considered in the model. The reactions were carried out under solvent-free conditions with epichlorohydrin as substrate and using 2.5% mol to 7% mol of non-metallic catalyst at temperatures between 50°C and 90°C. It was noticed that the reactions carried out at low concentrations of catalyst have an induction period (Figure 4-a) at the beginning of the process, though this was not apparent in reactions carried out at 7% mol of non-metallic catalyst. Additionally, at high temperatures the induction period was not observed even at low concentration of catalyst (Figure 4-b). Apparently at low temperatures and low catalyst concentrations, the activity of catalyst and its interaction rate in the system is slowest and we suppose according to the experimental data that the catalyst can be involved in the mechanism in a different form.

The conversion of epichlorohydrin was up to 90% at 70°C and 90°C after 480 min and 240 min of reaction respectively and the yields achieved were higher than 85%. The results suggest that one part of epichlorohydrin is converted into another product different from the cyclic carbonate. This aspect can have an impact in the modeling, and it must be considered.

On the other hand, we notice that the temperature has an accelerator role in the reaction and is more important than catalyst concentration. At high temperature the catalyst concentration does not have a meaningful impact in the consumption rate of epichlorohydrin.

Analysis of reaction mechanism

Organic cations can be used as catalyst for the reaction of carbon dioxide with epoxides. The most common reaction mechanism involves a nucleophilic attack of the catalyst which open the epoxide ring followed by the insertion of carbon dioxide to form the cyclic carbonate (Figure 5-a). In other way if the catalyst molecule is bulky, the attack can occur on carbon dioxide molecule followed by the interaction with epoxide and finally the formation of cyclic carbonate (Figure 5-b). In order to analyze the impact of reaction mechanism in the modeling, we compared the results obtained with both ways.

(A) (B)

Figure 5 -Reaction mechanism between CO2 and epoxides, with CO2 activation at the first step (A) or with epoxides activation at the first step (B) (Pescarmona and Taherimehr 2012)

Kinetic modeling

The main reactions leading to the formation of cyclic carbonates, such as carbon dioxide absorption, interaction between catalyst and carbon dioxide or catalyst and epoxide and the subsequent formation of cyclic carbonate are presented considering the two possible ways:

Way A: Catalyst/CO2 Way B: Catalyst/Epoxide Step1: 𝐶𝑂 2 (𝑔) ↔ 𝐶𝑂 2 (𝑙)
Step 2: 𝐶𝑂 2 (𝑙) + 𝐶𝐴𝑇 → 𝐼 1

Step 3:

𝐼 1 + 𝐸𝑃𝑂𝑋 → 𝐼 2 Step 4: 𝐼 2 → 𝐶𝐴𝑇 + 𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒 Step1: 𝐶𝑂 2 (𝑔) ↔ 𝐶𝑂 2 (𝑙)
Step 2: 𝐸𝑃𝑂𝑋 + 𝐶𝐴𝑇 → 𝐼 1

Step 3: 𝐼 1 + 𝐶𝑂 2 → 𝐼 2

Step 4: 𝐼 2 → 𝐶𝐴𝑇 + 𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒 We mentioned above that the steps 1 and 2 in the ways A and B change according to the interaction with the catalyst system, but the structure of mathematical model is the same. Thus, we present below the detailed mathematical model for the way A. Regardless of the way considered for the analysis the first order for each species was observed. According to the elementary reactions, the rate of formation of epichlorohydrin carbonate and intermediates and the rate of consumption of epichlorohydrin in a batch reactor with continuous feeding of carbon dioxide can be written as follows:

𝑑[𝐶𝑂 2 ] 𝑙 𝑑𝑡 = 𝑎 𝑖 𝑘 𝐿 𝐸 ( 𝑃 𝑖 𝐻𝑒 -[𝐶𝑂 2 ] 𝑙 ) -𝑘 1 [𝐶𝐴𝑇][𝐶𝑂 2 ] 𝑙 (1) 𝑑[𝐶𝐴𝑇] 𝑙 𝑑𝑡 = -𝑘 1 [𝐶𝐴𝑇][𝐶𝑂 2 ] l + 𝑘 3 [I 2 ]
(2)

𝑑[𝐼 1 ] 𝑙 𝑑𝑡 = 𝑘 1 [𝐶𝐴𝑇][𝐶𝑂 2 ] 𝑙 -𝑘 2 [I 1 ][𝐸𝑃𝑂𝑋] (3) 𝑑[𝐸𝑃𝑂𝑋] 𝑙 𝑑𝑡 = -𝑘 2 [I 1 ][𝐸𝑃𝑂𝑋] (4) 𝑑[𝐼 2 ] 𝑙 𝑑𝑡 = 𝑘 2 [I 1 ][𝐸𝑃𝑂𝑋] -𝑘 3 [I 2 ] (5) 𝑑[𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒] 𝑙 𝑑𝑡 = 𝑘 3 [I 2 ] (6)
Figure 6 shows the profiles of concentration for every specie involving on the detailed reaction mechanism which consider the catalyst interaction with carbon dioxide at the first step (Way A). Two cases are presented with and without the induction period at low (70°C) and high (90°C) temperature respectively.

We notice that the amount of epichlorohydrin carbonate increased gradually as the reaction proceeds. The attainment of equilibrium is produced after 480 min at 70°C and 300 min at 90°C. The results suggest that at low temperature, carbon dioxide saturates the liquid phase during the induction period. Then its concentration decreases once the acceleration phase begins. At high temperature this phenomena is present but is less observable because the reaction rate is faster than at low temperature. The intermediates concentrations are present in a small proportion to promote the reaction, we suppose that catalyst remains active along the reaction and only when the system is in the equilibrium the catalyst is recuperated.

(b) Figure 6 -Concentration profiles of species involving in the reaction of carbon dioxide with epoxides for detailed reaction mechanism -3 bar of CO2, 5% mol of non-metallic catalyst (a) 70°C (b) 90°C

Parameter estimation The numerical parameter estimation was based on the kinetic experiments where temperature and pressure were maintained constant. The model describing the reaction is composed of six differential equations where the kinetic rate constants should be estimated. During reactions, the intermediates concentrations were not observed neither quantified due to their low stability in the reactive system. In order to simplify the solution of the model, we have supposed that the production rate of intermediates is lower than the consumption rate. Otherwise a quasi-steady state (QSSA) could be assumed for the intermediates.

However the quasi steady state does not be applied on two intermediates. Since the system present an induction period, the intermediate formed in the first step with catalyst has an important role during the reaction and its concentration along time should be consider. This intermediate depends on the activation of carbon dioxide or epichlorohydrin (step 1). Therefore the quasi steady state can be only applied on intermediate 2.

In the quasi steady state the concentration of intermediate 2 is written in the form of an algebraic equation (Eq. 7):

[I 2 ] = 𝑘 2 [I 1 ][𝐸𝑃𝑂𝑋] 𝑘 3 (7)
The number of differential equations describing the system is reduced to five with an algebraic equation. The model predictions of the kinetic parameters were obtained by solving the ordinary differential equations (Eqs. 1-6) which describe the mass balance. A prediction of the components was obtained by applying a numerical differentiation method using the software Matlab ®. A genetic algorithm was used as optimization method to minimize the normalized root mean square error (NRMSE -Eq. 8):

𝑁𝑅𝑀𝑆𝐸 = √ ∑(𝑌 𝑖,𝑒𝑥𝑝 -𝑌 𝑖,𝑐𝑎𝑙𝑐 ) 2 𝑁 𝑌 𝑒𝑥𝑝,𝑚𝑎𝑥 -𝑌 𝑒𝑥𝑝,𝑚𝑖𝑛 (8)
Experimental concentrations 𝑌 𝑖,𝑒𝑥𝑝 were obtained from measurements and estimated concentrations 𝑌 𝑖,𝑐𝑎𝑙𝑐 were obtained from the models. N is the number of experimental points for each experiment and 𝑌 𝑒𝑥𝑝,𝑚𝑎𝑥 and 𝑌 𝑒𝑥𝑝,𝑚𝑖𝑛 are the maximal and minimal experimental concentration value respectively.

The curves of epoxide concentration as a function of time for the experimental points and estimated values from the model by using QSSA and the detailed model (way A) are presented in Figure 7. We compared the two simplified ways of catalyst activation, with the experimental data. It can be noticed that at low temperature only the detailed model can predict the induction period, but it presents some deviations at the equilibrium at high temperature (90°C). As we mentioned before, the kinetic results shown that one part of epichlorohydrin was not used to produce cyclic carbonate, therefore a secondary reaction can take place in the system modifying the equilibrium. The estimated values for the activation energy and the pre exponential factor are summarized in table 1.

CONCLUSION

The kinetics of synthesis of epichlorohydrin carbonate from carbon dioxide and epichlorohydrin in presence of a non-metallic catalyst was studied in the temperature range 50°C-90°C. Mass transfer effect was observed with the stirring speed, therefore the system depends on the hydrodynamic of the reactor. The temperature has an important role in the reaction rate, in agreement with the Arrhenius correlation. A kinetic model taking into account the mass transfer and the kinetic of the reaction was developed for the two activation ways. A simplification was introduced by using the quasi steady state on the intermediate 2. The detailed model describes accurately the induction period at low temperature. At high temperatures the models of the way A and the way B and the detailed model presents a similar behavior, however a deviation from the equilibrium is observed at the end of the reaction. We suppose that some phenomena in the equilibrium have not been taken into account in the model. Next work will consist in improving the modeling including the effect of the equilibria concentrations.

Figure 2 -

 2 Figure 2 -Schematic of stirred tank reactor GO -gas outlet; GI -Gas inlet; LI/LO -Liquid inlet/outlet; TStemperature sensor; PSpressure sensor; Qheater; Tsrtemperature sensor in reactor; Tsjtemperature sensor in jacket

Figure 3 -

 3 Figure 3 -Effect of temperature and stirring on the mass transfer coefficient of carbon dioxide into epoxides Effect of temperature and catalyst concentration

  Decrease in epichlorohydrin concentration with time at various concentration of non-metallic catalyst (a) 70°C (b) 90°C

  Fit of the detailed and simplified kinetic model with the experimental data at (a) 60°C (b) 70°C (c) 80°C (d) 90°C -3 bar of CO2 and 5% mol of catalyst.

Table 1 .

 1 Kinetic parameters estimated

		Pre-exponential factor	Activation Energy
	Reaction	A	(kJ/mol)
		Detailed method Way A
	𝐶𝑂 2 (𝑙) + 𝐶𝐴𝑇 → 𝐼 1	8.48E+17	160.66
	𝐼 1 + 𝐸𝑃𝑂𝑋 → 𝐼 2	4.09E+12	122.46
	𝐼 2 → 𝐶𝐴𝑇 + 𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒	7.31E+06	70.20
		QSSA Way A
	𝐶𝑂 2 (𝑙) + 𝐶𝐴𝑇 → 𝐼 1	1.35E+17	144.27
	𝐼 1 + 𝐸𝑃𝑂𝑋 → 𝐶𝐴𝑇 +
	𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒	3.81E+26	217.34
		QSSA Way B
	𝐸𝑃𝑂𝑋 + 𝐶𝐴𝑇 → 𝐼 1	8.37E+07	96.56
	𝐼 1 + 𝐶𝑂 2 → 𝐶𝐴𝑇 + 𝐶𝑎𝑟𝑏𝑜𝑛𝑎𝑡𝑒	1.28E+25	194.92


domain: math math.math-na sdv shs

text:
Maya De Buhan
email: maya.de-buhan@parisdescartes.fr

Chiara Nardoni
email: chiara.nardoni@upmc.fr

A facial reconstruction method based on new mesh deformation techniques

Keywords: forensic science, facial reconstruction, shape matching, elasticity, soft tissue deformation, finite elements

This article presents a new numerical method for facial reconstruction. The problem is the following: given a dry skull, reconstruct a virtual face that would help in the identification of the subject. The approach combines classical features as the use of a skulls/faces database and more original aspects: (i) an original shape matching method is used to link the unknown skull to the database templates; (ii) the final face is seen as an elastic 3d mask that is deformed and adapted onto the unknown skull. In this method the skull is considered as a whole surface and not restricted to some anatomical landmarks, allowing a dense description of the skull/face relationship. Also, the approach is fully automated. Various results are presented to show its efficiency.

Introduction

Facial reconstruction aims at recovering the facial appearance of an individual from the sole datum of the underlying skull. The facial reconstruction problem arises in various application fields like forensics, anthropology, archeology or history. In forensic science, facial reconstruction comes in the process of identification of deceased people. When all the usual methods of identification have failed and the skeletal remains are the sole data available for leading to a positive identification, facial reconstruction might be considered as an enhancing tool for 'recognition', producing a short list of candidates from which the individual may be identified by other endorsed methods of identification, see Wilkinson [START_REF] Wilkinson | Computerized forensic facial reconstruction: A review of current systems[END_REF]. In archeological investigations, facial reconstruction is employed with the purpose of identifying skeletal remains of famous people from the past.

The creation of the face from the skull is a procedure of approximation: from the observation of the cranium, one will not be able to recover a large amount of face features (eyes, hair, lips, ears). Moreover the facial likeness of a given individual changes considerably depending on factors like nutrition or aging. This flexibility may not be fully reflected on the subjacent skull. Any facial reconstruction tool is expected to account for the uncertainty related to the ill-posedness of the problem, no matter what method is employed (artistic, parametric, statistic, mechanical, etc.). From a mathematical point of view, this issue leads to at least two important difficulties: on the one hand, it raises the question of how to correctly characterize the solution, which might be a continuum spectrum of all the faces that are 'consistent' with a given skull rather than a single exemplar; on the other hand, it poses the problem of how to rigorously assess the accuracy of the result. Despite the intrinsic difficulty of the problem, the media are full of facial images that have been constructed on the basis of a single given skull. A fascinating survey of such cases can be found in the book of Prag and Neave [START_REF] Prag | Making Faces: Using Forensic and Archaeological Evidence[END_REF]. The work presented in this paper is part of the ongoing multi-disciplinary project ***.

The traditional facial reconstruction methods are based on manual procedures, producing 2d portraits or 3d sculptures. These methods basically consist of three common steps: (i) equip (a replica of) the raw skull with a sparse set of anatomical landmarks; (ii) apply an average soft tissue thickness to each skull landmark in order to estimate a corresponding landmark on the face; (iii) draw up or sculpt a face fitting the estimated landmarks. Most practitioners add a face muscle model in order to enrich the anatomical accuracy of the reconstruction, leading to the socalled Manchester method described in Wilkinson [START_REF] Wilkinson | Computerized forensic facial reconstruction: A review of current systems[END_REF]. The results obtained from forensic art are often quite plausible, as the medical artists may take anatomical, historical, archaeological or other types of expertise into account, giving the observer a feeling of coherence. However, the final result of a manual reconstruction depends on the subjectivity of the artist. Additionally, one single reconstruction requires several days of work of a well-experienced forensic artist, making impracticable the realization of multiple instances and feature variations. In order to alleviate these difficulties, several computer graphic software packages have been developed. These animation software packages use the same methodology as manual methods, allowing the expert to tune some modeling parameters and combine the human expertise with the flexibility of the software, as in Miyasaka et al. [START_REF] Miyasaka | The computer-aided facial reconstruction system[END_REF]. However, this approach does not eliminate the subjectivity in the reconstruction. During the last 30 years a large amount of work has been devoted to the conception of objective fully-automated methods. The common pipeline of modern facial reconstruction software is described in Claes et al. [START_REF] Claes | Computerized craniofacial reconstruction: Conceptual framework and review[END_REF]. First, an expert examines the unknown skull in order to determine anthropological parameters like age, gender and ethnicity. Then a virtual replica of the dry skull is produced and represented according to the modeling parameters. A craniofacial template encoding face, skull and soft tissue information is derived from a head database. Then an admissible geometric transformation drives the adaptation of the craniofacial template onto the unknown skull, according to the 'proximity' between the skulls. As a result, the template face is deformed onto the predicted face associated with the unknown skull, linking together information coming from both the database and the examination of the unknown skull. Finally a skin texture and hairiness are added to the reconstructed face.

The approach proposed in this article combines classical features -as the anthropological expertise for classifying the unknown skull and the use of a head database -with mathematical and computational skills as 3d geometric modeling. The contribution of this work is threefold. First, a method for generating a closed surface mesh model of the skull template is introduced. The method relies on an original mesh evolution technique. A template shape is iteratively deformed, producing a sequence of shapes that get 'closer and closer' to the source skull. Second, the elastic shape matching method first introduced in *** [5] is used. This technique is used to link skulls (or faces) with each other and learn about their similarities. This procedure is the heart of the reconstruction method and its most original feature. In particular it allows the method to be fully-automated by removing the need of landmarks. Third, the shape matching tool is combined with soft tissue deformation techniques from computational surgery for transporting the skull/face templates onto the unknown skull. The common denominator of the three parts is the use of elasticity equations for driving the shape deformation.

The paper is organized as follows: Section 1 deals with the acquisition and the construction of the database; Section 2 describes the representation of the skull templates; Section 3 presents the method for matching skull and face templates; Section 4 describes the process of reconstruction of the unknown face. In Section 5 some facial reconstructions are presented to prove the efficiency of the approach. Figure 1 summarizes the various steps of the procedure.

Data acquisition

Some recent numerical facial reconstruction methods rely on the a priori information contained in a database of coupled skull and face templates. The acquisition of both skull and face is accomplished by head CT scans of living subjects, allowing a good visualization of hard tissues. Standard segmentation tools and 3d reconstruction algorithms lead to the definition of dense surfaces of both skull and face from CT data. Unfortunately the invasiveness of this technique causes serious legal and ethical problems, preventing the constitution of a large database of healthy subjects. Due to this difficulty, several studies have proposed to exploit the relationship between soft and hard tissues by means of average soft tissue thickness measurements, for example Tilotta et al. [START_REF] Tilotta | A local technique based on vectorized surfaces for craniofacial reconstruction[END_REF], Hwang et al. [START_REF] Hwang | Facial Soft Tissue Thickness Database for Craniofacial Reconstruction in Korean Adults[END_REF] and De Greef et al. [START_REF] Greef | Large-scale in-vivo Caucasian facial soft tissue thickness database for craniofacial reconstruction[END_REF].

Soft tissue depth tables are usually used in combination with a large database of face templates as in Claes et al. [START_REF] Claes | Computerized craniofacial reconstruction: Conceptual framework and review[END_REF]. Facial templates can be acquired by non-invasive techniques such as stereophotogrammetry, Grewe et al. [START_REF] Grewe | Fast and Accurate Digital Morphometry of Facial Expressions[END_REF], making the constitution of a large database painless and easy. However, despite the acquisition of dense surface templates for describing the outer face, the average values of soft tissue thickness are systematically measured on a sparse set (< 53) of (manually positioned) anatomical landmarks. See Vezzetti et al. [START_REF] Vezzetti | Facial Landmarks for Forensic Skull-Based 3D Face Reconstruction: A Literature Review[END_REF] for a review of recent 3D landmark-based facial reconstruction methods. Since the manual measurement is time consuming and requires expertise in correctly identifying the landmarks, it is actually infeasible to extend these measurements to a dense distribution of points. See Kustar et al. [START_REF] Kustár | FACE R A 3D Database of 400 Living Individuals Full Head CT and Face Scans and Preliminary GMM Analysis for Craniofacial Reconstruction[END_REF] for a discussion about the use of sparse soft tissue measurements for facial reconstruction purposes. Several authors have claimed the importance of using a dense representation of the soft tissue information, for example by describing the intra-subject correlation in terms of the volume between the two boundary surfaces representing the face and the underlying skull, Nelson and Michael [START_REF] Nelson | The application of volume deformation to three dimensional facial reconstruction: A comparison with previous techniques[END_REF], Quatrehomme et al. [START_REF] Quatrehomme | A Fully Three-Dimensional Method for Facial Reconstruction Based on Deformable Models[END_REF]. The results presented in this study are based on a collection of 26 head CT scans of female healthy subjects aged between 20 and 40 years. The CT images used have been provided by the Statistical Facial Reconstruction project of Paris Descartes University presented in Tilotta et al. [START_REF] Tilotta | Construction and analysis of a head CT-scan database for craniofacial reconstruction[END_REF]. The 3d reconstruction from tomographic data described in this section, combined with an original wrapping method (see Section 2), allowed us to define closed surface models of both face and skull.

Image segmentation

Data segmentation consists in identifying the bone and soft tissues on stacked 2d gray-level images. This procedure has been carried out in a semi-automated way with the help of the software Amira developed by Stalling et al. [START_REF] Stalling | Amira: a Highly Interactive System for Visual Data Analysis[END_REF]. First, the CT slices are automatically presegmented using a multi-threshold technique. This step consists in partitioning the original images into subdomains whose boundaries are identified by given intensity values. The bone and soft tissues threshold values used are described in Tilotta et al. [START_REF] Tilotta | Construction and analysis of a head CT-scan database for craniofacial reconstruction[END_REF]. From the pre-segmentation step, two sets of binary images -respectively for bone and soft tissues-are obtained. By stacking these slices one can essentially detect the 3d structures. However the intensity-based segmentation is not enough to ensure a correct separation of the tissues of interest, due to the presence of noise on the data and artefacts occurring during the acquisition process. The binary images are then cleaned by removing the so-called islands (very small structures whose contours are defined by only a few pixels). These structures can be external to the tissues of interest (noise added during the acquisition process) or very thin internal structures (small bones inside the cranium). This action helps denoising the images. Moreover in most of the subjects, large artefacts due to dental filling are observed on the images. These defects need to be manually removed on each affected slice (see Figure 2).

Geometric mesh processing

From the segmented 2d slices, 3d mesh models of the skull and the face are generated by a marching cube algorithm, Stalling et al. [START_REF] Stalling | Amira: a Highly Interactive System for Visual Data Analysis[END_REF]. These initial meshes contain, in general, a prohibitive number of elements that are redundant and oversampled to correctly describe the geometry of the model. These dense meshes are simplified thanks to the surface remeshing tool mmgs, developed by Dapogny et al. [START_REF] Dapogny | Threedimensional adaptive domain remeshing, implicit domain meshing and applications to free and moving boundary problems[END_REF] by a local modification approach. The remeshing procedure aims at providing (i) a correct and accurate geometric approximation of the underlying 3d model (geometric mesh) and (ii) a computational mesh of high quality elements suitable for finite elements simulations (computational mesh). To ensure vicinity between the original and the remeshed triangulations, the remesher controls the discrete Hausdorff distance between the two sets of triangles. The element size of the simplified mesh is locally adapted to the surface curvatures, en-suring a correct approximation of the surface geometry. As illustrated in Figures 3 and4, this procedure also removes the 'staircase' artefacts due to the spatial discretization and connectivity.

Generation of surface model of the skull

The human skull is characterized by a complex structure, showing small details that are difficult to both acquire and handle numerically. Due to these difficulties, several authors opt for describing the skull through underlying anatomical or geometrical substructures. The most popular choice leads to the definition of a sparse set of anatomical landmarks, possibly coupled with a dense representation of the skull, see Claes et al. [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: Methodology and validation[END_REF], Vandermeulen et al. [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF] or Knyaz et al. [START_REF] Knyaz | Virtual face reconstruction based on 3D skull model[END_REF]. Some authors use automatically detected continuous crest-lines, Quatrehomme et al. [START_REF] Quatrehomme | A Fully Three-Dimensional Method for Facial Reconstruction Based on Deformable Models[END_REF]. The process of matching skull templates is then driven by the outlined feature structures, by requiring their best alignment.

The purpose of this section is to characterize the skull template in terms of a bounded domain, known by a closed mesh of its boundary. This issue is related to the more general problem of surface reconstruction from sample points. Surface reconstruction methods have been extensively investigated in the context of interface evolution via level set methods (see Claisse and Frey [START_REF] Claisse | A nonlinear PDE model for reconstructing a regular surface from sampled data using a level set formulation on triangular meshes[END_REF] and references therein) or via deformable surfaces, as in Duan et al. [START_REF] Duan | Shape Reconstruction from 3D and 2D Data Using PDE-Based Deformable Surfaces[END_REF] or Miller et al. [START_REF] Miller | Geometrically Deformed Models: A Method for Extracting Closed Geometric Models Form Volume Data[END_REF]. The specific task, already evoked in Jones [START_REF] Jones | Facial Reconstruction Using Volumetric Data[END_REF], Tilotta et al. [START_REF] Tilotta | A local technique based on vectorized surfaces for craniofacial reconstruction[END_REF], is to define a closed surface of the skull. In order to achieve this goal an original wrapping algorithm based on mesh deformation techniques is proposed.

First, a (possibly invalid) source triangulation T S of the skull geometry is generated by standard segmentation and 3d reconstruction tools (see Section 1). Then, a closed surface, for example a sphere, is iteratively deformed, producing a sequence of surfaces that are 'closer and closer' to the source triangulation T S . The final deformed surface partitions the ambient space into two sub-regions defining unambiguously an interior (resp. exterior) domain. This preliminary step is of crucial importance for matching shapes among the head database (see Section 3).

Presentation of the method

Let Γ 0 (resp. Γ 0 ) be the sphere of radius r > 0 (resp. R > r) and let Ω 0 be the volume between Γ 0 and Γ 0 . Suppose that Ω 0 encircles T S and is filled with a linear elastic material whose deformation is driven by the Elasticity equations recalled in the box below.

Elasticity equations

The displacement field of a shape Ω clamped at a part Γ D of its boundary Γ and deflated under the effect of internal pressure p ∈ H -1/2 (Γ N ) on 3 , where the stress tensor σ follows the Hooke's law:

Γ N = Γ \ Γ D is achieved as the unique so- lution u ∈ H 1 Γ D (Ω) 3 := {w ∈ H 1 (Ω) 3 , w = 0 on Γ D } of the following variational prob- lem: Ω σ(u) : ε(v) dx = Γ N pv • n ds, (2.1)
for all v ∈ H 1 Γ D (Ω)
σ(u) = 2µε(u) + λ tr(ε(u))I,
with λ, µ the Lamé coefficients and ε(u) = 1 2 (∇u + ∇u T ) the linearized strain tensor.

Hence, starting from Ω 0 deflated on its interior boundary Γ 0 , a sequence of shapes Ω k -with interior boundaries Γ k -is produced. Points of Ω k are advected according to Equation (2.1) until they intercept the triangulation T S . Doing so, at each step k the intersection between the vector displacement and the triangulation T S is checked. In case of multiple intersection points, the closest point is retained. Whenever a contact between the advecting shape and the source triangulation occurs, the first is clamped and forced not to cross the boundary of the latter.

The advecting sequence of internal boundaries Γ k gets closer and closer to the source mesh T S thanks to the strict inequality:

d H (Γ k+1 , T S ) < d H (Γ k , T S ),
where d H (•, •) is the Hausdorff distance between the two surfaces. Finally, for k sufficiently large, Γ k defines a closed boundary that wraps the source triangulation T S .

Numerical issues

From the numerical point of view, the solution of the Elasticity equations (2.1) is computed by the Finite Elements Method on a volume mesh, that is Ω 0 is filled with tetrahedra. For its part the boundary Γ 0 is discretized as a surface mesh T 0 . The iterative algorithm is performed in order to get a sequence T k of meshes with decreasing values of d H (T k , T S ). Given a fixed integer N , the procedure ends if all the points in T k have reached T S or if N iterations of the process occur without registering a new intersection. The latter condition deals with the potential presence of holes in T S .

Figures 5 and6 depict the wrapping of a skull triangulation. The parameters used for Equation (2.1) correspond to a very soft and compressible material. The procedure ends after 35 iterations with N = 5, running in a few minutes on a standard laptop computer. About 85% of points on the triangulation T 35 are clamped on the surface mesh T S .

Inter-subject shape matching

Shape morphing or matching arises in a wide variety of situations in areas ranging from biomedical engineering to computer graphics. Beyond the specific stakes to each particular application, the general issue is to find one transformation from a given 'template' shape Ω 0 into a 'target' Ω T (see Figure 7). Such a transformation may be used as a means to appraise how much Ω 0 and Ω T differ from one another -for instance in shape retrieval, classification or recognition -or to achieve physically the transformation from Ω 0 to Ω T (in shape registration, reconstruction, or shape simplification). See for instance Veltkamp [START_REF] Veltkamp | Shape matching: Similarity measures and algorithms[END_REF] and references therein for an overview of several related applications. In the facial reconstruction context, shape matching is a key ingredient for studying the shape database, making The problem is stated as follows: given a 'template' shape Ω 0 and a 'target' shape Ω T , the aim is to deforme the mesh of Ω 0 onto a computational mesh of Ω T .

To achieve this purpose, a method that has much in common with that of Bajcsy and Kovacic [START_REF] Bajcsy | Multiresolution elastic matching[END_REF] is proposed, borrowing techniques from shape optimization, and more generally optimal control. Under the assumption that Ω 0 and Ω T share the same topology, the desired transformation from Ω 0 to Ω T is obtained through a sequence of elastic displacements, which are obtained by minimizing an energy functional based on the distance between Ω 0 and Ω T . In this section, the mathematical framework and some numerical issues related to the method are briefly presented. Refer to *** [5] for more details and for 2d and 3d numerical examples.

Presentation of the method

The discrepancy between a reference shape Ω and a target shape Ω T is measured by the following functional J(Ω) of the domain:

J(Ω) = Ω d Ω T (x)dx, (3.1)
which involves the Euclidean signed distance function d Ω T to Ω T , defined as:

d Ω T (x) =    -d(x, ∂Ω T ) if x ∈ Ω T , 0 if x ∈ ∂Ω T , d(x, ∂Ω T ) if x ∈ c Ω T .
In the above formula, d(•, ∂Ω T ) denotes the usual Euclidean distance function to the boundary ∂Ω T .

In order to decrease the value of J(Ω), the domain Ω must expand in the regions of the ambient space where d Ω T is negative (that is, in the regions comprised in Ω T ), and to retract in those where it is positive. Note that the functional J(Ω) has a unique, global minimizer Ω = Ω T , and no extra local minimum point provided Ω T is connected. It is then expected that an iterative algorithm for minimizing J(Ω), starting from Ω 0 , will lead to an interesting way to transform Ω 0 into Ω T . This paves the way for an iterative algorithm, producing a sequence (Ω k ) k=0,... of shapes, which are 'closer and closer' to Ω T : at each step, Ω k is updated according to

Ω k+1 = (I + u)(Ω k ), (3.2)
where u is a suitable descent direction for J(Ω). Now, imagine that all the considered shapes Ω are filled with a linear elastic material. One can compute the unique solution u of the elasticity equations (2.1) where the pressure p is taken equal to -d Ω T . This vector field u is naturally a descent direction for J(Ω) since, by a classical calculation (see *** [5] for more details), the shape derivative of the function J(Ω) satisfies J (Ω)(u) ≤ 0.

Numerical issues

For numerical implementation, on the one hand the template shape Ω 0 is filled with tetrahedra and the Elasticity equations are solved on this volume mesh by the Finite Element Method. On the other hand the target shape Ω T is only supplied through its signed distance function, e.g. as a piecewise affine function on a large computational box.

Starting from the template shape Ω 0 , a gradient descent algorithm with adaptive step size is performed in order to get a sequence (Ω k ) k=0,... of domains with decreasing values of J(Ω k ). The algorithm stops when the step size is smaller than a fixed tolerance. The global mapping from Ω 0 to Ω T is easily recovered by the composition of the local displacements between each iteration.

In the proposed examples, the calculation of the signed distance function to Ω T is performed using the algorithm described in Dapogny and Frey [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. The error is computed as the L 2 -norm of the distance d Ω T calculated on the boundary of the final shape.

Figure 8 depicts the matching of two skull shapes. The convergence of the gradient descent procedure is obtained in 300 iterations and the final error equals 0.1 mm (much smaller than the minimal mesh size), revealing an excellent matching of Ω 300 with Ω T . Next, a face example is considered; see Figure 9. 400 iterations of the gradient descent algorithm have been performed to achieve convergence, running in a few minutes on a standard laptop computer. The error calculated on the boundary of the final shape Ω 400 is 0.2 mm (again, much smaller than the minimal mesh size).

Remark 3.1. The ears are not included in the face-matching process, since these structures are not linked with the underlying skull morphology. This procedure is simply done by defining patches on the face template and including only the interesting patches in the minimization of J(Ω).

Registration of the craniofacial template onto the unknown skull

The parametrization of skull, face and soft tissue information depends strongly on the nature of the database. Using CT scans of living subjects allows a dense representation of the craniofacial information. If closed surface models of both skull and face are available (see Section 2), then the craniofacial template may be described in terms of the 3d domain delimitated by the outer face and the inner skull boundaries, generating a 3d 'mask' that incorporates the soft tissue information. A single craniofacial template may be used for producing templates of all individuals in the database. The procedure relies on the shape matching technique (see Section 3), generating a set of vector fields, one for each individual, each of them allowing to deform the template shape onto a specific shape of the database.

Various methods have been proposed for deforming the craniofacial template onto the unknown skull. On the one hand, template deformation methods use generic mathematical transformations such as radial basis functions (Quatrehomme et al. [START_REF] Quatrehomme | A Fully Three-Dimensional Method for Facial Reconstruction Based on Deformable Models[END_REF]). Care has to be taken because no knowledge of facial anatomy is incorporated in these transformations leading sometimes to awkward looking faces. On the other hand face-specific methods (Claes et al. [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: Methodology and validation[END_REF], Paysan et al. [START_REF] Paysan | Face reconstruction from skull shapes and physical attributes[END_REF], Berar et al. [START_REF] Berar | 3d semi-landmark-based statistical face reconstruction[END_REF], Bai et al. [START_REF] Bai | Face Reconstruction from Skull Based on Least Squares Canonical Dependency Analysis[END_REF], Duan et al. [START_REF] Duan | 3D face reconstruction from skull by regression modeling in shape parameter spaces[END_REF]) use Principal Component Analysis on a facial database to generate a statistical shape model. The advantage is that faces are then deformed in a face-specific manner within statistical boundaries. The disadvantage is that the deformations are dependent on the database and can be too restrictive in case of small databases.

Here the method lies on the 'physical' deformation of the craniofacial template onto the unknown skull. This approach allows to overcome the drawbacks which affect existing template deformation methods. The mask is deformed according to a displacement field prescribed on the skull boundary, measuring the deformation of the skull template onto the unknown skull. Under the effect of boundary changes, the mask is allowed to deform as an elastic material, result- ing in a deformed face shape that is now adapted to the unknown skull. The use of elasticity for soft tissue deformation is largely customized in the field of computational maxillofacial surgery, Chabanas et al. [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF], Zachow et al. [START_REF] Zachow | Finite-Element Simulation for Soft Tissue Prediction[END_REF].

Craniofacial template generation

Let (S i ) i=0,••• ,N and (F i ) i=0,••• ,N be the collection of skulls and faces of the database. Thanks to the procedure described in Section 2, theses shapes can be at this point considered as closed surfaces. Then, the volume between S i and F i can be defined, generating a 3d mask M i (called craniofacial template). The craniofacial templates are first roughly aligned by an Iterative Closest Point algorithm. Then, a reference craniofacial template, say M 0 , is chosen and filled with an elastic material (see Figure 11). It is matched onto

all the M 1 , • • • , M N (with global displacements u 1 , • • • , u N )
thanks to the shape matching process described in Section 3, that is:

M i = (I + u i )(M 0 ) for all i = 1, • • • , N.
See Figure 13.

Note that this process generates meshes for all the craniofacial templates M 1 , • • • , M N that share the same number of elements (vertices, tetrahedra) and the same connectivity. The computation of average shapes therefore turns into trivial average of the vector fields. Any convex combination w of the form:

w = N i=1 α i u i , with α i ∈ [0, 1] and N i=1 α i = 1,
defines a new shape W through the mapping:

W = (I + w)(M 0 ).
See Figure 12 for an example of generation of new shapes as a convex combination of three templates.

Face-to-skull mapping via skull matching

Now let u 0 be the global displacement mapping the reference skull S 0 onto the target skull S T . The 3d mask M i is elastically deformed under the effect of the boundary changes (see Figure 14). The deformation of the template M i is achieved as the solution U i of an elastic problem similar to (2.1), in which the displacement condition U i = -u i + u 0 is imposed on the skull boundary and zero traction conditions on the face boundary. The soft tissue deformation step produces N new faces, each one of them obtained by linking a specific craniofacial template M i with the unknown skull. These facial items are then combined together to predict the most plausible face associated to the unknown skull.

Results and discussion

Use of facial reconstruction techniques in Forensic Science

An essential deal of the work of a forensic scientist regards the identification of post-mortem bodies. In the special case in which the facial remains are no longer available, the forensic scientist can resort to various techniques to achieve a positive identification of the dead: DNA assessment, dental record analysis, the identification of a tattoo, the serial number of a prosthesis, and other accepted methods of identification. All these techniques are based on a principle of comparison (with respect to a sample or to a database) and can indeed fail in absence of comparative elements. For instance, at the Forensic Medicine Institut in Paris, this situation shows up for 2 up to 5 cases per year. In this context, a facial reconstruction is then the last resort for reaching a positive identification. The key idea is to scan the skeletal remain to generate a virtual replica of the dry skull and to use a facial reconstruction method to produce a virtual estimation of the victim's face. Then, the produced output is compared with a pool of facial photos which contains the actual target and other individuals of similar age and ethnicity. Then, a group of volunteers is asked to choose which face of the pool is the most resemblant to the reconstruction (face pool assessment). In the context of a murder investigation, a call for witnesses can also be demanded. Indeed the reconstructed face should sufficiently resemble to the original victim's face to allow a recognition and enhance a positive identification. For examples of forensic cases which have employed facial reconstruction techniques see the book [START_REF] Prag | Making Faces: Using Forensic and Archaeological Evidence[END_REF] or the reviews [START_REF] Claes | Computerized craniofacial reconstruction: Conceptual framework and review[END_REF][START_REF] Wilkinson | Facial reconstruction: anatomical art or artistic anatomy ?[END_REF]. Actually in France, forensic facial reconstructions are produced using 2d reconstruction techniques and are rarely successful. Since a few years our multidisciplinary team works on the conception and development of a open-source software, based on the 3d deformation techniques presented in this paper, which will be used during criminal investigations in addition to other endorsed methods. The main concern of this project is to develop a software fully-automatic, usable even by non trained users, and allowing to produce with short time consuming a facial output (or a range of possible outputs) starting from the scan of the unknown skull.

Validation

The validation of a method for facial reconstruction is of tremendous importance to legitimate its use during a criminal investigation. The natural way to address this problem is the "leave-oneout scheme". One individual is removed from the database and the method is employed to reconstruct his/her face given the sole skull. Then the predicted face is compared with the available "unknown face". The discrepancy between the two shapes can be evaluated mathematically by computing distances between them (see below). However, since the final purpose of the method is a positive identification, a recognition test can be also employed for revealing the power of prediction of the method, as it is done in Claes et al. [START_REF] Claes | Computerized craniofacial reconstruction: Conceptual framework and review[END_REF]. A recognition test consists in showing the predicted face together with a sample of faces that contains the unknown face. Then the human volunteer indicates the face (or faces) that is (are) closest to the predicted one. The positive outcome will then correspond to the identification of the unknown face among the sample. Such a study for the method proposed in this article is under process and results will be published in a forthcoming paper.

Hereafter three examples of reconstruction are presented. According to the leave-one-out scheme, one individual is removed from the database and a prediction of his face is generated. Figures 15,16 and 17 show the results obtained by computing the mean of all the faces transported onto the unknown skull. Figure 18 shows the distribution of the error (signed distance function) over the surface for the three reconstructions. Remark that in the regions of interest, excluding ears, neck and eyes in particular, the error is less than 1 millimeter for most points. Also notice that in these regions, the thickness of the tissue is underestimated or overestimated. Using the information coming from the BMI of the individuals in the database could be a way to propose different reconstructed faces corresponding to different BMI of the unknown individual.

The discrepancy D between a surface Γ and a reference surface Γ 0 can also be evaluated by the following mean error:

D = 1 |Γ| Γ d 2 (x, Γ 0 )dx 1 2
, where |Γ| is the measure of Γ and d(•, Γ 0 ) is the Euclidean distance to Γ 0 . This error estimator is used for evaluating the vicinity between:

1. the skull templates of the database and the unknown skull;

2. the face templates of the database and the unknown face associated with the unknown skull;

3. the face templates after deformation onto the unknown skull and the unknown face associated with the unknown skull;

4. the face templates of the database and the reconstructed average face.

The minimal, maximal and mean values of D for the 25 templates of the database are reported in Table 1 in each of the third cases. In particular, the discrepancy between a face template and the unknown face is smaller after deformation for all the individuals. Thus, the elastic transformation used in the method seems to be a good tool to transport the faces close to the unknown face. Moreover, the discrepancy between the unknown face and the predicted face is smaller than the discrepancy with any individual in the database, so this measure D can be used for an automatic numerical identification.

Comparison with previous works

A comparison with previous facial reconstruction methods is a quite hard task. The reliability of two methods should be compared starting from the same database and the same test cases, to avoid the dependencies of the methods on the data. Such a comparison is out of the scope of this paper. Here we limit ourselves to compare the results of the test cases with the outcomes of some statistical facial reconstruction methods referenced in the literature [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: Methodology and validation[END_REF][START_REF] Tilotta | A local technique based on vectorized surfaces for craniofacial reconstruction[END_REF][START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF]. In Claes et al. [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: Methodology and validation[END_REF] a statistical method is developed for the purpose of estimating a facial outcome using landmarks. The experiment is carried out from a database of 118 individuals. The mean reconstruction error registered equals 1.4 millimeters, with a standard deviation of 1.04 millimeters. In Vandermeulen et al. [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF] a statistical method based on implicit surfaces was employed for the facial reconstruction test using a database of 20 individuals. The mean error registered equals 1.9 mm with a standard deviation of 1.7 mm. In Tilotta et al. [START_REF] Tilotta | A local technique based on vectorized surfaces for craniofacial reconstruction[END_REF] a local technique based on surface patches is used for statistically predicting some facial regions. The experiment is carried out using a database of 47 individuals. For both chin region and nasal region, the mean reconstruction equals 0.99 mm, ranging from 0.58 millimeters to 1.83 millimeters for the nose, and from 0.21 millimeters to 2.41 millimeters for the chin. The mean error registered with the method proposed in this paper is reported in Table 1.

Note that most of the error is concentrated on the ears, the neck, the eyes and lips regions. The distribution map of Figure 18 shows that for the three studied test cases the mean error over the remaining regions is smaller than 1 millimeter. The result is similar to the one reported in Tilotta et al. [START_REF] Tilotta | A local technique based on vectorized surfaces for craniofacial reconstruction[END_REF] (which only deals with the reconstruction of two facial regions) and improves the outcomes of Claes et al. [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: Methodology and validation[END_REF] and Vandermeulen et al. [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF], at least when considering the regions of interest. Also, note that all the experiments were conducted on a smaller database and with no need of landmarks, resulting in a significant improvement in the fastness and feasibleness of the overall process, and avoiding manual intervention.

Conclusion

The proposed reconstruction method lies on the 'physical' deformation of templates of coupled faces and skulls onto the unknown target skull.

In practice, the acquisition of full head scans of healthy subjects is still a difficult process. Most of the time one can have access to clinical data, meaning that the patients present morphological anomalies or that the scans are only partial (in the case of maxillo-facial examination for surgery purposes). The access to an adequate database of full heads of healthy subjects would enormously improve the final product of the method. Even if the experiments were carried out on a

Figure 1 -

 1 Figure 1 -From CT scans to densely matched surface templates (see Sections 1, 2 and 3).

Figure 2 -

 2 Figure 2 -Left: reference gray-level image. Middle: pre-segmentation of the bone using a threshold technique (red). Right: segmentation of the bone after denoising and artefact removal (red).

Figure 3 -

 3 Figure 3 -Left: 3d mesh of the skull after segmentation (2 179 332 triangles). Right: 3d mesh of the skull after geometric remeshing (15 6301 triangles).

Figure 4 -

 4 Figure 4 -Left: 3d mesh of the face after segmentation (1 286 542 triangles). Right: 3d mesh of the face after geometric remeshing (90 949 triangles).

Figure 5 -

 5 Figure 5 -Left: initial surface Γ 0 . Middle: advecting surface Γ 10 . Right: final surface Γ 35 .

Figure 6 - 8 Figure 7 -

 687 Figure 6 -Top: 3d model of the skull. Bottom: wrapped skull.

Figure 8 -

 8 Figure 8 -(a) Template shape Ω 0 and discrepancy with respect to the target shape. (b) Deformed shape Ω k for k = 300 and discrepancy w.r.t. the target shape. (c) Target shape.

Figure 9 -

 9 Figure 9 -(a) Template shape Ω 0 and discrepancy with respect to the target shape. (b) Deformed shape Ω k for k = 400 and discrepancy w.r.t. the target shape. (c) Target shape.

Figure 10 -

 10 Figure 10 -(a) Target shape Ω T as the zero level set of the signed distance function; (b) Adaptive remeshing of the computational box;

Figure 11 -

 11 Figure 11 -Generation of a 3d mask encoding soft tissue information.

Figure 12 - 3 Figure 13 -

 12313 Figure 12 -The corners of the triangle display three faces in the database. The top template is chosen as reference shape and matched onto the bottom right template and the bottom left template. The middle of the triangle and the mid-edges display the 'barycenter' faces computed by averaging the displacement fields.

Figure 14 -

 14 Figure 14 -Elastic deformation of the mask under the effect of skull changes. 14

Figure 15 -

 15 Figure 15 -Test case 1: (a) Original face. (b) Reconstructed face.

Figure 16 -

 16 Figure 16 -Test case 2: (a) Original face. (b) Reconstructed face.

Figure 17 -

 17 Figure 17 -Test case 3: (a) Original face. (b) Reconstructed face.

Figure 18 -

 18 Figure 18 -Discrepancy between the reconstructed and the original faces (a) for the test case 1 of Figure 15; (b) for the test case 2 of Figure 16; (c) for the test case 3 of Figure 17.

Acknowledgements

***

This work received a financial support by Idex Sorbonne Universités under the French funds "Investissement d'Avenir", reference ANR-11-IDEX-0004-02

small collection of 26 individuals, the preliminary results produced are very promising. The proposed method for shape matching allows an accurate registration. The method is simple to implement and does not require any a priori landmark correspondence, allowing for an automatic processing of the database.