Zeyun Shi

Pierre Alliez

Mathieu Desbrun

Hujun Bao

Jin Huang

Symmetry and Orbit Detection via Lie-Algebra Voting

: Adjoint invariant distance to orbit. Our Lie algebra voting approach to symmetry and orbit detection maps SE(3) transformations into points in a logarithmic space composed of a rotation part ω ∈ R 3 and a translation part u ∈ R 3 . The rotational orbit of the church and the translational orbit of the side railing (a) are mapped into collinear blue and red spheres respectively (a few transformations within these two orbits are marked with circled numbers to enhance comprehension). When the scene is centered, the two lines are orthogonal to each other and easy to distinguish (b). However, after a rigid translation of the scene, the rotational orbit now has u-values near the translation orbit points, making it impossible to automatically distinguish these two orbits using a Euclidean distance (d), while our adjoint invariant distance for orbit shows no discernible difference in results as evidenced by a binning of detected orbit sizes for both situations (e).

Introduction

Symmetry plays a central role in geometry. Whether they be in natural or man-made shapes or scenes, symmetric features represent visually defining properties that are easily recognizable, and can thus be exploited for semantic classification. Detecting these (global or local) symmetries is thus a core task in many applications such as modeling, reconstruction, and segmentation as it helps increasing robustness of further processing.

Previous work.

The importance of symmetry in geometry processing has resulted in a variety of approaches to detect near- † Corresponding author: hj@cad.zju.edu.cn symmetries in geometric datasets-a comprehensive review is provided in [START_REF] Wand | Symmetry in 3D geometry: Extraction and applications[END_REF]. Symmetry detection methods can be categorized as intrinsic or extrinsic. While intrinsic symmetry detection uses geodesic distances and considers isometric deformation as in [LTSW09, XZT * 09, JXCZ13], extrinsic symmetry detection typically looks for rigid transformations instead. We restrict our exposition to methods detecting partial, approximate and extrinsic symmetry in 3D meshes.

A number of methods have been proposed to identify the presence of global symmetries, based for instance on spectral analysis [BAK10, LCDF10], feature graph matching [BBW * 08, BBW * 09, THW * 14], descriptors [KFR04], or symmetry transforms [PSG * 06]. While a few of these methods can be used to detect local symmetries as well, finding partial symmetry within an input geometry is most efficiently achieved through voting [MGP06], by proceeding as follows. After an initial probing of the geometry that finds pairs of small, similar regions, a series of candidate transformations between each pair of neighborhoods are mapped as points in a high-dimensional space-typically R 7 to encode the Euler angles, translations, and scale factor of each transformation. Clustering in this transformation domain is then performed to find high-density regions, which are evidence of potential symmetries. Depending on the type of symmetry that one wishes to look for, various mappings to the transformation domain can be designed: any parameterization of the relevant symmetry group can be used, and clustering is performed in this parameter domain identified with R n , so that a notion of distance between transformations can be trivially evaluated using the Euclidean metric. This voting framework is intuitive and flexible, and has been shown to properly confirm a number of important symmetries.

Finding orbits, i.e., regular structures generated by combining k transformations, is more challenging, and as such, it has been significantly less studied. Translational orbits, both discrete and continuous, were discussed in the context of pattern-aware shape deformation [BWKS11]. Pauly et al. [PMW * 08] demonstrated that the voting framework can handle orbit detection too: with a proper choice of the map to the transformation domain, orbits are identifiable as lattice structures. Therefore, a more thorough analysis of the geometric configuration of samples in the transformation domain can detect more than just partial symmetries. However, specific maps are required to recognize specific orbits, stressing the importance of the choice of transformation domain.

Symmetry detection for urban scenes and building facades have also been sought after. Methods for this purpose often make assumptions on the orientation of the scene and of the facades. For example, in [START_REF] Kerber | Scalable symmetry detection for urban scenes[END_REF] an urban scene is considered to have consistent upward orientation, leading to a use of orientation histogram descriptors to detect symmetry candidates. In [ZXJ * 13], atomic elements of facade images are extracted via box abstraction, before being hierarchically decomposed by splitting and layering of facade structures. Instead, our approach detects local (even approximate) symmetries with no prior knowledge on the original position and orientation of the data to be robust to arbitrary inputs.

Current limitations. Despite the success of voting methods to find symmetries and orbits, we point out that this family of approaches suffers from a number of important shortcomings that limit the practical detection of shape symmetries. For instance, they rely on a customized mapping function from the space of transformations to a linear space based on which symmetries (or composition thereof) are looked for. Repeated mappings and clustering may thus be needed to detect a complete set of relevant symmetries. This lack of a unified transformation mapping also results in a need for various clustering methods to detect the presence of features such as concentration points or lattices; note that most of these clustering methods require parameter tuning of their own to perform efficiently and reliably. Finally, the notion of distance with which clustering is achieved is often coordinate system dependent, i.e., the symmetry and orbit detection can be significantly impacted by the translation, rotation, or scaling of the input geometry. This problem (detailed in Sec. 4.1) is particularly undesirable since what we call a symmetry in Euclidean space is, by very definition, frame and scale invariant. In practice, these limitations imply that one can only find symmetries reliably when one already knows they are present in the input geometry since manual tuning is necessary (see, e.g., Fig. 6). With all these issues, robust, "sight-unseen" symmetry discovery of geometric datasets is still a challenging task.

Contributions.

We propose a simple approach to (partial, local and global) symmetry and orbit detection for 3D datasets, which unifies all existing voting-based techniques into a single, versatile Liealgebra voting framework. Using the group structure of transformations that are combinations of rigid motions and scaling, we advocate the use of a logarithmic map to naturally encode symmetries and orbits as linear subspaces of the transformation domain. We also identify the importance of using an adjoint-invariant metric on the resulting (Lie algebra) space to enable robust, coordinate frame independent clustering of linear subspaces from a series of sampled transformations obtained through probing of the input geometry. Our group theoretic approach helps discover, systematically and automatically, relevant symmetry and orbit in arbitrary geometric datasets by removing many of the limitations of current methods.

Background

We begin our exposition by reviewing concepts and mathematical definitions which will be at the core of our contribution.

Transformations as Lie groups

Symmetry in shapes in 3D Euclidean space has been extensively studied since the antiquity. Most relevant to computer graphics are the notion of isometries, scaling, and reflections, as what we call shape is always defined up to an arbitrary translation, rotation, scale or even flip. Consequently, the transformations from the Euclidean space to itself that do not affect the shape of an object are indispensable tools to detect symmetries and orbits. The groups of all 3D rotations (SO(3)), rigid body motions (rotations and translations, denoted SE(3) ≡ SO(3) R 3 ), or even the more general case of similarities (rotation, translation and uniform (non-zero) scaling, denoted SIM(3)≡SE(3) R + ) have thus found countless applications in graphics (see, for instance, [START_REF] Huang | Linear combination of transformations[END_REF]). This latter group of similarities will be most relevant to our exposition, as the orientationreversing case of reflections will be treated algorithmically instead of geometrically (see Sec. 5.1). We note that similarity transformations in SIM(3) form not only a group (elements can be composed and always have an associated inverse), but a Lie group [START_REF] Murray R | A Mathematical Introduction to Robotic Manipulation[END_REF], i.e., it has a differential structure.

Lie algebra.

The tangent space at the identity of a Lie group canonically defines a Lie algebra, which we will denote as sim(3). Intuitively, the Lie algebra of the group of similarities represents infinitesimal similarity transformations, and thus offers a linear space where one can easily sum or subtract infinitesimal elements. A Lie group and its Lie algebra are strongly related, and one can perform calculations in one by first mapping to the other: the exponential (exp) and logarithm (log) maps can convert back and forth between Lie group and Lie algebra. Note that some restrictions apply in general; however, our work will only consider rotations limited to 2π, removing a number of ambiguities present in the general case.

Representation of elements.

A transformation T ∈ SIM(3) is typically represented as a 4×4 matrix of the form:

T = R t 0 w -1 (1)
where R∈ SO(3) is a 3×3 rotation matrix, t ∈ R 3 a translation, and w∈R is a scale. Applying this transformation to a geometric object S simply rotates S by R, translates it by t, then scales it by w. The matrix logarithm of this group element T maps it to an element of the Lie algebra which is also typically written out as a 4×4 matrix, this time expressed as [MLS94, Ead13]

log(T ) =     0 -ω 3 ω 2 u 1 ω 3 0 -ω 1 u 2 -ω 2 ω 1 0 u 3 0 0 0 -s.     (2)
The vector ω = (ω 1 , ω 2 , ω 3 ) can be understood as an infinitesimal rotation (where the direction of ω encodes the axis of rotation, while its norm is an angular speed between 0 and π), while u = (u 1 , u 2 , u 3 ) is akin to an infinitesimal translation and s to an infinitesimal change of scale factor. These 7 independent values in log(T ) can be grouped into a single vector log (T ) ≡ (ω t , u t , s) t . It should be noticed that in general, u = t since rotations and translations are intertwined in the calculation of the log matrix. The mapping T → log (T ) is often referred to as the logarithmic map. Note also that this logarithm of a transformation T has a clear geometric meaning, particularly when w = 1: in that case, s = 0, so ω and u represent the constant angular and linear velocity of a screw motion [vKC99] that brings the identity element of SIM(3) to T . With the scale added, one can think of the logarithm of a transformation T as a dilating screw motion, since the scale is also changing (i.e., growing or shrinking) along the screw motion.

Adjoint action. While a Lie algebra is the tangent space of its Lie group at the identity element, it is often important to work in the tangent space of an arbitrary transformation. The adjoint serves exactly this purpose, since it linearly transforms tangent vectors from one tangent space of the Lie group to another [START_REF] Murray R | A Mathematical Introduction to Robotic Manipulation[END_REF]. This adjoint map is easy to define from the notion of conjugation by an element F ∈ SE(3): a similarity T is conjugated by F through the composition FT F -1 . The adjoint map is the differential of conjugation at identity; that is, for an arbitrary element F of SE(3), the adjoint map Ad F with respect to F is the linear operator acting on Lie algebra elements as:

Ad F =   R 0 0 t × R R -t 0 0 1   , (3)
where t × is a skew-symmetric matrix satisfying t × d = t×d for an arbitrary vector d ∈ R 3 . Conjugation and adjoint will be especially useful when we talk about change of coordinate frames later. These group actions and geometric entities are thus not purely abstract: they, in fact, offer a coherent framework to manipulate 3D transformations in a principled manner.

Transformation voting strategy

Now that we have a general notion of shape congruence through the action of similarity

T = ∏ i G ri i
for arbitrary integers or reals r i (in continuous orbits), and the presence of such orbits signals the existence of repetitive patterns in the shape. Unlike a symmetry which is captured by a single transformation (e.g., a point in the transformation domain), a k-orbit is usually mapped to a k-dimensional linear structure (or lattice if r i ∈ N)-and can thus be detected by applying a robust shape detection approach such as RANSAC (repeated stochastic sampling consensus [FB81]) to the point cloud in the transformation domain.

In the remainder of this paper, we show that adopting a Lie algebra viewpoint when using the voting framework provides significantly more robust foundations for the discovery of symmetries and orbits.

Lie Algebra Domain

While previous works formulated various mappings from elements of SIM(3) to the transformation domain, we propose the use of a logarithmic map instead. This simple, easily computable map to the Lie algebra of similarity transformations satisfies key properties that will ease the process of symmetry and orbit detection through voting. It is also generic, as it can be used for most previouslystudied symmetries and orbits.

Logarithmic mapping

We introduced the logarithm map log in Sec. 2.1 which converts any element T of SIM(3) to a vector log T ∈ R 7 encoding the matrix logarithm of T . This operation offers a very convenient map to a general transformation domain. Just like the initial approach of Mitra et al.

[MGP06], it is seven dimensional; it is, however, a local diffeomorphism as well: transformations that are quite similar will be mapped to nearby points in R 7 . Such a basic requirement was not satisfied by the Euler angle based mapping of [MGP06]. Using this log-based map thus leverages the differential structure of Lie groups, since we directly use the Lie algebra representation of a transformation as its coordinates in the transformation domain.

Orbits as linear structures

Because the Lie algebra sim(3) is a linear space, we reap additional benefits as well. If F and G are two transformations of SIM(3) that commute, then the composition of F and G satisfies:

FG = exp(log(F) + log(G)).
Thus, if a k-orbit T is generated by transformations {G 1 , ..., G k } which commute with each other, then its logarithmic mapping to the transformation domain has the following property:

log (T ) = log ( ∏ i G ri i ) = ∑ i r i log (G i ). (4)
Consequently, a k-orbit always generates a k-dimensional linear subspace (line for k = 1, plane for k = 2, etc...) in the mapped domain. Our logarithmic mapping thus generalizes the orbit detection proposed in [PMW * 08] as we no longer have to construct a particular mapping to detect each type of k-orbit: our log-based transformation reveals the linear structure of all types of orbits in a unified manner. Note finally that our approach, much like previous work, cannot handle the case of non-commutative orbit; fortunately, most relevant orbits are commutative [PMW * 08].

Basic notion of distance

Finally, the log map has a last natural property that is particularly convenient to extract orbits. For each mapping to a transformation domain, a reliable distance measurement dist(F, G) between two transformations F and G must be formulated. One could consider the composition of F and G -1 as a notion of difference δ = FG -1 between the two, then use a matrix norm of δ to obtain a scalar measurement. However, this leads to inconsistencies since δ = δ -1 in general, which means that the distance between F and G is not equal to the distance between G and F.

Instead, the matrix logarithm does satisfy log(δ) =log(δ -1 ); we can thus define a consistent distance via the Euclidean norm (or any weighted variant) of the vector log (δ), such that:

dist 2 (F, G) = log (δ) t   α Id 0 0 0 β Id 0 0 0 γ   log (δ) ≡ log (δ) 2 E .
The parameters α, β and γ have to be fixed by the user as they reflect the arbitrariness of the importance ratio between a unit translation vector, a unit rotation and a unit scaling. Note that for α = β = γ = 1, the norm • E becomes the usual Euclidean norm of a vector in R 7 . The distance between F and G can be measured directly in the transformation domain as the (possibly weighted) Euclidean distance between the two points corresponding to F and G:

dist(F, G) = log (F) -log (G) E .
While this simple notion of distance in our transformation space is already substantially more reliable than previously proposed distances, we show next that another metric can be devised to further increase robustness of the symmetry discovery process.

Optimal Distances in Lie Algebra

Due to the shape congruence we discussed in Sec. 2.1, what constitutes a symmetry or orbit of an input object should not depend on its location in Euclidean space or on its scale. This invariance, which is key to robustness in shape analysis, has been surprisingly overlooked in previous work. Next, we discuss the specific mathematical requirements it imposes on the metric in our transformation domain, and deduce a practical notion of distance that vastly improves robustness over previous work.

Adjoint invariance

Consider an arbitrary rigid transformation F ∈SE(3). If one applies this transformation to an input geometry, and feeds this rotated and translated object to a symmetry detection algorithm, one should expect the output to be exactly the same as if F had not been applied: after all, applying the rigid F does not change the equivalence class in which the object is since its shape is entirely preserved. However, the transformations T i generated by the voting strategy described in Sec. 2.2 will be each transformed into their conjugation by F, i.e., they become of the form FT i F -1 . Alas, these conjugated transformations do not have the same pairwise distances as their non-conjugated counterparts in any of the transformation domains proposed in the past-and even the canonical metric • E defined in Sec. 3.3 is only invariant to rotations of the input. This is not just a mathematical wrinkle: this lack of invariance to similarity is a huge impediment to the discovery of obvious symmetries as demonstrated in Fig. 2, where a simple rotation and translation of a 3D model is enough to disrupt the proper extraction of trivial symmetries. Indeed, following the approach of [MGP06] (Sec. 2) and discarding the scale for simplicity, two sample points p i and p j on this model would find a rotation R i j to align their respective local frames and a translation t i j = p j -R i j p i . If the brick model is now translated globally by a vector T , then the new translation would become t i j = t i j +(I -R i j )T , which can be very different from t i j in general (even on flat surfaces if noise is present), preventing the detection of even simple structures. In order to enforce invariance to the placement and scale of the input, we instead need to find a special metric in the transformation domain that is adjoint invariant, i.e., that satisfies for any pair of transformations A and B:

log (A) -log (B) = Ad F log (A) -Ad F log (B) .
Intuitively, this is a Lie group extension of the idea of frame invariance for pairwise distances between points; but in our more complex setup, this is far from an easy requirement to satisfy. Not surprisingly, this type of requirement on the metric has been studied in robotics at the end of the '90s for the case of SE(3), see for instance [START_REF] Murray R | A Mathematical Introduction to Robotic Manipulation[END_REF]vKC99]. We can easily extend their findings to SIM(3): once we write out the requirements that an adjointinvariant metric must satisfy, we discover that only a 3-parameter family of quadratic forms M of the form:

M =   α Id 1 2 β Id 0 1 2 β Id 0 0 0 0 γ   (5)
is able to enforce this invariance. For these quadratic forms, one can trivially verify the adjoint-invariance property:

∀δ ≡ (ω t , u t , s) t ∈ sim(3), Ad F δ M = δ M = α ω, ω +β ω, u + γs 2 . (6)
Unfortunately, the resulting notion of distances (whose geodesics are screw motions [vKC99]) corresponds to only pseudo-Riemannian metrics: even after restricting the three parameters to be positive, this distance in the transformation domain can be negative, or even blind to translations when ω = 0. In other words, there is no perfect metric for our purpose. Despite this disappointing realization, one can still produce a useful metric that will be adjointinvariant over most of the transformation domain, which will be enough to render clustering significantly more robust.

Variational Distance

Although there is no perfect adjoint invariant metric on SIM(3), we propose a variationally-defined metric which is adjoint-invariant almost everywhere. We define distance • A based on the smallest norm over all adjoints of rigid body transformations, i.e., ∀δ ∈ sim(3),

δ A = min F Ad F δ 2 E , F ∈ SE(3). ( 7
)
That is, this distance should represent the smallest (weighted) Euclidean metric of δ over all possible coordinate frames defined by a coordinate system change F, which corresponds to the "shortest" screw motion from Id to δ. Expanding this definition via Eq. (3) leads to a simpler expression of the form:

Ad F δ 2 E = α ω 2 2 + β t × Rω + Ru -st 2 2 + γs 2 ,
The s component in δ here measures the difference of scaling between two sampled transformations, and is usually small. Thus, we approximate the above equation as

Ad F δ 2 E ≈ α ω 2 2 + β t × Rω + Ru 2 2 + γs 2 = α ω 2 +β u 2 +Ψ R,t (ω, u) +γ s 2 , where: Ψ R,t (ω, u) = ω T R T t T × t × Rω + 2ω T R T t T × Ru. The fact that R T R = Id and R T t × R = (R T t) × also implies Ψ R,t (ω, u) = Ψ Id,R T t (ω, u).
Thus, we can always assume R = Id without changing the value of the minimum over all F. Substituting R = Id into the optimization problem defining the norm, one gets:

δ A = min t α ω 2 + β t × ω + u 2 + γ s 2 . ( 8
)
This quadratic problem has a simple closed form solution. When ω = 0, u can be decomposed in the component u = ω/ ω , u ω parallel to ω and the remaining part u ⊥ = uu orthogonal to ω, thus picking t × ω = -u ⊥ leads to the solution. Thus, we get:

δ A = β u 2 + γ s 2 ω = 0 α ω 2 2 + β ω/ ω , u 2 + γ s 2 otherwise , (9)
Note that this formulation changes the term ω, u in Eq. ( 6), eliminating the issue of negative distances and returning a reasonable value when ω = 0. However, we traded positive-definiteness for continuity: this distance is now discontinuous at ω = 0. We regularize this distance by defining a vector ω increasingly aligned with u as the norm of ω approaches zero, and equal to the normalized vector ω/ ω otherwise, through:

ω = ω/( ω + ε)
for ε small (we use ε = 10 -6 in all our experiments). We can now define the final expression of this regularized distance as:

δ A ≡ α ω 2 2 +β((1-ω ) u + ω ω, u ) 2 + γs 2 . ( 10
)
This new distance definition combines the benefits of continuity and adjoint invariance (and consequently, of robustness) away from ω = 0. Equipped with this distance on sim(3), we can proceed as advocated in [MGP06] with a mean shift algorithm [START_REF] Cheng | Mean shift, mode seeking, and clustering[END_REF] to cluster samples in the transformation domain, where the sample mean m(δ) of a sample set S at δ is computed as

m(δ) = ∑ δs∈S K( δ -δs A ) δs / ∑ δs∈S K( δ -δs A ), (11)
where K used in our method is the unit Gaussian kernel. The resulting cluster centers will serve as reliable detectors of the significant symmetries present in the input.

Adjoint-invariant distance to orbit

The distance defined in Eq. (10) provides a measurement between mapped transformations, so that a voting strategy can be used to find the most significant symmetries. In order to detect orbits, however, we also need to provide a distance between a point in the transformation domain and an orbit. As orbits have a linear structure in our construction (see Sec. 3.2), a k-generator orbit can be represented by its k generators as a base matrix G of size 7×k. Any element in the orbit can be written as Gr with r ∈ R k ; thus an intuitive definition of point-to-orbit distance could be the minimal Euclidean distance from a point p to this linear structure, i.e., min

r p -Gr E . (12)
The k dimensional vector r can be viewed as the coefficients of the projection point of p in the orbit structure. From this quadratic optimization, we can directly derive an explicit formula for r, and compute the distance. However, this definition suffers from not being frame or scale invariant.

Extending our design of a distance between points, we can also define an adjoint-invariant point-to-orbit distance using:

dist G (p) = min F,r Ad F p -Ad F Gr E . (13)
Following the same reasoning as in Section 4.2, we find that the rotation R of F is also irrelevant to this distance definition: the above optimization problem can be solved by fixing R in F as Id, and solving in alternation for t and r via simple quadratic optimization.

Because the Hessian of t can be singular, we use Singular Value Decomposition to compute the Hessian's pseudo-inverse. Starting from r = 0 as the initial value, the optimization usually ends in less than 10 iterations, where convergence is tested through the gradient norm. After getting t and r, we use Eq. ( 10) to compute the distance to avoid the discontinuity issue.

We can now use RANSAC [FB81] to find the k-generator orbits present in the sampled transformations. Grid fitting can also be applied to each orbit as in [PMW * 08] if needed.

Results

Using the Lie algebra representation of transformations and our new adjoint-invariant distances, we can now reliably employ the voting framework to discover symmetries and orbits. In this section, we discuss how to efficiently implement our approach, and provide various results and comparisons with previous work.

Implementation Details

Our new transformation mapping and distances can be readily plugged into any existing voting framework. Our implementation mostly follows the pipeline advocated in [MGP06] (see their Fig. 2), which we briefly review to highlight our changes: 1. Sampling. A set P of sample points is uniformly selected from the model S; after pruning [MGP06], we get a reduced sample set P = {p i }.

Pairing.

A local frame is constructed for each point p i based on estimations of normal and principal curvature directions. From a random subset P ⊂ P with |P | = |P|/5, a transformation T i j is then computed for each pair p i ∈ P , p j ∈ P that aligns their local frames. Unlike Mitra et al.

[MGP06], who reject pairs if their distance in the signature space is too large, we compute the alignment error of neighboring patches surrounding p i , p j in terms of the average distance between the closest point pairs, and pairs with large alignment error are discarded. This procedure improves robustness to noise (see Fig. 13). 3. Clustering. The remaining set of transformations {T i j } is then clustered in logarithmic space using our adjoint-invariant distance metric. For symmetry detection, mean shift [START_REF] Cheng | Mean shift, mode seeking, and clustering[END_REF] clustering is performed, while RANSAC [FB81] is used to find orbits. 4. Extracting. As in [MGP06], we use each resulting cluster center as a representative transformation, and extract its associated symmetry region using a region growing algorithm.

We also point out a few important details that we used to further increase the accuracy and efficiency of symmetry and orbit detection.

Handling inverses. One simple, yet important issue that we identified is the ambiguity that may rise between a transformation and its inverse. As the pairs of samples that the initial sampling stage of the voting strategy generates are not ordered, one can arbitrarily find either a matching transformation T or its inverse depending on the order in which the samples are compared. While a large enough set of samples will create high densities around both transformations, this approach weakens the density of each transformation unnecessarily. Using our Lie algebraic representation of transformations, this problem can be addressed easily, since T maps to log (T ) while T -1 is mapped tolog (T ). Thus, the distance between two transformation A and B is further altered to read:

min( log (A) -log (B) A , log (A) + log (B) A ), (14)
to account for the inevitable presence of these doublets. This modified distance also avoids the typical issue due to the 2π periodicity of rotations: the distance between two rotations along the same axis by -π+ε and π-ε respectively will be small as expected.

Culling via metric estimates. Since symmetry and orbit detection relies on finding local accumulation of samples, any pair of transformations that are sufficiently dissimilar should be efficiently culled instead of being submitted to an exact distance computation. We obtain a speedup of a factor 3 by first checking the following lower bound on the distance defined in Eq. ( 10):

α ω 2 + γ s 2 ≤ δ A .
If this lower bound is beyond a given large threshold, we can safely skip further computations for this pair during clustering.

Dealing with reflections.

Recall that we did not consider reflections in our method so far, as their orientation-reversing effect creates complications in the mathematical definitions since reflection matrices have no real logarithm [START_REF] Culver | On the existence and uniqueness of the real logarithm of a matrix[END_REF]. We remedy this issue by applying a reflection to the whole model and adding the result to the input; more specifically: 1. Points are sampled from an input scene S and its YZ-plane reflection S = TYZS. (Note than any planar reflection of the original scene can be used instead: the purpose is simply to add one reflection, since all other reflections will be reachable via composition with rigid body transformations-which our method captures precisely thanks to our use of a logarithmic mapping.) 2. In the pairing step, we construct the set of transformations {T i j }

for which at least one point of p i , p j ∈ S ∪ S is in the original scene, i.e., we only deal with pairs satisfying {p i , p j } ∩ S = ∅ to avoid unnecessary computations. 3. In the clustering step, if more than half of the transformations in a cluster involves points coming from S, this cluster with transformation T is declared a reflection symmetry TYZT . 4. In the extracting step, for reflection symmetry TYZT , one of the seeds for patch growing p i , p j is in S, and after mapping it to S, the method proposed in [MGP06] can be directly used to extract the symmetry patches of S.

While this procedure is hidden from the user, this algorithmic trick is enough to deal with plane reflections with minimal changes to the code-see Fig. 3. However, it does double the computational cost of symmetry and orbit discovery. Note that a similar procedure can be applied to find point reflections.

Parameters. We used α = β = 1 and γ = 10 in all our experiments. The weight γ was set larger than α and β because we have found that the evaluation of the scale component of a transformation is typically not as reliable as the rotation and translation parts. The main reason is that the estimation of principal curvature uses a fixed neighborhood radius, which is suboptimal for point pairs with scaling symmetry [MGP06, PMW * 08]. The regularization parameter ε = 1e-6 is used in all our results, but any value in the range (T and T -1 are grouped together to increase reliability as we mentioned above). This pattern indicates that the symmetry close to the end of the orbit is harder to detect because of a smaller number of votes. Exploiting this characteristic helps us infer finer symmetries even for low sample counts.

Mean shift clustering.

An important parameter rarely discussed in symmetry detection is the kernel size σ used in mean shift clustering. We automatically set it through an initial calibration process: we randomly pick a small number of pairs of transformations and evaluate their variational distance • A to estimate the average distance h in the transformation set sampled from the object. We experimentally found that choosing a kernel size σ in the interval [0.01 h, 0.1 h] consistently leads to robust symmetry detection.

Comparisons of transformation spaces

In order to further demonstrate the advantages of our logarithmic mapping space over the usual R 7 space [MGP06], we compare the distribution of transformation samples on a series of model to provide a fair, parameter-free estimation of how a clustering algorithm is likely to find symmetries.

For a cluster C i of transformation samples (with center c i ), we first define its radial spatial distribution (inspired by [WS09]) as the function

τ i (r) = | f (c i , r) ∩ C i |/| f (c i , r) ∪ C i |,
where f (c i , r) represents the set of samples that are within a sphere centered at c i with radius r, and | • | represents the volume of a set. This notion of radial spatial distribution reveals the quality of a cluster, as shown in Fig. 4: when a cluster is clearly separated from other samples, its radial distribution will reach a plateau at its maximum value 1, from a minimum radius r min to a maximum radius rmax-see Fig. 4(a). Note that the ratio rmax/r min thus measures the compactness of the cluster. Instead, if a cluster overlaps other nearby clusters, then the maximum value that τ reaches indicates a separation distance of this cluster from other samples, and can thus be seen as a good measure of compactness-see Fig. 4(b).

In order to compare how the choice of transformation space affects symmetry detection, we need to generate a set of transformation samples for which we know what the results are supposed to be. We thus generate input scenes as follows:

1. For a given model M, we randomly generate three transformations T 1 , T 2 , T 3 , and apply them to M to get a new scene M composed of T 1 M, T 2 M and T 3 M. 2. We perform the sampling and pairing step of the traditional voting approach for M and get the transformation samples. 3. These transformation samples are then mapped to either our logarithmic, or the original R 7 transformation space.

By construction, the ground truth symmetries in this scene are known to be

T 2 T -1 1 , T 3 T -1 1 and T 3 T -1
2 (and their inverses, ignored for simplicity). Instead of performing the usual clustering which would be unfairly biased by the choice of (R 7 or variational) distance, we assign for each pair of transformations one of the three ground truth symmetries based on the one minimizing the local ICP error [BM92]-or we simply discard this transformation if the error is larger than 1% of the bounding box diagonal size. Note that this assignment is costly to compute, but has the advantage of being purely geometrical, and thus really indicate how intrinsically good a mapping is independent of the clustering method used on it. Finally, we measure the quality of the results using voting based on radial spatial distribution: a cluster (corresponding to one of the three symmetries) casts a vote for our logarithmic space or the R 7 space based on which embedding space leads to a larger max value of its radial spatial distribution-or for the embedding space with the larger ratio of rmax/r min if both spaces leads to the same maximum. This voting is done for 30 scenes that were randomly generated from a model. In our tests, the logarithmic space systematically outvotes the R 7 space by at least 50%, as demonstrated in the inset. Therefore, even if one ignores the advantages of our adjoint invariant clustering, the use of a logarithmic space is valuable in itself.

Adjoint invariance. Using radial spatial distributions, we can also demonstrate that our adjoint-invariant distance leads to improved robustness in the choice of kernel size used in mean shift clustering: as Fig. 5 illustrates, distributions change significantly when the scene is rotated and translated for the R 7 space mapping, while they remain unchanged in our method. This invariance to rigid transformations is crucial to robust symmetry and orbit detection.

2D visualization of transformation distributions. The role of adjoint invariance is illustrated in a simple example in Fig. 6, where a lamp model contains three identical flower-shape heads. As the visualization of a distribution of transformation samples in dimension 7 is difficult, we used multidimensional scaling (MDS) and multivariate kernel density estimation (MKDE) [START_REF] Venables | Modern Applied Statistics with S, Fourth ed[END_REF] to compute a 2D projection that best preserves the distances between samples in R 7 (a-d). Such a 2D visualization has the advantage of making the structure of this high-dimensional distribution visually obvious. While the distribution of transformations in our logarithm space is virtually unchanged for two different orientations and positions of the lamp (a-b), the original symmetry detection algorithm generates extremely different transformation distributions in R 7 (c-d). This behavior makes it impossible to pick a threshold for the clustering technique to decide which prominent clusters are potential symmetries. As a consequence, a lack of adjoint invariance leads to the inability to detect symmetries automatically and reliably. We point out that simply "normalizing" the position of the object or the scene by a global transformation does not help either as illustrated in Fig. 11, where multiple symmetries and orbits cannot be recovered with this procedure. In contrast, our metric • A performs as expected with no parameter tuning. Moreover, notice that these three symmetries (and their inverses) have cluster sizes that are clearly much larger than any other, so even without knowing the number or type of symmetries present in this example, the sudden change in cluster size (Fig. 6, bottom graph) is a clear indication that only three main symmetries are present. Our approach thus offers a sight-unseen detection of symmetries.

Ground truth comparison for a complex model

The necklace in Fig. 7 is composed of 32 components with identical links {C k }, k = 1, • • • , 32. The symmetry transformations between these components being known in advance, we can measure the accuracy of a symmetry detection algorithm by comparing the resulting clusters S = {S i } to the ground truth transformations G = {G i }, where each transformation G i maps a link C i -to another link C i + of the necklace. Ideally, the largest clusters in S should match the ground-truth symmetries G. We test this hypothesis for both our Lie-algebra voting and the R 7 voting of Mitra et al. [MGP06]. We Cluster centers corresponding to the correct transformations between the three heads are colored. While our Lie algebra approach is unchanged for different positions of the lamp, the original method differs widely based on position and orientation, with symmetries being sometimes co-located (see rectangle in (d)). We also show the distribution of cluster sizes (number of the transformations in each cluster, normalized by the largest one; bottom right), when a sudden change is visible with our Lie-algebra voting after the first six clusters-thus clearly identifying the three symmetries (and their inverses) as the only relevant ones. first find the "closest" transformation S j for each transformation G i ∈G via min j D(FG i F -1 , S j ), where F ∈ SE(3) is a random rigid transformation applied to the whole necklace to place the necklace into different positions, and D stands for the distance measurement (Euclidean distance for R 7 parameterization, adjoint invariant distance for our logarithm parameterization). We then define an error E(G i ) for each G i as the average pointwise matching error from C i -to C i + using S j (normalized by the bounding box diagonal of the model), since it measures how well we match each known symmetry. Finally, the average of all E(G i ) indicates how accurate the set of discovered symmetries S is. Fig. 7 shows the results (in log scale) of this test with mean shift clustering for 3 different positions (translation and rotation) of the necklace model for a large number of kernel sizes-a parameter that clustering heavily depends on. Our method produces identical results for all positions, while the R 7 approach exhibits very different behaviors. Moreover, our approach has better accuracy for most kernel sizes, or near equivalent accuracy for the best of the three positions. Finally, the R 7 approach fails to detect some of the symmetries G i in two out of the three positions. Note that the unit of kernel size on the graph in Fig. 7 is set to the minimal distance among the transformations in G, i.e. min Gi,G j∈G D(G i , G j ). The mean shift based symmetry detection error ∑ i E(G i ) as a function of cluster kernel size ratio is generally much smaller for our Lie-algebra voting approach compared to the R 7 method, and does not depend on the position of the model. Notice that for positions b and c, the R 7 method fails to detect some of ground truth symmetries no matter how the cluster kernel size is tuned.

Orbit detection

Linear structures. Fig. 8 demonstrates that our transformation domain exhibits clear linear structures when orbits are present, while previous R 7 representations do not share this property. The bumpy torus has a 6-fold rotational symmetry and many partial symmetries in the form of discrete rotations, while the Thai status model has a global 3-fold rotational orbit and the workpiece model has a clear translation orbit. The R 7 representation have no clear structures for these models, but our approach exhibits obvious (and relevant) structures.

Orbit detection. Our approach is able to discover both continuous and discrete orbits [PMW * 08,BWKS11]. Fig. 9 depicts the fire hydrant model containing a continuous rotational symmetry (highlighted by the dominant red line), the octopus model containing a 8-fold rotational orbit (clearly visible as 8 separate dots in the transformation domain), and the carter model with a 6-fold rotational orbit and a partial continuous rotational orbit (exhibiting a continuous line in the ω subspace, with 6 clear dots). Such clear features are thus easily detected and characterized by common clustering and linear subspace detection algorithms, without having to experimentally tune a variety of parameters. Robustness of the symmetry and orbit discovery process is thus much improved. The example in Fig. 1 contains a rotational orbit (represented by the blue dots) along the z axis and a translation orbit (represented by the red dots) along the x axis. When the scene is centered, orbits are easy to detect, even with a Euclidean distance, as the translation component of the rotational symmetry appears at the origin of the transformation domain. If the scene is translated along the y axis, however, the Euclidean distance can no longer distinguish the two orbits correctly. Our adjoint invariant orbit distance has no such issue, and all of our results were robustly handled without having to worry about whether the scene has to be relocated to detect symmetries and orbits.

Unified symmetry and orbit discovery

Detecting symmetries and orbits in a complex scene is a much harder task than detecting symmetries and orbits in a single model. Not only do the computational costs associated to both patch sampling and clustering increase, but the resulting complexity of the density field in the transformation domain can be quite high. However, because our approach satisfies the linear orbit property as well as an invariance to similarities, visualizing the transformation domain once basic clustering and RANSAC has been performed to hunt for symmetries and/or orbits is significantly easier. The scene in Fig. 11 includes various symmetries and orbits; without having to resort to many different mappings, all the SIM(3) symmetry and orbits are detected in a unified way. More importantly, if one uses a coordinate system dependent method instead, significant structures in the scene are missed because of their arbitrary relative locations.

Conclusions

We presented in this paper a new Lie algebra variant of the voting framework to discover, automatically and reliably, symmetries and orbits in geometric scenes. By leveraging the Lie group nature of symmetries and constructing a consistent metric that is fundamentally invariant to rigid body and scale change of the input models, we significantly improved the robustness and accuracy of the results. In particular, we introduced a logarithm mapping of transformations that unifies most previous symmetries and orbits, allowing for a significantly easier analysis of complex scenes since no manual tuning is required to discover symmetry and orbit. This unified Lie algebraic framework may also benefit common geometry processing tasks such as editing or reconstruction as it offers a principled way to identify, even without user interaction, partial symmetry in geometric data. In particular, we plan to investigate how symmetry and orbit detection can be used as a reliable prior in the context of reconstruction and abstraction of manmade geometry (where symmetries and orbits abound). Additionally, while we focused on finding simple geometric patterns (linear subspaces) in the voting space, detecting complex patterns is an interesting future work as it may allow detection of hierarchical and/or semantic structures in scenes.

Limitations. While we improved the reliability and robustness of the clustering step in the symmetry and orbit detection pipeline through a Lie-algebra voting and an adjoint-invariant distance, a few limitations of the original voting framework remain. First, our voting strategy is based on uniform sampling [MGP06], generating fewer points for symmetries involving small regions. Consequently, small symmetries have a smaller probability to be detected. For ex- ample, the 2π/18 symmetry of the corrugated teapot in Fig. 12 may be missed if too few transformation samples are used.

We showed in Fig. 6 that our method can automatically determine the number of significant symmetries without manually tuning the parameters by just tracking sudden changes in (normalized) cluster sizes. However, such clear, sudden changes may not be present in complex examples. For instance, our approach identifies Additionally, as our method depends on the accuracy of transformation estimates from sample pairs, very noisy data (Fig. 13) or symmetries with significantly different scales may incur estimation errors which will affect clustering quality. Similarly, while our experimental kernel size estimation for mean shift clustering (described in Section 5.1) works well in most of the cases, an adaptive strategy may be required to detect complex symmetries and orbits with many different scales.

Figure 2 :

 2 Figure 2: Lack of adjoint invariance. A lego brick has a 2generator translational orbit highlighted with the orange grid and yellow spheres (top left). While a "good" placement of this brick in space exhibits clear structures in R 7 (bottom left), simply rotating and translating the brick breaks this structure (right, two different views), rendering automatic orbit detection completely impossible.

Figure 3 :

 3 Figure 3: Reflections. Each point of these models is colored by its distance from the global reflection plane detected by our approach.

Figure 4 :

 4 Figure 4: Radial spatial distributions. When two clusters (blue and red, with their respective centers in green) are well separated, the radial spatial distribution τ of the blue cluster has a clear plateau with τ(r) = 1 (a). When two clusters overlap, the maximum of τ is lower than 1, but represents a good measure of compactness (b).

Figure 5 :

 5 Figure5: Invariance of radial spatial distribution. While a typical cluster's radial spatial distribution in our logarithmic space does not change if the scene is rotated and translated in space (blue; 3 different poses), the R 7 mapping space behaves unreliably (red).

Figure 6 :

 6 Figure 6: Adjoint invariance. The lamp model (bottom left) has 3 significant symmetries due to identical heads. For two different positions of the lamp in space, we plot both multidimensional scaling (MDS) and multivariate kernel density estimation (MKDE) [VR02] 2D visualizations of the resulting transformation distributions for our method ((a) and (b)) vs. the original R 7 method ((c) and (d)).Cluster centers corresponding to the correct transformations between the three heads are colored. While our Lie algebra approach is unchanged for different positions of the lamp, the original method differs widely based on position and orientation, with symmetries being sometimes co-located (see rectangle in (d)). We also show the distribution of cluster sizes (number of the transformations in each cluster, normalized by the largest one; bottom right), when a sudden change is visible with our Lie-algebra voting after the first six clusters-thus clearly identifying the three symmetries (and their inverses) as the only relevant ones.

Figure 7 :

 7 Figure 7: Ground truth comparison. (a) The necklace model. (b)The mean shift based symmetry detection error ∑ i E(G i ) as a function of cluster kernel size ratio is generally much smaller for our Lie-algebra voting approach compared to the R 7 method, and does not depend on the position of the model. Notice that for positions b and c, the R 7 method fails to detect some of ground truth symmetries no matter how the cluster kernel size is tuned.

Figure 8 :

 8 Figure 8: Linear structure of orbits. Our transformation domain (right, visualized through its ω or u coordinates for simplicity) always exhibits clear linear structures when orbits are present, while an Euler-angle based domain is unreliable. Transformation samples are colored by their alignment error (blue indicating small values). Cluster centers are shown as red dots.

cFigure 9 :

 9 Figure9: Rotational orbits. Our choice of Lie algebra coordinates renders rotational orbits particularly obvious in the transformation domain (middle), be they continuous (fire hydrant), discrete (octopus), or both (carter). We can then color the regions of the object exhibiting orbits (left) and/or visualize particular orbits (right).

Figure 10 :

 10 Figure10: Gallery. Our symmetry and orbit detection approach can find all the symmetries in these classical models with a single (logarithmic) mapping, and without the need to carefully move the model to a specific location.

Figure 11 :

 11 Figure11: Indoor scene. The dining room scene contains many symmetries and orbits; linear structures of orbits are very clear as illustrated by our visualization of ω and u components, despite the complexity of the scene. Bottom rows show the most significant orbits and some representative symmetries found in this scene using coordinate frames and lattices to indicate symmetries and orbits.

  992 significant symmetries of the necklace model from 45845 sampled transformations, but there is no clear demarcation that indicates that other clusters are simply small artifacts as the graph in the inset indicates. However, we significantly outperform the R 7 approach since it cannot even find all the symmetries for certain positions of the model.

c 2016 The Author(s) Author version

Acknowledgements. We would like to thank the anonymous reviewers for their valuable comments and suggestions. This work was partially supported by NSFC (No. 61522209, No. 61210007), the Fundamental Research Funds for the Central Universities (No. 2015XZZX004-19), the European Research Council (ERC Starting

Grant "Robust Geometry Processing" #257474), and the National Science Foundation's grant CCF-1011944.


Erhard Naess-Schmidt

Anna Tietze

Jakob Udby Blicher

Mikkel Petersen

Irene K Mikkelsen

Pierrick Coupé

José V Manjón

Simon Fristed Eskildsen

Automatic thalamus and hippocampus segmentation from MP2RAGE: comparison of publicly available methods and implications for DTI quantification

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

Introduction

The extensive use of magnetic resonance imaging (MRI) to investigate pathology in the brain entails identification of specific regions of interest (ROI) for quantitative analysis. Accurate manual tracing of deep brain structures, such as the thalamus and hippocampus, demands a high level of tracer expertise and preferably standardized segmentation protocols. Introducing automatic or semi-automatic techniques into postprocessing pipelines accelerate data analysis and offers reproducible and consistent decisions across datasets in large studies, which is crucial for obtaining reliable results [START_REF] Mulder | Hippocampal volume change measurement: quantitative assessment of the reproducibility of expert manual outlining and the automated methods FreeSurfer and FIRST[END_REF].

Several software solutions for automatic segmentation are publicly available. Frequently used software in clinical research include "Oxford Centre for Functional MRI of the Brain" (FMRIB) Software Library (FSL) (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/), FreeSurfer (http://surfer.nmr.mgh.harvard.edu), and Statistical Parametric Mapping (SPM) (http://www.fil.ion.ucl.ac.uk).

The segmentation techniques applied in FSL, FreeSurfer and SPM are model-based methods. In highly variable data, such as MRI of the brain, it may be difficult for the segmentation tools to model the ROIs with sufficient accuracy, even when the techniques are trained on representative datasets. To address this, multiatlas label fusion has been suggested and has demonstrated excellent segmentation abilities [START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF] [3] [START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF]. Label fusion relies on a representative image library with corresponding validated structure segmentations (atlases). Recently, multi-atlas segmentation techniques, such as patch-based segmentation, have become popular [5] [6]. Patch-based methods have the advantage of requiring a smaller training library compared to regular label fusion and is therefore relatively easy to implement in a local setting [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] [START_REF] Eskildsen | BEaST: Brain extraction based on nonlocal segmentation technique[END_REF]. Even though these improvements in segmentation algorithms have demonstrated highly accurate morphometric results 1 , most of the novel approaches are still not publicly available and therefore less used in clinical research. Moreover, the impact of segmentation accuracy on quantification of parameters from other imaging modalities, such as diffusion and perfusion MRI, is not well studied.

Quantitative diffusion tensor imaging (DTI) is widely used to investigate microstructural changes in tissue. In diseases that cause subtle microstructural changes, such as mild traumatic brain injury (mTBI), there is a need for sensitive biomarkers in clinically relevant areas of the brain. Thalamus and hippocampus are two deep brain structures, where previous DTI studies have shown microstructural changes linked to cognitive impairment [START_REF] Falangola | Non-Gaussian diffusion MRI assessment of brain microstructure in mild cognitive impairment and Alzheimer's disease[END_REF] [START_REF] Mitchell | Advances in Understanding Mechanisms of Thalamic Relays in Cognition and Behavior[END_REF], stress [START_REF] Vestergaard-Poulsen | Diffusion-weighted MRI and quantitative biophysical modeling of hippocampal neurite loss in chronic stress[END_REF] and headache [START_REF] Granziera | Structural abnormalities in the thalamus of migraineurs with aura: a multiparametric study at 3 T[END_REF]. Segmentation directly on the DTI maps is prone to inconsistency and bias, as DTI provides limited anatomical information. Unbiased and automatic studies rely on accurate T1-weighted (T1w) segmentation and co-registration for obtaining quantitative measurements within relevant brain regions. Thus, it is highly relevant to investigate the impact of automatic segmentation accuracy on these quantitative measures.

Patch-based segmentation methods [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] perform well on conventional T1w images, such as Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) [START_REF] Coupé | NeuroImage Simultaneous segmentation and grading of anatomical structures for patient ' s classi fi cation : Application to Alzheimer ' s disease[END_REF]. To the best of our knowledge, the accuracy of different automated segmentation methods has not yet been compared using T1w images from the recently proposed MP2RAGE sequence, which significantly reduces the intensity bias and provides superior grey matter (GM) to white matter (WM) contrast [START_REF] Marques | MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field[END_REF].

In this study, we compared the performance of a multi-atlas, patch-based segmentation method, as implemented in the online software platform volBrain (with two different training libraries), to three widely ap-plied methods implemented in FSL, FreeSurfer and SPM. We used manual segmentation as the gold standard and measured the segmentation accuracy of thalamus and hippocampus when imaged with MP2RAGE. Additionally, we applied the segmented masks of thalamus and hippocampus on co-registered fractional anisotropy (FA) and mean diffusivity (MD) maps for the purpose of evaluating the effect on the quantification of these diffusivity metrics.

Material and Methods

Participants

Twenty-two healthy subjects were included in the study (age range 19-40 years, 12 females). MP2RAGE images were acquired from all subjects, and DTI images in 10 subjects. All subjects were scanned on a Siemens Magnetom Skyra 3T MRI system with a 32-channel head coil. MP2RAGE parameters were acquired with TR=5s, TI 1 =0.7s, TI 2 =2.5s, α 1 =4°, α 2 =5°, a 3D sequence imaged at isotropic 1 mm 3 resolution (acquisition matrix: 240x256, 176 sagittal slices), Turbo factor of 176 as defined by Marques et al. [START_REF] Marques | MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field[END_REF]. DTI was acquired with 32 directions, 5 B0 acquisitions, TR=10.9s, TE=0.079, TI= 2.1s, imaged at isotropic 2.3 mm 3 resolution (acquisition matrix: 96x96, 38 axial slices), and inversion recovery based CSF suppression to reduce partial volume effects. Fig. 1 provides an overview of the methods and comparisons.

Pre-processing

MP2RAGE images were calculated as the ratio of two MPRAGE images, acquired with different inversion times [START_REF] Marques | MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field[END_REF]. As reported by Fujimoto and colleagues [START_REF] Fujimoto | Quantitative comparison of cortical surface reconstructions from MP2RAGE and multiecho MPRAGE data at 3 and 7T[END_REF], the amplified background noise in MP2RAGE images may introduce errors in FreeSurfer segmentations. To deal with this, we used a slightly different approach than Fujimoto and colleagues. We classified the intensities of the proton density-weighted image acquired during the second inversion recovery into four classes using a Fuzzy C-means algorithm [START_REF] Dudo | Pattern Classification, 2nd Editio[END_REF]. This robustly separates the background (1 class) from foreground (3 classes). To regain a "natural" background noise, we added the background noise of the second inversion recovery to the combined (flat) image. The modified image was used as input to all segmentation pipelines and no errors were detected. Despite inherent intensity normalization of the MP2RAGE images, all segmentation tools were run with intensity nonuniformity correction, as this was the default setting.

Diffusion data was eddy current and motion corrected using FSL, and EPI distortion correction was performed with ExploreDTI toolbox [START_REF] Leemans | ExploreDTI: a graphical toolbox for processing, analyzing, and visualizing diffusion MR data[END_REF]. During pre-processing the diffusion data was resampled to the space of the undistorted MP2RAGE image and FA and MD values were subsequently calculated with the Ex-ploreDTI pipeline.

Manual segmentation of ROIs

Thalamus and hippocampus from the 22 MP2RAGE images were manually segmented by an experienced neuroradiologist (EN: seven years of experience in neuroradiology) and a trained assistant (TA) using ITK-SNAP (www.itk-snap.org) [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF]. The thalami were first manually traced by EN in the axial plane using anatomical landmarks. Next an initial training phase of TA using the protocol outlined by Power et al. [START_REF] Power | Validation of a protocol for manual segmentation of the thalamus on magnetic resonance imaging scans[END_REF] and supervised by EN was established. TA then adjusted the thalami in all three principal planes according the protocol. Bilateral thalamus segmentation took 35-40 minutes per subject. The hippocampi were manually traced according to the EADC-ADNI segmentation protocol [START_REF] Boccardi | Delphi definition of the EADC-ADNI Harmonized Protocol for hippocampal segmentation on magnetic resonance[END_REF] by TA, initially supervised by EN. Segmentation of both hippocampi took 25-35 minutes per subject. As done in the EADC-ADNI protocol [START_REF] Boccardi | Delphi definition of the EADC-ADNI Harmonized Protocol for hippocampal segmentation on magnetic resonance[END_REF], all manual segmentations were performed in MNI space, where similar appearance of the nuclei is expected to improve tracing consistency and speed when using the segmentation protocols. The final segmentations were transformed back to scanner native space for comparison. Manual intra-operator reliability on hippocampus was tested 10 months after initial segmentations by TA on the ten subjects with DTI. A one hour training session reading the protocol and training on a separate subject was performed before segmentation of the ten subjects was carried out.

Automatic segmentation

The following provides a brief overview of the four processing methods volBrain, FSL, FreeSurfer and SPM with the applied settings.

volBrain: volBrain (http://volbrain.upv.es), which is an open access platform, is based on an advanced pipeline providing automatic segmentations of several brain structures [START_REF] Manjón | volBrain: An online MRI brain volumetry system[END_REF]. The version of volBrain used in the current setting involves an updated version of a recent patch-based method [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF]. We tested the segmentation method using two different libraries: 1) the default volBrain library (external) consisting of 50 conventional T1w images (MPRAGE and SPGR), and 2) our own manually segmented library (local) of 22 MP2RAGE images in a leave-one-out fashion. In both cases, the images were flipped across the mid-sagittal plane to artificially increase the library size as done in related work [START_REF] Eskildsen | BEaST: Brain extraction based on nonlocal segmentation technique[END_REF].

FSL: FSL is freely available for download online (http://fsl.fmrib.ox.ac.uk/fsl/ fslwiki/). The FMRIB's Integrated Registration & Segmentation Tool (FIRST), from FSL v5.0, was used to segment subcortical structures [START_REF] Patenaude | A Bayesian model of shape and appearance for subcortical brain segmentation[END_REF]. FSL FIRST is a model-based segmentation tool that uses training data from manually segmented images. In the current pipeline, we used the default setting of FIRST, which applies empirically optimized settings (numbers of modes and shape/boundary correction) for each structure (see http://fsl.fmrib.ox.ac.uk/fsl/fsl-4.1.9/first/index.html).

FreeSurfer: FreeSurfer image analysis suite (version 5.3), which is documented and freely available for download online (http://surfer.nmr.mgh.harvard.edu), was used in a default mode in the current setting. The segmentation maps are created using spatial intensity gradients across tissue classes and are therefore not simply reliant on absolute signal intensity [START_REF] Fischl | Whole brain segmentation: Automated labeling of neuroanatomical structures in the human brain[END_REF]. SPM: SPM is a MATLAB-based (MathWorks Inc.) freely available software, which can be downloaded online (http://www.fil.ion.ucl.ac.uk). Segmentation was performed with SPM12 and MATLAB R2015b by combining the unified segmentation tool with the neuroinformatics SPM template [START_REF] Ashburner | Unified segmentation[END_REF], which consists of multiple structures of the brain segmented in Montreal Neurological Institute (MNI) space.

Data and statistical analysis

The segmentations obtained from the four automatic methods were compared to the manual segmentations using volume, Dice similarity index (DSI), false positive rate (FPR), false negative rate (FNR), and Hausdorff distance estimated bilaterally. We report mean (M) and standard deviations (SD) and visualize data with boxplots.

DSI is defined as:

2C A+B and is the quotient of similarity ranging from zero to one. A and B are the number of voxels in segmentation A and segmentation B, respectively, and C is the number of voxels shared by the two segmentations. FPR and FNR were calculated respectively as number of false positive and false negative voxels as percentage of the total manually segmented number of voxels. Hausdorff distance, h, indicates the maximum distance error and is defined as the maximum distance, d, from the surface of segmentation A to the nearest point in the surface of segmentation B: h(A,B)= max aϵA {min bϵB {d(a,b)}}. Intra-rater reliability, volume, DSI, FPR, FNR and Hausdorff distance were analysed with two-way ANOVA and post-estimation was carried out, with a significance level of 0.05. Finally, FA and MD values were analysed using two-way ANOVA. Post-hoc analyses of FA and MD were carried out with a primary analysis of the automatic segmentations against the manual segmentation and secondary between the automatic methods. Results are presented at a significance level of 0.05 and in addition diffusivity results were reported with correction for multiple comparisons (60 test on diffusion metrics were carried out, which yields a Bonferroni corrected threshold of p=0.0008).

Results

Intra-operator reliability of manual segmentation

The 10 months intra-operator reliability test of hippocampus manual segmentation resulted in a mean volume difference of 3.1% (SD=4.9%), which were not significantly different (p>0.05). Mean DSI was 0.913 (SD=0.010) and mean FPR and FNR were respectively 10.4% (SD=3.4%) and 7.3% (SD=1.6%). Intraoperator Hausdorff distances ranged from 2.2 to 4.9 mm. For DTI metrics the overall model was significantly different for both FA (p<0.001) and MD (p<0.001). Post-estimation revealed a FA mean difference of 0.003 (SD=0.002) which was significantly different (p=0.003). MD obtained a mean difference of 0.006x10 -3 (SD=0.006x10 -3 ), which was significantly different (p=0.018). If Bonferroni corrected there is no significant difference between between the manual segmentations.

Thalamus and hippocampus volumes

Fig. 2 shows the volumes of thalamus and hippocampus for each of the segmentation methods. Overall the model was significantly different in both ROIs (p<0.001). There was no significant difference (p>0.05) in manual versus volBrain local, manual versus volBrain external and volBrain local versus volBrain external in thalamus, but all other comparisons for thalamus were significantly different (p<0.05). The hippocampus segmentations showed significantly higher volumes of volBrain external, FSL, FreeSurfer and SPM compared to the manual and volBrain local, and only FSL versus FreeSurfer and volBrain external versus SPM were not significantly different (p>0.05) from each other.

Manual versus automatic segmentation

Comparison of manual and automatic segmentation methods showed a substantial variation in DSI across the methods (see Fig. 3) and the overall model was significantly different (p<0.001) for both thalamus and hippocampus DSI, FPR and FNR. To maintain overview, only non-significant (P>0.05) p-values are marked in Fig. 3. All other p-values are significant (p≤0.001). DSI of the thalamus was significantly higher for volBrain local (M=0.913, SD=0.014) and volBrain external (M=0.868, SD=0.024) compared to FSL (M=0.806, SD=0.034), FreeSurfer (M=0.798, SD=0.049) and SPM (M=0.787, SD=0.031). FreeSurfer was not significantly different from FSL or SPM. FPR in the thalamus when segmented with FSL (M=41%) and SPM (M=42%) was significantly higher than the other segmentation methods. Over-segmentations are exemplified in Figs. 5A and5A, where the significantly lower FPR of volBrain local (M=9%) and external (M=14%) also can be observed. FreeSurfer FPR was significantly higher than volBrain and significantly lower than FSL and SPM. The mean FNR of the four methods ranged from 5% to 13% all being significantly different, except volBrain local versus SPM, volBrain external versus FreeSurfer, and FSL versus SPM

The DSI of the hippocampus demonstrated significantly different values between all methods, with volBrain local (M=0.892, SD=0.016) showing the best performance, followed by volBrain external (M=0.859, SD=0.014), FSL (M=0.808, SD=0.017), FreeSurfer (M=0.771, SD=0.023) and SPM (M=0.735, SD=0.038). A similar pattern was observed for FPRs, with volBrain local performing best (M=9%) followed by volBrain external (M=26%), FSL (M=36%), SPM (M=40%) and FreeSurfer performing worst (M=41%). FSL versus SPM and FreeSurfer versus SPM were the only methods which were not significantly different in FPR. Mean FNR ranged from 5% to 19% and all methods were significantly different, except volBrain local versus FreeSurfer.

In terms of Hausdorff distance, the overall model was significantly different in both thalamus and hippocampus (p<0.001). Fig. 4 shows the Hausdorff distances for the automatic hippocampus and thalamus segmentations with low distances indicating good performance. Post-estimation showed that all methods had significantly different Hausdorff distances (p<0.05) except volBrain external versus FSL in thalamus. The best performance was seen with volBrain local and the highest Hausdorff distances were measured with Free-Surfer in both thalamus and hippocampus.

Visual inspection of ROIs

Examples of manual segmentations and the corresponding automatic segmentations of the thalamus and hippocampus, overlaid on the T1w image and the FA map, are shown in Fig. 5 As illustrated, FreeSurfer, FSL and SPM generally over-segment the thalamus, especially the non-thalamic tissue near the border of the internal capsule (IC). volBrain external over-segments to a lesser extent and volBrain local demonstrated only subtle over-segmentation at the inferior and lateral border of the thalamus. The same pattern of over-segmentation is found in the hippocampus with more extensive over-segmentation by FSL, FreeSurfer, and SPM, but also slightly by volBrain external, compared to the manual (Fig. 5A and Fig. 5B). The over-segmentation of FSL, FreeSurfer, and SPM in the hippocampus is mainly restricted to the superior and the rostral part of the hippocampus in the transition to thalamus and fornix.

Diffusivity results -Thalamus

The model was overall significantly different in the diffusivity measurements for both FA (P<0.001) and MD (P<0.001) in thalamus. Fig. 6 shows mean FA and MD values in thalamus extracted from the six different segmentations and how MD values of all segmentations consistently change based on the segmentation method used, while FA values change less consistently. The volBrain local method provided the most accurate measurements compared to the manual segmentation.

Diffusivity results -Hippocampus

The model also provided overall significantly different results for FA (P<0.001) and MD (P<0.001) in hippocampus. Diffusivity results of the hippocampus are shown in Fig. 7. The figure illustrates the same consistent increase or decrease in MD between methods and subjects, but with different offsets and variation compared to the manual segmentation. FA showed a less consistent pattern. All automatic methods were significantly different from the manual values. When corrected for multiple comparisons, FA values for all methods stayed significantly different from the manual segmentation, except volBrain local and SPM and for MD all methods stayed significantly different except volBrain local.

Post-hoc analysis on diffusivity parameters between manual and automatic segmentation

Post-hoc analysis for thalamus and hippocampus are reported in Table 1.The post-hoc analysis of thalamus revealed that only volBrain local was not significantly different from the manual segmentation and obtained the lowest mean difference of M=-0.3% in FA and M=-0.1% in MD. The other methods obtained a higher mean difference, ranging from M=3% to M=9% in FA and M=1-3% in MD.

All methods obtained significantly different diffusivity parameters in the hippocampus when compared to the manual segmentation. The volBrain local demonstrated the most accurate result in the hippocampus, with a mean difference of M=-1% of FA and M=-0.5% of MD. If corrected for multiple comparisons volBrain local FA and MD were not significantly different form the manual and neither was the SPM result of FA.

Post-hoc analysis on diffusivity parameters between the automatic segmentation methods

Between-methods comparison revealed more variable results. For an overview see Table 3 with indication of corrected and un-corrected p-values for both thalamus and hippocampus. All methods, except FSL, Free-Surfer and SPM, were significantly different from each other, when measuring FA in the thalamus. When measuring MD in the thalamus, all five methods yielded significantly different results. For hippocampus FA measurements, only volBrain local stood out as different from all the other methods. Furthermore volBrain external was significantly different from FreeSurfer, while FSL versus FreeSurfer and SPM were also significantly different. For hippocampus MD, all methods were significantly different, except volBrain external versus FSL and FreeSurfer versus SPM.

Discussion

In this study, we evaluated the performance of a recent patch-based segmentation method [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] as implemented in volBrain [START_REF] Manjón | volBrain: An online MRI brain volumetry system[END_REF] and three widely used conventional methods as implemented in FSL [START_REF] Patenaude | A Bayesian model of shape and appearance for subcortical brain segmentation[END_REF], FreeSurfer [START_REF] Fischl | Whole brain segmentation: Automated labeling of neuroanatomical structures in the human brain[END_REF] and SPM [START_REF] Ashburner | Unified segmentation[END_REF]. Using MP2RAGE images, we tested the algorithms on two often investigated deep brain structures, the thalamus and the hippocampus. We found that the patch-based segmentation had the best overall accuracy. FreeSurfer, FSL and SPM all over-segmented the thalamus including non-thalamic tissue near the border of the IC and under-segmented in regions of the medial and lateral geniculate of the thalamus. In the segmentation of hippocampus, volBrain performed best followed by, FSL, FreeSurfer, and SPM. Moreover, we demonstrated that volBrain, based on a local library, was the only method, in which the diffusivity metrics of the thalamus did not differ significantly from the metrics obtained based on manual segmentation (Table 1). Analysis of hippocampus revealed that volBrain and SPM (although reporting low DSI) were not significantly different (Bonferroni corrected) from the manual method in terms of FA, and for MD only volBrain local was not significantly different, This demonstrates that segmentation accuracy impacts the obtained diffusivity results, and less accurate methods, such as FSL, FreesSurfer and SPM, do not produce consistent diffusivity results.

The accuracy of the patch-based segmentation method in our study is comparable to previous results on hippocampus segmentations using MPRAGE images [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] [START_REF] Tong | Segmentation of MR images via discriminative dictionary learning and sparse coding: Application to hippocampus labeling[END_REF]. A study by Patenaude et al. [START_REF] Patenaude | A Bayesian model of shape and appearance for subcortical brain segmentation[END_REF], using conventional T1w images and a leave-one-out comparison on its own library, found higher DSIs using FSL than found here. Patenaude and colleagues reported a mean DSI of 0.887 and 0.840 for the thalamus and hippocampus, respectively. This difference may reflect the importance of using coherent labelling protocols and similar imaging parameters within the template library. Patenaude et al. did, however, not reach the accuracy of the volBrain local segmentation in our study with DSI of 0.913 and 0.892 respectively. To compare the performance of the volBrain method with a training library different from MP2RAGE, we applied volBrain with an external training library consisting of MPRAGE and SPGR images. We found that volBrain still performed better than FSL, FreeSurfer and SPM (Fig. 3). The results of the intra-reliability test on hippo-campus further emphasize the advantage of automatic segmentation. We found a mean DSI of 0.915, which is consistent with previous findings by Frisoni et al. [START_REF] Frisoni | The EADC-ADNI harmonized protocol for manual hippocampal segmentation on magnetic resonance: Evidence of validity[END_REF] of DSI=0.89. This result is at an accuracy level of volBrain local. However, in contrast to manual segmentations, automatic methods are deterministic and yield consistent errors. Thus, automatic segmentation methods are more robust in a longitudinal setting. FSL, FreeSurfer and SPM over-segmented the structures with FPRs in the range of 15%-42%. This resulted in consistent inclusion of white matter in the segmented regions of thalamus and hippocampus (both grey matter structures) as qualitatively verified using FA maps (see Fig. 5A and5B). Regarding the volBrain method, no systematic over-or underestimation for thalamus was observed with neither local nor external libraries (FPR out-balanced FNR). Patches can capture texture similarities [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF], and this is perhaps why the patch-based method attains consistently high accuracy on both thalamus and hippocampus. volBrain local was unbiased for hippocampus, while volBrain external slightly over-segmented hippocampus (FPR M=26%). This was unexpected, because the two libraries were constructed using hippocampus masks segmented based on the same protocol (the EADC protocol), while the protocols differed for the thalamus libraries. For hippocampus, this may be explained by different interpretations (different operators) of the EADC-ADNI protocol in the segmentation procedure of the hippocampus, or by differences in contrast of the T1w images in the training libraries (MP2RAGE versus MPRAGE/SPGR).

The Hausdorff distance showed a stepwise increase between the manual and automatic methods, with mean values in the range 2-6 mm, the lowest being volBrain local followed by volBrain external, FSL, SPM and FreeSurfer (see Fig. 4). When considering the obtained FPR and FNR, the Hausdorff distance most likely reflects a maximum over-segmentation. However, evaluating the examples in Fig. 5A and 5B where the geniculate bodies of the thalamus are excluded (except for volBrain local), the distance may be due to undersegmentation in this specific region.

The intra-operator reliability test of manual hippocampus segmentation showed a consistent segmentation and no significant difference between volumes segmented with a time interval of 10 months. Our intraoperator DSI (M=0.913, SD=0.010) is in line with previous reports of manual hippocampus segmentation reliability (M=0.89, SD=0.01) [START_REF] Frisoni | The EADC-ADNI harmonized protocol for manual hippocampal segmentation on magnetic resonance: Evidence of validity[END_REF]. The DSI of repeated tracings reveal that manual segmentation of hippocampus has the same level of accuracy as between manual segmentation and the volBrain local method. The volBrain local method though has the advantage of being more consistent, faster and less costly when the library has been established [START_REF] Naess-Schmidt | Patch-Based Segmentation from MP2RAGE Images: Comparison to Conventional Techniques[END_REF].

The obtained segmentation accuracies are partly reflected in the analysis of the diffusivity metrics. The volBrain local method was the only method not yielding significantly different FA and MD results in the thalamus compared to the results obtained by manual segmentation with a mean difference of -1 % and -0.1 % in FA and MD respectively. The other methods yielded mean differences between 1 % and 9 % and were all significantly different in FA and MD compared to the values obtained with manual segmentation. This can be explained by the over-segmentation expanding into IC and the ventricular cerebral spinal fluid (CSF) (Fig. 5A and5B). In hippocampus the manual method was significantly different in FA compared to all methods (P<0.03). When correcting for multiple comparisons, FA values in both volBrain local and SPM and MD in volBrain local was not significantly different from manual measurements. The finding of SPM not being significantly different from the manual method, despite the inaccuracy of the SPM segmentation, can be explained by the segmentation expanding into both WM and CSF, which on average blur the FA differences as WM and CSF respectively represent higher and lower FA values. Higher FA values of volBrain external, FSL and FreeSurfer in the hippocampus can be explained by over-segmentation into areas at the transition to the thalamus and fornix. The difference between volBrain local and the other methods in the hippocampus segmentation is furthermore confirmed by the post-hoc analysis (Table 3), which shows that both volBrain external, FSL, FreeSurfer and SPM all significantly differ in FA from volBrain local estimates, but not from each other, if corrected for multiple comparisons. This is visualized in Fig. 7 by the offsets between volBrain local and volBrain external, FSL, FreeSurfer and SPM. The FA and MD results of the intra-operator reliability test in hippocampus showed a significant difference between the two segmentations. This difference is similar to that of the best automatic segmentation (volBrain local). In both cases the difference seems to be systematic (Fig. 7). However, this bias will be removed if using automatic methods in a longitudinal setting, as automatic methods are consistent and not prone to changing interpretations of the segmentation protocol. Finally, it should be noted that Bonferroni correction removed the significant differences between the manual segmentations.

Although the mean difference of FA and MD varies, all segmentation methods yielded consistent intersubject differences compared to the manual approach (Fig. 6 and7). This was most pronounced for MD results. Group comparisons may therefore relatively yield similar results, using the same method within the same study. However, in diseases and disorders with subtle structural changes where the influence of segmentation errors could blur the findings and result in reduced sensitivity, it is crucial to use the most accurate method to detect pathological changes. A study by Barbagallo et al. [26], found a significant difference in MD in the thalamus between amyotrophic lateral sclerosis (ALS) patients and controls (0.06 x 10 -3 , p=0.019) using FSL FIRST, but the FA difference of 0.01 was not significant (P>0.025). We speculate that such a result might have been significant if a more accurate segmentation method had been used. We found that volBrain local obtained the most accurate measurements compared to the manual segmentation (FA mean difference = -0.001), and all the other methods obtained a mean difference of FA higher than the 0.01 level obtained between groups in the study by Barbegello and colleagues. Although it is not directly comparable, the impact of using different methods (more or less accurate) in clinical studies should be investigated further.

The variation of our diffusivity measurements was considerably smaller compared to those reported in the study by Barbagello et al. The volBrain local method obtained SDs of FA values in the thalamus and hippocampus of respectively 0.008 and 0.007 and SDs of MD values of respectively 0.014 (x10 -3 ) and 0.017 (x10 - 3 ). In the Barbagello study, the corresponding SDs were 0.02 and 0.01 for FA and 0.05 (x10 -3 ) and 0.07 (x10 - 3 ) for MD. This could be due to the reliability of the MP2RAGE images, as pointed out in a recent study [START_REF] Okubo | MP2RAGE for deep gray matter measurement of the brain: A comparative study with MPRAGE[END_REF]. The MP2RAGE sequence is less influenced by B1 as well as M0 and T2*, improving the image contrast and sharpness, which makes it easier to discriminate between grey-and white matter structures [START_REF] Marques | MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field[END_REF]. Another reason for the higher FA and MD variation in the Barbagello study could be due to the ALS pathology.

In the present study, we executed the FSL pipeline with the empirical optimized default settings. Experimentation on optimizing the FSL modes and boundary correction adapted to the MP2RAGE images may improve the final segmentation of FSL. The influence of adapting the FSL segmentation tool, FIRST, to MP2RAGE remains to be investigated. However, Patenaude et al. [START_REF] Patenaude | A Bayesian model of shape and appearance for subcortical brain segmentation[END_REF], used an adapted and optimized setting to their T1w images (no specification of the T1 sequence) and did still not reach the DSI levels of the patch-based method in volBrain local. Furthermore it was not possible to use the same template library in volBrain and FSL, which would have been optimal for comparing the methods. Thus volBrain local has an advantage because it uses consistent training and testing data. This is similar to the Patenaude study [START_REF] Patenaude | A Bayesian model of shape and appearance for subcortical brain segmentation[END_REF], which showed good performance when using customized local settings. Patch-based methods have the advantage of requiring a relatively small library, which makes it feasible to implement and optimize locally.

State-of-the-art segmentation methods, like patch-based methods [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF], together with high image quality, as in MP2RAGE data, may lead to higher sensitivity in future studies of morphometry and of microstructural changes. It could also be of interest to apply the presented methods retrospectively to clinical studies and evaluate if statistical power and conclusions might be altered.

In conclusion, we have shown the potential of a recent, automatic patch-based segmentation method, volBrain, to provide more accurate thalamus and hippocampus segmentations in MP2RAGE images compared to conventional approaches. We have furthermore demonstrated that FA and MD values, extracted from co-registered DTI, deviated less from the reference of the manual segmentation, when using patchbased methods compared to the segmentations of FSL, FreeSurfer, and SPM. We have illustrated under-and particularly over-segmentations on T1w images and FA maps, especially for FSL, FreeSurfer and SPM. Finally, we propose that MP2RAGE images are more suitable for thalamus and hippocampus segmentation compared to conventional T1w images. 5 Examples of manual and automatic segmentations of thalamus and hippocampus presented in, A: the subject where volBrain local had the best performance, and B: where volBrain local had the worst performance. Upper two rows: thalamus in an axial view, overlaid on native T1 and coregistered FA images. Third row: 3D reconstructions of thalamus. The lower three rows contain similar visualizations for hippocampus segmentations. Left to right: manual, volBrain local, volBrain external, FSL, FreeSurfer, and SPM methods. Green areas indicate overlap between automatic methods and manual segmentation. Red indicates areas, which are included in the automatic, but not the manual method (false positives). Blue indicates areas, which are included by the manual, but not the automatic method (false negatives).

Fig. 1 Fig. 2

 12 Fig. 1 Overview. Left: Six segmentation methods, Middle: the two ROIs highlighted in red (upper: thalamus, lower: hippocampus) and right: three MRI comparisons (upper: T1 (DSI comparison), middle: FA, lower: MD)

Fig. 3

 3 Fig. 3 DSI, FPR and FNR for segmentations of the thalamus and the hippocampus using volBrain local (vBlocal) , volBrain external (vBext), FSL, FreeSurfer, and SPM compared to the manual "gold standard". Boxes indicate 25% and 75% percentile and the bold line the median. Whiskers indicate the most extreme point within 1.5 times the interquartile range. Data points outside this range are plotted individually. Horizontal bars indicate non-significant test for difference in DSI, FPR, and FNR. The remaining comparisons showed significant differences.

Fig. 4

 4 Fig.[START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF] Hausdorff distance of the automatic segmentations of hippocampus and thalamus compared to the manual segmentation. Boxes indicate 25% and 75% percentile and the bold line the median. Whiskers indicate the most extreme point within 1.5 times the interquartile range. Data points outside this range are plotted individually. The horizontal bar indicates a non-significant test for difference in Hausdorff distance. The remaining comparisons showed significant differences.

Fig.

  Fig.[START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] Examples of manual and automatic segmentations of thalamus and hippocampus presented in, A: the subject where volBrain local had the best performance, and B: where volBrain local had the worst performance. Upper two rows: thalamus in an axial view, overlaid on native T1 and coregistered FA images. Third row: 3D reconstructions of thalamus. The lower three rows contain similar visualizations for hippocampus segmentations. Left to right: manual, volBrain local, volBrain external, FSL, FreeSurfer, and SPM methods. Green areas indicate overlap between automatic methods and manual segmentation. Red indicates areas, which are included in the automatic, but not the manual method (false positives). Blue indicates areas, which are included by the manual, but not the automatic method (false negatives).

Fig. 6

 6 Fig.6FA and MD values for thalamus obtained by the six different segmentation methods.

Fig. 7

 7 Fig.7FA and MD values for hippocampus obtained by the six different segmentation methods and the second manual inter-rater segmentation (man2).

Table 1 Table 2

 12 FA and MD ( 1 x 10 -3 ) values from thalamus obtained by the six different segmentation methods. The values are expressed in: mean, standard deviation (SD), mean difference, 95% confidence interval (CI) and p-values ( * indicating significance at a 0.05 level and ** significance when Bonferroni corrected) Mean difference in FA and MD ( 1 x 10 -3 ) values of thalamus and hippocampus between the five automatic segmentation methods. Mean difference below the diagonal and p-values above ( * indicating significance at a 0.05 level and ** significance when Bonferroni corrected) 010 (-0.014; -0.006) <0.001 ** FSL 0.310 (0.008) 0.022 (0.016; 0.028) <0.001 ** 0.837 (0.015) 0.027 (0.023; 0.031) <0.001 ** FreeSurfer 0.316 (0.013) 0.029 (0.024; 0.035) <0.001 ** 0.819 (0.012) 0.007 (0.003; 0.012) <0.001 ** SPM 0.313 (0.008) -0.024 (-0.029;-0.020) <0.001 ** 0.845 (0.019) -0.035 (-0.045;-0.025) <0.001

See for example the recent MICCAI workshop on Multi-Atlas Labeling, https://masi.vuse.vanderbilt.edu/workshop2012/images/c/c8/MICCAI_2012_Workshop_v2.pdf

Acknowledgments

This work was funded in part by MINDLab UNIK initiative at Aarhus University, funded by the Danish Ministry of Science, Technology and Innovation, grant agreement number 09-065250, partly by the Spanish grant TIN2013-43457-R from the Ministerio de Economia competitividad and 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 (ANR-10-IDEX-03-02) by funding HL-DTI grant, Cluster of excellence CPU, LaBEX TRAIL (HR-DTI ANR-10-LABX-57) and the CNRS multidisciplinary project "Défi ImagIn.

Conflict of Interest

The Authors declare that they have no conflict of interest.

Ethical approval

"All procedures performed in the study involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. The study was a retrospective study. For this type of study formal consent is not required."