Alan Turing

Enzo Ferrante
email: ferrante.enzo@gmail.com

Nikos Paragios

Graph Based Slice-to-Volume Deformable Registration

Keywords: Slice-to-volume registration, graphical models, deformable registration, discrete optimization Slice-to-volume registration, modèles graphiques, recalage déformable, optimisation discrète

Deformable image registration plays a fundamental role in many clinical applications. In this paper we investigate the use of graphical models in the context of a particular type of image registration problem, known as slice-to-volume registration. We introduce a scalable, modular and flexible formulation that 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 pair-wise components) or exact definition of the data terms and the constraints on the plane selection. 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.

Introduction

Slice-to-volume deformable registration is an important problem in the community of medical image computing, which has received considerable attention during the last decade. In general terms, it seeks to determine the slice (corresponding to an arbitrary plane) from a given target volume that corresponds to the deformed version of a source 2D image. This slice is generally specified by a rigid transformation T . The source 2D image is deformed by a deformation field D towards improving the matching consistency between the deformed source image and the target slice.

Slice-to-volume registration is sometimes referred as 2D/3D registration, primarily due to dimension of the images involved in the registration process. Note that this term describes two different problems depending on the technology used to capture the 2D image: it might be a projective (e.g. x-ray) or sliced (e.g. ultrasound (US)) image. In this work we only focus on the latter case. Projective images have to be treated in a different way (basically a pixel in the 2D image does not correspond only to a voxel from the target volume, but to a projection of a set of them in certain perspective) and they are out of the scope of this paper. This is principally due to the fact that conventional image similarity terms cannot be used in the projective case. However, it should be noted that the proposed formulation with an appropriate definition of the matching and regularization cost could also accommodate a solution to this problem. We refer the reader to the comprehensive survey by [START_REF] Markelj | A review of 3d/2d registration methods for image-guided interventions[END_REF] for further information about this topic.

Motivation

A broad number of medical image computing applications benefit from slice-to-volume registration. One can cite, for example, image guided surgeries and therapies [FDB + 03], biopsies [XLF + 14], tracking of particular organs [GAV + 08] and minimally-invasive procedures [LZS + 13, HMG + 09]. In such a context, slice-to-volume registration is a key element for bringing 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 prior to the operation. During the procedure, 2D real time images are generated using different technologies (e.g. fluoroCT, US or interventional MRI slices). These intra-operative images refer to challenging acquisition constraints and inherit lower resolution and quality than the pre-operative ones. Moreover, tissue shift collapse as well as breathing and heart motion during the procedure, causes elastic deformation in the images. Non-rigid image registration is suitable to address this issue. The alignment of intra-operative images with pre-operative volumes augments the information that doctors have access to, and allows them to navigate the volumetric annotation while performing the operation.

Another interesting application is motion correction for image reconstruction. Here, the goal is to correct for misaligned slices when reconstructing a volume of a certain modality. A typical approach to solve it consists of 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 that introduced this idea in the context of functional MRI (fMRI) was presented by [KBB + 99]. More recently, applications of slice-to-volume registration to the same problem in different contexts like cardiac magnetic resonance (CMR) [CPN + 08], fetal images [SFC + 13] and diffusion tensor imaging (DTI) [JXC + 09] have shown promising results.

Although the goals of the motivating problems we have described are different, all of them require to perform (to some extent) slice-to-volume registration. In this work, we focus on the applications where we need to navigate a pre-operative volume using intra-operative images. However, the method we present is modular enough to be adapted to different image modalities and settings, and therefore can be applied to any of these problems.

Previous work

Several methods have been proposed during the recent years to deal with slice-to-volume registration. Some of them deal only with rigid registration, and therefore they cannot manage deformations due to tissue shift, breathing or heart motion. [START_REF] José Estépar | Towards real time 2d to 3d registration for ultrasound-guided endoscopic and laparoscopic procedures[END_REF], for example, proposed a method to register endoscopic and laparoscopic ultrasound images with pre-operative computed tomography volumes that potentially could work in real time. It is based on a new phase correlation technique called LEPART and it handles rigid registration.

[GAV + 08] tracks intra-operative MRI slices of prostate images with a pre-operative MRI volume. This monomodal registration is designed to provide patient tracking information for prostate biopsy performed under MR guidance, but is also constrained to rigid transformations. More recently, [START_REF] Aydin Eresen | Correlating 2d histological slice with 3d mri image volume using smart phone as an interactive tool for muscle study[END_REF] proposed a method that uses smart phone as a navigation tool for initial slice alignment followed by an overlap invariant mutual information-based refinement that estimates the rigid transformation.

Other methods tackle the challenging problem of non-rigid slice-to-volume registration using nonlinear models. Among these, there is a sub-category of approaches that uses several slices instead of a single one, in order to improve the quality of the results. Some examples are [OBH + 11] which uses a variational approach and [XLF + 14] who designed a two-step algorithm where initial rigid registration is followed by B-spline based deformable registration. Using several slices restricts the potential applications to the ones where more than one slice is available from the beginning. It also simplifies the problem by increasing the amount of available information. Our method performs slice-to-volume registration using a single input slice. Consequently, it can be adapted to a broader range of applications where just one slice is available at a time.

An important factor in medical image registration is the type of information that we use to determine correspondences between images, the so-called matching criterion. Methods can be classified according to this aspect in three main categories [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]: geometric, iconic and hybrid methods. In the first case, meaningful anatomical landmarks or features are extracted from the images. After that, the optimal transformation is obtained by mapping the corresponding landmarks. [START_REF] Dalvi | Fast feature based multi slice to volume registration using phase congruency[END_REF] applied this idea to solve slice-to-volume registration combining landmark extraction based on phase congruency and Gabor wavelets with Iterative Closest Point (ICP) optimization. One of the major disadvantages of this type of methods is that results are highly dependent on the quality of the extracted/determined landmarks. In the second case (iconic registration, also known as intensity based methods) [GKT + 08], image intensities are directly used to perform registration. A similarity measure is used to evaluate how similar the images are. Depending on the modalities we are trying to register, simple functions as sum of absolute differences (SAD) or more complex ones as mutual information (MI) might be applied on the intensity domain. [YWK + 10, FDB + 03] applied iconic registration using variations of the standard mutual information as similarity measure in the context of laparoscopic and iMRI-guided interventions, correspondingly. More Inria complex domain specific similarity measures have also been defined to deal with specific problems like [START_REF] Fuerst | Automatic ultrasound-mri registration for neurosurgery using the 2d and 3d lc2 metric[END_REF], where the authors define a metric that can be used to perform US to CT registration. Last but not least, hybrid methods combine both types of information in an effort to get the best of both worlds. Our framework is intensity based and modular with respect to the similarity measure. Therefore, it can be adapted to various image modalities by choosing different functions.

In this paper we extend our previous work presented in [START_REF] Ferrante | Non-rigid 2d-3d medical image registration using markov random fields[END_REF][START_REF] Ferrante | Slice-to-volume deformable registration: Efficient one shot consensus between plane selection and in-plane deformation[END_REF][START_REF] Ferrante | Implicit planar and inplane deformable mapping in medical images through high order graphs[END_REF] through the introduction of a single, mathematically rigorous and theoretically sound framework derived as a discrete labeling problem on a graphical model. Graphical models and discrete optimization are powerful formalisms that have been successfully used during the past years in the field of computer vision [START_REF] Wang | Markov random field modeling, inference & learning in computer vision & image understanding: A survey[END_REF]. In particular, rigid as well as non-rigid image registration have been formulated as a minimal cost graph problem where the nodes of the graph correspond to the deformation grid and the graph connectivity encodes regularization constraints. However, this technique has been applied mainly to mono-dimensional cases (2D-2D or 3D-3D). To the best of our knowledge, the only work that focuses on multidimensional image registration (apart of our previous articles that have been referenced at the beginning of this paragraph) using this type of techniques is [ZGK + 10]. However, it estimates only rigid transformations and works with projective images.

Contribution

This article contributes to enrich the standard graph-based deformable registration theory by extending it to the case of slice-to-volume registration. We present three different models to solve this challenging problem which vary in terms of graph topology, label space definition and energy construction. 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. We analyze the strong and weak points of every model and we perform comparative experiments. Validation is done using a monomodal MRI cardiac dataset and a multimodal brain dataset [MDMP + 12] including different inference methods.

Graph-based slice-to-volume deformable registration

An enormous variety of tasks in computer vision and medical image analysis can be expressed as discrete labeling problems. Low, mid and high-level vision tasks can be addressed within this framework. To this end, a visual perception task is addressed by specifying a task-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 [START_REF] Paragios | Discrete Visual Perception[END_REF].

In the context of graph-based discrete labeling problems, the model is composed by a graph G = V, E where vertices in V correspond to the variables while E is a neighborhood system (pair-wise & higher order cliques) that encodes the relationships among these variables. We also consider a discretized version of the search space that is represented by a discrete set of labels l ∈ L. The aim is to assign to every variable v ∈ V a label l v ∈ L.

Each time we choose to assign a label, say, l v1 to a variable v 1 , we are forced to pay a price according to the so-called energy function. This objective function is domain-specific and associates the observations to the model. It is formulated as the sum of singleton terms g v (l v ) (which depend only on one label l v ), pairwise terms f v1v2 (l v1 , l v2 ) (which depend on two variables l v1 , l v2 ) and high-order terms f v1...vn (l v i 1 , . . . , l

v |C i | i
) (which are associated to high-order cliques C i that depend on more than two variables). Our goal is then to choose a labeling which will allow us to recover the solution corresponding to the minimal value of the objective function. In other words, we want to choose a labeling that minimizes the sum of all the energy potentials, or equivalently the energy P(g, f ). This amounts to solving the following optimization problem:

argmin lp P(g, f ) = v∈V g p (l v ) + (v1,v2)∈E f v1v2 (l v1 , l v2 ) + Ci∈E f v1...vn (l v i 1 , . . . , l v |C i | i ), (1)
Performing parameter inference on this graphical model, could be an effective solution to a big variety of problems in computational medicine. Note that we make a distinction between singleton, pairwise and high-order terms, depending on the number of variables jointly interacting. It should be noted that most part of the graph-based vision models have explored mainly pairwise constraints (pairwise Conditional and Markov Random Field (CRF/MRF) models), because in these cases exact or approximate efficient inference of Maximum a Posteriori (MAP) solutions can be done. However, during the last few years, more and more high-order models and inference algorithms have been developed which offer higher modeling power and can lead to more accurate solutions of the problems [START_REF] Kohli | Higher-order models in computer vision[END_REF][START_REF] Komodakis | Mrf energy minimization and beyond via dual decomposition[END_REF]. Given such a general setting, let us now try to explore the expressive power of such models in the context of slice-to-volume deformable registration.

The task of slice-to-volume deformable registration can be expressed mathematically as follows. Given a 2D source image I and a 3D target volume J, we seek the 2D-2D in-plane local deformation field TD and the plane π[J] (i.e. a bi-dimensional slice from the volume J)

which in the most general case minimize the following objective function:

TD , π = argmin

T D ,π M(I • T D (x), π[J](x)) + R(T D , π), (2)
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. 2) 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) [RSH + 99] where each control point has local influence on the deformation and is used to approximate the dense deformation field TD . Having said that, alternative control point interpolation models can be used as well. Figure 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 [START_REF] Ferrante | Non-rigid 2d-3d medical image registration using markov random fields[END_REF]. 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 | Efficient mrf deformation model for non-rigid image matching[END_REF], we present a decoupled model where linear and deformable parameters are separated into two interconnected subgraphs which refer to lower dimensional label spaces [START_REF] Ferrante | Slice-to-volume deformable registration: Efficient one shot consensus between plane selection and in-plane deformation[END_REF]. 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, 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.

d 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 i = (φ, θ). (d) Displacement of the control points p i and p j when the corresponding labels l i = (d i , N i ) and l j = (d j , N j ) are applied. The planes π i and π j are those that contain the control points p i + d i , p j + d j and whose normals are N i , N j respectively.

Overparameterized approach

Let us consider an undirected pair-wise 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 2.a). The vertices v i ∈ V are moved by assigning them different labels u i ∈ L (where L corresponds to the label space) until an optimal position is found.

In order to deform the graph, we need to define a label space able to describe the inplane deformations and the plane selection variables. To this end, we consider a label space L that consists of 5-tuples l = (d x , d y , d z , φ, θ), where the first three parameters (d x , d y , d z ) define a displacement vector d i in the cartesian coordinate system (see Figure 2.b), and the angles (φ, θ) define a vector N i on a unit sphere, expressed using spherical coordinates (see Figure 2.c). Let us say we have a control point p i = (p xi , p yi , p zi ) and we assign the label l i = (d xi , d yi , d zi , φ i , θ i ) to this point. 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) 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 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 2.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) Inria Figure 3: Data term formulation for the overparameterized approach. The points x ∈ Ω 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. associated with each graph vertex and regularization terms F = {f ij (•, •)} (or pairwise potentials) associated with the edges. As we described in section 2, 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 | Deformable medical image registration: setting the state of the art with discrete methods[END_REF]. The minimization energy problem for the overparameterized formulaton is thus defined as:

P O (G, F ) = min i∈V g i (l i ) + γ (i,j)∈E f ij (l i , l j ), (4)
where l i , l j ∈ L are the labels assigned to the vertices v i , v j ∈ V respectively.

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 bi-dimensional 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 Ω 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):

g i (l i ) = Ωi δ(I(x), π i [J](x))dx.
(5) 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 ) = Ωi | I(x) -π i [J](x) | dx. (6)
In multimodal scenarios, where different modalities are compared (e.g. CT with Ultrasound 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 ) = - Ωi log p(I(x), π i [J](x)) p(I(x))p(π i [J](x)))
dx.

(7)

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 ). (8)
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 i , d j ) = || (p i + d i ) -(p j + d j ) || || (p o,i ) -(p o,j ) || . (9)
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 i , d j )) 2 + (1 -ψ i,j (d i , d j ) -1 ) 2 , (10)
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 4.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 . ( 11
)
F 2 is defined using this distance (Equation 11) and corresponds to the average of D πj (p i + d i ) and D πi (p j + d j ):

F 2i,j (l i , l j ) = 1 2 (D πj (p i + d i ) + D πi (p j + d j )). (12)
Inria 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(θ).

(

) 13
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 | Efficient mrf deformation model for non-rigid image matching[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 pair-wise 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 2.1 (see Figure 2.a) to improve the satisfaction of the desired geometrical constraints.

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 l I = (d x , d y ). 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 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 l P = (φ, θ, λ), where angles φ and θ define a vector N over a unit sphere, expressed through its spherical coordinates (see Figure 2.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 + 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 5 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, smoothness of the deformation and unique plane selection:

P D (I, P, D) = min α (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 ), (14)
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. Inria Figure 5: Data term formulation for the decoupled approach. It is similar to the formulation shown in Figure 3, but it combines labels from different label spaces. The points x ∈ Ω 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 in-plane deformation label l I = (d x , d y ) is applied to the corresponding imaginary grid point p k . Then, label l P = (N, λ) is used: the point is translated in the direction given by vector N as indicates scalar λ. Finally, the patch Ω k is sampled from plane π k with normal N , centered at the displaced point p k (in orange).

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 Ω 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 5):

e D i,j (l I i , l P j ) = Ω k δ(I(x), π k [J](x))dx. (15)
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 (l I i , l I j )) 2 + (1 -ψ i,j (l I i , l I j ) -1 ) 2 , ( 16
)
where ψ i,j is the distance defined in Equation 9.

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 12) point-to-plane distance:

e P i,j (l P i , l P j ) = 1 2 (D πj (p i ) + D πi (p j )). (17)
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 [START_REF] Glocker | Approximated curvature penalty in non-rigid registration using pairwise mrfs[END_REF], 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 twodimensional quasi-planar grid and they are displaced using 3D vectors 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, third-order 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 [START_REF] Kwon | Nonrigid image registration using dynamic higher-order mrf model[END_REF], 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 2.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 . The minimization energy problem becomes:

P H (D, R) = min (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 ), ( 18
)
where γ is a scaling factor and l i is a label associated with a displacement vector (d x , d y , d z ) and assigned to the node i.

The data term is defined over a disjoint set of triangular cliques, covering the entire 2D domain, as shown in Figure 6.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 Inria points of the clique. For a given similarity measure δ, the data term associated with the clique c is thus defined as:

D ijk (l i , l j , l k ) = Ω (l i ,l j ,l k ) δ(I(x), π d [J](x))dx, (19)
where x ∈ Ω (li,lj ,l k ) , and Ω (li,lj ,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 corresponding labels l i , l j , l k to the vertices.

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 6.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 [START_REF] Kwon | Nonrigid image registration using dynamic higher-order mrf model[END_REF]. Given three contiguous control points (p i , p j , p k ) and their corresponding displacement labels (l i , l j , l k ), it can be approximated as follows:

|| (p i + l i ) + (p k + l k ) -2 • (p j + 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 + l i ) + (p k + l k ) -2 • (p j + l j ) || d 2 2 , (20)
In-plane deformation smoothness is reinforced in the same manner as the previous mod- els -through a symmetric distance preserving approach. For the sake of clarity, we redefine Equation 10 as Ψ ij (l i , l j ) = (1 -ψ i,j (l i , l j )) 2 + (1 -ψ i,j (l i , 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 (l i , l j ) + Ψ jk (l j , l k ) 2 . ( 21
)
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 ), (22)
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 Inria 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 | Exact inference in multi-label crfs with higher order cliques[END_REF]).

In [START_REF] Ferrante | Non-rigid 2d-3d medical image registration using markov random fields[END_REF], the overparameterized approach was optimized using the FastPD algorithm [START_REF] Komodakis | Fast, approximately optimal solutions for single and dynamic mrfs[END_REF] while for the decoupled [START_REF] Ferrante | Slice-to-volume deformable registration: Efficient one shot consensus between plane selection and in-plane deformation[END_REF] and the higher order models [START_REF] Ferrante | Implicit planar and inplane deformable mapping in medical images through high order graphs[END_REF], 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 adapted 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 chose 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 [AKB + 12], a move-making algorithm which is a generalization of the classical Iterated Conditional Modes (ICM) [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] and has provided good approximations for several non-submodular models in different benchmarks . Both are approximate inference methods that can accommodate arbitrary energy functions, graph topologies and label spaces, and allow us to show how the three proposed approaches perform under different optimization strategies.

We have adopted the OpenGM2 library [KAH + 13] which implements both inference methods, and makes it possible to perform fair comparisons. It requires construction of a factor graph for every scheme (see Figure 7).

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] Frank R Kschischang | Factor graphs and the sum-product algorithm[END_REF]. Given a graphical model of any order G = V, E (like the ones described in this work), 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 our discrete labeling problem in the context of factor graphs is then given by:

E(x) = f ∈F ϕ f (l f 1 , ..., l f n ), (23)
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 . Figure 7 shows a comparison between the three models and the derivation of the corresponding factor graph in each case.

Experimental validation

We computed 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 one is a monomodal MRI heart dataset while the second one corresponds to multimodal US-MRI brain images. 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 inference methods, we used 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 used a maximum depth of 2 (for details about LF, we refer the reader to [KAH + 13]).

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

Monomodal dataset experiment

The monomodal dataset was 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 was generated from a temporal series of 3D heart MRI volumes M i as shown in Figure 8. For a given sequence in the heart dataset, every 2D slice I i was extracted from the corresponding volume M i at a position which is calculated as follows. Starting from a random initial translation T 0 = (T x0 , T y0 , T z0 ) and rotation R 0 = (R x0 , R y0 , R z0 ), 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 used σ 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 generated 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 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 Inria 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 used SAD as similarity measure given that 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 was used to initialize the registration of slice I i+1 .

Figure 10 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 approaches and optimization methods. The decoupled method outperforms the other two 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 used the segmentations being associated with the dataset. We computed 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 evaluated accuracy computing DICE coefficient, Hausdorff distance and Contour Mean Distance (CMD). We also provided 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 11 show the mean and standard deviation of the indicators at the three stages, for the three approaches and the two inference methods. It can be seen that results improve at each stage, achieving DICE coefficient of around 0.9 after deformation. Hausdorff distance and CMD decreased 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 10 for the rigid parameters). Figure 12 complements these results by showing DICE values disaggregated per sequence, while Figure 9 shows some qualitative results before, after rigid and after deformable registration. Finally, in terms of running time, Figure 13 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. Inria Figure 10: 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. Independently of the 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.

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 [MDMP + 12]. 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 Figure 11: 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 three methods (overparameterized, high-order and decoupled) using both inference strategies (LBP and LF).

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 Mutual Information as similarity measure instead of SAD.

Figure 11 summarizes the average DICE, specificity, sensibility, Hausdorff distance and Contour Mean Distance coefficients for all the series, while Figure 13 reports the running times. Figure 12 complements these results by showing DICE values disaggregated per sequence. While 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 Inria Figure 12: Final DICE (after deformation) comparison for every sequence (10 sequences in the heart dataset and 6 sequences in the brain dataset). Results are shown for all three approaches and two inference methods. In case of the heart dataset, decoupled method outperforms the other two in most part of the sequences. In the brain dataset, the highorder approach shows best performance in most cases. This is coherent with the aggregated results shown in Figure 11. 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 an extremely difficult task.

In all brain experiments we used initial grid size of 8mm, minimum patch size of 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 were chosen similarly as in the heart Figure 13: 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. experiments.

Comparative analysis

In this section, we aim at comparing different aspects of the three approaches we have presented 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 don't 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 14 shows a comparative analysis between the three approaches, using the two proposed inference methods, in terms of optimization time and final energy. Table 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 Inria Figure 14: 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. of a given label space by making smart decisions in terms of energy definition and graph topology.

But 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 that while the number of variables and edges augment linearly (it goes from N to 2N in case of variables, and from 2N -2 It results in better performance for the decoupled method as can be observed in Figure 13. 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 13, we achieve good experimental computation time because of the smaller label space.

Conclusion

We derived three new models from the standard graph-based deformable registration theory for slice-to-volume registration. We have shown promising results in a monomodal and a multimodal case, using different inference methods. 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 have presented 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 has been 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 [SOG + 10] could also be an interesting additional component of the registration criterion. Last but not least, domain/problem specific parameter learning [BGK + 13, KXP15] towards improving the proposed models could have a positive influence on the obtained results.

Acknowledgements

This research was partially supported by European Research Council Starting Grant Diocles (ERC-STG-259112). We thank Mihir Sahasrabudhe for proof-reading the paper, and Vivien Fecamp and Jorg Kappes for helpful discussions. Inria

Figure 1 :

 1 Figure 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 recon- structed 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 2 :

 2 Figure 2: (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

Figure 4 :

 4 Figure 4: (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 .

Figure 6 :

 6 Figure 6: 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.

Figure 7 :

 7 Figure 7: Factor graph derivation for the overparameterized (A), decoupled (B) and highorder (C) approaches. It shows the equivalence between cliques c i (unary, pairwise and high-order, depending on the model) and the corresponding factors f i . In red, we observe the cliques and factors associated with data terms, while in green and orange we represent those associated with regularization terms.

Figure 8 :

 8 Figure 8: 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).

Figure 9 :

 9 Figure 9: 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.

√

  N to 5N -4 √ N in case of pairwise edges), the number of labels decreases Parameter Overparameterized Decoupled High-order Label space |L D × L P | max(|L D |, |L P |) |L D |

Table 1 :

 1 Comparison among the three methods in terms of label space and graph topology. quadratically (from |L D × L P | to max(|L D |, |L P |)).

RR n°8803

Ferrante & Paragios

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Armando V D B Assis
email: armando.assis@outlook.com.br

Elementary Considerations on Prime Numbers and on the Riemann Hypothesis *

Keywords: Riemann hypothesis, prime number, number theory, disproof COMBINATORICS

The Riemann Hypothesis states the non trivial zeros of a mathematical function, the Riemann zeta function, are, all of them, points pertaining to a vertical line in the Argand-Gauss plane. Technically, this implies there exists an intrinsical order, a balance, within the structure of natural numbers, related to the building blocks, these the prime numbers. We prove the natural numbers are free to exist, i.e., the Riemann Hypothesis is false, which is the main purpose of this paper: to provide an elementary disproof of the Riemann Hypothesis. This version contains an appendix that should suffice for explanation on the cardinality of N, ℵ0, which is important here, emphasizing denumerable means countable [countably infinity], from which the scope of combinatorics is justified.

The Liouville function λ (n) [START_REF] Apostol | Introduction to analytic number theory[END_REF] depending on the variable n ∈ N = {1, 2, 3, 4, • • • } (in this paper, 0 / ∈ N) is given by:

λ (n) = (-1) ω(n) , (1)
where:

ω (n) = Number of prime factors of n, (2)
being these prime factors not necessarily distinct, counted with multiplicity. Hence, the image of λ (n) is {1, -1}. E.g.: λ (1) = 1, since 1 has not got any prime factors (ω (1) = 0); λ (2) = -1, since 2 has got just one prime factor (itself, since 2 is a prime number), from which ω (2) = 1; λ (20) = -1, since, counted with multiplicity, 20 has got 3 prime factors (20 = 2 × 2 × 5), from which ω (20) = 3. Of course, all the numbers n ∈ N, except the number 1, will have prime factors, being such quantity of prime factors, exactly ω (n), either odd or even. A prime number has got itself as its prime factor, hence λ (p) = -1, i.e., being p ∈ N prime: ω (p) = 1. By virtue of these considerations, one may be interested in a generic factorization with an even number of prime factors and also be interested in a factorization with an odd number of prime numbers. The presence of 1 as a factor in a given factorization does not matter, as well as for any quantity of the number 1 as a factor, except for the unique case of the factorization of the number 1, this latter being trivial and unique, once a consideration of 1 factor(s) in a given natural number n = 1 contributes nothing to ω (n), viz.:

ω (n) = ω (1 × 1 × • • • × 1 × • • • × 1 × n) , (3)
by the very elementary fact that we just count the quantity of prime factors of an n. Also, once a factorization is unique, except for the order of the prime factors of a given natural number, a given even (or odd) factorization with 2k (or 2k -1), k ∈ N, factors turns out to be unique under a combinatorial (combination) consideration with repetition with 2k (or 2k -1) elements. As a matter of fact and clarification, here, we should start to put these assertions under a more mathematically generalized sound. Let P N = {p 1 , p 2 , • • • , p N } be the set of all the first N prime numbers, so that P = P ∞ is the set of the prime numbers. One may consider 2k, k ∈ N, slots filled with any p l element, l ∈ {1, 2, • • • , N }, from the P N set, e.g.:

p 3 × p 5 × p 5 × p 2 × p 3 × p 3 × p 8 × • • • × p 2 (2k)-slots (4)
Of course, any permutation of (4) generates the very same configuration, generates the very same number due to the multiplication (× between slots). One should infer that repetition is allowed, viz., the elements in (4) do not need to be different. However, for an instantaneously fixed order of (4), a change in a given slot element (choosing a different one from P N ), with the elements within the remaining slots not changed, would lead to a new configuration for [START_REF] Borwein | The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike[END_REF]. A permutation of this latter new configuration does not change it. A given configuration represents a unique natural number with an even number of prime factors, since, in spite of permutation, its factorization (configuration) is unique. Similarly, one may consider 2k -1, k ∈ N, slots filled with any p l element, l ∈ {1, 2, • • • , N }, from the P N set, e.g.:

p 3 × p 5 × p 5 × p 3 × p 3 × p 8 × • • • × p 2 (2k-1)-slots (5)
Of course, any permutation of (5) generates the very same configuration, generates the very same number due to the multiplication (× between slots). One should infer that repetition is allowed, viz., the elements in [START_REF] Landau | Neuer Beweis der Gleichung ∞ k=1 µ (k) /k =[END_REF] do not need to be different. However, for an instantaneously fixed order of ( 5), a change in a given slot element (choosing a different one from P N ), with the elements within the remaining slots not changed, would lead to a new configuration for [START_REF] Landau | Neuer Beweis der Gleichung ∞ k=1 µ (k) /k =[END_REF]. A permutation of this latter new configuration does not change it. A given configuration represents a unique natural number with an odd number of prime factors, since, in spite of permutation, its factorization (configuration) is unique.

To consider the totality of numbers in N, one needs to impose N → ∞, N ∈ N, hence P N → P ∞ = P, and to consider the totality of 2k slots and to consider the totality of 2k-1 slots, these latter totalities by completely considering k (viz.: ∀k ∈ N).

Since a given configuration does not change with a permutation of its elements, we are dealing with a problem of combination [START_REF] Yaglom | Challenging Mathematical Problems with Elementary Solutions[END_REF] (the order does not matter). However, elements may be repeated, hence the combinatorics involved is neither C N,2k nor C N,2k-1 . To uniquely represent a configuration, we may define a convention. In fact, considering an hypothetical factorization with q ∈ N factors, e.g., as represented below:

p N × p 5 × p 1 × p 2 × p 3 × p 3 × p 8 × • • • × p 2 (q)-slots , (6)
it may be rearranged respecting the order of the indexes:

p 1 × p 2 × p 2 × p 3 × p 3 × p 5 × p 8 × • • • × p N (q)-slots , (7)
which is the very same factorization preserving its same elements, now written in increasing order of indexes, which, to uniquely be represented, may have each index increased by an amount exactly equal to the quantity of previous elements (which is a unique characteristic per element), also changing the label to x, viz.:

x 1 × x 3 × x 4 × x 6 × x 7 × x 10 × x 14 × • • • × x N +q-1 (q)-slots (8)
With this convention, each q-combination turns out to have got all of its elements with different indexes, with the maximum index occurring when a q-combination has got its last slot accupied by p N , from which we turn out to have a problem of simple q-combination of N + q -1 elements. Hence, the quantity Q N q of unique factorizations having q prime, with multiplicity allowed, factors taken from P N is:

Q N q = C N +q-1, q = (N + q -1) ! (N -1) ! q ! , (9)
and the quantity Q N q of unique factorizations having q prime factors turns out to be:

Q N q = lim N →∞ C N +q-1, q = lim N →∞ (N + q -1) ! (N -1) ! q ! = |N q | , (10)
the totality of distinct numbers belonging to N and having got q prime, with multiplicity allowed, factors. The totality of numbers belonging to N turns out to be:

Q N = 1 + ∞ q=1 Q N q = 1 + ∞ q=1 lim N →∞ C N +q-1, q = 1 + ∞ q=1 lim N →∞ (N + q -1) ! (N -1) ! q ! (11) = 1 + ∞ q=1 |N q | = |N| = ℵ 0 .
We will be interested in the infinite sum:

∞ n=1 λ (n) = λ (1) + ∞ n=2 λ (n) = 1 + ∞ n=2 λ (n) = 1 + ∞ k=1 Q N 2k - ∞ k=1 Q N 2k-1 , (12)
where:

Q N 2k = lim N →∞ C N +2k-1, 2k = lim N →∞ (N + 2k -1) ! (N -1) ! (2k) ! , (13)
and:

Q N 2k-1 = lim N →∞ C N +(2k-1)-1, 2k-1 = lim N →∞ (N + 2k -2) ! (N -1) ! (2k -1) ! , (14)
k ∈ N. Considering the partial difference:

D N k = Q N 2k -Q N 2k-1 , (15)
which, by virtue of the Eq. ( 9), leads to:

D N k = C N +2k-1, 2k -C N +(2k-1)-1, 2k-1 = (N + 2k -1) ! (N -1) ! (2k) ! - (N + 2k -2) ! (N -1) ! (2k -1) ! = (N + 2k -1) ! -2k(N + 2k -2) ! (N -1) ! (2k) ! = (N + 2k -1) (N + 2k -2) ! -2k(N + 2k -2) ! (N -1) ! (2k) ! = (N + 2k -1 -2k) (N + 2k -2) ! (N -1) ! (2k) ! = (N -1) (N + 2k -2) ! (N -1) ! (2k) ! = (N -1) (N + 2k -2) ! (N -1) (N -2) ! (2k) ! = (N + 2k -2) ! (N -2) ! (2k) ! = [(N -2) + 2k] ! (N -2) ! (2k) ! = 2k factors [(N -2) + 2k] [(N -2) + 2k -1] × • • • × [(N -2) + 1] (N -2) ! (N -2) ! (2k) ! = 2k factors [(N -2) + 2k] [(N -2) + 2k -1] × • • • × [(N -2) + 1] (2k) ! = 1 (2k) ! 2k j=1 [(N -2) + j] > d N k , (16)
where:

d N k = (N -2) 2k (2k) ! . (17)
Hence, from the Eqs. (15), ( 16) and (17):

∞ k=1 D N k = ∞ k=1 Q N 2k - ∞ k=1 Q N 2k-1 = ∞ k=1 1 (2k) ! 2k j=1 [(N -2) + j] > ∞ k=1 d N k = ∞ k=1 (N -2) 2k (2k) ! = cosh (N -2) -1. (18)
Now, consider the set N -{1}, hence N -{1} exhausts any prime factorization. This is so, since a given prime factorization represents one and only one element in N -{1} and a given element of N -{1} has got one and only one prime factorization. Hence, all the factorizations exhaust N -N -{1} is the entire set of possible prime factorizations and the entire set of possible prime factorizations is N -{1}. The arithmetic progression:

a n = 1 + (n -1)1 = n ⇔ a n = n, ( 19
)
with a n ∈ N -{1}, shows that the number of possible factorizations grows as n grows. Hence, since this is es-sentially the set N -{1}, the quantity of possible factorizations exhaustively grows [a given finite set of factorizations will have a factorization, say f , representing a greatest factorized number in this set and being oneto-one to an a n = f for this set, with the remaining factorizations (< a n )] as n → ∞. Hence, as discussed in the appendix [START_REF] See | [END_REF]:

n → 1 + ∞ k=1 lim N →∞ Q N 2k + ∞ k=1 lim N →∞ Q N 2k-1 , (20) viz.
:

n → 1 + ∞ k=1 Q N 2k + Q N 2k-1 = ∞. ( 21
)
By virtue of Eq. ( 21):

n 1 2 + → 1 + ∞ k=1 Q N 2k + Q N 2k-1 1 2 + = ∞. ( 22
)
Considering the partial sum:

s N k = Q N 2k + Q N 2k-1 , (23)
which, by virtue of the Eq. ( 9), leads to:

s N k = C N +2k-1, 2k + C N +(2k-1)-1, 2k-1 = (N + 2k -1) ! (N -1) ! (2k) ! + (N + 2k -2) ! (N -1) ! (2k -1) ! = (N + 2k -1) ! + 2k(N + 2k -2) ! (N -1) ! (2k) ! = (N + 2k -1) (N + 2k -2) ! + 2k(N + 2k -2) ! (N -1) ! (2k) ! = (N + 2k -1 + 2k) (N + 2k -2) ! (N -1) ! (2k) ! = (N + 4k -1) (N + 2k -2) ! (N -1) ! (2k) ! = (N + 4k -1) [(N -1) + (2k -1)] ! (N -1) ! (2k) ! = (N + 4k -1) 2k-1 factors [(N -1) + (2k -1)] {(N -1) + [(2k -1) -1]} × • • • × [(N -1) + 1] (N -1) ! (N -1) ! (2k) ! = (N + 4k -1) 2k-1 factors [(N -1) + (2k -1)] {(N -1) + [(2k -1) -1]} × • • • × [(N -1) + 1] (2k) ! = N + 4k -1 (2k) ! 2k-1 j=1 [(N -1) + j] < N + 4k -1 + 1 (2k) ! [(N -1) + 2k -1] 2k-1 < N + 4k (2k) ! [(N -1) + 2k -1 + 2k] 2k-1 < N + 4k (2k) ! (N + 4k -2 + 2) 2k-1 ∴ s N k < S N k , (24)
where:

S N k = N + 4k (2k) ! (N + 4k) 2k-1 = (N + 4k) 2k (2k) ! , (25)
remembering: k ∈ N. Furthermore:

S N k < Ξ N k , (26)
where:

Ξ N k = (N + N/α) 2k (2k) ! = 1 + 1 α 2k N 2k (2k) ! , (27)
∀ fixed α ∈ R * + . The Eq. ( 26) follows from the proposition:

4αk < N, ∀ α ∈ R * + . (28)
Since N is to exhaustively cover the set P ∞ = P, this latter being the entire set of prime numbers, viz., in:

P N = {p 1 , p 2 , p 3 , • • • , p N } = {2, 3, 5, • • • , p N } , (29)
N has been taken, defined, to exhaustively cover on demand the set of prime numbers, as we have defined from the beginning of this paper. Now, to conversely suppose the condition stated by the Eq. (28), i.e.:

∃ α ∈ R * + | 4αk ≥ N, (30)
one turns out to be, by hypothesis, considering an implied superior quota for the existence of prime numbers.

Suppose the Eq. ( 30) is correct. Hence, there exists, by virtue of the Eq. ( 30), an α 0 ∈ R * + such that the successive values of N to exhaustively cover the set of prime numbers never exceed 4α 0 k. Putting such 4α 0 k, with the obeyer number of the Eq. ( 30), a fortiori:

4α 0 k = [4α 0 k] + {4α 0 k} , (31)
with [4α 0 k] and {4α 0 k} being, respectively, the integer and the fractionary parts of 4α 0 k, one is led to a superior quota:

N = [4α 0 k] , (32)
since N ∈ N. By virtue of this superior quota for the quantity of prime numbers, there exist only finitely many primes

p 1 < p 2 < • • • < p N . Hence, let η = N f =1 p f > 2
, and consider η -1 ∈ N. Since η -1 is a product of primes, it turns out to have a common prime divisor p f with η, implying this common prime divisor p f divides η -(η -1) = 1: an absurd! Henceforth, Eq. ( 30) is an absurd, the proposition given by the Eq. ( 28) is correct ∀ fixed k ∈ N and, by virtue of the Eqs. ( 23), ( 24), ( 26) and ( 27), one turns out to be led to:

s N k = Q N 2k + Q N 2k-1 < Ξ N k = 1 (2k) ! 1 + 1 α N 2k .
(33) Back to the interest carried from the Eqs. (20), ( 21) and ( 22), now, we consider the consequence to the sum:

∞ k=1 s N k = ∞ k=1 Q N 2k + ∞ k=1 Q N 2k-1 , (34)
i.e.:

∞ k=1 s N k = ∞ k=1 Q N 2k + ∞ k=1 Q N 2k-1 < ∞ k=1 Ξ N k = ∞ k=1 1 (2k) ! 1 + 1 α N 2k = cosh 1 + 1 α N -1 (35)
From the reasonings we have carried throughout the march that has been led from that Eq. ( 12), accom-plished from its previous combinatorics, to the Eq. ( 35), one reaches, a fortiori, the following implied consequence:

lim n→∞ n l=1 λ (l) n 1 2 + = lim n→∞ 1 + n l=2 λ (l) n 1 2 + > lim N →∞ 1 + [cosh (N -2) -1] 1 + cosh 1 + 1 α N -1 1 2 + = lim N →∞ cosh (N -2) cosh 1 + 1 α N 1 2 + ∴ lim n→∞ n l=1 λ (l) n 1 2 + > 2 √ 2 lim N →∞ e (N -2) + e (2-N ) e (1+1/α)N + e -(1+1/α)N 1 2 + = 2 √ 2 lim
N →∞ e (N -2) 1 + e -2(N -2)

e (1+1/α)N 1 + e -2(1+1/α)N e (1+1/α)N 1 + e -2(1+1/α)N 1 + e -2(1+1/α)N

Since the Riemann Hypothesis [START_REF] Riemann | Ueber die Anzahl der Primzahlen unter einer gegebenen Grösse[END_REF] is true if and only if [START_REF] Borwein | The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike[END_REF][START_REF] Landau | Neuer Beweis der Gleichung ∞ k=1 µ (k) /k =[END_REF]:

lim n→∞ n l=1 λ (l) n 1 2 + = 0, ∀ fixed > 0, (37)
it follows that the convergence of the right-hand side of the Eq. (36), ∀ fixed > 0, is a necessary condition for the validity of the Riemann Hypothesis. It follows, then, that:

• If there exists some > 0 such that the right-hand side of the Eq. (36) diverges, then, the Riemann Hypothesis turns out to be false.

Supposing the Riemann Hypothesis is true, the condition:

1 -( + 1/2) (1 + 1/α) < 0 ( 38
)
holds ∀ fixed > 0. Hence:

( + 1/2) (1 + 1/α) > 1, (39)
∀ fixed > 0. But this latter condition is an absurd, since, for = 1/8 and α = 2, e.g.:

( + 1/2) (1 + 1/α)| (α, )=(2,1/8) = 15 16 < 1, (40)
implying the Riemann hypothesis turns out to be false.

1 2 +

 2 e N [1-( +1/2)(1+1/α)] .

ACKNOWLEDGMENTS

A.V.D.B.A is grateful to the Riemann Hypothesis. A.V.D.B.A is grateful to the Giants within the mathematical field, Human Beings and primes. A.V.D.B.A is deeply grateful to his Family.

APPENDIX

The Eq. (20) must be correct, otherwise the Fundamental Theorem of Arithmetic turns out to be contradicted, as one grasps. This theorem asserts the factorization is unique, except for the order of the prime factors being multiplied to generate a natural number. Hence, the ∞, the upper sum limit, is a symbol for the natural numerical limit of ℵ 0 . From this, follows the bijection between the natural numbers and their [unique] factorizations, with the degeneracy of a given factorization removed from the convention we have adopted in this paper, completely exhausting the members of both the sets, N and [their unique] factorizations, one-to-one. In fact, one can define the following function with domain N:

which is obviously invertible, with image N, having got, hence, a one-to-one correspondence, a bijection, where f -1 (n) is the inverse function of f (n). The cardinality of N (remember our definition of N, at the beginning of this paper, the positive integers) is ℵ 0 . By Eq. ( 41), it is crystalline the set of factorizations is the very same N, hence with the same cardinality, ℵ 0 . The sequence (

) tends to the infinity that is symbolized by the cardinal ℵ 0 . Hence, the cardinality of the factorizations must be ℵ 0 , otherwise there would exist some number without factorization and vice-versa, contradicting the Fundamental Theorem of Arithmetic.

Since N is denumerable, hence countable, the set of fac-torizations is [denumerable, hence countable] too. Counting, n must tend to as given by the Eq. ( 20), since n tends to the infinity that is symbolized by the cardinal ℵ 0 , asseverating, as well as the infinity quantity of all the factorizations, Eq. (20).


How to merge three different methods for information filtering ?

Introduction

Online reputation is a key information for public figures and companies in order to react to the public opinions and to anticipate them. Indeed, knowing what make their reputation good or bad allows them to make informed decision. For instance a company may make additional efforts on its call centers if it notices that its consumers are unsatisfied.

Monitoring online reputation of entities requires to be able to retrieve all opinions or reviews about them. Automatic approaches have then to deal with the noise generated by the recall-oriented retrieving techniques used. This noise is mainly the result of entity names ambiguity (e.g. jaguar which may refer to an animal or a car manufacturer). A classification step is required to filter out sources which do not actually mention the monitored entity. Topic detection is necessary to identify which matter is discussed in the source and finally the polarity of it has to be estimated (is the opinion positive, neutral or negative?). Each of these three issues is an open problem. Moreover, systems have to be able to process large amounts of incoming new documents in a short time to provide fresh feedbacks. Sources commonly used are news web sites, blogs, forums or more recently social networks such as Twitter.

We propose three approaches to filter tweets on whether or not they refer to a given entity. These approaches rely on tweets content and meta-data associated to them (timestamp, user names, . . . ) as well as on the information contained about the entity in a knowledge base. We also investigate if combining systems outputs with merging algorithms can improve the overall performances and if different strategies may be applied to promote a measure or another.

Each proposition is evaluated within the framework provided by the RepLab 2013 evaluation campaign and they all show competitive performances.

The remainder of this paper is organized as follows. Section 2 presents related works. Section 3 describes the proposed systems. Section 4 gives details about merging algorithms we used. Experiments are described in Section 5. In Section 6 we discuss the results before concluding in section 7.

Related Work

A decade ago, a TREC task called "Filtering" [START_REF] Robertson | The trec 2002 filtering track report[END_REF] had the following definition: finding documents relevant to a query in a stream of data. Effective approaches were in-spired by information retrieval techniques to score documents (Okapi (Robertson and al, 2002), Rocchio [START_REF] Schapire | Boosting and rocchio applied to text filtering[END_REF], ...).

In 2012, a new TREC task called "Knowledge Base Acceleration" (KBA) [START_REF] Frank | Building an entity-centric stream filtering test collection for trec[END_REF] started with a more entity centric definition: filtering a time-ordered corpus for documents that are related to a set of entities from Wikipedia. The best performing approach used one classifier (SVM) by entity tracked with features representing whether or not a term is in a document, regardless of its frequency [START_REF] Kjersten | The hltcoe approach to the trec 2012 kba track[END_REF]. Training data have however to be provided for each new entity "followed". Another successful approach capture intrinsic characteristics of related documents by relying on document centric features, entity profile related features and time features (Bonnefoy and al, 2013).

Recently information filtering on Twitter emerged. [START_REF] Lee | Mining spatio-temporal information on microblogging streams using a density-based online clustering method[END_REF] for instance followed the evolution of big and short terms events, like natural disasters, in real-time. RepLab 2012 Filtering task [START_REF] Amigò | Overview of RepLab 2012: Evaluating Online Reputation Management Systems[END_REF] follows the KBA 2012 definition but focus on Twitter as the source of incoming data (instead of news, blogs and forum posts). The submitted approaches rely on various sources of evidence like named entity recognition [START_REF] Villena-Román | DAEDALUS at RepLab 2012: Polarity Classification and Filtering on Twitter Data[END_REF], matches of terms between tweets and Wikipedia [START_REF] Younus | CIRGDISCO at RepLab2012 Filtering Task: A Two-Pass Approach for Company Name Disambiguation in Tweets[END_REF] or the importance of features specific to Twitter such as the presence of a user name in a tweet [START_REF] Peetz | From Sentiment to Reputation ILPS at RepLab 2012[END_REF] or the number of hashtags [START_REF] Chenlo | FBM-Yahoo! at RepLab[END_REF].

Merging metrics or methods used in natural language processing (NLP) and information retrieval can be seen, as shown in (Lamontagne and Abi-zeid, 2006), as a multi-criteria optimization problem: in particular, the ELECTRE methods [START_REF] Figueira | Multiple Criteria Decision Analysis: State of the Art Surveys[END_REF], which turned out to be efficient applied to industrial domains [START_REF] Gourion | Aide à la décision robuste pour la localisation d'un centre de traitement des déchets. Comparaison de méthodes d'analyse multicritères[END_REF], have been transposed to an NLP context (Carrillo and al, 2012) opted for a voting method to combine their runs with [START_REF] Chenlo | FBM-Yahoo! at RepLab[END_REF].

Methods

Cosine distance (TF-IDF-Gini)

The first approach consists in a supervised classification based on a cosine similarity. Vectors used to compute similarities are built using the Term Frequency-Inverse Document Frequency (TF-IDF) [START_REF] Salton | Term weighting approaches in automatic text retrieval[END_REF] and the Gini purity criterion [START_REF] Torres-Moreno | Opinion detection as a topic classification problem[END_REF].

Tweets are cleansed by removing hypertext links and punctuation marks, hashtags and @ before a user name. We have removed a set of tool-words and some entities ID. Terms are lower-cased. We generate a list of n-grams by using the Gini purity criterion.

We create terms (words or n-grams) models for both classes (related and unrelated tweets) and term frequencies are computed with the TF-IDF and Gini criterion. These models take into account the following meta-data: user id, entity id and language integrated as terms in the bag-of-terms of tweets.

A cosine similarity measures the distance between the bag-of-terms of a tweet and the whole bag built for each class and ranks tweets according to this measure.

KNN with discriminant features

The system tries to match each tweet in the test set with the K most similar tweets in the training set. Tweet similarity is computed using Jaccard measure on the bag-of-words discriminant representation of the tweets. As in section 3.1, each tweet is represented as a vector whose components are weighted according to TF-IDF and the Gini purity criterion. The process also takes into account tokens created from the meta-data (author, entity-id). The stoplist of section 3.1 has been used.

Adaptation KBA'12 system

For the KBA filtering task, a state-of-the-art approach consist in capturing intrinsic characteristics of highly relevant documents by mean of three types of features: document centric features, entity's profile features, and time features (Bonnefoy and al, 2013). Features are computed for each candidate document and, using a Random Forest classifier, used to determine if the document is related or not to a given entity.

Unlike previous approaches it doesn't require a new set of examples for each new entity. We want to measure the robustness of this approach by using it on another type of documents (i.e. tweets). No adjustments are made on it but tweets are however preprocessed: stop-words are deleted as well as @ before user names and hashtags are split. The classifier is trained on all related and unrelated examples for each type of entities (automotive, universities, banking and music/artists).

Merging algorithms

To improve the performances we use three ways of combining our systems outputs.

Linear combination of outputs score

N systems are available. For each tweet T of the test set, one system j associates each label L k with a confidence score s j (T, L k ) (j = 1, ..., N ). The output entity label L is chosen according to the following rule :

L = arg max k   N j=1 s j (T, L k )   (1)

ELECTRE I method

The goal of this method [START_REF] Roy | The outranking approach and the foundations of ELECTRE methods[END_REF] is to choose the best label from the entire set of labels ranked according to the different systems.

A relation S ⊂ L × L, denoted "over-ranking", is defined on the label set L: a label l over-ranks another label l ′ if l dominates l ′ on an "important" number of systems and if l ′ does not dominate "too much" l on the remaining systems.

More precisely, for each pair of labels (l, l ′ ), a concordance index c(l, l ′ ) is computed, corresponding to the proportion of systems where l dominates l ′ . l over-ranks l ′ if c(l, l ′ ) exceeds a concordance threshold, generally fixed around 2/3 and if l is not dominated by l ′ on the remaining systems above a veto threshold, which has been fixed here to v = 0.5.

The set of the best labels, possibly empty and denoted as the kernel of the relation S, consists in the labels which are not overanked by others. If there is no, or more than one, label in the kernel, this method is discarded and the merging algorithm described in the previous subsection, based on a linear combination of the scores, is applied.

PROMETHEE mono-criterion method

This method relies on a concordance matrix: for each pair of labels (l i , l j ), the matrix coefficient c ij corresponds to the concordance index c(l i , l j ) introduced in the previous subsection.

For each label l i , two sums are computed: s l (l i ) = j c ij and s c (l i ) = j c ji . s l (l i ) measures the tendency of l i to dominate the other labels, and s l (l i ) the tendency of l i to be dominated.

The final score of the label l i is the difference s l (l i )s c (l i ) and the dominant label is the one whose score is maximal.

Experiments

Replab 2013 Framework

The corpus is a bilingual (English and Spanish) collection of tweets containing the name of one of the 61 entities selected in four domains: automotive, banking, universities and music/artists. Tweets have been collected by querying the Twitter search engine1 . The dataset covers a period going from the 1 st of June 2012 to the 31 st of December 2012. 42,700 tweets have been provided for training purpose and 100,000 tweets for the evaluation. The training set is composed of the 700 first tweets retrieved for each entity. For each entity, at least 2,200 tweets have been collected.

Tweets, however, are not homogeneously distributed across the entities.

Systems are evaluated according to the following measures: Accuracy, Reliability and Sensitivity (Amigò and al, 2013). Reliability is defined as precision of binary relationships predicted by the system with respect to those that derive from the gold standard; and Sensitivity is similarly defined as recall of relationships. A F-measure is then used to combine both scores.

These measures are well adapted to the task but are really severe on unbalanced datasets.

Results

Table 1 shows results of our approaches against the official RepLab 2013 baseline and the median sys-tem among participants.

The baseline 2 is a supervised system that matches each tweet in the test set with the most similar tweet in the training set, and assumes that the annotations in the tweet from the training set are also valid for the tweet in the test set. Tweet similarity is computed using Jaccard distance and a straightforward bag-of-words representation of the tweets.

The method described in section 3.2 can be considered as an improved version of the baseline.

Two systems (KNN and KBA with a F-measure scores of respectively .381 and .341) have reached greater performances than the baseline on every measures. The confidence interval (.002 and .005 respectively for accuracy and F-Measure) computed following Polling Method [START_REF] Voorhees | Variations in Relevance Judgements and the Measurement of Retrieval Effectiveness[END_REF] shows that the difference between the systems is significant.

Merging strategies R-Elec (for ELECTRE) and R-LC (for Linear Combination) did not produce good selection rules since their performances remain lower than the best system taken alone.A natural merging strategy consisting in merging only the best systems on a development set gives better results (Naive LC and Naive Elec).

Moreover, a multi pass strategy (MPMS) merging systems in pair before considering merging all pairs improves Sensitivity and thus the F-Measure (.400) despite of a loss in term of accuracy and reliability. Finally, merging only the best (OTB) runs on each measure gives quite similar improvements.

These results show that using merging strategies to combine different systems lead to improvements, whatever the metric chosen. The key observation is that it is possible to pick a merging strategy according to the metrics we choose to focus on. A quite naive merging strategy (Naive) seems to result in a better precision (improvements in both Accuracy and Reliability). On the contrary, adopting a multi pass strategy (MPMS) allow to give a highest priority to recall in both classes (i.e. Sensitivity). Finally, if a compromise is preferred, we saw that promoting systems that did well on each measure (OTB) is a good option.

2 http://www.limosine-project.eu/events/replab2013

Conclusion

In this paper we presented some of the interesting features of the systems that we evaluated within the framework provided by RepLab 2013 as well as their performances. We proposed several combinations of them using different merging strategies in order to take benefit from the diversity of information offered by our systems. We also showed that these merging strategies have to be applied depending on the evaluation measures to offer in one hand the best results according to a specific measure or in the other hand to obtain a trade-off. Since a merging strategy cannot get the best score according to each metrics, we can accept a loss according to one metric if it has a real impact on the task official measures.

A more advanced view would be to apply a specific merger entity by entity, especially for unbalanced entities.

Table 1 :

 1 Results on the Filtering Task ordered according to the F-Measure.

	Approach	Accuracy	Reliability	Sensitivity	F-Measure
	MPMS	.899	.668	.367	.400
	OTB	.902	.651	.367	.386
	Naive LC	.904	.691	.364	.385
	Naive Elec	.903	.671	.363	.383
	k-NN	.890	.658	.357	.381
	KBA	.878	.619	.331	.341
	R-LC	.895	.680	.290	.313
	Baseline	.876	.461	.325	.312
	R-Elec	.892	.680	.281	.302
	Cosine	.834	.423	.331	.272
	RepLab Median	.826	.489	.286	..265

http://twitter.com/search