Mariano Tepper
email: mtepper@dc.uba.ar

Marta Mejail

Pablo Musé

Andrés Almansa
email: andres.almansa@telecom-paristech.fr

Andrés Almansa Boruvka

Meets Nearest

Neighbors 2011 Hal

Boruvka Meets Nearest Neighbors

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

Introduction

The computation of the minimum spanning tree (MST) is a classical problem in computer science. For an undirected weighted graph, it can be simply stated as finding a tree that covers all vertices, called a spanning tree, with minimum total edge cost. It is taught in every course of algorithms and data structure as an example where greedy strategies are successful and it is regarded as one of the first historical foundations of operations research.

The history of the MST problem up to 1985 was reviewed by [START_REF] Graham | On the history of the minimum spanning tree problem[END_REF]. Maybe the two most widely known algorithms to compute the MST are Prim's and Kruskal's [START_REF] Cormen | Introduction to Algorithms[END_REF]. There is a third classical algorithm by Boruvka [START_REF] Graham | On the history of the minimum spanning tree problem[END_REF] that mysteriously remained shadowed by the other two. This fact is emphasized by the fact that Boruvka's algorithm is also known as Sollin's algorithm, despite the fact that Sollin re-discovered it independently years later.

Under a restricted random-access computational model, [START_REF] Karger | A randomized linear-time algorithm to find minimum spanning trees[END_REF] published a randomized algorithm that runs on expected linear time. Up to date, [START_REF] Chazelle | A minimum spanning tree algorithm with inverse-ackermann type complexity[END_REF] is the fastest pointer-based algorithm to compute the MST. Table 1 summarizes the complexity of classical and state-of-the-art algorithms. [START_REF] Pettie | An optimal minimum spanning tree algorithm[END_REF] proposed an optimal theoretical algorithm which runs in time O(T * (m, n))

where n (respectively m) is the number of vertices (respectively edges) of the graph and T * is the minimum number of edge-weight comparisons needed to determine the solution.

The MST algorithm is particularly interesting for many data analysis tasks in computer vision and pattern recognition. A clear example is clustering, where the classical single-link hierarchical algorithm [START_REF] Jain | Algorithms for clustering data[END_REF] can be proved to be equivalent to computing the MST. In a seminal work, [START_REF] Zahn | Graph-theoretical methods for detecting and describing gestalt clusters[END_REF] studied the benefits of using the MST for clustering, which were lately checked in psychophysical experiments by [START_REF] Dry | The perceptual organization of point constellations[END_REF]. More recently, the MST received much attention maybe due to the growth in the size of clustering datasets [e.g. [START_REF] Carreira-Perpiñán | Proximity graphs for clustering and manifold learning[END_REF]Zemel, 2005, Felzenszwalb and[START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. The approximate MST (AMST), suboptimal but faster, also received attention for the same reasons [START_REF] Lai | Approximate minimum spanning tree clustering in high-dimensional space[END_REF].

We now slightly change the definition of the problem to a form more suitable for feature sets analysis (e.g. clustering). Then d and the pair (M, d) are said to be a metric on M and a metric space, respectively.

Definition 2 Given a metric space (M, d) and feature set X ⊆ M, the MST of X is defined as the MST of the weighted undirected graph G = (V, E) where each v i ∈ V is identified with a feature x i ∈ X, E = V ×V and the graph's weighting function ω :

E → R is defined as ω((v i , v j )) = d(x i , x j ). (1)
In other terms, the problem remains the same but the graph is now complete by definition. In such a context, and as the feature set gets larger, all the previous algorithms are worthless (as explained later in Section 1.1).

The problem is classically addressed by using metric spaces with exploitable specific characteristics, i.e. the Euclidean MST is contained in the Delaunay triangulation of X [START_REF] Eddy | Approximate single linkage cluster analysis of large data sets in high-dimensional spaces[END_REF]. Recent work has aimed at building an AMST [START_REF] Lai | Approximate minimum spanning tree clustering in high-dimensional space[END_REF] through a clever use of space-filling curves. The fractal nature of such curves imposes the use of a scale parameter, which is not easy to set automatically. [START_REF] Bentley | Fast algorithms for constructing minimal spanning trees in coordinate spaces[END_REF] and [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF] addressed the problem of using nearest neighbors search structures to compute the MST. The approach proved successful; moreover, using such structures allows in addition to compute the AMST in a natural and straightforward way. Both works are outdated and a revision in the light of novel nearest neighbors techniques and increasing computational power is much needed. More recently, [START_REF] Leibe | Efficient clustering and matching for object class recognition[END_REF] also use nearest neighbors techniques for hierarchical clustering using the average-link criterion. Although they improved the method's performance, their algorithm is not suitable for extremely large datasets.

In this work we address the MST problem without computing all distances in E. We build on Boruvka's approach, summarized in Algorithm 1, by an appropriate use of nearest neighbors search techniques.

Algorithm 1 Compute the minimum spanning tree of X using metric d

Require: X = / 0 Ensure: T = (X, E T ) is the minimum spanning tree of X using metric d

1:

E T ← / 0 2: while |E T | < |X| -1 do 3: E ′ ← / 0 4: for each connected component C of T do 5: (u m , v m ) ← arg min u∈C, v / ∈C d(u, v) 6: δ m ← d(u m , v m ) 7: E ′ ← E ′ ∪ {(u m , v m , δ m )} 8:
end for 9:

while E ′ = / 0 do 10:

(u m , v m , δ m ) ← arg min (u,v,δ )∈E ′ d 11: E ′ ← E ′ {(u m , v m , δ m )} 12: if E T ∪ {(u m , v m , d m )} does not contain cycles then 13: E T ← E T ∪ {(u m , v m , d m )} 14:
end if 15:

end while 16: end while

Complexity

Let us define n = |X|. The tasks to be performed are:

• creation of n singletons in O(n) time,
• n -1 unions of sets and • n 2 operations to find whether two nodes are in the same subtree or not.

The disjoint-sets data structure [START_REF] Tarjan | Worst-case analysis of set union algorithms[END_REF] gives extremely efficient operations on disjoint sets. The amortized complexity then is O(n 2 α(n 2 , n)) per iteration, where alpha is the extremely slow growing inverse Ackermann function, i.e. in practice α(n 2 , n) ≤ 4.

The number of connected components will be reduced by at least a factor of 2 in each iteration. At most log n iterations must be performed. Therefore, the complexity is O(n 2 α(n 2 , n) log n) which can be approximated without fear by O(n 2 log n).

The previous paragraphs contain an implicit fact: one must perform all n(n -1)/2 distance computations and thus a double-sided problem appears: in space : storing all n(n -1)/2 results for n ≥ 10 5 is prohibitive, in time : even if results are not stored, for n ≥ 10 5 the overall running-time is also prohibitive.

Keep in mind that, in modern pattern recognition applications, feature sets of 10 5 points or more are becoming common [START_REF] Frank | UCI machine learning repository[END_REF]. The rest of the paper is structured as follows. In Section 2 we propose a general approach to compute the MST using nearest neighbors search structures. Section 3 deals with search structures and in particular with a slight modification needed to compute the MST. Section 4 shows empirical results of the proposed approach on a synthetic dataset and Section 5 shows results on real image segmentation examples. Finally, some final remarks and future work are presented in Section 6.

A Nearest Neighbors Approach

Let us assume that we have a function NN d (A, b) that returns the nearest neighbor a ∈ A of b using metric d. We will discuss such functions and their implementation in Section 3.

The term arg min u∈C, v / ∈C d(u, v) in line 5 of Algorithm 1 can be straightforwardly expressed in terms of finding the nearest neighbor in the set V C:

u m = arg min u∈C d(u, NN d (V C, u)), ( 2
)
v m = NN d (V C, u m ). (3)
Let us define a constraint function ρ : X → {0, 1}. We propose to modify the function NN d (A, b) by adding an additional constraint ρ on the returned element. We denote it by NN d,ρ (A, b). In summary, it returns the nearest neighbor a ∈ A of b using metric d such that ρ(a) = 1. By setting

ρ(v) = (v / ∈ C) (4) we have NN d (V C, u) = NN d,ρ (V, u). ( 5
)
This kind of problem is sometimes referred to as Foreign Nearest Neighbors in the literature.

We are sure that the desired node v m is among the k nearest neighbors of u where k = |C| + 1. Therefore in the worst case, using a naive approach, NN d,ρ amounts to perform a k-nearest neighbors search and then a simple check among them by using ρ. Note that k is a dynamic (growing) quantity and it is not possible to fix it in advance. The problem is of a different nature than finding the MST in a constrained degree graph.

Rewriting Algorithm 1 in terms of these new elements results in Algorithm 2. Our work is similar to [START_REF] Bentley | Fast algorithms for constructing minimal spanning trees in coordinate spaces[END_REF]. They showed the pertinence of using nearest neighbors search structures to compute the MST in a Kruskal-like algorithm. Although they showed that the use of nearest neighbors search structures was fruitful for the computation of MSTs, their work received little attention.

Algorithm 2 Compute the minimum spanning tree of feature set X with metric d using nearest neighbors structures.

Require: X = / 0 Ensure: T = (V, E T ) is the minimum spanning tree of G using d

1:

E T ← / 0 2: while |E T | < |V | -1 do 3: E ′ ← / 0 4: for all connected components C of T do 5: u m ← arg min u∈C (u, NN d,ρ (V, u)) 6: v m ← NN d,ρ (V, u m ) 7: δ m ← d(u m , v m ) 8: E ′ ← E ′ ∪ {(u m , v m , δ m )} 9:
end for 10:

while

E ′ = / 0 do 11: (u m , v m , δ m ) ← arg min (u,v,δ )∈E ′ d 12: E ′ ← E ′ {(u m , v m , δ m )} 13: if E T ∪ {(u m , v m , d m )} does not contain cycles then 14: E T ← E T ∪ {(u m , v m , d m )} 15:
end if 16:

end while 17: end while Beside using of nearest neighbors, Bentley and Friedman also use priority queues to prune the number of nearest neighbors searches performed during the algorithm. First let us explain that Kruskal's algorithm is greedy: it creates a forest (i.e. a set of trees) where each isolated edge is a tree and gradually merges these trees by adding the smallest edge whose endpoints lie on different trees. They propose to store the nodes of the partial (i.e. already computed) forest, along with their foreign nearest neighbors, in a single and global priority queue where the priority of a node is the inverse of the distance to its foreign nearest neighbor. The use of a priority queue is indeed interesting in this context, as the next edge to add to the MST is at the top of the priority queue. The top of the queue is removed and the top-priority foreign nearest neighbors is added to the MST. This node is also added to the queue, after computing its foreign nearest neighbors. Additionally, the priority queue must be updated, since disjoint connected components are merged and some foreign nearest neighbors might not be foreigners anymore.

Note that it may not be necessary to update the entire priority queue. This is because the current priority of each of these nodes (the priority before the insertion in the MST) serves as an upper bound of its real priority (the priority after the insertion in the MST). The real priority of a node needs only to be computed when its current priority is on the top of the queue. Thus, by using a global priority queue Bentley and Friedman were able to speed up a Kruskal-like algorithm.

We already stated our interest in building on Boruvka's algorithm, hence a global priority queue is not suitable in this case. Alternatively, we propose to use several priority queues, one for each connected component in the partial MST. Each queue, only holds the nodes in its respective connected component, and their foreign nearest neighbor. After an insertion in the MST, two connected components are merged and their priority queues are also merged.

Note that the space complexity is still O(n). In the first iteration, there are n queues, each of length 1. In the second iteration there are roughly n/2 queues, each of length 2, and so on.

Algorithm 3 is the result of this modification. The operators top and pop return the top-priority element in the queue and remove it, respectively.

Approximate MST

We stated that our approach allows to compute approximate MSTs. Indeed, if we simply relax the search by finding the approximate nearest neighbors we end up with an approximate minimum spanning tree algorithm. Approximate nearest neighbors queries are much faster than exact ones, specially in high-dimensional spaces.

Typically, ANN d (X, u, η) ensures that, if the true nearest neighbor is at distance δ , the approximate nearest neighbor is at a distance lower than δ (1 + η). Note that AMSTs can also be obtained by using a probability bound on the nearest neighbor distance [START_REF] Toyama | Probably correct k-nearest neighbor search in high dimensions[END_REF]. [START_REF] Lai | Approximate minimum spanning tree clustering in high-dimensional space[END_REF] have previously studied AMSTs. Their approximation is obtained by using space-filling structures, i.e. Hilbert curves. Their work differs from ours in two central points. First, our algorithm allows to combine MSTs and AMSTs in a single framework, in which the only difference between them is a relaxation parameter. Their work is restricted to AMSTs. Second, Hilbert curves are fractal and the spacefilling accuracy follows an exponential scale. It relies on a scale parameter that has a non-intuitive meaning and which is difficult to chooseIt is not straightforward to set automatically a suitable scale for a given point set configuration. The relaxation parameter in our method has a clear interpretation and it is easy to monitor its effect.

Nearest Neighbors Search Structures

The problem of efficiently finding nearest neighbors has received much attention, as it is in the core of a great variety of problems in pattern recognition and classification. To name a few, non-parametric density estimation [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF], non-linear manifold learning [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and descriptors matching [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF].

A plethora of nearest neighbors search structures have been proposed over the years:

• Kd-trees and randomized Kd-trees [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF], [START_REF] Silpa-Anan | Optimised kd-trees for fast image descriptor matching[END_REF],

Algorithm 3 Compute the minimum spanning tree of X with metric d using nearest neighbors structures and multiple priority queues.

Require:

X = / 0 Ensure: T = (V, E T ) is the minimum spanning tree of G using d 1: for all v ∈ V do {each node is a connected component} 2: v m ← NN d,ρ (V, v) 3: d m ← d(v, v m ) 4:
add (v, v m ) to the empty priority queue Q {v} with priority d m 5: end for 6:

E T ← / 0 7: while |E T | < |V | -1 do 8: E ′ ← / 0 9: for all connected components C of T do 10: (v, v m ) ← top(Q C ) 11: while v m ∈ C do {not a foreign nearest neighbors} 12: pop(Q C ) 13: u m ← (v, NN d,ρ (V, v)) 14: δ m ← d(v, u m ) 15: add (v, u m ) to Q C with priority δ m 16: (v, v m ) ← top(Q C ) 17:
end while 18:

δ m ← d(v, v m ) 19: E ′ ← E ′ ∪ {(v, v m , δ m )} 20: end for 21: for all (u, v) ∈ E ′ do 22: C ← connected component of T containing u 23: C ′ ← connected component of T containing v 24: merge Q C and Q C ′ into Q C∪C ′ 25:
end for 26:

while • VP-trees and MVP-trees [START_REF] Yianilos | Data structures and algorithms for nearest neighbor search in general metric spaces[END_REF], [START_REF] Bozkaya | Distance-based indexing for high-dimensional metric spaces[END_REF]],

E ′ = / 0 do 27: (u m , v m , δ m ) ← arg min (u,v,δ )∈E ′ d 28: E ′ ← E ′ {(u m , v m , δ m )} 29: if E T ∪ {(u m , v m , d m )} does not contain cycles then 30: E T ← E T ∪ {(u m , v m ,
• M-trees and their extensions [START_REF] Ciaccia | Bulk loading the m-tree[END_REF]], [START_REF] Skopal | Revisiting m-tree building principles[END_REF], [START_REF] Skopal | Nearest neighbours search using the pm-tree[END_REF], [START_REF] Zhou | M + -tree: a new dynamical multidimensional index for metric spaces[END_REF]],

• hashing-based methods [START_REF] Andoni | Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions[END_REF],

• hierarchical k-means trees [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF] • list-of-clusters [START_REF] Chavez | An effective clustering algorithm to index high dimensional metric spaces[END_REF], [START_REF] Chávez | A compact space decomposition for effective metric indexing[END_REF].

The above list is, of course, non-exhaustive.

All metric search structures, except hashing-based methods, exploit, one way or another, the metric's triangular inequality to reduce the number of distance computations to be performed in a search. Hashing-based methods aim at finding hash functions which ensure that the probability of collision is much higher for objects that are close to each other than for those that are far apart.

List-of-clusters

Now we turn our attention to the list-of-clusters structure [START_REF] Chavez | An effective clustering algorithm to index high dimensional metric spaces[END_REF]Navarro, 2000, Chávez and[START_REF] Chávez | A compact space decomposition for effective metric indexing[END_REF]. It is reported to be very efficient and resistant to the intrinsic dimensionality of the data set. Our experiments, shown in Table 2, support this claim. It can also be implemented in primary and in secondary memory. Furthermore, it has the advantage of being easy to implement and understand, as it does not involve complex data structures. For this reason we choose list-of-clusters to explain the modifications needed to find foreign nearest neighbors in depth.

The core of the list-of-clusters construction algorithm is the choice of a center c ∈ X and a radius r c . The possible choices for them will be discussed later. Let us define

• internal elements:

I X,c,r c de f = {x ∈ X -{c}, d(c, x) ≤ r c },
• external elements:

E X,c,r c de f = {x ∈ X, d(c, x) > r c }.
A bucket is defined as the tuple (c, r c , I X,c,r c ). The process is repeated inside E X,c,r c recursively, producing at the end a list of buckets. This procedure is depicted in a recursion-free manner, in Algorithm 4. List-of-clusters, where only external elements are split, can be seen as a degenerated VP-tree, where external and internal elements are split.

Algorithm 4 Build a list-of-clusters from X using metric d

Require: X = / 0 Ensure: L X,d is the list of clusters of X L X,d ← / 0 while X = / 0 do select a center c ∈ X select a radius r c I ← {x ∈ X -{c}, d(c, x) ≤ r c } X ← X -I -{c} L X,d ← L X,d : (c, r c , I) end while
There are two decisions to make at the core of the building algorithm: the selection of the center c and the radius r c [START_REF] Chávez | A compact space decomposition for effective metric indexing[END_REF]. At iteration i, for selecting the center c i , we chose the element farthest to c i-1 in the remaining set. The objective of such a choice is to minimize the overlap between regions. By using partitions of fixed size, the radius r c is then easily deduced. Now we focus our attention to search itself. The original search algorithm in the list-of-clusters structure iterates through the list of buckets and performs exhaustive searches only when needed (determined by triangular inequalities). An exhaustive search occurs within a bucket's internal elements. It can be written in the following terms

x m ← arg min x∈I d(x, q) (6) d m ← d(x m , q). ( 7
)
The standard list-of-clusters search algorithm [START_REF] Chávez | A compact space decomposition for effective metric indexing[END_REF] is shown in Algorithm 5, with the introduction of the constraint ρ in the exhaustive search (lines 10 and 9). The operators head and tail return the first element in the list and remove it, respectively. Adding constraint ρ to any other search structure is similar. Algorithm 5 can be very easily modified to find approximate nearest neighbors. It is sufficient to rewrite lines 4 and 8 as follows

d c /r ≤ 1 + η (8) d c /(r c + r) ≤ 1 + η, (9)
where η is a relaxation parameter. Note that more complex nearest neighbors formulations can be used, such as probabilistic bounds [START_REF] Toyama | Probably correct k-nearest neighbor search in high dimensions[END_REF].

Algorithm 5 Search for the nearest neighbor p of q in the list-of-clusters L X,d with initial radius r and restriction ρ

Require: L X,d is not empty and r > 0 Ensure: p ∈ X is the nearest neighbor of q if ∃x ∈ X, d(x, q) ≤ r 1: repeat 2:

(c, r c , I) ← head(L X,d ) 3: d c ← d(c, q) 4: if d c ≤ r then 5: p ← c 6: r ← d c 7: end if 8: if d c ≤ r c + r then 9:
x m ← arg min

x∈I ρ(x)=1 d(x, q) 10: d m ← d(x, x m ) 11: if d m ≤ r then 12: p ← x m 13: r ← d m 14: end if 15: end if 16: L X,d ← tail(L X,d ) 17: until L X,d is not empty or d c ≤ r c -r

Experimental Results

As distance computations are the dominating speed factor, we measure performance and complexity as a function of them. We sample points from a uniform distribution in the unit hyper-cube. We tested with four different dimensionalities R 2 , R 5 , R 10 and R 20 . We compared the following methods:

Bvka: all distances are precomputed and stored in memory and then Algorithm 1 is performed. Bvka-O: Algorithm 2 where an online linear search is used to compute nearest neighbors. Bvka-LOC: Algorithm 2 where nearest neighbors are computed online by using the list-of-clusters search structure. Bvka-PQ-LOC: Algorithm 3 where nearest neighbors are computed online by using list-of-clusters search structure. Bvka-A η: Bvka-PQ-LOC modified to compute the AMST by using approximate nearest neighbors (see Equations 8 and 9) where η is the relaxation parameter.

A summary of these methods is presented in Table 3. Note that the reduced memory complexity of the algorithm warranties that we will be able to treat large datasets without "out of memory" issues.

Comparisons were made for relatively small feature sets (|X| ≤ 10 4 ) to be able to compare with a classical MST implementation. A summary of our results is shown in Figure 1. At a first approach and as expected, all algorithms have a polynomial dependency on the dataset size, since they are linear in a log-log plot. In fact, we are not aiming at improving the worst case complexity of Boruvka's algorithm, but at improving its expected performance.

Method

Solution Distances Space Search computed stored complexity speed 1 in a log-log scale. In low dimensions, Bvka-LOC is better than any classical algorithm while Bvka-PQ-LOC resists better the dimensionality increase.

Bvka MST n(n -1)/2 all O(n 2 ) - Bvka-O MST O(n 2 log n) none O(1) linear Bvka-LOC MST O(sn log n) none O(n) sub-linear Bvka-PQ-LOC MST O(sn log n) n -1 O(n) sub-linear Bvka-A η AMST O(sn log n) n -1 O(n) sub-linear
Results are indeed encouraging. Our method exhibits a very strong performance improvement in low dimensions, see the top row in Figure 1. Bvka-LOC and Bvka-PQ-LOC in both cases outperforms Bvka several orders of magnitude.

We can also notice a strong performance degradation of Bvka-LOC with the increase of dimensionality, see the bottom row in Figure 1. The only cause is the nearest neighbors search structure. It is a well known fact that the performance of nearest neighbors search structures tends to become linear in high-dimensions. In any case, our method is generic: any nearest neighbor structure can be used. Another structure may provide better results in high dimensions and we plan to explore these issues in future work.

Table 4 summarizes the results from Figure 1 by analyzing the slope of the different curves. The proposed approach lowers in practice the number of distance computations needed to solve the problem. The quadratic profiles of Bvka and Bvka-O are reduced to supralinear (e.g. n 1.6 approximately) by Bvka-LOC and Bvka-PQ-LOC. As stated, the latter shows a computational cost which is less sensitive to an increase in dimensionality.

We provide a simple example of the incidence of using the AMST, shown in Figure 2. We use X uniformly distributed on the square [0, 1] 2 and Euclidean distance. Computing the MST required 9613 distance computations with our algorithm, while taking 9155, 8705 and 7840 with η = 0.1, η = 0.2, η = 0.5 respectively. There is an important improvement in performance while the number of topology changes is small. Moreover, when carefully inspected, these changes are reasonable. It is a well known fact that (even little) jitter noise in the dataset greatly affects the topology of the MST [START_REF] Carreira-Perpiñán | Proximity graphs for clustering and manifold learning[END_REF]: computing the AMST can be seen as perturbing the dataset with such a noise. Usually η is chosen to be quite small, and its use has more meaning in large and high-dimensional datasets. In our toy example, keeping η small does not introduce changes in the topology of the tree. We exaggerated η to show actual topology changes.

A performance comparison between MSTs and AMSTs is shown in Figure 3. We use X uniformly distributed in the hyper-cube [0, 1] 20 and Euclidean distance. As argued before Bvka-LOC's performance tends to Bvka-O's in high-dimensions. Bvka-A greatly improves the performance: it is 1.7 and 1.62 times faster than Bvka-O and Bvka-LOC respectively when |X| = 10 4 .

Computing the MST for |X| = 10 5 is not possible with classical algorithms on standard computers, since approximately 5 • 10 9 distances must be computed and stored. This means more than 18.6 GB if we use 32 bits to store each computed distance. Using minimum memory (less than 20 MB), we were able to compute the MST using Euclidean distance, without, explicitly nor implicitly, exploiting the nature of the Euclidean space (i.e. without relying on Delaunay triangulations). Table 5 presents the resulting running times for all considered algorithms. Again, these results can be improved, as we did not perform any tuning of the list-of-clusters.

Moreover, since parallelization is straightforward, it can be exploited to boost the performance. Every iteration in both nearest neighbors-searching cycles, lines 4 to 9 in Algorithm 2 and lines 9 to 20 in Algorithm 3, can be run in parallel since it operates on a single disjoint connected component. We estimate that at least a 20X speedup can be achieved by using GPGPU.

Finally, more efficient search algorithms can be implemented for a given nearest neighbors structure that might increase the performance of the proposed algorithms, such as the best-bin-first or an optimized depth-first [START_REF] Samet | Depth-first k-nearest neighbor finding using the maxnearestdist estimator[END_REF]. [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] proposed a clustering method and used it for image segmentation. The algorithm is graph-based, building a forest of partial MSTs.

Application to Image Segmentation

The authors argue that:

"There are several possible ways of determining which feature points to connect by edges. We connect each point to a fixed number of nearest neighbors. Another possibility is to use all the Dimensions Bvka-PQ-LOC Bvka-A 0. neighbors within some fixed distance δ . In any event, it is desirable to avoid considering all O(n 2 ) pairs of feature points." [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] Felzenszwalb and Huttenlocher also test a version of their algorithm in which the graph layout is determined by the spatial neighborhood in the image. The previous citation is true in a context where all n(n -1)/2 distances are needed to compute the real MST.

On the one hand, it is clear that pruning the complete graph yields a faster algorithm. Nowadays a small resolution image has 640 × 480 = 307200 pixels. Using algorithms available in the literature, finding the MST is intractable.

On the other hand, from a conceptual point of view, pruning is dangerous. Results can be severely affected, as discussed by [START_REF] Fowlkes | How much does globalization help segmentation?[END_REF] and [START_REF] Cour | Spectral segmentation with multiscale graph decomposition[END_REF]. The use of k-nearest neighbors graphs or fixed radius graphs also brings a new complication: both parameters have to be carefully selected to yield a connected graph. For example, if k is fixed such that the resulting graph is not fully connected, an image where all pixels have the same color would end up being segmented! Based on the concepts of Section 2, we propose a new version of Felzenszwalb and Huttenlocher's algorithm. This version drops the original requirements of feature locality, becoming global. It exploits the fact that the proposed algorithm makes tractable the problem of computing the MST for a set of 10 5 ∼ 10 6 features.

The classical choice to partition a feature set in hierarchical clustering, is to globally thresholding the hierarchy at a fixed level [START_REF] Jain | Algorithms for clustering data[END_REF]. Felzenszwalb and Huttenlocher proposed a similar partitioning criterion but thresholds are locally adaptive, thus providing more realistic and useful partitions. This algorithm can also link disjoint nodes across the hierarchy, thus building a forest of partial MSTs, each of which may not be a subtree of the global MST. To illustrate the usefulness of the proposed MST algorithm, we use the splitting criterion by Felzenszwalb and Huttenlocher without allowing for such links across the hierarchy. In general, this modification does not involve a negative impact on the results but allows working with larger datasets.

Let again X be a feature set and d a suitable metric. Let T = (X, E T ) be the minimum spanning tree of X. We say

C ∈ X is a component if T C is a (connected) subtree of T with node set C. For a component C we define max(C) = max v i ,v j ∈C (v i ,v j )∈E T d(v i , v j ) + τ(C) (10)
The function τ can take many forms, but can be simply defined as τ(C) = s/|C|, where s acts as a scale parameter [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. Then, given two components C 1 and C 2 we define

min(C 1 ,C 2 ) = min (max(C 1 ), max(C 2 )) (11)
Felzenszwalb and Huttenlocher's algorithm is non other than Kruskal's with an additional criterion to avoid merging certain connected subcomponents in the MST. Our modified version, instead of constructing a MST, builds a minimum spanning forest. Given two disjoint components C 1 and C 2 they are only merged if

min(C 1 ,C 2 ) ≥ min v i ∈C 1 v j ∈C 2 d(v i , v j ) (12)
Each tree in the forest is finally detected as a cluster. The process is described in Algorithm 6.

Algorithm 6 Compute the clustering of feature set X using metric d

Require:

X = / 0 Ensure: S is a clustered partition of X 1: compute the MST T = (X, E T ) 2: n ← |X| 3: m ← |E T | 4: ẼT ← sort E T by non-decreasing edge weight. 5: S ← {{v 1 }, . . . , {v n }} 6: while ẼT = / 0 do 7: (v i , v j ) ← head( ẼT ) 8:
let C i and C j be the sets of S containing v i and v j respectively. 9:

if

(C i = C j ) ∧ (d(v i , v j ) ≤ min(C i ,C j )) then 10: S ← S -C i -C j ∪ C i ∪C j 11:
end if 12:

ẼT ← tail( ẼT ) 13: end while

In the following experiments, the color value in RGB space for each pixel is used as a feature in R 3 and Euclidean distance is used to construct the MST. As stated above, any metric could have been used without changing the algorithm hence, in this sense, the algorithm is parametric on the metric.

For our tests we compared three different versions of Felzenszwalb and Huttenlocher's algorithm: FH-Grid: Algorithm 6 using the image grid connectivity to prune the complete graph FH-Bvka-NN: Algorithm 6 using the complete graph and the proposed algorithm to compute the MST. FH-Bvka-A η: Algorithm 6 using the complete graph and the proposed algorithm to compute the AMST.

We tested these algorithms on images from the Berkeley segmentation dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. Since their size is 321 × 481 = 154401 pixels, computing the MST with any classical (or even with Chazelle's state-of-the-art algorithm) is not possible. The results on Figures 4, 5 and 6 (first three rows) show the importance of the role of global mechanisms in image segmentation.

Human visual system is able to capture and to use global characteristics. For example, we perceive a blue sky as a whole, even when a part appears through a hole in some object. Obviously, FH-Grid can not reproduce such behavior while global methods perform well.

By using the image connectivity, artifacts are introduced (clearly visible in all examples). Although this can be corrected by using locality directly on the color space, scale parameters are left free and must be correctly tuned to obtain satisfactory results.

Another issue can bee observed in the methods obtained with FH-Grid: the number of regions is badly overestimated. Accuracy in the localization of visually perceived region frontiers comes at the price of over-segmenting the image.

As a counter part, all global methods we tested have difficulties when segmenting images where several objects share common colors. This is also natural and is a widely adopted mechanism by animals in the from of camouflage. Animals copy colors and structures from their environment to avoid being detected by predators or preys. More careful and local inspection are necessary to detect them.

The last three rows in Figures 4, 5 and 6 illustrate the stability of using FH-Bvka-A with different levels of approximation compared to the exact FH-Bvka-NN: changing η maintains the global structures in the segmentation.

Final Remarks

The dominating factor when computing the MST of a feature set X is the number of distance computations to be performed. We presented a method for computing the MST based on a clever use of nearest neighbors search structures. It has O(n 2 ) and O(n) time and space complexities respectively. However, in practice it outperforms classical algorithms for large, and low dimensional, datasets.

The same algorithm with a slight modification can also be used to compute the AMST: instead of finding nearest neighbors, one finds approximate nearest neighbors. In high-dimensional datasets, we showed the performance increase that results from using AMSTs. Moreover, in our tests the AMST, as computed, has an stable behavior.

To show the pertinence of the proposed algorithm, we use it to improve a stateof-the-art image segmentation algorithm proposed by [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. Due to the memory and performance limitations of classical MST algorithms, Felzenszwalb and Huttenlocher's method is restricted to finding segmentations based on local connections or forced to rely on arbitrary connectivity parameters. Thanks to our input, a global segmentation can be obtained, without modifying the overall algorithm and without any extra parameter.

There are three conceptual main lines for future work. The first consists on performing an experimental evaluation of nearest neighbors search structures and their incidence on the performance of the proposed algorithm. This includes the evaluation of different criteria in list-of-clusters for selecting the centers and the radii. Second, we did not explore other search algorithms [START_REF] Samet | Depth-first k-nearest neighbor finding using the maxnearestdist estimator[END_REF] which may reduce the number of distance computations per query. Finally, when using approximate MSTs, the trade-off between enhanced speed and accuracy must be explored more carefully.

Last, from the implementation point of view, the proposed algorithms can be parallelized without any reformulation. Moreover, in list-of-clusters, the exhaustive search within a bucket can be implemented using vectorial processors as the bucket size is fixed.

Definition 1

 1 Given a set M and a function d: M × M → R such that, ∀x, y ∈ M, • d(x, y) ≥ 0 (non-negativity), • d(x, y) = 0 ⇔ x = y (identity of indiscernibles), • d(x, y) = d(y, x) (symmetry) and • d(x, z) ≤ d(x, y) + d(y, z) (triangle inequality).

FIG 3 .

 3 FIG 2. Comparison of the MST (using Bvka) vs the AMST (using Bvka-A η) for several levels of relaxation η.

FIG 4 .FIG 6 .

 46 FIG 4. Effect of using local, complete and approximate complete graphs. The oversegmentation produced by FH-Grid has been corrected by all the proposed methods. Notice also that the approximate versions of FH-Bvka-A show an stable behavior.

TABLE 1

 1 Major MST algorithms, their computational model and their time complexity (n and m are the number of nodes and the number of edges, respectively).Karger's and Chazelle's are amortized complexities. All algorithms have a spatial complexity of O(m) at least.

  Running-time comparison of search structures in the size and in the dimensionality of the dataset. The dataset is composed by uniformly distributed points. The values are the average of 1000 random queries and are expressed in milliseconds.

		Size Dimensions List-of-clusters kd-tree k-means tree	Linear
			2	0.04	0.22	0.1	2.76
			3	0.03	0.32	0.17	2.64
		10 4	4 5	0.08 0.04	0.4 0.61	0.45 1.11	2.59 2.56
			10	0.19	1.72	2.84	3.23
			20	0.8	5.62	4.16	4.48
			2	0.31	1.33	0.3	47.65
			3	0.49	2.04	0.91	52.84
		10 5	4 5	0.72 0.6	3.07 4.34	2.42 5.47	50.19 49.88
			10	2.7	14.51	53.52	59.26
			20	18.39	70.68	66.07	65.21
				TABLE 2
			d m )}
	31:	end if
	32:	end while
	33: end while

TABLE 3

 3 The methods compared in this work. s stands for average number of distance operations needed to complete a nearest neighbors search. The space required by the nearest neighbors search structure is O(n).FIG 1. Comparison in the number of distance computations as |X| grows. From left to right: top row, X ⊂ R 2 and X ⊂ R 5 ; bottom row, X ⊂ R 10 and X ⊂ R 20 . The radii in the list-of-clusters were chosen such that each bucket has |X|/2 internal elements. Both scales are logarithmic.

	Method	R 2	R 5	R 10	R 20
	Bvka	2	2	2	2
	Bvka-O	2.14 2.12 2.13 2.15
	Bvka-LOC	1.58 1.66 1.92 2.15
	Bvka-PQ-LOC	1.61	1.6	1.87 2.03
		TABLE 4
	Slopes of the different curves in Figure

Original image FH-Grid FH-Bvka-PQ-LOC FH-Bvka-A 0.1 FH-Bvka-A 0.2 FIG 5. Effect of using local, complete and approximate complete graphs. The oversegmentation produced by FH-Grid has been corrected by all the proposed methods. Notice also that the approximate versions of FH-Bvka-A show an stable behavior.

Acknowledgements

The authors wish to thanks Pablo Cancela for his fruitful comments and acknowledge financial support by CNES (R&T Echantillonnage Irregulier DCT / SI / MO -2010.001.4673), FREEDOM (ANR07-JCJC-0048-01), Callisto (ANR-09-CORD-003), ECOS Sud U06E01, ARFITEC (07 MATRH) and STIC Amsud (11STIC-01 -MMVP-SCV) and the Uruguayan Agency for Research and Innovation (ANII) under grant PR-POS-2008-003.


Joris Van Der Hoeven
email: vdhoeven@lix.polytechnique.fr

Efficient root counting for analytic functions on a disk *

Keywords: Root counting, reliable computation, Graeffe method A.M.S. subject classification: 68W25, 68W30, 65G20, 30B10, 42-04

In this note, we present a variant of an algorithm by Schönhage for counting the number of zeros of a complex polynomial in a disk. Our algorithm implements a few optimizations and also applies to more general analytic functions.

Introduction

Many algorithms have been proposed for the reliable computation of zeros of complex polynomials [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF][START_REF] Gourdon | Combinatoire, algorithmique et géométrie des polynômes[END_REF][START_REF] Kravanja | Computing the zeros of analytic functions[END_REF][START_REF] Pan | Univariate polynomials: nearly optimal algorithms for numerical factorization and root-finding[END_REF], assuming multiple precision arithmetic. A slightly less ambitious problem is to count the number of zeros in a given disk; see also [START_REF] Rump | Verification methods: Rigorous results using floating-point arithmetic[END_REF]Sections 13.2 and 13.3]. In this paper, we present a variant of Schönhage's method [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF], based on iterated Graeffe transforms, but with a few advantages.

We present our algorithm in the setting of ball arithmetic [vdH09], which is our preferred variant of interval arithmetic [START_REF] Moore | Interval Analysis[END_REF][START_REF] Alefeld | Introduction to interval analysis[END_REF][START_REF] Moore | Introduction to Interval Analysis[END_REF][START_REF] Rump | Verification methods: Rigorous results using floating-point arithmetic[END_REF]. In this framework, a large part of the burden of bound computations is taken away from our shoulders and moved to the underlying arithmetic. Moreover, the algorithm naturally applies for analytic functions, which are represented by an approximating polynomial and a bound for the error. Finally, during the iterated application of Graeffe transforms, some coefficients of the polynomial become very small. Ball arithmetic allows us to move such coefficients to a global error term and reduce the degree of the polynomial, thereby speeding up the overall algorithm.

Our algorithm is presented for analytic functions whose power series expansion is given explicitly up to a given order, together with an error bound for the tail. In practice, we are often given a program which may compute the expansion (and tail bound) up to any required order. If we are expecting a k-fold multiple root, then it is a good practice to compute the expansion up to order O(z τk ) for some 1 < τ < 2: this is usually only a constant times more expensive than the computation of an expansion up to order O(z k+1 ), but greatly improves the quality of the error bounds and the probability that our root counting algorithm will be successful.

Ball arithmetic

Basic principles

Let us briefly recall the principles behind ball arithmetic, while referring to [vdH09] for details. Given a normed vector space K, we will denote by K or B(K, R) the set of closed balls with centers in K and radii in R = {x ∈ R: x 0}. Given such a ball z ∈ B(K, R), we will denote its center by cen(z) and its radius by rad(z). Conversely, given z ∈ K and r ∈ R, we will denote by z + B(r) the closed ball with center z and radius r.

A continuous operation f : K r → K is said to lift into an operation f : K r → K on balls, which is usually also denoted by f , if the inclusion property

f (x 1 , , x r ) ∈ f (x 1 , , x r ) (1)
is satisfied for any x 1 , , x r ∈ K and x 1 ∈ x 1 , , x r ∈ x r . For instance, if K is a Banach algebra, then we may take

x + y = cen(x) + cen(y) + B(rad(x) + rad(y))

xy = cen(x) -cen(y) + B(rad(x) + rad(y))

x y = cen(x) cen(y) + B(rad(x) (|cen(y)| + rad(y)) + |cen(y)| rad(x)).

Similar formulas can be given for division and elementary functions.

Floating point arithmetic

In concrete machine computations, numbers are usually approximated by floating point numbers with a finite precision. Let R ˜= R p be the set of floating point numbers at a working precision of p bits. It is customary to include the infinities ±∞ in R ˜as well. The IEEE754 standard [START_REF]ANSI-IEEE Standard[END_REF] specifies how to perform basic arithmetic with floating point numbers in a predictable way, by specifying a rounding mode R ∈ {↓, ↑, } among "down", "up" and "nearest". A multiple precision implementation of this standard is available in the Mpfr library [START_REF] Hanrot | MPFR, a C library for multipleprecision floating-point computations with exact rounding[END_REF]. Given an operation f : R r → R, we will denote by f R : R ˜r → R ˜its approximation using floating pointing arithmetic with rounding mode R. This notation extends to the case when R and R ˜are replaced by their complexifications C and

C ˜= R ˜[i]. Setting C ˜= R ˜[i],
we will denote by C ˜or B(C ˜, R ˜) the set of closed balls in C with centers in C ˜and radii in R ˜ . In this case, we will also allow for balls with an infinite radius. A continuous operation f : C r → C is again said to lift to an operation f : C ˜r → C õn balls if (1) holds for any x 1 , , x r ∈ C ˜and x 1 ∈ x 1 , , x r ∈ x r . The formulas for the ring operations may now be adapted to

x + y = cen(x) + cen(y) + B(rad(x) + ↑ rad(y) + ↑ ǫ +,x,y ) x -y = cen(x) -cen(y) + B(rad(x) + ↑ rad(y) + ↑ ǫ -,x,y ) x y = cen(x) × cen(y) + B(rad(x) × ↑ (|cen(y)| + ↑ rad(y)) + ↑ |cen(y)| × ↑ rad(x) + ↑ ǫ ×,x,y ),
where ǫ +,x,y , ǫ -,x,y and ǫ ×,x,y are reliable bounds for the rounding errors induced by the corresponding floating point operations on the centers; see [vdH09] for more details. Given z ∈ C ˜, it will be convenient to denote by ⌊z ⌋ and ⌈z ⌉ certified lower and upper bounds for |z |.

Balls of polynomials and analytic functions

Let B = B(1) ⊆ C denote the closed unit disk and S ⊆ C the closed unit circle. Let A B and A S denote the rings of analytic functions on B resp. S. We define norms on A B and A S by

f = sup z ∈B |f (z)| (f ∈ A B ) f = sup z ∈S |f (z)| (f ∈ A S ).
By the maximum modulus principle, the inclusion A B A S preserves norms. Since our criterion for root counting will be based on a variant of Rouché's theorem, we will mainly consider analytic functions on S in what follows. In order to avoid confusion with B(r), we will denote by B S (r) the closed ball of radius r in A S :

B S (r) = {f ∈ A S : f r}.
In view of section 2.1, we may then consider the space B(C[z], R) of balls with centers in C[z] and radii in R . Any such ball can be written f = P + B S (r) for P ∈ C[z] and r ∈ R . For concrete machine computations, and in a similar way as in section 2.2, we may also consider the space

B(C ˜[z], R ˜) of balls with centers in C ˜[z] and radii in R ˜ . Given f ∈ B(C ˜[z],
R ˜), we will denote by ⌈⌈f ⌉⌉ a certified upper bound for f (z) on S. If f = P 0 + + P n z n + B S (r), then we may for instance take ⌈⌈f ⌉⌉ = ⌈P 0 ⌉ + ↑ + ↑ ⌈P n ⌉ + ↑ r.

Simplification

Consider a ball f = P + B S (r) ∈ B(C ˜[z], R ˜) with n = deg P 0 and let e P be the maximum of the exponents of the coefficients of P . Let q p be maximal such that 2 e P -q r 256 (n + 1) and take q = 0 if 256 (n + 1) 2 e P < r. Denote Z q = ±{0, , 2 q -1}. Truncation of the mantissas of the coefficients of P leads to a decomposition P = P head + P tail , where

P head ∈ Z q [i][z] 2 e P +1-q ⌈⌈P tail ⌉⌉ max r 64
, 2 e P +2-p . We define the simplification σ(f ) of f by

σ(f ) = P head + B S (⌈⌈P tail + B S (r)⌉⌉)
and notice that f ⊆ σ(f ).

Efficient multiplication

Assume now that we want to multiply two balls f = P + B S (r) and g = Q + B S (s) in an efficient way. Modulo simplification, we may assume without loss of generality that f = σ(f ) and g = σ(g). For some common q p, we thus have

P ∈ Z q [i][z] 2 e P +1-q Q ∈ Z q [i][z] 2 e Q +1-q
We may multiply the integer polynomials P * = P /2 e P +1-q and Q * /2 e Q +1-q using any fast classical algorithm, such as Kronecker substitution [START_REF] Pan | Polynomial and matrix computations[END_REF][START_REF] Zur Gathen | Modern Computer Algebra[END_REF]. When truncating (P * Q * ) 2 e P +e Q +2-2q back to a complex floating polynomial R ∈ C ˜[z], we have

⌈⌈R -(P * Q * ) 2 e P +e Q +2-2 q ⌉⌉ 2 e R +2-p .
Consequently, we may take

f g = R + B S (2 e R +2-p + ↑ r × ↑ ⌈⌈Q⌉⌉ + ↑ s × ↑ ⌈⌈P ⌉⌉ + ↑ r × ↑ s).

The algorithm

Given an analytic function f ∈ A S with no zeros on S, we will denote

κ(f ) = 1 2 p i S f ′ (z) f (z) dz. If f ∈ A B , then κ(f ) counts the number of zeros of f in the open unit disk. Now consider a ball f = P + B S (r) ∈ B(C ˜[z],
R ˜) of analytic functions. Whenever no f ∈ f admits a zero on S, then κ(f ) 6 κ(f ) does not depend on the choice of f ∈ f . The aim of this paper is to compute κ(f ) in this case. If some f ∈ f admit zeros which are too close to S, then the algorithm is allowed to fail.

Rouché's theorem

The method that we will use for certifying the number of roots relies on the following variant of Rouché's theorem.

Theorem 1. Let f (z) and g(z) be two analytic functions on S, such that

|f (z) -g(z)| < |g(z)|.
Then κ(f ) = κ(g).

Proof. We will use a similar argument as in [START_REF] Lang | Complex analysis[END_REF]page 158]. Let γ be the path on S which turns one time around the origin and let F = f /g. By our assumption, the path F • γ is contained in the open disk with center one and radius one. Since this disk does not contain the origin, we have [Lan76, Lemma 3, page 116]

W (F • γ , 0) 6 1 2 p i F •γ dz z = 0. But W (F • γ , 0) = 1 2 p i γ F ′ (z) F (z) dz = κ(f ) -κ(g), whence κ(f ) = κ(g).
Given a polynomial P = P 0 + + P n z n ∈ C ˜[z] and k ∈ N, let us denote by P k ¯the polynomial P 0 + + P k-1 z k-1 + P k+1 z k+1 + + P n z n .

Proposition 2. Consider a ball f = P + B S (r) ∈ B(C ˜[z], R ˜) and let k be the index for which ⌈P k ⌉ is maximal. If

⌈⌈P k ¯+ B(r)⌉⌉ < ⌊P k ⌋,

Root counting

Putting together the various pieces, we now have the following algorithm for root counting. Optionally, one may add an extra integer parameter in order to limit the number of recursive calls to, say, ⌈2 log 2 p⌉. Let k be such that

Algorithm root-count(f ) Input: f ∈ B(C ˜[z], R ˜) Output: κ(f )
|P k | is maximal If r |P k |, then abort If ⌈⌈P k ¯+ B S (r)⌉⌉ < ⌊P k ⌋, then return k + v Return root-count(Gr(P + B S (r)))
One nice property of this algorithm is that the degree of P often quickly decreases during successive recursive calls. More precisely, let z 1 , , z n be the zeros of P , ordered such that

|z 1 | < |z k | < 1 < |z k+1 | < < |z n |. Then Gr i (p) Gr i (p) k = (z -z k 2 i ) (z -z 1 2 i ) (1 -z k+1 -2 i z) (1 -z n -2 i z).
In the right hand size, all coefficients which smaller than 2 -p are discarded via simplification. Roughly speaking, for j < k, the j-th coefficient therefore only survives as long as (z j z k ) 2 i > 2 -p .

Alternative termination strategies

In our algorithm, satisfaction of the condition ⌈⌈P k ¯+ B S (r)⌉⌉ < ⌊P k ⌋ enabled us to produce a final certified root count. It may be possible to design quick certified root counts for other favourable cases. Any alternative root counting strategy may actually be tried before applying our algorithm, or even at every stage as a replacement of the default strategy based on the condition ⌈⌈P k ¯+ B S (r)⌉⌉ < ⌊P k ⌋.

One natural alternative root counting strategy is based on the evaluation of f at many points on S. Let N to be the smallest power of two which is larger than p n and take ω = e 2pi/N . Then we may efficiently evaluate P at 1, ω, , ω N -1 using the FFT. Furthermore, the distance between two successive powers of ω is bounded by 1/n. Now assume that

⌈⌈P ′ ⌉⌉ 2 n < ⌊P (ω k )⌋ -r (2)
for all k. Then it is guaranteed that f admits no zeros on S. Hence, κ(f ) = 1 2 p k=0 n-1 arg P (ω k+1 ) P (ω k ) .

When multiplying the number of points N by a constant K, we may replace 2 n by 2 K n in the denominator of (2). We may also consider the second order condition

⌈⌈P ′′ ⌉⌉ 4 n 2 + ⌈P ′ (ω k )⌉ 2 n < ⌊P (ω k )⌋ -r (3)

  or failure Let f 6 σ(f ) and write f = P + B S (r) If P = 0, then abort Let v 6 val(P ) and P 6 P z -v

* . This work has been supported by the ANR-09-JCJC-0098-01 MaGiX project, as well as a Digiteo 2009-36HD grant and Région Ile-de-France.

then κ(f ) = k.

Proof. Let f = P + ε with ε r. Setting g = f k z k , our assumption implies

for all z ∈ S. By theorem 1, we conclude that κ(f

Graeffe transforms

The second ingredient for our algorithm is the Graeffe transform, which we will use for analytic functions. Given f ∈ A S , we define its Graeffe transform g = Gr(f ) ∈ A S by

If f is actually a polynomial, then we notice that g has the same degree as f . The fundamental property of Graeffe transforms is:

Proposition 3. On a small annulus containing S, the roots of g are precisely the squares of the roots of f, when counting with multiplicities.

Proof. By continuity under small deformations, it suffices to consider the case when g has only simple zeros z near S. Now

for all z near S.

Let us show that the Graeffe transform can be lifted to an operation Gr:

for all f = P + B S (r) ∈ B(C ˜[z], R ˜) and f ∈ A S . If r = 0, then we may directly use the formula Gr(f ) = P odd 2 z -P even 2 , assuming that the squares are computed using the algorithm for ball multiplication in section 2.5. In general, we take Gr(f ) = Gr(P + B S (0)) + B S (2 (⌈⌈P odd ⌉⌉ + ⌈⌈P even ⌉⌉) r + 2 r 2 ) and claim that this definition satisfies the inclusion property. Indeed, given f ∈ f , there exists an ε ∈ B S (r) with f = P + ε. Hence, Gr(f ) = Gr(P ) + 2 P odd ε odd z + ε odd 2 z -2 P even ε evenε even 2 .

Since ε even ε and ε odd ε , we thus have Gr(f ) -Gr(P ) = 2 P odd ε odd z + ε odd 2 z -2 P even ε evenε even 2 2 P odd r + r 2 + 2 P even r + r 2 2 (⌈⌈P odd ⌉⌉ + ⌈⌈P even ⌉⌉) r + 2 r 2 .

This proves our claim.


Kaikai Shen

Pierrick Bourgeat

Jurgen Fripp

Fabrice Meriaudeau

Olivier Salvado

Detecting Hippocampal Shape Changes in Alzheimer's Disease using Statistical Shape Models

The hippocampus is affected at an early stage in the development of Alzheimer's disease (AD). Using brain Magnetic Resonance (MR) images, we can investigate the effect of AD on the morphology of the hippocampus. Statistical shape models (SSM) are usually used to describe and model the hippocampal shape variations among the population. We use the shape variation from SSM as features to classify AD from normal control cases (NC). Conventional SSM uses principal component analysis (PCA) to compute the modes of variations among the population. Although these modes are representative of variations within the training data, they are not necessarily discriminant on labelled data. In this study, a Hotelling's T 2 test is used to qualify the landmarks which can be used for PCA. The resulting variation modes are used as predictors of AD from NC. The discrimination ability of these predictors is evaluated in terms of their classification performances using support vector machines (SVM). Using only landmarks statistically discriminant between AD and NC in SSM showed a better separation between AD and NC.

Description of purpose

Early detection and diagnosis of Alzheimer's disease (AD) is a challenging task. Since the hippocampus is affected by atrophy in the earliest stage of disease, which may result in the reduction of volume and the change in the shape of hippocampus. Hippocampal volume has been previously used to classify AD from normal control (NC) subjects, as well as cases with mild cognitive impairment (MCI) [START_REF] Chupin | Anatomically constrained region deformation for the automated segmentation of the hippocampus and the amygdala[END_REF]. Shape information in the form of spherical harmonics (SH) has been used as features in the support vector machine (SVM) classification [START_REF] Gutman | Disease classification with hippocampal shape invariants[END_REF][START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF]. Statistical Shape Models (SSMs) have been used to model the variability in the hippocampal shapes among the population (e.g. [START_REF] Davies | Shape Discrimination in the Hippocampus Using an MDL Model[END_REF]). They usually rely on principal component analysis (PCA) to determine a lower dimensional subspace that accounts for the most variations. However, these modes of variations are not necessarily discriminant. In this study, we aim to improve the discrimination between AD and NC using shape information characterized by SSM. We use SSM to search for variations on the surface regions that are significantly discriminant between the two groups. The discriminant regions may be due to both volume and shape changes. We propose to use the morphological variation on these surface regions as variables to distinguish AD from NC cases. These variations are also correlated with the cognitive decline in AD.

Methods

Materials

The hippocampus segmentations used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (http://www.loni.ucla.edu/ADNI), which provides the hippocampal volumes segmented semi-automatically by SNT.

SSMs on discriminant points

An SSM is built upon a set of shape examples, in which each shape is represented by its n P landmark points. In this study, the hippocampal volumes in the training set were first registered and aligned via rigid transformations, followed by a groupwise optimization and fluid regularization on the shape image [START_REF] Pantazis | Statistical surface-based morphometry using a nonparametric approach[END_REF] to establish the correspondence across the training set. Once the correspondence across the training set is established, the shapes can be aligned by using Procrustes analysis, either through rigid or similarity transformations. The volume information is preserved in rigid transformations. In this case, the variation in the training data will be partly driven by the volumetric change of hippocampus due to tissue loss. If the shape surfaces are aligned via similarity transformations, the shapes are rescaled to normalize the hippocampal volume. This enables the SSM to be more specific to the change in the shape per se rather than the variation in the size of hippocampus. Given a training set {x i :i=1,2, … , n S } of n S shapes, in which the correspondence is established and the shapes are aligned rigidly or via similarity transformations, each of its shape containing n P landmarks can be represented as a 3n P -vector where p i (k) =(x i (k) , y i (k) , z i (k) ) is the position of the kth landmark point on the ith shape surface. For the kth landmark {p i (k) :i=1, 2, … , n S } in the SSM, a Hotelling's T 2 test can be performed [START_REF] Pantazis | Statistical surface-based morphometry using a nonparametric approach[END_REF] in order to assess its statistical significance in discriminating the NC group from the AD group. Selecting only the m(<n P ) landmarks that are significantly discriminant on the labelled data, i.e. those landmarks with p-value below a threshold, the shape is represented by the subset of landmarks which is a projection mapping to the regions more relevant to the pathology. A SSM can thus be built concerning only the regions identified by the statistical test. For shape vector x, the coefficients of variation modes on discriminant regions can be calculated as where $\widetilde{\mathbf{W}}$ is the matrix of eigenvectors describing the variation modes from the PCA performed on only significantly discriminant landmarks. Thus, we have four choices of building SSM on the given training set with established correspondence:

• MRA using all the hippocampal landmarks, with shapes aligned via rigid transformation,

• MSA using all the hippocampal landmarks, with shapes aligned via similarity transformation,

• MRD using hippocampal landmarks significant against AD, with shapes aligned via rigid transformation, • MSD using hippocampal landmarks significant against AD, with shapes aligned via similarity transformation.

Representation of the test data

For a given shape not in the training set, the correspondence between the SSM landmarks and the shape surface points needs first to be established, so that the given shape can be represented in the same vector space as the training data. The SSM is deformed to fit the smoothed target surface, minimizing L1 distance metric between the SSM generated surface and the target. Without the assumption of correspondences between two shape vectors x and y, the L1 distance d S between their surfaces S(x) and S(y) can be defined as the sum of the Euclidean distance from each point in x to its closest point on y and from each point y to x. Thus we can fit the SSM to the target surface S(y) by the optimization of parameters using Powell's algorithm where T is a similarity transformation with 7 degrees of freedom, and b parameterizes the deformation of the SSM. For the each landmark generated by SSM, we can find the closest point on the target surface as the corresponding landmark point. By aligning the shape vector of correspondences to the shapes in the SSM through Procrustes analysis, we have the n P -D shape vector representing the surface S(y) in the same space as the shapes in the SSM, ready for further analysis.

Classification method

The shape descriptors from SSM are evaluated in terms of their performance when being used as features in classification algorithms. SVMs are widely used in solving the classification problem. It usually maps the feature space to higher dimension via a kernel function, and finds the optimal hyperplane with the largest margin separating the classes. We use the classifier trained by the SVM to test the discrimination ability of the features. In this study, the radial basis function was used as the kernel of SVM. The features of each shape x for classification are the coefficient variation modes. In order to avoid modelling noise, less significant components produced by PCA are not included in the feature set. A subset of features is selected using random forest by minimizing the out-of-bag (OOB) error [START_REF] Diaz-Uriarte | Gene selection and classification of microarray data using random forest[END_REF].

Results

In the experiment, hippocampi on both the right and left sides of 60 NC and 60 AD subjects from ADNI data were used as a training set to built the SSM. Hotelling's T 2 test was performed on each SSM between the AD and NC groups, with the resulting significance maps shown in Fig. 1. Landmarks with p < 0.01 68.4 79.5 84.6 85.5 A PCA was performed on the landmarks with p < 0.05 and p < 0.01. For the purpose of comparison, conventional PCA on all the landmarks was also performed. The first 15 principal components in each of the SSM were used for feature selection. These components explained approximately 90% of variations in the training set. In addition to the shape variations, the hippocampal volume normalized by TIV was also used as independent features. The resulting SSMs were evaluated on a testing set of 117 subjects (39 AD and 78 NC), also from ADNI. The results of the accuracy on the testing set are shown in Table 1. When using all the points to perform the PCA, the features extracted from the SSM built using a rigid transform reflect largely the volume differences between the hippocampi, whereas the ones extracted from the SSM built using a similarity transform are a better reflection of local shape variations. This is effect results in the MRA yielding better classification results than MSA (77.8% vs 72.6%) as volume alone is a good discriminant (83.8%). Adding volume to the features extracted using MRA and MSA increased both of their accuracy, but to a much less extent for MRA (+6.0%) compared to MSA (+12.8%). Restricting the PCA to discriminant points on shapes aligned via similarity transformation produced better discrimination between NC and AD than using all the points. Using the shape variation on discriminant points extracted from MSD yielded better accuracy than the global variation incorporating both shape and volume characterised by MRA. The best classification results were obtained when volume was added as an additional feature to the shape features obtained using the SSM built using a similarity transform on the discriminant points (MSD).

New and breakthrough work to be presented

• Model the shape variation between the AD and NC population of hippocampus on the regions affected by atrophy identified by statistical tests;

• Improving the disease classification performance by using the shape descriptors extracted from atrophy affected areas on the hippocampal surface.

Conclusions

The shape of the hippocampus can provide valuable information for the diagnosis of AD. The variation modes of the hippocampus among the population as modelled by the SSM can be used to classify AD against NC. The conventional PCA in SSM is performed on all the landmark points on the shape which is a good representation of the original shape data in a lower dimensional subspace, while it might be not discriminant between two groups. By applying a statistical test on the landmark points in the SSM, we can identify the regions on the hippocampal surface which are more discriminant between AD and NC groups. The PCA performed on this subset produced variation modes which were used as features for the classification between these two groups.

Fig. 1 :

 1 Fig. 1: Signicance map by Hotelling's T 2 test, performed on rigid and similarity aligned SSMs. Top: superior view; bottom: inferior view.

Table 1 .

 1 Results of SVM accuracy on a test set of 117 cases (%).

	With Volumety


Soumya Ghose

Arnau Oliver

Robert Martí

Xavier Lladó

Jordi Freixenet

Joan C Vilanova

Fabrice Meriaudeau

Prostate Segmentation with Local Binary Patterns Guided Active Appearance Models

Keywords: Prostate Segmentation, Active Appearance Model, Texture, Local Binary Patterns

Real-time fusion of Magnetic Resonance (MR) and Trans Rectal Ultra Sound (TRUS) images aid in the localization of malignant tissues in TRUS guided prostate biopsy. Registration performed on segmented contours of the prostate reduces computational complexity and improves the multimodal registration accuracy. However, accurate and computationally efficient segmentation of the prostate in TRUS images could be challenging in the presence of heterogeneous intensity distribution inside the prostate gland, and other imaging artifacts like speckle noise, shadow regions and low Signal to Noise Ratio (SNR). In this work, we propose to enhance the texture features of the prostate region using Local Binary Patterns (LBP) for the propagation of a shape and appearance based statistical model to segment the prostate in a multi-resolution framework. A parametric model of the propagating contour is derived from Principal Component Analysis (PCA) of the prior shape and texture information of the prostate from the training data. The estimated parameters are then modified with the prior knowledge of the optimization space to achieve an optimal segmentation. The proposed method achieves a mean Dice Similarity Coefficient (DSC) value of 0.94±0.01 and a mean segmentation time of 0.6±0.02 seconds when validated with 70 TRUS images of 7 datasets in a leave-one-patient-out validation framework. Our method performs computationally efficient and accurate prostate segmentation in the presence of intensity heterogeneities and imaging artifacts.

INTRODUCTION

More than 670,000 people are diagnosed with prostate cancer every year worldwide 1 accounting for over 32,000 deaths in Noth America. [START_REF]Prostate Cancer Statistics -Key Facts[END_REF] In clinical practice, TRUS guided needle biopsy is performed to confirm the presence of cancerous tissues in prostate. However, localization of malignant tissues in TRUS is difficult due to the low soft tissue contrast in TRUS images. Multi-core prostate biopsy is performed from 6 to 12 cores to improve positive detection of malignant tissue. [START_REF] Rodriguez | Multi-core prostatic biopsy[END_REF] MR images provide higher contrast for soft tissues of the prostate that allows a better localization of cancerous tissues. However, interventional MRI guided biopsy is expensive and complicated. The localization of malignant tissues could be improved by real-time registration of a 2D TRUS video sequence with a 3D pre acquired MRI during TRUS guided biopsy. [START_REF] Xu | Real-time MRI-TRUS Fusion for Guidance of Targeted Prostate Biopsies[END_REF] Registration performed on segmented contour of TRUS and MR images will aid in designing computationally efficient and accurate registration procedures. [START_REF] Yan | Optimal Search Guided by Partial Active Shape Model for Prostate Segmentation in TRUS Images[END_REF] Thus, 2D segmentation of the prostate from TRUS video sequence is necessary for the process. However, automatic or semi-automatic accurate segmentation of the prostate in TRUS image is a challenging task due to low SNR, intensity heterogeneities inside the prostate gland and presence of imaging artifacts like the speckle noise and shadow regions. Low SNR produces spurious and broken edges with traditional edge detection filters and intensity heterogeneities inside the prostate gland inhibit the design of a global descriptor for the prostate. Added to these shadow artifacts, speckle noise and micro-calcifications is a hindrance to accurate segmentation for the segmentation algorithms. The challenges involved in accurate prostate segmentation in TRUS images are illustrated in Fig. 1. To address the challenges of prostate segmentation in TRUS images Shen et al. [START_REF] Shen | Segmentation of Prostate Boundaries from Ultrasound Images Using Statistical Shape Model[END_REF] and Betrouni et al. [START_REF] Betrouni | Segmentation of Abdominal Ultrasound Images of the Prostate Using A priori Information and an Adapted Noise Filter[END_REF] proposed to use prior prostate shape information in their models. Prior shape information made the models robust to imaging artifacts and detected true prostate edges in low SNR. Yan et al. [START_REF] Yan | Optimal Search Guided by Partial Active Shape Model for Prostate Segmentation in TRUS Images[END_REF] proposed to incorporate prior shape information of the prostate in a partial active shape model framework. The model selectively discarded salient points for building the shape model to improve segmentation accuracy in presence of shadow artifacts in TRUS images. Gong et al. [START_REF] Gong | Parametric Shape Modeling Using Deformable Superellipses for Prostate Segmentation[END_REF] proposed to use deformable super ellipse to segment the prostate. Shape constraints of the model was effective in achieving impressive segmentation results.It has been demonstrated by Zhan et al. [START_REF] Zhan | Deformable Segmentation of 3D Ultrasound Prostate Images Using Statistical Texture Matching Method[END_REF] that incorporating shape and intensity prior information in prostate segmentation methods improve the segmentation accuracy. Cootes et al. [START_REF] Cootes | Active Appearance Models[END_REF] provided an efficient framework for combining shape and intensity prior in their Active Appearance Model (AAM). AAM framework was used by Medina et al. [START_REF] Medina | A 2D Active Appearance Model For Prostate Segmentation in Ultrasound Images[END_REF] to segment prostate in 2D TRUS images with an overlap ratio of 96%. However, time complexity involved with AAM is high making the model unsuitable for real-time procedures. [START_REF] Wolstenholme | Wavelet Compression of Active Appearance Models[END_REF] Larsen et al. [START_REF] Larsen | Texture Enhanced Appearance Models[END_REF] proposed to use truncated wavelet coefficients in AAM framework to improve on segmentation accuracy of facial images and achieve reduction in segmentation time.

In this paper, we propose a novel AAM that is propagated by texture features of Local Binary Patterns (LBP) [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] transformed images in a multi-resolution framework. The motivation behind such an approach is to provide precise and concise texture information for the propagation of the AAM to improve on computational time and segmentation accuracy. The performance of our method is validated with 70 TRUS images of 7 datasets in a leave-one-patient-out validation framework. Experimental results show that our method is unaffected by low SNR, intensity heterogeneities and shadow and speckle noise. The key contributions of this work are:

• Enhancement of texture information of the prostate region from local binary patterns. [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] • Use of the texture information in training and propagation of an active appearance model that improves on computational time and segmentation accuracy when compared to the traditional active appearance model. [START_REF] Cootes | Active Appearance Models[END_REF] To the best of our knowledge this is the first attempt to use enhanced texture features of the prostate from LBP in training and propagation of AAM. The rest of the paper is organized as follows. The texture driven AAM is formulated in Section 2 followed by quantitative and qualitative evaluation of our method in Section 3. We finally draw conclusion in Section 4. The working schema of the model is illustrated in Fig. 2 2

. METHOD

In our method the traditional AAM is adapted for the incorporation of texture information extracted from the LBP transformed image. Traditional AAM is presented here, followed by a comprehensive discussion about enhancement of texture information using LBP.

Active Appearance Model

Generalized Procrustes Analysis of the Point Distribution Model (PDM) [START_REF] Cootes | The Use of Active Shape Model for Locating Structures in Medical Images[END_REF] built from manually segmented contours aligns the PDM in a common reference frame. Principal Component Analysis (PCA) of the aligned PDMs is used to identify the principal modes of shape variations. Similarly, PCA of the intensity distribution warped into correspondence using a piece wise affine warp and sampled from shape free reference is used to identify the principal modes of intensity variations. Principle modes of shape and intensity variations are combined and a third PCA ensures that primary components of combined model variations are identified. On completion of building the combined model, an unseen test image is segmented exploiting the prior knowledge of the nature of the optimization space. [START_REF] Larsen | Texture Enhanced Appearance Models[END_REF] The model could be formalized in the following manner. In eq. 2 let E {s} and E {t} represent the shape and intensity models where s and t are the shape and the intensities of the corresponding training images, s and t denote the mean shape and intensity respectively, then Φ s and Φ t contain the first p eigenvectors of the estimated joint dispersion matrix of shape and intensity and θ represents the corresponding eigenvalues.

E {s} = s + Φ s θ (1)
E {t} = t + Φ t θ
In addition to the parameters θ, four parameters, two translations, rotation and scale are represented by ψ. In order to infer the parameters of θ and ψ of a previously unseen image, a Gaussian error model between model and pixel intensities is assumed. [START_REF] Larsen | Texture Enhanced Appearance Models[END_REF] Furthermore, a linear relationship between changes in parameters and difference between model and image pixel intensities ∆t is assumed as shown in eq.2.

∆t = X ∆ψ ∆θ (2)
X is estimated from weighted averaging over perturbation of model parameters and training examples. Eq.2 is solved in a least square manner fitting error as shown in eq.3.

ψ θ = (X T X) -1 X T δt (3)
The problem is computationally expensive. To reduce computational burden we propose to use a truncated basis for the representation of the pixel intensities achieved with LBP. However, this will introduce the additional time requirement of transformation of the image into a new representation. Since the transformation is based on sparse matrix, the computational burden can be considerably reduced as stated by Larsen et al. [START_REF] Larsen | Texture Enhanced Appearance Models[END_REF]

Local Binary Patterns

Ojala et al. [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] proposed LBPs as computationally efficient discriminative transformations for texture analysis. Given a pixel in the image, LBP is computed by thresholding the neighborhood of the pixel. LBP of a pixel is defined as in eq. 4

LBP P,R = P -1 p=0 s (g p -g c ) 2 P , s (x) = 1, x ≥ 0 0, x < 0 (4)
where, g c is the gray value of the centered pixel, g p are the gray values of the neighborhood pixels, P is the number of neighbors and R is the radius of the neighborhood. Signed difference g pg c is unaffected by change in mean luminance that provides gray scale invariance to the operator. Rotational invariance (RI) is defined in eq. 5.

LBP ri P,R = min {ROR (LBP P,R,i ) |i = 0, 1....., P -1}|

where ROR (x, i} performs a circular bit-wise right shift on the P -bit number x, i times to assign a unique identifier to each rotation invariant local binary pattern. Ojala et. al. observed some fundamental patterns depicting the texture properties of an image under study. They named these patterns as Uniform circular patterns (U) of LBP transformed image. Uniform circular pattern is defined as in eq. 6 and 7

LBP riu2 P,R = P -1 p=0 s (g p -g c ) if U (LBP P,R ) ≤ 2 P + 1 otherwise (6)
where

U (LBP P,R ) = |s (g p-1 -g c ) -s (g 0 -g c )| + P -1 p=1 |s (g p -g c ) -s (g p-1 -g c )| (7)
Superscript riu2 refers to rotational invariant "unif orm" patterns. Where "unif orm" patterns are given by eq. 7. Examples of LBP transformed images are given in Fig. 3 We have used LBP, uniform LBP and rotational invariant uniform LBP to produce LBP transformed representation of the image and used the PCA of the texture of the prostate region to propagate our AAM in a multi-resolution framework. The LBP transformed representation of the image is a truncated representation of the original intensity and hence fitting of the new image after LBP transformation is computationally efficient. Multi-resolution fitting improves segmentation accuracy and reduces the computational time of the process. The multi-resolution fitting of the model is illustrated in Fig.

RESULTS

We have validated the accuracy and robustness of our approach with 70 TRUS images of prostate in leave-onepatient-out evaluation strategy. The images are of resolution 538 × 418 obtained from 7 prostate datasets. Our method was implemented in Matlab 7 on an Intel Quad Core Q9550 processor of 2.83 Ghz processor speed and 8 GB RAM. We have used most of the popular prostate segmentation evaluation metrics like Dice Similarity Coefficient (DSC), 95% Hausdorff Distance (HD), Mean Absolute Distance (MAD), Maximum Distance (MaxD), specificity, sensitivity, accuracy and mean segmentation time to evaluate our method. We have used two different values for the number of sample points (P = 16, 8) and three different radii (R = 1, 2, 8) for the LBP and presented the best results in Table 1. LBP 16,2 corresponds to the traditional LBP as given in eq. 4, U (LBP 16,2 ) represents the uniform LBP as given in eq. 7 , and LBP riu 16,8 represents rotational invariant uniform LBP as given in eq. 6. Subscript values denote neighbors and radius respectively. The results are given in Table 1 In Table 1 we observe that; in general LBP and different variants of LBP guided AAMs perform better than the traditional AAM comparing the distance metrics HD and MAD. This could be attributed to the fact that gray scale invariance of the LBP operator helps gray-scale normalization of AAM, suppressing intensity heterogeneities, speckle noise and micro-calcifications in the process. LBP enhance underlying texture information of the prostate region. PCA of the LBP enhanced texture information reduces noise and helps in identifying principle modes of texture variation inside the prostate region. Surprisingly, computational complexities involved with traditional LBP method is less compared to traditional AAM considering the fact that the image had to be transformed into LBP map before fitting. On closer inspection it was revealed that average time for transforming 70 images into LBP map was 0.025 seconds and it produced an accurate representation of the prostate texture thereby reducing differences with the mean model. This in turn reduces the time for computationally expensive procedure of fitting a new image. The performance of Uniform LBP is better compared to the other variants of LBP. This is due to the fact that uniform LBP produces a uniform profile of the prostate region that is different from the neighboring tissue textures and therefore facilitates in improving the accuracy of segmentation. However, the procedure is computationally expensive compared to other variants of LBP. It is also observed that the performance of rotational invariant uniform LBP LBP riu 16,8 is almost comparable to traditional AAM. The robustness of the proposed method against low SNR, intensity heterogeneities, shadow artifacts, speckle noise and micro-calcification is illustrated in Fig. 5. On initialization, a section of the mean model is located in a region of shadow artifact. The model successfully avoids the artifact and segments the prostate (black contour) with an accuracy of 97%.

The mean model is initialized by clicking on the center of the prostate judged from visual inspection. However, the model could be initialized automatically with the assumption that the prostate is located at the center of the TRUS image. A more sophisticated approach would be to use an initial rough classification to initialize the mean model. Low standard deviation value associated with sensitivity, specificity and accuracy indicate that the final segmentation result is not affected by the manual initialization. To further validate our claim ten images were selected randomly and the mean model was initialized for 10 times in ±5 pixels around the center for uniform LBP U (LBP 16,2 ). Mean accuracy of these initializations is 96±0.006%. Low standard deviation associated with accuracy metric indicates that the process is indifferent to manual initialization as long as the initialization is close to the prostate center. To provide qualitative result of our method we present a subset of results obtained using uniform LBP operator U (LBP 16,2 ) in Fig. 6 Comparison of different prostate segmentation methodologies is difficult in absence of public datasets and standardized evaluation metrics. In addition, the methods are developed with a wide variety of algorithms with specific application requirements. Nevertheless, to have an overall qualitative estimate of the functioning of our method we have compared our uniform LBP AAM with some of the state-of-the-art works in the literature (see Table 2) for TRUS images. Please note '-' in Table 2 means information is not available.

The average DSC value of our method with 0.95±0.01 is comparable to mean DSC values achieved by Betrouni et al. [START_REF] Betrouni | Segmentation of Abdominal Ultrasound Images of the Prostate Using A priori Information and an Adapted Noise Filter[END_REF] and Ladak et al. [START_REF] Ladak | Prostate Segmentation from 2D Ultrasound Images[END_REF] and inferior to the area overlap measure of Abolmaesumi et al. [START_REF] Abolmaesumi | Segmentation of Prostate Contours from Ultrasound Images[END_REF] of 98%. However,

CONCLUSION AND FUTURE WORKS

A novel approach of AAM propagation using LBP enhanced texture profile with the goal of segmenting the prostate in 2D TRUS images has been proposed. Our approach is accurate, computationally efficient and robust to low SNR, intensity heterogeneities of the prostate tissues, shadow artifacts, speckle noise and microcalcification in TRUS images. While the proposed method is validated with prostate mid gland images, effectiveness of the method against base and apical slices is yet to be validated. Computational time of the process is fast but not suitable for real time segmentation of the prostate. We would like to explore the possibilities of optimization of our code and use CUDA platform to achieve necessary parallelization and hardware acceleration to achieve real-time segmentation of the prostate. We would like to explore the possibilities of automatic model initialization and extend the model for 3D segmentation of the prostate.

Figure 1 .

 1 Figure 1. Prostate segmentation challenges , A=Low SNR, B=Micro Calcification, C=Intensity heterogeneity inside prostate, D=Shadow Artifacts, E=Speckle Noise.

Figure 2 .

 2 Figure 2. Schematic representation of our approach. The final segmentation is given in red contour and ground truth in green. Abbreviations used AAM = Active Appearance Model, PDM = Point Distribution Model, GPA = Generalized Procrustes Analysis, LBP = Local Binary Patterns, RIU = Rotational Invariant Uniform.

4

 4

Figure 3 .Figure 4 .

 34 Figure 3. Variants of LBP transformed image(a) 2D TRUS image of the prostate, (b) LBP transformed image, (c) Uniform LBP transformed image, radius=2, samples=16 (d) Rotational invariant uniform LBP transformed image.Radius=2 and Samples=16.

Figure 5 .

 5 Figure 5. (a)Artifacts in TRUS image of the prostate, A=Low SNR, B=Micro Calcification, C=Intensity heterogeneity inside prostate, D=Shadow Artifacts, E=Speckle Noise. (b) Manual initialization of the mean model, (c) Final segmentation result.

Figure 6 .

 6 Figure 6. Qualitative segmentation results. The white contour gives the ground truth and the black contour gives the obtained result.

Table 2 .

 2 Prostate Segmentation Evaluation: A Qualitative Comparison with the state-of-the-art algorithms for TRUS it is to be noted that we have used more images (24 images) compared to 6 images of Abolmaesumi et al. Our MAD value of 1.50±0.41 mm is comparable to Gong et al.[START_REF] Gong | Parametric Shape Modeling Using Deformable Superellipses for Prostate Segmentation[END_REF] average distance value of 1.36±0.6 mm and to Yan et al.[START_REF] Yan | Optimal Search Guided by Partial Active Shape Model for Prostate Segmentation in TRUS Images[END_REF] MAD of 1.79±0.95 mm. From these observations we can safely conclude that qualitatively our method performs well in overlap and contour accuracy measures over some of the works in literature.

	Images
	Reference	Year Area Overlap Accuracy	Contour Accuracy	Datasets
	Yan 5	2009 -	Average MAD 1.79±0.95 mm	10 datasets
	Abolmaesumi 16 2004 Average Area overlap 98%	-	6 images
	Gong 8	2004 -	Mean error 1.36±0.6 mm	125 images
	Betrouni 7	2004 Average Area overlap 93±0.9%	Average distance 3.77±1.3 pixels 10 images
	Shen 6	2003 Average Area overlap error 3.98±0.97% Average distance 3.2±0.87 pixels 8 images
	Ladak 17	2000 Average Area accuracy 90.1±3.2%	Average MAD 4.4±1.8 pixels	117 images
	Our	Average DSC 0.95±0.01	Average MAD 1.50±0.41 mm	7 datasets/70
	Method			images

ACKNOWLEDGMENTS

This research has been funded by VALTEC 08-1-0039 of Generalitat de Catalunya, Spain and Conseil Régional de Bourgogne, France.