Slim Medimegh

Jean-Yves Pierron

Jean-Pierre Gallois

Frédéric Boulanger
email: frederic.boulanger@lri.fr

A New Approach of Qualitative Simulation for the Validation of Hybrid Systems

Keywords: hybrid systems, qualitative simulation, verification, validation

Hybrid systems are specified in a heterogeneous form, with discrete and continuous parts. Simulating such systems requires precise data and computational power in order to synchronize continuous changes and discrete transitions. However, in the early design stages, the lack of precision about some parameters forbids such simulations. Qualitative simulation consists in computing only some qualities of the behavior, without computing the exact values. We present here a new approach for the qualitative simulation of hybrid systems, which relies on an abstract model of the state variables and of the qualitative aspects of their evolution. This model is analyzed by the Diversity symbolic execution engine to build the qualitative behaviors of the system.

Introduction

Embedded software has become essential in most industrial sectors: energy, transport, telecommunications, healthcare, etc. To ensure a high level of reliability, it is essential to carry out, in the early phases of the development cycle, an analysis of the behavior of the system when it interacts with its environment. This environment usually involves various business knowledge (mechanical, hydraulics, electronics) and it is important to describe each of them in an appropriate formalism. The whole system (the software and its environment) is specified in a heterogeneous manner and contains discrete and continuous parts. The simulation of such hybrid systems requires precise data and computational power in order to detect changes in the continuous values and to synchronize them with discrete transitions. However, during the early stages of the design, the exact value of some parameters is not known yet, while it is already necessary to check the possible behaviors of the system to make some design decisions.

Qualitative simulation can help here. When the behavior of the continuous state variables is described by differential equations, it partitions their domain of variation into discrete domains of qualitative behavior (increasing, decreasing, or constant) according to the sign of their first derivative. Symbolic execution can then produce a tree of abstract behaviors, which can be combined with the model of the discrete part of the system. When the differential equations are not available, qualitative simulation can be performed using a qualitative model which describes only the causal links between the qualitative behavior of the variables and their derivatives, in the form of a state machine.

In this article, we present a new approach of qualitative simulation, which relies on a qualitative model of the differential equations of the continuous part of the system, taking into account the first and second derivatives of the state variables, and on a model of execution implemented in the Diversity tool to compute the qualitative behavior of the hybrid system.

Hybrid Systems

Hybrid systems are dynamic systems that explicitly and simultaneously combine continuous and discrete behaviors. The modeling of hybrid dynamical systems can be made using various formalisms, among which hybrid automata, which are defined by a set of continuous variables and states, and discrete transitions that include guards and assignments to these variables. The evolution of the continuous variables can be described by different kinds of differential equations.

Polyhedral hybrid automata have linear differential equations and guards. Rectangular hybrid automata have constant differential equations and guards that are comparisons to constants. Timed hybrid automata have differential equations and guards involving time [START_REF] Gallois | Qualitative simulation and validation of complex hybrid systems[END_REF], and can be handled by Uppaal [START_REF] Bengtsson | Automated verification of an audio-control protocol using uppaal[END_REF], Kronos [START_REF] Bozga | Kronos: A model-checking tool for real-time systems[END_REF] and recently TiAMo [START_REF] Bouyer | Symbolic optimal reachability in weighted timed automata[END_REF]. Some tools and methods deal with other types of hybrid systems: Hytech [START_REF] Henzinger | Hytech: A model checker for hybrid systems[END_REF] for rectangular hybrid automata with linear guards on the transitions, and CheckMate for non linear continuous evolutions with linear guards [START_REF] Chutinan | Computational techniques for hybrid system verification[END_REF].

Qualitative Simulation

Qualitative simulation is a re- The purpose of qualitative simulation is to develop a representation that makes it possible to reason about the behavior of the physical system without the quantitative information. Figure 1 shows the kind of qualitative changes that we are interested in.

For hybrid systems, the discrete part of the system is not influenced by the qualitative abstraction process. However, continuous variables are discretized in order to consider only their qualitative changes. Therefore, continuous transitions become discrete transitions, and the resulting system is entirely discrete and can be treated by traditional methods (property verification, test generation) for the verification of discrete systems.

Related Work

Many methods have been developed for qualitative reasoning, for instance:

-The component approach in EnVision by Kleer and Brown, builds a qualitative model of the system from its components. It takes into consideration the qualitative sign of the variables. It determines the qualitative states and the transitions between these states [START_REF] Kleer | A qualitative physics based on confluences[END_REF]. -The Qualitative Process Theory developed by K. D. Forbus models objects and processes which cause changes to the objects. Qualities of objects are defined by inequalities about quantities associated to these objects, and the qualitative behavior of a system is obtained by reasoning on the processes [START_REF] Forbus | Qualitative process theory[END_REF]. -Qualitative Simulation, in which the qualitative model of the system is a differential qualitative equation which is an abstraction of an equivalence class of ordinary differential equations. It relies on the continuity properties of the variables to build a sequence of qualitative states describing the behavior of the system [START_REF] Kuipers | Qualitative simulation[END_REF]. Many tools have been developed for qualitative analysis: KAM [START_REF] Yip | KAM : a system for intelligently guiding numerical experimentation by computer[END_REF] which deals with nonlinear system with two freedom degrees, Maps [START_REF] Zhao | Computational dynamics: Modeling and visualizing trajectory flows in phase space[END_REF] which deals with second and third order systems, PSX2NL [START_REF] Nishida | Artificial intelligence research in japan grammatical description of behaviors of ordinary differential equations in two-dimensional phase space[END_REF] which deals with ordinary differential equations, POINCARE [START_REF] Sacks | Automatic analysis of one-parameter planar ordinary differential equations by intelligent numeric simulation[END_REF] which also deals with ordinary differential equations with one parameter. -Qualitative simulation for polynomial continuous evolutions developed by Tiwari [START_REF] Tiwari | Series of abstractions for hybrid automata[END_REF]: the system is modeled as a discrete system whose guards are polynomial, and these guards can be evaluated by the QEPCAD tool.

Our approach is similar to the Tiwari method, but based on qualitative values for avoiding the evaluation of polynomial guards with QEPCAD.

Qualitative Simulation with Diversity

We have implemented our approach, which focuses on first and second derivatives, in Diversity, a symbolic execution engine developed at CEA LIST in order to prove safety properties and to automatically generate tests from state machines models of complex systems [START_REF] Gallois | Le test structurel pour la vérification de spécifications de systèmes industriels: L'outil agatha[END_REF][START_REF] Rapin | Behavioural unfolding of formal specifications based on communicating automata[END_REF]. To analyze complex systems that include continuous parts in hybrid automata, CEA LIST had to extend Diversity with qualitative methods [START_REF] Gallois | Qualitative simulation and validation of complex hybrid systems[END_REF]. For this, the state machine of hybrid models has continuous variables, which must be abstracted by a discretization method.

Discretization process

Each continuous variable is represented by a state machine which encodes its qualitative state: increasing, decreasing, constant etc. The guard of the transitions depends on the sign of the derivative of the variable, which is constrained by the differential equations. These equations are analyzed to determine under which conditions some qualitative changes are possible or not. When the change is not possible, the corresponding execution branch is cut. When the change is possible, the branch is tagged with the condition of the change.

When the structure of the equations is not suitable for processing by a solver, or to avoid the cost of solving the differential equations, it is possible to perform qualitative simulation without differential equations, by using only a qualitative model of these equations. The results are of course less precise, but this approach can be used at very early stages of the design.

From this model, Diversity generates a tree of behaviors, which can be analyzed for reachability properties or cycles. Execution traces can be compared with the result of quantitative simulations, and the compliance of the numerical simulation to the qualitative simulation can be established.

A New Approach for Qualitative Simulation

In order to improve the qualitative simulation without differential equations in Diversity, we developed a model of execution which constraints the qualitative evolution of the state variables, and computes their qualitative behavior.

Model of Execution

In this model of execution, we consider only continuous changes of the state variables and their first derivative. The second derivative, which is the result of input forces, may not be continuous (but this choice is not definitive). For instance, the value of a state variable cannot change from Negative to Positive without being Null in between. Moreover, the value cannot change from Negative to Null unless the first derivative is Positive, and the same is true for the first derivative with regard to the second derivative.

These constraints of continuity and deriv- ability form a model of execution which reflects the physical nature of the phenomenon that we are modeling. It is modeled using state machines as illustrated in Figure 2 for the control of the value of a state variable.

A similar automaton controls the changes of the first derivative with respect to the second derivative. These state machines are nondeterministic. For instance, in the < 0 state, if the first derivative is Positive, we can either stay in the current state, or go to the = 0 state.

Computing Qualitative Behaviors

In order to compute the qualitative behavior of a state variable, we observe the changes of the values attached to this variable. For instance, when the second derivative is Negative and the first derivative is Null, we have a Maximum. However, the current value of a state variable and the current value of its first and second derivatives are not enough to make a difference between a maximum and a start of decrease after a constant stage. We therefore model each continuous state variable in the system by five values: the current (x) and previous (x -1 ) values of the variable, the current ( ẋ) and previous ( ẋ-1 ) values of its first derivative, and its second derivative (ẍ).

Taking into account only ẋ-1 , ẋ and ẍ, we identified 13 qualitative states, which cover our analysis goals, without the need for higher order derivatives:

Constant when ẋ-1 , ẋ and ẍ are null, FlexStartIncrease when ẋ-1 = 0, ẋ = 0 and ẍ > 0. The first derivative is not positive yet, but the dynamics is transitioning toward an increase. FlexStartDecrease when ẋ-1 = 0, ẋ = 0 and ẍ < 0. Similar to the previous case when transitioning toward a decrease. StartIncrease when ẋ-1 = 0 and ẋ > 0. The first derivative has just become positive. Notice that there is no condition on the second derivative, which may have come back to 0. StartDecrease when ẋ-1 = 0 and ẋ < 0. Increase when ẋ-1 > 0 and ẋ > 0. The first derivative is positive. Decrease when ẋ-1 < 0 and ẋ < 0. The first derivative is negative. Maximum when ẋ-1 > 0, ẋ = 0 and ẍ < 0. The three conditions are necessary to distinguish this case from other qualitative states such as inflection points. Minimum when ẋ-1 < 0, ẋ = 0 and ẍ > 0. The three conditions are necessary to differentiate this case from other qualitative states such as inflection points. FlexIncrease when ẋ-1 > 0, ẋ = 0 and ẍ > 0. Inflection point during an increase. FlexDecrease when ẋ-1 < 0, ẋ = 0 and ẍ < 0. Inflection point during a decrease. StopIncrease when ẋ-1 > 0, ẋ = 0 and ẍ = 0. The state variable reaches a plateau at the end of an increase. StopDecrease when ẋ-1 < 0, ẋ = 0 and ẍ = 0. The state variable reaches a plateau at the end of a decrease.

By additionally taking into account x -1 and x, it is possible to distinguish sub-cases in these qualitative states, for instance "reaching a null maximum" when reaching a maximum with x -1 < 0 and x = 0.

There are impossible cases. Two are due to the continuity of the variable: the conjunction of x -1 < 0 and x > 0, and the conjunction of x -1 > 0 and x < 0 are impossible. Two others are due to the continuity of the first derivative (no angular points). Ten others are due to the fact that each variable change should be consistent with the previous value of the first derivative. For instance, with x -1 = 0 and ẋ-1 = 0, we necessarily have x = 0, because starting from 0 with a null first derivative, the variable cannot change. Other cases are similar.

Illustrative Example

This technique was applied to the simple example of a bouncing ball. In this example, we have only one state variable, the height of the ball above the ground, noted z. The usual way to model the bouncing ball is to consider that it has only one state, in which it is in free fall, with z = -g and g = 9.81m.s -2 . The bounce is modeled with a single transition, which is triggered when the ball hits the ground (z = 0), and reverses the speed of the ball with a damping factor 0 < c ≤ 1. This model and its qualitative abstraction are shown in Figure 3. The quantitative behavior of this model, as obtained in the Ptolemy II environment [START_REF] Eker | Taming heterogeneity -the ptolemy approach[END_REF], is illustrated in Figure 4. This kind of model is an abstraction of what really happens, because if the ball were really changing its speed instantaneously, there would be an exchange of a finite amount of energy in zero time between the ball and its environment, which corresponds to an infinite power.

z = -g z = 0 / ż = -c. ż z > 0 ż = 0 z < 0 z = 0 / ż = - ż z > 0 ż = 0

Qualitative Simulation of the Raw Hybrid Model

When performing the qualitative simulation of this unphysical model using our model of execution, Diversity finds a dead-lock: This dead-lock happens because the model of the ball makes the first derivative change from negative to positive without going through zero, locking the state machine which models the first derivative in its "negative" state.

Flex Start Decrease z ż-1 ż z-1 z < 0 = 0 = 0 > 0 > 0 Start Decrease z ż-1 ż z-1 z < 0 = 0 < 0 > 0 > 0 Decrease and bounce z ż-1 ż z-1 z < 0 < 0 > 0 > 0 = 0 Dead-lock z ż-1 ż z-1 z < 0 < 0 > 0 = 0 < 0
Our model of execution enforces the continuity and derivability of physical states, but the model of the ball violates these properties. However, we have to handle such models because that is the way engineers model systems, and quantitative simulation tools are able to handle them as shown in Figure 4.

Adjusting the Execution Model

A first solution is to adjust our execution model in order to allow such "unphysical" transitions. With this model of execution, the qualitative state changes to "Increase" after the bounce, then "Maximum", then to "Start Decrease" before cycling through the "Decrease" state, as shown in Figure 6. Therefore, the ball was able to bounce, and the maximum of its height was detected, but not the minimum that was reached during the bounce. The issue here is that the detection of the qualitative states relies on the continuity and derivability of the state variable. The discontinuous change of the velocity of the ball prevents the detection of the minimum.

Flex Start Decrease z ż-1 ż z-1 z < 0 = 0 = 0 > 0 > 0 Start Decrease z ż-1 ż z-1 z < 0 = 0 < 0 > 0 > 0 Decrease (and bounce) z ż-1 ż z-1 z < 0 < 0 > 0 > 0 = 0 Increase z ż-1 ż z-1 z < 0 > 0 > 0 = 0 > 0 Maximum z ż-1 ż z-1 z < 0 > 0 = 0 > 0 > 0 Start Decrease z ż-1 ż z-1 z < 0 = 0 < 0 > 0 > 0

Adjusting the model

Another solution is to adjust the model of the ball. However, we said that we must be able to handle models as engineers create them. In the case of the bouncing ball, the model can be automatically adjusted by replacing the bouncing transition by a series of transitions. The algorithm for this is as follows:

z < 0 z > 0 z > 0 z > 0 ż = 0 z = 0 ż = 0 ż > 0
changing the derivative from < 0 to > 0 is illegal. The only legal path is to go from < 0 to = 0 and then from = 0 to > 0, so we replace the illegal transition by the sequence of these two transitions. changing the derivative from < 0 to 0 and from 0 to > 0 requires that the second derivative be positive, so we make the second derivative positive during the bounce (the second derivative can be discontinuous).

The resulting state machine is shown in Figure 7, where the additional states have a gray background. The first additional state, reached when z = 0, sets the second derivative to "Positive" because this is required to make the first derivative change from "Negative" to "Null". When the first derivative becomes null, the second additional state is reached, and when it becomes positive, as requested by the initial un-physical transition, we can go back to the free fall state of the ball, with a reversed velocity.

With this adjusted model of the bouncing ball and the physical model of qualitative simulation, we obtain the qualitative behavior shown in Figure 8. The qualitative simulation now correctly detects a minimum when the ball hits the ground and bounces, and a maximum when the ball reaches its highest position after having bounced. However, the original model was adjusted to match the model of computation of the qualitative behavior. How is the qualitative behavior we obtained related to the behavior of the original model?

Flex Start Decrease z ż-1 ż z-1 z < 0 = 0 = 0 > 0 > 0 Start Decrease z ż-1 ż z-1 z < 0 = 0 < 0 > 0 > 0 Decrease (hit ground) z ż-1 ż z-1 z > 0 < 0 < 0 > 0 = 0 Minimum (bounce) z ż-1 ż z-1 z > 0 < 0 = 0 = 0 < 0 Start Increase z ż-1 ż z-1 z < 0 = 0 > 0 = 0 < 0 Increase z ż-1 ż z-1 z < 0 > 0 > 0 < 0 = 0 Maximum z ż-1 ż z-1 z < 0 > 0 = 0 = 0 > 0 Start Decrease z ż-1 ż z-1 z < 0 = 0 < 0 > 0 > 0

Discussion about the Qualitative Model of Execution

In this qualitative simulation, we observe an unexpected artifact: the height of the ball above the ground becomes negative during the bounce. This comes from the fact that in the adjusted model, the first derivative takes some time to become null when the ball hits the ground, so the height goes from null to negative. Since we adjusted the model to make it more physical, we have to check if this artifact has some physical meaning.

The condition z = 0 for triggering the bounce corresponds to the surface of the ball hitting the ground. Therefore, the z state variable represents the height of the center of the ball above the ground. z becomes null at the moment when the ground starts pushing the ball upwards. Therefore, the z > 0 action we added to adjust the model represents the result of this pushing. This makes the ball slow down and stop ( ż goes from negative to null), but since ż is still negative, z is still decreasing, so it goes from null to negative. For the physical system, this corresponds to the deformation of the ball, which converts kinetic energy into elastic tension in the ball (and heat if the collision is not elastic). When the ball is deformed against the ground, its center is below the position where it stood when the ball hit the ground, so z is negative.

We consider that the states and transitions added to adjust the model of the ball to the model of computation of the qualitative behavior are only reconstructing the physical behavior of the system that was abstracted in the original model. In the example of the bouncing ball, this reconstruction was possible in a single and systematic manner. However, this preliminary work did not allow us to verify that such adjustments are unique, or that there exists a minimal one for more complex systems.

When illustrating our approach with the example of the bouncing ball, we have only shown one of the qualitative behavior found by the Diversity tool. Among the other behaviors, some are physically possible. For instance, after bouncing on the ground, the ball can go up indefinitely, without reaching a maximum. This is possible if the ball reaches the escape velocity, but since we ignore the exact value of the velocity of the ball, we cannot tell, using qualitative simulation, whether the ball will reach a maximum or fly upwards forever.

Another possible behavior is more troubling: the ball can fall forever towards the ground without reaching it. This is indeed possible in our qualitative model of computation because it is only first order. We qualify only the variations of the state variable, not the variations of its derivative. Therefore, when the ball is falling, we only know that its height is decreasing, not that its velocity is increasing in magnitude. Without this information, it is possible that the ball goes slower and slower towards the ground, never reaching it, and the good thing is that the tool finds this behavior. We therefore plan to build a more complex, second order model of computation which will eliminate such physically impossible behaviors from the qualitative simulation.

Conclusion and future work

In this work, we have established a qualitative model of computation based on the first and the second derivative of the state variables to describe the behavior of hybrid systems without differential equations (with only a qualitative model of the equations). Our approach has been applied to the typical hybrid system example of the bouncing ball. We have shown that we can find a qualitative path that allows us to describe the whole behavior of the ball, respecting the continuity and derivability of the first and the second derivative. Using only the first two derivatives is enough to obtain the qualitative behavior of the variables in terms of variation, minimum and maximum. The tree of qualitative behaviors generated by Diversity can be used to prove reachability properties or to detect cycling behaviors.

To continue this work, we plan to enrich the model of computation with second order qualities (qualitative variations of the first derivative) in order to automatically eliminate some unphysical behaviors from the simulation results. Another track to explore is dealing with several bounded state variables, with constraints on their qualitative values or the qualitative value of their derivatives.

The main goal of this approach is to allow the qualitative analysis of hybrid models during the early stages of the design of a system. The abstraction process used to transform an hybrid automata into a qualitative discrete model, and the automatic transformation of this model into one which does obey the constraints of "physicality" encoded in our model of execution, are still very experimental. They can be considered as a kind of abstraction glue to combine the continuous and the discrete models, and its semantics needs a more precise definition to improve the confidence in the results of the qualitative simulation.

Fig. 2 .

 2 Fig. 2. Qualitative changes

Fig. 3 .

 3 Fig. 3. Hybrid automaton (left) and qualitative (right) automaton of the bouncing ball

Fig. 4 .

 4 Fig. 4. Quantitative behavior of the bouncing ball

Fig. 5 .

 5 Fig. 5. Qualitative simulation of the raw model

Fig. 6 .

 6 Fig. 6. Qualitative simulation with discontinuous transitions

Fig. 7 .

 7 Fig. 7. Adjusted Model for the Bouncing Ball

Fig. 8 .

 8 Fig. 8. Qualitative simulation with adjusted model


Jose E Romero

Pierrick Coupé

Rémi Giraud

Vinh--Thong Ta

Vladimir S Fonov

Mallar Min Tae Park

M Mallar Chakravarty

Aristotle N Voineskos

Jose V Manjón

Min Tae

CERES: A new cerebellum lobule segmentation method

Keywords: MRI, cerebellum lobule segmentation, non--local multi--atlas patch--based label fusion, optimized patchmatch

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The human cerebellum is a neuroanatomical structure within the human brain located below the cerebrum and connected to the brainstem through the cerebellar peduncles. Although it represents a small percentage of the total intracranial volume, about 10%, it plays a key role in motor coordination and learning [START_REF] Manto | Handbook of the Cerebellum and Cerebellar Disorders[END_REF]. The cerebellum acts as a modulator for the motor cortex signals to perform fine and precise movements through the integration of the inputs signals from the somatosensory system and cerebellar cortex. It has been shown that cerebellar damage leads to important motor dysfunctions such as limb ataxia, balance alterations and other deficits in motor coordination [START_REF] Kase | Cerebellar infarction. Clinical and anatomic observations in 66 cases[END_REF][START_REF] Compston | Multiple sclerosis[END_REF][START_REF] Davie | Persistent functional deficit in multiple sclerosis and autosomal dominant cerebellar ataxia is associated with axon loss[END_REF][START_REF] Klockgether | The clinical diagnosis of autosomal dominant[END_REF]. Recent neuroimaging studies indicate that cerebellum is also involved in many cognitive functions [START_REF] Stoodley | The cerebellum and cognition: evidence from functional imaging studies[END_REF] such as attention, emotion or behavior [START_REF] Timmann | The human cerebellum contributes to motor, emotional and cognitive associative learning[END_REF].

Despite the cerebellum being involved in many brain functions, it has been widely understudied.

Lately, several publications have focused on the study of the cerebellum and its implication in several neuropsychiatric pathologies. For example, patients with schizophrenia showed a reduced vermis [START_REF] Okugawa | Selective reduction of the posterior superior vermis in men with chronic schizophrenia[END_REF] and subjects with Alzheimer's disease presented smaller posterior lobes [START_REF] Thomann | The cerebellum in mild cognitive impairment and Alzheimer's disease --a structural MRI study[END_REF]. Nevertheless, there exist publications with inconclusive results regarding the role of the cerebellum in different neuropsychiatric disorders [START_REF] Nenadic | Distinct pattern of brain structural deficits in subsyndromes of schizophrenia delineated by psychopathology[END_REF][START_REF] James | Cerebellar, prefrontal cortex, and thalamic volumes over two time points in adolescent--onset schizophrenia[END_REF]. This could be explained by the fact that accurately segmenting the cerebellum is quite challenging due to its complex anatomy.

Cerebellar anatomy consists of a white matter tree structure located behind the pons. It is divided into two hemispheres (left and right), with each white matter branch surrounded by a layer of grey matter that creates folds called foliations. These grey matter folds are denominated cerebellum lobules. The size, position and number of (visible) lobules is highly variable between subjects which makes the segmentation process challenging.

Manual segmentation (delineation) by an expert requires a huge effort since it is a time consuming and tedious process. Delineation by non--expert raters may be faster but still impractical for studies applied to a large cohort of subjects. Further, both alternatives require dealing with inter and intra--rater variability. A few semi--automatic methods have been developed offering large structure parcellation as the brainstem, white matter and grey matter [START_REF] Weier | Evaluation of a new approach for semi--automatic segmentation of the cerebellum in patients with multiple sclerosis[END_REF]. These methods intend to reduce the time required by a human rater while keeping the reliability of manual delineation. However they are still time consuming.

An alternative to the manual segmentation is using multi--atlas based segmentation [START_REF] Aljabar | Multi--atlas based segmentation of brain images: Atlas selection and its effect on accuracy[END_REF] which consists of transferring and combining multiple manual segmentations to the case to be segmented. One of the first automated algorithms specifically developed for cerebellum segmentation is SUIT [START_REF] Diedrichsen | A spatially unbiased atlas template of the human cerebellum[END_REF]. In this method, a set of manually labeled images was used to build a probabilistic atlas that was then warped into the target image space. However, SUIT is not able to capture all the anatomical variability due to the use of a single atlas approach.

In addition, the method developed by Diedrichsen is semi--automatic since it includes a cerebellum extraction step that is suggested to be manually corrected for optimal results. Fully automatic methods appeared in the last years. [START_REF] Bogovic | Automated segmentation of the cerebellar lobules using boundary specific classification and evolution[END_REF] published a method called ACCLAIM using deformable models and random forest boundary classification. Recently Yang et al. (Yang et al., 2015) proposed a segmentation method combining multi--atlas label fusion and tissue/boundary classification in a graph--cut segmentation framework to address the high variability of the cerebellum anatomy achieving good results at the expense of a high computational burden (execution time of several hours per case). Also recently, a new approach based on multi--atlas segmentation called MAGeT Brain [START_REF] Chakravarty | Performing label--fusion--based segmentation using multiple automatically generated templates[END_REF][START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF] has been proposed for cerebellum lobule segmentation. This new approach aims at solving the problem of having very few exemplar cases (manual segmentations) of a given segmentation protocol. It uses the estimation of a large set of non--linear deformations to map the example cases to a set of unlabeled templates increasing the number of candidate segmentation to apply a final majority vote based label fusion. MAGeT Brain demonstrated better results than single atlas and classical multi--atlas label fusion using majority voting methods at the expense of a high computational cost in terms of both memory and execution time. Finally, Weier et al. published another multi--atlas segmentation method named RASCAL [START_REF] Weier | Rapid Automatic Segmentation of the Human Cerebellum and its Lobules (RASCAL) --Implementation and Application of the Patch--based Label--fusion Technique With a Template Library to Segment the Human Cerebellum[END_REF] based on a multi-atlas patch--based label fusion. It utilizes a library of manually segmented cases that is non--linearly registered into the target image space to perform the labeling process.

In this paper, we introduce a new pipeline for cerebellum lobule segmentation that is based on an adaptation of a recently proposed segmentation method called Optimized PatchMatch Label fusion (OPAL) [START_REF] Ta | Optimized patchMatch for near real time and accurate label fusion[END_REF][START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF]. This method consists of a multi--atlas patch-based [START_REF] Rousseau | A supervised patch--based approach for human brain labeling[END_REF] segmentation with a non--local label fusion technique that produces fast and accurate segmentations using a library of manually segmented cases. The original method developed for hippocampus segmentation (Coupe et al., 2011) has been adapted to segment the cerebellum anatomy using a library of non--linearly registered cases instead of using only an affine registration as in the original method. A post--processing step has been added to enforce regularity of the different lobule labels. Also, given the small size of the template library (only 5 cases) we extended the number of templates of the library by adding automatically segmented cases.

Materials and methods

In the following sections, the data used in this paper and the proposed method details are described. We have called the proposed method CERES (for CEREbellum Segmentation).

Image data

In this work, two datasets have been used. A High Resolution (HR) set of 5 MR T1--weighted images has been used to construct the library of manually labeled cases. A second dataset of 20 standard resolution images has been used to evaluate the proposed method and to compare it with other related approaches. The HR images are publicly available at the CoBrALab website (http://cobralab.ca/atlases/Cerebellum.html). Both the HR images used as input and the manually labeled validation dataset are the same as those used in [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF].

HR dataset

These cerebellar atlases were created based on HR MR images from 5 healthy volunteers (2 male, 3 female, aged 29-57). We used T1--weighted images that were acquired on a GE Discovery MR 750 3 T system (General Electric, Waukesha, WI) using an 8--channel head coil. High--resolution T1-weighted images were acquired using the 3D inversion--prepared fast spoiled gradient--recalled echo acquisition, FSPGR--BRAVO, in a scan time of ~ 20 min with the following parameters: TE/TR = 4.3 ms/9.2 ms, TI = 650 ms, α = 8°, 2NEX, FOV = 22 cm, slice thickness = 0.6 mm, 384 × 384 in-plane steps for an approximate isotropic resolution of 0.6 mm dimension voxels. Reconstruction filters, ZIPX2 and ZIP512, were also used resulting in isotropic 0.3 mm 3 voxels with a final matrix size of 489x734x503 voxels [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF]. These HR MR images were manually segmented by two expert raters that allowed for both inter--and intra--rater comparisons of segmentations to validate the consistency of the manual segmentation protocol (for details of the segmentation protocol please see [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF]). The cerebellum was manually delineated into 26 structures: White matter and Lobules I--II, III, IV, V, VI, Crus I, Crus II, VIIb, VIIIa, VIIIb, IX, and X considering left and right hemispheres as shown in Figure 1. One of the authors of this manuscript, Min Tae M. Park, was one of the manual raters of these data.

Standard resolution dataset

MR images from 10 healthy control and 10 patients with schizophrenia (15 male, 5 female, average 30.1 years age, range 25-35 years) were used to create a validation set. Imaging data were acquired using an 8--channel head coil on a 1.5 T GE Echospeed system (General Electric Medical Systems, Milwaukee, WI), which permits maximum gradient amplitudes of 40 mT/m. Axial inversion recovery prepared spoiled gradient recall images were acquired: echo time (TE) = 5.3 ms, repetition time (TR) = 12.3 ms, time to inversion (TI) = 300 ms, flip angle = 20°, number of excitations (NEX) = 1 (124 contiguous images, 1.5 mm thickness, 0.78 x 0.78 mm in plane resolution). These images were linearly registered to the MNI space using MRITotal software [START_REF] Collins | Automatic 3D Inter--Subject Registration of MR Volumetric Data in Standardized Talairach Space[END_REF] to perform the manual labeling thus resulting in a 1 mm 3 image resolution and a matrix size of 181x217x181 voxels. The whole cerebellum is manually segmented on these 20 cases dataset and the left cerebellum was segmented on 10 of the cases while in the other 10 the right cerebellum was segmented (24 lobules: I--II, III, IV, V, VI, Crus I, Crus II, VIIb, VIIIa, VIIIb, IX, and X). The definition of the lobules is the same as the atlases and it is shown in Figure 1. One of the authors of this manuscript, Min Tae M. Park, was one of the manual raters of these data.

Preprocessing

All the images are preprocessed to enhance their quality and to transform them to a common intensity and geometric space. The steps included in this preprocessing are the following:

Denoising

The Spatially Adaptive Non--local Means Filter is applied to reduce the noise in the images. This filter was chosen because it is able to automatically deal with both stationary and spatially varying noise levels (Manjón et al., 2010).

Inhomogeneity correction in native space

The N4 bias field correction [START_REF] Tustison | N4ITK: improved N3 bias correction[END_REF] is applied to correct intensity inhomogeneities across the images. This step is applied to improve the registration process but not used as a final bias correction method as will be described later.

MNI affine registration

All the images have to be registered to a common coordinate space so that the anatomy within the case to be segmented and the library cases is consistent. To this end, the images are linearly registered (affine transform) to the Montreal Neurological Institute (MNI) space using the MNI152 template. This is done using the Advanced Normalization Tools (ANTs) [START_REF] Avants | Advanced normalization tools (ANTS)[END_REF].

Inhomogeneity correction in MNI space

N4 method is used to remove the inhomogeneity of the images in the MNI space [START_REF] Tustison | N4ITK: improved N3 bias correction[END_REF]. We use the MNI125 intracranial mask to better correct the brain region. We found experimentally that the application of N4 method in MNI space has a superior performance compared to its application on native space.

Cerebellum area cropping

A cropping step is applied to reduce the computational burden of the method limiting the processing to the cerebellum area. A cropping 3D window was calculated from the limits of the average of all the atlas labels of the library in the MNI space (including a margin of 5 voxels in all 3 dimensions).

Low dimensional non--linear registration estimation

To achieve a better cerebellum anatomic matching between the target image and the library of manually labeled templates, the cropped area of both (target and templates) is non--linearly registered to the cropped MNI152 atlas (the non--linear transformation is estimated but not applied). Although the segmentation process could be performed at this non--linear MNI152 space, the inverse transformation needed to map back the segmentation result on the original linear MNI space would reduce segmentation quality due to the introduction of interpolation artifacts.

Applying successive transformations implies performing successive interpolations. Instead, we segment the target case in the linear MNI space avoiding the use of a non--linear inverse transformation on the result. To this end, a customized subject--specific library is created by concatenating the forward non--linear transformation to the MNI152 atlas of each library case and the inverse non--linear transformation of the subject to be segmented. By doing this we are applying the transformation once, so only one interpolation is needed. Although the custom library has interpolation errors these are minimized during the label fusion step.

The registration step has a huge impact on segmentation accuracy. We tested 3 different resolution levels (low=4 mm 3 , medium=2 mm 3 and high=1 mm 3 ). In general, a high dimensional transformation (1 mm 3 ) is expected to better match the reference anatomy. However, this is highly dependent on the data topology and we found that the low dimensional registration (4 mm 3 ) provided the best segmentation results (being also the most efficient in terms of computational load). The non--linear registration was performed using ANTs tool [START_REF] Avants | Advanced normalization tools (ANTS)[END_REF] with cross correlation similarity measure.

Intensity normalization

As previously commented, it is necessary to normalize the images in order to obtain the same intensity values across subjects for cerebellar tissues (white matter and grey matter). This is required to perform the patch--based segmentation that finds matches between the target image and the library of templates using an intensity--based similarity metric. To this purpose, a local intensity normalization based on regions of interest (ROIs) is applied in a similar manner to the method proposed by [START_REF] Asman | Muti--Atlas Segmentation using Non--Local STAPLE[END_REF] in the multi--atlas segmentation challenge of MICCAI 2012.

We use the average MNI152 atlas as the reference image that is mapped to the target space by applying the previously computed non--linear inverse transformation. Then the majority voting segmentation obtained from the library is used to define several ROIs (one for each label). Next, the median intensity for each ROI is calculated for both target and atlas and a second order polynomial function is fitted to map the corresponding median intensities. Finally, the intensities of the target case are normalized to the atlas intensities by evaluating the fitted function at every target voxel intensity value. Note that we use median intensity values instead of mean values because it was found to be more robust to the presence of outliers.

Library construction

To construct the library of templates, the HR dataset mentioned in Section 2.1 was preprocessed as described in Section 2.2. Label maps delineated on native space were linearly mapped to the MNI space. Then, images and labels were left--right flipped to double the library size (from 5 to 10 labeled cases). The fact that the cerebellum is nearly symmetric allows the segmentation process to take benefit from this step as long as we are increasing the number of anatomical references (patterns) locally. Images and labels were then cropped to the cerebellum region. Afterwards, the non--linear transformation of these cropped images to the cropped MNI atlas was estimated. These transformations are used to non--linearly register the entire library to the space of the case to be segmented in a similar way as done in [START_REF] Weier | Rapid Automatic Segmentation of the Human Cerebellum and its Lobules (RASCAL) --Implementation and Application of the Patch--based Label--fusion Technique With a Template Library to Segment the Human Cerebellum[END_REF]. Finally, the images were intensity normalized. Note that the resulting library at MNI space has a 1 mm 3 resolution (original HR data had a 0.3 mm 3 voxel dimensions).

Segmentation pipeline description

Preprocessing

Any new unlabeled case is preprocessed as previously described by applying the denoising, the linear registration to the MNI space, the inhomogeneity correction, the cropping, the non--linear registration and the local intensity normalization. Finally, the estimated non--linear transformation is used to create the subject--specific library. After this process, all the library templates are non-linearly registered to the target case requiring only one non--linear transformation estimation instead of N non--linear transformations (where N is the number of subjects in the library).

Labeling

Non--local patch--based label fusion

Non--local patch--based label fusion is a multi--atlas segmentation technique that combines multiple reference atlas segmentations to produce a new segmentation that fits the anatomy of a new unlabeled case [START_REF] Coupé | Patch--based Segmentation using Expert Priors: Application to Hippocampus and Ventricle Segmentation[END_REF]. In this strategy, the final label for each image voxel is estimated by performing a weighted label fusion of a set of sample patches within a specific search area centered on the target voxel for all N cases in the library. The estimation is computed as follows:

𝑣 𝑥 ! = 𝑤(𝑥 ! , 𝑥 !,! )𝑦 !,! !∈! ! ! !!! 𝑤(𝑥 ! , 𝑥 !,! ) !∈! ! ! !!! (1)
where V i represents the search area, N is the number of subjects in the template library, y s,j is a possible label from the voxel x s,j at the position j in the subject s and w(x i , x s,j ) is the patch similarity defined as:

𝑤 𝑥 ! , 𝑥 !,! = 𝑒𝑥𝑝 ! ! ! ! !!(! !,! ) ! ! ! !
(2)

where P(x i ) is the patch centered at x i , P(x s,j ) the patch centered at x j in the templates and ||.|| 2 is the normalized L2 norm (normalized by the number of elements) calculated from the distance between each pair of voxels from both patches P(x i ) and P(x s,j ). h is a normalization parameter that is estimated from the minimum of all patch distances within the search area.

Although the non--local label fusion is a powerful segmentation technique, one of its drawbacks is its high computational cost (mainly due to the patch similarity estimation) that limits its use on large images and big libraries. Recently, a fast patch matching techniques have been proposed that significantly reduces the computational burden of non--local methods [START_REF] Ta | Optimized patchMatch for near real time and accurate label fusion[END_REF][START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF].

The PatchMatch algorithm

The PatchMatch (PM) algorithm [START_REF] Barnes | PatchMatch: a randomized correspondence algorithm for structural image editing[END_REF] is an efficient strategy to find patch correspondences between two images (originally 2D). It is based on the hypothesis that given a patch p(i,j) ∈ an image A that matches a patch p(i',j') ∈ an image B, then the adjacent patches of p(i,j) will probably match the adjacent patches of p(i',j').

This algorithm is implemented in three steps: initialization, propagation and random search. First, in the initialization, a random patch from B is assigned to each patch in A. Second, adjacent patches are tested to improve the current match by following the concept of adjacency. Third, with the aim of avoiding local minima, random samples are selected within a neighborhood and compared to the current match. The second and third steps are repeated iteratively to improve the correspondences. This method quickly converges and gives an approximate nearest neighbor field (ANN) that is a set of correspondences between images A and B.

Optimized PatchMatch (OPAL)

A 3D extension of PM to deal with MR images was recently proposed by [START_REF] Ta | Optimized patchMatch for near real time and accurate label fusion[END_REF]. This method, called OPAL, applies PM to 3D MR images and finds patch correspondences between a given image S and a library of L templates. One of the benefits of the PM algorithm is that its complexity only depends on the size of the image A and not on the size of the library L. This advantage allows computation of patch correspondences between an image S and a large library L = {T 1 ,...,T N } without using any template preselection at an almost constant computational cost (discarding library loading time). Compared to PM method, OPAL has a constrained random initialization since typically both target and library are located in the same geometric space (MNI space) and therefore good matches for the target patches are expected to be close in the library space. For further details, see the original OPAL paper [START_REF] Ta | Optimized patchMatch for near real time and accurate label fusion[END_REF] and its recent extension paper [START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF].

Multi--scale OPAL label fusion

In patch--based segmentation, the patch size is a key parameter that is strongly related to the structure to be segmented and the image resolution. This can be empirically optimized for a specific dataset. It can be seen in the literature that multi--scale approaches improve segmentation results [START_REF] Manjón | Nonlocal Intracranial Cavity Extraction[END_REF]. In OPAL [START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF] independent and simultaneous multi-scale and multi--feature ANN fields were computed. In this paper, a multi--scale approach was also used but only voxel intensities are considered as a feature. Finally, several different ANNs are computed for different scales resulting in different label probability maps and have to be combined. In OPAL, the final estimation is obtained by simply averaging these maps by using a late fusion [START_REF] Snoek | Early versus late fusion in semantic video analysis[END_REF]. All maps were equally contributing to the final estimation. In this paper, two scales are used and a different label dependent adaptive weighting scheme is proposed to fuse these maps (eq. 3).

𝑝(𝑙)

= 𝛼 𝑙 𝑝 ! (𝑙) + (1 -𝛼(𝑙))𝑝 ! (𝑙) (3)
where p 1 (l) is the probability map of scale 1 for label l, p 2 (l) the probability map of scale 2 for label l, p(l) is the final probability map for label l and 𝛼(𝑙) ∈ [0,1] is the mixing coefficient for label l optimized using a gradient descent technique. For this optimization, we minimized the segmentation error (i.e., maximizing DICE coefficient). For this process, since the training dataset was very small (N=5) we used the validation dataset using a leave--one--out approach to avoid bias in the method comparison section. Specifically, we performed 20 estimation rounds (removing one case each time). Details are given in the experiments section.

Label regularization

As a consequence of the reduced size of the template library, the segmentations produced by the labeling method may present slightly irregular edges. To reduce this effect we propose a post-processing step to regularize the label maps by convolving them with a smoothing kernel of 5x5x5 voxels. Then, the label with the maximum weight is assigned to each voxel. This penalizes narrow structures being coherent with the layered definition of the cerebellum lobules. To prevent small lobules from being eliminated, the original (not convolved) maps are combined with the corresponding smoothed one using a mixing parameterβ:

𝑝 ! 𝑙 = 𝑝 𝑙 + 𝑙 • 𝐻 𝑝(𝑙) (4)
where 𝑝 ! (𝑙) represents the final weight of label 𝑙, H is the convolution kernel and β(𝑙) is the regularization mixing coefficient of label l optimized in an orthogonal manner performing a dichotomic search to find the optimal value for each label. We performed the optimization using the validation dataset following a leave--one--out approach as done for the estimation of the alpha coefficients.

Figure 2 presents an overview of the whole pipeline showing the workflow for a new case to be segmented.

Library extension

To mitigate the effect of the reduced size of the library of manually segmented templates, an automatic library extension strategy was used. To this purpose, 100 randomly selected unlabeled images from healthy subjects were used. The images belong to the IXI --Information eXtraction from Images (EPSRC GR/S21533/02) dataset which is publicly available at its website (http://biomedic.doc.ic.ac.uk/brain--development).

The library extension strategy consisted of applying the described segmentation pipeline to the new set of images and to add these images and the corresponding segmentations to the library as done in some previous works [START_REF] Eskildsen | BEaST: brain extraction based on nonlocal segmentation technique[END_REF][START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF].

Experiments and results

Several experiments were performed to estimate optimal method parameters and to compare the performance of the proposed method in terms of accuracy and computational time with other related state--of--the--art approaches.

Method parameters

As OPAL [START_REF] Giraud | An Optimized PatchMatch for multi--scale and multi--feature label fusion[END_REF] has been developed and evaluated for hippocampus segmentation, an optimization was required to adapt a set of parameters in order to obtain optimal cerebellum segmentation results. To measure the accuracy we used the DICE [START_REF] Zijdenbos | Morphometric analysis of white matter lesions in MR images: method and validation[END_REF] coefficient for the 13 structures of the cerebellum (left and right labels have been combined into a single one). Results for the following experiments are calculated over the validation set. Note that since only half cerebellum was manually labeled, results were obtained only comparing the segmentation results with the labeled area (discarding the other half cerebellum). In all the experiments we set the patch sizes to 3x3x3 and 9x9x9 voxels, respectively. The restricted search area was set to 7x7x7 voxels. The number of PMs was set to 32 and the number of iterations of OPAL to 4.

Library construction

We mentioned that CERES works using a preprocessed library of manually segmented cases. In this preprocess the images are located in the MNI space by linear registration. Then the images are cropped and registered again to better fit the cerebellum anatomy. This second registration has a significant impact in the segmentation. We tested three different approaches to estimate the impact of the different transformation options. The first one consists on performing a linear registration in the cropped MNI space. In the second approach a non--linear transformation is estimated to register the training templates to the cropped MNI space. Finally, in the third one, we generate a subject specific library by nonlinearly registering all the training templates to the cropped image of the subject in the MNI space. In table 1, the results of these three different options are presented.

As can be noted, the subject specific library construction option was the best one for all the structures. In this experiment, we applied OPAL method without any of the proposed optimizations to measure just the influence of the library construction approach.

Muli--scale mixing coefficient

As explained before, OPAL obtains the matches for every image patch at different scales and then performs a late fusion. We compared the uniform averaging scheme (α=0.5) as done in original OPAL paper with the proposed adaptive method. Table 2 shows the results for the two compared options. Note that we optimized these coefficients 20 times using leave--one--out cross validation.

Results in table 2 show the average dice of the 20 cases using its corresponding set of coefficients (estimated from the 19 remaining cases of the validation set, i.e. excluding the case under study). In Tab. 3, for sake of simplicity we only show the average of the 20 different estimations. In following experiments, the 20 different sets of alpha coefficients will be used for validation to prevent the validation cases to be included in its own segmentation process.

Impact of label regularization

We performed a similar experiment to measure the effect of the label regularization on the segmentation results. Table 4 shows the results of the proposed regularization method computed over the 20 subjects within the validation dataset. Note that we optimized these coefficients 20 times using leave--one--out cross validation as previously described and the final DICE coefficients shown in table 4 were computed using the corresponding parameter set from the leave--one--out estimation. Also, in the following experiments, the 20 different sets of beta coefficients will be used to prevent the validation cases to be included in its own segmentation process.

Despite the fact that the global improvement is limited, an important increase in lobule I--II, which is a difficult structure to segment given its high variability within subjects, can be seen due to the regularization. The optimal values of α coefficient were 18 for lobule I--II and 4 for the rest of structures including the background label. To more specifically evaluate the impact on the label contour, we calculated the Hausdorff distance [START_REF] Beauchemin | On the Hausdorff distance used for the evaluation of segmentation results[END_REF] for each label. Despite the differences not being significant, regularization tends to produce lower distances.

Impact of extended library size

In this experiment we evaluated the impact of the number of subjects from the extended library used to perform the labeling step. The extended library included the original 10 manually labeled templates and 100 automatically segmented images from IXI dataset. From the library, the N most similar cases (in the mean squared error sense) to the target are selected. As can be noticed in Figure 3 the plateau is reached around 50 cases. Note that we have a noticeable oscillation for values higher than 50. As explained before, PatchMatch algorithm involves a random search step that is performed iteratively. In contrast to an exhaustive search PatchMatch converges faster because it doesn't explores the whole library. Thus we obtain an approximated nearest neighbor field. This implies that libraries with different sizes will have a different set of visited patches even if the majority of the templates are the same. For this reason, library sizes larger than the optimal (50) present oscillation around the optimal result. The use of this library extension enabled to improve the mean DICE coefficient from 0.7702 (N=10) to 0.7729 (N=50). Therefore, we set N=50 as the default parameter of the method.

Method comparison

The proposed method (with all the optimizations and the extended library) was compared with 3 state--of--the--art methods: SUIT [START_REF] Diedrichsen | A spatially unbiased atlas template of the human cerebellum[END_REF], MAGeT [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF][START_REF] Chakravarty | Performing label--fusion--based segmentation using multiple automatically generated templates[END_REF] and RASCAL [START_REF] Weier | Rapid Automatic Segmentation of the Human Cerebellum and its Lobules (RASCAL) --Implementation and Application of the Patch--based Label--fusion Technique With a Template Library to Segment the Human Cerebellum[END_REF]. The different compared methods were applied to segment the validation dataset to evaluate their performance. In Table 5 CERES method is compared with MAGeT and RASCAL methods for all 12 structures and the whole cerebellum. As can be noticed, CERES outperformed both method overall (only MAGeT had a slightly better result for lobule VI). As SUIT works with a different set of labels, the output segmentations of MAGeT, RASCAL and CERES were adapted by fusing lobules I--II, III and IV for a meaningful comparison as done in [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF]. Lobules Crus I and IX were not taken into account due to differences in their definition. We processed the validation images using SUIT as described in http://www.icn.ucl.ac.uk/motorcontrol/imaging/suit_function.htm without applying any manual correction. In Table 6, DICE coefficient values for the different compared methods are shown. CERES provided the best results overall (MAGeT results were better for lobule I--IV, Vi and X although differences were not significant). In Figure 4, some example images from the different method segmentations are shown.

Execution time

One of the main advantages of our method is its efficiency as it produces competitive results in terms of accuracy with a reduced execution time. The slowest method in the comparison was MAGeT (approximately 4 hours per case) due to the need of performing many non--linear registrations. RASCAL takes around 110 minutes to provide the segmentation result. Finally, SUIT (without manual correction) and CERES require around 5 minutes to produce the segmentation.

Figure 5 shows the computational cost of the different parts of the proposed method. As can be noticed the most time consuming part is the non--linear registration estimation. Note also the segmentation part took only 1 minute thanks to the efficiency of the OPAL method, while the preprocessing remains the most time consuming step of the proposed pipeline.

CERES, this time is highly reduced by using OPAL making the patch matching process no longer the bottleneck of the segmentation procedure.

Another important feature of CERES is its adaptability. As it uses a library of manually segmented templates as a knowledge base it can learn new anatomies by adding new cases to the library. In our experiments, the main limitation of CERES is the small size of the manually labeled template library (only 5 cases). In previous multi--atlas label fusion approaches it has been shown that a library of at least around 20 cases is required in order to achieve optimal segmentation results. Therefore, it is expected that CERES performance will be improved when using a larger library. This will be done by adding more manually segmented cases as long as they become available.

Finally, in order to make CERES pipeline fully and freely accessible to the scientific community, we have integrated this pipeline as a part of our online volBrain platform (http://volbrain.upv.es). At the time of writing this paper, CERES pipeline has been already applied to more than 3000 cases (in around 5 months) with a failure rate lower than 1%.

Figure 1 :

 1 Figure 1: Example of one HR case and its corresponding manual segmentation.

Figure 2 :

 2 Figure 2: Segmentation pipeline overview.

Figure 3 :

 3 Figure 3: Evolution of DICE coefficient by number of cases from the extended library used for segmentation over the validation set.

Figure 4 :

 4 Figure 4: Manual, SUIT, MAGeT, RASCAL and CERES segmentation examples over the cases with the best, median and worst DICE for CERES. Corresponding DICE is provided for each segmentation. Note how SUIT tends to produce oversegmentation and always fills the fissures between lobules. Also note the irregular boundary between lobules I--IV(red) and V(yellow) produced by RASCAL in all three cases. Finally, in the worst case, it can be seen that lobule VIIIA(purple) in MAGeT and RASCAL segmentations is barely visible.

Table 1 :

 1 Mean and standard deviation DICE results for CERES segmentation over the validation set.Columns show the average values and standard deviation for the lobules and the whole cerebellum. The values are calculated in the native space. Best results marked in bold.

	Lobule	Linear reg.	Non--linear reg.	Subject specific
	AVERAGE	0.7207 ± 0.0574	0.7549 ± 0.0452	0.7595 ± 0.0474
	Lobule I--II	0.4022 ± 0.2157	0.4165 ± 0.1993	0.4226 ± 0.2110
	Lobule III	0.6644 ± 0.2039	0.7044 ± 0.1728	0.7068 ± 0.1714
	Lobule IV	0.6879 ± 0.1625	0.7156 ± 0.1301	0.7213 ± 0.1374
	Lobule V	0.7263 ± 0.1424	0.7448 ± 0.1315	0.7514 ± 0.1376
	Lobule VI	0.8392 ± 0.0437	0.8620 ± 0.0322	0.8667 ± 0.0322
	Lobule Crus I	0.8678 ± 0.0297	0.8937 ± 0.0174	0.8970 ± 0.0161
	Lobule Crus II	0.7564 ± 0.0668	0.7973 ± 0.0545	0.8024 ± 0.0553
	Lobule VIIB	0.6076 ± 0.1354	0.6667 ± 0.1107	0.6742 ± 0.1107
	Lobule VIIIA	0.7277 ± 0.0705	0.7783 ± 0.0691	0.7830 ± 0.0692
	Lobule VIIIB	0.8012 ± 0.0541	0.8456 ± 0.0383	0.8499 ± 0.0358
	Lobule IX	0.8623 ± 0.0383	0.8813 ± 0.0322	0.8849 ± 0.0313
	Lobule X	0.7053 ± 0.0830	0.7529 ± 0.0438	0.7540 ± 0.0473
	Cerebellum	0.9274 ± 0.0107	0.9346 ± 0.0091	0.9358 ± 0.0091

Table 2 :

 2 Mean DICE results for CERES segmentation over the validation set. Columns show the average values for the different lobules and the whole cerebellum for the two compared approaches. Best results marked in bold.

	Lobule	α = 0.5	α adaptive (LOOCV)
	AVERAGE	0.7595	0.7640
	Lobule I--II	0.4226	0.4364
	Lobule III	0.7068	0.7237
	Lobule IV	0.7213	0.7281
	Lobule V	0.7514	0.7517
	Lobule VI	0.8667	0.8673

Table 3 :

 3 Values for the adaptive coefficient optimized using a gradient descent technique.

	Label	α coefficient	Label	α coefficient
	Lobule I--II	0.82	Lobule VIIB	0.59
	Lobule III	0.81	Lobule VIIIA	0.59
	Lobule IV	0.71	Lobule VIIIB	0.61
	Lobule V	0.59	Lobule IX	0.60
	Lobule VI	0.50	Lobule X	0.79
	Lobule Crus I	0.81	WM	0.51
	Lobule Crus II	0.60	Background	0.82

Table 4 :

 4 Mean DICE and Hausdorff distance results for CERES segmentation over the validation set. Columns show the average values for the lobules and the whole cerebellum. The values are calculated in the native space. Best results are marked in bold.

	Structure	No regularization		Regularization
		DICE	Hausdorff	DICE	Hausdorff

Table 5 :

 5 Mean DICE coefficient and standard deviation values forMAGeT, RASCAL and CERES.

	Structure	MAGeT	RASCAL	CERES	Intra--rater manual
	AVERAGE	0.7320 ± 0.0568	0.6890 ± 0.0524	0.7729 ± 0.0427 †	0.833
	Lobule I--II	0.3960 ± 0.1424	0.3260 ± 0.2178	0.5201 ± 0.1555* †	0.639
	Lobule III	0.6800 ± 0.1741	0.6379 ± 0.2165	0.7213 ± 0.1572	0.751
	Lobule IV	0.6980 ± 0.1440	0.6627 ± 0.1611	0.7271 ± 0.1346	0.818
	Lobule V	0.7320 ± 0.1398	0.6666 ± 0.1560	0.7561 ± 0.1332 †	0.881
	Lobule VI	0.8710 ± 0.0359	0.7969 ± 0.0523	0.8695 ± 0.0316 †	0.912
	Lobule Crus I	0.8870 ± 0.0257	0.8383 ± 0.0351	0.9007 ± 0.0152 †	0.904
	Lobule Crus II	0.7780 ± 0.0679	0.7340 ± 0.0667	0.8096 ± 0.0569 †	0.900
	Lobule VIIB	0.5990 ± 0.1487	0.5820 ± 0.1137	0.6850 ± 0.1205* †	0.863
	Lobule VIIIA	0.7300 ± 0.0934	0.6757 ± 0.1426	0.7926 ± 0.0759* †	0.860

Intra--rater accuracy for the manual segmentations from the validation set is also shown for evaluation purposes. Best results are marked in bold (significant differences (p<0.05) are marked with * for

MAGeT and † for RASCAL comparison)

.

Table 6 :

 6 Mean and standard deviation for DICE coefficient values for SUIT, MAGeT, RASCAL and CERES for the SUIT segmentation protocol. Best results are marked in bold (significant differences (p<0.05) are marked with * for SUIT, † forMAGeT and ‡ for RASCAL comparison).

	Structure	SUIT	MAGeT	RASCAL	CERES
	AVERAGE	0.7097 ± 0.0689	0.7639 ± 0.0792	0.7140 ± 0.0487	0.7888 ± 0.0409 * ‡
	Lobule I--IV	0.7435 ± 0.0880	0.8055 ± 0.0964	0.7703 ± 0.1107	0.7898 ± 0.1021 *
	Lobule V	0.6598 ± 0.1026	0.7429 ± 0.1369	0.6730 ± 0.1560	0.7561 ± 0.1332 * ‡
	Lobule VI	0.7800 ± 0.0543	0.8762 ± 0.0365	0.7994 ± 0.0523	0.8695 ± 0.0316 * ‡
	Lobule Crus II	0.7430 ± 0.0631	0.7787 ± 0.0678	0.7300 ± 0.0667	0.8096 ± 0.0569 * ‡
	Lobule VIIB	0.5701 ± 0.1572	0.6013 ± 0.1476	0.5761 ± 0.1137	0.6850 ± 0.1205 * ‡
	Lobule VIIIA	0.7134 ± 0.0996	0.7330 ± 0.0928	0.6701 ± 0.1426	0.7926 ± 0.0759 * † ‡
	Lobule VIIIB	0.7721 ± 0.0596	0.8012 ± 0.0607	0.7654 ± 0.0931	0.8533 ± 0.0390 * † ‡
	Lobule X	0.6955 ± 0.0512	0.7721 ± 0.0475	0.7275 ± 0.0680	0.7548 ± 0.0469 *

Acknowledgements

We want to thank the IXI --Information eXtraction from Images (EPSRC GR/S21533/02) datasets promoters for making available a valuable resource to the scientific community. This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (HL--MRI ANR--10--IDEX--03--02), Cluster of excellence CPU, TRAIL (HR--DTI ANR--10--LABX--57) and the CNRS multidisciplinary project "Défi ImagIn". This research was also supported by the Spanish grant TIN2013--43457--R from the Ministerio de Economia y competitividad.

Discussion

We have presented a new pipeline for cerebellum lobule segmentation called CERES that works in a fully automatic manner and that is able to provide accurate results in a low computation time. We have shown that CERES produces state--of--the--art results outperforming current cerebellum segmentation methods such as SUIT, MAGeT and RASCAL. Main contributions of the presented method are: 1) the use of an efficient custom made non--linear library for each case to be segmented 2) an adaptive multi--scale implementation and 3) a new label regularization step.

In terms of accuracy, CERES obtained the best results for most of the labels. Each of the improvements proposed in the method allowed to consistently leverage the results being the custom library construction the most significant one (0.7207 to 0.7595). The adaptive multi--scale mixing and the label regularization steps allowed to further boost the results to obtain 0.7729 final DICE coefficient. More interestingly, CERES almost reaches intra--rater accuracy for the whole cerebellum segmentation (CERES = 0.9377 vs. human = 0.941) [START_REF] Park | Derivation of high--resolutionMRI atlases of the human cerebellumat 3 T and segmentation using multiple automatically generated templates[END_REF]. It is worth noting that the high accuracy obtained by CERES was achieved using only 5 manually labeled templates. It is expected that results with larger manually labeled template libraries will provided significantly better results by covering a larger variability of cerebellum shapes and conditions (including pathological ones).

From the efficiency point of view, CERES was found to produce the segmentation results in a very short time (5 minutes per case) being the most efficient of the compared methods. This can be explained by two main reasons. First, multi--atlas label fusion techniques usually require the processing of multiple non--linear registrations which is quite time consuming. CERES only requires a single non--linear registration since it applies the N non--linear mapping from a set of pre-computed transformation to the MNI space. Furthermore, this single non--linear registration being low--dimensional and applied to a cropped volume is much faster estimated that high dimensional registration typically used in multi--atlas methods. Second, one of the most computationally demanding steps in non--local multi--atlas label fusion methods is the patch matching process. In