Stavros Alchatzidis
email: stavros.alchatzidis@ecp.fr

Aristeidis Sotiras
email: aristeidis.sotiras@uphs.upenn.edu

Evangelia I Zacharaki
email: evangelia.zacharaki@centralesupelec.fr

Nikos Paragios
email: nikos.paragios@ecp.fr

A Discrete MRF Framework for Integrated Multi-Atlas Registration and Segmentation

Keywords: multi-atlas segmentation, medical imaging, Markov Random Fields, discrete optimization

Multi-atlas segmentation has emerged in recent years as a simple yet powerful approach in medical image segmentation. It commonly comprises two steps: i) a series of pairwise registrations that establish correspondences between a query image and a number of atlases, and ii) the fusion of the available segmentation hypotheses towards labeling objects of interest. In this paper, we introduce a novel approach that solves simultaneously for the underlying segmentation labels and the multi-atlas registration. The proposed approach is formulated as a pairwise Markov Random Field, where registration and segmentation nodes are coupled towards simultaneously recovering all atlas deformations and labeling the query image. The coupling is achieved by promoting the consistency between selected deformed atlas segmentations and the estimated

Introduction

Segmentation, or the process of assigning voxels to distinct anatomical regions or tissue types, is a fundamental task in medical image analysis. The accurate delineation of anatomical structures is the cornerstone of quantitative analysis that aims, among other, to understand normal and diseased anatomical variability. The increasing availability and size of high resolution imaging data along with the widespread adoption of imaging across clinical and research practice further underline the importance of reliable, efficient and accurate image segmentation. While manual segmentation is the gold standard, the complexity of the task, as well as the high time requirement and the associated cost, make it prohibitive on a large scale. In such a setting, automatic image segmentation provides a valuable alternative.

As a consequence, important research efforts have been focused on developing automatic segmentation algorithms. Among the wealth of the developed techniques, segmentation via registration [START_REF] Rohlfing | Quo Vadis, Atlas-Based Segmentation?[END_REF] stands as a unique example in medical image processing. Registration is used to map the grayscale image of an atlas to the query image, while the estimated deformation is subsequently used to propagate the available labels and provide an estimate of the segmentation. However, a single atlas is limited with respect to its ability to capture the variability of the population.

Extending registration-based segmentation by incorporating multiple atlases has thus emerged as a natural extension towards tackling the aforementioned limitation [START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF]. Multi-atlas segmentation, fueled by the maturity of the available registration techniques and advancements in computer hardware that partly alleviate its high computational cost, has gained significant popularity and found numerous applications in medical image analysis. For example, one may cite full brain segmentation [START_REF] Akhondi-Asl | Simultaneous Truth and Performance Level Estimation Through Fusion of Probabilistic Segmentations[END_REF][START_REF] Artaechevarria | Combination Strategies in Multi-Atlas Image Segmentation: Application to Brain MR Data[END_REF][START_REF] Asman | Robust Statistical Label Fusion Through Consensus Level, Labeler Accuracy, and Truth Estimation (COLLATE)[END_REF][START_REF] Iglesias | A unified framework for cross-modality multi-atlas segmentation of brain MRI[END_REF][START_REF] Rousseau | A Supervised Patch-Based Approach for Human Brain Labeling[END_REF], skull-stripping [START_REF] Doshi | Multi-atlas skull-stripping[END_REF], hippocampal [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF][START_REF] Wang | Multi-Atlas Segmentation with Joint Label Fusion. Pattern Analysis and Machine Intelligence[END_REF], prostate [START_REF] Dowling | Fast Automatic Multi-atlas Segmentation of the Prostate from 3D MR Images[END_REF][START_REF] Klein | Automatic segmentation of the prostate in 3D MR images by atlas matching using localized mutual information[END_REF] and heart [START_REF] Isgum | Multi-Atlas-Based Segmentation With Local Decision Fusion—Application to Cardiac and Aortic Segmentation in CT Scans[END_REF][START_REF] Van Rikxoort | Adaptive local multi-atlas segmentation: Application to the heart and the caudate nucleus[END_REF] delineation.

Multi-atlas segmentation methods produce state of the art results in various settings by essentially applying a common pipeline. The pipeline comprises two basic steps. First, rigid or non-rigid registration is performed in order to align the atlas images to the query one and propagate their labels forming a set of candidate segmentations. Subsequently, the derived multiple segmentation hypotheses are fused to produce the final labeling of the query image.

One of the most popular fusion strategies is the majority voting one [START_REF] Rohlfing | Quo Vadis, Atlas-Based Segmentation?[END_REF]. Since early approaches in multiatlas segmentation, majority voting has been established as an intuitive, simple and robust baseline method. Many efforts have concentrated on developing more sophisticated fusion algorithms by incorporating additional information in terms of either local appearance [START_REF] Artaechevarria | Combination Strategies in Multi-Atlas Image Segmentation: Application to Brain MR Data[END_REF], neighborhood information [START_REF] Fonov | Multi-atlas labeling with population-specific template and non-local patch-based label fusion[END_REF][START_REF] Asman | Non-local statistical label fusion for multi-atlas segmentation[END_REF], or by exploiting the statistical correlation of inter-atlas errors to calculate optimal weights [START_REF] Wang | Multi-Atlas Segmentation with Joint Label Fusion. Pattern Analysis and Machine Intelligence[END_REF].

Alternative fusion techniques adopt probabilistic models to derive per voxel weights. Modeling the segmentation candidates as noisy observations of the true segmentation and using Expectation-Maximization estimation of the per atlas confusion matrix has been proposed in [START_REF] Warfield | Validation of Image Segmentation and Expert Quality with an Expectation-Maximization Algorithm[END_REF][START_REF] Rohlfing | Expectation Maximization Strategies for Multi-atlas Multi-label Segmentation[END_REF]. This method has been extended in many ways, for example including a smoothness prior [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF], using appearance-based ranking [START_REF] Cardoso | STEPS: Similarity and Truth Estimation for Propagated Segmentations and its application to hippocampal segmentation and brain parcelation[END_REF], incorporating non-local intensity proposals [START_REF] Asman | Non-local statistical label fusion for multi-atlas segmentation[END_REF] and performing iterative atlas selection [START_REF] Langerak | Label Fusion in Atlas-Based Segmentation Using a Selective and Iterative Method for Performance Level Estimation (SIMPLE)[END_REF]. In another related approach, a generative probabilistic model of label fusion was proposed in [START_REF] Sabuncu | A Generative Model for Image Segmentation Based on Label Fusion[END_REF].

A common principle behind most multi-atlas segmentation approaches is that registering the atlas images to produce candidate solutions and segmenting the query image are treated separately, in two independent steps. However, registration could benefit from taking into account the underlying segmentation towards establishing more accurate correspondences. Thus, ap-proaches that treat registration and segmentation through fusion as an inter-weaved process have recently appeared.

In [START_REF] Asman | Robust GM/WM Segmentation of the Spinal Cord with Iterative Non-local Statistical Fusion[END_REF] the authors segmented white and gray matter areas of the spinal cord by iteratively, affinely registering atlas masks to an estimated target segmentation, and fusing the registered masks and images to get a new estimate of the target segmentation. In a similar context, a probabilistic generative model based on the Demons registration framework was proposed by [START_REF] Iglesias | A unified framework for cross-modality multi-atlas segmentation of brain MRI[END_REF]. The model coupled registration by discouraging atlas to target transformations that disagree with atlas to precomputed mean template to target transformations. The model also imposed a distance transform based prior on the target mask, a membership field that corresponds voxels to atlases and a mixture of gaussians appearance model linking anatomical regions to target intensity probabilities. An approximate variational EM scheme was used to find the most probable model parameters. A similar work was presented in [START_REF] Tang | Bayesian parameter estimation and segmentation in the multi-atlas random orbit model[END_REF], where the Large Deformation Diffeomorphic Metric Mapping registration framework was employed.

Apart from the multi-atlas process our method can be seen in the same perspective as classical joint segmentation/registration methods. A first joint approach for voxel-based registration and segmentation is proposed in [START_REF] Wyatt | Noble: MAP MRF joint segmentation and registration of medical images[END_REF]. The joint class histogram is used to drive the Mutual Information metric for the registration step and a Gaussian mixture model is inferred to drive the segmentation step. Validation is carried out over rigid segmentations and 2D images. Application in medical images as well as a generative model to link the two processes can be found in [START_REF] Ashburner | Unified segmentation[END_REF] where intensity uniformity as well as spatial priors are introduced in the segmentation model. Use of MRFs can be seen again in the works of Xiaohua et al. first on MRI data [START_REF] Xiaohua | Simultaneous Segmentation and Registration for Medical Image[END_REF] and then in a formulation for contrast enhanced breast MRI scans [START_REF] Xiaohua | Simultaneous Segmentation and Registration of Contrast-Enhanced Breast MRI[END_REF]. Recent works include lung [START_REF] Xue | Joint registration and segmentation of serial lung CT images for imageguided lung cancer diagnosis and therapy[END_REF] and prostate segmentation [START_REF] Ou | Non-rigid registration between histological and MR images of the prostate: A joint segmentation and registration framework[END_REF], as well as brain lesions identification in MRI data [START_REF] Seghier | Lesion identification using unified segmentation-normalisation models and fuzzy clustering[END_REF][START_REF] Parisot | Joint Tumor Segmentation and Dense Deformable Registration of Brain MR Images[END_REF].

We complemented previous approaches for integrated registration and label fusion segmentation in [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] by introducing a discrete formulation based on Markov Random Field theory. Latent variables include the displacements of the grid nodes of a B-Spline transformation model as well as voxel segmentation variables. Segmentation additionally takes into account class likelihoods produced by a discriminatively trained classifier. Constraints were imposed by taking into account how congruent the proposed segmentations are with respect to the proposals of the rest of the atlases as well as the classifier produced likelihoods. Registration and segmentation variables are coupled by using an appearance-based weighting similar to the one used in local fusion strate-gies. As a consequence, as votes are weighted, atlases that do not match well locally will have a minor contribution to the inferred segmentation mask, resulting in a local "soft" atlas selection scheme.

In this work, we build upon the work presented in [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] by extending the formulation, providing a more detailed description of the method and reporting a more extensive validation setting. In this work, we estimate membership fields by introducing a local atlas selection scheme. This scheme, explicitly models variables for selecting parts of each atlas labeling by comparing them directly to the estimated underlying segmentation. As a consequence, membership field images are directly produced by optimizing the model. The fact that the atlas selection is achieved by comparing segmentations, facilitates the use of images of different modalities to be part of the dataset.

The remainder of this paper is organized as follows.

In Section 2 we formulate the problem in the continuous domain, while in Section 3 we present the decomposition of the objective function toward a mapping into a discrete graph structure. Experimental validation is discussed in Section 4, while Section 5 concludes the paper.

Problem Formulation

We consider a dataset of N annotated images A = {A 0 , . . . , A N -1 }. Each image comes with a corresponding segmentation mask where the anatomical regions of interest have been annotated, forming the set S = {S 0 , . . . , S N -1 }. Each voxel in the segmentation image is assigned to a segmentation label corresponding to one of M anatomical classes, S i (x) ∈ {0, . . . , M -1}. In this paper, we refer to an atlas as the aggregation of an intensity image A i and its corresponding segmentation mask S i . Moreover, we consider that an image I is given as input to be segmented into anatomical regions. Henceforth, we are going to interchangeably refer to this image as either target or query image. The output of the proposed algorithms comprises: i) a set of membership field images F = {F 0 , . . . , F N -1 }, F i (x) ∈ {0, 1} denoting if an atlas influences a point x in the query image; ii) the segmentation mask S I corresponding to the target image; and iii) a set of deformation fields D = {D 0 , . . . , D N -1 }, where D i denotes the deformation field mapping of A i to I.

Method Outline

The goal of the proposed method is to simultaneously solve for the parameters of the final segmentation S I of the query image and the set of deformation fields D. Hence, the agreement of the deformed segmentation masks with the underlying estimated segmentation may be taken into account during the estimation of the deformation fields, leading to more accurate correspondences and consequently, improved segmentation. Our basic premise is that by allowing the two problems to interact, the quality of the respective solutions will rise due to the additional available information.

Furthermore, we aim to improve the final segmentation S I by taking into account class specific appearance priors. The motivation behind incorporating prior segmentation probabilities lies in that fact that image registration is often trapped in local minima when trying to match areas of high anatomical variability (e.g., brain cortex). In such a setting, appearance information constitutes an alternative, more reliable cue that can robustly guide segmentation [START_REF] Sdika | Combining atlas based segmentation and intensity classification with nearest neighbor transform and accuracy weighted vote[END_REF].

Continuous Energy Formulation

We formulate the problem as an energy minimization one. The proposed energy consists of three components: i) a registration component comprising a matching term (M ), that quantifies the level of alignment between each atlas and the query image, and a regularization term (R d ) that enforces the smoothness of the deformation field; ii) a segmentation component comprising an appearance prior term (S P ), that measures the loglikelihood of the segmentation with respect to the probabilities (π) learned during a training phase. iii) a coupling term (C) that takes into account the labeling that is proposed by the atlases, over the domain indicated by F, and encourages their agreement with the estimated segmentation S I . Finally, a regularization term (R f ) is imposed on the membership fields allowing the atlases to influence the derived segmentation in a smooth spatially varying fashion.

The energy has the following form:

E(D, A, S, S I ) = M (D, A, I) + R d (D) Registration + S P (S I , π) Segmentation + C(D, F, I, S, S I ) + R f (F) Coupling (1)
The first two energy terms correspond to the standard energy that is commonly minimized in multi-atlas segmentation frameworks, while the third term is common in segmentation frameworks. The fourth term introduces the main novelty of this work, i.e. the coupling between the segmentation and the multi-atlas registration. Let us now detail each term of the previous energy.

Registration

The registration component aims to align all atlases to the query image. As such, it is a generalization of standard pairwise image registration.

Matching criterion Given any dissimilarity intensitybased criterion ρ, the matching term takes the following form:

M (D, A, I) = N -1 i=0 Ω ρ(I, A i • D i (x))dx. ( 2
)
This term is the summation of the independently evaluated dissimilarity criteria between all atlases and the query image I.

Deformation smoothness Image registration is an illposed problem according to Hadamard's definition. In order to account for this fact, regularization is necessary. Typical regularization settings penalize the nonsmoothness of deformation fields. Thus, given a smoothness inducing function ψ, the regularization term takes the following form:

R d (D) = N -1 i=0 Ω ψ(D i (x))dx. (3)
In other words, this term evaluates the smoothness of all deformation fields mapping from an atlas to the target image and sums the independent evaluations.

Transformation model In this work, the popular Free Form Deformations (FFDs) transformation model [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF][START_REF] Kybic | Fast parametric elastic image registration[END_REF] is used. Free Form Deformations parametrize the transformation D j (x) by a linear combination of K control points:

D j (x) = x + K-1 i=0 ω i (x)φ i , (4)
where φ i is the displacement of control point i and ω i (x) is an interpolation or weighting function that determines the influence of the control point i to the image point x. In the current approach, we use N uniformly distributed grids of control points superimposed over the image domain (one corresponding to each atlas) and cubic B-splines as the weighting functions.

Segmentation

Segmentation Prior The quality of the segmentation hypotheses that are provided by warping the given atlases is conditioned upon the quality of the registration. Image registration is often overwhelmed when trying to establish correspondences between highly variable anatomical regions leading to inaccurate results.

In such cases, one can exploit additional cues to enhance segmentation estimation. Local appearance provides complementary information that can be incorporated into segmentation. Assuming a probability function on the candidate labeling of the form π x (l) where l ∈ {0, . . . , M -1}, we aim to penalize all segmentations that go against prior information:

S P (S I , π) = Ω -log(π x (S I (x)))dx. ( 5
)
Such probabilities can be learned using any modern classification method.

Coupling

In this work, we close the circuit between multi-atlas registration and label fusion. We allow segmentation to influence registration by imposing agreement constraints between the estimation of the underlying segmentation and the atlas hypotheses. As a consequence, improved segmentation accuracy may be achieved through the refinement of the registration result. The above is modeled by means of introducing a third term in the energy. This term penalizes in a controlled manner deformations that lead to disagreement between the hypotheses and the estimated segmentation. Specifically the penalties are not defined in advance but are updated based on the segmentation consistency and smoothness of membership fields at each iteration.

C(D, F, I, S, S I ) = N -1 i=0 Ω F i (x)κ(S i • D i (x), S I (x))dx, (6)
where κ(a, b) is a penalizing disagreement function, for which we assume κ(a, a) = 0.

Membership field smoothness

We assume that only a spatially varying subset of the dataset is pertinent to deducing the correct segmentation mask at any point. Dealing with images, an important natural prior is spatial smoothness:

R f (F) = N -1 i=0 Ω ψ(F i (x))dx, ( 7
)
where ψ is again a smoothness inducing function.

Markov Random Field Formulation

We use Markov Random Field (MRF) theory to formulate the above minimization problem in a discrete context. The problem is represented by a graph G = (V, E), where V denotes the set of nodes that encode the latent variables, and E the set of edges that encode the interactions between the variables. The graph is associated with an energy of the form:

E M RF (l) = p∈V g p (l p ) + (p,q)∈E f pq (l p , l q ), (8)
where random variables p take values from a discrete set of solutions L, g p (l p ) measures the cost of assigning a value l p to the variable p and f pq (l p , l q ) is a pairwise function that determines the cost of assigning different values l p and l q to the variables p and q.

Graph Structure

The constructed graph should encode the multi-atlas registration, the segmentation and the constraints that integrate the two problems. Let us now detail how the graph is constructed to achieve this.

Multi-Atlas Registration

Let us recall that the deformation model is parametrized by N deformation grids. This is encoded in the MRF graph G by a set of N isomorphic grid graphs G D = {G D0 , . . . , G D N -1 }. For every control point in the deformation grid that is superimposed onto image A i , there is a node p i ∈ V Di that represents its displacement. Since grids are isomorphic, p indexes a common control point position, while i indexes the grid. The edge system of each grid E Di is created by assuming a regular 6-connectivity scheme. The label set L D for this set of variables is a quantized version of the displacement space. A label assignment

l d pi ∈ L D (with p i ∈ V Di ) is equivalent to displacing the control point p i by displacement d p i .
Segmentation An additional graph G S = (V S ) is employed to model segmentation. A node p s ∈ V S encodes a random variable and corresponds to a voxel in the target image whose position is indexed by the subscript s.

We should also emphasize the fact that the nodes that form the segmentation graph are not connected to one another. The set of possible solutions L S represents the set of anatomical regions augmented by the background label. We refer to a potential anatomical label in L S by l s .

Coupling Integrating segmentation and multi-atlas registration is achieved by coupling segmentation and deformation graphs. The set of edges E C connects nodes of V S with nodes of V D . In order to create the coupling edge system, we connect every node p ∈ V D with nodes of V S that correspond to voxels that are influenced by a displacement of p.

Local Atlas Selection

We parametrize membership fields by taking into account the spatial support of the deformation nodes. Voxels within the support of a control point share the same membership state ([0, 1]). To model this, we augment the label set of the deformation nodes by considering the Cartesian product between the deformation label set L D and a binary selection label set L E = {0, 1}. Thus, for a node p a label l d indexes a pair (d p , e d p ), where e d p ∈ L E . A node p is selected when e d p = 1, otherwise it is deselected. If p is deselected, it will not penalize inconsistent candidate segmentations and it will not be influenced by them. In the following section, only the relevant part of the pair (d p , e d p ) appears in the right hand side of the equations.

MRF Energy

The continuous energy in Eq. 1 is mapped to a discrete MRF energy of the form in Eq. 8. In short, we map i) the matching term M (Eq. 2) to the unary potentials of the deformation variables (Eq. 9), ii) the deformation smoothness penalty term R d (Eq. 3) to pairwise potentials between deformation variables (Eq. 10), and iii) the coupling penalty C (Eq. 6) to one pairwise potential between registration and segmentation variables and one unary potential over deformation variables (first and second part of the right hand side of Eq. 12, correspondingly), as well as iv) the membership field smoothness penalty term R f to an additional pairwise potential between deformation nodes (Eq. 13).

Multi-Atlas Registration

Multi-atlas registration is performed by registering in a pairwise fashion all atlases to the target image. Formulating pairwise registration in a discrete setting has been shown in [START_REF] Glocker | Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods*[END_REF]. For completeness reasons, we briefly discuss here how the matching term M and the regularization term R of Eq. 1 are mapped to unary and pairwise potentials.

As far as the matching term is concerned, we are interested in quantifying how well the assignment of a displacement label l d pi ∈ L D to a node p i ∈ V Di aligns atlas A i to the target image. This is measured by the following unary potential:

g M pi (l d pi ) = Ω ωpi (x)ρ(A i • D pi , I(x))dx. ( 9
)
D pi is the transformation induced by the movement of the control point p in the i-th deformation grid by the displacement l d pi . The weighting function ωpi determines the contribution of the point x to the unary potential of the control point p. This function is similar to the ω weighting functions used in the FFD deformation model (Eq. 4).

Regarding the regularization term, [START_REF] Glocker | Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods*[END_REF] shows that it can be efficiently modeled by pairwise potentials. A discrete approximation of the gradient of the spatial transformation can be computed by taking the vector difference between the displacements of neighboring nodes that belong to the same deformation grid:

f R piqi (l d pi , l d qi ) = d pi -d qi , (10)
where d pi is the displacement applied to node p in the i-th deformation grid, indexed by l d pi .

Segmentation

In order to assign a class label to every voxel of the target image, we take into account learned appearance model for every class. The appearance model is encoded in the form of a probability distribution π(l) and can be naturally incorporated in the MRF model by setting the unary potentials of the segmentation grid for every label to the negative log-probability of the respective class:

g SP q S (l s q S ) = -log(π(l s q S )).

Integrated Segmentation and Multi-Atlas Registration

We want to encourage the agreement between the estimated segmentation and the warped segmentation masks. Thus, we penalize control point displacements of grid G Di that result in warping the segmentation mask of the corresponding atlas i in a fashion that does not agree with our final segmentation:

f C piq S (l d pi , l s q S ) = ωpi (s) • Ind(S i • D pi (s), l s q S ), (11)
where p i belongs to the grid G Di , q S belongs to G S and Ind(x, y) = 0 when x = y 1 when x = y.

Algorithm 1 Agreement percentages estimation algorithm.

1: Define a p i = Ω ωp i (x)(1-Ind( Ŝ(x),S i (x))dx) Local Atlas Selection Due to the sparse way we model selection variables, any candidate deformation d p corresponds to a segmentation mask that agrees at a certain percentage to the segmentation variables within the support of p. In order to model local atlas selection, we need to determine a set of such percentages below which a candidate segmentation is unacceptably incongruent with the consensus segmentation, and thus should lead to a node being deselected. We call these percentages agreement percentages and index them with a p . Note that the agreement percentage does not depend on the grid index i but only on the inter-grid control point index p. Given an agreement percentage a p , atlas selection may be enforced by introducing an additional unary cost for deformation nodes. Thus the total coupling cost becomes

Ω ωp i (x)
f A S piq S (l d pi , l s q S ) = e d pi • f C piq S (l d pi , l s q S ) + (1 -e d pi ) • (1 -a p ), (12)
where e d pi is equal to 1 when p i is enabled, and 0 otherwise. As a consequence, nodes for which the (pairwise) segmentation cost f C is very high when warping an atlas i, are disabled for this atlas. When a node p i is disabled the cost to be paid is 1 -a p regardless of the level of disagreement of the deformation with its corresponding segmentation nodes. Thus, a disabled deformation node will not affect segmentation variables and conversely will not be affected by them.

Agreement estimation Agreement percentages are estimated by locally comparing each segmentation mask with an estimate of the consensus segmentation, followed by a clustering scheme to arrive at robust per control point estimates. This process is summarized in Algorithm 1. First, agreement percentages per control point are computed (step 1). Then a 2-means clustering algorithm is used to separate them into inliers and outliers (steps 3-10). Finally, a variance reduction criterion is used to determine the validity of the clustering into two sets and a p is defined accordingly (steps [START_REF] Dowling | Fast Automatic Multi-atlas Segmentation of the Prostate from 3D MR Images[END_REF][START_REF] Iglesias | A unified framework for cross-modality multi-atlas segmentation of brain MRI[END_REF][START_REF] Fonov | Multi-atlas labeling with population-specific template and non-local patch-based label fusion[END_REF][START_REF] Glocker | Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods*[END_REF][START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF].

Local Atlas Selection Smoothness Membership fields are encoded over deformation nodes. We want to enforce smoothness over the fields to achieve concise regions of influence for every atlas:

f S pi,qi (l d pi , l d qi ) = Ind(e d pi , e d qi ). ( 13
)

MRF energy parameterization

To conclude the outline of the discrete energy, we summarize the terms along with the parameters controlling each term's weight:

E M RF (l) = g M pi (l d pi ) + λf R piqi (l d pi , l d qi ) + σg SP q S (l s q S ) +αf A S piq S (l d pi , l s q S ) + βf S pi,qi (l d pi , l d qi ), (14)
λ controls the deformation field smoothness, σ encodes the prior likelihood weight, α specifies how much registration influences segmentation and reversely how much segmentation affects registration, and β regulates the smoothness of membership fields.

MRF Optimization through Dual Decomposition

DD-MRF [START_REF] Komodakis | MRF energy minimization and beyond via dual decomposition[END_REF] has been introduced as a framework for MRF optimization, offering global optimality guarantees. Its flexibility in terms of possible energy types, its ability to report the quality of the final solution as well as its optimality guarantees are the merits we considered in opting for its use. DD-MRF works by receiving as input a decomposition of the initial graph (primal problem) into subgraphs (dual problems). It initializes the costs of the dual problems using the costs of the primal problem. It then proceeds by iteratively finding a global optimum for each subproblem, compare the subproblem solutions and update their costs.

In short, at each iteration: i) dual subproblems are solved in an optimal fashion (usually by dynamic programming); ii) the dual energy, defined as the sum of the energies of the optimal dual solutions, is computed; iii) a solution to the primal is inferred by the multitude of, possibly conflicting, dual solutions; iv) the difference between the energy of the primal and the sum of the duals is computed (referred to as primal-dual gap); v) the primal dual gap is used to update subproblem costs. This way agreement is induced between the next subproblem solutions leading iteratively to a coherent, globally optimal solution. The way dual costs are updated guarantees that the euclidean distance of the current solution to the set of the globally optimum solutions will decrease monotonically.

Subproblem decomposition

In order to optimize a graph G, DD-MRF requires an input set of subproblems SD = {SD i , • • • SD n }, such that i SD i = G. Thus, we decompose the problem into a series of subproblems that can be exactly optimized through dynamic programming. Deformation grid subgraphs G Di are decomposed into sets C Di where every element of C Di corresponds to a chain subproblem.

The coupling term is divided into tree subproblems. There is one such subproblem per deformation node consisting of the deformation node p, the segmentation nodes within its support to which it is connected, as well as the edges connecting them. We call such a tree subproblem T p and T Di = {T p |p ∈ G Di } denotes the tree subproblems that include all the nodes of a grid subgraph G Di . Thus, the subproblem decomposition for the problem at hand is given by the set SD = {T D0 , . . . , T D N , C D0 , . . . , C D N }.

Local atlas selection model optimization However, a direct one-shot optimization of the above decomposition is not viable for a problem of this size and type. The inclusion of the switching potential that is encoded in the coupling cost makes optimization very hard. This has a crippling effect on both the time required for convergence and the quality of the resulting solution. In order to facilitate the solution of the problem, we further decompose the problem by breaking SD into N parts which we will optimize sequentially using DD-MRF, while updating the remaining parts with the solutions of the parts that have been already solved.

This iterative procedure estimates the deformation and selection fields for each atlas sequentially. For a given atlas i we project the coupling potentials of every other atlas by considering them stationary (zero vector displacement label for every deformation node). We then proceed by solving the single atlas decomposition over the updated segmentation potentials.

The algorithm consists of the steps shown in Alg. 2. The process starts by computing the potential functions for all candidate solutions. It then iterates over the atlases, projecting the updated constraints to the segmentation grid (step 3), optimizing the corresponding part of the decomposition to infer deformation variables

% deselected nodes

Fig. 1 Left: Closeness ranking of atlases to target IDs for leave-one-out cross-validation. Lighter boxes signify closeness of an atlas to a target. Right: Percentage of deselected nodes per atlas with respect to target ID as returned by our method when using all 17 atlases. The reported percentages are much more significant than what they appear, considering that 30%-40% of control points have a support that is essentially background. See Fig. 5 for a visual representation. Elements in the diagonals should be ignored as naturally no atlas was used to segment itself.

and the best current segmentation labeling (step 4), updating deformation fields and selection images (step 5) and the corresponding constraints (step 6). The final segmentation is the one produced by the last iteration.

Algorithm 2 Sequential model optimization.

1: Define the potentials as described in previous sections. 2: for each atlas k do 3: Define g P roj q S (l s q S ) = i∈N \k f C p i q S (0, l s q S ) over G S 4:

Optimize using DD-MRF the subproblem set

{T D k , C D k } 5:
Update the deformation field D k and membership field image F k from optimal labeling of G D k . 6:

Recalculate f C p k q S (0, l s q S ) 7: end for

Validation

We have validated our method using leave-one-out crosscorrelation over publicly data available on the Internet Brain Segmentation Repository (IBSR). We specifically use the skull stripped version of the dataset provided in [START_REF] Rohlfing | Image similarity and tissue overlaps as surrogates for image registration accuracy: widely used but unreliable[END_REF]. The dataset consists of 18 T1-weighted MR Images with 1.5mm slice thickness. Images and masks have been linearly registered and cropped to 145 × 158 × 123 from their initial resolution of 256 × 256 × 128. We report mean results over labels that were annotated in more than half of the images of the dataset. Results are reported over 14 symmetric (left, right hemisphere) annotations (Thalamus, Caudate, Putamen, Pallidum, Hippocampus, Amygdala, Accumbens, Lateral Ventricle, Inferior Lateral Ventricle, Cerebral White Matter, Cerebral Cortex, Cerebellum White Matter, Cerebellum Cortex, vessel) and 3 non-symmetric (3rd Ventricle, 4th Ventricle, Brain Stem).

Given a target image, all other images have been registered using affine registration and compared using Normalized Cross Correlation (NCC) to create a closeness ranking (see Fig. 1 (left) ). Experiments using less than the full dataset ( N < 17 atlases), refer to this ranking to select the N closest atlases to the target image. This ranking is also used by the sequential optimization detailed in Sec. 3.3 to determine the sequence of the atlases, from the closest to the furthest.

Registration parameters

The same iterative procedure as [START_REF] Glocker | Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods*[END_REF] has been used in order to best cope with the computational efficiency vs accuracy trade-off. Dense deformation fields are produced by control point displacements using Cubic B-Spline functions. For all registrations we use two levels of deformation control points with spacing 7mm and 3.5mm. We use the deformation fields produced by the coarse level to initialize the finer one. For deformation nodes, candidate displacements are uniformly sampled over each axis, 12 per axis for a total of 37. For each control point level we iterate 5 times, reducing with each iteration the area of candidate displacements by a factor of 0.66. Normalized Cross Correlation (NCC) is the dissimilarity function (ρ) used for the matching criterion in all experiments. Segmentation Likelihoods Training has been conducted by removing the target image from the training set. The following configuration has been used to learn prior per voxel likelihoods. For each training set image, we sampled up to 150 samples belonging to each label and then sampled in a spatially uniform manner another 15000 samples avoiding duplicates. Three types of features were used: i) median, entropy, standard deviation, kurtosis and skewness sample statistics on a 3D patch with sides of 5, 7 and 9 voxels; ii) Gabor features using 6 per axis orientations and 3 scales; iii) HOG3D features on a 11 × 11 × 11 patch broken up to 8 (2 per dimension) subpatches using 4 orientations; and iv) normalized voxel positions and distance from the center voxel. The total size of the feature vector was 258. We used the Random Forest framework [START_REF] Breiman | Random Forests[END_REF] to discriminatively learn local per voxel probabilities for our target image.

We used 200 trees of maximum depth 20 and opted to stop splitting at 20 samples per node.

Optimization convergence Concerning DD-MRF convergence, we stop iteration when the mean change rate of the primal dual gap computed over the last 10 iterations reaches the 10% of the mean change rate of the primal dual gap computed over all iterations. This criterion leads to a uniform optimization quality over all possible target images, number of atlases and models.

Methods We compare our method against: i) pairwise registration fused using majority voting referred in the figures as Majority, ii) pairwise registration fused using appearance-based locally weighted fusion. We use an exponential weighting function over a local patch: e -ρN CCs(x) where ρ is a slope parameter, N CC s (x) the NCC comparison of an atlas and a target patch around Parameter estimation Parameters were estimated by maximizing the Dice coefficient when segmenting 4 handpicked target images by using 4 different datasets that consist of 7 atlases. First, the parameter λ was found for the pairwise registrations (λ = 0.01). For the appearancebased locally weighted fusion the slope parameter was found (ρ = 3) and the patch size was set (s = 9). Keeping λ the same, the parameters σ and α were set for [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] (σ = 1.5, α = 0.3). Finally σ and α and β were set for the proposed method (σ = 2.5, α = 0.1, β = 0.025). The parameter introduced in Alg. 1 was set to 3

N 2
where N the number of atlases in the dataset.

Results

In this section, we validate the proposed method through qualitative and quantitative comparisons with competing methods. Moreover, we opt to shed light upon the mechanics of our method and in particular, how local atlas selection affects the segmentation framework.

When we refer to a specific atlas id x out of the dataset we will use the notation A x . In terms of quantitative results the complementarity of overlap and surface distance measures has been recently made evident in [START_REF] Rohlfing | Image similarity and tissue overlaps as surrogates for image registration accuracy: widely used but unreliable[END_REF]. Overlap is measured by the Dice coefficient that is defined as Dice(A, B) = In Fig. 2 we present comparison results of segmentations produced using all 17 available atlases. The proposed method is consistently better than competing methods although [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] seems to perform in some cases (e.g A 2 ) particularly well. As far as Dice coefficient is concerned, the results are statistically significant. The proposed method outperforms both majority voting (paired t-test, P < 10 -5 ) and local appearance-based voting (paired t-test, P < 10 -5 ). Moreover, it performs better than [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] (paired t-test, P < 10 -2 ) despite not exploiting appearance similarity information to weigh the coupling term. Regarding the SMSD metric, the proposed method outperforms majority voting (paired ttest, P < 10 -3 ), local appearance-based voting (paired t-test, P < 10 -2 ) and [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] (paired t-test, P < 10 -2 ). In Table 1 mean Dice overlap, SMSD are presented for the methods considered.

To qualitatively appraise the quality of the obtained results, we show in Fig. 4 sagittal views of the ground truth segmentation as well as the results obtained with all examined methods. We note that the proposed method has retrieved more details in the highlighted areas. In Fig. 5, we show the ground truth segmentation and the corresponding views of segmentation proposals by the most and least deselected atlases along with their respective membership fields. The two atlases appear to propose very different candidate segmentations. The segmentation that is proposed by the most deselected atlas, differs in many locations when compared to the ground truth. The proposed method correctly estimates most areas where erroneous proposals have been made and correctly disables the membership field. Moreover, correctly estimating the membership field assures that atlases with bad matching will not act as noise in the coupled registration/segmentation scheme.

A close look at Fig. 1 (right) adds additional insight into the atlas selection mechanism. The figure shows us how the algorithm prunes atlases in a concise manner. For example, we see how A 9 is treated in average as an outlier, while A 0 is consistently part of the set of atlases that drives the method. Additionally, we can see a clear clustering pattern. A 0 -A 4 and A 12 -A 17 form a group in which atlases seem to segment one another, while the rest of the dataset seems to have a much more uniform node selection pattern.

Lastly, in Fig. 3 we address the issue of scaling with respect to the number of employed atlases. As far as the Dice coefficient is concerned (Fig. 3 (left)), [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] outperforms the proposed method when few atlases are used. This is due to the fact that the definition of the thresholds in Alg.1 is a noisy procedure.

As the number of used atlases increases, the mean quality of the dataset decreases. As consequence, our method starts to outperform [START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF] because of its ability to extract good quality matchings from the newly added atlases, without allowing them to affect already established matchings. An analogous pattern can be seen in Fig. 3 (right) for the case of SMSD.

Discussion

In this paper, we presented a method that integrates registration and segmentation fusion in a pairwise MRF framework. The proposed approach allows registration parameters to be updated based on the segmentation estimation while membership fields estimation assures that bad matchings will not deteriorate segmentation quality. The experimental results demonstrate the potential of the proposed approach.

Fig. 2

 2 Fig. 2 Top: Symmetric mean surface distance with respect to the ID of the target image. Bottom: Dice coefficient with respect to the ID of the target image. All 17 available atlases have been used to segment a target image. Majority refers to majority voting after pairwise registration, Local refers to the appearance-based locally weighted voting, Soft Selection refers to [2], while Hard Selection corresponds to the proposed method.

  2 * |A∩B| |A|+|B| over annotations A and B. Surface distance is measured by the Symmetric Mean Surface Distance (SMSD) defined as follows: If  is the set of surface voxels of volume A and d(A, B) = i∈A min j∈B e(i, j) |A| with e(i, j) being the euclidean distance, then SM SD(A, B) = d( Â, B) + d( B, Â) 2 .

Fig. 3

 3 Fig. 3 Left: Mean Dice coefficient with respect to the number of atlases used. Right: SMSD with respect to the number of atlases used.Majority refers to majority voting after pairwise registration, Local corresponds to appearance-based locally weighted voting, Soft Selection refers to [2], while Hard Selection corresponds to the proposed method.

Fig. 4

 4 Fig.4Sagittal views of the ground truth segmentation and the result produced by the proposed method, majority voting, local appearance-based voting and[START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF]. For the segmentation shown all 17 atlases have been used. Ellipses highlight areas where the proposed method has recovered more details than the competing ones.

Fig. 5

 5 Fig.5Ground truth segmentation and pairs of warped atlas segmentation mask and membership field for the atlas with the least (6% of the total number of control points) and most (11% of the total number of control points) deselected control points. For the segmentation shown all 17 atlases have been used. Note that the areas where segmentation errors occur (highlighted by red ellipses) have been correctly deselected.

Table 1

 1 Mean Dice and SMSD using 17 atlases for all four methods. point x and s the diameter of the patch. We refer to this method in the figures as Local, iii)[START_REF] Alchatzidis | Discrete Multi Atlas Segmentation using Agreement Constraints[END_REF], which is dubbed as Soft Selection because it couples registration and segmentation through the use of appearance-based local soft selection. The proposed method is termed as Hard Selection in the figures.

	Method Dice	SMSD (mm)
	Hard S. 0.789 (0.018) 1.019 (0.114)
	Soft S.	0.786 (0.019)	1.059 (0.105)
	Major.	0.775 (0.018)	1.072 (0.138)
	Local	0.779 (0.017)	1.052 (0.122)

Acknowledgements

This research was partially supported by European Research Council Grant Diocles (ERC-STG-259112).


Innovative IPM for Winter Wheat-based Rotations: First Results of ex post Sustainability Assessment of Cropping Systems Tested at INRA (France

Within the context of the PURE project (WP2), innovative IPM cropping systems were designed for winter wheat-based rotations in the Paris basin area, at INRA in France. We used a three-step prototyping method to design the cropping systems: (1) crop successions and agricultural practices were defined for each system, (2) the prototypes were ex ante assessed and, (3) the most promising systems were tested in field trials and ex post Here, we presented results of the ex post sustainability assessment after one complete rotation.

METHODOLOGY

Three cropping systems were designed according to a gradient of pesticide use intensity: (1) current agricultural practices (C.S.) with a conventional use of pesticides, (2) intermediate level of IPM (I.S.) with a reduction in pesticide use and (3) advanced level of IPM (A.S.) where no pesticides are allowed (for more details, see the poster of Colnenne-David et al., 2015).

Ex post assessment was performed with DEXiPM (Pelzer et al., 2012) on these systems after one complete rotation.

METHODOLOGY

Three cropping systems were designed according to a gradient of pesticide use intensity: (1) current agricultural practices (C.S.) with a conventional use of pesticides, (2) intermediate level of IPM (I.S.) with a reduction in pesticide use and (3) advanced level of IPM (A.S.) where no pesticides are allowed (for more details, see the poster of Colnenne-David et al., 2015).

Ex post sustainability assessment was performed with DEXiPM (Pelzer et al., 2012) on these systems after one complete rotation.

Colnenne-David C. 1,2 , Grandeau G. 1,2 , Tanneau V. 1,2 , Bénézit M. 3 , Angevin F. 3 , Lefèvre L. 1,2 , Doré T. 2,1 1 INRA, UMR 211, 78850 Thiverval-Grignon, France (caroline.colnenne@grignon.inra.fr), 2 AgroParisTech, UMR 211, 78850 Thiverval-Grignon, France, 3 INRA, UAR 1240, 78850 Thiverval-Grignon, France

RESULTS

After one complete rotation, all cropping systems achieve a "high" score (4/5) in terms of overall sustainability. However, result is obtained by very combinations of performances on the three sustainability pillars: economical pillar is higher in the C.S., environmental pillar is higher in the A.S. The social sustainability has remained medium (3/5) for all systems.

There is a clear hierarchy between these systems in terms of the environmental sustainability, which can be explained by various uses of pesticides. TFI are as follows: C.S. = 4.7 < I.S. = 1.8 < A.S. = 0.0 In the C.S., high level of pesticide uses led to a decline in all subcomponents of the environmental sustainability. On the contrary, in the A.S. (i.e. without any pesticide), environmental performances are judged very high.

DISCUSSION -CONCLUSION

Combining innovative strategies in the I.S. and A.S. (i.e. high diversity of species sown in the rotation, the use of resistant variety mixtures, high seed density and delayed sowing dates) led to a decrease of pesticide applications and to an improvement of the environmental performances.

Performance results show t hat it is difficult to meet various objectives. In the C.S., economical performances are high (i.e. high yields, and good produce quality) while environmental performances are judged as medium (i.e. high pesticide use). In the A.S., the performances are the opposite of the C.S. results and are explained by medium yields, with low produce quality, and no pesticide use.

Main results of ex ante and ex post sustainability assessments of the three cropping systems are close to each other. Therefore, DEXiPM seems to be a relevant tool to perform initial assessments required during the innovative cropping system design processes.

The PURE project is supported by the European commission through the seventh framework program (FP7/2011-2014) under contract number 265865