Camille Couprie
email: c.couprie@esiee.fr

Xavier Bresson
email: xbresson@cityu.edu.hk

Laurent Najman
email: l.najman@esiee.fr

Hugues Talbot
email: h.talbot@esiee.fr

Leo Grady
email: leo.grady@siemens.com

Surface reconstruction using Power Watershed

Keywords: optimization, point measurements, Graph cuts, total variation

à 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

This paper develops a watershed-based algorithm providing a global optimal solution to the surface reconstruction problem from a set of scattered points. Surface fitting is a challenging problem when dealing with data containing sparse noise, gaps, and outliers. The set of points may be for example acquired by several scans of an object (range scanning). In this context, regularization-based methods have been shown to be robust when the points are lacking connectivity, ordering information, and may be contaminated by noise. While there exist numerous explicit surface extraction techniques that estimate the exact positions of surface points, in this work we will focus on implicit surface representation. Implicit surfaces may be represented by level sets (e.g. [START_REF] Zhao | Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method[END_REF]) or binary partitions (e.g. [START_REF] Lempitsky | Global Optimization for Shape Fitting[END_REF]).

Local methods for surface reconstruction including the MPU method [START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF] are sensitive to noise, as shown in the experiments of [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF]. Among the recent global approaches, the Poisson method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] and FFT-method [START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF] are more robust to noise, however they require orientation information.

The method of Jalba and Roerdink [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF] makes it possible to avoid having to estimate orientation information. This is achieved by computing approximations of Coulomb potentials in a grid as an input for the convection method of [START_REF] Zhao | Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method[END_REF]. We propose a different approach that allows us to perform a global optimization of the surface reconstruction problem without requiring the computation of such field as Coulomb potentials.

A generic optimization-based regularization formulation for shape fitting minimizes the total variation functional weighted by the distance function from the set of points P . More generally, given two positive numbers p and q, we consider an object indicator partition u solution of min u∈[0,1] Ω w(z) p |∇u(z)| q dz subject to u(z) = 0 ∀z ∈ Ω in , and u(z) = 1 ∀z ∈ Ω out , [START_REF]Stanford 3D scanning repository[END_REF] where Ω in is the set of labels inside the surface and Ω out is the set of labels outside the surface. The weight function w is defined at every point z in a grid as w(z) = d P (z), where d P (z) is the distance map from the points. When p is finite and q = 1, Eq. ( 1) leads to a binary solution u [START_REF] Strang | Maximum flows through a domain[END_REF]. A solution can be deduced in the discrete setting using e.g. the network flow technique [START_REF] Ford | Maximal flow through a network[END_REF], also known as Graph cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. Augmenting path max flow implementations are fast and efficient in 2D but memory consuming as the connectivity increases, for instance in 3D. Lempitzky and Boykov [START_REF] Lempitsky | Global Optimization for Shape Fitting[END_REF] have overcome this problem by limiting the size of the search for a solution in the grid while still guaranteeing a global optimum. However their solution is based on restrictive assumptions assuming that the data points are provided with an estimate of the surface orientation. Furthermore, at low resolution, results exhibit metrication artifacts and look blocky, so a high resolution is essential for getting smooth results using the Graph cuts method.

We propose a watershed-based approach in order to provide a way to quickly obtain smooth surfaces at a high resolution without any need to pre-estimate the surface orientation. Recently, Couprie et al. introduced the power watershed method [START_REF] Couprie | Power Watersheds: A Unifying Graph Based Optimization Framework[END_REF][START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF], which can be seen as an anisotropic discretization of (1) with p → ∞ and q = 2. Although this technique was introduced in the context of image segmentation, the authors described how the method could be used as an optimization method for various functionals in [START_REF] Couprie | Anisotropic Diffusion Using Power Watersheds[END_REF]. In the present paper we show that the power watershed method of Couprie et al. is well-suited to address the surface reconstruction problem. The idea of using watersheds for surface reconstruction from a set of points is quite natural. It can be seen as a an extension of the classical "coffee bean" segmentation example, where a watershed is applied on a filtered distance function to separate overlapping convex objects [START_REF] Beucher | Micromorph reference manual, applications and solutions[END_REF].

For further comparison, we also propose in this paper to examine a wellknown weighted isotropic discretization of (1) (for q = 1) known as "Total Variation" (TV). We mention that there exist other TV-based approaches developed to overcome metrication artefacts, e.g. [START_REF] Ye | A fast variational method for surface reconstruction from sets of scattered points[END_REF][START_REF] Goldstein | Geometric applications of the split Bregman method: Segmentation and surface reconstruction[END_REF]. However, due to space constraints, further comparison with these techniques will be the subject of future work. Finally, we will also show that our algorithm compares favorably with existing surface reconstruction algorithms [START_REF] Amenta | The power crust[END_REF][START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF][START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF][START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF].

Method

For solving the surface reconstruction problem, we first place the points cloud onto a regular grid. Our method aims to label the nodes of the grid as an indicator of the object to reconstruct. Since the power watershed is defined on a graph, we begin by casting the surface reconstruction problem formulation in discrete terms.

A graph consists of a pair G = (V, E) with vertices v ∈ V and edges e ∈ E ⊆ V × V with cardinalities n = |V | and m = |E|. An edge, e, spanning two vertices, v i and v j , is denoted by e ij . We define an edge set corresponding to a 4 or 8-connected lattice (or 6-, 18-or 26-connected in 3D). A weighted graph assigns a (typically non-negative and real) value to each edge called a weight. The weight of an edge e ij is denoted by w(e ij ) or w ij . A plateau is a maximal set of connected edges with identical weight.

Given foreground F and background B node values (also called seeds), and p, q two real positive values, the energy presented for binary segmentation in [START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF] is a discretization of (1) given by min

x eij ∈E w p ij |x i -x j | q s.t. x(F ) = 1, x(B) = 0. (2)
In this energy, x is a labeling indicating the foreground and background membership. The edge values w ij can be interpreted as weights enforcing a regularization of the contours, such that any (usually unwanted) high-frequency content is penalized in x. The definition of the weights for surface reconstruction from a set of points P is based on the construction of a discrete distance map d P , in a grid bounding the set of points.

w ij = min (d P (i), d P (j)) , (3)
where d P (i) is the discrete Euclidean distance between the node i and the set of points P . We recall that exact Euclidean discrete distance map may be obtained in linear time using the algorithm of Hirata [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF], and that high-quality ordered algorithms also exist [START_REF] Ragnemalm | The euclidean distance transform in arbitrary dimensions[END_REF] . The background seeds may simply correspond to the frame, or bounding box of the lattice. The foreground seeds can be given by the maxima of the distance function that are not connected to the frame. The distance map may be previously filtered, for example using an attribute filter to obtain more robust markers [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Najman | Building the component tree in quasi-linear time[END_REF]. As we illustrate in the remainder, the energy defined in (2) essentially forces x to remain smooth within the object, while allowing it to vary quickly close to point clusters near the boundary of the object. The data constraints enforce fidelity of x to a specified configuration, taking the values zero and one as the reconstructed object indicator. Observe that the values of x may not necessarily be binary when the value of q is strictly greater than one, which is a positive point for the surface reconstruction problem as we will further explain in this work.

The different values of p and q lead to different algorithms for optimizing the energy. When the power of the weight, p, is finite, and the exponent q = 1, we recover the Graph cuts energy which can be optimized by a max flow algorithm. When q = 2, we obtain a combinatorial Dirichlet problem also known as the Random walker problem [START_REF] Grady | Random walks for image segmentation[END_REF]. As described in [START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF][START_REF] Couprie | Power Watersheds: A Unifying Graph Based Optimization Framework[END_REF], when the exponent p tends toward infinity, the cut obtained when minimizing the energy is a watershed cut [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF], which has been shown to be equivalent to Maximum Spanning Forests [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] (MSF). Furthermore, an algorithm is presented to optimize the unique watershed that optimizes the energy for q = 2 and p → ∞. The power watershed (PW) algorithm is recalled in Algorithm 1.

Algorithm 1: Power watersheds algorithm p → ∞, q = 2 Data: A weighted graph G(V, E) comprising known labels x(B), x(F ).

Result: A labeling x solution of (2).

while any node has an unknown label do Find a maximal subgraph S ∈ G composed of edges of maximal weight; if S contains any nodes with known x then Find x S minimizing (2) for q = 2 on the subset S; Consider all x S values produced by this operation as known; else Merge all of the nodes in S into a single node, such that when the value of x for this merged node becomes known, all merged nodes are assigned the same value of x and considered known;

This set of parameters q = 2 and p → ∞ is particularly interesting :

1. The power watershed algorithm has a well-defined behavior in the absence or lack of weight information (presence of plateaus). An example is shown at Figure [START_REF]Stanford 3D scanning repository[END_REF]. 2. The worst-case complexity of the power watershed algorithm in the case p → ∞ is given by the cost of optimizing (2) for the given q. In the best-case scenario (all weights have unique values), the power watershed algorithm has the same asymptotic complexity as the algorithm used for a MSF computation (quasi-linear) (See [START_REF] Chazelle | A minimum spanning tree algorithm with inverse-ackermann type complexity[END_REF] for more details). In practical applications where the plateaus have a size less than some fixed value K, then the complexity of the power watershed algorithm matches the quasi-linear complexity of the standard watershed algorithm. Although the PW algorithm is fast, it can be further accelerated for the specific case of weight defined according to a distance map [START_REF] Beucher | Micromorph reference manual, applications and solutions[END_REF]. Given the foreground and background seeds, we define a narrow band S as the set obtained by thresholding the distance map d P with the the smallest threshold T S such that the connected components are divided between at least an interior (foreground) and an exterior (background). We then compute the PW only on this incomplete distance map, which saves both time and memory. In practice, it is possible to avoid computing an exact distance map on the full grid, using for instance an ordered algorithm propagating from the point cloud [START_REF] Ragnemalm | The euclidean distance transform in arbitrary dimensions[END_REF].

Recall that the exterior seed is connected to the frame of the image, so this is a simple unambiguous connectivity criterion. Applying PW on S is guaranteed to provide the same (unique) global optimizer of the energy as the solution on the full grid, because the connectivity criterion ensures this computation yields a Jordan curve (2D) or surface (3D) separating foreground and background seeds and passing through only already-computed distance values. Performing the PW computation on the full map would not change this result because no distance weight would be lowered, and so no new surface with a lower weight could be found on this full map. Conversely, thresholding the distance map at a lower value than T S would change the result, as any Jordan curve or surface separating foreground and background seeds and computed on this map would necessarily cross some nodes where the distance map had not been computed. Note this does not mean that this thresholding criterion is necessary and sufficient for optimal computation and least memory usage, as an adaptive threshold depending on the local point density could be used instead. However, we leave it as future work to find a better criterion than the simple global threshold.

3 Results and comparative evaluation

Comparison with Graph cuts and Total Variation

We now demonstrate the performance of the power watershed algorithm for surface reconstruction with respect to two graph-based methods discretizing the energy defined in equation ( 1), namely the weighted total variation (TV) and the Graph cuts (GC) method.

Our first experiment consists of finding a contour fitting sparse and noisy dots in a 2D plane. Figure 2 compares the result of TV minimization, the Graph cuts result, a maximum spanning forest result, and the result obtained with power watershed algorithm (q = 2, p → ∞). We observe that all resulting contours are excluding outliers. The Graph cuts results demonstrate that this algorithm is less sensitive to noise, but the contours are blocky because the obtained object indicators are binary. We note that a post-processing step for producing a smooth isosurface from such a binary object reconstruction has recently been proposed by Lempitzky [START_REF] Lempitsky | Surface extraction from binary volumes with higher-order smoothness[END_REF]. The TV contour is the smoothest one, which may be an unwanted effect for rendering details in surfaces. Thus, the TV method requires us to adapt the parameter p in the exponent of the weights to the desired smoothness of the surface. The maximum spanning forest result passes through most points and the contour looks noisy. The power watershed result demonstrates good performance for fitting the dots, while avoiding both blocky contours and noise. This good performance is due to the presence of interwoven plateaus in the distance map. During the execution of the algorithm, the Random walker (RW) algorithm is called several times around the dots, resulting in a smooth output x. An isocontour or isosurface computation at the 0.5 level is thus providing smooth contours compared to the binary results obtained with Graph cuts or maximum spanning forests.

Figure 3 shows surfaces reconstructed from noisy scanned dot sets using the power watershed (PW) algorithm. We used the coordinates of points acquired from scans of several 3D shapes (bunny, Buddha) from the Stanford database available online [START_REF]Stanford 3D scanning repository[END_REF] In our experiments, we embedded those points in 6-connected grids. Quantitative comparisons for the fitting quality are difficult because not all data points are required to be part of the surface. However, we show in Figure 4 that the power watershed method is producing reasonable results even if the point cloud is corrupted by Gaussian noise. We also compared our results to the results obtained using TV minimization and Graph cuts at Figure 6. We can observe that the surface obtained with the Graph cuts method is quite blocky. Using the same rendering method to render the output x minimizing the power watershed energy, the power watershed algorithm obtains a smoother surface showing significantly more details. Figure 5 shows that the power watershed method performs well even with a small amount of seeds. In contrast, the total variation method requires a large amount of seeds placed close to the searched surface.

Computation times and memory requirements

In our comparisons, we used the C++ software library of Lemptitzky and Boykov available online which implements the touch-expand algorithm that minimizes the GC energy. For that purpose, the touch-expand algorithm calls a max flow algorithm on a partial graph, a band around the points which is extended when the solution touches the boundary of the band. We also compare to a fast TVbased solver often called split-Bregman algorithm [START_REF] Goldstein | The split Bregman method for ℓ1-regularized problems[END_REF] but identical in principle to the Alternating Direction Method of Multipliers [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Setzer | Split Bregman algorithm, douglas-rachford splitting and frame shrinkage[END_REF]. This TV solver is implemented in C and may be called from Matlab. Finally, we also implemented the PW method in C. Computation times of the compared algorithms for the bunny and Buddha data sets are provided on Table 1. The three methods have large memory requirements when applied to the full grid. However, the graph cuts and power watershed banded methods are less memory intensive. For the bunny dataset, the touch-expand algorithm optimizing the Graph cuts energy only needs to allocate 3.6% of the full grid size. The solution space for the bunny dataset is much larger for our power watershed algorithm and reaches 33% of the full grid, because some points are widely spaced out. Thus, on this image, the touch expand algorithm is faster than the power watershed method. However, points clouds may require a more extensive expansion of the research area for the Graph cut touch expand method. For example, in the Buddha dataset, 8.75%

of the grid size is needed for computing the Graph cuts solution. In such cases, the power watershed approach is much faster. The TV approach has currently no guarantee to produce a global optimum when called on a banded graph, so the TV algorithm was implemented only on a full grid. On the Buddha dataset for example, given a 207 × 505 × 207 grid, the TV-based method requires 3G of RAM. 1. Timing experiments (in seconds) using an PC with a 3GHz Intel dual-core processor and 2G of RAM. The memory requirements are given in Mega Octets, and the computation times in seconds. Note that the seeds used are different for the three methods. For TV and PW, the seeds are imposed as hard constraints and are located far from the point cloud (see Fig. 5). The GC touch expand method uses soft constraints seeds computed from the normals given in input with the point cloud. We observe that Graph cuts are fast when the point cloud configuration does not require to expand the research area too much, like in the Bunny case (3.6%). However in the Buddha example, there is a need for a larger expansion, where the real complexity of the Graph cut algorithm become visible.

Point

Comparison with some other approaches

Following the quantitative information given in [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF], we compare the performances of power watershed for the Stanford bunny reconstruction with different methods [START_REF] Amenta | The power crust[END_REF][START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF][START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF][START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF]. All these methods, including the power watershed, reconstruct a surface close to the Bunny set of dots, with an error comprised between 2 × 10 -4 and 6 × 10 -4 . In terms of computation time, the Power Crust method [START_REF] Amenta | The power crust[END_REF] is about 10 times slower, and the Poisson method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] 3 times slower than our PW approach. Although the FFT [START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF], MPU [START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF] and Hoppe et al [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF] methods are fast, the FFT method suffers from large memory requirements limiting the grid resolution, and Hoppe et al and MPU methods produces artifacts in the presence of noise (See [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF]). The method of Jalba and Roerdink [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF], based on Coulomb potentials, uses 4 times less memory than our power watershed method, but is 5 to 10 times slower on a CPU and is still slower using a GPU. Furthermore our PW is more flexible in the choice of the markers. In our experiments, the amount of markers is not very large as shown in Fig. 5, but a strategy using larger seeds could be employed to reduce the size of the solution space and the computation time.

Conclusion

The power watershed method can be used to efficiently produce surfaces fitting noisy measurements. Contrary to standard watershed algorithms and the Graph cuts approach, the unique solution provided by the power watershed is not binary, resulting in the reconstruction of both smooth and detailed surfaces. In addition, this method is fast and not limited by large memory requirements, when using a restriction of the solution space. Finally, in comparison with other methods, our power watershed is robust to seed placement, and requires fewer parameters to be set. In practice, when using the power watershed method, close-fitting markers are not as mandatory as in other methods. Future work will follow several directions: memory and computation times improvements are still achievable, in particular by improving the touch-expand idea used in the Graph cuts optimization and adapting it for power watershed. The power watershed energy could also be modified to add some priors to the reconstructed surface, such as local orientation. Finally, we hope to demonstrate the efficiency of the power watershed technique for solving related problems such as multiview reconstruction.

Fig. 1 .

 1 Fig. 1. (a) Three dots in a 4 × 5 lattice. (b) Associated lattice weighted by a distance map from the dots according to (3) and squared, with Foreground and Background seeds. Note that only the ordering of the weights counts in the power watershed algorithm, thus the squared distance map can be used directly and produces the same solution as if the weights in (3) were used. (c) First steps of the power watershed algorithm to optimize (2) in the case q = 2 and p → ∞. Nodes having a maximum weight are merged. (d) A plateau of weight 2 (in green) including different seeded nodes is encountered. The Random walker algorithm is applied to label the nodes on the plateau. (e,f) New plateaus of weight 1 and 0 are encountered, the Random walker algorithm is applied, (g) Final labeling x solution of (2). The isocontour is represented in red.

Fig. 2 .

 2 Fig.2. Comparison of surface reconstruction from a set of points in 2D, using the total variation method, Graph cuts, a maximum spanning forest algorithm and finally the power watershed algorithm.

Fig. 3 . 7 × 10 - 3

 37103 Fig. 3. Power watershed results obtained from points clouds. (a) Total size of scans : 362272 points, Grid size : 234×297×301. (b) Number of scans used : 341072 points, Grid size : 275 × 276 × 668 (c) Total size of scans : 2748318 points, Grid size 382 × 270 × 171.

Fig. 4 .

 4 Fig. 4. Power watershed results obtained from noisy points clouds, corrupted by Gaussian noise of variance σ. The error was computed as the average distance between the obtained isosurface points to the original point cloud. The error is given in percentages of the diagonal of the bounding box of the data points.

Fig. 5 .

 5 Fig. 5. Robustness to seed quantity: this figure shows slices of solutions obtained for reconstructing the bunny surface, obtained on a grid of size 251 × 248 × 195, using different seeding strategies. The interior seed is colored in red, and the background seed in blue, and resulting isolines in green. (a,b,c) : Results obtained with TV minimization, and (d) Result obtained with Power Watershed.

Fig. 6 .

 6 Fig. 6. Grid size : 234 × 297 × 301. (a) Total variation minimization result, (b) Graph cuts result, (c) power watershed result. Isosurfaces at 0.5 have been extracted on all results, and were downsampled by 2 to render the surfaces.

  Table

		Method	Grid	Grid	Peak	Time	Time /
	cloud		size	used memory		Nb voxels
	Bunny	GC	234 × 257 × 301 3.6%	212	9	1.4 × 10 -5
	Bunny	PW	234 × 257 × 301 31%	1180	51 9.1 × 10 -6
	Bunny	TV	195 × 248 × 251 100%	900	122 1.0 × 10 -5
	Buddha	GC	276 × 671 × 277 8.7%	850	253 5.6 × 10 -5
	Buddha	PW	276 × 671 × 277 16%	1500	52 6.3 × 10 -6


Joël Van Neerbos
email: joelvneerbos@gmail.com

Laurent Najman
email: l.najman@esiee.fr

Michael H F Wilkinson
email: m.h.f.wilkinson@rug.nl

Towards a Parallel Topological Watershed: First Results

In this paper we present a parallel algorithm for the topological watershed, suitable for a shared memory parallel architecture. On a 24-core machine an average speed-up of about 11 was obtained. The method opens up possibilities for segmentation of gigapixel images such as found in remote sensing routinely.

Introduction

The watershed transformation is a popular tool for segmenting grayscale images, introduced by S. Beucher and C. Lantuéjoul [START_REF] Beucher | Use of watersheds in contour detection[END_REF]. It can be used to segment an image into regions with similar gray values. Due to ever increasing image sizes, several parallel algorithms have become available for different watershed paradigms [START_REF] Roerdink | The Watershed Transform: Definitions, Algorithms and Parallelization Strategies[END_REF]. In contrast, no parallel algorithms for the topological watershed [START_REF] Bertrand | On topological watersheds[END_REF][START_REF] Couprie | Topological grayscale watershed transformation[END_REF] are available. In the topological watershed framework, some of the grayscale information from the original image is preserved, which may be useful for further processing, such as reconnection of corrupted contours. Also, this grayscale information can be used to determine the significance of watershed lines [START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF].

In this paper we present a parallel algorithm for the topological watershed on shared-memory parallel machines. The algorithm is based on parallellizing each of the stages of the sequential algorithm [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF], and including a multi-pass stage such as used in most parallel algorithms for the regular watershed [START_REF] Roerdink | The Watershed Transform: Definitions, Algorithms and Parallelization Strategies[END_REF].

We will first describe the topological watershed briefly, after which we describe the parallel algorithm. Timing results are presented in Section 4, followed by the conclusions.

Topological Watersheds

The topological watershed [START_REF] Bertrand | On topological watersheds[END_REF][START_REF] Couprie | Topological grayscale watershed transformation[END_REF] was introduced to include grey level information in the end result, in such a way that the significance of each watershed line is preserved. Specifically, the topological watershed preserves the pass values [START_REF] Bertrand | On topological watersheds[END_REF][START_REF] Najman | Watersheds, mosaics and the emergence paradigm[END_REF], i.e. the highest altitude of the lowest path between any two minima.

The pass value of two points in the image is related to a concept of separation of the points. The points p and q are said to be k-separated if the following conditions apply: -There exists a path from p to q with maximum value k -1 -There exists no path from p to q with a maximum value lower than k -1 -Both p and q have a value lower than k -1 A path that satisfies the first two conditions for some k is called a lowest path in this paper. If a lowest path from p to q contains no value that is higher than both p and q, then p and q are not separated, but are linked. Separation is illustrated in Fig. 1. The top path in Fig. 1(b) connects two 5-separated pixels. The left path connects pixels that are not separated (but linked). The right path connects two 3-separated pixels.

In contrast to most watershed algorithms, topological watershed does not produce a binary image, but a gray-scale image. Intuitively, grey values of the pixels are obtained as follows:

-All pixels in a basin have the same gray value, namely the value of the minimum from the input image that is contained within the basin. -The values of the pixels on the watershed lines are as low as possible, without changing the separation relations between the basins. If two pixels from different basins were k-separated in the input image, they should still be k-separated in the topological watershed of the image.

The generic algorithm for computing a topological watershed of an image F proceeds by iteratively lowering some point satisfying a condition of destructibility until there is no more such points. More precisely, a point x is said to be W-destructible for F (where W stands for Watershed) if its altitude can be lowered by one without changing the number of connected components of the level set

F [k] = {p ∈ E; F (p) < k}, with k = F (x). A map G is called a W-thinning of F if it
may be obtained from F by iteratively selecting a W-destructible point and lowering it by one. A topological watershed of F is a W-thinning of F which contains no W-destructible point.

Both the watershed and a topological watershed of the image from Fig. 1(a) are shown in Fig. 2. Fig. 1(a)(c) shows that the separated pixel pairs from Figure 1(b) are still respectively 5-separated and 3-separated. The third pair is not shown as it was not separated.

Parallel Implementation

To obtain a fast algorithm, we need to lower a W-destructible point not by 1 as in the generic algorithm, but as much as possible. Following the theoretical results in [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF], the lowest value to which a point can be lowered can be computed from the Min-tree C(F ). This Min-tree, or component tree, can be obtained using the parallel algorithm from [START_REF] Wilkinson | Concurrent computation of attribute filters using shared memory parallel machines[END_REF]. Furthermore, the Min-tree allows one to easily check if two points are separated. More precisely, the separation is related to the altitude of the Lowest Common Ancestor (LCA) between two nodes of C(F ).

For efficiency, we need to perform preprocessing so that the LCA can be found in constant time. This is done using the algorithm in [START_REF] Schieber | On finding lowest common ancestors: simplification and parallelization[END_REF]. After this we need to lower all W-destructible pixels in the image to the value of their watershed basin. A pixel is W-destructible if its value can be lowered without linking two local minima into a single basin.

A key element of the algorithm is the function W-Destructible, which determines to which grey level a pixel can be lowered in the output. This uses the LCA value of neighbouring components of a the Min-tree. Because the function W-Destructible is a local function, we can parallelize the topological watershed algorithm for n threads simply by dividing the image into n tiles and assigning one tile to each thread. An example division for a 2D image is illustrated in Figure 3(a) and a division for a 3D volume is shown in Figure 3(d).

Problems arise when examining border pixels, because their neighbours in other tiles can be changed at any time by their assigned threads, producing incorrect results. We can solve these problems by letting each thread process its tile in different stages, and synchronizing all threads after each stage. Each tile is divided into sub-tiles, and a different sub-tile is processed in each stage. The tiles are divided in such a way that no two adjacent sub-tiles need to be processed at the same time.

Whether or not two sub-tiles are adjacent depends on the connectivity, as shown in Fig. 3. Fig. 3(a) shows an image divided into 12 tiles for 12 threads. With 4-connectivity, each tile should be divided into 2 sub-tiles, as shown in Fig. 3(b). The dark sub-tiles represent the sub-tiles that are processed in the first A single application of the multi-stage algorithm suggested above is however insufficient. For example, a certain pixel x may need to be lowered to the value of some local minimum y to obtain a topological watershed of the input image. However, if pixels x and y are part of different sub-tiles, it is possible that pixel x will not get the value of y the first time its sub-tile is processed. Multiple iterations may be needed to obtain the desired result, as in the case of regular watersheds [START_REF] Roerdink | The Watershed Transform: Definitions, Algorithms and Parallelization Strategies[END_REF].

The sequential topological watershed algorithm [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF] processes all pixels in an order that is determined by a priority queue, built as preprocessing stage. In the parallel implementation, this is done when the pixels are processed for the first time, but visiting each pixel in every later iteration is not necessary. Instead, we keep track of the pixels that have been changed, and only add those pixels to the priority queue that are adjacent to pixels in other sub-tiles that changed recently. For this purpose we will use the binary map pxChanged, that will store for each pixel whether or not its value has changed recently. The use of this map is illustrated in Fig. 4.

The parallel algorithm is distributed over two procedures: the procedure InitializeQueue , which corresponds to the first part of the sequential algorithm, and the procedure TopologicalWatershedTile. , which corresponds to the second part. Both procedures are given below. As its input, the InitializeQueue procedure needs the image F , the mintree C(F ) and the corresponding component map Ψ that associates to each point the node it belongs to in C(F ). Additionally, it also needs the sub-tile T in which it should operate, and the current state of the pxChanged map. If the pixel turns out to be W-destructible, it is added to the priority queue with its priority set to the level to which the pixel may be lowered. Also, this new level is stored in the map K and a pointer to the component to which the pixel may be added is stored in the map H. The pixel itself is not lowered yet. Thus, the output consists of the priority queue L, the maps K and H, and the updated pxChanged map. The queue L and the maps K and H are all local, but the map pxChanged is global, and may be read and modified by other threads while this procedure is being executed. However, each thread will only write in the part of the pxChanged map that corresponds to its current sub-tile, and will only read in adjacent sub-tiles that are processed in a different stage, so no conflicts emerge.

The algorithm starts by initializing the priority queue L in line 01. It then proceeds by setting the pxChanged map to false for every pixel in sub-tile T . If the procedure is run during the first iteration, then lines 04 to 08 are executed. In the first iteration, L is initialized exactly like in the original sequential algorithm, apart from the fact that only the pixels within T are processed instead of all pixels in F . If the procedure is called after the first iteration, lines 09 to 18 are executed. In these lines the algorithm checks all border pixels for changed neighbours, and tests the pixels for W-destructibility if any changed neighbours are found. If a pixel turns out to be W-destructible, it is added to the priority queue, and the maps K and H are updated as before.

The second procedure, TopologicalWatershedTile, needs the same input as the procedure InitializeQueue. The output consists of the updated image F , the update pxChanged map, and the binary variable anyChanges. For All border pixels p of T Do 11.

addP ← false 12.

For All neighbours q of p Do 13.

If pxChanged[q] = true Then addP ← true 14.

If addP = true Then 15.

c ← W-Destructible(F, p, C(F ), Ψ ) 16. If c = ∅ Then 17. i ← level of c; Li ← Li ∪ {p} 18. K(p) ← i; H(p) ← pointer to c
anyChanges is used to quickly determine if any changes have occurred in the sub-tile during the execution of this procedure. The procedure TopologicalWatershedTile starts by initializing the value of anyChanges. The function call on line 02 produces an initialized priority queue L, as well as initialized maps K and H. The rest of the function is mostly the same as the second part of the sequential algorithm. Line 08 is added, where the pxChanged map and the anyChanges variable are updated. Also, a new restriction is added to line 09, saying that only neighbours of p that lie within the sub-tile T should be added to the priority queue.

With the procedure TopologicalWatershedTile implemented, we can now define the main parallel algorithm: the procedure ParallelTW. As its input, it needs the pixel mapping F , the Min tree C(F ) and the corresponding component map Ψ . Because each thread will run this procedure independently, the global map pxChanged needs to be provided to each thread as well. However, no initial values need to be stored in it. Additionally, each thread is provided its identifier id. The first thread gets id value 0, the second gets value 1 and so on. The output of the procedure is the updated map F , that now contains the topological watershed of the input image.

The algorithm keeps looping until a topological watershed of the input image is found. In each iteration, the loop starting on line 02 is executed once by each thread. Line 04 starts a loop that visits all stages. The number of stages is equal to the number of sub-tiles assigned to each thread, as shown in Fig. 3. Line 05 then determines the location and dimensions of the sub-tile to be processed by the current thread in the current stage. Some examples of the sub-tiles that While ∃p ∈ L k Do 05.

L k = L k \{p} 06. If K(p) = k Then 07. F (p) ← k; Ψ (p) ← H(p) 08.
pxChanged[p] ← true; anyChanges ← true 09.

For All neighbours q of p within T , with k < F (q) Do 10.

c ← W-Destructible(F, q, C(F ), Ψ ) 11.

If c = ∅ Then K(q) ← ∞ 12.

Else 13.

i ← level of c 14.

If K(q) = i Then 15.

Li ← Li ∪ {q}; K(q) ← i 16.

H(q) ← pointer to c should be processed by each thread in the first stage are displayed in Fig. 3. In the other stages the sub-tiles that are processed should have a similar pattern, always assuring that no two adjacent sub-tiles are processed at the same time.

The topological watershed of the tile is then computed on line 06. If the algorithm returns that there have been some changes, then this is stored for the current thread in the (global) anyChangesThr array. After this, a standard barrier function is called, that just waits until all threads have reached this barrier and then lets all thread continue. This is done to ensure that no thread will start with the next stage until all threads are done with the current one. When all threads have finished processing all their sub-tiles, the first thread will check if there have been any changes in any of the threads. If there have not been any changes at all, a topological watershed has been found and all threads will terminate. Otherwise, each thread will go to the next iteration by starting again with the main loop from line 02.

Min-Tree Compression

Wilkinson et al. [START_REF] Wilkinson | Concurrent computation of attribute filters using shared memory parallel machines[END_REF] described how to parallelize the computation of a min-tree.

Because the algorithm from [START_REF] Wilkinson | Concurrent computation of attribute filters using shared memory parallel machines[END_REF] also deals with features of the min-tree that we won't use, its implementation is simplified somewhat. Basically, the sequential algorithm is parallelized by letting multiple threads each compute the min-tree of a different part of the input image, and merging the min-trees of the parts afterwards. This algorithm uses a representation of the min-tree which is as large as the image or volume itself, i.e. each pixel or voxel is a node, containing a pointer to its parent. Only those nodes which have a parent with grey level greater than its own are relevant to either filtering or the watershed computation. For All stages s Do 05.

T ← current sub-tile, based on id and s 06.

TopologicalWatershedTile(F, C(F ), Ψ , T , pxChanged) 07.

If anyChanges = true Then anyChangesThr[id] ← true 08.

Barrier() 09.

If id = 0 Then 10.

anyChangesAtAll ← false 11.

For All threads t Do 12.

If anyChangesThr[t] = true Then 13.

anyChangesAtAll ← true 14.

Barrier() 15.

If anyChangesAtAll = false Then done ← true

These nodes are called level roots. When the min-tree is built in parallel, not all pixels of a given min-tree node directly point to the level root. This means finding a level root is costly. Because the LCA algorithm inspecs the level roots often, we need to compress the tree in the sense that all parent pointers always point to a level root, yielding what is refered to as a canonical representation of the tree [START_REF] Najman | Building the component tree in quasi-linear time[END_REF]. This reduces the computation time of the LCA algorithm. The algorithm is shown in Alg. 4.

The procedure CompressTree needs a map F , a component tree C(F ) and a corresponding component map Ψ as its input, as well as the identifier id of the thread that executes it. The output of the procedure consists of the compressed component tree and map. Function LevelRoot2is used by CompressTree to obtain the level roots, without ever writing in memory sections not assigned to the current thread.

Algorithm 4

The Min-tree compression algorithm.

Function LevelRoot2 (Input c, F , C(F )) 01. If c = root(C(F )) ∨ F (c) = F (parent(c)) Then 02. Return c 03. Else 04. Return LevelRoot2(parent(c), C(F ), F ) Procedure CompressTree (Input F , C(F ), Ψ , id; Output C(F ), Ψ ) 01. For All components c ∈ segment id of C(F ) Do 02. parent(c) ← LevelRoot2(parent(c), F , C(F )) 03. For All pixels p ∈ segment id of F Do 04. Ψ (p) ← LevelRoot2(Ψ (p), F , C(F )) (a) (b) (c) (d)

Results

The implementation was tested on the following four different input images: (i) a 4000 × 4000 satellite image of an airfield, (ii) a 4000 × 4000 image in which each pixel has a random gray value (iii) a 1000 × 1000 image with a spiral shaped plateau, and (iv) an angiogram with dimensions 256 × 256 × 128. The spiral image was chosen as a (near) worst-case scenario, in the same way as in [START_REF] Roerdink | The Watershed Transform: Definitions, Algorithms and Parallelization Strategies[END_REF]. The random noise image was another extreme case. The images are displayed in Fig. 5. All data sets were tested on a 4-socket, 6-core per sockets, opteron-based machine with 128 GB of memory, using 1 to 24 threads, and various connectivities (4 and 8 in 2D, 6 and 26 in 3D). The memory is divided into banks assigned equally to each processor sockets, but is accessible as shared memory amongst all cores. Accessing memory of another porcessor socket does incur a slight speed penalty. In all cases, each thread was initially assigned an equal slice of the image or volume during the min-tree construction phase, as in [START_REF] Wilkinson | Concurrent computation of attribute filters using shared memory parallel machines[END_REF], and an equal tile as described in section in the remainder of the algorithm as detailed above.

The previously existing sequential algorithm by Couprie et al. [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF] and the newly implemented parallel algorithm were run on the same machine with the same input (the satellite image, with 4-connectivity using 1 thread) to compare the two versions on a single core.

On the satellite image, we obtained a wall-clock time of 86.36 s at 4 connectivity, for a single thread. This dropped to 17.14 s at 8 threads, 9.71 s at 16 threads, and 7.17 s at 24 threads. At 8 connectivity, all times rose, and the timings were 135.17 s at 1 thread, decreasing in a similar fashion to 12.04 s at 24 threads. The noise image (at 4 connectivity) showed a very similar pattern, decreasing from 112.83 s on 1 thread to 11.95 s on 24 threads. The influence of connectivity was most profound in the 3D case. At 6 connectivity, timings run from 73.18 s at 1 thread to 5.7 s at 24 threads. At 18 connectivity these figures rise to 214.71 s and 16.62 s respectivily. increasing further to 359.22 s and 33.37 s at 26 connectivity. As expected, the (smaller) spiral image behaved very differently, with timings of 4.29 s on 1 thread, 1.71 s on 8, 1.58 s on 16, and increasing again at 24 threads to 1.99 s. This is due to the increase in the number of iterations of the topological watershed algorithm.

The overall behaviour is shown in Fig. 6. In Fig. 6(a) the diagonal line shows the ideal speedup, the line at the bottom shows the speedup for the spiral image. The line through the gray area shows the speedup for the satellite image with 4-connectivity, which is computed for each number of threads individually. The grey area represents the results of the all tests, excluding for the spiral image. The top and bottom of the grey area are defined as the average speedup of their stage plus and minus the standard deviation, respectively. The diagonal lines show the ideal speedup, where the speedup is equal to the number of threads.

In Fig. 6(b) we show the contributions of the four stages to the total wallclock time of the algorithm. The bottom layer represents the time consumed by the construction of the min-tree, the second layer from the bottom represents the time taken to compress the min-tree, the third layer shows the time it takes to perform the preprocessing for the LCA, and finally the remaining layer on top represents the time consumed by the final stage which produces the topological watershed.

The individual speedups are shown in Figure 7. The speedup of each stage is represented by a black line and a grey area, where the line shows the speedup of that stage when computing the satellite image with 4-connectivity. The corresponding grey area represents the results of the remaining tests as in Fig. 6, again excluding the spiral image.

The construction of the component tree and the final stage that computes the topological watershed parallelize quite well, while the tree compression and the LCA preprocessing parallelize rather poorly. Fortunately, the tree compression takes up only a small percentage of the total wall-clock time, even when 24 threads are used. However, this percentage will probably increase when more threads are used. The LCA has a larger impact on the total wall-clock time of the algorithm, and will cause the speedup of the total algorithm to decrease more severely when using a larger number of threads.

Note however, that all speedups reported are relative to the wall-clock times of the parallel algorithm using one thread. In practice, the previously existing sequential implementation performs about 1.5 times faster than parallel algorithm when only one thread is used. On two threads, the current algorithm has the edge, ableit by a small margin. Further optimizations in the parallel algorithm, for example by using a better tree compression method, may reduce this difference.

Conclusions

This paper described a way to parallelize the computation of the topological watershed. An implementation that was created according to this description showed that a reasonably good speedup could be achieved while using up to 24 threads, and the trend in the results suggests that even better speedups may be achieved when more than 24 threads are used.

The parallel implementation may be improved further by creating a more compact min-tree in the min-tree construction or compressing it more in the tree compression stage. For example, each node in the min-tree that has only one child can be merged with that child, setting all pixels belonging to its component to the gray level of its child, and setting the parent of its child to be the parent of the node itself. This simple change could significantly reduce the LCA preprocessing time, and may also improve the wall-clock time of the last stage where the topological watershed is computed.

Another improvement could be to use a more sophisticated method of parallel list ranking. Better performing algorithms already exist, the algorithm in this paper was used for simplicity reasons only. A list ranking algorithm that parallelizes better and has higher speedups, may significantly reduce the total wall-clock time of the algorithm, especially when using larger numbers of threads.

Furthermore, the image is now divided into tiles which are each processed by their own thread. This division is only based on the image dimensions and the number of threads, not on the contents of the image. Taking into account the contents of the image, maybe in combination with the min-tree and component map, while dividing the image among the threads may reduce the communications needed between the different threads, which could lead to a faster algorithm. We can also try other ways of tiling the image, following [START_REF] Matas | Parallel Algorithm for Concurrent Computation of Connected Component Tree[END_REF].

In short, there is still room for improvement in the parallel algorithm proposed in this paper, but in its current form it can already be used to greatly speed up the computation of the topological watershed, compared to the previously existing sequential algorithm.

Fig. 1 .

 1 Fig. 1. Separation. (a) a digital grey-scale image; (b) lowest paths between three pairs of pixels.

Fig. 2 .

 2 Fig. 2. Watershed and topological watershed: (a) watershed of Fig. 1(a); (b) topological watershed of Fig. 1(a); (c) lowest paths between the separated pairs of pixels from Fig. 1(b).

Fig. 3 .

 3 Fig. 3. Tiling: (a) 2D image; (b) tiling for 4-connectivity; (c) same for 8-connectivity; (d) 3D volume; (e) tiling for 6-connectivity; (f) same for 18-connectivity and (g) 26connectivity;

  (e), (f) and (g). Again, the dark sub-tiles represent the sub-tiles to be processed in the first stage.

Fig. 4 .

 4 Fig.4. An example of a pxChanged map in an iteration after the first. The white pixels in the four adjacent sub-tiles are pixels that are marked as 'changed' in the pxChanged map. The pixels in the current sub tile that are adjacent to such a changed pixel, marked with x in this figure, are added to the priority queue of the sequential algorithm.

The variable Algorithm 1

 1 Queue Initialization Procedure InitializeQueue (Input F, C(F ), Ψ , T , pxChanged; Output L, K, H, pxChanged) 01. For k From kmin To kmax -1 Do L k ← ∅ 02. For All pixels p ∈ T Do pxChanged[p] ← false 03. If first iteration Then 04. For All pixels p ∈ T Do 05. c ← W-Destructible(F, p, C(F ), Ψ ) 06. If c = ∅ Then 07. i ← level of c; Li ← Li ∪ {p} 08. K(p) ← i; H(p) ← pointer to c 09. Else 10.

Algorithm 2

 2 The algorithm for a single tile Procedure TopologicalWatershedTile (Input F, C(F ), Ψ , T , pxChanged; Output F , pxChanged, anyChanges) 01. anyChanges ← false 02. InitializeQueue(F, C(F ), Ψ , T , pxChanged) 03. For k From kmin To kmax -1 Do 04.

Algorithm 3

 3 The parallel topological watershed algorithm Procedure ParallelTW (Input F, C(F ), Ψ , pxChanged, id; Output F ) 01. done ← false 02. While not done Do 03. anyChangesThr[id] ← false 04.

Fig. 5 .

 5 Fig. 5. The four test input images. (a) shows a satellite image, (b) an image with random pixel values, (c) an image with a spiral-shaped plateau and (d) an angiogram, which is a 3D volume.

Fig. 6 .

 6 Fig. 6. The overall behaviour of the algorithm: (a) The total speedups in the performed tests, relative to the wall-clock time for a single thread with the same input; (b) the contributions of the four stages to the total wall-clock time of the algorithm, when computing the topological watershed for the satellite image with 4-connectivity.

Fig. 7 .

 7 Fig. 7. The speedups of the four stages of the algorithm: (a) Top: min-tree construction, bottom: tree compression; (b) Top: topological watershed computation, bottom: LCA preprocessing.