The steepest watershed : from graphs to images

The steepest watershed : from graphs to images

Fernand Meyer
Centre de Morphologie Mathématique
Mines-ParisTech, Paris, France
fernand.meyer@mines-paristech.fr
Abstract

The watershed transform is a powerful and popular tool for segmenting objects whose contours appear as crest lines on a gradient image : it associates to a topographic surface a partition into catchment basins, defined as attraction zones of a drop of water falling on the relief and following a line of steepest descent. To each regional minimum corresponds a catchment basin. Points from where several distinct minima may be reached are problematic as it is not clear to which catchment basin they should be assigned. Such points belong to watershed zones, which may be thick. Watershed zones are empty if for each point, there exists a unique steepest path towards a unique minimum. Unfortunately, the classical watershed algorithm accept too many steep trajectories, as they use too small neighborhoods for estimating their steepness. In order to nevertheless produce a unique partition they do arbitrary choices, out of control of the user. Finally, their shortsidedness results in unprecise localisation of the contours. We propose an algorithm without myopia, which considers the total length of a trajectory for estimating its steepness ; more precisely, a lexicographic order relation of infinite depth is defined for comparing non ascending paths and chosing the steepest. For the sake of generality, we consider topographic surfaces defined on node weighted graphs. This allows to easily adapt the algorithms to images defined on any type of grids in any number of dimensions. The graphs are pruned in order to eliminate all downwards trajectories which are not the steepest. An iterative algorithm with simple neighborhood operations performs the pruning and constructs the catchment basins. The algorithm is then adapted to gray tone images. The neighborhood relations of each pixel are determined by the grid structure and are fixed ; the directions of the lowest neighbors of each pixel are encoded as a binary number. Like that, the graph may be recorded as an image. A pair of adaptative erosions and dilations prune the graph and extend the catchment basins. As a result one obtains a precise detection of the catchment basin and a graph of the steepest trajectories. A last iterative algorithm allows to follow these downwards trajectories in order to detect particular structures such as rivers or thalweg lines of the topographic surface.

1 Introduction

The watershed line or divide line of a topographic surface is the boundary separating its catchment basins. A drop of water falling on this surface glides along a line of steepest descent until it is captured by a regional minimum. A catchment basin is the attraction zone of a minimum. Catchment basins may overlap and the overlapping zone is precisely a divide line, as a drop of water falling on it may glide towards several distinct minima. Any gray tone image may be considered as a topographic surface, where the altitude is proportional to the gray-tone. Consider in particular the gradient of an image to segment. In absence of noise or texture, the inside of the objects and of the background appears as minima and the contours appear as crest lines separating the minima. Each object, appearing in the gradient image as the catchment basin of a minimum, is easily extracted by the watershed transform. Each minimum gives birth to a catchment basin ; with two many meaningless minima, over-segmentation occurs. Marker based segmentation regularizes the gradient image ; it consists in flooding the topographic surface in order to keep only one regional minimum for each object of interest [2]. More generally, closing reconstructions or floodings regularize the gradient as they completely fill some catchment basins up to their lowest pass point. As a result these catchment basins are absorbed by neighboring basins, yielding a coarser segmentation. To a series of increasing floodings will be associated a hierarchy [16]. The filling of lakes and subsequent absorption of their catchment basins by neighboring regions may be ordered according some geometric criteria, such as the depth of the lakes [9], [18], [19], their area or their volume [22],[23]. They may also be filled in an interactive mode and be a building block for interactive segmentation [26].

The watershed being a powerful tool for segmentation has been victim of its success : a number of definitions and algorithms have been published, claiming to construct a watershed line or catchment basins, although they obviously are not equivalent. It is out of the scope of the present paper to give an exhaustive bibliography of the concept of watershed. Distinct algorithms published in the literature will produce different results. But even the same algorithm may produce distinct results if one changes the scanning order of the image. Jos. Roerdink gives a review of the most popular watershed definitions and implementations [21]. An historical analysis with many references on how the watershed idea has been developed, triggered both by theoretical considerations and by technological possibilities for its implementation may be found in [17].

We propose an algorithm for the following situation : a gray tone image is given, represented on a grid or on a graph and we desire delineating its catchment basin. Moreover we desire constructing a partition into catchment basins, often at the price of arbitrary divisions of the zones where they overlap. We will compare the performances between algorithms dealing with the same situation. These algorithms do not cover the situation where a topographic surface is constructed by a shortest path algorithms and represents for each pixel its distance to predefined seeds as in [14], [7]. The Voronoï tessellation associated to this distance also constitutes a partition ; furthermore, if one applies the shortest path algorithm by Dijkstra, one may assign to each node one ancestor, through which its distance to the nearest neighbor has been constructed. The algorithm produces a minimum spanning forest.

Further we do not consider the algorithms which construct a thin watershed line separating catchment basins. In a hierarchy one goes from a fine to a coarse partition by merging adjacent regions. This operation is immediate if one deals with partitions : one assigns to all regions to be merged the same label. It is however more problematic if the contour is materialized between the regions and paradoxical situations may be met if one does not carefully chose the graph representing the images [5]. This is an additional reason why to prefer watershed zones without boundaries between regions ; furthermore representing contours wastes space in the image and makes it impossible to segment adjacent small structures.

As a conclusion, we consider here the watershed transforms associating to a topographic surface a partition into catchment basins, defined either by water streaming down or flowing up. The first definition defines the catchment basins as the domain of attraction of a drop of water falling on the surface and gliding downwards to reach a minimum. The second floods the domain from sources placed at the minima ; lakes containing distinct minima meet along the divide line. The progression of the rain or of the flood follows lines of steepest descent. Whether a drop of water glides down a surface or a flood progresses upwards finally has no importance, the result will be the same if they follow identical trajectories. Among all non increasing trajectories between a node and a minimum, the algorithm has to chose the steepest ; if several equivalent trajectories connect a node with distinct minima, this node belongs to two overlapping catchment basins ; if the target is a partition, we have to divide the overlapping zone between the two adjacent basins. It is not clear and often out of control of the user, how this division is made by the algorithms published in the literature. The resulting partition may dramatically vary if one applies one type of division rules or another. The same algorithm may also give different results if one simply changes the scanning order of the image. This makes the positioning of the contours is neither stable nor precise.

Ideally, if each pixel had only one downwards trajectory towards a unique minimum, catchment basins would not overlap. In the present work we do not reach complete unicity but get as close as possible to this ideal situation. We propose a method which selects the steepest trajectories possible. We propose an algorithm without myopia, which considers the total length of a trajectory for estimating its steepness. It turns out that for natural images these trajectories often are unique.

For the sake of generality we develop the algorithm for node weighted graphs. The result may then easily be transposed to images defined on arbitrary grids and any number of dimensions.

The outline of the paper is the following. We first consider topographic surfaces defined on node weighted graphs. The graphs are pruned in order to eliminate all downwards trajectories which are not the steepest. An iterative algorithm with simple neighborhood operations performs the pruning and constructs the catchment basins. The algorithm is then adapted to gray tone images. The graph structure itself is encoded as an image thanks to the fixed neighborhood structure of grids. A pair of adaptive erosions and dilations prune the graph and extend the catchment basins. As a result one obtains a precise detection of the catchment basins and a graph of the steepest trajectories. A last iterative algorithm allows to follow selected downwards trajectories in order to detect particular structures such as rivers or thalweg lines of the topographic surface.

2 The watershed on weighted graphs

2.1 Weighted graphs

A non oriented graph is a collection of vertices or nodes and of edges ; an edge being a pair of vertices [1],[8]. Each node is weighted with a real number .

The subgraph spanning a subset is defined as , where are the edges linking two nodes of

The partial graph associated to a subset of edges has the same nodes as but has less edges ; it is defined by

A gray tone image may be considered as a graph, in which each pixel becomes a node, with a weight representing its gray tone. Neighboring pixels are linked by an edge.

A subgraph of a node weighted graph is a flat zone, if any two nodes of are connected by a path where all nodes have the same altitude.

A subgraph of a graph is a regional minimum if is a flat zone and all neighboring nodes have a higher altitude.

2.2 Drainage graphs and their catchment basins

We associate to an oriented graph , called drainage graph, encoding all possible directions a drop of water may follow when falling on a given node outside a regional minimum. If is one of the lowest neighboring nodes of in , an arc is created in from towards If and both do not belong to regional minima and have the same weights, then two arcs are created, one from to and one from to

A drainage path between two nodes and is a sequence of nodes such that two successive nodes and are linked by an arc.

Lemma 1

From each node outside a regional minimum starts a drainage path to a regional minimum

Proof. Any node outside a regional minimum has a lower neighboring node, if it does not belong to a plateau. Otherwise it belongs to a plateau, containing somewhere a node with a lower neighboring node, as the plateau is not a regional minimum. The plateau being connected, there exists a path of constant altitude in the plateau between and . Following this path, it is possible, starting at node to reach the lower node This shows that for each node there exists an oriented path to a lower neighboring node. Taking this new node as starting node, a still lower node may be reached. The process may be repeated until a regional minimum is reached.   

Thanks to this lemma, it is possible to define the catchment basins of the minima.

The catchment basin of a minimum is the set of all nodes from which starts an oriented path towards in

The watershed zones are the nodes belonging to more than one catchment basin.

If a node has several lowest neighbors, each arc towards one of them is the first arc of a drainage path. Each lowest neighbor in turn may have several second lowest neighbors, multiplying the number of drainage paths with the same origin ; if two such paths link a node with two distinct minima, then belongs to a watershed zone.

The restricted catchment basins are the nodes which belong to only one catchment basin : it is the difference between the catchment basin of a minimum and the union of catchment basins of all other minima. From a node in a restricted catchment basin, there exists a unique non ascending path towards a unique regional minimum.

2.3 Steepest drainage paths

For some nodes in a drainage graph there exist several oriented paths towards regional minima. If a node has several lower neighbors, each arc towards one of them may be the first arc of such a path. Each lowest neighbor in turn may have several lowest neighbors, multiplying the choices of drainage paths. This is the reason why thick watershed zones appear, as several oriented paths, starting from the same origin node may reach distinct regional minima.

The number of drainage paths with the same origin is reduced if we consider not only the first neighbors of each node, but larger neighborhoods. Suppose that a node has two lowest neighbors and ; each arc and is the first arc of an oriented path towards a minimum. Suppose now that the lowest neighbor of is lower than the lowest neighbor of indicating that the path passing through is steeper than the path passing through : both have identical weights on the first two nodes, but the third node differenciates them. More generally, two paths may have identical weights on the first nodes and distinct nodes on the node A lexicographic order relation will be defined for comparing the steepness of non ascending paths. steepness of a non ascending path will be based on the lexicographic order. The weights of the nodes along an oriented drainage path, starting from a given node and following the path downwards, are by construction a series of never increasing weights. A lexicographic order may be defined for comparing paths: a path is steeper than a path if or if there exists an index such that for we have and

The number of distinct paths with the same origin towards distinct minima sharing exactly the same list of non increasing weights is extremely low and often reduced to one if one considers natural images. This shows that if one considers only the paths with maximal steepness in a drainage graph, the size of the watershed zones is strongly reduced. They will be empty, if there exists only one path with maximal steepness linking each node with a regional minimum.

2.4 Pruning drainage graphs

If a path of a drainage graph has a maximal steepness, then any sub-path obtained by skipping the first nodes is still a path of maximal steepness. If another path of origin were steeper, then concatenating and would produce a path steeper than This lemma has a consequence: if is not the first arc of a steepest path with origin then no other steepest path will pass through this edge. As a matter of fact, if were a steepest pass passing through the sub-path of starting at would be a steepest path, which is not the case.

The previous analysis shows that if we cut all arcs of a drainage graph which are not the highest edge of a steepest path, we obtain a partial graph with exactly the same steepest paths as the original one. For cutting the highest edge of each non steepest path, the naive approach would compare the paths with the same origin two by two, follow each one downwards until two edges are found with distinct weights : the path leading to the highest edge would be the pass with a lower steepness and its initial edge could be pruned.

A more clever pruning algorithm relies on the erosion on the drainage graph, which lets the isolated regional minima nodes unchanged and assigns to each other node the weight of its lowest neighbors to which it is linked by an arc. Consider two paths and with identical weights for the first nodes and verifying indicating that is steeper than Eroding the graph once assigns to the node the weight and to the nodes and the weights With successive erosions, the weights of the nodes along the paths and file past the two first nodes. The erosion assigns to the node the weight and to the nodes and respectively the weights and which are different. Hence, after erosions, the graph is not a drainage graph anymore, as the edge links with which is not one of its lowest neighbor. This edge should thus be pruned in order to obtain a partial graph which still is a drainage graph. This cutting achieves the intended pruning : the first edge of a non steepest path is cut.

Figure 1: Construction of the catchment basins of a node weighted graph

This leads to the following algorithm for obtaining the steepest drainage graph, starting with an edge weighted graph :
* construct a drainage graph by linking by an arc each node of outside the regional minima with each of its lowest neighbors. Detect, label the regional minima and replace their inside arcs by edges. Figure 1_B represents the drainage graph of fig.1_A.
* Repeat until stability : a) erode the graph and cut all arcs linking a node to another which is not one of its lowest neighbors ; b) if is an oriented arc from to and if holds a label whereas has no label, then is assigned the label of , the arc between and is replaced by an edge and all arcs with origin are suppressed. If points to several nodes with distinct labels, one of them is chosen and the same procedure applied.
Applied on fig.1_B the erosion and pruning produces fig.1_C ; the label propagation produces fig.1_D . The next iteration produces fig. 1_E and 1_F. After the last erosion, pruning and final label propagation, all nodes are labeled and the steepest drainage graph obtained in fig. 1_G. The final partition is superimposed on the initial graph in fig.1_H.

3 The watershed on images

3.1 Representing an oriented graph for images.

In order to transpose to images the preceding algorithm initially defined for graphs, one has to find a representation of the drainage graph itself. The nodes are simply the pixels of the image to which the watershed algorithm is applied. The nodes hold three types of valuations and will be represented on three images. First, a gray tone image holding the initial distribution of gray tones and its evolution as the algorithm proceeds. The second image holds the labels of the minima and of the catchment basins as they expand. The last image is more original as it has to encode the drainage graph itself.

The grids on which images are represented have a regular structure, where each node has the same number of neighbors, in identical positions. Numbering the neighbors according to their direction allows to represent the neighborhood relation of each pixel with a binary number, where each bit encodes for one direction. The n-th bit is set to 1 if and only if there exists an arrow between the central point and its n-th neighbor ; such oriented arrows having as origin a node are called out-arrows of . Figure 2 shows the numbering of the directions in a hexagonal raster and the corresponding bit planes (on the right). The bottom image gives an example of encoding a particular neighborhood configuration : 2 + 1 + 32 = 35. This representation has been introduced by F. Maisonneuve in his seminal work on watersheds [11].

Figure 2: The encoding of the directions of the neighbors in the hexagonal raster, and the weights of the corresponding bits in the binary representation of the arrows. An example with three arrows with weights 2, 1 and 32 is represented by the binary number 100011, i.e. the decimal number 35.

The algorithm creates and updates 3 images: 1) the gray tone image itself, 2) an image of labels representing the regional minima and the catchment basins in construction, 3) the encoding of the out-arrows representing the arcs of the drainage graph.  Fig. 3 presents the three images. Fig. 3.1 represents the gray tone image. Fig. 3.2 represents the image of out-arrows encoding the drainage graph. Fig. 3.3 represents a label at the central position, represented as a colored dot.

Figure 3: The three images used for constructing the catchment basins : a gray tone image, an image of arrows and an image of labels. 

Fig. 4.1 presents a gray tone image. Fig. 4.2 combines the three images: the gray tone value for each node, a colored dot representing the label of the minima and the initial out-arrows encoding the arcs.

Figure 4: A gray tone image, the arrows representing its drainage graph and the labels of the regional minima. 

3.2 An adaptive erosion and dilation, guided by the arrows

The successive prunings of arrows are easily implemented with two neighborhood transformations, the first being an adaptive erosion for propagating the gray tone values and pruning the arrows, the second being an adaptive dilation for propagating the labels of the regional minima as they are extended.

3.2.1 An adaptive erosion

We define an adaptive erosion and combined pruning on the arrowed image. It uses and updates both the gray tone image as the arrows image. If a pixel is without arrows, it is left unchanged. Otherwise, it is replaced by the lowest of its arrowed neighbors and only the arrows towards these neighbors are kept. In figure 5 a pixel has two lowest neighbors, but only one of them is arrowed. Hence only the arrow towards this node is kept.

Figure 5: An example of adaptive erosion and pruning. 

In figure 6 the two lowest neighbors are not arrowed : they are discarded and only the lowest arrowed neighbor is taken into consideration : its value is assigned to the central pixel and only the arrow towards it is kept.

Figure 6: Adaptive erosion and pruning

When no arrow is present, the central pixel is left unchanged as in figure 7.

Figure 7: In the absence of any arrow, the central pixel is left unchanged.

3.3 An adaptive dilation

The adaptive dilation propagates the labels of the regional minima. It modifies both the image of labels and of arrows. It is guided by the drainage graph.

Recall that the labels are represented by strictly positive values, the pixels without labels having the value in the labeled image. Furthermore, the algorithm cares for the fact that labeled pixels have no arrows. This is true at initialization, when the regional minima get their labels. It is also true during the expansion of the catchment basins, as a pixel loses its arrows as soon it gets a label.

The label propagation is done by an adaptive dilation of the label images guided by the arrows image. One considers the pixels without labels (a pixel with a label has no out-arrows) ; such a pixel gets the highest label of its arrowed neighbors. In the absence of arrowed and labeled neighbors, this value is 0. The operation is thus an adaptive dilation. Furthermore, every time a pixel gets a label, its arrows are suppressed, as it now belongs to a catchment basin. This produces an additional pruning of the drainage graph. Below we illustrate the combination of the adaptive erosion and dilation in a number of situations

Remark 2

1) Neighboring nodes of the same catchment basins are not linked by an edge ; they are identified by the label they hold.
2) In case where a pixel without label has 2 or more out-arrows towards distinct labeled pixels, only one of them has to be chosen in order to get a partition. Chosing the highest of them constitutes an adaptive dilation. . As a matter of fact, we may chose arbitrarily one of them. This is the only place where a choice takes place in the algorithm. It divides the catchment zones and produces a partition. Such situations are nevertheless rare, as we propagate the labels along the trajectories whose steepness is estimated by taking into account their total length. The necessity of a choice appears only in the case where two trajectories have exactly the same distribution of node weights, from top to bottom.

3.4 Combination of the adaptive erosion and dilation : illustration

The following figure shows how the adaptive erosion and dilation are used in sequence. The erosion assigns to the central pixel the value of its lowest arrowed neighbor ; only the arrow towards this neighbor is kept. And when this arrow points towards a labeled neighbor, this label is propagated to the central pixel by the adaptive dilation, and its arrows are suppressed. Figures 8,9,10,11 present how, in different neighborhood configurations, the combined adaptative erosion and dilation erode the gray tone, prune the arrows and propate the labels.

Figure 8: Adaptative erosion, pruning and guided dilation of the labels 
Figure 9: Adaptative erosion, pruning and guided dilation of the labels 
Figure 10: In the case where two or more distinct labels are present in the direction of arrows, the highest of them, or one, chosen arbitrarily, is propagated
Figure 11: Case where a labeled pixel is present in the neighborhood of the central pixel, but as it is not arrowed, it is not propagated to the central pixel.

3.5 The complete watershed algorithm

Figure 12 shows that the algorithm can be initialized with arrows in all directions. After the first combined adaptive erosion and dilation, only the arrows towards the lowest neighbors are kept and the labels propagated accordingly. After the initialization phase, the combined erosion of gray tones, pruning of arrows and dilation of labels is iteratively applied until the labels cover the total domain. Convergence is attained after 4 iterations for figure 12.

Figure 12: Initial gray tone distribution, labeling of the regional minima and arrowing in all directions, followed by 4 iterations of a combined adaptive erosion and dilation.

By recording the arrows existing for each pixel, just before it gets its label, one obtains the final and steepest drainage graph, where all prunings have been done as illustrated in fig. 13. The trajectories of a drop of water falling on the surface are extremely selective and narrow. We will use this selectivity in the last part of the paper for following lines of steepest descent and thalweg lines on a topographic surface.

Figure 13: Final catchment basins, and recording of the last arrow distribution of each node before it gets its label. 

3.6 Successive steps of the pruning

During the successive adaptive erosions, no choice is ever made. A choice may appear necessary when two when two distinct arrowed and labeled pixels are present in the neighborhood of a pixel. The adaptive dilation choses the highest of them. Other rules may be introduced, as for instance a random choice. The occurence of such choices is rare in natural images, as they only happen when two distinct drainage paths linking a node with two distinct minima has exactly the same distribution of weights.

Figure 14: Successive prunings of the arrows. 

3.7 Complexity

After the initial detection and labeling of the minima, the algorithm needs a number of iterations equal to the largest distance between a pixel in a catchment basin and its regional minimum. Each iteration consists in the combination of the adaptive erosion and dilation.

3.8 The problem of the plateaus

The plateaus pose a particular problem to all watershed algorithms which only consider local neighborhoods. As a matter of fact, a drop of water falling on a plateau has no clear direction for reaching the nearest regional minimum. The classical solution consists in constructing a geodesic distance transform to the lowest neighbors of the plateau and to follow the steepest descent line on this function. Fig. 15.1.1 shows a topographic surface containing several plateaus with value The geodesic distance within each plateau to its lower boundary is illustrated in fig. 15.1.2. This produces a topographic surface on which we may compute the drainage graph, as illustrated in fig. 15.2.2. By comparison, the steepest drainage graph produced by our algorithm is more selective and has less arrows It is illustrated in fig. 15.2.1. Furthermore no special treatment is required for dealing with the plateaus : they are treated as any other part of the topographic surface.

Remark 3

Some watershed algorithms use arrows for the construction of the watershed. F.Maisonneuve is the first to assign arrows to all pixels with lower neighbors, and then iteratively completing the arrowing inside plateaus [11]. F. Lemonnier, in a hardware implementation of the watershed, constructs separately the arrows of the drainage graph and those of the geodesic distance within the plateaus, before regrouping both and completing the watershed construction [10].

Figure 15: Top: On the left a topographic surface and on the right the geodesic distance to its lower boundary on each plateau. Bottom: On the left the arrows of the steepest drainage graph and on the right the arrows associated to the distance function. 

4 Illustration on a real image

The figures 16 and 17 illustrate the method on a real image. They contain in the top row a gray tone image and its gradient. The bottom row contains on the left the labeled regional minima and on the right the associated catchment basin. In fig. 16.2.2 they hold the same labels as the minima they contain. Fig. 17.2.2 is a mosaic image where each catchment basin takes the mean gray tone of the initial image in this region. Fig. 18 represents the final drainage graph and shows the image of the arrows encoded in false color.

Figure 16: Initial image, gradient, labeled regional minima and labeled catchment basins. 
Figure 17: Initial image, gradient, labeled regional minima and as final result the catchment basins containing each the mean gray tone of the initial image. 
Figure 18: Central part of the arrow image of the steepest flooding graph represented in false color. 

4.1 Contrast with the flooding algorithms : the watershed of a digital elevation map

Fig. 19.1.1 presents a digital elevation map of an existing landscape : each gray tone represents an altitude. Due to sensing errors, there are spurious regional minima in the topographic surface. As the rivers leave the image in the direction of the see, the only minima which make sense are those on the boundary of the image. A marker image is produced, equal to the relief on the boundary of the image and to elsewhere. The highest flooding of the relief under this mask has all its minima touching the boundary (see fig. 19.1.2 ); all other spurious minima due to sensing errors have been suppressed. lexicographic watershed produces the partition in catchment basins illustrated by fig. 19.1.3. Its catchments are the real catchment basins of the topographic surface, containing each a river, appearing as a thalweg line.

Figure 19: The gray tone image represents a digital elevation map. The central image shows all regional minima touching the boundary of the image. The right image presents the catchment basins.

As announced earlier, the labels are propagated along lines of steepest descent ; at each iteration they progress one step further. When they reach the top of a drainage path, they stop. This may be clearly seen in fig. 20 showing the extension of the catchment basins after 20, 40, 60, 80, 100 and 120 iterations of the elementary step of the algorithm (adaptive erosion of gray tones, dilation of labels and pruning of arrows). The construction of the catchment basins touching the lower boundary is achieved after a low number of iterations as they are small. They stop growing, independently of the adjacent catchment basins. The adjacent catchment basin, associated to the largest river reaches them after many more iterations.

This is in absolute contrast with the flooding algorithms which construct the catchment basins as attraction zones of the regional minima. For such algorithms, a catchment basin stops growing only when it arrives in contact with another catchment basin. The most well known such algorithms implement the topographic distance, as recalled below.

4.1.1 The flooding algorithms

Steepest descent lines are followed both by a drop of water gliding downwards or by an increasing flood, which invades a topographic surface. The great majority of watershed algorithms mimic the flooding of a topographic surface. The regional minima are the sources, the level of the flood is uniform : as it increases, the new lakes are created and old ones expand. The algorithm cares that two lakes containing distinct minima do not merge. As a matter of fact, the floodings progresses according lines of steepest descent of the surface furthermore, these lines can be modeled as the geodesics of the topographic distance function [20],[15]. The catchment basins represent the Voronoï tessellation of the regional minima for the topographic distance to the minima, whose altitude has been set to 0. The following classical algorithm is derived from this definition. It uses a hierarchical queue which correctly organizes the flooding in the presence of plateaus [13]: a hierarchical queue is a series of first in first out queues, with a priority order between the queues. The pixels are put in the queue corresponding to their altitude ; pixels with a lower altitude having a higher priority than pixels with a higher altitude. Furthermore, pixels of the same queue are treated on a first in first out basis : like that the pixels of plateaus of uniform altitude are treated in the order of increasing distances to the boundary of the plateaus. The algorithm constructing the zones of influence of the minima for the topographic distance is the following:

Label the nodes of the minima and introduce them in a hierarchical queue HQ each with a priority equal to their weight. 

As long as the HQ is not empty, extract the node with the highest priority from the HQ:

For each unlabeled neighboring (on the flooding graph) node of

*

* put in the queue with priority

Remark 4

The hierarchical queue structure does most of the job for a proper flooding. Lower pixels are flooded before higher ones, due to the hierarchy between the queues ; pixels near the boundary of a plateau are treated before more inwards pixels, due to the first in first out administration of each queue. This management of flooding saves the day for the short-sighted watershed algorithms, which accept too many downstream trajectories, resulting in thick watershed zones. The hierarchical queue structure (or other similar structures) helps making a proper division of these thick watershed zones. The formulation of the algorithm is deceptive, as it gives an elegant mathematical definition of catchment basins, but the fact that they largely overlap is hidden ; the illustrations look convincing, however, the credit for this quality is to be given more to the data structure for implementing them as to the way they are defined. The algorithms with the highest myopia are those based on the flooding ultrametric distance [12],[4].

Figure 20: Extension of the catchment basins after 20, 40, 60, 80, 100 and 120 iterations of the combined adaptive erosion of gray tones, pruning and dilation of labels.
Figure 21: Construction of the catchment basins with an algorithm based on uniform flooding

The grows of catchment basins for both algorithms may now be compared. Fig. 20 shows the extension of the catchment basins after 20, 40, 60, 80, 100 and 120 iterations of the combined adaptive erosion of gray tones, pruning and dilation of labels. One remarks a striking difference with the classical algorithm for constructing the watershed which is based on the simulation of a flooding, where a flood starting from the regional minima grows with a uniform altitude and progressively invades the topographic surface. Fig. 21 precisely shows the progressive construction of the catchment basin based on the flooding distance; the unflooded part appears in dark blue, the flooded parts appears as colored labels. The successive levels of flooding represented are 75, 105, 135, 170, 205 and 240 (gray-tones on a scale ). The flooding algorithm is a greedy shortest distance algorithm based on the topographic distance. Each catchment basin stops growing everywhere it meets another catchment basin. With the steepest path algorithm on the contrary, the catchment basins are dilated at each iteration by a dilation of size one and they stop growing when they have reached their full extension, have they reached another catchment basin or not. As a consequence, if we suppress the label of a minimum at initialization, the catchment basin of this minimum will remain empty.

5 Downstream trajectories of a drop of water

5.1 The algorithm for downstream propagation

The successive prunings of a drainage graph leaves a minimal graph containing the steepest downstream trajectories for each node. Given a few labeled starting points, it is possible to follow the downstream trajectories, simply by following the downstream arrows. Fig. 22 illustrates the operator which is used. If a pixel has a label (fig. 22.1) and an arrow towards a neighboring pixel (fig. 22.2) then the label is propagated from to (fig. 22.3). If has already a label, the maximum of both labels is chosen. This expansion is repeated until stability.

Figure 22: Downstream following of a drop of water: a labeled node is expanded in the directions of its arrows to its neighboring nodes.

In fig. 23.1 two starting points are chosen, and assigned distinct labels, represented respectively by a red and a green dot. After downwards propagation following the arrows, both trajectories appear as labeled nodes in fig. 23.2, each of them reaching a regional minimum.

Figure 23: Two starting nodes are chosen and assigned a label. On the right figure, the downwards trajectories are illustrated. 

This method is applied on the same DEM image in fig. 24.1 after the construction of the steepest drainage graph. A number of positions are chosen by hand in in fig. 24.2, where a drop of water will start its downwards trajectory, highlighted in the image on the right. The drop of water duly glides downwards until it meets and follows a river and ultimately reaches the boundary of the image in the direction of the sea, as illustrated in fig. 24.3. Each river keeps the label of its highest source and keeps this label unless it meets another label which is higher. The connected components colored in red and touching the boundary of the image represent the regional minima of the image.

Figure 24: After construction of the arrows of the DEM, a number of starting points are chosen and labeled (central image) and the downstream trajectories of a drop of water falling on these points highlighted in the right image.

5.2 Application to the detection of fibers, cracks, thalweg lines.

Fig. 25 contains two spirals, intricated one in another, ending with a dark regional minimum. The background is brighter than the spiral. A red dot and a green dot have been put on the other extremity of both spirals. The regional minimum has been marked by a blue dot. The spirals are plateaus with a uniform gray-tone, with one lower boundary in the regional minimum. The steepest descent trajectories taking their origin in the red and green tots are represented in fig.26. They start from the red and green dots and follow the stripes until they reach the minimum.

Figure 25: Left : Two spirals on a bright background ending in a dark minimumRight : A red and and a green dot mark the extremities of the spiral. A blue dot the minimum. 

The spiral stripes in fig. 26 are plateaus of constant altitude, without internal structure for centering the trajectories, which appear at some places as large as the stripes in the right image. For a better centering of the trajectories, we have transformed the binary image into a gray tone image by constructing a geodesic distance to the lower borders of the plateaus. The arrows of the steepest drainage graph within a small square crossing a stripe are shown in fig. 27. The stripes are well centered and the trajectory along their thalweg well delineated. On this new relief, the detected trajectories are much thinner (fig.28)

5.2.1 Comparison with shortest path algorithms

The classical algorithms for following and highlighting thin and elongated dark structures (let us call it fiber) in a noisy bright background rely on shortest path algorithms (for instance cracks in a porous medium, hairs, glass fibers, vessels in 2 or 3 dimensions etc.). As the gray tone along the fiber is darker than in the background, the integral of gray tone along the fiber is smaller than along a path with the same length lying in the brighter background. This integral may be seen as a weighted distance transform. The method consists in computing the weighed distance along the fiber, first in one direction, then in the opposite direction ; the sum of both distance transforms is then minimal along the fiber. The method has been first proposed for detecting cracks in porous material [24],[25]. The method works well if the cracks or fibers or more or less rectilinear.

However, if the fiber is tortuous, like the spirals above, then a shortest path algorithm between the regional minimum and the red dot (resp. green) will find shortcuts and will not follow the spiral. The classical solution consists in a stepwise progression : one progresses along the fiber on a short distance, so as to remain on it and initiates a new progression at the arrival point. Like that, little jump after little jump, one progresses. The length of the jump depends upon the contrast between fiber and background [6].

On the contrary, using the drainage graph does nor present this weakness : the trajectories are delineated and based on a long range lexicographic distance which makes it possible to follow them completely from top to down, until a regional minimum is met. Contrarily to the shortest path algorithm, the algorithm has not to be used stepwise, does not depend upon the contrast between foreground and background. It is also able to follow several fibers at the same time. However, it is certainly more sensitive to noise or missing data than the shortest path algorithms.

Figure 26: The spiral stripes have a constant altitude, that is they are plateaus, without internal structure for centering the trajectories, which appear at some places as large as the stripes. 
Figure 27: The arrows of the steepest drainage graph after constructing the geodesic distance of each node of a plateau to the lower border of the plateau. 
Figure 28: Replacing each strip by a distance transform to its lowest boundary permits to center the downwards trajectories inside the stripes.

6 Discussion

The present work is the first implementation and validation of a more ambitious work (under publication) where we consider node or edge weighted graphs and show that they are in fact equivalent from the point of view of the watershed construction. We show how to complete the missing weights, such that the weight of a node is the smallest weight of its adjacent edges, i.e. an erosion from edges to nodes, and the weight of an edge the highest weight of its extremities, i.e. a dilation from nodes to edges. Erosion and dilation being adjunct from each other, the associated opening is particularly interesting : the edges invariant by this opening are all edges which are the lowest edge of one of their extremities. This indicates that a flood coming from this higher extremity may pass through this edge.

We next define two erosions, assigning to a node weight(resp. edge) the minimal weight of its adjacent nodes (resp. edges). After eroding edges and nodes, a new graph is created which is not necessarily invariant by the opening We prune it by suppressing all edges which are not invariant by Repeating this operation, erosion/pruning until stability, produces a minimal graph, where only the steepest paths remain.

In the present work also we use successive erosions of the drainage graph and prunings so as to get a drainage graph containing only the steepest paths. The difference lies in the implementation of the algorithm. Here we have no possibility to assign weights to edges and we have to use oriented arcs for representing the drainage graph. As the neighborhood structure is fixed, the repartition of arrows having a pixel as extremity may be encoded as a binary number of length equal to the number of neighbors.

The present algorithm could also be applied to an edge weighted graph after transforming it into a node weighted graph, obtained by assigning to each node a weight equal to the smallest weight of its adjacent edges.

7 Conclusion

We have presented a particularly simple parallel algorithm for the construction of the watershed of gray tone images. As it is defined for node weighted graphs, it may be adapted for any type of grid and any number of dimensions. It extracts from a node weighted graph the smallest drainage graph possible, without any arbitrary choices. Arbitrary choices only remain in the propagation of labels, in the case where a given node is linked with two distinct minima by 2 drainage paths with exactly the same weights in the same order.

The algorithm is based on simple parallel neighborhood operators, which are particularly suitable for an implementation in special hardware or in graphics processors.

Many variants of this algorithm can be thought of:

  • It may be optimized in different ways. For instance as soon a pixel has only one arrow, its label can arrive only from this direction. In terms of graphs, the corresponding edge could be contracted. In terms of images, one could connect such pixels with a ”union-find” type of algorithm.

  • The algorithm could also be used as a pre-processing step, with only a small number of iterations performed, in order to significantly prune the graph. The arrows may then be inverted (if arrows in an image arrows after inversion of the arrows). Like that any classical algorithm based on hierarchical queues would consider propagating the labels only in the direction of the inverted arrows.

  • It may also be used as a preprocessing step in some applications where a high speed is required : apply a fixed number of steps of pruning, and then suppress all possibilities of choice, by keeping only one arrow for each node. The algorithm concludes by extracting and labeling the connected components. (A.Bieniek, in [3], keeps an arrow only between each node and one of its lowest neighbors and then labels the resulting graph.)

The lexicographic watershed selects a minimal set of downstream directions. Furthermore, the selectivity of the pruning allows to construct long range downstream trajectories, even in the presence of high tortuosity.

—————————————————

References

  • [1] C. Berge. Graphs. Amsterdam: North Holland, 1985.
  • [2] S. Beucher and C. Lantuéjoul. Use of watersheds in contour detection. In Proc. Int. Workshop Image Processing, Real-Time Edge and Motion Detection/Estimation, 1979.
  • [3] Moga A. Bieniek, A. An efficient watershed algorithm based on connected components. Pattern Recognition, 33(6):907–916, 2000.
  • [4] J. Cousty, G. Bertrand, L. Najman, and M. Couprie. Watershed cuts: Minimum spanning forests and the drop of water principle. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 31(8):1362 –1374, aug. 2009.
  • [5] Jean Cousty, Michel Couprie, Laurent Najman, and Gilles Bertrand. Grayscale watersheds on perfect fusion graphs. In Ralf Reulke, Ulrich Eckardt, Boris Flach, Uwe Knauer, and Konrad Polthier, editors, Combinatorial Image Analysis, volume 4040 of Lecture Notes in Computer Science, pages 60–73. Springer Berlin / Heidelberg, 2006.
  • [6] Thomas Deschamps and Laurent D. Cohen. Fast extraction of minimal paths in 3d images and applications to virtual endoscopy. Medical Image Analysis, 5(4):281 – 299, 2001. ¡ce:title¿Soft Tissue Deformation¡/ce:title¿.
  • [7] A.X. Falcao, J. Stolfi, and R. de Alencar Lotufo. The image foresting transform: theory, algorithms, and applications. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 26(1):19 – 29, jan 2004.
  • [8] M. Gondran and M. Minoux. Graphes et Algorithmes. Eyrolles, 1995.
  • [9] M. Grimaud. New measure of contrast : dynamics. Image Algebra and Morphological Processing III, San Diego CA, Proc. SPIE, 1992.
  • [10] F. Lemonnier. Architecture Electronique Dédiée aux Algorithmes Rapides de Segmentation Basés sur la Morphologie Mathématique. PhD thesis, E.N.S. des Mines de Paris, 1996.
  • [11] F. Maisonneuve. Sur le partage des eaux. Technical Report-Centre of Mathematical Morphology-Mines-ParisTech, 1982.
  • [12] B. Marcotegui and S. Beucher. Fast implementation of waterfalls based on graphs. ISMM05 : Mathematical Morphology and its applications to Signal Processing, pages 177–186, 2005.
  • [13] F. Meyer. Un algorithme optimal de ligne de partage des eaux. In Proceedings Congrès AFCET, Lyon-Villeurbanne, pages 847–857, 1991.
  • [14] F. Meyer. Color image segmentation. In IEE Fourth International Conference on Image Processing and its Applications, pages 303–306, April 1992.
  • [15] F. Meyer. Topographic distance and watershed lines. Signal Processing, pages 113–125, 1994.
  • [16] Fernand Meyer. Flooding and segmentation. In John Goutsias, Luc Vincent, and Dan S. Bloomberg, editors, Mathematical Morphology and its Applications to Image and Signal Processing, volume 18 of Computational Imaging and Vision, pages 189–198. 2002.
  • [17] Fernand Meyer. The watershed concept and its use in segmentation : a brief history. (arXiv:1202.0216v1), 2012.
  • [18] L. Najman. Morphologie Mathématique: de la segmentation d’images à l’analyse multivoque. PhD thesis, Université Paris-Dauphine, 1994.
  • [19] L. Najman. Geodesic saliency of watershed edges and hierarchical segmentation. IEEE Trans. Pattern Anal. Machine Intell, 16(3):175–182, 1996.
  • [20] Laurent Najman and Michel Schmitt. Watershed of a continuous function. Signal Processing, 38(1):99 – 112, 1994. Mathematical Morphology and its Applications to Signal Processing.
  • [21] Jos B. T. M. Roerdink and Arnold Meijster. The watershed transform: Definitions, algorithms and parallelization strategies. Fundamenta Informaticae, 41:187–228, 2001.
  • [22] C. Vachier. Extraction de Caractéristiques, Segmentation d’Image et Morphologie Mathématique. PhD thesis, E.N.S. des Mines de Paris, 1995.
  • [23] C. Vachier and F. Meyer. Extinction values: A new measurement of persistence. In I. Pitas, editor, 1995 IEEE Workshop on Nonlinera Signal and Image Processing, pages 254–257, 1995.
  • [24] Jeulin D. Vincent, L. Minimal paths and crack propagation simulations. Acta Stereologica, 8(2 II):487–494, 1989.
  • [25] Luc Vincent. Minimal path algorithms for the robust detection of linear features in gray images. In Proceedings of the fourth international symposium on Mathematical morphology and its applications to image and signal processing, ISMM ’98, pages 331–338, Norwell, MA, USA, 1998. Kluwer Academic Publishers.
  • [26] F. Zanoguera, B. Marcotegui, and F. Meyer. A segmentation pyramid for the interactive segmentation of 3-d images and video sequences. In John Goutsias, Luc Vincent, and Dan S. Bloomberg, editors, Mathematical Morphology and its Applications to Image and Signal Processing, volume 18 of Computational Imaging and Vision, pages 223–232. Springer US, 2002.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
240210
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description