Correlated Noise. CRF Conditional Random Field. DAG

Keywords: Spatial Bayes Network, STMT Spatial Triplet Markov Tree, UGM Undirected Graphical Model, UN Uncorrelated Noise AUC Area Under Curve, Ba Background, BHC Beam Hardening Correction, BIC Bayesian Information Criterion Convolutional Neural Network, CT Computed Tomography, CVD Cardiovascular Disease Magnetic Resonance Imaging, NC Nodular Calcification, PET Positron Emission Tomography, PR Precision Recall, pyIS pyImSegm, RNN Recurrent Neural Network, ROC Receiver Operating Characteristic Stochastic Expectation Maximization, SH Specimen Holder, SPD Symmetric Positive Definite, SPE Stochastic Parameter Estimation, ST Soft Tissue, T-Gibbs Tempered-Gibbs, UD Upward Downward, VI Variational Inference

In the morning all the birds are chirping and everything is going on. I grab a coffee and get to work. Then when I get tired, I take off. If I run out of ideas, I go for a run, and things fall into place.

-Bernd Heinrich, interviewed by Runner's World. Je te tiens tout d'abord à remercier Isabelle Bloch et Cédric Richard d'avoir accepté d'être les rapporteurs de mon travail. Je remercie également Wojciech Pieczynski d'avoir présidé mon jury de thèse. Merci à Mauro Gargiulo (Professeur à l'Université de Bologne) d'avoir honoré le jury de sa présence, en qualité de membre invité. Un grand merci à mes directeurs de thèse, Christophe Collet et Nabil Chakfé, et à mon encadrant de thèse, Emmanuel Monfrini. Depuis le début, ils m'ont fait confiance, m'ont encouragé et soutenu. Merci de m'avoir ouvert les portes des laboratoires ICube et Geprovas pour ces trois intenses années de thèse qui viennent de s'écouler... Merci à Jean-Baptiste Courbot, qui m'a transmis une partie de ses connaissances et de son expérience d'enseignant-chercheur, en donnant beaucoup de son temps et de son énergie.

Merci à ma famille. En particulier, la présence sans faille de mes deux mamies et de mes parents, depuis le tout début, a été essentielle à la réussite de cette thèse. Merci à ma soeur, qui démontre que le succès peut prendre une infinité de formes.

Finalement, je veux dire merci à toutes les personnes que j'ai croisées pendant ces dernières années, je suis certain d'avoir appris quelque chose de chacune de nos rencontres. C'était un plaisir de passer du temps à vos côtés et d'avancer, dans la thèse ou dans toutes les autres aventures qu'offre la vie. Donner une liste de noms me semble risqué et grandement insatisfaisant ; j'ai confiance dans le fait que vous vous reconnaîtrez...

List of Notations

Graph Theory

G A graph. E
Set of edges of a graph.

S

Set of vertices of a graph. N s Neighborhood of site s. P(s) Set of parents of s. s - A parent vertice of s.

S

Set of vertices with at least a parent vertice. s s s s s s s ++ The set of all descendants of s. s s s s s s s --The set of all ascendants of s.

xiii

Résumé étendu

Introduction

Les maladies cardiovasculaires sont la première cause de mortalité dans le monde. Les coûts de ces maladies pour les sociétés sont humain et économique, ils sont déjà importants (chiffrés en plusieurs centaines de milliards de dollars) et ne cessent de croître 1 2 . Le vieillissement et l'augmentation de la population mondiale sont des facteurs contribuant à l'augmentation de ces pathologies. Pour y faire face, la recherche est poussée à mieux connaître ces maladies et développer de nouveaux traitements plus efficaces.

Dans ce contexte, les approches modernes en chirurgie vasculaires sont basées sur des opérations mini-invasives guidées par l'image et sur l'implantation de biomatériaux (dont l'élément le plus connu est le stent). On parle de chirurgie endovasculaire. Ces nouveaux types d'intervention offrent une alternative à la chirurgie ouverte classique. Cependant, la chirurgie endovasculaire, déjà largement pratiquée, soulève des questions auxquelles la recherche n'a actuellement que des réponses trop partielles. Ces questions concernent aussi bien le comportement mécanique des biomatériaux dans le corps du patient et leur implantation par le chirurgien, que la correction d'images scanner illisibles et l'utilisation per-opératoire de ces images.

Les problématiques présentées ci-dessus relèvent d'une recherche fortement pluri-disciplinaire. Dans cette thèse, nous présentons des contributions aux traitement des images médicales. Ces contributions ont pour but de permettre la création de nouveaux outils qui doivent servir à d'autres chercheurs (biomécaniciens, chirurgiens, spécialistes du textile, radiologues, histopathologistes, ...) pour, in fine, mieux connaître les maladies cardiovasculaires, leur traitement et leur prévention. Nos contributions sont liées à une principale application : celle de la segmentation des images médicales.

Les images que nous traitons dans la thèse sont de type Computed Tomography (CT) ou micro-Computed Tomography (mCT). Les principaux défis dans la segmentation des images médicales en chirurgie vasculaire concernent la complexité et les dégradations subies par images à rayons-X lorsque celles-ci contiennent un biomatériau métallique, la disponibilité limitée des données et la taille des données à traiter. Dans la section suivante, nous décrivons plus en détails ces problèmes, en mettant en avant les méthodes graphiques probabilistes utilisées pendant la thèse pour les traiter. La des données de chirurgie vasculaire que nous souhaitons pouvoir segmenter.

Méthodes

Dans les modèles graphiques probabilites appliqués à la segmentation des images, des sommets du graphe sont associés à des variables aléatoires qui représentent l'image observée, et d'autres sommets sont associés à des variables aléatoires qui correspondent aux classes dans l'image segmentée. Les arêtes du graphe représentent une relation de dépendance entre les deux variables reliées.

Les arêtes peuvent être orientées ou non-orientées ce qui modifie la formulation probabiliste du modèle [START_REF] Mourad | A survey on latent tree models and applications[END_REF].

Les modèles graphiques probabilistes les plus connus pour la segmentation sont ceux de la famille des modèles de Markov cachés [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF]. Parmi les modèles de Markov cachés les plus populaires, nous notons les chaînes et arbres de Markov [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF] [START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF] (modèles graphiques orientés) et les champs de Markov (Kato and Zerubia, 2012) (modèle graphique non-orienté). Nous étudions largement ces modèles et leurs extensions : les modèles de Markov cachés couples et triplets [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] [START_REF] Lanchantin | Unsupervised segmentation of randomly switching data hidden with non-Gaussian correlated noise[END_REF] [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF]. Les modèles probabilistes markoviens offrent un moyen simple et intuitif d'introduire de la dépendance entre les pixels de l'image, ils sont une réponse pertinente aux problématiques posées par les images étudiées pendant la thèse.

Nous passons maintenant en revue les principales notions théoriques sur les modèles probabilistes graphiques, à la lumière des problématiques propres à notre application.

Modèles probabilistes et segmentation d'images bruitées

Une problématique particulière à laquelle nous devons faire face est celle des artéfacts. Lors de l'acquisition des images à rayons-X, les biomatériaux mé-talliques présents dans le corps du patient vont intéragir avec les rayons-X ce qui va créer, à la sortie de l'algorithme de reconstruction CT, de forts artéfacts. Ces derniers empêchent de discerner facilement l'anatomie environnante. D'un point du vue du traitement du signal les artéfacts peuvent être modélisés comme du bruit spatialement corrélé. Ce type bruit peut être naturellement modélisé dans certains modèles graphiques probabilistes comme les nouveaux modèles de Markov couples et triplets [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF] qui sont présentés dans la thèse. En parallèle, les corrélations spatiales sont classiquement étudiées dans les modèles de type Gaussian Markov Random Fields (Rue and Held, 2005) (GMRF) qui sont également étudiés dans notre travail.

Modèles probabilistes et disponibilité des données

Les biomatériaux sont récents et leur étude par l'imagerie après explantation n'a débutée que récemment. Ainsi, la quantité de données est encore relativement restreinte et les données sont, pour la plupart, non-annotées 3 . Dans de tels cas, les approches non-supervisées sont privilégiées, voire nécessaires. Les modèles probabilistes graphiques sont souvent utilisés car ils bénéficient d'algorithmes relativement efficaces dans les cas où les données sont rares ou manquantes. En contexte non-supervisé, les paramètres des distributions de probabilité des variables aléatoires du modèle graphique doivent être estimés à l'aide seulement de l'observation. Dans ce contexte, les approches de type Expectation-Maximization [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] (EM) qui ont pour objectif de maximiser la vraisemblance des données sont les plus connues. Plusieurs versions de ces algorithmes existent et elles sont déclinables pour les modèles graphiques orientés et non-orientés [START_REF] Celeux | On Stochastic Versions of the EM Algorithm[END_REF] [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]. Nous nous intéressons, dans notre travail, à ce type d'approches non-supervisées et leur procédure d'estimation de paramètres.

Le cas où une base de données annotées est disponible (approche supervisée) est aussi étudié dans la thèse. La segmentation des images bidimensionnelles est faite par un réseau de neurones convolutionnel [START_REF] Chetan | Deep neural network models for computational histopathology: A survey[END_REF]. Ces approches qui relèvent de l'apprentissage profond sont introduites dans la thèse. En effet, elles donnent les meilleurs résultats dans un grand nombre de problèmes de segmentation supervisée et sont devenues incontournables dans le domaine. En revanche, l'étape ultérieure qui correspond à la reconstruction d'images segmentées tridimensionnelles est plus complexe. Elle nécessite des approches plus fines, notamment des associations de modèles d'apprentissage profond et de modèles probabilistes graphiques, que nous étudions également dans le manuscrit [START_REF] Ben-Cohen | Fully convolutional network for liver segmentation and lesions detection[END_REF] [START_REF] Kamnitsas | Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation[END_REF] [START_REF] Alexey | Deep sequential segmentation of organs in volumetric medical scans[END_REF].

Modèles probabilistes et le coût de l'inférence

Les images de chirurgie vasculaire de type mCT (qui offrent une bien meilleure résolution spatiale que les images CT) sont de très grande taille. Les tâches d'inférence dans les modèles graphiques probabilistes associés à des images de grande taille peuvent être très coûteuses en temps de calcul, voire même infaisables.

De nombreuses approches sont alors fondées sur une approximation de l'étape d'inférence. Par exemple, la technique de l'inférence variationnelle approche la distribution de probabilité pour laquelle l'inférence est coûteuse par une distribution simplifiée. Les recherches sur cette approche sont aujourd'hui très actives (C. Zhang et al., 2018) et l'inférence variationnelle est étudiée en détail dans notre travail pour des modèles graphiques orientés et non-orientés.

Les champs de Markov cachés sont les modèles graphiques probabilistes les plus populaires pour la segmentation non-supervisée d'images médicales. Or l'inférence dans les modèles graphiques non-orientés (auxquels appartiennent les champs de Markov) ne peut se faire, sauf dans des cas particuliers, de manière directe. L'inférence dans ces modèles repose sur des calculs indirects (qui peuvent être approximants ou non) comme l'approche itérative de l'échantillonneur de Gibbs (S. Geman and D. Geman, 1984). Les chaînes de Markov cachées peuvent alors être une alternative judicieuse car leur structure intrinsèquement unidimensionelle et les possibilités d'inférence exacte offrent des coûts calculatoires faibles relativement aux approches classiques d'inférence dans les champs. Ces propriétés sont rencontrées dans [START_REF] Bricq | Unifying framework for multimodal brain MRI segmentation based on Hidden Markov Chains[END_REF] [START_REF] Courbot | 2-step robust vertebra segmentation[END_REF] et vues dans cette thèse également.

Contributions

Après avoir introduit les principales notions sur les modèles graphiques probabilistes, nous donnons ici les principales contributions de la thèse.

Le modèle Gaussian Pairwise Markov Fields

Nous présentons un nouveau modèle de champs de Markov couples cachés, nommé Gaussian Pairwise Markov Fields (GPMF), qui permet la segmentation non-supervisée d'images corrompues par du bruit spatialement corrélé modélisé par un GMRF (Rue and Held, 2005). Nous étudions, théoriquement et expérimentalement, dans quelle mesure ce modèle est une généralisation du modèle de champ de Markov caché classique. Nous proposons également un algorithme stochastique itératif pour l'estimation non-supervisée des paramètres. Le modèle est utilisé pour segmenter la matière organique dans des images mCT touchées par des artéfacts de stents. La Figure 0.3 illustre cet axe de recherche.

Le modèle Spatial Triplet Markov Tree

Dans cette thèse, nous développons également un nouveau modèle d'arbre de Markov triplet nommé Spatial Triplet Markov Tree (STMT) basé sur (Courbot, Monfrini, et al., 2018). Ce modèle intègre des variables aléatoires auxiliaires afin d'enrichir les possibilités de modélisations. Le modèle STMT est une généralisation des arbres de Markov caché classique.

De plus nous étudions ses liens avec les champs de Markov cachés classiques. Le nouveau modèle semble en effet exhiber des corrélations semblables aux champs mais offre des possibilités d'inférence exacte grâce à la structure d'arbre. Une telle propriété a pour principal avantage d'éviter le recours à des procédures itératives souvent approximantes, plus longues et induisant des pertes de précision. Nous étudions les corrélations à l'intérieur du nouveau modèle STMT avec la technique de l'inférence variationnelle à variables aléatoires auxiliaires. Le modèle STMT est aussi utilisé dans le cadre de la segmentation de calcifications sur des images mCT. La

Inférence variationnelle dans des modèles plus complexes

Nous étudions des modèles plus généraux dans le sens où plus de dépendances directes sont modélisées entre les variables aléatoires.

Pour les modèles graphiques orientés, alors que tous les modèles de Markov cachés classiques sont associés à des graphes ne présentant ni cycles ni semicycles, nous proposons un nouveau modèle qui contient des semi-cycles appelé Spatial Bayes Network. La représentation graphique de ce modèle est donnée Figure 3.9. L'introduction de semi-cycles rend l'inférence beaucoup plus complexe et nous résolvons ce problème en ayant recours à l'inférence variationnelle.

Pour les modèles graphiques non-orientés, il est possible de complexifier le modèle en ajoutant des dépendances conditionnelles entre les variables qui agrandissent les voisinages. Le modèle des fully-connected Conditional Random Fields est un modèle probabiliste populaire en segmentation d'images qui est associé à un graphe totalement connecté (voir Figure 0.5b). Nous étudions à nouveau plusieurs techniques d'inférence variationnelle afin de pouvoir mener l'inférence dans ce modèle.

Segmentation de stents

Dans les images CT, les stents peuvent apparaître fortement déformés à cause des artéfacts qui les entourent. Nous proposons une modélisation basée sur les chaînes de Markov cachées et un modèle de bruit particulier qui permet de segmenter finement le stent dans un environnement complexe (tel que celui de la Figure 0.2). Notre approche opère de manière non-supervisée et est entièrement automatique grâce à l'algorithme d'estimation des paramètres que nous développons également. Une illustrations des résultats de cet algorithme est donnée dans la Figure 0.6.

Segmentation histologique tridimensionnelle d'artères

Dans une dernière partie, nous mettons au point un protocole pour la création de la première base de données d'images annotées d'artères explantées touchées par l'athérosclérose. Ensuite, nous développons un réseau de neurones convolutionnel de type U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] afin de procéder à la segmentation bidimensionnelle histologique des images. L'objectif de ce projet unique est de permettre une première analyse histologique de l'artère uniquement à l'aide du scanner. Une reconstruction tridimensionnelle propre est obtenue grâce à un posttraitement basé sur une inférence variationnelle dans un modèle de champs aléatoires conditionnels [START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF]. Nous étudions également une amélioration de cette dernière procédure d'inférence variationnelle utilisant les chaînes de Markov. Un exemple de segmentation tridimensionnelle d'une artère est visible dans la Figure 0.7.

Conclusion

Les modèles développés au cours de la thèse apportent des contributions sur des enjeux cruciaux de la modélisation des images par modèles graphiques probabilistes tout en répondant à des problématiques modernes de la chirurgie vasculaire.

Introduction

Advances in vascular surgery and current issues

Cardiovascular Diseases

Cardiovascular Diseases (CVDs) are a major cause of mortality in the world, and particularly in developed countries, with an increasing human and monetary cost for societies. In 2016, approximately 17.9 million people died worldwide from CVDs, which represents 31% of all deaths 4 . In 2015, 49 million people (9.6% of the population) in the European Union where living with a CVD, the total costs are estimated to e210 billion with e111 billion of direct health care costs and e99 billion indirect costs (productivity loss and informal care of people) 5 . A major cause of the rise of CVDs is also the increasing lifetime: this causes new diseases, related to old age, to be seen more widely and needing to be treated. But CVDs are also often linked with life hygiene: smoking, sedentarity and nutrition are habits with direct effects on health in general and more particularly on the cardiovascular system 6 .

CVDs include a number of heart and blood vessel conditions. Among the most prevalent conditions, one can note ischemic heart disease, stroke and high blood pressure. These conditions often relate to another condition called atherosclerosis. Atherosclerosis is a disease of the arterial wall (not restricted to coronary arteries). It consists in the formation of plaques, called atheromatous plaques, which mostly contain lipids, macrophage cells, connective tissues and calcium (Rafieian-Kopaei et al., 2014). The deposition of calcium contributes to the plaque hardening and the sclerosis of the artery, but it is also what makes the plaque clearly visible on X-ray images. Therefore, one commonly refers to the atheromatous plaques as calcifications since the latter are the most easily identifiable elements of the plaques. Those plaques may grow and/or break which disturbs or, in the most severe cases, stops the blood flow (total occlusion of the artery, rupture of the artery, etc.). These situations are medical emergencies [START_REF] Delia | An overview of the evolution of the atherosclerotic plaque: from fatty streak to plaque rupture and thrombosis[END_REF]. However, in some cases, atherosclerotic plaques can remain stable and/or regress [START_REF] Dave | Plaque regression and plaque stabilisation in cardiovascular diseases[END_REF]. Much is yet to discover about this phenomenon which is summarized in Figure 0.8.

In this thesis, our research will be applied to pathologies affecting the arteries of the lower limbs, as well as their surgeries and their treatments, which are described in the next section. More precisely, we will focus on the superficial femoral artery and the popliteal artery which are located in the human body in Figure 0.9. The choice of the femoropopliteal arterial segment was made because it is still the most challenging area for treatment, it is also the most often affected segment in the arterial tree.

Remark: Strictly speaking, cardiovascular refers to the heart and its main arteries and veins. Hence, we will use the adjective vascular to refer to the arteries of the lower limb. However, the problematics and pathologies that affect the vessels are identical or very similar.

Remark: Veins will not be mentioned in this thesis since they do no develop atherosclerosis (G. [START_REF] Saul | Physical fitness, dynamic extra-arterial pressures, and the pathogenesis and distribution of atherosclerosis[END_REF].

Vascular surgery: open and endovascular approaches

Vascular surgery began with the first attempt to control a bleeding vessel. (S. G. [START_REF] Friedman | The elements of statistical learning[END_REF] The traditional approach in vascular surgery, often called open surgery, relies in opening the patient body up to the affected vessel. The vessel is then repaired or replaced. This approach can lead to high-risk interventions and might not be tolerated by already weak patients such as old patients or patients with a chronic illness. In such cases, the latter have increased chances of complications, including death, following the intervention.

A more recent approach is endovascular surgery, in which new interventions and protocols are regularly created since the last two decades of the 20 th century. In this interventional technique, the vessel is reached and treated from the inside. With an incision in a peripheral artery, the surgical tools and treatments are brought to the desired location, by the inside of the blood vessel, following the patient's vascular tree. The treatment often includes a biomaterial. Biomaterials are introduced in the next section. Figure 0.9.: Arteries of the lower limbs (anterior view). From [START_REF] Openstax | [END_REF].

Endovascular interventions rely heavily on imaging to see inside the patient which is not opened: the terms image-guided surgery and mini-invasive are then often associated with endovascular surgery. This surgical technique is often a solution proposed to patients who cannot undergo open surgery. The risk of medical complications and the patient's traumatism following an endovascular procedure is reduced [START_REF] Kirt | Has open vascular surgery disappeared?[END_REF]. While endovascular techniques have classically been linked with a higher monetary cost than open surgery procedures, the gradually decreasing cost of the treatment and of the patient's length of stay becomes an additional asset for endovascular procedures [START_REF] Stanke | Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources[END_REF] [START_REF] Clayton | Financial viability of endovascular aortic repair in the modern era: a single center experience[END_REF].

Note that interventional techniques, either from endovascular or open surgery, are constantly improved and updated. A complete reference to endovascular history, techniques and biomaterials is available in (Jing et al., 2018).

Endovascular surgery is a modern approach which raises new questions on biomaterials and new challenges on the per-operative use of imagery during surgical interventions and at all the other stages of the patient's treatment. Those questions are the motivations behind this thesis. Biomaterials are found in many medical fields: knee or hip prosthesis in orthopedic surgery, intraocular lens in ocular surgery, heart valve in cardiac surgery, etc. We focus here in biomaterials in vascular surgery. Following (Jing et al., 2018), we list the most common biomaterials:

Biomaterials and related problematics

• Stents are metallic biomaterials widely used to solve problems such as stenosis or dissection of the artery. Once they have been brought to the location of the lesion via the intravascular path, they are released. They are either self-expandable or expanded with a surgical balloon. A stent is illustrated in Figure 0.10a.

• Stentgrafts are distinguished from the (bare) stents described above by the covering (most of the time in polytetrafluoroethylene or polyethylene terephthalate 7 ) added over the stent. Stentgrafts are used mainly in aneurysm repairs or aortic dissections for their capacity to seal the vessel wall. The Figure 0.10b depicts a stentgraft.

• In the most severe cases of lesions, such as complete stenosis or occlusion, the vessel can be replaced with vascular protheses. They may originate from an autogenous vein (the ideal case for biocompatibility) or from an artificial prosthesis which can again be made of PTFE or PET. Note that, as opposed to stents and endografts, an implantation of a vascular prosthesis involves an open-surgery procedure. However they offer better performance and stability than stentgrafts in terms of biocompatibility, fracture resistance and ability to seal the vessel wall. A vascular prosthesis can be seen in Figure 0.10c.

To treat a given lesion, the choice of the model and size of the biomaterial depends mainly on the patient's imaging analysis and the surgeon's experience. Remark: In this thesis, we will mainly work with medical images of explants. In vascular surgery, explanting is the act of removing an implanted biomaterial or a part of a blood vessel. Note that when a biomaterial is explanted, it is almost always encapsulated in some biological material. An explant can be seen in Figure 0.11.

Problematics

Biomaterials in vascular surgery are recent treatments. While these devices benefit from constant improvements and are increasingly implanted, many questions remain unanswered:

• The mechanical behavior of the implanted stents are almost unknown.

Simulations, or in vitro experiments, only partially reflect the reality of the constraints applied on a biomaterial implanted in a living body. In vivo mechanical data are still very much ignored and biomaterial manufacturers would greatly take advantage of such data. For example, the optimal stent design is still much debated [START_REF] Raffort | Atherosclerosis: process, indicators, risk factors and new hopes[END_REF].

• The choice of the endovascular procedure is still very dependent on the experience of the surgeon. Objective mathematical tools and databases to guide clinicians are still lacking [START_REF] Ohana | Detailed crosssectional study of 60 superficial femoral artery occlusions: morphological quantitative analysis can lead to a new classification[END_REF].

• In many cases, the cause(s) of failure of the vascular biomaterials are ignored or unclear. Why did the stentgraft covering tear? Why did the stent fracture? This highlights the need of new analysis tools for explanted biomaterials to investigate biomaterial failures [START_REF] Chakfé | What do we learn from explant analysis programs?[END_REF]) [START_REF] Lejay | Explanted vascular and endovascular graft analysis: where do we stand and what should we do?[END_REF].

• Although the whole branch of endovascular surgery relies heavily on images, those medical images remain widely unused. Indeed we could benefit from deeper, large-scale and automated analyses, carried by new tools that researchers in image processing could provide [START_REF] Raffort | Atherosclerosis: process, indicators, risk factors and new hopes[END_REF]. The perspectives of image processing for vascular surgery are the Introduction central motivation for this thesis and are discussed in depth later in this introduction.

We have only listed a small sample of questions arising today in the field of vascular surgery. Yet, addressing these questions already requires joint work from research teams from many fields (biomechanics, image processing, vascular surgery, textile sciences, artificial intelligence, histolopathologists, etc.). However, this translational research is also dependent on the availability of data (explanted biomaterials, explanted arteries, patient's clinical images, patient's clinical data, antecedents and outcomes, etc.). Since this research involves human subjects on a new topic, the data is scarce and sometimes unavailable due to legal issues. Availability of the data is a crucial point in the field of image processing as we will see later.

The Geprovas laboratory

The Geprovas11 laboratory (Groupe Européen de Recherche sur les Prothèses liées à la Chirurgie Vasculaire), located in Strasbourg, France, has been founded in 1993 by Professor Nabil Chakfé. It has become a worldwide actor in the field of biomaterials, with an original and unique expertise for devices from vascular surgery. The Geprovas is organized around four activities: the collect and analysis of explants, the innovative research on biomaterials, the medical education program and the clinical analysis program. One of the main goals of Geprovas is to remove barriers and motivate translational research to help answer questions in medical research such as the questions we listed in the previous section.

The work presented in this thesis is initiated and supported by the explant analysis program of the Geprovas. In this context, our work benefited from the wide knowledge of the Geprovas experts as well as from its unique database of explants from vascular surgery.

Medical image processing and new challenges

In the previous section, we mentioned medical imaging and its key role in the modern approaches in vascular surgery. In this second part of the introducion, we present medical imaging in more details. We start with very general definitions and gradually shift our focus towards the mathematical methods linked with image analysis studied in this thesis. The goal is to develop approaches which are automated. Automated analysis is indeed the modern way to process images relying on the power of computers to support medical practice and research.

Medical image processing

Main concepts and vocabulary

A bidimensional (respectively tridimensional) image is a set of values arranged on a rectangular (respectively cubic) grid12 . For numerical images, the grid is discretized. If the values are scalars, the image is called grayscale, if the values are vectors, the image is a color image in some color space [START_REF] Fieguth | Statistical image processing and multidimensional modeling[END_REF].

The most common medical image processing operations are image segmentation, which consists in dividing an image into non-overlapping regions with homogeneous properties, image registration, which deals with associating the objects/features in one image with those in one or more other images, and image denoising, which aims at estimating a noise-reduced image given an image corrupted by noise [START_REF] Fieguth | Statistical image processing and multidimensional modeling[END_REF]. Image segmentation is the task studied in this thesis and all the methods we develop can be linked with the final goal of segmentation.

A numerical image has to be constructed with some devices, in medical image processing, it is related to the notion of image modality.

Medical imaging modalities

The observed image that we have to process can be of different types in medical research. Those are called the imaging modalities: they are acquired with different devices relying on different physical phenomena to form the image. Thus, the imaging of the same physical object will result in different images with different properties according to the modality. Following [START_REF] Suetens | An introduction to conditional random fields[END_REF], the main medical imaging modalities for diagnosis and treatment13 are now listed:

• X-ray Computed Tomography (CT) is an imaging technique based on Xrays. X-rays are electromagnetic waves with wavelengths around 10 -10 m which are attenued differently according to the matter it interacts with. The computation of the attenuations undergone by X-rays at certain locations of the space are at the foundation of X-ray images. Computed Tomography (CT) refers to the technique of the computations of the attenuations. X-rays have been used in medical imaging starting from the end of the 19 th century. Figure 0.12a shows a typical CT scan image from vascular surgery.

• Magnetic Resonance Imaging (MRI) is an imaging modality introduced in the medical field in the 1970s. In MRI, the reconstructed image illustrates the magnetic properties of the object of interest. In normal operating conditions, MRI is safer for patients than CT imaging or PET imaging (see next item) because MRI does not rely on ionizing radiations.

• Positron Emission Tomography (PET) is a technique from the field of nuclear medical imaging. It uses a tracer molecule that will be involved in a metabolic process of the body. The tracer molecule is injected to the patient and can be localized because some identified radioactive atoms of the molecule emit γ-rays (wavelength below 10 -11 m). This way, an abnormal metabolism can be detected. PET is used clinically since the 2000s.

Note that some modalities describe the anatomy, this is called structural imaging and this is the case for most CT and MRI performed. Others depict a function, this is called functional imaging; PET is a functional imaging modality. In more specific contexts, a modality can be ambivalent and relate to both structural and functional imaging.

In this thesis, we will only develop approaches for X-ray images. Indeed, X-ray images are more common than any other modalities in vascular surgery notably because current MRI devices can be unsafe for patients with a metallic implant. Indeed, MRI might induce stent dislodgement, heating and important image artifacts [START_REF] Jabehdar Maralani | The one hundred layers tiramisu: Fully convolutional densenets for semantic segmentation[END_REF]. X-ray images have also the advantage of being a totally non-invasive technique14 , the body of the patients does not need to be touched.

Our work will be more particularly focused on the potential assets of microcomputed tomography (mCT or microCT) X-ray images used for the research in vascular diseases. The formation of mCT images follows the same principles as CT X-ray images described above but they are suited for the imaging of small objects (animals15 , stones, wood, etc.). Notably, mCT modality deals particularly well with the imaging of explanted arteries. Those images offer a much greater spatial resolution than classic CT scans since pixel sizes are on the scale of the micrometer [START_REF] Flannery | Estimation of space deformation model for non-stationary random functions[END_REF]. The first report of a microtomographic image dates back to 1982 in (J. C. [START_REF] Jim | U-Net: deep learning for cell counting, detection, and morphometry[END_REF]. Figure 0.12b shows a mCT image of an arterial cross-section.

Remark: We will also mention optical microscopy, which is an image modality that uses light to image a very thin section of an object of interest, which often cannot be analyzed or seen by the naked eyes [START_REF] Michael | Simulated annealing: From basics to applications[END_REF]. Figure 0.12c describes an arterial cross-section seen in optical microscopy. The work in this thesis will not deal with a direct processing such images, however microscopic images will be essential in the project described in Chapter 5.

Main challenges

In this section, we describe the difficulties linked with medical image processing that motivate the approaches we developed during the thesis.

Metallic biomaterials and artifacts

We have discussed above that metallic stents (and metallic biomaterials in general) are an issue for MRI scanning but some problems also arise in CT scans. Due to high variations of attenuations in metallic parts according to the X-ray energy level, the CT reconstruction algorithms introduce artifacts in the reconstructed images. The artifacts are problematic since they hide the underlying anatomy. Figure 0.13 illustrates stent metallic artifacts in vascular X-ray scans. This issue is becoming increasingly important because of the development of biomaterials (orthopedic prostheses, stents, pacemakers, etc.).

To answer this problem, many approaches to Metal Artifact Reduction (MAR) have been developed in the literature. See (H. S. [START_REF] Hyoung | Computed tomographic beam-hardening artefacts: mathematical characterization and analysis[END_REF] for an introduction on stent metallic artifacts and MAR.

Availability of the data

In some medical applications data are scarce while in other fields data abound. This variable then plays a role on the processing approach to follow and on the problem complexity. In vascular surgery, biomaterials were usually destroyed after explantation or patient's death and relatively little research could be carried. As stated earlier, the goal of the explant analysis program of the Geprovas is precisely to collect such data, and make it available in a database that can be used in subsequent studies. An image analysis task in which no images can be used to learn the algorithm parameters, i.e., where the algorithm has to work with the observed image alone, is called unsupervised. In the opposite case, it is called a supervised task. Both approaches will be seen in the thesis. Note that a supervised problem requires a human interaction at some stage, for example, to add annotations to the data.

Remark: Collections of biological samples used for scientific investigation, or biobanks, are initiatives that are slowly gaining popularity [START_REF] V Pw Scholtes | Biobanking in atherosclerotic disease, opportunities and pitfalls[END_REF].

Computational costs

Image processing can rapidly become a computationally intensive task, especially for medical images. We encounter this problem in our work dealing with mCT images which are big, high resolution images. In other problematics, when a database with hundreds of elements must be processed, the treatments must be wisely chosen. In this thesis, some of the works study the problem of computational complexity.

Remark: Our developments are however not constrained by speed considerations. Our primary data are images of explanted materials, thus the time constraint is not comparable at all to developments linked with videos or realtime analyses such as [START_REF] Innocent Nwoye | Weakly supervised convolutional LSTM approach for tool tracking in laparoscopic videos[END_REF].

Image processing and probabilistic modeling

In this thesis, the main contributions will be focused on probabilistic graphical models applied to image segmentation, particularly Hidden Markov Models (HMM) [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF] and their extensions to Pairwise and Triplet Markov Models [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] [START_REF] Lanchantin | Unsupervised segmentation of randomly switching data hidden with non-Gaussian correlated noise[END_REF]) [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF]. Markovian models offer a simple and intuitive way to introduce dependencies between the pixels of the images and thus, they can be a relevant answer to the problematics previously presented.

For example, it will be seen that stent metallic artifacts (illustrated in Figure 0.13) can be modeled by spatially correlated noise: a pixel value at a location is highly dependent on the neighboring pixel values. Such a noise is naturally taken into account in Pairwise and Triplet Markov Models. Such correlations are also studied in models based on Gaussian Markov Random Field (Rue and Held, 2005) (GMRF). In this thesis, these models are then studied as possible improvements to MAR algorithms and as a way to perform segmentation in images degraded with spatially correlated noise.

A second asset to probabilistic models is that they are relatively efficient in unsupervised problems. Following the popular Expectation Maximization [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] (EM) algorithm, there is a rich literature of learning algorithms in unsupervised context for graphical models [START_REF] Celeux | On Stochastic Versions of the EM Algorithm[END_REF] [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]. They are then competitive methods for problems with scarce and/or missing data.

Probabilistic models also benefit from a very active research community focusing on optimization methods. Efficient approaches exist to approximate computations that cannot be done exactly. See for example (C. Zhang et al., 2018) which deals with Variational Inference (VI). They also offer much modelization power with possibly auxiliary random variables (Lanchantin et al., Introduction 2011) and can then fit many complex real life problems (Courbot, Monfrini, et al., 2018).

Probabilistic graphical models will be introduced in depth in Chapter 1.

Outline of the thesis

This thesis is composed of five chapters and their appendices. Chapter 1 consists in a literature review of the main concepts and developments around probabilistic modeling with graphical models with a special focus on image segmentation. This chapter sets up all the definitions and notations that will be used more in depth in subsequent chapters. Chapters 2, 3 and 4 present the core research of this thesis dealing with the improvement and development of new probabilistic models for image segmentation. Chapter 2 proposes a new model of Markov Fields, called Gaussian Pairwise Markov Fields (GPMF), for the unsupervised segmentation of regions that are corrupted with long-spatially correlated noise. The model is tested on both artificial and real datasets. Chapter 3 describes a new model of Markov Trees, called Spatial Triplet Markov Tree (STMT), aiming at being a deterministic counterpart to Markov Fields, which are also used in the context of unsupervised image segmentation. Theoretical results and applications on real data are provided. STMTs make use of auxiliary random variables, the enhanced correlations they introduce are largely studied. Chapter 4 proposes a review on Deep Learning for medical imaging and describes the development of a Convolutional Neural Network for semantic segmentation. This segmentation is subsequently improved by a new approximate inference approach for a model of Conditional Random Fields that is also presented in Chapter 4. Chapter 5 illustrates the applications to vascular surgery of the probabilistic graphical models developed along the thesis.

The appendix of each chapter gives the full details of calculations and additional observations.

List of publications

During this thesis, several journal and conferences articles have been written. They are listed below according to the part of the thesis they deal with. For completeness, the other publications from the author of the thesis are also listed even if their topic is not directly related to the content of this thesis.

First we list the publications primarily directed to the signal processing community:

• Chapter 2:

-Under review: Unsupervised Segmentation with Gaussian Pairwise Markov Fields, H. - (Gangloff, Courbot, et al., 2019) • Chapter 5:

- (Gangloff, Monfrini, Collet, et al., 2020): Unsupervised segmentation of stents corrupted by artifacts in medical X-ray images, H.

Gangloff, E. Monfrini, C. Collet, N. Chakfé, International Conference on Image Processing Theory, Tools and Applications, 2020.

- (Gangloff, Monfrini, Collet, et al., 2019): Segmentation de stents dans des données médicales à rayons-X corrompues par les artéfacts, H. Gangloff, E. Monfrini, C. Collet, N. Chakfe, Colloque GRETSI, 2019.

• Not in the thesis:

- [START_REF] Gorynin | Performance comparison across hidden, pairwise and triplet Markov models' estimators[END_REF] - [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF]: Assessing the segmentation performance of pairwise and triplet Markov models, I. Gorynin, H. Gangloff, E. Monfrini, W. Pieczynski, Signal Processing, 2018.

- [START_REF] Gangloff | Improved Centerline Tracking for new descriptors of atherosclerotic aortas[END_REF] We then present the publications primarily directed to the vascular surgery community: E) is formed by a set of vertices S, that we also refer to as sites or nodes, and by a set of edges E. We have E ⊂ S × S, since, in general, not all vertices are connected by an edge. The set of vertices connected to a vertice s ∈ G is called the neighborhood of s and is denoted N s . Note that we have s / ∈ N s . A clique is a subset of S with the property that every two elements are connected by an edge. A clique is then a fully-connected subset. A maximal clique of G is clique of G which cannot accept any more vertice from G without breaking the clique property. A tree is a connected graph without any cycle. In a graph G, the edges can either be directed, undirected or a mix of both. This leads to different families of probabilistic models which we discuss in the next sections. In the graphical representation of the probabilistic graphical models we will use arrows to represent directed edges. The absence of arrow represents an undirected edge. These first definitions are illustrated in Figure 1.1.

A graph G = (S,
In a probabilistic graphical model, each vertice s of G is associated with a random variable whose name then contains the name of the vertice as subscript, for example, X s .

Elements of probability theory

Let X be a random variable from a probability to space, equipped with a probability measure p, to a measurable space (E, E). If E is a finite or countable set, then X is said to be a discrete random variable. If E is an uncountable 1.2. Graphical modeling set, then X is said to be a continuous random variable. In both cases, we are interested in measuring the probability of X taking certain values, i.e., we want to evaluate the quantities of the type P X (A) = p(X ∈ A), ∀A ∈ E. Such computations involve P X another probability measure. P X is called the law of X which is more precisely called distribution in the discrete case and density in the continuous case. The definitions are similar for random vectors. When considering the law of random vectors made of both continuous and discrete random variables, we may use indistinctly the term distribution or density. If P X is the law of a random variable X, we may write X ∼ P X . If X is a discrete random variable, we denote by x its realization, ∀x ∈ E. When A is a singleton, i.e., A = x, ∀A ∈ E, we have p(X ∈ A) = p({X = x}). When there is no ambiguity we write p({X = x}) = p(X = x) or even p({X = x}) = p(x)1 to denote the probability of a realization x of X. Similarly, if X is a continuous random variable, in case of no ambiguity, we might use p(x) to refer to the density of X. If X ∼ P X , we denote the expectation of X by E[X], which can also be written, with a slight notation abuse, E x∼p(x) [x].

The conditioning of a probability law on the realization of some other random variable(s) will be written classically with a vertical bar |. For example, using the notation shortcuts mentioned above, the law of X given the realization y of random variable Y will be written p(x|y).

All our definitions extend to the multivariate case: random variables then become random vectors. Some of the random vectors we will work with are more precisely stochastic processes. However, stochastic processes and their properties are out of the scope of this thesis and little will be said about them.

Directed Graphical Models

Directed Graphical Models (DGMs), also known as Beliefs Networks, have directed edges represented by an arrow between two vertices. Let (s, s ′ ) ∈ S 2 , if there is a directed edge from s to s ′ then s is called the father node of s ′ , conversely, s ′ is called the child node of s. A node s may have several fathers, which are refered to as the set of parents of s, denoted P(s). A node for which the parent set is empty is refered to as a root node. Furthermore we denote S the set of vertices which have at least one father. Figure 1.2 illustrates the new notions we have just introduced.

The notion of Directed Acyclic Graph (DAG), also known as Bayesian Networks, refers to graphs where all the edges are directed and where the graphs do not contain directed cycles. The directed edges of a DGM leads to a partial ordering of its vertices (Wainwright and [START_REF] Martin | Graphical models, exponential families, and variational inference[END_REF]. A directed tree is a DAG which does not contain any directed cycles and not any semi cycles. A semi cycle is the DAG that has been obtained after reversing some arrow directions of a directed cycle.

In this thesis we will focus on several kinds of directed rooted trees, called arborescences. In such graphs the root vertice is unique, it is then the only vertice with no edge pointing towards it, and all the other edges are directed away from it. Note that in arborescences, the fact that there is only one root node implies that all the other non-root nodes have exactly one father. When all the vertices of an arborescence have one son (except for the last vertice), we obtain a directed chain. DGMs are specified by local conditional probabilities which are given for every random variables associated to every site: ∀s ∈ S the associated conditional probability is p(x s |x x x x x x x P(s) ). If r is a root node, the associated conditional probability is p(x r ). The joint probability distribution defining a DGM is obtained, using the definition of conditional probabilities, by multiplying the local conditional probabilities:

p(x s1 , x s2 , . . . , x s N ) = r∈S\ S p(x r ) s∈S p(x s |x x x x x x x P(s) ), (1.1)
where s 1 , . . . , s N are the N vertices of S.

Sampling in a DGM is performed by ancestral sampling: starting from the root nodes we follow the partial order over the nodes, and we sample each time the associated random variable using the local conditional probability [START_REF] Goodfellow | Deep learning[END_REF].

Graphical modeling

Conditional independence of random variables given other random variables (or Markov property) is an important property when dealing with graphical models. The procedure of d-separation can determine whether there is conditional independence [START_REF] Mourad | A survey on latent tree models and applications[END_REF]. A particular result of this procedure that will be underlying many of our studies is for arborescences: let A, B and C be subsets of S. ∀a ∈ A, ∀c ∈ C, the random variables {X a } a∈A and {X c } c∈C are said to be conditionally independent given B if every path (without considering the direction of the edges) between a and c goes through a vertice from B.

Undirected Graphical Models

Undirected Graphical Models (UGMs), also known as Markov Random Fields (MRFs), have undirected edges. They are defined differently from DGMs. We call potential functions, or factors or unnormalized probability distributions, a real non-negative function φ C ({x s } s∈C ) = φ C (x x x x x x x C ) associated with a clique C. Potential functions describe the interactions between random variables. We then have that an UGM is defined by a set of random variables X X X X X X X which admits for joint distribution:

p(x x x x x x x) = 1 Z C∈C φ C (x x x x x x x C ), (1.2)
where C is the set of all cliques of G and Z is the partition function or normalization constant, defined as:

Z = x x x x x x x C∈C φ C (x x x x x x x C ). (1.
3)

The partition function ensures that the distribution sums to 1. As a special parametrization, a probability distribution which takes the form of Equation 1.4 is known as a Gibbs distribution:

p(x x x x x x x) = 1 Z exp (-E(x x x x x x x)/T ) = 1 Z exp - 1 T C∈C ψ C (x x x x x x x C ) (1.4)
where E is called the energy function, ψ C are some potential functions and T is a parameter called the temperature (with T > 0). All the MRFs seen in this thesis will follow a Gibbs distribution. An equivalent definition for UGMs is possible due to the Hammersley-Clifford theorem [START_REF] Hammersley | Markov field on finite graphs and lattices[END_REF]. A set of random variables X X X X X X X, defined such that:

p(x x x x x x x) > 0, for all set of realizations X X X X X X X = x x x x x x x, p(x s |x x x x x x x S\{s} ) = p(x s |x Ns ), ∀s ∈ S, (1.5)
forms an UGM or MRF, with respect to the neighborhood N . The MRF definition of Equation 1.5 is called the full conditional definition. The Hammersley-Clifford theorem proves that the definitions of Equations 1.4 and 1.5 are equivalent. It then appears that the conditional independence property (Markov property) for UGMs is simpler than for DGMs. A random variable at site s,

Chapter 1. Main concepts in probabilistic modeling

given the realizations of the random variables located in its neighborhood N s , is conditionally independent from all the remaining random variables (i.e., from all the vertices in S \ N s ). Sampling realizations from UGMs is not as straightforward as for DGMs. Except from few exceptions where sampling can be done without approximations [START_REF] Stanke | Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources[END_REF], Markov Chain Monte Carlo (MCMC) approaches are the most popular way for this task. The Gibbs sampler is a widely used algorithm to sample from UGMs. It belongs to the MCMC approaches and is a special case of the Metropolis-Hastings algorithm (S. Geman and D. Geman, 1984). The Gibbs sampler relies on the full conditional probabilities and generates a Markov chain of samples which is proved to converge to the joint distribution. The Gibbs sampler is interesting when samples from the joint distribution cannot be directly drawn, as it is the case in general. The algorithm is given in Algorithm A.1.

Extensions of the Gibbs sampler algorithm are an active research topic, notable examples are the chromatic Gibbs sampler [START_REF] Gonzalez | Parallel gibbs sampling: From colored fields to thin junction trees[END_REF] or the blocked Gibbs sampler [START_REF] Brown | Sampling strategies for fast updating of Gaussian Markov random fields[END_REF]. Another extension for the Gibbs sampler consists in improving the exploration of the modes of the joint distribution using tempered transitions [START_REF] Ruslan | Learning in Markov random fields using tempered transitions[END_REF]. Lastly, in the particular case of Gaussian Markov Random Fields studied in Chapter 2, sampling can be done in a single iteration of a fully-blocked Gibbs sampler or a single sampling from a standard normal distribution followed by Fourier-based transformations (Rue and Held, 2005).

Remark: The notions of trees and chains also exist for UGMs, they will not be mentioned in this chapter for brevity.

Remark: In practice, many models, such as those studied in this thesis, have both directed and undirected edges. In such cases, the associated probability distribution, the parameter estimation process and the inference process for these mixed models will straightforwardly follow either from the DGM theory or UGM theory.

Remark: Note that Factor Graphs [START_REF] Frey | Extending factor graphs so as to unify directed and undirected graphical models[END_REF] are a popular kind of graphical models that aims at unifying UGMs and DGMs. Factor Graphs are however out of the scope of this thesis.

Inference in probabilistic models

Inference in a probabilistic model can refer to several closely related computational tasks:

• Computing a marginal distribution p(x x x x x x x A ), A ⊂ S.

• Computing a conditional distribution p(x x x x x x x A |x x x x x x x B ), A ⊂ S, B ⊂ S and A = B.

Inference in probabilistic models

• Computing a function of a probability distribution. Two important inference problems are:

-Marginal inference2 : ∀s ∈ S, xs = argmax xs p(x s ).

-Maximum A Posteriori (MAP) inference: x x x x x x x = argmax x x x x x x x p(x x x x x x x). In Section 1.5.2, these two inference problems will be studied in a Bayesian context.

Inference rapidly becomes a computationally intractable task for general graphs, hence the need for approximate computations in order to relieve the computational burden. Exact computations are however possible in DAGs and is still the subject of ongoing research for more general graphs. We use the term exact, or direct, for an inference process which benefits from exact computations, and approximate, or indirect, for an inference process which requires approximate computations. In the next sections, we first review key concepts for exact inference, then for approximate inference.

Exact inference

One of the basics of exact inference is the Variable Eliminitation algorithm which enables computing isolated marginals in arbitrary DGMs or UGMs (N. L. Zhang and Poole, 1994) [START_REF] Michael | Graphical models[END_REF]. However in most of the cases we are interested in several marginals and running many times the VE algorithm becomes computationally prohibitive and a waste of resources (most of the operations are the same throughout the instances of the algorithm).

In order to compute several marginals, the principle of the Variable Elimination algorithm is reformulated into a message-passing algorithm called Sum Product algorithm (or Belief Propagation algorithm with summations) [START_REF] Mark S Pearce | Radiation exposure from CT scans in childhood and subsequent risk of leukaemia and brain tumours: a retrospective cohort study[END_REF] [START_REF] Matthew | Structured belief propagation for NLP[END_REF]. The Sum Product algorithm is the basis for exact marginal inference in undirected acyclic graphs. In the case of MAP inference in undirected acyclic graphs, the Max Product algorithm (or Belief Propagation with maximizations) [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF] is the basic approach.

In the case of marginal inference applied to arborescences, the Sum Product algorithm becomes equivalent to the popular Forward Backward (FB) algorithm [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF], generalized in the Upward Downward (UD) algorithm (J.-B. Durand andGonçalves, 2001) (Laferté et al., 2000). The FB algorithm is the foundation to many extensions following its principles. The FB and the UD algorithms are frequently used for inference in particular arborescences: the Markov Chains and the Markov Trees that will be presented in Section 1.6 and studied in this thesis. The FB and UD algorithms are presented respectively in Algorithm A.2 and in Algorithm 3.1. In the case of MAP inference in arborescences, the Viterbi algorithm [START_REF] Viterbi | Error bounds for convolutional codes and an asymptotically optimum decoding algorithm[END_REF]) is a popular approach (a version of the Max Product algorithm dedicated to arborescences). The Viterbi algorithm has been generalized to Markov Trees in [START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF].

Approximate inference

Approximate inference is the general approach for all the cases but those mentioned in the previous section. The first approach we refer to is the direct extension of the BP algorithm to arbitrary graphs: the Loopy Belief Propagation algorithm [START_REF] Matthew | Structured belief propagation for NLP[END_REF]. Despite being a popular approach to approximate inference, the convergence of the Loopy Belief Propagation algorithm is not guaranteed: the message massing procedure can never terminate on arbitrary graph structures.

Variational Inference has been introduced more recently. It recasts approximate inference into a deterministic optimization problem (Michael I Jordan et al., 1999) (Lauritzen and[START_REF] Steffen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF]) (L. K. [START_REF] Lawrence K Saul | Mean field theory for sigmoid belief networks[END_REF]. In this context, a probability distribution p in which inference is complex is approximated by a variational distribution q in which inference is easy (often exact using the methods from Section 1.3.1). The optimization problem consists in finding q which minimizes the reverse Kullback-Leibler (KL) divergence between p and q. More on Variational Inference will be seen in Chapters 3 and 4. We can note that it is a very active research topic (C. Zhang et al., 2018).

Let us now mention popular approaches for approximate inference in UGMs defined by a Gibbs distribution. Well-known approaches for the computation of the mode of a density are the simulated annealing approaches [START_REF] Kirkpatrick | High resolution simulation of nonstationary Gaussian random fields[END_REF]. Serial Gibbs Simulated Annealing (SA) (S. Geman and D. Geman, 1984) is a probabilistic approach based on samplings from a Gibbs distribution (Equation 1.4) with a varying temperature to ease the exploration of the modes of the distribution. The algorithm is theoretically guaranteed to converge towards the mode of the distribution, but putting SA into practice can be cumbersome and one can end up with an approximation of the optimum we want to find. SA is described in Algorithm A.3. The literature on simulated annealing approaches is large. [START_REF] Michael | Simulated annealing: From basics to applications[END_REF] offers a good introduction on the topic. The Gibbs sampler presented in Section 1.2.3, as well as other sampling approaches, are then often part of an inference process. Additional comments on inference will be made in Section 1.5.2.

Remark: There are other popular approaches from different subfields to answer inference problems, in particular for MAP inference which seems more studied for image segmentation. For example, two other main approaches are based on linear programming [START_REF] Koller | Machine learning methods for histopathological image analysis[END_REF] or graph-cut [START_REF] Jörg | Multicuts and perturb & MAP for probabilistic graph clustering[END_REF]. However such approaches are out of scope of our work.

Parameter estimation

Either local conditional probabilities of DGMs or potential functions of UGMs include some parameters that need to be estimated. This section is devoted to the estimation of those parameters.

Remark: This section on parameter estimation does not cover the topic 1.4. Parameter estimation of learning the structure of the graphical models (Drton and Maathuis, 2017) which is a closely related topic. We also do not mention how parameters can be estimated in a fully-Bayesian statistics context (see Section 1.5.2).

Maximum Likelihood estimation

In a graph G, if the realizations of all the random variables are available, or observed, then we say that we possess the complete data. If the realizations of some random variables are not available, or hidden, then we say that we have incomplete data.

Let X X X X X X X the vector of the hidden variables of the model and

Y Y Y Y Y Y Y the vector of the observed variables. Z Z Z Z Z Z Z = (X X X X X X X, Y
Because of the logarithm that ties together the parameters, maximizing Equation 1.9 is much harder than maximizing Equation 1.8. The Expectation Maximization (EM) algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is an iterative approach to perform such a maximization in the case of hidden data. It alternates between an E-step in charge of inferring the missing values and a M-step which performs the actual maximization with the completed data. The EM algorithm is guaranteed to converge to a point of null gradient (local maximum or saddle point).

Note that the E-step is a step of probabilistic inference, and the methods seen in Section 1.3 need to be used. The EM algorithm is given in Algorithm A.4.

Parameter estimation for UGMs

In general, the parameter estimation techniques for DGMs cannot be applied to UGMs for mainly one reason: the partition function is intractable for UGMs and computations must always be approximate. The algorithms we review here are stochastic approximations of the previous algorithms or modified versions of the ML estimation.

Remark: In some exceptions the partition function is tractable (see for example [START_REF] Suetens | An introduction to conditional random fields[END_REF])), then one can apply the techniques of Section 1.4.2 to UGMs.

Probabilistic models for image segmentation

Stochastic Maximum Likelihood [START_REF] Younes | On the convergence of Markovian stochastic algorithms with rapidly decreasing ergodicity rates[END_REF], also discovered as Persistent Contrastive Divergence [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF].

Identity 1.11 can also be used to approximate the untractable summation which appears in the case of ML estimation with hidden variables (see Equation 1.9). This then extends the Stochastic Maximum Likelihood and Persistent Contrastive Divergence algorithms to the case of incomplete data.

Pseudolikelihood

One way to avoid the problem of the partition function is to change the objective function so that there is no more partition function to compute. This concept is known as the pseudolikelihood approach.

There exists many different pseudolikelihoods. A popular choice introduced in (J. [START_REF] Besag | Nearest-neighbour systems and the auto-logistic model for binary data[END_REF] is the pseudolikelihood where the joint distribution is a simple product of the factors: PL(θ θ θ θ θ θ θ; y y y y y y y) = log c∈C φ C (y y y y y y y C ).

(1.12)

The difficulty involved by the partition function is now discarded and one may continue the optimization with the approaches from Section 1.4.2.

Alternative forms for the EM algorithm

In the EM algorithm introduced in Section 1.4.2. Both the inference (E) and maximization (M) step can be complex for UGMs for all the reasons we have given so far. Moreover the other well-known limitations of the EM algorithm, such as its dependence to initial parameters and slow convergence, has motivated the development of stochastic versions of the EM algorithm [START_REF] Celeux | On Stochastic Versions of the EM Algorithm[END_REF]. Among them, the Stochastic Expectation Maximization (SEM) [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF] algorithm will be studied in this thesis. Its principle is to replace the E-step by sampled realizations of the hidden variables to complete the data. The algorithm is not proved to converge, however it empirically behaves well in practice. Algorithm A.5 describes the SEM algorithm.

When the inference step (E) is intractable, another common idea is to approximate this step with variational approximations of the EM algorithm [START_REF] Michael I Jordan | An introduction to variational methods for graphical models[END_REF].

Remark: Stochastic versions of the EM algorithm can also be applied to DGMs.

1.5. Probabilistic models for image segmentation 1.5.1. Definitions and context In this section we give the general setting in which probabilistic graphical models can be used for image segmentation.

Chapter 1. Main concepts in probabilistic modeling

For each pixel of the image, the associated graph will contain two vertices, both at site s. An observed variable Y s is associated with the first vertice and its realization is taken equal to the pixel value at site s. A hidden variable X s is associated with the second vertice, X s has value in the set of classes. The total number of sites is also the number of pixels and is taken equal to N . Then, more formally, X X X X X X X is a hidden process with value in

Ω Ω Ω Ω Ω Ω Ω N the set of classes of the image, X X X X X X X forms the segmented image. Y Y Y Y Y Y Y
is the observed process with value in R N or Z N depending on the pixel values.

The setting we have just described is the foundation for many works in image segmentation. In particular, Hidden Markov Models described in Section 1.6 are popular for image segmentation.

Bayesian image segmentation

In this section we focus our development on Bayesian inference applied to the context of image segmentation. We are interested in the posterior distribution of X X X X X X X given the observed process Y Y Y Y Y Y Y . Using Bayes theorem:

p(x x x x x x x|y y y y y y y) ∝ p(x x x x x x x)p(y y y y y y y|x x x x x x x), (1.13) where p(x x x x x x x) is called the prior distribution and p(y y y y y y y|x x x x x x x) is the conditional likelihood distribution.

Once the posterior distribution is computed, a decision needs to be taken to produce the actual segmentation. This is done in the context of Bayes decision theory. We first need to define a loss function L which measures the compatibility of an estimation with the hidden truth. In the case of image segmentation, L measures the compatibility between the estimated segmentation x x x x x x x and the hidden class image x x x x x x x. The optimal estimation, having observed Y Y Y Y Y Y Y = y y y y y y y, is defined to be that minimizing the posterior expected loss, L(x x x x x x x|y y y y y y y):

L(x x x x x x x|y y y y y y y) = E x x x x x x x∼p(x x x x x x x|y y y y y y y) [L(x x x x x x x, x x x x x x x)] = x x x x x x x p(x x x x x x x|y y y y y y y)L(x x x x x x x, x x x x x x x).

(1.14)

One common loss function, not restricted to image processing, is the zero-one loss function defined by:

L 1 (x x x x x x x, x x x x x x x) = 0, if x x x x x x x = x x x x x x x, 1, otherwise.
(1.15)

It can be shown that the estimation which minimizes L(x x x x x x x|y y y y y y y) when L = L 1 , is the Maximum A Posteriori (MAP) (S. Geman and D. Geman, 1984) Remark: The Bayesian framework of estimation (centered around the use of prior distributions) still holds in a more general setting than image segmentation. In particular, prior distributions can be used over the model parameters themselves in the case of fully Bayesian parameter estimation, which is not covered in this thesis. See [START_REF] Ghahramani | Estimation of generalized multisensor hidden Markov chains and unsupervised image segmentation[END_REF] for an introduction and [START_REF] Morris | Fully Bayesian image segmentation-an engineering perspective[END_REF] [START_REF] Dobigeon | Hierarchical Bayesian sparse image reconstruction with application to MRFM[END_REF] [START_REF] Vacar | Unsupervised joint deconvolution and segmentation method for textured images: a Bayesian approach and an advanced sampling algorithm[END_REF] for study examples.

Discriminative and generative models

We conclude this section by noting the difference between discriminative and generative modeling (Minka, 2005) (Sutton and[START_REF] Suetens | An introduction to conditional random fields[END_REF].

In the case of image segmentation, computing the posterior distribution p(x x x x x x x|y y y y y y y) is the end goal in order to effectively segment the image. There are two ways of doing so. On the one hand, one can directly model the posterior distribution p(x x x x x x x|y y y y y y y), this is called discriminative modeling. Chapter 4 of the thesis deals with such a modeling approach. On the other hand, one can model the joint distribution p(x x x x x x x, y y y y y y y) first, which is, by Bayes rule, the product of a prior p(x x x x x x x) and a conditional likelihood p(y y y y y y y|x x x x x x x). The posterior is then obtained with Equation 1.13. This is called generative modeling and Chapters 2 and 3 study generative models.

While it is hard to predict which type of modeling will perform best on a given dataset, it is commonly admitted that generative modeling is more suited to unlabeled datasets or small datasets where modeling a prior can be beneficial and have a regularizing effect. However, on larger annotated datasets, discriminative models tend to perform more accurately (Ng and [START_REF] Radford | On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes[END_REF].

Hidden Markov Models

Hidden Markov Models (HMMs) [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF] are the most popular type of probabilistic graphical models. HMMs are generative models which relate to the Bayesian formulation with a prior and a conditional likelihood presented in Section 1.5.2. The image processing context is generalized here. HMMs are a central element of this thesis, and this section is dedicated to these models and some of their extensions.

Chapter 1. Main concepts in probabilistic modeling

Main families of HMMs Hidden Markov Chains

A Hidden Markov Chain (HMC) [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF]) is a probabilistic model particularly suited for sequence modeling. The associated graph G is here a directed chain, called, in this context, Markov Chain (MC). Let X X X X X X X = (X 1 , X 2 , . . . , X N ) be a random vector or process. More precisely we have that X

X X X X X X is a MC if: p(x s |x s-1 , . . . , x 1 ) = p(x s |x s-1 ), ∀s ∈ {2, . . . , N }.
(1.19) Equation 1.19 means that the realization of the random variable at site s, given the preceding realizations in the chain, depends only on the realization of the random variable at site s -14 . Using the conditioning property, one can show that the joint distribution of the MC X X X X X X X is then: .20) Note that in our work we will consider MCs that are homogeneous: p(x s |x s-1 ) is the same ∀s ∈ {2, . . . , N }, see [START_REF] Brémaud | Discrete Probability Models and Methods: Probability on Graphs and Trees, Markov Chains and Random Fields, Entropy and Coding[END_REF] for more details. A MC is then defined by a root distribution and transition distribution.

p(x x x x x x x) = p(x 1 )p(x 2 |x 1 ) . . . p(x N |x N -1 ). ( 1
In the HMC model, X X X X X X X then is a MC whose variables are hidden in which we want to perform an inference task (p(x x x x x x x) is the prior). It is associated to an observed process Y Y Y Y Y Y Y whose variables are taken, in the most common case, independent given x x x x x x x (p(y y y y y y y|x x x x x x x) is the conditional likelihood) and depending only of the hidden realization at the same site. This is the hypothesis of the Independent Noise

(IN). The joint distribution (X X X X X X X, Y Y Y Y Y Y Y ) of such HMC-IN is then: p(x x x x x x x, y y y y y y y) = p(x 1 ) s∈{2,...,N } p(x s |x s-1 ) s∈{1,...,N } p(y s |x s ). (1.21)
In HMCs, inference is directly computed with the Forward Backward (FB) algorithm. In an unsupervised context, using the EM algorithm with FB at the E-step is known as the Baum Welch algorithm [START_REF] Jelinek | Design of a linguistic statistical decoder for the recognition of continuous speech[END_REF]. Figure 1.4 illustrates a HMC with Independent Noise.

The HMC has been used in a wide variety of contexts, famous ones are gene prediction [START_REF] Stanke | Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources[END_REF], speech processing [START_REF] Toda | Histopathologic characterization of peripheral arteries in subjects with abundant risk factors: correlating imaging with pathology[END_REF] or stock option 1.6. Hidden Markov Models forecasting (Gupta and Dhingra, 2012). In Chapter 5 of this thesis, we develop a model of HMC for image processing.

Hidden Markov Trees

We focus on two particular arborescence graphs which are known as Markov Trees (MTs) [START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF]:

• Dyadic MT: all the vertices of the arborescence have two sons (except on the last layer).

• Quadtree MT: all the vertices of the arborescence have four sons (except on the last layer).

Let S = {S 1 , . . . , S L } be the layer-wise subdivision of the set of vertices. S 1 is the first layer containing only one vertice, the root r. Dyadic or Quadtree MTs then have similar factorizations to MCs. Let X X X X X X X be a MT, then, ∀s ∈ S:

p(x s |x x x x x x x s ′ ∈S,s ′ =s ) = p(x s |x P(s) ).
(1.22) Equation 1.22 is the direct extension of Equation 1.19 to general arborescences.

Recall that in arborescences, the set of parent of a vertice is either empty or contains only one vertice denoted s -. The joint distribution of a MT X X X X X X X is then:

p(x x x x x x x) = p(x r ) s∈ S p(x s |x s -).
(1.23)

Hidden Markov Trees (HMTs) are similarly defined as HMCs by a root distribution and a transition distribution. In a HMT, the hidden process X X X X X X X has MT joint distribution. The conditional likelihood distribution implying the observed process Y Y Y Y Y Y Y can be freely formed, but in the most popular case, it is taken as independent given X X X X X X X realizations and depending only on realizations of X X X X X X X at the same site (case of Independent Noise as in Section 1.6.1). The joint distribution

(X X X X X X X, Y Y Y Y Y Y Y ) for such a HMT-IN is then: p(x x x x x x x, y y y y y y y) = p(x r ) s∈ S p(x s |x s -) s∈S p(y s |x s ).
(1.24) HMTs and derived models, known as Latent Tree Models, are popular in multi-resolution analysis [START_REF] Kinebuchi | Image interpolation using wavelet based hidden Markov trees[END_REF] [START_REF] Durand | Statistical inference for hidden Markov tree models and application to wavelet trees[END_REF] and in phylogenetic analysis [START_REF] Mourad | A survey on latent tree models and applications[END_REF] [START_REF] Zwiernik | Semialgebraic statistics and latent tree models[END_REF]. To a lesser extent, HMTs have also been used for image processing [START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF] [START_REF] Hanzouli | PECT/CT image denoising and segmentation based on a multi observation and multi scale Markov tree model[END_REF]. Models based on HMTs are presented in Chapter 3 of this thesis.

(J.-B.

Hidden Markov Fields

Hidden Markov Fields (HMFs) (S. Z. [START_REF] Li | Markov random field modeling in image analysis[END_REF] are UGMs but they are constructed as the other previously described HMMs. In HMFs, the hidden process

Chapter 1. Main concepts in probabilistic modeling Y s X s X r
Y Y Y Y Y Y .
X X X X X X X has the joint distribution of a Markov Field (Section 1.2.3). As an example, using the same conditional likelihood as before, i.e, an Independent Noise, we

get the classical joint distribution (X X X X X X X, Y Y Y Y Y Y Y
) of a HMF with Potts prior (J. Besag, 1986)5 and independent Gaussian noise:

p(x x x x x x x, y y y y y y y) = 1 Z exp     s∈S s ′ ∈ Ns∪{s} βδ x s ′ xs + s∈S log √ 2πσ xs - (y s -µ xs ) 2 2σ 2 xs     .
(1.25) In the last equation, δ is the Kronecker function and β, µ and σ are the model parameters (the parameters will be seen in depth in Chapter 2). The graphical model of a HMF with Independent Noise can be seen in Figure 1.6.

From the pioneering works of (S. Geman and D. Geman, 1984) and [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF], HMFs have been used a lot in image processing fields such 1.7. Conclusion as image segmentation [START_REF] Schmitt | A survey on image data augmentation for deep learning[END_REF] (Mignotte et al., 1999) (Kato and[START_REF] Jörg | Multicuts and perturb & MAP for probabilistic graph clustering[END_REF] or computer vision [START_REF] Wang | Markov random field modeling, inference & learning in computer vision & image understanding: A survey[END_REF].

Pairwise and Triplet extensions

Pairwise HMMs make the hypothesis that the joint distribution

(X X X X X X X, Y Y Y Y Y Y Y ) of the HMM is a Markov Chain/Markov Tree/Markov Field.
For example, considering chains, making the Pairwise Markov Chain assumption leads us to the general joint distribution of a Pairwise Markov Chain:

p(x x x x x x x, y y y y y y y) = p(x 1 , y 1 ) s∈{2,...,N } p(x s , y s |x s-1 , y s-1 ) (1.26)
It appears that this equation is a strict generalization of the classic HMC (Equation 1.21). This means that more correlations between random variables can be modeled. In particular, X X X X X X X is not restricted to be a Markov process anymore. Similar observations can be made in the case of Pairwise Markov Trees and Pairwise Markov Fields.

Triplet HMMs go one step further by introducing a third auxiliary random process

V V V V V V V with values in Λ N .
We then make the assumption that the triplet

(X X X X X X X, V V V V V V V , Y Y Y Y Y Y Y ) is a Markov Chain/Markov Tree/Markov Field.
Again, in the case of the chains, a Triplet Markov Chain has for joint distribution:

p(x x x x x x x, v v v v v v v, y y y y y y y) = p(x 1 , v 1 , y 1 ) s∈{2,...,N } p(x s , v s , y s |x s-1 , v s-1 , y s-1 ) (1.27)
Again, the greater generality of TMCs over HMCs (Equation 1.21) but also over PMCs (Equation 1.27) is to be noted. In particular, neither X

X X X X X X, V V V V V V V , (Y Y Y Y Y Y Y , V V V V V V V ), (X X X X X X X, Y Y Y Y Y Y Y ) nor (X X X X X X X, V V V V V V V
) is necessarily a Markov process anymore. New pairwise models are developed in Chapter 2 of the thesis, while developments around triplet models can be found in Chapter 3. In the literature, applications of pairwise and triplet HMMs are much less numerous, but they exhibit each time better performances than the classical HMMs [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] [START_REF] Courbot | Pairwise Markov fields for segmentation in astronomical hyperspectral images[END_REF] [START_REF] Hanzouli-Ben | A framework based on hidden Markov trees for multimodal PET/CT image co-segmentation[END_REF] [START_REF] Gorynin | Pairwise Markov models for stock index forecasting[END_REF]) (Courbot, Monfrini, et al., 2018) [START_REF] Lanchantin | Unsupervised segmentation of randomly switching data hidden with non-Gaussian correlated noise[END_REF] [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF].

Remark: Auxiliary random variables are central to deep probabilistic models. Triplet Markov Trees and Spatial Bayes Networks developed in Chapter 3 can be seen as belonging to this family of models. Other deep probabilistic models will not be studied in this thesis and we refer the reader to [START_REF] Goodfellow | Deep learning[END_REF] for an in-depth overview.

Conclusion

Probabilistic graphical models are powerful approaches to express direct dependencies within a set of random variables capable of modeling a complex Chapter 1. Main concepts in probabilistic modeling phenomenon. Using graphical models is a way to control the expressiveness of the direct dependencies between random variables and the tractability of the models thanks to the conditional independencies. A successful model is based on two important tasks: the estimation of the parameters of the probability distribution and the inference task. Such tasks have been introduced in this first chapter. As we have seen, it is little to say that there is a wide variety of approaches to solve these challenges and put into practice probabilistic models.

Our presentation also focused on probabilistic graphical models which are very propular in the signal processing field: the Hidden Markov Models. Indeed, the latter offer a good compromise between the dependencies they introduce and their computational performances. HMMs are a family of models which is central in this thesis.

Many of our contributions presented in the next chapters are new probabilistic models that enrich the dependencies between the random variables. Our goal is to model more complex phenomena to address problems from real world applications.

Introduction

Probabilistic models for spatially correlated random variables

This chapter studies the task of unsupervised image segmentation, in the context of Bayesian image segmentation presented in Section 1.5.2, when strong spatially correlated noise corrupts the image. In this case, classical approaches reach their limits and new dedicated models need to be considered to improve the accuracy of the segmentation. We will focus on Pairwise Markov Fields (PMFs) presented in Section 1.6.2. Such models relax the Markovian assumption on the hidden process, which enables the modeling of more complex correlations, while keeping Bayesian inference easily available. Figure 2.1 depicts the undirected graphical models of the classical HMF model and of several PMF models.

We introduce Gaussian Random Fields (GRFs) and Gaussian Markov Random Fields (GMRFs) (Rue and Held, 2005) which are powerful probabilistic models capable of dealing with a large variety of correlated random variables and especially with long-range spatial correlations between pixels in images. GRFs and GMRFs are Undirected Graphical Models (UGMs) defined with re-

spect to a graph G. Let X X X X X X X = (X 1 , . . . , X n ) T , n < ∞, X X X X X X X is a GRF with mean µ µ µ µ µ µ µ (a (n × 1) vector) and Symmetric Positive-Definite (SPD) covariance matrix Σ (a (n × n) matrix), if,
and only if, its density is a multivariate Gaussian function:

p(x x x x x x x) = (2π) -n 2 det (Σ) -1
While the graphical representation of a GRF is a fully connected graph, a conditional independence assumption (Markov property) is made to define a GMRF. The latter is defined as a GRF with respect to G but in which a vertice s is only connected to the vertices of a subset N s ⊂ S called the neighborhood of s. The associated random variable, X s , is then independent of X s ′ , ∀s ′ ∈ S \N s , given the realizations of

X s ′′ , ∀s ′′ ∈ N s . In such case, Q s ′ ,s ′′ = 0 with Q = Σ -1 ,
where Q is called the precision matrix.

Then formally, the density of a GMRF with mean µ µ µ µ µ µ µ and precision matrix Q is given by:

p(x x x x x x x) = (2π) -n 2 det Q -1 -1 2 exp - 1 2 (x x x x x x x -µ µ µ µ µ µ µ)Q(x x x x x x x -µ µ µ µ µ µ µ) . (2.2)
The conditional independence assumption is a way to alleviate the computational cost of the probabilistic model. Indeed, smaller neighborhoods (few edges in G) lead to a sparse Q which makes computations based on Equation 2.2 much easier. Appendix B recalls the main definitions and properties of GRFs and GMRFs.

Researches have been carried for regression and classification using conditional GMRF, i.e., these studies consider a discriminative probabilistic model [START_REF] Radosavljevic | Continuous Conditional Random Fields for Regression in Remote Sensing[END_REF] [START_REF] Vemulapalli | Gaussian conditional random field network for semantic segmentation[END_REF] [START_REF] Petrović | Gaussian Conditional Random Fields for Classification[END_REF].

Introduction

However, to the best of our knowledge, there exists no literature on generative probabilistic models using GMRF for image segmentation. This chapter aims at proposing such a new probabilistic modelization. Indeed, generative models might be advantageous when no training dataset is available, for example, in the case of unsupervised image segmentation.

We introduce a new model, called Gaussian Pairwise Markov Fields (GPMF), which belongs to the PMF family. It combines a generalization of the classical HMF model and the ability to model strongly correlated variables as a GMRF, while preserving tractability. First, we show that the PMF hypothesis is a natural way to answer the problem of modeling correlated noise and introducing long correlations by using the GMRF model as conditional likelihood. Secondly, the new model introduces depencies between the observed and latent processes which generalize the usual Markovianity hypothesis of the hidden process. Besides, we handle the unsupervised image segmentation problem which requires a crucial step of parameter estimation. To this end, we propose a stochastic parameter estimation algorithm for the PMF models. In the case of unsupervised image segmentation, the GPMF model performs better than other classical unsupervised approaches.

Remark: Pairwise Markov Field models have been defined independently in [START_REF] Dimitrova | Graphical Models Over Heterogeneous Domains and for Multilevel Networks[END_REF]) and (Y. [START_REF] Park | Learning the network structure of heterogeneous data via pairwise exponential Markov random fields[END_REF]. However, those works mostly focus on structure learning without latent variables (Drton and Maathuis, 2017), which significantly modifies the work hypothesis. To the best of our knowledge, there is no work on structure learning in Pairwise Markov Fields with latent variables. Thus, despite similar model names, the content of our chapter deals with a different problem since latent variables are central to our work.

Link with image processing in vascular surgery

In the context of mCT X-rays scans of human arteries that contain a metallic stent. These images are very noisy because of the strong artifacts caused by the interactions between the X-rays and the metallic stent. To improve the design of biomaterials and their implantation in the human body, it is crucial to analyze, in situ, the biomaterial when it fails and when it needs to be explanted from the patient body. In such analyses, we need to precisely segment, the stent, the organic material and the background despite the strong noise in the image. The scarcity of such data requires the use of an unsupervised approach as developed in this chapter. In this image processing problem, the stent artifacts are modeled as a correlated noise. An automatization of the segmentation process could help processing more data and create enhanced inputs for biomechanical studies (such as 3D meshes) carried to increase the knowledge about the vascular diseases.

In Section 2.5.4, we present the application of the new probabilitic model dedicated to handle correlated noise in the task of unsupervised segmentation of mCT images in presence of metallic artifacts.

Chapter 2. Gaussian Pairwise Markov Fields X X X X X X X Y Y Y Y Y Y Y (a) HMF-IN X X X X X X X Y Y Y Y Y Y Y (b) PMF X X X X X X X Y Y Y Y Y Y Y (c) PMF-UN X X X X X X X Y Y Y Y Y Y Y (d) HMF-CN
(d). Within the PMF models, (X X X X X X X, Y Y Y Y Y Y Y )
is the Markov process. GPMFs belongs to the PMF family, as well as PMFs-UN and HMFs-CN which are intermediate cases between the HMF and PMF (with full direct dependencies) models in terms of direct dependencies. PMFs-UN and HMFs-CN will be introduced and studied in Section 2.3. We note the numerous correlations that can be introduced between random variables when the pairwise assumption holds.

Gaussian Pairwise Markov Fields

2.2.1. Model definition X X X X X X X = (X 1 , . . . , X N ) is a discrete-valued random vector with values in Ω N , with Ω = (ω 0 , . . . , ω K-1 ). Y Y Y Y Y Y Y = (Y 1 , . . . , Y N ) is a real-valued random vector with values in R N . (X X X X X X X, Y Y Y Y Y Y Y
) is a stationary Markov process on a graph, whose vertices are indexed by S such that |S| = N . The Markov process is defined with respect to the neighborhood N and we have N = (N X ∪ N Y ) ⊂ S. Indeed, in general, neighborhoods are different if we consider the hidden or observed variables, thus, a neighborhood can be decomposed on a part linked to the hidden variables (N X ) and another to the observed variables (N Y ).

We use the notation p to refer to factors (or unnormalized probability distributions, see Section 1.2.3). The classical Hidden Markov Field model with Independent Noise (HMF-IN) is defined by the joint distribution of Equation 1.25 that we rewrite as:

p(x x x x x x x, y y y y y y y) ∝ s∈S p(x s |x x x x x x x N X s )p(y s |x s ). (2.
3)

The undirected graphical model of a HMF-IN is given in Figure 2.1a.

The PMF family of models generalizes the HMF-IN model as mentioned in Section 1.6.2. Recall that in the PMF family, the assumption of (

X X X X X X X, Y Y Y Y Y Y Y ) being a Markov field is made. Note that this implies that X X X X X X X given some realizations 2.2. Gaussian Pairwise Markov Fields y y y y y y y of Y Y Y Y Y Y Y is a Markov field.
A PMF is defined by the distribution:

p(x x x x x x x, y y y y y y y) ∝ s∈S p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ), (2.4)
The greater generality of the pairwise fields hypothesis enables the modeling of more complex correlations between variables. Indeed, neither X

X X X X X X nor Y Y Y Y Y Y Y is necessarily Markovian.
Following the notations of Section 1.5, the X X X X X X X are the hidden variables, and the Y Y Y Y Y Y Y are the observed variables. The new GPMF model then offers new modelization properties, notably, it is capable of handling Gaussian spatially correlated noise. Moreover, in the context of Bayesian image segmentation, an estimation of the hidden truth will be done under the MAP criterion (Equation 1.16) and the MPM criterion (Equation 1.18). The MAP estimator is here classicaly estimated with the Simulated Annealing (Algorithm A.3). The MPM will be approximated by the Marroquin algorithm (Algorithm A.6) [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF] using the local expression of the posterior Markov field:

p(x s |x x x x x x x N X s , y y y y y y y) = p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) x ′ s p(x ′ s , y s |x x x x x x x N X s , y y y y y y y N Y s )
.

(2.5)

Description of the GPMF distribution

In our new Gaussian Pairwise Markov Field (GPMF) model we define the factor of Equation 2.4 by, ∀s ∈ S:

p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) = exp -2 s ′ ∈N X s V (x s , y s ) + 1 2 Q s,s ȳ2 s + s ′ ∈N Y s Q s,s ′ ȳs ȳs ′ , (2.6)
where we denote ȳs = y sµ xs and the potential function V is given by:

V (x s , y s ) = -δ x s ′ xs β 1 - 1 2 (ȳ s -ȳs ′ ) 2 , (2.7) with δ x s ′
xs the Kronecker delta function. Equation 2.6 is parametrized by a granularity parameter β, spatially varying means µ xs , and a precision matrix Q. These model parameters are detailed in Section 2.2.4.

Remark:

The potential function of Equation 2.7 can be seen as an enhanced Potts potential (S. Z. [START_REF] Li | Markov random field modeling in image analysis[END_REF]) (Pereyra et al., 2013) where the granularity coefficient β is adjusted to the observations available in the neighborhood of each site.

The GPMF conditional likelihood

The local expression of the Pairwise Markov Field of Equation 2.4 derives from the following joint Gibbs distribution (see Appendix C.1 for a proof):

p(x x x x x x x, y y y y y y y) = 1 Z exp -E(x x x x x x x, y y y y y y y) , (2.8)
where E(x x x x x x x, y y y y y y y) is the joint energy and Z is the normalization constant. It can be shown that:

E(x x x x x x x, y y y y y y y) = s∈S s ′ ∈ N X s ∪{s} -δ x s ′ xs β 1 - 1 2 (ȳ s -ȳs ′ ) 2 + s∈S s ′ ∈ N Y s ∪{s} 1 2 Q s,s ′ ȳs ȳs ′ -β 1 - 1 2 (ȳ s -ȳs ) 2 , = s∈S s ′ ∈ N X s ∪{s} -δ x s ′ xs β + s∈S s ′ ∈ N X s ∪{s} 1 2 δ x s ′ xs β(ȳ s -ȳs ′ ) 2 + s∈S s ′ ∈ N Y s ∪{s} 1 2 Q s,s ′ ȳs ȳs ′ -β.
(2.9)

Using standard calculus rules, one can show that:

s∈S s ′ ∈ N X s ∪{s} 1 2 δ x s ′ xs β(ȳ s -ȳs ′ ) 2 =
1 2 ȳ y y y y y y T P ȳ y y y y y y, (2.10)

where P is a matrix with elements:

P s,s ′ =      2 s ′ ∈N X s δ x s ′ xs β, if s = s ′ , -2δ x s ′ xs β, if s ′ ∈ N X s , 0 otherwise. (2.11)
We have l = |N X | + 2. Note that, P approximates a diagonally banded matrix for two reasons. First, beyond a certain range (the N X neighborhood), all the terms in P are null, starting from the diagonal. Second, among the remaining terms, some are null because x s = x s ′ . Hence, with the restriction β ∈ R * + , P has the property of diagonal dominance which makes P a SPD matrix. If Q is also a SPD matrix, R = P + Q is a SPD matrix1 . The energy can then be

Gaussian Pairwise Markov Fields

written:

E(x x x x x x x, y y y y y y y) = (s,s ′ )∈S 2 1 2 P s,s ′ ȳs ȳs ′ + (s,s ′ )∈S 2 1 2 Q s,s ′ ȳs ȳs ′ - s∈S s ′ ∈ N X s ∪{s} δ x s ′ xs β, = (s,s ′ )∈S 2 1 2 R s,s ′ ȳ ȳs ′ - s∈S s ′ ∈ N X s ∪{s} δ x s ′ xs β.
(2.12)

The last term of the energy gets canceled between the denominator and the numerator in Equation ( 3). Using the result of the integral of the multivariate Gaussian, we finally have:

p(y y y y y y y|x x x x x x x) = 1 (2π) N det(R -1 )
exp -1 2 ȳ y y y y y y T Rȳ y y y y y y , (2.13) which is the density of a Gaussian Markov Random Field (GMRF) with nonstationary mean µ µ µ µ µ µ µ x x x x x x x and precision matrix R = P + Q (Cressie and Verzelen, 2008) (Rue and Held, 2005). R corresponds to the precision matrix Q which is perturbed by P .

While we need to restrict β ∈ R * + , typical values of β are much smaller than 1. Hence, in practice, the non-zero entries in P are much smaller than the entries of Q at the same site. Then, one assumes R = P + Q ≈ Q. Since R is the precision matrix of a GMRF, it is sparse and so is Q. This approximation is the foundation for the parameter estimation procedure we propose, see Section 2.4.

The model parameters

We now detail the model parameters introduced in the previous section. The non-stationary mean vector µ µ µ µ µ µ µ, associated with the GMRF, is dependent on x x x x x x x, the realizations of the field X X X X X X X, such that ∀s ∈ S,

µ s = µ xs ∈ R. β ∈ R *
+ is a coefficient whose role is similar to the granularity parameter in a classic Potts model (Pereyra et al., 2013).

Q is a SPD matrix which is approximated to be the precision matrix of the GMRF defined by the conditional likelihood, i.e., R ≈ Q, as explained in Section 2.2.3. This approximation leads to a much simpler parameter estimation procedure.

In order to use computationally efficient spectral methods for matrix manipulations, we also make a periodic boundary assumption. A detailed presentation of the periodic boundary assumption and its consequences is given in Appendix B.2. Then Q is a block-circulant matrix with circulant blocks whose inverse is the covariance matrix Σ = Q -1 . Σ is defined as a stationary covariance matrix with variance σ ∈ R * + , and associated to an exponential correlation function with decay r ∈ R * + defined in Equation B.

(X X X X X X X, Y Y Y Y Y Y Y ), which leads to the Markovianity of Y Y Y Y Y Y Y given a realization of X X X X X X X
, implies that Q is a sparse matrix. This leads to efficient computations of the GPMF model equations.

Remark: The variance of the GMRF is stationary, as opposed to the mean of the GMRF. The simulation of the GMRF must be carried through its conditional equations (Rue and Held, 2005) [START_REF] Brown | Sampling strategies for fast updating of Gaussian Markov random fields[END_REF] which are updated at each iteration since the neighbor values change. In this context it is known that introducing a non-stationary variance is very complex [START_REF] Fuglstad | Exploring a new class of non-stationary spatial Gaussian random fields with varying local anisotropy[END_REF] and is out of scope of this thesis.

Related models: the PMF model family

In this section we continue to consider

(X X X X X X X, Y Y Y Y Y Y Y )
as a Markov process and we present several models related to GPMF, and show that GPMF generalizes them. Using the conditioning formula we know that:

p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) =p(x s |x x x x x x x N X s , y y y y y y y N Y s )p(y s |x x x x x x x N X s , y y y y y y y N Y s , x s ). (2.14)
Thus, the factor of Equation 2.4 can be split according to Equation 2.14 and, after normalization, a PMF has a distribution that can be written as:

p(x x x x x x x, y y y y y y y) ∝ s∈S p(x s |x x x x x x x N X s , y y y y y y y N Y s )p(y s |x x x x x x x N X s , y y y y y y y N Y s , x s ). (2.15)
The last equation is a way to look at PMFs which highlights their greater generality and the more complex direct dependencies between variables than in HMFs (compare with Equation 2.

3).

Then we propose new models based on Equation 2.14, by new factors with different conditional independence assumptions. These new intermediate models will be less general than the GPMF model but could also perform well at a smaller computational cost. We first present some new factors and then we use them to define models within the PMF model family with Equation 2.14.

Let us recall the Potts prior for Markov fields (Kato and Zerubia, 2012), it is based on the following local conditional probabilities:

p(x s |x x x x x x x N X s ) = exp   -2 s ′ ∈N X s V P (x s , x s ′ )   , (2.16)
where V P is the Potts potential function (J. [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] involving of a granularity coefficient β:

V P (x s , x s ′ ) = -δ x s ′ xs β.
(2.17)

Let us also recall the independent Gaussian conditional likelihood (Kato and Zerubia, 2012) which leads to the local conditional probabilities:

p(y s |x s ) = exp -ln( √ 2πσ) - ȳ2 s 2σ 2 .
(2.18)

Related models: the PMF model family

We also have the strict multivariate generalization of Equation 2.18 to the conditional GMRF likelihood in its local form (Rue and Held, 2005) [START_REF] Brown | Sampling strategies for fast updating of Gaussian Markov random fields[END_REF]:

p(y s |x s , x x x x x x x N X s , y y y y y y y N Y s ) = exp -ln 2πV s - (y s -M s ) 2 2V s , (2.19)
with

M s = µ xs -V s s ′ ∈S s ′ =s Q s,s ′ ȳs ′ and V s = Q -1 s,s . (2.20)
Let us introduce a generalized version of the Potts potential taking into account the spatial context:

p(x s |y s , x x x x x x x N X s , y y y y y y y N Y s ) = exp -2 s ′ ∈N X s -δ x s ′ xs β 1 - 1 2 (ȳ s -ȳs ′ ) 2 .
(2.21)

Given the local conditional probabilities stated so far, three models belonging to the PMF family can be designed: [START_REF] Courbot | Pairwise Markov fields for segmentation in astronomical hyperspectral images[END_REF].

• Potts-Gaussian Markov Field (P-GMRF)2 , with local distributions given by Equations 2.16 and 2.19, which belongs to the HMF-CN family (depicted in Figure 2.1d). We are not aware of similar models in the literature.

Table 2.1 summarizes the local distributions of these models. Note that GPMF, PMF-UN and P-GMRF all introduce more direct dependencies than the P-IN model. PMF-UN enhances the local distributions for the hidden variables and P-GMRF for the observed variables. Besides, note that all these models are submodels of the GPMF in the sense that they have direct dependencies that are ignored with respect to the GPMF direct dependencies. The submodels all ignore some direct dependencies that the pairwise factorization of Equation 2.4 integrates.

E ′ (x x x x x x x, y y y y y y y) = s∈S s ′ ∈ N X s -δ x s ′ xs β + s∈S s ′ ∈N Y s 1 2 Q s,

Parameter estimation

Stochastic Parameter Estimation

In this section, we develop an algorithm for the unsupervised parameter estimation task in the GPMF model. Without loss of generality, in the following development, we consider K = 2, i.e. Ω = {ω 0 , ω 1 }. Therefore the model is described with 5 parameters. Let θ θ θ θ θ θ θ ∈ Θ be the vector of parameters, then:

θ θ θ θ θ θ θ = {µ ω0 , µ ω1 , β, σ, r} {µ 0 , µ 1 , β, σ, r}, with Θ = R 2 × (R * + ) 3 . (2.22)
We develop a variation of the Stochastic Expectation Maximization algorithm [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF], which we call Stochastic Parameter Estimation (SPE). We first give the statistical estimators of the parameters given the complete data (x x x x x x x, y y y y y y y).

Generalization of the Linear Least Square estimator for GPMF

A generalization3 of the approach of (Derin and H. Elliott, 1987) is established to retrieve the parameter β in the PMF family, using the Linear Least Square (LLS) estimator (J. [START_REF] Friedman | The elements of statistical learning[END_REF] and a completed pair of realizations of (x x x x x x x, y y y y y y y). The derivation is done for the GPMF model and is similar for the other models.

First note that, ∀s ∈ S: p(x s , x x x x x x x N X s |y y y y y y y) p(x x x x x x x N X s |y y y y y y y) = p(x s |x x x x x x x N X s , y y y y y y y),

= p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) x ′ s ∈Ω p(x ′ s , y s |x x x x x x x N X s , y y y y y y y N Y s ) , = p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) x ′ s ∈Ω p(x ′ s , y s |x x x x x x x N X s , y y y y y y y N Y s )
.

(2.23)

Here the second equality has been seen in Equation 2.5. Then we have, ∀s ∈ S:

p(x s , y s |x x x x x x x N X s , y y y y y y y N Y s ) p(x s , x x x x x x x N X s |y y y y y y y) = x ′ s p(x ′ s , y s |x x x x x x x N X s , y y y y y y y N Y s ) p(x x x x x x x N X s |y y y y y y y) , (2.24)
where we can make the same key observation as in (Derin and H. Elliott, 1987): the right-hand side of the last equation is independent of the realization x s ∈ Ω.

Then so is the left-hand side. Then, ∀s ∈ S, ∀(x s , x ′ s ) ∈ Ω (2.25)

Now, taking the exponential on each side, and using the expression of Equation 2.6, ∀s ∈ S, ∀(x s , x ′ s ) ∈ Ω 2 :

ln p(x s , x x x x x x x N X s |y y y y y y y) p(x ′ s , x x x x x x x N X s |y y y y y y y) + 1 2 Q s,s ȳ2 s + t∈N Y s Q s,t ȳs ȳt - 1 2 Q s,s (ȳ ′ s ) 2 - t∈N Y s Q t,s ȳ′ s ȳt = 2 t∈N X s V (x s , y s ) -2 t∈N X s V (x ′ s , y s ), (2.26)
where we define ȳ′ s = (y sµ x ′ s ) which gives:

ln p(x s , x x x x x x x N X s |y y y y y y y) p(x ′ s , x x x x x x x N X s |y y y y y y y) + 1 2 Q s,s ȳ2 s + t∈N Y s Q s,t ȳs ȳt - 1 2 Q s,s (ȳ ′ s ) 2 - t∈N Y s Q t,s ȳ′ s ȳt = β 2 t∈N X s δ xt xs 1 - 1 2 (ȳ s -ȳt ) 2 - 2 t∈N X s δ xt x ′ s 1 - 1 2 (ȳ ′ s -ȳt ) 2 .
(2.27)

The last equation can be written for each site s, ∀(x s , x ′ s ) ∈ Ω 2 . All these equations can be put in the form:

a s = ln p(x s , x x x x x x x N X s |y y y y y y y) p(x ′ s , x x x x x x x N X s |y y y y y y y) + 1 2 Q s,s ȳ2 s + t∈N Y s Q s,t ȳs ȳt - 1 2 Q s,s (ȳ ′ s ) 2 - t∈N Y s Q t,s ȳ′ s ȳt , (2.29) Chapter 2. Gaussian Pairwise Markov Fields and b s = t∈N X s 2δ xt xs -2δ xt x ′ s + (µ x ′ s -µ xs )(-ȳ ′ s -ȳs + 2ȳ t ) .
(2.30)

with ȳ′ s = (y sµ x ′ s ). The probabilities p(x s , x x x x x x x N X s |y y y y y y y), ∀x s ∈ Ω, denote the probability of encountering a certain configuration of hidden variables at site s and its neighboring sites. These probabilities are independent of s and are estimated using the frequency estimator. The LLS estimator then gives that: .

(2.33)

Let us now write an estimator of the range of the exponential correlation function, r, given the completed data. As given in (Cressie, 1992)[Section 2.4], the correlogram of the field can be estimated by:

Ĉ(d) = 1 |D(d)| (s,s ′ )∈N (d) 1 σ2 (y s -μxs )(y s ′ -μx s ′ ), (2.34)
where D(d) is the set of pixel pairs whose Manhattan distance on the torus is d ∈ N + , i.e.:

D(d) = {(s, s ′ ) ∈ S 2 : min(|s 1 -s ′ 1 |, L 1 -|s 1 -s ′ 1 |)+ min(|s 2 -s ′ 2 |, L 2 -|s 2 -s ′ 2 |) = d}. (2.35)
where L 1 and L 2 are the lengths of the sides of the image we treat. The exponential correlation function has the form u = e -d r , where u is the correlation value, d the distance (between 0 and N ) and r the range to estimate. There is no constant in front of the exponential since the GMRF has been standardized, so we will estimate r ∈ R + * by fitting the exponential correlation function to the data points Ĉ(d) [START_REF] Ottar | Nonparametric spatial covariance functions: estimation and testing[END_REF]. Using the LLS estimator again we get:

Experiments and Results

In the context of latent variables, we only possess the observations and not the complete data. Thus the SPE algorithm successively repeats two steps. The first consists in simulating an estimation of the hidden layer, x x x x x x x, which is then used to form the completed data (x x x x x x x, y y y y y y y), to approximate the complete data 4 . The second step uses the parameter estimators defined previously over the completed data. Algorithm 2.1 details the steps of the procedure.

SPE is a variant of the Stochastic Expectation Maximization (SEM) [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF]. SPE differs from SEM because β and r are not estimated using a ML estimator. In SPE, at each iteration we sample from the posterior law to get the completed data, using Gibbs sampling, run during K(t) iterations. We set K(t) = t as proposed in (Carreira-Perpiñan and [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF]. In practice we stopped the SPE algorithm after 30 iterations, such a number was enough for the estimators to converge to a value. Data: θ θ θ θ θ θ θ 0 = {β 0 , µ 0 0 , µ 0 1 , σ 0 , r 0 }, the initial set of parameters, y y y y y y y, the observations, x x x x x x x 0 , the initial configuration for Gibbs sampler. Result: θ θ θ θ θ θ θ = { β, μ0 , μ1 , σ, r}, the estimated parameters. t ← 1 while convergence is not attained do /* Posterior sampling with a Gibbs sampler initialized at x x x x x x x t-1 and run during K(t) steps */ x x x x x x x t ∼ p(x x x x x x x|y y y y y y y; θ θ θ θ θ θ θ t-1 ) /* Estimation */ LLS estimator for β t (Eq.2.31). ML estimator for µ t 0 and µ t 1 (Eq.2.32). ML estimator for σ t (Eq.2.33). Estimation via correlogram for r t (Eq.2.36).

θ θ θ θ θ θ θ t ← {β t , µ t 0 , µ t 1 , σ t , r t } t ← t + 1 end Remark:
The same estimation procedure can be used to estimate the parameters of the other related models described in Section 2.3.

Experiments and Results

This section illustrates the models from the PMF family in the practical task of image segmentation in situations of various complexity. These models are also compared to other methods from the literature of unsupervised image segmentation.

In this section we consider real-world images but also semi-real images. The latter are real binary images artificially corrupted by correlated noise, in order to test the relevance and tractability of the new model and its counterparts. The images used come from the "1070-Binary Shape Database"5 . We want to evaluate the capacity of our model to handle correlated noise, the main purpose of this work, but also to test the parameter estimation procedures with the proposed SPE algorithm (Algorithm 2.1). Thus, from the natural binary image we construct an observation with an additive correlated Gaussian noise over the image. Such a noise is modeled by a GRF. Beforehand, we discuss the efficiency of sampling from the probabilistic models.

Improved sampling with Tempered-Gibbs sampler

In this section we study more in depth the problematic of sampling from the new distributions. Indeed, it is known that the classic Gibbs sampler procedure (or other MCMC-based sampling approaches) suffers from a poor exploration of the probability distribution. The algorithm is moreover dependent on the initialization of the Gibbs sampler. Running several Gibbs sampler with different random initializations might not be satisfactory. These issues are discussed in [START_REF] Radford | On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes[END_REF] which also introduces parallel tempering. The idea of parallel tempering is to run in parallel several Gibbs samplers at different temperatures (in the same meaning as in Simulated Annealing (S. Geman and D. Geman, 1984) presented in Algorithm A.3). Samples from the Gibbs samplers at high temperature can swap and become the current state of a Gibbs sampler of a lower temperature. At high temperature the distribution is less severely peaked and the Gibbs sampler might explore more easily this distribution. This idea at the origin of parallel tempering is represented in Figure 2.2. This approach is called the Tempered Gibbs sampler (T-Gibbs sampler).

In our case we use the parallel tempering approach to improve the sampling procedures from the P-GMRF and GPMF models. We develop a similar methodology as [START_REF] Cho | Parallel tempering is efficient for learning restricted Boltzmann machines[END_REF] where the temperature factor was used to rescale some precise terms of the energy function rather than the whole energy. In the algorithm that we propose, the energy is tempered (multiplied by a positive scalar inferior to one) so that, at high temperature the probability distribution of the complex model tends towards the distribution of the P-IN model. We follow this idea since sampling is easier in the P-IN model. This uses the fact that the PMF models are more general than the P-IN model. To do so we modify the potentials and make them dependent on a temperature parameter that we explicitly add in the notation:

Experiments and Results

T1

< T2 < T3 = ∞ T 1 T 2 T 3 Figure 2.2.:
The parallel tempering approach: as the temperature increases, the probability distribution (of the family Equation 1.4) is flattened. It ultimately tends towards an uniform probability distribution as the temperature tends towards infinity. In a less multimodal distribution, the Gibbs sampler might more easily explore the distribution and produce interesting samples.

• Equations 2.6 and 2.7 become:

pT (x s , y s |x x x x x x x N X s , y y y y y y y N X s ) = exp -2 s ′ ∈N X s V T (x s , y s ) + 1 2 Q s,s ȳ2 s + 1 T s ′ ∈N Y s Q s,s ′ ȳs ȳs ′ , (2.37)
where

V T (x s , y s ) = -δ x s ′ xs β 1 - 1 2T (ȳ s -ȳs ′ ) 2 .
(2.38)

• Equation 2.19 becomes:

pT (y s |x s , x x x x x x x N X s , y y y y y y y N Y s ) = exp -ln 2πV s - (y s -M T,s ) 2 2V s , (2.39) with M T,s = µ xs - V s T s ′ ∈S s ′ =s Q s,s ′ ȳs ′ and V s = Q -1 s,s . (2.40)
The tempered versions of the potentials we have just defined are used to sample from the GPMF and P-GMRF model, in the T-Gibbs sampler, according to Table 2

.2.
Finally, Algorithm 2.2 details the T-Gibbs sampler. In this algorithm, let K be the number of Gibbs samplers that we run in parallel. We then have K

Model

Equations used in T-Gibbs P-GMRF Eq. 2.16 and Eq. 2.39 GPMF Eq.2.37 Table 2.2.: Equations used to sample from the models with the T-Gibbs.

associated and ordered temperatures starting from the temperature at which the probabilistic model is defined; which is 1 in our case. Two chains with successive temperatures can swap with probability:

α k (x x x x x x x i+1 k , x x x x x x x i+1 k+1 ) = min exp ( 1 T k - 1 T k+1 )(E(x x x x x x x i+1 k ) -E(x x x x x x x i+1 k+1 )) , 1 . (2.41)
In our experiments, the set of temperatures is fixed, for all the PMF models, to 16 linearly spaced temperatures ranging from 1 to 20 (Salakhutdinov, 2009).

Algorithm 2.2: Tempered Gibbs sampler

Data: {T k } k=K k=1 , a set of ordered temperatures, {x x x x x x x 1 k } k=K k=1 , initial
k ∈ {K -1, . . . , 1} do x x x x x x x i+1 k ← x x x x x x x i+1 k+1 with probability α k (x x x x x x x i+1 k , x x x x x x x i+1 k+1 ) end i ← i + 1 end
The T-Gibbs sampler can replace the classic Gibbs sampler in the Marroquin algorithm for MPM computation [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF].

Remark:

The tempered Gibbs sampler is clearly much more computationally intensive than the classical Gibbs sampler. This constitutes an important current drawback despite the improved performance that we will study in the next section.

Remark: An additional experiment characterizing the T-Gibbs sampler is available in Appendix C.2.

Experiments and Results

Supervised segmentation of semi-real images with the PMF models

In this section the task of image segmentation is performed with the models from the PMF family in the context of supervised segmentation. The pair of complete data (x x x x x x x, y y y y y y y) is then available to estimate θ θ θ θ θ θ θ from it (with the estimators described in Section 2.4) before estimating x x x x x x x using y y y y y y y and θ θ θ θ θ θ θ. The PMF models are here tested in the case of supervised image segmentation according to the two segmentation algorithms MAP, MPM. The MPM will be computed either with the Marroquin's algorithm [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF] or the T-Gibbs algorithm. We compute the average error rate for each PMF model in the segmentation of a series of images from the dataset with varying noise levels. Recall that the noise is modeled by the realization of a GRF for the semi-real image formation. In the first case, the noise parameter that is varying is ∆µ = |µ 1µ 0 |, in the second case, the range r of the noise is the variable. The experiments are summed up in Figures 2. First of all, we note the overall reduced error rate thanks to the new pairwise models over the classical P-IN model. The GPMF model leading to the best segmentations for all noise levels. This behavior illustrates the capacity for the new model to take into account the correlated noise.

Note that in the case of varying ranges, for the two segmentation criteria (especially for the MAP), the GPMF and P-GMRF models perform worse than the P-IN model at the smallest correlation ranges. Two reasons might explain the phenomenon. First, the model definitions: when r → 0, the GPMF model does not totally tend towards the P-IN model, for example. Second, convergence issues with the Gibbs sampler can happen: the error rates vary importantly between the three algorithms used to compute the criteria6 . The MPM computed with T-Gibbs exhibits the best error rates for GPMF and P-GMRF (as expected from an improved sampling procedure). Notably, in the critical zone where r < 2, the averaged error rates are increased by several points for these two models. This issue is, however, unlikely to appear in practical cases when the noise is sure to be correlated, such as our applications. Indeed, we here highlight a theoretical drawback linked with the limits of stochastic Gibbs sampling.

Unsupervised segmentation on semi-real images

In this section, we address the problem of unsupervised image segmentation. We use the semi-real images previously introduced but consider only the observations Y Y Y Y Y Y Y = y y y y y y y. We here consider the GPMF model and the MPM criterion computed with Marroquin's algorithm: it is the best compromise between segmentation performance and aomputation time. To compare the GPMF model, we used three approaches developed for unsupervised image segmentation: • The KMeans clustering7 algorithm [START_REF] Arthur | Statistical inference for probabilistic functions of finite state Markov chains[END_REF].

• The pyImSegm8 (pyIS) segmentation algorithm proposed in [START_REF] Blei | Supervised and unsupervised segmentation using superpixels, model estimation, and graph cut[END_REF]. The core of this technique relies on a Markov Random Field energy minimization problem using super-pixel based and graph-cut based approaches.

• Based on the recent developments of [START_REF] Rubel | Additive Spatially Correlated Noise Suppression by Robust Block Matching and Adaptive 3D Filtering[END_REF] for correlated noise reduction (which were not available online), we propose to combine the Block-Matching 3D filter9 from [START_REF] Noel Cressie | Conditional-mean least-squares fitting of Gaussian Markov random fields to Gaussian fields[END_REF] and the graph-cut algorithm10 from [START_REF] Blei | Supervised and unsupervised segmentation using superpixels, model estimation, and graph cut[END_REF]). This approach is then the combination of a popular noise reduction technique and a popular graph-cut based segmentation. We call this approach BM3D+GC.

Let us now compare the unsupervised segmentations obtained with the GPMF model with the results given by other models. Figure 2.5 illustrates the segmentation performance of the models for a varying noise level. Figure 2.6 depicts some segmentations of images from the database. We notice, in all cases, the superiority of the new probabilistic model GPMF in its capacity to handle the correlated noise, improving the segmentation error and giving stable results. The GPMF model always performs best or equally best at all noise levels. In the best scenarii, the GPMF model increases the average error rate by about four points. The BM3D+GC and pyIS approaches are more unstable and perform worse when dealing with an image corrupted with correlated noise, despite manually tuned hyperparameters. Finally, note that, as expected, the overall error rates are higher in this unsupervised segmentation experiment than in the supervised segmentation experiment of Section 2.5.2.

On real world images

We now present the model in a real world application from the medical field, where unsupervised image segmentation needs to be done on strongly spatially corrupted data. The context of these experiments has been presented in Section 2.1. To the best of our knowledge, there is no counterpart of our model in the medical literature involving spatially-correlated noise for stent segmentation in medical images. We treat 512 × 512-pixel 2D images.

The goal of the experiment is to segment precisely the organic material from the background of the images. The two classes we wish to distinguish are visible by the naked eye, but they are very corrupted by the artifacts and For the ∆µ row, the other parameters of the GRF noise were fixed to σ = 0.5 and r = 2. For the r row, they were fixed to µ 0 = 0, µ 1 = 0.6 and σ = 0.5. The dashed red lines highlight the situations when the parameter configuration is identical in both graphs. As expected, in such cases, we notice that the models have identical performances relative to each other. are challenging to segment using automatic unsupervised algorithms. Such a processing would help identify the nature of the elements in contact with the stent, to precisely understand the possible interactions of the stent that could have lead to its failure. A dedicated model needs to be created since the artifacts hide the underlying material. Figure 2.7 depicts the results obtained by three methods described from the previous section (BM3D+GC, P-IN and GPMF). The results of the pyIS method were essentially similar to the results of BM3D+GC and are omitted here for brevity, as well as the results form the KMeans algorithm which are omitted because of the poor performance of the algorithm studied in Section 2.5.3. The P-IN and the BM3D+GC results are particularly prone to misclassifications because of the artifacts.

While the overall error rate is in favor of the GPMF model, it does not truly reflects the capacity of the model to resolve correlated noise and to offer a proper segmentation. Therefore we consider the organic material as the true class and the background as the false class. We then compute the False Negative (FN) and False Positive (FP) percentages of pixels in the areas of interest, around the stent where the correlations are the strongest. Those data are given in Table 2.3. We can see that the GPMF model best captures the correlated noise and resolves much of the stent artifact problem: it offers the best compromise between the FN and FP scores. In addition the results obtained with the GPMF model were also judged very satisfying by the pathologists.

The GPMF model and its application to real world images is further discussed in Chapter 5.

Conclusion

In this chapter we proposed a new probabilistic model and two new submodels which belong to the PMF family. We have seen that the PMF assumption can lead to modelization embedding much more complex correlations between random variables than the modelizations based on HMFs. This was used to construct a generative latent variable model with a GMRF based likelihood. By proposing a new parameter estimation procedure, the GPMF model then offered the best results, in the unsupervised segmentation of images corrupted with strongly correlated noise, when compared with other classical unsupervised segmentation techniques. : Unsupervised segmentations of organic material in corrupted Xrays images. For these 2 cases, which illustrate different explanted stents, we have: the real image y y y y y y y, the ground truth image x x x x x x x, the BM3D-GC segmentation, the P-IN segmentation and the GPMF segmentation. The segmentation criterion for the probabilistic models is the MPM. The stent parts (brightest pixels in y y y y y y y) were segmented beforehand by a thresholding technique and then considered as image borders. They appear in red on the segmented images and did not take part in the segmentation. The ground truth could be provided by experts since a histological analysis is available. Error rates with respect to the ground truth on the whole image appear in parenthesis.

Conclusion

As an application of this model to real world data, we were able to solve the problem of organic material segmentation in an image where strongly correlated noise is caused by metallic biomaterials. The GPMF model deals well with correlated noise and limits the number of misclassifications, as opposed to methods unsensitive to the spatial context of pixels.

Chapter 3.

Spatial Triplet Markov Trees

Introduction

This chapter introduces probabilistic models that belong to the Directed Graphical Model (DGM) family. We will study both direct and approximate inference. We have seen in the previous chapters some difficulties that arise when inference must be approximate. That is why, the main focus of this chapter is to develop a new model called Spatial Triplet Markov Tree (STMT) in which we elaborate enhanced correlations with the help of auxiliary random variables, but we preserve the possibility of performing the inference directly.

The enhancement of correlations in the Markov Tree (MT) structure are local, spatial, in the same sense of the correlations found in classical Markov Random Fields models (such as the models from the previous chapter). The STMT model is then particularly suited for image segmentation, we show that it performs better than the classical MT model. The presentation of the new model and its application to image segmentation is done in the first part of this chapter.

In the second part of the chapter, Variational Inference (VI) is developed more in depth. VI [START_REF] Steffen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF]) is an approach to perform (approximate) inference with a complex probability distribution where exact computations are not available. VI recasts the inference problem as an optimization problem which relies on a simpler distribution, the variational distribution, used to approximate the complex one. Additional approaches to VI consist in constructing a variational distribution with a certain complexity in itself (direct dependencies between random variables) which might, moreover, involve auxiliary random variables. Such developments might lead to further improvements in the inference (Agakov and [START_REF] Felix | An auxiliary variational method[END_REF]). We will apply these principles and we will also use VI as a way to better understand the STMT model.

Markov Tree models

This section focuses on particular MT models. The purpose of this section is to present the models that will be used in the context of Bayesian image segmentation (Section 1.5.2). We continue to denote X X X X X X X as the hidden process and Y Y Y Y Y Y Y as the observed process. We focus our work on dyadic and quadtree MTs whose layers of sites are S = {S1 , . . . , S L }.

Hidden Markov Trees

Recall the introduction of HMTs made in Section 1.6.1. The joint distribution of a HMT with Independent Noise (HMT-IN) is given in Equation 1.24. Figure 1.5 depicts the graphical model of a dyadic HMT-IN.

As mentioned earlier, the inference tasks can be done with exact computations in HMTs. This can be done with the Upward Downward (UD) algorithm 1 [START_REF] Monfrini | Image and signal restoration using pairwise Markov trees[END_REF] , ∀x s ∈ Ω, ∀s ∈ S, where s s s s s s s ++ denotes the set of all the descendants of s (vertices that can be reached from s by following the directed edges) [START_REF] Pieczynski | Arbres de Markov couple[END_REF].

Remark: Algorithm 3.1 is written for the HMT-IN model [START_REF] Monfrini | Image and signal restoration using pairwise Markov trees[END_REF][START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF] where observed variables are associated to hidden variables on the last layer only (as drawn in Figure 1.5). This algorithm is easily adapted if there is no, more or less observed variables (J.-B. [START_REF] Durand | Computational methods for hidden Markov tree models-An application to wavelet trees[END_REF].

Spatial Triplet Markov Trees

Spatial Triplet Markov Trees (STMTs) is a model of MTs based on a triplet Markov process that we have defined in Section 1.6.2. Along with the previously Backward algorithm used for Hidden Markov Chains, or as a principled message passing schedule in the Belief Propagation algorithm.

Chapter 3. Spatial Triplet Markov Trees introduced processes X X X X X X X and Y Y Y Y Y Y Y , let V V V V V V V be a discrete-valued random process with values in Λ N . V V V V V V V
plays the role of an auxiliary process. In the Bayesian approach to image segmentation that we follow, V (3.5) Equation 3.5 reveals (similarly to Equation 1.27 in the case of chains) the triplet assumption:

V V V V V V is part of the hidden process, along with X X X X X X X. Let T T T T T T T be a process such that T T T T T T T = (X X X X X X X, V V V V V V V , Y Y Y Y Y Y Y ).
T T T T T T T = (X X X X X X X, V V V V V V V , Y Y Y Y Y Y Y ) is a Markov process.
A generalized version of the UD algorithm (Algorithm 3.1) can be proposed for STMTs: it suffices to replace the occurences of the X X X X X X X process by the new hidden process

(X X X X X X X, V V V V V V V ).
Then, to recover information on the original process X X X X X X X, one needs to marginalize over the auxiliary random variables, ∀s ∈ S:

p(x s ) = vs∈Λ p(x s , V s = v s ). (3.6)
The rest of this section is dedicated to the definition of the STMT model in its dyadic and quadtree versions2 . The STMT model takes advantage of the triplet structure for exact inference with UD. It also takes advantage of auxiliary random variables to introduce much richer correlations between random variables as compared to the classical HMT model. More precisely, the model aims at introducing local, spatial correlations similarly to MRFs (Section 1.6.1). However, as already mentioned, the model preserves direct computations of the inference tasks, as opposed to MRFs.

Remark:

The STMT model will be described in a general form containing the visible process Y Y Y Y Y Y Y only at the last resolution and the Independent Noise hypothesis will always hold.

The dyadic model

We now describe the construction of the dyadic STMT model whose graphical representation is given in Figure 3.1 (in case of Independent Noise). We consider that the observations are only associated with vertices of the lowest layer, with the hypothesis of the independent noise.

A node at site s ∈ S \ S L is a triplet (X s , V ← s , V → s ); we then assume that we have the factorization:

p(t t t t t t t s |t t t t t t t s-) = p(x s , v v v v v v v s |x s -, v v v v v v v s -), = p(x s |x s -, v v v v v v v s -)p(v ← s |x s -, v v v v v v v s -)p(v → s |x s -, v v v v v v v s -). (3.7) However, a node at site s ∈ S L is a quadruplet (X s , V ← s , V → s , Y s ); we then 3.2. Markov Tree models (V ← r , X r , V → r ) (a)
← → s - L ← → s L R ← → s R X V
) = p(x s , v v v v v v v s , y s |x s -, v v v v v v v s -, y s -), = p(x s , v v v v v v v s |x s -, v v v v v v v s -, y s -)p(y s |x s , v v v v v v v s , x s -, v v v v v v v s -, y s -), = p(x s , v v v v v v v s |x s -, v v v v v v v s -)p(y s |x s ), = p(x s |x s -, v v v v v v v s -)p(v ← s |x s -, v v v v v v v s -)p(v → s |x s -, v v v v v v v s -)p(y s |x s ). (3.8)
To further define the transition laws of each of the variables, we define the notion of inner and outer variables for the V variables. We define that, within

Left (L) (resp. Right (R)) sons, V ← s L (resp. V → s R ) is an outer variable and V → s L (resp. V ← s L
) is an inner variable. Figure 3.2 illustrates these concepts for a particular node.

We now detail Equations 3.7 and 3.8, according to the variable type (left, right, inner or outer ). For X s sons:

p(x L s |x s -, v v v v v v v s -) = p(x L s |X s -= x s -, V s -= v ← s -), p(x R s |x s -, v v v v v v v s -) = p(x R s |X s -= x s -, V s -= v → s -), (3.9)
for inner V s sons:

p(v → s L |x s -, v v v v v v v s -) = p(v → s L |x s -= x s -, V s -= v → s -), p(v ← s R |x s -, v v v v v v v s -) = p(v ← s R |X s -= x s -, V s -= v ← s -), (3.10)
and for outer V s sons:

p(v ← s L |x s -, v v v v v v v s -) = p(v ← s L |X s -= v ← s -, V s -= x s -), p(v → s R |x s -, v v v v v v v s -) = p(v → s R |X s -= v → s -, V s -= x s -).
(3.11)

Remark: The conditioning of the V s variables can be seen as being the same conditioning the nearest X neighbor of X s on a same resolution. The triplet tree framework then provides a way to simulate X variables conditionally to the realizations of V variables which behave similarly to the neighboring X variables. This is notable since such links are strictly impossible in classical

Markov Tree models

Type of son at s Inner V s sons at s Outer I V s sons at s

Outer II V s sons at s NW V → s , V ց s , V ↓ s V ← s , V տ s , V ↑ s V ր s , V ւ s NE V ↓ s , V ւ s ,V ← s V ↑ s , V ր s , V → s V տ s , V ց s SE V ← s , V տ s , V ↑ s V → s , V ց s , V ↓ s V ր s , V ւ s SW V ↑ s , V ր s , V → s V ← s , V ↓ s , V ւ s V տ s , V ց s Table 3.
1.: Definitions of the inner, outer I and outer II V s sons for the quadtree model.

MTs.

Finally for the root node

s ∈ S 1 , p(x s , v v v v v v v s ) is chosen as a discrete distribution on Ω × Λ 2 .

The quadtree model

We now describe the construction of the quadtree STMT. Again, we consider that the observations are only associated with vertices of the lowest layer, with the hypothesis of the independent noise.

Each father s -has four sons s N W , s N E , s SE and s SW , where NW, NE, SE and SW stand respectly for North West, North East, South East and South West. For example, at site s N W , the random variables are X s N W and

V V V V V V V s N W . The composition of the octuplet V V V V V V V s N W is again spatially described by V V V V V V V s N W = (V ← s N W , V տ s N W , V ↑ s N W , V ր s N W , V → s N W , V ց s N W , V ↓ s N W , V ւ s N W ).
Similarly to the dyadic case we define inner and outer V variables for each type of father node. However the outer type is now divided in two categories outer I and outer II. The definitions are given in Table 3

, V V V V V V V s ).
The transition distributions are similar to Equation 3.7 (when s ∈ S \ S L ) and Equation 3.8 (when s ∈ S L , case with observed random variables3 ), however p

(x s , v v v v v v v s |x s -, v v v v v v v s -
) is now assumed to factorize as:

p(x s , v v v v v v v s |x s -, v v v v v v v s -) =p(x s |x s -, v v v v v v v s -)p(v ← s |x s -, v v v v v v v s -)p(v տ s |x s -, v v v v v v v s -) p(v ↑ s |x s -, v v v v v v v s -)p(v ր s |x s -, v v v v v v v s -)p(v → s |x s -, v v v v v v v s -) p(v ց s |x s -, v v v v v v v s -)p(v ↓ s |x s -, v v v v v v v s -)p(v ւ s |x s -, v v v v v v v s -).
(3.12)

This factorization is independent of the type of the father. However, we will now specify the conditioning variables in these transition laws and the latter will be dependent on the type of the father as in the dyadic case. Thus, we will propagate the homogeneity of the realizations at one resolution.

Chapter 3. Spatial Triplet Markov Trees

Each father s -has four sons s N W , s N E , s SE and s SW , with, for the X s sons:

         p(x s N W |x s -, v v v v v v v s -) = p(x s N W |x s -= x s -, V V V V V V V s -= (v ← s -, v տ s -, v ↑ s -)), p(x sNE |x s -, v v v v v v v s -) = p(x sNE |x s -= x s -, V V V V V V V s -= (v ↑ s -, v ր s -, v → s -)), p(x s SE |x s -, v v v v v v v s -) = p(x s SE |x s -= x s -, V V V V V V V s -= (v → s -, v ց s -, v ↓ s -)), p(x s SW |x s -, v v v v v v v s -) = p(x s SW |x s -= x s -, V V V V V V V s -= (v ↓ s -, v ւ s -, v ← s -)).
(3.13)

We now detail the transition laws for the inner V s sons at site s N W :

     p(v → s N W |x s -, v v v v v v v s -) = p(v → s N W |X s -= x s -, V V V V V V V s -= (v ↑ s -, v ր s -, v → s-)), p(v ց s N W |x s -, v v v v v v v s -) = p(v ց s N W |X s -= x s -, V V V V V V V s -= (v → s -, v ց s -, v ↓ s -)), p(v ↓ s N W |x s -, v v v v v v v s -) = p(v ↓ s N W |X s -= x s -, V V V V V V V s -= (v ↓ s -, v ւ s -, v ← s -)).
(3.14)

The same remark as for Equation 3.10 needs to be made for Equation 3.14: given a father s -and two sons s, s ′ , ∀s ∈ {N W, N E, SE, SW }, the inner V s sons take the same conditioning as the nearest X s ′ (neighbor of X s ) on the same resolution. Since inner V s will always have the same conditioning as a X s ′ whose father is also at site s -, we have s = s ′ . Outer I V s sons take the same conditioning as the nearest X s ′ but under a switching condition concerning the conditioning: the V s -with the same direction as V s is switched with X s -. This gives the transition laws of the outer I V s sons at site s N W :

     p(v ← s N W |x s -, v v v v v v v s -) = p(v ← s N W |X s -= v ← s -, V V V V V V V s -= (x s -, v տ s -, v ↑ s -)), p(v տ s N W |x s -, v v v v v v v s -) = p(v տ s N W |X s -= v տ s -, V V V V V V V s -= (v ← s -, x s -, v ↑ s )), p(v ↑ s N W |x s -, v v v v v v v s -) = p(v ↑ s N W |X s -= v ↑ s -, V V V V V V V s -= (v ← s -, v տ s -, x s -)).
(3.15)

Outer II V s sons take the same conditioning as the nearest X s ′ , for s and s ′ two nodes which are sons of s -(leading again to s = s ′ ). Moreover we also have a switching condition: the only V s -variable which plays a role in the conditioning of both X s and X s ′ is switched with X s -. This gives the transition laws for the outer II V s sons at site s N W :

p(v ր s N W |x s -, v v v v v v v s -) = p(v ր s N W |X s -= v ↑ s -, V V V V V V V s -= (x s -, v ր s -, v → s -)), p(v ւ s N W |x s -, v v v v v v v s -) = p(v ւ s N W |X s -= v ← s -, V V V V V V V s -= (x s -, v ↓ s -, v ւ s -)), (3.16)
To illustrate Equations 3.13, 3.14, 3.15 and 3.16 (and the other transition laws that have not been written), Figure 3.3 illustrates a transition between a node s -and its 4 sons s N W , s N E , s SE and s N W with the spatial context.

Finally for the root node s

∈ S 1 , p(x s , v v v v v v v s ) is chosen as a discrete distribution on Ω × Λ 8 .
Remark: In all the thesis, when dealing with image segmentation, the quadtree version of the STMT model (as well as the HMT model) will be used.

Image segmentation

← տ ↑ ր → ց ↓ ւ s - N W ← տ ↑ ր → ց ↓ ւ s N W N E ← տ ↑ ր → ց ↓ ւ s N E SW ← տ ↑ ր → ց ↓ ւ s SW SE ← տ ↑ ր → ց ↓ ւ s SE X V

Image segmentation

The models introduced in the previous section will be tested in the context of image segmentation. In this section we present the family of functions used to define the transitions in HMTs and in STMTs, as well as a method to estimate the parameters of these distributions in an unsupervised context.

The Potts-like transition distributions

Definition

Recall the Potts potential function (J. [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] defined in Equation 2.16. Let us define conditional probabilities in arborescences that follow this function.

Then, in the HMT model, we can form a transition distribution that we call Potts-like distribution:

p(x s |x s -) ∝ exp(αδ x s - xs ), (3.17)
where α ∈ R + is a parameter which regulates the similarity between two linked realizations of the hidden process in the HMT model. The higher α the more probable is the fact that the realizations will be equal. A Potts-like distribution for the STMT model is defined by:

p(x s |x s -, v v v v v v v s -) ∝ exp   αδ x s - xs + v s -∈v v v v v v v s - βδ v s - xs   , (3.18)
where α ∈ R + and β ∈ R + are parameters with the same role as before. In all the following work, the transitions of the HMT and STMT models will use the Potts-like transition distributions. Moreover, the transitions will be considered homogeneous: ∀s ∈ S, they are Potts-like transitions equally parametrized.

Remark: Note that when β = 0, Equations 3.17 and 3.18 are equal and the HMT and STMT models thus become equivalent.

Linear Least-Square estimator of the Potts-like transition parameter

We now establish an estimator for the α and the β parameters based on the Linear Least-Square (LLS) estimator first proposed in (Derin and H. Elliott, 1987). The method requires the complete data, hidden and observed, which correspond to the pair (X X X X X X X = x

x x x x x x, Y Y Y Y Y Y Y = y y y y y y y) for HMT and the triplet (X X X X X X X = x x x x x x x, V V V V V V V = v v v v v v v, Y Y Y Y Y Y Y = y y
y y y y y) for STMT. Note that the derivation of the LLS estimator for the granularity coefficient of the GPMF model (Section 2.4.1) is close to the equations that we now establish.

In the case of the α estimation for HMTs, we can write, ∀s ∈ S: Taking the natural logarithm on both sides, plugging the potential from Equation 3.17, using the definition of the independent Gaussian noise and rearranging, we finally get:

p(x s ,
α(δ x s - xs -δ x s - x ′ s ) = ln p(x s , x s -|y y y y y y y) p(x ′ s , x s -|y y y y y y y) -ln 1 2πσ 2 xs + (y s -µ xs ) 2 2σ 2 xs + ln 1 2πσ 2 x ′ s + (y s -µ x ′ s ) 2 2σ 2 xs .
(3.23)

Image segmentation

The probabilities of the type p(x s , x s -|y y y y y y y), ∀(x s , x s -) ∈ Ω 2 , are estimated with the frequency estimator. Then the same LLS procedure as described in Section 2.4.1 can be used to find α.

Then, following strictly the same procedure for the STMT model, starting from, ∀s ∈ S:

p(x s , x s -, v v v v v v v s -|y y y y y y y) p(x s -, v v v v v v v s -|y y y y y y y) = p(x s |x s -, v v v v v v v s -, y y y y y y y), (3.24)
we can write, in the case of the independent Gaussian noise, ∀s ∈ S,

∀(x s , x ′ s ) ∈ Ω 2 : α(δ x s - xs -δ x s - x ′ s ) + β   v s -∈v v v v v v v s - δ v s - xs -δ v s - x ′ s   = ln p(x s , x s -, v v v v v v v s -|y y y y y y y) p(x ′ s , x s -, v v v v v v v s -|y y y y y y y) - ln 1 2πσ 2 xs + (y s -µ xs ) 2 2σ 2 xs + ln 1 2πσ 2 x ′ s + (y s -µ x ′ s ) 2 2σ 2 xs . (3.25)
Again, the probabilities of the type p(x s , x s -, v v v v v v v s -|y y y y y y y) are estimated with the frequency etimator. The LLS procedure presented in Section 2.4.1 is adapted to handle two unknown variables and to estimate α and β.

Iterative Parameter Estimation for Trees

This section describes fully the Iterative Parameter Estimation for Trees (IPET) algorithm which will be used to estimate the parameters in unsupervised problems of image segmentation. This procedure is a deterministic iterative procedure close to the SEM algorithm [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF]. Moreover, we describe the IPET method for the STMT model, its adaptation for the HMT model is straightforward.

Without loss of generality, in the following developments we consider a two class segmentation problem, i.e., Ω = {ω 0 , ω 1 }, and the case of Gaussian independent noise parametrized by the class of the realization of the underlying hidden variable. This noise model adds four parameters to the modelizations:

µ 0 µ ω0 , µ 1 µ ω1 , σ 0 σ ω0 , σ 1 σ ω1 .
Finally for the STMT model, this leads to six parameters to estimate. Let θ θ θ θ θ θ θ ∈ Θ be the vector of parameters, then: favor of STMT). Note the strict equivalence between the models when β = 0, this fact already mentioned earlier is now confirmed empirically. Furthermore we can see that as β increases, i.e., the importance of the realization of the random variables which are neighbors to the father increases, the advantage of the STMT model increases. This was to be expected since the HMT-IN model cannot take into account such spatial and local correlations. However, a bigger α, i.e., the dependence on the realization of the father random variable increases, leads to a smaller advantage for the STMT model. This is due to the fact that, with a bigger α, the importance of the neighbors is lowered in the Potts-like transitions. In such cases the HMT-IN model is more relevant.

θ θ θ θ θ θ θ = {α, β, µ 0 , µ 1 , σ 0 , σ 1 }, with Θ = (R + ) 2 × (R) 2 × (R * + ) 3 . ( 3
In the second experiment, let us fix σ = σ 0 = σ 1 and ∆µ = |µ 0µ 1 |. Then σ and ∆µ are the varying parameters in the STMT samples. We plot in Figure 3.5 the quantity err HM T -INerr ST M T as a function of ∆µ and σ. We can see that for all the noise levels, the error rates are in favor of the STMT model. Moreover we note that the smaller the noise levels, the more similar the performances of both models.

Semi-real image experiment

In this section, we work in the context of unsupervised image segmentation on semi-real images from the "1070-Binary Shape Database" 5 already used in Chapter 2. We aim at highlighting the apparent capacity for the STMT model to take into account the spatial context. The observations y y y y y y y in this section are then the binary images from the data corrupted with an additive Gaussian noise. The noise level will be varying during our experiment.

The parameters of the HMT and STMT models are computed unsupervisedly with the IPET algorithm (Algorithm 3.2). The tree models are then compared to the results of a Hidden Markov Field with Independent Noise (HMF-IN). This probabilistic model is indeed a reference to model local spatial correlations and the motivation behind the STMT model. In the binary classification task we consider (Ω = {ω 0 , ω 1 }), the HMF-IN model is defined with a prior factor taken as a Potts factor (J. [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] with bias:

p(x s |x x x x x x x N X s ) = exp   -δ ω0 xs α -2 s ′ ∈N X s -δ x s ′ xs β   , (3.29)
where α ∈ R is the bias parameter6 and β ∈ R * + is the granularity parameter, and by the conditional likelihood factor: Figure 3.7 illustrates, for a varying noise level, the error rates in the unsupervised segmentation task for each model. This experiment clearly attests the advantage of the STMT model over the HMT model in unsupervised segmentation, the STMT model always performs equivalently or better than the HMT model. This result confirms the conclusions made on simulated data in Section 3.3.3. The HMFs remain the best performing model for σ ≤ 1.0. However, it then performs worse than STMTs at high noise levels, this might be due to the failures of the stochastic algorithms that we use for the models. Indeed, HMFs rely on indirect stochastic inference and parameter estimation which might fail at such noise levels. It is then notable that the STMT improvement is obtained while keeping a direct and analytical inference process. Figure 3.8 reports selected graphical examples of the experiment described above.

p(y s |x s ) = exp -ln( √ 2πσ) - ȳ2 s 2σ 2 . ( 3
Experiments of atherosclerotic calcification segmentation are given in Chapter 5.

Spatial Triplet Markov Trees for Auxiliary Variational Inference in Spatial Bayes Networks

This section now offers a more theoretical study of the distribution of the original hidden process X X X X X X X within the STMT model. In order to understand the much more complex and richer correlations that the STMT model introduces,

STMTs for Auxiliary Variational Inference in SBNs

we first define a new model of Bayesian network (or Directed Acyclic Graph (DAG)) which differs from the previous MT models. We call this new model Spatial Bayes Network (SBN) and we will study its relationship to the STMT model. More precisely, as will see in the next section about the model definition, SBNs exhibit directly the spatial correlations that we want to model (but cannot directly model without breaking the Markovian assumption on the process) in the STMT construction. This comes at the price of an inference process that cannot be done through direct computations. However, since the direct computations are available within the STMT model, the pair of models SBN/STMT can be used in the context of Variational Inference (VI) as we will see in Section 3.4.2.

Remark:

The following development focuses on a study of the hidden and auxiliary processes, therefore, there will be no mention of a visible process, and the latter will not be drawn in the graphical representations. This can also be interpreted as working in graphical models where all the variables are visible.

Spatial Bayes Networks

The SBN model is a DAG based on the MT model but SBNs contain semi cycles. SBNs are not arborescences anymore and the computation for the inference step must be approximate. The additional edges have the purpose to better capture local correlations between random variables.

The dyadic model

In order to distinguish between the two sons of a father node s -in a dyadic tree: let s L and s R be respectively the left and right son of s -. We also define s ← (resp. s → ) as the left (resp. right) neighboring node of s. Let v : S → S be a mapping from a node to a neighboring node of its father, such that:

v : s → (s -) ← if s is a left node, (s -) → if s is a right node. (3.31)
The SBN model (illustrated in Figure 3.9) has the following distribution:

p(x x x x x x x) =p(x r ) s∈ S p(x s |x s -, x v(s) ). (3.32)
The quadtree model

The development above, dedicated to dyadic trees, can be straightforwardly extended to quadtrees, to form the quadtree SBN. There are 4 types of nodes (except from the root), which we call North West (NW), North East (NE),

X r X s - X (s -) ← X (s -) → X s Figure 3
.9.: Graphical model corresponding to a dyadic SBN. As an illustration, a node s, its father s -, the left neighbor of its father (s -) ← , the right neighbor of its father (s -) → and the root node r have been annotated. South East (SE) and South West (SW). Now the v function is redefined as a mapping from S to S3 with:

v : s →          ((s -) ← , (s -) տ , (s -) ↑ ) if (s -) NW node, ((s -) ↑ , (s -) ր , (s -) → ) if (s -) NE node, ((s -) → , (s -) ց , (s -) ↓ ) if (s -) SE node, ((s -) ↓ , (s -) ւ , (s -) ← ) if (s -) SW node.
The joint distribution p(x x x x x x x) of the quadtree SBN uses this newly defined v function. It is then written as:

p(x x x x x x x) =p(x r ) s∈ S p(x s |x s -, x x x x x x x v(s) ). (3.33)
For clarity, we only illustrate locally the quadtree SBN model. Figure 3.10 depicts the relations between a father vertice and all its sons.

Exploring the model similarities by sampling

We propose to explore the capacity of the STMT, SBN and MT models to take into account the spatial context by clamped sampling in these models. Recall that in all this section, we do not consider any visible process (Y Y Y Y Y Y Y ). Clamped sampling in these layered models means that we fix a layer to a given realization and sample the next generation(s) by ancestral sampling. We fix a realization to be that of a MRF with a Potts potential since it is our reference for modeling local spatial interactions. Figure 3.11 illustrates the experiment and provides an averaged difference rate. From this score, and also the two visual examples of the figure, we can say that the SBN and STMT models seem to better capture the spatial context than the MT model whose realizations are farther away from the MRF reference. Indeed, the difference rate between the MT realizations and MRF reference is more than twice the difference rate between the STMT realizations and the MRF reference.

This experiment provides a first insight on the similarity between STMT and SBN. The next section and the rest of the chapter study this similarity more in depth in the context of Variational Inference.

Variational inference

Variational Inference (VI) [START_REF] Steffen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF] [START_REF] Michael I Jordan | An introduction to variational methods for graphical models[END_REF] (Murphy, 2012) [START_REF] Blei | Supervised and unsupervised segmentation using superpixels, model estimation, and graph cut[END_REF]) (C. Zhang et al., 2018) is an approach to perform inference in probabilistic models where direct computations are not feasible. In the general context the goal is to find a variational distribution q(x x x x x x x) which approximates well a posterior p(x x x x x x x|y y y y y y y) in which computations for inference need to be approximated. To do so, VI casts an optimization problem which aims at minimizing the Kullback Leibler (KL) divergence between the two distributions; it is denoted KL(q(x x x x x x x)||p(x x x x x x x|y y y y y y y)). By definition:

KL(q(x x x x x x x)||p(x x x x x x x|y y y y y y y)) = E x x x x x x x∼q(x x x x x x x) [log q(x x x x x x x)] -E x x x x x x x∼q(x x
x x x x x) [log p(x x x x x x x|y y y y y y y)].

(3.34)

The variables of this optimization problem are the factors of the q distribution. The structure of q is chosen so that inference in q is easier than in p. Most of the time, q is simple enough for inference tasks in q to be directly computable. Therefore, at the end of the optimization problem, when q is fitted to p in the sense of the KL divergence, the complex inference in p will be approximated by an easier inference in q. It can be shown [START_REF] Blei | Supervised and unsupervised segmentation using superpixels, model estimation, and graph cut[END_REF] that minimizing Equation 3.34 is equivalent to maximizing the Evidential Lower BOund (ELBO) quantity defined as:

ELBO(q) = log p(y y y y y y y) -KL(q(x x x x x x x)||p(x x x x x x x|y y y y y y y)),

= E x x x x x x x∼q(x x x x x x x) [log p(x x x x x x x, y y y y y y y)] -E x x x x x x x∼q(x x x x x x x) [log q(x x x x x x x)]. (3.35)
Remark: As stated above, the rest of this section will not consider any visible process Y Y Y Y Y Y Y , hence the conditioning on y y y y y y y in the equations, such as Equation 3.34, will drop. Note however that all the results could be readily extended to integrate a visible process. Then, we have that:

ELBO(q) = E x x x x x x x∼q(x x x x x x x) [log p(x x x x x x x)] -E x x x x x x x∼q(x x
x x x x x) [log q(x x x x x x x)], = -KL(q(x x x x x x x)||p(x x x x x x x)). (3.36) Therefore, in what follows, we will refer to the optimization problem as maximizing the opposite of the KL divergence:

-KL(q(x x x x x x x)||p(x x x x x x x)) = E x x x x x x x∼q(x x x x x x x) [log p(x x x x x x x)] -E x x x x x x x∼q(x x x x x x x) [log q(x x x x x x x)]. ( 3

.37)

Remark: The rest of this section will consider the MT and SBN models in their dyadic version. All the results can also be written for the quadtree versions of the models.

In our case, p(x x x x x x x) is the joint distribution of a dyadic SBN given in Equation 3.32. In the context of structured VI (introduced in the next paragraph), we will reparametrize and work with another representation of Equation 3.32. Indeed, we need to use a concise notation embedding the notion of clusters of

STMTs for Auxiliary Variational Inference in SBNs

variables for each of the terms of the product. The notations are adapted from (C. [START_REF] Besag | Nearest-neighbour systems and the auto-logistic model for binary data[END_REF]. We can write Equation 3.32 as: p(x x x x x x x) = ds p ds (3.38) where d s represents the cluster of variables (x s , x s -, x v(s -) ) (note that these clusters overlap) and p ds is a shortcut for p(d s ). In the following we develop the necessary material to approximate Equation 3.32 by a variational distribution:

1. with totally independent random variables (Section 3.4.3).

2. with a MT structure (Section 3.4.4).

3. with a STMT structure (Section 3.4.5).

The first option is the most popular and is called Mean Field VI [START_REF] Michael I Jordan | An introduction to variational methods for graphical models[END_REF]. The two other options give the variational distribution a structure, it is called structured VI (Wiegerinck, 2000) (C. Bishop and Winn, 2003) (Z. Ghahramani and M. I. Jordan, 1996) [START_REF] Olariu | Modified variational Bayes EM estimation of hidden Markov tree model of cell lineages[END_REF]. It consists in using a structured variational distribution, but which is still simpler than the distribution p, in order to better approximate the correlations in the original intricate distribution. Structured VI can lead to dramatic increases in performances thanks to the improved modeling of the correlations. In the structured VI scenarii of this study, we take advantage of the MT and STMT structures and of the direct inference computations. Algorithmically, the steps of the VI procedure for a variational distribution q can be summarized:

1. Initialize all the factors of q.

2. Update each factor of q with the expression that minimizes the KL divergence.

3. Check convergence and repeat step 2 if needed.

The convergence can be assessed, e.g., by monitoring the values of the cost function (the KL divergence) or monitoring stationarity in the estimated variational parameters.

Remark: Once the terms of the variational distribution are learnt, a final inference task still needs to be done. For example, one could compute the MPM of the variational distribution, which would be used to approximate the MPM of the original intractable distribution.

Mean Field Variational Inference in SBNs

Let q be the variational distribution. If q is subject to the Mean Field (MF) assumption, all the random variables are independent:

q(x x x x x x x) = i q i (x i ). (3.39)
The quantity to maximize (Equation 3.37) becomes in this case:

-KL(q(x x x x x x x)||p(x x x x x x x)) = E x x x x x x x∼q(x x x x x x x) [log p(x x x x x x x)log q(x x x x x x x)],

= x x x x x x x i q i (x i ) log p(x x x x x x x) - k log q k (x k ) , (3.40)
where p is the joint dyadic SBN distribution 3.32. MF is the most popular form of VI. The derivation of the update equations to solve the MF optimization problem is done in Appendix D.1.

Markov Tree Variational Inference in SBNs

Let r be the variational distribution. r is here defined with the structure of a MT, it then follows a distribution of the form of Equation 1.23. Using a MT structure as a posterior to approximate the complex posterior p is a compromise between additional structure that might better reflect the correlations in p and tractability. Indeed, once the posterior transitions are learnt through the VI optimization, the posterior marginals are easily computable with the UD algorithm.

We then reparametrize the MT expression in terms of clusters:

r(x x x x x x x) = cs r cs , (3.41)
where c s represents the cluster of variables (x s , x s -). Note that the clusters overlap.

The quantity to maximize (Equation 3.37) becomes in this section:

-KL(r(x x x x x x x)||p(x x x x x x x)) = E x x x x x x x∼r(x x x x x x x) [log p(x x x x x x x)log r(x x x x x x x)],

= x x x x x x x cs r cs log p(x x x x x x x) - c k log r c k . (3.42)
The derivation of the update equations for this optimization problem is given in Appendix D.2.

Auxiliary variable Variational Inference

We now study VI with auxiliary random variables whose first occurence in the literature is in (Agakov and [START_REF] Felix | An auxiliary variational method[END_REF]. Let V V V V V V V be the auxiliary process, both the target distribution, p, and the variational distribution, which we call here t, will be augmented with variables from V

V V V V V V 7 .
Remark: As stated in (Agakov and [START_REF] Felix | An auxiliary variational method[END_REF], in this thesis, we augment p with the condition that p(x x x

x x x x, v v v v v v v) = p(x x x x x x x)p(v v v v v v v|x x x x x x x)
. This ensures an easy 3.4. STMTs for Auxiliary Variational Inference in SBNs marginalization of the V V V V V V V process in order to compare the distributions of the X X X X X X X process. In Section 3.4.7 we show how to integrate auxiliary random variables in p, following the previous joint distribution decomposition when p is a SBN.

Then, considering again temporarily observed random variables for generality, the new KL divergence with auxiliary random variables to minimize is:

KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v|y y y y y y y)) = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(x x x x x x x, v v v v v v v)] -E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x x x x x x, v v v v v v v|y y y y y y y)], = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(x x x x x x x, v v v v v v v)] -E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x
x x x x x|y y y y y y y)]-

E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(v v v v v v v|x x x x x x
x, y y y y y y y)].

(3.43)

Following (Kingma, 2017), we give the relation with the previous objective (without auxiliary variables), i.e., the KL divergence of Equation 3.34:

KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v|y y y y y y y)) = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(x x x x x x x, v v v v v v v)] -E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x
x x x x x|y y y y y y y)]-

E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(v v v v v v v|x x x x x x x, y y y y y y y)] = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(x x x x x x x)] + E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(v v v v v v v|x x x x x x x)]- E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x x x x x x|y y y y y y y)] -E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(v v v v v v v|x x x x x x x, y y y y y y y)] = KL(t(x x x x x x x)||p(x x x x x x x|y y y y y y y)) + KL(t(v v v v v v v|x x x x x x x)||p(v v v v v v v|x x x x x x x, y y y y y y y)) ≥0 (non-negativity)
, ≥ KL(t(x x x x x x x)||p(x x x x x x x|y y y y y y y)). (3.44) This allows us to draw a conclusion similar to (Kingma, 2017): the use of auxiliary random variables leads to the minimization of an objective function that is worse as concerns the original hidden variables of interest. However, the use of auxiliary variables offers much more flexibility in the modelization and enables the practitioner to add rich correlations between the variables of interest. The use of a well structured distribution also reduces local optima, which often helps since VI leads, in general, to a non convex optimization problem (C. Zhang et al., 2018). Those advantages outweight, in practice, the fact that we are optimizing an objective function with a worse lower bound.

Remark:

Note that, in Equation 3.44, we also used the factorization

t(x x x x x x x, v v v v v v v) = t(x x x x x x x)t(v v v v v v v|x x x x x x x).
However, in the case of a STMT variational distribution t, an explicit expression t(x x x x x x x) as not been established yet and will be part of further research on the topic. Therefore Equation 3.44 constitutes here a theoretical insight but cannot be used directly.

Remark:

The conditioning on y y y y y y y that appears in Equations 3.45 and 3.44 will disappear in the following development. Then, for reasons seen in the previous section, we will refer to the optimization problem as maximizing the opposite of the KL divergence with auxiliary random variables:

-KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v)) = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x x x x x x, v v v v v v v)] -E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log t(x x x x x x x, v v v v v v v)].
(3.45)

STMT auxiliary variable Variational Inference in SBNs

In this section we propose a STMT distribution t as the variational distribution for auxiliary variable VI in SBNs. We do not include the observed process in the modelization. Then we have: 3.46) which corresponds to Equation 3.7 without the observed process.

t(x x x x x x x, v v v v v v v) = s∈S t(x s , v v v v v v v s |x s -, v v v v v v v s -), = s∈S t(x s |x s -, v v v v v v v 1 s -)t(v ← s |x s , v v v v v v v 2 s -)t(v → s |x s , v v v v v v v 3 s -), (
V V V V V V V 1 s -, V V V V V V V 2 s - and V V V V V V V 3 s -are different subsets of V V V V V V V s -.
Recall that at the end of the VI procedure, exact marginal computation can be done, again with a generalized UD algorithm.

For the cluster parametrization, we have:

t(x x x x x x x, v v v v v v v) = cs t cs c ′ s t c ′ s c ′′ s t c ′′ s , (3.47) with c s = (x s , x s -, v v v v v v v 1 s -), c ′ s = (v ← s , x s -, v v v v v v v 2 s -) and c ′′ s = (v → s , x s -, v v v v v v v 3 s -).
As previously mentioned, the clusters overlap.

Remark: We separated the auxiliary process, at each site, using their Left or Right attribute. However this choice is arbitrary and it is only done here for clarity. It can be omitted without modifying the optimization problem.

We want to maximize Equation 3.45, which becomes, with p the distribution of a SBN augmented with auxiliary random variables8 :

-KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v)) = E (x x x x x x x,v v v v v v v)∼t(x x x x x x x,v v v v v v v) [log p(x x x x x x x, v v v v v v v) -log t(x x x x x x x, v v v v v v v)], = x x x x x x x,v v v v v v v cs t cs c ′ s t c ′ s c ′′ s t c ′′ s log p(x x x x x x x)- c k log t c k + c ′ k log t c ′ k + c ′′ k log t c ′′ k . (3.48)
The derivations of the update equations are given in Appendix D.3.

STMTs for Auxiliary Variational Inference in SBNs

Experiments & Results

Small SBN example

We now consider variational inference on the small SBN network given in Figure 3.12a. Similar experiments have been conducted in the same context, to evaluate a VI approximation, for example in [START_REF] Steffen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF]. Due to its small size the SBN of Figure 3.12a represents a slightly modified probability distribution from that the SBN used up to now. Indeed, we needed to treat in a specific fashion the root node to induce SBN-like correlations on a network with three layers only. It follows that p has the following distribution: The developments of the previous sections are used to conduct VI with MT and STMT. Note that for the STMT VI, we also need to provide the SBN with auxiliary nodes [START_REF] Felix | An auxiliary variational method[END_REF]. We need to keep the property that p

p(a, a ← , a → , b, c, d, e, f, g) =p(a)p(a ← )p(a → )p(b|a, a ← )p(c|a, a → ) p(d|b)p(e|b, c)p(f |c, b)p(G|c). ( 3
(x x x x x x x, v v v v v v v) = p(x x x x x x x)p(v v v v v v v|x x x x x x x), where x x x x x x x = {a, a ← , a → , b, c, d, e, f, g} and v v v v v v v = {b ← , b → , c ← , c → , d ← , d → , e ← , e → , f ← , f → , g ← , g → } in order to ensure that p(x x
x x x x x) (Equation 3.49) is the same between the three VI procedures. In STMT VI, p is then defined so that the auxiliary random variables have the same structure of dependencies in p as in t, the STMT variational distribution. We then have:

p(x x x x x x x, v v v v v v v) =p(x x x x x x x)p(b ← |b, a → )p(b → |a, a → )p(c ← |a, a ← )p(c → |a, a ← ) p(d ← |c → , b ← )p(d → |b ← , c → )p(e ← |b → , c ← )p(e → |c ← , b → ) p(g ← |c → , b ← )p(g → |b ← , c → )p(f → |b → , c ← )p(f ← |c ← , b → ), (3.50) with                    p(b|a, a ← ) = p(b ← |a, a ← ) = p(c ← |a, a ← ), p(c|a, a → ) = p(b → |a, a → ) = p(c → |a, a → ), p(d|b) = p(d ← |b ← ) = p(e ← |b ← ), p(e|b, c) = p(d → |b ← , b → ) = p(f ← |c ← , b → ), p(f |c, b) = p(g ← |c → , c ← ) = p(e → |b → , c ← ), p(g|c) = p(f → |c → ) = p(g → |c → ).
(3.51)

The model p with auxiliary nodes is described in Figure 3.12b. The random variables are chosen with values in {0, 1}. Our goal is to estimate the true marginals p(x = 1), ∀x ∈ x x x x x x x of the SBN of Figure 3.12a. In this synthetic example those true marginals are given as well as the transitions in p (we do not cover the parameter estimation problem). Moreover, the transitions of Equation 3.49 which also totally define Equation 3.50 are taken randomly.

A A ← A → B C D E F G (a) A A ← A → B B ← B → C C ← C → D D ← D → E E ← E → F F ← F → G G ← G → (b)
A A ← A → B C D E F G (a) A A ← A → B B D E F G (b) A A ← A → B B ← B → C C ← C → D D ← D → E E ← E → F F ← F → G G ← G → (c)

Results

We define a marginal error for a variational distribution q and a random variable x by e q (x) = p(x = 1)q(x = 1). These errors are computed and stored over 1000 different SBNs p whose transitions are randomly chosen. Figure 3.14 depicts the values of the KL divergence of the three VI procedures. We see a rapid convergence which is related to the small size of the considered SBN. We then choose to stop the VI process after 30 iterations. The value for MF VI and MT VI are comparable and are shown on the same graph. The minimization is better in the case of MT VI, the structured VI; this is reflected in Figure 3.15 which illustrates the goodness of the estimated marginals. The boxplots of Figure 3.15 illustrates the interest of the STMT structure to approximate the marginals of SBN: this VI procedure exhibits in all cases the smallest error with respect to the true marginals. In absolute value, this error is remarkably very small. We can also conclude that progressively adding structure increases the quality of the approximation since MF VI performs worse than MT VI which performs in turn worse than STMT VI. Note that the left/right symmetry of the SBN can be found in the behavior of the error rates: D is similar to G, E is similar to F , and so on.

Conclusion

This chapter was dedicated to the development, study and application of a new probabilistic model, the STMT model. This model combines the benefits of the triplet structure and auxiliary random variables modeling to introduce rich correlations between random variables while preserving exact inference computations. More precisely, we have seen that the correlations modeled are close to the spatial correlations of MRFs, making the new model particularly suited for image processing. In particular, the STMT model outperforms the classical HMT model.

Conclusion

In the second part of this chapter, the developement of the new SBN model and the use of the STMT/SBN pair of models in the context of VI offers an even better understanding of the correlations among the STMT model. The correlations of the X X X X X X X process in SBNs are experimentally found to be very close to the correlations of X X X X X X X in STMTs. Hence the richness of the STMT model.

Introduction

Supervised segmentation

In this chapter we study some supervised segmentation approaches. To this extent, the developed methodologies will differ from those of the previous chapters which focused on unsupervised segmentation methods. Deep Learning (DL) [START_REF] Goodfellow | Deep learning[END_REF] will also be introduced in this chapter; it is the core segmentation approach for the practical problem that we will address. Subsequent to the DL approach, in order to refine the results, we rely on a discriminative UGM model called Conditional Random Fields (CRFs) [START_REF] Suetens | An introduction to conditional random fields[END_REF].

Remark: The content of this chapter has been developed after the creation of an annotated dataset of mCT images during the course of the thesis. Thanks to this valuable dataset we were able to study the complementarity of DL approaches with probabilistic graphical models. As we will see in this chapter, state of the art approaches in medical image segmentation often rely on mixing both of these approaches.

The purpose of the dataset is to gather knowledge to answer the problem of histological classification of the biological components on the mCT images of femoropopliteal arteries (see Figure 0.9 for the location of these arteries in the lower limbs). The end goal of the project is to perform a first automatic 3D histological analysis of the explant by just using a mCT image. Such an algorithm could save time, money and help process a growing number of biological data. Figure 4.1 depicts the dataset. To the best of our knowledge, there exists no similar study in the literature.

The whole project was created from scratch at the Geprovas laboratory since 2018: starting from the explantations of the arteries from human bodies to the output of the segmentation algorithms. The project is presented in two parts in this thesis. The rest of this chapter focuses on the mathematical aspects of this study, while details about the medical background of the project, the protocols and the full results are given in Chapter 5.

Deep Learning approaches for histopathological image analysis

During the last decade, the interest of the researchers towards DL approaches have skyrocketed, particularly in the field of image segmentation and analysis. The mathematical basics of DL is given in Section 4.2.1. Our study deals with the application of DL to medical image processing, and more precisely, histopathological image processing. Reviews of image processing techniques in this field are published on a regular basis, and more frequently since the beginning of the DL era which was popularized in [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. One can compare a pre-DL era review of histopathological image analysis such In the context of histopathological image segmentation with DL, Convolutional Neural Networks (CNNs) are the reference approach which yields state of the art results (Srinidhi et al., 2019) (Haque and[START_REF] Hanzouli-Ben | A framework based on hidden Markov trees for multimodal PET/CT image co-segmentation[END_REF]. The seminal paper [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] demonstrated that these approaches can learn features in the data that can outperfom expert modelizations from researchers which were traditionally proposed in machine learning. The mathematical foundations and a description of the CNN we use in our study are given in Section 4.2.

3D data segmentation with deep learning approaches

To conclude this introduction, we now focus on a particular point which is crucial to our application and for which CNNs currently reach their limits: the problem of volumetric semantic segmentation. Most often in the literature, CNNs work hand-in-hand with another machine learning approach in order to perform the volumetric segmentation. The final goal of this study is also to reconstruct a 3D volume using the 2D annotations that we possess. This is unsurprisingly a common topic of interest in medical imaging. In general, the two main difficulties in this problem are, first, that the 3D annotations are almost inexistent, mainly because they are too complicated to make for human experts. Second, working with fully 3D CNNs is also computationally prohibitive. Several approaches built on 2D CNNs or related models have been proposed to provide the final 3D segmentations, we review them now.

The first type of approaches is to train directly a CNN to perform 3D convolutions to reconstruct a final 3D volume. The methodology can use a single fully 3D CNN (Çiçek et al., 2016) or use the 3D convolutions in hybrid approaches making use of previous results from 2D CNNs [START_REF] Yang | Multicomponent deformable models coupled with 2D-3D U-Net for automated probabilistic segmentation of cardiac walls and blood[END_REF]) (X. [START_REF] Li | H-DenseUNet: hybrid densely connected UNet for liver and tumor segmentation from CT volumes[END_REF]. Then because of the computational cost of 3D convolutions, one can note 2.5D approaches which aim at integrating spatial data with strong restrictions. For example, one can find propositions to consider successively small 3D chunks of consecutive slices [START_REF] Ben-Cohen | Fully convolutional network for liver segmentation and lesions detection[END_REF]. One also finds the idea of combining segmentations from the three axes of the data cube to produce a segmentation that can take into account more spatial context than a single 2D segmentation [START_REF] Zhao | 3D brain tumor segmentation through integrating multiple 2D FCNNs[END_REF]. Another type of popular approaches are based on Recurrent Neural Networks (RNNs) [START_REF] Goodfellow | Deep learning[END_REF] which process iteratively the slices through the datacube in successive passes to provide a segmentation taking into account all the axes [START_REF] Alexey | Deep sequential segmentation of organs in volumetric medical scans[END_REF] [START_REF] Cai | Improving deep pancreas segmentation in CT and MRI images via recurrent neural contextual learning and direct loss function[END_REF]. The final methodology that has been popularized during the last years and that we will follow is the 3D refinement with fully-connected Gaussian CRFs [START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF] of 2D segmentations produced by a CNN. The first occurence of such a work can be found in [START_REF] Kamnitsas | Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation[END_REF]. In Section 4.3, we present an approach built on a 3D refinement of the 2D dimensional segmentations with the fully-connected Gaussian CRFs.

Convolutional Neural Networks

Remark: A mathematical equivalence between computations in RNNs and Mean Field inference computations in CRFs has been demonstrated [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF]. This thus brings these two approaches closer.

Convolutional Neural Networks

The Deep-Learning approach

DL is based on the paradigm that computers learn by themselves from the data the features that optimally answer the problem [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF]. A DL model is composed of many layers that transform the input data up to an output layer which realizes the required task. Among DL models in the context of unsupervised learning (not treated in this thesis) one can cite Restricted Boltzmann Machines [START_REF] Hinton | A practical guide to training restricted Boltzmann machines[END_REF], Deep Belief Networks (Bengio et al., 2007), or Variational Autoencoders (Kingma and[START_REF] Diederik | Auto-encoding variational bayes[END_REF]. CNNs [START_REF] Fukushima | Neocognitron: A self-organizing neural network model for a mechanism of visual pattern recognition[END_REF] [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] are deep models trained in a supervised setting. CNNs are the most popular deep models for medical image processing, they are described in 4.2.2. Another supervised technique for this task are the RNN [START_REF] Hochreiter | Long short-term memory[END_REF] approach which is gaining popularity [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF]. RNNs are nonetheless out of the scope of this thesis.

Network definition

The U-Net network

In this work we construct a CNN based on the U-Net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. Since its conception, the U-Net network has yielded state of the art results in the field of hispathological slice analysis and has become a reference approach [START_REF] Jim | U-Net: deep learning for cell counting, detection, and morphometry[END_REF] upon which many other architectures were proposed [START_REF] Jabehdar Maralani | The one hundred layers tiramisu: Fully convolutional densenets for semantic segmentation[END_REF] [START_REF] Oktay | Attention u-net: Learning where to look for the pancreas[END_REF] [START_REF] Zhou | Unet++: A nested u-net architecture for medical image segmentation[END_REF]. The U-Net architecture is depicted in Figure 4.2.

The U-Net network retains the key elements of a CNN, particularly in the first part of the network, the contractive path or encoder. Convolution layers are the foundations of the network, they are parametrized by weights and biases and play the role of filters that extract features from the data. Each convolution layer is followed by a non-linear activation layer and a max-pooling layer. The latter successively reduces the dimension of the previous filtering operation, or feature map. This dimension reduction enables, in particular, to learn higher-level features (at different depths), which are likely to be more invariant to each training sample. Thus, the CNN predictions are more prone to generalize on unseen data. The expansive path, or decoder, is the second part of the network specific to network architectures from the U-Net family. In the decoder we find upsampling and concatenation operations which combine, respectively, the contextual, higher-order, features acquired in the encoder part with a precise higher resolution image. The network ends with a feature map of the size of the original image where all the pixels are assigned a score of The network presented here has input images of size 256 × 256 (which is adjustable). Note that the Batch Normalization layer is not present in the original article [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. Following the contractive path, the layers lose width and height but gain depth: we build a more condensed and higherlevel representation of the data.

Convolutional Neural Networks

belonging to a certain class.

Remark: In our U-Net modelization we added Batch Normalization layers [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] after each convolution layer.

U-Net training and evaluation

The U-Net network is trained in a supervised manner by solving an optimization problem over the parameters (weights and biases) of the network with the help of the back-propagation algorithm [START_REF] David E Rumelhart | Learning representations by back-propagating errors[END_REF]. The loss function of the optimization problem varies depending on the problem. The Dice loss function is a popular approach in the case of imbalanced dataset1 such as ours [START_REF] Sadegh | Tversky loss function for image segmentation using 3D fully convolutional deep networks[END_REF] [START_REF] Sudre | Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations[END_REF]. The Dice score D is defined by:

D(A, B) = 2|A ∩ B| |A| + |B| , (4.1)
where A and B are two images. We have the property that ∀A, B, D(A, B) ∈ [0, 1]. In the case where A and B are two segmentation maps for a same semantic class, a perfect segmentation, i.e., a perfect superposition of A and B, yields D(A, B) = 1. In the worst case, i.e., when A and B have no common pixel, we have D(A, B) = 0. From the Dice score we form the Dice loss D as

∀A, B, D(A, B) = 1 -D(A, B) ∈ [0, 1].
The Dice metric has the advantage of being insensitive to the absolute number of pixels in a class, hence its popularity in case of imbalanced data set. For a multi-class segmentation problem, we consider a generalized Dice score, the mean Dice score which is defined by the mean of the Dice scores [START_REF] Fidon | Generalised wasserstein dice score for imbalanced multi-class segmentation using holistic convolutional networks[END_REF]. For a L-label problem, where the segmentation maps to evaluate are organized as A A A A A A A = (A 1 , . . . , A L ) and B B B B B B B = (B 1 , . . . , B L ), the mean Dice score is defined by:

mD(A A A A A A A, B B B B B B B) = 1 L l∈L D(A l , B l ). (4.2) Since ∀A l , B l , D(A l , B l ) ∈ [0, 1], we have that mD(A A A A A A A, B B B B B B B) ∈ [0, 1]. The mean Dice loss is defined by ∀A A A A A A A, B B B B B B B, m D(A A A A A A A, B B B B B B B) = 1 -mD(A A A A A A A, B B B B B B B
). The latter is the loss function we use in the rest of this thesis.

The training is done on the training set which is composed of original slices but also augmented slices which come from data augmentation techniques (Shorten and Khoshgoftaar, 2019) (rotation, rescaling and gamma correction) applied to the original slices. Medical applications rely heavily on data augmentation techniques since medical data is often scarce, hard to gather and complex to annotate [START_REF] Morra | Artificial Intelligence in Medical Imaging: From Theory to Clinical Practice[END_REF]. A small portion2 of the training set, the validation set, is kept apart and is not used in the back-propagation algorithm. It serves as a means to tune hyperparameters of the neural network such as the number of iterations, or epoches, in the learning process.

Once the training is achieved the network is evaluated on a third dataset, the test set, which have been totally held out during the training process not to introduce any biases in the performance evaluation. The score obtained on the test set then reflects the generalization capability of the network to segment unseen data.

Fully-connected Conditional Random Fields

Model definition

In this section we define the CRF model from [START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF] that we use in a post-processing step to reconstruct segmented 3D volumes. We call this model fully-connected CRF (fcCRF). Recall the introduction on discriminative models from Section 1.5.3. The fcCRF model is a discriminative model that we then define with a Gibbs distribution:

p(x x
x x x x x|y y y y y y y) = 1 Z exp (-E(x x x x x x x|y y y y y y y)) ,

where Z is a normalization constant. The energy function E(x x x x x x x|y y y y y y y) is composed of unary and pairwise factors:

E(x x x x x x x|y y y y y y y) = s∈S ψ u (x s ) + (s,s ′ )∈S 2 ψ p (x s , x s ′ ). (4.4)
To avoid complicated expressions, the dependency of the factors on the observations Y Y Y Y Y Y Y = y y y y y y y will not explicitly appear. The unary potentials ψ u (x s ), ∀s ∈ S, ∀x s ∈ Ω, are taken as the label map outputs from the CNN classifier described in the previous section. The pairwise potentials are defined, ∀(s, s ′ ) ∈ S 2 , ∀(x s , x s ′ ) ∈ Ω 2 , by:

ψ p (x s , x ′ s ) = (1 -δ x s ′ xs ) 2 m=1 w m k m (f f f f f f f s , f f f f f f f s ′ ), = (1 -δ x s ′ xs ) w 1 exp - |s -s ′ | 2 2θ 2 α - |y s -y s ′ | 2 2θ 2 β k1(f f f f f f fs,f f f f f f f s ′ ) + w 2 exp - |s -s ′ | 2 2θ 2 γ k2(f f f f f f fs,f f f f f f f s ′ ) . (4.5)
f f f f f f f s , ∀s ∈ S
, is a feature vector whose components are the location, s, and the observed value at this location, y s . The weights w 1 and w 2 multiply, respectively, an appearance kernel, k 1 , parameterized by θ α and θ β , and a smoothness kernel, k 2 parametrized by θ γ . More precisely, the appearance kernel is the kernel used in bilateral filters [START_REF] Paris | Bilateral filtering: Theory and applications[END_REF] which is the 4.3. Fully-connected Conditional Random Fields product of a Gaussian parametrized by the spatial distance and a Gaussian parametrized on the pixel intensity distance. This results in a filter capable of smoothing while preserving edges.

Optimized Mean Field Variational Inference for fcCRFs

The Mean Field Variational Inference algorithm

In this section, we study an approximate inference procedure for estimating marginals in the fcCRF model. The fully-connected property of the networks forbids any kind of exact computations and requires several approximations for the problem to be tractable. The standard approach is to perform Mean Field (MF) Variational Inference (VI) in the fcCRF model. We already introduced MF VI in Section 3.4.3. Here, observed random variables need to be considered. We minimize the quantity KL(q(x x x x x x x)||p(x x x x x x x|y y y y y y y)), or, equivalently, maximize its opposite, with p(x x x x x x x|y y y y y y y) given in Equation 4.3 and q(x x x x x x x) = i q i (x i ). By following similar steps as in Appendix D.1 but integrating visible random variables, and using the energy from Equation 4.4, we find that3 , ∀j ∈ S, ∀x j ∈ Ω:

q j (x j ) = 1 Z j exp   -ψ u (x j ) - x ′ ∈Ω (1 -δ x ′ xj ) 2 m=1 w m i =j k m (f f f f f f f j , f f f f f f f i )q i (x ′ )   , (4.6
) where Z j is a normalization constant. The segmentation is then performed following the MPM criterion such that, ∀j ∈ S:

xMF V I j = argmax x∈Ω q j (x j ). (4.7)
The particular ordering of the summations in Equation 4.6 enables a MF algorithm divided in three steps which can be further optimized. The MF algorithm is given in Algorithm 4.1.

Remark:

The update of all the q j terms in Algorithm 4.1 can be done in a parallel fashion.

Optimization via approximate convolutions

In this section, we study the optimization proposed in [START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF] for an efficient MF inference in fcCRFs. In terms of computational performances, the bottleneck of Algorithm 4.1 is the message passing step. Algorithm 4.1: MF VI for fcCRFs [START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF] Data: w 1 , w 2 , θ α , θ β , θ γ , the model parameters, y y y y y y y, the observations. Result: q(x x x x x x x), the MF approximation of p(x x x x x x x|y y y y y y y). /* Initialize q */ q j (x j ) ← exp (-ψ u (x j )) , ∀j ∈ S, ∀x j ∈ Ω while convergence is not attained do for j ∈ S do

/* Message passing */ for x ′ ∈ Ω do for m ∈ {1, 2} do qj,m (x ′ ) ← i =j k m (f f f f f f f j , f f f f f f f i )q i (x ′ ) (Optimized in Section 4.3.2) end end for x j ∈ Ω do /* Compatibility transform */ qj (x j ) ← x ′ ∈Ω (1 -δ x ′ xj ) 2 m=1 w m qj,m (x ′ ) /* Local update */ q j (x j ) ← exp (-ψ u (x j ) -qj (x j )) end Normalize q j end end 106 4.

Markov Chain Variational Inference for fcCRFs

We first start by rewriting the message passing step in terms of a convolution. For m = 1 we have:

qj,1 (x ′ ) = i =j k 1 (f f f f f f f j , f f f f f f f i )q i (x ′ ), = i∈S k 1 (f f f f f f f j , f f f f f f f i )q i (x ′ ) -q j (x ′ ), = [ Ġθα,θ β * q(x ′ )](f f f f f f f j ) -q j (x ′ ), (4.8)
where * is the convolution operator, Ġθα,θ β is a Gaussian kernel in the augmented space and q is the q distribution in the augmented space. The notion of augmented space is needed to interpret the appearance kernel as a convolution. Details about the augmented space are given in Appendix E.1. For m = 2 we have:

qj,2 (x ′ ) = i =j k 2 (f f f f f f f j , f f f f f f f i )q i (x ′ ), = i∈S k 2 (f f f f f f f j , f f f f f f f i )q i (x ′ ) -q j (x ′ ), = [G θγ * q(x ′ )](f f f f f f f j ) -q j (x ′ ), (4.9)
where G θγ is a Gaussian kernel of parameter θ γ . In this case of the smoothing kernel, the interpretation in terms of a convolution is straightforward. For tractability of the MF VI algorithm, the methodology of (Paris and F. [START_REF] Olariu | Modified variational Bayes EM estimation of hidden Markov tree model of cell lineages[END_REF] is used to approximate the convolutions. The observation is that a convolution acts as a low-pass filter, hence its result can be approximated by performing a convolution on a downsampled version of the filtered object. The result is then upsampled again. In practice, the downsampling step is taken equal to the kernel standard-deviation and the downsampling is done via averaging.

Remark: A further optimization consists in performing the multidimensional Gaussian filtering independently on each axis due to the separability of the Gaussian kernel. Note also that the parameters θ α and θ γ can be vary on each axis such that we have θ θ θ θ θ θ θ α = (θ α,x , θ α,y , θ α,z ) and θ θ θ θ θ θ θ γ = (θ γ,x , θ γ,y , θ γ,z ).

Markov Chain Variational Inference for fcCRFs

Chapter 3 illustrated the dramatic increase in performances in structured VI procedures over MF VI procedures. That is why in this section we propose to study the gain in the classical MF VI for fcCRFs with a Markov Chain (MC) structured VI inference approach. This is an original approach that we propose; to the best of our knowledge there is no study in the literature that deals with structured VI in fcCRFs.

Markov Chains for image processing

While MCs are intrinsically unidimensional data structures, they have been numerous examples of their use to process data of higher dimension [START_REF] Bricq | Unifying framework for multimodal brain MRI segmentation based on Hidden Markov Chains[END_REF][START_REF] Yahiaoui | Markov Chains for unsupervised segmentation of degraded NIR iris images for person recognition[END_REF] (Giordana and Pieczynski, 1997) (Lanchantin et al., 2011) [START_REF] Fjortoft | Unsupervised classification of radar images using hidden Markov chains and hidden Markov random fields[END_REF]. The advantage of the MC structure resides in the small computational cost and the exact inference procedures that are available. The key step is to scan the higher dimensional space with a well designed path. The most common approach is to scan the space with space filling curves such as the Hilbert-Peano scan [START_REF] Sagan | Space-filling curves[END_REF]. This is used in a 2D context in (Giordana and Pieczynski, 1997) and in a 3D context in [START_REF] Bricq | Unifying framework for multimodal brain MRI segmentation based on Hidden Markov Chains[END_REF]. However the path could be adapted, according to the specific application, for a better statistical segmentation [START_REF] Courbot | 2-step robust vertebra segmentation[END_REF] (Meriem [START_REF] Yahiaoui | Markov Chains for unsupervised segmentation of degraded NIR iris images for person recognition[END_REF].

Scanning the data with Markov Chains

We developed the MC VI approach for inference in fcCRFs to answer the strong and limiting assumption of fully independent random variables in the variational distribution of the MF approach. Indeed, with the goal in mind to refine 3D segmentations produced by a simple stacking of 2D slices segmented by the CNN, it seems important to take as much spatial context as possible in the post-processing step. We introduced the concept of a MC scan in Section 4.4.1.

In this section, we define another way to use MCs in the processing of higher dimensional data which is based on parallel and independent MCs. Indeed, for computational reasons, it is prohibitive to perform a Peano scan such as defined in [START_REF] Bricq | Unifying framework for multimodal brain MRI segmentation based on Hidden Markov Chains[END_REF] for 3D volumes inside a VI procedure.

The variational distribution that we propose is made of independent MCs that goes through the data cube following one of the cube axes. This follows an idea also found in the Factorial Hidden Markov Models introduced in (Z. [START_REF] Ghahramani | Factorial hidden Markov models[END_REF]. There is, however, no reason to privilege an axis over another, we then propose to perform parallel MC VIs along each axis and then to merge their result. More precisely, for a 3D image with dimensions M × N × K, respectively the height, width and depth. Then, we propose to perform six parallel MC VIs, one in each direction of each axis. This is illustrated in Figure 4.3.

For a given MC VI procedure, e.g., in the cases top-down or down-top (N ×K MCs of length M ), the variational distribution is (recall the definition of a MC given in Equation 1.20):

l(x x x x x x x) = N ×K n=1 l n 1 (x n 1 ) M m=2 l n m (x n m |x n m-1 ). (4.10)
The quantity to maximize, -KL(l(x x x x x x x)||p(x x x x x x x|y y y y y y y)), then becomes in this section:

-KL(l(x x x x x x x)||p(x x x x x x x|y y y y y y y)) = E x The details of the derivation of the update equations for the MC VI approach are given in Appendix E.2. Once the update equations are found we use the same approach as shown in Section 4.3.2, that is, we interpret the bottleneck equation (summation over all the sites) as a convolution. Algorithm 4.2 summarizes the MC VI for fcCRFs. The update of all the terms of the variational distribution in Algorithm 4.2 can be done in a parallel fashion; this is the foundation of the practical MC VI methodology we present in Section 4.4.2.

Performing the MC VI procedures yields six posterior marginal probability distributions. Let these distributions be ∀s ∈ S, ∀x s ∈ Ω, l bf s (x s ), l f b s (x s ), l td s (x s ), l dt s (x s ), l lr s (x s ) and l rl s (x s ), associated, respectively to the front-back, back-front, top-down, down-top, left-right and right-left MC VI procedure. Then the final posterior marginal probability distribution results from taking the product of each MC VI, at each site and for each class. That is, the final marginal probability distribution is defined by, ∀s ∈ S, ∀x s ∈ Ω:

ls (x s ) = l bf s (x s )l f b s (x s )l td s (x s )l dt s (x s )l lr s (x s )l rl s (x s ). (4.12)
From l we perform the final segmentation following the MPM criterion such that, ∀s ∈ S:

xMCV I s = argmax x∈Ω ls (x s ). (4.13)
Remark: Other merging rules were tested such as a majority voting for each pixel based on the segmentations provided by each MC VI, or taking the maximal probability between each MC VI, at each site and for each class. Those criteria performed worse than Equation 4.12 and their results will not be reported.

Remark:

The approach we propose is straightforwardly transposed in the context of 2D images by performing four parallel MC VIs (top-down, down-Algorithm 4.2: MC VI for fcCRFs (example of top-down or down-top cases)

Data: w 1 , w 2 , θ α , θ β , θ γ , the model parameters, y y y y y y y, the observations. Result: l(x x x x x x x), the MC approximation of p(x x x x x x x|y y y y y y y).

/* Initialize l */ l n m (x n m |x n m-1 ) ← exp (-ψ u (x n m )) , ∀n ∈ {1, . . . , N × K}, ∀m ∈ {1, . . . , M }, ∀(x n m , x n m-1 ) ∈ Ω 2 Forward Backward (Algorithm A.2) to get the marginals of l while convergence is not attained do for n ∈ {1, . . . , N × K} do for m ∈ {1, . . . , M } do for x ′ ∈ Ω do for r ∈ {1, 2} do ln m,r (x ′ ) ← (m ′ ,n ′ ) =(m,n) k r (f f f f f f f n m , f f f f f f f n ′ m ′ ) x ′′ ∈Ω l n ′ m ′ -1 (x ′ |x ′′ )l n ′ m ′ -1 (x ′′ ) (Optimized in Section 4.3.2 by recognizing a convolution) end end for (x n m , x n m-1 ) ∈ Ω 2 do ln m (x n m |x n m-1 ) ← x ′ ∈Ω (1 -δ x ′ x n m ) 2 r=1 w r ln m,r (x ′ )+ x n m-1 ∈Ω l n m-1 (x n m-1 |x n m-1 )l n m-2 (x n m-2 )× log l n m-1 (x n m-1 |x n m-2 ) l n m (x n m |x n m-1 ) ← exp -ψ u (x n m ) -ln m (x n m |x n m-1 ) end Normalize l n m end end Forward Backward (Algorithm A.
2) to get the marginals of l end

Experiments and Results

top, left-right and right-left). The final posterior marginal distribution is also computed by adapting Equation 4.12 to a product of four terms.

Experiments and Results

Segmentation via Deep Learning

We here describe the experiment of histological image segmentation and its results with the CNN described in Section 4.2.

Experimental set-up

The medical context surrounding this application is described in Chapter 5, the focus in this section is on the mathematical aspects of the problem.

The annotated dataset is made of 91 original slices of the type shown in Figure 4.1. We consider a 6-class segmentation problem; the mCT image pixels are then assigned to one of the 6 labels: soft tissue, fatty tissue, sheet calcification, nodular calcification, specimen holder and background. The class acronyms are respectively ST, FT, SC, NC, SH, Ba. Among the 91 annotated slices, 10 are held apart from the beginning of the experiment, they form the test set. The remaining 81 original annotated slices form the training set. Then the train and test sets undergo data augmentation algorithms (rescaling, rotation and gamma correction) [START_REF] Schmitt | A survey on image data augmentation for deep learning[END_REF]. This raises the train and test set to respectively 1620 and 200 elements. Eventually the train set is split into the actual train set and a validation set whose size is set to 10% of the former training set. The network we use is the same as in Figure 4.2 but with inputs of size 512 × 512 (images are cropped to fit this size). The network is then trained with respect to the mean Dice loss function. The experiments are coded in Python3 using the Tensorflow4 library.

Results

Figure 4.4 illustrates the evolution of the mean Dice loss on the train and validation set during training. We can see that, after 50 iterations, the loss is stabilized and the training stops. Table 4.1 gives the Dice score for each class on the test set. We can see some important discrepancies in the results between the classes. These discrepancies reflect the intrinsic complexity of the dataset. Indeed, the latter presents high class imbalances, high similarities between classes and notable differences between the image dynamic ranges. Note however that one key goal was to distinguish, on mCT images, the two types of calcifications namely, SC and NC. This is an important histopathological aspect. While such a task seems extremely complex for the naked eye, our DL approach can, to some extent, answer this problem. Future research might further improve this point. An example of a 2D prediction made by the CNN on a slice from the test set is given in More results and in-depth discussions about this semantic segmentation problem and some histopathological interpretations are available in Chapter 5.

Experimental comparisons of MF VI and MC VI on semi-real images

Before applying the fcCRF VI procedures on real data, we evaluate the relative performance of MF VI and MC VI on synthetic data. We consider the case of unsupervised grayscale two-class image segmentation corrupted by additive Gaussian noise. The set of hidden classes is Ω = {ω 1 , ω 2 }. The additive Gaussian noise is parametrized by a mean and a standard deviation for each class: (µ ω , σ ω ), ∀ω ∈ Ω. A first segmentation is performed using a Gaussian Mixture Model (GMM) [START_REF] Mourad | A survey on latent tree models and applications[END_REF]. The posterior probabilities of the GMM are used for the initialization of the unary potentials of the fcCRF model, for each pixel and for each class. We define as err GM M , err M F V I-P e and err M CV I-P a the error rates in the segmentation with, respectively, the GMM, the fcCRF model with MF VI, the fcCRF model with MC VI with Peano scan and the fcCRF model MC VI with parallel MCs. In all cases we fix the number of iterations in the VI procedures to 50. In the following, we provide experiments offering a general overview of the behavior of the models.

In Figure 4.6 we plot the error rates of the approaches for the unsupervised segmentation of images. First we emphasize the fact the MC VI approach appears as more efficient for inference in fcCRFs than the classical MF VI approach both with the parallel MCs and the Peano scan implementations. In particular, in the case of error rate over MF VI, and provides a better or equivalent error rate for all noise levels. We also note that the GMM approach is not competitive against the three others, its main drawback lies in its inability to take into account any spatial context. The fcCRF approach appears as an efficient post-processing technique which is built upon the probability map given by the GMM method.

The fcCRF approach significantly improves the GMM segmentation with more than 10 point gains in segmentation error. Figure 4.6 provides evidence of the similar performances of the parallel MCs and Peano scan implementations of the MC VI procedures in terms of error rate. However they are not equivalent in terms of computational cost. Indeed, in Figure 4.7 we plot the time taken for each procedures as a function of the size of the image5 . The Peano scan approach suffers mainly from the fact that it cannot be parallelized and must process very long MCs. This experiment discards the use of the Peano scan in favour of our proposed parallel MCs scan.

In Figure 4.8 we plot the surface whose points are defined by err M F V Ierr M CV I-P a where the varying parameters are the kernel weights w 1 and w 26 . Except for the pathologic case w 1 = 0, the surface error is always above 0, suggesting that the MC VI (with parallel MCs) approach that we propose always gives better results that the classical MF VI approach for inference in fcCRFs.

Finally, in Figure 4.9 we illustrate graphically the results on particular images. The images come from the "1070-Binary Shape Database"7 .

We conclude this section by noting the evidence of the advantage in favour of the MC VI procedures (both with Peano scan and parallel MCs) for fcCRFs in the case of image segmentation. While we did not work on finding the op-

µ ω1 = 0, µ ω2 = 1, θ α = θ β = θ γ = 1.
The new MC VI methodologies appear to be the best performing methods for all noise levels. The results are here averaged on 10 randomly chosen images from the dataset. The results are averaged over 10 simulations.

timal set of parameters for the models, we performed several experiments that exhibited better results for the MC VI for any non-pathological parameter set.

Similarly to the conclusion of Chapter 3, the structured VI approach seems to better take into account the spatial context thanks to the MCs, leading to better segmentation rates.

Remark: In the following, all references to MC VI will refer to the parallel MCs version of MC VI.

Post-processing with fcCRFs

The rest of this chapter is dedicated to improving the CNN segmentation with a post-processing method based on fcCRFs. We here compare the results in the segmentations provided by the CNN (Section 4.2.2), by a fcCRF postprocessing based on MF VI (Section 4.3.2) and by a fcCRF post-processing based on the new MC VI (Section 4.4).

MF VI with fcCRFs for segmented 3D volumes

In this section we employ the fcCRFs with MF VI to segment 3D data: the whole artery can be reconstructed in 3D. Figure 4.10 illustrates the resulting segmentation on a mCT from a test artery. The post-processing enables taking into account the third axis to reconstruct a volume with spatial context from all the dimensions. More segmentation are provided in Chapter 5.

The soft tissue class misclassifications

We focus on the problem of correcting spurious classifications from the CNN which erratically misclassifies the soft tissue class into the background class leading to an inhomogeneous segmented volume. This problem is similar to [START_REF] Kamnitsas | Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation[END_REF], which also uses the fcCRFs to a similar end. The problem is described in Figure 4.11. (Kamnitsas et al., 2017) uses fcCRFs and MF VI to process the CNN output. They rely on the smoothing capabilities and the modeling of the 3D context provided by the fcCRFs to solve the misclassifications. Illustration of such process is done in Figure 4.12a. However, MF VI does not seem suitable in our case: it results in the loss of the small calcification pixels which are essential in our application. Those unsatisfactory results made us consider a new methodology based on MC VI. It is described in the next section.

Using MC VI with fcCRFs to recover the missing soft tissue class

In this section we describe the processing based on MC VI. MC VI offers the additional possibility to refine the segmentation with known constraints from the application. Indeed, the VI procedure produces a non-stationary transition matrix between the states in Ω. One then has an understable view of the model that has been learnt before the final inference in the variational distribution with the Forward Backward algorithm. Since a MC transition matrix is easily comprehensive, we also propose to adjust the probabilities that are learnt in order to correct learnt transitions that does not fit the desired application. We want to penalize the predictions made in favour of the background class in the areas were the soft tissue class is missing. This can be done thanks to the transition matrix of the MCs that are learnt during the MC VIs processes. A similar idea of matrix adjustment to exhibit class relationships to improve classifications predictions can be found in [START_REF] Fidon | Generalised wasserstein dice score for imbalanced multi-class segmentation using holistic convolutional networks[END_REF] for example.

We start by detecting the areas where the soft tissue is missing by using a binary closing operator. Then, within these areas denoted A, at each iteration of the MC VI, the probability of transitioning from any class to the background class is lowered so that, for the transition matrix of any MC VI, with ω bck being the background class:

∀ω ∈ Ω \ {ω bck }, ∀s ∈ A, l s (ω bck |ω) = min ω ′ ∈Ω l s (ω ′ |ω). (4.14)
Such an adjustment improved greatly the segmentations enabling the soft tissue class to be recovered thanks to the spatial context brought by each VI. An example of MC VI is given in Figure 4.12b (single slice illustrated) and in Figure 4.13c (3D segmented volume). We note how the MC VI methodology leads to the most refined results. It is also computationally efficient since both the MC VI procedures and the processing of each MC inside each of the procedure are parallelized. The computational bottleneck of MC VI remains the same as in MF VI, namely, the bilateral filtering step.

Remark: Such an adjustment in the MC VI procedure does not possess an equivalent adjustment within the MF VI procedure. The closest experiment : KL values (up to a constant) for the VI procedures as a function of the iterations regarding the experiment from Figure 4.12. We can see that all the procedures converge in a few iterations.

that was carried consisted in a direct modification of the estimated marginal probabilities q s (x s ), ∀s ∈ S, in the same way as done in Equation 4.14, which resulted in bad performances.

Remark: In all this section, as suggested by [START_REF] Kamnitsas | Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation[END_REF], the parameters were manually tuned. Further developments might consider an optimization procedure for automatic fcCRF parameter selection.

Convergence of the VI procedures

Finally we give an insight about the convergence of all the VI approach. This is in done in Figure 4.14. We can see that in all cases the convergence is very quick to happen. As discussed in the supplementary material of (Krähenbühl and Koltun, 2011), the KL values (Equation 4.11) are computed up to a constant because the partition function of p(x x x x x x x|y y y y y y y) cannot be computed. The partition function can be discarded since it appears as an additive term which only depends on the observations y y y y y y y and the CRF parameters which are constant accross the VIs.

Conclusion

In this chapter we treated an unique dataset of mCT of atherosclerotic arteries, whose histological ground truths have been annotated by an expert during the course of the thesis. That is why we studied supervised segmentation approaches and more precisely, deep learning approaches. We also studied the complementary of deep learning and probabilistic approaches.

We first processed the data and trained a CNN to segment any new mCT data of atherosclerotic data into six histologically meaningful classes. We then 4.6. Conclusion studied the problem of using the CNN outputs to reconstruct a 3D segmented volume. We used the fcCRF probabilistic model with both its classical MF VI procedure and a new structured MC VI procedure. The latter seems to enable a much better use of the spatial context than the former, leading to better results, both on synthetic and real data. Thus, we brought new solutions to the problem of volumetric semantic segmentation introduced at the beginning of the chapter.

More medical context, interpretations and results, as well as illustrations of this study are available in Chapter 5.

Introduction

The previous chapters were dedicated to the mathematical presentations of new probabilistic models motivated by new challenges in the field of image processing for medical imaging. In this chapter, the focus is on the application of some of the models discussed previously on real problems and real data collected in the field of vascular surgery.

The first section describes the development of a Hidden Markov Chain model for the precise segmentation of stents corrupted by artifacts in X-ray images. We then present additional results of the applications of the Gaussian Pairwise Random Field model (Chapter 2) and of the Spatial Triplet Markov Tree model (Chapter 3) for the segmentation of organic material and calcifications in Xray images. The third section completes the presentation and evaluation of the U-Net Convolutional Neural Network (Chapter 4) with more background on the histologic problem, on the dataset construction as well as additional illustrations.

Unsupervised segmentation of stents corrupted by artifacts

Context and motivation

In this section, we describe a statistical approach based on the Hidden Markov Chain (HMC) model, introduced in Section 1.6.1, to finely segment stents in medical images corrupted by strong artifacts.

Medical context

The in vivo behavior of the stent is far from being fully understood [START_REF] Lejay | Explanted vascular and endovascular graft analysis: where do we stand and what should we do?[END_REF]. However, the amount of medical images is constantly growing, offering new opportunities to learn from clinical cases. The interactions between the stent struts and the calcifications are believed to be a major cause of failure of the treatments based on the implantation of a stent. The study of such interactions by combining information of both CT and mCT scans is a new approach to the problem, it requires the development of new imaging tools which can finely segment the metallic stent components.

The most complex images to process often include a broken stent in a calcified environment. Images are also often corrupted by strong artifacts. Figure 5.1 illustrates the typical images from which one needs to segment the stent. It appears that the task of unsupervised segmentation in X-ray scans is a complex problem. Notably, artifact and calcification pixels are close to the stent pixels both in intensity, geometry and localization, which tends to produce stent False-Positive classifications with traditional segmentation methods as we will see in Section 5.2.3. Some stent, calcification and artifact components are indicated by, respectively, red, blue and green arrows. We note the spatial resolution differences between CT scans and mCT scans.

Unsupervised segmentation of stents corrupted by artifacts

Related literature

To the best of our knowledge, there exists no dedicated method to confidently segment stents in such images in an unsupervised fashion1 . As a consequence, people still use simple approaches such as manual thresholding which fails in complex cases. Thus, such cases are discarded from clinical datasets despite their medical interest [START_REF] Perrin | Patient-specific simulation of endovascular repair surgery with tortuous aneurysms requiring flexible stent-grafts[END_REF].

Relatively few works deal directly with the topic. Works such as [START_REF] Klein | Automatic segmentation of the wire frame of stent grafts from CT data[END_REF] or [START_REF] Langs | Learning deformation and structure simultaneously: In situ endograft deformation analysis[END_REF] that also address stent segmentation are not suitable for us since they do not consider the stent in a calcified environment. In our case, the calcifications and the stent are so close in appearance and sometimes in geometry that dedicated methods need to be developed to distinguish both classes during segmentation. Moreover, we have to discard approaches such as [START_REF] Demirci | 3d stent recovery from one x-ray projection[END_REF] since they are stent specific and we need our method to be independent on the stent model to process our data.

However, the problem of stent segmentation is somehow very close to the problem of blood vessel segmentation (see [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF] for a review on the topic which concerns traditional approaches). Notably, some multi-scale filters have met success in the task of blood vessel segmentation [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] [START_REF] Merveille | Tubular Structure Filtering by Ranking Orientation Responses of Path Operators[END_REF], these methods will be used for comparisons in the result section. During the last years, Deep Learning approaches are performing

A HMC dedicated to handling strong artifacts

The HMC path We already mentioned in Section 4.4.1 the importance of the transformation of higher dimensional data into one dimensional data, so that the latter can be processed by a MC.

In this section we propose a new kind of path dedicated to handling the artifacts. Each connected component (of the stent class from the initial segmentation) of each slice is visited using a snail path (snail paths are described in (M. [START_REF] Yahiaoui | Implementation of unsupervised statistical methods for low-quality iris segmentation[END_REF]). The sense of rotation is alternated between each snail path. This way, physical continuity in the path is simulated, gathering all the stent parts in a 1D sequence which allows efficient computation. Such a path is particularly suited to get a smoothed and regular surface for tubular structures [START_REF] Courbot | Vertebra segmentation based on two-step refinement[END_REF]. Moreover, in the HMC model one has to learn the transition probabilities between the classes. A small morphologic dilation is then performed over the result of the initial segmentation (see the end of this section) for the snail path to be a little bigger than the connected components isolated during the initial segmentation. Thus, we make sure that 5.2. Unsupervised segmentation of stents corrupted by artifacts all the voxels belonging to the stent class are part of the sequence of pixels to classify.

Figure 5.2 illustrates our choice for the path of the MC through the slices of our data cubes.

The HMC statistical model

The segmentation problem presented in this section aims at classifying pixels into two classes, stent and rest, in an unsupervised fashion. Recall the context of Bayesian image segmentation introduced in Section 1.5.2. Let X X X X X X X be the vector of hidden random variables of the class image taking their value in If the hypothesis of independent Gaussian noise is made then the parameters that we need to learn are the initial probabilities of the MC, the transition probabilities of the MC, the mean and the standard deviation of the Gaussian distribution associated with each class. The parameters are learnt unsupervisedly with the SEM algorithm (Algorithm A.5) similarly to the parameter estimation steps described in Chapters 2 and 3. The statistical inference step is performed by the FB algorithm (Algorithm A.2) followed by the MPM algorithm (Algorithm A.6).

Ω = {stent, rest}. Let Y Y Y Y Y Y Y

Improving the noise model to handle strong artifacts

Classically, such as in [START_REF] Courbot | Vertebra segmentation based on two-step refinement[END_REF], the noise model in a HMC, also called conditional likelihood, is chosen to be a mixture of Gaussian distributions (MoG). Using developments on generalized mixture models in Hidden Markov Models such as [START_REF] Delignon | Estimation of generalized mixtures and its application in image segmentation[END_REF] and [START_REF] Pieczynski | Estimation of generalized mixture in the case of correlated sensors[END_REF], we propose to work with a noise model involving mixture of exponential distributions (MoE) to improve the results. MoEs have already been explored for Hidden Markov Trees in [START_REF] Monfrini | Estimation de mélanges généralisés dans les arbres de Markov cachés, application à la segmentation des images de cartons d'orgue de barbarie[END_REF]. Indeed, the "all-or-nothing" behavior enabled by using the exponential distribution seemed particularly suited to handle the strongest artifacts. The exponential probability density function is given by:

f (x; λ, δ) = λ exp (-λ (x -δ)) x > δ, 0 x ≤ δ, (5.1)
where (λ, δ) ∈ (R + * , R) are the new parameters to estimate for each class in the SEM procedure. Moreover, a Bayesian Information Criterion (BIC) [START_REF] Wit | All models are wrong...': an introduction to model uncertainty[END_REF] score comparison was conducted and we showed that a MoE fitted better the empirical distribution (obtained after the coarse segmentation) for the stent class and the rest class than a MoG. It is described in Appendix F.1.

The complete segmentation approach

The HMC needs to be initialized by an initial segmentation, we propose that its initialization correspond to the response of a Frangi filter [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] followed by a region growing algorithm and a watershed algorithm (Haidekker, 2011). This forms the Preprocessing part of Figure 5.3. The Frangi filter is first performed and we keep only voxels whose probability to belong to a tubular structure is above a certain threshold (dependent on the nature of the image we want to segment). A region-growing algorithm is performed in 3D to select only the biggest connected component (the pixels from the stent class). A watershed algorithm is performed on every slice to separate connected components that should be disjoint but still appear stuck together because of the artifacts. Remark: A step of upsampling is needed for images of very low resolution such as CT scans (see, e.g, Figure 5.1b). The idea is that the stent parts on the image should be increased up to a certain width (we used 10 pixels in our applications) so that the successive filters used for segmentation show responses.

Remark: The Frangi filter is the most important part of the Preprocessing algorithm. The subsequent HMC algorithms can be seen as a statistical refinement of the Frangi filter response, increasing the robustness of the latter and making it fully automatic.

Results

We begin by illustrating segmentations by our fine segmentation method. We start by considering the refinement introduced by using a noise model involving a mixture of exponential distributions, which better reflects the nature of the artifacts in the images. Indeed, as depicted in Figure 5.4 on several cases with very strong artifacts, the recovered stent structure is much thinner and much more circular as the original metallic structure of the stent. Figure 5.5 illustrates some more segmented 3D meshes of broken stents.

As stated in the introduction, we can apply our new metal segmentation method to improve Metal Artifact Reduction (MAR) procedures. Beam Hardening Correction (BHC) is built upon the Linear Interpolation technique and they form the most classical approaches to MAR (Verburg and Seco, 2012) (Y. [START_REF] Lianwen | Convolutional neural network based metal artifact reduction in x-ray computed tomography[END_REF]. We show, with Figure 5.6, on various CT scans, restored images with reduced artifacts using our fine statistical segmentation approach followed by BHC. Results are compared to a manual segmentation, to a seg-

Unsupervised segmentation of stents corrupted by artifacts

mentation based on Frangi filter2 [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] and to a segmentation based on the RORPO filter 2 [START_REF] Merveille | Tubular Structure Filtering by Ranking Orientation Responses of Path Operators[END_REF], all three followed by BHC. Note that all methods operate in 3D but 2D slices are presented for a better visual assessment of the results. Such results highlight the importance of the metallic segmentation step in MAR problems.

The absence of ground truth forces us to make a qualitative assessment of the segmentations. Our algorithm is the most likely to offer a precise refinement over the stent pixels as well as avoiding False-Positive classifications resulting from calcifications or artifacts. The precise segmentation that we propose reduces the remaining artifacts at the end of the MAR procedure.

Remark: It is notable that the results from the three methods used for comparisons depend on a final step of manual thresholding. On the contrary, the proposed HMC does not need such manual intervention since the statistical refinement is fully automatic, avoiding a time consuming step for the expert which is also subject to the operator subjectivity. In this section we provide some more illustrations on real world images of segmentations using the GPMF model. We discuss the current performances of the model and explore its current limits. Recall the context of Section 2.5.4, the goal is to segment the organic biomaterial in mCTs of stented arteries corrupted by spatially correlated noise. In this context, Figures 5.7 A closer examination of the data seems to reveal that a bad performance of the GPMF is mainly caused by a smoothing effect that entails a quite substancial loss of details in the segmentation. The stationarity assumption might explain the problem, indeed, the model might reflect poorly the reality in zones where there is no or very few spatial correlation.

We might argue that the most interesting zones are the ones directly around the stents. Indeed, an important biomechanical question that initiated this work is whether the stent is in contact with calcifications. Such kind of interactions are believed to be of crucial importance in many cases of the failure of the treatment such as stent fracture. Hopefully, at such locations of the images the GPMF model is performing well in all cases: see the FN/FP rates computed in Table 5.1.

Remark: It also seems that the most complex cases for the GPMF model are the ones where the artifacts are the strongest, or more spread on the image. However this last remark has not currently been explored since the literature seems to lack a good measure of the metallic artifacts on an image.

Unsupervised calcification segmentation with STMTs

In this section, we work in the context of unsupervised image segmentation with STMTs. We want to compare, on real images, their relative performance with the HMF and HMT models. The goal of the experiment is to segment the atherosclerotic calcifications on a mCT image: the brightest parts of the mCT images (see Figure 5.9). We also add additive independent Gaussian noise to form more diverse observations Y Y Y Y Y Y Y = y y y y y y y. Such a construction reflects situations where the mCT image would be of less good quality. Our objective is to evaluate, on real images, the capacity of STMTs to take the spatial context into account spatial context while offering fully deterministic inference computations. That is why the segmentation of atherosclerotic calcifications is the chosen objective: they are very homogeneous and rounded shapes which is expected to favor modelizations which can take into account the spatial context.

In this experiment we compare the HMF-IN model defined in Equations 3.29 and 3.30, the (quadtree) HMT-IN model, defined in Equation 1.24, and the (quadtree) STMT with Independent Noise, defined in Equations 3.5. The parameters are estimated unsupervisedly, with the IPET algorithm (Algorithm 3.2), as described in Section 2.3 for HMFs-IN and Section 3.3.2 for STMTs and HMTs-IN. Figures 5.10,5.11 and 5.12 depict such an experiment. The HMF-IN model seems to better segment the proposed shapes. However, for the higher noise levels, as stated in Chapter 3, we also encounter the issues arising from the stochasticity of the HMF-IN inference and parameter estimation procedures. Indeed, the results can differ a lot between two runs of the same experiment.

We also see that STMTs seem more prone to capture the spatial correlations than the HMT-IN model. As a result they offer better segmentations than the HMT-IN model as the noise level increases. At the lowest noise levels, STMTs and HMTs-IN perform similarly. When slightly better error rates are in favor of the HMT-IN model versus the STMT model, as in the first column of Figure 5.12, we find a contradiction with the results on synthetic data of Chapter 3. We might argue that this is due to the real world data we treat: the noise model might not be the best modelization for the image. The LLS estimation of the parameters might also be a limiting factor.

Histological segmentation of microCT with Deep Learning

In this section we provide more context, details and results to the original project presented in Chapter 4.

Motivation

Our goal is to study peripheral obstructive artery diseases in femoropopliteal arteries. They are common affections of the blood vessels which are still not well understood [START_REF] Yahagi | Pathophysiology of native coronary, vein graft, and in-stent atherosclerosis[END_REF]) (Torii et al., 2019). On the one hand, microCT has been proven to be a valuable tool for the recognition of the different atherosclerotic components [START_REF] Jinnouchi | Micro-Computed Tomography Demonstration of Multiple Plaque Ruptures in a Single Individual Presenting With Sudden Cardiac Death[END_REF]. On the other hand, histology, which is the gold standard to analyze and study the atherosclerotic process is a destructive technique and requires a skilled expert. Therefore, in this project we aim for an automatic histologic segmentation of microCT images based on artificial intelligence and, more precisely, Deep Learning (DL) which has become the state of the art approach for biomedical semantic image segmentation. With the algorithm we developed, the classic histologic process might be replaced by an automatic and non-invasive segmentation for a first analysis. Such an algorithm will save time and enable to process more data in a shorter time, since the histopathologists will be guided towards the most interesting biological materials. Eventually, the algorithm would help develop knowledge of the atherosclerotic process that would lead to better and more personalized treatments.

Protocol and database construction

In this section we give the main steps of the protocol that we set up for the construction of the database that was used as learning data for the DL algorithm.

The steps are illustrated with images in Figure 5.14. In all the descriptions of this work, the colormap used for the histologic segmentation is given in Figure 5.15. The steps of the proposed protocol are:

1. Retrieval of femoropopliteal arteries from amputated legs from patients who had undergone transfemoral amputations3 . These explanted arteries were collected thanks to the Geprovas collaborative retrieval program. Figure 5.13 shows explanted arteries from the study.

2. The microCT 3D images of the arteries were acquired at the CVPath Institute, Inc, (Gaithersburg, MD, USA) using a Nikon X-Tek XT H 225ST4 .

3. Histology was then performed on the specimens as described in (Torii et al., 2019). Some histologic slices are depicted in Figure 5.14a.

4. Co-registration was subsequently performed manually between the mi-croCT images and the histologic slices obtained during the two steps described above. This was performed with VGSTUDIOMAX 3.05 and ImageJ6 . The result of this step consists in pairs of data: the microCT 2D image with its histologic ground truth. Some co-registered microCT images are illustrated in Figure 5.14b.

5. An expert then annotated the microCT images using the histologic ground truths in the GIMP software 7 . It was decided to segment the data into 11 classes. It has then been decided with the histopathologists that 6 classes were of interest to develop a first version of the algorithm. We now list the 6 classes, with the subclasses that possibly compose them:

• soft tissue (ST ): soft tissue, formaldehyde, thrombus, fibrous plaque.

• fatty tissue (FT ): fatty tissue, lipid pool.

• sheet calcification (SC ).

• nodular calcification (NC ).

• specimen holder (SH ).

• background (Ba): background, true lumen. In the list above, we also specify the acronym that might be used to refer to the class. This step results in the creation of data pairs (microCT image and the corresponding class image) that are the data the DL algorithm will be trained on. Some annotated slices are illustrated in Figure 5.14c.

The subsequent steps of the project are described in Chapter 4. Those steps deal with the set-up and training of a Convolutional Neural Network (CNN) to achieve the automatic segmentation goal we described in Section 5.4.1. We now focus on the analysis of the resulting segmentations provided by the CNN.

Remark:

The choice to reduce the number of classes from 11 to 6 has been made for various reasons. It is essentially a compromise between providing a segmentation meaningful from a histologic viewpoint but also to set up a solvable problem from a DL viewpoint. Indeed, this project represents exploratory work where data are rare, highly imbalanced and the some classes exhibit high similarity in images.

Histological segmentation of microCT with Deep Learning

Results

2D segmentation on the test set Table 4.1 references the Dice scores on the test set. These results indicate several degrees of difficulties in the automatic recognition of the classes. First of all, the soft tissue, background and specimen holder classes are well segmented (Dice scores > 0.86). Secondly, the CNN is capable of responding quite well to a crucial histologic analysis: the differenciation of the calcification classes, sheet calcification and nodular calcification, whose Dice scores are, respectively, 0.85 and 0.64. This ability is remarkable since the differenciation is almost impossible by the naked eye, even for a trained expert. Finally, we note that the fatty tissue class is the most complex to segment (Dice score of 0.41). This can be explained by the fact that this class is really similar in appearance to the soft tissue class. It is also under represented in the training data, which might be another cause for the high number of misclassifications.

To go further in the analysis, we give in Table 5.2 the multiclass confusion matrix associated with the segmentations of the CNN on the test set. We can then be more precise about the comments we made on the Dice scores. Looking at the SC and NC rows of the table we see that the main source of errors in the segmentation of the calcifications comes from a confusion between the two calcification types. Interestingly, by reading the FT row, we get that the reason for the bad FT score is that this class is most of the time misclassified as the ST class.

By considering the segmentation of each class as a binary One versus the rest classification problem we can compare the segmentation of each class using the Precision Recall (PR) curve and the Area Under Curve (AUC). A bigger AUC reveals a better classifier for a given class. PR curves and AUCs are given in Figure 5.16. This analysis exhibits the same hierarchy between classifiers as the Dice scores discussed above. Note that because of the high class imbalance the Receiver Operating Characteristic (ROC) curves are not exploitable, they are biased by the huge number of negatives from the Ba class and cannot be properly compared. See (Murphy, 2012) for more details on PR and ROC curves.

To conclude this section on the 2D segmentations, Figure 5.17 depicts some segmented slices by the CNN that illustrate the detailed discussion we give above. The illustrations reflect the complexity to correctly classify pixels from the FT class and the confusion that sometimes arises when classifying pixels into one of the calcification classes.

3D segmentation on test femoropopliteal arteries

We now provide in Figure 5.18, some examples of the complete segmentation of an artery by the CNN after fcCRF and MF VI post-processing (see Chapter 4). Such results are directly understandable by the histopathologists who can then analyze the explant at a histologic level with a microCT acquisition only. The expert can then decide how to go further in their analysis.

Chapter 5. Applications to vascular surgery

Conclusion

In this chapter we offered an applicative perspective on the probabilistic models described along the thesis. Our goal was to show to what extent they could be the foundation of new tools used in the field of vascular surgery to develop and enhance the medical research as well as the clinical workflow. The results in this chapter are founded on CT, mCT and histological images. This highlights the importance of the imagery and the processing of the produced images in vascular surgery. The key role of images is expected to grow with time.

Conclusions and Perspectives

Gaussian Pairwise Markov Fields

Family of the probabilistic model studied: Hidden Markov Fields Model, innovation: Gaussian Pairwise Markov Fields Issue that has been addressed: Unsupervised segmentation of images corrupted with spatially correlated noise Perspective that have been opened in medical research: Fine automated analysis of the biological environment of explanted stents

The GPMF model introduced in Chapter 2 is a way to confidently segment organic biomaterial in a scan corrupted by long-spatially correlated noise. We showed that GPMFs handle the stent artifacts and can offer the best compromise to limit the number of False Positive and False Negative pixels, as compared to other classical approach of segmentation.

We can give two perspectives based on this work. The first would be to introduce a spatial non-stationarity for the range and the variance, i.e., to allow these parameters to vary according to the location on the images. This would indeed better reflect the correlated noise observed on the images since its strenght is not the same everywhere. One could then expect further improvements in the segmentations by, for example, reducing the undesired smoothing effect that can happen. One way to study this issue is to complexify the correlation function so that it integrates non-stationarities [START_REF] Flannery | Estimation of space deformation model for non-stationary random functions[END_REF][START_REF] Kirkpatrick | High resolution simulation of nonstationary Gaussian random fields[END_REF][START_REF] Innocent Nwoye | Weakly supervised convolutional LSTM approach for tool tracking in laparoscopic videos[END_REF]. Another approach might consider building a probabilistic graphical model which can switch between submodels with different parameters (Courbot, Monfrini, et al., 2018) [START_REF] Vacar | Unsupervised joint deconvolution and segmentation method for textured images: a Bayesian approach and an advanced sampling algorithm[END_REF]. Another perspective is based on the construction of an annotated dataset (described in Section 5.4.2) of the images we want to segment8 . This would enable us to construct a discriminative version of the GPMF model related to other models of CRFs recently developed [START_REF] Radosavljevic | Continuous Conditional Random Fields for Regression in Remote Sensing[END_REF] (Krähenbühl and Koltun, 2011) [START_REF] Vemulapalli | Gaussian conditional random field network for semantic segmentation[END_REF] [START_REF] Petrović | Gaussian Conditional Random Fields for Classification[END_REF]. Such a discriminative formulation of the model would potentially improve the resulting segmentations (Ng and [START_REF] Radford | On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes[END_REF] and offer a computational advantage for the inference process.

Stent segmentation in X-ray images

Family of the probabilistic model studied: Hidden Markov Chains Model, innovation: Hidden Markov Chain with specific noise model and path Issue that has been addressed: Unsupervised segmentation of deformed elements without prior information Perspective that have been opened in medical research: Automated analysis of stent deformations and fractures

The fine segmentation of stents in CT and mCT images developed in Chapter 5 faces the problem of artifacts which are a common type of corruption in image processing of biomaterials. We proposed a statistical model dedicated to handling these artifacts and restore the original smooth and tubular stent shape.

On the signal processing aspect, the model we developed could be improved with a noise model which explicitly takes into account the correlated noise such as Pairwise or Triplet Markov models [START_REF] Gorynin | Assessing the segmentation performance of pairwise and triplet Markov models[END_REF]. One could also algorithmically improve the Markov Chain path, which is currently based on a slice-by-slice processing, by extending the path computation algorithm so that it works in full 3D. The problem of stent segmentation could also be approached in a supervised fashion if a database of annotated images were developed. A first Deep Learning algorithm for stent segmentation has been proposed recently for X-ray images in [START_REF] Breininger | Intraoperative stent segmentation in X-ray fluoroscopy for endovascular aortic repair[END_REF].

On the medical aspect, being able to confidently segment stents in artifacts paves the way to further automated stent analyses. The main challenge would consist in using the extracted 3D stent as a basis for numerical biomechanical analyses carried to understand the cause of the stent failures and improve the treatments [START_REF] Langs | Learning deformation and structure simultaneously: In situ endograft deformation analysis[END_REF][START_REF] Maaike A Koenrades | Quantitative Stent Graft Motion in ECG Gated CT by Image Registration and Segmentation: In Vitro Validation and Preliminary Clinical Results[END_REF]. Other types of automated processing over the 3D segmented stent could consist in an automated stent model recognition or stent fracture localization. Such algorithms would enhance the Geprovas workflow of explant analysis.

Appendix A.

Main algorithms in probabilistic modeling

Algorithm A.1: Gibbs sampler (S. Geman and D. Geman, 1984) to sample from an UGM using the full conditionals.

Data: p(x s |x x x x x x x Ns ), ∀s ∈ S, the full conditional equations of an UGM distribution p(x x x x x x x). Result: Series of samples x x x x x x x 0 , x x x x x x x 1 , . . . , x x x x x x x n . n ← 1 Initialize x x x x x x x 0 while convergence is not attained do /* Initialize the new samples */ x x x x x x x n ← x x x x x x x n-1 /* Update at each site */ for s ∈ S do

x n s ∼ p(x n s |x x x x x x x n Ns ) end n ← n + 1 end Algorithm A.3: Simulated Annealing via serial Gibbs sampling (S. Geman and D. Geman, 1984). This algorithm allows only changes in the states towards a state of lower energy. We then easily get stuck in local minima. A practical implementation usually relies on several runs of the algorithm starting from different initial configurations. In Equation 1.4, we make the temperature parameter explicit so that: p T (x x x x x x x|y y y y y y y) = Data: p(x s |x x x x x x x Ns ), ∀s ∈ S the full conditional equations of an UGM distribution p(x x x x x x x), K, the number of Gibbs samples. Result: x x x x x x x the MPM estimator of p(x x x x x x x). while convergence is not attained do for k ∈ {1, . . . , K} do x x x x x x x k ← last sample from Gibbs sampling with p(x x x x x x x) end Frequency estimator of p(x s = ω), ∀s ∈ S, ∀ω ∈ Ω using the realizations (x x x x x x x 1 , . . . , x

x x x x x x K ) xMP M s = argmax ω∈Ω p(x s = ω) end B.

Spectral methods for GMRFs

with N × N blocks, if, and only if, it can be written as: A block-circulant matrix is fully specified by its base or one block column or one block row.

C =      C 0 C 1 C 2 . . . C N -1 C N -1 C 0 C 1 . . . C N -2 . . . . . . . . . . . . C 1 C 2 C 3 . . . C 0      , (

B.2.2. Spectral properties

The following properties linked with the Fourier transform are available for circulant and block-circulant matrices. We will recall the ones for block-circulant matrices, hence we will use the bidimensional Fourier transform.

Let C be a block-circulant matrix with base C b , then we have that C -1 has base (C -1 ) b with:

(C -1 ) b = IDFT2(DFT2(C b ) • (-1)). (B.

B.2.3. Torus assumption for GRFs and GMRFs

The torus assumption, also known as periodic boundary assumption supposes that the boundaries of the grid (formed by the image) are joined together as on a torus. Figure B.1 illustrates the process. The main consequence2 of the torus assumption for a stationary GRF/GMRF is that its covariance matrix is then block-circulant and each block is circulant itself. As a consequence, all the spectral properties for circulant matrices presented in Section B.2.2 can be used. The most used property in our study is the computationally efficient matrix inversion given in Equation B.8. On toruses, one has to modify the distance used. For example, on a two dimensional torus, the Euclidean distance becomes, for a a a a a a a = (x 1 , y 1

                      
Σ (0,0),(0,0) . . . Σ (0,0),(0,lj ) Σ (0,0), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (0,0), (1,lj ) . . . Σ (0,0),(li,0) . . . Σ (0,0),(li,lj ) . . . . . . . . . Σ (0,lj ),(0,0) . . . Σ (0,lj ),(0,lj ) Σ (0,lj ), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (0,lj ), (1,lj ) . . . Σ (0,lj ),(li,0) . . . Σ (0,lj ),(li,lj ) Σ (1,0),(0,0) . . . Σ (1,0),(0,lj ) Σ (1,0), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (1,0), (1,lj ) . . . Σ (1,0),(li,0) . . . Σ (1,0),(li,lj ) . . . . . . . . .

Σ (1,lj ),(0,0) . . . Σ (1,lj ),(0,lj ) Σ (1,lj ), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (1,lj ), (1,lj ) . . . Σ (1,lj ),(li,0) . . . Σ (1,lj ),(li,lj ) . . . . . . . . . . . . . . . . . . Σ (li,0),(0,0) . . . Σ (li,0),(0,lj ) Σ (li,0), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (li,0), (1,lj ) . . . Σ (li,0),(li,0) . . . Σ (li,0),(li,lj ) . . . . . . . . . Σ (li,lj ),(0,0) . . . Σ (li,lj ),(0,lj ) Σ (li,lj ), [START_REF]T-Gibbs complementary experiment data[END_REF]0) . . . Σ (li,lj ), (1,lj ) . . . Σ (li,lj ),(li,0) . . . Σ (li,lj ),(li,lj ) : Building the covariance matrix of an image: the blue arrow highlights the order in which we will consider the neighbors of x 0,0 in the image and their order of appearance in the covariance matrix.

                       . . . Σ ( 0
indexing relation is then, for a l i × l j two dimensional grid:

Cov(x i,j , x i ′ ,j ′ ) = Σ (i,j),(i ′ ,j ′ ) , ∀(i, i ′ ) ∈ {0, . . . , l i } 2 , ∀(j, j ′ ) ∈ {0, . . . , l j } 2 . (B.12)

B.3. Application to GMRF simulation

In this section we illustrate GMRFs by drawing samples from the models in the 2D case using the Fourier based formulas of Section B.2.2. We start by constructing a valid covariance matrix for the field. This can be done by specifying the variance of the field, choosing a correlation function c, parametrized by a range r, and following Algorithm B.1. In Algorithm B.2, we give the algorithm that can be followed to draw samples from the zero-mean Gaussian random field with base precision matrix q q q q q q q. As an illustration, Figure B.3 depicts samples from different stationary GM-RFs. For complete introductions to GRF simulation, one can read [START_REF] Pichot | Algorithms for stationary Gaussian random field generation[END_REF], [START_REF] Brown | Sampling strategies for fast updating of Gaussian Markov random fields[END_REF] or [START_REF] Dirk | Spatial process generation[END_REF].

Appendix C.

Complements on GPMF

C.1. Derivation of the single site equations

In this section we give details on the derivation of the single site equations for the GPMF model, starting from the joint distribution of (X X X X X X X, Y Y Y Y Y Y Y ) given under its energy form in Equation 2.9. The process can be applied for all the processes with a similar energy definition, in particular the other models from the GPMF family.

We start using by the definition the conditional probabilities, for a fixed site s ∈ S we have: (C.1)

Note that in the denominator, we sum over all the realizations of X s and integrate over the domain of Y s , s is fixed; this summation and integration do not happen at other sites. Now, using Equation 2.9 and developing all the terms, we have

p(x s , y s |x N X s , y N Y s ) = exp - s∈S s ′ ∈N X s V (x s , x ′ s ) . . . xs R exp - s∈S s ′ ∈N X s V (x s , x s ′ ) . . . • • • + 1 2 s∈S s ′ ∈S Q s,s ′ (y s -µ xs )(y s ′ -µ x s ′ ) • • • + 1 2 s∈S s ′ ∈S Q s,s ′ (y s -µ xs )(y s ′ -µ x s ′ ) dy s N D (C.2)
Focus on the denominator D, we can extract all the terms that do not imply the fixed site s from the sum and the integral: (C.4)

D = exp - s ′ ∈S s ′ =s
The denominator is now a constant and can be discarded. We focus on the expression: Appendix D.

p(x s , y s |x N X s , y N Y s ) ∝ exp - (s ′ ,s ′′ )∈S×N X s ′ s ′ =s∨s ′′ =s V (x s ′ , x s ′′ ) + 1 2 (s ′ ,s ′′ )∈S 2 s ′ =s∨s ′′ =s Q s ′ ,s ′′ (y s ′ -µ x s ′ )(y s ′′ -µ x s ′′ ) .

Complements on Variational Inference

D.1. Mean Field Variational Inference in SBNs

In the section we derive the parameter update equations that are needed to solve the maximization of the opposite of KL divergence of Equation 3.40. We first isolate terms involving q j from the other:

-KL(q(x x x x x x x)||p(x x x x x x x)) = xj q j (x j )

x x x x x x x\xj i =j q i (x i ) log p(x x x x x x x) xj q j (x j )

x x x x x x x\xj i =j

q i (x i )   k =j
log q k (x s ) + log q j (x j )   .

(D.1)

To maximize -KL(q||p) with the constraints ∀i ∈ S, xi q i (x i ) = 1, we introduce Lagrangian multipliers λ i , ∀i ∈ S, (for each of the constraints) and we define the new target to maximize:

-KL(q(x x x x x x x)||p(x x x x x x x)) = -KL(q(x x x x x x x)||p(x x x x x x x))

+ i λ i xi q i (x i ) -1 . (D.2)
We now consider the functional derivative of -KL(q(x x x x x x x)||p(x x x x x x x)) with respect to q j , and then the derivative with respect to λ j , both derivatives ∀j ∈ S. Term by term we have: q j (x j )

∂ ∂q j    xj q j (x j )
x x x x x x x\xj i =j q i (x i ) log q j (x j )

   =

x x x x x x x\xj i =j q i (x i ) (log q j (x j ) + 1) , = log q j (x j ) + 1, (D.5) (D.6) where Equation D.3 depends on q j through p(x x x x x x x), Equation D.4 depends on q j but Equations D.5 and D.6 do not depend on q j . The second derivative gives:

∂ ∂q j i λ i xi q i (x i ) -1 = i λ i ,
∂ ∂λ j
-KL(q(x x x x x x x)||p(x x x x x x x)) = xj q j (x j ) -1.

(D.7)

We set the derivatives to 0, and we then need to solve the following system for q j , ∀j ∈ S:

             x x
x x x x x\xj i =j q i (x i ) log p(x x x x x x x)k =j log q k (x s )log q j (x j ) -1 + i λ i = 0, xj q j (x j ) -1 = 0.

(D.8)

With elementary manipulations we get:

              
q j (x j ) = exp i λ i -1 exp x x x x x x x\xj i =j q i (x i ) log p(x x x x x x x)k =j log q k (x s ) , xj q j (x j ) = 1.

(D.9)

By integrating the first equation (i.e. summing on all the x j ∈ Ω) and injecting the second in the first we get: (D.12) Plugging Equation D.12 in the first line of Equation D.9, ∀j ∈ S, ∀x j ∈ Ω: Now we remark that all the terms that do not involve x j get simplified. We can start by simplifying the second term of the exponential: xi) [log p(x x x x x x x)] , (D.15) with Z j a normalization constant. We now turn to the precise case of our study: p is the distribution of a SBN (using Equation 3.38). Equation D.15 can be further simplified, ∀j ∈ S: qi(xi) [log p(x x x x x x x)] , where D j is the set of clusters of variables containing variable x j . The expectation distributes over the terms of D j , let us now make explicit one of these expectation terms. For a particular d j = (x s 1 , x s 2 , x s 3 ), we have that one of the three variables is x j . Without loss of generality, if x s 1 = x j : xi) log p dj =

q j (x j ) = 1 Z exp   x x
q j = 1 Z j exp E {xi} i =j ∼ i =j qi(
q j (x j ) = 1 Z j exp E {xi} i =j ∼ i =j
E {xi} i =j ∼ i =j qi(
x s 2 ,x s 3 q s 2 (x s 2 )q s 3 (x s 3 ) log p (xj ,x s 2 ,x s 3 ) . (D.17 where we used the fact that c k r c k = c k =cj r c k + r cj . We want to maximize this quantity with the constraints xi = r(x i |x i -) = xi r ci = 1, ∀i ∈ S. Therefore, similarly to Section D.1, we introduce Lagrangian multipliers λ i , ∀i ∈ S, and the quantity to maximize becomes:

-KL(r(x x x x x x x)||p(x x x x x x x)) = -KL(r(x x x x x x x)||p(x x x x x x x)) + -KL(r(x x x x x x x)||p(x x x x x x x)) =E {rc i } c i =c j ∼ c i =c j rc i log p(x x x x x x x)log r cj -1 with Z j a normalization constant. The expectation appearing in Equation D.28 requires to sample from the joint law x x x x x x x \ {x j , x j -} given x j and x j -. Such a sampling can be done when r is a MT because of the straightforward and sparse decomposition. Note also that, in fact, we only need to sample the variables in each cluster (upon which the expectation is computed, i.e. d j or b j ). Then X X X X X X X does not need to be fully sampled if we recompute the local marginal distributions and the local transition distributions between the MT VI procedure. Note that further developments would consider rewriting Equation D.28 similarly to Equation D.17.

To summarize the steps of the MT-structured variational inference in SBNs:

1. Initialize r for all factors.

D.3. STMT auxiliary variable Variational Inference in SBNs

2. ∀j ∈ S, update r cj with the expression that maximizes Equation 3.42, i.e., Equation D.28.

3. Check convergence and repeat step 2 if needed.

Remark: Because r is not fully factorized, we observe, in the update Equation D.28, as opposed to the MF variational inference, a dependency on variables of r itself. Moreover, for a given update of r cj , the variables of p and r which play a role are all in the Markov blanket of x j .

D.3. STMT auxiliary variable Variational Inference in SBNs

As opposed to the previous sections, we now need to derive the update equations of three variational parameters t ci , t c ′ i and t c ′′ i , ∀i ∈ S, to solve the maximization of Equation 3.45.

Update for t cj

We have, by isolating a particular t cj :

-KL(t||p) = xj ,x j -,v n(j -) (D.29)

t cj x x x x x x x\{xj ,x j -} v v v v v v v\{v n(j -) } ci =cj t ci c ′ s t c ′ s c ′′ s t c ′′ s log p(x x x x x x x, v v v v v v v)-   c k log t c k + c ′ k log t c ′ k + c ′′ k log t c ′′ k     , = xj ,x j -,v n(j -) t cj x x x x x x x\{xj ,x j -} v v v v v v v\{v n(j -) } ci =cj t ci c ′ s t c ′ s c ′′ s t c ′′ s log p(x x x x x x x, v v v v v v v) (A) - xj ,x j -,v n(j -) t cj x x x x x x x\{xj ,x j -} v v v v v v v\{v n(j -) } ci =cj t ci c ′ s t c ′ s c ′′ s t c ′′ s c k =cj log t c k (B) - xj ,x j -,v n(j -) t cj x x x x x x x\{xj ,x j -} v v v v v v v\{v n(j -) } ci =cj t ci c ′ s t c ′ s c ′′ s t c ′′ s c ′ k log t c ′ k (C) - xj ,x j -,v n(j -) t cj x x x x x x x\{xj ,x j -} v v v v v v v\{v n(j -) } ci =cj t ci c ′ s t c ′ s c ′′ s t c ′′
As previously done, we rewrite this objective to maximize introducing a Lagrange multipliers λ i , ∀i ∈ S, to fulfill the constraints, ∀i ∈ S:

xi,v ← i ,v → i t(x i |x i -, v v v v v v v i -)t(v ← i |x i -, v v v v v v v i -)t(v → i |x i -, v v v v v v v i -) = xi,v ← i ,v → i t ci t c ′ i t c ′′ i = 1.
(D.30) We then want to maximize:

-KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v)) = -KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v)) + i λ i   xi,v ← i ,v → i t ci t c ′ i t c ′′ i -1   F , = A -B -C -D -E + F.
(D.31)

Plugging Equation D.29 into Equation D.31, the functional derivatives of each term, with respect to t cj , are:

∂ ∂t cj {A} = E c i =c j tc i c ′ i t c ′ i c ′′ i t c ′′ i [log p(x x x x x x x, v v v v v v v)], ∂ ∂t cj {B} = E c i =c j tc i c ′ i t c ′ i c ′′ i t c ′′ i   c k =cj log t c k   , ∂ ∂t cj {C} = E c i =c j tc i c ′ i t c ′ i c ′′ i t c ′′ i   c ′ k log t c ′ k   , ∂ ∂t cj {D} = E c i =c j tc i c ′ i t c ′ i c ′′ i t c ′′ i   c ′′ k log t c ′′ k   , ∂ ∂t cj {E} = E c i =c j tc i c ′ i t c ′ i c ′′ i t c ′′ i [log t cj + 1], = log t cj + 1, ∂ ∂t cj {F } = i λ i .
(D.32)

The derivative of Equation D.31 with respect to λ j , ∀j ∈ S, is:

∂ ∂λ j -KL(t(x x x x x x x, v v v v v v v)||p(x x x x x x x, v v v v v v v)) = xj ,v ← j ,v → j t cj t c ′ j t c ′′ j . (D.33)
We set Equations D.32 and D.33 equal to 0 and, following the same manipulations as done from Equation D.8 to Equation D.13, we can find the expression of the variational factor t cj minimizing the divergence, ∀j ∈ S: D.3. STMT auxiliary variable Variational Inference in SBNs (D.34) where Z and Z ′ are different normalization constants. Indeed some simplifications were made on on the terms of the sums which do not involve c j . Following the same developments as for other VI procedures, we can exhibit easily the fact that many more terms which do not involve c j can be simplified. Indeed, since STMTs still have a relatively sparse and highly structured associated graph, by carefully selecting the subsets of the variables involved in the expectations, we can straightforwardly simplify the last expression as previously done.

t cj = 1 Z exp   E c i =c j tc i c ′ i t c ′ i c ′′ s t c ′′ s log p(x x x x x x x, v v v v v v v) - c k =cj log t c k - c ′ k log t c ′ k - c ′′ k log t c ′′ k   , = 1 Z ′ exp   E c i =c j tc i c ′ i t c ′ i c ′′ s t c ′′ s log p(x x x x x x x) - c k =cj log t c k - c ′ k log t c ′ k - c ′′ k log t c ′′ k   ,

Update for t c ′ j

Similarly we find: (D.35) where Z and Z ′ are different normalization constants. Again, further simplifications can be made straightforwardly.

t c ′ j = 1 Z exp E c i tc i c ′ i =c j t c ′ i c ′′ i t c ′′ i log p(x x x x x x x, v v v v v v v) - c k log t c k - c ′ k =c ′ j log t c ′ k - c ′′ k log t c ′′ k   , = 1 Z ′ exp E c i tc i c ′ i =c j t c ′ i c ′′ i t c ′′ i log p(v v v v v v v|x x x x x x x) - c k log t c k - c ′ k =c ′ j log t c ′ k - c ′′ k log t c ′′ k   ,
l n m (x n m |x n m-1 ) . . . l 2 (x n 2 |x n 1 )l n 1 (x n 1 ),

(E.8)
and, for a fixed x n m-1 :

x n m-2 ∈Ω l n m-1 (x n m-1 |x n m-2 )l n m-2 (x n m-2 ) ≈ x x x x x x x (m-1) --∈Ω |(m-1) --| l n m-1 (x n m-1 |x n m-2 ) . . . l n 2 (x n 2 |x n 1 )l n 1 (x n 1 )
.

(E.9)

These approximations can be done at each iteration of the update algorithm using the Forward Backward (FB) algorithm (Algorithm A.2)3 to compute the

The manipulation performed in Equation E.10 was already seen in Section 4.3.2. It enables the reinterpretation of the equation as a convolution.

Remark: For completeness, let us give Equation E.6 when formulated in terms of clusters of variables as we did in all the other VI derivations of this thesis. We have, ∀j ∈ S: Constant terms with respect to c j are successively simplified. We have, from the condition on c k , that the only possibility for c k to have a non-empty intersection with c j is that c k = (c j -, c (j -) -). This also leads to Equation E.7. Note that the last line of Equation E.11, without the last summation on log l c k , ∀c k = c j , falls back on the derivation of the MF updates given in Equation 4.6.

  Bold lower case letters represent vectors. A Capital letters represent matrices. ∝ Equality up to a constant factor. δ Kronecker function. Indicator function. | • | Cardinal of a set or absolute value of scalar. . A distance function (Euclidean, ℓ1, ...). det(•) Determinant of a matrix. • T Transpose operator on vectors or matrices. × Product (explicit notation) or Cartesian product. • Element-wise power. ⊙ Element-wise multiplication. DFT2(•) Bidimensional Discrete Fourier Transform. IDFT2(•) Bidimensional Inverse Fourier Transform. ℜ(•) Real part of complex number. ℑ(•) Imaginary part of complex number. ∪ Set union operator. ∩ Set intersection operator. ∧ Logical and operator. ∨ Logical or operator.

  Figure 0.2.: Coupe de CT scan typique des données à traiter. Le stent se trouve dans un environnement complexe, entouré de calcifications et d'artéfacts. Notons qu'un bruit corrélé est une modélisation pertinente de ce phénomène.

Figure 0

 0 Figure 0.3.: Segmentation non-supervisée de la matière organique. De gauche à droite : le mCT observé, la segmentation par champ de Markov caché classique, la segmentation par le nouveau modèle GPMF. La résolution des artéfacts par le modèle GPMF est bien meilleure, elle réduit le nombre de faux-positifs et faux-négatifs.

FigureFigure 0 . 5 .

 05 Figure 0.5.: Exemples de modèles probabilistes graphiques plus généraux.Gauche : Spatial Bayes Network. Droite : fully-connected Conditional Random Field. Les ronds blancs représentent une variable cachée et les carrés gris représentent une variable observée.

Figure 0

 0 Figure 0.6.: Segmentation d'un stent explanté par notre approche.

Figure 0

 0 Figure 0.7.: Segmentation tridimensionnelle d'une artère pathologique. Haut : Image rayons-X originale. Bas : Notre segmentation histologique.

Figure 0

 0 Figure 0.8.: Summary of the atherosclerotic process. From (CNX, 2016).

Biomaterials

[Figure 0

 0 Figure 0.10.: Common biomaterials in vascular surgery: (a) a stent manufactured by Cook 8 , (b) a stentgraft manufactured by Bard 9 , (c) a vascular prosthesis manufactured by Gore 10 .

Figure 0 .

 0 Figure 0.11.: Explanted stent from a superficial femoral artery. The stent itself can be seen, it is encapsulated in a biological material (artery).

  (a) CT: slice of a CT scan which depicts an abdominal crosssection of a human body. A stent implanted in the aorta is visible at the center of the image. Artifacts are also visible, see Figure 0.13 for more details. (b) mCT: cross-section of a stented and explanted superficial femoral artery.(c) Microscopy: histological crosssection of an explanted, stented and totally occluded superficial femoral artery.

Figure 0 .

 0 Figure 0.12.: Illustrations of the CT, mCT and microscopy medical imaging modalities with examples from vascular surgery.

SFigure 1 . 1 .

 11 Figure 1.1.: Examples and notations of graph theory I.

SFigure 1 . 2 .Figure 1 . 3 .

 1213 Figure 1.2.: Examples and notations of graph theory II.

Figure 1 .

 1 3 shows examples of the newly defined graphs.

X 1 Y 1 Figure 1 . 4 .

 1114 Figure 1.4.: Graphical model corresponding to a HMC with Independent Noise (HMC-IN) for N = 5. X 1 is the root node. The white circled vertices are associated with the hidden process X X X X X X X. The gray squared vertices correspond to the variables from the observed process Y Y Y Y Y Y Y .

Figure 1 .

 1 Figure 1.5 shows the graphical model of a dyadic HMT-IN.HMTs and derived models, known as Latent Tree Models, are popular in multi-resolution analysis[START_REF] Kinebuchi | Image interpolation using wavelet based hidden Markov trees[END_REF] (J.-B.[START_REF] Durand | Statistical inference for hidden Markov tree models and application to wavelet trees[END_REF] and in phylogenetic analysis[START_REF] Mourad | A survey on latent tree models and applications[END_REF] [START_REF] Zwiernik | Semialgebraic statistics and latent tree models[END_REF]. To a lesser extent, HMTs have also been used for image processing[START_REF] Salomé | Coregistration of peripheral atherosclerotic plaques assessed by conventional CTangiography, micro-CT and histology in CLTI patients[END_REF] [START_REF] Hanzouli | PECT/CT image denoising and segmentation based on a multi observation and multi scale Markov tree model[END_REF]. Models based on HMTs are presented in Chapter 3 of this thesis.

Figure 1 . 5 . 1 X 1 Figure 1

 15111 Figure 1.5.: Graphical model corresponding to a dyadic HMT with Independent Noise (HMT-IN) for L = 4. X r is the root node. The white circled vertices are associated with the hidden process X X X X X X X. The gray squared vertices correspond to the variables from the observed process Y Y Y Y Y Y Y .

Figure 2 . 1 .

 21 Figure 2.1.: Graphical Models of the classical Hidden Markov Field with Independent Noise (HMF-IN) (a), the Pairwise Markov Field (PMF) with full direct dependencies (b), the Pairwise Markov Field-Uncorrelated Noise (PMF-UN) (c), and the Hidden Markov Field-Correlated Noise (HMF-CN)(d). Within the PMF models, (X X X X X X X, Y Y Y Y Y Y Y )is the Markov process. GPMFs belongs to the PMF family, as well as PMFs-UN and HMFs-CN which are intermediate cases between the HMF and PMF (with full direct dependencies) models in terms of direct dependencies. PMFs-UN and HMFs-CN will be introduced and studied in Section 2.3. We note the numerous correlations that can be introduced between random variables when the pairwise assumption holds.

  b b are real vectors with N elements with generic term, ∀s ∈ S, ∀(x s , x ′ s ) ∈ Ω 2 :

  β = (a a a a a a a T a a a a a a a) -1 a a a a a a a T b b b b b b b.(2.31)Estimator for the other parametersThe Maximum Likelihood (ML) estimator is used to estimate µ 0 , µ 1 , σ. The expressions are, for i ∈ {0, 1}:

  e e T = O . . . N and w w w w w w w T = log Ĉ(0) . . . log Ĉ(N ) .

Algorithm 2 . 1 :

 21 The Stochastic Parameter Estimation (SPE) procedure to train the GPMF model.

  3

  and 2.4 and we now describe and interpret the results.

Figure 2 Figure 2

 22 Figure 2.3.: Error rate in function of varying ∆µ parameters and of several estimators in the supervised segmentation of the 'dude' images of the dataset. Each row, from top to bottom, respectively corresponds to the MAP, MPM (Marroquin's algorithm) and MPM (T-Gibbs algorithm) segmentation criteria. The other parameters of the GRF noise were fixed to σ = 0.5 and r = 2. The dashed red lines are common between the graphs of Figure 2.4, for each row respectively.

Figure 2 . 5 .

 25 Figure 2.5.: MPM (Marroquin's algorithm) unsupervised segmentation over the 'dude' images of the dataset. For the ∆µ row, the other parameters of the GRF noise were fixed to σ = 0.5 and r = 2. For the r row, they were fixed to µ 0 = 0, µ 1 = 0.6 and σ = 0.5. The dashed red lines highlight the situations when the parameter configuration is identical in both graphs. As expected, in such cases, we notice that the models have identical performances relative to each other.

Figure 2

 2 Figure 2.6.: Unsupervised segmentation of images from the dataset with the new model and 4 other models. The segmentation criterion is MPM with Marroquin's algorithm. The error rates with respect to the ground truth appear below each sample. For the 3 images, the noise parameters are µ 0 = 0, µ 1 = 2, σ = 1, r = 3.

Figure 2 . 7 .

 27 Figure 2.7.: Unsupervised segmentations of organic material in corrupted Xrays images. For these 2 cases, which illustrate different explanted stents, we have: the real image y y y y y y y, the ground truth image x x x x x x x, the BM3D-GC segmentation, the P-IN segmentation and the GPMF segmentation. The segmentation criterion for the probabilistic models is the MPM. The stent parts (brightest pixels in y y y y y y y) were segmented beforehand by a thresholding technique and then considered as image borders. They appear in red on the segmented images and did not take part in the segmentation. The ground truth could be provided by experts since a histological analysis is available. Error rates with respect to the ground truth on the whole image appear in parenthesis.

3. 4 . 1 .

 41 Spatial Bayes Networks . . . . . . . . . . . . . .

3. 4 . 2 .

 42 Variational inference . . . . . . . . . . . . . . . 3.4.3. Mean Field Variational Inference in SBNs . . . .3.4.4. Markov Tree Variational Inference in SBNs . . . .3.4.5. Auxiliary variable Variational Inference . . . . . . 3.4.6. STMT auxiliary variable Variational Inference in SBNs . . . . . . . . . . . . . . . . . . . . . . . 3.4.7. Experiments & Results . . . . . . . . . . . . . . 3.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . .

   Goncalves, et al., 3.2. Markov Tree models Algorithm 3.1: Upward-Downward algorithm Data: y y y y y y y, a realization of the observed process, p(x s |x s -), ∀s ∈ S, transitions of the hidden process, p(x r ), r ∈ S 1 , distribution at the root vertice. Result: p(x s |x s -, y y y y y y y), ∀s ∈ S, posterior transitions, p(x s |y y y y y y y), ∀s ∈ S, posterior marginals. /* Compute recursively β(x s ) = p(y y y y y y y s s s s s s s ++ |x s ), ∀x s ∈ Ω, ∀s ∈ S */ β(x s ) = p(y s |x s ), if s ∈ S N , β(x s ) = t∈s s s s s s s + xt β(x t )p(x t |x s ) otherwise.(3.1) /* Compute the remaining posterior transitions β(x s ), ∀(x s , x s -) ∈ Ω 2 , ∀s ∈ S */ p(x s |x s-, y y y y y y y) = β(x s )p(x s |x s-) xs β(x s )p(x s |x s-) , (3.2) /* Compute the posterior marginal at root, ∀x r ∈ Ω */ p(x r |y y y y y y y) = β(x r )p(x r ) xr β(x r )p(x r ) the remaining posterior marginals, ∀x s ∈ Ω, ∀s ∈ S */ p(x s |y y y y y y y) = x s - p(x s -|y y y y y y y)p(x s |x s -enables structured parallel retrieval of the posterior transition distributions (p(x s |x s -, y y y y y y y), ∀s ∈ S) and of the posterior marginal distributions (p(x s |y y y y y y y), ∀s ∈ S). Algorithm 3.1 describes the UD algorithm which makes use of an intermediate quantity β(x s ) = p(y y y y y y y s s s s s s s ++ |x s )

Figure 3 . 1 .

 31 Figure 3.1.: Graphical model corresponding to a dyadic STMT with Independent Noise. (a) depicts all the correlations, (b) gives a condensed view highlighting the preserved MT structure. In (b), a node T T T T T T T s , its father T T T T T T T s -and the root node T T T T T T T r have been annotated. The shadowed squares correspond to the variables from the observed process.

Figure 3 . 3 .

 33 Figure 3.3.: A father node s -linked to its 4 sons (s N W , s N E , s SE , s SW ) in the STMT model. Each node is composed of 9 random variables X s and V s . The directionality of the V s is specified as well as their type (inner, outer I or outer II).

  .26) We use the Maximum Likelihood (ML) estimators µ 0 and µ 1 , ∀i ∈ {0, are also used for σ 0 and σ 1 , ∀i ∈ {0, 1}: σi = 1 s∈S {xs=i} s∈S {xs=i} (y sμxs )

Figure 3

 3 Figure 3.4.: err HM T -INerr ST M T as a function of α and β. The other parameters were fixed to µ 0 = 0, µ 1 = 0.6 and σ 0 = σ 1 = 0.5. We note a segmentation gain (in terms of error rates) raising up to 4 points in favour of the STMT model.

Figure 3

 3 Figure 3.5.: err HM T -INerr ST M T as a function of ∆µ and σ. We note a segmentation gain (in terms of error rates) raising up to 14 points.

Figure 3 . 8 .

 38 Figure 3.8.: Unsupervised segmentation of semi real images with the HMF, HMT and STMT models. Ground truths are on the first row. The percentages indicate the segmentation error rate.

Figure 3 .

 3 Figure 3.10.: Local illustrations of the correlations within the quadtree SBN model: all the sons of a father vertice.

Figure 3 .

 3 Figure 3.11.: Clamped samplings of the original X X X X X X X process (last layer only) from the SBN, STMT and MT models. The first row represents the MRF reference simulations. The other rows represent the clamped samplings. The average difference rate between the realizations of each model and the MRF reference realization is computed over 100 simulations (from which only two are graphically depicted).

  .49) We are interested in retrieving the marginals in the SBN using VI. We successively consider three VI techniques: MF VI (Figure 3.13a), MT VI (Figure 3.13b) and STMT VI (with auxiliary nodes) (Figure 3.13c).

Figure 3 .

 3 Figure 3.12.: Target distributions in the VI procedures: (a) SBN, (b) Adapted SBN with auxiliary nodes.

Figure 3 .Figure 3 .

 33 Figure 3.13.: Variational distributions for the VI procedure: (a) MF VI, (b) MT VI and (c) STMT VI.

GFigure 3 .

 3 Figure 3.15.: Boxplots of the dispersion of the error around the true marginal, for each vertice, for the three VI procedures. We show the results over 1000 different experiments (different p transitions randomly chosen).

Figure 4 .

 4 Figure 4.1.: A specimen from the created dataset: (a) the histological ground truth (optical microscopy), (b) the mCT scan and (c) the annotated mCT scan. Recall that the histological ground truth is not used by the methods described in this chapter. It is only needed for the human expert to correctly annotate the mCT image. The labels on the annotated image are soft tissue, fatty tissue, sheet calcification, nodular calcification, thrombus, specimen holder and background, respectively in blue, green, purple, pink, red, lime and cyan (see Figure 5.15). In the segmentation problem described in Section 4.5.1, the thrombus and soft tissue classes are merged.

Figure 4 . 2 .

 42 Figure 4.2.: The U-Net architecture. Data flow from left to right and first gothrough the contractive path (encoder ) and then the expansive path (decoder ). The network presented here has input images of size 256 × 256 (which is adjustable). Note that the Batch Normalization layer is not present in the original article[START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. Following the contractive path, the layers lose width and height but gain depth: we build a more condensed and higherlevel representation of the data.

Figure 4 . 3 .

 43 Figure 4.3.: The different MC VIs in the image data cube that can be run in parallel. (a) depicts the MC VIs of types front-back and back-front, each involving M × N MCs of length K in both directions of the third dimension. (b) depicts the MC VIs of types top-down and down-top (N × K MCs of length M in both directions of the first dimension). (c) depicts the MC VIs of types left-right and rightleft (M × K MCs of length N in both directions of the second dimension).

Figure 4

 4 Figure 4.4 illustrates the evolution of the mean Dice loss on the train and validation set during training. We can see that, after 50 iterations, the loss is stabilized and the training stops. Table4.1 gives the Dice score for each class on the test set. We can see some important discrepancies in the results between the classes. These discrepancies reflect the intrinsic complexity of the dataset. Indeed, the latter presents high class imbalances, high similarities between classes and notable differences between the image dynamic ranges. Note however that one key goal was to distinguish, on mCT images, the two types of calcifications namely, SC and NC. This is an important histopathological aspect. While such a task seems extremely complex for the naked eye, our DL approach can, to some extent, answer this problem. Future research might further improve this point. An example of a 2D prediction made by the CNN on a slice from the test set is given in Figure4.5.

Figure 4 . 4 .

 44 Figure 4.4.: Mean Dice loss values on the train and validation set as a function of the epoches. Recall that the training is stopped according to the behavior of the validation loss. ST FT SC NC SH Ba Dice 0.94 0.41 0.85 0.64 0.86 0.99Table 4.1.: Dice scores on test set for each class.

Figure 4 .

 4 6b, MC VI offers up to a 8 point gain in the 4.5. Experiments and Results (a) MicroCT slice (b) CNN predictions (c) Ground truth

Figure 4 . 5 .

 45 Figure 4.5.: Example of 2D prediction by the CNN on a sample from the test set. It is the same sample as illustrated in Figure 4.1 where, as noted, the classes soft tissue and thrombus have been merged. The colormap is given in Figure 5.15.

  err M F V I err M CV I-P e err M CV I-P a (a) w1 = w2 = 1 w1 = 1, w2 = 2

Figure 4 .

 4 Figure 4.6.: Error rates in unsupervised segmentation of the two VI procedures as a function of the noise level σ = σ ω1 = σ ω2 . The other parameters are fixed toµ ω1 = 0, µ ω2 = 1, θ α = θ β = θ γ = 1.The new MC VI methodologies appear to be the best performing methods for all noise levels. The results are here averaged on 10 randomly chosen images from the dataset.

Figure 4 . 7 .

 47 Figure 4.7.: Timing of the VI procedures as a function of the size of the image (a size n means an image of dimensions 2 n × 2 n ).

Figure 4 . 8 .

 48 Figure 4.8.: err M F V Ierr M CV I as a function of w 1 and w 2 . The other parameters areµ ω1 = 0, µ ω2 = 1, σ ω1 = σ ω2 = 1.5, θ α = θ β = θ γ = 1.The results are averaged over 10 simulations.

Figure 4 .

 4 Figure 4.9.: Unsupervised segmentation of semi-real images with the two VI procedures on variable images and noise levels σ = σ ω1 = σ ω2 . We have w 1 = 1, w 2 = 2, µ ω1 = 0, µ ω2 = 1, and θ α = θ β = θ γ = 1.The parallel MCs approach was taken for the MC VI.

Figure 4 .

 4 Figure 4.10.: Example of a 3D prediction by the CNN on a mCT from a test artery. We note that the irregularities found in the fatty tissue class in green are to be expected since this class exhibit the worst Dice score. The support class has been made invisible for clarity and the opacity of the soft tissue class has been lowered. The colormap from Figure 5.15 is used.

Figure 4 .

 4 Figure 4.11.: Illustration of spurious classifications CNN of the soft tissue class. While clearly visible by the naked eye, the soft tissue class is only partially segmented on this slice. A major cause of this problem might be the lack of context from the depth axis in the slice-byslice segmentation performed by the CNN. The color code from Figure 5.15 is used.

4. 5 .

 5 Experiments and Results(a) x x x x x x x M F V I (b) x x x x x x x M CV I

Figure 4 .

 4 Figure 4.12.: Illustration of fcCRF post-processings in 3D of the image in Figure 4.11. In (a), the case of the MF VI, we can see that the spatial context offers some improvement to solve the misclassification from the soft tissue class. However this comes at the price of high kernel variances for smoothing which leads to loss of details in the calcifications parts which are the most important components of the image. However in (b), the case of the MC VI with modified transition matrix, we can see that our proposed approach offers the best improvements while preserving the most the details within the segmented calcifications. None of the approaches could recover the missing calcification in the segmentation. For both VIs the parameters were set to w 1 = 2, w 2 = 2, θ α,x = θ α,y = θ α,z = 4, θ γ,x = θ γ,y = θ γ,z = θ β = 2. The color code from Figure 5.15 is used.

Figure 4 .Figure 4 .

 44 Figure 4.13.: Illustration of fcCRF post-processings in 3D (full volumes). The same parameters and colormap as in Figure 5.15 are used. The smoothest segmentation is obtained in (c) with the proposed MC VI approach thanks to its ability to better correct the spurious CNN classifications.

Figure 5 .

 5 Figure 5.1.: Example of input data (2D views), we see the complex environment in which the stent lies: artifacts and calcifications are notable.Some stent, calcification and artifact components are indicated by, respectively, red, blue and green arrows. We note the spatial resolution differences between CT scans and mCT scans.

Figure 5 .

 5 Figure 5.2.: The HMC path, to transform 3D data into 1D data, illustrated on 3 successive slices.

  be the vector of random variables representing the observed grayscale image (the X-ray scan), the elements of Y Y Y Y Y Y Y take their value in R.

Figure 5 . 3 .

 53 Figure 5.3.: Diagram of the fine stent segmentation step.

Figure 5 .

 5 Figure 5.3 summarizes the complete segmentation approach.

Figure 5 . 4 .

 54 Figure 5.4.: 3D segmentation of stents via manual thresholding (Thresh.), MoG approach and MoE approach.

Figure 5 .

 5 Figure 5.5.: 3D views of segmented mCT and CT stent images from the database.

Figure 5

 5 Figure 5.6.: MAR with the BHC algorithm when using different stent segmentation techniques. The last row, which is the one including our segmentation technique, shows restored images less affected by the artifacts. This can be seen especially in the center region of the stent. Artifacts may be hard to see on the original slices because of the image dynamic range. The scales used in the Frangi filter are {3, 6, 9}. The scales used in the RORPO filter are {100, 150, 200}.

  and 5.8 illustrate some more segmentations with the classical HMF model, the Potts-Independent Noise (P-IN) model, and the new GPMF model. It is notable that the overall performance of the GPMF model seems worst in Figure 5.8 than in Figure 5.7. Indeed, the latter segmentations are examples of the current limitations of the GPMF model that we associate with the noise stationarity assumption made in the model.

Figure 5 . 7 .Figure 5 . 8 .Figure 5

 57585 Figure 5.7.: P-IN and GPMF segmentations of organic material in mCT on more examples. The blue areas correspond to the areas of FN/FP computations. The indicated percentage represents the error rate in the segmentation computed over the whole image.

Figure 5 .

 5 Figure 5.10.: Case A. Unsupervised segmentations of atherosclerotic calcifications. The percentages below the segmentations indicate the error rates with the ground truth. The ground truth is shown in Figure 5.9a.

Figure 5 .

 5 Figure 5.11.: Case B. Unsupervised segmentations of atherosclerotic calcifications. The percentages below the segmentations indicate the error rates with the ground truth. The ground truth is shown in Figure 5.9b.

Figure 5 .

 5 Figure 5.12.: Case C. Unsupervised segmentations of atherosclerotic calcifications. The percentages below the segmentations indicate the error rates with the ground truth. The ground truth is shown in Figure 5.9c.

Figure 5 .

 5 Figure 5.13.: Pictures of two explanted arteries after amputation, the first step of our protocol.

  (a) Three histologic slices resulting from Step 3. (b) Three microCT images correlated with their histologic truth after Step 4. (c) Three expert annotated microCT images obtained after Step 5.

Figure 5 .Figure 5 .

 55 Figure 5.14.: Illustrations of the different steps of the protocol on three examples.

Figure 5 .

 5 Figure 5.16.: Precision-Recall curves and their AUCs obtained by the CNN for each class on the test set.

  B.6) where C i is a circulant n × n matrix with base c c c c c c c i . The base of C is the n × N matrix C

8 )

 8 Let C and D be two block-circulant matrices with bases C b and D b , then we have CD = IDFT2(DFT2(C b ) ⊙ DFT2(D b )). (B.9) Let C be a block-circulant matrix with base C b , then Λ = DFT2(C b ) (B.10) is the matrix of the eigenvalues of C.

  (a) The image representation with a matrix

  Figure B.2.: Building the covariance matrix of an image: the blue arrow high-lights the order in which we will consider the neighbors of x 0,0 in the image and their order of appearance in the covariance matrix.

  p(x s , y s |x N X s , y N Y s ) = p(x s , y s |x S\{s} , y S\{s} ),xs∈Ω R p(x s , y s , x S\{s} , y S\{s} )dy s x x = (x s , x S\{s} ), y y y y y y y = (y s , y S\{s} ).

(C. 5 )

 5 Before further simplifications, note the fact that to get Equations C.3, C.4 and C.5 we used the following equality on sets, for a fixed s ∈ S:{(s ′ , s ′′ ) : s ′ ∈ S ∧ s ′′ ∈ S} ={(s ′ , s ′′ ) : s ′ ∈ S ∧ s ′′ ∈ S ∧ (s ′ = s ∧ s ′′ = s)} ∪ {(s ′ , s ′′ ) : s ′ ∈ S ∧ s ′′ ∈ S ∧ (s ′ = s ∨ s ′′ = s)}.Error rate as a function the iterations of Gibbs and T-Gibbs sampler in the segmentation of the observed images of Figure C.1 with the P-GMRF model. The purple vertical lines indicate iterations explored in (b), where we represent, from left to right, Iteration 0, Intermediate iteration and Iteration 30. Illustration of some iterations given in (a) to highlight the different behaviors of the Gibbs and T-Gibbs samplers.

Figure C. 2 .

 2 Figure C.2.: Supervised segmentation experiment with the P-GMRF model, to compare the T-Gibbs and Gibbs samplers in typical cases. The final segmentation after 30 iterations of both samplers is better in the case of the T-Gibbs sampler. In a first typical scenario (found in the dude12occ4 case) the T-Gibbs sampler avoids a total loss of the image details and ends up with a much better error rate. The second case (dog45 case) illustrates, in the case of the T-Gibbs, a faster convergence and a final result which conserves more details despite a similar error rate to the Gibbs sampler.

  E {xi} i =j ∼ i =j qi(xi) ds log p ds , i =j ∼ i =j qi(xi)

)D

  Using this result in the general expression of the variational parameter, ∀j ∈ x x\{xj ,x j -} ci =cjr ci log r cj x x\{xj ,x j -} ci =cj r ci c k =cj log r c k ,(D.19)

  to find the functional derivative of Equation D.20. Plugging Equation D.19 in Equation D.20, we exhibit the functional derivatives with respect to each r cj for each of the terms: x x\{xj ,x j -} ci =cj r ci log p(x x x\{xj ,x j -} ci =cj r ci log r cj x x\{xj ,x j -} ci =cj r ci (log r cj + 1the functional derivative of Equation D.20 is, ∀j ∈ S: ∂ ∂r cj

-

  of Equation D.20, with respect to λ j , is, ∀j ∈ S:∂ ∂λ j -KL(r(x x x x x x x)||p(x x x x x x x)) = xj r cj -1. (D.26)We set Equations D.25 and D.26 equal to 0 and, following the same manipulations as done from Equation D.8 to Equation D.13, we can find the expression of the variational factor r cj minimizing the divergence, ∀j ∈ S:log r cj = 1 Z exp   E {rc i } c i =c j ∼ c i =c j rc i a normalization constant.This expression can be simplified because many of the terms do not depend on c j . Let D j be a set whose elements are clusters of variables d j containing x j , B j a set of clusters of variables b j also containing x j but with the condition b j = c j . Equation D.27 simplifies in: r cj = r(x j |x j -i } c i =c j ∼ c i =c j rc i

  u (x j ) + E c i =c j lc ii =j ψ p (x j , x i )c k =cj log l c k   , with c k s.t. c k ∩ c j = ∅ and c k = c j ,(E.11)

F. 1 .

 1 Figure F.1.: Histograms for the realizations of the random variables for the stent and rest classes, on a typical data slice, after the preprocessing segmentation. The density values are divided by the total number of pixel in the slice on the three histograms.

  3. The Euclidean distance Chapter 2. Gaussian Pairwise Markov Fields on the torus (whose dimensions are the size of the image) is then used as defined in Equation B.11. The assumption of Markovianity for

Table 2 .

 2 s ′ ȳs ȳs ′ . 1.: Factorization of the related intermediate models.

	Model P-IN	p(x s , y s |x x x x x x x Ns y y y y y y y Ns ) factorizes using Eq. 2.16 and Eq. 2.18
	PMF-UN	Eq. 2.21 and Eq. 2.18
	P-GMRF	Eq. 2.16 and Eq. 2.19

  states of the parallel Markov chains.

	Result: x x x x x x x i 1 , the sample at the model temperature.
	i ← 1 while not converged do /* Gibbs sampler for each chain */
	for k ∈ {K, . . . , 1} do x x x x x x x i+1 k drawn by Gibbs sampler from p T k (x x x x x x x) with initialization at x x x x x x x i k
	end
	/* Test to swap chains */
	for

Table 2 .

 2 3.: FN and FP rates computed in the blue areas, for each model, for each case of Figure 2.7.

  Contents 3.1. Introduction . . . . . . . . . . . . . . . . . . . . . . 3.2. Markov Tree models . . . . . . . . . . . . . . . . . 3.2.1. Hidden Markov Trees . . . . . . . . . . . . . . .

3.2.2. Spatial Triplet Markov Trees . . . . . . . . . . . 3.3. Image segmentation . . . . . . . . . . . . . . . . . 3.3.1. The Potts-like transition distributions . . . . . . 3.3.2. Iterative Parameter Estimation for Trees . . . . . 3.3.3. Experiments and Results . . . . . . . . . . . . . 3.4. STMTs for Auxiliary Variational Inference in SBNs

  with outer inner Figure 3.2.: Details of the STMT construction: a father node s -linked to its 2 sons (s L , s R ). The directionality of the V s is specified as well as their type (inner or outer).

	assume that we have the factorization:
	p(t t t t t t t s |t t t t t t t s-

  .1. At a site s ∈ S, should it be North West or North East or South East or South West, we find the couple (X s

  Unsupervised segmentation of semi-real images with HMTs, STMTs and HMFs. The figure depicts the error rate of the three models with a varying additive Gaussian noise level. It appears that the STMT model always performs equally or better than the HMT model, offering up to a 5 point improvement in the error rates with respect to the HMT model. The HMF model is the best performing model offering up to a 10 point gain over the STMT model.

		3.4. STMTs for Auxiliary Variational Inference in SBNs
	Error rate to truth	0.00 0.20 0.40
		0.2	0.4	0.6	0.8	1	1.2	1.4
					σ
			HMF		HMT		STMT
	Figure 3.6.: Error rate in unsupervised segmentation
	Figure 3.7.:
								.30)
	The description and parameter estimation for a HMF-IN model without bias
	has been developed in Section 2.3. The extension to estimate the bias parameter
	α is straightforward by considering a two-unknown-variable LLS as already
	mentioned in Section 3.3.1.

  Chapter 4. Image segmentation with Deep Learning and Conditional Random Fields

	(256, 256, 1)	(256, 256, 8)	(256, 256, 8)					(256, 256, 16)	(256, 256, 8)	(256, 256, 8)	(256, 256, 6)
	Input							Output
			(128, 128, 8)	(128, 128, 16)	(128, 128, 16)			(128, 128, 32)	(128, 128, 16)	(128, 128, 16)
					(64, 64, 16)	(64, 64, 32)	(64, 64, 32)	(64, 64, 64)	(64, 64, 32)	(64, 64, 32)
							(32, 32, 32)	(32, 32, 64)	(32, 32, 64)	(32, 32, 128)	(32, 32, 64)	(32, 32, 64)	2D Convolution +Batch Normalization +Relu Activation 2D Max Pooling
								2D Up Convolution
								(16, 16, 64)	(16, 16, 128)	(16, 16, 128)	Copy +Concatenate on 3 rd axis height and depth of a layer (x, y, z) Triplets indicate width,

Table 5 .

 5 1.: FN and FP rates computed in the blue areas for each model, for each case ofFigures 5.7 and 5.8.

	5.3. Unsupervised segmentation of the organic
	material and calcifications
	5.3.1. Organic material segmentation with GPMFs

Table 5 .

 5 2.: Multiclass confusion matrix of the CNN results of the test set. Each cell has been normalized by the total number of pixels in the test set so that the sum of all the cells gives 1. The numbers in the cells must be multiplied by ×10 -2 .

	5.4. Histological segmentation of microCT with Deep Learning
	❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ Truth Pred. ST	FT	SC	NC	SH	Ba
	ST	9.94 0.13	0.07	0.02	0.02	0.2
	FT	0.50 0.31 0.002 0.0006 0.002	0.004
	SC	0.06 0.007 0.78	0.11	0.0009 0.0005
	NC	0.05	0	0.17	0.25	0.001	0.002
	SH	0.15 0.34 0.008 0.003	2.73	0.32
	Ba	0.17 0.06 0.003 0.001	0.19	83.4

  (s ′ ,s ′′ )∈S 2 s ′ =s∧s ′′ =s Q s ′ ,s ′′ (y s ′µ x s ′ )(y s ′′µ x s ′′ ) (s ′ ,s ′′ )∈S 2 s ′ =s∨s ′′ =s Q s ′ ,s ′′ (y s ′µ x s ′ )(y s ′′µ x s ′′ ) dy s .The first exponential term in D now gets simplified with terms from the numerator N . What remains then in N is the same terms as in the second exponential in D. So we can write:p(x s , y s |x N X s , y N Y ′ ,s ′′ )∈S 2 s ′ =s∨s ′′ =s Q s ′ ,s ′′ (y s ′µ x s ′ )(y s ′′µ x s ′′ ) ′ ,s ′′ )∈S 2 s ′ =s∨s ′′ =s Q s ′ ,s ′′ (y s ′µ x s ′ )(y s ′′µ x s ′′ ) dy s .

						V (x s ′ , x s ′′ )
					s ′′ ∈N X s ′ s ′′ =s
	+	1 2
							(C.3)
	×	xs	R	exp -	s ′ =s∨s ′′ =s s ′ (s ′ ,s ′′ )∈S×N X	V (x s ′ , x s ′′ )
	+	1 2
			• • • + 1 2 (s • • • + 1 2 (s

s ) = exp -(s ′ ,s ′′ )∈S×N X s ′ s ′ =s∨s ′′ =s V (x s ′ , x s ′′ ) . . . xs R exp -(s ′ ,s ′′ )∈S×N X s ′ s ′ =s∨s ′′ =s V (x s ′ , x s ′′ ) . . .

Les travaux de cette thèse utilisent la base de données d'explants du laboratoire Geprovas (https://geprovas.org) qui a vu le jour grâce à son programme de collecte et d'analyse d'explants.

Polytetrafluoroethylene (PTFE) and Polyethylene Terephthalate (PET) are thermoplastics. Note that their biocompatibility with the human body is debated[START_REF] Chakfé | Editor's Choice-European Society for Vascular Surgery (ESVS) 2020 Clinical Practice Guidelines on the Management of Vascular Graft and Endograft Infections[END_REF].

https://www.cookmedical.com/products/224e3666-308f-4244-8695-6fd23bbd671c/

https://www.crbard.com/Peripheral-Vascular/en-US/Products/ FLUENCY-Plus-Endovascular-Stent-Graft

https://www.goremedical.com/me/products/vgstretch

https://geprovas.org

The definition straightforwardly extends to higher dimensions.

The use of these modalities for medical diagnosis and treatment defines the medical field of radiology.

Provided we tolerate the ionizing radiations which are problematic for some patients[START_REF] Mark S Pearce | Radiation exposure from CT scans in childhood and subsequent risk of leukaemia and brain tumours: a retrospective cohort study[END_REF] but also clinicians[START_REF] Brun | Radiation protection for surgeons and anesthetists: practices and knowledge before and after training[END_REF].

15 Some illustrations of the possibilities offered by mCTs in the imaging of small animals are depicted in (Roque-Torres, 2020).

Note that this last notation shortcut then confuses p with P X but no ambiguity is raised in the work of this thesis.

Also called Maximum Posterior Mode (MPM).

Broyden-Fletcher-Goldfarb-Shanno, the names of the main contributors to the algorithm.

s -1 is the father node of s in the associated gaphical representation, seeFigure 1.4.

The Potts and Ising (J. E.[START_REF] Besag | Nearest-neighbour systems and the auto-logistic model for binary data[END_REF] models are equivalent when the hidden stochastic process is binary valued.

exp -1 2 (x x x x x x xµ µ µ µ µ µ µ)Σ -1 (x x x x x x xµ µ µ µ µ µ µ) .(2.1)

The conditions on β and Q for the SPDness of R are respected, see Section

2.2.4.

Following Appendix C.1, it can be shown that the P-GMRF model arises from the joint energy E ′ :

This constitutes a generalization of the approach since we derive the estimators starting from the posterior distributions as opposed to the prior distribution as done in the original article.

This step of inference is here carried stochastically since an exact analytical expression has not been established.

https://vision.lems.brown.edu/content/available-software-and-databases

MPM has already been recognized as more stable in practical applications that the MAP estimator[START_REF] Courbot | Pairwise Markov fields for segmentation in astronomical hyperspectral images[END_REF].

Implementation from https://opencv.org.

Implementation from https://github.com/Borda/pyImSegm.

Implementation from http://www.cs.tut.fi/~foi/GCF-BM3D/.

Implementation from https://github.com/shaibagon/GCMex/.

The UD algorithm can be interpreted in two ways. Either as a generalization of the Forward

This section is a new model which builds upon and redefines the only occurence of the STMT prior to our work, found in(Courbot, Monfrini, et al., 2018).

In such case the Y process is taken into account and we consider the probability of a decuplet (Xs, V V V V V V V s, Ys), ∀s ∈ S L .

Equations from sites from other resolutions could be derived and integrated to the LLS however we found out that they do not contribute significantly to the final solution.

Strictly speaking, there should be a bias parameter for each of the class. However, in this binary classification context, the proposed parametrization can be seen as forcing the bias parameter for class ω 1 to be 0. α is then the bias parameter for class ω 0 and it is estimated under this constraint.

Note that if p is not augmented with auxiliary random variables, certain terms of t involving the distributions of auxiliary random variables cannot be updated.

In Section 3.4.7, we show how p can integrate auxiliary random variables in the SBN case.

In the sense that some semantic classes appear much less often than others.

Typically 10% of the training set.

The full demonstration can also be found in the supplementary material of[START_REF] Krähenbühl | Efficient inference in fully connected crfs with gaussian edge potentials[END_REF].

https://www.tensorflow.org

Note, moreover, that the use of the Peano scan can only be done on images whose dimensions are a power of two, thus making this approach less practical.

 6 We empirically found that these parameters are more critically related to the segmentation performance than the kernel standard

deviations. 7 https://vision.lems.brown.edu/content/available-software-and-databases

When this project was carried at the beginning the thesis, we had very few and only unannotated data at our disposal.

This filter response over the data cube is manually segmented afterwards, with a global threshold.

We started our study with six explanted arteries.

https://www.nikonmetrology.com/en-gb/product/xt-h-255-st

https://www.volumegraphics.com/

https://imagej.nih.gov/ij/

https://www.gimp.org

This dataset was not available at the time of the work on the GPMF model.

Potentially taken on the torus, see Section B.2.3.

Under an appropriate ordering of the indices, see Section B.2.4.

Under an appropriate ordering of the indices, see Section B.2.4.

Typically I = [0,

255] for grayscale images, or another integer interval depending on the image format. I can also be vectorial for color images.

Note that in this case, we use the FB algorithm without any observations, this is equivalent to taking all the conditional likelihood terms of the equations of Algorithm A.2 as equal to 1.

Acknowledgements

θ θ θ θ θ t ←= {α t , β t , µ t 0 , µ t 1 , σ t 0 , σ t 1 } t ← t + 1 end

Experiments and Results

Synthetic data experiment

In this first experiment, we assess experimentally the strict generalization of STMTs over HMTs-IN in the case of the Potts-like potential. We sample realizations of the STMT model and segment the observed images with both models. The segmentation criterion is the MPM and the parameters used for the segmentation are the parameters used for the simulation, the true parameters. We are in the context of supervised segmentation. Again the hypothesis of the independent Gaussian noise parametrized by the underlying class (Section 3.3.2) is done for both models. Let us denote the segmentation error rate of the HMT-IN model, resp. STMT model, err HM T -IN , resp. err ST M T .

In the first experiment, α and β are the varying parameters in the STMT samples. We plot the quantity err HM T -INerr ST M T as a function of α and β in Figure 3.4. We first notice that the values taken by the surface are always positive which means that the STMT segmentations are always better or equal to the HMT-IN segmentations (the maximal value is a 4-point improvement in Chapter 4.

Image segmentation with Deep Learning and Conditional Random Fields

Contents Remark: A comparison with the results we obtained on semi-real images in Section 2.5.2 corroborates the fact that the stationarity assumption in the GPMF model is the most limiting factor. Indeed, in the semi-real experiment, the noise stationarity assumption holds and we do not observe any of the problems that we observe here. The STMT model from Chapter 3 illustrates the interest of auxiliary random variables in probabilistic models. We explored to what extent this model could be a close alternative to MRF models with the capability of performing exact inference. This could be an asset over an approximate inference algorithm when dealing with large data such as medical images. For a better understanding of STMTs, we showed its closeness to the SBN model in the context of auxiliary random variable VI.

Further research on the proximity of STMTs and SBNs might follow the direction of deriving theoretical results quantifying the relations between the two models. One might also consider extending our auxiliary random variable VI algorithm to quadtrees to process images. Such an algorithm is expected to be very efficient based on the preliminary results that we obtained in the case of dyadic trees.

Conclusions and Perspectives

Histologic segmentation of atherosclerotic arteries with Deep Learning

Family of the probabilistic model studied: Conditional Random Fields Model, innovation: Improved inference in Conditional Markov Random Fields Issue that has been addressed: Improving the inference in a probabilistic model for post-processing the bidimensional segmentations of a Convolutional Neural Network to reconstruct a tridimensional segmentation Perspective that have been opened in medical research: A first histological analysis of the vascular segment using the X-ray image only Thanks to the U-Net approach and the fcCRF post-processing described in Chapter 4, we reconstructed tridimensional arteries with annotated histologic elements. This study confirms the important role that microCT images could have in vascular research.

The main perspectives in this work is first to increase the score in the segmentation of each of the classes (particularly, in the current implementation of our model, the fatty tissue class) but also to integrate more classes and offer a more detailed histologic segmentation. As an example, it would be interesting to treat arteries containing a stent. The neural network approach, combined with a pre-or post-processing approach, could additionally segment the stent in the images as well as to handle the artifacts such images would contain. To this end, one could consider exploring recent advances in CNN, such as more recent versions of the U-Net network [START_REF] Zahangir Alom | Recurrent residual convolutional neural network based on u-net (r2u-net) for medical image segmentation[END_REF] [START_REF] Oktay | Attention u-net: Learning where to look for the pancreas[END_REF]. One could also study the interest of including shape priors within the CNN (Nosrati and Hamarneh, 2016) [START_REF] Oktay | Anatomically constrained neural networks (ACNNs): application to cardiac image enhancement and segmentation[END_REF][START_REF] Ravishankar | Learning and incorporating shape models for semantic segmentation[END_REF]. Shape priors could be elaborated jointly with the experts collaborating on the project. Algorithm A.2: Forward Backward (FB) algorithm, rescaled version from [START_REF] Pierre | Baum's forward-backward algorithm revisited[END_REF]. The original unscaled algorithm from [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF]) is prone to underflows; it lacks the rescaling by the term κ s . Note that in the original unscaled version of FB we have α(x s ) = p(x s , y 1 , . . . , y s ) and β(x s ) = p(y s+1 , . . . , y N |x s ), ∀s ∈ S. Data: y y y y y y y, a realization of the observed process, p(x s |x s -), ∀s ∈ S, transitions of the hidden process, p(x s |x s -), ∀s ∈ S, p(x r ), distribution at the root vertice. Result: p(x s |y y y y y y y), ∀s ∈ S, the posterior marginals. /* Compute the rescaled forward probabilities α * (x s ) = p(x s |y), ∀s ∈ S, the recursion is defined by */

/* Compute the rescaled backward probabilities β * (x s ) = p(ys+1,...,yN |xs) p(xs+1,...,x N |y1,...,ys) . The recursion is defined by */ Appendix B.

Complements on GMRFs

B.1. Stationary GRFs and GMRFs

In this section we introduce GRFs and GMRFs, the latter is derived from GRFs with an important assumption that is made about the covariance matrix. First, we need to define the notion of positive-definite function. A positivedefinite function of a real variable is a complex-valued function f : R → C such that for any numbers y 1 , . . . , y n , the n × n matrix A = (a i,j ) n i,j=1 , a i,j = f (y iy j ) is positive semi-definite. In our application, f : R → R and the resulting A will be a real SPD.

A GRF or GMRF is called stationary if the following conditions hold:

• the mean vector is constant (does not depend on the vertice s, ∀s ∈ S).

• the covariance between two random variables in the field only depends on the norm (ℓ1, Euclidian...) between those two points, such that, for some positive-definite function f :

In Equation B.1, f is also called the covariance function. In the case of stationary GRFs or GMRFs, the covariance function is directly linked to a correlation function c such that:

where σ 2 ∈ R + is a constant variance (imposed for stationary GRFs or GMRFs [START_REF] Abrahamsem | A review of Gaussian random fields and correlation functions[END_REF]). Note that for unitary variance, the correlation function and the covariance function are equal.

By definition, we require the covariance function to be a positive-definite function, then Equation B.2 shows that a correlation function c is also a positive-definite function. Let h be a distance function (Euclidean or ℓ1 for example 1 ). Some classical correlation functions are, ∀(s, s ′ ) ∈ S 2 :

• the exponential correlation function:

• the Gaussian correlation function:

The parameter r is called the correlation length, decay or range. Note that the Gaussian and exponential correlation functions can both be derived from the Matérn correlation function [START_REF] Abrahamsem | A review of Gaussian random fields and correlation functions[END_REF], a more general correlation function. Moreover, these correlation functions are given under their isotropic form. An example of anisotropic Gaussian or exponential correlation function would consist in setting different correlation lengths according to the direction. This is considered for example in (Dietrich and Newsam, 1997).

In this thesis, to perform segmentation with the GPMF model, we consider isotropic stationary covariance functions but with non-stationary mean vector. This is a case of non-stationarity much easier to handle than a non-stationary covariance matrix [START_REF] Fuglstad | Exploring a new class of non-stationary spatial Gaussian random fields with varying local anisotropy[END_REF].

Remark: The modeling of non-stationary GRFs and GMRFs is an active research area out of scope of this thesis. One can see [START_REF] Gelfand | Handbook of spatial statistics[END_REF][Chapter 10] or (Risser, 2016) for an overview.

B.2. Spectral methods for GMRFs

This section is a condensed view of the development made in (Rue and Held, 2005) to establish computationally efficient formulas for GMRF manipulation, which is essential for the tractability of the algorithms developed in Chapter 2. These formulas are based on the Fourier transform and on a periodic boundary assumption, also called torus assumption. The latter assumption leads to a special form of the covariance matrix. We start by reviewing key elements of algebra over these matrices.

B.2.1. Circulant matrices

We here introduce circulant matrices and their related spectral properties. A n × n matrix C is circulant if, and only if, it has the form:

All the MRF models of Chapter 2 are developed under the torus assumption.

Remark: Another approach to using the spectral properties of circulant matrices without making the torus assumption is via circulant embedding. The method is based on the observation that, in the example of the 2D case, the covariance matrix of the GMRF will be block-Toeplitz with Toeplitz blocks 3 . This matrix is then embedded in a block-circulant with circulant block matrices and the spectral methods can be used. A presentation of this method in the 1D case is available in (Dietrich and Newsam, 1997) where the method was introduced and in (Kroese and Botev, 2013) for the 2D case.

B.2.4. Indices ordering

As noted before, for a Toeplitz or circulant matrix structure to appear for the covariance matrix, one needs to use an appropriate indexing of the matrix.

Focusing on covariance matrices of two dimensional GRF/GMRF, an entry in the covariance matrix will be stored in row major order, the mathematical Algorithm B.1: Construction of q q q q q q q Result: q q q q q q q, the base of the precision matrix for the GMRF. for i = {1 . . . l i } do for j = {1 . . . l j } do q q q q q q q i,j ← σc (d torus ((0, 0), (i, j)), r) end end Algorithm B.2: Sampling a GMRF with via the Fourier transform properties Data: q q q q q q q, the base of the precision matrix of the model. Result: x x x x x x x, a sample of the GMRF. Sample z z z z z z z, a 0-mean complex field with independent elements distributed according to: ℜ(z i ) ∼ N (0, 1), ℑ(z i ) ∼ N (0, 1) /* Matrix of the eigenvalues */ Λ ← DFT2(q q q q q q q) /* Perform the sampling */

C.2. T-Gibbs sampler complementary experiment

For a fixed s ∈ S, it can also be rewritten:

Now it is important to note the following relation on cardinality, caused by the symmetry of the elements. For a fixed s ∈ S:

The same equalities slightly modified can be written for the set involving the set of the neighbors (recall that, by definition, the set of neighbors of the site s does not contain s), for a fixed s ∈ S:

with the following relation on cardinality:

By definition, in Equations C.9 and C.10, both s ′ and s ′′ can not simultaneously be equal to s. Finally, using the preceding remarks on sets, we can rewrite Equation C.5 in a simplified form, always for s ∈ S fixed:

(C.11)

Remark: Note that Equations C.10 and C.11 only consider the current case where all the cliques are of size 2.

C.2. T-Gibbs sampler complementary experiment

In this section, we consider the P-GMRF model which tends to produce, for particular initializations, samples that worsen. This can be stabilized using the T-Gibbs sampler and the parametrization described in Section 2.5.1. We consider the supervised segmentation of the dude12occ4 and dog45 images from the dataset with the P-GMRF model. S, ∀x j ∈ Ω:

We emphasize that the choice of x s 1 = x j is purely arbitrary and made to illustrate the fact that one of the variables in d j is x j , ∀d j ∈ D j .

Remark: Due to the sparsity of the connections in SBNs, the summations induced by the expectations can be computed exactly. However, in denser graphical models, some additional approximations must be introduced for the tractability of these expectations. These approximations consider a modified objective function which constitutes an upper bound to the KL divergence to minimize. Such procedures to obtain a coarser but more tractable objective function often consist in the simplification (e.g. linearization) of the term log p(x x x x x x x). See for example [START_REF] Wang | Markov random field modeling, inference & learning in computer vision & image understanding: A survey[END_REF] or [START_REF] Depraetere | A comparison of variational approximations for fast inference in mixed logit models[END_REF].

Finally the MF inference in SBN is done according to the following steps:

1. Initialize q for all factors.

2. ∀j ∈ S, update q j with the expression that maximizes Equation 3.40, i.e., Equation D.18.

3. Check convergence and repeat step 2 if needed.

Following this procedure, the divergence is proved to decrease at each update of q j , ∀j ∈ S (Wiegerinck, 2000).

D.2. Markov Tree Variational Inference in SBNs

This section gives details to solve the optimization problem given in Equation 3.42. We first isolate one of the factors, r cj of the variational distribution r. This leads to splitting sums and products according to the different clusters c i , ∀i ∈ S. It follows:

Similarly we find:

where Z is a normalization constant. Again, further simplifications can be made straightforwardly.

Remark: We can conclude this section recalling the remark at the end of the section on MT VI: in the final update equations for STMT VI we see more dependencies on the distribution t itself than in Equation D.28. Indeed this comes from the fact that t, which follows a STMT distribution, has more direct dependencies between its random variables than r which follows a MT distribution.

Appendix E.

Complements on fcCRFs

E.1. Augmented space for the bilateral kernel computations

In this section we demonstrate how the term 4.8) can be interpreted as a convolution so that it can be approximated with the approach seen in Section 4.3.2. Recall that k 1 is defined in Equation 4.5 and is called the bilateral kernel. It is the product of a Gaussian on the spatial distance and a Gaussian on the pixel value difference [START_REF] Paris | Bilateral filtering: Theory and applications[END_REF]. The concept of augmented space has been defined in [START_REF] Chakfé | Editor's Choice-European Society for Vascular Surgery (ESVS) 2020 Clinical Practice Guidelines on the Management of Vascular Graft and Endograft Infections[END_REF] based on (Paris and F. Durand, 2006). We now study how the derivations from the latter adapt to our case.

We have, ∀x ′ ∈ Ω, ∀j ∈ S:

where G θα (|ss ′ |) and G θ β (|y sy s ′ |) are the two Gaussian kernels respectively on the spatial distance and pixel value difference. The augmented space technique now appears by summing over an auxiliary variable defined in the space of the pixel intensities I 1 . The so-called augmented space is the product space S × I. It follows, with δ the Kronecker function:

We can see that Ġθαθ β = G θα (|j -i|)G θ β (|y j -ξ|) defines a Gaussian kernel on S × I. Then, ∀x ′ ∈ Ω, define the augmented space variational distribution q(x ′ ) with, ∀(i, ξ) ∈ S × I:

It follows that we can finally write Equation E.1 as a convolution in the augmented space:

As an illustration, if I is an unidimensional space (grayscale image) and S the set of vertices of a d-dimensional image, the augmented space variational distribution q(x ′ ),∀x ′ ∈ Ω has d + 1 finite dimensions. The data cube it forms is sparse, it is null everywhere except at one point of the augmented axis I, for each site i of the image in S, where it takes the value q i (x ′ ). We then have:

(E.5)

E.2. Markov Chain Variational Inference in fcCRFs

In this section we derive the equations to solve the optimization problem of Equation 4.11. In what follows, N is the number of parallel MCs and M is their length. The initial steps of the derivation are similar to the Markov Tree VI update equations in Appendix D.2. However, we need to consider the visible random variables. We get, ∀n ∈ {1, . . . , N }, ∀m ∈ {2, . . . , M } 2 and ∀(x n m , x n m-1 ) ∈ Ω 2 : marginal distributions that are used. These approximations are obvious for computational reasons.

We now focus on the term A of Equation E.7 which requires a summation on all the sites: Appendix F.

Complements on the applications to vascular surgery F.1. Fine stent segmentation: exponential mixture models

In this section we show that a mixture of exponential (MoE) distributions fits better the observed realizations of random variables (the image) that we want to segment. Works such as [START_REF] Delignon | Estimation of generalized mixtures and its application in image segmentation[END_REF] [START_REF] Pieczynski | Estimation of generalized mixture in the case of correlated sensors[END_REF] and [START_REF] Monfrini | Estimation de mélanges généralisés dans les arbres de Markov cachés, application à la segmentation des images de cartons d'orgue de barbarie[END_REF] involve the use of the exponential distribution in Hidden Markov Models and motivated the work we present here.

We will use the BIC score [START_REF] Wit | All models are wrong...': an introduction to model uncertainty[END_REF] as a way to find the best fitting mixture to our data. The BIC score for a model is defined as:

List of Tables

List of Algorithms A.4. Expectation Maximization (Dempster et al., 1977) . . . . . . . . A.5. Stochastic Expectation Maximization [START_REF] Miguel A Carreira-Perpiñan | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF] . . . . . . . A.6. Marroquin algorithm for the MPM computation [START_REF] Loeliger | Probabilistic solution of ill-posed problems in computational vision[END_REF] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1. Construction of q q q q q q q . . . . . . . . . . . . . . . . . . . . . . . . . . B.2. Sampling a GMRF with via the Fourier transform properties . .