Keywords: Manual segmentation, semi-automatic segmentation, Masseter muscle

La Segmentation semi-automatique

On appelle segmentation, l'identification sur les images des voxels (pixels en deux dimensions) appartenant à un même objet d'intérêt. Dans notre cas trois objets : le crâne, le visage et le muscle masséter.

La segmentation du crâne et du visage peut être effectuée rapidement à l'aide des niveaux de gris. En effet, il est facile de réaliser automatiquement la segmentation des éléments qui ont un fort -La gestion d'une base réduite de formes en 3D totalement automatique, -Une nouvelle méthode anatomique de reconstruction faciale, -Un outil de validation et de comparaison des méthodes de reconstruction, -Une méthode de segmentation semi-automatique.

Résumé et état de l'art

Ce manuscrit se divise en 3 parties distinctes, chacune composée de 2 à 3 articles soumis ou en cours de soumission. Ci-dessous se trouve un descriptif détaillé de chacun de ces articles, ainsi que l'état de l'art et les résultats qui lui sont associés.

Partie I -Travail préliminaire

Cette première partie regroupe un travail préliminaire sur les relations anatomiques entre le crâne et le muscle masséter, qui nous ont permis de faire l'hypothèse selon laquelle il est possible de déduire la forme du masséter uniquement à partir de celle du crâne. Et le travail sur la gestion d'une base réduite de formes 3D, qui nous permet effectivement de faire les reconstructions du masséter.

Corrélation entre céphalométrie et muscle masseter

Ce travail a été effectué en collaboration avec Adeline Kerbrat, Patrick Goudot et Thomas Schouman.

Bibliographie

, on peut en visualiser plusieurs exemples fréquemment utilisés. FIGURE 13 -Exemple de différents landmarks céphalométriques Les études céphalométriques des dernières décennies ont permis de classifier différents phénotypes [195] de longueur faciale [75] ou de classes dentaires [7, 107], qui ont aidé les cliniciens à Application pratique de cette méthode semi-automatique Ce travail a été effectué en collaboration avec Camille Cardot, Damiano Lombardi, Thomas Schouman et Wided Siai.

Le propos de cet article est de comparer l'utilisation de notre méthode semi-automatique, exposée ci-dessus, à la méthode de référence manuelle actuelle, par des praticiens. Il s'agit donc entre autres d'une remise en contexte pratique de cette méthode.

). FIGURE 26 -Les deux logiciels de segmentation utilisés : à droite, Amira, et à gauche, la méthode semi-automatique -A new anatomical method of facial reconstruction, -A tool for validating and comparing reconstruction methods, -A semi-automatic segmentation method.

Summary and state of the art

This manuscript is divided into 3 distinct parts, each composed of 2 to 3 articles submitted or in the process of being submitted. Below is a detailed description of each of these articles, as well as the state of the art and the results associated with it.

Part I -Preliminary Work

This first part includes preliminary work on the anatomical relationships between the skull and the masseter muscle, which allowed us to make the hypothesis that it is possible to deduce the shape of the masseter only according to that of the skull. And the work on the management of a reduced 3D shapse base, which effectively allows us to propose the reconstructions of the masseter.

Correlation between cephalometry and muscle masseter

This work was conducted in collaboration with Adeline Kerbrat, Patrick Goudot and Thomas Schouman.

Introduction

Materiels and Methods

anthropologiques connus pour leur constance évolutive, nommés landmarks [30]. Ces mesures peuvent par exemple caractériser des épaisseurs de tissus mous, ou d'autres marqueurs décrivant la variabilité de ces tissus [START_REF] Peter | Computerized craniofacial reconstruction : conceptual framework and review[END_REF][START_REF]Computerized forensic facial reconstruction[END_REF]. Toutefois bien que l'identification de ces landmarks soit rendue possible par un expert, leurs propriétés géométriques sont floues et ne permettent pas de les identifier automatiquement via un ordinateur.

Le projet FaciLe se démarque de ces méthodes par deux aspects :

en basant son approche de la reconstruction faciale tridimensionnelle non pas sur les statistiques mais sur une méthode déterministe, utilisant des équations aux dérivées partielles (EDP), de l'analyse numérique et de la modélisation géométrique 3D. -en rendant l'entièreté du processus complètement automatique et sans intervention humaine, après acquisition des données du crâne.

Objectif principal de cette thèse

La thématique de la reconstruction faciale est le « fil rouge » de cette thèse, qui décrit l'ensemble de la réflexion scientifique : de l'acquisition et du traitement des données, aux résultats en sciences légales. L'objectif est de développer une méthode complémentaire à celle développée dans le cadre de la thèse de Chiara Nardoni [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], valorisant l'ensemble du projet, tout en apportant des précisions anatomiques dans le processus. Différents problèmes relatifs à cette nouvelle méthodologie, ont été soulevés, occasionnant durant ce travail des études parallèles pour le compléter. Il s'agit par exemple de l'automatisation complète du processus, ou encore la segmentation de tissus mous dans des images médicales.

Démarche scientifique : une application des mathématiques

Les mathématiques appliquées, branche des mathématiques, s'intéressent à l'étude et l'application de modèles à d'autres disciplines, telle que la physique, la chimie, la biologie, l'économie, ou encore la médecine. Elles regroupent toutefois un grand nombre de sous-domaines, en partie indépendants des applications elles-mêmes, comme l'analyse numérique, le contrôle, l'optimisation, la théorie des jeux, ou encore les probabilités.

C'est dans cette vaste discipline que s'inscrit cette thèse de mathématiques appliquées dédiée à la reconstruction faciale tridimensionnelle pour la médecine légale. Différents aspects du problème nous amèneront à utiliser des outils mathématiques variés, tels que les équations différentielles, la réduction de dimensionnalité, les éléments finis... FIGURE 3 -Reconstruction faciale tridimentionnelle pour un couple de forme 'départ' et 'cible' W 0 et W T dans R d (d = 2, 3) donné, à trouver un moyen de transformer W 0 en W T . Ce type de méthode joue un rôle important dans des applications en modélisation géométrique et en analyse de formes [START_REF] Remco | Shape matching : Similarity measures and algorithms[END_REF]. Beaucoup de chercheurs des communautés mathématique et infographique s'y sont intéressés [START_REF] Ruzena | Multiresolution elastic matching[END_REF][START_REF] Mirza | Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms[END_REF][START_REF] Serge | Shape matching and object recognition using shape contexts[END_REF]67,[START_REF] Bruno | A numerical algorithm for L2 semi-discrete optimal transport in 3D[END_REF][START_REF] Konrad | A linear elastic force optimization model for shape matching[END_REF]. Une telle déformation peut servir à mesurer à quel point les deux formes diffèrent ou aussi à simplement accomplir la transformation d'un objet en un autre. Ici nous mesurons la différence entre la forme de 'départ' et celle 'cible' selon la fonctionnelle de forme suivante :

J(W) = Z W d W T (x)dx,
impliquant la fonction distance signée d W T de ∂ W T . La fonction distance signée est une fonction « level set » d'un domaine : qui pour chaque point x de l'espace retourne la distance Euclidienne d(x, ∂ W T ) au bord ∂ W T , le signe associé à cette distance est positif si le point se situe à l'intérieur du domaine, négatif si non. Le « mapping » global entre W 0 et W T est obtenu par un procédé itératif qui minimise J(E (W 0 )) (avec E la deformation élastique), inférant à chaque étape une direction de descente, via un algorithme de gradient, pour la fonctionnelle J. Les conditions optimales associées à ce problème d'optimisation de forme sont calculées en utilisant le cadre de la dérivée de forme d'Hadamard [3,[START_REF] Jean | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]96]. En pratique, l'espace admissible des déformations est restreint par une paramétrisation élastique du déplacement. Les équations de l'élasticité ont été implémentées avec la méthode des éléments finis. FIGURE 6 -Exemple 2D simple de morphing Cette première phase du travail a été principalement portée par la thèse de Chiara Nardoni [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], encadrée par Maya de Buhan et Pascal Frey. Nous avons effectué certains de ces travaux ensemble : particulièrement le morphing dit surfacique et sa comparaison avec celui dit volumique (détaillé dans le Chapitre 6). A l'issue de cette thèse dont j'ai repris le travail sur la reconstruction faciale, j'ai développé une méthode complémentaire présentée ci-dessous, et détaillée dans le Chapitre 7, se concentrant sur l'objectif principal de cette thèse, la reconstruction faciale.

Anatomie de la face

Toutefois, cette méthode élastique ne tient pas compte pas de l'anatomie sous jacente du visage, et donc, plus précisément, des composantes fonctionnelles du masque de tissus mous évoqué ci-dessus. Des articles comme celui de [START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF], nous rappellent pourtant l'importance des connaissances anatomiques, accumulées dans diverses disciplines, sur le visage et dont le domaine de la reconstruction faciale ne prend que très peu en compte et de manière non systématique. Nous nous sommes donc tournés vers nos collègues anthropologues et chirurgiens, participant au projet, afin d'améliorer ce premier modèle, en prenant mieux en compte l'anatomie du visage.

Introduction -En français

Le muscle masséter

Le visage est principalement, en terme de volume, constitué de muscles et de graisse, recouverts de peau. La graisse étant une composante dont nous ne retrouvons pas d'information sur le crâne, il a été décidé de la considérer à la fin du processus comme une variable à ajuster une fois la forme du visage identifiée. On peut en revanche facilement identifier jusqu'à 100 muscles sur un visage. Deux d'entre eux, plus volumineux que les autres, ont attiré notre attention : le muscle masséter (Figure 7 Du fait de leur importance fonctionnelle et morphologique respective, nous avons décidé de nous concentrer plus particulièrement sur le masséter du fait de son impact sur la forme du visage. Comme on peut le voir sur la Figure 7, le muscle masséter se situe au niveau de la joue, en avant de l'oreille, le long de la branche de la mandibule. Alors que le muscle temporal, qui se situe lui au niveau des tempes, au dessus de l'arcade zygomatique, est souvent plus fin et couvert par les cheveux, ayant donc dans notre problème un impact moins important dans un premier temps.

Base de données et réduction de dimensionnalité

Il nous faut alors créer un modèle de ce muscle masséter. Afin de le décrire au mieux, aussi simplement que possible, nous utilisons une méthode de réduction de complexité, l'analyse en composantes principales. Décrire la forme du masséter ne nécessite alors plus les 5000 coordonnées 3D des points qui constituaient son maillage, mais seulement n = 5 à 20 coefficients bien choisis, selon le niveau de précision souhaité (Figure 8) : Une fois la reconstruction du masséter (et à terme, du temporal) effectuée, la méthode de déformation élastique exposée dans la section précédente est appliquée, cette fois-ci non pas à partir d'un masque crâne/visage mais crâne+muscles/visage. Cela intégrant des informations plus précises dans ce procédé, tenant compte de l'anatomie, et mettant l'accent sur le lien entre forme et fonction, idée fondatrice de cette thèse.

Acquisition des données

L'ensemble des travaux présentés ci-dessus a été rendu possible grâce à la mise en place d'une base de données associant crâne, visage et muscles pour un même individu comme mentionné auparavant. Cette base de données est constituée d'images médicales 3D qui permettent l'extraction de toutes les informations nécessaires à l'application de cette méthode. Toutes les imageries scanner utilisées dans ce travail de thèse proviennent du service de Chirurgie Maxillo-Faciale de l'Hopital de La Pitié Salpétrière à Paris. Cette base de données ne recoupe en rien celle utilisée par Chiara Nardoni avant moi, qui avait obtenu ses données d'un autre projet, sans lien académique avec le projet FaciLe, sur la reconstruction faciale menée par Françoise Tilotta et Yves Rozenholc [START_REF] Françoise | Construction and analysis of a head CT-scan database for craniofacial reconstruction[END_REF][START_REF] Tilotta-Yasukawa | Contribution à la reconstitution faciale en médecine légale : Étude préliminaire d'une méthode tridimensionnelle à partir de l'imagerie numérique[END_REF]. Il semblait plus pertinant à cette phase du projet FaciLe de constituer notre propre base de donnée.

Scanner : précision et difficultés

Le choix du type d'imagerie n'est pas anodin. En effet bien que l'IRM présente une meilleure précision en ce qui concerne les tissus mous, la donnée d'entrée de notre problème reste le crâne seul. La meilleure imagerie 3D pour visualiser les os avec précision est, à l'heure actuelle, la tomodensitométrie ou encore appelé scanner, dont on peut voir un exemple Figure 9. Malgré les avantages liés à l'ensemble des informations connues par les images de scanner, il faut garder à l'esprit que de nombreux problèmes intrinsèques aux images altèrent l'acquisition et la visualisation des données en trois dimensions, comme c'est le cas des artefacts créés par le métal utilisé pour soigner les dents (braquettes, couronnes...).

En plus de ces difficultés, certaines contraintes, venant cette fois du modèle, compliquent d'autant plus le pré-traitement des données, par exemple : ne pas tenir compte de la colonne vertébrale dans notre modèle de crâne 3D. Ce problème lié à la segmentation d'images médicales est récurrent dans bien des domaines qui ont besoin d'utiliser comme pré-requis des bases de données 3D d'organes. Une des principales limites est le temps consacré à l'aquisition des données, qui ne sont obtenues qu'en petite quantité, et qui limite d'autant plus la qualité des résultats présentés dans ces différents travaux. C'est dans ce contexte où le peu de données a entravé nos résultats, que nous avons été amenés à développer une application de segmentation semi-automatique pour les tissus mous. A l'aide de la base réduite décrite précédemment, nous connaissons a priori la forme de l'objet à identifier sur les images, dans notre cas le muscle masséter. Il suffit d'identifier les coefficients associés à l'objet d'intérêt, via la résolution d'un système prenant en entrée les coordonnées de quelques points sélectionnés sur l'image x ( j) ⇤ , j = 1, . . . , m., avec m le nombre de points sélectionnés :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m.

Cette méthode est indépendante des niveaux de gris de l'image, et peut donc être appliquée sur tous objets et toutes imageries.

Validation des méthodes développées

L'ensemble des méthodes présentées ont été validées d'un point de vue mathématique, mais pour certifier de leur efficacité pratique dans les disciplines où elles ont pour objectif d'être utilisées, nous avons mis en place des validations expérimentales.

Reconstruction et médecine légale

Pour la reconstruction d'un visage à partir d'un crâne, la façon la plus simple de certifier si la méthode est une réussite, est de vérifier si un individu est en mesure de reconnaître ce visage. C'est pourquoi nous avons mis en place une application de sondage : un utilisateur à qui on présente une photo doit choisir parmi 5 propositions, le visage, décrit par un maillage surfacique, qui lui semble être celui de la personne sur la photo (Figure 11).

Cette opération est réalisée en double aveugle et randomisée pour 4 catégories : la surface 3D du visage segmentée directement, la surface 3D de la reconstruction élastique simple, la surface 3D de la reconstruction avec le masséter segmenté directement et la surface 3D de la reconstruction avec le masséter lui-même reconstruit. Cette application permet donc de comparer les résultats de différentes méthodes sur un même jeu de données. En effet, les données utilisées en reconstruction faciale sont souvent extraites de données médicales qui ne sont pas cessibles à d'autres équipes, limitant les comparaisons possibles sur les méthodes à proprement parler.

Segmentation

Dans le cas de la segmentation semi-automatique, l'application développée dans cet objectif étant l'utilisation par des cliniciens ou des chercheurs, nous avons mis en place une étude pour comparer le temps et l'erreur de segmentation avec cette méthode contre la segmentation manuelle de deux praticiennes.

La segmentation manuelle demande d'être précis sur l'ensemble des coupes de scanner où se situe le masséter, alors que la segmentation semi-automatique ne nécessite d'identifier qu'une cinquantaine de points, tel que présenté dans la Figure 12.

Principales contributions

Les principales contributions de cette thèse ont été : retrouve déjà cette idée dans la littérature avec [START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF], toutefois notre approche a pour but de compléter leurs résultats en s'intéressant à des morphologies plus variées et à des mesures souvent peu mises en valeur.

Dans le cadre de ce travail et parmi les très nombreuses mesures céphalométriques qui existent, nous en avons sélectionné cinq qui permettent de compléter au mieux la littérature déjà existante sur le sujet [START_REF] Ariji | Three-dimensional morphology of the masseter muscle in patients with mandibular prognathism[END_REF][START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF]. Il s'agit des angles SNA, SNB, ANB, FMA et entre le plan Bispinal et Mandibulaire, que l'on peut visualiser sur la Figure 16.

Pour ce qui est des mesures effectuées sur le masséter, la communauté s'accorde à dire que les paramètres les plus à propos pour le décrire sont sa hauteur, sa largeur, son épaisseur et son volume [START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF].

FIGURE [START_REF] Siqi | Multi-scale structured CNN with label consistency for brain MR image segmentation[END_REF] -Les mesures céphalométriques que nous avons utilisées pour notre étude L'étude a été menée sur 51 patients du service de chirurgie maxillo-faciale de l'hôpital de La Pitié Salpêtrière, à Paris. Nous avons identifié plus en détaillé les corrélations présentes entre la forme du crâne et celle du muscle masséter.

Le FMA est dans notre étude le paramètre de mesure angulaire osseuse, qui présente plus de corrélation avec les mesures massétérines. En effet, il est inversement corrélé aux 4 mesures du masséter testées (hauteur, largeur, épaisseur et volume). Ce qui signifie qu'un muscle masséter long, large, épais et volumineux est associé à un phénotype dit de « face courte » avec un angle mandibulaire fermé. Ce résultat est similaire à ceux trouvés dans [START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF][START_REF] Kiliaridis | Masseter muscle thickness measured by ultrasonography and its relation to facial morphology[END_REF].

Toutefois peu d'études abordent ce problème selon la position antéropostérieure du maxillaire, ici nous montrons une corrélation du SNA avec l'épaisseur et le volume du masséter et du SNB avec sa hauteur et son volume, contrairement à [START_REF] Munutsugu Kubota | Maxillofacial morphology and masseter muscle thickness in adults[END_REF].

Du fait de la présence de ces corrélations significatives, ces conclusions nous ont donc permis de poser comme hypothèse, pour notre méthode de reconstruction faciale, qu'il serait possible de déduire à partir de la forme des os du crâne, celle du muscle masséter.

Processus ingéniérique de la constitution d'une base réduite de formes 3D

Ce travail a été effectué en collaboration avec Yvon Maday et Chiara Nardoni.

Il s'agit ici d'un processus complet de traitement d'une base de données de formes 3D. Toutes les méthodes mathématiques utilisées ont déjà été présentées individuellement dans la littérature. Cependant, elles sont ici assemblées et ajustées pour répondre aux caractéristiques de nos applications : reconstruction faciale et segmentation semi-automatique. Il est toutefois tout à fait envisageable d'utiliser le même procédé de traitement et de gestion d'une base réduite de données 3D, dans des objectifs et applications différents tels que les jeux vidéos, la biologie moléculaire, etc.

Les étapes

Plusieurs étapes utilisant des outils mathématiques variés sont nécessaires afin de pouvoir utiliser notre base de données de forme.

1. Acquisition des données : dans notre contexte, il s'agit de la segmentation des images médicales. La seconde partie de ce manuscrit est entièrement consacrée à ce problème épineux. La bibliographie et la présentation des méthodes mises en place sont proposées dans la section suivante. L'ensemble des données de scanner collectées pour créer cette base de formes provient du service de Chirurgie Maxillo-Faciale de l'hôpital de La Pitié Salpêtrière, à Paris. Ils ont été acquis avec deux types de fenêtres : tissus mous et os (Figure 17). [32] comme les « Marching Cubes », introduit par [START_REF] William | Marching cubes : A high resolution 3D surface construction algorithm[END_REF]. Ce dernier permet de créer une triangulation surfacique à partir d'un volume de voxels (pixels en 3D), toutefois ils sont aussi connus pour laisser des artéfacts et des éléments invalides dans le maillage, pour des simulations réalisées aves des éléments finis. Il est donc nécessaire d'y ajouter une étape de remaillage, permettant de retrouver une géométrie correcte et de contrôler la qualité du maillage souhaité. Pour cela nous avons utilisé mmg [START_REF] Charles | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF] Cécile | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF], qui nous permet de contrôler la distance de Hausdorff entre le maillage d'entrée et le maillage de sortie, ainsi que d'autres paramètres comme la longueur maximale et minimale des arêtes des triangles du maillage. Les résultats de l'algorithme de Marching Cubes et celui du remaillage sont donnés dans la Figure 18. Alignement : Souvent sous-estimé, l'alignement cohérent des objets de la base de données entre eux, peut avoir un fort impact sur la suite du processus. Effectivement, si les données sont acquises de plusieurs sources qui n'ont pas les mêmes paramètres d'acquisition, les objets peuvent différer à un seul mouvement rigide et/ou une homothétie près. Il faut donc corriger les translations, rotations et homothéties qui influent sur les données. La première phase consiste à centrer les objets W en l'origine du repère à l'aide du calcul de leur barycentre x G respectif défini par :

x G = 1 |W| Z W xdx avec x 2 W et |W| le volume d'W.
On obtient après discrétisation sur le maillage de l'objet (détaillé dans le chapitre 2, page 49) :

x G = 1 2|W| Â i x 2
i n i ds i avec n la normale extérieure à la surface de l'objet W. Si nécessaire, c'est à ce moment qu'une mise à l'échelle peut être mise en place. Pendant la deuxième phase, nous calculons les axes d'inertie pour pouvoir les aligner sur ceux d'un objet de référence, grâce à la matrice d'inertie I W :

I W = Z W (|x| 2 Id x ⌦x)dx
avec Id la matrice identité, et x ⌦x la matrice x •x T . Ces deux étapes permettent d'identifier tous les objets de la Figure 19 comme identiques en forme, si on aligne les axes d'inertie et que l'on compare les maillages.

FIGURE 19 -Cinq muscles masséter identiques à un mouvement rigide près

Il est important de noter que notre méthode s'attache particulièrement à ne pas sélectionner de points spécifiques à mettre en correspondance, même pour l'alignement, car nous souhaitions un processus totalement automatique. Cela est rendu possible par la méthode décrite précédemment combinée aux étapes 4 et 5 (warping et morphing) décrites ci-dessous. 4. Warping (facultatif) : Cette étape n'est nécessaire que lorsque l'objet dont on veut créer une base réduite de forme n'est pas topologiquement fermé. L'algorithme qui consiste à utiliser les équations de l'élasticité pour déformer un maillage template fermé sur un objet ouvert, a été développé dans [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF]. L'idée principale est présentée dans la Figure 20. 5. Morphing : Le principe de cette étape est de déformer une forme de référence W 0 , dans notre cas une ellipsoïde, vers une forme objectif W k , ici chaque muscle masséter de notre base de données. Cela nous permet de considérer tous les objets de la base comme des déplacements u k d'un même maillage de référence, ce qui rend toute comparaison entre eux bien plus facile à mettre en place. Cela permet aussi de s'affranchir de la mise en correspondance de « landmarks », à identifier sur tous les objets de la base de données. Un chapitre complet de cette thèse est dédié à l'explication détaillée de cette méthode et de sa bibliographie dans la troisième partie.

Bibliographie

La segmentation d'images médicales consiste à identifier et labelliser les voxels (pixels 3D) de l'ensemble de l'image ou d'un objet d'intérêt. Dans notre cas, nous nous sommes intéressés au muscle masséter dans le contexte d'une application spécifique : la reconstruction faciale.

Cependant, les méthodes de segmentation automatiques ou semi-automatiques ont un plus vaste champ d'application. Elles sont par exemple très largement utilisées dans le « diagnostic assisté par ordinateur », qui est devenu un des sujets principaux de recherche en imagerie médicale et en radiologie [START_REF] Kunio | Overview on research and development of computer-aided diagnostic schemes[END_REF][START_REF] Kunio | Current status and future potential of computer-aided diagnosis in medical imaging[END_REF][START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Bradley | Computer-aided detection and diagnosis at the start of the third millennium[END_REF][START_REF] Ronald | Road maps for advancement of radiologic computer-aided detection in the 21st century[END_REF] dans les dernières décennies. Dans les années 80, le concept de diagnostic automatique ou assisté par ordinateur était déjà connu depuis des études datant des années 60 [START_REF] Hc Becker | Digital computer determination of a medical diagnostic index directly from chest X-ray images[END_REF][START_REF] Gwilym | Computer diagnosis of primary bone tumors : A preliminary report[END_REF][START_REF] Phillip | Automated computer analysis of radiographic images[END_REF], mais ces premières tentatives n'étaient pas un réel succès. En effet, la détection automatique d'organes ou de lésions s'est révélée extrèmement difficile.

De nos jours, les techniques se sont de plus en plus nombreuses, en partie grâce au développement technologique des ordinateurs et des méthodes d'acquisition. Certains tissus, dont le contraste avec le tissu voisin est suffisamment important, comme les os, sont facilement détectables automatiquement [START_REF] Yan | A new accurate and precise 3-D segmentation method for skeletal structures in volumetric CT data[END_REF][START_REF] Shafry Mohd | 3D bones segmentation based on CT images visualization[END_REF][START_REF] Jing | Fast segmentation of bone in CT images using 3D adaptive thresholding[END_REF]. Il en va de même pour les organes majeurs du corps humain et les pathologies qui leur sont associées, comme le cerveau et le coeur [START_REF] Siqi | Multi-scale structured CNN with label consistency for brain MR image segmentation[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF][START_REF] Thomas | Handbook of medical imaging[END_REF], qui ont toujours suscité un grand intérêt pour les soignants. Les méthodes d'utilisation communes pour ces types d'applications sont les méthodes dites de seuillage ou de détection des bords. La première consiste à regrouper dans un même label les pixels qui ont un niveau de gris proche ou plus exactement contenu par une ou plusieurs valeurs seuils [START_REF] Ramesh | Thresholding based on histogram approximation[END_REF]. La seconde, qui s'appuie aussi sur la détection des niveaux de gris, cherche les frontières des régions anatomiques en analysant les transitions rapides de ces niveaux, et assemble les voxels homogènes dans un même label [START_REF] Antonio | Evaluation of methods for ridge and valley detection[END_REF]. Ces deux méthodes sont fondées sur la détection des niveaux de gris, c'est ce qui permet, entre autres, aux experts de différencier deux régions anatomiques et de réaliser une segmentation manuelle [START_REF] Dzung | Current methods in medical image segmentation[END_REF]. Ce type de stratégie, mimant la connaissance des experts ont été les premières à être mises en place [START_REF] Kunio | Diagnostic imaging over the last 50 years : research and development in medical imaging science and technology[END_REF][START_REF] Kunio | Computer-aided diagnosis in radiology : potential and pitfalls[END_REF]. Différentes classifications de ces méthodes existent, on peut par exemple citer [START_REF] Daniel | Medical image segmentation : Methods and software[END_REF].

Toutefois, plus récemment, d'autres types de méthodes, fondées au moins en partie sur les statistiques, se sont progressivement développées, comme les atlas guidés, les modèles déformables ou encore l'utilisation massive de l'apprentissage. Elles se différencient entre autres par la quantité de connaissances a priori injectée dans la méthode [START_REF] Carlos | Advances in auto-segmentation[END_REF].

Les atlas statistiques utilisent comme référence une ou plusieurs structures d'intérêt déjà segmentées, appelées atlas, comme connaissance préalable pour de nouvelles segmentations [START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF][START_REF] Torsten | Quo vadis, atlas-based segmentation ?[END_REF]. La segmentation d'une nouvelle image repose alors sur un alignement et une déformation qui permet de définir la transformation optimale entre l'atlas et la nouvelle image. Différents types d'algorithmes existent pour optimiser le déplacement [START_REF] Arish | Auto-segmentation of normal and target structures in head and neck CT images : A feature-driven model-based approach[END_REF][START_REF] Thirion | Image matching as a diffusion process : an analogy with Maxwell's demons[END_REF], ou l'utilisation de l'atlas lui-même [2,[START_REF] Olivier | Efficient selection of the most similar image in a database for critical structures segmentation[END_REF][START_REF] Jinzhong | Atlas ranking and selection for automatic segmentation of the esophagus from CT scans[END_REF]. Ces méthodes de segmentation se sont révélées efficaces lors des grands concours de segmentation, organisés par « the American Association of Physicists in Medicine », ou « the Medical Image Computing and Computer-Assisted Intervention Society ».

Les modèles déformables nécessitent quand à eux plus d'informations a priori. Lorsque ces dernières sont disponibles, les variations caractéristiques de la forme de la structure d'intérêt peuvent être utilisées pour former des modèles statistiques de forme. Cela limite les résultats obtenus des segmentations, à des formes anatomiquement plausibles [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF]. Toutefois, ces méthodes sont aussi moins flexibles pour la raison précédemment citée. Par ailleurs, la taille et le contenu des données d'entraînement limitent aussi la performance des résultats. Cette méthode reste cependant fréquemment utilisée pour les structures anatomiques pelviennes [START_REF] Daniel | Model-based segmentation of medical imagery by matching distributions[END_REF][START_REF] Hans | A 3D statistical shape model of the pelvic bone for segmentation[END_REF][START_REF] Heiko | Automatic Segmentation of the Pelvic Bones from CT Data Based on a Statistical Shape Model[END_REF].

L'apprentissage statistique, appliqué à la segmentation d'images, regroupe à lui tout seul un large panel de méthode en plein essor, telle que le support vecteur machine, les arbres de décisions, le deep learning ou encore les réseaux de neurones [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF]. Ces derniers sont déjà retrouvés dans de multiples applications comme la segmentation de tumeur dans le cerveau, ou dans le foie, etc [START_REF] Jinsa | Lung cancer classification using neural networks for CT images[END_REF][START_REF] Mala | Neural network based texture analysis of CT images for fatty and cirrhosis liver classification[END_REF][START_REF] Malay | Cardiac arrhythmia classification using neural networks with selected features[END_REF][START_REF] Saritha | Classification of MRI brain images using combined wavelet entropy based spider web plots and probabilistic neural network[END_REF]. Les résultats obtenus par ces algorithmes sont prometteurs dans le cadre des applications pour lesquels ils sont développés. Mais il n'est pas toujours évident de les étendre à tous les types d'images et cela nécessite souvent un jeu de données important au vu de la variabilité humaine sous-jacente aux données médicales [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF].

La méthode que nous proposons s'appuie sur ces dernières techniques (modèle déformable, atlas statistique et apprentissage) et emprunte certains outils à chacune d'entre elles, pour les combiner de manière à réaliser notre objectif. Elle est comparable à [START_REF] Florian | Shape-aware surface reconstruction from sparse 3D point-clouds[END_REF], mais avec une vision déterministe du problème, et une application pratique dans des cas réels de segmentation par des praticiens. En effet le plus souvent les méthodes de segmentation sont développées dans un but spécifique, et il est important de noter qu'il n'y a, pour le moment tout du moins, pas de méthode universelle de segmentation. Chaque algorithme répond à un ensemble unique de contraintes et d'ojectifs, rendant la méthode mise en place particulièrement intéressante pour une application précise et parfois moins pour d'autres. Les objectifs plus précis que nous nous sommes donnés pour le cadre de ce travail sont : -une méthode rapide, -utilisable par des cliniciens, -plus spécifiquement pour la segmentation des tissus mous, tels que des muscles, -et dont les limites sont parfois difficiles à percevoir sur scanner.

Nous nous sommes limités ici à explorer les grandes méthodes existantes, cependant pour plus de détails, il est possible de se référer à [START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF][START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF][START_REF] Geert | A survey on deep learning in medical image analysis[END_REF].

Introduction du concept

Ce travail a été effectué en collaboration avec Damiano Lombardi et Yvon Maday. Le bref article qui en est extrait est présenté dans l'appencice A.

Nous présentons ici une méthode de segmentation d'imagerie médicale semi-automatique fondée sur un modèle, minimisant l'intervention d'un expert. A l'aide d'une base réduite de formes 3D (constituée selon la méthode exposée dans le chapitre précédent, et donc représentées comme des déformations élastiques), nous pouvons approximer la forme de l'objet à identifier sur les images, comme un muscle sur un scanner. Il suffit d'identifier les coefficients associés à l'objet d'intérêt, via la résolution d'un système linéaire prenant en entrée les coordonnées de quelques points sélectionnés sur l'image. Cette méthode est indépendante des niveaux de gris de l'image, et peut donc être appliquée sur tous objets et toutes imageries.

Soit E ⇢ R 3 notre forme de référence, x un point de E et x 2 F ⇢ R 3 avec F le bord de l'objet segmenté. On peut extraire de l'ACP de l'ensemble des objects de notre base de données, représenté par des déplacements u de la forme de référence sur ces objets, les modes f i . On obtient alors :

u ⇡ u n = n  i=1 a i j i .
On sélectionne un certain nombre de points x ( j) sur notre image, et on peut alors écrire que :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m.
En choisissant les x ( j) ⇤ tels qu'ils soient les points de E les plus proches de x ( j) , on obtient le système linéaire suivant à résoudre :

8 > < > : Fa = d, F jik = j i,k (x ( j) ⇤ ), d j,k = x ( j) k x ( j) ⇤k , (1)
Un exemple implémenté en 2D est proposé Figure 22. À gauche de la figure, on peut voir la forme de référence qui est transportée vers différentes formes cibles, constituant ainsi notre base de données de formes 3D représentée par des champs de déplacement. À droite, la reconstruction (en rouge) d'une image générée aléatoirement (en bleu) par la résolution d'un système linéaire, permettant d'identifier les coefficients de notre base réduite.

avec (µ 1 , µ 2 , µ 3 ) 2 [0, 1] 3 , b = 1.5, et k = l(m + 1).
Il suffit d'un simple graphe, Figure 23, pour voir l'utilité de cette paramétrisation sur la décroissance des valeurs propres après l'ACP, par rapport à une ACP faite directement sur les harmoniques sphériques. Les nombreuses difficultées rencontrées lors des segmentations manuelles sont celles exposées dans la bibliographie ci-dessus. Cela rend le travail extrêmement long et fastidieux, puisqu'il faut en moyenne 1h45 pour segmenter un seul muscle. Pour la segmentation semi-automatique, bien que certaines complications liées intrinsèquement à l'exercice de la segmentation étaient toujours présentes, la rapidité et la simplicité du processus, rend l'ensemble de l'acquisition plus agréable et plus rapide (en moyenne 9m30 par muscle). Une difficulté associée plus spécifiquement à cette méthode a tout de même perturbé les opérateurs : il s'agit de l'impossibilité de visualiser en temps réel le modèle 3D de la segmentation. Cela reste toutefois bien plus rapide et moins pénible pour un praticien, et on peut sans mal espérer améliorer l'algorithme nous permettant de calculer l'alignement, qui est le plus chronophage.

( W ! R x 7 ! F (x 1 , x 2 , x 3 ). (3)
Par ailleurs une brève étude de reproductibilité inter-opérateur a été menée. Que ce soit en terme de temps de segmentation ou d'erreur de segmentation, notre méthode semi-automatique de segmentation est reproductible d'un opérateur à l'autre. Tandis que la segmentation manuelle n'est, en terme de temps, pas reproductible. Il y a donc un gain dans la fiabilité des résultats avec notre méthode.

Partie III -Reconstruction faciale

La dernière partie présente d'abord le « morphing », outil dont nous nous sommes déjà servis dans le processus permettant de créer une base réduite de formes 3D, et que nous allons maintenant aussi appliquer à la reconstruction faciale, dans le second chapitre de cette partie.

Le morphing

Ce travail a été effectué en collaboration avec Maya de Buhan et Chiara Nardoni, il a en partie déjà été présenté dans les travaux de thèse de cette dernière, et va donner suite à une publication commune.

L'objectif général du morphing est de trouver une transformation d'une forme de référence donnée vers un objet cible. Ici nous proposons de définir cette transformation comme une déformation élastique obtenue en minimisant une fonctionnelle fondée sur la distance entre les deux formes. Pour accomplir cet objectif, nous utilisons des méthodes d'optimisation de forme. L'optimisation de forme consiste à optimiser une fonctionnelle de coût F(W) dans un ensemble de forme admissible O ad . Un problème générique d'optimisation de forme peut s'écrire : min

W2O ad F(W), ou max W2O ad F(W).
Les applications de ce type de problème sont vastes et amènent des problèmes aussi bien théoriques que numériques à résoudre : la définition du domaine admissible, la notion de dérivation dans cet espace, les propriétés géométriques du domaine, etc. De nombreuses notions de dérivation par rapport au domaine sont décrites dans la littérature, nous ferons ici référence à la méthode de variation de bord d'Hadamard [3,96,[START_REF] Francois | Quelques résultats sur le contrôle par un domaine géométrique[END_REF].

Notre méthode a plus en commun avec Simon et al. [START_REF] Konrad | Elasticity-based matching by minimising the symmetric difference of shapes[END_REF][START_REF] Konrad | A linear elastic force optimization model for shape matching[END_REF] car elle propose un modèle de déformation inspiré de la physique, et avec Bacjcy et al. [START_REF] Ruzena | Multiresolution Elastic Matching[END_REF] qui emprunte des techniques au contrôle optimal et à l'optimisation de forme. Pour une bibliographie détaillée de toutes ces méthodes, il est possible de se référer à [START_REF] Remco | Shape matching : Similarity measures and algorithms[END_REF]. Nous présentons ici deux fonctionnelles de forme : une « volumique » J et une « surfacique » P, puis à l'aide de différentes simulations numériques, nous comparons leurs résultats.

Une fonctionnelle volumique

Nous mesurons la différence entre la forme de référence W et la cible W T par la fonctionnelle de forme suivante :

J(W) = Z W d W T (x)dx,
qui implique d W T la fonction distance signée à W T , définie comme suit :

8x 2 R d , d W T (x) = 8 < : d(x, ∂ W T ) if x 2 W T , 0 if x 2 ∂ W T , d(x, ∂ W T ) if x 2 c W T .
On peut alors calculer la dérivée de forme de J :

8q 2 W 1,• (R d , R d ), J 0 (W)(q ) = Z ∂ W d W T q • n ds. avec q : R d 7 ! R d un petit difféomorphisme.
Puis trouver la direction de descente associée, et à l'aide d'un algorithme itératif qui minimise J. On obtient alors une séquence de forme de plus en plus proche de W T . Cet algorithme a été implémenté avec la méthode des éléments finis en 2D et en 3D.

On peut voir avec la Figure 27, la déformation d'un disque en un lapin, obtenu avec la minimisation de J : Nous avons aussi comparé la fonctionnelle J avec une fonctionnelle plus classique de surface P, telle que :

P(W 0 ) = Z ∂ W 0 (d(s, ∂ W T )) 2 ds,
qui dépend de la distance euclidienne d(., ∂ W T ) à ∂ W T , évaluée au bord de W 0 . La fonctionnelle P est plus intuitive, puisqu'elle essaie d'aligner deux formes en fonction de leurs bords et a l'avantage d'être bien définie sur une surface et pas une forme. Cependant, en pratique, cette fonctionnelle s'avère numériquement plus difficile à mettre en place. En effet, la dérivée de forme de P implique de calculer la courbure moyenne au bord :

8q 2 W 1,• (R d , R d ), P 0 (W)(q ) = Z ∂ W ✓ ∂ d 2 ∂ n + kd 2 ◆ q • n ds, avec k la courbure moyenne de ∂ W.
Par ailleurs, pour minimiser P, un point de ∂ W doit se déplacer vers le point le plus proche de ∂ W T , ce qui peut générer une ambiguité, menant à des singularités numériques, empêchant la convergence globlale de l'algorithme.

On peut voir sur la Figure 28, un premier exemple où la fonctionnelle surfacique parvient au même résultat que la fonctionnelle volumique (Figure 6), mais sur un exemple un peu plus compliqué, des singularités apparaissent. La reconstruction faciale a pour but d'estimer l'apparence faciale d'un individu à partir de la seule connaissance de son crâne. Cela se base sur l'hypothèse qu'il existe une importante corrélation entre la structure osseuse et celle des tissus mous [START_REF] Isabel | Past, present, and future of craniofacial superimposition : Literature and international surveys[END_REF], ce que nous avons aussi confirmé dans le premier chapitre de cette thèse. Toutefois, ce processus reste une approximation de la forme du visage, en effet les caractères secondaires comme les cheveux, la couleur de peau ou des yeux ne sont pas identifiés via ces méthodes. Par ailleurs, d'autres paramètres qui eux influencent la forme même du visage ne sont pas toujours faciles à identifier sur le crâne comme l'âge ou le sexe, voire pas du tout comme la corpulence. Certaines études proposent alors un ensemble de reconstructions possibles pour pallier ce problème [START_REF] Peter | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths : methodology and validation[END_REF][START_REF] Sachio | The computer-aided facial reconstruction system[END_REF]. D'un point de vue mathématique, on dit que le problème est mal posé, ce qui va engendrer deux types de difficultés : la caractérisation de la solution qui n'est pas unique, celle du degré de précision obtenu via ces reconstructions, ainsi que la stabilité.

De nos jours, la reconstruction faciale a divers objectifs de la médecine à l'anthropologie. Dans le projet FaciLe, l'objectif est la reconstruction d'un visage en médecine légale afin d'identifier la personne décédée.

Une longue histoire L'histoire de la reconstruction faciale est assez fascinante. En effet, on en retrouve des traces très tôt dans les civilisations humaines, dans les pratiques liées aux ancêtres ou les rites mortuaires. En 1930, une collection de crânes humains plâtrés (Figure 29), datant d'il y a 9000 ans, est retrouvée sous le sol d'une maison néolithique à Jericho [START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF]. FIGURE 29 -Crânes plâtrés de Jericho (7000 av JC) préservé au British Museum de Londres (à gauche) et Jordan Archeological Museum d'Amman (à droite). Images extraites de [START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF] Dans ce contexte la signification symbolique est prépondérante sur l'anatomie. Cependant, on peut considérer ces visages comme une des premières tentatives de reconstruction faciale. On retrouve par la suite ce type de pratique dans de nombreuses civilisations. On peut citer l'Egypte antique et les masques mortuaires pharaoniques par exemple.

Ce n'est toutefois qu'à partir du 15 ème siècle, en Europe, que les connaissances en anatomie se développent, en particulier grâce à l'essor de la dissection jusque là interdite. Ce qui va permettre progressivement un gain de connaissances puis d'intérêt pour la reconstruction faciale dans un cadre plus judiciaire. Ce n'est seulement qu'à la fin du 19 ème siècle, voir le début du 20 ème , que la reconstruction faciale entre à part entière dans les sciences légales. Deux écoles de référence se développent alors, l'école Russe et l'école Américaine. L'anthropologue russe Mikhail Gerasimov (1907-1970) pensait que la forme des muscles était très liée aux structures osseuses sous-jacentes [START_REF] Mikhail | The face finder[END_REF][START_REF] Gerasimov | The reconstruction of the face from the basic structure of the skull[END_REF], et qu'il fallait pour faire une reconstruction faciale correcte avoir des connaissances anatomiques et artistiques très poussées. L'école Américaine, dont on peut citer l'anthropologue Wilton Krogman comme pionnier [START_REF] Marion | The Human Skeleton in Forensic[END_REF], fonde sa méthode non pas sur la reproduction de détails anatomiques mais sur des tables de mesures d'épaisseur de tissus mous. Dans les années 80, on voit apparaître une troisième école, qui s'inspire des deux précédentes, et qui est développée par Richard Neave et Carolin Wilkinson. Cette méthode traditionnelle se fonde sur des procédures manuelles telles que : des reconstructions 2D sous forme de portraits [START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF][START_REF] Karen | Forensic art and illustration[END_REF], ou 3D sous forme de bustes (Figure 30). FIGURE 30 -Méthode manuelle 2D (en haut) et 3D (en bas) de reconstruction faciale. Images extraites de [START_REF]Facial reconstruction-anatomical art or artistic anatomy ?[END_REF] Ces méthodes manuelles, bien qu'obtenant souvent des résultats plausibles, sont rapidement critiquées, principalement pour leur manque de standardisation et leur trop grande subjectivité liée à l'artiste anatomiste lui-même [START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF]. Par ailleurs ce type de reconstruction pouvait prendre jusqu'à plusieurs jours de travail à un expert avant d'obtenir un visage. L'apparition des ordinateurs, des méthodes numériques et de visualisation donne une nouvelle opportunité à la reconstruction faciale. Celle de l'automatisation du processus, permettant un gain de temps et une moindre subjectivité des méthodes. De nos jours, la plupart des méthodes digitales existantes associent modèles statistiques et bases de données.

Une méthode qui tient compte de l'anatomie C'est dans ce vaste contexte, que les méthodes développées au sein du projet FaciLe prennent racines. Notre méthode de reconstruction faciale, s'appuie sur le travail préliminaire de Chiara Nardoni et Maya de Buhan au sein du projet.

La méthode que nous proposons s'inspire largement de l'histoire de la reconstruction faciale, en prêtant une attention particulière à l'anatomie musculaire, qui a toutefois été délaissée dans les méthodes automatiques actuelles [START_REF] Carl | The human masseter muscle and its biological correlates : a review of published data pertinent to face prediction[END_REF]. Elle complète la méthode originale développée dans [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], qui se détachait déjà des méthodes statistiques en se fondant sur la déformation élastique d'un masque de parties molles.

Dans un premier temps, et avant de pouvoir mettre en place toute reconstruction, il nous faut constituer une base de données 3D avec des crânes, des visages et des muscles masséters associés pour une même personne. La constitution de cette base se fait comme indiqué dans le chapitre 2 de cette thèse, cela nous permet de relier entre eux l'ensemble des éléments de la base (Figure 31). FIGURE 31 -Chaque crâne de la base est relié à un crâne de référence après l'étape du « morphing », et chaque visage au visage associé à ce crâne de référence, par la même méthode. Il en va de même pour les muscles, et ainsi chaque élément de notre base de données est relié. Image extraite de [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF] Un dernière étape est nécessaire avant de pouvoir faire la reconstruction d'un visage à partir d'un crâne inconnu, celle de créer, pour chaque couple crâne/visage de la base, un masque qui représente les tissus mous (Figure 32). C'est à partir de ces masques que la reconstruction d'un nouveau visage est généré. On applique la déformation élastique d'un crâne de la base, identifié comme proche du crâne inconnu, vers ce crâne inconnu, au masque de tissu mous associé au crâne de la base choisie.

Il est maintenant possible de mettre en place une reconstruction. A l'aide de la base réduite de forme 3D mise en place pour le muscle masséter, nous commençons par ajouter sur le crâne inconnu la reconstruction de ce muscle. Puis le crâne est fermé, à l'aide de la technique de « warping » décrit dans le chapitre 2 (Figure 33). Il ne reste plus que l'étape de déformation d'un masque de la base de données vers notre crâne inconnu, qui intégre déjà l'information supplémentaire contenue dans le muscle masséter. Nous obtenons ainsi un visage pour notre crâne inconnu, Figure 34. Pour cela nous avons développé une application, permettant à des utilisateurs d'identifier à partir de deux photos d'un même individu (face et profil), le visage 3D qui leur paraît correspondre. Quatre types de tests sont proposés, à l'aveugle, où l'opérateur doit identifier à partir d'une photo :

1. des vrais visages segmentés en 3D, 2. des reconstructions effectuées avec la méthode de [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], 3. des reconstructions effectuées avec la méthode exposée précédemment mais où le vrai muscle masséter a été ajouté, 4. des reconstructions effectuées avec la méthode exposée précédemment avec le masséter reconstruit. Une capture d'écran de l'application est visible Figure 11.

Annexes

Je présente dans ces annexes :

la courte note, extraite du Chapitre 3 de cette thèse, publiée dans les Comptes Rendus Mathématique de l'Académie des Sciences, avec Damiano Lombardi et Yvon Maday, le « Read Me » (notice explicative) du logiciel de reconstruction faciale mis en place pendant ma thèse qui permet de rendre l'ensemble du processus totalement automatique, -une application pratique de reconstruction faciale qui a permis des simulations en anthropologie, en collaboration avec Martin Friess.

identified via a computer.

The FaciLe project differs from these methods in two main ways :

-by basing the proposed approach to three-dimensional facial reconstruction not on statistics but on a deterministic method, using partial differential equations (PDEs), numerical analysis and 3D geometric modelling. -by making the entire process completely automatic and without human intervention, after the acquisition of skull data.

Main objective of this thesis

The theme of facial reconstruction is the common thread of this thesis, which describes the whole scientific process : from data acquisition and processing to results in the forensic sciences. The goal is to develop a method complementary to the one developed in Chiara Nardoni's thesis [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], enhancing the overall project, while providing anatomical details in the process. Various problems relating to this new methodology were raised, causing parallel studies to be carried out during this work. These include, for example, the complete automation of the process, or the segmentation of soft tissue in medical images.

Scientific approach : an application of mathematics

Applied mathematics, a branch of mathematics, focuses on the study and application of models to other disciplines, such as physics, chemistry, biology, economics, or medicine. However, they include a large number of sub-domains, partly independent of the applications themselves, such as numerical analysis, control, optimization, game theory and probability.

It is in this vast field that this thesis of applied mathematics dedicated to three-dimensional facial reconstruction for forensic medicine is part of. Different aspects of the problem will lead us to use various mathematical tools, such as differential equations, dimensionality reduction, finite elements...

FIGURE 37 -Three-dimensional facial reconstruction

Modeling and face

A model

One of the first steps in applied mathematics is the modelling of the problem posed. Modelling here refers to the entire process by which mathematics can be used in an experiment-based or observational science to explain, describe or predict a phenomenon. A mathematical model is therefore the identification of ad'hoc quantities that describe well the state under investigation and a set of equation on theses quantites in order to apply mathematical tools, techniques and theories to those. one shape to another. Here we measure the difference between the 'starting' and 'target' shapes according to the following functional shape :

J(W) = Z W d W T (x)dx,
involving the signed distance function d W T to ∂ W T . The signed distance function is a « level set » function of a domain : which for each point x of the space returns the Euclidean distance d(x, ∂ W T ) to the edge ∂ W T , the sign associated with this distance is positive if the point is located inside the domain, negative if not. The overall « mapping » between W 0 and W T is obtained by an iterative process that minimizes J(E (W 0 )) (with E the elastic deformation), inferring at each step a direction of descent, via a gradient algorithm, for the functional J. The optimal conditions associated with this shape optimization problem are calculated using the Hadamard shape derivative framework [3,[START_REF] Jean | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]96]. In practice, the permissible deformation space is restricted by elastic displacement parameterization. The elasticity equations were implemented using the finite element method.

FIGURE 40 -Simple 2D morphing example

This first phase of the work was mainly carried out by Chiara Nardoni's thesis [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], supervised by Maya de Buhan and Pascal Frey. We have carried out some of these studies together with her : particularly the so-called boundary morphing and its comparison with the so-called volumic morphing (detailed in Chapter 6). At the end of this thesis, I have developed a complementary method presented below (and detailed in Chapter 7), focussing on the main aim of this thesis the facial reconstruction.

Anatomy of the Face

However, this elastic method does not take into account the underlying anatomy of the face, and therefore, more precisely, the functional components of the soft tissue mask mentioned above. Articles such as the one in [START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF], however, remind us of the importance of anatomical knowledge, accumulated in various disciplines, on the face and of which the field of facial reconstruction takes little account and in a non systematic way. We therefore turned to our anthropological and surgical colleagues involved in the project to improve this first model, taking better into account the anatomy of the face.

The masseter muscle

The face is mainly, in terms of volume, made up of muscles and fat, covered with skin. Since fat is a component for which there is no information on the skull, it was decided to consider it at the end of the process as a variable to be adjusted once the shape of the face has been identified.

On the other hand, it is easy to identify up to 100 muscles on a face. Two of them, larger than the others, caught our attention : the masseter muscle (Figure 41 Due to its respective functional and morphological importance, we have decided to focus more specifically on the masseter because of its impact on the shape of the face. As can be seen in Figure 41, the masseter muscle is located at the cheek, in front of the ear, along the branch of the mandible. While the temporal muscle, which is located at the temples, above the zygomatic arch, is often thinner and covered by the hair, thus having a less significant impact on our problem at first.

Database and reduction of dimensionality

We then need to create a model of this masseter muscle. In order to best describe it, as simply as possible, we use a method of complexity reduction, a principal component analysis (PCA). Describing the shape of the masseter no longer requires the 5, 000 3D coordinates of the points that constituted its mesh, but only n = 5 to 20 well chosen coefficients, depending on the desired level of precision (Figure 42) : Once the reconstruction of the masseter (and possibly, of the temporal) has been completed, the elastic deformation method exposed in the previous section is applied, this time not from a skull/face mask but from a skull+muscles/face mask. This incorporates more precise information into this process, taking into account anatomy, and emphasizing the link between form and function, the founding idea of this thesis.

Data acquisition

All the work presented above was made possible thanks to the implementation of a database combining skull, face and muscles for the same individual as mentioned above. This database consists of 3D medical images that allow the extraction of all the information necessary for the application of this method. All the scanner data collected to create this shape database comes from the Maxillo-Facial Surgery department of La Pitié Salpêtrière Hospital in Paris. This database does not in any way overlap with the one used by Chiara Nardoni before me, who had obtained her data from another project, with no academic link to the FaciLe project, on the facial reconstruction carried out by Françoise Tilotta and Yves Rozenholc [START_REF] Françoise | Construction and analysis of a head CT-scan database for craniofacial reconstruction[END_REF][START_REF] Tilotta-Yasukawa | Contribution à la reconstitution faciale en médecine légale : Étude préliminaire d'une méthode tridimensionnelle à partir de l'imagerie numérique[END_REF]. Indeed it seemed more relevant to this phase of the FaciLe project to build our own database.

CT-sacn : precision and difficulties

The choice of the type of imaging is not insignificant. Indeed, although MRI has a better precision for soft tissue, the input to our problem is the skull alone. The best 3D imaging to precisely visualize bones is currently the CT scanner, an example of which can be seen in Figure 43. Despite the advantages associated with all the information known by scanner images, it should be kept in mind that many problems inherent to images affect the acquisition and visualization of data in three dimensions, as is the case with artifacts created by the metal used to treat teeth (brackets, crowns...).

In addition to these difficulties, some constraints, this time coming from the model, further complicate the pre-processing of the data, for example : not taking the spine into account in our 3D skull model.

Semi-automatic Segmentation

Segmentation is the identification on images of voxels (three-dimensional pixels) belonging to the same object of interest. In our case three objects : the skull, the face and the masseter muscle.

Segmentation of the skull and face can be done quickly using grey levels. Indeed, it is easy to automatically segment elements that have a high contrast. The results shown in Figure 44 This problem of segmentation of medical images is found in many research fields, which need to use 3D organ databases as a prerequisite. A main limitation is the time devoted to the acquisition of data, which is therefore obtained in limited quantities, and which limits the quality of the results presented in these various studies. It is in this context, where the limitation of data has compromised our results, that we have proposed and developed a semi-automatic segmentation application for soft tissues. Using the reduced basis described above, we know a priori the shape of the object to be identified on the images, in our case the masseter muscle. It is only necessary to identify the coefficients associated with the object of interest, via the resolution of a system taking as input the coordinates of a few points selected on the image x ( j) ⇤ , j = 1, . . . , m., with m the nomber of selected points :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m.
This method is independent of the grayscale of the image, and can therefore be applied to all objects and images.

Validation of the methods developed

All the methods presented have been validated from a mathematical point of view. We also have set up experimental validations, in order to certify their practical effectiveness in the disciplines where they are to be used.

Reconstruction and Forensic Medicine

For the reconstruction of a face from a skull, the easiest way to certify that the method is successful is to check if an individual is able to recognize this face. This is why we have set up a survey application : a user to whom a photo is presented must choose among 5 proposals, the face, described by a surface mesh, which seems to him to be the one of the person in the photo (Figure 45).

This operation is performed in double blind and randomized for 4 categories : the 3D surface of the true face directly segmented, the 3D surface of the simple elastic reconstruction, the 3D surface of the reconstruction with the directly segmented masseter and the 3D surface of the reconstruction with the reconstructed masseter itself.

FIGURE 45 -Facial reconstruction validation application

This application therefore allows to compare the results of different methods on the same dataset. Indeed, the data used in facial reconstruction are often extracted from medical data that are not transferable to other teams, limiting possible comparisons on the methods themselves.

Segmentation

In the case of semi-automatic segmentation, the application developed for this purpose being the use by clinicians or researchers, we have set up a study to compare the time and error of segmentation with this method against the manual segmentation of two practitioners.

Manual segmentation requires accurency on all scanner sections where the masseter is located, while semi-automatic segmentation requires only about 50 points to be identified, as shown in Figure 46.

Main Contributions

The main contributions of this thesis are : -Management of a reduced 3D shapes database fully automatic, however our approach aims to complete their results by focusing on more varied morphologies and measures that are often undervalued.

Within the framework of this work and among the very numerous cephalometric measurements that exist, we have selected five that best complement the existing literature on the subject [START_REF] Ariji | Three-dimensional morphology of the masseter muscle in patients with mandibular prognathism[END_REF][START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF]. These are the angles SNA, SNB, ANB, ANB, FMA and between the Bispinal and Mandibular plane, which can be seen in Figure 50.

As for the measurements taken on the masseter, the community agrees that the most appropriate parameters to describe it are its height, width, thickness and volume [START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF].

FIGURE 50 -The cephalometric measurements we used for our study

The study was conducted on 51 patients in the maxillofacial surgery department of La Pitié Salpêtrière Hospital in Paris. We have identified significant correlation between the shape of the skull and the shape of the masseter muscle.

In our study, the FMA is the bone angle measurement parameter, which is more correlated with masseter measurements. Indeed, it is inversely correlated to the 4 masseter measurements tested (height, width, thickness and volume). This means that a long, wide, thick and voluminous masseter muscle is associated with a phenotype called « short face » with a closed mandibular angle. This result is similar to those found in [START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF][START_REF] Kiliaridis | Masseter muscle thickness measured by ultrasonography and its relation to facial morphology[END_REF].

However, few studies address this problem according to the anterior posterior position of the maxilla, here we show a correlation of the SNA with the thickness and volume of the masseter and the SNB with its height and volume, unlike [START_REF] Munutsugu Kubota | Maxillofacial morphology and masseter muscle thickness in adults[END_REF].

Because of the presence of these significant correlations, these conclusions have therefore allowed us to hypothesize, for our facial reconstruction method, that it would be possible to deduce the shape of the masseter muscle from the skull bones.

Engineering process of building a reduced base of 3D shapes

This work was conducted in collaboration with Yvon Maday and Chiara Nardoni. This is a complete methodology for processing a 3D shape database. All the mathematical methods used have already been individually presented in the literature. However, they are here assembled and adjusted to meet the characteristics of our applications : facial reconstruction and semi-automatic segmentation. However, it is entirely possible to use the same process for processing and managing a reduced 3D database, for different purposes and applications such as video games, molecular biology, etc.

The steps involved

Several steps using various mathematical tools are necessary in order to be able to use our 3D shape database.

1. Data acquisition : in our context, it is the segmentation of medical images. The second part of this manuscript is entirely devoted to this thorny problem. The bibliography and presentation of the methods implemented are proposed in the following section. All the scanner data collected to create this shape database comes from the Maxillo-Facial Surgery department of La Pitié Salpêtrière Hospital in Paris. They were acquired with two types of windows : soft tissue and bone (Figure 51). . The latter allows to create a surface triangulation from a volume of voxels (3D pixels), however they are also known to leave artifacts and invalid elements in the mesh for finite element simulations. It is therefore necessary to add a remeshing step, allowing to find a correct geometry and to control the quality of the desired mesh size. For that we used mmg [START_REF] Charles | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF] Cécile | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF], which allows us to control the Hausdorf distance between the input mesh and the output mesh, as well as other parameters such as the maximum and minimum length of the edges of the mesh triangles. The result of the Marching Cube algorithm and the remeshing algorithm is given in Figure 52.

3. Alignment : often underestimated, the coherent alignment of database objects can have a strong impact on the rest of the process. Indeed, if the data are acquired from several sources FIGURE 52 -Mesh of a masseter muscle after the marching cube algorithm on the left and after remeshing on the right that do not have the same acquisition parameters, the objects may differ from only a rigid movement and/or homothety. It is therefore necessary to correct the translations, rotations and homoteties that affect the data. The first phase consists in centering the objects W in the origin of the reference system using the calculation of their respective barycenter x G defined by :

x G = 1 |W| Z W xdx
with x 2 W and |W| the volume of W. After discretization on the mesh of the object (detailed in chapter 2, page 49) :

x G = 1 2|W| Â i x 2 i n i ds i
with n the external normal on the surface of the object W.

If necessary, this is when a scaling can be done. During the second phase, we calculate the axes of inertia to be able to align them with those of a reference object, using the inertia matrix I W :

I W = Z W (|x| 2 Id x ⌦x)dx
with Id the identity matrix, and x ⌦x the matrix x •x T . These two steps allow to identify all the objects in Figure 53 as identical in shape, if the axes of inertia are aligned and the meshes are compared. It is important to note that our method is particularly attentive not to select specific points to be mapped, even for alignment, because we wanted a fully automatic process. This is made possible by the method described above combined with steps 4 and 5 (warping and morphing) described below.

4. Warping (optional) : This step is only necessary when the object to be created as a reduced shape base is not topologically closed. The algorithm which consists in using the elasticity equations to deform a closed template mesh on an open object, was developed in [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF]. The main idea is presented in Figure 54. 5. Morphing : The idea of this step is to deform a reference form W 0 , in our case an ellipsoid, into an objective form W k , here each muscle masseter from our database. This allows us to consider all the objects in the database as u k displacements of the same reference mesh, which makes any comparison between them much easier to set up. It also eliminates the need to match landmarks, to label on all objects in the database. A complete chapter of this thesis is dedicated to the detailed explanation of this method and its bibliography in the third part.

6. Reduction dimension : In order to extract the main information from our database, and to best represent the shapes it contains as simply as possible, we perform a redcution dimmension : a principal component anaylsis (PCA)or a greedy algorithm. For this purpose we calculate the modes that allow us to represent all the elements of the database by a few coefficients.

the context of a specific application : facial reconstruction. However, automatic or semi-automatic segmentation methods have a wider scope of application. For example, they are widely used in « computer-assisted diagnosis », which has become one of the main research topics in medical imaging and radiology [START_REF] Kunio | Overview on research and development of computer-aided diagnostic schemes[END_REF][START_REF] Kunio | Current status and future potential of computer-aided diagnosis in medical imaging[END_REF][START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Bradley | Computer-aided detection and diagnosis at the start of the third millennium[END_REF][START_REF] Ronald | Road maps for advancement of radiologic computer-aided detection in the 21st century[END_REF] in recent decades. In the 1980s, the concept of automatic or computer-assisted diagnosis was already known from studies dating back to the 1960s : [START_REF] Hc Becker | Digital computer determination of a medical diagnostic index directly from chest X-ray images[END_REF][START_REF] Gwilym | Computer diagnosis of primary bone tumors : A preliminary report[END_REF][START_REF] Phillip | Automated computer analysis of radiographic images[END_REF], but these first attempts were not a real success. Indeed, the automatic detection of organs or lesions has proven extremely difficult.

Nowadays, techniques are more and more numerous, partly thanks to the technological development of computers and acquisition methods. Some tissues, whose contrast with the neighbouring tissue is sufficiently important, such as bones, are easily detected automatically [START_REF] Yan | A new accurate and precise 3-D segmentation method for skeletal structures in volumetric CT data[END_REF][START_REF] Shafry Mohd | 3D bones segmentation based on CT images visualization[END_REF][START_REF] Jing | Fast segmentation of bone in CT images using 3D adaptive thresholding[END_REF]. The same is true for the major organs of the human body and the pathologies associated with them, such as the brain and heart [START_REF] Siqi | Multi-scale structured CNN with label consistency for brain MR image segmentation[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF][START_REF] Thomas | Handbook of medical imaging[END_REF], which have always aroused great interest among carers. The common methods of use for these types of applications are the so-called thresholding or edge detection methods. The first consists in grouping in the same label the pixels which have a close grey level or more exactly contained by one or more threshold values [START_REF] Ramesh | Thresholding based on histogram approximation[END_REF]. The second, which also relies on the detection of grey levels, looks for the boundaries of anatomical regions by analyzing the rapid transitions of these levels, and assembles the homogeneous voxels into a single label [START_REF] Antonio | Evaluation of methods for ridge and valley detection[END_REF]. These two methods are based on the detection of grey levels, which allows experts to differentiate two anatomical regions and perform manual segmentation [START_REF] Dzung | Current methods in medical image segmentation[END_REF]. This type of strategy, mimicking the knowledge of experts, is the first to be implemented. Different classifications of these methods exist, for example, we can quote [START_REF] Daniel | Medical image segmentation : Methods and software[END_REF].

However, more recently, other types of methods have gradually developed, such as deformable models, guided atlases or the massive use of machine learning. They differ, among other things, by the amount of knowledge a priori injected into the method [START_REF] Carlos | Advances in auto-segmentation[END_REF].

Statistical atlases use as a reference one or more already segmented structures used as prior knowledge for new segmentations [START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF][START_REF] Torsten | Quo vadis, atlas-based segmentation ?[END_REF]. The segmentation of a new image is then based on a registration and deformation that allows to define an optimal transformation between the atlas and the new image. Different types of algorithms exist to optimize the deformation [START_REF] Arish | Auto-segmentation of normal and target structures in head and neck CT images : A feature-driven model-based approach[END_REF][START_REF] Thirion | Image matching as a diffusion process : an analogy with Maxwell's demons[END_REF], or the use of the atlas itself [2,[START_REF] Olivier | Efficient selection of the most similar image in a database for critical structures segmentation[END_REF][START_REF] Jinzhong | Atlas ranking and selection for automatic segmentation of the esophagus from CT scans[END_REF]. These segmentation methods have proven effective in major segmentation competitions organized by « the American Association of Physicists in Medicine », or « the Medical Image Computing and Computer-Assisted Intervention Society ».

Deformable models require more information a priori. When the latter are available, characteristic variations in the interest structure shape can be used to statistical shape models. This limits the results obtained from segmentations to anatomically plausible shapes [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF]. However, these methods are also less flexible for the reason mentioned above. In addition, the size and content of the training data also limits the performance of the results. However, this method is still frequently used for pelvic anatomical structures [START_REF] Daniel | Model-based segmentation of medical imagery by matching distributions[END_REF][START_REF] Hans | A 3D statistical shape model of the pelvic bone for segmentation[END_REF][START_REF] Heiko | Automatic Segmentation of the Pelvic Bones from CT Data Based on a Statistical Shape Model[END_REF].

Statistical learning, applied to image segmentation, groups together a wide range of methods in full swing, such as machine vector support, decision trees, deep learning or neural networks [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF]. These are already found in multiple applications such as tumor segmentation in the brain, or in the liver, etc. [START_REF] Jinsa | Lung cancer classification using neural networks for CT images[END_REF][START_REF] Mala | Neural network based texture analysis of CT images for fatty and cirrhosis liver classification[END_REF][START_REF] Malay | Cardiac arrhythmia classification using neural networks with selected features[END_REF][START_REF] Saritha | Classification of MRI brain images using combined wavelet entropy based spider web plots and probabilistic neural network[END_REF]. The results obtained by these algorithms are promising for the applications for which they are developed. But it is not always obvious to extend them to all types of images and often requires a large dataset in view of the human variability underlying the medical data [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF].

The method we propose is based more on these last techniques and borrows some tools from each of them, to combine them in order to achieve our objective. It is comparable to [START_REF] Florian | Shape-aware surface reconstruction from sparse 3D point-clouds[END_REF], but with a deterministic view of the problem, and a practical application in real cases of segmentation by practitioners. Indeed, most often segmentation methods are developed for a specific purpose, and it is important to note that there is, at least for the moment, no universal segmentation method. Each algorithm responds to a unique set of constraints and objectives, making the method implemented particularly interesting for a specific application and sometimes less so for others. The more specific objectives that we have set for the framework of this work are :

-a fast method, -suitable for clinicians, -more specifically for the segmentation of soft tissues, such as muscles, -and whose limits are sometimes difficult to perceive on a scanner.

We have limited ourselves here to explore the main existing methods, however for more details, it is possible to refer to [START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF][START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF][START_REF] Geert | A survey on deep learning in medical image analysis[END_REF].

Introduction of the concept

This work was conducted in collaboration with Damiano Lombardi and Yvon Maday. The brief article that is extracted from this work is presented in appencice A.

We present here a model-based semi-automatic medical imaging segmentation method, minimizing the intervention of an expert. Using a reduced basis of 3D shapes ( constituted according to the method described in the previous chapter, and therefore represented as elastic deformations), we can approximate the shape of the object to be identified on the images, like a muscle on a scanner. It is enough to identify the coefficients associated with the object of interest, via the resolution of a linear system taking as input the coordinates of some selected points on the image. This method is independent of the grayscale of the image, and can therefore be applied to any object and any image.

Let E ⇢ R 3 be our reference shape, x a point of E and x 2 F ⇢ R 3 with F the boundary of the segmented object. We can extract from the PCA of all the objects in our database, represented by displacements u of the reference shape on these objects, the modes f i . We then obtain :

u ⇡ u n = n  i=1 a i j i .
We select a certain number of points x ( j) on our image, and we can then write that :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m.
By choosing the boldsymbolx

( j)
⇤ such that they are the closest points of E to x ( j) , we obtain the following linear system to solve :

8 > < > : Fa = d, F jik = j i,k (x ( j) ⇤ ), d j,k = x ( j) k x ( j) ⇤k , (4)
An example implemented in 2D is proposed Figure 56. To the left of the figure, we can see the reference shape that is transported to different target shapes, thus constituting our database of 3D shapes represented by displacement fields. On the right, is the reconstruction (in red) of an image randomly generated (in blue) by the resolution of a linear system, allowing us to identify the coefficients of our reduced base. In order to test our 3D concept, briefly explained in the previous 2D section, we decided to use the following parameterization, applied to spherical harmonics, to create 100 3D shapes as a basis for reconstruction :

8 < :

a Real k = p µ 3 (l+1) b cos(2pk p µ 1 ) cos(2pk p µ 2 ) a Imag k = p µ 3 (l+1) b sin(2pk p µ 1 ) sin(2pk p µ 2 ) with (µ 1 , µ 2 , µ 3 ) 2 [0, 1] 3 , b = 1.5, and k = l(m + 1).
A simple graph, Figure 57, is sufficient to see the utility of this parameterization on the decay of eigenvalues after the PCA, compared to an PCA made directly on spherical harmonics.

FIGURE 57 -Decrease of eigenvalues after PCA, with (in blue) and without (in black) parameterization

A certain number of points are randomly selected on a new shape, randomly generated using spherical harmonics, which we wish to reconstruct. After resolution of the linear system 4 allowing to identify the coefficients a of the PCA : it is possible to reconstruct by the linear combination of these coefficients and the PCA modes, the shape of the interest. We can then compare the error made between our reconstruction and the real form, to the projection error due to the fact that we only use a limited number of modes and that we therefore make an approximation. In this context, the two errors found are of the same order of magnitude and differ little.

Application to masseter muscle segmentation

This work, still in progress, was conducted in collaboration with Damiano Lombardi and Yvon Maday. The purpose of this article is to compare the use of our semi-automatic method, described above, with the current manual reference method by practitioners. It is therefore, among other things, a practical contextualization of this method.

Two operators with anatomical knowledge of the muscle masseter region, but no extensive expertise in radiology or segmentation, each segmented 20 masseters : first manually, then semiautomatically, at least one month later.

For manual segmentation, we used the Amira software [START_REF] Detlev | Amira : A highly interactive system for visual data analysis[END_REF]. For semi-automatic segmentation, a graphical interface that can be used by clinicians has been implemented (Figure 60). The many difficulties encountered during manual segmentation are those described in the bibliography above. This makes the work extremely long and tedious, since it takes on average 1h45 to segment a single muscle. For semi-automatic segmentation, although some complications intrinsically linked to the exercise of segmentation were still present, the speed and simplicity of the process makes the whole acquisition more pleasant and faster (on average 9m30 per muscle). A difficulty associated more specifically with this method has nevertheless disturbed the operators : it is the impossibility of visualizing in real time the 3D model of the segmentation. However, it remains much faster and less painful for a practitioner, and one can easily hope to improve the algorithm allowing us to calculate the alignment, which is the most time-consuming.

In addition, a brief inter-operator reproducibility study was conducted. Either in terms of segmentation time or segmentation error, our semi-automatic segmentation method is reproducible from one operator to another. While manual segmentation, in terms of time, is not reproducible. There is therefore a gain in the reliability of the results with our method.

Partie III -Facial Reconstruction

The last part first presents the « morphing », a tool that we have already used in the process of creating a reduced base of 3D shapes, and that we will now also apply to facial reconstruction, in the second chapter of this part.

Morphing Algorithm

This work was conducted in collaboration with Maya de Buhan and Chiara Nardoni, it has already been partly presented in the latter's thesis work, and will follow up on a joint publication.

The general objective of morphing is to find a transformation from a given reference shape to a target object. Here we propose to define this transformation as an elastic deformation obtained by minimizing a functional one based on the distance between the two shapes. To achieve this objective, we use shape optimization methods. Shape optimization consists in optimizing a functional with a cost of F(W) in an eligible shape set O ad . A generic shape optimization problem can be written : min

W2O ad F(W), ou max W2O ad F(W).
The applications of this type of problem are vast and lead to both theoretical and numerical problems to be solved : the definition of the admissible domain, the notion of derivation in this space, the geometric properties of the domain, etc. Many notions of derivation in relation to the domain are described in the literature, we will refer here to Hadamard's method of edge variation [3,96,[START_REF] Francois | Quelques résultats sur le contrôle par un domaine géométrique[END_REF].

Our method has more in common with Simon et al. [START_REF] Konrad | Elasticity-based matching by minimising the symmetric difference of shapes[END_REF][START_REF] Konrad | A linear elastic force optimization model for shape matching[END_REF] because it proposes a deformation model inspired by physics, and with Bajcsy et al. [START_REF] Ruzena | Multiresolution Elastic Matching[END_REF] which borrows techniques from optimal control and shape optimization. For a detailed bibliography of all these methods, we refer to [START_REF] Remco | Shape matching : Similarity measures and algorithms[END_REF]. We present here two functional forms : a « volumetric » J and a « surface » P, then with the help of different numerical simulations, we compare their results.

A volumetric functionnal

We measure the difference between the reference form W and the target W T by the following functional form :

J(W) = Z W d W T (x)dx,
which implies d W T the signed distance function at W T , defined as follows :

8x 2 R d , d W T (x) = 8 < : d(x, ∂ W T ) if x 2 W T , 0 if x 2 ∂ W T , d(x, ∂ W T ) if x 2 c W T .
We can then computed the shape derivative of J :

8q 2 W 1,• (R d , R d ), J 0 (W)(q ) = Z ∂ W d W T q • n ds.
with q : R d 7 ! R d a small diffeomorphism.

Then find the associated descent's direction, and using an iterative algorithm that minimizes J. We then obtain a shape sequence that is closer and closer to W T . This algorithm has been implemented with the 2D and 3D finite element method.

We can see with Figure 61, the deformation of a disc into a rabbit, obtained with the minimization of J : We also compared the J functional with a more traditional P boundary functional, such as :

P(W 0 ) = Z ∂ W 0 (d(s, ∂ W T )) 2 ds,
which depends on the Euclidean distance d(., ∂ W T ) to ∂ W T , evaluated at the edge of W 0 . The P functional is more intuitive, since it tries to align two shapes according to their edges and has the advantage of being well defined on a surface and not a shape. However, in practice, this functional is numerically more difficult to implement. Indeed, the shape derivative of P implies calculating the average curvature at the edge :

8q 2 W 1,• (R d , R d ), P 0 (W)(q ) = Z ∂ W ✓ ∂ d 2 ∂ n + kd 2 ◆ q • n ds,
with k the mean curvature of ∂ W.

In addition, to minimize P, a point of ∂ W must move to the nearest point of ∂ W T , which can generate ambiguity, leading to numerical singularities, preventing the algorithm's global convergence.

We can see on Figure 62, a first example where the functional surface achieves the same result as the functional volume (Figure 40), but on a more complicated test, some singularities appear.

In view of these results, it is the first functional that has been chosen for the facial reconstruction application, of which morphing is a key ingredient. FIGURE 62 -2D morphing of a disc into a rabbit correlation between bone structure and soft tissue structure [START_REF] Isabel | Past, present, and future of craniofacial superimposition : Literature and international surveys[END_REF], which we also confirmed in the first chapter of this thesis. However, this process remains an approximation of the shape of the face, as secondary characteristics such as hair, skin colour or eyes are not identified by these methods. In addition, other parameters that influence the shape of the face itself are not always easy to identify on the skull, such as age or sex, or even not at all, such as stature. Some studies then propose a set of possible reconstructions to overcome this problem [START_REF] Peter | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths : methodology and validation[END_REF][START_REF] Sachio | The computer-aided facial reconstruction system[END_REF].

From a mathematical point of view, it is said that the problem is not well posed, which will generate three types of difficulties : the characterization of the solution, which is not unique, the characterization of the degree of precision obtained through these reconstructions and also the stability of the proposed solutions.

Nowadays, facial reconstruction has various objectives from medicine to antropology. In the FaciLe project, the objective is to reconstruct a face in forensic medicine in order to identify the deceased person.

A long history

The history of facial reconstruction is quite fascinating. Indeed, traces of it can be found very early in human civilizations, in practices related to ancestors or mortuary rites. In 1930, a collection of flattened human skulls (Figure 63), dating back 9000 years, was found under the ground of a neolithic house in Jericho [START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF]. FIGURE 63 -Flat skulls of Jericho (7000 BC) preserved at the British Museum in London (left) and the Jordan Archeological Museum in Amman (right). Images from [START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF] In this context the symbolic meaning is predominant over anatomy. However these faces can be considered as one of the first attempts at facial reconstruction. This type of practice is later found in many civilizations, such as ancient Egypt and pharaonic death masks, for example. However, it was not until the 15th century, in Europe, that anatomical knowledge developed, particularly thanks to the rise of dissection, which had until then been prohibited. This will gradually lead to a gain in knowledge and interest in facial reconstruction in a more judicial context. It is only at the end of the 19th century, or even the beginning of the 20th, that facial reconstruction becomes a full part of the legal sciences. Two reference schools were then developed, the Russian school and the American school.

The Russian anthropologist Mikhail Gerasimov (1907-1970) believed that the shape of the muscles was closely linked to the underlying bone structures [START_REF] Mikhail | The face finder[END_REF][START_REF] Gerasimov | The reconstruction of the face from the basic structure of the skull[END_REF], and that it was necessary to have very advanced anatomical and artistic knowledge to make a correct facial reconstruction. The American school, whose anthropologist Wilton Krogman can be cited as a pioneer [START_REF] Marion | The Human Skeleton in Forensic[END_REF], bases its method not on the reproduction of anatomical details but on tables measuring soft tissue thickness. In the 1980s, a third school appeared, inspired by the two previous ones, and developed by Richard Neave and Carolin Wilkinson. This traditional method is based on manual procedures such as : 2D reconstructions in the form of portraits [START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF][START_REF] Karen | Forensic art and illustration[END_REF], or 3D in the form of a bust (Figure 64). These manual methods, although often producing credible results, are quickly criticized, mainly for their lack of standardization and their too great subjectivity related to the anatomist artist himself [START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF]. In addition, this type of reconstruction could take up to several days for an expert to work before obtaining a face.

The emergence of computers, digital and visualization methods gives a new opportunity to facial reconstruction. The one of process automatization, allowing time savings and less subjectivity of methods. Nowadays, most existing digital methods combine statistical models and databases.

A method that takes anatomy into account

It is in this vast context that the methods developed within the FaciLe project take root. Our facial reconstruction method is based on the preliminary work of Chiara Nardoni and Maya de Buhan within the project.

The method we propose is largely inspired by the history of facial reconstruction, with particular attention to muscle anatomy, which has however been neglected in the current automatic methods [START_REF] Carl | The human masseter muscle and its biological correlates : a review of published data pertinent to face prediction[END_REF]. It complements the original method developed in [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], which already differed from statistical methods based on the elastic deformation of a soft tissues mask.

First, and before we can set up any reconstruction, we need to build a 3D database with skulls, faces and masseter muscles associated for the same person. The constitution of this database is done as indicated in chapter 2 of this thesis, this allows us to link all the elements of the database together (Figure 65). FIGURE 65 -Each skull on the base is connected to a reference skull after the « morphing » step, and each face to the face associated with this reference skull, by the same method. The same applies to muscles, so every element of our database is linked. Image from [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF] A last step is necessary before a face can be reconstructed from an unknown skull, that of creating, for each skull/face pair of the base, a mask that represents the soft tissues (Figure 66). It is from these masks that the reconstruction of a new face is generated. The elastic deformation of a skull from the base, identified as close to the unknown skull, is applied to this unknown skull, using a soft tissue mask associated with the skull of the chosen base. It should be kept in mind that it is possible to implement the same method for other facial elements and that the masseter muscle is a first illustration of this process.

All that remains is the step of deforming a mask from the database to our unknown skull, which already integrates the additional information contained in the masseter muscle. This gives us a face for our unknown skull, Figure 68. FIGURE 68 -Reconstruction of a face from its skull (left), to be compared with the real 3D face (right)

Validation

Facial reconstruction methods are particularly difficult to evaluate, compare and validate. Indeed, each team which proposes a solution to this problem has worked on a limited set of data, often of medical origin, and which therefore does not have the possibility to be transferred to other teams for comparison. Moreover, the ultimate objective of facial reconstructions in forensic medicine is the recognition and identification of a deceased person. It therefore seemed important and interesting to us to set up a reproducible validation method that would allow us to understand and compare these methods.

For this purpose we have developed an application, allowing users to identify from two pictures of the same individual (face and profile), the 3D face that seems to correspond to them. Four types of tests are proposed, blind, where the operator must identify from a picture :

1. real faces segmented in 3D, 2. reconstructions performed using the [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF] method, 3. reconstructions carried out with the method described above but where the real masseter muscle has been added, 4. reconstructions performed with the method previously described with the reconstructed masseter. A screenshot of the application is shown in Figure 45.

Appendices

In these appendices, I detail : -the short note, extracted from Chapter 3 of this thesis, published in the Comptes Rendus Mathématiques de l'Académie des Sciences, with Damiano Lombardi and Yvon Maday, -the "Read Me" (explanatory leaflet) of the facial reconstruction software implemented during my thesis, which makes it possible to make the entire process fully automatic, -a practical application of facial reconstruction that have allowed simulations in anthropology, in collaboration with Martin Friess. Results : We found numerous correlations between cephalometric measurements and masseter muscle dimensions in dysmorphic patients. The most strongly correlated parameters was the FMA, with a Pearson's correlation coefficient between -0.440 and -0.340 for all masseter dimensions and always a p-value < 0.01. Some other correlation were proved between the SNA and thickness (p=0.021) and volume (p=0.005) of the masseter, and the SNB and height (p=0.002) and volume (p=0.009) of the masseter.

Conclusion :

This study confirmed previously reported correlations between the masseter dimensions and facial vertical proportions in a sample of patients exhibiting a wide range of maxillomandibular dysmorphoses, but also demonstrated some correlation with the sagittal proportions of the face.

Key Words

Cephalometric analysis, Masseter muscle.

Introduction

Human craniofacial skeleton exhibits a wide range of variations [START_REF] Beni | The pattern of craniofacial associations : A morphological and methodological correlation and factor analysis study on young male adults[END_REF]. Skull anatomy can be described with different phenotypes mainly defined by vertical and sagittal (antero-posterior) parameters. These parameters allow categorization of most phenotypes into one of the following : short, normal or long face for vertical parameters [START_REF] Tr Gebeck | Analysis-concepts and values. Part I[END_REF], retruded, normal, or protruded for sagittal parameters [7,107]. We can describe the different facial typologies according to three categories based on the vertical dimension : long, balanced, and short-face [START_REF] Opdebeeck | The short face syndrome[END_REF][START_REF] Stephen | The long face syndrome : vertical maxillary excess[END_REF], and three categories based on the Francfort-mandibular angle : hyperdivergent (>30˚), normodivergent (20˚to 30˚), and hypodivergent (<20˚).

Because they insert onto facial bones, notably the jaws, masticatory muscles influence facial growth patterns, and eventually facial typology. Several authors studied the link between masticatory muscle contraction force, facial typology, and their descriptive anatomy [START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF][START_REF] Ingervall | Relation between facial morphology and activity of the masticatory muscles. An electromyographic and radiographic cephalometric investigation[END_REF][START_REF] Tanya | Human mandibular shape is associated with masticatory muscle force[END_REF]. The masseter is the largest muscle of the face [START_REF] Yu | Botulinum toxin A for lower facial contouring : a prospective study[END_REF]. It is strongly related to the shape of the mandible, due to its major insertions on the mandibular angle. Masseter volume and cross-sectional area (CSA) are significantly associated with vertical dimensions of the face [START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF]. It has also been demonstrated that masseter thickness presented the most significant correlation with skeletal dimensions in all directions [START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF]. Masseter muscle volume greatly influences facial morphology and is documented to be strongly related to skeletal anatomy [START_REF] Carl | The human masseter muscle and its biological correlates : a review of published data pertinent to face prediction[END_REF][START_REF] Viviana | The Predictability from Skull Morphology of Temporalis and Masseter Muscle Cross-Sectional Areas in Humans[END_REF]. However, most of these studies used samples of healthy subjects, with regards to jaw discrepancy, or very restrictive definitions of facial typology (e.g. absolute dimensions of the face). Differences were found between the masseter morphology of patients exhibiting mandibular prognathism when compared to healthy subjects but it is unclear how these correlations apply to patients with significant jaw discrepancy.

The aim of this study was to assess the validity of previously reported correlations and to look for other possible correlations between specific masseter shape parameters and facial typology as assessed by cephalometric analysis in a sample of patients representing a wide range of facial phenotypes, including patients with significant jaw discrepancy.

Masseter Measurements

Muscle measurements were performed on axial and coronal reconstructions of the soft tissue window acquisition using Horos R 2.1.1 (Horos Project, GNU Lesser General Public License, Version 3 (LGPL-3.0)). Muscle measurements were performed on axial and coronal reconstructed sections obtained from the soft tissue window acquisition using Horos R 2.1.1 (Horos Project, GNU Lesser General Public License, Version 3 (LGPL-3.0)). Both the right and left masseters were measured in a first series of 20 patients. As no significant differences were found between both sides in this preliminary series, only the right masseter was measured for the following 31 patients and we used the right masseter measures of the whole sample for further analysis of the correlations.

The CT-scans were first reoriented by aligning the axial views onto the Francfort horizontal plane and the sagittal views onto the midsaggital plane. We chose 3 variables of interest to describe the masseter muscles : width (anteroposterior distance in a horizontal plane), thickness (maximal distance between masseter lateral edge and outer cortex of the mandibular ramus), and height (maximal vertical distance between zygomatic arch and mandibular basilar border). Width and thickness were measured using an axial slice through the lingula. Height was measured using a coronal view aligned on masseter maximum height, Figure 1.2.

We performed manual slice-by-slice segmentation of the masseter using Amira R 2019-05 (Zuse Institute Berlin, Berlin, Germany) [START_REF] Detlev | Amira : A highly interactive system for visual data analysis[END_REF] (Figure 1.3 (b)). The volume of the 3D mesh generated from the segmentation was computed using a simple algorithm that we devolpped in C langage, using the 3D tetrahedra of the mesh.

Statistical Analysis

All statistical analyses were done with the statistical package of R R (3.5.3, 2019-03-11, The R Foundation for Statistical Computing). The mean value, standard deviation (SD) and range of each measurement were calculated. We checked that the data distribution followed a normal law with a Shapiro-Wilk test. Pearson's correlation coefficients between masseter variables and cephalometric variables were calculated as well as data plots with linear and non linear regressions. The significance threshold used for all tests was p = 0.05.

Discussion

This study showed that the morphology of the masseter muscle is strongly related to skull cephalometric parameters in a sample of patients exhibiting a wide range of maxillo-mandibular dysmorphoses, the most strongly correlated being the FMA.

We found numerous correlations between cephalometric measurements and masseter muscle dimensions. Most studies in the literature assessed the relationship between muscle thickness and facial height. In our study, the FMA angle was significantly inversely correlated with masseter length, height, width and volume. Moreover, the MMA was also significantly inversely correlated with muscle length, height and width. These results are similar to those found in the literature [START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF][START_REF] Kiliaridis | Masseter muscle thickness measured by ultrasonography and its relation to facial morphology[END_REF]. According to the literature and to our study, thick, long and large masseter muscles are associated with short-face phenotypes characterized by closed gonial angles (hypodivergent), while short, thin masseters are related with long-face phenotypes characterized by open gonial angles (hyperdivergent) [START_REF] Philip | Masseter muscle volume measured using ultrasonography and its relationship with facial morphology[END_REF][START_REF] Munutsugu Kubota | Maxillofacial morphology and masseter muscle thickness in adults[END_REF].

In our findings, the SNA angle (describing the relative horizontal position of the maxilla to the cranial base) was significantly correlated with the masseter muscle volume and thickness. This means that the more protruded maxilla, the thicker the masseter is. The SNB angle (horizontal position of the mandible to the cranial base) was significantly correlated with the masseter muscle volume and height. Again, the more protruded mandible, the thicker and longer the masseter is. These results differ from the literature : Kubota et al [START_REF] Munutsugu Kubota | Maxillofacial morphology and masseter muscle thickness in adults[END_REF] found no associationbetween the masseter muscle dimensions and SNA or SNB. Similarly, Zepa et al did not find any significant difference between the masseter muscle dimensions of Class II and Class III patients [START_REF] Katrina | Three-dimensional evaluation of musculoskeletal in Class II and Class III patients[END_REF]. The mechanism that could explain the SNA angle correlation with masseter dimensions in our findings is unclear. The masseter muscle inserts onto the mandibular angle and the zygomatic arch. On the lingual side of the mandibular ramus lies the medial pterygoid muscle. It occupies the vast majority of the infratemporal fossa, which is located right behind the maxilla. Along with the masseter, this muscle forms the pterygomasseteric sling around the mandibular angle. Because of this continuity and of their identical function, their dimensions are related. Indeed, Weijs et al. evaluated masseter muscle and medial pterygoid muscle morphology of 16 volunteers. They found a strong correlation between the CSAs of the masseter muscle and the medial pterygoid muscle [START_REF] Wa Weijs | Physiological cross-section of the human jaw muscles[END_REF]. Therefore, we hypothesize that this correlation between a more anteriorly positioned maxilla and a thicker masseter could be the result of a bulky pterygoid muscle pushing the maxilla forward.

The ANB angle characterizes the horizontal relationship of the maxilla to the mandible. In our study, this angle was significantly inversely correlated only with muscle height. The greater distance between the maxilla and the mandible, the shorter masseter. However, the relevance of such a correlation based on the total range of variations of ANB, including positive and negative values which correspond to two opposite phenotypes, is questionable. A subgroup analysis based on the sign of the ANB value could be more informative but this was not possible in our series because of the limited sample size.

Becht and al. [START_REF] Michael | Evaluation of masseter muscle morphology in different types of malocclusions using cone beam computed tomography[END_REF] compared masseter muscle length and orientation in a population of 150 patients divided in subgroups according to facial height, anteroposterior skeletal pattern (ANB angle) and symmetry : no significant differences were found between the length of masseter muscle in the three anterior-posterior groups. However, the orientation of the muscle differed according to anterio-posterior facial pattern : the orientation of the masseter was more vertical in Class III compared to Classes I and the orientation of the masseter was more horizontal in Class II compared to Class I. ANB is probably related to the orientation of the masseter muscle rather than to the muscle's height.

In addition to the correlation with craniofacial phenotype, masseter dimensions also depend to some extent on genetics and muscle tissue fiber type composition. The vertical dimension of the face mostly is strongly correlated with muscular fiber composition ; a reverse correlation has been demonstrated between masseter's height and the percentage of type II muscular fibers [START_REF] Sciote | Masseter function and skeletal malocclusion[END_REF]. However, it is still unclear whether the percentage of fiber types in the muscle is the cause or the consequence of the facial type.

Computed tomography allows qualitative and quantitative analysis of muscle morphology [START_REF] Stavros | Ultrasonographic thickness of the masseter muscle in growing individuals with unilateral crossbite[END_REF][START_REF] Munutsugu Kubota | Maxillofacial morphology and masseter muscle thickness in adults[END_REF]. Although the muscular structures are slightly less clearly visible than with Magnetic Resonance Imaging (MRI), measurements obtained on CT-scans or MRI are comparable [START_REF] Kiliaridis | Masseter muscle thickness measured by ultrasonography and its relation to facial morphology[END_REF]. Moreover, CT-scans allow evaluate both craniofacial skeletal morphology and masseter muscle shape. MRI is particularly suitable to analyze muscular tissue. It enables in vivo measurement of many anatomical parameters of the masseter, such as CSA [START_REF] Van Spronsen | Comparison of jaw-muscle bite-force cross-sections obtained by means of magnetic resonance imaging and high-resolution CT scanning[END_REF][START_REF] Van Spronsen | Relationships between jaw muscle cross-sections and craniofacial morphology in normal adults, studied with magnetic resonance imaging[END_REF][START_REF] Wa Weijs | Relationships between masticatory muscle cross-section and skull shape[END_REF] and volume [START_REF] Gionhaku | Relationship between jaw muscle volume and craniofacial form[END_REF][START_REF] Hannam | Relationships between the size and spatial morphology of human masseter and medial pterygoid muscles, the craniofacial skeleton, and jaw biomechanics[END_REF][START_REF] Kitai | Human masticatory muscle volume and zygomatico-mandibular form in adults with mandibular prognathism[END_REF]. However, it is less suitable for bone tissue assessment and cephalometric analysis based on MRI is not validated. Ultrasonography allows in vivo measurement of the CSA and volume of the muscle, but the technique requires a trained practitioner. In our study, we also controlled the reproducibility bias risk by splitting bone tissue and soft tissue measurements, which were done by 2 distinct operators.

In this study, we investigated muscular shape and volume in a sample of patients exhibiting facial dysmorphosis eligible to surgical treatment. In addition to specific post-operative physical therapy, surgical correction of the jaw discrepancy changes the muscle functional pattern. Therefore, it would be interesting to investigate how these functional changes may translate into dimensional changes in the masseter muscle after surgery, as Westermark et al. hypothesized [START_REF] Anders | Three-dimensional osteotomy planning in maxillofacial surgery including soft tissue prediction[END_REF], and to identify a muscular pattern so as to predict relapse.

Although we were able to demonstrate significant correlations, increasing the sample size of our study would allow increase the power of the correlations. This would be especially interesting to investigate separately subgroups based on ANB angle as defined by above or below the norm rather than considering it as a continuous variable. Indeed, increased or reduced ANB angles describe very different facial types which may exhibit specific patterns of correlation with masseter variables. We could thus better describe the link between the morphology of the skull on that of the muscles, which could help predict the shape of one from the other. This could be useful for facial morphology reconstruction or simulation or for computer-assisted surgery applications. A second limitation that can be exposed is the absence of a measurement reproducibility study for the masseter dimensions. This would possibly make it easier to establish strict anatomical criteria that would be comparable from one study to another.

Conclusion

As a conclusion, we showed in this study that the morphology of the masseter muscle is strongly related to the vertical proportions of the skull and, to some extent, to the sagittal, in a sample of patients exhibiting a wide range of maxillo-mandibular dysmorphoses, suggesting that it is possible to deduce the shape of the masseter from that of the skull.

Introduction

The processing of medical images has always been and remains a difficult task. Different mathematical methods have been put in place to extract as much data as possible as easily as possible. Depending on the desired applications, each of these methods has its own advantages and disadvantages. Our objective in this chapter is to present a method for creating a 3D database of organs that is easy to use for different applications, of which we will detail two examples in Parts II and III, respectively medical image segmentation and facial reconstruction for forensic medicine.

Modeling anatomical elements such as muscles or bones in 3D is now increasingly useful for multiple applications such as reconstruction, simulation or even surgery [64,[START_REF] Mullermuller | Morphometric Analysis of Human Bone Biopsies : A Quantitative Structural Comparison of Histological Sections and Micro-Computed Tomography-Comparison of high-resolution computed tomography and histologic sections[END_REF][START_REF] Ranjanomennahary | Comparison of radiograph-based texture analysis and bone mineral density with three-dimensional microarchitecture of trabecular bone[END_REF]. Building a database of anatomical organs is the first step in this type of research. This therefore represents a major challenge for physicians and can have significant implications for clinical practice [START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Stefanie | A preclinical micro-computed tomography database including 3D whole body organ segmentations[END_REF]. Indeed, this type of database can be used to improve diagnosis or prognotic processes [START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Kim | Medical information service system based on human 3D anatomical model[END_REF][START_REF] Nobutaka | BodyParts3D : 3D structure database for anatomical concepts[END_REF], or even facilitate the identification of anatomical structures for students in training [START_REF] Sohyun | The impacts of three-dimensional anatomical atlas on learning anatomy[END_REF].

The set of all variations of a given anatomical element constitutes a huge sampling (based on the total human population, ie about 7.5 billion) of what can be considered as a manifold, M ie a continuous set that varies in a rather regular way. The complexity of this manifold informs on the diversity of the given element, which, most of the times, is rather low. In order to mathematically apprehend this complexity, the notion of Kolmogorov n-width is rather helpful. It is defined by

d n (M , X) = inf X n ⇢X sup f 2M inf g2X n || f g||
where X denotes the way differences between two elements in M is measured (it is associated to a Hilbert space that contains M ). The Kolmogorov n-width measure to which extend all elements of M can be represented as a linear combination of very few basic structures. This mathematical representation is associated in the medical literature to anatomical atlases that is based on a rich enough sample that represents quite well the variability of the given anatomical element.

Several problems emerge from this definition of anatomical atlases for the construction of a digital 3D atlas [START_REF] Michel | An algorithmic overview of surface registration techniques for medical imaging[END_REF][START_REF] Montagnat | A review of deformable surfaces : topology, geometry and deformation[END_REF] :

1. how to identify those characteristic structures that best describe stability and anatomical variability ? 2. how to match all the objects in the database in order to compare them as well as possible ? 3. how to model a new structure using this database ? Currently most of the construction methods of these atlases are largely automatic, so it is more interesting to look at a clearer classification [START_REF] Maxime | Modèles Statistiques de la forme d'organes du corps humain. Application à la recontruction faciale[END_REF]. There are anatomical atlases of different structures : brain, heart, skull, vertebrae, liver, mandible, lungs ; which can be classified according to whether these organs are rigid (bone), elastic (such as brain or liver), or elastic and articulated (such as heart or face) [START_REF] Darrell | Brain Atlases : Their Development and Role in Functional Inference[END_REF]. Here we will be particularly interested in the representation chosen for the elements of the atlas.

An anatomical structure can be represented by different common generic and invariable characteristics : set of points, lines, surfaces or volumes. However, there may be a difference between the representation chosen for the statistical atlas and that used by the mapping method. Finally, some atlases accumulate representations : points and lines, lines and surfaces [START_REF] Stefano | Representation, analysis, and recognition of 3D humans : A survey[END_REF].

The representation by characteristic points, or landmarks, is the representation traditionally used for the anatomy of the head. It is still widely used because the characteristic points allow a simple statistical modeling of the shape to be characterized. There are many cephalometric atlases that have been in place for several decades in 2D and 3D [28, [START_REF] Delaire | L'analyse architecturale et structurale cranio-faciale (de profil). Principes théoriques. Quelques exemples d'emploi en chirurgie maxillo-faciale[END_REF][START_REF] Alexander | In vivo reliability of 3D cephalometric landmark determination on magnetic resonance imaging : a feasibility study[END_REF][START_REF] Marécaux | A new 3D cephalometric analysis for planning in computer aided orthognatic surgery[END_REF][START_REF] Manuel | A 3D cephalometric protocol for the accurate quantification of the craniofacial symmetry and facial growth[END_REF]. The automatic detection of cephalometric points has also been the subject of numerous studies, such as those of [START_REF] Timothy | Use of active shape models for locating structures in medical images[END_REF][START_REF] Tim | An evaluation of active shape models for the automatic identification of cephalometric landmarks[END_REF] and more recently of [START_REF] Min | Automatic 3D cephalometric annotation system using shadowed 2D image-based machine learning[END_REF]. An example is provided in Figure 2.1. FIGURE 2.1 -2D cephalometric method, Images extracted from [START_REF] Tim | An evaluation of active shape models for the automatic identification of cephalometric landmarks[END_REF] Landmarks play a fundamental role in computational anthropology and are defined as discrete anatomical loci that can be recognized as homologous in all specimens under investigation. In morphometric practices the sparse set of anatomical landmarks are used to build a statistical shape model and study variations. However, we have decided to drop these landmarks because of the difficulties of acquisition, manual or automatic, but also because these points are only really reliable on the bone, and more difficult to transpose on the skin, and more particularly on the « intermediate » soft tissues like muscles. Some authors have already developed landmark-free, fully automated methods for shape analysis. An example with application in paleontologyis presented in [START_REF] Benjamin | A landmark-free method for three-dimensional shape analysis[END_REF]. Our choice was therefore based on surface characteristics.

Many atlases or statistical shape models use a representation of the shape by one or more surfaces modeled by meshes. The hypothesis necessary for the statistical modelling step is the following : each vertex of the mesh at the end of the mapping step is a point which, if not characteristic (no biological correspondence), represents the same place on each copy. These points are then called semi-characteristic points, or semi-landmarks, and were introduced by Bookstein [30]. Our objective here is closer to that of [START_REF] Paul | Detection, visualization and animation of abnormal anatomic structure with a deformable probabilistic brain atlas based on random vector field transformations[END_REF], whose aim is not to achieve a set of variations around this form but to quantify an individual's anatomical deviation from the copies in the database. The methods used are however very different.

With the development of data science and visualization, especially for 3D objects, it has become increasingly important to consider a way to synthesize information. It also allows us to better represent the data. A good way to represent a family of complex object with only a few coefficients, is to build a reduced base [START_REF] Harold | Analysis of a complex of statistical variables into principal components[END_REF]112], with a so-called Principal Component Analysis (PCA), most of the time. The first analysis of this type, applied to the processing of shapes, is presented in [START_REF] Timothy | Combining point distribution models with shape models based on finite element analysis[END_REF][START_REF] Timothy | Active shape models-their training and application[END_REF]. But in order to be able to performe this kind of analysis, a Procruste Analysis like first discribe in [START_REF] John | The Procrustes program : Producing direct rotation to test a hypothesized factor structure[END_REF], with the aim to performed rigid registration, is needed.

In this chapter we present a method to build a database of 3D shapes, easily representative and useful for modeling. We are developing a threaded example of 3D shape representation in medical imaging, in order to create a 3D anatomical atlas [START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF]. We choose to present our methodology on specific example, the skull, the face and the masseter muscle, which is the most voluminous muscle of the face, used for mastication.

The first part presents the pre-processing of the data, where we detail the segmentation, the meshing and the alignement process. The second shows how to describe, compare and therefore model each object in the database. The last part describes an application to facial reconstruction and segmentation, note that it is necessary to keep in mind that this procedure can be applied to many other applications, such as structural classification which is a basic task in molecular biology, or in animation cinema and video game creation, etc. [START_REF] Benjamin | Feature-based similarity search in 3D object databases[END_REF].

Preprocessing

Data Acquisition : Segmentation of CT-Scan

To process data in three dimensions, the first step is to acquire the raw data. This step depends on the type of data needed for a specific problem. In our case, a good way to obtain a 3D representation of anatomical elements is to use medical imaging such as MRI or CT-Scan. The data we collected are CT scans, from the oral surgery department of La Pitié Salpêtrière Hospital, Paris, France.

Segmentation is a necessary step to move from a cross-sectional image to a 3D mesh representation, it corresponds to the process of assigning a label to each voxel (3D pixel) of an image. This step can be done in different ways : manually, semi-automatically or automatically. The choice of method depends on the type of image and the object of interest. In our example, the segmentation of the masseter was first manually segmented by anatomical experts, then semi-automatically using the [START_REF] Maday | Fast semi-automatic segmentation based on reduced basis[END_REF] method. In the Figure 2.2, we can visualize an axial section of a CT-Scan, with the selection of the pink masseter muscle, performed manually. Other objects of interest are easier to segment automatically, such as skull bones or face, on the same scanners. However, most of the time, it is necessary to manually modify the result obtained automatically, using multi-thesholding method, due to the occurrence of artefacts, as is often the case with teeth (Figure 2.3). In addition, in our specific case, other adjustments may be of particular interest, such as the removal of the spine (Figure 2.4), or segmentation of the mandible separately. The generation of a 3D triangular mesh from discrete volumetric data is generally called isosurface extraction [32]. The standard tool used for isosurface extraction is the « Marching Cubes » algorithm, introduced in [START_REF] William | Marching cubes : A high resolution 3D surface construction algorithm[END_REF]. The algorithm takes as input a cubic lattice T with an implicit function f defined at each T vertex. In our context, the lattice is given by the set of voxels in which the discrete volumetric data are placed. The binary function f is defined to be positive at a given vertex if it belongs to the considered label, negative otherwise. First, when a T edge is composed of two opposite sign ends, a mesh vertex is created along the edge. Secondly, the previously defined mesh vertices are connected by triangles.

Despite the simplicity of this classic method of « Marching Cubes », it is well known that sometimes it can produce topologically ambiguous configurations or artificial holes, leading to invalid triangulation. To address these problems, several variants have been developed [START_REF] Akio | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[END_REF][START_REF] Leif | Feature sensitive surface extraction from volume data[END_REF][START_REF] Sergey | Approximation of isosurface in the marching cube : ambiguity problem[END_REF][START_REF] Claudio | A modified look-up table for implicit disambiguation of marching cubes[END_REF][START_REF] Gregory | The asymptotic decider : resolving the ambiguity in marching cubes[END_REF][START_REF] Rüdiger | Real-time exploration of regular volume data by adaptive reconstruction of isosurfaces[END_REF][START_REF] Van Gelder | Octrees for faster isosurface generation[END_REF]. For this work, we used the algorithm available in the Amira software [START_REF] Detlev | Amira : A highly interactive system for visual data analysis[END_REF]. In Figure 2.5, we show some example of isosurface-cell intersection in the Marching Cubes algorithm. The yellow (respectively orange) dots correspond to positive (respectively negative) values of f . The orange zone corresponds to the reconstructed surface patch. The whole isosurface is constructed by gluing the local patches. The first line of the Figure 2. 5 (a) show some examples of configurations from the look-up table of Marching Cubes algorithm [START_REF] William | Marching cubes : A high resolution 3D surface construction algorithm[END_REF]. The second line (b) present a configuration from the look-up table of Marching Cubes algorithm [START_REF] William | Marching cubes : A high resolution 3D surface construction algorithm[END_REF] presenting a topological ambiguity. The grey zones correspond to two possible connections associated to the same configuration. At last, the bottom line (c) Topological incoherence generated from the look-up table of Marching Cubes original algorithm [START_REF] William | Marching cubes : A high resolution 3D surface construction algorithm[END_REF] (left) and repairing (right).

Mesh and remesh

To represent a 3D object, we thus use a triangular surface mesh. As pointed out in the previous section, triangulations obtained from « Marching Cubes » procedure contain, in general, a lot of ill-shaped elements which are redundant and oversampled to correctly describe the geometry of the model for computation. We are thus confronted with a problem of surface remeshing. The quality of a triangulation is strongly related to the targeted application and can involve vertex sampling, regularity, complexity, size of the elements or shape. Very often, the initial triangulation provided as a support for the geometry knowledge is not appropriate for a numerical computation due to very poor quality elements. It is therefore important to assure that each item of the database has a valid mesh.

Generally, the remeshing procedure aims at providing a correct and accurate geometric approximation of the underlying 3D model, and a computational mesh of high quality elements suitable for finite elements simulations. We discuss here only some basic ideas, see [5] for futher references.

Software like mmg [START_REF] Charles | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF] Cécile | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF] for the remeshing process can be used with the appropriate parameters, in order to obtain a valid mesh with a computable number of triangles. The parameters we considered were :

the minimum (h min ) and maximum (h max ) size of an element that can be used for remeshing, -the Hausdorff distance (h hausd ), which controls the proximity authorized between the original mesh and its remeshing, -the mesh gradation (h grad ), i.e. the size variation between neighboring elements.

One of the most important parameter is the Hausdorff distance, which defines the tolerance we accept between the input and the output mesh.

Let (X, d) be a complete metric space, and K(X) all the compacts set of X. Let A and B be compacts set of K(X), we use d A (x) = inf a2A d(a, x) the distance at A. We define their Hausdorf distance by :

∂ (A, B) = max(sup a2A d B (a), sup b2B d A (b)).
According to the inferred geometric information and the input defined by the user parameters the method constructs an isotropic metric size map, to which the final mesh is adapted. Figure 2.6 illustrates the remeshing process and influence of the Hausdorff distance parameter over the remeshing procedure. This value is chosen according to the tolerance allowed in the segmentation process. As illustrated in this same figure, the remeshing procedure also removes the « staircases » artefacts due to the spatial discretization and connectivity. FIGURE 2.6 -A part of a 3D mesh of a masserter muscle, extracted from Amira Software with the marching cube algorithm, on the right ; and remesh with mmg, on the left, with a small Hausdorff distance parameters (top), or with a poor Hausdorff distance parameters (bottom).

Alignment

A recurring problem in areas that use three-dimensional object databases, the alignment step is really important to « standardise » data and to be able to provide a better model. Indeed to compare our data we first need to realize that there is no identical data up to a rigid movement and an homothety.

A rigid movement is defined as a transformation that, when acting on a vector v, produces a transformed vector T (v) of the form :

T (v) = Rv + L where R T = R 1 is a rotation, ( 2.1)
with L a vector representing a translation of the origin.

An homothety can be defined as a transformation of an affine space determined by a point S that refers to its center and a nonzero number l that refers to its ratio, which gives

8M, M 7 ! S + l ! SM. (2.2)
These transformations : rotation, translation and homothety are called similarities, and do not modify the shape of the object. Thus, once we have eliminated, the uncertainity generated by those similarities, we are certain to compare the shape of each object of the database. Therefore in practice, to be able to compare the shape of two 3D objects, we first need to center, scale and give them the right orientation accordingly.

Centering the objects at the origin We call « centering the object », the translation of the object to the origin point Õ = (0, 0, 0). In order to do that, we choose to compute the barycentric coordinates. And subtract them from the coordinates of all the object's points.

We write x a point of R 3 , with components (x 1 , x 2 , x 3 ). Let xG be the barycenter of the object of interest W, we write (x G 1 , x G 2 , x G 3 ) its coordinates. We have :

xG = 1 |W| Z W xdx, (2.3)
with x 2 W and |W| the volume of W.

With the consideration of the first component of xG , and we note that x 1 = div((

x 2 1
2 , 0, 0)), then from the divergence-formula we have

R W div(ũ)dx = R ∂ W ũ.
ñds with ñ being the exterior normal to the surface of W, with components (n 1 , n 2 , n 3 ), allows to transform a volumetric integral to a surfacic one. We have :

x G 1 = 1 2|W| Z ∂ W x 2 1 n 1 ds . (2.4)
By discretizing on our mesh, and we get :

x G 1 = 1 2|W| Â i x 2 i,1 n i,1 ds i , (2.5)
with ds i the area of the triangle of interest.

We proceed similarly for x G 2 and x G 3 . We thus obtain the barycentric coordinates of W, which we can subtract from the coordinates of all the points x of W to center up our object at the origin.

Rotation : alignement of inertia axes

In a second step, we wish to orient all our objects accordingly. The goal is to be able to identify as a unique object two surface-meshes which would only differ from a rotation.

Let W be a domain containing the origin Õ as barycenter. The angular momentum of W in relation to a direction axis d when the solid is rotating around d is given by :

M( d) = Z W x ^w(x)( d ^x)dx (2.6)
with w being the angular velocity, d ^x being the direction vector of the displacement of the particle x (instentaneous direction). This is a generalisation of the formula giving the angular velocity of a sum of particles. Therefore :

M( d) = wI W d, (2.7)
where I W is the inertia matrix.

I W = Z W (|x| 2 Id x ⌦x)dx, (2.8)
with Id the identity matrix, and x ⌦x the matrix x •x T .

Proposition 2.2.1 The inertia matrix I W is symmetric, positive definite, therefore diagonalisable in an orthogonal basis.

Along the components, the inertia matrix is written as :

I W 3,3 = Z W (x 2 1 + x 2 2 + x 2 3 ) x 2 3 dx 1 dx 2 dx 3 , (2.9)
= Z W (x 2 1 + x 2 2 ) dx 1 x 2 x 3 .
(2.10)

Then we obtain :

I W 1,1 = Z W (x 2 2 + x 2 3 ) dx 1 dx 2 dx 3 , (2.11)
I W 2,2 = Z W (x 2 1 + x 2 3 ) dx 1 dx 2 dx 3 , (2.12)
(2.13)

And :

I W 1,2 = Z W x 1 x 2 dx 1 dx 2 dx 3 , (2.14)
I W 1,3 = Z W x 1 x 3 dx 1 dx 2 dx 3 , (2.15)
I W 2,3 = Z W x 2 x 3 dx 1 dx 2 dx 3 . (2.16)
(2.17)

Computation of the coefficients of I W : we transform the expressions below using the Green formula, similarly as the barycenter.

For instance :

I W 1,2 = Z W div 0 @ 0 0 x 1 x 2 x 3 , 1 A dx 1 dx 2 dx 3 (2.18) = Z ∂ W x 1 x 2 x 3 n 3 dx 1 dx 2 dx 3 . (2.19)
The same holds for the other components.

After discretization, as for the barycenter, we obtain for (2.19) :

I W x i,1 x i,2 = Â i x i,1 x i,2 x i,3 n i,3 ds i (2.20)
Using these discretized formula, we implemented this procedure on a simple case first (Figure 2.7). The same 3D brick has been rotated in three different orientations on the right of first line of this figure and, after using our algorithm, they are all set back the same direction as the brick of the left of the figure, on the bottom line.

In Figure 2.8 we give another exemple, more specific to our problem and a less obvious than the previous one : the masseter muscle. This whole process is quite useful to detect rigid motion differences between two objects and align them. However in some cases, depending on the goal, we need more accuracy on a specific part of the object. For example on our CT-Scan, we are interested in aligning more specifically the eyes, the nose, or/and the mouth more than the whole head, especially when the acquisition prameters or the position of the patients are not all the time the same, hence affecting the data. In this case, we can add a last step, presented below.

Iterative closest point (ICP)

This algorithm (ICP) is employed to minimise the difference between two clouds of points. One is the fixed target, while the other one (the source) is aligned to best match the reference. This algorithm iteratively transforms (through translations and rotations) the clouds of points of the source to minimise the distance between two clouds i.e. the sum of the squared differences between the coordinates of the matched pairs is minimised.

This method was first introduced by [START_REF] Paul | Method for registration of 3-D shapes[END_REF]. As suggested in the section above, we can use the ICP algorithm to obtain a better accuracy : indeed for the same « face » example we can only select the points of the region of interest, and make the others follow the same transformation.

However, if the first two steps are not performed before, i.e. the translation and the rotation, this algorithm may fail. Indeed, a first good initialization is essential to the ICP algorithm for an pertinent result. Most of the time this is done via landmarks selection (semantically identifiable point on the object), but in our facial reconstruction and segmentation applications, it was precisely important to get rid of it. This is one of the interesting points of this fully automatic method of processing 3D data for a reduced shape basis : the matching of shapes is not done by a selection of landmarks.

Function-to-shape representation

In order to make the link between the items of our database and to be able to summarize information, we need a way to compare them. The preprocessing of the data was needed to eliminate identical cases which differ only from size or alignement, and obtain a smaller Kolmogorov dimension of the object family. This being done, in order to find a way to summarize the information, we need comparable meshes with the same number of points, triangles, and with the same topology for each object of the database.

The idea being to be able to match all objects of the same type in the database, such as skulls for example. As we have decided not to select landmarks, we will define semi-landmarks, via a dense mesh used as a deformed template on all objects in the database. For that we propose as a model an elastic shape matching from a reference shape W 0 to each object of the database W T j . A supplementary step is needed before if the object is not closed. This case does not appear in the face or muscle case, but for creating a database of skull one would need to perform this first warping step.

Warping (if necessary)

As mentioned above, if the object of interest is not closed, this step is mandatory before performing the elastic shape matching. A relevant example of this problem, in our case, is the skull (Figure 2.9). FIGURE 2.9 -Warping a skull from a sphere Warping, as described in [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], is a tool allowing to define a closed surface model from an unorganized source model. In this context, we focus on the reconstruction of a connected shape starting from the discrete data of a source triangulation. Such a triangulation could present one or more undesirable particularities : invalid, not closed, or topologically intersected, etc. Here we use an algorithm based on mesh deformation techniques, where a closed surface mesh model is iteratively deformed, producing a sequence of surfaces which are "closer and closer" to the source triangulation, in the sense of the Hausdorff distance.

Morphing

This step really enables us to describe and put in correspondance every object of the database. Once performed, we will no longer consider a database of 3D mesh objects but of displacement a function with the same number of components for every object.

The aim of a morphing algorithm is to transport a template to a source object by minimising a shape functional, as described in more details [START_REF] Buhan | Volumic and Boundary shape matching : comparison of two optimization methods[END_REF][START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF], and also in the Chapter 6 of this thesis.

For reasons detailled in Chapter 6, this algorithm need a volumic mesh to be processed. Indeed, we have shown that surfaced morphing is not very reliable for complex shapes and is difficult to implement in 3D, so choosing volume morphing is simpler and faster. Thus we add a brief step of 3D meshing of the surface shape database. We can briefly described the used method as follow.

Let W 0 and W T ⇢ R d be respectively our ellipsoid template E (Figure 2.10) and one "target" masseter of our data T i . As they share the same topology, our purpose is to map W 0 onto W T by minimising an energy functional which measures the discrepancy between the shapes. In order to achieve our goal, we introduce the following functional J(W) :

J(W) = Z W d W T dx, (2.21)
where W = F(W 0 ) the ensemble of eligible shapes for F the elastic transformation, and d W T is the signed distance function to W T , defined as :

8x 2 R d , d W T (x) = 8 < : d(x, ∂ W T ) if x 2 W T , 0 if x 2 ∂ W T , d(x, ∂ W T ) if x 2 c W T . (2.22)
In order to decrease the value of J(W), the domain W must expand in the regions of the ambient space R d where d W T is negative, and to retract in those where it is positive. By deforming the template shape as above we can always obtain better values for J until the deformed template shape reaches the target shape W T .

Using a classical calculation (see [3,96]), the functional (2.21) is shape differentiable at any masseter shape W 2 O ad , with O ad a set of admissible shape, and the associated shape derivatives reads :

8q 2 W 1,• (R d , R d ), J 0 (W)(q ) = Z ∂ W d W T q • n ds. (2.23)
Then a descent direction is inferred from previous formulas and an iterative algorithm devoted to the minimisation of J produces a sequence of shapes which are 'closer and closer' to W T .

First we fill all the considered shapes W with a linear elastic material. Then, we assume that any such shape W contains a given subset w b W on which it is clamped. We now obtain a descent direction for J(W) as the unique solution u W belonging to

H 1 w (W) d := v 2 H 1 (W) d , v = 0 2 w of the variational problem : 8v 2 H 1 w (W) d , Z W s (u W ) : e(v) dx = J 0 (W)(v) = Z ∂ W d W T v • n ds. (2.24)
where e(v) = 1 2 (-v + -v T ) is the linearised strain tensor and s is the associated stress tensor via the Hooke's law with Lamé coefficients l , µ : s (u) = 2µe(u) + l tr(e(u))I.

(2.25)

Figure 2.10 shows our ellipsoid template E which is deformed on 3 target shape T 1 , T 2 , T 3 , and we obtain three result shapes F 1 , F 2 , F 3 described by their displacement function u 1 , u 2 , u 3 . FIGURE 2.10 -Ellipsoid map with elastic morphing to three masseter muscles At the end of this step, each of the shapes in our database is no longer described by its original mesh obtained by a Marching Cube algorithm, but by the mesh of a reference deformed by elastic morphing. The high density of the template mesh and the prior alignment of the shapes between them allows us to consider each shape in point-to-point correspondence. It is important to note that at this stage, only the surface mesh is retained, the inner tetrahedrons being added only to facilitate the numerical calculation of the elastic morphing.

Parametrization : Dimension Reduction

We can now consider each object of the database as a surfacic displacement u of the reference mesh to each of them. To extract the information of this displacement database, we choose to performe a parametrization. At first with a PCA method, because it is a well known method easy to implement. Then a second time to be more coherent, and keep all the process deterministic, we present a greedy method.

Principal Components Analysis

PCA is a method used to reduce the complexity of the data, if it is in space with a small Kolmogorov dimension. This allows for a large number (in our case we have acquired 126 snapshots) of similar functions u k to be represented with just a handful of well-chosen linearly combined components.

It involves computing the variance-covariance matrix A with entries a k`= hu k , u `i (where h., i stands for an appropriate scalar product in the space of anatomical displacement items u). Then we extract the eigenvectors and eigenvalues of A sorted in decreasing order, and construct the reduced basis :

j `= 126 Â k=1 v l k u k , (2.26)
where v k are the components of the `the eigenvectors of A, where we retain only the first few ones : `= 1, . . . , n, n << 126.

These modes allow to approximate, by linear combination, all the other object of the database :

a k 1 j 1 + a k 2 j 2 + ... + a k n j n , (2.27)
where a k `= Âhuk,j`i Âhj`,j`i is a good approximation of each u k with very few modes (n of the order of 5). Thus an approximation of an unknown muscle can be proposed by determining its coefficients a `, `= 1, . . . , n.

At first we used the L 2 scalar product to compute the matrix A, and we obtained the presented decrease of eigenvalues, Figure 2.11. FIGURE 2.11 -Decrease of the 10 first eigenvalues for our masseter muscles example A possible improvement, to expect a better representation of the eigenvalues, would be to use the H 1 scalar product. Nervertheless it is necessary to compute the gradient in this case, and it needs to be done during the resolution of the elastic morphing.

Once the whole process is performed, from the data acquisition to the PCA, we only need a few coefficients (n between 5 and 15 depending the precision needed and the application) to describe an object like the masseter muscle for exemple, instead of 5000 3D coordinates of the mesh points.

Table 2.1 presents the discrepancies between the different reconstructions and the target morphed masseter, using 5, 10 and 20 coefficients. On the left, we can see the error computed with all the element in the database used for performing the PCA. On the right, the PCA modes were computed without the masseter reconstructed, for which we compute the error.

The visualisation of those errors is shown in the next secrion called « 3D Visualisation ». From the series of snapshots of masseters M i , i = 1, . . . , N, that we have, we want to extract a different basis, better suited to our futur goal.

We first choose a « representative masseter » as being the first element in the basis, we name this masseter M ⇤ 1 . The choice of M ⇤ 1 can be random, or it can be the object that best represents all the others (i.e. it minimizes the projection error on each member of the database).

We then approximate each masseter M i in the vetorial space spanned by M 0

M i ' l i 1 M ⇤ 1
with a coefficient l i 1 fitted so that the points that are selected on M i are best approximated (in a mean square sense) by

l i 1 M ⇤ 1 .

Let us define by M ⇤

2 the masseter that is the worse approximated (everywhere, and not only in the reliability region) by l i 1 M ⇤ 1 . We have now a two dimensional approximation by combining M ⇤ 1 and M ⇤ We proceed recursively by defining, at step k an approximation

M i ' k  j=1 l i j M ⇤ j
where the coefficients l i j are fitted so that M i are best approximated (in a mean square sense) by  k j=1 l i j M ⇤ j . and M ⇤ k+1 is chosen as being the masseter that is the worse approximated by  k j=1 l i j M ⇤ j .

Table 2.2 presents the discrepancies between the different reconstructions and the target morphed masseter, using 5, 10 and 20 coefficients. On the left, we can see the error computed with all the element in the database used for performing the Greedy algorithm. This explained why the minimum of the error is 0. On the right, the Greedy modes were computed without the masseter reconstructed, for which we compute the error. The visualisation of those errors is shown in the next secrion called « 3D Visualisation ».

3D Visualisation

We can see on Figures 2.12 and 2.13, the respective reconstructions of two masseters, one being in the database we used for the reconstruction of the modes, the next one not being part of this database. For each of these figures, the top line corresponds to the reconstruction made using the PCA modes and the bottom line is made using the Greedy modes. -we can see that the results are better when the object to be reconstructed is itself in the database -the POD method seem to prove more efficient for the « best reconstruction » possible here, -using 20 coefficients of 126 ellement in the database, most of the error is under 0.5 -1mm in the Figure 2.12, and under 1 mm for the Figure 2.13.

From this easy description of an anatomical object, different applications can be developped, in the next section, we briefly introduce two of them.

Application of reduces basis of 3D shapes 2.4.1 Shape prediction of anatomical element

It is possible to learn the coefficients a `, `= 1, . . . , n, from a training set, which links two objects of the database. In our example of the masseter, we know from the previous chapter, that the shape of this muscle depends a lot of the skull shape. So if we have 2 reduced basis of 3D shapes, one for the muscle and one for the skull, computed from the same CT-scan data set, we can learn the coefficient a of the muscle from the ones of the skull. And so, create muscles from the only knowledge of the skull.

m k L ! a k i
, where m k are measurments of the known object or their own a coefficients extracted from the previous process, L a machine learning process (like a neural network), and a k i the coefficients we want to obtain by learning.

This application is detailled in a full section of the Chapter 7, for facial reconstruction.

Segmentation and reconstruction

Segmentation of soft tissus

As already mentioned before, segmentation of soft tissues is a difficult process. But the ability of describing an object like a muscle with only few coefficients to identify, can make the whole process easier. Indeed by selecting only a few points on medical images, we can resolve a linear problem which gives us these a coefficients, and let us reconstruct the expected shape on the image. All this process is described in details in the next part of this thesis.

Facial Reconstruction

The process of facial reconstruction involves giving an approximation of the face of an unknown person from the only knowledge of its skull shape. An adequate database for Facial reconstruction must contain couples of faces and skulls coming from an healthy individual together with anthropological parameters like age, sex, Body Mass Index (BMI). Since computer science emerge in forensic sciences, this imply a database of skulls, faces and even muscles 3D shapes linked to each other.

Our method can be a really good way to treat those datas and can also be useful for the facial reconstruction application itself. Once the database is acquired as described before, it becomes possible with only some mesurements on the skull (or just with the 3D shape of the skull bones) to reconstruct muscles and/or faces.

All this process is described in details in the last part of this thesis, in the Chapter 7.

Conclusion

This chapter presents a way to describe a 3D shapes database with a reduced basis, after performing some preprocessing step, starting from segmentation. Even if most of the tools presented in this chapter are usually used in this kind of process to creat anatomical atlas, some points nevertheless deserve a particular interest. First all the method is devoted to be fully automatic, and this without selecting landmarks. In addition, the whole procedure is based on a purely deterministic approach, which allows us to build a database of dense three-dimensional surface meshes. Unlike most of the « statistical atlases » described in the literature, which are based as their name suggests on statistical, and sometimes probabilistic, models. The PCA, proposed in this chapter, being thought of as a reduction of complexity, rather than as a statistical regularization, could also be a greedy algorithm.

To conclude we presente two different applications, where this approach is efficient, which are detailled in the two next part of this thesis. Semi-automatic segmentation method in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Second Part : Segmentation (3D pixels) of the entire image or object of interest [68]. This has been for years, the effort of specialists and it is quite time consuming. In Figure 3.1, we try to illustrate this both definition : voxel and segmentation. Beyond the specifics task of each particular application, the general issue is to find a way to identify automatically an anatomical region of interest in order to speed up the process. Indeed automatic or semi-automatic segmentation methods have a wider scope of application. For example, they are widely used in « computer-assisted diagnosis », which has become one of the main research topics in medical imaging and radiology [START_REF] Kunio | Overview on research and development of computer-aided diagnostic schemes[END_REF][START_REF] Kunio | Current status and future potential of computer-aided diagnosis in medical imaging[END_REF][START_REF] Kunio | Computer-aided diagnosis in medical imaging : historical review, current status and future potential[END_REF][START_REF] Bradley | Computer-aided detection and diagnosis at the start of the third millennium[END_REF][START_REF] Ronald | Road maps for advancement of radiologic computer-aided detection in the 21st century[END_REF] in recent decades. In the 1980s, the concept of automatic or computer-assisted diagnosis was already known from studies dating back to the 1960s : [START_REF] Hc Becker | Digital computer determination of a medical diagnostic index directly from chest X-ray images[END_REF][START_REF] Gwilym | Computer diagnosis of primary bone tumors : A preliminary report[END_REF][START_REF] Phillip | Automated computer analysis of radiographic images[END_REF], but these first attempts were not a real success. Indeed, the automatic detection of organs or lesions has proven extremely difficult.

Nowadays, techniques are more and more numerous, partly thanks to the technological development of computers and acquisition methods. When the contrast is high enough (e.g. for bones, or between bones and tissues) automatic segmentation is quite efficient and various methods are proposed and implemented [START_REF] Yan | A new accurate and precise 3-D segmentation method for skeletal structures in volumetric CT data[END_REF][START_REF] Shafry Mohd | 3D bones segmentation based on CT images visualization[END_REF][START_REF] Jing | Fast segmentation of bone in CT images using 3D adaptive thresholding[END_REF]]. One automatic segmentation example is shown in Figure 3.2 (b). Likewise for anatomical major organs like the brain or the heart [START_REF] Siqi | Multi-scale structured CNN with label consistency for brain MR image segmentation[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF][START_REF] Thomas | Handbook of medical imaging[END_REF], which have always aroused great interest among carers. On the contrary, when the contrast is not high enough (e.g. segmentation between different soft tissues on CT scans) the most common imaging modalities fail and there are very few approaches available. One manual segmentation exemple of soft tissues is show, Figure 3.2 (c).

The common methods of use for these types of applications are the so-called thresholding or edge detection methods. The first consists in grouping in the same label the pixels which have a close grey level or more exactly contained by one or more threshold values [START_REF] Ramesh | Thresholding based on histogram approximation[END_REF]. The second, which also relies on the detection of grey levels, looks for the boundaries of anatomical regions by analyzing the rapid transitions of these levels, and assembles the homogeneous voxels into a single label [START_REF] Antonio | Evaluation of methods for ridge and valley detection[END_REF]. These two methods are based on the detection of grey levels, which allows experts to differentiate two anatomical regions and perform manual segmentation [START_REF] Dzung | Current methods in medical image segmentation[END_REF]. This type of strategy, mimicking the knowledge of experts, is the first to be implemented. They are simple and effective approaches to segment anatomical details featured by high contrast, as, for instance, the bones (that appear as light objects on dark background). However it is more challenging for the detection of objects that do not have a marked gray scale signature, as the muscles. Different classifications of these methods exist, for example, we can quote [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF][START_REF] Geert | A survey on deep learning in medical image analysis[END_REF][START_REF] Daniel | Medical image segmentation : Methods and software[END_REF]. Therefore, more recently, other types of methods have gradually developed, such as deformable models, guided atlases or the massive use of machine learning. They differ, among other things, by the amount of knowledge a priori injected into the method [START_REF] Carlos | Advances in auto-segmentation[END_REF].

Statistical atlases use as a reference one or more already segmented structures used as prior knowledge for new segmentations [START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF][START_REF] Torsten | Quo vadis, atlas-based segmentation ?[END_REF]. The segmentation of a new image is then based on a registration and deformation that allows to define an optimal transformation between the atlas and the new image. Different types of algorithms exist to optimize the deformation [6,[START_REF] Arish | Auto-segmentation of normal and target structures in head and neck CT images : A feature-driven model-based approach[END_REF][START_REF] Thirion | Image matching as a diffusion process : an analogy with Maxwell's demons[END_REF], or the use of the atlas itself [2,[START_REF] Olivier | Efficient selection of the most similar image in a database for critical structures segmentation[END_REF][START_REF] Jinzhong | Atlas ranking and selection for automatic segmentation of the esophagus from CT scans[END_REF]. These segmentation methods have proven effective in major segmentation competitions organized by « the American Association of Physicists in Medicine », or « the Medical Image Computing and Computer-Assisted Intervention Society ».

Deformable models require more information a priori. When the latter are available, characteristic variations in the interest structure shape can be used to statistical shape models. This limits the results obtained from segmentations to anatomically plausible shapes [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF]. However, these methods are also less flexible for the reason mentioned above. In addition, the size and content of the training data also limits the performance of the results. This method is still frequently used for pelvic anatomical structures [START_REF] Daniel | Model-based segmentation of medical imagery by matching distributions[END_REF][START_REF] Hans | A 3D statistical shape model of the pelvic bone for segmentation[END_REF][START_REF] Heiko | Automatic Segmentation of the Pelvic Bones from CT Data Based on a Statistical Shape Model[END_REF].

Statistical learning, applied to image segmentation, groups together a wide range of methods in full swing, such as machine vector support, decision trees, deep learning or neural networks [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF][START_REF] Alexander | Accurate automated volumetry of cartilage of the knee using convolutional neural networks : data from the osteoarthritis initiative[END_REF]. These are already found in multiple applications such as tumor segmentation in the brain, or in the liver, etc [START_REF] Jinsa | Lung cancer classification using neural networks for CT images[END_REF][START_REF] Mala | Neural network based texture analysis of CT images for fatty and cirrhosis liver classification[END_REF][START_REF] Malay | Cardiac arrhythmia classification using neural networks with selected features[END_REF][START_REF] Saritha | Classification of MRI brain images using combined wavelet entropy based spider web plots and probabilistic neural network[END_REF]. The results obtained by these algorithms are promising for the applications for which they are developed. But it is not always obvious to extend them to all types of images and often requires a large dataset in view of the human variability underlying the medical data [START_REF] Manaswini | A Survey on Applications of Machine Learning Techniques for Medical Image Segmentation[END_REF].

It is important to keep in mind that most often segmentation methods are developed for a specific purpose, and it is important to note that there is, at least for the moment, no universal segmentation method. Each algorithm responds to a unique set of constraints and objectives, making the method implemented particularly interesting for a specific application and sometimes less so for others.

We propose here to develop a model based semi-automated method to help segmentation, which have most in common with the three last types of methods cited above. In particular, the shape of the object of interest is represented by an elastic deformation field of a reference configuration. A reduced basis representation of the object of interest is introduced, that exploits an a posteriori knowledge of a database of available segmentations. Given an image to be segmented, the elastic deformation field can be approximated as a linear combination of a small number of modes. This provides, in turn, a powerful regularisation method and a significant speed-up of the procedure since only few apropriate points need to be selected by the user in order to define the field. We can find a similar idea in [START_REF] Florian | Shape-aware surface reconstruction from sparse 3D point-clouds[END_REF], but with a probabilistic point of view, and where all the modes of the PCA are used as a regularisation and not as a dimentional reduction. Moreover, PCA is here a quick way to make a first test which will be later, in chapter 4, replaced by a greedy algorithm, making the process entirely deterministic.

The proposed method is complementary to the ones exposed before, first because it requires some preliminary segmentations to « build » the reduced basis, and second because it can be used on indifferently on CT and MR image for any type of object, even the ones difficult to identify like soft tissues on a CT. This work remains a preliminary presentation of the concept in 2D and 3D, we wish to develop. Some details will be better explained in the next chapter, with more realistic and complex examples. While improving the points of weakness that could be highlighted from these first very simplified results.

The chapter is structured as follows : in Section 3.2, the problem and the notation are presented, in the next Section 3.3 some 2D numerical experiments are detailed, then Section 3.4 presents first 3D numerical results.

Notation and problem setting.

Most medical imaging modalities produces a gray-scale image, i.e. a function defined over a subset of R 2 (or R 3 for 3D imaging), f : R 2 ! [0, 1] associating a scalar value to each point of the space (denoted later on by x = (x 1 , x 2 ) (or by x = (x 1 , x 2 , x 3 ) in 3D)). After a manual segmentation of a region of interest is performed, it is possible to model this region by using an elastic morphing that maps a reference object : e.g. an ellipse in R 2 (or an ellipsoid in R 3 ) onto the boundary of the two-dimensional (or three-dimensional) segmented region. In the present work, for the sake of simplicity, we will restrict to two-dimensional configurations. Having a collection of segmented objects (belonging to a coherent family) is equivalent to having a set of 2D elastic displacement fields mapping a same (reference) ellipse into the different segmented objects.

Let us assume that the reference configuration is an ellipse E ⇢ R 2 ; a point x 2 ∂ E is mapped through a displacement field into a point x 2 ∂ F ⇢ R 2 , where F is the boundary of the segmented object. It holds :

(

Id + u : E ! F x 7 ! x = x + u(x ), (3.1)
and u 2 H1 (E).

From an acquired database 1 consisting of N 2 N ⇤ such elastic displacements : D = u (1) , . . . , u N , let us assume that a set of n 2 N ⇤ Proper Orthogonal Decomposition (POD) modes is extracted :

j i : E ! F, i = 1, . . . , n. (3.2)
When a new image to be segmented is available, the unknown elastic displacement mapping the reference configuration onto the image can be approximated by using the POD basis as follows :

u ⇡ u n = n  i=1 a i j i . (3.3)
The goal is to find the coefficients a i in such a way that the reference ellipse, mapped by the displacement field Id + u n is as close as possible to the actual image, hence providing an already good guess for the segmentation. In order to identify the coefficients a i a semi-automated procedure is devised. The image does not correspond, in general, to the reference configuration (this would be the case in which a i = 0, 8i). A number of m 2 N ⇤ points x ( j) 2 ∂ F is selected in which there is a discrepancy between the actual image and the reference configuration. Each of these points corresponds to an equation for the displacement field : there exists a x ( j) ⇤ 2 ∂ E such that :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m. (3.4)
This is a non-linear system of equations in the variables a i , x ( j) ⇤ ; however, if the points x ( j) ⇤ are known, this system is a linear system of equations for the coefficients a i . In the present work, the points x ( j) ⇤ are chosen as the ones minimising the distance to the point x ( j) 2 :

x ( j) ⇤ = arg min x 2E kx ( j) x k 2 `2 , j = 1, . . . , m. (3.5)
This corresponds to a regularisation for the solution of the non-linear system Eq.(3.4). The linear system to be solved reads :

Fa = d, (3.6)
F jik = j i,k (x ( j) ⇤ ), (3.7)
d j,k = x ( j) k x ( j) ⇤k , (3.8)
where F is a tensor, assumed to be of maximal rank (remember that it is constructed from POD vectorial modes). Of course, as soon as m n, the system is solved in a least-square sense, and admits a unique solution. When discretizing and solving this system, the following choices are made. The minimisation step described by Eq.(3.5) is solved by using a nearest neighbor algorithm that, among the points of the discretised ∂ E selects the one which is closest to a given x ( j) . Then, the vector system in Eq.(3.6) is unfolded and solved by using a rank revealing QR decomposition.

Numerical experiments in 2D 3.3.1 Reference configuration

In this section, some numerical experiments are shown on a simple 2D test case to assess the performances of the proposed semi-automated segmentation approach. The reference configuration, depicted in Figure 3.3, is a unit circle, divided into two semicircles.

FIGURE 3.3 -Reference configuration

2. Remark that this point might not be unique.

Generation of the database set

The images are defined by deforming the reference configuration by a parameter dependant field, generated as follows. The displacement of each point on the reference circle is normal to the circle and constructed by summing the set of N h 2 N ⇤ harmonics (eigenfunctions of the Laplace-Beltrami operator), scaled with a parameter dependant coefficient : let the undeformed configuration points be denoted by x 2 R 2 , r = kx k 2  1, J 2 [0, 2p] (i.e. x 1 = r cos(J ), x 2 = r sin(J ), the unit vector in the radial direction be denoted by êr and the displacement providing the deformed image configuration be x 2 R 2 :

x = X(x ) = x + r N h  j=1 g j ( j + 1) b ! êr , (3.9)
where g j = g j (µ), have an amplitude which decreases as a power law of the harmonics numbers :

( j +

Singular Values Decrease

In the present work K = 100 harmonics are used to generate the displacement fields. A set of N = 100 uniformly distributed in [0, 1] 2 samples of the parameter µ are taken, and, for each, the radial displacement is constructed. The POD basis was computed. The singular value decay is shown in Figure 3.6, exhibiting a faster decay when compared to the one typical of the spherical harmonics for a given regularity.

FIGURE 3.6 -Decrease of singular values with POD parametrisation

Remark that the displacement is purely radial and that the center of the reference circle r = 0 does not move. This corresponds to the case where the various anatomic regions that have been segmented are « aligned » and properly scaled. We are aware that in practice this point is crucial to achieve good outcomes and is not developed in this chapter, but in the next one. Furthermore, given this way of building a database of elastic displacements, the POD modes reduce to a set of vector displacements in the radial direction whose intensity is given by the value of the harmonics. The discretised shape is built by discretising uniformly the interval J using N J = 512 points.

Application and results

The image (the deformation of the reference shape) is cut by N p 2 N ⇤ horizontal planes, represented by an green dotted line in Figure 3.7, uniformly distributed on the image, to mimic a medical exam like CT or MRI which are sliced. Each plane intersects the only outer boundary of the object in two points (m = 2N p ), each of which corresponding to one of the halves of the deformed circle. These points, represented by an orange star in Figure 3.7, provide an equation for the displacement field to be approximated by using the POD modes. Several numerical experiments are made to assess the performances of the proposed reconstruction method, by varying the number of modes and the number of cutting planes used. In particular, the number of modes used varied from n = 6 to n = 20 and, for each choice the number of planes was varied from the minimal number of planes needed to have a square system n  m to m = 2n. In Figure 3.8 twos example of reconstruction is proposed, when the images were produced by using N h = 100 harmonics, n = 20, m = 80. The result is quite precise, all the details of the image are well represented, whether it is the outer edge where points are selected or the edge of the inner hemicycles where no points are selected.

Reconstruction error

The reconstruction errors are shown in Figure 3.9 for a population of 100 target displacements (that do not belong to the database used to construct the POD modes). The errors are measured in `2 norm with respect to the discretized images. The average error and standard deviation are shown when the minimal number of planes and the maximum number of planes are used. In general, the errors are compatible with the ones considered as compatible with realistic applications. When the minimum number of planes is used a stability issue can be seen when a large number of modes is used and the standard deviation is quite large. On the contrary, when N p = 2n (orange), the lack of stability is delayed, the errors decrease and the standard deviation is much smaller, reaching a minimum for n = 17.

First 3D numerical exploration

In order to better understand the complexity of the situation in 3D, we chose to also work on a simulated problem. We therefore generated a set of shapes using spherical harmonics to represent our dataset. Then we extracted the modes with a POD, and reconstructed these shapes according to the basic methods explained in the previous section.

Spherical harmonics -Definition in 3D

Spherical harmonics are special functions defined over the unit sphere. They can be derived as eigenfunctions of the Laplace-Beltrami operator on the sphere. They can be particularly useful to solve rotation invariant problems. They are a complete set of orthonormal functions of L 2 (S 2 ), and thus may be used to represent functions defined on the surface of a sphere, just as circular functions (sinuses and cosines) are used to represent functions on a circle via Fourier series.

The harmonic polynomials P(x, y, z) of degree l form a vector space of dimension 2l + 1, and can be expressed in spherical coordinates (r, q , f ) :

r l Y l,m (q , f ), with l  m  +l, where q 2 [0, p] and f 2 [0, 2p].
The cartisian coordinates of a point can be retrieved by the change of coordinates :

8 < : x = r sin f cos q y = r sin f sin q z = r cos q
In the present work, we will consider the unit sphere and hence r = 1.

A real basis of spherical harmonics can be defined in terms of their complex analogues by setting :

Y l,m = 8 > < > : p 2( 1) m Im[Y |m| l ] if m < 0, Y 0 l if m = 0, p 2

Parameterization of shapes

In order to test the proposed method, we define two families of shapes, constructed by using spherical harmonics, in a similary way to the 2D method presented in section 3.3.2.

Let x 2 S 2 be a point of the unit sphere, which is the reference (undeformed) configuration. A shape is defined as the result of a mapping of the unit sphere. Let x 2 F ⇢ R 3 be a point of the mapped surface. It holds :

( Id + u : S 2 ! F x 7 ! x = x + u(x ), (3.11)
where u is the displacement field, defined as follows :

u = L  l=1 l  m= l a R lm Y R lm + a I lm Y I lm Let L 2 N.
For the first familly, the coefficients are defined as follows :

8 < : a R lm = g R lm (1+l) b a I lm = g I lm (1+l) b (3.12)
with g R,I lm are sample of a uniformly distributed random variable in [ 1 2 , 1 2 ], and b the regularisation .

R The shape family that is constructed here is considered large enough, because all the generated coefficients are independent.

An exemple of a shape generated by this paramatrization with a great regularisation (b = 2), is shown in 3.11.

FIGURE 3.11 -Example of shape generated with spherical harmonic functions, from 2 point of view

Aiming at studying shapes featured by common patterns (as one would expect from anatomical parts), we introduce a second family of shapes, in which the coefficients defining the displacement are no longer independent. Instead, they depends upon a small number of parameters.

Let µ = [µ 1 , µ 2 , µ 3 ] 2 [0, 1] 3 .
The coefficients of the real part and the imaginary part were chosen as below : An exemple of a shape generated by this second paramatrization (b = 1.5) is shown in Figure 3.12. FIGURE 3.12 -Exemple of shape generate with a parametrized spherical haromonic function Once those two types of shape generated, we perform a Proper Orthogonal Decomposition (POD), by computing the matrix C i j , its eigenvalues l i and its eigenvectors v (i) j :

C i j = 1 N N Â k=1 N Â h=1 u i h,k u j h,k
with u i h,k being the displasment at point i between the initial sphere and the shape generated by the spherical harmonics.

Thus we compare the decrease of the eigenvalues of the two types of shape in Figure 3.13, in blue the parameterized one and in black the non-parameterized. We can see that the blue curve decreases quicker at first and even during the plate phase, because it represents a parametric family of shapes with a small Kolmogorov dimension (see Chapter 2), just as it is the case of our example on masseter muscles. FIGURE 3.13 -Exemple of shape generate with spherical haromonics function

Reconstruction of a new shape

Thus, from the acquired database consisting of N 2 N ⇤ such displacements : D = u (1) , . . . , u (N) let us assume that a set of n 2 N ⇤ Proper Orthogonal Decomposition (POD) modes is extracted :

j i = N Â j=1 v (i) j u ( j)
We are now able to approximate, by using the POD basis, the unknown displacement mapping from the sphere reference configuration onto the object F as in the 3.3 equation :

u ⇡ u n = n  i=1 a i j i
with n << N, the maximal number of objects consider in the database, here N = 100.

The goal is now to find the coefficients a i , in such a way that the reference sphere, mapped by the displacement field Id + u n is as close as possible to the actual object. In order to identify the coefficients a i , we used the method described in the previous section. Where a number of m 2 N ⇤ , m n points x ( j) 2 F is randomly selected in the actual object F, which allows us to build a linear system that is easy to solve by using a rank revealing QR decomposition.

Analytic error

To simulate the segmentation situation and have a better idea of what can be expected from an error point of view, we performed an analysis of the error of the reconstructed shape, and compared it with the error of projection.

The simulation was performed with these specific parameters :

-for the generation of the hamornic spherical basis, l = 10 which means 66 modes were used to construct N = 100 random cases with the parametrisation of the section before, -for the reconstruction of those cases, differents types of reconstructions were tested with L = 3, 4, 5 and different regularity b = 0.5, 1, 1.5, we compare the projection error with the reconstruction error with different numbers of selected points : 40, 80, 120, 160, 240, 320, 480, -the reference sphere had 20002 points.

We obtain the results presented in Table 3.1 for the projection error and in Table 3 Comparing both these tables, we can conclude that for a regular enough shape (b 1.0), and mostly thanks to the parameterization introduced in the previous section, our reconstruction error is quite similar to the projection one using 200 randomly selected points, which represent 1% of the data knowledge. In the present work a semi-automated reconstruction method is proposed to help the segmentation of anatomical details in the cases in which the contrast is not enough to use fully automated techniques. The method is based on the use of a reduced basis to approximate an elastic displacement mapping a reference configuration into the actual shapes of the anatomical details of interest. The method uses a nearest neighbors optimisation procedure coupled to a least square approach to solve the problem. A 2D and a 3D test cases are proposed to mimic realistic applications and assess the performances of the method. A set of synthetic generated images was used and the error of the numerical experiments showed that the proposed approach can lead to a viable tool. In this respect, further investigations are mandatory both from the theoretical and applied point of view. In particular a systematic error analysis will be carried out and more realistic tests will be performed, in 3D configuration and by taking the noise and the other experimental error sources into account.

Especially in 3D, several challenges are added to those previously described, especially in the case of a real application. The first one is the constitution of the manual segmentation database : time-consuming and requiring the involvement of an anatomist expert at first. However, the idea is that we will soon be able to propose a first phantom, which will require only a few corrections and will save a lot of time while continuing to feed the database. All these aspects are discussed in detail in the following chapter.

Introduction

In medical imaging, segmentation of anatomical parts is still a chalenge. Indeed, it is a recurring and time consuming problem, as seen in the previous chapter.

Segmentation is the process of dividing a medical image into distinct anatomical regions [68]. For 3D medical images, this amounts to attribute a label to each voxel (3D pixel) of the image. The gold standards have been, for years, manual segmentations performed by anatomical experts, but on top of the unavoidable errors, it is a time consuming process. The general issue is thus to find a way to identify automatically an anatomical region of interest in order to speed up the process, and eventually standardize as much as possible the generated error. The different existing methods and the associated literature are presented in more detail in the previous chapter.

We propose here to develop a model based semi-automatic method to help segmentation, according to the same general principle as previously exposed in the chapter 3. In particular, the shape of the object of interest is represented by an elastic deformation field of a reference configuration. A reduced basis representation of the object of interest is introduced, that exploits an a posteriori knowledge of a database of available segmentations. Given an image to be segmented, the elastic deformation field can be approximated as a linear combination of a small number of modes. This provides, in turn, a powerful regularisation and a significant speed-up of the procedure needs since few apropriate points need to be selected by the user in order to define the field.

This problem and the associated method were first set in last chapter where we proposed a proof of concept. We present here the same problem, extended and addapted in 3D, with a practical application of the masseter muscle segmentation on CT-Scan. However, due to various problems encountered in the 3D generalisation of the techniques, for a realistic application, used for the very simple previous examples, some steps of our segmentation algorithm are significantly modified.

The chapter is structured as follows : a first part about the modes generation according 2 methods (POD and Greedy algorithm). In a second part, we present an overview of the new method adapted for a 3D realistic problem. Finally, a third part which present the application of the latter for the particular case of the masseter muscle on CT-Scan.

Modes in the reference configuration

The first step of our method involves, as referred on several occasions, computing a reduced basis in which all the elements of our dataset are represented by only a few parameters. All the stages of the constitution of our database are explained in Chapter 2. The construction of the modes can be summarized as follows.

By exploiting the morphing of a reference shape in 3D, the ellipsoid presented in the Chapter 2, it was possible to define the elastic displacements u (Figure 2.10) and from them, the modes.

Two extraction methods of the modes are possible in our application : 1. POD modes, quick to implemented and computed. These are the ones we used for the proof of concept and for the numerical exploration in the previous chapter. 2. Greedy modes, were all masseters are represented by a linear combinaison of the ones of the database. A specificity is to be noted here : we obtained these modes by taking into account that for our application, the knowledge of the masseters shape is not uniform on the image. The points that will be selected (detail in the step 1 of the methods), due to the nature of the image, won't be uniformly distributed on the image, but restricted to a certain region. Which can present an important advantage for the reconstruction part, even if that kind of algorythm is slower at first.

A set of n 2 N ⇤ modes is retained. For the present application, we choose to use n = 20. And each masseter can now be represented as :

n  i=1 a i j i . (4.1)
with a i the coefficients and j i the modes POD or compute with our greedy algorithm.

It is important to specify that due to an artificial translation present for some snapshots in the database, the number of couples (image-segmentation) allowing the implementation of the whole process (including the reconstruction with projection on the image, i.e. the visualisation of the segmentation) is here reduced to 45 masseters (against 126 in our complete database).

In this context, we were all the more confronted with a weakness of this type of model : the need for a consistent database. To overcome this problem, we have implemented an automatic shape generator tool, which from existing shapes and with the help of interpolation and spherical harmonics, allows us to multiply the dataset. We also add in this process small rotation and translation to favour a better registration. Due to the fact that this part of the method is still under development, details of its design and results will not be presented in the following section.

Overview of the method for 3D segmentation in a real case

As a first step, and in order to clarify the structure of the segmentation algorithm developed here and the modifications made in respect of Chapter 3, we begin in this section by presenting an overview of the different steps and their sequence.

Notation

In what follows : -N p 2 N ⇤ : number of points selected by the user on the image, -N s 2 N ⇤ : number of snapshots in the database, -N 2 N ⇤ : number of points of the reconstructed mesh, n 2 N ⇤ : number of modes.

Pipeline

The pipeline we have set up is directly linked to the segmentation application : patient-specific geometric transformations must be derived solely from the Dicom image, and geometric transformations that are not linked to the image are common to all patients. This ensures that, for a new patient, we can produce a muscle reconstruction without the need for any additional information (as is the case in realistic applications).

Let the image domain be the subset W ⇢ R 3 . The image is considered as a function F :

F : ( W ! R x 7 ! F (x 1 , x 2 , x 3 ). (4.2)
We therefore present here the different steps necessary for the semi-automatic segmentation of a new masseter muscle, from the acquisition of the dicom to the projection of the final segmentation on the dicom images or in 3 dimensions.

-Step 1 -An interface that allows the operator to select a set of points describing the muscle and the associated bounding box.

A bounding box W m ⇢ W is determined on the image, such that the muscle is contained in it. Moreover, let N p 2 N ⇤ , a set of points p (i) 2 W m are chosen, i = 1, ..., N p .

The most important thing in this step is to make sure that the whole muscle is contained in the box, or more precisely that all the selected points are contained in the bounding box.

-Step 2 -A coordinate change from dicom images to our reduced basis space [0, 1] 3 .

An affine transformation is built, such that :

X : ( W m ! [0, 1] 3 x 7 ! x . (4.3)
We have x = X (x), and, since the transformation is invertible, x = X 1 (x ).

The points are transported into [0, 1] 3 , that we call the reference configuration, where the modes are defined and where the reconstruction takes place.

-Step 3 -The reconstruction of the 3D masseter from the selected points, i.e. the identification of the best coefficient a to fit those points in our reduced basis.

By exploiting the points, a mesh consisting of N 2 N ⇤ degrees of freedom is reconstructed :

v (i) , i = 1, ..., N .
These two last step differ significantly from the method described in Chapter 3. All the details and needed examples are given in the next section.

-Step 4 -Visualisation and validation of our segmentation on the dicom images.

The mesh is transported into the image configuration by using u (i) = X 1 (v (i) ) and a phantom is computed based on them, which is : F :

( W ! R x 7 ! F (x 1 , x 2 , x 3 ). (4.4)
This allows the user to check the segmentation and modify it if it is not fitting.

All of these steps, briefly named here, are detailed in the following section.

4.4 Application in a real case : the masseter segmentation

Graphical interface development

All the interface that is presented in this part has been implemented during this thesis, using python language and the package « OpenCv ». the masseter. The main instruction given to the user is to make sure that the whole muscle is located inside the selected box, even if it seems a little large at first.

In practice it is necessary to select two rectangles : one on the axial section and one on the frontal section, as can be seen in the example in the Collecting this data for all 45 Dicoms-Segmention couples in our database was the first step in setting up our reconstruction method.

Coordinates Change

Once the points collected, we need to send them to our reference configuration thanks to the bounding box selected at the same time, and the reversible application X .

In Figure 4.4, we can visualize the point cloud selected in the image reference, on the left, and the masseter we want to reconstruct in the unit box, reference configuration, where the modes was computed, on the right. When a point is selected on a dicom image, a local value (i, j, k) based on the pixel numbering is attributed to it. We need to change this coordinates to the cartesian system x, y, z, to be able to perform our analysis. This is made possible because the slice thickness information, voxel spacing, and number of slices s are recorded in the dicom files.

In Figure 4.5, we can visualize a schema of a head that summarizes the orientation of the axes x, y, z according to the image coordinates system. The Cartesian coordinates of the points are deduced from this, according to the section plane (axial, coronal or sagittal) in which they have been selected. Once the Cartesian coordinates have been computed, we must align the acquired data with those in our database, i.e. the space where the POD was performed. This registration step has a significant impact on the results, and should not be underestimated.

As the reference point of the initial bounding box, the point G (x min , y min , z min ) is chosen, where x min is the minimum x-coordinate of the selected points in the different sections, projected into the Cartesian coordinate system. It is from this point G that the box is centred in zero : 8 > < > : that we apply to all the points of the box, such as :

x t = x x G y t = y y G z t = z z G .
X a = hC, X t i L 2 , (4.7)
with X t vector (x t , y t , z t ) and X c vector (x , y c , z c ).

To ensure that all the transported masseters of our database are in our reference space, the unit box, [0.1] 3 , we applied a discount factor l identical to all the masseters, when we construct modes (Section 4.3). We therefore apply the same reduction coefficient to the selected points on the image :

X = l X c + 1 l 2 , (4.8)
with in our case l = 0.7, the highest coefficient for which the 45 masseters in our database are contained in the unit box.

The result of this registration procedure can be observe in the Figure 4.6, where the points cloud is now almost superimposed on the masseter mesh. This simple 2D illustration does not allow for easy visualization of the distance of the points (in red) to the masseter (in blue) in 3D. However, if we display the rotation of this figure in space (on a computer screen), we can see that all the points are close to the original object.

Reconstruction

To be able to perfom our reconstruction and as suggested in a previous section, we need to identified the points x ( j) ⇤ of our reference ellipsoid E by minimizing the distance to the selected point. Here we presented a potential solution to this problem.

The novelty comes in this reconstruction phase. In the previous chapter, given a point q (i) = X (p (i) ), the first step consisted in assessing which was the degree of freedom of the ellipsoid that was more likely to be the one we need to deform to get close to the point. After having perfomed this task, the problem becomes linear. The problem is that the result could be strongly affected by the correspondence point-dof (that takes care of the non-linearity).

In the present approach, we try to perform the reconstruction by solving the non-linear problem. To this end, a functional is introduced to measure the discrepancy between the points to be matched and the reconstructed mesh.

Let v j 2 [0, 1] 3 , j = 1, ..., N , and v ( j) = x ( j)

`+ Â n i=1 j i j a i be the j-th reconstructed mesh point. Let Q = {q (1) , ..., q N p } be the selected points cloud and the reconstructed mesh be V = {v (1) , ..., v N }.

Then :

J : Q ⇥V ! R + , (4.9) J = 1 2 N p  i=1 min j21,....,N |q (i) v ( j) | 2 . (4.10)
The meaning is : when we have to evaluate the distance from the point q (i) we consider the mesh point which is the closest to it, we compute the Euclidean distance squared, and then we sum over all the contributions.

The points v ( j) are function of the coefficients a 2 R n . The reconstruction amounts to solve :

a ⇤ = arg inf a2R n J (a). (4.11)
This problem is known in the literature as min-sum-min because of the structure of the cost functional. To solve this problem we have been inspired by the method developed by A. Demianov in [START_REF] Alexei | On the solution of min-sum-min problems[END_REF]. However, it was necessary to adapt this algorithm to our specific case. Indeed, the examples presented in the article remain quite simple and in 2D, therefore less difficult and expensive to implement. We have thus initially chosen to leave aside the combinatorial part which is exposed in this paper. The present approach to solve our min-sum-min problem reads as follows :

1. The initial guess is found by looking, among the elements in the database, to the one which is the closest to the points cloud. Let A = {a (1) , ..., a (N s ) } be the set of the coefficients resulting from the best reconstruction of the database snapshots, in the above introduced sense (minimising J ). 2. When solving for it, a correspondence « point-dof » is found, namely :

j (i) ⇤ = arg inf j2{N } |q (i) v ( j) | 2 . (4.13)
That makes it possible to define the set J = { j(i), i = 1, ..., N}. The set V | J , whose cardinality is V | J = N p is the restriction of the mesh to the points closest to Q.

3. Then we se the set V | J to get rid of spurious translations and linear transformations (intermediate registration step). Update the sets J and V J .

4.

A least square problem is solved :

a ⇤ = arg inf a ||q x ` F| J a|| 2 l 2,Np . (4.14)
5. And finally recompute J, update V J and repeat from step 4 of this procedure.

In step 5, the demonstration of convergence in this type of min-sum-min problem is not trivial at this stage of the method's development. It would therefore be necessary to define more precisely a stop criterion. Furthermore, it is also possible to add a translation/rotation stage between these different steps to optimise alignment. However, it is also necessary to determine when this stage would be the most judicious and whether it is more interesting to repeat it several times or not.

Visualisation

Once all these steps have been completed, the central operations of the algorithm are carried out. What we want to do now is to show to the user the result of its semi-automatic segmentation. In order to visualise our reconstructions in 3D or to project on the dicom images, we need to reverse the registration process that sent them to the reference configuration. This is easily achieved by reversing the transformation used in section 4. It is important to keep in mind that the manual segmentation obtained here (in yellow, on the left of Figure 4.7), is possibly not perfect, as it was performed by a single operator with a developed knowledge of mandibular anatomy but not of segmentation. The differences observed are therefore difficult to interpret with this 3D representation. Indeed, this 3D visualization does not allow us to judge the quality of the underlying segmentation. This is why we present in Figure 4.8, the projection of this reconstruction directly on the image. Here we can visualise the boundary of the edge of the masseter segmentation on an axial section, by the manual method (in red), our reconstruction using the POD modes (in green), and the one with the greedy algorithm modes (in blue). FIGURE 4.8 -Projection of the reconstruction on the original dicom images : on the left our recontruction with the POD modes, on the middle the manual segmentation, on the right our reconstruciton with the greedy modes. (Please note that the quality of the images presented above is not optimal, this is not related to the results themselves, but more to a lack of control of the viewing software. However, they do allow us to observe a first result which is a good basis for discussing the method).

One of the main comments that can be made from these first results is that our method, whether it is carried out with the POD or Greedy method, achieves one of our main objectives : to unravel a first segmentation from little information ask the operator (50 points and a bounding box). Moreover, the Figure 4.8 shows that the manual segmentation extends slightly over the bone, which is not the case for our reconstructions. This can be explained by different reasons.

First of all, this is a very specific example : a single axial section for a single masseter in our database. We therefore cannot generalize to all cases, but this allows us all the same to make some hypotheses to understand how this could have happened even occasionally.

The operator who performed this manual segmentation probably did not make such a mistake when selecting regions of interest in the slices, but more likely later when using 3D regularization tools that allow to obtain a more coherent 3D shape from an anatomical point of view. This therefore amounts to saying that these are the stages of post-processing of the 3D object which would induce an error in this tedious and time-consuming segmentation.

In the case of our reconstructions, the method we have developed forces the final shape to stay close to the selected points. However, as explained previously, the user who selected the point cloud given as input to our algorithm, clearly visualizes the change in density between the muscle and the bone, and has very probably selected points quite precisely at this level. This therefore allows us not to make any mistakes at this level.

Anyway and in order to be able to discuss this point in more detail, it is important in the rest of this work to do a more in-depth analysis of the error that was made during our reconstruction. For this, additional work to obtain the gold standard of segmentation at the present time : manual segmentations of different experienced operators for the same masseter, is necessary. This will allow us a more objective comparison using methods currently available in the literature, such as the « dice coefficient », in order to further improve our method still under development.

Discussion and Conclusion

This method, which is still under development, meets several of our initial criteria, better than the method initially proposed in Chapter 3. However, it is still under development and there is still much that needs to be improved. These include the following :

1. Regularisation of least squares : The coefficients a 2 R n are known for all the elements of the database. A possibility would be to constrain a to belong to a certain sub-set [a (min) 1

, a

(max) 1 ] x ... x [a (min) n , a (max) n
] for a new reconstruction.

2. The discrepancy measure for the min-sum-min problem seems to be reasonable. Are there other ways, with easier functionals to be minimised ? A test was done by using exponential kernels adopted in computational geometry, but the result was depending from the kernel length (free parameter) in a strong way. A more detailed study of the literature or the use of this parameter could improve this point.

3. Enrich the database by considering random expansions based on Laplace harmonics to perturb the reference ellipsoid (we can ensure the fields to be still in H 1/2 (d |W `)), as suggested in the section 4.2. This would add a bit of noise to the basis but would also allow it to be increased considerably, which would be particularly useful when the base is really small. It is also possible to account for rigid displacements using this method by adding a small random translation and/or rotation to the deformation. This would improve the alignment and correspondence of the points, between those selected on the images and those of the mesh of the reconstructed masseter.

4. Can we devise a procedure to work directly with the phantom ? This would make it possible to show the overall improvement of the reconstruction in real time at each selected point, and would also potentially save computing time.

5. The ultimate goal of this work is to be able to completely automate the process. This saves time and hopefully reduces human error. Some ideas are being considered to get rid of the manual selection of points.

Introduction

3D reconstructions of the facial skeleton obtained from a CT-scan are widely used in clinical practice to analyze facial morphology disorders or lesions of the facial skeleton. More recently, specific digital tools were developed to perform surgical planning and simulation of the morphological result. Any surgical treatment can be planned based on available preoperative diagnostic data [START_REF] Tahir | A shape prior-based MRF model for 3D masseter muscle segmentation[END_REF][START_REF] Tahir | Using a flexibility constrained 3D statistical shape model for robust MRF-based segmentation[END_REF]. To obtain such virtual model, the hard tissue -bone tissue -represented in the CT scan must be identified to generate the corresponding 3D object through a segmentation process. Bone segmentation is quite easy to perform because its density is radically different from that of adjacent tissues.

However, when a more comprehensive model is needed, soft tissues should also be segmented, such as for a dimentional study of a muscle (cf. Chapter 1) or a facial reconstruction methods (cf. Chapter 7). This implies the skin surface and some of the facial muscles that can impact facial morphology outcome. The masseters are the muscles contributing the most to facial morphology at rest (cf. Chapter 1). Other facial muscles mostly determine facial movements but have no significant contribution to facial volumes. The segmentation of the skin surface is also easy to obtain. However, the segmentation of the muscles from CT-scan images is more challenging as the surrounding tissues have similar density with blurred margins between them. The task is even more arduous when the CT-scan contains artifacts, noise or inhomogeneities. This prevents the automatisation of segmentation. Yet the current clinical practice -manual segmentation [START_REF] Tahir | A shape prior-based MRF model for 3D masseter muscle segmentation[END_REF] requires tedious visual identification and manual highlighting of muscles. The process is very time-consuming and subject to operator variability. Accuracy is limited, and it currently requires the input of expert operators [START_REF] Tahir | A shape prior-based MRF model for 3D masseter muscle segmentation[END_REF]. A possible improvement can involve automatisation of this process for segmenting target tissues [START_REF] Yousef | The effect of labeled/unlabeled prior information for masseter segmentation[END_REF].

The purpose of this study is to validate the semi-automatic segmentation method, introduced in the two previous chapters, by comparing it to manual segmentation, for clinicien uses.

Materiels and Methods

Medical Image Data

CT-scan data were collected from a Maxillo-facial surgery Department database. All included patients were informed their data would be collected anonymously for scientific purpose and that they could refuse the scientific use of their data. 26 patients (19 female patients/7 male patients) aging from 16 to 47 years were included.

The selection criteria we considered were : -Preoperative scans exhibiting no bone implants to avoid artifacts.

-No edentulous patients.

-A field of acquisition including the whole masseter muscles.

All CT-scans were obtained using a Discovery CT750HD system (GE Healthcare, Chicago, Illinois, United-States). Slice thickness was 1.25 mm and the pixel spacing was (0.46 mm, 0.46 mm).

Segmentation process

We call here segmentation the process of attributing a label to each voxel (3D pixel) of an image to generate a corresponding 3D object.

A first case was used as training by the two operators with the two methods, to get used to the softwares. After this first training, the segmentation of 40 masseters from 26 CT-scans, was performed. Each operator processed 20 masters segmentation in a seperate database. We needed to limit ourself to this number of segmentation due to the time required for the manual segmentation. First, manual segmentations were done, then the semi-automatic one was performed by the same operator on the same data. The time required for the two segmentations was recorded for comparison.

Manual segmentation

The manual segmentation was performed using Amira R [START_REF] Detlev | Amira : A highly interactive system for visual data analysis[END_REF] (develloped at the ZIB, Berlin, version 2019.05) as shown in Figure 5.1. Each CT-scan was reconstructed in coronal, sagittal and axial planes. Manual segmentation was obtained by delineating the contours of the masseter muscle in axial and coronal views, on every slices. Various selection tools offered by the software have been used for this purpose, such as « region growing », or « selective brusheing ». In a final step, the continuity and relevance of the traces were re-verified. To identification of the cranial and caudal limits of the masseter muscle, and so to define the upper and lower limits, we used the coronal slices. The cranial limit of the masseter was marked at contact with the zygomatic arch. Then double checking was made between the coronal and axial sections in order to check the coherence of the segmentation.

After this first step of segmentation, a specific protocol was used in Amira R to exclude irregularities and inconsistencies in the 3D objects.

The protocol included the following actions : -« Remove Island » : to delete isolated parts of less than 15 voxels.

-« Fill Holes » : this function allows to automatically incorporate all the voxels situated in the muscle body and which have not been selected by the operator. -« Smooth Label » : smooth the 3D general contour, with a tolerance set to 3.

At last, the segmentation was double-checked by a third expert.

In Table 5.1, we show the time of segmentation by operator for the two methods.

Operator Semi-automatic segmentation was very significantly faster than manual, (p=5.61e-23). The inter-operator time was significantly different with the manual method p = 6.67e-5, whereas with the semi-automatic one no difference were found p = 0.07.

The descriptives statistics of the error between the manual and semi-automatic methods are presented for each operator in the Table 5 The inter-operator error (computed for the semi-automatic segmentation) were found not significantly different between the two operators, with p = 0.79.

Discussion

We have presented in this paper a proof-of-concept study focused on the masseter muscle from CT-scan images. The masseter is one example among many structures that can be segmented using a semi-automatic process and the type of image can vary also. Semi-automatic segmentation can be of use in highly diverse clinical applications/settings.

FIGURE 5.3 -1-Artifacts, 2-Masseter layout hardly discernible due to artifacts

The general process of the segmentation reveals some major difficulties : First, images may contain artifacts (Figure 5.3), often caused by implants or dental crowns and fillings, which can change the intensity of a voxel so that its classification becomes uncertain [START_REF] Kathrin | Comparison between manual and semi-automatic segmentation of nasal cavity and paranasal sinuses from CT images[END_REF].

Secondly, because of the blurred boundaries between the masseter and its adjacent anatomical structures, it was sometimes difficult, even for experts, to determine with precision where the muscle boundaries are.

Consequently, it is difficult to achieve a perfect segmentation, due to the occurrence of undersegmentation or over-segmentation. This means, in the first case, pixels belonging to different objects are classified as belonging to the same object ; and on the contrary, in the second case, voxels belonging to the same object are classified as being of different objects [START_REF] Benjamin | Feature-based similarity search in 3D object databases[END_REF][START_REF] Tahir | A shape prior-based MRF model for 3D masseter muscle segmentation[END_REF][START_REF] Tahir | Using a flexibility constrained 3D statistical shape model for robust MRF-based segmentation[END_REF][START_REF] Kathrin | Comparison between manual and semi-automatic segmentation of nasal cavity and paranasal sinuses from CT images[END_REF].

We encountered several difficulties in manual segmentation, among which : -Lateral and medial limits were very difficult to identify on the coronal and sagittal slices, as such we mostly used the axial slice (Figure 5.4). Many techniques proposed in the literature are semi-automated or even fully automated, leading to highly consistent segmentations across multiple users [START_REF] John | The semiotics of medical image Segmentation[END_REF], but not on every type of anatomical objects or medical images. In terms of medical ethics and law, the responsibility for errors of diagnosis may lie with the medical practitioner and, therefore, the correction of the automatic segmentation errors is required to avoid negligence [START_REF] Sara | Semi-automatic segmentation for prostate interventions[END_REF], hence the interest of semi-automatic segmentation.

The impossibility to visualize in real time the 3D shape of the masseters as in manual segmentation has been the main difficulty encountered during the semi-automatic segmentation. Indeed the operator often verify their slice segmentation by checking the 3D anatomy. This can also be considered as an improvement as it decreases the operator subjectivity in the segmentation process. The anatomy information being provided by the database with the semi-automatic method and not by the operator anymore.

It has to be noticed, however, that the 2 operators had segmented the masseters manually at first. It made the semi-automatic technique, on the same muscle, easier even if a delay was respected between the 2 segmentations. Nevertheless, as mentioned in the methods section, this method requires a first database of the muscle of interest to be able to work, and this database is mostly done manually.

This study shows a significant gain of time when using the semi-automatic method. We also evidenced that the semi-automatic technique leads to better inter-operator « reproducibility ». The errors made by the semi-automatic segmentation tools developed is in average 2.85 mm (standard deviation 0.9 mm). It is important to keep in mind that, this method is still under development and some improvement can be expected, for the registration time, the error and the reproductibility study, which is really important in our case since the shape approximation depends on the amount and location of the points selected.

The semi-automatic technique uses a 3D muscle database and a model-based computer process with the goal to decrease the time of segmentation, with similar precision, as well as making segmentation more accessible.

Finally, it must be kept in mind that the segmentation times are relative to the operator's level of experience. A well-trained operator, with a better knowledge of segmentation processes, would probably have reduced segmentation times. However, it can be reasonably assumed that this would probably be the case for both methods presented. The semi-automatic method would then keep its advantage.

Conclusion

This chapter presents a new, simpler and less time-consuming technique for segmenting soft tissue from CT data sets, than the usual manual segmentation. Muscle segmentation is a difficult task due to its resemblance to surrounding soft tissue. However, this new method shows that by incorporating prior knowledge of the shape (with a database and a compute model method), clinically acceptable results can be obtained through a semi-automatic segmentation software.

Since this method is still under development, it would be interesting to test a number of additional tools to determine how the result would be affected. We can cite : -a complete inter-and intra-operator reproducibility study, -a guest phatom, on 2D images of the anatomic object, as the points selection is made for the semi-automatic method, a 3D visualization, of the anatomic object, as the points are selected for the semi-automatic method, -the option of manually editing the anatomical object obtained semi-automatically, to allow the enrichment of a database under construction. This would speed up the initial process, while allowing to continue acquiring new shapes for the dataset.

INTRODUCTION PART III

This last part contains two chapters, about the primary objective of this thesis work : facial reconstruction and the morphing method used to obtain it.

In chapter 6, the morphing method is detailed, which constitutes one of the main elements of the facial reconstruction method detailed later on. This work is a follow-up to the one presented by Maya de Buhan and Chiara Nardoni in [START_REF] De | An optimization method for elastic shape matching[END_REF], and was carried out jointly with both of them. Chiara Nardoni has already presented it in her own thesis [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF]. In particular, my contribution to the development of the boundary method, its comparison with the volumic method, and 2D and 3D numerical simulations. This last chapter (chapter 7 of the thesis) is the outcome for the FaciLe project, dedicated to facial reconstruction for forensic medicine, and presents an improvement of the method proposed by Maya de Buhan and Chiara Nardoni in [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], taking into account the underlying anatomy of the face. All the process is fully automatic, which allowed us to propose an application to validate our method compared to the one mentioned above via the recognition from a picture of different reconstructions by a population of more than 100 people.

*the nice figures presented above are the results of unfortunate numerical simulations carried out during the work on Chapter 6 and 7, and they also illustrate some of the difficulties of numerical implementation work. criterion depending on the squared Euclidean distance to the boundary of the target shape. Then, a novel volume criterion depending on the signed distance function to the target shape is used. Both deformation models are driven by a computational mesh and have been implemented in a finite elements setting. We present numerical examples in two and three dimensions to compare their efficiency as well as their drawbacks.

Key Words

Shape optimization, Shape matching, Level-set function, Finite elements, Mesh deformation.

Introduction

This article is devoted to the problem of 2d and 3d shape matching (or shape morphing). The main issue is to find a transformation from a given 'template' shape W 0 into a 'target' shape W T . Such a transformation may be used to compare the two shapes W 0 and W T and find their similarities which is a key ingredient in shape retrieval, classification, recognition or shape tracking. One can also be interested in computing the physical deformation of W 0 into W T to reconstruct or simplify the target shape, or simply in transporting quantities of interest as textures from W 0 to W T .

In our approach, the desired transformation is realized as a sequence of elastic displacements, which are obtained by minimizing an energy functional based on the distance between W 0 and W T . Indeed, we introduce two functionals devoted to measure the discrepancy between shapes : the first one uses a boundary criterion that depends on the square of the Euclidean distance between the boundaries of the two shapes, whereas the second one is based on a volume criterion depending on the signed distance function (that will be defined in (6.8)) computed on the target shape. The latter has been already briefly presented in [START_REF] De | An optimization method for elastic shape matching[END_REF]. In the current article, we complete the theoretical framework, give more details about the numerical implementation, show new 2d and 3d examples and, above all, propose a comparison of the original volume criterion with the more intuitive boundary one to prove its superiority.

Our method has much in common with [START_REF] Konrad | Elasticity-based matching by minimising the symmetric difference of shapes[END_REF][START_REF] Konrad | A linear elastic force optimization model for shape matching[END_REF] as it proposes a physically inspired deformation model and also with [START_REF] Ruzena | Multiresolution Elastic Matching[END_REF], borrowing techniques from optimal control and shape optimization. For a general overview on shape matching methods, we refer to the survey [START_REF] Remco | Shape matching : Similarity measures and algorithms[END_REF]. Let us only mention here a few techniques.

In the field of computational anatomy, a series of articles (see e.g. [START_REF] Mirza | Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms[END_REF]67,[START_REF] Ulf | Computational anatomy : An emerging discipline[END_REF]) have suggested to describe the deformation of W 0 onto W T via a diffeomorphic map. The restraint to diffeomorphims prevents discontinuous or self intersecting deformations. This approach leads to the so-called Large deformation Diffeomorphic Metric Mapping (LDDMM) method which has received a lot of interest in the last years, especially for dealing with medical data. The LDDMM approach considers shapes embedded in a metric space and describes the diffeomorphism between two items W 0 and W T as the dynamic flow of a velocity field v. The search for v is then casted as an optimal control problem. The criterion to be minimized is the sum of a matching functional -representing the cost of the flow which transports W 0 onto W T -and a regularization term, guaranteeing the smoothness of the deformation. The structure of metric space allows the definition of similarity measures between shapes, in terms of geodesics in the shape space. The LDDMM formulation produces a global deformation of the ambient space, which in particular matches the shapes of interest into each other. In [START_REF] Kamal | Variational shape matching for shape classification and retrieval[END_REF] the energy functional associated to a LDDMM problem is used for defining a measure of similarity between two shapes, defined as the geodesic path in the space of diffeomorphism by minimizing an energy functional. The LDDMM framework has been used for landmarks-based matching using the spline interpolation theory [START_REF] Sarang | Landmark matching via large deformation diffeomorphisms[END_REF], unlabelled point-sets [START_REF] Joan | Diffeomorphic matching of distributions : A new approach for unlabelled point-sets and sub-manifolds matching[END_REF], 3d stacks of density images [START_REF] Mirza | Variational and computational methods for flows of diffeomorphisms in image matching and growth in computational anatomy[END_REF][START_REF] Mirza | Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms[END_REF], curves [START_REF] Joan | Modeling planar shape variation via Hamiltonian flows of curves[END_REF][START_REF] Joan | Large deformation diffeomorphic metric curve mapping[END_REF] or surfaces [START_REF] Joan | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Vaillant | Surface matching via currents[END_REF]. In [START_REF] Vaillant | Surface matching via currents[END_REF] a theoretical framework is provided for encoding surfaces embedded in R 3 as currents ; the optimal control problem is defined consequently as matching of currents. The resulting mapping is used to study features of organs, detect anomalies, etc.

Conformal geometry theory has been used for translating a 3d shape matching problem into a 2d one by conformal parametrization of the boundary, via harmonic maps, conformal maps, least square conformal maps [START_REF] Xianfeng | Matching 3d shapes using 2d conformal representations[END_REF][START_REF] Ming | Optimized conformal surface registration with shape-based landmark matching[END_REF][START_REF] Sen | 3D surface matching and recognition using conformal geometry[END_REF]. In a series of articles devoted to shape matching and retrieval, 3d shape matching is achieved by extraction and comparison of topological features, namely Reeb graphs [START_REF] Masaki | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF], or others topological indicators. However, a Reeb graph is computed from the analysis of the isosurfaces of a real-valued function defined in a manifold ; for many shapes the extraction of such a function is not a trivial task.

More recently, in the field of Computer Graphics, the optimal transport point of view has been used to displace an input tetrahedral mesh onto a given object [START_REF] Bruno | A numerical algorithm for L2 semi-discrete optimal transport in 3D[END_REF][START_REF] Zhengyu | Optimal mass transport for shape matching and comparison[END_REF] and to register 2d images [START_REF] Steven | Optimal mass transport for registration and warping[END_REF]. In the optimal transport point of view a measure (derived on a density function supported by a computational mesh or defined as a sum of Dirac masses) is displaced into another by minimizing the celebrated Monge-Kantorovitch optimal mass transportation cost. Methods from optimal transport have the merit to be insensitive to a rigid-body motion of the target shape, which is of tremendous importance for searching shape similarity, and to allow topology changes.

Several shape matching methods use landmarks or anchor points for which the correspondence is supposed known and the deformation is guided by this information. Depending on the application, these feature points are defined manually (for example in the case of anatomical landmarks) or by automatic extraction, for instance as critical points of a well-suited scalar function defined on the shape. In [START_REF] Serge | Shape matching and object recognition using shape contexts[END_REF] for example, the authors start by distributing sample points on the contour of both shapes, that will be matched according to their 'shape context', namely a shape descriptor defined as an histogram of the relative position of the sampled points with respect to a given point on the shape. Our method does not make use of any preprocessing feature extraction. In the proposed problem the 'template' shape W 0 is numerically described by means of a computational mesh whereas the 'target' shape W T is known only via its boundary. The sole restrictive assumption is that both shapes share the same topology.

The paper is organized as follows. We first introduce the shape optimization framework in 6.2. Both methods, boundary and volumic, are presented in 6.3. Eventually, in 6.4, we compare the results obtained with the two functionals thanks to several numerical examples.

A shape optimization framework

Generally speaking, shape optimization is about optimizing a cost functional F(W) -or objective functional -over a set of admissible shapes O ad . A generic shape optimization problem reads as : min

W2O ad F(W), or max W2O ad F(W).
When the optimization problem is cast over a set O ad of shapes, some additional difficulties come into play, from both theoretical and numerical points of view. First, one needs a way for describing admissible variations over the space O ad and a notion of differentiation over the space O ad , which is needed for inferring optimality conditions. Moreover, in most problems the dependency of the objective functional on the domain is highly non-trivial ; the functional F may depend on geometric properties of the domain (i.e. curvature, normals, thickness) but also on a state u W -achieved as the solution of a PDE posed on W -and its derivatives. In this case we have to deal with the derivation of the state u W under consideration with respect to the domain.

Shape sensitivity analysis : Hadamard's shape derivative

Several notions of differentiation with respect to the domain are available in the literature. Here we refer to the Hadamard's boundary variation method (see e.g. [3,96,[START_REF] Francois | Quelques résultats sur le contrôle par un domaine géométrique[END_REF]). The general idea is to FIGURE 6.1 -Variation W q of a shape W according to Hadamard's method. establish a correspondence between variations of the domain and diffeomorphisms of R d (d = 2, 3). Thus, variations of a given shape W are considered under the form :

W q = (I + q )(W),
where q : R d 7 ! R d is a 'small' vector field. Indeed, each admissible variation W q of W is parametrized in terms of a transformation of the form I + q , which remains 'close' to the identity. The admissible vector field q is sought among the Banach space W 1,• (R d , R d ) of bounded and Lipschitz functions q : R d 7 ! R d endowed with the norm :

||q || W 1,• (R d ,R d ) := ||q || L • (R d ) d +||-q || L • (R d ) d⇥d , 8q 2 W 1,• (R d , R d ).
The convenience of the choice of this space follows from the following result (from [3]). Lemma 6.2.1 For any

q 2 W 1,• (R d , R d ) such that ||q || W 1,• (R d ,R d )  1, the application I + q is a Lipschitz homeomophism with Lipschitz inverse.
Using the correspondence between vector fields and domains, the Hadamard shape derivative is defined in terms of the classical derivative over a functional space. Definition 6.2.1 A function F(W) of the domain is said to be shape differentiable at W if the mapping q 7 ! F(W q ), from W 1,• (R d , R d ) into R, is Fréchet differentiable at q = 0. The associated Fréchet differential is denoted as q 7 ! F 0 (W)(q ) and called the shape derivative of F. The following expansion then holds :

F(W q ) = F(W) + F 0 (W)(q ) + o(q ), where |o(q )| ||q || W 1,• (R d ,R d ) q !0 ! 0.
The form q 7 ! F 0 (W)(q ) (at least at first order) does not depend on the tangential component of q on ∂ W. More precisely we have the following statement (see [3]). Proposition 6.2.2 Let W be bounded, Lipschitz domain in R d . Let F be a functional differentiable at W. If q • n = 0, n standing for the outer normal to ∂ W, then we have :

8q 2 W 1,• (R d , R d ), F 0 (W)(q ) = 0.
Moreover, for a large class of functionals F, the shape derivative F 0 takes the following form

8q 2 W 1,• (R d , R d ), F 0 (W)(q ) = Z ∂ W f q • n ds, (6.1)
where f is a scalar field defined on ∂ W, depending on the problem at stake. The choice of W 1,• (R d , R d ) as underlying Banach space guarantees that the application I + q is a Lipschitz homeomophism with Lipschitz inverse. However, in order to compute the shape derivative of a surface functional, we will sometimes require the vector field q to enjoy more regularity that the one inherited from the space

W 1,• (R d , R d ).
In such a case we will consider the space

C 1,• (R d , R d ) := W 1,• (R d , R d ) \ C 1 (R d , R d
) endowed with the norm :

8q 2 C 1,• (R d , R d ), ||q || C 1,• (R d ,R d ) := k  l=0 sup a2N d ,|a|=l || ∂ q ∂ x a || L • (R d ) d .
Note that an equivalent definition of shape derivative holds over the space C 1,• (R d , R d ).

General procedure for the minimization of a shape functional

From the expression of the shape derivative (6.1), we deduce that the vector field

f n, (6.2)
is a natural descent direction for F. Indeed, if q = t f n with t > 0 small enough, it guaranties that :

F(W q ) = F(W) t Z ∂ W f 2 ds + o(t 2 ) < F(W).
This paves the way for an iterative algorithm, producing a sequence (W k ) k=0,... of shapes, with decreasing values of F(W k ) : at each step, W k is updated according to

W k+1 = (I + q k )(W k ), where q k is (an extension to W k of) t f n. (6.3)
From the discrete point of view, the procedure generates a sequence of meshes (T k ) k=0,... . At each iteration the mesh T k is updated by pushing its vertices according to the descent direction :

8x vertex of T k , x 7 ! x + q k (x). (6.4)
The formal procedure described above boils down to deforming a shape W in the negative direction of the L 2 (∂ W) d gradient of the differential q 7 ! F 0 (W)(q ) :

F 0 (W)(q ) = h f , q • ni L 2 (∂ W) d . (6.5)
Unfortunately, the descent direction associated to the vector field (6.2) reveals unsuited for at least two reasons :

-the vector field featured in (6.2) is defined only on the boundary of W ; it has therefore to be extended to W as a whole so that it can be a guide for displacing the vertices of T ; if no particular attention is paid to this extension, the extended displacement field may impose an important stretching in W, making the motion of the vertices of T via (6.4) impossible to achieve without invalidating the mesh. To overcome both these difficulties we rely on a procedure of extension and regularization of the descent direction (6.2) by using the gradient of q 7 ! F 0 (W)(q ) associated to another inner product. This velocity extension -regularization issue is quite classical in shape optimization (see [START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] De | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF] and references therein), and can be thought of as an efficient preconditioning of the naive procedure (6.3). The general strategy is described in the following. Suppose that the space of admissible deformations is endowed with a structure of Hilbert space V . Let h•, •i V be the inner scalar product in V . Therefore, we can consider the following linear form l(•) acting on elements of V :

8v 2 V, l(v) = F 0 (W)(v).
Now let a V (•, •) be a suitable symmetric, V -elliptic bilinear form acting on elements of V . The Lax-Milgram theorem ensures the existence of a unique u W 2 V solution of the variational problem :

8v 2 V, a V (u W , v) = l(v).
The vector field u W is a descent direction for F since :

F 0 (W)( u W ) = a V (u W , u W ) < 0.
A suitable choice in our context is to use the extension inherited from a particular mechanical system, that of linear elasticity. See 6.3.2 for an example of such a procedure.

6.3 Two functionals for shape matching 6.3.1 Shape matching without constraint on the displacement Let W 0 and W T ⇢ R d be respectively 'template' and 'target' shapes. Under the assumption that they share the same topology, our purpose is to map W 0 onto W T by minimizing an energy functional which measure the discrepancy between the shapes. In order to achieve our goal, we introduce two shape functionals. At first, consider the following boundary functional P(W) :

P(W) = Z ∂ W d 2 (s, ∂ W T ) ds, (6.6)
where d(•, ∂ W T ) denotes the usual Euclidean distance function to ∂ W T . This functional is based on a classical paradigm -aligning two shapes by aligning their boundaries -and is used in a wide variety of situations in the context of non rigid registration (see [START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF] for example). Then, consider the following functional J(W) of the domain :

J(W) = Z W d W T dx, (6.7)
which involves the signed distance function d W T to W T , defined by :

8x 2 R d , d W T (x) = 8 < : d(x, ∂ W T ) if x 2 W T , 0 if x 2 ∂ W T , d(x, ∂ W T ) if x 2 c W T . (6.8)
Note that the functional J is defined on the whole domain of interest W. In order to decrease the value of J(W), the domain W must expand in the regions of the ambient space R d where d W T is negative (that is, in the regions comprised in W T ), and to retract in those where it is positive (see 6.2). By deforming the template shape as above we can always obtain better values for J until the deformed template shape reaches the target shape W T . Therefore, J(W) has a unique, global minimizer W = W T if we assumed that W T 2 O ad . A functional similar to (6.7) has also been used in [START_REF] Charles | GEOMETRICAL SHAPE OPTIMIZATION IN FLUID ME-CHANICS USING FREEFEM++[END_REF] for shape optimization purposes.

By a classical calculation (see [3,96]), the functionals (6.6) and (6.7) are shape differentiable at any (enough) regular shape W, and the associated shape derivatives reads respectively :

8q 2 C 1,• (R d , R d ), P 0 (W)(q ) = Z ∂ W ✓ ∂ d 2 ∂ n + kd 2 ◆ q • n ds, (6.9)
where k is the mean curvature of ∂ W, and

8q 2 W 1,• (R d , R d ), J 0 (W)(q ) = Z ∂ W d W T q • n ds.
(6.10) Following 6.2.1, a descent direction is inferred from formulas (6.9) and (6.10). An iterative algorithm devoted to the minimization of P (or J) will produce a sequence of shapes which are 'closer and closer' to W T . FIGURE 6.2 -Deformation of the template shape decreasing the value of J.

Shape matching with extension and regularization of the descent direction

As recalled in 6.2.1, it is crucial for numerical purposes to use the gradient of the mapping q 7 ! P 0 (W)(q ) associated to another inner product instead of the one featured in (6.5). Therefore, assume that all the considered shapes W are filled with a linear elastic material. Also, imagine that any such shape W contains a given subset w b W on which it is clamped. We now obtain a descent direction for P(W) as the unique solution u W belonging to H 1 w (W) d := {v 2 H 1 (W) d , v = 0 in w} of the linearized elasticity system :

8v 2 H 1 w (W) d , Z W s (u W ) : e(v) dx = P 0 (W)(v) = Z ∂ W ✓ ∂ d 2 ∂ n + kd 2 ◆ v • n ds, (6.11)
where e(v) = 1 2 (-v + -v T ) is the linearized strain tensor and s is the associated stress tensor via the Hooke's law with Lamé coefficients l , µ : s (u) = 2µe(u) + l tr(e(u))I. This vector field u W is a descent direction for P(W) since P 0 (W)(u W )  0 thanks to (6.11), and its advantages over the 'natural' deformation field

⇣ ∂ d 2 ∂ n + kd 2 ⌘
n are twofold : 1. u W is defined on the whole shape W ; owing to the regularizing effect of elliptic equations, it is intrinsically smoother (see for instance [42]) ; 2. Owing to the mechanical features of elastic displacements (notably their 'rigidity'), it is expected that u W will be more amenable to the displacement of the mesh T into a valid mesh via (6.4) ; see e.g. [START_REF] Baker | Mesh Movement and Metamorphosis[END_REF] for an example of use of elastic displacements in the context of mesh displacement. Furthermore, we can infer a descent direction for J(W) as the unique solution u W belonging to H 1 w (W) d of the variational problem :

8v 2 H 1 w (W) d , Z W s (u W ) : e(v) dx = J 0 (W)(v) = Z ∂ W d W T v • n ds. (6.12)
R From the numerical point of view, the choice of a subset w corresponds to a global alignment of shapes (cf. 6.3). This restriction we used to guarantee the well-posedness of Problems (6.11) and (6.12) could be replaced by adding a 0 th -order term. Note also that by construction the boundaries ∂ w and ∂ W do not intersect, preventing the loss of regularity at the interface between boundary conditions of different natures.

Numerical issues and examples

As far as the numerical setting is concerned, the template shape W 0 is discretized as a simplicial mesh (filled with triangles in 2d and tetrahedra in 3d). The computational mesh is non uniform : it is refined in the vicinity of the boundaries according to a curvature based sizing function and coarsened in the interior of the domain. This has proven to prevent severe distorsion/tangling of the elements avoiding the need to remesh the domain. In all the 2d examples, the chosen template mesh T 0 is a disk that has about 18 000 triangles. The set w chosen for aligning W 0 and W T is a small disk located in the interior of both shapes. Starting from the template shape W 0 we perform a standard gradient descent algorithm with adaptive step size in order to get a sequence of pairs (W k , T k ) of domains and their corresponding meshes with decreasing values of the functional. The general algorithm is the following : Algorithme 6.1 Start from an initial shape W 0 , for k = 0, . . . , convergence 1. compute the value of the shape functional ; 2. compute the state u W k of the considered elasticity system on W k by a P 1 finite elements discretization ; 3. perform a line search for getting the optimal step size t k so as to guarantee the decrease of the shape functional ; 4. advect the shape W k according to this displacement field, so as to get

W k+1 = (I + t k u W k )(W k ).

⌅

The algorithm stops when the step size is smaller than a fixed tolerance e. The global mapping from W 0 to W T is easily recovered by the composition of the different displacements obtained at each iteration. Note that by moving the mesh at each iteration we are able to devise a deformation between two shapes which are not close from one another by using a linear elastic model, which accounts in principle for 'small' displacements. Note that that we can compose the deformations from the template W 0 onto two targets W T ,W 0 T in order to obtain the mapping from W T to W 0 T and conversely.

Boundary functional P

In the boundary case, the computation of the shape gradient P 0 (W k ) requires at each iteration k the computation of : -the normal derivative of the distance function d at each vertex x i of the boundary ∂ W k by the formula :

∂ d 2 ∂ n = -(d 2 ) • n(x i ) = 2(x i p W T (x i )) • n(x i ),
where p W T (x i ) stands for the projection of x i on ∂ W T ; -the discrete mean curvature k of the boundary ∂ W k . In 2d we follow the approach described in [START_REF] Charles | GEOMETRICAL SHAPE OPTIMIZATION IN FLUID ME-CHANICS USING FREEFEM++[END_REF]. Let {x i } i=0,...,n be the list of points lying on ∂ W k . Let us suppose that the edges of the boundary have been sorted in such a way that for each vertex x i the vertex x i 1 (resp. x i+1 ) denotes the vertex placed before (resp. after) x i , according to a counterclockwise orientation. In such a situation the normal vector at point x i is estimated as the orthogonal of the (normalized) vector ! x i x i 1 . Then the discrete radius of curvature at vertex x i is given by the formula :

r i = 1 4 ✓ ! x i x i 1 • ! x i x i 1 n(x i ) • ! x i x i 1 + ! x i x i+1 • ! x i x i+1 n(x i ) • ! x i x i+1
◆ .

The discrete curvature k i at x i is taken as the inverse of the radius of curvature when the latter is not vanishing, and is set to zero otherwise. Let us now present several examples of 2d shape matching using the boundary functional. First, consider the test cases of 6.4. The algorithm converges in few iterations and is able to perfectly recover the target shape, even in regions of high curvature. Then, the same template is deformed onto the targets represented in 6.7. In these cases, the algorithm fails to converge, invalidating the mesh after few iterations. The occurrence of singularities is due to a lack of differentiability of the distance function when the set of projections of a point over the boundary of the target shape is not a singleton. In such case the shape derivative is not well defined since we have as many descent directions as projection points. This is an important drawback of the proposed functional.

Volume functional J

For the volume functional, the computation of the shape gradient J 0 (W k ) does not demand the evaluation of the discrete curvature or the normal derivative of the distance function. The only information required about the target shape W T is the datum of its signed distance function e.g. a P 1 piecewise affine function on the fixed mesh of a large computational domain D, which can be defined on a possibly non-conforming mesh (e.g. showing small gaps, overlapping entities, etc.). The calculation of the signed distance function to W T is performed using the algorithm [START_REF] Charles | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. An example of signed distance function to a target shape is shown in 6.6-(a). Let us present several 2d numerical examples to show the performances of the method. Both target and template meshes are embedded in a unit computational box with dimensions [0, 1] 2 . At first, we consider the test cases depicted in 6.5 that result in an excellent recovery of W T . In 6.8, note that the algorithm is also able to recover isolines of high curvature resulting in an excellent matching, even for template and target shapes very far from one another. Also notice that the use of linear elasticity equations makes the deformation smooth preserving the quality of the mesh under deformation. 6.6 shows the contour superposition between the deformed template mesh and the contour of the target shape model under consideration. Note that the template shape is never remeshed during the process. This allows to keep a point-to-point correspondence between the template mesh and the final mesh of the target, which may be useful in applications. However, in the example depicted in 6.8-(g)-(i), the matching in the extremity of the tail could be improved by adding new boundary points in this area.

Next, we consider a 3d example ; see 6.9. Both the target and the template meshes are embedded in a unit computational box D = [0, 1] 3 . The shapes W 0 and W T are aligned by choosing a small ball w in W 0 \ W T . The template mesh T 0 has about 9 000 triangles, and 1 300 iterations of the gradient descent algorithm have been performed to achieve convergence for a tolerance e = 10 6 , running in a few minutes on a standard laptop computer. The L 2 norm of the distance between the boundaries of the final shape W 1 300 and the target shape W T is 5.04 ⇥ 10 4 (much smaller than the minimal mesh size). Eventually, in 6.10 an ellipsoid template mesh is mapped onto a shape of hippocampus that was obtained from segmentation of medical images.

Comparison of the proposed methods

The meshes obtained with both functionals are valid, and as we never remesh the template during the process, it is easy to compare the results and the associated errors. The convergence analysis relative to all the considered 2d examples is detailed in Table 6.1. The quantitative analysis of the error is performed by computing the discrepancy between the boundaries of the last deformed template W k and the target shape W T .

For the 2d examples depicted on 6.4 and 6.5, we obtain similar results using the two proposed functionals. The main advantage of the volume functional is the ability to converge in all the proposed examples, with a little number of iterations. The boundary functional has shown good performances and even more accuracy in few examples (see 6.4 generating singularities. Note that the signed distance function to W T is not differentiable at points of the space having multiple projections over ∂ W T . However this lack of differentiability does not create ambiguities in the volume case because the expression of the shape derivative depends only on the value of the signed distance function and not on the projection points. Hence, the volume functional based on the paradigm 'ask to W 0 to fill all the space occupied by W T ' has proven to be more efficient in most situations with respect the one based on the alignment of boundaries.

In [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF] where the volumic method is used to efficiently infer correspondences in a database of 3d shapes coming from segmentation of medical images. The inferred correspondence is used to drive a new numerical method for facial reconstruction, which is a challenging problem in forensic science. Note also that the volume functional, although much lesser sensitive that its boundary counterpart, is not independent to the initial alignment of both shapes : a rigid body motion applied to the initial shape could lead to a very different mapping between template and target. Future works could compare combine the proposed method with techniques which are unsensitives to the initial alignment, notably the ones based on optimal transport theory. Indeed, the discrepancy between shapes could be measured by means of the Wasserstein distance (that, precisely, involved in optimal transport).

Introduction

Facial reconstruction (FR) is the process of deducing, from the sole knowledge of a dry skull, the shape of a deceased person's face [START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF]. It is still a challenge today with many different applications that can be medical, legal, or even historical.

The work on this paper began with a forensic science goal to identify bodies that could not otherwise be recognized. An essential part of this process is therefore the identification of the deceased's face. This provides a means of police investigation in some cases [146,[START_REF] Valentine | Forensic facial identification : a practical guide to best practice[END_REF] and valuable evidence in otherwise unsolvable cases [108,[START_REF] Iscan | Forensic Analysis of the Skull : Craniofacial Analysis, Reconstruction, and Identification[END_REF][START_REF] Taylor | Skull protection and preparation for reconstruction[END_REF]. Often the good results obtained draw the attention of the media and the public to the reconstructed image of a face [START_REF] Valentine | Forensic facial identification : Theory and practice of identification from eyewitnesses, composites and CCTV[END_REF][START_REF] Christopher | Craniofacial identification[END_REF]. However, the recognition of other elements such as personal clothing or accessories (glasses, hats, necklaces) displayed with facial recognition triggers recognition rather than the estimated facial morphology itself [START_REF] Mikhail | The face finder[END_REF][START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF][START_REF] Carl | Facial approximation-from facial reconstruction synonym to face prediction paradigm[END_REF]. Subsequently, the usefulness of these methods was verified in practice and often used in high-profile cases [START_REF] William | Use of facial approximation techniques in identification of Green River serial murder victims[END_REF]108,[START_REF] Iscan | Forensic Analysis of the Skull : Craniofacial Analysis, Reconstruction, and Identification[END_REF][START_REF] Carl | Facial approximation-from facial reconstruction synonym to face prediction paradigm[END_REF][START_REF] Taylor | Skull protection and preparation for reconstruction[END_REF][START_REF] Valentine | Forensic facial identification : a practical guide to best practice[END_REF]. However, it is beneficial to maximize their abilities to generate recognition from facial morphology. It is useful to use science to improve or test methods in the laboratory and in the field [START_REF] Carl | An overview of the latest developments in facial imaging[END_REF].

The FR problem is thorny because it is really ill-posed. Indeed, the solution is not unique : a skull does not correspond to a unique face [START_REF] Dennis | Probabilistic joint face-skull modelling for facial reconstruction[END_REF]. Some parameters such as ageing or corpulence are respectively difficult or impossible to extrapolate from the knowledge of a dry skull alone [42,[START_REF] Dennis | Probabilistic joint face-skull modelling for facial reconstruction[END_REF][START_REF] Sachio | The computer-aided facial reconstruction system[END_REF]. These difficulties make the FR a research subject that is still active today.

However, many methods have already been developed to try to solve the FR problem. The first use of facial approximation in a forensic case was reported in 1916 [START_REF] Author | Police find house used for murders[END_REF][START_REF] Author | Reconstruct body to solve murders[END_REF], but the method had been studied in the context of anthropological and anatomical research since 1898 [START_REF] Hv | Die Leistungsfahigkeit physiognomischer Rekonstruktionsversuche auf Grundlage des Schädels[END_REF][START_REF] Büchly | Die Persistenz der Rassen und die reconstruction der Physiognomie prahistorischer Schädel[END_REF]. Traditional methods are most often carried out manually in 2D [START_REF] Evenhouse | Computer-aided forensic facial reconstruction[END_REF][START_REF] Stephan | Building Faces from Dry Skulls : Are They Recognized Above Chance Rates ?[END_REF][START_REF] Karen | Two-Dimensional Facial Reconstruction from the Skull[END_REF][START_REF] Taylor | The wisdom of bones : facial approximation on the skull. Sous la direction de Marks MK CLEMENT JG[END_REF][START_REF] Douglas | Computer-Assisted Facial Reproduction[END_REF] or 3D [START_REF] Pat | Facial sculpture on the skull for identification[END_REF][START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF][START_REF] Stephan | Building Faces from Dry Skulls : Are They Recognized Above Chance Rates ?[END_REF][START_REF]Forensic Facial Reconstruction[END_REF] (portraits or sculptures). These methods essentially consist of three common steps [START_REF]Computerized forensic facial reconstruction[END_REF] :

1. equip the raw skull with a series of sparse anatomical landmarks, 2. apply an average thickness of soft tissue to each skull mark to estimate a corresponding mark on the face, 3. draw (2D) or sculpt (3D) a face corresponding to the estimated landmarks. Most of the time these traditional methods are performed by an artist who is an expert in facial anatomy. Although the results are often very plausible, the subjectivity of practitioners is inevitable [START_REF] Stephan | Anthropological facial 'reconstruction' -recognizing the fallacies, 'unembracing' the errors, and realizing method limits[END_REF]. These methods are also time-consuming and can take from a few days to a week to complete a reconstruction. These are two criteria, subjectivity and long implementation time, that we are trying to get rid of.

In addition, the blooming up of numerical methods over the past few decades has greatly affected many disciplines. Numerical simulations for FR have therefore also been developed since early 1990s [START_REF] Evenhouse | Computer-aided forensic facial reconstruction[END_REF][START_REF] Douglas | Computer-Assisted Facial Reproduction[END_REF][START_REF] Vanezis | Application of 3-D computer graphics for facial reconstruction and comparison with sculpting techniques[END_REF][START_REF] Vanezis | Facial reconstruction using 3-D computer graphics[END_REF], and through automatisation processes, have made it possible to overcome subjectivity and really accelerate the whole process. At this stage, numerical reconstructions are most often carried out with a statistical model and a database. The database can be either a 3D shaped database that associates a face with a skull, or a soft tissue thickness database with a given landmark, or a combination of the two. Here two types of approaches can be described : 1-computerized facial modelling using retroactive devices and 3D scanning of the skull [129,[START_REF]Computerized forensic facial reconstruction[END_REF][START_REF] Christopher | Craniofacial identification[END_REF][START_REF] Wilkinson | A blind accuracy assessment of computer-modeled forensic facial reconstruction using computed tomography data from live subjects[END_REF], and 2-computerized facial reconstruction using more complex and computerized 3D routines [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF][START_REF] Peter | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths : methodology and validation[END_REF][START_REF] Peter | Computerized craniofacial reconstruction : conceptual framework and review[END_REF][START_REF] Pierre | Anthropological Facial Approximation in Three Dimensions (AFA3D) : Computer-Assisted Estimation of the Facial Morphology Using Geometric Morphometrics[END_REF][START_REF] Parks | Preliminary performance assessment of computer automated facial approximations using computed tomography scans of living individuals[END_REF][START_REF] Gérald | A fully three-dimensional method for facial reconstruction based on deformable models[END_REF][START_REF] Subsol | Automatic 3D facial reconstruction by feature-based registration of a reference head. Sous la direction d[END_REF][START_REF] Tu | Towards a statistical basis for facial deformation modes in reconstruction[END_REF][START_REF] Tu | Face reconstructions using flesh deformation modes[END_REF][START_REF] Turner | A novel method of automated skull registration for forensic facial approximation[END_REF][START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF][START_REF] Vandermeulen | Automated facial reconstruction[END_REF].

Several recent examples of these methods, which have been made and tested in the last 20 years, can be given [START_REF] Carl | Facial approximation-from facial reconstruction synonym to face prediction paradigm[END_REF].

In all of these approaches, there is near-universal reliance on measurements of the thickness of the soft tissue of the face [START_REF] Peter | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths : methodology and validation[END_REF][START_REF] John | Making Faces Using Forensic and Archaeological Evidence[END_REF]129,[START_REF] Stephan | Beyond the Sphere of the English Facial Approximation Literature : Ramifications of German Papers on Western Method Concepts*[END_REF][START_REF] Carl | Facial approximation-from facial reconstruction synonym to face prediction paradigm[END_REF][START_REF] Herbert | On Gerasimov's plastic facial reconstruction technique : new insights to facilitate repeatability[END_REF], which renders the Russian-American classification entirely obsolete [START_REF] Stephan | Beyond the Sphere of the English Facial Approximation Literature : Ramifications of German Papers on Western Method Concepts*[END_REF][START_REF] Carl | Facial approximation-from facial reconstruction synonym to face prediction paradigm[END_REF][START_REF] Herbert | On Gerasimov's plastic facial reconstruction technique : new insights to facilitate repeatability[END_REF].

Starting with the simple methodology of placing soft tissue layers on the dry skull before deforming a facial mesh towards these estimated facial landmarks [START_REF] Peter | Computerized craniofacial reconstruction : conceptual framework and review[END_REF]. An other possibility, explored in [START_REF] Thomas | A method for automatic forensic facial reconstruction based on dense statistics of soft tissue thickness[END_REF], consists in using, in addition to these measurements, a face statistical shape model, allowing to obtain a global estimation, taking more into account the possible variations. However [START_REF] Carl | Facial soft tissue depths in craniofacial identification (part I) : an analytical review of the published adult data[END_REF] shows that, whatever the methods used, soft tissue thickness measurements are linked to a significant error (of about 10%).

Using this idea of local approximation initially, via statistical shape models, some studies uses a least square regression to obtain a smooth face [START_REF] Qingqiong | A regional method for craniofacial reconstruction based on coordinate adjustments and a new fusion strategy[END_REF][START_REF] Yongli | Craniofacial reconstruction based on a hierarchical dense deformable model[END_REF]. This idea of local approximation also includes an analysis of regional facial features, such as the position of eyeballs in orbits [START_REF] Craig | Commentary on : Stephan CN. Facial approximation : globe projection guideline falsified by exophthalmometry literature[END_REF]65,[START_REF] Pierre | Anatomical Placement of the Human Eyeball in the Orbit-Validation Using CT Scans of Living Adults and Prediction for Facial Approximation*[END_REF][START_REF] Kim | Three-dimensional prediction of the human eyeball and canthi for craniofacial reconstruction using cone-beam computed tomography[END_REF][START_REF] Parks | Automated facial recognition and candidate list rank change of computer generated facial approximations generated with multiple eye orb positions[END_REF][START_REF] Stephan | Facial Approximation : Globe Projection Guideline Falsified by Exophthalmometry Literature[END_REF][START_REF] Stephan | The Placement of the Human Eyeball and Canthi in Craniofacial Identification[END_REF][START_REF] Stephan | Further Evidence on the Anatomical Placement of the Human Eyeball for Facial Approximation and Craniofacial Superimposition[END_REF][START_REF] Wilkinson | Measurement of Eyeball Protrusion and Its Application in Facial Reconstruction[END_REF], the shape of the nose [4,[START_REF] Pat | Facial sculpture on the skull for identification[END_REF][START_REF] Mala | Pronasale position : an appraisal of two recently proposed methods for predicting nasal projection in facial reconstruction[END_REF][START_REF] Af Ridel | Skeletal dimensions as predictors for the shape of the nose in a South African sample : a cone-beam computed tomography (CBCT) study[END_REF][START_REF] Christopher Rynn | Appraisal of traditional and recently proposed relationships between the hard and soft dimensions of the nose in profile[END_REF][START_REF] Stephan | Predicting Mouth Width from Inter-Canine Width-A 75% Rule[END_REF][START_REF] Strapasson | Forensic facial reconstruction : relationship between the alar cartilage and piriform aperture[END_REF][START_REF] Herbert | On Gerasimov's plastic facial reconstruction technique : new insights to facilitate repeatability[END_REF][START_REF] Utsuno | Establishment of a prediction method for the mid-facial region of unknown human Mongoloid skeletal remains[END_REF], the width of the mouth [START_REF] Murphy | Mouth width prediction in craniofacial identification : cadaver tests of four recent methods, including two techniques for edentulous skulls[END_REF][START_REF] Wilkinson | The relationship between the soft tissues and the skeletal detail of the mouth[END_REF], the shape and size of the lips [START_REF] Mala | Vertical lip position and thickness in facial reconstruction : a validation of commonly used methods for predicting the position and size of lips[END_REF][START_REF] Wilkinson | The relationship between the soft tissues and the skeletal detail of the mouth[END_REF], the morphology of the ears [1,[START_REF] Guyomarc'h Et | The validity of ear prediction guidelines used in facial approximation[END_REF], the morphology and thickness of the mastication muscles [START_REF] Carl | The human masseter muscle and its biological correlates : a review of published data pertinent to face prediction[END_REF].

We can also mention the methods used to identify the parameters that characterize the mapping between a skull and its face via a database correlating several pairs of skull and associated face. First [START_REF] Wuyang | A PCA-Based method for determining craniofacial relationship and sexual dimorphism of facial shapes[END_REF], which uses a statistical form model applied to a large CT-scan database. Others, such as [66,[START_REF] Pascal | A 3D face model for pose and illumination invariant face recognition[END_REF][START_REF] Sarah | A spatially-dense regression study of facial form and tissue depth : towards an interactive tool for craniofacial reconstruction[END_REF], suggest a last step that take into account corpulence.

Finally, a last type of approach consists in statistically simulating the database [START_REF] Dennis | Probabilistic joint face-skull modelling for facial reconstruction[END_REF], also allowing to present not a unique face but a set of possibilities. A large majority of these methods are based on statistical or probabilistic models.

The method presented in this article takes advantage of these recent contributions, as well as the history of the FR, but attempts to distinguish itself by a deterministic approach to the problem. Our method consists of adding the facial muscles one by one, then the fat and skin, as the artists did but using mathematical modelling. This requires building parametric models of facial muscles to be able to adapt them to the underlying facial skeleton. For the moment, we have focused on the masseter muscle, which is the largest muscle of the face and influences its shape in a prominent way [START_REF] Yu | Botulinum toxin A for lower facial contouring : a prospective study[END_REF]. Then we applied the elastic deformation method of [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF]. We also propose here, an application to validate the process from the point of view of recognition, which also allows us to compare our methods with others.

Anatomy

The starting points of our approach are anatomical and medical. The main idea is to integrate this knowledge into the mathematical model. Indeed, the article by Stephan et al. [START_REF] Carl | The human masseter muscle and its biological correlates : a review of published data pertinent to face prediction[END_REF] shows us that although the anatomy of the muscle of the face, is well defined in the literature, few recently published FR methods take this information into account.

In this first part we are therefore interested in the anatomy of the masseter muscle, and the possible repercussions of its shape on that of the face.

Link between the skull and the face

The anatomy of the skull largely determines the shape of the face. Different studies show, for example, the link between long face and short face and mandibular angle [START_REF] Opdebeeck | The short face syndrome[END_REF][START_REF] Stephen | The long face syndrome : vertical maxillary excess[END_REF]. Others have a link between the angular measurements of the skull and masticatory muscles such as the masseter [START_REF] Ariji | Three-dimensional morphology of the masseter muscle in patients with mandibular prognathism[END_REF][START_REF] Faouzi | Relationship between dimensions of muscles of mastication (masseter and lateral pterygoid) and skeletal dimensions : study of 40 cases[END_REF][START_REF] Boom | A comparison of human jaw muscle cross-sectional area and volume in long-and short-face subjects, using MRI[END_REF][START_REF] Girard | Significance of the Frankfort-mandibular plane angle to prosthodontics[END_REF]. This is also what we contributed to in the first chapter of this thesis, by including new features.

The important link between bone morphology and muscles morphology described in the literature was the first step in the development of this FR method. The originality is to do this anatomical reconstruction numerically and not manually, as [START_REF] Mikhail | The face finder[END_REF][START_REF] Gerasimov | The reconstruction of the face from the basic structure of the skull[END_REF].

The masseter muscle

The human head has nearly a hundred muscles. A large majority of them are so-called mimic muscles. However, the most voluminous muscles are those of the masticatory system, which is composed of four main muscles : the masseter, the temporal, the internal pterygoid and the lateral pterygoid.

The masseter muscle is one of the most important muscles of the face, as musch by it is its size, its location as its function. Indeed, it has the largest volume of muscle of the face, and therefore it participates to the actual shape of the face [START_REF] Yu | Botulinum toxin A for lower facial contouring : a prospective study[END_REF]. The masseter is an elevator muscle of the mandible, thus its role in mastication is also very important : it is one of the most powerful muscles in the mandible. Another pertinent muscle candidate, superficial and quite large, would have been the temporal muscle, but because of its location, it is often hidden by the hair, and in this area of the face it does not seem to exceed half of the soft tissue thickness [START_REF] Stephan | The superficial temporal fat pad and its ramifications for temporalis muscle construction in facial approximation[END_REF]. This is why we chose to study first the influence of the masseter muscle on the morphology of the face. Keeping in mind to add other muscles to our facial reconstructions afterwards, such as the temporal muscle.

Many anatomical studies describe masseter muscles such as [START_REF] Bonnefoy | Anatomie descriptive et fonctionnelle de l'articulation temporo-mandibulaire[END_REF]. Stratified dissections and anatomical sections in different spatial planes showed that the muscle masseter had a typically pennate structure consisting of a succession of alternating musculoaponeurotic layers. The muscle has three clearly differentiated parts : the superficial, intermediate and deep beams. The masseter insertions are on the gognion and zygomatic bone [START_REF] Busato | Management and control of isotonic contraction generated stress : evaluation of masseter muscle deformation pattern by means of ecography[END_REF].

Facial reconstruction process

We call here « reconstruction », the process of giving a face to a skull, for forensic investigation [START_REF] Iscan | Forensic Analysis of the Skull : Craniofacial Analysis, Reconstruction, and Identification[END_REF]. We have to keep in mind that it is an approximation process, which will gave us a guess of the face associated to the unknown skull : one will not be able to recover a large amount of face features (eyes, hair, lips and ears). The main parameters of interest for us is therefore the shape of the face.

The first step of our methods is to reconstruct from mesurements on the skull the muscles like the masseter. And in a second step we add a mask of the remaining soft tissus with a general reconstruction process, based on elastic deformation process [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF].

From the 3D meshes obtained after segmentation we therefore apply the method developed in chapter 2 of this thesis, and the following steps are implemented :

1. elastic deformation of an ellipsoid of reference to each muscle, in order to obtain a mesh having the same characteristics (number of points, triangles, etc.) whatever the masseter muscle considered (this elastic deformation is obtained using the volumis morphing algorithm, described in the Chapter 6), 2. Principal Component Analysis of all the displacements of the reference ellipsoid mesh to all the masseter muscles of our database (in the way described in the Chapter 2) to get the basic shapes, 3. reconstruction of a masseter muscle from only some measurements made on the skull with a linear combinaison fo the basic shapes. First of all, we can present the case of Figure 7.4, where we use Reduced Basis (RB) to reconstruct a muscle whose true shape is known but which was not used to compute the principal components. The three muscles on the left are represented using respectively 5, 10 and 20 coefficients, while the muscle on the right is our target. With 5 coefficients, we can see that the overall shape of the muscle is already captured. However, it takes between 10 and 20 coefficients to obtain a more accurate result. In consideration of our purpose in this chapter (and as opposed to those required in Part II of this thesis), we have subsequently decided to consider a masseter by its approximation to the first 5 coefficients, the main information being contained within.

The next step is to link the first five RB coefficients of the masseters muscles with the shape of their associated skull. Here we propose and compare two approaches to represent the skull :

A -anthropological measurements on the skull, a method that is similar to what can be found in existing methods, or B -to use a parametric model to describe the mandible (bone underlying the masseter), allowing a global automatisation of the process.

A -Measurements on the skull

We first made various measurements on the skull, nine in total, between anatomical landmarks. In Figure 7.5 we can visualize these different landmarks :

1. Glabella : most prominent point in the median sagittal plane between the supraorbital ridges, 2. Orbital : most inferior point on margin of orbit, 3. Gnathion : point located perpendicular on mandibular symphysis midway between pogonion and menton, 4. Point F : intersection between the dental plane and most anterior part of the mandible branch, In a second step, in order to avoid having to perform these various tedious measures, we took a direct interest in the shape of the mandible. We used the same process, described in Chapter 2 of this thesis, for mandibles as for masseters, in order to describe them using only a few coefficients, extracted from the RB. To best match the shapes, it was necessary to « cut » the mandibles in half, each half referring to the masseter muscles on the corresponding side. The warping step, described as optional for masseters, proved to be necessary and very practical here. This allowed to consider a mandible as a closed 3D surface, see Figure 7.7.

After the implementation of this whole process, we also obtained a parametric model for the mandibles, each described by the first ten RB coefficients of all the mandibles in the database. We now have two ways to describe the shape of the skull, and more specifically that of the mandible : 9 anthropological measurements or the ten first RB coefficients of the mandibles. The first method consisting in linking measurements on the skull and the shape of the masseter muscle has been largely legitimized by the literature but also by the first chapter of this thesis. However which is not the case with the second approach, where we want to link the coefficients of the mandibular RB with those of the masseter. Table 7.2 summarizes the excellent correlations, 0.8986 on average (p-value < 0.05, with a Pearson correlation test), found between the first five coefficients of the two RBs. This allows us to validate our approach conceptually, before moving on to the next step. With the two sets of information characterizing the shape of the skull and more particularly the mandible, it is then possible, using a machine learning algorithm, such as neural networks, to deduce the first 5 RB coefficients of the masseters. This allows us to find the shape of a muscle on a dry skull to identify. In Figure 7.8, we present the reconstructions of masseter muscles after learning, on two examples. The two masseters on the left are represented after deduction of the learning algorithm of the first 5 RB coefficients, the first according to the 9 cephalometric measurements described above, the second, on the left, according to the mandibular RB coefficients. The middle masseter, in blue, is our target. After examining some examples, it seems quite clear that learning through cephalometric measurements is more effective. Indeed, in Figure X, the color maps show a lower error for reconstruction via measurements (left) compared to that via mandibular RB (middle). However, this remains only a first observation because when we made this first analysis we had only 68 masseter couples and measurements at our disposal, and 30 masseter couples and mandibles segmented separately (and thus allowing us to make a RB). So very little data for effective learning. It is in response to this problem of lack of data that we have begun work on semi-automatic segmentation presented in Part II.

In the context of facial reconstruction and in view of the correct results found by learning through measurements made on the skull, we first used this reconstruction of the masseters. We keep in mind that the second method (mandibular RB), in the light of the significant correlation found, may be more efficient, or at least as effective, with more data in the long term. This would allow full automation of the process.

Results presentation : combinaison with a general method of facial reconstruction

Our method has mostly two steps. The first one is to add to the dry skull, the muscles we reconstruct from knowledge of the skull, measurements or 3D shape (Figure 7.9). This is what we presented in the section above and in more details described in the second chapter of this thesis. This allows local anatomical knowledge to be integrated into the complete facial reconstruction process, this is what is new in comparison to [START_REF] Chiara | Mesh deformation strategies in shape optimization. Application to forensic facial reconstruction[END_REF] and which is original because it has not yet been developed, numerically and automatically, to our knowledge. The second step is to deform a mask of the other soft tissues (Figure 7.10) with the methods of [START_REF] Buhan | Volumic and Boundary shape matching : comparison of two optimization methods[END_REF][START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], but instead to consider as input the dry skull, we consider the skull augmented with the masseter muscles reconstruction. To be able to carry out this second step, it is necessary to go through a « warping » step described in chapter 2, which allows us to describe our skull + masseter group as a 3D surface mesh. We can visualize on Figure 7.11, the same skull and its masseter as that of Figure 7.9, but after the « warping » process. FIGURE 7.11 -The same example skull as in Figure 7.9, but after using the warping algorithm to close it Once all the skulls in the database have been closed, with the masseter muscle, segmented or reconstructed, it is possible to map each skull/face pair of the database, via a deformed skull/face reference pair, using an elastic morphing described in the previous chapter, to each of these so-called pairs. Thus, when identifying a new skull, it is sufficient to match it with the reference skull/face pair, to be able to deform any face from the database to this new skull.

In view of the limited data, all the reconstructions were carried out via a "Leave-One-Out" process. This allowed us to perform 24 face reconstructions in a simulated « real situation », since during the reconstruction of an individual A, the skull/face pair of this individual was removed from the database, before performing the reconstruction.

However, there is one last point to be clarified : when an unknown skull must be identified, which skull from the database is chosen to reconstruct its face, via elastic deformation of the skull/associated face mask ? In our case we decided to select the 3 skulls from the database, the closest to the unknown skull (according to their L 2 scalar product) , then to average the three faces obtained. On this point, different methods could be proposed and compared in order to obtain better results.

We can present one exemple of facial reconstruction performed by this new method in the Figure 7.12 : FIGURE 7.12 -One example of reconstruction with our method on the left versus the 3D segmented face on the right With this only Figure 7.12, it is really difficult to appraise the accurancy of the method, with the naked eye. In Figure 7.13, we can more easily compare the discrepancy of thoses recontructions with the real face by adding a color map, which express the error made between the surface reconstruction and the 3D segmented surface of the face. The Figure 7.13 presents, on the left, the reconstruction of an individu of our dataset made with the methode of [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF] without the masseter muscle (front and side) ; on the middle, the reconstruction performed with the segmented masseter muscle ; and on the right, the reconstruction performed with the reconstructed masseter muscle (using the Chapter 2 method). After comparison a similar result is obtain, exepted in the region of the anatomical details ad-up of the masseter, where a better accurancy can be shown for both of the reconstruction including a masseter muscle. Moreover, these two reconstructions performed with the masseter muscle differ little from each other. This leads us to believe that our muscle reconstruction is precise enough for this application. The Figure 7.13 is only an example of reconstruction. We used these three methods on a set of 24 skulls, to reconstruct their associated faces, with a leave-one-out method. The descriptive statistics of the errors found for these different applications are expressed in the Table 7. These results seem completely disproportionate compared to the visualization of Figure 7.13, where the error is less than 5mm on the majority of the face. Several remarks expressed below may be the cause of this.

R

-The ears are not included in the face-matching process, since these structures are not linked with the underlying skull morphology. So the errors made on them has to not be take into account. -The neck is taken into account in the shape-matching process but it is not a good indicator of the mistake made. Indeed, depending on the acquisition modalities of each scanner, the images do not all stop at the same level, leading to a large variability that is not representative of the face.

To illustrate this, we present below, in Figure 7.14, a second visualization with a wider color scale between 0, 2 and 10mm against 0, 1 and 5mm. This allows us to confirm that the majority of the error on the face is between 0 and 2mm. To test these hypotheses in more detail, and to be able to make a more interesting and accurate comparison between the three methods, one idea would be to select an area of interest and calculate the associated errors.

This procedure is often used in facial reconstruction and would allow us to set aside areas we know we are not able to reconstruct with the methods used (ears and neck), and compare these results with the existing published literature. It could also allow us to study in detail a particular anatomical area such as the area where the masseter is located. These first results, which are still encouraging, show the interest of using anatomical knowledge for the facial reconstruction process, even with numerical tools. But this kind of validation, often carried out in the field of facial reconstruction, is purely theoretical. Thus, we also develop in the following section a more realistic method of facial reconstruction recognition.

Validation

In order to show that thoses methods of reconstruction can be useful in forensic science, we have built a validation test, which compares the ability of a large population to identify, from a picture, a 3D volume of the face of the same individual and their reconstructions with and without the muscles. On the Figure 7.15, we show an example of reconstruction, from our database, performed in three different ways : without any muscles, with the real 3D segmented masseter and with the reconstruted masseter. This is the recognition of those four different faces, from a picture of the same individual, that we want to evaluat. To be able to determine if our methods can improve the recognition of someone's face in comparison of other methods, we also need to show that people are able to identify a reconstruction face as well as with the segmented 3D representation, and later on better with the masseter muscle than without. After research in the literature, there is to our knowledge only few practical validation method tested for the recognition of facial reconstructions in forensic sciences. Most of the time on a small reconstruction sample (about ten),with few test participants (less than 100) and in 2D. The results are mixed and tend to show little interest in facial reconstruction for forensic identification [START_REF] Stephan | Assessing facial approximation accuracy : How do resemblance ratings of disparate faces compare to recognition tests ?[END_REF][START_REF] Stephan | Building Faces from Dry Skulls : Are They Recognized Above Chance Rates ?[END_REF].

We have thus concived and realized a survey. The survey is constructed as follow : a front and side view of an individual is presented to the user, who must choose one among five 3D face representations, one of which is the one he recognizes as the individual of the pictures. Every participants who answers to the suvey have 10 questions to validate.

For each of these 10 questions, the five 3D faces are selected from one of these four categories (hidden from the participant) :

-the real 3D segmented face -the reconstruction 3D face without the masseter muscle -the reconstruction 3D face with the 3D segmented masseter muscle -the reconstruction 3D face with the masseter muscle reconstructed

The selection of the categories, chosen for a question, is made randomly.

For now, 122 persons use this test application, and they all performed the test on 10 cases, randomly chosen from the 24 individuals of our database. We are aware of the too few participants who answered the survey for now, but those results are a really interesseting first attempt to perfome a validation of the recognition process and a good way to compare different methods [START_REF] Stephan | The reproducibility of facial approximation accuracy results generated from photo-spread tests[END_REF].

On the Figure 7.16 and 7.17, the application and how to use it on a example case are shown (at the end of the chapter).

Results

In this section we first present the Table 7.4 of the normalized descriptive statistics of our survey. The signification of each colums is detailed above :

-RealSkin = the real 3D segmented face -SkullOnly = the reconstruction 3D face without the masseter muscle -RealMass = the reconstruction 3D face with the 3D segmented masseter muscle -RBMass = the reconstruction 3D face with the masseter muscle reconstructed In average a 3D real face is recognized in 79.6% of cases, a 3D reconstruction 45.4% of cases, while chance would be one among five (20%).

With a basic t-test we can show that the real face is with any doubt best regonized with respect to the reconstruction from pictures, but not perfectly either. And we found no significantly differences between the three reconstructions preposed, even if the RBMass is slightly larger. This can be attributed to different interpretations :

1. the accurancy of the method and the recognition process are two differents things, and people do not seem to look at the masseter region to recognize someone, 2. the only few people we tested is not enough to demonstrat a real difference.

Discussion and Conclusion

This article show a complementary method to [START_REF] De | A facial reconstruction method based on new mesh deformation techniques[END_REF], with more accurancy in terms of distance to the real face, with slight improvement in terms of recognition. A way to even improve our results would be to -add the reconstruction of other muscles than the only masseter, -add a parametric way to control the fat tissus, in order to better represent the corpulence, -increase the number of individuals in the database, particularly to diversify morphologies as much as possible, -look for a more subtle way to define the skulls closest to the unknown skulls, or to combine the faces obtained.

It also propose a way to compare the results of differents facial reconstruction studies, with a protocol validation easy to implement. This has not been done extensively in the field of facial reconstruction before. To complete this first approach, we could consider completing this test with two other tests : -submit 5 pictures of individuals in order to identify a face reconstructed in 3D, this would allow a more realistic situation to be simulated, compare for a selected picture, which of the two reconstructions, with or without muscle masseter, is identified as the best. In addition to that an interesting study would also be to compare our validation method with the result obtained by automatic face recognition by computer, a currently prolific field of research, and also with the more classical statistic facial reconstruction [START_REF] Parks | Biometric correspondence between reface computerized facial approximations and CT-derived ground truth skin surface models objectively examined using an automated facial recognition system[END_REF][START_REF] Parks | Recognizability of computer-generated facial approximations in an automated facial recognition context for potential use in unidentified persons data repositories : Optimally and operationally modeled conditions[END_REF]. -enlarging the database is the first thing we need to do to improve our results, which are produced on a small database. Two possibilities are then possible : to acquire new "real" data manually segmented (at least partially), which would allow a better consideration of the different existing morphologies ; but also a diversification of the shapes already present in the database, generated with the help of perturbation of these shapes by spherical harmonics ; -the whole part dedicated to the data alignment between them is central in the 3D objects management procedures, as discussed in this thesis. It is therefore quite possible to imagine different methods improving this aspect in the work presented in chapters 3 and 4, as the impact on the results could be significant. One proposal would be to integrate rigid movements to the variations added by the spherical harmonics used to enlarge the database, as suggested above ; -in order to be as systematic as possible a more detailed comparison of the 3D results presented in Chapter 4 with the pre-existing literature results, either from a numerical error point of view or from a 3D visualization point of view.

Finally as far as morphing is concerned, as presented in the Chapter 6, other aspects can still be explored : -the use of a faster method than the gradient descent used, in order to accelerate the process ; -the combination of the two functionalities presented : the volume one, fast and efficient in many situations, and the surface one, less robust but more accurate when the target shape has singularities. This would allow to better take advantage of the advantages of each of the two methods ; -improve the so-called surface method, which does not succeed systematically, by including a remeshing step in some iterations, thus avoiding the algorithm to get stuck due to a lack of points.

Although there are still several possible improvements for each of these methods, as detailed above for example, they have all already been implemented on different applications and real cases, detailed below. This allows feedback on the potential and effectiveness of each, as well as the possibility of critical thinking to further improve them.

1. the semi-automatic segmentation, which uses the 3D shape database of the masseter muscle to model the object to be segmented. This method under development can still be improved, but it already speeds up the first part of the process. 2. facial reconstruction with anatomical precision on the masseter muscle, a tool we used to reconstruct the face of the New Caledonian warlord Atai, presented in the introduction (Figure 7.18). In this method it would be possible to add other muscles, to integrate a nose reconstruction method to obtain more details, or to allow reconstruction on altered skulls (without teeth or with a damaged mandible, cf. Appendice B).

As regards to facial reconstruction, one aspect that was left out during this thesis is the reconstruction of secondary characters. In the Figure 7.18, a possible interpretation, designed by Sebastien Le Gall (ISCD, Sorbonne Université, Paris), of Atai's face is presented, based on known historical information about him. However, it would be interesting to try to systematize this approach using work on the genetic determination of secondary traits (skin, eyes and hair colours, hair type, etc.) for exemple.

The clinical context has also not been put aside, since on the one hand, thanks to an analysis between cephalometry and the masseter muscle size, we were able to establish and confirm the anatomical hypotheses on which all the tools and applications developed are based (Chapter 1). On These validation methods are still under development and can be further improved, especially with regard to face recognition (Chapter 7). Indeed, they are for the moment dependent only on visual recognition by subjects participating in the study, whereas many facial recognition software already exist. It would therefore be interesting to validate these recognition applications and compare the results.

A.1 Introduction

One of the recurring and time consuming challenges in medical imaging is the segmentation of parts of interest. Segmentation is the process of dividing a medical image into distinct anatomical regions [68]. This has been for years the effort of specialists and it is quite time consuming. Beyond the specifics takes to each particular application, the general issue is to find a way to identify automatically an anatomical region of interest in order to speed up the process. When the contrast is high enough (e.g. for bones, or between bones and tissues) automatic segmentation is quite efficient now and various methods are proposed and implemented [START_REF] Yan | A new accurate and precise 3-D segmentation method for skeletal structures in volumetric CT data[END_REF][START_REF] Shafry Mohd | 3D bones segmentation based on CT images visualization[END_REF][START_REF] Jing | Fast segmentation of bone in CT images using 3D adaptive thresholding[END_REF]. Likewise for anatomical major organs like the brain or the heart [START_REF] Siqi | Multi-scale structured CNN with label consistency for brain MR image segmentation[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF][START_REF] Thomas | Handbook of medical imaging[END_REF]. On the contrary, when the contrast is not high enough (e.g. segmentation between different soft tissues on CT scans) the most common imaging modalities fail and there are very few approaches available.

Understandably enough, most of the current segmentation methods are based on the gray scale which can let experts easily differentiate 2 objects [START_REF] Dzung | Current methods in medical image segmentation[END_REF]. Several classifications of those methods exits, we can cite [START_REF] Daniel | Medical image segmentation : Methods and software[END_REF] for an overview. We limit ourselves to mentioning here a few approaches, for some more detail surveys we can look at [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF][START_REF] Geert | A survey on deep learning in medical image analysis[END_REF], for example.

Thresholding and edge detection are two very common approach used for image segmentation [START_REF] Antonio | Evaluation of methods for ridge and valley detection[END_REF][START_REF] Ramesh | Thresholding based on histogram approximation[END_REF]. They are simple and effective approaches to segment anatomical details featured by high contrast, as, for instance, the bones (that appear as light objects on dark background). However it is more challenging for the detection of objects that do not have a marked gray scale signature, as the muscles. We can also cite the region based method. Three other types of process can be exposed here : the deformable models [START_REF] Tobias | Statistical shape models for 3D medical image segmentation : a review[END_REF], the atlas guided [START_REF] Eugenio | Multi-atlas segmentation of biomedical images : a survey[END_REF] and all the approaches dealing with artificial intelligence, which has recently blooming up in the field of medical image segmentation [START_REF] Neeraj | Automated medical image segmentation techniques[END_REF][START_REF] Daniel | Medical image segmentation : Methods and software[END_REF].

We propose here to develop a model based semi-automated method to help segmentation, which have most in common with the three last types of methods cited above. In particular, the shape of the object of interest is represented by an elastic deformation field of a reference configuration. A reduced basis representation of the object of interest is introduced, that exploits an a posteriori knowledge of a database of available segmentations. Given an image to be segmented, the elastic deformation field can be approximated as a linear combination of a small number of modes. This provides, in turn, a powerful regularisation and a significant speed-up of the procedure needs since few apropriate points need to be selected by the user in order to define the field. We can find a similar idea in [START_REF] Florian | Shape-aware surface reconstruction from sparse 3D point-clouds[END_REF], but with a probabilistic point of view, and where all the modes of the PCA are used as a regularisation and not as a dimentional reduction.

The proposed method is complementary of the ones exposed before, first because it requires some preliminary segmentations to "build" the reduced basis, and second because it can be used on indifferently on CT and MR image for any type of object, even the difficult one to identify like soft tissues on a CT.

The paper is structured as follows : in Section A.2 the problem and the notation are presented. Then, in Section A.3 some numerical experiments are detailed.

A.2 Notation and problem setting.

Most medical imaging produces a gray-scale image, i.e. a function defined over a subset of R2 (or R 3 for 3D imaging), f : R 2 ! [0, 1] associating a scalar value to each point of the space (denoted later on by x = (x 1 , x 2 ) (or by by x = (x 1 , x 2 , x 3 ) in 3D)). After a manual segmentation of a region of interest is performed, it is possible to model this region by using an elastic morphing that maps a reference object : e.g. an ellipse in R 2 (or an ellipsoid in R 3 ) onto the boundary of the two-dimensional (or three-dimensional) segmented region. In the present work, for the sake of simplicity, we will restrict to two-dimensional configurations. Having a collection of segmented objects (belonging to a coherent family) is equivalent to having a set of 2D elastic displacement fields mapping a same (reference) ellipse into the different segmented objects.

Let us assume that the reference configuration is an ellipse E ⇢ R 2 ; a point x 2 ∂ E is mapped through a displacement field into a point x 2 ∂ F ⇢ R 2 , where F is the boundary of the segmented object. It holds :

( Id + u : E ! F x 7 ! x = x + u(x ), (A.1) and u 2 H1 (E).

From an acquired database consisting of N 2 N ⇤ such elastic displacements : D = u (1) , . . . , u N 1 let us assume that a set of n 2 N ⇤ Proper Orthogonal Decomposition (POD) modes is extracted :

j i : E ! F, i = 1, . . . , n. (A.2)
When a new image to be segmented is available, the unknown elastic displacement mapping the reference configuration onto the image can be approximated by using the POD basis as follows :

u ⇡ u n = n  i=1 a i j i . (A.
3)

The goal is to find the coefficients a i in such a way that the reference ellipse, mapped by the displacement field Id + u n is as close as possible to the actual image, hence providing an already good guess for the segmentation. In order to identify the coefficients a i a semi-automated procedure is devised. The image does not correspond, in general, to the reference configuration (this would be the case in which a i = 0, 8i). A number of m 2 N ⇤ points x ( j) 2 ∂ F is selected in which there is a discrepancy between the actual image and the reference configuration. Each of these points correspond to an equation for the displacement field : there exists a x ( j) ⇤ 2 ∂ E such that :

x ( j) = x ( j) ⇤ + n  i=1 a i j i (x ( j) ⇤ ), j = 1, . . . , m. (A.4) This is a non-linear system of equations in the variables a i , x ( j) ⇤ ; however, if the points x ( j) ⇤ are known, this system is a linear system of equations for the coefficients a i . In the present work, the points x ( j) ⇤ are chosen as the ones minimising the distance to the point x ( j) 2 :

x ( j) ⇤ = arg min

x 2E kx ( j) x k 2 `2 , j = 1, . . . , m. (A.5)

This corresponds to a regularisation for the solution of the non-linear system Eq.(A.4). The linear system to be solved reads :

Fa = d, (A.6) F jik = j i,k (x ( j)
⇤ ), (A.7)

d j,k = x ( j) k x ( j) ⇤k , (A.8)
where F is a tensor, assumed to be of maximal rank (remember that it is constructed from POD vectorial modes). Of course, as soon as m n, the system is solved in a least-square sense, and admits a unique solution. When discretising and solving this system, the following choices are made. The minimisation step described by Eq.(A.5) is solved by using a nearest neighbor algorithm that, among the points of the discretised ∂ E selects the one which is closest to a given x ( j) . Then, the vector system in Eq.(A.6) is unfolded and then solved by using a rank revealing QR decomposition.

A.3 Numerical Experiments

In this section, some numerical experiments are shown on a simple 2D test case to assess the performances of the proposed semi-automated segmentation approach. The reference configuration, depicted in Figure 1.a), is a unit circle, divided into two semicircles. The images are defined by deforming this configuration by a parameter dependant field, generated as follows : the displacement of each point on the reference circle is normal to the circle and constructed by summing the set of N h 2 N ⇤ harmonics (eigenfunctions of the Laplace-Beltrami operator), scaled with a parameter dependant coefficient : let the undeformed configuration points be denoted by x 2 R 2 , r = kx k 2  1, J 2 [0, 2p] (i.e. x 1 = r cos(J ), x 2 = r sin(J ), the unit vector in the radial direction be denoted by êr and the displacement providing the deformed image configuration be x 2 R 2 :

x = X(x ) = x + r where a j = a j (µ), have an amplitude which decreases as a power law of the harmonics numbers : ( j + 1) b with b = 1.5, and µ is a parameter that identifies each image. Consider the following radial displacement, depending upon the parameter µ 2 [0, 1] 2 . The harmonics coefficients are chosen according to the following expression : Remark that the displacement is purely radial and that the center of the reference circle r = 0 does not move. This corresponds to the case where the various anatomic regions that have been segmented are "aligned" and properly scaled. Furthermore, given this way of building a database of elastic displacements, the POD modes reduce to a set of vector displacements in radial direction whose intensity is given by the value of the harmonics. The discretised shape is built by discretising uniformly the interval J using N J = 512 points. The image (the deformation of the reference shape) is cut by N p 2 N ⇤ horizontal planes, uniformly distributed on the image, to mimic a medical exam like CT or RMI which are sliced. Each plane intersects the outer boundary of the object in two points (m = 2N p ), each of which corresponding to one of the halves of the deformed circle. These points provide an equation for the displacement field to be approximated by using the POD modes.

Several numerical experiments are made to assess the performances of the proposed reconstruction method, by varying the number of modes and the number of cutting planes used. In particular, the number of modes used varied from n = 6 to n = 20 and, for each choice the number of planes was varied from the minimal number of planes needed to have a square system n  m to m = 2n. In Figure A.1.b) an example of reconstruction is proposed, when the image was produced by using N h = 100 harmonics, n = 20, m = 80. The result is quite precise, all the details of the image are well represented, whether it is the outer edge where points are selected or the edge of the inner hemicycles where no points are selected.

The reconstruction error are shown in Figure A.2.b for a population of 100 target displacements (that do not belong to the database used to construct the POD modes). The errors are measured in `2 norm with respect to the discretised images. The average error and the standard deviation are shown when the minimal number of planes and the maximum number of planes are used. In general, the errors are compatible with the ones considered as compatible with realistic applications. When the minimum number of planes is used a stability issue can be seen when a large number of modes is used and the standard deviation is quite large. On the contrary, when N p = 2n (orange), the lack of stability is delayed, the errors decrease and the standard deviation is much smaller, reaching a minimum for n = 17.

A.4 Conclusions and Perspectives.

In the present work a semi-automated reconstruction method is proposed to help the segmentation of anatomical details in the cases in which the contrast is not enough to use fully automated techniques. The method is based on the use of a reduced basis to approximate an elastic displacement mapping a reference configuration into the actual shapes of the anatomical details of interest. The method uses a nearest neighbors optimisation procedure coupled to a least square approach to solve the problem. A 2D test case is proposed to mimic realistic applications and assess the performances of the method. A set of synthetic generated images was used and the error of the numerical experiments showed that the proposed approach can lead to a viable tool. In this respect, further investigations are mandatory both from the theoretical and applied point of view. In particular a systematic error analysis will be carried out and more realistic tests will be performed, in 3D configuration and by taking the noise and the other experimental error sources into account.

Especially in 3D, several challenges are added to those previously described, especially in the case of a real application. The first one is the constitution of the manual segmentation database : time-consuming and requiring the involvement of an anatomist expert at first. However, the idea is that we will soon be able to propose a first phantom, which will require only a few corrections and will save a lot of time while continuing to feed the database. Moreover, the mapping of the points would remain assured by the same idea as in 2D from a deformed reference configuration. However, instead of using a nearest neighbor algorithm, we would opt to define a probability density of each point in the point cloud of the corresponding dataset. In a second step, it is envisaged to replace the POD by a glutton algorithm, to increase robustness.

B.1.1 Github

All the code develloped during my PhD and needed for the process are availlable on the github of the ISCD Tools Box. It is necessary to install the whole « Commons » library and « facilepython ». To understand in detail the installation of these libraries, it is possible to refer to the « Read Me » that corresponds to them on github. The one concerning « facile-python » is also detailed below.

B.1.2 Python Librairies

The scripts are written in python 3, of which you must have an installed version. Anaconda is the most convenient to install and avoids having to re-download libraries, but classic python 3 can be perfectly sufficient and takes up less disk space.

The additional python libraries to install are : -scipy -numpy -scikit-image They can be installed with : pip3 install scikit-image numpy scipy ou conda3 install scikit-image numpy scipy

B.1.3 Executables

You must also have access to the following executables (their path must be specified correctly in the lib-exe.py file, so that they can be invoked as they are on the command line) :

-tetgen -mmg3D -mmgs -blender -warping -morphing -medit All of those are in open access on the web, either on the ISCD GitHub or on the website of the official developper.

B.1.4 Installation testing

In the scripts folder, from « facile-python », are present all the necessary python files, which can be used as independent executables (see the Section A.5).

« tests.py » is the file that allows you to run the code on test data. This allows you to test that everything is installed and that the different sripts used during the reconstruction are working.

Il se lance ainsi : python3 scripts/tests.py -d tests/ [-m] Note that the -m option will open each result with medit, useful to check that the signed distances are correctly calculated, that the references for the masks are the right ones, or that the meshes are not degenerated. lib-exe.execute("mmgs-O3 %s -o %s -hausd %f -hgrad %d" % (maillage.mesh, sortie.mesh, 0.01, 1.23))

B.4.3 Run several functions or scripts in parallel

When it is necessary to handle many cases, and therefore try to benefit from parallelization, it is best to encapsulate the actions to be performed in parallel in a function taking a single argument (this is what is done in the « database.py » file), and to use the « parallel » function of the lib-exe file : lib-exe.parallel(myFunction, myArguments) This function will loop through "myArguments" (which must therefore be a list), and will launch the "myFunction" function on each of the elements of "myArguments", in parallel.

Note that the function used may contain calls to programs executed on command lines.

B.4.4 Handle and edit a step

The « databse.py » script performs a set of function calls on files contained in folders.

To add a new step (or modify an existing one) :

-we create a function that takes as argument the name of the file to convert, -in this function, we run the corresponding script with the arguments created from the file name, to which we have added the path of the relevant folders, -we create a list of all the files that we want to process in the corresponding folder, -only files whose names are not already in the file output folder after processing are kept in this list, -to test or perform only a few steps, it is best, for example, to use the test() function of the database.py file by following the example presented.

the conclusion below. We can appreciate the ease with which these tests were set up : fully automatic, with the possibility of running them in parallel, it took less than 12 hours to obtain these results.

In addition, the stability of the simulations provided by the software is also very interesting : the overall shape of the face remains the same in all the simulations, showing that the change of a bone element alone does not affect the entire face.

Moreover, local changes in the bone shape (teeth, nose) also only locally affect the shape of the face, confirming the overall stability but also the fact that the software really takes into account the underlying elements of the bone skeleton. This last very positive point in the case where the skull is complete allows us to attest to a certain precision of the method. However, in cases where the reconstruction of some bones is necessary for the process because they are poorly preserved (as is often the case in anthropology), a small modification can lead to a significant local modification of the shape of the face. It is therefore important to keep this aspect of the problem in mind. To solve this problem, numerical methods of automatic reconstruction of the shape, and not manual, are entirely possible, and remains an ongoing research approach.

It is important to note that the entire project in which the facial reconstruction software presented here is available on forensic medicine. Indeed, the works presented in this appendix are to be carried out from an anthropological point of view. The two disciplines having of course common points, but not the same objectives. Forensic medicine focuses on the identification of a deceased person, while anthropology makes hypotheses about the faces of the entire human lineage.

  numéro 1) et le muscle temporal (Figure 7 numéro 2).

FIGURE 7 -

 7 FIGURE 7 -Situation anatomique (1) du muscle masséter, (2) du muscle temporal

FIGURE 8 -

 8 FIGURE 8 -Reconstruction d'un muscle masséter via ACP avec 5 coefficients, 10, et 20, face au vrai muscle

FIGURE 9 -

 9 FIGURE 9 -Image de scanner (a) fenêtre d'aquisition "tissus mous", et (b) "osseuse"

FIGURE 10 -

 10 FIGURE 10 -Segmentation (a) du crâne, (b) du visage, et (c) des muscles masséters

FIGURE 11 -

 11 FIGURE 11 -Application de validation de la reconstruction faciale

FIGURE 12 -

 12 FIGURE 12 -Identification de points à la bordure du masséter avec l'application semi-automatique de segmentation

FIGURE 17 -

 17 FIGURE 17 -Les deux fenêtres d'acquisition des images scanner (a) tissus mous, et (b) os

FIGURE 18 -

 18 FIGURE 18 -Maillage d'un muscle masséter après l'algorithme du marching cubes à gauche et après remaillage à droite

FIGURE 20 -

 20 FIGURE 20 -Description schématique du warping : (a) Calcul du champ de déplacement (flèche noire) sur le template W k (en gris) entourant le modèle source (en rouge), (b) Identification des points d'intersection entre le champ de déplacement et le modèle source (points en blanc), et (c) Actualisation de la forme selon le déplacement calculé. Schéma extrait de [158].

FIGURE 22 -µ 3 (µ 3 (

 2233 FIGURE 22 -Template non déformé à gauche et une déformation et sa reconstruction à droite

FIGURE 23 -

 23 FIGURE 23 -Décroissance des valeurs propres après ACP, avec (en bleu) et sans (en noir) paramétrisation

FIGURE 25 -

 25 FIGURE 25 -Représentation 3D d'un crâne et de son visage associé : (a) et la superposition des masséters segmentés manuellement et semi-automatiquement, (b) et du muscle masséter reconstruit après segmentation semi-automatique uniquement

FIGURE 27 -

 27 FIGURE 27 -Morphing 2D d'un disque en un lapin

FIGURE 28 -

 28 FIGURE 28 -Morphing 2D d'un disque en un lapin

FIGURE 32 -

 32 FIGURE 32 -Assemblage des deux surfaces crâne et visage d'un même individu (à gauche), maillage tétrahédrique de l'espace entre le crâne et le visage, appelé masque (à droite).

FIGURE 33 -

 33 FIGURE 33 -Exemple de crânes avec masséters reconstruits (à gauche) puis fermés par warping (à droite), en dessous section médiale de ces mêmes crânes

FIGURE 34 -

 34 FIGURE 34 -Reconstruction d'un visage à partir de son crâne (à gauche), à comparer avec le vrai visage 3D (à droite).

  number 1) and the temporal muscle (Figure 41 number 2).

FIGURE 41 -

 41 FIGURE 41 -Anatomical situation (1) of the masseter muscle, (2) of the temporal muscle

  f j , with a coefficient applying to the j-ième principale componentf , for a given individual k.

FIGURE 42 -

 42 FIGURE 42 -Reconstruction of a masseter muscle via PCA with 5 coefficients, 10, and 20, facing the real muscle

FIGURE 43 -

 43 FIGURE 43 -Scanner images (a) "soft tissue" acquisition window, and (b) "bone" acquisition window

  (a) and (b)

FIGURE 44 -

 44 FIGURE 44 -Segmentation (a) of the skull, (b) of the face, and (c) of the masseter muscles

FIGURE 46 -

 46 FIGURE 46 -Points identification at the edge of the masseter with semi-automatic segmentation application

FIGURE 51 -

 51 FIGURE 51 -The two scanner image acquisition windows (a) soft tissue, and (b) bone

FIGURE 53 -

 53 FIGURE 53 -Five masseter muscles identical with the exception of a rigid movement

FIGURE 54 -

 54 FIGURE 54 -Schematic description of the warping algorithm : (a) Computation of the displacement field (black arrow) on the template W k (in grey) surrounding the source model (in red), (b) Identification of the intersection points between the displacement field and the source model (white points), and (c) Update of the shape according to the calculated displacement. Schema from [158].

FIGURE 56 -

 56 FIGURE 56 -Template undeformed on the left and deformation and reconstruction on the right

FIGURE 59 -

 59 FIGURE 59 -3D representation of a skull and its associated face : (a) and the superposition of manually and semi-automatically segmented masseters, (b) and the reconstructed masseter muscle after semi-automatic segmentation only

FIGURE 60 -

 60 FIGURE 60 -The two segmentation software used : on the right, Amira, and on the left, the semi-automatic method

FIGURE 61 -

 61 FIGURE 61 -2D morphing of a disc into a rabbit

FIGURE 64 -

 64 FIGURE 64 -Manual 2D (top) and 3D (bottom) facial reconstruction method. Images from[START_REF]Facial reconstruction-anatomical art or artistic anatomy ?[END_REF]

FIGURE 66 -

 66 FIGURE 66 -Assembly of the two skull and face surfaces of the same individual (left), tetrahedral mesh of the space between the skull and the face, called a mask (right)

FIGURE 67 -

 67 FIGURE 67 -Example of skulls with reconstructed masseters (left) then closed by warping (right), below medial section of the same skulls

FIGURE 1 . 2 -

 12 FIGURE 1.2 -Measurement of the masseter dimensions on (a) frontal and (b) axial views.

TABLE 1 . 3 -

 13 Correlation coefficients between masseter and skeletal measurements (ns : Non-Significant)In Figure1.4, the scatter plots of masseter measurements that are significantly correlated to angular bones dimensions are highlighted in blue.

FIGURE 1 . 4 -

 14 FIGURE 1.4 -Data plots showing significant correlation, highlighted in blue. The continuous straight lines represent linear regressions, which explain the 20% (R 2 ) average variability of the data, while the discontinuous straight lines represent nonlinear regression.

FIGURE 2 . 2 -

 22 FIGURE 2.2 -Segmentation of the masseter muscle, in pink, on a head CT-Scan

FIGURE 2 . 3 -

 23 FIGURE 2.3 -Automatic skull segmentation in blue, with thresholding method, in the presence of dental artifacts

  , 2.3 or 2.4. From all these 2D images that are manualy segmented, we can obtain a 3D representation of an object of interest, here the skull, the face and the masseter muscle.

FIGURE 2 . 5 -

 25 FIGURE 2.5 -Some example of isosurface-cell intersection in the Marching Cubes algorithm.

FIGURE 2 . 7 -

 27 FIGURE 2.7 -Alignment of a 3D pavement after using this method

FIGURE 2 . 8 -

 28 FIGURE 2.8 -Alignment of 2 masseters with an ellipsoid template : on the left colomn the raw muscles after remeshing ; on the right after the alignement process on the elipsoid template (in the middle).

2

 2

FIGURE 2 . 12 -

 212 FIGURE 2.12 -One object in our reduced basis reconstructed with 5, 10 and 20 coefficients with the POD method on the top and with the Greedy one on the bottom. On the right the control

3

  Proof-of-Concept in 2D . . . . . . . . . . . . . 101

4

 4

FIGURE 3 . 1 -

 31 FIGURE 3.1 -On the left, we illustrate the definition of a scanner image as a 3D matrix whose unitary element is the voxel ; and on the right is shown the selection of one of these voxels in red, which corresponds to the segmentation process.

FIGURE 3 . 2 -

 32 FIGURE 3.2 -Segmentation examples : (a) an axial CT-Scan image of a head, (b) an automatique segmentation of the skull, and (c) a manual segmentation of the masseter muscle.

1 )

 1 b with b = 1.5, and µ is a parameter that identifies each image. The level of differente level of regulation can be illustrated by the Figure 3.4.

FIGURE 3 . 4 - 2 ⇥

 342 FIGURE 3.4 -Illustration of the regularisation coefficient (a) b = 1, (b) b = 1.5, (c) b = 2.Consider the following radial displacement, depending upon the parameter µ 2 [0, 1] 2 . The harmonics coefficients are chosen according to the following expression :g j = A j 2⇥ sin (2p j p µ 1 ) + cos (2p j p µ 2 ) ⇤ ,

FIGURE 3 . 5 -

 35 FIGURE 3.5 -Examples of automatically generated shapes in our 2D database

FIGURE 3 . 7 -

 37 FIGURE 3.7 -Examples of automatically generated shapes in our 2D database

FIGURE 3 . 8 -

 38 FIGURE 3.8 -Visualisation of two results (a) of a slightly distorted image and (b) of much more distorted image

FIGURE 3 . 9 -

 39 FIGURE 3.9 -Test case presented in Section 3.3 : Error mean and standard deviation

( 1 )

 1 FIGURE 3.10 -Visual representations of the first few real spherical harmonics

µ 3 (µ 3 (

 33 l+1) b cos(2pl(m + 1) p µ 1 ) cos(2pl(m + 1) p µ 2 ) l+1) b sin(2pl(m + 1) p µ 1 ) sin(2pl(m + 1) p µ 2 ) (3.13)

TABLE 3 . 2 -

 32 Mean reconstruction error 3.5 Conclusions and Perspectives.

  Figure 4.3.

FIGURE 4 . 3 -

 43 FIGURE 4.3 -An example of a bounding box selection on a frontal (left) and an axial slice (right)

FIGURE 4 . 4 -

 44 FIGURE 4.4 -(a) Points cloud selected in the image coordinates system, (b) the masseter mesh associated in the unit box

FIGURE 4 . 5 -

 45 FIGURE 4.5 -Schema of the head of our 3D image data

( 4 . 5 )

 45 with (x t , y t , z t ) the points coordinates in the translated bounding box.We then build the change of orientation matrix C :

FIGURE 4 . 6 -

 46 FIGURE 4.6 -Cloud points superimposed on the masseter mesh associated in the unit box

5 . 2 .

 52 The first results obtained, based on a dataset of 45 masseters to build the modes, are visualised on the Figure 4.7. In this figure we see the segmented (left) and reconstructed (right) 3D masseter, placed on the corresponding skull face couple. As well as their superimposition, for better comparison, in the middle of the figure.

FIGURE 4 . 7 -

 47 FIGURE 4.7 -Reconstruction in 3D of a masseter projected in its initial space : on the left the manualy segmented one, on the right our reconstruction, on the middle a comparison of them both.

FIGURE 5 . 1 -

 51 FIGURE 5.1 -Segmentation of the masseter from the 3 scanner cross-sectional planes : A : The final 3D form, B : The axial plane, C : The coronal plane, D : The sagittal plane. 1 : The separation of the two traces shows a segmentation of two different beams of the masseter. 2 : the segmentation of the sagittal plane deduced from the other two planes.

FIGURE 5 . 4 -FIGURE 5 . 5 -

 5455 FIGURE 5.4 -The posterolateral contour of the masseter is traced first on the axial (1) and then coronal section (2), followed by a check on the sagittal section(3)

0 TFIGURE 6 . 3 -

 063 FIGURE 6.3 -Target and template shapes sharing a common fixed subset w.

FIGURE 6 . 4 -

 64 FIGURE 6.4 -Simple 2d examples for the boundary functional. Boundary of the target shape (blue contour) and computational mesh of the template shape W k at iteration k.

FIGURE 6 . 5 -

 65 FIGURE 6.5 -Simple 2d examples for the volume functional. Boundary of the target shape (blue contour) and computational mesh of the template shape W k at iteration k.

FIGURE 6 . 6 -

 66 FIGURE 6.6 -(a) Isolines of the signed distance function. (b)-(c)-(d) Discrepancy between W T and the deformed template shape for the test cases of 6.8.

FIGURE 6 . 7 -FIGURE 6 . 8 -

 6768 FIGURE 6.7 -Three 2d examples for the boundary functional. Boundary of the target shape (blue contour) and computational mesh of the template shape W k at iteration k. The use of the surface functional (6.6) generates singularities.

FIGURE 6 . 9 -FIGURE 6 . 10 -

 69610 FIGURE 6.9 -An example in 3d with the volume functional : the mapping on the Igea model. (a) Template shape W 0 . (b) Deformed shape W k for k = 20. (c) Deformed shape W k for k = 1 500. (d) Isovalues of the signed distance function to the target shape W T . (e) Discrepancy between W T and W 1 500 .

FIGURE 7 . 1 -

 71 FIGURE 7.1 -Anatomical location of the masseter muscle and its three layers : (a) a schema, (b) a picturs of a dissection we performed (1-the superficial layer, 2-the vertical medium layer 3-the deep layer of the masseter)

FIGURE 7 . 4 -

 74 FIGURE 7.4 -RB reconstruction of a masseter muscle : from (a) 5 coefficients, (b) 10 coefficients, (c) 20 coefficients and (d) the real muscle

5 .

 5 Point H : intersection between the dental plane and most posterior part of the mandible branch, 6. Gonion : most posterior inferior point on angle of mandible, 7. Coronion : tip of the coronoid process of the mandible, 8. Condylion Superior of the Temporo-mandibular articulation : most superior point on the condyle of mandible, 9. Porion : most superior point of outline of external auditory meatus, 10. Point E : intersection between the dental plane and the middle of most exterior part of the mandible branch, 11. Point G : intersection between the dental plane and the middle of most interior part of the mandible branch.

FIGURE 7 . 5 -

 75 FIGURE 7.5 -Landmarks use for the measurments made on the skull

FIGURE 7 . 7 -

 77 FIGURE 7.7 -(a and b) Mandible mesh model issued from segmentation of medical images, (d and c) Closed surface model of a half mandibule generated with our warping algorithm

FIGURE 7 . 8 -

 78 FIGURE 7.8 -Two examples of reconstructed masseter via machine learning, on the left with the measurements, on the right with the RB of the mandibule, and on the middle the target shape.

FIGURE 7 . 9 -

 79 FIGURE 7.9 -Example of a skull of our database with (a) no muscle, (b) the real segmented masseter muscle, and (c) our masseter muscle reconstruction

FIGURE 7 . 10 -

 710 FIGURE 7.10 -An example of a mask which include the knowledge of the masseter muscle

FIGURE 7 . 13 -

 713 FIGURE 7.13 -Color map of the errors made on an example of reconstruction without (on the left), with the segmented masseter muscle (on the middle), and with the reconstructed masseter muscle (on the right).

  3. The signification of each colums is detailed above :-SkullOnly = the reconstruction 3D face without the masseter muscle -RealMass = the reconstruction 3D face with the 3D segmented masseter muscle -RBMass = the reconstruction 3D face with the masseter muscle reconstructed

FIGURE 7 . 14 -

 714 FIGURE 7.14 -Color map of the errors made on an example of reconstruction, different color scale, without (on the left), with the segmented masseter muscle (on the middle), and with the reconstructed masseter muscle (on the right).

FIGURE 7 . 15 -

 715 FIGURE 7.15 -One example reconstruction with the three tested cases with (a) only the skull, (b) the real segmented masseter, (c) the reconstructed masseter, and (d) the 3D segmented face

FIGURE 7 . 16 - 1 FIGURE 7 . 17 -

 7161717 FIGURE 7.16 -Application of validation -1

FIGURE 7 .

 7 FIGURE 7.18 -(Top) 3D Skull then facial reconstruction of Atai from the front and profile, (Bottom) on the left a bronze mortuary casting of the face of Atai, and on the middle and the right a suggested face with secondary features, for the reconstruction of Atai, by Sebastien Le Gall

FIGURE A. 1 -

 1 FIGURE A.1 -Test case presented in Section A.3 : (a) Reference configuration (b) an image and its reconstruction.

a j = A j 2 ⇥ sin (2p j p µ 1 )

 21 + cos (2p j p µ 2 ) ⇤ , In the present work K = 100 harmonics was used to generate the displacement fields. A set of N = 100 uniformly distributed in [0, 1] 2 samples of the parameter µ are taken, and, for each, the radial displacement is constructed. The POD basis was computed. The singular value decay is shown in Figure A.2.a, exhibiting a faster decay if compared to the one typical of the spherical harmonics for a given regularity.

FIGURE A. 2 -

 2 FIGURE A.2 -Test case presented in Section A.3 : (a) Decreasing of singular values with parametrisation (b) Error mean and standard deviation.

FIGURE C. 6 -

 6 FIGURE C.6 -Comparison of the position of the mandible between simulations 3-a and -b

TABLE 2 .

 2 1 -Discrepancies between reconstructions and target, using PCA algorithme

			Leave-All-in	Leave-One-Out
		5 coef	10 coef 20 coef 5 coef	10 coef 20 coef
	Mean 0.95495 0.50475 0.30059 1.01208 0.55879 0.37108
	SD	0.27560 0.10551 0.04472 0.30112 0.13366 0.08206
	Min	0.47886 0.31587 0.20424 0.49172 0.34276 0.22627
	Max	2.04961 0.86739 0.45761 2.10330 1.04464 0.69908
	Greedy Algorithm

TABLE 2 .

 2

	2 -Discrepancies between reconstructions and target, using Greedy algorithme
			Leave-All-in	Leave-One-Out
		5 coef	10 coef 20 coef 5 coef	10 coef 20 coef
	Mean 1.23136 0.66722 0.39767 1.30032 0.75220 0.49460
	SD	0.42973 0.24113 0.18634 0.38255 0.19295 0.10824
	Min	0.00000 0.00000 0.00000 0.52633 0.33344 0.24367
	Max	2.27286 1.15129 0.65106 2.47522 1.44458 0.99699

TABLE 5 .

 5 .2.

		Operator 1	Operator 2
	Mean	2.8 mm	2.9 mm
	Max	4.4 mm	3.9 mm
	min	0.8 mm	0.7 mm
	SD	1.0 mm	0.8 mm

2 -Descriptive statistic of the error performed with the semi automatic segmentation

TABLE 6 .

 6 for instance) but it is very sensitive to the initial alignment of both target and template shapes, failing to converge in most cases and 1 -Comparison between the boundary and volume functionals. Computation error and number of iterations needed to reach convergence for the 2d cases.

		Boundary functional		Volume functional
	Test case • • • 1. Flower 2. Lip 3. Rabbit 4. Cat 5. Monkey	L 2 mean 4.8 ⇥ 10 4 2.9 ⇥ 10 5 Do not converge Do not converge Do not converge	Iterations 400 500 1 000 680 300	L 2 mean 6.4 ⇥ 10 4 1.9 ⇥ 10 4 1.0 ⇥ 10 3 5.7 ⇥ 10 4 1.2 ⇥ 10 3	Iterations 500 2 000 500 1 500 2 300

Table 7 . 1

 71

	:

TABLE 7 .

 7

1 -Description of the measurements made on the skull

FIGURE 7.6 -Exemple of mesurements made in 3D B -A parametric model of the mandibule

TABLE 7 .

 7 2 -Correlation between the 5 coefficients of the mandibule and the masseter RB, with cMass the masseter coefficients and cMand the mandibule coefficients, both in descending order

		cMass5	cMass4	cMass3	cMass2	cMass1
	cMand5	0.9844	0.9923	0.8725	0.9075	0.9444
	cMand4	0.9802	0.9956	0.8855	0.9178	0.9342
	cMand3	0.9723	0.9930	0.8791	0.9068	0.9209
	cMand2	0.9749	0.9929	0.8763	0.9066	0.9266
	cMand1	0.9801	0.9894	0.8613	0.8974	0.9404
	Machine learning results

TABLE 7 .

 7

		Skull Only	RealMass	RBMass
	Mean	11.74 mm	11.85 mm	11.77 mm
	Sd	14.09 mm	14.09 mm	13.97 mm
	Min	72.66 mm	72.63 mm	72.49 mm
	Max	0.086 mm	0.086 mm	3.37 mm

3 -Description of reconstrution errors

  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . 199 List of Figures . . . . . . . . . . . . . . . . . . . . . . . 205 Bibliography . . . . . . . . . . . . . . . . . . . . . . . 207

		Introduction
	A.2	Notation and problem setting.
	A.3	Numerical Experiments
	A.4	Conclusions and Perspectives.
	B	Annexe Facial Reconstruction Software
		185
	B.1	Installation
	B.2	Creation of the database
	B.3	Facial reconstruction from a new skull
	B.4	How to add step ?
	B.5	Independent use of the scripts
	C	Annexe Anthropological Application 193
	C.1	Context
	C.2	The test case
	C.3	Numerical results
	C.4	Conclusion
		Books
		Articles
		Phd Thesis
		Misc

This database is assumed to be large enough to represent well the full coherent family we are interested in.

This database is assumed to be large enough to represent well the full coherent family we are interested in

Remark that this point might not be unique.

Acknowledgments

The authors would like to kindly thank P. Frey, Y. Maday, C. Dapogny and the partners of project FaciLe for fruitful discussions.

All the tested scripts, mentionned above, are in fact called in a « global » file allowing to create a reduced 3D database, named « database.py », which manages the different steps described in chapter 2 (data download, scaling, alignment, warping, morphing, etc).

This script recovers the initial files present in the « mesh » folder, then performs the required transformations one after the other, writing each result in a folder corresponding to the step performed.

The script is executed as follows : python3 scripts/database.py -d DATA-DIRECTORY -t TEMPLATES-DIRECTORY with :

-DATA-DIRECTORY : folder in which all the files will be processed, and where the different step results will be stored in separate folders created by this script. Note that the root folder « DATA-DIRECTORY » must be created beforehand. -TEMPLATES-DIRECTORY : folder containing the .mesh files used as templates for the different steps of the algorithm (sphere used for skull warping, reference elipsoid for masseter alignment and morphing, etc.).

B.3 Facial reconstruction from a new skull

Once the reduced database of 3D shapes generated by the script presented above, there is only one quick step left to rebuild a new face.

You have to run the script « reconstruct.py », this way : python3 scripts/reconstruct.py -i UNKNOWN-SKULL -d DATA-DIRECTORY -t TEMPLATES-DIRECTORY with : UNKNOWN-SKULL the mesh of the skull whose face we want to reconstruct.

B.4 How to add step ? B.4.1 Script python stand-alone

To create a new step, the idea is to create, as much as possible, a new python script that can be executed in "stand-alone" mode with a few arguments on the command line.

Some guidelines for a new script :

-manage arguments on the command line with « argparse », -check that the arguments are correct (that the files exist for example), -use mesh = lib-msh.Mesh("fichier.mesh") to open a new file .mesh edit « mesh.verts », « mesh.tris », « mesh.scalars », « mesh.tets » which correspond to the mesh elements, -use mesh.write("sortie.mesh") to write a file .mesh incorporate an external program (such as warping or morphing, code in C and previously compile).

B.4.2 Incorporate an external program

To incorporate an external program into a python script, the most convenient way is to use the « lib-exe.execute() » function of the « lib-exe.py » file. For example : lib-exe.execute("mmgs-O3 maillage.mesh -o output.mesh -hausd 0.01 -hgrad 1. 23") or, in a more practical form to use variables as program arguments :

carve.py

Creates a rough shell of a .mesh file with a « Space Carving » technique, which has the effect of voxelizing the object. May be more appropriate than a convex envelope to prepare the warping and obtain a manifold surface after remeshing for example.

python3 carve.py -i INPUT -o OUTPUT [-r RESOLUTION]

INPUT : .mesh file to carve OUTPUT : output .mesh file RESOLUTION : resolution of the voxel lattice to consider, default at 51. Set a smaller value if the points are scattered, you can set a larger value for a very dense mesh size.

convert.py

Converts a .obj or .stl file to .mesh using meshlabserver to avoid binary file problems. .mesh or .sol file to write, If the format is .mesh, the file and its associated .sol will be created, if .sol, only the solution.

download.py Downloads all .stl and .mesh files from a folder on an FTP server, and renames them according to their original name ("003-OsB.obj" will be downloaded as "003-bone.obj"). python3 download.py -a IP -u USER -p PWD -i FTPFOLDER -o LOCALFOLDER IP : IP address of the FTP server USER : username on the server PWD : associated password FTPFOLDER : path of the folder on the FTP server from which to download the files LOCALFOLDER : local folder where to download the files fill.py Fills a surface with a volumetric mesh and an icosphere in the center. This step is necessary to be able to perform the morphing.

python3 fill.py -i INPUT -o OUTPUT [-c cx cy cz] [-r RADIUS]

INPUT : .mesh file for the surface to be filled OUTPUT : .mesh output file cx, cy, cz : coordinates of the center of the icosphere to be included in the surface (default 0,0,0) RADIUS : radius of the icosphere (default 0.1)

icp.py

Computes the result of the transformation by ICP of a source mesh to a target mesh. The two meshes must already be « relatively » close to each other to obtain good results.

python3 icp.py -s SOURCE -t TARGET -m MATRIX[-mIts ITERATIONS][-mPts MAXPTS][-tol TOLERANCE]

SOURCE : .mesh file of the source (object to be transformed) TARGET : target .mesh file (file to transform to) MATRIX : .txt file in which to write the obtained transformation matrix ITERATIONS : maximum number of iterations, default 200 MAXPTS : maximum number of points on which to calculate the PCI (practice for large mesh sizes), default 5000 TOLERANCE : Convergence tolerance, default 1e-4

levelset.py

Extract the surface corresponding to a specified value from levelset (i.e. the reference surface 10 following extraction with mmg3d) for a tetrahedral .mesh file containing signed distance information.

python3 levelset.py -i INPUT -o OUTPUT[-d HAUSDORFF][-l LEVELSET]

INPUT : input .mesh file, must match with a mshdist output OUTPUT : .output mesh, will only contain the surface HAUSDORFF : Hausdorff distance, default at 0.01 LEVELSET : value for the levelset, default at 0 mask.py Generates a volumetric mask from an inner and outer surface.

python3 mask.py -i INTERIOR -e EXTERIOR -o OUTPUT [-t TEMPLATE]

INTERIOR : inner surface EXTERIOR : outer surface OUTPUT : output file TEMPLATE : .mesh file which will allow to write an assigned .sol inside the mask (Dirichlet) corresponding to the displacement field between the inner surface and this file. Note that the interior and exterior surfaces must be manifold, and therefore have no intersections. For the writing of Dirichlet's conditions, the .sol to be written should actually contain the morphing results.

merge.py

Merges several .meshe files into one. python3 merge.py -i FILES -o OUTPUT FILES : list of .mesh files to merge, supports wildcards OUTPUT : .mesh merged file

morph.py

Initiates the morphing of a template mesh to a signed distance.

python3 morph.py -s SOL-FILE -t TEMPLATE -o OUTPUT SOL-FILE : .sol file target (signed distance) for the morphing TEMPLATE : .mesh template file to deform towards the signed distance OUTPUT : .mesh morphed file

pca.py

Initiates a PCA to reconstruct an « unknown » mesh based on a set of meshes. Several actions can be performed with this script : compute the coefficients only, or rebuild a mesh in an existing database using n components. By default, the script computes the PCA on the input meshes, reconstructs the unknown mesh using the first 15 components. INPUT : .mesh file representing the surface to which we want to compute the distance. OUTPUT : .mesh file that will contain the result VOLUME : .mesh file containing tetrahedra and on which the signed distance will be calculated.

Note that by default, a volumetric mesh enclosing the object is created and adapted to the target object. This then allows mshdist to compute a signed distance on a mesh whose elements are smaller around the surface. If a volume is provided, it is used as a reference mesh, but is still remeshed to fit the surface. It would therefore be necessary to modify the script to directly calculate the signed distance on the volume provided, which would actually be equivalent to simply running mshdist.

split.py

Cuts a .mesh file in its middle on the X axis, keeping the part of the positive or negative x's. Based on blender to make a "clean" cut, this script is different from the others in that it must be executed via blender from the command line.

blender -background -python split.py --i INPUT -o OUTPUT -x X INPUT : .input mesh to be cut in its center (must therefore be centered in X) OUTPUT : output file to write X : +1 (to keep the x-positive part and make a symmetry) or -1 (left part, no other transformations)

Be careful, the files thus created are not centered in 0,0,0 !

transform.py

Can be used to : scale to 1 (the maximum size of the object will be 1) and center an .mesh to 0.5,0.5,0.5 , simply scale and/or move a .mesh file, apply a transformation matrix to a mesh (e. g. from the ICP calculation).

python transform.py -i INPUT -o OUTPUT [-s sx sy sz] [-t tx ty tz] [-c] [-m MATRIX]

INPUT : Input .mesh file OUTPUT : transformed .mesh file sx, sy, sz : scale factor according to x, y and z tx, ty, tz : translation according to x, y and z c : Boolean flag to be used alone to scale 1 and translate an .mesh into 0.5,0.5,0.5 MATRIX : .txt file containing a transformation matrix to be applied to the object warp.py Performs the warping of a template sphere on another object, with the objective of closing this object. python3 morph.py -i FILE -t TEMPLATE -o OUTPUT FILE : .mesh file corresponding to the object to warp TEMPLATE : .mesh template file to be deformed (a sphere) OUTPUT : .mesh warped file