Keywords: volume mapping and context specific methods deformable image registration, markov random field, slice-to-volume registration, graphical models, discrete optimization i . . . . . . . . . . . . . . .

Titre : Recalage déformable à base de graphes: mise en correspondance coupe-vers-volume et méthodes contextuelles Mots clés : recalage déformable, champs aléatoires de Markov, recalage slice-to-volume, modèles graphiques, optimisation discrète Résumé : Les méthodes de recalage d'images, qui ont pour but l'alignement de deux ou plusieurs images dans un même système de coordonnées, sont parmi les algorithmes les plus anciens et les plus utilisés en vision par ordinateur. Les méthodes de recalage servent à établir des correspondances entre des images (prises à des moments différents, par différents senseurs ou avec différentes perspectives), lesquelles ne sont pas évidentes pour l'oeil humain. Un type particulier d'algorithme de recalage, connu comme « les méthodes de recalage déformables à l'aide de modèles graphiques » est devenu de plus en plus populaire ces dernières années, grâce à sa robustesse, sa scalabilité, son efficacité et sa simplicité théorique. La gamme des problèmes auxquels ce type d'algorithme peut être adapté est particulièrement vaste. Dans ce travail de thèse, nous proposons plusieurs extensions à la théorie de recalage déformable à l'aide de modèles graphiques, en explorant de nouvelles applications et en développant des contributions méthodologiques originales.

Notre première contribution est une extension du cadre du recalage à l'aide de graphes, en abordant le problème très complexe du recalage Les deux autres contributions inclues dans ce travail de thèse, sont liées au problème de l'intégration de l'information sémantique dans la procédure de recalage (indépendamment de la dimensionnalité des images). Actuellement, la plupart des méthodes comprennent une seule fonction métrique pour expliquer la similarité entre l'image source et l'image cible. Nous soutenons que l'intégration des informations sémantiques pour guider la procédure de recalage pourra encore améliorer la précision des résultats, en particulier en présence d'étiquettes sémantiques faisant du recalage un problème spécifique adapté à chaque domaine.

Nous considérons un premier scénario en proposant un classificateur pour inférer des cartes de probabilité pour les différentes structures anatomiques dans les images d'entrée. Notre méthode vise à recaler et segmenter un ensemble d'images d'entrée simultanément, en intégrant cette information dans la formulation de l'énergie. L'idée principale est d'utiliser ces cartes estimées des étiquettes sémantiques (fournie par un classificateur arbitraire) comme un substitut pour les données non-étiquettées, et les combiner avec le recalage déformable pour améliorer l'alignement ainsi que la segmentation. Notre dernière contribution vise également à intégrer l'information sémantique pour la procédure de recalage, mais dans un scénario différent. Dans ce cas, au lieu de supposer que nous avons des classificateurs arbitraires pré-entraînés à notre disposition, nous considérons un ensemble d'annotations précis (vérité terrain) pour une variété de structures anatomiques. Nous présentons une contribution méthodologique qui vise à l'apprentissage des critères correspondants au contexte spécifique comme une agrégation des mesures de similarité standard à partir des données annotées, en utilisant une adaptation de l'algorithme « Latent Structured Support Vector Machine ».

Image registration methods, which aim at aligning two or more images into one coordinate system, are among the oldest and most widely used algorithms in computer vision. Registration methods serve to establish correspondence relationships among images (captured at different times, from different sensors or from different viewpoints) which are not obvious for the human eye. A particular type of registration algorithm, known as graphbased deformable registration methods, has become popular during the last decade given its robustness, scalability, efficiency and theoretical simplicity. The range of problems to which it can be adapted is particularly broad. In this thesis, we propose several extensions to the graph-based deformable registration theory, by exploring new application scenarios and developing novel methodological contributions.

Our first contribution is an extension of the graph-based deformable registration framework, dealing with the challenging sliceto-volume registration problem. Slice-tovolume registration aims at registering a 2D image within a 3D volume, i.e. we seek a mapping function which optimally maps a tomographic slice to the 3D coordinate space of a given volume. We introduce a scalable, modular and flexible formulation accommodating low-rank and high order terms, which simultaneously selects the plane and estimates the in-plane deformation through a single shot optimization approach. The proposed framework is instantiated into different variants based on different graph topology, label space definition and energy construction. Simulated and real-data in the context of ultrasound and magnetic resonance registration (where both framework instantiations as well as different optimization strategies are considered) demonstrate the potentials of our method.

The other two contributions included in this thesis are related to how semantic information can be encompassed within the registration process (independently of the dimensionality of the images). Currently, most of the methods rely on a single metric function explaining the similarity between the source and target images. We argue that incorporating semantic information to guide the registration process will further improve the accuracy of the results, particularly in the presence of semantic labels making the registration a domain specific problem.

We consider a first scenario where we are given a classifier inferring probability maps for different anatomical structures in the input images. Our method seeks to simultaneously register and segment a set of input images, incorporating this information within the energy formulation. The main idea is to use these estimated maps of semantic labels (provided by an arbitrary classifier) as a surrogate for unlabeled data, and combine them with population deformable registration to improve both alignment and segmentation.

Our last contribution also aims at incorporating semantic information to the registration process, but in a different scenario. In this case, instead of supposing that we have pre-trained arbitrary classifiers at our disposal, we are given a set of accurate ground truth annotations for a variety of anatomical structures. We present a methodological contribution that aims at learning context specific matching criteria as an aggregation of standard similarity measures from the aforementioned annotated data, using an adapted version of the latent structured support vector machine (LSSVM) framework.

First and foremost, I would like to thank Prof. Nikos Paragios for trusting and giving me the opportunity to do my PhD in the Center for Visual Computing of CentraleSupelec. His ideas are the basis of this thesis. His insights and advice together with our discussions guided my journey through the long path towards my PhD. But more importantly, his support after occasional failures, his jokes and constant stimulation allowed me to keep on going and do not give up in the difficult moments.

I am extremely grateful to the members of the jury. To Julia Schabel and Sandy Wells for reviewing this manuscript, and Pierre Jannin, Nicholas Ayache and Christos Davatzikos for reading this thesis and accepting to be part of the jury. It is a great honor to have them in my jury.

Many thanks to Prof. Silvio Savaresse for receiving me in Stanford during three amazing months, to Roland Angst for his valuable advises and discussions during this period, and to all the guys from CVGL lab for the nice moments we shared in the Bay Area.

I would like to thank the French State for the financial support for this thesis, and the Argentinian State for allowing me to study during 6 years in our public universities before moving to France, and for keeping the Argentinian House in Paris (where I lived during three amazing years) alive during more than eight decades.

The fantastic four years I spent in France would have not been the same without all my friends, which are too many to be named but too important to be omitted. The new ones that I met in the lab and in the Argentinian House at Cite Universitaire, and those that were always there supporting me from Argentina. Without them, my days in Paris would not be as happy. Special thanks to Puneet, Evgenios, Ishan and Aline, for being the best office mates one could ever ask for, to Natalia and Carine for always helping me with the paper work, to Vivien for the translation of the thesis abstract to French and to Fede for his help with the thesis writing until the last minute.

Finally, I would like to thank my family for their moral support; specially to my parents for their daily call, which made me feel in my hometown when I was actually thousands of kilometers away from home. To my Dad, for his never-ending optimism towards me and his wise words, always pronounced in the right moment. And to my Mum, because perhaps unconsciously, she was the one who transferred me her passion for knowledge, teaching and science. The same passion that brought me into the amazing world of research.

List of Figures

Introduction

Image registration algorithms aim to align two or more images into the same coordinate system and they are widely used in computer vision [Szeliski 2006]. Inspired by Horn and Shunk [Horn 1980] and the work of Lucas and Kanade 1 in the eighties [START_REF] Lucas | [END_REF]], the research community of computer vision has made major efforts towards developing more accurate and efficient strategies. Registration methods seek to establish correspondence relationships among images (captured at different times, from different sensors or from different viewpoints [Brown 1992]) which are not obvious for the human eye. The important image variations in conjunction with their diverse origin, make image registration one of the most challenging tasks in the field. Consequently, we observe the absence of a general algorithm suited to solve the general registration problem [Zitová 2003], and research has still to be done in order to address it. The literature in the area is vast and, every year, a significant number of new papers relevant to this problem are published in the most prestigious conferences and journals. Graph-based deformable registration methods [START_REF] Glocker | [END_REF]] have become popular during the last decade due to the fact that they exhibit robustness, scalability, efficiency and theoretical simplicity. The range of problems to which it can be adapted is broad. In this thesis, we propose several extensions to the graph-based deformable registration theory, by exploring new application scenarios and introducing novel methodological contributions. Before presenting our work, let us briefly motivate the research on image registration.

Motivation

Image registration is a central component for a large number of visual perception problems. The classical optical flow estimation problem [Fortun 2015], for example, can be formulated as a registration problem where we aim at recovering the displacement of intensity patterns in two different images. Motion analysis, video restoration, vehicle navigation systems and video surveillance, just to name a few, are cases in which optical flow estimation and image registration are critical. Another application involves remote sensing images. Accurate registration algorithms are essential in supporting Earth scientists as they mosaic satellite images and track changes of the planet's surface over time for environmental, political and fundamental science studies [Le [START_REF] Moigne | [END_REF]. In the context of medical image analysis (main 4 Chapter 1. Introduction interest of this thesis), image registration methods are ubiquitous: diagnostic, planning, surgical and radiotherapy procedures are a few examples where these methods are essential.

The information provided by different anatomical imaging modalities (such as computed tomography (CT), magnetic resonance imaging (MRI), ultrasound (US), X-ray, video sequences obtained by laparoscopy or endoscopy, etc.) as well as functional modalities (single-photon emission computed tomography (SPECT), positron emission tomography (PET), functional MRI (fMRI), etc.) is of complementary nature. Therefore, proper integration or fusion of data obtained from the independent images is desirable [Maintz 1998]. Image registration plays a key role in this task, since it allows bringing different image modalities into spatial correspondence. Longitudinal studies [Csapo 2012] -where images are gathered for the same subjects repeatedly over time-require image registration to monitor tumor growth, post surgery control, etc. Motion correction when acquiring 3D images of moving objects (such as fetal brain imaging [START_REF] Rousseau | [END_REF], Jiang 2007a] or fMRI reconstruction [Kim 1999]) is another problem tackled by image registration, where acquired images can once registered compensate the motion that occurred during the their acquisition. Another example is image guided surgery, where pre-operative 3D volumes and intra-operative 2D images need to be fused to guide surgeons during medical interventions. Image registration is fundamental in this process: it allows doctors to navigate 3D pre-operative high-resolution annotated data using low-resolution 2D images acquired in real-time during surgery. Several examples can be cited in this case such as image guided surgeries and therapies [START_REF] Fei | [END_REF]], biopsies [Xu 2014a], radio frequency ablation [Xu 2013], tracking of particular organs [START_REF] Gill | [END_REF]] and minimally-invasive procedures [START_REF] Liao | [END_REF], Huang 2009].

Graph-based Deformable Registration

Image registration methods can be classified according to different criteria, such as the transformation model (used to establish spatial correspondences), the matching criterion (means of defining similarities between images) and the optimization method (used to minimize the aforementioned matching criterion and thus infer the optimal transformation model). Given the importance of this problem, a large number of techniques has been proposed. We refer to the work by [Maintz 1998,Sotiras 2013] presenting a comprehensive study about medical image registration.

An active research approach to address deformable registration during the last years, has been the use of discrete methods and graphical models [START_REF] Glocker | [END_REF]]. The term deformable is used to denote the fact that the observed images are associated through a nonlinear dense transformation, or a spatially varying deformation model [Sotiras 2013]. Discrete methods and graphical models [Wang 2013, Paragios 2014] are powerful formalisms that have been successfully applied to many vision tasks including image registration. Casting image registration as a labeling problem offers a simple yet robust and efficient way to model this problem. In such a formulation, the graphical model consists of nodes corresponding to a deformation grid and the connectivity is used to impose regularization constraints. These nodes are labeled with displacement vectors corresponding to a discrete search space of allowable transformations.

Such a formulation, introduced in [START_REF] Glocker | [END_REF]] is modular with respect to the matching criterion and can be efficiently optimized through different discrete optimization methods (see Section 2.4). Glocker proposed a formulation to deal with registration of images of the same dimension (i.e. 3D to 3D or 2D to 2D registration). Numerous extensions to this work were then introduced, proving the robustness of the framework. Just to name a few, let us consider: (i) the work of [Sotiras 2009] extending this framework to perform population registration (i.e. to identify an homology between more than two images at the same time), (ii) the work from the same authors coupling geometric and iconic registration method [Sotiras 2010a], where the matching criterion is modified so that it accounts for landmark correspondences and intensity based similarity measures, (iii) multi-atlas segmentation using graph-based deformable registration was proposed by [Alchatzidis 2014a[START_REF] Alchatzidis | [END_REF], while in (iv) [Fecamp 2015a,Fecamp 2015b] authors focused on how linear and deformable parameters can be decoupled when performing graph based registration. In another line of research, regularization constraints based on learned deformation priors were incorporated to the framework [Glocker 2009b], while an adaptive co-registration and segmentation of brain tumors where the sampling resolution of the label space is adapted according to the image content is proposed in [START_REF] Parisot | [END_REF]]. Last but not least, a modified version of the initial work was also applied to detect temporal changes in remote sensing images [START_REF] Vakalopoulou | [END_REF]].

The number of extensions proposed to the graph-based deformable registration framework, demonstrates how flexible and powerful it can be. Let us now present the main ideas developed in this thesis, which contemplate new application scenarios and introduce a number of methodological contributions.

Contributions

We begin this thesis by including, in chapter 2, an initial overview concerning discrete models in computer vision. We derive the standard semantic segmentation and image registration problems within this framework, and we present the inference algorithms used in this dissertation to optimize such models.

The first class of contributions is presented in chapter 3. We introduce a generic graph-based deformable registration framework, to deal with the challenging sliceto-volume registration problem. As we mentioned, the initial work by [START_REF] Glocker | [END_REF]] proposes a formulation for deformable registration of images in dimensional correspondence. However, in slice-to-volume registration, we aim at registering a 2D image within a 3D volume, i.e. we seek a mapping function which optimally maps a tomographic slice I to the 3D coordinate space of a volume J. This problem is substantially different from the 3D-3D or 2D-2D case; therefore, the standard registration approach cannot be directly extrapolated to the slice-to-volume scenario. We instantiate this generic framework into three different models to solve this challenging problem. Our aim is to demonstrate how flexible and powerful the graph theory is in terms of expressive potential of the modeling process, while solving a new problem using graphical models. This general formulation can be expressed through different discrete labeling problems on a graph by changing its topology, label space and energy terms.

First, we propose the so-called overparameterized method, which combines linear and deformable parameters within a coupled formulation on a 5-dimensional label space [Ferrante 2013]. We adopt a metric free locally over-parametrized graphical model seeking simultaneously a linear plane transformation and an in-plane dense deformation. Image similarity is encoded in singleton terms, while geometric linear consistency of the solution (common/single plane) and in-plane deformations smoothness are modeled in a pair-wise term. The main advantage of such a model is the simplicity provided by its pairwise structure, while the main disadvantage is the dimensionality of the label space which makes inference computationally inefficient and approximate (limited sampling of search space).

Motivated by the work of [START_REF] Shekhovtsov | [END_REF]], we present an alternative model referred as the decoupled approach. In this case, linear and deformable parameters are separated into two interconnected subgraphs which refer to lower dimensional label spaces [Ferrante 2015b]. The main contribution is the ability to decouple the plane selection and the in-plane deformation parts of the transformation towards reducing the complexity of the model, while being able to obtain simultaneously the solution for both of them. To this end, the plane selection process is expressed as a local graph-labeling problem endowed with planarity satisfaction constraints, which is then directly linked with the deformable part through the data registration likelihoods. The resulting model is modular with respect to the image metric, can cope with arbitrary in-plane regularization terms and inherits excellent properties in terms of computational efficiency.

The last slice-to-volume registration model presented in this work is the highorder approach [Ferrante 2015a]. Differently from the strategy adopted for the decoupled model, here we aim at reducing the dimensionality of the label space by augmenting the order of the graphical model. Using third-order cliques which exploit the expression power of this type of variable interactions, we simplify the label space to standard 3D displacement vectors. Geometrical consistency (unique plane selection) and deformation smoothness (in-plane deformations) as well as image similarity are encoded in different third order cliques. Such a model provides better satisfaction of the global deformation constraints at the expense of quite challenging inference.

The second class of contributions is presented in chapter 4. The objective is to introduce semantic information within the registration process (independently of the dimensionality of the images). The notion of semantic segmentation refers to either partial or full labeling of the source/target volumes. Currently, most of the methods rely on a single metric function that explains the similarity between the source and target images. We argue that incorporating semantic information to guide the registration process, will further improve the accuracy of the results, in particular, in the presence of semantic labels which then make the registration a domain specific problem.

In section 4.2, we consider the scenario where we are given a classifier inferring probability maps for the anatomical structures being present in the images. The method seeks to simultaneously register and segment a set of input images, incorporating this information within the energy formulation. The central idea is to use these estimated maps of semantic labels (provided by an arbitrary classifier) as a surrogate for unlabeled data, and combine them with deformable registration to improve both alignment and segmentation. To this end, images are deformed towards a common space where consensus with respect to image similarities (according to a given arbitrary metric) and semantic labels is achieved. We employ an approach sharing conceptual similarities with α-expansion [START_REF] Boykov | [END_REF]], combined with efficient linear programming optimization methods. Promising results demonstrate the potential of our method on two different datasets, containing annotations of challenging brain structures.

Our last contribution -presented in section 4.3-also exploits semantic information in the registration process with a rather different objective. In this case, instead of supposing that we have more or less accurate classifiers at our disposal, we are given a set of ground truth annotations for different anatomical structures. We present a method that aims at learning context-specific matching criteria as a weighted aggregation of standard similarity measures from the aforementioned annotated data. As we will show in section 2.3, graph-based registration is modular with respect to the matching criterion. In other words, alternative similarity measures can be used depending on the type of images we are registering. Even so, the definition of this metric function is among the most critical and complicated aspects of the registration process. To alleviate this task, we show how different content specific metrics can be combined towards defining more efficient matching criteria. To this end, we adapted the standard latent structured support vector machine (lssvm) framework [START_REF] Yu | [END_REF]] so that optimal combinations of standard similarity measures can be learned from examples.

We conclude this work in chapter 5, presenting some general conclusions and giving further ideas and future research directions.

Chapter 2

Discrete Models in Computer Vision

In this chapter we provide an overview as it concerns discrete models in computer vision. First, we describe how classical discrete optimization and graphical models theory can be used to cast vision tasks as a graph labeling problem, giving a simple example where the semantic segmentation problem is solved within this framework. Second, following [START_REF] Glocker | [END_REF][START_REF] Glocker | [END_REF], we focus on how deformable image registration can be formulated using discrete models. In an effort to make this thesis self-contained, we also describe some standard discrete inference algorithms, applied in various parts of our work to perform the optimization task. We conclude with a small discussion about the relevance of discrete methods when modeling vision tasks.

Graphical Models in Computer Vision

Graphical models have been part of the computer vision field for more than three decades. In late nineties, they started to gain more attention in the community, thanks to the advent of efficient inference methods. They are also known as probabilistic graphical models (PGMs) [START_REF] Koller | [END_REF]], since they encode probability distributions over high-dimensional spaces, and can be efficiently manipulated using ideas from discrete optimization and computer sciences. They enjoy a combination of simplicity -in mathematical terms-and efficiency that makes them the choice of preference for a wide spectrum of vision problems. Low, mid and high-level vision tasks can be modeled within this framework. Visual perception is addressed by defining a specific parametric model, associating it to the available observations (images) through an objective function and optimizing the model parameters given both, the objective and the observations [Paragios 2014]. Markov Random Fields (MRF) are a particular class of graphical models, ubiquitous in computer vision research [Wang 2013]. This powerful modeling tool corresponds to an undirected graph G encoding a probability distribution governed by the local independence assumption (also known as local Markov or markovian property). Each node of this graph represents a variable from the probability distribution. Dependence among the variables is modeled using the graph neighbor system, meaning that two variables are dependent of each other if and only if there is an edge between them. This dependency can be extended from order two to any order, through the concept of clique. A clique is defined as set of fully connected nodes; therefore, all the variables represented by nodes in a given clique depend on each other. The Markov property imposes that a node is independent of any other node given all its neighbors. This holds for any distribution that can be represented by a MRF.

Formally, a MRF is an undirected graph G = V, C , where V is a set of nodes (every node i ∈ V is associated to a variable x i ∈ X ) and C is the set of cliques defining the graph connectivity. The Hammersley-Clifford [START_REF] Hammersley | [END_REF]] theorem states that the probability distribution of any MRF can be factorized over its cliques as follows:

p(x) = 1 Z c∈C ψ c (x c ), x = {x c }, (2.1)
where Z is the so-called partition function (used to normalize the probability distribution so that x p(x) = 1), ψ c is the positive potential function associated to clique c and x c is a possible configuration (assignment) for the variables x i being part of the clique. Note that this definition is general enough to accommodate discrete and continuous variables. In this work, we focus on the discrete case, where variables x i associated to nodes i ∈ V take values l from a discrete label space L. This is known as a discrete MRF where x is a labeling (or assignment or realization) that imputes a label l i ∈ L to every variable x i ∈ X , x c ⊂ x is a subset of x that assigns a label to every variable in the clique c and ψ c is a clique potential mapping a positive real value to every possible configuration x c .

Graphical Models in Computer Vision

In this framework, nodes represent the variables of interest while domain specific/prior information is expressed through the clique potentials. Once the objective function is defined, the optimal solution can be obtained by performing inference via Maximum a posteriori (MAP) estimation (i.e. the solution x will be the one that maximizes the probability function p):

x = arg max x∈X p(x) = arg max x∈X 1 Z c∈C ψ c (x c ). (2.2)
A more convenient representation of the joint distribution p(x) can be obtained, so that the objective function is defined as a sum of clique potentials instead of a product. To this end, let us define a new type of potentials φ such that φ c (x c ) = -log ψ c (x c ). Therefore, the probability distribution p(x) is redefined in terms of φ c as:

p(x) = 1 Z exp{- c∈C φ c (x c )}, x = {x c }. (2.3)
Let us denote E(x; Ψ), the energy of a MRF, as the sum of the clique potentials Φ = {φ c }, ∀c ∈ C given by:

E(x; Φ) = c∈C φ c (x c ).
(2.4)

We reformulate the aforementioned MAP inference problem in a discrete MRF from equation 2.2 as an equivalent minimization problem where the objective function is the MRF energy E(x; Φ) as follows:

x = arg min x∈X E(x; Φ) = arg min x∈X c∈C φ c (x c ) (2.5)
Note that this equivalent energy minimization problem, does not compute the actual probability function p(x). Therefore, the partition function Z (frequently defined as the sum of all the potentials evaluated in all possible labelings) is not required. This is desirable since computing Z as a normalizing factor for p(x) requires the evaluation of the potentials for all possible labelings, usually intractable considering its definition.

The underlying graph can involve cliques of different sizes. The size of the maximum clique indicates the order of the MRF, which is one of its most distinctive properties. MRFs of order two are known as pairwise, whereas those whose order is higher than two are referred as high-order. In the same way, the energy terms of equation 2.5 involving two variables are called pairwise terms, while those with more than two variables are the high-order terms. The terms associated to only one variable are the unary terms. It is worth mentioning that the vast majority of the graph-based vision methods exploit pairwise models. The reason is that exact or approximate efficient inference solutions can be obtained in these cases. The quality of the inference depends on several facts such as the cardinality of the label space or different mathematical properties of the energy function [START_REF] Kolmogorov | [END_REF]]. Nonetheless, during the last few years, novel high-order models and inference algorithms have been developed leading to more accurate solutions [Kohli 2012[START_REF] Komodakis | [END_REF], Ishikawa 2009]. In chapter 3, we explore the expressive power of such models in the context slice-to-volume deformable registration.

The cardinality of the label space defines either a binary problem (two labels) or a multi-label problem (more than two labels). Kolmogorov and coworkers [Kolmogorov 2004] showed that, for a particular type of binary problems where the energy fulfills a condition known as submodularity, the global minimum can be found in polynomial time using graph-cuts [START_REF] Boykov | [END_REF]].

An alternative structure frequently used to represent graphical models is factor graph. Some inference libraries, such as OpenGM2 [Kappes 2013] (used in this thesis), model the MAP problem using factor graphs since they explicitly describe the factorization of the joint distribution in the graph. It leads to a clear conceptualization of the underlying factorization. For the sake of completeness, we include in the appendix (section A.1) a brief description of this representation.

Several approaches were proposed to perform MAP inference on graphs, depending on the structure, order and energy of the discrete MRF. We refer the reader to the comprehensive survey [Kappes 2013] where a comparative study of modern inference techniques is presented. For demonstration purposes, in the remaining of this chapter we present two generic formulations of graphical models towards addressing two fundamental problems in computer vision and medical imaging. The semantic segmentation problem is used as an introductory and simple example. Then, image registration (the problem that we tackle in this thesis) is formulated within the same framework.

Semantic Segmentation as a Discrete Labeling Problem

Semantic segmentation is a typical example illustrating how a vision task can be solved through MAP inference on a discrete MRF. Given an image I(i) : Ω ∈ 2 → , we would like to assign a semantic class to every pixel i ∈ Ω. Let us consider the case where we have a map indicating the probability of every pixel to be part of a given semantic class l ∈ L (label). This map could be obtained using any classifier. This probability map is represented by function P i (l) : L → [0, 1], which is a prior of the semantic classes, where i indicates a pixel from the original image.

Let G = V, E be a graph representing a MRF, where nodes in V are associated to variables (pixels) and E is a standard 4-connected neighborhood system defining a two-dimensional lattice (grid). In such a simple 4-connected neighborhood system each node is connected to its adjacent neighbors in left, right, top, and bottom directions. Variables x i associated to nodes i ∈ V can take labels l i from a label space L. A labeling S = {l i | i ∈ V} assigns one label to every variable. We define the energy E(S; Φ) consisting of unary potentials V i ∈ Φ and pairwise potentials V ij ∈ Φ such that it is minimum when S corresponds to the best labeling.

Deformable Image Registration as a Discrete Labeling Problem 13

Unary terms V i (l i ) = -log(P i (l i )) assign low energy to high probability values. Pairwise terms encode an spatial homogeneity constraint by simply encouraging neighbor variables to take the same semantic label V i,j (l i , l j ) = [l i = l j ]. This type of pairwise term is known as Potts model, and it is defined as:

[l i = l j ] =    1, if l i = l j 0, otherwise.
(2.6)

Finally, the energy minimization problem is:

S * = arg min E(S; Φ) = arg min i∈V V i (l i ) + λ (i,j)∈E V i,j (l i , l j ).
(2.7)

S * represents the optimal label assignment. Note that this energy is a metric in the space of labels L. According to [START_REF] Boykov | [END_REF]], V is called a metric in the space of labels L if, given any three labels α, β, γ ∈ L it satisfies:

V (α, β) = 0 ⇔ α = β, (2.8) V (α, β) = V (β, α) ≥ 0, (2.9) V (α, β) ≤ V (α, γ) + V (γ, β).
(2.10) Thus, if the energy fulfills these conditions, it is guaranteed that by using αexpansion technique we can find a solution Ŝ whose energy lies within a factor of 2 with respect to the optimal energy (i.e. E( Ŝ) ≤ 2.E(S * )). α-expansion is a well known move-making technique to perform approximate inference using graphcuts [START_REF] Boykov | [END_REF]]. It proved to be accurate in a broad range of vision problems. We refer the reader to [START_REF] Boykov | [END_REF]] for a complete discussion on energy minimization using α-expansion.

Deformable Image Registration as a Discrete Labeling Problem

In the previous section, we described how semantic segmentation can be modeled through MAP inference on a discrete MRF. This problem is inherently discrete since the solution is given by a finite set of discrete labels indicating the class that corresponds to every pixel of the input image. However, in other cases like deformable image registration, the nature of the problem is not discrete. Thus, modeling it through discrete methods is not an obvious task. In this section, following [START_REF] Glocker | [END_REF], Glocker 2009c[START_REF] Glocker | [END_REF], we describe how deformable image registration can be intelligently formulated as a discrete labeling problem. Deformable image registration is an NP-complete problem [START_REF] Keysers | [END_REF]] and consists in recovering a transformation T : d → d mapping a source image I Figure 2.1: (a) Overlapping between the source (I) and target (J) images, before registration. As it can be observed, the aligment bewteen the images is not perfect. (b) Overlapping between the deformed source image and the target, after registration. The deformable registration algorithm recovered a deformation field T(x) that was applied to the source image, in order to align them.

to a target image J (I, J : Ω ⊂ d → ), so that optimal alignment is achieved1 (the warped image is denoted by I • T ). Figure 2.1 illustrates this idea. This can be mathematically formulated as follows: (2.11) where M represents the data term (similarity between the warped source image I • T and the target image J) and R the regularization term (geometric consistency/realistic deformations). Here T indicates the optimal transformation. T is usually defined using a dense deformation field D : d → d (relative with respect to the current position), mapping every position x ∈ Ω from image I to its corresponding position T (x) in image J:

T = arg min T M(I, J; T ) + R(T ),
T (x) = x + D(x).
(2.12)

In such a setting, the dense deformation field D has as many parameters as pixels in the image (one d-dimensional displacement vector for every x ∈ Ω). Such dimensionality of the graph results on high complexity which should be avoided (e.g. for a standard volumetric image of 256 × 256 × 256 voxels, the degrees of freedom of our parameteric model would be 256 3 = 16777216 3-dimensional vectors). Therefore, a sparse representation of the deformation field is adopted. It is a reparametrization which uses k control points distributed uniformily within

Deformable Image Registration as a Discrete Labeling Problem 15

the image domain to represent the complete deformation field. An interpolation strategy acting on the displacement of these control points is used to reconstruct the final solution. The simplification introduced by such a model is significant: we can interpolate the displacement vectors that correspond to every image position x by simply manipulating a few control points. Given a set of control point displacements D = {d 1 ..d k }, the dense deformation field D(x) used in equation 2.12 is reparameterized as: (2.13) where ω i are weighting (or interpolation) factors determining the influence of the control point in the interpolated final position (the closer the control point is to the current position x, the higher the influence of the displacement d i ). A common choice to define these weighting factors is the use of B-splines, resulting in the well-known Free Form Deformation (FFD) model. FFDs were popularized in the medical image registration community by [Rueckert 1999]. They adopt a regular grid as parametric model. Every control point contributes locally to the interpolation function. Other interpolation models like Thin-Plate Splines (see section 3.1.3.2) could be adopted. However, the local support for deformations and inherent computational efficiency make FFDs a good choice to deal with deformable image registration in a discrete scenario.

D(x) = k i=1 ω i (x)d i ,
Once the parametric model is chosen, let us proceed with a formal definition of the corresponding graphical model and the energy function associated with it. Let us consider an undirected graph G = V, E , where nodes are associated to variables representing control points, and E is a regular grid (2D or 3D depending on the dimensionality of the images) with a 4-connected neighborhood system in case of 2D registration, and a 6-connected neighborhood system in case of 3D registration (where edges connecting nodes along the z axis are added to the standard 4-connected system). The edges (i, j) ∈ E are used to impose smoothness on the deformation field and render the problem well posed.

The label space L consists of labels l i ∈ L associated to d-dimensional displacement vector d l i ∈ d (i.e. the integer label l i could be seen as an index value for the displacement vector d l i ). Several strategies can be used to discretize the space. Two common approaches are the so-called dense and sparse sampling (illustrated in figure 2.2). When we assign a label l i to a node i, the control point p i is moved according to the displacement vector d l i . A labeling Γ assigns one label l i to every node i ∈ V. Recall that every discrete label l i is associated in our framework to a d-dimensional displacement vectors d l i . Thus, we define D Γ as the set formed by the displacement vectors associated to the labeling Γ. Hence, by considering different labelings Γ, we explore a discretized version of the naturally continuous space of deformation fields (see example in Figure 2.3. For example, let us consider a simple 2D label space formed by 5 labels L = {0, 1, 2, 3, 4}. We assign to every label a vector of length 1 in each direction (according to the sparse sampling shown in figure 2.2), together with the null vector. Therefore, the vectors associated to every label would be:

d 0 = (0, 0), d 1 = (1, 0), d 2 = (0, 1), d 3 = (-1, 0), d 4 = (0, -1).
If we are dealing with a graph of 6 nodes, then the labeling Γ = {0, 3, 1, 1, 3, 2} is interpreted as the following control point based representation of the displacement field

D Γ = {d 0 , d 3 , d 1 , d 1 , d 3 , d 2 }.
The search across this space is guided by the minimization of the MRF energy defined as:

E(Γ; G, F ) = i∈V g i (l i ) + λ (i,j)∈E f i,j (l i , l j ), (2.14)
where Γ is a labeling assigning one label l i to every node i ∈ V, G are the unary potentials playing the role of data (or matching) terms and F are the pairwise terms acting as regularizers of the estimated deformation field. The MAP inference problem (a discretized version of the continuous formulation from equation 2.11, where G = {g i (•)} i∈V is associated to data terms M and F = {f ij (•, •)} (i,j)∈E to regularization terms R) is therefore defined as:

Γ = arg min Γ E(Γ; G, F ) = arg min Γ i∈V g i (l i ) + λ (i,j)∈E f i,j (l i , l j ).
(2.15)

The final dense deformation field D is interpolated from the set of displacement vectors D Γ associated to the parameters Γ, as previously explained.

The data term measures the matching quality (according to the similarity measure) between the deformed source image I and the target image J, given the displacement vectors. In this thesis, we only consider iconic matching criteria (i.e. image intensities based functions). However, other matching criteria could be considered (we refer the reader to the review by [Sotiras 2013] for complete discussion about this topic in the context of deformable medical image registration). Every pixel (or voxel) will contribute to the data term of a given control point, according to its distance. The control points have only a limited local image support. We consider a patch-based similarity measure where the region of support for control point p i is given by the image patch Ω i (note that, most frequently, the complete image domain Ω is covered by the union of patches Ω i , i.e. Ω = i Ω i ). Therefore, the unary potentials are defined as:

g i (l i ) = Ω i ω i (x)δ(I • T d l i (x), J(x))dx, (2.16)
where δ is an iconic similarity measure. Examples of standard similarity measures are the sum of absolute differences (SAD) and the sum of squared differences. In the context of different image modalities, the intensity correlation is not obvious and usually more complex functions based on statistical measures such as mutual information (MI) or cross-correlation (CC) may be adopted (for further information about iconic similarity measures, we refer the reader to section 3.1.2.1). The data term is the driving force of the optimization process. However, even if the matching criterion is optimally minimized, the resulting deformation field could be erroneous from geometrically view point. The regularization term plays a key role in this process, since it imposes geometrical and smoothness constraints towards anatomically realistic solutions. We define it as a pairwise term acting on the two labels (displacement vectors) assigned to neighbor nodes. Using pairwise terms, we can encode different smoothness constraints, such as the standard Pott's model (penalizes deviations in the displacement assignment), truncated absolute or quadratic difference (penalize the magnitude of the vectors difference up to a bound), approximated curvature penalty (approximates a second order derivative of the displacement field) or the distance preserving approach (penalizes changes in the distances between neighboring control points with respect to their initial position) [Glocker 2009a]. Let us illustrate, for example, the distance preserving approach, defined as the ratio between the current position of the control points p i , p j and their original position p o,i , p o,j :

Chapter 2. Discrete Models in Computer Vision

f i,j (l i , l j ) = || (p i + d l i ) -(p j + d l j ) || || (p o,i ) -(p o,j ) || .
(2.17)

As mentioned, alternative definitions could be considered for this term, depending on different constraints and domain specific restrictions.

A common strategy to improve the accuracy of the method and increase its capture range consists in using a pyramidal approach. It allows a progressive search space exploration, while keeping low complexity. The idea is to start with displacement vectors with important distance and progressively reduce their spread, centering them around the current solution. It is clear that increasing the number of labels to infinity will converge to a continuous formulation. Since this is intractable from computational perspective [START_REF] Glocker | [END_REF]], the idea behind this approach is to iteratively compose several deformation fields, while refining the grid of nodes and the label space, so that big as well as small local deformations are captured. At every iteration we consider a fixed label space and constant number of grid nodes, we perform MAP inference for equation 2.15 and we compose the resulting deformation field with the previous one.

MAP Estimation on Graphical Models

The computer vision tasks that were presented in the previous sections (i.e. semantic segmentation and deformable image registration) were formulated as MAP inference problems. In this section, we focus on how these MAP problems can be solved through discrete optimization.

Discrete optimization of MRFs is, in general, an NP-hard problem [START_REF] Shimony | Finding MAPs for belief networks is NPhard[END_REF]]. However, in special cases, it can benefit from very efficient solutions. The trivial brute force algorithm (i.e. trying all possible combination of labels for each and every variable) has an exponential complexity that makes such an approach unsuitable. More efficient algorithms have been developed during the last two decades which boosted the use of graphical models in the field of computer vision. They can be classified in three main categories according to [Kappes 2013 ii. Message passing methods, in which messages are calculated and propagated between nodes in a graph. This propagation can be seen as a re-parametrization of the original problem aiming to establish special properties in the re-weighted function that makes inference easier. Examples are the standard Loopy Belief Propagation (LBP) [Murphy 1999] and Three Re-weighted Belief Propagation (TRBP) [START_REF] Wainwright | [END_REF]].

iii. Max-flow and move-making algorithms make use of the well know max-flow min-cut [Boykov 2004] algorithm from graph theory, which can optimally solve some instances of discrete energies. These methods are usually combined with greedy strategies that iteratively minimize over the label space by solving a sequence of max-flow min-cut sub problems. Examples are α-expansion, αβ swap [START_REF] Boykov | [END_REF]] and FastPD [Komodakis 2007a] algorithms. Simpler movemaking algorithms not using max-flow, but still applying the strategy of starting with an initial labeling and iteratively moving to a better one until a convergence criterion is met, are the classical Iterated conditional modes (ICM) [Besag 1986] and its generalization Lazy Flipper [START_REF][END_REF]].

In the spirit of making this thesis self contained, in this coming section we describe three MAP inference algorithms used in our work. Loopy Belief Propagation and Lazy Flipper are used in chapter 3, since they can be applied to arbitrary graph topologies, label spaces and energy terms. In that chapter, we consider three different approaches to solve slice-to-volume deformable registration, varying in terms of graph topology, label space definition and energy construction. Therefore, we adopted these two algorithms which can be used to perform MAP inference in all these models. FastPD is used in chapter 4 where we propose context-specific pairwise deformable registration models which can be efficiently optimized through this method. In what follows, we briefly describe these three approaches.

Loopy Belief Propagation

Message passing algorithms estimate a solution by iteratively passing messages around the variables of the random field. These messages m ij (sent from a node i to a node j) are actually vectors of size | L | (cardinality of the label space), where every scalar entry represents what node i thinks about assigning label l to the node j (every entry is therefore m ij (l j )). Given, for example, a pairwise random field with unary potentials φ i and pairwise potentials φ ij , the message m ij (l j ) is computed as:

m ij (l j ) = min l i ∈L (φ i (l i ) + φ ij (l i , l j ) + r∈N (i)\{j} m rj (l i )),
(2.18)

where N (i) is the set of neighbors of i. Once a node (variable) i receives all the messages from its neighbors, it computes its beliefs (also vectors of size | L |) in a label l i as follows:

bel i (l i ) = φ i (l i ) + j∈N (i) m i j(li). (2.19)
The messages are iteratively passed from one node to its neighbors until no change occurs from one iteration to the next one. When convergence is achieved, the MAP labeling is obtained for every node i as the label l i that minimizes the belief bel i (l i ).

Note that both, messages and beliefs computed for a given node i, depend on the messages received from its neighbors j ∈ N (i). Therefore, we need a way to initialize the message calculation process. If the graph that underlies the MRF is a tree, this process is initialized in the roots since messages for these nodes can be calculated considering just the potentials (the set N (i) \ {j} from equation 2.18 is empty). In this case, at convergence, the solution is guaranteed to be optimal for arbitrary energies. If the structure is not a tree, messages are passed in arbitrary order, and the algorithm is not guaranteed to converge in a finite number of iterations. LBP has shown good experimental performance in empirical studies [Murphy 1999] with arbitrary graphs and it was successfully applied in this thesis to optimize different type of graphical models in 3.

Lazy Flipper

Lazy Flipper is a movie-making algorithm proposed by [START_REF][END_REF]]. It is a generalization of the well-known ICM [Besag 1986] which explores (exhaustively or not) the search space.

The idea is to start from an arbitrary initial assignment and perform successive flips of variables that reduce the local (and consequently the global) energy to be minimized. A greedy strategy is adopted to explore the space of solutions: as soon as a flip reducing the energy is found, the current configuration is updated accordingly. In a first stage, only one variable is flipped at a time (in the same way that ICM). However, once a configuration is found whose energy can no longer be reduced by flips of one variable, a new stage starts where all subsets of two connected variables (i.e. variables that are linked by an edge in the graph) are considered. This strategy is applied, considering sets of maximum size k. This parameter controls the search depth. For k = 1, it specializes to ICM. For bigger values of k a tradeoff between approximation quality and runtime is established, which in the limit converges to an exhaustive search over only the connected subgraphs (intractable in most of the cases). The algorithm is implemented using two efficient data structures: the connected subgraph tree (CS-tree, used to efficiently and uniquely enumerate the subgraphs so that only connected subsets of variables are considered) and a tag list that prevents the repeated assessment of unsuccessful flips.

The solutions obtained using Lazy Flipper are guaranteed to be optimal within a given Hamming distance in the label space. This move-making algorithm can be applied to problems modeled with arbitrary graphs, energies and label spaces. That is the reason why we choose such a basic but robust approach to compare different types of models in chapter 3. We refer the reader to [B Andres, JH Kappes 2010, Andres 2012] for a complete description and discussion about this method.

FastPD

FastPD is a discrete optimization algorithm based on principles from linear programming and primal dual strategies, which at the same time generalizes α-expansion [START_REF] Komodakis | [END_REF]]. Therefore, it may be included in categories (i) and (iii) from the classification scheme proposed at the beginning of this section. One of the main advantages of FastPD is its modularity/scalability, since it deals with a much wider class of problems than α-expansion, being an order of magnitude faster while providing the same optimality guarantees when performing metric labeling [Komodakis 2007b].

FastPD solves a series of max-flow min-cut problems on a graph. In that sense, it is similar to α-expansion which also performs MAP inference on multi-label problems by solving successive binary max-flow min-cut problems. The main difference between these approaches is the construction of the graph where max-flow min-cut algorithm is applied. α-expansion constructs the binary problem by restricting the label space, so that the only options for a given variable are to remain in its current assignment, or to take a label α (which varies in every iteration). Instead, FastPD constructs these binary problems by performing a Linear Programming Relaxation (LPR) of the integer program that represents the discrete MRF formulation. It builds upon principles drawn from the duality theory of linear programming, applying the well known primal-dual schema to the relaxed version of the MRF integer programming formulation.

For the sake of clarity, let us explicitly show how we can cast a standard MAP inference problem on a pairwise MRF as an integer program [Komodakis 2005]. Given, for example, the MAP formulation for the image registration problem presented in equation 2.15, we define the following integer program:

min i∈V α∈L g i (α)x i (α) + i∈V j∈N (i) αβ∈L f ij (α, β)x ij (α, β), (2.20) s.t. α∈L x i (α) = 1 ∀i ∈ V, (2.21) α∈L x ij (α, β) = x j (β) ∀β ∈ L, ∀i ∈ V ∧ j ∈ N (i), (2.22) β∈L x ij (α, β) = x i (α) ∀α ∈ L, ∀i ∈ V ∧ j ∈ N (i),
(2.23)

x i (•), x ij (•, •) ∈ {0, 1} (2.24)
In this formulation x i (•) and x ij (•, •) are indicators for the discrete labels l i assigned to each node i, used to linearize the MRF energy. A binary variable x i (α) is equal to 1 when node i chooses label α. The same holds for the pairwise case x ij (•, •). Equation 2.21 guarantees that every variable is labeled with only one label at a time, while equations 2.22 and 2.23 keep consistency between variables x i (•) and x j (•). This integer formulation is still an NP-hard problem, given the integrality constraint expressed in equation 2.24. Therefore, this integer program becomes a standard linear program when this constraint is relaxed, by allowing variables x i (•) and x ij (•, •) to take continuous positive values.

The aforementioned linear program task is used by FastPD as the primal problem, while the dual of the resulting LP is used as the FastPD dual problem. Then, the primal-dual schema is iteratively applied. At every iteration, a max-flow mincut problem is solved, improving both primal optimality and dual feasibility.

Given that FastPD can be applied to a wide range of pairwise energies (it only requires the pairwise potentials to be positive) providing good balance between efficiency and accuracy, it is used in chapter 4 where both models introducing context specific information to the registration process are conceived as pairwise MRFs.

Why Discrete Methods?

Discrete methods have several advantages when compared with continuous approaches for image registration. First, discrete algorithms are inherently gradientfree, while most part of continuous methods require the objective function to be differentiable. Gradient-free methods do not require computation of the energy derivative. Therefore, it may be applied to any complex energy function (allowing the user to define its own similarity measures in case of registration problems). The only requirement is that this function must be evaluable in a variety of possible discrete labelings. Second, continuous methods are quite often prone to be stuck in local minima when the functions are not convex. In case of discrete methods, even complicated functions could potentially be optimized using large neighbor search methods. The main limitation is the discretization of the continuous space; however, as suggested by [Glocker 2010], 'the optimality is bounded by the discretization, but with intelligent refinement strategy the accuracy of continuous methods can be achieved'. Third, parallel architectures can be used to perform non-sequential tasks required by several discrete algorithms (such as message calculation in LBP) leading to more efficient implementations. Fourth, by using a discrete label space we can explicitly control its range and resolution (it can be useful to introduce prior information, as it will be shown in chapter 3), while in continuous models it is not clear how this type of information can be used to constraint the solution. Last but not least, discrete frameworks such as discrete MRF provide a modular and principled way to combine prior knowledge with data likelihood, what makes it applicable to a wide range of vision tasks [Wang 2013].

Chapter 3

Slice-to-Volume Registration

In this chapter we investigate the use of graphical models in the context of slice-tovolume registration. We start with a comprehensive literature review, that presents most of the methods that have been proposed to deal with this challenging problem. Then, we introduce a scalable, modular and flexible formulation which can accommodate low-rank and high-order terms, that simultaneously selects the plane and estimates the in-plane deformation through a single shot optimization approach. The proposed framework is instantiated into different variants seeking either a compromise between computational efficiency (soft plane selection constraints and approximate definition of the data similarity terms through pairwise components) or exact definition of the data terms and the constraints on the plane selection. Simulated and real data in the context of US and MRI registration (where both framework instantiations as well as different optimization strategies are considered) demonstrate the potentials of our method.

Literature Review

Image registration is the process of aligning and combining data coming from more than one image source into a unique coordinate system. During the last decades, image registration has become one of the pillars of computer vision and medical imaging. In medical imaging, the problem of slice-to-volume registration has received further attention during the last decade. In this case, instead of registering images which same dimension, we seek to determine the slice (corresponding to an arbitrary plane) from a given 3D volume that corresponds to an input 2D image.

Several applications requiring slice-to-volume mapping have emerged and pushed the community towards developing more accurate and efficient strategies. Such medical imaging tasks can be classified in two main categories: those related to image fusion for image guided interventions; and those related to motion correction and volume reconstruction. In the first category, pre-operative 3D images and intra-operative 2D images need to be fused to guide surgeons during medical interventions. Slice-to-volume registration plays a key role in this process, allowing the physicians to navigate 3D pre-operative high-resolution annotated data using low-resolution 2D images acquired in real-time during surgery. In the second category, the goal is to correct for misaligned slices when reconstructing a volume of a certain modality. A typical approach to solve this task consists in mapping individual slices onto a reference volume in order to correct for inter-slice misalignment. Again, the development of accurate slice-to-volume registration algorithms is crucial to successfully tackle this problem.

Slice-to-volume registration is also known as 2D/3D registration, due to the dimension of the images involved in the registration process. However, this term is ambiguous since it describes two different problems depending on the technology used to capture the 2D image: it may be a projective (e.g. x-ray) or sliced/tomographic (e.g. US) image. Even if both problems share similarities in terms of image dimensionality, every formulation requires a different strategy to estimate the solution. The lack of perspective and different image geometry [START_REF] Birkfellner | [END_REF]] inherent to both modalities, make it necessary to come up with distinct strategies to solve these registration problems. Moreover, a pixel in any 2D projective image does not correspond to only one voxel from the volume (this is the case for slice-to-volume), but to a projection of a set of them in certain perspective. Therefore, the type of functions used to measure similarities between the images is necessarily different in every case. While most of the projective 2D/3D image registration methods require to bring the images into dimensional correspondence (by different strategies like projection, back-projection or reconstruction [START_REF] Markelj | [END_REF], in case of slice-tovolume registration, pixels from the 2D image can be directly compared with voxels from the volume. In this introduction we focus on the latter case, while a more comprehensive overview about projective 2D to 3D image registration is presented in [START_REF] Markelj | [END_REF] Slice-to-volume registration could be considered as an extreme case of 3D-3D registration, where one of the 3D images contains only one slice. Even if theoreti-cally true, 3D-3D registration methods can not be extrapolated in a straightforward way to the slice-to-volume scenario. This holds, particularly, for registration methods based on image information, since the descriptors used to quantify similarities between images, normally assume that the amount of information available from both images is balanced. The fact that a single slice (or even a few sparse slices) provides less information than an entire volume, should be explicitly considered in the problem formulation. Moreover, specific geometrical constraints like planarity satisfaction and in-plane deformation restrictions, arise in the case of slice-to-volume registration, which are not applicable in the setting of dimensional correspondence.

In this introduction, we discuss the literature of this emerging field, offering a comprehensive survey of the articles based on slice-to-volume registration, proposing a categorical study of the related algorithms according to a taxonomy and analyzing the advantages and disadvantages of each category. We propose a taxonomy based on previous surveys on image registration [START_REF] Markelj | [END_REF], Sotiras 2013], and adapted to the particular slice-to-volume case. We classified the related work according to several principles: (i) matching criterion (section 3.1.2), which specifies the strategies to define similarities between the images; (ii) transformation model (section 3.1.3), indicating the nature of the models used to bring images into spatial correspondence; (iii) optimization (section 3.1.4), differentiating the approaches according to the strategy used to infer the optimal transformation model; (iv) number of slices (section 3.1.5), which splits the methods in two groups, according to whether they require a single or multiple source image slices; and finally the (v) applications (section 3.1.6), where we identify the main clinical scenarios that have motivated research in the field.

Definition of Slice-to-Volume Registration

Let us start by giving a formal definition of slice-to-volume registration. Given a 2D image I and a 3D volume J, we seek a mapping function Θ which optimally aligns the tomographic slice I to the 3D coordinate space of J, through the minimization of the following objective function:

Θ = arg min Θ M(I, J; Θ) + R(Θ), (3.1)
where M represents the image similarity term (the so called matching criterion) and R the regularization term. Note that this mapping may be rigid or non-rigid, depending on whether we allow image I (or its corresponding reformatted slice from J) to be deformed or not. If we estimate only a rigid mapping (i.e. we calculate a 6 degrees of freedom rigid transformation or even a more restrictive one), we name the problem rigid slice-to-volume registration. In case we also infer some sort of deformation model or we consider more expressive linear transformations (such as affine transformations), we call it non-rigid registration. We refer the reader to section 3.1.3 for more information about different transformation models.

Matching criterion M measures the similarity between the 2D image and its corresponding mapping (slice) to the 3D volume. Usually, it is defined using inten-sity information or salient structures from I and J. A complete discussion about matching criteria in the context of slice-to-volume registration (where we also consider methods that do not use any image information but rely on external sensors) is presented in section 3.1.2.

The regularization term R imposes constraints on the solution that can be used to render the problem well posed. It also may encode geometric properties on the extended (plane selection and plane deformation in case of non-rigid registration) transformation model. The choice of regularizer depends on the transformation model. While simple models like rigid body transformations can be explicitly estimated even without regularizer, the term R becomes crucial in more complex non-rigid scenarios to ensure realistic results.

We aim at optimizing the energy defined in equation 3.1, by choosing the best Θ that aligns the 2D and 3D images. Depending on the variables we are trying to infer, and their optimality guarantees, they can be classified in different categories. A full study of this topic is included in section 3.1.4.

The general definition given in equation 3.1 considers a single slice as the input to the registration process. However, for the sake of completeness, in this survey we also consider the so called multi slice-to-volume registration approaches, for which several (but sparse) slices are registered to a full 3D volume. It is worth noting that we do not include methods that perform a prior 3D reconstruction from the input 2D slices before registration (see, for example, the work by [START_REF] Arbel | [END_REF][START_REF] Arbel | [END_REF]), since these methods reduce the problem to the classic 3D-3D scenario, which is not within the scope of our work. On the contrary, we consider those methods that directly register the subset of 2D slices (which can be orthogonal, parallel, arbitrary or even without an a-priori known spatial relation), rather than reconstructing a volume. Discussion and classification of the approaches according to this criterion is presented in section 3.1.5.

Matching Criterion

The matching criterion (also known as (dis)similarity measure, merit function or distance function) quantifies the level of alignment between the images, and it is typically used to guide the optimization process of the transformation model. Depending on the nature of information exploited in the matching process, registration methods can be classified as iconic -we use voxel intensities to quantify similarity-, geometric -we use a sparse set of salient image locations to guide the registrationor hybrid methods -we combine both strategies-. In the particular case of slice-tovolume registration, there are also some approaches which instead of using image information, they rely on other non-image technologies; we will refer to them as sensor based methods.

Slice-to-volume registration excluding medical imaging-based methods is mainly performed using two different technologies: optical (OTS) and electromagnetic (EMTS) tracking systems [Birkfellner 2008]. Optical systems employ different types of cameras (such as infrared (IR), standard RGB or laser cameras) to track markers which help to identify the current position of the objects. Electromagnetic positioning systems perform tracking based on a system of transmitter, sensors and processing unit that localizes the position and orientation of a target object by measuring electromagnetic field properties. Both methods inherit constraints on the nature of transformation that they can estimate, since only rigid transformations can be obtained through these technologies. Moreover, accurate initial calibration is usually required at the beginning, due to the fact that positioning is done through the estimation of the relative displacement from the previous one and therefore errors can be propagated and accumulated through all the process. Instead, registration algorithms exploiting image based information -iconic or geometric-can deal with elastic anatomical changes and they are less sensitive to initial errors (since it is simple to correct them iteratively during the registration process). On the negative side, such methods can be more sensitive to image noise, which is frequently present in intra-operative, real-time and low-quality modalities, normally corresponding to the input 2D image. In addition, the amount of information of an image slice is sparse when compared to a volumetric image, resulting on ambiguities in terms of image matching that render the registration problem ill-posed. Different strategies have been developed to deal with these problems depending on the choice of the matching criterion.

Image registration can be monomodal -when the slice and the volume are captured with the same type of image technology-or multimodal -when slice and volume refer to different modalities, e.g. US slice and MRI or CT volume-. In the former case, the task of measuring the similarity between the images is simpler, since pixel/voxel intensity values corresponding to the same anatomical structure are highly correlated, or even identical, in both images. Therefore, iconic methods may perform better since the main issue associated to them -the difficulty to explain image similarities using pixel/voxel correspondences-is already solved. In case of multimodality, where the relation between pixel intensities is not obvious, there are two major alternatives: to continue using the iconic matching criterion but defining more complex similarity measures, or to adopt a geometric or sensor based strategy which appears to be more robust when dealing with different image modalities.

Iconic

Iconic matching criteria are defined using image intensity information. Similarity between images is measured using functions that act on the pixel/voxel intensity level. Standard signal processing tools, information theoretical approaches or even similarity measures defined for particular image modalities can be considered. The challenge lies in describing both images on a common space where they can be compared, in particular in the context of multimodal setting, where voxel intensities corresponding to the same anatomical or functional structures are dissimilar. There are two desirable properties sought at the definition of any iconic similarity measure: (i) to be convex, since it simplifies the optimization process and (ii) to be discriminative, in the sense that it assigns distinct values to different tissues or Chapter 3. Slice-to-Volume Registration anatomical structures.

Typically, in the context of slice-to-volume registration, an iconic matching criterion δ(I 1 , I 2 ) :

Ω 1 × Ω 2 →
is defined to quantify the similarity between two slices (or patches, according to whether we specify a global or local function). Such a similarity measure varies depending on the modalities we are trying to register. In a monomodal scenario, simple similarity measures such as the sum of absolute differences (SAD) [Ferrante 2013, Ferrante 2015a, Ferrante 2015b[START_REF] Leung | [END_REF], sum of square sifferences (SSD) [Fogtmann 2014,Heldmann 2009[START_REF] Leung | [END_REF][START_REF] Marami | [END_REF], Miao 2014a, Osechinskiy 2011b, Po Su 2013, Seshamani 2013, Xu 2008, Yu 2008, Yu 2011] or even mean of square differences [Fogtmann 2014[START_REF] Gholipour | [END_REF][START_REF] Gholipour | [END_REF],Honal 2010,Kim 2008b,Kim 2010a,Tadayyon 2010a,Tadayyon 2010b] of the intensity values can be used. In vector notation, SAD can be also seen as the L 1 norm of vectorized image, whereas SSD would correspond to the L 2 norm. These metrics assume that there is a straightforward correspondence between the intensity values in both images, which is not necessarily the case.

More complex metrics, exploiting statistical properties of the observed intensity values in both images, have also been proposed. These methods, on top of handling identity transformations, can cope with piece-wise linear relationships between the intensities in the images to be registered. In these cases, image pixels in both images are seen as entries of two random vectors X and Y. Crosscorrelation (CC) is a well known function widely used in the fields of signal processing and statistics, also applied in several slice-to-volume registration studies [START_REF] Birkfellner | [END_REF][START_REF] Elen | Automatic 3-D Breath-Hold Related Motion Correction of Dynamic Multislice MRI[END_REF][START_REF] Fei | [END_REF],Fei 2003a,Fei 2003b,Fei 2004b,Frühwald 2009,Hummel 2008,Jiang 2007a,Jiang 2007b,Jiang 2009,Kainz 2015,Kim 2005,Miao 2014a,Noble 2005[START_REF] Osechinskiy | [END_REF], Osechinskiy 2011b, Xu 2008, Yan 2012[START_REF] Zarow | [END_REF]. CC measures the correlation between the entries of X and Y. It is simple to compute and, more importantly, invariant to shifts and scaling in the intensity domain. Another metric is the correlation ratio (CR), which has shown promising results even in multimodal image registration [Roche 1998]. It measures functional dependencies between X and Y, taking values between 0 (no functional dependence) and 1 (purely deterministic dependence). CR is intrinsically asymmetrical, since the two variables (images) do not play the same role in the functional relationship. In other words, unlike CC, CR offers different values depending on the order that images were considered. CR has been used as an iconic criterion for slice-to-volume registration in [START_REF] Marami | [END_REF], Osechinskiy 2011b, Smolíková-Wachowiak 2005].

Information theoretic similarity measures are usually the choice of preference for multimodal registration, and slice-to-volume registration is not an exception. The most popular is mutual information (MI), which measures the statistical dependence or information redundancy between the image intensities of corresponding distributions in both images, that is assumed to be maximal if the images are geometrically aligned [Maes 1997]. It requires an estimation of joint and marginal probability density functions (PDFs) of the intensities in every image. Given the information sparse nature of slice-to-volume registration when compared to the volume-to-volume scenario, the estimation of these joint PDFs for every slice -especially in slices of low im-age resolution/number of samples-is a hard task and may redound to poor MI-based registration results. An alternative approach to improve MI-based slice-to-volume registration was developed by incorporating informative PDF priors in the context of fMRI time-series registration [Bhagalia 2009]. First, it was shown that slices located near the middle of the head scans give more reliable PDFs and MI estimations because they refer to a richer information space than the end slices (top or bottom). End-slices registration is then guided by a joint PDF prior based on intensity counts from registered center-slices. Alternatively, a better MI calculation in sliceto-volume registration can be achieved by using a PDF estimate that retains as much information about voxel intensities as possible from the higher resolution anatomical data set, when registering 2D MR scout scan to a complete 3D MR brain volume [Chandler 2004]. MI is a widely used similarity measure for slice-to-volume registration, adopted in an important number of methods in the last decades (e.g. [START_REF] Birkfellner | [END_REF][START_REF] Brooks | [END_REF], Eresen 2014[START_REF] Fei | [END_REF], Ferrante 2015b, Fogtmann 2014[START_REF] Gill | [END_REF], Huang 2009, Kim 1999[START_REF] Museyko | [END_REF], Park 2004[START_REF] Rousseau | [END_REF], Seshamani 2013, Smolíková-Wachowiak 2005, Tadayyon 2010a, Xiao 2011, Xu 2014a, Yeo 2004, Zakkaroff 2012[START_REF] Zarow | [END_REF]). One of the main drawbacks of MI, is that it varies when the overlapping area between the images changes, i.e. it is not invariant to changes in the overlap region throughout registration. It could happen that while estimating the transformation model, some potential solutions lie out of the volume. In such cases, an overlap invariant function would be of choice. To this end, a modified version of MI, the normalized mutual information (NMI), can be applied, which is simply the ratio of the sum of the marginal entropies and the joint entropy [Studholme 1999]. Another advantage of NMI with respect to MI is its range: it conveniently takes values between 0 and 1. NMI has been used as well for slice-to-volume registration [START_REF] Chandler | [END_REF][START_REF] Elen | Automatic 3-D Breath-Hold Related Motion Correction of Dynamic Multislice MRI[END_REF][START_REF] Gefen | [END_REF], Hummel 2008, Jiang 2007a, Kainz 2015, Kuklisova-Murgasova 2012[START_REF] Leung | [END_REF][START_REF] Marami | [END_REF], Miao 2014a[START_REF] Rousseau | [END_REF].

Using prior knowledge like segmentation masks during the iconic registration process can be useful. These approaches, also known as region-based methods, employ intensities information or statistics to describe a pre-segmented region. Chan-Vese metric [Chan 2001] for instance, aims to minimize the intensity variances on the regions inside and outside a given segmentation contour. Nir and coworkers [Nir 2011[START_REF] Nir | [END_REF]] applied this matching criterion to the problem of aligning multiple slices of histological images to in vivo MR images of the prostate.

Border information is another low level visual cue that was exploited by iconic methods. It is usually determined from the intensity gradient of the images, which gives an idea of the image structure defined by intensity changes, independently of their actual value. MI as well as SAD or SSD can be applied on top of the gradient magnitudes of both images. In [START_REF] Brooks | [END_REF]], for example, MI between the gradient magnitudes of an US slice and MRI volume is used. [Po Su 2013] applies SSD on both image intensities and gradients, combining them in a unique similarity measure. In [Xu 2014b], the matching criterion is defined using the normalized gradient field of the images, while CC of both intensity and gradient Chapter 3. Slice-to-Volume Registration magnitude was adopted in [Xu 2008].

More robust similarity measures have been introduced more recently to deal with problems in medical image registration, and in the slice-to-volume case in particular. Remarkable contributions have been made by [Wein 2008a, Fuerst 2014] to this field. They propose different similarity measures based on the simulation of US images from MRI and from CT, which can deal with these challenging multimodal registration problems. In [Wein 2007, Wein 2008a], novel methods for simulation of ultrasonic effects from CT data are presented, together with a new similarity measure entitled Linear Correlation of Linear Combination (LC 2 ), which is invariant to missing simulation details, yielding smooth properties and a global optimal corresponding to the correct alignment. Since they simulate US imaging effects with respect to the probe geometry, the original B-mode scan planes of the sweep are used instead of 3D reconstruction, making it suitable for multi slice-to-volume registration. LC 2 was used by [Fuerst 2014], although without the simulation process: the similarity measure is defined by locally matching US intensities to both MRI intensity and gradient magnitude.

Another robust similarity measure, the modality independent neighborhood descriptor (MIND) was proposed by [Heinrich 2012] for multimodal rigid and deformable registration. It is based on the concept of local self-similarity at the patch level in order to create a multi-dimensional descriptor through ranking of the local intensity distribution of the two images, therefore providing a very good representation of the local shape of an image feature. It can be computed in a dense fashion across all the pixels (or voxels) of the images; once it is computed, the SSD of the MIND representations can be used as a similarity measure. Thanks to its point-wise (pixel or voxel-wise) calculation nature, it can be adapted to almost any registration algorithm. An extension to MIND, named Self Similarity Context (SSC) is also estimated using patch-based self-similarities [Heinrich 2013]. In [Cifor 2013b], MIND is successfully used in a multi slice-to-volume registration framework to align untracked freehand 2D US sweeps to CT volumes.

A last family of robust similarity measures was introduced by [Rivaz 2014c, Rivaz 2014b, Rivaz 2014a] in the context of mutlimodal US/MRI image registration. In [Rivaz 2014c], a metric called Contextual Conditioned Mutual Information (Co-CoMI) was proposed. The metric aims at tackling one of the main drawbacks of classic MI based methods, that is taking into account the intensity values of corresponding pixels and not of neighbor. Images are treated as "bag of words" and consequently contextual information is ignored. CoCoMI overcomes this limitation by conditioning the MI estimation on contextual information. In [Rivaz 2014b], Self Similarity α-MI (SeSaMI) -another MI based matching criterion-is proposed. α-MI is usually calculated on multiple features like intensities and their gradients, as opposed to standard MI which is usually calculated on intensities only. SeSaMI combines this multi-feature α-MI formulation with self-similarities in a kNN α-MI registration framework by penalizing clusters (i.e. the nearest neighbors) that are not self-similar. Finally, in [Rivaz 2014a], a CR similarity measure is introduced. Robust Patch Based Correlation Radio (RaPTOR) computes local CR values on small patches and adds them to form a global cost function. Authors claim a property that makes suitable such methods: their metric is invariant to important spatial intensity inhomogenity, what is especially useful when dealing with US images due to wave attenuation, shadowing and enhancement artifacts. One of the main advantages of these Rivaz's metrics (CoCoMI, SeSaMI and RaPTOR) is that their gradient can be derived analytically, and therefore the cost function can be efficiently optimized using stochastic gradient descent methods.

Geometric

Geometric registration finds correspondences between meaningful anatomical locations or salient landmarks [START_REF] Joshi | [END_REF][START_REF] Glocker | [END_REF]]. These methods aim at minimizing an energy function that, for a given transformation, measures the discrepancy between the key-points detected in both, the 2D slice and the volumetric image. Simplicity of the registration process once the landmarks are appropriately determined, no sensitivity to initializations and a wider capture range in terms of deformation are the main strengths of such approaches. On the other hand, the landmark detection and matching processes are not that trivial, and errors on their position compromise the accuracy of the registration process. Moreover, due to the sparsity of the key-points, the quality of the deformations may become insufficient (due to the limited support on the interpolation).

The early work by [Gourdon 1994] presents a geometric method to perform slice-to-volume registration between a curve and a surface. In this work, Gourdon and Ayache exploited the knowledge about the differential properties computed on both, the curve and the surface, to constrain the rigid matching problem. The most relevant contribution of this work is their discussion about how differential constraints can be used to register rigidly a curve on a surface. However, they used the basic Marching Cubes algorithm [Lorensen 1987] to extract them from simulated and medical data. The extraction of these structures remains an arduous task for images with low resolution, a scenario often valid in slice-to-volume registration.

Extracting distinctive features becomes even more complicated for medical images (opposed to natural ones), since the first ones are usually not as discriminative (lack of texture) as the last ones. Invariance to scaling, rotation and changes in illumination or brightness constitute useful properties for methods seeking to extract salient points. Highly distinctive features (in both spatial and frequency domains) simplify the matching task since it is likely that they are correctly matched. Classical examples of such descriptors successfully applied in different computer vision tasks are the Scale-invariant feature transform (SIFT) [Lowe 1999], the Harris detector [START_REF] Harris | [END_REF]], the Histogram of Oriented Gradients (HOG) [START_REF] Dalal | [END_REF]], Speeded Up Robust Features (SURF) [START_REF] Bay | [END_REF]], etc. More recently, features learned using deep learning have been successfully applied for scenarios involving massive amount of annotated data for training [Long 2014]. These features are used to extract the salient landmarks that will guide the geometric registration process. Once the landmarks are established in both slice and volumetric images, there are two choices. First, we can perform a matching process which results in establishing correspondences between pairs of points from the two sets. Once the matching process is finished, correspondences can be used to estimate the desired -rigid or non-rigid-transformation models. Alternatively, methods directly estimating the transformation model without inferring any correspondence could also be applied. Sometimes, both correspondences and transformation model can be inferred at the same time, like in case of Iterative Closest Point (ICP) [Besl 1992] algorithm. ICP is an algorithm subsequently improving the matching of point pairs. It minimizes the sum of geometric distances between the transformed set of source image landmarks and the closest detected landmarks in the target one. It is a simple and fast method that follows the closest neighbor principle to perform the matching task, which often converges to local minimum though.

[ [START_REF] Dalvi | [END_REF]] proposed a slice-to-volume registration approach which uses ICP. First, it extracts phase congruency information from the slices/volume using oriented 2D Gabor wavelets. Then, using local non maximum suppression, a robust and accurate set of feature points is automatically obtained, which are subsequently matched by ICP inferring a rigid body transformation. [START_REF] Nir | [END_REF]] presented a particle filtering framework using a geometric matching criteria where the closest point principle is applied (they use the SSD between pairs of closest points in the two sets). In their proposal, the inferred transformation model involves 9-DOF (affine) matrices instead of 6-DOF transformations (rigid body). Nir's method was applied by [Moradi 2013] as a pre-processing step in a tumor detection learning based framework, to register US to whole-mount histopathology references of prostate images, mapping the location of tumors to the US image domain. Another method that uses ICP was proposed by [START_REF] Yavariabdi | [END_REF]] for deformable registration. They register 2D transvaginal US (TVUS) images with 3D MRI volumes to localize endometrial implants. They use contour to surface correspondences through a novel variational one-step deformable ICP method, finding a smooth deformation field while establishing point correspondences automatically. The main drawback of this approach is that it relies on the user to segment the organs. Manual pelvic organ segmentation is a laborious and time consuming task, and consequently the applicability of the approach in real scenarios is limited. Note that in the three papers discussed in this paragraph, different types of transformation models where estimated based on the closest point principle, showing the flexibility of such approach.

Signed distance functions (or signed distance maps) can be used to avoid the landmark matching step when performing non-rigid registration of shapes and points. In this case, the landmarks or shapes are assigned to zero distance, while the rest of the pixels of the image are labeled with the distance to the nearest geometrical primitive (landmark or boundary). Once the distance map is created, the optimal transformation model can be estimated by means of standard iconic registration (e.g. using SSD on the distance maps) [Paragios 2003, Xiaolei Huang 2006, Taron 2009]. In [START_REF] Boer | [END_REF]], signed distance functions are used to register slices of histological images with a pre-reconstructed 3D model. Authors proposed to use a naive brute force approach for optimization, which incurs in extremely high computational cost.

Sensor Based Methods

Sensor based systems are an alternative to image-based methods (even if images are often used for their initialization). Slice position and orientation are continuously updated from the current one. Such information can be used to extract the corresponding slice from a given volume assuming an initial correspondence. Optical (OTS) and electromagnetic (ETMS) tracking systems are the most widely used technologies to perform this task. OTS systems determine in real-time the position of an object by tracking the positions of either active or passive markers attached to the object. It requires a line-of-sight to be maintained between the tracking device and the instrument to be tracked. This might cause an inconvenience to the physician -especially in image guided interventions-during his work, resulting highly inconvenient. EMTS systems usually consist of three basic components: the electromagnetic field generator, a system control unit that interfaces with a PC, and tracked sensor coils together with their interfaces to the system control unit. By measuring the behavior of each coil, the position and orientation of the object can be determined. EMTS do not inherit the same constraints as the line-of-sight requirements but it is sensitive to distortion from nearby metal sources and exhibits limited accuracy compared to optical tracking [Birkfellner 2008]. Moreover, neither OTS nor EMTS can deal with elastic deformations between volume and slice. In general, one can claim that sensor based methods are usually more reliable (when they are calibrated correctly) than image based methods since they are not affected by ambiguities that arise during the image interpretation process.

In [Gholipour 2011], a three-dimensional magnetic field sensor is used to track the motion of a subject during MRI scanning. It allows estimating the location of the slices for volume reconstruction. The rigid 6-DOF transformation of the sensor in three dimensions (6-DOF) is calculated in real-time using the native gradient fields of the MRI scanner. The relative 3D location of each slice is computed through the sensor motion parameters at the time of slice acquisition. Other works exploiting EMTS systems can be found here [Hummel 2008, Olesch 2011a, Olesch 2011b, Xu 2008].

Optical tracking systems can also replace classical image based slice-to-volume registration algorithms. In [Schulz 2014] a slice-by-slice motion correction for fMRI image reconstruction (see section 3.1.6.2 for a complete discussion about this problem) is achieved thanks to an optical tracking system. Authors proposed to use three tracking cameras with embedded image processing, to track the position of the three optical markers attached to the skull using goggles. Tracking information is then used to replace the classical slice-to-volume registration step [Kim 1999] necessary to account for motion correction during fMRI image reconstruction. In [Bao 2005], optical tracking is used to register laparoscopic US to CT images of a phantom liver, where an an infrared camera was used to localize the US probe. The advantage of infrared cameras is that the position sensor does not pick up interference from reflections and ambient light. However, the line-of-sight requirement still holds and therefore sensor must not be blocked. Other papers using OTS to track the position of image slices when registering them with volumetric images can be found here [Huang 2009, Penney 2006[START_REF] Penney | [END_REF], Schulz 2014].

More recently [Eresen 2014], smartphone tracking technology was considered as a navigation tool to initialize a slice-to-volume registration process between a histological 2D slice and a MR volume. The Inertial Measurement Unit (IMU) of the smartphone (also available as standalone components) is used to define the orientation of the slice. Given the orientation and position data of an IMU sensor, one can then interpolate the corresponding slice from the MR volume. This interactive alignment is applied to determine initial orientation for the 2D slice, and is refined using an iconic MI based registration process optimized via brute force. IMU tracking systems use a combination of accelerometers and gyroscopes to measure acceleration and angular velocity, respectively. Since acceleration is the second derivative of position with respect to time, and angular velocity is the first derivative, angular changes integrated over time from a known starting position yield translation and rotation (i.e. a 6-DOF transformation) [Birkfellner 2008]. These sensors are cheap and widely available. However, a major issue of using IMUs for tracking refers to accumulation of errors (either systematic or statistical), leading to degradation estimation over time. Kalman filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]] can be used to deal with this type of issues and have been often adopted from the scientific community.

Note that, as described in the next section, several papers that propose sensorbased methods, are actually a combination of this technology with some imagebased registration technique, resulting in what we call a hybrid method. [Sotiras 2010b] states that the hybrid matching criteria take advantage of both, iconic and geometric approaches, in an effort to get the best of both worlds. This category can be extended by including sensor based technologies as it can also be combined with iconic or geometric methods. [Mitra 2012b, Mitra 2012a] proposed a hybrid slice-to-volume registration approach, combining geometric and iconic matching criteria in a probabilistic framework to register transrectal US (TRUS) with MR images of the prostate. The geometric component is based on shape-context [Belongie 2002] representations of the segmented prostate contours; Bhattacharyya distance [Bhattachayya 1943] between the shape-context histograms of the two shapes is used to find point to point correspondence in every axial MR image. The Chi-square distances between the TRUS slice and each of the MR slices are calculated and used to determine the matching slice. Once the TRUS-MR slice pair with the minimum Chi-square distance is determined, it is used to retrieve a 2D rigid transformation (in-plane rotation and translation) between them. This transformation is applied to the remaining MR slices to ensure similar 2D in-plane rigid alignment with the 2D TRUS slice.

Hybrid

The iconic step is performed by measuring the similarity (NMI and CC are used as metric) between the TRUS image and every rotated axial slice of the MRI. Finally, shape and image similarity measures are transformed into probability density functions and mapped into a statistical similartity framework towards retrieving the MR slice that resembles the TRUS image. Two issues limit the applicability of such an approach in clinical practice. First, the registration is determined by the manual segmentations of the prostate on both images. Second, the method assumes that the TRUS slices are parallel to the MR axial planes; therefore, the result of the registration does not consider any out-of-plane rotation, which could occur in a realistic scenario.

Tracking information (coming from a EMTS or OTS system) can be combined in a hybrid approach with iconic or geometric criteria to perform slice-to-volume registration. In this case, the tracking signal is used to calculate the additive position of a new 2D slice with respect to the previous one, bringing them to the same coordinate system. Iconic or geometric registration is used to calculate the initial transformation that relates this unique coordinate system with the volumetric image. This approach was used by [START_REF] Penney | [END_REF]] to register intra-operative US images to a pre-operative MR volume. The algorithm was extended in [Penney 2006] to deal with 2D US to volumetric CT registration. [START_REF] San | Towards real time 2D to 3D registration for ultrasound-guided endoscopic and laparoscopic procedures[END_REF] refines the rigid tracking information provided by an EM system attached to a laparoscopic or endoscopic probe. The refinement is made in terms of translation parameters, using an edge-based iconic registration method combined with a phase correlation technique. In [Xu 2008], EM tracking and intra-operative iconic image registration are used to superimpose MRI data on TRUS images of the prostate. To this end, a three-step algorithm is applied, with the intermediate one corresponding to slice-to-volume registration. At the beginning of the TRUS procedure, the operator performs a 2D tracked axial sweep from the prostate's base to its apex, which is used to produce a volumetric US reconstruction of the prostate volume. The reconstructed US volume and the MRI are manually aligned. Thus, 2D TRUS slices are registered to the US reconstructed volume (reducing the MRI/US multimodality problem to US/US monomodality) while the estimated transformation can still be used to recover the corresponding MRI slice. During the intervention, out-of-plane motion compensation is achieved using intermittent multi slice-to-volume registrations between nearly real-time 2D US images and the 3D US image. This rigid body registration is based on minimizing SSD with a Simplex algorithm. The final step refers to a 2D-2D registration of US (using CC and gradient based similarity measures) and it seeks to account for in-plane misalignment. In this workflow, EM tracking information is used to simplify the problem from multi to monomodality registration, and combined with iconic slice-to-volume registration to account for motion compensation. Other hybrid methods that combine tracking information with iconic or geometric registration approaches are [START_REF] Sun | Image guidance of intracardiac ultrasound with fusion of pre-operative images[END_REF], Yan 2012] In [Ghanavati 2010] a method that combines the information from a tracked freehand 2D US transducer with an iconic matching criterion based on US image simulation was proposed, in the context of Total Hip Replacement (THR) surgeries.

They use a statistical shape model (SSM) of the pelvis, constructed from several CT images. They developed a multi-slice to volume registration method, to register multiple 2D US slices to a statistical atlas of the pelvis. The mean shape of the atlas is first registered rigidly to the freehand 2D US images, based on an US simulation method [Wein 2008a]. Then, the atlas is instantiated by measuring the iconic similarity between the actual US and the US images simulated from the instance of the atlas. This two-step algorithm estimating the rigid transformation and followed by the optimization of the deformable parameters, was turned into a one-shot optimized approach in [Ghanavati 2011].

In [Wein 2008b], a position sensor attached to the patient's skin was used to extract a scalar surrogate measurement, which represents an anterior-posterior translation used to detect and compensate for respiratory movements of the liver. This information is combined with an iconic criterion to estimate an affine 12-DOF transformation model which maps pre-operative plans and imaging into the interventional scenario. [Huang 2009] proposed a hybrid slice-to-volume registration strategy exploiting multiple temporal frames. A multimodality image navigation system was introduced that integrates 2D US images with pre-operative cardiac CT volumes, using electrocardigram (ECG) information, optical tracking and iconic MI based registration. The ECG is used to synchronize the US images to the corresponding pre-operative dynamic 3D CT, depending on the cardiac phase indicated by the ECG. Spatial information given by the optical tracker is used to produce a near-optimal starting point for every slice, that is refined by maximizing the MI similarity measure. Such an approach works in real-time, making it applicable to real interventional scenarios. A different combination of sensors to solve the registration problem was proposed in [Hummel 2008], where EM and OTS sensors are combined with iconic registration to register endoscopic US to CT data.

Transformation Model

Transformation models explain the relation between the slice and the volume being registered, and are the outcome of the registration process. They are often classified according to their degrees of freedom. Rigid transformations deal with global rotations and translations, while deformable models -the complex case-can produce local in-plane and out-of-plane deformations. The richness of the deformation model is proportional to the number of parameters we need to specify it, and therefore the trade-off is to be found between the model complexity and power of expression.

Rigid

The simplest transformation model accounts for rotation and translation parameters. It is usually expressed as a 6 degrees of freedom (6-DOF) transformation Θ R composed by 3 rotation and 3 translation parameters. Such a basic model is the most common choice in the literature for slice-to-volume registration. Rigid transformations are expressive enough to explain simple slice-to-volume relations. They can deal with in-plane and out-of-plane translations and rotations. Clinical scenarios that do not inherit image distortion -like simple inter-slice motion correction [Jiang 2007b, Kim 2008b[START_REF] Rousseau | [END_REF], Smolíková-Wachowiak 2005] or basic nature image guided surgeries [START_REF] San | Towards real time 2D to 3D registration for ultrasound-guided endoscopic and laparoscopic procedures[END_REF][START_REF] Birkfellner | [END_REF][START_REF] Gill | [END_REF]]can be modeled with rigid transformations.

When out-of-plane motion is avoided, even simpler models can be used. [Zakkaroff 2012] proposed to recover in-plane slice rotations in cardiac MR series, using the stack alignment transform. In-plane translation along X and Y and rotation around a user-supplied center of rotation for the individual slices were parameterized independently. It also includes a parameter for global translation along the Z direction. Hence, it combines the individual slice transforms into a unique space of parameters -that contains 3-DOF per slice instead of 6-so that all of them can be optimized simultaneously. This type of models has been largely criticized suggesting that (especially in the context of volume reconstruction) the exact corresponding slice can only be found through a slice-to-volume registration which considers outof-plane rotation and translations as well. In [Xiao 2011], authors considered a slice-to-slice matching initial step, where every 2D image from a histological volume is matched to one axial MRI slice, via a group-wise MI based comparison. It implies a restricted transformation model involving simple slice correspondences (1-DOF). Subsequently, they correct the out-of-plane misalignment applying 2D-2D affine registration for each pair of matched slices, and a posterior 3D-3D affine registration between the histology pseudo-volume (reconstructed after transforming every 2D slice) and the MRI volume. Restricted rigid body transformations can be a convenient initialization component of a complete slice-to-volume registration pipeline.

6-DOF rigid body transformations are part of nearly all slice-to-volume registration algorithms. Literature seeking deformable registration [Xu 2014a], often initially employs rigid alignment to account for big range displacements. The standard way to estimate 6-DOF rigid transformations, consists in minimizing an energy functional (based on an iconic or geometric matching criterion) often with a continuous optimization algorithm (see section 3.1.4.1) where the search space is part of the Euclidean group SE(3) of rigid transformations. Sensors as described in section 3.1.2.3 can be used as well for rigid alignment as proposed in a number of papers [Chandler 2008b[START_REF] Elen | Automatic 3-D Breath-Hold Related Motion Correction of Dynamic Multislice MRI[END_REF][START_REF] Fei | [END_REF], Fei 2003a, Frühwald 2009, Fuerst 2014[START_REF] Gholipour | [END_REF],Huang 2009,Jiang 2007a,Kainz 2015,Kim 1999,Kim 2010a,Nir 2011[START_REF] Penney | [END_REF][START_REF] Rousseau | [END_REF],Smolíková-Wachowiak 2005,Xu 2014b,Yeo 2004,Yu 2011].

Non-Rigid

In this review, we call non-rigid transformations to those models which perform -at some extent-changes in the structure of the images. These changes vary from simple operations -which can still be modeled using linear transformations, such as Chapter 3. Slice-to-Volume Registration scale, flip or shear-to more complex models that produce local deformations.

Linear Models

Linear transformations are the first step towards non-rigid registration. In [Wein 2008a], authors suggested that affine models can eliminate most of the large-scale deformations when registering intra-operative 2D US slices to CT images. These deformations are particularly observed between CT and US exams in different respiratory configurations, and could be handled with affine transformations. However, in a different scenario where histological slices of the prostate are registered to a MRI volume, Nir and coworkers [START_REF] Nir | [END_REF]] concluded that the affine transformation model was unable to capture deformations, especially in the rectum area. This is due to the fact that local deformations cannot be expressed from global affine models. In such cases, deformable registration is required (further described in the next section).

Recently, affine models were used for registering 2D slices and 3D images of the bones [Hoerth 2015]. Their proposal includes a novel way to initiate the registration between 2D slices and 3D images of the bones. The method uses the Generalized Hough Transform (GHT) [Ballard 1981] to identify suitable starting positions. They create a template version of the 2D slice by thresholding its gradient vector field, used to explore the 3D volume trough a GHT shape-detection process, providing a set of initial positions. Standard affine registration process is then used to deform images according to the initial configurations of the GHT space and updated using a MI based criterion. Such a principle could be considered to handle alternative transformation models. Other papers that estimate affine models to solve slice-to-volume registration are [START_REF] Gefen | [END_REF][START_REF] Micu | [END_REF][START_REF] Museyko | [END_REF], Wein 2008b, Xiao 2011].

Deformable Models

Elastic deformations are powerful transformations with strong adoption in slice-tovolume mapping. The expressive power of these models depends on the number of parameters used to define them. An extensive description of different types and categories of deformable models used for non-rigid registration is presented in [Sotiras 2013] and [Holden 2008]. Here, we focus on models applied for slice-tovolume mapping.

Thin-Plate Splines (TPS) are frequently used to generate a dense deformation field from a sparse set of control points. These methods involve a set of control points that can be located in arbitrary positions, which are usually obtained by detecting salient structures. Radial basis function (RBF) -where the value at any interpolation point is given as a function of its distance from the control pointsare centered at the control points and combined with affine terms to define an interpolation function. TPS minimize a bending energy based on this interpolation function, which gives a closed-form solution whose uniqueness is guaranteed in most cases. TPS can be decomposed into an affine and a local component. In [Osechinskiy 2011a, Osechinskiy 2011b], TPS are used to parametrize a smooth 3D deformation of a 2D surface (slice). Control points are placed in a regular grid on the 2D image domain, and a 3D warp is defined using three independent TPS functions. Similarly, in [Miao 2014a], several regular 2D grids of control points define one TPS deformation model per slice in a multi slice-to-volume registration scenario. [Kim 2000] also applied TPS to a single slice-to-volume registration problem, and compared two variants of the elastic mapping (involving different number of control points to support the TPS model) with a rigid body registration algorithm. Authors extracted two main conclusions: (i) TPS based registration outperformed rigid body registration, at least in their multi slice fMRI scenario where local deformations are encountered due to the local field induced deformations or localized out-of-plane motion artifacts; and (ii) the number and location of the control points have a significant impact on the final results.

Free-form deformation (FFD) models are also widely applied for medical image registration. Originated in computer graphics [Sederberg 1986], they became popular thanks to the seminal work by [Rueckert 1999]. In this model, the weighting functions are cubic B-splines. The control points, have limited local support (in contrast to TPS where the control points influence the complete domain during interpolation) and are uniformly distributed over the image domain in a grid-like manner [START_REF] Glocker | [END_REF]]. FFDs were presented in section 2.3 of this thesis since they have been used when casting deformable image registration as a discrete problem [START_REF] Glocker | [END_REF]]. In [Ferrante 2013, Ferrante 2015a, Ferrante 2015b], FFDs were used in a graph-based discrete optimization framework to perform slice-to-volume deformable mapping. In this model, a 2D grid-like graph encodes at the same time the plane position (rigid body transformation) and the in-plane deformation of a slice with respect to its corresponding position at the 3D. The in-plane deformation is obtained through a FFD interpolation, where the nodes of the graph play the role of control points. Similarly, but in a continuous setting, [Osechinskiy 2011b] optimizes the position of the control points defined on a 2D grid, with unknown variables corresponding to 3D displacements. In a different scenario, [Fuerst 2014] proposed the use of a 3D grid based FFD within the bounding box of a tracked US sweep. Therefore, they estimated a single 3D deformation field out of a 3D grid, which can be used to deform all the slices contained in the US sweep. The main advantage of the FFD over TPS, is that it does not require solving a linear system for computing weights from displacements. However, FFD imposes constraints in terms of modeling, since the control points must be placed in a regular grid and must enclose the domain boundaries, whereas for TPS, control points can -in principle-be arbitrarily placed in the domain [Osechinskiy 2011b]. Other interesting works combining Bsplines to model a deformable transformation in the context of slice-to-volume registration can be found here [START_REF] Brooks | [END_REF],Honal 2010[START_REF] Museyko | [END_REF],Po Su 2013,Rivaz 2014c, Rivaz 2014b, Rivaz 2014a, Tadayyon 2010a[START_REF] Tadayyon | [END_REF], Xu 2014a].

An alternative model based on finite-element (FE) meshes has been used in [START_REF] Marami | [END_REF]] to model a dynamic linear elastic deformation field. It imposes a regularization constraint on permissible volume transformations based on an iconic image similarity criterion. The advantage of FE models is that they allow incorporating the dynamic behavior of tissue deformation into the registration process, through physically meaningful constraints.

Optimization Method

Optimization methods aim to determine the instance of the transformation model that minimizes a function based on the aforementioned matching criterion (see section 3.1.2). Depending on the nature of the variables being involved, those methods can be classified to continuous or discrete. The continuous approaches exploit the entire space of parameters, while the discrete ones a discretized/quantized version of the admissible solutions.

Numerous problems in computer vision and medical imaging are inherently discrete (like semantic segmentation); however, this is not the case of slice-to-volume image registration, where the search space is continuous. Most of the published methods about slice-to-volume registration adopt a continuous formulation. Nevertheless, recent works on image registration have focused on discrete formulations [START_REF] Glocker | [END_REF][START_REF] Glocker | [END_REF], both for projective and sliced 2D/3D registration [Zikic 2010b, Zikic 2010a, Ferrante 2013, Ferrante 2015a, Ferrante 2015b]. The next subsections present existing work and discuss the limitations of continuous and discrete inference methods in the context of slice-to-volume registration. We also consider a third category of heuristic methods, which are independent of the nature of the variables and can be applied to a wider range of problems at the expense of not providing optimality warranty. These methods are usually applied when finding an optimal solution is impossible or impractical given the nature of the problem objective function.

Continuous

Continuous optimization algorithms are usually iterative methods. They infer the best value for a set of parameters by iteratively updating them. A common mathematical formulation for this strategy is given by:

Θ t+1 = Θ t + ωd t t = 0, 1, 2, 3.... (3.2)
where Θ is the vector of parameters, d t is the search direction at iteration t and ω is the step size or gain factor. The search direction can be calculated using different strategies.

Continuous optimization methods can be classified according to different criteria (e.g. convex or concave minimization, solvers for linear or non-linear functions, smooth or non smooth problems, etc.). Here we simply classify them depending on whether they perform gradient-or non-gradient-based optimization, i.e., whether they exploit first (or higher order) derivatives of the energy function to compute the search direction, or they rely on other strategies. An interesting comparative anal-ysis where methods coming from both classes are used to solve slice-to-volume registration of histological slices with MR images is presented by [Osechinskiy 2011b].

Gradient-based Methods

Gradient-driven methods use the derivative of the objective function to guide the optimization process. In case of convex and differentiable functions, these methods are endowed with optimality guarantees. Otherwise, convergence to local minima is possible, and they could be sensible to initialization. Their main drawback is the requirement of analytical derivation or numerical estimation of the energy function derivatives, reducing their applicability since it is usually complicated.

Gradient descent is the simplest strategy in this category, where the search direction d t is given by the negative gradient of the energy function. It refers to the standard continuous optimization method, and it has been widely applied to the problem of slice-to-volume registration [START_REF] Rousseau | [END_REF][START_REF] Rousseau | [END_REF], Bhagalia 2009, Huang 2009, Tadayyon 2010a[START_REF] Tadayyon | [END_REF], Kim 2010b[START_REF] Marami | [END_REF], Yu 2011, Zakkaroff 2012, Miao 2014a, Xu 2014a].

Conjugate gradient methods use conjugate directions instead of the local gradient to estimate d t . Energy function with the shape of a long and narrow valley, can be optimized using fewer steps than standard gradient descent approach, resulting in faster convergence. [START_REF] Elen | Automatic 3-D Breath-Hold Related Motion Correction of Dynamic Multislice MRI[END_REF], Osechinskiy 2011b] have applied this strategy to estimate rigid and non-rigid slice-to-volume mapping functions, respectively.

Quasi-Newton optimization strategies are used when the Jacobian or Hessian of the energy function cannot be calculated or it is too expensive to be computed. In this case, the search direction at time t is computed based on an estimation of the Hessian, calculated using information provided by the previous iteration t -1. Quasi-Newton optimization was applied to perform deformable registration of brain slices to MR images in [Kim 2005], using several similarity measures. It was also used to register 2D intra-operative US images with pre-operative volumes in [START_REF] Brooks | [END_REF]]. One of the most popular quasi-Newton methods is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm used in [Honal 2010] to correct for breathing motion artifacts during MRI acquisitions. A limited-memory version of BFGS (L-BFGS) is particularly suited for problems involving large numbers of variables and was applied to solve multi-slice to volume registration in [Fogtmann 2014, Xu 2014a].

Gauss-Newton methods can deal with non-linear least squares functions. An approximation of the Hessian matrix is used which, once combined with the gradient, provides a good estimate of the search direction for such functions. Thus, it only requires computing first order derivatives (in contrast to the standard Newton method where the actual Hessian matrix must be computed). In [START_REF] Heldmann | [END_REF], Olesch 2011a, Olesch 2011b], Gauss-Newton method is used for multimodal registration, where 2D US slices of the liver were registered to a pre-operative CT volume. Their matching criterion computes the Sum of Squared Differences (SSD) on the segmentation of the liver vessel structure in both, 3D and 2D slices.

SSD is a least square function and therefore, Gauss-Newton method can be efficiently applied. This method was also used in [START_REF] Lasowski | [END_REF]], solving an Iteratively Reweighter Least Squares (IRLS) process in order to estimate a rigid transformation between fluoroCT and CT images. The Levenberg-Marquardt algorithm is an alternative approach to solve non-linear least squares problems, that requires also only first order derivatives. For well-behaved smooth functions, this algorithm can take more time than standard Gauss-Newton method. However, the method is more robust than the standard Gauss-Newton, meaning that it could find a solution even with bad initializations. To the best of our knowledge, this optimization algorithm was only applied to estimate rigid body transformations; however, it has been used to optimize several iconic matching criteria in a variety of domains such as reconstruction of 3D cell images [Yu 2011], endoscopic interventions [Hummel 2008] and fetal brain imaging [Kim 2008b, Kim 2010a].

Non-Gradient based Methods

Derivative free methods eliminate the differentiability condition of the objective function. They are applicable to a wider range of functions, including noisy, nondifferentiable or even unknown functions, where we only have a black-box that returns an output value given a set of input parameters. These are cases that frequently arise when dealing with medical image registration.

The simplest non-gradient based method that has been used to solve slice-tovolume registration is the local search or best neighbor method. In this case, each parameter to be estimated is perturbed in turn using the stepsize ω, and the value of the objective function is calculated. The search direction d t is then estimated as the one that produced the greatest improvement in the objective function. The main drawback of this approach is its dependency from the choice of the initial stepsize ω. A big stepsize can result in an algorithm moving outside its capture range. On the other hand, for stepsizes that are too small, the optimization may become trapped in a local optimum [START_REF] Penney | [END_REF]]. Such a greedy algorithm has been applied to estimate rigid body and affine transformations mapping slice-tovolume, which have less DOF than deformable models. Different image modalities like US to MRI [START_REF] Penney | [END_REF]] or CT [Penney 2006[START_REF] Sun | Image guidance of intracardiac ultrasound with fusion of pre-operative images[END_REF]] images, as well as fluoroCT images [START_REF] Micu | [END_REF]] were considered in the clinical setting.

Nelder-Mead [Nelder 1965] (also known as downhill simplex method or amoeba method) is the most popular derivative-free method. It relies on the notion of simplex (a n + 1 polytope living in a n-dimensional space) to explore the space of solutions in a systematic way. At every iteration, the method constructs a simplex over the search surface, and the objective function is evaluated on its vertices. The algorithm moves across the surface by replacing, at every iteration, the worst vertex of the current set by a point reflected through the centroid of the remaining n points. The method can converge to a local optimum for objective functions that are smooth and unimodal. However, it exhibits more robust behavior for complex parameter space compared to standard gradient-based methods, providing a good compromise between robustness and convergence time [START_REF] Leung | [END_REF]]. It has been widely used for slice-to-volume registration applications [START_REF] Fei | [END_REF], Fei 2003a, Fei 2003b, Fei 2004b, Frühwald 2009[START_REF] Gill | [END_REF],Hummel 2008,Kim 1999,Kim 2000[START_REF] Leung | [END_REF],Noble 2005, Osechinskiy 2011b,Park 2004,Wein 2008a,Xiao 2011,Xu 2014b[START_REF] Yeo | [END_REF],Yeo 2004].

Another popular non-gradient based method is the well-known Powell's method. It explores the search space by performing bi-directional searches along N different vectors. Usually, these vectors initially refer to the canoninc directions. Then, the search directions are updated using linear combinations of the earlier ones. The algorithm iterates until no significant improvement is made. [Smolíková-Wachowiak 2005] applied Powell's method to register two dimensional cardiac images to preoperative 3D images. Authors claim that similarly to all local techniques, Powell's method converge to local minima, but it is generally robust and accurate, and exhibits fast convergence. [Wein 2008b] takes advantage of the fact that Powell's method performs line search in specific directions. They apply Principal Component Analysis (PCA) on the 12 parameters of an affine transformation, reducing the search space to the 3 most significant PCA modes. Therefore, Powell direction search is initilized with the most significant PCA modes, which assures good performance and robustness. Other papers applying Powell's method to our problem can be found here [Fei 2003a[START_REF] Gholipour | [END_REF][START_REF] Gholipour | [END_REF],Jiang 2007a,Jiang 2007b, Jiang 2009[START_REF] Micu | [END_REF], Osechinskiy 2011b].

Discrete

In a discrete scenario, slice-to-volume registration is expressed as a discrete labeling problem over a graphical model. This is, indeed, the approach considered in this thesis. For the sake of completeness, in this section we are briefly reviewing our own publications, which are (to the best of our knowledge) the only discrete based slice-to-volume registration algorithms proposed to date.

Graphical models and discrete optimization are powerful formalisms that have been successfully used during the past years in the field of computer vision [Wang 2013]. Image registration can be formulated as a minimal cost graph problem where the nodes of the graph correspond to the control points of a deformation grid and the graph connectivity imposes regularization constraints. Even if this technique has been applied mainly to mono-dimensional cases (2D-2D or 3D-3D) [START_REF] Glocker | [END_REF]], our works were focused on extending this theory to the case of slice-to-volume registration [Ferrante 2013, Ferrante 2015a, Ferrante 2015b]. In section 2.3, we showed how monodimensional deformable image registration can be formulated within the discrete MRFs framework, as proposed by [START_REF] Glocker | [END_REF], Glocker 2009c[START_REF] Glocker | [END_REF]]. We refer to [Wang 2013] for a complete discussion about discrete methods in computer vision. [Ferrante 2013] addressed deformable slice-to-volume registration using graphical models and discrete optimization in the so-called overparameterized approach (see section 3.2.1). Grid nodes encode at the same time the plane location and the in-plane deformations. A pairwise model which combines linear and deformable parameters within a coupled formulation on a 5-dimensional label space was used. FastPD was adopted as the optimization algorithm given the properties of the energy function (pairwise terms and non-submodularity). Motivated by the work of [START_REF] Shekhovtsov | [END_REF]], in an effort to reduce the dimensionality of the label space, [Ferrante 2015b] presented a different model (the so called decoupled model, see section 3.2.2 for a complete description) where linear and deformable parameters are now separated into two interconnected subgraphs which refer to lower dimensional label spaces. It reduces the dimensionality of the label space by increasing the number of edges and vertices, while keeping a pairwise graph that was optimized using loopy belief propagation. Finally, in [Ferrante 2015a], a high-order approach is presented (see section 3.2.3), where the label space dimensionality reduction is achieved by augmenting the order of the graphical model, using third-order cliques which exploit the expression power of this type of variable interactions. Such a model provides better satisfaction of the global deformation constraints at the expense of complex inference. Loopy belief propagation was also used as the optimization algorithm in this case.

Miscellaneous

Independently of whether the variables are continuous or discrete, different strategies can be used to explore the space of solutions. Heuristic or metaheuristic can be used to obtain acceptable solutions especially in cases where we deal with nonlinear, non-convex or black-box optimization problems, even if we know that no optimality guarantee is provided. In some cases (e.g. initializing a more complex registration procedure) having an approximately good solution is enough; in these scenarios, different strategies can be envisioned.

Evolutionary algorithms are among the most popular bio-inspired metaheuristics. In evolutionary computation, the concept of biological evolution is used to explore the search space. The individuals of a given population are the candidate solutions to our problem and they evolve according to different laws (such as mutation, recombination and selection). Evaluation of the quality for a given solution is performed through a fitness function, which is in fact our objective function. The idea is that, over the generation sequence, individuals with better and better fitness are generated, leading to a good solution. These algorithms are usually stochastic. If this is the case, new candidate solutions might be sampled according to a multivariate normal distribution in the space of parameters. One evolutionary algorithm applied in slice-to-volume registration problems is the covariance matrix adaptation evolutionary strategy (CMA-ES) [Ghanavati 2010,Tadayyon 2010a[START_REF] Tadayyon | [END_REF]]. This method updates the covariance matrix of the aforementioned distribution, by estimating a second order model of the underlying objective function similar to the approximation of the inverse Hessian matrix. Another case based on genetic algorithms to solve our problem of interest was proposed by [Gefen 2005[START_REF] Gefen | [END_REF]]. Different optimization strategies may be combined as well depending on the parameter space. In [Tadayyon 2010a[START_REF] Tadayyon | [END_REF], for example, the CMA-ES was applied to optimize the translation parameters of the rigid transformation model, while rotation and deformable ones were estimated using a standard gradient descent approach. In their application to prostate MR images, the CMA-ES was not able to optimize a 6-DOF search space as it diverged on rotations regardless of scaling.

Another popular metaheuristic is the well known simulated annealing (SA) method. It mimics the physical process that metal atoms suffer when the material is heated and then slowly cooled. In order to avoid local minima, the algorithm explores new directions which lower the objective function, but also, with a certain probability, those that raise the objective. This probability decreases with the number of iterations, so that it stays away from local minima in early iterations and is able to explore globally for better solutions. [START_REF] Birkfellner | [END_REF]] used this method to register fluroCT slice to CT volumes optimizing several standard iconic matching criteria, while [Cifor 2013b] applied it to a multi-slice to volume registration case where a robust modality-independent similarity measure was optimized.

Number of Slices

In this review we include single and multi slice-to-volume registration methods. While in the first case the estimated slice-to-volume mapping function considers just a single 2D slice, in the second case several slices are mapped to the 3D volume.

A trivial extension of any single slice-to-volume registration method to the multislice case could be achieved by simply applying the algorithm to every input slice independently. It would allow parallelizing the process -since every registration could be performed in parallel-but, at the same time, contextual information would be lost since no relation among the input slices would be considered. Several methods of this type have been proposed. We can cite [Fei 2003a, Ferrante 2015a, Osechinskiy 2011b[START_REF] Yavariabdi | [END_REF] among others.

Different strategies can be adopted to register multiple slices to a given volume. Here we only review methods that perform a consistent transformation to the set of input slices, but still conceive the process as slice-based without seeing them as a complete and unique volume. Methods that reconstruct a new volume using the input slices, and then perform a standard 3D-3D registration process, are not taken into account. Note that, in some works like [START_REF] Brooks | [END_REF], Yan 2012, Xu 2014a], experiments have been performed to compare the performance of multi slice-tovolume versus volume-to-volume registration. In [Xu 2014a], the comparison was done in the context of MR prostate images, in order to show that multi slice-tovolume registration is sufficient in capturing prostate motion intra-operatively. The conclusion was that multi slice-to-volume registration was able to produce results that were close enough to volume-to-volume registration. This is an encouraging result, considering the data computation reduction archived when using only a few slices instead of a complete volume. [START_REF] Brooks | [END_REF]] and [Yan 2012] arrived at similar conclusions in different medical scenarios. This consensus indicates the reliability of using sparse slices instead of the full volume when possible. Several advantages are associated to this strategy: (i) the computational requirements are lower thanks to the reduction in the amount of data to be processed (only a few slices vs a complete volume), (ii) the omission of the potentially complex reconstruction step, (iii) greater adaptability of the similarity measures and (iv) easier parallelization [START_REF] Brooks | [END_REF]].

In the presence of multiple slices, multi slice-to-volume procedure -as opposed to single slice procedure-greatly increases the quality of the matching leading to more robust registration methods [Yu 2011]. This is even more interesting when the relative position between the slices is known. In [START_REF] Chandler | [END_REF], Chandler 2008b], for example, multi slice-to-volume registration was used to correct for misaligned cardiac anatomy in Short Axis (SA) images by registering stacks of two parallel slices (which are supposed to be aligned between them) to a high-resolution 3D MR axial cardiac volume. In other cases, different configurations like orthogonal slices [START_REF] Gill | [END_REF][START_REF] Leung | [END_REF][START_REF] Tadayyon | [END_REF],Miao 2014b,Xu 2014a], slices aligned in arbitrary positions (obtained using a tracking system) [Wein 2008a,Heldmann 2009[START_REF] Heldmann | [END_REF],Olesch 2011a,Olesch 2011b,Cifor 2013b] or even temporal sequences of slices [Miao 2014a], are registered to a volume. Another case is [Yu 2008,Yu 2011], where a bi-protocol was proposed to reconstruct a microscopic volume of a cell. The bi-protocol is composed by two sets of multiple slices imaging the same cell, which are captured with two different geometries. Both sets of slices are then registered to each other and a final volume is reconstructed using a simple interpolation strategy.

Let us recall that multi slice-to-volume registration involves lower computational complexity compared to 3D-3D registration methods. Furthermore, when compared to single slice methods, multi slice-to-volume registration certainly improves the robustness of the registration process by augmenting the image support. Dealing with several slices requires more computing power than the single slice case; however, on occasions like freehand US sweeps, slices contain redundant information that could be avoided. In that sense, as the complexity of the scheme is proportional to the number of input slices, [Wein 2008a] proposed a strategy which starts out by selecting only a few key-slices. These are used to estimate a rigid or affine transformation model, mapping the US slices to a CT volume. Since neighboring frames of the freehand US sweep contain overlapping information, only a few key slices are selected. They chose the most informative slices by selecting those with the highest image entropy. [Olesch 2011a] applied the same key-slice selection technique to the variational deformable registration framework proposed by [Heldmann 2009]. To this end, well distributed slices throughout the volume were considered, containing meaningful information in terms of entropy. In another work, [Olesch 2011b] extended this key-slice selection technique by introducing a focused registration strategy that only considers slices which are in a given region of interest (ROI).

In [Jiang 2007b], a multi slice-to-volume registration method -entitled Snapshot magnetic resonance imaging (MRI) with Volume Reconstruction (SVR)-is proposed to deal with reconstruction of MRI of moving subjectss. After imaging the target object of interest repeatedly (producing multiple overlapped stacks of slices of the same moving object) to guarantee sufficient sampling, one of the stacks is chosen.

Then, an iterative mutli slice-to-volume registration strategy is applied where the stacks are subsequently registered (in a rigid way) to the reference. Once all the stacks are registered, they are averaged to produce a new reference. Then, the number of temporally contiguous slices in every stack is reduced so that the data is divided into subpackages that are temporally contiguous although the slices in each subpackage may not be spatially adjacent. The process is then repeated: the subpackages are registered to the average image, which is again updated, and the number of slice per stack is once more decreased; this process is repeated until each slice is treated in isolation (reducing to the single slice-to-volume case) and the final average volume is reconstructed. In [Jiang 2007a, Jiang 2009] the idea was extended to allow the reconstruction of diffusion tensor images in the same scenario. In the same work, [Jiang 2007b] tested the benefit of using both, parallel and orthogonally acquired slices with prospectively acquired data and simulated cases as well. More recently, a fast multi-GPU accelerated framework to perform SVR was presented [Kainz 2015]. They proposed a fully parallel SVR approach for the reconstruction of high-resolution volumetric data from motion corrupted stacks of images. Parallelization is performed at two levels: (i) at the slice level, multiple slices are treated separately for large parts of the reconstruction process and (ii) since pixel/voxel based operations are independent of each other, they can be also executed in parallel. Authors claimed that their approach is five to ten times faster than the fastest currently available multi-CPU frameworks. [START_REF] Rousseau | [END_REF][START_REF] Rousseau | [END_REF]] proposed the registration of multiple sets of orthogonal 2D MRI slices into a high resolution MRI volume. The work of [Jiang 2007b] shares similarities some with them where, instead of reducing the number of slices per subpackage until each slice is treated in isolation, a two-step approach is performed. In this setting, a global alignment of the low resolution images was performed first and then each slice of the low resolution images was registered to the reconstructed high resolution volume. Only orthogonal sets of slices were considered for reconstruction.

The work by [Kim 2008b,Kim 2010a], built on top of [START_REF] Rousseau | [END_REF][START_REF] Rousseau | [END_REF]], proposed a new approach entitled Slice Intersection Motion Correction (SIMC). It considers the registration process directly in terms of the intersections of each pair of slices in the stacks, avoiding the intermediate volume estimation process. Independent per-slice rigid transformations are estimated by minimizing a global energy function defined by the sum of dissimilarity measures of all intersection profiles between any two orthogonal slices. SIMC can be used in this scenario, considering intersecting lines of voxels instead of standard patch or global based calculation. [Kim 2008b, Kim 2010a] applied SIMC to the problem of reconstructing MR images of a moving human fetal brain. Since then, several extensions to this work have been published, like for example the work presented in [Kim 2010b] expanding the method to the case of Diffusion Tensor Images and [Kio Kim 2010, Kim 2011] which modified it in order to account for bias field inconsistency correction on fetal brain MR images.

An approach based on particle filtering (PF) was considered in [Nir 2011[START_REF] Nir | [END_REF][START_REF] Nir | [END_REF]] to deal with multi slice-to-volume registration of histological images to MR or US volumes. PF framework derives an optimal estimation of the parameters in a Bayesian fashion, tackling two of the main issues that arise from multimodal registration, such as susceptibility to initialization and optimal solutions.

Applications

A broad number of clinical settings and applications can benefit from slice-to-volume registration. In this review we classify them in two main categories.

Image Fusion and Image Guided Interventions (IGI)

Several medical procedures such as image guided surgeries and therapies [START_REF] Fei | [END_REF]], biopsies [Xu 2014a], radio frequency ablation [Xu 2013], tracking of particular organs [START_REF] Gill | [END_REF]] and minimally-invasive procedures [START_REF] Liao | [END_REF], Huang 2009] belong to this category. In this context, slice-to-volume registration brings high resolution annotated data into the operating room. Generally, pre-operative 3D images such as computed tomography (CT) or magnetic resonance images (MRI) are acquired for diagnosis and manually annotated by expert physicians. During the surgical procedure, 2D real time images are generated using different technologies (e.g. flu-oroCT, US or interventional MRI slices). The alignment of intra-operative images with pre-operative volumes augments the information that physicians have access to, and allows them to navigate the volumetric annotation while performing the operation. These intra-operative images inherit lower resolution and quality than the pre-operative ones. Moreover, tissue shift collapse as well as breathing and heart motion during the procedure, cause elastic deformation in the images, what makes slice-to-volume registration an extremely challenging task. A statistically significant improvement in alignment has been demonstrated when comparing automatic methods to manual (human) results, showing the importance of automatic slice-to-volume registration algorithms in the context of image fusion and IGI [Frühwald 2009].

The pioneering work of [START_REF] Fei | [END_REF][START_REF] Fei | [END_REF], Fei 2003a, Fei 2004a, Fei 2004b] introduced iconic slice-to-volume registration to the problem of image fusion in the context of image guided surgeries. The motivation was that low-resolution Single Photon Emission Computed Tomography (SPECT) can be brought to operating by pre-registering it with a high-resolution MRI volume, which could be subsequently fused with live-time iMRI. That is how, by registering the high-resolution MR image with live-time iMRI acquisitions, Fei and coworkers could map the functional data and high-resolution anatomic information to live-time iMRI images for improved tumor targeting during thermal ablation. [START_REF] Birkfellner | [END_REF]] used slice-to-volume registration to fuse 2D fluoroCT with volumetric CT, which is a well know tool for image-guided biopsies in interventional radiology. In this case, the pre-interventional diagnostic high resolution CT with contrast agent is used to localize a lesion in the liver. However, during the inter-vention, the lesion is no longer visible. Thus, localizing the slice of the CT that corresponds to the intra-operative fluoroCT allows doctors to find the lesions during the biopsy. This approach only considers rigid transformations. However, interventional procedures like radio frequency ablation (RFA) or image-guided biopsies, which use fluoroCT as image guiding technology, are performed while the patient is breathing continuously. Therefore, deformations should also be taken into account when registering with the pre-operative static CT image. The influence of such deformations and the reliability of performing non-rigid registration in such scenario was discussed in [START_REF] Micu | [END_REF][START_REF] Lasowski | [END_REF]. It was claimed that a 2D-3D non-rigid registration solution -based on the single low quality fluoroCT-cannot be precise as required to performed medical procedures. This is mainly due to the poor support in terms of liver anatomical features (mainly vessels) provided by the fluoroCT slices. They proposed to overcome this limitation by providing an adaptive visualization [START_REF] Lasowski | [END_REF]] of the volume area surrounding the minimum estimated pose. This approach copes with the uncertainty in estimating the deformation and brings more information than a single registered slice. Their method performs rigid slice-to-volume registration, and includes views of the CT-Volume determined along flat directions of the out-of-plane motion parameters next to the minimum pose. Laparoscopic and endoscopic interventional procedures also exploit slice-to-volume registration. [START_REF] San | Towards real time 2D to 3D registration for ultrasound-guided endoscopic and laparoscopic procedures[END_REF], proposed a method to register endoscopic and laparoscopic US images with pre-operative computed tomography volumes in real time. It is based on a new phase correlation technique called LEPART accounting for rigid registration. Other methods applying slice-to-volume registration on minimally invasive procedures can be found here [START_REF] Heldmann | [END_REF], Bao 2005].

Motion Correction / Volume Reconstruction

The second category is motion correction and volume reconstruction. Here, the goal is to correct misaligned slices when reconstructing a volume of a certain modality. A typical approach to solve this task consists in mapping individual slices within a volume onto another reference volume in order to correct the inter-slice misalignment. The popular map-slice-to-volume (MSV) method introduced this idea [ Kim 1999]. More recently, applications of slice-to-volume registration to the same problem in different contexts like cardiac magnetic resonance (CMR) [Chandler 2008a[START_REF] Elen | Automatic 3-D Breath-Hold Related Motion Correction of Dynamic Multislice MRI[END_REF], fetal images [Seshamani 2013] and diffusion tensor imaging (DTI) [Jiang 2009] have shown promising results. In these problems, it is usually assumed that a single slice is coherent, in the sense that spatial inconsistency only happens at the inter-slice level.

Slice motion correction in the context of volume reconstruction, typically involves iterative registration of slices to a target volume. The target volume may be an anatomical reference or it could be reconstructed at each iteration using current estimates of slice motion, considering all possible views of the subject. [Kim 1999] introduced the map-slice-to-volume (MSV) approach for inter-slice motion correction in time series of fMRI image. In such cases, head motion rep-resents the major source of error in measuring intensity changes related to given stimuli in fMRI time series [Yeo 2004]. The aim of MSV is to retrospectively remap slices that are shifted by head motion to their spatially correct location using an anatomical MRI volume as reference. The MSV approach estimates 6-DOF rigid transformations independently for every fMRI slice, by minimizing the energy based on mutual information, using a Nelder-Mead downhill simplex method for optimization. This method was presented as an alternative to the slice-stack approach. Instead of considering slice-wise registration, it assumes stacks of slices being already registered among them, and ignores the inter-slice motion inherent to multi slice echo planar imaging (EPI) acquisition sequence (since each slice is excited at a sequential time interval). The MSV showed better performance than previous volume-to-volume registration methods, but the reliability of the estimated position parameters for the end cap slices was low due to the limited information support of smaller regions (less textured area and more background). [Park 2004] proposed to overcome this limitation using Joint Mapping of Slices into Volume (JMSV). JMSV is a multi-slice registration method that jointly estimates a rigid body transformations per slice, while penalizing the implied acceleration in the motion trajectory of the subject -i.e. abrupt changes in the motion parameters of sequentially acquired slices. Other extensions to the standard MSV include (i) accounting for deformable registration [Kim 2000] through TPS transformation models; (ii) improving motion correction capability of MSV with concurrent iterative field-corrected reconstruction [Yeo 2004, Yeo 2008]; (iii) extensive evaluation of the activation detection performance of MSV and effects of temporal filtering of motion parameter estimates [START_REF] Yeo | [END_REF]] and spin saturation effect [Kim 2008a]; (iv) improving the MI matching criterion estimation in the low-complexity end-slices (slices near top or bottom of the head scans, presenting poorer information) by incorporating joint probability density functions of image intensities estimated from successfully registered center-slices in the same time-series.

After the fundamental work of [Kim 1999] and its extensions (most of them related to fMRI image reconstruction), another problem requiring to correct for inter-slice motion started to attract attention from the community: fetal brain MR imaging. In this case, there is no anatomical reference available to be used as a target volume; therefore, a reference volume is calculated during the registration process using the current estimate of the slices. [Jiang 2007a] and [START_REF] Rousseau | [END_REF]] established the basis of a new family of methods which rely on iterations of sliceto-volume registration and scattered data interpolation (SDI) to perform superresolution reconstruction of moving objects -in particular, fetal brain MRI-. Both approaches share a similar and iterative slice-to-volume registration scheme, but differ on the SDI method: while [START_REF] Rousseau | [END_REF]] use Gaussian kernel-based SDI, [Jiang 2007a] relied on a regular grid of control-points that acts as a cubic Bspline to perform SDI. In a more recent work, [START_REF] Gholipour | [END_REF][START_REF] Gholipour | [END_REF] criticized SDI-based approaches stating that none of these techniques can guarantee the convergence of the reconstruction to, at least, a local optimal solution. Thus, they proposed the use of a maximum likelihood (ML) error minimization approach, which guarantees the convergence of the reconstructed image to match the motioncorrected slices through steepest descent error minimization. An extension to this method was presented in [Gholipour 2011], using a real-time sensor based tracking (i.e. a non-image based approach) to estimate an initial position of every slice, which is then refined through a retrospective slice-to-volume registration approach based on ML error minimization.

Graph-based Slice-to-Volume Deformable Registration

In the previous section, we presented a literature review about slice-to-volume registration of biomedical images. In the current section, we introduce a novel sliceto-volume registration framework, which models this challenging task as a discrete labeling problem. As mentioned, our work is pioneer in using discrete methods to solve this problem. We already discussed how image registration can (in general) be formulated using MRFs and discrete inference in section 2.3. Here we extend this approach to the particular case of slice-to-volume registration, and we propose three different derivations that can solve this problem. In section 3.1.1 we give a general definition of slice-to-volume registration, which is reflected in equation 3.1. This definition contains a generic transformation model Θ that maps the 2D input image I to the corresponding slice from volume J. In our non-rigid discrete setting, that function includes a 2D-2D in-plane local deformation field TD and a plane π[J] (i.e. a rigid transformation that specifies a bi-dimensional slice from the volume J). The general optimization problem defined in equation 3.1 is therefore reformulated as:

TD , π = arg min T D ,π M(I • T D (x), π[J](x)) + R(T D , π), (3.3)
where M represents the data similarity term and R the regularization term. The data term M measures the matching quality between the deformed 2D source image and the corresponding 3D slice. The regularization term R imposes certain constraints on the solution that can be used to render the problem well posed. It also imposes certain expected geometric properties on the extended (plane selection and plane deformation) deformation field. The plane π, that minimizes the equation, indicates the location of the 3D volume slice that best matches the deformed source image. The deformation field TD represents the in-plane deformations that must be applied to the source image in order to minimize the energy function.

The fundamental idea behind our approaches is quite intuitive: we aim at deforming a planar 2D grid in the 3D space, which encodes both the deformation field TD and the plane π at the same time. This grid is super-imposed to the 2D source image and consists of control points that jointly represent the in-plane deformation and the current position of the 2D image into the 3D volume. The source image is positioned within the volume by applying different displacement vectors with respect to the control points of the superimposed grid. These displacements are chosen such that a given energy (see Eq. 3.3) is minimized to best fit the matching criterion M. Since they can be moved without any restriction, geometric constraints are imposed through the regularization term R in order to keep a smooth deformation field and a planar grid. Given that we impose a soft planar constraint, the resulting grid is approximately planar. Therefore, we reconstruct the final solution by projecting all the points into a regression plane which is estimated out of the current position of the points. The rigid transformation that indicates the position of the regression plane is considered as π. Finally, the projected grid is interpreted as a 2D Free Form Deformation model (FFD) [Rueckert 1999] where each control point has local influence on the deformation and is used to approximate the dense deformation field TD (see Figure 3.2). Having said that, alternative control point interpolation models can be used as well. Figure 3.1 illustrates this workflow.

This general formulation can be expressed through different discrete labeling problems on a graph by changing its topology, the label space definition and the energy terms. As we mentioned, in this work we propose three different approaches to derive slice-to-volume registration as a discrete graph labeling problem. First, we propose the so-called overparameterized method, which combines linear and deformable parameters within a coupled formulation on a 5-dimensional label space [Ferrante 2013]. The main advantage of such a model is the simplicity provided by its pairwise structure, while the main disadvantage is the dimensionality of the label space which makes inference computationally inefficient and approximate (limited sampling of search space). Motivated by the work of [START_REF] Shekhovtsov | [END_REF]], we present a decoupled model where linear and deformable parameters are separated into two interconnected subgraphs which refer to lower dimensional label Figure 3.2: Recovering the rigid and deformable parameters. Since the resulting grid after inference is approximately planar, we recover the final solution by projecting all the points into a regression plane which is estimated out of the current position of the points. The rigid transformation π indicates the position of the regression plane. The projected grid is interpreted as a 2D FFD, used to approximate the dense deformation field TD .

spaces [Ferrante 2015b]. It allows us to reduce the dimensionality of the label space by increasing the number of edges and vertices, while keeping a pairwise graph. Finally, in the high-order approach [Ferrante 2015a], we achieve this dimensionality reduction by augmenting the order of the graphical model, using third-order cliques which exploits the expression power of this type of variable interactions. Such a model provides better satisfaction of the global deformation constraints at the expense of quite challenging inference.

Overparameterized Approach

Let us consider an undirected pairwise graph G O = V, E super-imposed to the 2D image domain with a set of nodes V and a set of cliques E. The nodes V (a regular lattice) are interpreted as control points of the bi-dimensional quasi-planar grid that we defined in the previous section. The set of edges E is formed by regular 4-neighbors grid connections and some extra edges introduced to improve the propagation of the geometrical constraints (see figure 3.3.a). The vertices v i ∈ V c) Unit vectors in spherical coordinates corresponding to the last two coordinates of a label from the overparameterized approach N l i = (φ, θ). (d) Displacement of the control points p i and p j when the corresponding labels l i associated to (d l i , N l i ) and l j associated to (d l j , N l j ) are applied. The planes π i and π j are those that contain the control points p i + d l i , p j + d l j and whose normals are N l i , N l j respectively. are moved by assigning them different labels l i ∈ L (where L corresponds to the label space) until an optimal position is found. Γ O represents a complete labeling that assigns one label to each variable (control point).

In order to deform the graph, we define a label space describing the inplane deformations and the plane selection variables. To this end, we consider a label space L where every label l i is associated to 5-tuples (d x , d y , d z , φ, θ), where the first three parameters (d x , d y , d z ) define a displacement vector d l i in the cartesian coordinate system (see figure 3.3.b), and the angles (φ, θ) define a vector N l i on a unit sphere, expressed using spherical coordinates (see figure 3.3.c). Let us say we have a control point p i = (p xi , p yi , p zi ) and we assign to this point the label l i , associated to the 5-tuple (d xi , d yi , d zi , φ i , θ i ). So, the new point position p i after assigning the label is calculated using the displacement vector as given by the following equation:

p i = (p xi + d xi , p yi + d yi , p zi + d zi ).
(3.4)

Additionally, we define a plane π i containing the displaced control point p i and whose unit normal vector (expressed in spherical coordinates and with constant radius r = 1) is N l i = (φ i , θ i ). One of the most important constraints to be considered is that our transformed graph should have a quasi-planar structure, i.e. it should be similar to a plane; the plane π i associated with every control point p i is used by the energy term to take into account this constraint. Figure 3.3.d shows how to interpret the labels for two given points p i and p j .

The energy to be optimized is formed by data terms G = {g i (•)} (or unary potentials) associated with each graph vertex and regularization terms F = {f ij (•, •)} (or pairwise potentials) associated with the edges. As it was mentioned, the first ones are typically used for encoding some sort of data likelihood, whereas the later ones act as regularizers and thus play an important role in obtaining high-quality results [START_REF] Glocker | [END_REF]]. The energy associated to the overparameterized formulation is defined as:

P O (Γ O ; G, F ) = i∈V g i (l i ) + γ (i,j)∈E f ij (l i , l j ), (3.5)
where l i , l j ∈ L are the labels assigned to the vertices v i , v j ∈ V respectively. Then, we can define slice-to-volume registration as a MAP inference problem using the overparameterized formulation as: ΓO = arg min

Γ O P O (Γ O ; G, F ), (3.6)
where ΓO is the optimal labeling used to recover the rigid transformation and dense deformation field as explained in the previous section.

The formulation of the unary potentials that we propose is independent of the similarity measure. It is calculated for each control point given any intensity based metric δ capable of measuring the similarity between two 2D images (e.g sum of absolute differences, mutual information, normalized cross correlation). This calculation is done for each control point p i , using its associated plane π i in the target image J and the source 2D image I. An oriented patch Ω l i over the plane π i (centered at p i ) is extracted from the volume J, so that the metric δ can be calculated between that patch and the corresponding area from the source 2D image (see Figure 3.4):

g i (l i ) = Ω l i δ(I(x), π i [J](x))dx.
(3.7)

One of the simplest and commonly used similarity measures is the Sum of Absolute Differences (SAD) of the pixel intensity values. It is useful in the monomodal scenario, where two images of the same modality are compared and, therefore, the grey intensity level itself is discriminant enough to determine how related are the two images. Its formulation in our framework is:

g SAD,i (l i ) = Ω l i | I(x) -π i [J](x) | dx.
(3.8)

In multimodal scenarios, where different modalities are compared (e.g. CT with US images), statistical similarity measures such as Mutual Information (MI) are generally used since we can not assume that corresponding objects have the same intensities in the two images. MI is defined using the joint intensity distribution p(i, j) and the marginal intensity distribution p(i) and p(j) of the images as:

g MI ,i (l i ) = - Ω l i log p(I(x), π i [J](x)) p(I(x))p(π i [J](x)))
dx.

(3.9) As we can see in the previous examples, our framework can encode any local similarity measure defined over two two-dimensional images.

Let us now proceed with the definition of the regularization term. Generally, these terms are used to impose smoothness on the displacement field. In our formulation, the pairwise potentials are defined using a linear combination of two terms: the first (F 1 ) controls the grid deformation assuming that it is a plane, whereas the second (F 2 ) maintains the plane structure of the mesh. They are weighted by a coefficient α as indicates the following equation:

f ij (l i , l j ) = αF 1i,j (l i , l j ) + (1 -α)F 2i,j (l i , l j ).
(3.10)

The in-plane deformation is controlled using a distance preserving approach: it tries to preserve the original distance between the control points of the grid. Since this metric is based on the euclidean distance between the points, it assumes that they are coplanar. We use a distance that is symmetric, based on the ratio between the current position of the control points p i , p j and their original position p o,i , p o,j :

ψ i,j (d l i , d l j ) = || (p i + d l i ) -(p j + d l j ) || || (p o,i ) -(p o,j ) || . (3.11)
Once we have defined ψ ij , the regularizer should fulfill two conditions: (i) it has to be symmetric with respect to the displacement of the points, i.e. it must penalize equally whenever the control points are closer or more distant; (ii) the energy has to be zero when the points are preserving distances and monotonically increasing with respect to the violation of the constraint. The following regularization term fulfills both conditions for a couple of nodes i, j ∈ V labeled with labels l i , l j :

F 1i,j (l i , l j ) = (1 -ψ i,j (d l i , d l j )) 2 + (1 -ψ i,j (d l i , d l j ) -1 ) 2 ,
(3.12) The plane preservation term is based on the average distance between a given control point and the plane defined from the neighboring ones (see Figure 3.5.b). The aim is to maintain the quasi-planar structure of the grid. Given that the distance between a point and a plane is zero when the point lies on the plane, this term will be minimum when both of the control points are on the same plane.

The distance between a point p = (p x , p y , p z ) and a plane π defined by the normal vector N = (n x , n y , n z ) and the point q = (q x , q y , q z ) is calculated as:

D π (p) = | n x (p x -q x ) + n y (p y -q y ) + n z (p z -q z ) | n 2 x + n 2 y + n 2 z . (3.13)
F 2 is defined using this distance (equation 3.13) and corresponds to the average of D π j (p i + d l i ) and D π i (p j + d l j ):

F 2i,j (l i , l j ) = 1 2 (D π j (p i + d l i ) + D π i (p j + d l j )). (3.14)
Recall that normal vectors in our label space are expressed using spherical coordinates with a fixed radius r = 1 (unit sphere). However, the formulation that we presented uses cartesian coordinates. Therefore, the mapping from one space to another is done as follows:

x = r sin(θ) cos(φ), y = s sin(θ) sin(φ), z = r cos(θ).

(3.15)

Note that such pairwise terms are non submodular since we include the current position of the points (which can be arbitrary) in their formulation and therefore the submodularity constraint is not fulfilled. In this context, even if there is no energy bounding that guarantees certain quality for the solution of the optimization problem, good empirical solutions are feasible since we are in a pairwise scenario. Still, two issues do arise: (i) high dimensionality of the label space and consequently high computational cost, (ii) insufficient sampling of the search space and therefore suboptimal solutions. In order to address these issues while maintaining the pairwise nature of the methods, we propose the decoupled method inspired by [START_REF] Shekhovtsov | [END_REF]]. We consider decoupling the label space into two different ones and redefining the topology of the graph, so that we can still capture rigid plane displacements and in-plane deformation.

Decoupled Approach

We propose to overcome the limitations of the overparameterized method by decoupling every node of the previous approach in two different ones: one modeling the in-plane deformation and another the position of the plane. This is somewhat analogous to creating two separated graphs of the same size and topology corresponding to different random variables and label spaces. Once spaces have been decoupled, different sampling strategies can be used for them. Another advantage of this approach is that we can define distinct regularization terms for edges connecting deformation nodes or plane position nodes. It allows to regularize in a different way the deformation and the plane position, imposing alternative geometrical constraints for every case.

Since data term computation requires the exact location of the node, both position and deformation labels are necessary. Both graphs can thus be connected through a pairwise edge between every pair of corresponding nodes. Therefore, new pairwise potentials are associated with these edges in order to encode the matching measure.

Formally, the decoupled formulation consists of an undirected pairwise graph G D = V, E with a set of nodes V = V I ∪ V P and a set of cliques E = E I ∪ E P ∪ E D . V I and V P have the same cardinality and 4-neighbor grid structure. Nodes in V I are labeled with labels that model in-plane deformation, while labels used in V P model the plane position. Edges from E I and E P correspond to classical grid connections for nodes in V I and V P respectively; they are associated with regularization terms. Edges in E D link every node from V I with its corresponding node from V P , creating a graph with a three dimensional structure; those terms encode the matching similarity measure. Note that E I and E P can be extended with the same type of extra edges defined in section 3.2.1 (see figure 3.3.a) to improve the satisfaction of the desired geometrical constraints. The proposed graphical model is depicted in figure 3.6 and corresponds to a pairwise MRF, without unary potentials.

We define two different label spaces, one associated with V I and one associated with V P . The first label space, L I , is a bidimensional space that models in-plane deformation using displacement vectors d l I = (d x , d y ) associated to labels l I . The second label space, L P , indicates the plane in which the corresponding control point is located and consists of labels l P representing different planes. In order to specify the plane and the orientation of the grid on it, we consider an orthonormal basis Figure 3.6: Structure of the decoupled graph. The green nodes are included in V I while the orange ones in V P , modeling in-plane deformations and plane position respectively. Edges connecting V I nodes are part of E I and those connecting V P nodes are part of E P ; they are associated with regularization terms. Dotted lines represent cliques in E D that encode the matching similarity measure. Using this information we can reconstruct an imaginary grid similar to the one presented in section 3.2.1. In the image we can appreciate how we associate two nodes of the graph with one control point of the grid. acting on a reference point in this plane. Using this information, we can reconstruct the position of the control points of the grid. The planes parametrization is given by vectors (φ, θ, λ) associated to labels l P , where angles φ and θ define a vector N over a unit sphere, expressed through its spherical coordinates (see figure 3.3.c). This value, together with parameter λ, defines the position of the plane associated with the given control point. This is an important advantage of our method: we could use prior knowledge to improve the way we explore the plane space, just by changing the plane space sampling method.

As it concerns the considered plane sampling method, the final position of every control point p k of the grid is determined using the pairwise term between two graph nodes (v I k ∈ V I and v P k ∈ V P ) and their respective labels (l I k ∈ L I and l P k ∈ L P ). Imagine we have a plane π k with normal vector N that contains the displaced control point p k +d l I k . Parameter λ indicates the magnitude of the translation we apply to π k in the direction given by N in order to determine the plane's final position (see figure 3.7 for a complete explanation). Given that we can associate different planes to different control points (by assigning them different labels l P ), we need to impose constraints that will force the final solution to refer to a unique plane.

The energy that guides the optimization process involves three different pairwise terms, which encode data consistency between the source and the target, smooth-Chapter 3. Slice-to-Volume Registration ness of the deformation and unique plane selection:

P D (Γ D ; I, P, D) = α (i,j)∈E I e I i,j (l I i , l I j ) + β (i,j)∈E P e P i,j (l P i , l P j ) + (i,j)∈E D e D i,j (l I i , l P j ),
(3.16) where α, β are scaling factors, e I i,j ∈ I are in-plane deformation regularizers (associated to edges in E I ), e P i,j ∈ P are plane consistency constraints (associated with edges in E P ) and e D i,j ∈ D are data terms (associated with edges in E D ). l I i , l P i are labels from label spaces L I and L P respectively, while Γ D is a complete labeling that assigns one label to every variable. Thus, we aim at solving the following MAP inference problem: ΓD = arg min (3.17)

Note that, in this case, we reconstruct the final position of the control points p k as a combination of the labels l I k and l P k assigned by the optimal labeling ΓD . Once the grid is reconstructed, we can apply the same strategy (projecting the control points to the regression plane) to extract both, the ridig transformation and the dense deformation field, as explained at the beggining of section 3.2.

The data term is defined for every control point of the imaginary grid p k using the information provided by two associated graph nodes. It is encoded in the pairwise term e D ∈ E D . To this end, we extract an oriented patch Ω l I k ,l P k over the plane π k (centered at p k ) from the volume J, so that the similarity measure δ can be calculated between that patch and the corresponding area over the source 2D image (see figure 3.7):

e D i,j (l I i , l P j ) = Ω l I i ,l P j δ(I(x), π k [J](x))dx. (3.18)
We define two different regularization terms. The first controls the in-plane deformation; it is defined on V I and corresponds to a symmetric distance preserving penalty:

e I i,j (l I i , l I j ) = (1 -ψ i,j (d l I i , d l I j )) 2 + (1 -ψ i,j (d l I i , d l I j ) -1 ) 2 , (3.19)
where ψ i,j is the distance defined in equation 3.11.

The second term penalizes inconsistencies in terms of plane selection, and is defined on V P . We use the earlier defined (at is concerns the overparameterized model, in equation 3.14) point-to-plane distance:

e P i,j (l P i , l P j ) = 1 2 (D π j (p i ) + D π i (p j )). (3.20)
where p i and p j are the positions after applying label l P i , l P j to p i , p j respectively. Note that these terms are similar to the ones of the former approach. However, there is an important difference regarding the parameters they use. In case of the overparameterized approach, parameters are always 5-dimensional labels. In the current approach, parameters are at most 3-dimensional, thus reducing the complexity of the optimization process while also allowing a denser sampling of the solution space. Conventional pairwise inference algorithms could be used to optimize the objective function corresponding to the previously defined decoupled model. Such a model offers a good compromise between expression power and computational efficiency. However, the pairwise nature of such an approach introduces limited expression power in terms of energy potentials. The smoothness (regularization) terms with second order cliques are not invariant to linear transformations such as rotation and scaling [Glocker 2009c], while being approximate in the sense that plane consistency is imposed in a rather soft manner. These concerns could be partially addressed through a higher order formulation acting directly on the displacements of the 2D grid with 3D deformation labels. Furthermore, the data term is just a local approximation of the real matching score between the deformed source 2D image and the corresponding target plane; by introducing high-order terms we could define it more accurately.

High-Order Approach

The new formulation consists of an undirected graph G H = V, E with a set of nodes V and a set of third-order potentials E = E D ∪ E R . The nodes are control points of our two-dimensional quasi-planar grid and they are displaced using 3D vectors d l i = (d x , d y , d z ) associated to labels l i ∈ L H . We define two types of cliques in E. Cliques in E D are triplets of vertices with a triangular shape and they are associated with data terms. Those in E R , are collinear horizontal and vertical third-order cliques associated with regularization terms.

Unlike the previous methods, which require extra labels to explicitly model the plane selection, high-order potentials explicitly encode them. Furthermore, thirdorder triangular cliques can also explicitly encode data terms, since the corresponding plane can be precisely determined using the position of these 3 vertices. We use triplets of collinear points for regularization terms. According to [Kwon 2008], this allows us to encode a smoothness prior based on the discrete approximation of the second-order derivatives using only the vertices' position. Therefore, we define a simple three dimensional label space of displacement vectors which is sampled as shown in figure 3.3.b.

The energy to be minimized consists of data terms D ijk associated with triangular triplets of graph vertices (i, j, k) ∈ E D and regularization terms R ijk associated with collinear horizontal and vertical triplets (i, j, k) ∈ E R . It is defined in the following equation:

P H (Γ H ; D, R) = (i,j,k)∈E D D ijk (l i , l j , l k ) + γ (i,j,k)∈E R R ijk (l i , l j , l k ), (3.21)
where γ is a scaling factor and l i is a label associated to a displacement vector d l i = (d x , d y , d z ) assigned to the node i. Γ H is a labeling that assigns one label to every variable used to reconstruct the final rigid transformation and dense deformation field. Slice-to-volume registration is thus posed as a MAP inference problem on a high-order graph following: ΓH = arg min

Γ H P H (Γ O ; D, R). (3.22)
The data term is defined over a disjoint set of triangular cliques, covering the entire 2D domain, as shown in figure 3.8.a. Its formulation is independent of the similarity measure δ and it is calculated for each clique c = (i, j, k) ∈ E D using the source 2D image I and the corresponding plane π d [J] extracted from the target volume J, defined by the three control points of the clique. For a given similarity measure δ, the data term associated with the clique c is thus defined as: (3.23) where x ∈ Ω l i ,l j ,l k , and Ω l i ,l j ,l k corresponds to the triangular area defined by the control points of clique c = (i, j, k) over the plane π d [J], after applying the corre- sponding labels l i , l j , l k to the vertices.

D ijk (l i , l j , l k ) = Ω l i ,l j ,l k δ(I(x), π d [J](x))dx,
Smoothness and plane consistency are also imposed using higher order cliques. We define a clique for every set of three collinear and contiguous grid nodes (in horizontal and vertical directions as depicts figure 3.8.b). We also introduce extra cliques formed by nodes that are collinear but not contiguous. The aim is to propagate the regularization so that the planar structure is conserved. The regularization term, as previously, seeks to satisfy the plane structure of the grid and the smoothness nature of the in-plane deformations.

Planar consistency can be easily enforced by propagating a null second-derivative constraint among collinear triplets of points. In fact, a null second-derivative for these cliques does not impose just a planarity constraint but it also aims at regularizing the grid structure. Thanks to the third-order cliques, we can accurately approximate a discrete version of the second-order derivative [Kwon 2008]. Given three contiguous control points p i , p j , p k and their corresponding displacement vec- tors d l i ,d l j ,d l k associated to the labels l i , l j , l k , it can be approximated as follows:

|| (p i + d l i ) + (p k + d l k ) -2 • (p j + d l j ) ||.
Based on this idea, we define the following energy term that is proportional to the second derivative, and normalized with the original distance between the control points, d:

R A ijk (l i , l j , l k ) = || (p i + d l i ) + (p k + d l k ) -2 • (p j + d l j ) || d 2 2 , (3.24)
In-plane deformation smoothness is reinforced in the same manner as the previous models, through a symmetric distance preserving approach. For the sake of clarity, we redefine equation 3.12 as

Ψ ij (l i , l j ) = (1 -ψ i,j (d l i , d l j )) 2 + (1 - Chapter 3. Slice-to-Volume Registration ψ i,j (d l i , d l j ) -1
) 2 , and we apply it to both pairs of contiguous points that form the clique (i, j, k):

R B ijk (l i , l j , l k ) = Ψ ij (d l i , d l j ) + Ψ jk (d l j , d l k ) 2 . (3.25)
The equation that regularizes the grid is a weighted combination of both terms R A ijk and R B ijk :

R ijk (l i , l j , l k ) = (1 -α)R A ijk (l i , l j , l k ) + αR B ijk (l i , l j , l k ), (3.26)
where α represents a weighting factor used to calibrate the regularization term.

Results and Discussion

Let us now proceed with a systematic evaluation of the proposed methods. One of the main aspects shared across methods is the inference algorithms used to produce the desired solution.

Inference Methods

Depending on their cardinality and regularity, objective functions can be optimized using a variety of discrete optimization algorithms which offer different guaranties.

It must be noted that the regularization terms presented in our three models are non submodular, since we include the current position of the points (which can be arbitrary) in their formulation. Therefore, submodularity constraint is fulfilled neither in the pairwise nor in the high-order terms (for a clear definition of submodularity in pairwise and high-order energies, we refer the reader to the work of [START_REF] Ramalingam | [END_REF]).

In [Ferrante 2013], the overparameterized approach was optimized using the FastPD algorithm [Komodakis 2007b] while for the decoupled [Ferrante 2015b] and the higher order models [Ferrante 2015a], we consider loopy belief propagation networks. For the sake of fairness, in order to improve the confidence of the comparison among the three methods, in this work we reimplement it to be optimized with the same algorithms. Therefore, results in this work can not be directly compared with our previous works.

Given the variety of models presented in this work, we choose two different inference methods that can deal with arbitrary graph topologies and clique orders, coming from two standard inference algorithm classes: (i) Loopy Belief Propagation (LBP), a well know message passing algorithm that has been extensively used in the literature; and (ii) the Lazy Flipper (LF) by [START_REF][END_REF]], a movemaking algorithm which is a generalization of the classical Iterated Conditional Modes (ICM) [Besag 1986] and has provided good approximations for several nonsubmodular models in different benchmarks . Both are approximate inference methods that can accommodate arbitrary energy functions, graph topologies and label Figure 3.9: Factor graph derivation and labels spaces corresponding to the overparameterized, decoupled and high-order approaches. It shows the equivalence between cliques c i in the first column (unary, pairwise and high-order, depending on the model) and the corresponding factors f i in the second column. In red, we observe the cliques and factors associated with data terms, while in green and orange we represent those associated with regularization terms. In the third column, we include a figure representing the label space associated to every model (orange vectors and planes are associated to different labels). Note that the overparameterized approach is defined as the Cartesian product between the displacements and plane selection labels, while in the decoupled approach these label spaces are independent.

spaces, and allow us to show how the three proposed approaches perform under different optimization strategies. We refer the reader to sections 2.4.1 and 2.4.2 for a brief description about LBP and LF, respectively.

We adopt the OpenGM2 library [Kappes 2013] which implements both inference methods, and makes it possible to perform fair comparisons. It requires construction of a factor graph for every scheme (we refer the reader to the appendix in section A.1 for a complete definition of factor graphs). Figure 3.9 shows a comparison between the three models and the derivation of the corresponding factor graph in each case. return bestEnergy, G 16: end procedure

Experimental Validation

We compute results on two different datasets for the three methods, using the two inference algorithms (LBP and LF) in order to validate both the resulting 2D-2D deformation field and the final plane estimation. The first dataset is a monomodal MRI heart dataset while the second one corresponds to multimodal US-MRI brain images.

We also run the same experiments using a continuous approach to estimate a rigid transformation model that serves as baseline for comparison. The continuous optimization is performed using the Simplex algorithm [Nelder 1965], which is a gradient-free method that has been successfully used to deal with continuous formulations of slice-to-volume registration, as discussed in section 3.1.4.1. We run the continuous optimization until convergence, which was always achieved in a few seconds. We optimized a global energy where the similarity measure was computed for the complete image, since no local deformation model is considered.

In the three methods, we adopt a pyramidal approach, using different grid resolution levels, from coarse to fine spacing between the control points. For each grid resolution, some iterations of the registration algorithm are performed, choosing the optimal set for each one and updating the control point positions with this information. During the inner iterations of one grid level, the size of the displacement vectors that form the deformation label space as well as the parameter variation of the plane label space are reduced in order to improve the search space sampling. A pseudocode of the algorithm is shown in Algorithm 1. For every registration case, we run the inference algorithm several times (more precisely, the inference method is executed a number of times equal to the product between grid refinement levels and label refinement levels). For a single execution of both discrete inference methods, we use the same compound stopping criterion based on the energy gap between iterations and maximum running time. The algorithms run until the energy improvement between two iterations is smaller than a fraction of the current energy (we use ε = 0.01%) or until a timeout is reached (we use 60 seconds as timeout).

For LF we use a maximum depth of 2 (for details about LF, we refer the reader to section 2.4.2 or to [Kappes 2013]).

In the following subsections we describe the datasets and present quantitative results.

Monomodal Dataset Experiment

The monomodal dataset is derived from a temporal series of 3D heart MRI volumes. It consists of 10 sequences of 19 MRI slices which have to be registered with an initial volume. The slices are extracted from random positions in the volumes while satisfying spatio-temporal consistency. The ground truth associated with this dataset refers to the rigid transformation used to extract every 2D slice of every sequence (it is used to validate the plane estimation or rigid registration) and a segmentation mask of the left endocardium, that can be used to validate the quality of the estimated deformation field.

The dataset is generated from a temporal series of 3D heart MRI volumes M i as shown in figure 3.10. For a given sequence in the heart dataset, every 2D slice I i is extracted from the corresponding volume M i at a position which is calculated as follows. Starting from a random initial translation T 0 = (T x 0 , T y 0 , T z 0 ) and rotation R 0 = (R x 0 , R y 0 , R z 0 ), we extract the first 2D slice I 0 from the initial volume M 0 . Then, gaussian noise is added to every parameter of the transformation in order to generate the position of the next slice at the next volume. We use σ r = 3 • as rotation and σ t = 5mm as translation parameters. Those parameters generate maximum distances of about 25mm between the current and the succeeding plane. In this way, we generate 2D sequences that correspond to trajectories inside the volumes. Since the initial 3D series consists of temporally spaced volumes of the heart, there are local deformations between them due to the heartbeat; therefore, extracted slices are also deformed.

The resolution of the MRI volume is 192 × 192 × 11 voxels and the voxel size is 1.25mm × 1.25mm × 8mm. The slices of the 2D sequences are 120 × 120 pixels with a voxel size of 1.25mm × 1.25mm.

Experiments for the 3 methods were performed using equivalent configurations. In all of them we used 3 grid refinement levels, 4 steps of label refinement per grid level, initial grid size of 40mm and minimum patch size (for similarity measure calculation) of 20mm. In case of the overparameterized approach we used α = 0.8, γ = 1 and 342 labels; for the decoupled approach we used α = 0.8, β = 0.2, 25 Figure 3.10: Heart dataset construction. Given a series of 3D MRI volumes of a beating heart (A), we extract ten different random trajectories (B). Every trajectory is composed of twenty different positions from which we extract the 2D slices (C). labels in the 2D deformation space and 91 in the plane selection space; and finally, for the high-order approach we used α = 0.5, γ = 1.10 and 19 labels. Parameters α, β, γ were chosen using cross-validation. The number of labels in every label space was chosen to make the search spaces as similar as possible.

Results are reported (for every approach and every inference method) for 10 sequences of 19 images, giving a total of 190 registration cases. We use SAD as similarity measure since we are dealing with monomodal registration. The idea is to register every 2D slice I i (which plays the role of an intra-operative image) to the same initial volume M 0 (which acts as the pre-operative image). The resulting position of the slice I i is used to initialize the registration of slice I i+1 .

Figure 3.13 shows results in terms of rigid transformation estimation. We measured the distance between the transformation parameters, and reported the average of the 190 registration cases. It resulted in less than 0.02rad (1.14 • ) for rotation and less than 1.5mm for translation parameters in all the discrete approaches and optimization methods. The decoupled method outperforms the others by orders of magnitude in terms of reduction of the standard deviation and the mean error.

To measure the influence of the deformation in the final results, we use the dataset annotations. We compute statistics for the segmentation overlapping at three different stages: before registration (i.e. between the source image and the target volume slice corresponding to the initial transformation), after rigid registration (i.e. between the source image and the target volume slice corresponding to the estimated transformation) and after deformable registration (i.e. between the deformed source image and the target volume slice corresponding to the estimated transformation). We evaluate accuracy computing DICE coefficient, Hausdorff distance and Contour Mean Distance (CMD). We also provide sensitivity (which measures how many pixels from the reference image are correctly segmented in test image) and specificity (which measures how many pixels outside the reference image are correctly excluded from the test image) coefficients to complete the analysis. Results presented in figure 3.14 show the mean and standard deviation of the indicators at the three stages, for the three approaches and the two inference methods. Observe that results improve at each stage, achieving DICE coefficient of around 0.9 after deformation. Hausdorff distance and CMD decrease at each stage until a total reduction of around 66%. Decoupled method still outperforms the others after deformation in all the indicators, and presents a substantial improvement in terms of standard deviation reduction with respect to them (it is consistent with the results we showed in figure 3.13 for the rigid parameters). All the three discrete methods outperform the continuous rigid baseline. Figure 3.15 complements these Figure 3.12: 12 registration cases of the same sequence, before and after deformable registration. The overlapping images (in light blue we show the source image and in red the target) showed before registration corresponds to the source image and a slice taken from the volume at the initial position. The overlapping after registration corresponds to the deformed source image and the slice taken from the volume at the estimated plane position.

results by showing DICE values disaggregated per sequence, while figures 3.11 and 3.12 shows some qualitative results before and after registration.

Finally, in terms of running time, figure 3.16 presents the average value for the three approaches and the two inference methods, together with the distribution with respect to data cost computation and optimization time. As we can see, the decoupled method again outperforms the other two when inference is performed using LBP. We run all the experiments (brain and heart datasets) on an Intel Xeon W3670 with 6 Cores, 64bits and 16GB of RAM.

Multimodal Experiment

Another dataset was used to test our approaches on multimodal image registration. The dataset consists of a preoperative brain MRI volume (voxel size of 0.5mm × 0.5mm × 0.5mm and resolution of 394 × 466 × 378 voxels) and 6 series of 9 US images extracted from the patient 01 of the database MNI BITE presented in [Mercier 2012]. The intra-operative US images were acquired using the prototype neuronavigation system IBIS NeuroNav. We generated 6 different sequences of 9 2D US images of the brain ventricles, with resolution around 161 × 126 pixels and pixel size of 0.3mm×0.3mm. The brain ventricles were manually segmented in both modalities. The estimated position of the slice n was used to initialize the registration process of slice n + 1. Slice 0 was initialized in a position near the ground truth using the rigid transformation provided together with the dataset. We computed statistics as we did in the previous experiment, but in this case based on the overlap between ventricle segmentations. Since we registered different modalities, we used Figure 3.13: Rigid transformation estimation error for the heart dataset. We measured the distance for every one of the 6 rigid parameters, for the three approaches using LF and LBP as inference methods and for the continuous rigid approach. The discrete methods outperform the results obtained using the baseline. Independently of the discrete inference method, the decoupled approach outperforms the other two in terms of average and standard deviation of the estimated error, for all the 6 parameters.

Mutual Information as similarity measure instead of SAD. 3.14: Segmentation overlapping statistics computed before, after rigid and after deformable registration for both datasets (190 registration cases for the heart dataset and 54 for the brain dataset). In the case of deformable registration, the source segmentation mask was deformed using the estimated deformation field. Results are reported for the continuous approach that only estimates rigid parameters using the Simplex method, and for three discrete methods (overparameterized, high-order and decoupled) using both inference strategies (LBP and LF).

the decoupled method does better in terms of computational time (independently of the inference method), the high-order method achieves better results in terms of segmentation statistics. It must be noted that, in this case, we are dealing with a more complex problem than in the case of monomodal registration; consequently, the increment obtained in terms of accuracy for both, rigid and deformable registration, is smaller. Given that we are dealing with highly challenging images of low resolution being heavily corrupted from speckle, those results are extremely promising. It is known to the medical imaging community that explaining correspondences between different modalities is a difficult task (see section 3.1.2).

In all brain experiments we use initial grid size of 8mm, minimum patch size of Figure 3.15: Final DICE (after deformation) comparison for every sequence (10 sequences in the heart dataset and 6 sequences in the brain dataset). Results are reported for the continuous approach that only estimates rigid parameters using the Simplex method, and for three discrete methods (overparameterized, high-order and decoupled) using both inference strategies (LBP and LF). In case of the heart dataset, decoupled method outperforms the other two in most part of the sequences. In the brain dataset, the high-order approach shows better performance in most cases. This is coherent with the aggregated results shown in Figure 3.14. 13mm, 16 bins to construct the mutual information histograms, 3 grid levels and 4 steps of label refinement per grid level. In case of the overparameterized approach we used α = 0.9, γ = 0.1 and 342 labels; for the decoupled approach α = 0.015, β = 0.135, 25 labels in the 2D deformation space and 91 in the plane selection space; finally, for the high-order approach α = 0.7, γ = 0.05 and 19 labels. Parameters are chosen similarly as in the heart experiments.

Comparative Analysis

In this section, we aim at comparing different aspects of the three approaches we present in this paper, namely label spaces, graph topology and computational time. Without loss of generality, some assumptions are made regarding the models. First, we consider only square grids where N is the number of control points and consequently √ N is the number of nodes per side. Second, for the sake of simplicity we As we can observe, data cost computation represents a bigger portion of the total time in the brain dataset than in the heart dataset. This is due to the similarity measure: while in the monomodal case (heart) we use a simple SAD, in the multimodal case (brain) we need a more complex measure like mutual information. Note that data cost computation time remains constant when we vary the inference method (with small fluctuations due to operating system routines which ran during the experiment) but not across different models. Note that we do not compare with the rigid estimation using Simplex methods, since it estimates only 6 parameters (in a few seconds) while the discrete methods deal with hundreds.

do not consider the extra cliques introduced to improve the geometrical constraints propagation, since they are contemplated as an alternative strategy which may or may not be adopted.

Figure 3.17 shows a comparative analysis between the three approaches, using the two proposed inference methods, in terms of optimization time and final energy. Table 3.1 presents a compendium of the most critical parameters related to the proposed methods. Let us start with the label spaces. We divide them into two types: displacement space (L D ) and plane selection space (L P ). The first one contains the displacement vectors (2D or 3D, depending on the model) applied to the control points, while the second one contains the set of planes that can be chosen. In terms of cardinality of the label spaces, the overparameterized approach has the highest complexity, given by the cartesian product between the displacements and all the possible planes, |L D × L P |. The decoupled model is dominated by the maximum of the cardinality of both label spaces, max(|L D |, |L P |). Finally, for the high-order model it depends only on |L D | since it is not necessary anymore to explicitly model which planes can be chosen -the triangles defined by the triplets of points describe a plane (and even more, a patch on this plane) by themselves. It clearly illustrates how we can reduce the complexity of a given label space by Figure 3.17: Comparison between total optimization time and final energy using two different optimizers (LBP corresponds to circles and LF to crosses). Results are shown for the overparameterized approach (in blue), the high-order approach (in orange) and the decoupled approach (in green). The gray lines connect data points corresponding to the same registration case. As it can be observed with respect to the final energy, both methods are equivalent in general (without considering the outliers). However, there are more important differences in terms of computational time. In the high-order approach, where the label space is small, LF outperforms LBP since convergence is achieved in a few seconds, independently of the dataset. For bigger label spaces (like decoupled and overparameterized approaches), LBP converges faster in case of the heart dataset, where SAD is used as similarity measure and therefore the energy is smooth. The last case is when we use MI as similarity measure (brain dataset) and we have big label spaces: there is no clear pattern in this case. Note that we do not compare with the rigid estimation using Simplex methods, since it estimates only 6 parameters (in a few seconds) while the discrete methods deal with hundreds. making smart decisions in terms of energy definition and graph topology.

However, there is always a trade-off. This strong reduction in the size of the label space, has an effect on other parameters like number of cliques and number of variables. In case of the decoupled model, the main advantage is related to the fact Figure 3.18: Comparison between total optimization time and results accuracy (measured using DICE coefficient) using two different optimizers (LBP corresponds to circles and LF to crosses). Results are shown for the overparameterized approach (in blue), the high-order approach (in orange) and the decoupled approach (in green). The gray lines connect data points corresponding to the same registration case. Note that these results are consistent with those shown in figure 3.17. Indeed, one can observe that graphs in figure 3.17 are essentially a flipped version (over the X axis) of graphs included in this figure. This evidences a high correlation between low energy values and high accuracy of the results, proving that the energy is appropriately modeled.

that while the number of variables and edges augment linearly (it goes from N to 2N in case of variables, and from 2N -2 √ N to 5N -4 √ N in case of pairwise edges), the number of labels decreases quadratically (from

|L D × L P | to max(|L D |, |L P |)).
It results in better performance for the decoupled method as can be observed in figure 3.16. A consequence of the third-order cliques in the high-order method is higher computation costs. Even then, judging from the running times reported in figure 3.16, we achieve good experimental computation time because of the smaller label space.

Finally, we include a comparison in terms of memory footprints (see table 3.2) among the three methods, using two different optimizers. We reported the maximum amount of memory that a process consumed while running one registration case for the heart dataset. As expected, the overparemterized model requires more

Parameter

Overparameterized Decoupled High-order Label space

|L D × L P | max(|L D |, |L P |) |L D | # variables N 2N N # 1st order cliques N - - # 2nd order cliques 2N -2 √ N 5N -4 √ N - # 3rd order cliques - - 4N -6 √ N + 2
Table 3.1: Comparison among the three methods in terms of label space and graph topology.

memory than the other two approaches. Results also suggest that LF is more efficient in terms of memory consumption than LBP, since given the same graphical model, LF always outperforms LBP in terms of memory consumption.

Overparameterized Decoupled High-order LBP 362 MB 37.9 MB 65.4 MB LF 342 MB 34.9 MB 44.7 MB Table 3.2: Memory footprint comparison among the three methods, using two different optimizers. The reported value corresponds to the maximum amount of memory that the process used while running in every case.

Conclusions and Further Ideas

We derive three new models from the standard graph-based deformable registration theory for slice-to-volume registration. We show promising results in a monomodal and a multimodal case, using different inference methods, and we compare them with a baseline rigid approach were inference is performed using continuous optimization. The proposed framework inherits the advantages of graph-based registration theory: modularity with respect to the similarity measure, flexibility to incorporate new types of prior knowledge within the registration process (through new energy terms) and scalability given by its parallelization potential.

The three methods we present aim at optimizing different types of energy functions in order to get both, rigid and deformable transformations that can be applied independently, according to the problem we are trying to solve. An extensive evaluation in terms of different statistical indicators is presented, together with a comparative analysis of the algorithmic and computational complexity of each model. This work constitutes a clear example of the modeling power of graphical models, and it pushes the limits of the state-of-the-art by showing how a new problem can be solved not just in one, but in three different ways.

Numerous future developments built upon the proposed framework can be imagined. Alternative optimization methods and in particular second order methods in the context of higher order inference could improve the quality of the obtained solution while decreasing the computational complexity. The integration of geometric information (landmark correspondences) combined with iconic similarity measures [Sotiras 2010b] could also be an interesting additional component of the registration criterion. Moreover, sensor based technologies (like the ones described in section 3.1.2.3) could be combined with our iconic method to improve its stability when image occlusions or abrupt changes occur. Domain/problem specific parameter learning [Baudin 2013,Komodakis 2015] towards improving the proposed models could have a positive influence on the obtained results. The learning framework that will be presented in section 4.3 follows this direction.

The applications explored in this chapter are mainly related to pre and intraoperative image fusion. However, as it was mentioned in section 3.1.6, other scenarios related to motion correction and volume reconstruction require slice-to-volume registration methods. It would be interesting to explore how our approaches perform in these cases.

The three models we proposed aim at finding the optimal deformation field that relates a 2D slice to a 3D volume. From this deformation field, we extract rigid and deformable parameters as explained in section 3.2. However, it was mentioned that in some specific problems, rigid models are enough since deformations are not considered (see section 3.1.3.1). If we would like to apply our method to these cases, unnecessary efforts would be done to infer the deformation field and then extract the rigid parameters. Therefore, rigid (only) slice-to-volume registration using discrete models constitutes an open research direction. In this regard, the work of [Zikic 2010b] about linear image registration using MRF and discrete methods could be extended to the case of rigid slice-to-volume registration.

Literature Review

Medical images are n-dimensional vectors of intensity values. Context specific/semantic information can augment the expressive power of such images. This information may be encoded in different ways such as data tags (e.g age, gender, weight, disease, etc), anatomical or functional annotations, landmarks, etc. In this chapter, we discuss how such context specific or semantic information can be exploited to improve registration. Various approaches were envisioned to incorporate priors into the registration process. Without being exhaustive, let us present the most commonly adopted approaches towards incorporating context information to the registration process: (i) knowledge-based transformations, where the information is encoded within the deformation model (ii) modality-specific similarity measures, adopted when we know a priori the type of images we want to register (iii) metric learning, where we aim at learning more accurate matching criteria from examples and (iv) segmentationaware strategies, which directly incorporate segmentation priors to the registration process.

Knowledge-based Transformations

Prior knowledge as it concerns the elastic mapping within deformable registration can be achieved either by learning a high dimensional statistical model of deformations or by exploiting biomechanical/biophysical properties of tissues [Sotiras 2013]. In the first case, statistical deformation models are used to capture information about deformation across a population of subjects. Principal Component Analysis (PCA) was used in several studies (e.g [Wang 2000, Tang 2009]) to represent such statistical models. Then, registration can be performed estimating the linear combination of the PCA modes of variation that optimize a given criterion. In [Wouters 2006] the PCA parameters were estimated by minimizing the mutual information similarity criterion. The advantage of such models is that they constrain the set of admissible deformations to be a linear combination of a basis learned from examples. Such methods are manifold constrained since only "observable" mappings can be reproduced. Manifold enhanced methods is an alternative that seek to minimize the distance of the mapping from the space of admissible solutions. One can refer to the work presented in [Glocker 2009b] where the deformation was expressed through a graph-based control-driven method and relative constraints on the deformation of control points were learned. Such constraints were used in the form of costs during the inference of the deformations in new data.

The second alternative, exploiting the same principle, uses biomechanical or biophysical inspired models. Building such models is a challenging task since not only the material properties but also the geometry and boundary conditions of the organs must be considered. Tumor growth models, for example, have been successfully used to guide the estimation of the deformation field [START_REF] Clatz | [END_REF]] through such biphysical constraints. Breast simulation based on Finite Elements Model (FEM) was used by [START_REF] Del Palomar | A finite element model to accurately predict real deformations of the breast[END_REF] to accurately predict real deformations of the breast. Other organs like the prostate have also been considered. In [START_REF][END_REF] a biomechanical model was used to produce admissible simulations of the deformations generated by the TRUS probe insertion. This information is then considered when performing registration between pre-operative and intra-operative prostate images.

The advantage of such approaches relates to the DOF used to represent the transformation. The aforementioned models decrease the number of transformation parameters when compared to the standard dense deformation field. Thus, exploring the search space becomes simpler, leading to improved efficiency and lower computational complexity at the expense of decreasing the expressive power of the mapping through constraints on the admissible transformations.

Modality-Specific Similarity Measure

The use of context-specific metric functions is an alternative constraint to the deformation priors. Clinical settings/applications usually involve known image modalities to be registered. Therefore, an effort was invested to design modality specific similarity measures adapted to them. The simplest approach consists in simulating the target modality from the source (through either a physical simulation system or a machine learning approach), reducing the problem to a monomodal scenario. Simulating one modality from another can be achieved considering the physical properties of the imaging device [Sotiras 2013]. In [Wein 2007, Wein 2008a] a novel method to generate simulated US images from CT data was proposed. They introduce a new similarity measure entitled Linear Correlation of Linear Combination (LC 2 ), which is invariant to missing simulation details, yielding smooth properties and a global maximum at the correct alignment. The same LC 2 was used in [Fuerst 2014] for registering US to MRI images. The similarity measure was defined by locally matching US intensities to both MRI intensity and gradient magnitude. These techniques exploit knowledge relative to specific types of image modalities. Therefore, they achieve accurate results when the underlying physical principles are known and allow the simulation of new modalities. However, such methods fail when generalizing to arbitrary image technologies.

Distance Metric Learning

Metric learning is an alternative way to incorporate context specific and prior information to the similarity measure. It seeks to determine a mapping that will bring samples from similar classes close (corresponding anatomical structures in the different modalities) and far apart the data from different classes [Boukouvalas 2011]. This concept can be applied to different fields using distance metrics, such as classification, clustering and retrieval tasks. A complete survey about distance metric learning was presented in [Yang 2006]. From this general definition, the meaning of metric learning in the context of image registration can be simply interpreted as learning a domain specific matching criterion that allows the comparison of any two given image modalities. Hence, instead of using single universal similarity measures such as MI, SAD or NCC, we learn specific matching criteria optimally adapted for a given task.

In [Yang 2006] methods were classified as supervised or unsupervised. Supervised methods require a set of annotated data (i.e. some sort of similarity information, either in the form of a full scalar distance between samples, a ranking or some proximity of some samples to others or just a separation between samples that are deemed similar and dissimilar [Michel 2014]). Metric learning consists in determining a mapping that will bring all data points within the same classes close, while separating all the data points from different classes. In unsupervised methods (referred as manifold learning approaches in [Yang 2006]) the main idea is to learn an underlying low-dimensional manifold where geometric relationships (e.g. distance) between most of the observed data are preserved. Then, we can map the original data into this simpler representation endowed with a metric that correctly represents the data similarity [Bronstein 2010].

Several image registration methods using supervised metric learning can be found in the recent literature. [START_REF] Lee | [END_REF]] derived a rigid multi-modality registration algorithm where the similarity measure is learned in a discriminative manner, such that the target and correctly deformed source image receive high similarity scores. The learning was performed using a max-margin approach known as structured prediction. The structured prediction approach of [START_REF] Lee | [END_REF]] uses pairs of pre-aligned images as training data. They assume that the similarity of two images decomposes into the similarities of local regions (patches). Therefore, the learning is performed at the patch level. In a different but still supervised setting, [Bronstein 2010] proposed the use of sensitive hashing to learn a multi-modality distance metric, that can be applied to data coming from two different spaces. The idea is to embed the input data from two arbitrary spaces into the binary space where Hamming distance is minimum for corresponding patches and high for the non corresponding ones. This mapping is expressed as a binary classification problem and can be efficiently learned using boosting algorithms. The authors applied this method to perform multimodal non-rigid registration. The training set involved pairs of perfectly aligned images and consisted on a collection of positive and negative pairs of patches. Such a metric learning approach was plugged into the standard graphbased deformable registration framework in [Michel 2011]. In [Toga 2008], instead of requiring pre-aligned images, the training set consisted of non-aligned images with manually annotated patch pairs (landmarks). Coarse-to-fine strategies were adopted where a global similarity measure was learned first, followed by a fine similarity measure that can capture correspondences at the higher resolution image level. Learning was performed by selecting and fusing a large number of image features. More recently, approaches based on deep learning have started to gain popularity. An interesting discussion about this topic was presented in [Long 2014], suggesting that features learned using convnets are at least as useful (and sometimes considerably more useful) than conventional ones when performing alignment of elements of the same class.

Unsupervised learning methods for image registration have also been studied. [START_REF] Wachinger | [END_REF]] proposed to apply manifold learning (through Laplacian eigenmaps) to learn structural representations of multimodal images. The idea was to calculate dense descriptors that represent the structural information of image patches which do not depend on the intensity values of the images but on the structures in the patch. Once the original images were projected to the learned manifold, standard L1 or L2 norm were used for the registration of the structural images. In [Ou 2011] the notion of mutual saliency was considered to locally adapt the metric function involving a weighted sum over a large space of features. A different approach based on unsupervised deep learning was proposed by [START_REF] Wu | [END_REF]]. The composition of a basis of filters was learned to effectively represent the observed image patches. These features were then incorporated in state-of-the-art registration methods and tested on MR images of the brain.

Segmentation-Aware Strategies

Segmentation labels can be an interesting alternative to encode prior information. Given an image I (2D or 3D) and a finite set of semantic classes C, a hard segmentation mask assigns a label c ∈ C to every pixel (voxel) of I. Soft segmentation masks (or probability maps) indicate how likely a given pixel (voxel) is to be associated with a class c ∈ C. We will call segmentation-aware strategies to the methods that incorporate -at some extent-these segmentation maps (hard or soft) directly to the registration process.

A popular way to incorporate segmentation information into the registration process is by performing both tasks at the same time. The idea is to exploit the inter-dependencies between the two problems to improve both, registration and segmentation results. This is known as co-registration and segmentation. A unified segmentation framework integrating tissue classification and image registration within the same generative model was proposed in [Ashburner 2005]. Their model is based on a mixture of Gaussians (MOG) incorporating smooth intensity variation, non-linear registration and soft segmentations. Another example is given in [START_REF] Mahapatra | [END_REF], where segmentation information was integrated to improve MRF-based elastic image registration. They formulated the data term of the MRF energy as a combination of image intensity and mutual dependence of registration and segmentation information. Hence, the labels were defined as the joint occurrence of displacement fields and segmentation class probability. The main disadvantage of this approach is the manual intervention required to define a mask around the object of interest. The area within this mask was considered to get the intensity distributions of the object, which was then used to determine the mutual dependence of the segmentation mask and the displacement vectors. [START_REF] Wyatt | [END_REF]] also showed that one can obtain more accurate and robust results by seeking a joint solution to registration and segmentation. They used a MRF-based approach as well, inducing spatial regularization into the registration process through a class-based entropy measure. Gaussian Mixture Models were coupled with an MRF prior to perform the segmentation task, while the registration relied on the segmentation by minimizing the joint class histogram between both images. In [Parisot 2014], an elegant formulation seeking to take advantage of concurrent segmentation and registration of brain tumor images was proposed. Such a framework used a data-driven adaptive graphical model endowed with uncertainties and was used to adapt the registration task to the presence of the tumor increasing also the segmentation quality. One can also refer to early works like [Yezzi 2003], where a geometric, variational framework that uses active contours to simultaneously segment and register features from two or more images was proposed. They minimize a set of energy functionals, one for each image. Similar work in the context of computer vision and motion analysis and tracking was presented in [Paragios 2005].

Another alternative considering segmentation masks within the registration process is the atlas-based method. In the context of medical image analysis, an atlas is usually defined as the combination of an intensity image (template) and its corresponding semantic annotations (the atlas labels). These annotations can be expressed in the form of discrete segmentation labels or probability maps (probabilistic atlas). The image segmentation problem can be solved using either a single [Cabezas 2011] or multiple atlases [Iglesias 2015]. The idea is to register the atlas(es) to the target image and, then, transfer the labels using the obtained deformation field. In the multi-atlas setting, a label fusion strategy must be considered to generate the final (consensus) segmentation mask. Several works have been presented where simultaneous registration and segmentation are performed, leading to improvements in both tasks by exploiting the inter-dependencies between the two problems. [Gass 2014] introduced a novel MRF-based approach for segmenting medical images while simultaneously registering an atlas non-rigidly. Even if they modeled both type of variables (deformation vectors and segmentation classes) they did it in a sum-label way, i.e. the solution label space is the union of segmentation and registration label spaces and therefore the total number of labels is the sum of the number of registration and segmentation labels (instead of the product as in [START_REF] Mahapatra | [END_REF]). This strategy reduces the complexity of the label space by using a two-layer model where variables in every layer take values from different spaces. The solution of this pairwise model was efficiently found using α-expansion [START_REF] Boykov | [END_REF]]. A variational formulation that simultaneously register an atlas shape to a novel brain image and segments it at the same time was presented in [Wang 2005]. They solved a coupled set of nonlinear PDEs using efficient numerical schemes. The main disadvantage of this approach is that it is limited to binary segmentation problems.

Methods mentioned in the previous paragraph make use of a single atlas. However, multi-atlas approaches where prior information is incorporated to the registration process have also been considered. In [START_REF] Heckemann | [END_REF]], a hierarchical multi-atlas approach was proposed. It introduced tissue-class information to improve image registration and, consequently, segmentation results. The modification used only tissue classification information (learned offline) to guide the three first hierarchy registration levels (using cross-correlation matching criterion to align the tissue probability maps). The resulting transformation was then used as a starting point for detailed registration using the gray intensity values and NMI as the similarity measure. Another multi-atlas framework that incorporates soft segmentation priors was proposed in [Alchatzidis 2014a]. The motivation was based on the idea that prior segmentation probabilities can help registration to avoid local minima, specially in areas of high anatomical variability like brain cortex. They used a MRF-based formulation where each atlas is associated with a deformation field (registration variables), while the target image is associated with a segmentation map (segmentation variables). Prior per voxel probabilities, produced through learning of local features, are taken into account in the energy formulation. Intensity values, smoothness constraints and a coupling term imposing agreement between the segmentation label and the segmentation hypotheses produced from the deformed atlases are considered as well in the energy formulation.

Introduction

In recent years, multi-atlas segmentation methods have become one of the most widely used segmentation techniques for biomedical applications. Multi-atlas based segmentation principle employs an annotated dataset of atlases (images with their corresponding ground truth segmentation), in order to segment an unseen image (target). The atlases are registered to the target image and the resulting deformed candidate segmentations are fused, generating the final segmentation mask. Normally, these approaches allow to segment one target image at a time. In this work, we propose to boost the spatial transformation estimation process by integrating learned soft priors on the structures of interest. To this end, we propose a population deformable registration framework [Sotiras 2009] that, given a set of input volumes, creates an optimal (unknown) reference space through their mutual deformation. We achieve this through a discrete groupwise co-registration process endowed with segmentation consistency constraints which makes possible to segment not one, but multiple images at the same time. This gives rise to a ground truth-free, data-driven registration process that leads to more accurate results than standard multi-volume segmentation strategies.

Multi-atlas segmentation approaches involve two stages: registration and label fusion. Different non-linear approaches have been proposed in the literature with regards to registration [Sotiras 2013]. Methods differ in terms of the deformation model (physical models, interpolation theory or based on specific knowledge of the domain), matching criterion (geometric, iconic or hybrid) and optimization technique (continuous or discrete).

Label fusion is a post processing step, that merges the candidate segmentations obtained after wrapping the atlases to the target image. Such strategies may be simple, like for example majority voting, or more complex involving probabilistic frameworks [START_REF] Rohlfing | [END_REF][START_REF] Sabuncu | A Generative Model for Image Segmentation Based on Label Fusion[END_REF]]. An extensive review of multi-atlas seg-mentation methods can be found in [Iglesias 2015]. A comprehensive evaluation of the accuracy of the atlas-based brain structure segmentation with several different publicly available registration methods was provided in [START_REF] Klein | [END_REF]].

Our method bears some similaries with a standard multi-atlas segmentation algorithm, but targets a different objective. The segmentation results are a byproduct of the population registration where context labels are used to improve the registration performance. Moreover, we aim at segmenting all input volumes simultaneously, while standard approaches do it independently. In the literature, concept similarities with the proposed framework can be found in [START_REF] Heckemann | [END_REF]] and [Alchatzidis 2014a]. These methods use probabilistic priors obtained through a pre-trained classifier to improve segmentation and registration results. Here, opposite to [START_REF] Heckemann | [END_REF]], we do not perform independent prior and intensity based registration steps, but we do it using a compound matching criterion that considers both at the same time. Furthermore, in contrast to the work of [Alchatzidis 2014a] where segmentation variables are explicitly modeled, here we only model registration variables, reducing the number of parameters to be estimated and making the approach tractable and scalable. Moreover, both methods aim at segmenting a single target image, whereas ours considers a population of images simultaneously.

Given a set of unseen images, our method infers the atlases on-the-fly by using previously learned classifiers, and incorporates this information to the energy formulation. We evaluate the effectiveness of our approach on the task of segmenting challenging sub-cortical structures in two brain image datasets. We compare the performance of the proposed method with the standard pairwise atlas-based segmentation, and we evalute it using accurate as well as weak priors. Comparison with ground truth driven co-registration and segmentation is also included for the sake of completeness.

Prior-based Co-Registration and Segmentation

The co-registration and segmentation algorithm is formulated as an energy minimization problem. The input is the set of 3D images D I = {I 1 , I 2 , ..., I N } (with The first term, E I (T ; D I , δ), imposes compactness on appearance of the registered population at the pixel level. It depends on a similarity measure δ and the image intensity of every volume in the input set D I . It is defined as:

I i : Ω ⊂ R 3 → R)
E I (T ; D I , δ) = Ω δ(I 1 • T 1 (x), I 2 • T 2 (x), . . . , I N • T N (x))dx, (4.2)
where δ is a function measuring the coherence between equivalent pixels (or patches centered at given positions, depending on the formulation) in every deformed image.

It is nothing more than the notion of similarity measure, which is extended so that it accounts for more than one image at the same time. Such a term can be either higher order (explicit modeling of global similarity between the deformed volumes), or lower order (decomposition of the global term into a sum of local terms).

The second term, E S (T ; D S , η), measures the degree of agreement among the deformed priors. It depends on the segmentation masks D S and it uses a function η that quantifies this agreement. Hamming distance of the thresholded probability maps or more complex functions could be used in this case. The term is defined as follows:

E S (T ; D S , η) = c∈C Ω η(S 1 • T 1 (x, c), S 2 • T 2 (x, c), . . . , S N • T N (x, c))dx. (4.3)
Note that we sum over all possible semantic labels c ∈ C. If we have priors which are more reliable than others, we could weight them according to the confidence level.

Last but not least, we define a regularization term E R that introduces geometric or anatomical constraints to the deformation fields. In medical image registration, we usually assume that the deformation field is smooth. Different type of regularizers R could be used, which are commonly defined as a convex function acting on the gradient of the deformation field. We apply this constraint independently to every deformation field T i as follows:

E R (T ; R) = N i=1 Ω R(T i (x))dx.
(4.4)

The estimation of the desired transformations can be formulated as a MAP inference problem where we aim at minimizing the following:

T = arg min T E(T ; D I , D S , δ, η, R).
(4.5)

Once the optimal deformation fields T are estimated, we can reconstruct the final segmentation masks DS . To this end, we first propagate every original segmentation prior from D S to the common space, generating warped segmentation masks S i • Ti . Then, say we want to estimate the final segmentation mask for volume k, then we back propagate all the warped segmentation masks S i • Ti from the common space, to the space of image k by using the inverse of the deformation Figure 4.1: Co-registration and segmentation workflow. The input consists of images (I i ) and their corresponding segmentation likelihoods S i . After running the prior co-registration algorithm (which considers the priors in the energy formulation through the semantic consistency term), the optimal deformation fields are used to deform the maximum-a-posteriori of the segmentation likelihoods, which are considered as segmentation hypothesis. A label fusion strategy is then applied to generate the final segmentation mask. The example in this figure illustrates the process to generate the final segmentation for image I 1 . The same strategy is followed to generate the rest of the segmentation.

field. This results in the following equation with respect to the segmentation of the k example:

Ŝk = fuse(S 1 • T1 • T -1 k , S 2 • T2 • T -1 k , . . . , S k , . . . , S N • TN • T -1 k ). (4.6)
In an abuse of notation, S • T indicates that the deformation field T is applied to the initial domain with images intensities being replaced with semantic classes c ∈ C.

The method is modular with respect to the fusion strategy applied to the backpropagated priors to compute the final segmentation mask. We used a simple majority voting, assigning to every voxel the class c ∈ {0, ..., C-1} with the highest number of votes after back-propagation; although simple, this scheme proved to be fairly effective. Figure 4.1 illustrates the complete process. Let us now proceed with a formal definition of the approach in the context of graphical models.

Graphical models and discrete optimization techniques can address quite effectively this type of vision problems. However, note that the energy to be minimized in the MAP problem (defined in equation 4.1) involves high-order dependencies in the data terms E I and E S since all the deformation fields must be considered at the same time. A proper derivation of a discrete MRF using this formulation will produce high-order terms which are difficult to optimize and significantly increase the computational complexity of the method. Alternative models can be derived for this MAP problem, by doing a smarter choice in terms of graph topology, discrete energy formulation and label space. We propose to follow the strategy adopted by [Sotiras 2009], that involved an approximation of the global terms. It consists of assuming that, when estimating a given deformation field T k , the rest of the images do not move within the current iteration. This is somehow analogous to the expansion moves of the well known α-expansion [START_REF] Boykov | [END_REF]], in the sense that we follow the philosophy of the 'very large neighborhood search' techniques [Ahuja 2002]: we start with a feasible solution (in our case, the identity deformation fields) and iteratively try to obtain a better solution that improves the solution of a volumes towards all "static" volumes.

Iterative Prior-based Co-Registration and Segmentation

The co-registration and segmentation algorithm presented in the previous section is reformulated as an iterative process that solves several energy minimization problems defined on independent Markov Random Fields. The step-by-step procedure is described in Algorithm 2. The energy terms E I , E S and E R (intensity matching, semantic prior consistency and regularization terms, respectively) are replaced by E I , E S and E R which consider a single transformation T k at a time (at a given iteration), instead of considering the complete set T . The co-registration process is performed iteratively: only one source image is chosen at each step and registered considering all other images in the set (simultaneously) as targets. To this end, we optimize an objective function E based on a compound criterion where image similarities (E I ), learned priors (E S ) and deformation field regularization (E R ) are considered.

The matching component E I uses again an iconic similarity measure δ. It could be Sum of Absolute Differences (SAD), Mutual Information (MI), Normalized Cross Correlation (NCC) or any other metric that captures similarity between the two images. It is defined as follows:

E I T k ; Ī, I k = N i=0,i =k Ω δ( I i , I k • T k (x))dx. (4.7)
Note that we refer to I i as the latest deformed version of every original input image I i . Given C semantic classes C = {0, ..., C -1}, the segmentation prior component E S takes the following form:

E S T k ; S, S k = c∈C N i=0,i =k Ω η( S i , S k • T k (x, c))dx, (4.8)
where η is the same function measuring the degree of agreement among the deformed segmentation masks (used to define E S in equation 4. Register I k to the other images in the set Ī = {Ii} i∈{1,...,N },i =k , while optimizing E:

Tk = arg min T k EI T k ; Ī, I k + ES T k ; S, S k + ER (T k ) (4.10) 7:
Update DI and DS as:

I k ← I k • T k and S k ← S k • T k 8:
Update deformation field Tk ← Tk • Tk 9:

until All members of DI have been selected as target once 10:

until No more changes in the transformation Ti or the maximum of iterations is reached 11:

Obtain T = T1, T2, ..., TN

12:

for each image I k ∈ DI do 13:

for each segmentation prior Si ∈ DS do 14:

Deform Si to the native space of I k :

S i = Si • T -1 k (4.
E R (T k ) = (N -1) Ω R (T k (x)) dx, (4.9)
where we scale the regularization term by (N -1) for normalization purposes. This iterative process is repeated until convergence, meaning that all images are aligned in a common pose space. The same back-propagation and label fusion strategy defined in equation 4.6 is used to produce the final segmentation masks.

Discrete Formulation

The non-rigid registration between each target image I k and the rest of the images Ī = {I i } i∈{1,...,N },i =k is formulated as a discrete energy minimization problem. We consider Free Form Deformation (FFD) [Rueckert 1999] as our discrete transformation model, which parametrizes the deformation fields T k as a linear combination of K control points. These are regularly placed forming a uniformly distributed grid.

We cast the energy minimization task as a discrete optimization problem of a pairwise discrete MRF. As usual, it is modeled using an undirected graph G = V, X superimposed to the 3D image I k , where V and X denote the nodes and the edges of the graph, respectively. The nodes are interpreted as random variables, while the edges encode the interaction between the variables. The random variables of the MRF model correspond to the displacements of the control points and edges encode dependencies between them. The solution space that can be explored around a control point is given by labels l living in a label space L, where every label l ∈ L is associated to a displacement vector d l . Thus, in the context of discrete optimization, the energy of equation 4.10 becomes:

E MRF (Γ; G, F ) = p∈V g p (l p ) + λ (p,q)∈X
f pq (l p , l q ), (4.12)

where l p , l q are the labels assigned to the nodes p, q ∈ V, λ is a scaling factor and Γ is a labeling that assigns a label l i to every variable in the MRF. The data term g p (l p ) ∈ G is formulated as the combination of the matching and the segmentation prior term:

g p (l p ) = Ωp δ( I i , I k • T d lp k (x))dx + β c∈C Ωp η( S i , S k • T d lp k (x, c))dx, (4.13)
where T

d lp k
is the transformation induced by the movement of the control point p when displacement vector d lp (associated to label l p ) is applied. β is a coefficient, which determines the influence of segmentation priors on the optimization problem. The data term is defined for an image patch Ω p centered at control point p, as it was discussed in section 2.3. Note that for β = 0, the problem reduces to the standard approach; the only difference in this case is that it considers all images in the input set, instead of considering just one target image like in the pairwise case.

The pairwise potentials are defined as the regularization of the deformation field. A discrete approximation of the gradient of the spatial transformation can be computed by taking the difference between the displacements of the neighbouring nodes:

f pq (l p , l q ) = d lp -d lq .
(4.14)

The method is modular with respect to the regularization term. Therefore, different type of regularizers could be considered depending on the problem.

In order to infer the best possible labeling, we employ the Fast-PD [START_REF] Komodakis | [END_REF]] algorithm described in section 2.4.3. It is an efficient discrete optimization method based on linear programming relaxation, that has shown promising results when applied to multi-label problems with similar types of energies.

Experiments and Results

We evaluate the performance of our approach on two medical image datasets by segmenting subcortical structures of the brain. These structures are located below the cerebral cortex and refer to complex shapes difficult to segment. We tested using two different types of priors: reliable priors generated using a convolutional neural network (CNN) and weaker segmentations obtained through random forest (RF) classification.

The first dataset is the publicly available Internet Brain Segmentation Repository (IBSR) [Rohlfing 2012]. It includes 18 annotated 3D T1-weighted MR scans with slice thickness of around 1.3 mm. We use the set of 16 subcortical annota-Figure 4.4: Standard pairwise multi-atlas registration approach. Given a set of atlases (image + ground truth segmentation mask) and a target image to segment, the pairwise multi-atlas registration strategy consists in registering (individually) every atlas image to the target. Once all of them are registered, the resulting deformation fields are used to warp the ground truth segmentation masks to the target space. These deformed masks are considered as segmentation hypotheses. A label fusion step is then applied to generate the final segmentation mask considering the hypotheses. tions, including left and right lateral ventricle, thalamus, caudate, putamen, pallidum, hippocampus, amygdala, and accumbens. The second dataset involves patients with Rolandic Epilepsy (RE), and contains 17 children with epilepsy and 18 matched healthy individuals. T1-weighted magnetic resonance images (MRI) scans were performed by a 3 T scanner (Philips Acheiva) with an in-plane resolution of 256 × 256 and slice thickness of 1 mm. The left and right putamen structures were manually annotated by an expert.

Pairwise Multi-Atlas Baseline. As a baseline, in the first experiment, we implement the standard multi-atlas segmentation based on pairwise registration. In the pairwise approach, all atlases are independently registered to the target image, then the ground truth annotations are deformed to generate the segmentation on the target space (see Figure 4.4). Pairwise registration has been conducted in a discrete fashion using the standard formulation presented in section 2.3, without any prior guidance. Figures 4.2-4.5 show the values for Dice overlap coefficient, Haussdorff distance (HD) and contour mean distance (CMD) between the output segmentation of the pairwise approach and the ground truth on IBSR dataset. Figure 4.6 illustrates the same measurements for the RE dataset. In both cases our approach using CNN priors (Coreg+CNN) performs better than the pairwise segmentation baseline. It should be noted though that the pairwise baseline has a clear advantage with respect to the proposed approach: the segmentation masks that are deformed and used to generate the final segmentation result are the actual ground truth. Instead, the prior based methods (Coreg+CNN and Coreg+RF) do not use at any point the ground truth masks; they only use the estimated probability maps.

CNN Priors (Coreg+CNN)

We train a CNN classifier using the formulation presented in Appendix A.2. To this end, we split the IBSR dataset into three sets. Each time, we use two of the sets as training data and the third set as test data. Similarly, we split RE into two subsets of equal size, using one for training and one for testing, each time. We train on both datasets for 35 epochs starting with a learning rate of 0.01 and dropping it at a logarithmic rate until 0.0001. We train using standard SGD with a momentum of 0.9 and a softmax loss. For all our experiments we use MATLAB and the deep learning library MatConvNet [START_REF] Vedaldi | MatConvNet-convolutional neural networks for MATLAB[END_REF]]1 .

RF Priors (Coreg+RF)

In a second experiment, we compare the performance of the proposed segmentation method using CNN priors, with an approach based on less accurate priors obtained using Random Forest (RF) classification. The RFbased probability maps are computed according to [Alchatzidis 2014b], and the same co-registration process, as explained in 4.2.2, is applied to produce the final ), we achieve results close to the ground truth after co-registration and segmentation in the majority of the cases. Moreover, when we use weak priors like (Coreg+RF), we improve the initial results given by the thresholded priors, specially in terms of CMD and HD.

showed that by using less accurate priors (like those learned using RF), we also obtain good results after co-registration-segmentation as shown in Figures 4.3, 4.2, 4.5 and 4.6. There, we observe significant improvements between RF and Coreg+RF, especially in terms of Hausdorff and contour mean distances. This is due to our formulation, which includes a novel energy term that enables any prior to help during registration. Our method may be used in combination with weak priors learned from weak annotations, such as bounding boxes, requiring less time to be produced than accurate segmentation masks, yet delivering good results.

Our method has considerable advantages over standard multi-atlas segmentation. Given a set of target volumes, multi-atlas segmentation would repeatedly register a set of ground-truth masks and perform label fusion individually for each test volume. Contrary to that, we compute the segmentation probabilities once, and then segment all the volumes "simultaneously". If numerous ground-truth masks were to be used for the registration step, our method benefits from substantial computational gains, as complexity depends only on the number of input volumes we want to segment. In case of large datasets, it removes the burden of selecting an appropriate ground-truth subset to perform multi-atlas segmentation more efficiently; one simply has to compute the probability masks on the input volumes. Another difference is that, in classical multi-atlas, appearance features are used to compute the deformation fields between source and target during the registration process, while no prior on the target image is available. We exploit more sophisticated, learned features to drive the co-registration process. These features are computed for all volumes involved, and are directly related to the desired final output.

Deformable Registration through Learning of Context-Specific Metric Aggregation

In this section we describe the second contribution of this chapter: a novel weakly supervised discriminative algorithm for learning context specific registration metrics as a linear combination of conventional similarity measures. The advantage of combining conventional metrics rather than performing actual metric learning is computational efficiency. Such metrics have been well optimized over the past two decades and therefore both their strengths and limitations are known while efficient implementation strategies have been developed for them as well. The challenge remains to find the optimal relative weighting (or parameters) of different metrics forming the data term. Hand-tuning these parameters would result in sub optimal solutions and quickly become infeasible as the number of metrics increases. Furthermore such hand-crafted combination can only happen at global scale (entire volume) and therefore will not be able to account for the different tissue properties of the observed anatomy. In order to alleviate this problem, we propose a weakly supervised learning algorithm for estimating these parameters locally, conditioned to the data semantic classes. The objective function of our formulation turns out to be a special case of non-convex function, known as the difference of convex function. We use the well known concave convex procedure to obtain the minima or the saddle point of the optimization problem. We show the efficacy of our approach on three highly challenging datasets in the field of medical imaging, for different anatomical structures and image modalities.

Introduction

Deformable registration involves the definition of a similarity criterion that, once endowed with a deformation model and a smoothness constraint, determines (through an optimization strategy) the optimal transformation to align two given images. The definition of this metric function is among the most critical components of the registration process. In that spirit, recent efforts have been made towards learning them from data. A discussion about metric learning for image registration was presented at the beginning of this chapter. As mentioned, machine learning methods, like support vector machines [START_REF] Lee | [END_REF]] or similarity sensitive hashing [Bronstein 2010] have been considered to find mappings between different image modalities. However, these methods rely on the existence of exact massive correspondences at learning which are usually difficult to obtain. Instead, we focus on cases where the ground-truth is given in the form of segmentation masks, which is more common in real scenarios than availability of dense deformation fields or pre-aligned images.

In a number of applications, segmentation maps might be available in one of the volumes to be registered (e.g. adaptive radiotherapy, patient follow up exams or multi-atlas segmentation). Human anatomy is domain specific (deformations as well) and depending on the organ to be imaged one can expect different ade-

Deformable Registration through Learning of Context-Specific Metric Aggregation 103

quacy between the local structure and different mono or multimodal metrics (which somehow reflect the observed tissue properties). Hence, the available segmentation maps could be used to locally suggest efficient similarity measures depending on the organ. The aim of this section is to propose a novel theoretical framework where domain specific optimal combination of metrics is derived according to the clinical task / observed images. The degrees of freedom of the learning procedure are the class weight combination of the different metrics.

It has been empirically shown that the data terms in the energy function have great influence on the accuracy of the solution to the registration problem. Normally, a data term is the function of a metric such as mutual information, sum of absolute differences, normalized cross correlation, etc. A particular data term is thus chosen based on the application. We argue that instead of using one metric, the locally adaptive and content-specific combination of different metrics can further improve the accuracy of the registration task, in particular, in the presence of semantic labels which then make the registration a domain specific problem.

The idea of combining different similarity measures (not necessarily considering them as domain or content-specific) has already been explored. In [START_REF] Cifor | [END_REF]] it was shown that multichannel registration produces more robust registration results when compared to using the features individually. In this case, they used gray intensity value, phase congruency and local phase as features. In a posterior work, [Cifor 2013a] proposed a methodology that does not require to explicitly weight the features, by estimating different deformation fields from each feature independently, and then composing them into a final diffeomorphic transformation. Such a strategy produces multiple deformation models (as many as the number of metrics) that in general are locally inconsistent. Therefore, their combination will not be trivial and, in the general case, not anatomical meaningful.

In deformable registration, it is crucial to choose the right relative weighting between the different metrics and the pairwise smoothness term. One naive way would be to choose these relative weights by cross validation (or hand tuning) over the parameters. Clearly, this approach quickly becomes infeasible as the number of metrics increases while it can be applied only at global scale and cannot accommodate local anatomical differences. In order to circumvent the above mentioned problem, we would like to learn the relative weights from a given training dataset using a learning framework. We propose a novel discriminative learning framework, based on the well known structured support vector machines (ssvm) [Taskar 2003[START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF], to learn the relative weights (or the parameters). The ssvm and its extension to latent models lssvm [START_REF] Yu | [END_REF]] have received considerable attention in the recent years for parameter learning. Few recent examples of such frameworks are [Baudin 2013, Kumar 2015].

Our focus is mainly on the 3D to 3D deformable registration problem where the input and the output images are 3D volumes. However, the same framework can be trivially used for other registration problems as well, as 2D to 2D, or even for slice-to-volume registration. One of the key issues we face in all these problems is that the ground truth deformations are not known. This leads us to adopt a weakly supervised learning framework where we treat the ground truth deformations as the latent variables. We model the latent variable imputation problem as the deformable registration problem with additional constraints. These constraints ensure that when the latent deformations are applied to the source image, the deformed source image is maximally aligned with the target image. The alignment accuracy is measured based on a loss function.

Our learning framework, similar to the lssvm, is a special case of non-convex optimization problems, known as the difference of convex functions. The local optimum or the saddle point of such non-convex function can be obtained efficiently using the well known cccp algorithm [Yuille 2003]. We demonstrate the efficacy of our framework using three challenging medical imaging datasets.

Multi-Metric Deformable Registration

Let us assume that we are given a source 3D volume (or image) I, source 3D segmentation mask S I , and the target 3D volume (or image) J. The size of the segmentation mask is the same as that of the corresponding image. The segmentation mask is formed by the elements (or voxels) s k ∈ C, where C is the set of classes. Without loss of generality, we assume that the elements in the class set C are the discrete variables starting from one.

The problem is formulated using a MRF that consists of a regular grid graph G = V, E , where V is the set of nodes and E is the set of edges. Each node i ∈ V corresponds to a control point p i . Each control point p i is allowed to move in the 3D space, therefore, can be assigned a label l i ∈ L associated to a displacement vector d l i ∈ 3 (in figure 2.2 we illustrate how this label space is built as a discretization of the 3D space). Notice that each 3D displacement vector is a tuple defined as d l i = (d x , d y , d z ), where d x , d y , and d z are the displacements in the x, y, and z directions, respectively. A quantized version of the deformation field D Γ is associated to a labeling Γ of the graph G. In another words, the labeling Γ ∈ L |V | assigns to a each node i ∈ V a label l i associated to a displacement vector d l i . Note that D Γ is a control point based representation of a dense deformation field, which is interpreted as a FFD. We will call DΓ to the dense deformation field obtained as an FFD interpolation of the control point based representation D Γ . We denote the control point p i as the new control point when the displacement d l i is applied to the original control point p i . Let us define a patch ΩI l i on the source image I centered at the displaced control point p i (after applying label l i ). Similarly, we define Ω J i as the patch on the target image J centered at the original control point p i , and ΩS I l i as the patch on the input segmentation mask centered at the displaced control point p i .

Given the definition of these variables, let us define the unary feature vector corresponding to the i th node for a given label l i as U i (l i ; I, J) = (u 1 ( ΩI

l i , Ω J i ), • • • , u n ( ΩI l i , Ω J i )) ∈ R n ,
where n is the number of metrics and u j ( ΩI l i , Ω J i ) is the unary feature corresponding to the j th metric evaluated using the patches ΩI l i and Ω J i . For example, in case the j th metric is the mutual information (MI) then the unary feature In the example, we can observe that the liver (in yellow) and the kidney (in violet) are the dominant classes for the two highlighted control points, respectively. Note that source image (a) is annotated with semantic labels, while the target (b) is a simple intensity image.

u M I ( ΩI l i , Ω J i )
is the mutual information between the patches ΩI l i and Ω J i . In case of single metric, then n is equal to 1. Recall that we have |C| number of classes. Therefore, given a weight matrix W ∈ R n×|C| , where W (i, j) denotes the weight for the i th metric corresponding to the class j, the unary potential for the i th node for a given label l i is computed as follows:

Ūi (l i ; I, J, S I , W ) = w(c) U i (l i ; I, J) ∈ R, (4.15)
where, w(c) ∈ R n is the cth column of the weight matrix W and c is the most dominant class in the patch on the source segmentation mask ΩS I l i . The dominant class c is obtained as follows:

c = arg max c∈C f ( ΩS I l i , c), (4.16)
where, f ( ΩS I l i , c) is the number of voxels of class c in the patch ΩS I l i . Notice that one can use other criteria to find the dominant class.

The pairwise clique potential between the control points p i and p j is defined as f (l i , l j ), where f (•, •) is the L 1 norm between the two displacement vectors d l i , d l j corresponding to the input labels (other regularizers could be used as well, as it was mentioned in section 2.3). Under this setting, the multi-class energy function corresponding to the deformable registration task is defined as:

E(Γ; I, J, S I , W ) = i∈V Ūi (l i ; I, J, S I , W ) + (i,j)∈E f (l i , l j ).
(4.17) Therefore, we aim at finding the optimal labeling Γ (associated to the quantized version of the deformation field D Γ) by solving the following problem: Γ = arg min Γ∈L |V | E(Γ; I, J, S I , W ).

(4.18)

Approximate Inference Algorithm Similar to [START_REF] Glocker | [END_REF]], we adopt a pyramidal approach that allows us to refine the search space at every level and, at the same time, capture a big range of deformations. We use the well known FastPD [Komodakis 2007a] (see section 2.4.3) as the inference algorithm at every level of the pyramid. Notice that the energy function from equation 4.17 is defined over the nodes and the edges of the sparse graph G which represents the deformation field, and not over the dense voxels and the neighbourhood system defined over the input image I. The reason being that the input images are too big and thus can not be optimized efficiently. Once we obtain the optimal labeling Γ, we estimate the dense deformation field DΓ from its quantized representation D Γ using the FFD interpolation model [Rueckert 1999] in order to warp the input image, as it was already mentioned.

Learning the parameters

In the previous section we assumed that the weight matrix W is given to us. However, this assumption becomes unrealistic quickly as the number of metrics and the classes increases. In order to circumvent this problem, we propose an algorithm to learn the weights using a given dataset. Our algorithm is based on the well known latent structured svm framework [Taskar 2003[START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF][START_REF] Yu | [END_REF]] which optimizes an upperbound on the empirical risk. Instead of learning the complete weight matrix at once, we learn the weights for each class c ∈ C individually. From now onwards, the weight vector w c denotes a particular column of the weight matrix W , which represents the weights corresponding to a particular class. We use the words 'parameters' and 'weights' interchangeably. In what follows we talk about the learning algorithm in details.

Preliminaries

Dataset We consider a dataset

D = {(x i , y i )} i=1,••• ,N . Each x i is a pair repre- sented as x i = (I i , J i )
, where I i is the source volume (or the source image) and J i is the target volume (or the target image). Similarly, each y i is a pair represented as y i = (S I i , S J i ), where S I i and S J i are the segmentation masks for the source and target images, respectively. The size of each segmentation mask is the same as that of the corresponding images. As stated earlier, the segmentation mask is formed by the elements (or voxels) s k ∈ C, where C is the set of classes.

The Loss Function The loss function ∆(S I , S J ) ∈ R ≥0 evaluates the similarity between the two segmentation masks S I and S J . Higher value of ∆(., .) implies

Deformable Registration through Learning of Context-Specific Metric Aggregation 107

stronger dissimilarity between the segmentations. Since our final evaluation is based on the dice coeffcient, therefore, we use a dice based loss function. Note that we follow the empirical risk minimization strategy, in which we should optimize the risk over which the evaluation is performed. The loss is therefore defined as:

∆(S I , S J ) = 1 -DICE(S I , S J ). (4.19)
We approximate the dice between the segmentation masks as defined below:

DICE(S I , S J ) = 2 i∈V |φ(p I i ) ∩ φ(p J i )| |φ(p I i )| + |φ(p J i )| , (4.20)
where, φ(p I i ) and φ(p J i ) are the patches at the control point p i on the segmentation masks S I and S J , respectively. The function |.| represents the cardinality of a given set. The above approximation of the dice makes it decomposable over the nodes of the graph G. As will be discussed shortly, this decomposition allows us to train our algorithm very efficiently.

Joint Feature Map Given the parameters w c for a particular class, the labeling Γ, and the input tuple x, the mutli-class energy function (4.17) can be trivially converted into class based energy function as follows:

E c (Γ; x, w) = w c i∈V U i (l i ; x) + w p (i,j)∈E f (l i , l j ). (4.21)
where, w p ∈ R ≥0 is the parameter associated with the pairwise term. Let us denote the parameter vector w ∈ R n+1 as the concatenation of w c and w p . The energy function (4.21) is linear in w and can be written as:

E c (Γ; x, w) = w Ψ(Γ; x). (4.22)
where, Ψ(Γ; x) ∈ R n+1 is the joint feature map defined as:

Ψ(Γ; x) =          i∈V U 1 i (l i ; x) i∈V U 2 i (l i ; x) . . . i∈V U n i (l i ; x) (i,j)∈E f (l i , l j )          (4.23)
Notice that the energy function (4.21) does not depend on the source segmentation mask S I . Source segmentation masks in the energy function (4.17) are used to obtain the dominant class using the equation (4.16), which in this case is not required. However, we will shortly see that the source segmentation mask S I plays a crucial role in the learning algorithm to compute the loss function.

Latent Variables Ideally, the dataset D must contain the ground truth deformation labeling Γ corresponding to the source image I in order to compute the energy term defined in the equation (4.21). Since annotating the dataset with the ground truth deformation is non-trivial, we use them as the latent variables in our algorithm. As will be seen shortly, we impute these deformations using the given dataset ensuring that the loss (as defined in the equation (4.19)) between the source segmentation mask when deformed using the imputed deformation field, and the target segmentation mask is minimized. In what follows, we will refer to this step indistinctly as 'latent assignment step' or 'segmentation consistent registration'.

The Objective Function

Given the dataset D, we would like to learn the parameter vector w such that minimizing the energy function (4.21) leads to a deformation field which when applied to the source segmentation mask gives minimum loss with respect to the target segmentation mask. Let us denote S • DΓ the deformed segmentation when the dense deformation DΓ is applied to the segmentation mask S. Therefore, ideally, we would like to learn w such that:

w * = arg min w 1 N i ∆(S I i • DΓ , S J i ). (4.24)
where, Γ = arg min Γ E(Γ; x i , w). The above objective function is the empirical risk minimization based formulation. However, the objective function is highly non-convex in w, therefore, minimizing it directly makes the algorithm sensitive in terms of convergence to bad local minima. In order to circumvent this problem, we optimize a regularized upper bound on the loss as follows:

min w,{ξ i } 1 2 ||w|| 2 + α||w -w 0 || 2 + C N i ξ i , (4.25) s.t. min Γ,∆(S I i • DΓ ,S J i )=0 w Ψ(x i , Γ) ≤ w Ψ(x i , Γ) -∆(S I i • DΓ , S J i ) + ξ i , ∀ Γ, ∀i (4.26) w p ≥ 0, ξ i ≥ 0, ∀i. (4.27)
The above objective function minimizes an upper bound on the dice based loss function denoted as the variable ξ i , known as the slack. The first term in the objective function ||w|| 2 , is the regularization term used to avoid overfitting. The effect of the regularization term is controlled by the hyper-parameter C. The second term is the proximity term. This ensures that the learned w is close to the initialization w 0 . The effect of the proximity term can be controlled by the hyperparameter α.

Intuitively, for a given input sample, the constraints of the above objective function tries to enforce the condition that the energy corresponding to the best possible deformation field (with minimum loss, additional constraint to enforce coherence) must always be less than the energy corresponding to any other deformation field with a margin proportional to the loss. Notice that, the term

Deformable Registration through Learning of Context-Specific Metric Aggregation 109

Algorithm 3 The cccp Algorithm. 1: D, w 0 , C, α, η, the tolerance ε.

2: t = 0, w t = w 0 .

3: repeat 4:

Segmentation consistent registration. For a given w t , impute the latent variables Γi for each sample by solving the problem (4.29).

5:

Obtain the updated parameters w t+1 by solving the convex optimization problem (4.30). The most violated constraint can be found by solving the problem (4.31).

6: t = t + 1 7: until The objective function of the problem (4.25) does not decrease more than ε.

min Γ,∆(S I i • DΓ ,S J i )=0 w Ψ(x i , Γ
) in equation 4.26, makes the problem non convex. Shortly we will see how to upperbound this term, also known as the latent variable imputation step, in order to make the problem convex.

The Learning Algorithm

The objective function (4.25) that optimizes an upper bound on the empirical risk, is non-convex. Hence, it can not be optimized efficiently to obtain the optimal set of parameters. However, it can be shown that the objective function is a non-convex function that can be re-written as the difference of convex functions, which can be seen as the sum of the convex and the concave functions. For such family of non-convex functions, the well known cccp algorithm [Yuille 2003]) can be used to obtain efficiently a local minima or a saddle point (the cutting plane algorithm can produce the optimal solution, but in practice this is not possible due to the computational constraints associated with the dimension of the problem). Broadly speaking, the cccp algorithm consist of three steps -(1) upperbounding the concave part at a given w, which leads to an affine function in w; (2) optimizing the resultant convex function (sum of convex and affine functions is convex); (3) repeating the above steps until the objective can not be decreased beyond a given tolerance of ε.

The cccp algorithm for the optimization of the objective function (4.25) is shown in the Algorithm 3 (a detailed version of the algorithm is presented in the appendix A.3). The first step of upperbounding the concave functions (Algorithm 3, Line 4) is the same as the latent imputation step, which we call the segmentation consistent registration problem. The second step is the optimization of the resultant convex problem (Algorithm 3, Line 5), which, in this case, is the optimization of the ssvm. The optimization leads to updating the parameters. We use the well known cutting plane algorithm [Joachims 2009] for this purpose. In what follows we discuss these steps in detail.

Segmentation Consistent Registration

As already discussed, the ground truth deformation is not known a priori. Thus, in this step, we generate the best possible ground truth deformation field at a given w. This is same as the latent imputation step of the cccp algorithm. Recall that we are interested in learning the parameters w such that the upper bound on the loss function, defined in equation (4.19), is minimized. This leads us to formulate the latent imputation step as an inference problem with additional constraints. These additional constraints ensure that the imputed deformation field deforms the input image such that the loss between the deformed input image and the target image is minimized. Mathematically, for a given parameter vector w, the latent deformation is imputed by solving the following problem: Γi = arg min

Γ∈L |V | ,∆(S I i • DΓ ,S J i )=0 w Ψ(x i , Γ). (4.28)
The above problem is hard to solve and may not have unique solution. Thus, we solve the relaxed version of the above problem as defined below:

Γi = arg min

Γ∈L |V | w Ψ(x i , Γ) + η∆(S I i • DΓ , S J i ). (4.29)
where, η controls the relaxation trade-off parameter. Since the loss function used is decomposable, the above problem is equivalent to the inference of the deformable registration with simple modifications on the unary potentials. Thus, it can be solved efficiently using the FastPD inference algorithm discussed in the section 2.4.3.

Updating the Parameters Once the latent variables have been imputed or the concave functions have been upperbounded, the resultant objective function can be written as:

min w,{ξ i } 1 2 ||w|| 2 + α||w -w 0 || 2 + C N i ξ i , s.t. w Ψ(x i , Γi ) ≤ w Ψ(x i , Γ) -∆(S I i • DΓ , S J i ) + ξ i , ∀ Γ, ∀i w p ≥ 0, ξ i ≥ 0, ∀i. (4.30)
where, Γi is the labeling associated to the quantized latent deformation field D Γi imputed by solving the problem (4.29). Intuitively, the above objective function tries to learn the parameters w such that the energy corresponding to the imputed deformation field is always less than the energy of any other deformation field with a margin proportional to the loss function with some positive slack. Notice that the above objective function has exponential number of constraints, one for each possible labeling Γ ∈ L |V | . In order to alleviate this problem we use the cutting plane algorithm [Joachims 2009]. Let us briefly talk about the idea behind the cutting plane algorithm. For a given w, each labeling Γ gives a slack. Therefore, instead of minimizing all the slacks for a particular sample at once, we rather find the labeling that leads to the maximum value of the slack and store this in a set known

Deformable Registration through Learning of Context-Specific Metric Aggregation 111

as the working set. This is known as finding the most violated constraint. Therefore, instead of using exponentially many constraints, we now optimize our algorithm over the constraints stored in the working set. This process is repeated till no constraints can be added to the working set. The main ingredient of the above discussed cutting plane algorithm is finding the most violated constraint. As discussed earlier, the most violated constraint for the i th sample is the deformation field associated to the labeling that maximizes the slack corresponding to this sample. Rearranging the terms in the constraints of the objective function (4.30) to obtain the slack, ignoring the constant term w Ψ(x i , Γi ), and maximizing it with respect to the possible deformations (which is equivalent to minimizing the negative of it), leads to the following problem solving which gives the most violated constraint: Γi = arg min

Γ∈L |V | w Ψ(x i , Γ) -∆(S I i • DΓ , S J i ) . (4.31)
Since the loss function is decomposable, again, the above problem is equivalent to the deformable registration with modifications on the unary potentials. Thus, it can be solved efficiently using the FastPD based approximate inference algorithm discussed in the section 2.4.3.

Prediction

Once we obtain the learned parameters w c for each class c ∈ C using the Algorithm 3, we form the matrix W where each column of the matrix represents the learned parameter for a specific class. This learned matrix is then used to solve the registration problem defined in the equation 4.17 using the approximate inference algorithm in a new setting where semantic labels are available in the source image but not in the target (like for example, the scenario of atlas-based segmentation).

Experiments and Results

We evaluated our registration algorithm and the learning framework on three different medical datasets -(1) RT Parotids, (2) RT Abdominal, and (3) IBSR. These datasets involve several anatomical structures, different image modalities, and inter/intra patient images, which makes the deformable registration task on these dataset highly challenging. In all the experiments, we cross validate the hyper parameters C and α, and use η = 50. We use four different metrics in all the experiments: (1) sum of absolute differences (sad), (2) mutual information (mi), (3) normalized cross correlation (ncc), and (4) discrete wavelet coefficients (dwt). In all the experiments (single and multi-metric) we used the same set of parameters for the pyramidal approach based inference algorithm (discussed in the section 4.3.2). These parameters are as follows: 2 pyramid levels, 5 refinement steps per pyramid level, 125 labels, and distance between control points of 25mm in the finer level.

The running time for each registration case is around 12 seconds. For the training, we initialized w 0 with the hand tuned values for each metric: w 0 = (0.1, 10, 10, 10), From second to sixth column, we observe the overlapping between the warped source (in green) and the target (in red) segmentation masks, for the multiweight algorithm (MW) and for the single metric algorithm using sum of absolute differences (SAD), mutual information (MI), normalized cross correlation (NCC) and discrete wavelet transform (DWT) as similarity measure. This is coherent with the numerical results reported in Figures 2 and 3 from the main paper.

for sad, mi, ncc, and dwt, respectively. Below we give details about the different datasets and discuss the results obtained in each of them.

RT Parotids

The first dataset (RT Parotids) contains 8 CT volumes of head, obtained from 4 different patients, 2 volumes per patient. The volumes are captured in two different stages of a radiotherapy treatment in order to estimate the radiation dose. Right and left parotid glands were segmented by the specialists in every volume. The dimensions of the volumes are 56 × 62 × 53 voxels with a physical spacing of 3.45mm, 3.45mm, and 4mm, in x, y, and z axes, respectively. We generated 8 pairs of source and target volumes using the given dataset. Notice that, while generating the source and target pairs, we did not mix the volumes coming from different patients. We splitted the dataset into train and test, and

RT Abdominal

The second dataset (RT Abdominal) contains 5 CT volumes of abdomen for a particular patient captured with a time window of about 7 days during a radiotherapy treatment. Three organs have been manually segmented by the specialists: (1) sigmoid, (2) rectum, and (3) bladder. The dimensions of the volumes are 90 × 60 × 80 voxels with a physical spacing of 3.67mm, 3.67mm, and 4mm, in x, y, and z axes, respectively (there are small variations depending on the volume). We generated a train dataset of 6 pairs and test dataset of 4 pairs. The results on the test dataset are shown in the Figure 4.12, while qualitative results can be found in image 4.9.

IBSR The third dataset (IBSR) is the well known Internet Brain Segmentation Repository dataset, which consists of 18 brain mri volumes, coming from different patients. Segmentations of three different brain structures are provided: white mater (WM), gray mater (GM), and cerebrospinal fluid (CSF). We used a down- sampled version of the dataset to reduce the computation. The dimension of the volumes are 64 × 64 × 64 voxels with a physical spacing of 3.75mm, 3.75mm, and 3mm in x, y, and z axes, respectively. To perform the experiments, we divided the 18 volumes in 2 folds of 9 volumes on each fold. This gave a total of 72 pairs per fold. We used an stochastic approach for the learning process, where we sample 10 different pairs from the training set, and we tested on the 72 pairs of the other fold. We run this experiment 3 times per fold, giving a total of 6 different experiments, with 72 testing samples and 10 training samples randomly chosen. The results on the test dataset are shown in the Figure 4.13.

As evident from the Figures 4.11, 4.12 and 4.13, the linear combination of similarity measures weighted using the learned coefficients outperforms the single metric based registration. In all the cases the dice for the multi-metric is higher than the dice for the single metric based registration, or it is as good as the best dice obtained using the single metric registration (please refer to the Figure 4.12, 'Parotr' to see the case in which the multi-metric is at least as good as the best obtained using the single metric). The results for the 'Sigmoid' organ in the Figure 4.12 show that in some cases the multi-metric based registration can significantly outperform the single metric based registration. Table 4.3.4 shows the average dice value per organ for the three datasets, considering the single and multi-metric approaches. The last column indicates the average dice improvement that our proposed method provides when compared with the single metric version. We achieved maximum average improvement of 0.033, 0.082 and 0.037 in terms of dice coefficient for RT Parotids, RT Abdominal and IBSR respectively. to each dataset. In figures 4.9 and 4.10 we include some qualitative results on the three challenging datasets in order to highlight the effects of learning the weights of different metrics. In the first one (figure 4.9), we present the overlapping of the segmentation masks in different views for one registration case from RT Abdominal and RT Parotids datasets, using single and multi-metric approaches. The observed results are coherent with the numerical results reported in figures 4.11,4.12. We observe that multiweight algorithm gives a better fit between the deformed and ground truth structures than the rest of the single similarity measures, which are over segmenting most of the structures showing a poorer registration performance.

In the second graph (figure 4.10), we include results for the IBSR dataset. According to the scale in the bottom part of the image, extreme values (which mean high differences between the images) correspond to blue and red colors, while green indicates no difference in terms of intensity. Note how most of the big differences observed in the first column (before registration) are reduced in the multiweight algorithm, while some of them (specially in the peripheral area of the head) remain when using single metrics.

Discussion

The main contribution presented in this section is a novel and sophisticated framework for learning aggregations of image similarities in the context of deformable image registration. We also proposed a multi-metric MRF based image registration algorithm that incorporates such metric aggregations by weighting different similarity measures depending on the anatomical regions. We showed that associating different similarity criteria to anatomical regions yields results superior to the classic single metric approach. In order to learn this mapping in real scenarios where ground truth is generally given in the form of segmentation masks, we proposed to conceive deformation fields as latent variables and implement our problem using the lssvm framework. The main limitation of our method is the need of segmentation masks for the source images in testing time. However, different real scenarios like radiation therapy or atlas-based segmentation methods fulfill this condition and can be improved through this technique.

Note that, at prediction (testing) time, the segmentation mask is simply considered to determine the metrics weights combination per control node (as indicated in equation 4.16). The segmentation labels are not used explicitly at testing time to guide the registration process that is purely image based. This is different from the approach presented in section 4.2, where we discussed how segmentation priors can be directly incorporated to registration process in a population-registration scenario. However, as it was mentioned, the co-registration and segmentation method requires a classifier able generate priors for both, source and target images. In the multi-metric approach proposed in this section, segmentation masks are only required (at testing time) for the source image and are only used to choose the best learned metric aggregation. The idea could be further extended to unlabeled data (as it concerns the source image) where the dominant label class per control node is the output of such a classification/learning method.

Deformable Registration through Learning of Context-Specific Metric Aggregation 117

Dataset Organ SAD MI NCC DWT MW Average dice increment for MW Parotl 0,756 0,760 0,750 0,757 0,788 0,033 RT Parotids Parotr 0,813 0,798 0,783 0,774 0,811 0,019 Bladder 0,661 0,643 0,662 0,652 0,736 0,082 Sigmoid 0,429 0,423 0,432 0,426 0,497 0,070 RT Abdominal Rectum 0,613 0,606 0,620 0,617 0,660 0,046 CSF 0,447 0,520 0,543 0,527 0,546 0,037 GM 0,712 0,725 0,735 0,734 0,761 0,035 IBSR WM 0,629 0,658 0,669 0,661 0,682 0,028 Table 4.2: Average dice value per organ, for the single and multi-metric approaches, are reported for the three datasets. The last column indicates the average dice improvement that our proposed method provides when compared with the single metric approaches. We can observe improvements of a maximum of 8% points in terms of dice coefficient.

Conclusions and further ideas

In this chapter we explored two different strategies to incorporate prior knowledge to the registration process. In section 4.2 we proposed a novel method for co-registration and segmentation of a population of images, guided by semantic label likelihoods. We evaluated two different types of priors: probability maps obtained from a deep CNN and weaker annotations obtained with standard RF. Our approach aimed to: i) infer deformations that are anatomically plausible (smoothness constraints); ii) establish visual consistencies between all volumes according to an arbitrary metric; iii) enforce segmentation consistencies among all volumes, according to the predicted likelihoods. Experimental evaluation on conventionally adopted benchmarks, as well as on an additional clinical problem, demonstrates the effectiveness of our approach. Our experiments also show the importance of reliable segmentation priors. Label likelihoods extracted with a deep CNN outperform alternative methods and can replace ground truth annotations in co-registration with minimum loss in performance. In this context, future research directions include formulating the problem using high-order models to allow simultaneous estimation of all the deformation fields and to avoid the iterative approach. The use of decomposed optimization methods allowing efficient, one-shot optimization [Komodakis 2011] is under investigation.

Our second contribution in this chapter was a novel, generic, context-specific registration algorithm that is based on the principle of learning the combination of different metrics. Rather than employing highly customized solutions that suffer from scalability, portability and modularity, our approach relies on producing class-dependent metrics as linear combination of widely known and conventional monomodal and multimodal metrics. Consequently, the proposed registration method is scalable and modular and can be adjusted in any setting by simply changing the linear weights. These weights could be learned trivially offline using partially annotated data or could easily integrate new metrics. Extensive experimental validation on various challenging datasets demonstrated the potentials of the proposed method. The registration weights were learned by minimizing an upperbound on the dice based loss function. Dice is conventionally used but does not offer a very convincing picture as it concerns registration performance. The integration of alternative accuracy measures such as the Hausdorff distance between surfaces or even real geometric distances for anatomical landmarks could further enhance the performance and the robustness of the method. The use of alternative parameter learning mechanisms [Komodakis 2015] is also another interesting approach to explore.

Last but not least, an important research direction that involves the two contributions is testing the proposed methods on other real clinical problems, where co-registration and segmentation as well as context specific metrics are of great importance, such as adaptive radiotherapy.

Chapter 5

Conclusions and future directions

Contributions

In this thesis we propose several contributions to the theory of graph-based deformable registration. The first one (presented in chapter 3), focuses on extending it to a new application scenario not yet explored. We present three different models to deal with the challenging problem of slice-to-volume registration. We define low-rank and high-order graphical models, that simultaneously select the plane and estimates the in-plane deformation through a single shot optimization approach. We discuss the peculiarities of each model in a comparative study, emphasizing the robustness of this theory and illustrating the always underlying compromise between model complexity and power of expression. While exploring different approaches, we transfer the complexity from the label space to the topology of the graph, and we also experiment with high-order terms. Different types of discrete optimizers are used to perform parameter inference on two different datasets of medical images. We also compare the quality of our results with a standard rigid approach based on Simplex optimization of a global similarity measure. The results show that, even in a problem which is continuous by nature like slice-to-volume deformable registration, we can achieve equivalent performance by defining the label space in a smart way and exploring it using iterative large neighborhood search strategies.

Two different strategies incorporating prior information to the standard graphbased registration process are proposed in chapter 4. Both approaches demonstrate that semantic and context specific information can boost the registration algorithms by further improving the accuracy of the results. The first one (described in section 4.2) shows how prior information about the anatomical structures of a given set of images, can be encompassed within the energy function of a discrete population registration algorithm. We seek at performing joint registration/segmentation of the input images while exploiting context information provided by pre-trained classifiers. We experiment with both reliable as well as weaker priors. The first were obtained using a CNN classifier, while for the second we used random forests. The experimental validation shows that by using accurate segmentation priors we can obtain results close to the ground truth. Moreover, we also show that less accurate priors will also lead to acceptable results after performing co-registrationsegmentation. The second contribution (discussed in section 4.3), adapts the well known lssvm framework to learn content specific aggregations of metrics from anno-tated data. We implemented this idea within a multi-metric version of the standard graph-based deformable registration algorithm, and we evaluate its performance on three datasets involving several organs, different image modalities, and inter/intra patient images, which make the deformable registration task highly challenging. In all these cases, the multi-metric approach produced more robust registration results than using the features individually. One of the key issues we faced when designing this approach, is that the ground truth deformations are not known a priori. Instead, annotations usually come in the form of hard segmentation or probability masks. This led us to adopt a weakly supervised learning framework where we treat the ground truth deformations as the latent variables. In such a way, our proposal is applicable to a wider range of problems than those methods requiring ground truth deformations.

Future directions

The contributions presented in this dissertation push the limits of the state-of-theart methods concerning medical image registration. However, they are only small steps towards improving this challenging task in different directions where research has still to be done. To conclude this work, let us discuss some envisioned future research directions.

Boosting context specific methods

The advent of accurate classifiers efficiently augmenting the information provided by standard image modalities, rises a number of interesting questions regarding image registration, which are only partially answered in this thesis. Can this information be used to further improve registration results? Would it be possible to directly learn the deformation fields from segmentation masks, instead of conceiving them as latent variables in this process? Could we benefit from the results of the coregistration and segmentation process to boost the classifiers used to create the initial priors? All these open questions define new directions in which research can be conducted.

Combining our two context specific strategies might be possible. If we think about the co-registration and segmentation strategy, considering that we have the inferred semantic labels at our disposal, the proposed approach can support classspecific combinations of metrics that could be learned offline using the learning approach from section 4.3. Once the metric aggregations are learned, incorporating a multi-metric data term within the co-registration-segmentation approach would be simple since the inferred priors could be used as part of equation 4.16 to find the dominant classes. In such a way, having a classifier alleviates the main limitations of our multi-metric registration algorithm, i.e. the requirement of segmentation masks for the source image at testing time. The long-term goal would be to study the gains of combining different metrics per class and using them as content-adaptive potentials in the registration/segmentation energy.

Future directions 121

Developing novel high-order models

Other interesting open problem is related to the use of high-order graphs and inference methods especially adapted to these cases. As it was shown, the power of expression of graphical models grows together with their complexity. However, different types of high-order models can be efficiently optimized [Komodakis 2009,Dokania 2015]. In section 3.2.3 we present a high-order model to solve slice-to-volume registration which is optimized using standard inference methods, whereas in section 4.2 we formulate prior based co-registration and segmentation as a high-order problem, which is then approximated through a pairwise formulation. In both cases more accurate results could be obtained if proper high-order inference adapted to the particular type of energies would be considered. Accordingly, dual decomposition methods [Komodakis 2009] could be contemplated as an alternative solution.

Extending and improving graph-based slice-to-volume registration

Numerous future developments built upon the proposed slice-to-volume registration framework can be envisioned. New hybrid data terms integrating geometric information, iconic similarity measures [Sotiras 2010b] and sensor based technologies (like the ones described in section 3.1.2.3) may lead to improvements in terms of accuracy and stability when image occlusions or abrupt changes occur. Domain/problem specific parameter learning [Baudin 2013, Komodakis 2015] towards improving the proposed models could have a positive influence on the obtained results. The learning framework presented in section 4.3 could be simply adapted to perform this task.

Slice-to-volume applications explored in this thesis are mainly related to pre and intra-operative image fusion. However, as it was mentioned in section 3.1.6, other scenarios like motion correction and volume reconstruction require slice-tovolume registration methods. It would be interesting to explore how our approaches performs in these cases. Moreover, a single slice at a time is considered in our formulation. However, as discussed in section 3.1.5, employing a multi slice-to-volume procedure -as opposed to single slice procedure-certainly improves the robustness of the registration process by augmenting the image support. This is even more interesting when the relative position between the slices is known a priori (e.g. the works of [START_REF] Gill | [END_REF], Heldmann 2009, Cifor 2013b, Yu 2011]). Our single slice-tovolume registration framework could be directly extended to the multi slice scenario by performing independent registrations for each slice. However, the prior knowledge about the relative positions or interaction dynamics among the slices as well as the wider image support would be ignored. A possible extension considering one grid-graph per slice (modeling in-plane deformations and plane selection parameters) and inter-graph connections (encoding inter-slice dependencies) could do the job.

Considering discrete rigid slice-to-volume registration where only 6-DOF are estimated, could also be of interest. The three models we proposed aim at perform-software libraries, have rekindled the interest in deep learning and have made it possible to harness their power. Krizhevsky et al. published [Krizhevsky 2012] results demonstrating clear superiority of deep architectures over hand-crafted features or shallow networks, for the task of image classification. Since then, CNNs have helped set new performance records for many other tasks; object detection [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF]], texture recognition [START_REF] Cimpoi | Deep filter banks for texture recognition and segmentation[END_REF]] and object semantic segmentation [Chen 2014] just to name a few.

The CNN used in this work is similar in spirit to [START_REF] Prasoon | Deep feature learning for knee cartilage segmentation using a triplanar convolutional neural network[END_REF]], but with some notable differences. In [START_REF] Prasoon | Deep feature learning for knee cartilage segmentation using a triplanar convolutional neural network[END_REF]] the authors train one CNN for each of the three orthogonal views of MRI scans, for knee cartilage segmentation, with the loss being computed on the concatenated outputs of the three networks. The inputs to each CNN are 28 × 28 image patches and the output is a softmax probability of the central pixel belonging to the tibial articular cartilage. In contrast, our method operates on full 2D image slices, exploiting context information to accurately segment regions of interest in the brain. In addition, we use fully convolutional CNNs [START_REF] Long | [END_REF]] to construct dense segmentation maps for the whole image, instead of classifying individual patches. Furthermore, our method handles multiple class labels instead of delivering a foreground-background segmentation, and it does that efficiently, performing a single forward pass in 5ms.

We build a CNN for semantic segmentation, using ideas that were recently introduced for the Deeplab model by Chen et al [Chen 2014]. Deeplab is a fully convolutional neural network that computes a dense labelling of the input image. Compared to other systems, such as the one proposed by Long et al [START_REF] Long | [END_REF]], Deeplab avoids excessive downsampling of the feature maps by skipping the subsampling steps in the last convolutional layers. In order to keep the same effective receptive field size, the authors introduce holes in a spirit similar to the "atrous" (with holes) algorithm [Mallat 1999]. For more details, we refer to [Chen 2014].

In our case we are interested in computing accurate priors for the location of sub-cortical regions in the human brain. Employing the hole trick, allows us to retain a reasonable output resolution, increasing the probability of capturing these fine structures. For training, the straightforward way would be to simply fine-tune the last three convolutional layers in the original Deeplab network, starting from the VGG-16 initialization for the remaining weights. However, given the very different nature of natural RGB images and MR image data (RGB vs grayscale, varying vs black background), we decided to train a fully convolutional network from scratch.

Our choice is not free of technical hurdles. Medical image datasets are usually smaller in size compared to natural image datasets, and high-quality ground truth annotations for segmentation are hard and costly to obtain. In our case, we only have a few 3D scans at our disposal, which increases the risk of overfitting. As we mentioned, the repeated pooling and sub-sampling steps that are applied to the input image, decrease the output resolution, making it difficult to detect and segment finer structures in the human brain. We attempt to overcome these hurdles through a series of design choices for the network. First, we use a shallow network, compared to VGG-16/Deeplab, composed of five pairs of convolutional/max pooling

A.3 cccp algorithm (detailed version)

Here we include a detailed version of the cccp algorithm used to solve the learning problem of section 4.3 (see algorithm 3).

Algorithm 4

The cccp algorithm (detailed version).

1: D, w 0 , C, α, η, the tolerance ε. 2: t ← 0.

3: w t ← w 0 4: repeat 5:

For a given w t , impute the latent variables Γi for each sample by solving the problem: Γi = arg min Γ∈L |V | w t Ψ(x i , Γ) + η∆(S I i • DΓ , S J i ) .

6:

Initialize the constraint set for each sample: W i ← ∅, ∀i. w Ψ(x i , Γ) -∆(S I i • D Γi , S J i ) .

9:

Update constraint set if Γi is sufficiently violated.

W i ← W i ∪ Γi , ∀i. , vol. 21, no. 11, pages 977-1000, 2003. (Cited on page 3.)

  d'une tranche avec un volume. Le recalage d'une tranche avec un volume est le recalage 2D dans un volume 3D, comme par exemple le mapping d'une tranche tomographique dans un système de coordonnées 3D d'un volume en particulier. Nos avons proposé une formulation scalable, modulaire et flexible pour accommoder des termes d'ordre élevé et de rang bas, qui peut sélectionner le plan et estimer la déformation dans le plan de manière simultanée par une seule approche d'optimisation. Le cadre proposé est instancié en différentes variantes, basés sur différentes topologies du graph, définitions de l'espace des étiquettes et constructions de l'énergie. Le potentiel de notre méthode a été démontré sur des données réelles ainsi que des données simulées dans le cadre d'une résonance magnétique d'ultrason (où le cadre d'installation et les stratégies d'optimisation ont été considérés).

2. 1

 1 Overlapping between the source and target images before and after deformable registration. . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Different type of sampling strategies for the discrete label space used for image registration . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Deformable registration as a discrete graph labeling problem. . . . . 3.1 Basic workflow to perform slice-to-volume registration based on graphical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Recovering the rigid and deformable parameters. . . . . . . . . . . . 3.3 Connectivity structure of the graph for a grid and overparameterized label space interpretation . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Data term formulation for the slice-to-volume overparameterized approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Regularization term formulation for the slice-to-volume overparameterized approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Structure of the graphical model in the slice-to-volume decoupled approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Data term formulation for the slice-to-volume decoupled approach . 3.8 Different types of cliques used in the high-order formulation of sliceto-volume registration . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Factor graph derivation and labels spaces corresponding to the overparameterized, decoupled and high-order approaches. . . . . . . . . . 3.10 Heart dataset construction . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Slices extracted from three different sequences of the heart dataset before and after registration . . . . . . . . . . . . . . . . . . . . . . . 3.12 12 registration cases of the same heart sequence, before and after deformable registration . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Rigid transformation estimation error for the heart dataset . . . . . 3.14 Segmentation overlapping statistics computed before, after rigid and after deformable registration for heart and brain datasets . . . . . . 3.15 Final DICE (after deformation) comparison for every sequence (10 sequences in the heart dataset and 6 sequences in the brain dataset). 3.16 Average running time for the three slice-to-volume registration approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Comparison between total optimization time and final energy using two different optimizers . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Comparison between total optimization time and results accuracy (measured using DICE coefficient) using two different optimizers . . 3.19 Visual results for one slice from four of the six brain sequences (each row correspond to a different sequence) . . . . . . . . . . . . . . . . Chapter 1

Figure 2 . 2 :

 22 Figure 2.2: Different type of sampling strategies for the discrete label space used for image registration: sparse (a) and dense (b). Note that every label in this formulation corresponds to a displacement vector (in orange).

Figure 2 . 3 :

 23 Figure 2.3: Deformable registration is formulated as a discrete labeling problem on a grid-graph, superimposed to the source 2D image. The nodes of the graph correspond to the variables of a Markov Random Field, and model the control points of a FFD. The space of plausible deformations is explored by choosing different labelings.

  of Slice-to-Volume Registration . . . . . . . . . . . 25

Figure 3 . 1 :

 31 Figure 3.1: Basic workflow to perform slice-to-volume registration based on graphical models. (1) A 2D input image I and a 3D target volume J are given as input data. (2) A grid is superimposed to image I. The process is initialized using a 6-DOF rigid transformation T 0 that specifies the initial position of the grid within the volume J. (3) The grid is deformed by optimizing an energy function. (4) The plane π and the deformation field TD are reconstructed from the final state of the optimized grid. (5) TD is used to deform image I, and it is provided as output together with the corresponding slice π[J].

Figure 3

 3 Figure 3.3: (a) Connectivity structure of the graph for a grid of size 5x5. The gray edges are standard 4-neighbor connections while the orange ones correspond to the extra cliques introduced to improve the geometrical constraints propagation. (b) Displacement vectors corresponding to the first three elements of a label from the overparameterized approach d l i = (d x , d y , d z ). (c) Unit vectors in spherical coordinates corresponding to the last two coordinates of a label from the overparameterized approach N l i = (φ, θ). (d) Displacement of the control points p i and p j when the corresponding labels l i associated to (d l i , N l i ) and l j associated to (d l j , N l j ) are applied. The planes π i and π j are those that contain the control points p i + d l i , p j + d l j and whose normals are N l i , N l j respectively.

Figure 3

 3 Figure 3.4: Data term formulation for the overparameterized approach. The points x ∈ Ω l i are used to calculate the unary potential. π[J](x) returns the intensity of the point in the 2D slice corresponding to the plane π i in the 3D image, whereas I(x) returns the 2D image intensity. δ represents the similarity measure.

Figure 3

 3 Figure 3.5: (a) In plane regularization term: the dotted line represents the distance used in F 1 , i.e. the distance between the points assuming they are coplanar. (b) Plane structure regularization term: the dotted line represents the distance between one of the control points and the plane corresponding to the other one. This information is used to compute the term F 2 .

  Γ D P D (Γ D ; I, P, D).

Figure 3

 3 Figure 3.7: Data term formulation for the decoupled approach. It is similar to the formulation shown in figure 3.4, but it combines labels from different label spaces. The points x ∈ Ω l I k ,l P k are used to calculate the unary potential. π k [J](x) returns the intensity of the point in the 2D slice corresponding to the plane π k in the 3D image, whereas I(x) returns the 2D image intensity. δ represents the similarity measure. In order to compute the final position of the sampled patch in the volume, the displacement d l I k = (d x , d y ) associated to the in-plane deformation label l I k is applied to the corresponding imaginary grid point p k . Then, label tuple (N, λ) associated to the plane selection label l P k is used: the point is translated in the direction given by vector N as indicates scalar λ. Finally, the patch Ω l I k ,l P k is sampled from plane π k with normal N , centered at the displaced point p k (in orange).

Figure 3

 3 Figure 3.8: Different types of cliques used in the formulation. (a) Example of a triangular clique used for data term computation. The green patch Ω corresponds to the clique (i, j, k) and it is used to calculate the data term. (b) Examples of vertical (i 1 , j 1 , k 1 ) and horizontal (i 2 , j 2 , k 2 ) collinear third-order cliques used to regularize the grid structure.

  Pseudo-code corresponding to the pyramidal approach based sliceto-volume registration 1: procedure 2D3DRegistration(I: Source,J:Target,T 0 :Initial guess)

Figure 3 .

 3 Figure 3.11: Slices extracted from three different sequences of the heart dataset before and after registration. Input slice (a) is initialized in a position specified by a rigid transformation within the volume, whose slice corresponds to (b). After deformable registration, a deformation field (c) is estimated. (d) shows the difference between initial images (a) and (b), while (e) shows the difference between (a) and the corresponding slice extracted after rigid registration. Finally, (f) corresponds to the results after deformable registration (i.e., the difference between the deformed version of slice (a) and the slice corresponding to the estimated transformation). Red indicates bigger differences between the images. Note how these values are changing before (d), after rigid (e) and after deformable (f) registration.

Figure 3 .

 3 Figure 3.14 summarizes the average DICE, specificity, sensibility, Hausdorff distance and Contour Mean Distance coefficients for all the series, while figure 3.16 reports the running times. Figure 3.15 complements these results by showing DICE values disaggregated per sequence. Figure 3.19 show some qualitative results. While

Figure

  Figure3.14: Segmentation overlapping statistics computed before, after rigid and after deformable registration for both datasets (190 registration cases for the heart dataset and 54 for the brain dataset). In the case of deformable registration, the source segmentation mask was deformed using the estimated deformation field. Results are reported for the continuous approach that only estimates rigid parameters using the Simplex method, and for three discrete methods (overparameterized, high-order and decoupled) using both inference strategies (LBP and LF).

Figure 3 .

 3 Figure3.16: Average running time expressed in seconds, for one registration case, for the three approaches running on the heart dataset (a) and the brain dataset (b), using LF and LBP. Blue part corresponds to data cost computation while orange part corresponds to the optimization time. As we can observe, data cost computation represents a bigger portion of the total time in the brain dataset than in the heart dataset. This is due to the similarity measure: while in the monomodal case (heart) we use a simple SAD, in the multimodal case (brain) we need a more complex measure like mutual information. Note that data cost computation time remains constant when we vary the inference method (with small fluctuations due to operating system routines which ran during the experiment) but not across different models. Note that we do not compare with the rigid estimation using Simplex methods, since it estimates only 6 parameters (in a few seconds) while the discrete methods deal with hundreds.

Figure 3 .

 3 Figure 3.19: Results for one slice from four of the six brain sequences (each row correspond to a different sequence). (a) Source 2D US image. (b) Slice extracted from the MRI corresponding to the initial position of the plane. (c) Deformed source image overlapped with the estimated deformation field. (d) Blending between initial images (US and corresponding MRI slice). (e) Blending between final images (Deformed US image and estimated MRI slice). (f) Overlapping between initial segmentations. (g) Overlapping between segmentations after registration.

  and their corresponding segmentation priors (soft or hard) D S = {S 1 , S 2 , ..., S N } which are associated to the possible segmentation classes c ∈ C = {0, ..., C-1} as S i : Ω × C → {0, 1}. The desired output is a set of final multi-label segmentation masks DS = Ŝ1 , Ŝ2 , ..., ŜN together with the deformation fields T = T1 , T2 , ..., TN which warp every image to a common coordinate space. We optimize an objective function E based on a compound matching criterion where both image similarities and learned priors are considered as well as deformation smoothness for the deformation fields. Let us now, without loss of generality, proceed with a concept definition of our approach. The objective function that could be used to solve such a problem involves three terms: E(T ; D I , D S , δ, η, R) = E I (T ; D I , δ) + E S (T ; D S , η) + E R (T ; R) .

Figure 4

 4 Figure 4.2: Average Dice coefficient of subcortical structures in IBSR dataset.

Figure 4 . 3 :

 43 Figure 4.3: Average Hausdorff distance of subcortical structures in IBSR dataset.

Figure 4 . 5 :

 45 Figure 4.5: Average contour mean distance on sub-cortical structures of IBSR dataset.

Figure 4

 4 Figure 4.6: Average Dice coefficient, Hausdorff distance, and contour mean distance for the left and right putamen structures in RE dataset.

Figure 4 . 7 :

 47 Figure 4.7: Segmentation results for three different views (axial, coronal, and sagittal). Coreg-CNN and Coreg-GT (oracle) generate more accurate segmentations compared to pairwise and RF-based approaches. Coreg-CNN can be used as a reliable substitute for ground truth annotations in multi-atlas co-registrationsegmentation (better viewed in color and magnified).

Figure 4

 4 Figure 4.8: Context-specific metric aggregations. The multi-metric deformable registration algorithm uses a different aggregation of metrics depending on the context. It is determined by the dominant class in the corresponding image support area.In the example, we can observe that the liver (in yellow) and the kidney (in violet) are the dominant classes for the two highlighted control points, respectively. Note that source image (a) is annotated with semantic labels, while the target (b) is a simple intensity image.

Figure 4 . 9 :

 49 Figure4.9: Overlapping of the segmentation masks in different views for one registration case from RT Abdominal (first and second rows) and RT Parotids (third and fourth rows) datasets. The first column corresponds to the overlapping before registration between the source (in blue) and target (in red) segmentation masks of the different anatomical structures of both datasets. From second to sixth column, we observe the overlapping between the warped source (in green) and the target (in red) segmentation masks, for the multiweight algorithm (MW) and for the single metric algorithm using sum of absolute differences (SAD), mutual information (MI), normalized cross correlation (NCC) and discrete wavelet transform (DWT) as similarity measure. This is coherent with the numerical results reported in Figures 2 and 3 from the main paper.

Figure 4

 4 Figure 4.10: Qualitative results for one slice of one registration case from IBSR dataset.In this case, since showing overlapped structures in the same image is too ambiguous given that the segmentation masks almost cover the complete image, we are showing the intensity difference between the two volumes. This is possible since images are coming from the same modality and they are normalized. The first column shows the difference of the original volumes before registration. From second to sixth column we observe the difference between the warped source and the target images, for the multiweight algorithm (MW) and the single metric algorithm using sum of absolute differences (SAD), mutual information (MI), normalized cross crorrelation (NCC) and discrete wavelet transform (DWT) as similarity measure.

Figure 4 .

 4 Figure 4.11: Results for the RT parotids dataset for the single-metric registration (sad, mi, ncc, dwt) and the multi-metric registration (mw). The weights for the multi-metric registration are learned using the framework proposed in this work. 'Parotl' and 'Parotr' are the left and the right organs. The red square is the mean and the red bar is the median. It is evident from the results that using the learned linear combination of the metrics outperforms the single-metric based registration.

Figure 4 .

 4 Figure4.12: Results for the RT abdominal dataset for the single-metric registration (sad, mi, ncc, dwt) and the multi-metric registration (mw). The weights for the multi-metric registration are learned using the framework proposed in this work. 'Bladder', 'Sigmoid', and 'Rectum' are the three organs on the dataset. The red square is the mean and the red bar is the median. It is evident from the results that using the learned linear combination of the metrics outperforms the single-metric based registration.

Figure 4 .

 4 Figure4.13: Results for the IBSR dataset for the single-metric registration (sad, mi, ncc, dwt) and the multi-metric registration (mw). The weights for the multimetric registration are learned using the framework proposed in this work. 'CSF', 'Grey Mater', and 'White Mater' are different structures in the brain. The red square is the mean and the red bar is the median. It is evident from the results that using the learned linear combination of the metrics outperforms the singlemetric based deformable registration.

Figure 4 .

 4 Figure 4.14: The top row represents the sample slices from three different volumes of the RT Parotids dataset. The middle row represents the sample slices of the RT Abdominal dataset, and the last row represents the sample slices from the IBSR dataset.

  Obtain the most violated constraint (compute Γi for each sample):Γi = arg min Γ∈L |V |

10 2.2 Semantic Segmentation as a Discrete Labeling Problem . . 12 2.3 Deformable Image Registration as a Discrete Labeling Prob- lem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 MAP Estimation on Graphical Models . . . . . . . . . . . . 18 2

Contents 2.1 Graphical Models in Computer Vision . . . . . . . . . . . . . .4.1 Loopy Belief Propagation . . . . . . . . . . . . . . . . . . . . 19 2.4.2 Lazy Flipper . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.3 FastPD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Why Discrete Methods? . . . . . . . . . . . . . . . . . . . . .

22 10 Chapter 2. Discrete Models in Computer Vision

Graph-based Slice-to-Volume Deformable Registration . . . 51

	3.1.2 Matching Criterion . . . . . . . . . . . . . . . . . . . . . . . .	26
	3.1.3 Transformation Model . . . . . . . . . . . . . . . . . . . . . .	36
	3.1.4 Optimization Method . . . . . . . . . . . . . . . . . . . . . .	40
	3.1.5 Number of Slices . . . . . . . . . . . . . . . . . . . . . . . . .	45
	3.1.6 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . .	48
	3.2 3.2.1 Overparameterized Approach . . . . . . . . . . . . . . . . . .	53
	3.2.2 Decoupled Approach . . . . . . . . . . . . . . . . . . . . . . .	58
	3.2.3 High-Order Approach . . . . . . . . . . . . . . . . . . . . . .	62
	3.3

Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 64

	3.3.1 Inference Methods . . . . . . . . . . . . . . . . . . . . . . . .	64
	3.3.2 Experimental Validation . . . . . . . . . . . . . . . . . . . . .	66
	3.3.3 Comparative Analysis . . . . . . . . . . . . . . . . . . . . . .

Chapter 3. Slice-to-Volume Registration

2. Prior-based Co-Registration and Segmentation 93 Algorithm 2

  3) and S i is the latest 4.Iterative Co-registration Segmentation algorithm 1: procedure ICS(DI : {I1, I2, ..., IN }, DS : {S1, S2, ..., SN })

	2:	Initialize the deformation fields T1, T2, ..., TN as null (identity) deformation fields
	3:	repeat
	4:	repeat
	5:	Randomly select one image from DI as temporal target image I k
	6:

  version of the original input mask S i . S is defined as {S i } i =k and S k is the segmentation prior for the target image.Finally, the regularization term E R from equation 4.4 is redefined as E R following:

				11)
	15:	end for
	16:	Apply label fusion (e.g., Majority Voting) on {S i} i∈{1,...,N } to obtain Ŝk
	17:	end for
	18:	Output: T =	T1, T2, ..., TN and DS =	Ŝ1, Ŝ2, ..., ŜN
	19: end procedure

deformed

Table 4 .

 4 1: Comparison of the average and standard deviation of the Dice coefficient, HD and CMD on the ISBR and RE datasets. As we can observe, when we use accurate priors (Coreg+CNN

  Γi ∈ W i , ∀i w p ≥ 0, ξ i ≥ 0, ∀i. No working set W i can be further updated.Update the parameters: w t ← w 14: until Objective of the problem does not decrease more than ε. Imaging. International Society for Optics and Photonics, 2010. (Cited on pages 40, 79 and 122.) [Zitová 2003] Barbara Zitová and Jan Flusser. Image registration methods: A survey. Image and Vision Computing

	10:	Solve the following optimization problem to obtain w:
		min w,{ξ i }	1 2	||w|| 2 + α||w -w 0 || 2 +	C N i	ξ i ,
		s.t.
	11: until 12: t ← t + 1
	13:

w Ψ(x i , Γi ) ≤ w Ψ(x i , Γi ) -∆(S I i • DΓ , S J i ) + ξ i , ∀

Note that we are assuming that both images have the same dimension d (i.e. we are performing monodimensional deformable registration, e.g. using

2D or 3D images). When we mix different dimensionalities (for example, we register 2D to 3D images), the formulation of the problem must be different. In the chapter 3 of this thesis, we present an extension to this standard formulation, adapted to the case where we aim at registering 2D slices within a 3D volume.

3.4 Conclusions andFurther Ideas . . . . . . . . . . . . . . . . . 78

The source code of the classifier can be downloaded from https://github.com/tsogkas/ brainseg

Acknowledgments

Chapter 4

Context Specific Registration Methods

In this chapter we explore how semantic information can be used to improve the registration process. Existing methods for medical image registration rely on a single metric function explaining the similarity between the source and target images. We argue that incorporating semantic information to guide the registration process, will further improve the accuracy of the results. To this end, we consider two different scenarios. In the first scenario, we are given a set of unlabeled input images and pre-trained classifiers. These classifiers can estimate semantic probability maps indicating how likely every voxel corresponds to a given anatomical structure. This information is incorporated within a co-registration and segmentation algorithm that jointly maps all the input volumes to a common reference space while generating a semantic segmentation mask per image. In the second scenario, we consider the case where hard semantic segmentation masks are available for a training dataset. We use them towards learning context-specific aggregations of standard similarity measures. The learned matching criterion is integrated to the standard graph based image registration formulation, resulting in a multi-metric algorithm endowed with a spatially varying similarity metric function conditioned on the anatomical structures.

Prior-based Co-Registration and Segmentation

In this section we propose a modular and scalable framework for dense co-registration and co-segmentation that incorporates (soft) prior information. The central idea of this segmentation-aware method is to use estimated maps of semantic labels (provided by an arbitrary classifier) as a surrogate for unlabelled data, and combine them with (population) deformable registration to improve both alignment and segmentation. The method seeks to deform all volumes towards an arbitrary reference space where the transformed images exhibit consensus with respect to image similarities (according to some arbitrary metric) and semantic labels. We employ an approach sharing conceptual similarities with α-expansion, combined with efficient linear programming.

Chapter 4. Context Specific Registration Methods

segmentation. Note that the quality of these priors is not as good as the one of the CNN priors. andRF) Figures 4.2-4.6 show the average Dice coefficient, HD and CMD for the co-registration-segmentation method based on RF priors. To further demonstrate the interest of using the iterative co-registration process on top of CNN/RF priors, we also report the results of the simple segmentation methods based on taking the maximum probability given by CNN and RF without co-registration. In this setting, given a CNN/RF prior, the segmentation class of every voxel is simply chosen as the class with the highest probability. These results show that the CNN priors (without any registration process) are more accurate than the weak RF priors.

Thresholded Priors (CNN

Ground Truth Based Oracle (Coreg+GT) So far we have highlighted the positive effect of using priors in our proposed co-registration-segmentation method. As mentioned before, the main idea is to use the probability maps as a surrogate for the ground truth annotations to guide the co-registration process. In order to assess the maximum potential of our approach, we implement an oracle variant that can provide us with an upper-bound for the performance of our method. Hence, we employ the ground truth segmentation masks for all 3D volumes, except for the target image, for which we keep the priors computed by the CNN (this could be considered as an almost-ideal case in which our classifier can generate extremely accurate probability maps). We refer to this method as Coreg+GT in the figures.

Unsurprisingly, the results of Coreg+GT after running Algorithm 2 outperform all other combinations (See Figures 4.2-4.6). Nonetheless, performance of the Coreg+CNN variant is close to Coreg+GT in most cases. This evidence solidifies our original claim, that reliable priors can act as a practical substitute for gold standard annotations for multi-atlas segmentation.

In Figure 4.7 we show some qualitative results, while in 4.2.4 we report the average Dice coefficient, HD and CMD for all structures in the IBSR and RE datasets. The oracle combination Coreg+GT performs best, as expected, with Coreg+CNN a close second, outperforming the other approaches. However, Coreg+CNN has the additional advantage of using learned priors instead of ground truth annotations, which are not always available.

Discussion

In this section, we propose a generic co-registration-segmentation method that can benefit from (soft) segmentation priors and produce results which are close to those generated using ground-truth. Our method is not a multi-atlas segmentation algorithm, since the segmentation results are a by-product of the co-registration where context labels are used to improve the registration performance. We have shown that, by using accurate segmentation priors (like CNN based segmentation maps) we can obtain such results close to the ground truth. More importantly, we also Chapter 5. Conclusions and future directions ing deformable slice-to-volume registration finding the optimal dense deformation field. From this deformation field, we extract rigid and deformable parameters as explained in section 3.2. However, we mentioned that in some specific problems, rigid models are enough since deformations are not considered (see section 3.1.3.1). Unnecessary efforts are invested in exploring an immense parameter space when applying our method to these cases, where only 6-DOF should be considered. Therefore, rigid (only) slice-to-volume registration using discrete models constitutes an open research direction. In this regard, the work of [Zikic 2010b] about linear image registration using MRF and discrete methods could be extended to the case of rigid slice-to-volume registration.

Appendix

A.1 Factor Graphs

In this appendix, we give a brief description of the factor graph data structure, which was used in this thesis to implement the three graphical models proposed in Chapter 2 that aim at solving the challenging slice-to-volume registration problem.

A factor graph G is a bipartite graph that factorizes a given global energy function, expressing which variables are arguments of which local functions [START_REF] Kschischang | [END_REF]. Given a graphical model of any order G = V, E , we can derive a factor graph G = V , F , E . Here, V is the set of variable nodes formed by the nodes of G, F is a the set of all the factors f ∈ F (where every f is associated to one clique G), and the set E ⊂ V × F defines the relation between the nodes and the factors. Every factor f has a function ϕ f : V n → R associated with it, that might correspond to one of the data or regularization terms defined in previous sections. The energy function of a discrete labeling problem like the one presented in 2.4 in the context of factor graphs is then given by:

where x corresponds to a given labeling for the complete graph and l f 1 ...l f n are labels given to the variables in the neighborhood (or scope) of the factor f .

A.2 Convolutional Neural Networks for Brain Segmentation

In this appendix, we give a short description of the Convolutional Neural Network (CNN) based classifier used to learn the segmentation priors applied in section 4.2. This classifier is part of our recent segmentation method that was presented in [START_REF] Shakeri | [END_REF]].

In the past few years, deep learning techniques, and particularly CNNs, have rapidly become the tool of choice for tackling challenging computer vision tasks. CNNs were popularized by Lecun, after delivering state of the art results on handwritten digit recognition [START_REF] Lecun | [END_REF]]. However, they fell out of favor in the following years, mostly due to hardware and training data limitations. Nowadays, the availability of large-scale datasets (e.g. ImageNet), powerful GPUs and appropriate , where C is the number of class labels. We use 1-pixel stride for all convolutional layers and 0.5 activation probability for all dropout layers. In A.1 we report the layer types and parameter values in detail. At test time, our network processes 2D slices from a 3D volume. The output is a three-dimensional array of probability maps (one for each class/brain structure), obtained via a softmax operation. We can readily obtain a brain segmentation at this stage, by simply resizing the output probability array to the input image dimensions using bilinear interpolation and assigning at each pixel the label with the highest probability. As we discuss in section 4.2, these segmentations can be used by the co-registration-segmentation algorithm to guide the deformation process.