Watersnakes: energy-driven watershed segmentation

330 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 25, NO. 3, MARCH 2003 Watersnakes: Energy-Driven Watershed Segmentation Hieu Tat Nguyen, Marcel Worring, and Rein van den Boomgaard Abstract—The watershed algorithm from mathematical morphology is powerful for segmentation. However, it does not allow incorporation of a priori information as segmentation methods that are based on energy minimization. In particular, there is no control of the smoothness of the segmentation result. In this paper, we show how to represent watershed segmentation as an energy minimization problem using the distance-based definition of the watershed line. A priori considerations about smoothness can then be imposed by adding the contour length to the energy function. This leads to a new segmentation method called watersnakes, integrating the strengths of watershed segmentation and energy based segmentation. Experimental results show that, when the original watershed segmentation has noisy boundaries or wrong limbs attached to the object of interest, the proposed method overcomes those drawbacks and yields a better segmentation. Index Terms—Watershed segmentation, energy-based segmentation, topographical distance, snakes. æ 1 S INTRODUCTION EGMENTATION is a fundamental problem in image analysis. It should yield a partitioning of the image into disjoint regions, uniform according to some feature such as gray value, color, or texture. The segmentation process can rely both on the uniformity of the feature within the regions or on edge evidence. In both cases, the result should be a balance between adherence to the possibly noisy and incomplete data and smooth segmentation results suited for further analysis. There are two major approaches in segmentation: energybased and watershed-based. In the first approach, the segmentation is obtained as result of the minimization of a so-called energy function. The second approach is mainly based on the watershed algorithm from mathematical morphology. This paper investigates the relation between these approaches. A large number of segmentation methods in literature use the first approach [1], [2], [3], [4], [5]. To balance data adherence and smoothness, the energy is composed of a data-driven term and a regularization term. The purpose of the second term is to impose a priori knowledge, usually smoothness of the region contour, on the segmentation result. The data-driven terms used in literature can be classified into two classes: contour-based (or snake-based) and region-based. Snake-based methods, as originally defined by Kass and Witkin [4], use gradient information as input data. A contour, balancing maximal edge evidence with minimal curvature and curve length is found. The contour representation used in this method, may encounter problems with topological changes, like self-intersections [6], during optimization. To avoid these problems, Osher and Sethian proposed a level set approach [6]. The contour is represented as a level curve of a . The authors are with the Intelligent Sensory Information Systems Group, University of Amsterdam, Faculty of Science, Kruislaan 403, NL-1098 SJ, Amsterdam, The Netherlands. E-mail: {tat, worring, rein}@science.uva.nl. Manuscript received 15 Nov. 2000; revised 8 Sept. 2001; accepted 14 Aug. 2002. Recommended for acceptance by L. Vincent. For information on obtaining reprints of this article, please send e-mail to: tpami@computer.org, and reference IEEECS Log Number 113145. 3D surface or boundary of a growing region. As the representation is intrinsic, topological changes can be handled. The level set method has also been used in the geodesic model of Caselles et al. [5] and the bubble model of Tek and Kimia [7]. Other common issues that need to be addressed in the snake-based methods are: the decision when to stop moving the contour, overcoming gaps between broken edges, and the selection of the initial contours. The geodesic model deals well with gaps, while the bubble model offers a solution for initial contour selection by randomly initializing a large number of seeds that grow and merge afterwards. In both models, the criterion to stop contour motion, however, depends on a predefined parameter. As a consequence, the contour may converge to nonsignificant edges. The region-based methods use a global energy function, computed for the entire area of the regions, rather than their boundaries only. The earliest models are the Bayesian segmentation of Geman and Geman [8] and the one proposed in Mumford and Shah [1]. In [3], Leclerc suggests a segmentation approach based on the minimum description length criterion. In [2], Zhu and Yuille considered this criterion in the continuous domain and proposed a segmentation method named region competition. In all existing region-based methods, the uniformity in a region is described using a parameterized probability distribution of intensity. The energy yields a trade-off between how well the models describe the region and the smoothness of the contours. Since parameters of the distributions are unknown, the algorithm needs to go back and forth between the determination of the region parameters and the determination of the region labels. The methods, therefore, are complex and computationally expensive. Apparently, the watershed algorithm from mathematical morphology takes a very different approach, compared to the energy-based methods described. The input is a relief function representing edge evidence, where the morphological gradient is a common choice for computing such a relief. By viewing this function as a mountain landscape, object boundaries are determined as watershed lines. In [9], [10], the watershed algorithm is implemented via region growing, where seeds are the regional minima of the 0162-8828/03/$17.00 ß 2003 IEEE Published by the IEEE Computer Society Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION relief. To prevent oversegmentation due to a possibly large number of minima in the original edge evidence function, the watershed algorithm from selected markers is usually employed [9], [11]. This technique selects markers from significant minima and modifies the relief accordingly such that nonsignificant minima are filled up. The markers are usually extracted as low gradient zones in the image, smoothed by a morphological filter such as opening by reconstruction [12]. The watershed algorithm has been proven powerful for contour detection [13], [14] as well as image segmentation [9], [11]. Compared to the other methods mentioned, the watershed has several advantages, including the proper handling of gaps and the placement of boundaries at the most significant edges (see Section 4). It is difficult, however, to impose a priori information, in particular smoothness, on the watershed lines. As the algorithm relies purely on the edge evidence, the boundaries between regions may be noisy. A simple postprocessing smoothing of the watershed line does not take into account the observed information and, therefore, important corners with high edge evidence may be lost. So, the question is whether we can represent the watershed segmentation as the result of the minimization of an energy function. If so, a priori information can be imposed by adding appropriate regularization terms to the energy. This is achieved in this paper, leading to a new method called watersnakes. The paper is organized as follows: Section 2 provides the necessary background on watershed segmentation. In Section 3, we derive the cost function, whose minimization is equivalent to the watershed segmentation. We then add the contour length to this function to produce a segmentation with smooth boundaries. Section 4 derives properties of watersnakes, including a detailed comparison with the above energy-based methods. Section 5 develops algorithms for the implementation. Finally, Section 6 shows experimental results. 2 DEFINITION AND PROPERTIES OF THE WATERSHED Several distance-based definitions of the watershed have been proposed by Meyer in [15], Najman and Schmitt in [13], and Preteux in [16]. Although the definitions are slightly different, all papers show that, with an appropriate measure of the distance traveled over the relief, the watershed line is the set of points equidistant to the regional minima of the relief function. The definitions are valid for both the continuous and the discrete case. In this section, we briefly consider the definitions and results needed for this paper. 2.1 The Topographical Distance The heart of the distance-based definition of the watershed is the concept of topographical distance. Suppose the watershed is defined for a given relief function fðxÞ : X ! 7 IR on some domain X  IR2 . Definition 1. For any smooth function f, the topographical distance between two points x and y is a spatial distance weighted with the gradient norm jrfj: 331 Lðx; yÞ ¼ inf > e yŠ 2½x Z jrfð ðsÞÞjds; ð1Þ where ½xe>yŠ denotes the set of all possible paths from x to y. As indicated, the above definition is valid for smooth functions only. When the function is nonsmooth, the topographical distance is defined in a more complicated way [15]. For the case where f is defined on a digital grid: X  ZZ2 , a discretized version of (1) has been proposed by Meyer in [15]: L~ðx; yÞ ¼ min  n X rðtiÿ1 ; ti Þ:distðtiÿ1 ; ti Þ; ð2Þ i¼2 where distðÞ denotes the Chamfer distance [17]. The minimum is taken over all possible paths of adjacent pixels  ¼ t1 ; t2 ; . . . ; tn from x to y, where t1 ¼ x, tn ¼ y, and ti ; tiþ1 are in the same Chamfer neighborhood. Here, rðtiÿ1 ; ti Þ is an approximation of the gradient norm. Based on the work by Meyer in [15], we propose to compute rðtiÿ1 ; ti Þ as follows: 8 if fðtiÿ1 Þ < fðti Þ < LSðti Þ if fðtiÿ1 Þ > fðti Þ ð3Þ rðtiÿ1 ; ti Þ ¼ LSðtiÿ1 Þ : minfLSðti Þ; LSðtiÿ1 Þg if fðtiÿ1 Þ ¼ fðti Þ; where LSðxÞ is the lower slope at x [15]. The definition is identical to the one of Meyer in [15] except when fðtiÿ1 Þ ¼ fðti Þ. As we have shown in [18], (3) is advantageous over the definition in the reference in cases where the latter leads to the shifting of the watershed line. 2.2 The Watershed Line Now, suppose the function f has a finite number of distinct regional minima [9], denoted M1 ; . . . ; MK . Let i be the level of f on Mi . For each regional minimum, we consider the corresponding topographical distance transform: Li ðxÞ ¼ Lðx; Mi Þ ¼ inf Lðx; yÞ: y 2 Mi ð4Þ For each point, K values of distance to the regional minima are obtained. The catchment basins and the watershed line are then defined as: Definition 2. The catchment basin CBi of a regional minimum Mi is defined by: CBi ¼ fx 2 X j 8j 6¼ i; 1  j  K : i þ Li ðxÞ < j þ Lj ðxÞg: ð5Þ Definition 3. The watershed line of the function f is the set of points not belonging to any catchment basin: [ ð6Þ W SðfÞ ¼ X n CBðMi Þ: i In case all regional minima are at the same level, the watershed is the topographical SKIZ (skeleton by zones of influence [19]) of the regional minima. It is important to note that the function f can be fully reconstructed from the distances to the regional minima as follows: [20], chapter 5], [21]: fðxÞ ¼ min f 1iK i þ Li ðxÞg: Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. ð7Þ 332 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 25, NO. 3, MARCH 2003 In [16], CðÞ is called the connection cost of . Here, we show that it can be used to define the reconstructed relief. Proposition 1. g1 ðxÞ can be interpreted as the maximal level the water from the markers has to reach before flooding x: g1 ðxÞ ¼ ð10Þ min CðÞ; > e Sm Š 2½x where ½x e>Sm Š denotes the set of paths, thus linking x and Sm . Fig. 1. Illustration of the watershed line and topographical distance in the one-dimensional case. As illustrated in Fig. 1, i þ Li ðxÞ equals fðxÞ within the catchment basin CBi . When the pixel moves beyond the catchment basin, for the 1D case i þ Li ðxÞ is the reflection of fðxÞ with respect to the horizontal line passing through the watershed point. For the 2D case, a similar behavior of i þ Li ðxÞ can be observed, except that the reflection is far more complicated as the watershed line has varying height. The watershed line may be thick, i.e., have a nonzero area. It may also have so-called barbs which are branches of zero area with an end point [13]. Several approaches have been proposed to define a thin watershed line [11], [16], [22]. These methods add the geodesic distance to the topographical distance. This makes the definition much more complicated. In addition, we notice that the geodesic SKIZ itself may also be thick (see Appendix A). As a consequence, the referenced approaches yield a thinner watershed line, compared to the one in Definition 3, but they cannot guarantee a watershed line of zero area. 2.3 The Watershed from Selected Minima In practice, to prevent oversegmentation, the watershed line is constructed from a given set of regions called markers. A modification of the relief function is required [11]. First, all points in the markers are given a lowest value, say, 0. A new relief is constructed by the recursive conditional erosion: gnþ1 ¼ maxff; "gn g; ð8Þ where g0 ðxÞ is the function, which has the value 0 on the markers and 1 otherwise, and "gn denotes the erosion of gn with a minimum disk. The watershed transformation is then applied to g1 which has the markers as its sole minima. The classical watershed line of g1 is a subset of the watershed line of the original function f with the most significant edges retained. We add another interpretation of g1 that has a generic property with respect to (8). Let Sm be the set of markers. Given a path  from x to Sm . Let CðÞ ¼ max fðzÞ: z2 ð9Þ The proof is given in [18]. Another point of view on g1 , which also does not require recursion, has recently been given by Meyer in [23]. We have described in this section the distance-based definition of the watershed line. We have modified the definition of the topographical distance in a discrete grid [15]. The modified definition guarantees the placement of the watershed line at ridgelines of the relief and the possibility to recover the relief from the topographical distance to regional minima. We also give a new and more generic definition for the reconstructed relief in the watershed algorithm from markers. This new definition does not require the recursive erosion as in the traditional approach [11]. The definition of the watershed, based on the topographical distance, is needed in the remaining part of the paper for the investigation of the relation between the watershed and the energy minimization. 3 WATERSHED AS A MINIMIZATION PROBLEM In this section, we establish the relation between watershed and energy-based minimization. 3.1 The Energy Function Unless stated otherwise, the notations in this subsection are used for both the continuous and the discrete case. Let us first define segmentation as a partition of the image space: Definition 4. A partition of the image space X is a set of connected regions 1 ; . . . ; K such that aÞ K [ i ¼ X; i¼1 bÞ 8i 6¼ j : cÞ Intð iÞ i ¼ \ j ¼ ;; ð11Þ i; where IntðRÞ denotes the interior of set R and R its closure. The last condition is meaningful only in the continuous case. It ensures that the boundaries of the regions do not have barbs. This is equivalent to saying that each boundary segment must separate at least two different regions. Definition 5. A partition segmentation if 8x 2 i : i 1; . . . ; K þ Li ðxÞ ¼ min f 1jK is called a watershed j þ Lj ðxÞg: ð12Þ Considering (7), the right-hand side of (12) actually equals fðxÞ. It follows from Definition 5 that each region of a watershed segmentation contains one and only one catchment basin and the borders of the regions lie within the possibly thick watershed line (see Fig. 2). Note that the inverse statement is Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION 333 In conclusion, Theorem 2 states the equivalence of the watershed to energy minimization. Fig. 2. (a) Illustration of a watershed segmentation. (b) Illustration of a partition that is not a watershed segmentation. true only in case of a thin watershed line. When the watershed is thick, one can find a partition for which the region boundaries lie in the watershed region but (12) does not hold. An important question is whether a watershed segmentation actually exists. The existence is obvious in case of thin watershed lines. However, when the watershed region is thick, it is not clear how to assign pixels in the watershed region to the catchment basins such that Definitions 4 and 5 are satisfied. To this end, we have proven the following. Proposition 2. Every catchment basin CBi is a connected set. See [18] for a proof. This proposition is then needed for the proof of the following theorem. Theorem 1. For a smooth function f, defined on an open, bounded, and connected domain X , there always exists at least one watershed segmentation. Proof. See Appendix A. The main idea is to assign all pixels in the watershed region to one catchment basin unless (12) or the requirements of a partition are violated. u t Although we prefer to accept the fact that the watershed line can be thick, the proof of Theorem 1 actually specifies a way to construct a watershed line, which has zero area, no barbs, and at the same time, satisfies (12). As noted in the previous section, none of the existing approaches can guarantee a watershed segmentation with an ideal thin watershed line. Given Theorem 1, we can now take a next step by proving: 1; . . . ; Theorem 2. A partition function: Eð 1; . . . ; KÞ ¼ K ZZ X i¼1 K f minimizes the following i þ Li ðxÞgdx ð13Þ i if and only if it is a watershed segmentation. The proof is given in Appendix B. Looking at (13), we see that in case i ¼ 0, EðÞ is actually the sum of the volumes of the topographical distance function over the regions. We remark that, in fact, SKIZ with respect to any distance measure can be obtained by minimizing a cost function similar to (13). In our case, the minimization of E yields the SKIZ with respect to topographical distance, which is the watershed line. We call the function E energy as this is the common term for a cost function whose minimization yields a desired segmentation. 3.2 Watershed and PDEs—Imposing Smoothness A PDE form of the watershed line can be obtained by considering the minimization of the energy E in (13) along the steepest descent direction. RR f iþ The differentiation of the area functional i Li ðxÞgdx with respect to a deformation of the border @ i @ of i results in the equation @x ni , where x@ @t ¼ f i þ Li ðx@ Þg~ is a point on @ i and ~ ni is the unit normal vector of @ i pointing inwards into i [2]. Thus, the value i þ Li plays the role of a “force” compressing the region i from its boundary. A boundary point adjacent to two regions i and ni and ð j þ Lj Þ~ nj . Since j is under two forces ð i þ Li Þ~ ~ ni , the motion is given by nj ¼ ÿ~ @x@ ¼ fð @t i ÿ jÞ þ ðLi ÿ Lj Þg~ ni : ð14Þ Thus, x@ moves as long as it still resides inside one of the 6 0 and stops catchment basins where ð i ÿ j Þ þ ðLi ÿ Lj Þ ¼ when this value equals zero. In the end, the boundary resides in the watershed line. Having represented the watershed segmentation as an energy minimization, the smoothness of the region boundary is now achieved by adding the boundary length to the energy function 0 1 Z ZZ K X B C dsA; ð15Þ f i þ Li ðxÞgdx þ Eð 1 ; . . . ; K Þ ¼ @ i¼1 @ i i where is the weighting coefficient. We call the segmentation method based on the minimization of this energy function watersnake in order to distinguish it from the original watershed. In this case, the compressing force acting on a point x@ at the border of a ni , where i denotes the region i is ð i þ Li þ i Þ~ curvature of @ i at x@ . When two regions i and j are neighbors, we have i ¼ ÿj ¼  on their common boundary. The motion of the boundary is then caused by the vector sum of the two compressing forces @x@  ¼ ð @t i ÿ jÞ þ ðLi ÿ Lj Þ þ 2  ~ ni : ð16Þ Some other PDE-based formulations for the watershed line have been used in literature. In [7], Tek and Kimia notice that the watershed algorithm on a relief function f 1~ @ can be implemented by the evolution equation @x @t ¼ f n. In [24], Maragos and Butt give another equation: @x@ 1 ~ @t ¼ jrfj n. Both these equations describe the evolution of marker contours before the markers meet, i.e., when the contours reside within the corresponding catchment basins only. Our PDE (14), in contrast, describes the motion of an arbitrary contour toward the watershed line. The formulation in [7] and [24] can be used for computing the original watershed line. However, in case the evolution equation also contains the regularization term, these approaches are not appropriate. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. 334 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 25, NO. 3, MARCH 2003 Fig. 3. Classification of segmentation methods. 4 COMPARISON OF WATERSNAKES WITH ENERGY-BASED SEGMENTATION METHODS In this section, we compare watersnakes with other energybased segmentation methods. It is clear from (15) that the energy function of watersnakes belong to the region-based group (see Fig. 3). In particular, it can be compared to the energy used by Zhu and Yuille in [2]. The latter energy is defined as follows: 1 0 Z Z Z K XB C ð17Þ dsA; ÿ log P ðIx jai Þdx þ min @ i¼1 i @ i where P ðIx jai Þ denotes the probability density of the measured intensity at pixel x under the assumption that the pixel belongs to the region i and ai denotes the parameter vector of this distribution. It can be observed that the energy function of watersnake in (15) and the energy function of the region competition method in (17) are similar except that the topographical distance i þ Li ðxÞ in (15) replaces the term ÿ log P ðIx jai Þ in (17). Both these terms measure the homogeneity within region i but the topographical distance has the major advantage that it does not contain any parameters. The watersnake does not have to go back and forth between updating region parameters and updating pixel’s labels. The minimization procedure is therefore simpler than those in the region competition and other methods using region-based energy. The watersnake model has also advantageous properties over the snake-based methods: 1. While the snake-based methods [4], [5], [7] might converge to weak edges, depending on the tuning of parameters, the watershed line always corresponds to the most significant edges between the markers. When the smoothing term is present, the watersnake deviates from the watershed line to reduce the local curvature until ð i ÿ j Þ þ ðLi ÿ Lj Þ is compensated with 2 . The smoothing effect is inversely proportional to the sharpness of edges since the value of ð i ÿ j Þ þ ðLi ÿ Lj Þ increases faster for sharp edges when the contour moves away from the watershed line. 2. In case evidence for edges between the markers is low, the resulting contours still lie between the markers. The watersnake is therefore able to fill the gaps between broken edges. In conclusion, the watersnake is a fusion of energy-based segmentation and morphology-based segmentation and combines the advantages of these two approaches. 5 IMPLEMENTATION There are two possible approaches for implementing the minimization of the energy (15) in the discrete domain. The first approach takes the motion (16) as its starting point. Discretization of this equation requires a discrete estimator of curvature. The second approach discretizes the energy in (15) and defines a method for minimizing it. This approach requires a discrete estimator of contour length. These two approaches are investigated in detail in the two sections that follow. 5.1 Algorithm Based on Region Growing For the discretization of PDEs of contour motion depending on curvature, the well-known level set method of Osher and Sethian in [6] can be applied. This method is appropriate for solving (16) in case of two markers. The method would, however, grow expensive in case of multiple markers. We have to keep track of every boundary between every pair of potentially adjacent regions each with its own PDE. At the same time, we would have to take care that the topology of the entire partition is preserved. A simpler approach is to modify the region growing step of the original watershed algorithm such that the region borders are kept smooth during growth. Starting from a set of markers, the algorithm iteratively adds to the markers the pixels connected to their border. For the pixel selection, the growing process takes into account the value of the compressing force i þ Li þ i at each boundary pixel. The part of the border, where this force is weakest, is Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION 335 Although these coefficients are obtained with circles, they also yield a good estimate for straight boundaries, as indicated in Fig. 5. Note that the above contour length estimator will not be accurate for small-size regions like a single pixels since it is derived for disks of larger radius. 5.3 Algorithm Based on Energy Discretization In this section, we propose an algorithm that aims to minimize the following discrete version of the energy (15) Fig. 4. For the indicated pixel: n4 ðxÞ ¼ 2; n8 ðxÞ ¼ 2; nk ðxÞ ¼ 5. E~ð propagated in the outward direction first. The curvature i at a pixel is computed from the local boundary pixels [25]. Every time a pixel is added to the markers, the curvature at the boundary pixels in the neighborhood is recomputed. A similar idea has been proposed in the segmentation methods of [26] and [27]. In [26], the criterion for adding a pixel to a marker is based on both contrast and contour complexity. The latter term is defined as the number of contour points added to the marker contour when the pixel in consideration is added to the marker. This approach is then extended in [27] for region growing applied to flat zones. Note that, although the intuitive approach of keeping region contours smooth during growth yield smooth resulting boundaries in many practical situations, it does not guarantee to yield a solution of (16). As a consequence, we do not consider this algorithm any further in the paper. We now concentrate on (15). 5.2 Intermezzo: Local Contour Length Estimate Before describing in the next section an algorithm minimizing the energy (15), we need a contour length estimator. We estimate the contour length using the number of pairs of neighboring pixels. In this way, the contour length can be re-estimated locally when the contour moves. We use the Chamfer neighborhood. Given pixel x 2 i , let n4 ðxÞ; n8 ðxÞ, and nk ðxÞ be the number of pixels which are in the four, eight, and “knight-move” connected neighborhood of x but outside i (see Fig. 4). The perimeter of i is then estimated as: Z X ½w4 n4 ðxÞ þ w8 n8 ðxÞ þ wk nk ðxފ; ð18Þ ds ¼ @ x2 i i where w4 , w8 , and wk are the weighting coefficients. The motivation of (18) is the observation that contour length is proportional to the number of neighbor pairs crossing the contour. As will be seen, the advantage of this approach is that we can calculate how much region perimeter increases when a pixel is added to the region border by looking at local information only. To determine the values of the weighting coefficients, we do an experiment to find the relation between n4 , n8 , nk , and the theoretical perimeter of circles and rectangles. The results are shown in Fig. 5. The results for circles confirm the linear relationship between the perimeter and n4 , n8 , and nk , namely, 2r  0:79n4  0:56n8  0:18nk , where r is the circle radius. Taking the final length as the average of the three individual estimates, we set w4 ¼ 0:26; w8 ¼ 0:19; wk ¼ 0:06: ð19Þ 1; . . . ; KÞ ¼ K Xn X i¼1 x2 i i þ L~i ðxÞþ o ½w4 n4 ðxÞ þ w8 n8 ðxÞ þ wk nk ðxފ : ð20Þ The smoothing term used in (20) is somewhat similar to that in the Bayesian model [3], [28]. This term is nonzero only near the region boundary. The algorithm needs a partition of the image as starting point. For example, the original watershed segmentation can be used. Then, an exchange of boundary pixels between the regions is performed such that the energy (20) is minimized. Here, we define boundary pixels as ones having at least one neighbor with a different label. Suppose a boundary pixel x is currently assigned to i and it is adjacent to region j . Let n4 ðx; jÞ; n8 ðx; jÞ, and nk ðx; jÞ be the number of pixels with label j, respectively, in the 4, 8, and “knight-move”-connected neighborhood of x. If we change the pixel assignment from label i to label j, n4 ðx; jÞ pairs of 4-connected pixels disappear, while n4 ðx; iÞ new pairs are created. With the same observation for the other two kinds of connectivity, it can be verified that the change of energy (20) is
E~ðx; i ! jÞ ¼ ð j ÿ i Þ þ ðL~j ðxÞ ÿ L~i ðxÞÞ n þ 2 w4 ½n4 ðx; iÞ ÿ n4 ðx; jފ þ w8 ½n8 ðx; iÞ ÿ n8 ðx; jފ ð21Þ o þ wk ½nk ðx; iÞ ÿ nk ðx; jފ : Relating this to (16), note that the last term in the right-hand side of (21) can be regarded as a measure for the curvature (but a inaccurate one). Obviously, reassignment of x should occur if
E~ðx; i ! jÞ < 0. The value of
E~ is called the stability of x [28]. Junction points may have several reassignment possibilities. In this case, the stability of a junction pixel is defined as the minimal value of
E~. For each reassignment, the pixel with lowest stability is selected. Since reassignments can only reduce the energy E~, the convergence of the algorithm is guaranteed. As noted, (21) uses a crude measure of curvature. This measure takes into account only the local information at the pixel considered and, therefore, smoothing occurs only at a small scale. To smooth corner structures at larger scales, the algorithm should run in multiscale mode. The minimization result at a low-resolution level is used as the initial state for the minimization at the next higher resolution level. One advantage of multiscale smoothing is that resulting boundaries are not jagged as often happens in methods based on a single scale [2]. In the latter approach, only global corners are smoothed but local corners remain unaffected, causing jaggedness. Another advantage of the multiscale mode is the avoidance of undesired local minima of the energy function. As we have proven, the original watershed segmentation gives the Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. 336 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 25, NO. 3, MARCH 2003 Fig. 5. Results of the experiment for the determination of w4 , w8 , and wk . (a) The numbers n4 , n8 , and nk are counted at the border of digital disks of radius ranging from 5 to 200 and plotted as function of the radius. The plot shows a linear relationship between the numbers of neighbor pairs and the perimeter of the disks: 2r  0:79n4  0:56n8  0:18nk . (b) The numbers n4 , n8 , and nk are counted at the border of a digital rectangle of size 210  194, the orientation of which is then rotated with respect to the grid. The plots show the ratios (in percentage) of 0:79n4 , 0:56n8 , and 0:18nk , respectively, to the theoretical perimeter of the rectangle as function of the rotation angle. (c) The same as in (b) but for the combination ð0:79n4 þ 0:56n8 Þ=2. (d) The same as in (b) but for the combination of ð0:79n4 þ 0:56n8 þ 0:18nk Þ=3. global minimum of the energy function (13). However, due to adding the contour length as in (15) and (20), the energy function may have many local minima. As the minimization procedure finds a local minimum only, the result depends on the choice of the initial segmentation. The coarse-to-fine strategy prevents the algorithm from getting trapped in an undesired local minimum by providing a good initial segmentation, obtained at lower resolution levels. The algorithm is summarized below: It consists of two stages. The result of Stage 1 is identical to the original watershed segmentation as in [15], [22]. Since Stage 1 increases regions based on a minimum distance criterion, it is also similar to the algorithm in [29] except for a different distance measure. The major contribution of our algorithm, therefore, is Stage 2 where the segmentation result of Stage 1 is smoothed by minimizing energy (20). ALGORITHM STAGE 1 1. Compute K distance transforms Li ðxÞ ( see [15], [22] ). 2. Initialize the regions i from the regional minima Mi . 3. From unassigned pixels on the outer boundary of the regions, select one with minimal value of i þ Li ðxÞ, where i is the label of the adjacent region. Assign the selected pixel to region i . 4. Iterate Step 3 until all pixels are assigned. STAGE 2 1. Start from a low-resolution level by subsampling the label image as well as the K distance transforms. 2. For every boundary pixel compute the stability according to (21). 3. Select the boundary pixel with lowest stability. When this value is negative, perform the reassignment. 4. Recompute the stability of the pixels in the chamfer neighborhood of the reassigned pixel according to (21). 5. Iterate Steps 3 and 4 until no reassignment is possible. 6. Repeat Steps 2-5 for higher resolution levels. Let us give more details on the downscaling process in Stage 2. The input for the minimization at the upper or coarse resolution level l is obtained from the previous level l ÿ 1 by subsampling the label image and the K distance maps by a factor of two in both dimensions. To be specific, a label of a pixel x at level l is selected from the labels of the four corresponding pixels x1 , x2 , x3 , and x4 at the finer resolution l ÿ 1 according the majority criterion. The K topographical distance maps could be recomputed for each level. However, since this recomputation is time consuming, we use a simplified approach in which the topographical distance Li ðxÞ at level l is obtained by averaging of the values of Li ðx1 Þ, Li ðx2 Þ, Li ðx3 Þ, and Li ðx4 Þ at level l ÿ 1. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION 337 Fig. 6. (a) A brain image. (b) The relief computed from morphological gradient. (c) The markers extracted. (d) The result of the original watershed segmentation, shown for comparison. Since the final segmentation is obtained at the finest resolution, the result is always a local minimum of the energy (20). Note also that the algorithm is guaranteed to converge since the energy function always decreases. 6 RESULTS In this section, we show illustrations of the performance of watersnakes in segmentation using the algorithms described. Before the segmentation started, images were smoothed by the opening by reconstruction filter [12]. For the smoothed images, the morphological gradient rm ¼ I ÿ "I was computed, where I and "I, respectively, denote the gray-value dilation and erosion of I with the minimum disk. For color images, we defined the gradient as the maximum of the morphological gradients computed for the three RGB channels. The markers were extracted as low-gradient zones, where the gradient is below 0.5 times of the estimated standard deviation of the gradient over the entire image. Let  0 if x 2 markers ð22Þ hðxÞ ¼ rm ðxÞ otherwise and gðxÞ ¼  0 1 if x 2 markers otherwise: ð23Þ The relief f was obtained by the reconstruction by erosion of g with mask h. The reconstruction algorithm can be found in [30]. Fig. 6a shows an original brain slice image of size 375  371. Fig. 6b shows the relief image obtained. The markers extracted are shown in Fig. 6c. The result of the original immersion-based watershed algorithm [10], [9] is shown in Fig. 6d for comparison. The segmentation results obtained with the energy discretization-based algorithm, presented in Section 5.3, are shown in Fig. 7 for different values of the smoothing coefficient . The minimization was performed at three scales. At the lowest resolution level, the image was reduced by a factor of four in both dimensions. Although the algorithm used a simple measure for curvature, in multiscale mode, the smoothing result is quite satisfactory. More segmentation results are shown in Figs. 8 and 9. To focus the reader attention, in all examples except for the top row of Fig. 8, we show the results for the object of interest only. The result in the top row of Fig. 8 can be compared to that in [2], which aimed to segment the same image. Our segmentations have less regions. Moreover, the regions are smooth and the most prominent regions in the image. As observed in the results, the watersnake algorithm has two advantages over the watershed. The first one is the ability to impose smoothness on the segmentation result. In Figs. 6d and 8a, the original watershed line is rather noisy and jagged. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. 338 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 25, NO. 3, MARCH 2003 Fig. 7. The segmentation result of the watersnake algorithm based on energy discretization with (a) ¼ 10, (b) ¼ 50, (c) ¼ 100, and (d) ¼ 150. Note, in comparison with the original watershed segmentation in Fig. 6d, that the results in figure are smoother, but still identify the main objects. Fig. 8. Segmentation by (a) watershed and (b) watersnake ( ¼ 50). In the bottom row, the result is shown for the object of interest only. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION 339 Fig. 9. Segmentation by (a) watershed and (b) watersnake ( ¼ 50). The results are shown for the object of interest only. These contours are usually not what the user desires. Smooth contours as in the results of Figs. 7 and 8b are important for further investigation of object shape. The second important advantage of the watersnake algorithm is that it prevents the formation of wrong limbs, the unwanted regions attached to the object of interest. These limbs often occur in the watershed segmentation due to a leak at the object boundary. Examples are shown in the bottom row of Fig. 8a and in Fig. 9a. As observed in Figs. 8b and 9b, the watersnake algorithm has significantly shortened the wrong limbs, resulting in much better segmentation compared to the left column. Apparently, smoothness of the resulting boundary is always at the cost of losing adherence to edges at some corners. Examples are some protrusions in Figs. 7c and 7d or the foot of the cat in the top row of Fig. 9b, which have been smoothed over. We still cannot give an optimal scheme for the determination of the coefficient since this issue is subjective and depending on the application. Nevertheless, the proposed method is advantageous over the traditional watershed and postsmoothing of the result. As noted in Section 4, boundaries in our results are not smoothed equally but rather depending on the sharpness of local edges, i.e., abruptness of intensity change. Since important edges are usually sharp, they tolerate the smoothing better than weak edges. This can also be observed in Fig. 7. Pay attention to the indentation on the top of the brightest region in the middle. Despite high curvature, this indentation remains unaffected in Figs. 7a and 7b, while protrusions at other weaker edges have been smoothed over. The intactness of the true limbs of the tree in Fig. 8b is another example. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. 340 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, Fig. 10. (a) An example of a thick geodesic SKIZ. The bold lines illustrate the shortest geodesic paths from x to the seeds A1 and A2 . (b) A partition that has a boundary of zero area. 7 CONCLUSION This paper has established a connection between the wellknown watershed segmentation from mathematical morphology and energy-based segmentation methods. Using the distance-based approach for the definition of the watershed line, we have derived an energy function whose minimization is equivalent to the watershed segmentation. Using this point of view, a priori knowledge on boundary smoothness or shape can now be taken into account. In particular, we obtain smooth watersheds by adding the length of the region boundary to the energy function. It yields a new segmentation method called watersnake that is implemented by a multiscale algorithm. Significant improvements in smoothness of segmentation results, using the new method, have been confirmed by experiments. Apart from the incorporation of a priori knowledge into the watershed segmentation, the results of the paper can also be used for the combination of the watershed lines obtained on different relief functions. The minimization of a linear combination of the energy functions yields a natural trade-off in segmentation between different types of edge indicating functions. In a forthcoming paper, we are exploring this for tracking regions in an image sequence. APPENDIX A PROOF OF THEOREM 1 In this appendix, we prove the existence of at least one watershed segmentation according to Definition 5. When the watershed line is thick, the assignment of the watershed pixels into the catchment basins is ambiguous. This is known as the plateau problem. As noted, several methods [10], [16], [11], [22] use the geodesic distance in plateau to solve this ambiguity. Unlike the euclidean SKIZ [19], the geodesic SKIZ may be thick (see Fig. 10). Therefore, none of the above methods can guarantee an ideally thin watershed line that has a zero area (or equivalently, empty interior). In essence, the proof presented below specifies a way to construct a watershed segmentation that has boundaries of zero area. Our solution to the plateau problem is to assign all pixels in the plateau to one catchment basin unless the constraints of connectedness and regularity of the partition are violated. If it is desired that the watershed line lies in the middle of the plateau, one can also use the geodesic distance to divide the VOL. 25, NO. 3, MARCH 2003 Fig. 11. Illustration of the watershed region in case of three catchment basins. plateau first and then use the proposed approach to divide the remaining geodesic SKIZ that may be thick. We shall call a set A regular if IntðAÞ ¼ A. Such a set does not have barbs although it may contain barbs of other sets inside it. If all regions of a partition are regular and, furthermore, the border of the union set X does not have barbs, then the border of every region does not have barbs either. We use the following proposition that has been proven in [18]. Proposition 3. Let S be an open set. Given K seeds: A1 ; . . . ; AK which are disjoint subsets of S, where every Ai , i 2 f1::Kg, is connected and regular. We assume further that every pixel in S is connected to one of Ai . Then, we can partition S into K disjoint regions R1 ; . . . ; RK such that every region Ri , i 2 f1::Kg is connected and regular, and, furthermore, Ai  Ri . Using this proposition, we now can prove Theorem 1. Proof of Theorem 1. Given x 2 X , we use I ðxÞ to denote the set of indices such that 8i 2 I ðxÞ : i þ Li ðxÞ ¼ min f 1jK j þ Lj ðxÞg: ð24Þ Let ðxÞ be the number of indices in I ðxÞ and ð25Þ S ðnÞ ¼ fx 2 X jðxÞ  ng: SK ð1Þ ðKÞ ¼ X . It can be Note that S ¼ i¼1 CBi and S verified that S ðnÞ is an open set and, furthermore, every pixel in S ðnÞ is connected to one of the regional minima (using a similar proof as in Proposition 2). To illustrate why we need sets S ðnÞ , let us consider the segmentation with three catchment basins, i.e., K ¼ 3. The watershed, being the set of points with equal distances to the two nearest catchment basins, is then the union of the following sets (see Fig. 11): W12 ¼ fx j 1 þ L1 ðxÞ ¼ 2 þ L2 ðxÞ; þ L1 ðxÞ < 3 þ L3 ðxÞg ¼ fx j 2 þ L2 ðxÞ ¼ 3 þ L3 ðxÞ; 2 þ L2 ðxÞ < 1 þ L1 ðxÞg 1 W23 W13 ¼ fx j 1 W123 ¼ fx j 1 þ L1 ðxÞ ¼ 3 þ L3 ðxÞ; 1 þ L1 ðxÞ < 2 þ L2 ðxÞg þ L1 ðxÞ ¼ 2 þ L2 ðxÞ ¼ 3 þ L3 ðxÞg: To complete the segmentation, the regions of the catchment basins need to be extended to fill up the Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. NGUYEN ET AL.: WATERSNAKES: ENERGY-DRIVEN WATERSHED SEGMENTATION watershed region. Note, however, that points in W12 cannot be assigned to the catchment basin CB3 because this would violate (12) in our watershed definition. For that reason, we do not construct a partition of the whole image immediately. Instead, we grow the catchment basins within set S ð2Þ first. By our definition: 341 We consider the two following cases: 1. S ð2Þ ¼ fxjðxÞ  2g ¼ CB1 [ CB2 [ CB3 [ W12 [ W23 [ W13 : In this way, points in W12 will not be assigned to CB3 since W12 and CB3 are disconnected, as will be proven below. The final segmentation is constructed by assigning points in the remaining region W123 to one of the regions obtained from the partition of S ð2Þ . We are now back to the main proof with a general K. In order to prove the existence of a watershed segmentation, we construct a series of partitions of S ðnÞ ðnÞ ðnÞ S ðnÞ ¼ R1 [ . . . [ RK ; n ¼ 1; . . . ; K ð26Þ APPENDIX B PROOF OF THEOREM 2 ð27Þ Proof. In this appendix, we prove the equivalence of the watershed to the minimization of energy (13). From (7), it follows that i þ Li ðxÞ  fðxÞ. Therefore, for any partition 1 ; . . . ; K we have: with ð1Þ Ri ¼ CBi ; i ¼ 1; . . . ; K: 2. There exists an index ‘ 2 I ðzm Þ and ‘ 62 I ðxÞ. This implies ‘ þ L‘ ðzm Þ ¼ min0jK f j þ Lj ðxÞg. Since this equality also holds at the limit z , we have ‘ 2 I ðz Þ. On the other hand, I ðxÞ  I ðz Þ. It follows that I ðxÞ [ f‘g  I ðz Þ and, therefore, ðz Þ, the number of indices in I ðz Þ is larger than that in I ðxÞ: ðz Þ > n. This contradicts to the fact that is contained in S ðnÞ . I ðzm Þ is a subset of I ðxÞ and ðzm Þ < n. This implies zm 2 S ðnÿ1Þ . Since i 62 I ðxÞ ) i 62 I ðzm Þ and, thereðnÿ1Þ ðnÿ1Þ ðnÞ fore, zm 62 Ri . It follows that zm 2 Rj  Rj ðnÞ for some j 6¼ i. At the same time, zm 2 Ri , this ðnÞ ðnÞ contradicts to the fact that Ri \ Rj ¼ ;. u t ðnÞ ðnÞ The regions R1 ; . . . ; RK are obtained by applying ðnÿ1Þ ðnÿ1Þ ; . . . ; RK as seeds. Then, Proposition 3 for S ðnÞ with R1 ðKÞ ðKÞ the final regions fR1 ; . . . ; RK g constitute a disjoint ðKÞ partition of X , where each region Ri is connected and Eð KÞ ¼ 1; . . . ; i¼1 regular. ðKÞ It remains to prove that fRi g is a watershed segmentation. We shall prove that for every n ¼ 1; . . . ; K, 8x 2 ðnÞ Ri : i þ Li ðxÞ ¼ min f 0jK j þ Lj ðxÞg: ð28Þ We employ the principle of mathematical induction. Since ð1Þ Ri ¼ CBi , (28) holds for n ¼ 1 by the definition of the catchment basins. Assuming that (28) holds for n ÿ 1, we shall prove that it also holds for n. ðnÞ Let x 2 Ri . We assume conversely that i þ Li ðxÞ > min f 1jK j þ Lj ðxÞg: That also means i 62 I ðxÞ. As assumed, (28) holds for n ÿ 1, it follows that ðnÿ1Þ ðnÞ ðnÞ . Furthermore, since 8j 6¼ i : Ri \ Rj ¼ ; x 62 Ri ðnÿ1Þ ðnÞ ðnÞ ðnÿ1Þ  Rj , we also have Ri \ Rj ¼ ;. From and Rj ðnÞ ðnÿ1Þ x 2 Ri , it follows that x does not belong to any Rj , j 6¼ i. It follows that ðxÞ > n ÿ 1. ðnÞ At the same time, as assumed, x 2 Ri  S ðnÞ and from the definition of S ðnÞ , we have ðxÞ  n. Combining ðxÞ > n ÿ 1 and ðxÞ  n, we conclude ðxÞ ¼ n. ðnÞ ðnÞ Since Ri is connected, one can find a path in Ri from x to CBi , where ð0Þ ¼ x and ð1Þ 2 CBi . Observe that the sets I ð ð0ÞÞ and I ð ð1ÞÞ are different. Therefore, if we travel along from x to CBi , a change in I ð ðtÞÞ occurs somewhere. Let z be the first point where I ð ðtÞÞ changes. In formulas, z ¼ ðt Þ, where t ¼ supft0 j80  t < t0 : I ð ðtÞÞ  I ðxÞg. By this definition, one can find a sequence zm ¼ ðtm Þ, where tm > t and limm!1 zm ¼ z such that I ðzm Þ 6 IðxÞ. K ZZ X  f i þ Li ðxÞgdx i K ZZ X i¼1 fðxÞdx ¼ i ZZ ð29Þ fðxÞdx: x2X For a watershed segmentation, the function E is minimized since (29) becomes an equality. On the other hand, if (29) is an equality, then ZZ f i þ Li ðxÞ ÿ fðxÞgdx ¼ 0: 8i : i Since i þ Li ðxÞ ÿ fðxÞ  0, it follows that ZZ f i þ Li ðxÞ ÿ fðxÞgdx ¼ 0 U for any set U  i . That means: 8x 2 Int i : i þ Li ðxÞ ¼ fðxÞ and as both Li ðxÞ and fðxÞ are continuous, this equality holds also on Intð i Þ ¼ i : 8x 2 This implies i 1; . . . ; : i K þ Li ðxÞ ¼ fðxÞ: ð30Þ is a watershed segmentation. t u ACKNOWLEDGMENTS The authors would like to thank Professor Arnold Smeulders for his valuable comments on earlier versions of this paper. They would also like to thank the anonymous reviewers for their helpful comments and suggestions that have improved the paper. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply. 342 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] D. Mumford and J. Shah, “Optimal Approximation by Piecewise Smooth Functions and Associated Variational Problems,” Comm. Pure and Applied Math., vol. 42, pp. 577-684, 1989. S.C. Zhu and A. Yuille, “Region Competition: Unifying Snakes, Region Growing, and Bayes/MDL for Multiband Image Segmentation,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 18, no. 9, pp. 884-900, Sept. 1996. Y.G. Leclerc, “Constructing Simple Stable Descriptions for Image Partitioning,” Int’l J. Computer Vision, vol. 3, no. 1, pp. 73-102, 1989. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active Contour Models,” Int’l J. Computer Vision, vol. 1, no. 4, pp. 321-331, 1987. V. Caselles, R. Kimmel, and G. Sapiro, “Geodesic Active Contours,” Int’l J. Computer Vision, vol. 22, no. 1, pp. 61-79, 1997. S.J. Osher and J.A. Sethian, “Fronts Propagating with Curvature Dependent Speed: Algorithms Based on Halmiton-Jacobi Formulations,” J. Computational Physics, vol. 72, pp. 12-49, 1988. H. Tek and B.B. Kimia, “Image Segmentation by ReactionDiffusion Bubbles,” Proc. Int’l Conf. Computer Vision, ICCV ’95, pp. 156-162, 1995. S. Geman and D. Geman, “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 6, no. 6, pp. 721-741, Nov. 1984. S. Beucher and F. Meyer, “The Morphological Approach of Segmentation: The Watershed Transformation,” Mathematical Morphology in Image Processing, E. Dougherty, ed., chapter 12, pp. 43-481, New York: Marcel Dekker, 1992. L. Vincent and P. Soille, “Watershed in Digital Spaces: An Efficient Algorithm Based on Immersion Simulations,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 13, no. 6, pp. 583-598, June 1991. F. Meyer and S. Beucher, “Morphological Segmentation,” J. Visual Comm. and Image Representation, vol. 1, no. 1, pp. 21-46, Sept. 1990. P. Salembier, “Morphological Multiscale Segmentation for Image Coding,” Signal Processing, vol. 38, no. 3, pp. 359-386, Sept. 1994. L. Najman and M. Schmitt, “Watershed for a Continuous Function,” Signal Processing, vol. 38, no. 1, pp. 99-112, July 1994. L. Najman and M. Schmitt, “Geodesic Saliency of Watershed Contours and Hierarchical Segmentation,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 18, no. 12, pp. 1163-1173, Dec. 1996. F. Meyer, “Topographic Distance and Watershed Lines,” Signal Processing, vol. 38, no. 1, pp. 113-125, July 1994. F. Preteux, “On a Distance Function Approach for Gray-Level Mathematical Morphology,” Mathematical Morphology in Image Processing, E. Dougherty, ed., pp. 323-350, New York: Marcel Dekker, 1992. G. Borgefors, “Distance Transforms in Digital Images,” Computer Vision, Graphics, and Image Processing, vol. 34, pp. 344-371, 1986. H.T. Nguyen, M. Worring, and R. van den Boomgaard, “Watersnakes: Energy-Driven Watershed Segmentation,” Technical Report 12, Intelligent Sensory Information Systems Group, Univ. of Amsterdam, 2000, available at http://www.science.uva.nl/ research/reports-isis/2000/ISISreport12.ps. J. Serra, Image Analysis and Mathematical Morphology. New York: Academic Press, 1982. P.L. Lions, Generalized Solutions of Hamilton-Jacobi Equations. Boston: Pitman, 1982. P.L. Lions, E. Rouy, and A. Tourin, “Shape-from-Shading, Viscosity Solutions and Edges,” Numerische Mathematik, vol. 64, pp. 323-353, 1993. J.B.T.M. Roerdink and A. Meijster, “The Watershed Transform: Definitions, Algorithms and Parallelization Strategies,” Fundamenta Informaticae, vol. 41, pp. 187-228, 2000. F. Meyer, “Flooding and Segmentation,” Mathematical Morphology and Its Application to Image and Signal Processing, J. Goutsias, L. Vincent, and D.S. Bloomberg, eds., pp. 189-198, Kluwer, 2000. P. Maragos and M.A. Butt, “Advances in Differential Morphology: Image Segmentation via Eikonal PDE and Curve Evolution and Reconstruction via Constrained Dilation Flow,” Mathematical Morphology and Its Application to Image and Signal Processing, H. Heijmans and J. Roerdink, eds., pp. 167-174, 1998. D.G. Lowe, “Organization of Smooth Image Curves at Multiple Scales,” Proc. Int’l Conf. Computer Vision, pp. 558-567, 1988. VOL. 25, NO. 3, MARCH 2003 [26] M. Pardas and P. Salembier, “Time-Recursive Segmentation of Image Sequences,” Proc. European Assoc. Signal, Speech, and Image Processing, pp. 18-21, 1994. [27] B. Marcotegui and F. Meyer, “Bottom Up Segmentation of Image Sequences for Coding,” Annals of Telecommunications, vol. 52, no. 78, pp. 397-407, 1997. [28] P.B. Chou and C.M. Brown, “The Theory and Practice of Bayesian Image Labeling,” Int’l J. Computer Vision, vol. 4, no. 3, pp. 185-210, 1990. [29] R. Adams and L. Bischof, “Seeded Region Growing,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 16, no. 6, pp. 641-647, June 1994. [30] L. Vincent, “Morphological Gray Scale Reconstruction in Image Analysis: Applications and Efficient Algorithms,” IEEE Trans. Image Processing, vol. 2, pp. 176-201, 1993. Hieu Tat Nguyen received the Eng and MSc degrees in computer technology from the University “Lvivska Polytechnica” of Lviv, Ukraine in 1994. He received the PhD degree in computer science from the University of Amsterdam, the Netherlands, in 2001. He is currently a postdoctoral fellow with the Intelligent Sensory Information Systems group at the University of Amsterdam. His current research interests include image sequence analysis, object tracking, object recognition, active learning, mathematical morphology, and content-based image retrieval. Marcel Worring received the masters degree in computer science (honors) from the Free University, Amsterdam. His PhD thesis was on digital image analysis and was obtained from the University of Amsterdam, The Netherlands, in 1993. He became an assistant professor in 1995 in the Intelligent Sensory Information Systems group at the University of Amsterdam. His current interests are in multimedia information analysis, in particular, document and video analysis. He has been a visiting researcher in the Department of Diagnostic Imaging at Yale University (1992) and at the Visual Computing Lab at the University of California, San Diego (1998). Rein van den Boomgaard graduated from Delft University in 1988 and received the PhD degree from the University of Amsterdam in 1992. He is an assistant professor in the Intelligent Sensory Information Systems (ISIS) group at the University of Amsterdam. Currently, he is working in the field of computer vision and image processing with research interests in color vision, scalespace theory, and mathematical morphology. Since 1999, he has also been an independent consultant for companies needing advice and feasibility studies in the areas of image processing and computer vision. . For more information on this or any other computing topic, please visit our Digital Library at http://computer.org/publications/dlib. Authorized licensed use limited to: IEEE Xplore. Downloaded on December 29, 2008 at 11:51 from IEEE Xplore. Restrictions apply.