Serge Beucher

Fast implementation of large erosions and dilations in Mamba

Introduction

Large structuring elements are used in many morphological transformations. They are needed when computing size distributions, in filtering operators (and especially in alternate sequential filters), for regularised gradients and large top-hat transforms, for residual operators (skeletons by openings, critical balls, ultimate openings and so on). Moreover, in some applications, images are large and even very large. Remind that Mamba is able to process images containing up to 16.10 6 pixels, value which corresponds for instance to a 4096x4096 pixels image. So, it is not unusual, when dealing with these images, to need very large erosions, dilations, openings or closings. Performing an ultimate opening on such an image may require the use of opening sizes up to 2000! Most of the time, large size operators use iterations of an elementary one. These implementations are mainly hardware implementations based on GPU, parallel processors, pipe-line processors [START_REF] Clienti | Architectures flot de données dédiées au traitement d'images par Morphologie Mathématique -Doctorat Morphologie mathématique[END_REF], [START_REF] Brambor | Algorithmes de la morphologie mathématique pour les architectures orientées flux -Doctorat Morphologie Mathématique[END_REF]. These implementations are often very specific, they need a big amount of hardware ressources (memory, calculation power). Moreover, they are not very flexible. For example, they do not allow to use easily more complex structuring elements as octogons on the square grid or dodecagons on the hexagonal one. But it is a matter of fact that these shapes become compulsory when large sizes are at stake. Square shapes in particular are too coarse (and ugly!) to be used efficiently.

Software implementations are even fewer. Their main drawback is that they most often need, to be effective, a specific representation of the image to be processed (for instance, a specific coding of the intercepts for a binary image). These approaches are very clumsy, and they work only in very simple contexts. Moreover, all these implementations are, most of the time, unable or at least inefficient to cope with edge effects which inevitably occur. Some hardware implementations define a thick edge around the image window. This, indeed, avoids edge biases but needs a great amount of memory and multiply by a factor 4 to 8 the computation time. Other implementations do not take care of these errors and accept them, arguing that they are not so important. This way of doing is obviously the worst solution and should not be agreed.

General description

The solution proposed in this document is based on a fast decomposition of operations by doublets of points and on the use of these operations to generate various structuring elements: segments, hexagons, squares, dodecagons and octogons. We also paid attention to propose solutions where edge effects have been addressed and also where the two definition modes (euclidean or geodesic) for these operators are possible.

The efficiency of the implementation described below depends on the availability of fast operators for image shifting and for dilations and erosions by a doublet of points. Fast operator means that the computation time is supposed to be independant of the size of the operation. For instance, the computation time for an image shifting of size 3 should be the same (or approximately the same) as for a shifting of size 300... Fortunateley, these operators are available in Mamba. They are named shift, supFarNeighbor and infFarNeighbor. See the Mamba Reference Manual [START_REF] Beucher | Mamba Image Library Python Reference[END_REF] for a more detailed description of these transforms. These operators are written in C and wrapped into Mamba. The only (slight) restriction in their use is that the shift or doublets of points directions must be defined among the main directions of the grid in use.

Operations with large segments

The classical way of performing a dilation by a segment of size n in a given direction consists in iterating n size-1 dilations. So, if your segment is large, the transformation will take a rather long time. However, it is very easy to get the same result faster by using dilations by doublets of points. This algorithm seems to have been used for the first time in [START_REF] Sternberg | Computer Architectures Specialized for Mathematical Morphology -Workshop on Algorithmically Specialized Computer Organizations[END_REF] (although we are not sure of it, more information about this is welcome...) and in a hardware architecture implementation.

Let us illustrate this with an example. Suppose that we want to dilate a function or a set by an horizontal segment of size 21. Instead of applying 21 successive size 1 dilations, let us perform the following sequence of dilations by a doublet of points: a dilation of size 6, followed on the previous result by a dilation of size 8, followed by a dilation of size 4, then 2 and finally 1. As illustrated in Figure 1, we obtain finally a dilation by a segment of size 21 but with only 5 steps instead of 21. We have potentially multiply by 4 the computation speed. How to decompose the initial size n into a suitable sequence as shown before? In fact, when n is equal to 2 i -1 (i being any positive integer value), the decomposition is straigthforward: the successive powers of 2, 2 k , with 0 < k < i are used. We can easily see that these dilations produce at the end a dilation by a segment of size n, provided that we apply these successive dilations in decreasing order (largest one first, size 1 dilation ending necessarily the sequence). When the size is not equal to but greater than 2 i -1, the decomposition is not complicated: we start by substracting 2 i_ 1 from n, 2 i -1 being the largest value lower than n. The value n-2 i 1 corresponds to the size of the first dilation by a doublet of points. Note that n-2 i 1 is always lower than or equal to 2 i -1.

Then, we decompose the number 2 i -1 as explained above. If the order of operations is respected (and especially if the first one is the dilation of size n-2 i 1 ), the final transformation is a dilation by a segment of size n. If you start from a set made of a single point, at the end of these successive operations, you will obtain a segment of size n without holes or missing points. Note also that no edge effects appear: if the initial point is at a distance less then n from the edge, the final result is the intersection of the segment of size n with the field of analysis (Figure 1b).

Compared to the classical implementation, this one is dramatically faster especially when the size n of the linear structuring element is large. If 2 i-1 n2 i , the speed increase is equal to n/i, that is n / int(log 2 (n)). In fact, this speed gain is lower because the supFarNeighbor operation in Mamba is a little bit slower than supNeighbor, dilation by an elementary segment. However, for large sizes, this gain is impressive: a dilation by a size 300 segment is about 30 times faster! Erosions by large segments can be performed in the same way, thanks to the availability of inFarNeighbor in Mamba. Note also that these transformations can be equally performed in euclidean or geodesic modes (the edge configuration can be defined in the operators arguments list).

These two operators are implemented in largeLinearErode and largeLinearDilate in the Mamba erodilLarge.py module (see annex).

Operations with large squares

Performing rapidly erosions and dilations with large squares is straightforward. Indeed, a square belongs to the Steiner class of polygons, which means that erosions and dilations with this structuring element can be decomposed into successive operations with segments. Each segment corresponds to a linear part of its boundary (or more precisely its half-boundary). Therefore, a square dilation can be obtained by four successive dilations with segments in the directions 0, π/2, π and 3π/2 (Figure 2). No edge effect is likely to happen for pixels close to the edge because points sent outside the image window by one segment dilation (these points are lost) have no effect on the following dilations which do not need to reinsert them inside the window.

A similar implementation is used for erosions. Theses two operators are named largeSquareErode and largeSquareDilate in erodilLarge.py (see annex).

Operations with large hexagons

Hexagons belong also to the Steiner class. An hexagonal dilation can be performed (theoretically) by three successive linear segmentations in the directions 0, π /3 and 2π /3.

Unfortunately, in practice, this approach does not work properly. The reason? Edge effects... These edge effects have already been encountered when defining the basic hexagonal dilation (size 1 hexagon), see [START_REF] Beucher | Algorithmic description of erosions and dilations in Mamba -To be[END_REF]. It is due to the fact that, for pixels near the edge of the image, a dilation in a given direction may propagate a point outside the image window. However, if this dilation is not the last one, the next dilation using another direction should reintroduce inside the image a point resulting of the dilation of the previous pixel, that is not possible because this previous pixel is unfortunately lost (see Figure 3). In order to cope with this problem, the solution, for the elementary hexagon, was to perform a supremum of dilations by segments in the six directions of the hexagonal grid. However, this solution is no longer possible with size n (n > 1) hexagons because it would correspond to dilations by « snowflakes » or asterisks. Thus, we are faced to a dilemma: on the one hand, the Steiner decomposition is compulsory as it is the only way to generate dilations by large hexagons with segments and, on the other hand, we know that this approach leads to unacceptable edge effects. To solve this dilemma, a mixed solution is proposed here, based on suprema of partial dilations obtained by composition of dilations by segments. Indeed, when performing an hexagonal dilation by concatenating linear dilations, different orders are possible: we can use direction 1 first, followed by direction 3 and direction 5. Or we can use another order. Six different orders are available, each one leading to the same result, except if the pixel to be dilated is near an edge. In this case, some orders may produce incorrect results if edge effects appear as illustrated at Figure 3a above.

Therefore, if we perform two dilations, using different directions orders, we may expect that the supremum of these two dilations will produce a result where edge effects have been removed. For instance, if we use successive directions 4, 6 and 2 in the above example, the union of both results produces the awaited hexagonal dilation. In fact, the second dilation alone gives a right result, however using both dilations produces a satisfactory result, the point being near the left or right edges (Figure 3b). However, we are not sure that this combination will give a correct result whatever the position of the pixel in the image. In fact, in this example, it is not the case. If the pixel to be dilated is near the top left corner of the image, the result is wrong (Figure 3c).

In fact, this approach does not work properly. Whatever the combinations which are used, it is always possible to exhibit counter-examples of pixels which are placed in such a position that edge effects occur.

Figure 4: Result of two successive dilations in directions 1 then 3 (left) and in two directions 3 then 1 (right). The first one produces an important edge effect.

Therefore, to be sure to obtain a correct result for all pixels in the image, it is necessary, not only to perform the supremum of partial dilations but also to combine the directions of the linear dilations in such a way that possible edge effects occuring with some combinations of dilations will be compensated by the others. To insure this, the following combination has been used (Figure 4): we perform a dilation by a segment in direction 1 followed by a dilation in direction 3. Then, we perform on the initial image a dilation by a segment in direction 3 followed by a dilation in direction 1. The supremum of these two dilations is computed. Then, the linear dilation of this supremum in direction 5 is performed. Using two dilations in directions 1 and 3 by simply inverting their order allows to remove possible edge effects which could occur when either direction 1 or direction 3 dilation goes over the edge.

Figure 5: First step of the algorithm (left), use of directions 1, 3 and 5. Second step (middle), use of directions 4, 6 and 2. Note that edge effects do not occur because directions 4 and 6 have been used twice. Supremum of the two previous results (right).

It may happen that both directions produce edge effects (see Figure 5 and examples above). This is why the entire procedure is applied again on the initial image using directions 6, 4 and 2 this time. Finally, the supremum of these two operations is computed to obtain the final result without edge effects.

Nevertheless, this procedure assumes that, for at least one direction of linear dilation, the result of this dilation is always contained inside the the image window (Figure 6). If it is not the case, this procedure cannot avoid edge effects leading to a wrong result. This configuration occurs when the hexagon size n is larger then half the window dimensions (width, height). The only solution to avoid this flaw consists in splitting the operation into successive dilations of size i in the following way: It may seem odd to consider structuring elements larger than the size of the image. However, this situation may happen, when dealing with residual transforms for instance.

• Perform a dilation of size i = min(width, height)/2 if n > i, of size n if not • Compute n = (n -i) • Repeat until n0
The same procedure applies, mutatis mutandis, to large hexagonal erosions.

Both operations are implemented in the erodilLarge.py module. They are named respectively largeHexagonalErode and largeHexagonalDilate. They work on binary or greytone images. The edge of the image can be defined in the operations arguments list.

Operations with large octogons

Large octogons are more interesting than squares on the square grid because they are more isotropic. We already know that an octogonal operation (dilation or erosion) can be achieved by concatenating a square one followed (or preceded) by an operation with a conjugate square or diamond [START_REF] Beucher | Algorithmic description of erosions and dilations in Mamba -To be[END_REF]. In order to obtain an isotropic polygon of size n, the respective sizes of the square and diamond operations are given by:

n 1 =  2 12
×n for the square size

n 2 =n-n 1 = n 12
for the diamond size that is approximately:

n 1 ≃0.41421 n n 2 =n-n 1
If the size of images is less than 10000 pixels in height or width, this approximation is sufficient.

We already know how to perform rapidly transformations with large squares. Regarding conjugate squares (or diamonds), these structuring elements are also Steiner polygons. However, although they could be obtained by linear dilations in directions multiples of π/4, this way of doing arises two problems. The first problem comes from edge effects which have already been described with hexagons. They have the same origin. The second problem is due to the fact that the linear dilations used to generate conjugate squares are defined on the conjugate grid. This leads to dilations with holes as some points are never concerned (Figure 7). So, to avoid edge effects, an approach based on suprema of dilations by doublets of points is applied to obtain dilations by diamonds. To get a size n diamond, the size is decomposed in the same way as it was done for linear structuring elements. For each value i of the decomposition, a dilation by a "cross" structuring element (Figure 8) is performed.

This structuring element is obtained by the union of dilations by doublets of points in the four horizontal and vertical directions. Dilations are concatenated for all the values of the size decomposition (Figure 8b).

The structuring element generated by this procedure is not a true conjugate square as some holes may appear inside it as shown at Figure 8. Moreover, these holes may sometimes be larger than one pixel wide, especially near the edges. Therefore, this procedure cannot be used to obtain unbiased dilations or erosions by conjugate squares. But there is an interesting characteristic of this transformation which allows to use it to produce octogonal operations: there is no point missing on the boundary of the structuring element. This can be proved easily when the size n of the structuring element is equal to 2 i -1. This size can be decomposed into the successive decreasing values: Each point on the boundary is then generated by the dilation by a size 1 "cross" structuring element.

2 i-1 , 2 i-2 , ..., 2 , 1
The size 1 cross center point is itself generated by a dilation by a size 2 cross and so on (see Figure 9). It is then possible to define a path of horizontal or vertical dilations, all of them inside the image window, from each boundary point towards the initial central point.

More precisely, if the initial central point is at coordinates (0, 0) and the boundary point at coordinates (x, y) with xy=n=∑ k=0 i -1 2 k , it is always possible whatever x and y to write:

x= ∑ k ∈I 1 2 k and y= ∑ k ∈I 2 2 k
with I 1 ≠I 2 and I 1 ∪I 2 ={0,i -1}

When n is not equal to but larger than 2 i -1, we can perform a first dilation of size n 1 = n -2 i + 1 followed by the decomposition of the value 2 i -1 described above. But as n 1 ≤2 i -1 , there is always an overlap of the boundaries produced by the size 2 i -1 dilations of the cross of size n 1 (Figure 10). For any size of the dilation, the boundary of the diamond structuring element is complete. So, although it is not possible to obtain a true conjugate square by these means, it is sufficient to get dilations by octogons. Indeed, we saw that the respective sizes n 1 for the square and n 2 for the conjugate square to get an octogonal transformation of size n are given by:

n 1 =  2 12 ×n and n 2 =n-n 1 = n 12
So, if we start from an initial point and if we perform first a dilation by a "sparse" diamond of size n 2 as described above, we generate boundary points at a maximum distance n 2 from this center point (indeed, we generate many more points...). Now, if each of these points together with the original center point are dilated by a square of size n 1 , all these squares cover the interior of the previous hollow conjugate square because n 1 n 2 , leading at the end to the generation of an octogon of size n (Figure 11). This implementation, both for dilations and erosions, has been realised in largeOctogonalDilate and largeOctogonalErode operators. We may wonder why this approach has not been used with hexagonal structuring elements (that is performing iterations of dilations by a structuring element made of the summits and the center of an hexagon). We let the reader verify that this approach is not very convenient as it produces important edge effects.

Operations with large dodecagons

We saw previously that performing fast operations with large hexagons was rather tricky if edge effects must be avoided.These edge effects are mainly due to the use of directions non parallel to the edges. They do not appear with squares because only horizontal and vertical directions are used for the dilations or erosions and because no outside point is needed to generate inside ones (each point can be generated by linear dilation from another inside point). We saw also that these dilations by horizontal and vertical doublets of points allow to build hollow conjugate squares which can be used for octogonal operations. Dodecagonal dilations and erosions being performed by concatenating hexagonal operations and operations by conjugate hexagons, let us describe how faster operations with these latter structuring elements can be achieved.

We explained previously how complex operations with hexagons can be to avoid edge effects. We are faced to the same difficulties with conjugate hexagons. Therefore, the same approach as the one used for conjugate squares has been implemented for conjugate hexagons. This approach generates also sparse or hollow conjugate hexagons. However, these structuring elements are sufficient to perform unbiased dodecagonal erosions and dilations.

Let us describe the algorithm for dilations (performing erosions is straightforward). The first idea consists in using successive dilations by the following structuring element made of the center and the summits of a conjugate hexagon of size i (Figure 12):

Figure 12: Initial structuring element used to compute sparse conjugate hexagons (left). This structuring element is obtained by combining shifts (in green) and supFarNeighbor operators (in red) on the hexagonal grid with supFarNeighbor transforms (in blue) on the square grid. Note that the size of this structuring element is 2).

Each size i is obtained by the size decomposition already used above. This structuring element is built by a combination of shifts and of dilations by doublets of points as illustrated at Figure 12. However, edge effects are likely to happen: in directions 1 or 4 first, where the shift operation may push a point outside the image window if the center point is close to the vertical edges. To avoid this bias, one could perform the shift operation twice, in the right and left directions, before computing the corresponding dilation by a doublet. However, this is rather complex and a simple trick allows to cope with this problem. The trick consists in performing a dilation by a doublet of points in the vertical directions 1 or 5 on the square grid! Edge effects may also appear in other directions as illustrated at Figure 13b. In this case, they come from the proximity of the center point to the horizontal edges. If we consider two successive dilations, the fact that a point has been pushed outside the image window in the first dilation prevents the generation of a necessary point in the following operation (Figure 13c). In order to eliminate this problem, the previous structuring element can be amended by adding two points which correspond to the horizontal summits of the hexagon inscribed in the conjugate hexagon (Figure 13d). So, this point will contribute to the generation of the point which was missing in the previous implementation (Figure 13e). These two points are added by a dilation by doublets of points in directions 3 and 5 (hexagonal grid). If i is the size of the conjugate hexagon (size expressed in number of steps to get it from the elementary conjugate hexagon), the size j of the inscribed hexagon is given by j = 3 i / 2. note that, if i is an odd value, j is given by j=int(3i / 2). In this case, the added point is always inside the conjugate hexagon. This procedure builds a hollow conjugate hexagon without edge effects thanks to the horizontal and vertical dilations by doublets. Moreover, as for the conjugate square, no point is missing on the boundary (the proof is similar to the one used with conjugate squares). Finally, the respective sizes n 1 and n 2 needed for the dilations by an hexagon and a conjugate hexagon to get a size n dodecagon are given by:

n 1 =32-3n≃0.4642×n n 2 = n-n 1 / 2
Now, to cover entirely the conjugate hexagon by dilating its summits and its center by an hexagon in order to fill the holes which remain in it, the minimal size n' of the hexagonal dilation must be n' = n 2 . The actual size of the hexagonal dilation is n 1 =3 n 2 . Therefore, n1 > n2 which insures the covering and the generation of an unbiased dodecagon of size n (Figure 13f). The computation of erosions and dilations by large dodecagons is realised in erodilLarge.py by largeDodecagonalErode and largeDodecagonalDilate.

Performances

The performances of these operators can be coarsely estimated by comparing the number of elementary operations needed to achieve the transformation with the classical approach and with these new implementations.

The transformations with linear structuring elements need the same number of elementary image operators per loop. As this number of loops is lower in the new algorithm than in the classical one, the increase of the computation speed is noticeable from size 2 transform (in fact, it is not exactly the case as the computation time of supFarNeighbor and infFarNeighbor operators is higher than it is with supNeighbor and infNeighbor transforms).

Large square operations are always faster than classical ones because the latter ones do not use the Steiner decomposition of the square contrary to the large ones, so that the number of elementary operations in the first implementation is twice the number of such operators in the second one.

If we compare the complexities of the hexagonal classical transformations and of the large hexagonal ones, we can assume that the number of elementary operations to achieve a classical operation of size n is approximately equal to (7 n + 1). The large hexagonal operators are more complex and need a large number of elementary operators per loop. However, this number of loops is given by int(log2(n)+1) when the size of transformation is n. The number of elementary operators is therefore equal to 10 int(log 2 (n) + 1) + 14.

The following table indicates the approximate number of operators for hexagonal transforms in the range of sizes [0, 10] and for a size equal to 100. We see that, from size 7, the computation speed for the new implementation overtakes the classical one. Finally, the curves above (Figure 14) illustrate the speed comparisons between the large hexagonal erosion and the classical hexagonal erosion. We see a dramatic enhancement of the performances obtained by this implementation so that, for size larger than 20, the computation time is almost independant of the size.

Conclusions

We have shown that, from a simple idea (splitting a linear transformation into successive operations by doublets of points), a correct implementation needs an undoubted effort to be efficient and unbiased. But the gain in speed is worth this effort. This complexity explains surely why this approach has been seldom used especially for hardware implementations. It is a matter of fact that a transformation letting edge effects appear is not acceptable.

The implementation which has been achieved in the Mamba library could still be enhanced by coding some operators directly in C instead of Python. However, the present realisation is a compromise between the ease of use and the never ending race to better performances.

Annex:

Mamba erodilLarge.py module script

""" This module provides a set of functions performing erosions and dilations with large structuring elements. They are built with special shift operators written in C, together with special 'infFarNeighbor' and 'supFarNeighbor' functions.

"

Figure 1 :

 1 Figure 1: (a) successive dilations of an initial point (top line) by doublets of points of respective sizes 6, 8, 4, 2 and 1. In red, points added at each step. (b) Verification that no edge effect may occur.

Figure 2 :

 2 Figure 2: Successive dilations of a point by segments in directions 0, π/2, π and 3π/2 to obtain a square (Steiner polygon).

Figure 3 :

 3 (a) Dilation of a point near the edge by three successive dilations by segments in directions 1, 3 and 5. (b) Same operation with segments in successive directions 4, 6 and 2. (c) The union of these two results do not prevent edge effects when the initial point is in the upper left corner.

Figure 6 :

 6 Figure6: When the size of the structuring element is higher than half the size of the window (here, size is 12 whereas window height is 11), edge effects are unavoidable.

Figure 7 :

 7 Figure 7: Dilations by segments defined on the conjugate grid. The result is a diamond also defined on the conjugate grid (hollow diamond).

Figure 8 :

 8 Figure 8: "Cross" structuring element of size 4 (left) and result of successive dilations of a point by this structuring element (sizes 4, 2 and 1) producing a sparse conjugate square of size 7 (right).

Figure 9 :

 9 Figure 9: The green pixel on the diamond boundary at coordinates (5, 2) is generated by a size 1 dilation of the red pixel, itself obtained by a size 2 dilation of the blue one. The blue pixel has been generated by a size 4 dilation of the initial yellow center point. These points are always propagated inside the image window (refer to the possible position of an edge indicated by an horizontal line in the figure).

Figure 10 :

 10 Figure 10: Boundary overlapping: the boundaries of the size 7 diamonds generated by the red and blue pixels overlap, which insures that no boundary point of the final size 10 diamond is missing.

Figure 11 :

 11 Figure 11: Size 10 octogon obtained by a size 6 sparse conjugate square dilation followed by a size 4 square dilation. We see that the dilations of the center point and of the boundaries of the conjugate square are sufficient to get an octogon without holes.

Figure 13 :

 13 Figure 13: (a) Concatenation of dilations by the structuring element described above (figure 12). (b) Edge effects, 2 points are missing. (c) Missed points should have been generated by points which fell outside the image window at the previous step. (d) Modified structuring element, 2 horizontal points are added. (e) Edge effects have been corrected: the previously missing points are generated by the two added points (red points near the edge). (f) Dodecagons obtained by an hexagonal dilation of these sparse conjugate hexagons.

Figure 14 :

 14 Figure 14: Respective computation times for classical hexagonal erosions (in blue) and for large hexagonal ones (in red). The top chart compares binary operations and the bottom one, 8-bits erosions. Times are given in milliseconds. The image size is 768x576 pixels. These tests have been made with a Core 2 Duo 8400 Intel processor, 3 GHz clock, 4 Go RAM memory running on Linux Fedora 13, 64 bits, SSE2 instructions set enabled.

  , imOut, 6, i, grid=mamba.HEXAGONAL, edge=edge) mamba.supFarNeighbor(prov1, imOut, 5, j, grid=mamba.SQUARE, edge=edge) mamba.shift(prov1, prov2, 5, i, val, grid=mamba.HEXAGONAL) mamba.supFarNeighbor(prov2, imOut, 1, i, grid=mamba.HEXAGONAL, edge=edge) mamba.supFarNeighbor(prov2, imOut, 3, i, grid=mamba.HEXAGONAL, edge=edge) j = 3*i/2 mamba.supFarNeighbor(prov1, imOut, 2, j, grid=mamba.HEXAGONAL, edge=edge) mamba.supFarNeighbor(prov1, imOut, 5, j, grid=mamba.HEXAGONAL, edge=edge) def largeDodecagonalErode(imIn, imOut, size, edge=mamba.FILLED):""" Erosion by large dodecacagons (hexagonal grid). Basically, it is the same operation as the previous one where classical erosions have been replaced by erosions by large structuring elements, and where a "partial" erosion by a conjugate hexagon is used. Dilation by large dodecacagons (hexagonal grid). Basically, it is the same operation as the previous one where classical dilations have been replaced by dilations by large structuring elements, and where a "partial" dilation by a conjugate hexagon is used. Erosion by a large diamond (losange) on square grid. This diamond is not completely filled. It is for internal use only.""" prov = mamba.imageMb(imIn) mamba.copy(imIn, imOut) for i in _sizeSplit(size): mamba.copy(imOut, prov) mamba.infFarNeighbor(prov, imOut, 1, i, grid=mamba.SQUARE, edge=edge) mamba.infFarNeighbor(prov, imOut, 3, i, grid=mamba.SQUARE, edge=edge) mamba.infFarNeighbor(prov, imOut, 5, i, grid=mamba.SQUARE, edge=edge) mamba.infFarNeighbor(prov, imOut, 7, i, grid=mamba.SQUARE, edge=edge) def largeOctogonalErode(imIn, imOut, size, edge=mamba.FILLED): """ Erosion by a large octogon (square grid). This operation uses erosions by large squares and large diamonds previously defined.

	an enhanced shift operator).
	""" prov1 = mamba.imageMb(imIn) # Large dodecagons
	prov2 = mamba.imageMb(imIn) def _sparseConjugateHexagonErode(imIn, imOut, size, edge=mamba.FILLED):
	mamba.copy(imIn, imOut) sizemax = min(imIn.getSize())/2 """
	for i in _sizeSplit(size): # if size larger than sizemax, the operation must be iterated to prevent edge effects. Erosion by a conjugate hexagon. The structuring element used by this operation
	mamba.supFarNeighbor(imOut, imOut, dir, i, grid=grid, edge=edge) n = size is not complete. Some holes appear inside the structuring element. Therefore, this
	mamba.copy(imIn, imOut) operation should not be used to obtain true conjugate hexagons dilations (for
	# Large hexagons while n > 0: internal use only).
	"" def largeHexagonalErode(imIn, imOut, size, edge=mamba.FILLED): s = min(n, sizemax) """
	# Contributor: Serge BEUCHER, August 1, 2010 """ largeLinearDilate(imOut, prov1, 6, s, grid=mamba.HEXAGONAL, edge=edge)
	Erosion by large hexagons using erosions by large segments and the Steiner largeLinearDilate(prov1, prov1, 4, s, grid=mamba.HEXAGONAL, edge=edge) prov1 = mamba.imageMb(imIn)
	import mamba decomposition property of the hexagon. largeLinearDilate(imOut, prov2, 4, s, grid=mamba.HEXAGONAL, edge=edge) prov2 = mamba.imageMb(imIn)
	import mambaComposed as mC def _sizeSplit(size): """ This internal function splits the size of the structuring element into a list of successive and decreasing sizes (except the first one). Successive erosions or dilations by double points produce an erosion or dilation by a segment of length 'size'. """ sizeList=[] incr=1 while size>incr: sizeList.append(incr) size=size-incr incr=2*incr sizeList.append(size) sizeList.reverse() return sizeList # Elementary operators for large structuring elements def largeLinearErode(imIn, imOut, dir, size, grid=mamba.DEFAULT_GRID, edge=mamba.FILLED): """ Erosion by a large segment in direction dir in a reduced number of iterations. Uses the erosions by doublets of points (supposed to be faster, thanks to an enhanced shift operator). """ mamba.copy(imIn, imOut) for i in _sizeSplit(size): mamba.infFarNeighbor(imOut, imOut, dir, i, grid=grid, edge=edge) def largeLinearDilate(imIn, imOut, dir, size, grid=mamba.DEFAULT_GRID, edge=mamba.EMPTY): """ Edge effects are corrected by an erosion with transposed decompositions followed by inf operations (see documentation for further details). This operator is quite complex to avoid edge effects. """ prov1 = mamba.imageMb(imIn) prov2 = mamba.imageMb(imIn) sizemax = min(imIn.getSize())/2 # if size larger than sizemax, the operation must be iterated to prevent edge effects. n = size mamba.copy(imIn, imOut) while n > 0: s = min(n, sizemax) largeLinearErode(imOut, prov1, 6, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(prov1, prov1, 4, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(imOut, prov2, 4, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(prov2, prov2, 6, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(prov1, prov2, prov1, "inf") largeLinearErode(prov1, prov2, 2, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(imOut, prov1, 1, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(prov1, prov1, 3, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(imOut, imOut, 3, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearErode(imOut, imOut, 1, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(prov1, imOut, prov1, "inf") largeLinearErode(prov1, imOut, 5, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(imOut, prov2, imOut, "inf") n = n -s def largeHexagonalDilate(imIn, imOut, size, edge=mamba.EMPTY): """ Dilation by large hexagons using dilations by large segments and the Steiner decomposition property of the hexagon. Edge effects are corrected by dilations with transposed decompositions followed by sup operators. largeLinearDilate(prov2, prov2, 6, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(prov1, prov2, prov1, "sup") largeLinearDilate(prov1, prov2, 2, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearDilate(imOut, prov1, 1, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearDilate(prov1, prov1, 3, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearDilate(imOut, imOut, 3, s, grid=mamba.HEXAGONAL, edge=edge) largeLinearDilate(imOut, imOut, 1, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(prov1, imOut, prov1, "sup") largeLinearDilate(prov1, imOut, 5, s, grid=mamba.HEXAGONAL, edge=edge) mamba.logic(imOut, prov2, imOut, "sup") n = n -s # Large squares def largeSquareErode(imIn, imOut, size, edge=mamba.FILLED): """ Erosion by large squares using erosions by large segments and the Steiner decomposition property of the square. No edge effects are likely to happen with a square structuring element. """ largeLinearErode(imIn, imOut, 1, size, grid=mamba.SQUARE, edge=edge) largeLinearErode(imOut, imOut, 3, size, grid=mamba.SQUARE, edge=edge) largeLinearErode(imOut, imOut, 5, size, grid=mamba.SQUARE, edge=edge) largeLinearErode(imOut, imOut, 7, size, grid=mamba.SQUARE, edge=edge) def largeSquareDilate(imIn, imOut, size, edge=mamba.EMPTY): """ Dilation by large squares using dilations by large segments and the Steiner decomposition property of the square. No edge effects are likely to happen with a square structuring element. """ largeLinearDilate(imIn, imOut, 1, size, grid=mamba.SQUARE, edge=edge) largeLinearDilate(imOut, imOut, 3, size, grid=mamba.SQUARE, edge=edge) mamba.copy(imIn, imOut) val = mamba.computeMaxRange(imIn)[1]*int(edge==mamba.FILLED) for i in _sizeSplit(size): mamba.copy(imOut, prov1) j = 2*i mamba.infFarNeighbor(prov1, imOut, 1, j, grid=mamba.SQUARE, edge=edge) mamba.shift(prov1, prov2, 2, i, val, grid=mamba.HEXAGONAL) mamba.infFarNeighbor(prov2, imOut, 4, i, grid=mamba.HEXAGONAL, edge=edge) mamba.infFarNeighbor(prov2, imOut, 6, i, grid=mamba.HEXAGONAL, edge=edge) mamba.infFarNeighbor(prov1, imOut, 5, j, grid=mamba.SQUARE, edge=edge) mamba.shift(prov1, prov2, 5, i, val, grid=mamba.HEXAGONAL) mamba.infFarNeighbor(prov2, imOut, 1, i, grid=mamba.HEXAGONAL, edge=edge) mamba.infFarNeighbor(prov2, imOut, 3, i, grid=mamba.HEXAGONAL, edge=edge) j = 3*i/2 mamba.infFarNeighbor(prov1, imOut, 2, j, grid=mamba.HEXAGONAL, edge=edge) mamba.infFarNeighbor(prov1, imOut, 5, j, grid=mamba.HEXAGONAL, edge=edge) def _sparseConjugateHexagonDilate(imIn, imOut, size, edge=mamba.EMPTY): """ Dilation by a conjugate hexagon. The structuring element used by this operation is not complete. Some holes appear inside the structuring element. Therefore, this operation should not be used to obtain true conjugate hexagons dilations (for internal use only). """ prov1 = mamba.imageMb(imIn) prov2 = mamba.imageMb(imIn) mamba.copy(imIn, imOut) val = mamba.computeMaxRange(imIn)[1]*int(edge!=mamba.EMPTY) for i in _sizeSplit(size): mamba.copy(imOut, prov1) j = 2*i mamba.supFarNeighbor(prov1, imOut, 1, j, grid=mamba.SQUARE, edge=edge) mamba.shift(prov1, prov2, 2, i, val, grid=mamba.HEXAGONAL) mamba.supFarNeighbor(prov2, imOut, 4, i, grid=mamba.HEXAGONAL, edge=edge) n1 = int(0.4641*size) n1 += abs(n1 % 2 -size % 2) n2 =(size -n1)/2 _sparseConjugateHexagonErode(imIn, imOut, n2, edge=edge) largeHexagonalErode(imOut, imOut, n1, edge=edge) def largeDodecagonalDilate(imIn, imOut, size, edge=mamba.EMPTY): """ """ n1 = int(0.4641*size) n1 += abs(n1 % 2 -size % 2) n2 =(size -n1)/2 _sparseConjugateHexagonDilate(imIn, imOut, n2, edge=edge) largeHexagonalDilate(imOut, imOut, n1, edge=edge) """ Dilation by a large diamond (losange) on square grid. This diamond is not completely filled. It is for internal use only. """ prov = mamba.imageMb(imIn) mamba.copy(imIn, imOut) for i in _sizeSplit(size): mamba.copy(imOut, prov) mamba.supFarNeighbor(prov, imOut, 1, i, grid=mamba.SQUARE, edge=edge) mamba.supFarNeighbor(prov, imOut, 3, i, grid=mamba.SQUARE, edge=edge) mamba.supFarNeighbor(prov, imOut, 5, i, grid=mamba.SQUARE, edge=edge) """ n1 = int(0.41421*size + 0.5) n2 = size -n1 largeSquareErode(imIn, imOut, n1, edge=edge) _sparseDiamondErode(imOut, imOut, n2, edge=edge) def largeOctogonalDilate(imIn, imOut, size, edge=mamba.EMPTY): """ Dilation by a large octogon (square grid). This operation uses dilations by large squares and large diamonds previously defined. """ n1 = int(0.41421*size + 0.5) n2 = size -n1 largeSquareDilate(imIn, imOut, n1, edge=edge) mamba.supFarNeighbor(prov2""" def _sparseDiamondDilate(imIn, imOut, size, edge=mamba.EMPTY): _sparseDiamondDilate(imOut, imOut, n2, edge=edge)
	Dilation by a large segment in direction dir in a reduced number of iterations. This operator is quite complex to avoid edge effects. largeLinearDilate(imOut, imOut, 5, size, grid=mamba.SQUARE, edge=edge)
	Uses the dilations by doublets of points (supposed to be faster, thanks to """ largeLinearDilate(imOut, imOut, 7, size, grid=mamba.SQUARE, edge=edge)

mamba.supFarNeighbor(prov, imOut, 7, i, grid=mamba.SQUARE, edge=edge) def _sparseDiamondErode(imIn, imOut, size, edge=mamba.FILLED):

"""

Acknowledgements

The author wishes to thank Nicolas Beucher for his kind and benevolent help in the achievement of this implementation. Nicolas wrote in particular the C and Python codes for the shift, supFarNeighbor and infFarNeighbor operators. He also debugged my ackward Python scripts, tested them and evaluated the performances of this implementation.


Nicolas Beucher

Serge Beucher

Hierarchical Queues: general description and implementation in MAMBA Image library

Introduction

Hierarchical queues (HQ) are extremely efficient structures for fast computation of morphological transformations such as watershed lines or geodesic reconstructions.

The recent development of MAMBA, new morphological image processing library, led the authors of this paper to review the definition and implementation of the watersheds and geodesic reconstructions algorithms based on these hierarchical queues. This library incorporates very fast implementations of these operators, based on HQ. These implementations increase deeply the performances of MAMBA in real-time applications.

Indeed, the documents describing HQ are quite succinct [START_REF] Beucher Serge | The Morphological approach of segmentation: the watershed transformation[END_REF][START_REF] Meyer Fernand | Procédé de traitement d'images par files d'attente hiérarchisées[END_REF]. Unfortunately, these descriptions are far from being accurate and complete. The HQ implementation as it is described in these documents will be used as a starting point for the description of its implementation in the MAMBA library. Regarding the watershed transform, two implementations are available: the first one where only the catchment basins are generated and the second one where these catchment basins are separated by the watershed lines. However, only the first implementation is described in the above mentioned documents. The second one is not addressed although it is more difficult to achieve. The two implementations however are available in some software libraries [START_REF] Rivest Jean-François | Xlim3D, an image analysis tool[END_REF][START_REF] Morph-M | Image Processing Software Specialized in Mathematical Morphology[END_REF].

The HQ implementation as it is realised in these libraries will be used as a starting point for the description of its implementation in the MAMBA library. Some of its peculiarities and defects will be recalled before indicating which enhancements have been achieved in MAMBA.

Regarding geodesic reconstruction, the description of HQ implementations in the above documents exhibits a peculiar feature which makes this algorithm to be tautological. We shall come back to this characteristic and we shall describe an implementation that avoids this problem. A hierarchical queue is the assembly of N simple queues (Figure 1 (a)). A queue is also called FIFO ("First In First Out") register. The queues contain tokens which are introduced and extracted according to a mechanism which will be described below. In each queue, tokens are extracted in the same chronological order as they were introduced. Each token represents a pixel of the image. The only information carried by a token is the coordinates (x, y) of its corresponding pixel in the image. There is no obligation for a pixel to be represented by only one token. Each queue has a single priority level. This priority level is generally linked to the grey values of the pixels, but it does not necessarily correspond to them directly. By convention, the lower the priority number, the higher the priority of the queue. Therefore, queue 0 has the highest priority (this convention comes from the fact that, in the watershed implementation which was the first HQ-based implementation, pixels of lower grey values are processed first). The number N corresponds to the number of grey levels in the image (256 grey levels in MAMBA library).

All queues are always open at their top: at any time, a token can be inserted in the corresponding priority queue. However, only the highest priority queue can be emptied and cleared: the token which is extracted from the highest priority queue is the one which arrived first in the queue (Figure 1 (b)). Finally, once the highest priority queue is empty, it is removed and the next queue can begin to empty. A lower priority queue cannot be processed if all the higher priority queues are not emptied. Conversely, when a queue is empty and is removed, it cannot be created again.

Figure 2: The four components of the structure: original image, label image, status image and HQ Three images are associated with a HQ (Figure 2). The first one is the image to be treated, f. This image defines the number of queues and the tokens priorities. The second image, called result image and denoted g both contains the initialisation state of the process and the result of the transformation when achieved. The status image s contains various tags (or flags) which control the status and the use of each pixel (and therefore of its corresponding token(s)) of the image.

Watershed transformation implementation by hierarchical queues (HQ)

Let us see how this general architecture is used to build watershed transformations (WTS).

Generating a WTS without catchment basins boundaries

Let us describe first how a WTS is realised with a HQ in the simplest case where only the catchment basins are built. This implementation is available in MAMBA (basinSegment).

Initialisation of the process

The image g is called label image in the WTS implementation. It is initialized with values corresponding to the labels of the different sources of flooding (minima or markers). Each source of flooding (connected component) is assigned a positive integer value (a label). It is not necessary that the different labels be consecutive values. You may also assign the same label to different connected components. In this case, these connected components will be considered as a single flooding source. Pixels which are not labelled are assigned a zero value in g.

The HQ is initialized by storing tokens corresponding to the pixels labelled in image g in their respective queues. In the case of the watershed transform (WTS), this priority level of each queue corresponds to the grey level of the token (pixel). The pixels of lower grey level have the highest priority (queues are numbered according to the grey levels, queue 0 has the highest priority). Remind that these tokens carry one single information: the coordinates of the corresponding pixel in f or g. Therefore, in the sequel, when there is no ambiguity, a pixel or its corresponding token will be given the same naming. The order of introduction of tokens in a queue is arbitrary. It generally corresponds to the image scanning order.

The status image s uses two tags:

• "CANDIDATE" indicates that the pixel has not been processed (not pushed on the queue) yet.

• "QUEUED" indicates that the pixel has been queued.

The status image his initialized with all pixels set to "CANDIDATE" except those corresponding to the labelled pixels in g which are set to "QUEUED".

HQ processing

The HQ operation is controlled by the following algorithm (pseudo-code) applied to each token: -Mark y a s QUEUED i n t h e s t a t u s image s .

} } }

At the end of the process, each catchment basin is tagged with a strictly positive label corresponding to the label of the source. The label image thus indicates, during the process, which catchment basin each pixel belongs to.

Example and discussion

This operation is illustrated in a simple case (one-dimensional).

Initial state: in this example, the flooding sources correspond to the minima of f. Tokens (b), (c) and (g) are pushed on the HQ in their respective priority queues. Image g is initialised with two labels (1 and 2). The status image contains CANDIDATE tags except for source pixels which are marked QUEUED.

First step of the process: token (b) is extracted, its only CANDIDATE neighbor (a) takes the same label in g, is pushed on the HQ (priority 1) and is tagged QUEUED in s. Same treatment is applied to token (c). Priority 0 queue is then empty and is removed.

Second step: token (g) is processed, it neighbors (f ) and (h) are placed on the HQ and tagged QUEUED. Token (a) has no CANDIDATE neighbor.

Next step: tokens (d) and (h) are extracted. Token (d) sends (e) on the queue. (h) pushes (i) on the current priority queue. So (i) is processed and pushes (j). The process ends up when (f ), (j) and (k) are processed.

Let us recall the importance of the rule for managing disappeared queues. Sometimes, indeed, lower priority tokens appear when the corresponding stack no longer exists. In this case, the token is inserted into the current stack. This situation occurs frequently with marker-controlled WTS, but not only in this case. Figure 3 shows the use of a HQ when initialisation is performed with any kind of markers and not with minima of f. Indeed, in this case, tokens appear, of higher priority than the current queue. Everything happens as if the grey value of the corresponding pixel was raised to the value corresponding to the current height of the flooding, which can be considered, somehow, as a change of homotopy "on the fly" [START_REF] Beucher | Segmentation Tools in Mathematical Morphology[END_REF]. However, this rule is also essential in the case of the construction of a simple WTS where it can happen that points with a height less than the current flood have not been taken into account yet (especially in the buttonholes structures, frequently encountered in real images, see [START_REF] Beucher | Segmentation Tools in Mathematical Morphology[END_REF]).

Strengths and weaknesses

This implementation is quite simple and therefore it is very fast. In MAMBA, the original image f is always a 8-bit (greyscale) image. The label image g is a 24-bit image and the status image s a 8-bit image. In fact, g and s are sharing respectively the three lower bytes and the higher byte of a 32-bit image. The depth of g allows to use up to 2 24 -1 markers (16,777,215) in the computation of the WTS. If a marker is assigned to each pixel of an image (a rather clumsy hypothesis by the way...), this means that the algorithm could manage images containing this number of pixels (typically 4096x4096 images for instance). Fortunately, in most applications, we need a much less number of markers. The computation time is proportional to the size of image f which must be a greyscale image. It is possible, however, to use this implementation to compute a watershed transformation on a 32-bit image (see section 6.1.1). remember also that this implementation is not isotropic as tokens queueing introduces a processing order which has important consequences on the result.

Watershed algorithm with real boundaries (watershed lines)

When catchment boundaries are produced, the watershed algorithm implementation is more complex. The status image and its tags are more sophisticated. The label image may also be more complex.

Note that, although this implementation has been widely used in past and recent libraries, to the best of our knowledge, no formal description of this algorithm exists in the literature. Therefore, this presentation will give us the opportunity to highlight some particularities in the initialisation process and in the algorithm operation.

Initialisation of the process

A label image g and a status image s are still used. The first one contains the catchment basins labels as previously and the other one tags indicating the status of image points. This status may be one of the following:

• "CANDIDATE", if the point is candidate to be inserted in the HQ.

• "QUEUED", if the point has already been inserted in the HQ.

• "RG LABELLED", if the point received a label corresponding to a catchment basin.

• "WTS LABELLED", if it belongs to the watershed line (WTS).

Compared to the previous algorithm, two more tags have been added. Their purpose is to specifically handle the status of the watershed points.

The label image contains strictly positive label values corresponding to the various catchment basins in progress and the value 0 corresponding to points not assigned yet. Associated with that image, there is a HQ.

The HQ and the status image initialisations are more complex than in the previous algorithm. The initialisation is performed by placing the initial marker points on the queue and by storing their corresponding labels in the label image. They are also tagged "RG LABELLED" in the status image. When these tokens come out of the queue, they are, by definition, already marked "RG LABELLED" and their neighbors can only have a "CANDIDATE", "WTS LABELLED" (surprising, but possible ...) or "RG LABELLED" status but with a label identical to the extracted point (in this version of the WTS with boundaries, original markers cannot be adjacent). All candidate neighbors must therefore be pushed on the HQ thus ensuring the spreading of the flood. The label of these points could possibly be replaced with a label. . . of same value! They will never become WTS points.

HQ processing, pseudo-code

Once initialised, the operation on this HQ is the following:

As l o n g a s t h e HQ i s n o t empty , do : { -E x t r a c t a t o k e n x from t h e HQ ( from t h e a c t i v e h i g h e r p r i o r i t y queue )

-For e a c h o f i t s n e i g h b o r i n g p o i n t s : { -I f t h e n e i g h b o r p o i n t i s "CANDIDATE" , i t s c o o r d i n a t e s a r e s t o r e d i n a b u f f e r . T h i s b u f f e r s i z e i s e q u a l t o t h e number o f n e i g h b o r p o i n t s ( 8 i n s q u a r e g r i d , 6 i n h e x a g o n a l g r i d ) and t h e r e f o r e c o n t a i n s p o t e n t i a l l y l i a b l e t o f l o o d i n g n e i g h b o r s . -O t h e r w i s e , i f t h e n e i g h b o r p o i n t i s "RG LABELLED" , t h e n :

{ -I f t h e e x t r a c t e d p i x e l i s "RG LABELLED" and i f i t s l a b e l i s d i f f e r e n t from t h e l a b e l o f t h e n e i g h b o r p o i n t , t h e e x t r a c t e d p i x e l i s t h e r e f o r e a WTS p o i n t and i t i s marked a s s u c h ( "WTS LABELLED ) . -O t h e r w i s e , t h e e x t r a c t e d p i x e l i s marked "RG LABELLED" and i t i n h e r i t s t h e l a b e l o f i t s n e i g h b o r i n t h e l a b e l image . } } -I f a t t h e end , when a l l n e i g h b o r s have been t a k e n i n t o a c c o u n t , t h e e x t r a c t e d p o i n t i s marked "RG LABELLED" , t h e p o i n t s o f t h e b u f f e r a r e i n s e r t e d i n t o t h e HQ and a r e marked "QUEUED" ( s o a s n ot t o be i n s e r t e d more than o n c e on t h e HQ) . }

Let us explain briefly this algorithm. Any token leaving the queue should have its status determined definitely: it is either "RG LABELLED" or "WTS LABELLED". But this status is determined only when all its neighboring points have been examined. Depending on the final status, its "CANDIDATE" neighbors can then be pushed or not to the HQ: if the pixel is "RG LABELLED", it is legitimate for its "CANDIDATE" neighbors to be pushed on the queue because they are adjacent to a flooded point and then, the flood will spread to these points. On the contrary, if the point belongs to a watershed line (it is "WTS LABELLED"), it cannot propagate flooding as it is not flooded itself. Its "CANDIDATE" neighbors must remain "CANDIDATE". This is acknowledged by storing its "CANDIDATE" neighbors in a temporary buffer and by pushing them to the HQ once it is certain that the central point is "RG LABELLED" and not "WTS LABELLED".

Note that the order of examination of the neighbors is not imposed and it depends on the implementation. In the following examples, we shall suppose that this order is clockwise starting from the upper right neighbor.

Note also that, in this algorithm, all the neighbors are processed even when the status of the extracted token has been fully determined by a previous neighbor.

Example of processing

Here is an example of the algorithm applied to the image in Figure 4, defined on an hexagonal grid (the same algorithm applies to a square grid image without any problem). The numerical values correspond to the value of each pixel (here the contrast is artificially increased). Pixels that are not markers (markers all have the value 2) are marked with a letter. The image depicted in Figure 5 is a buttonhole (see [START_REF] Beucher | Segmentation Tools in Mathematical Morphology[END_REF]). Lips of the buttonhole are indicated in green. Its entry point is point (n). The awaited watershed lines at the end of the process are in red.

Initialisation step. At left, the label image with labels in orange, blue and purple. These pixels correspond to the initial three markers (the colors correspond to the values of these labels). These pixels are pushed on the HQ during initialisation. They are all pushed in priority 2 queue (they have not all been represented) and are "RG LABELLED" in the status image (middle image, in green). The black dots in the status image are not yet included in the HQ ("CANDIDATE" points). At right, the corresponding status of the queue after initialisation.

Treatment of priority 2 queue. Each pixel of the queue is surrounded either by pixels of same label or by "CANDIDATE" pixels. Therefore, each token in the priority 2 queue pushes, when extracted, its candidate neighbors on the HQ. These points are labelled "QUEUED" in the status image (in grey). The label image is not modified.

Situation when priority 3 queue has been entirely processed: tokens (o) and (p) were extracted. (o) found point (n) as a candidate ("CANDIDATE") to the inclusion in the HQ. This point has been labelled "QUEUED" (grey), inserted on the HQ in the priority queue 4 corresponding to its value and, during the checking of all its neighbors, the status of point (o) itself was modified. It took the purple label. Finally (p) was extracted, it produced no insertion of a neighboring point (no "CANDIDATE" black neighbor) and it took the blue label.

Treatment of priority 4 queue. Tokens (q) (r) and (n) are extracted from the queue. They become WTS points. However, although (n) has point (k) as a candidate, it is not included in the HQ because point (n) is not a point marked "RG LABELLED" (it is not flooded). The queue processing ends (no points are inserted). The priority 5 queue is empty, the process continues with the priority 6 queue.

Treatment of priority 6 queue. Tokens (i), (j) (l) and (m) are extracted in this order. (i) is labelled (purple label) and, being labelled "RG LABELLED", it pushes on the HQ tokens (h) and (k). The priority 5 queue being disappeared, (h) and (k) are placed on the current queue.

(h) and (k) are extracted from the queue and labelled as WTS points. The rest of the process is similar and will produce the same result as above.

Processing the last queue (priority 7). Points (a), (c), (d), (e) and (f ) become WTS points. Point (b) has only two types of neighbors: WTS points and a single labelled point. It takes the label of this latter point (orange label). Note that the final result is different from the result expected above. It is however fully valid. The difference comes from the fact that points have not been processed in the same order in the HQ.

Strengths and weaknesses

This algorithm extends the principles established in the first algorithm (which computes only the basins) to allow the computation of the watershed lines. It remains fast (although not as fast as the first one) and straighforward.

Another important property of this implementation is that the WTS is idempotent (which was not true for some WTS implementations, with or whitout HQ, which can be found in some mathematical morphology software libraries).

Idempotence is a very important property of a watershed transformation. It means that, if we take the result of a watershed transform (that is the valued watershed function, see [START_REF] Beucher | Segmentation Tools in Mathematical Morphology[END_REF]) and if we apply again a watershed transformation on this function, the result is unchanged. This property is not always fulfilled in some watershed HQ implementations and a shift of the watershed lines may be observed. This defect is very annoying especially when the watershed transformation is used in hierarchical segmentations. However, we can guarantee that the watershed implementation in MAMBA is idempotent.

However, a remaining problem concerns points that are not labelled (because they are surrounded by WTS points). These points are never inserted in the HQ. Therefore, they are still labelled "CANDIDATE" at the end of the procedure. So, it is necessary to modifiy this tag (to change it to "WTS LABELLED") to obtain the final WTS result. This is done through a simple look-up table applied on all the pixels of the status image. Note that, as, in some applications, the result of the watershed transform is not used as such but after a threshold to get the binary image of the watershed lines, this final cleaning could be merged with the thresholding through an appropriate look-up transform. However, for the sake of simplicity, this cleaning procedure has been kept in the MAMBA library.

Example of configuration inducing a "thick" WTS (top left). For the sake of simplicity, label and status images have been mixed (color conventions are the same). At right, the colored dots correspond to different catchment basins and associated labels. The grey points are the points "QUEUED" stored on the HQ. Point (a) is not yet inserted. The extraction of grey points stored on the queue will cause the insertion of point (a) in the current priority queue (priority 6) because it is "CANDIDATE". All extracted points will be considered as WTS points because they are all surrounded by points with at least two different labels. Point (a) is then removed from the queue and it is not labelled because it is only surrounded by WTS points.

Summary and discussion

Both watershed transform implementations in MAMBA (basinSegment and watershedSegment) are efficient and straightforward. The initialisation process is simple and the HQ operation is in accordance with the general description given in section 2. Moreover, the watershed transform in MAMBA is idempotent. It is also fast and its computation time is more or less proportional to the number of pixels in the image as it can be shown in Figure 6. These values have been obtained on a 2.66 GHz processor running on Windows XP with different sizes of images. Only the watershed or catchment basins computation times have been calculated. The markers labelling has not been taken into account. As expected, the basinSegment transformation is faster than the watershedSegment one.

Reconstruction algorithms

HQ can also be used efficiently to perform geodesic reconstructions by dilations or erosions. A geodesic reconstruction is obtained by successive geodesic dilations of an image g, the marker image, under an image f, the mask image, until idempotence is reached (no evolution of the transform). Using HQ to perform this transformation allows a dramatic reduction of its computation time. The dual reconstruction is obtained by achieving successive erosions of image g above image f until idempotence.

Dual reconstruction by HQ

Let us start with the dual reconstruction using a HQ, where the priority management in the queues is similar to the one which is used in the WTS operators. We refer again to [START_REF] Beucher Serge | The Morphological approach of segmentation: the watershed transformation[END_REF] for the description of the algorithm where only this dual reconstruction (erosion-reconstruction) is described. Let us recall it before addressing the issue raised by the initialisation proposed in this description.

Initialisation

Dual reconstruction R * f (g) of a function f by a function g with g ≥ f, uses, as the WTS, a status image s. This status image (different from image g) will contain the various labels (or tags) indicating the status of the different points of image g during the reconstruction process. The label image g (also called marker image) contains the image to be rebuilt. This image will be modified and will hold the rebuilt image at the end of the reconstruction. Image f (also called mask image) remains unchanged. Three labels are possible for the points of the label image:

• "FINAL" indicates that the point has received its final value in the reconstruction.

• "QUEUED" indicates that the point has already been inserted in the HQ, but that it has not taken its final value yet.

• "CANDIDATE" means that the point has not been taken into account (this label corresponds normally to the value 0).

HQ handling in the geodesic reconstruction is very similar to its operation in the WTS, at least when the watershed lines are not drawn (catchment basins only). In particular, the priority of each queue is directly related to the grey values of image g, as it was the case in the WTS implementation. This means that pixels (and tokens) of lower grey value are processed first in this algorithm.

However, in the algorithm described in [START_REF] Beucher Serge | The Morphological approach of segmentation: the watershed transformation[END_REF], the HQ initialisation presents a strange singularity. The HQ initialization is made by storing at least one point of each regional minimum of function g. But nothing is said about how to extract those points belonging to regional minima of g. This is all the more disturbing since the geodesic reconstruction is conventionally used to extract extrema (minima and maxima) of an image. With this initialisation, the geodesic reconstruction by HQ seems tautological: to realise it, we must first detect points belonging to the minima of image g. But, usually (and in Mamba in particular), these points are detected by an reconstruction operator! Everything is happening as if, to perform a reconstruction using a HQ, we were obliged to first perform a mimina extraction by means of a geodesic reconstruction realised by another algorithm, a recursive oerator for instance (this operator exists in MAMBA but it is slower than the reconstruction by HQ, see section 4.3).

Fortunately, this level of complexity is not necessary. Indeed, it is essential that at least one point belonging to each minimum of g be included in the initial set of tokens. But we can add to this set any number of extra points without changing in any way the final result. In particular, we can initiate the process with all the points of g, which eliminates the screening step described above.

This ensues from a well-known property of the reconstruction (direct or dual). You can replace image g by any image g' (g' ≥ f ) provided that minima of g and g' are identical. One can even further relax the condition by considering functions g and g' such that the intersection of the connected components of their minima is pairwise non-empty (Figure 7). Figure 7: The two images g and g' provide identical dual reconstructions as their minima are superposed Another important point is to ensure that g ≥ f. As it is not always the case, this condition is enforced by using the image max(g(x), f (x)) in the initialisation phase which is therefore the following:

The first token (c) is removed from the queue of highest priority. It is marked "FINAL" (red) in the label image. The neighbor (b) of (c) is "CANDIDATE", it is stored in the HQ in the stack of priority max(g(c), f(b)) = 0, value also taken by g(b). The point (b) is marked "QUEUED" (green) in the label image. The same treatment is applied to point (d), neighbor of (c). We note that (b) and (d) are inserted in the queue whereas they were already present (on higher priority queues). We will see that this feature does not affect the result.

Token (b) exits from the queue, it is marked "FINAL". Point (a), "CANDIDATE" neighbor of (b) takes the value max(g(b), f(a)) = 0 and is inserted into the priority queue 0. It is marked "QUEUED" in the label image. Same treatment is applied for point (d), marked "FINAL" that pushes point (e) on the stack (marked "QUEUED") with priority g(e) = max(g(d), f(e)) = 4.

Token (a) leaves the stack (it is marked "FINAL"). Point (a) had no neighbors labelled "CANDIDATE", no points are inserted and the stack of priority 0 is empty.

The priority 1 queue is empty. Point (a) exits, it is already marked "FINAL" and it has no "CANDIDATE" neighbor. It cannot have one since it has been processed already and its neighbors have therefore already been taken into account. This shows that to insert a point several times in the HQ has no adverse consequences because, when it is reinserted, it is always in a lesser or equal priority queue than the queue it was previously. Point (j) is extracted. It takes the "FINAL" label and it pushes ponts (i) and (k) on the stack with the same priority. The image g and the label image are modified accordingly. We see that the point (k) is reinserted into the queue with the same priority.

Token (k) is extracted. This token was inserted in the initialisation phase, but it does not matter. We can already guess that the second token corresponding to (k) and located in the same queue will simply be removed when its turn comes. (k) is marked "FINAL" and has no "CANDIDATE" neighbors. Point (i) is also marked "FINAL" and pushes its "CANDIDATE" neighbor (h) on the priority 2 ((h) is marked "QUEUED"). g(h) takes the value 2. The second occurrence of (k) is simply removed.

Token (b) is exited. Token (h) comes out and pushes its neighbor (g) on the current priority queue. g(g) is amended.

Point (g) is processed and pushes its neighbor (f ) on the current priority queue. g (f ) is modified.

Point (f ) is processed and does not push any point in the HQ (no neighboring point is "CANDIDATE").

End of the process. Tokens (g), (h) and (i) are extracted. (e) is marked "FINAL". Finally, (d), (e) and (f ) are removed from the priority queue 5.

Strengths and weaknesses

The previous example shows that there is no problem to initialise the HQ with all the points of the original image g. In fact, this approach amounts to generate a recursive process in the label image and in the HQ since tokens are reintroduced thus replacing the first insertions and somehow cancelling them. Note that this algorithm requires a HQ larger than the number of image points. This size can theoretically be up to twice the size of the image (if a single minimum is present, which could, associated to a very specific configuration, lead to the reinsertion of all other points). We will see below how the HQ has been designed in MAMBA in order to simplify its management and to cope with this problem.

Due to the insertion, another consequence is that the processing time is not proportional to the number of pixels in the image, although the reinserted pixels processing is very fast in the second pass (since it is sufficient to unstack them). We could reduce the number of reinsertions and, thereby, reduce the HQ size and processing time by initialising the process with a reduced set of tokens as already explained above. However, this always needs some pre-processing (in order to generate this reduced set) which at the end decreases the performances of the overall process.

Reconstruction by HQ

The direct reconstruction R f (g) by dilation can also be performed with a HQ, mutatis mutandis. Indeed, the reconstruction R f (g) of image g under f is equal to R f (g) = m -(R * m-f (m -g) where m can be any revolving value (see [START_REF] Beucher | Geodesy and geodesic transformations[END_REF]). If m is equal to 255, (m -f) and (m -g) represent the inverted images. Therefore, the algorithm design for the dual reconstruction must be modified accordingly.

Initialisation

The initialisation step is identical for the reconstruction to the initialisation realised for the dual reconstruction. The HQ is loaded with tokens representing the pixels of min(g, f ) to ensure that g ≤ f. The only modification lies on the priority of the queues. The priority order must be reverted. Pixels of image g of highest values will be processed first. This means that the priority 0 queue (highest priority queue considering the conventional priority numbering) contains tokens corresponding to pixels with a grey value equal to 255.

HQ operation

The pseudo-code describing the HQ operation is also almost identical. Queues processing starts with highest priority queues which now contain tokens of highest grey values. The rule for changing g values becomes:

new v a l u e o f g ( y ) = min ( g ( x ) , f ( y ) )

This corresponds to the propagation of a geodesic dilation starting from the highest values of g.

Discussion: reconstruction by HQ versus its directional implementation

MAMBA offers two possibilities to perform a reconstruction (or its dual operation) on a 8-bit image, the directional reconstruction (build function in mambaComposed.geodesy.py module) based on the buildNeighbor operator and the hierarchical reconstruction (hierarBuild function in mamba.py module). The first version has the advantage of supporting all MAMBA image depths. It is built with optimised operators (buildNeighbor is a wrapper of MB BldNb8 that uses SSE2 instructions). Its main drawback is that its processing time is not proportional to the image size but to the complexity of the reconstruction. The HQ implementation does not use SSE2 instructions or specific optimizations but its characteristics make its running time almost proportional to the number of pixels in the image and not to the complexity of its reconstruction.

An example of complex image to be reconstructed is given in Figure 8. A single white dot at top left of the image is chosen as starting point for the reconstruction. The higher the number of zigzags, the higher the complexity of the reconstruction. If this result validates the theory, however, it is hardly representative of a real situation. To compare the two approaches in a realistic situation, we conducted performance tests on image minima detection using dual reconstructions. The results are shown in Figure 10. The control parameter is the depth of the minima.

These results are clearly to the benefit of the HQ-based algorithm.

MAMBA HQ implementation

The implementation of HQ algorithms in MAMBA is structured around tokens (MB Token structure in the C code). A token is a simple structure which points only to the next token, thus allowing to chain them together in threaded lists. A list then corresponds to a queue and becomes a simple structure indicating the first and last tokens of the queue (thus allowing to stack another token on its top). Tokens are created in a 2-dimensional array of the same size as the processed image which allows to link a token with its corresponding pixel through its position in the tokens array (Figure 11). The hierarchical queues algorithms were implemented in MAMBA for watershed transformations with two different functions. The first one, basinSegment, performs the watershed transformation without generating the watershed lines. The corresponding C Implementation (MB Basins function) is the simplest of those using a HQ in MAMBA. The second function, watershedSegment, generates the catchment basins and the watershed lines (basins boundaries). Its C implementation (MB Watershed function) is rather complex to ensure its idempotence. For both algorithms, the number of tokens stacked in the queue cannot exceed the number of pixels in the image. Moreover, there is at any time a single occurrence of a given pixel, i.e. of its corresponding token in the HQ.

Regarding the implementations of geodesic reconstructions by HQ in MAMBA, two functions were created, hierarBuild and hierarDualBuild, which perform respectively a reconstruction and its dual operation. C implementations (MB HierarBld and MB HierarDualBld functions) are relatively simple, although they introduce some tricks to solve the specific problems coming with these reconstruction algorithms by HQ (such as initialisation with all image points). In this case, some pixels can be introduced a second time in the HQ even though the initial introduction has not been processed. So we sometimes need two tokens per pixel. To be able to safely use these two sets of tokens, the tokens array is composed of two parts. The first one is dedicated to the HQ initialisation, whilst the queue processing and the neighbors tokens insertion are using the second part of the tokens array (Figure 12). During treatment, it is possible to have two tokens corresponding to the same pixel in the queue. We must therefore control, for each token extracted from the queue, that its associated pixel has not already been processed. Another feature of the reconstruction by HQ is that the result is a 8-bit image. So, contrary to the watershed transform where the label image and the result were stored in different byte planes of a single 32-bit image, we cannot mix the result and the label image in reconstruction by HQ algorithms. Therefore, the label image is defined in a separate table.

Improvements and extensions of HQ

This section is devoted to the description of some extensions and improvements which can be applied to the HQ operations described previously to address some specific domains. Among them, we shall focus our interest on 32-bit images and 3D images. Regarding 32-bit images, extensions of HQ watershed transforms and HQ geodesic reconstructions will be described. Concerning 3D images, only general ideas aiming at extending HQ structures to this class of images will be given.

Hierarchical queues for 32-bit images

Using the HQ algorithms on 32-bit images is at first challenged by a problem of memory requirement. Indeed, on 8-bit images, the range of pixel values is 0 to 255. Creating a HQ of 256 queues presents no difficulty on any modern computer or even handheld devices (smartphone,. . . ). 32-bit images however can potentially have pixel values range between 0 and 2 32 -1. Creating a HQ with more than 4 billions queues is simply impossible.

The solution is to split the processing of 32-bit images slice by slice using efficiently the HQ structure designed for 8-bit images. 32-bit images are cut into different slices of equal thickness. Each slice is then processed in the HQ structure designed for 8-bit images. Finally, the results of the processing of these different slices are stacked again to produce the expected transformation on the initial 32-bit image. Obviously, some appropriate rules must be enforced to obtain unbiased results. The computation time of the transformations is generally proportional to the range of values found in the 32-bit images. This time is not much longer than the computation time awaited with a direct 32-bit implementation because this implementation would also require to operate a hierarchical list proportional to this range.

Of course, this approach could be directly implemented into a low-level C function that would directly process 32-bit images thus simplifying the process. This would also allow more flexibility and speed (for example the thickness of the slices could be increase to 512, 1024,. . . instead of the "forced" 256).

Watershed transforms for 32-bit images

Let us describe first the extension to 32-bit images of the watershed transform where only catchment basins are generated. This algorithm has been made possible thanks to the fact that it is possible, in MAMBA watershed implementation, to stop the flooding process at any altitude.

Figure 13 Figure 13: Principle of realisation of a watershed transform by slicing of a 32-bit image with a 8-bit HQ Note the importance of stopping the flooding process at level 254. Indeed, if the flooding was continued to level 255, it would propagate on the flat zones added at this level. Therefore, the catchment basins would be adjacent and the entire process would end at the first step. Moreover, these catchment basins would be deeply biased as the topography above this altitude would not be taken into account.

A similar approach can be used to compute the watershed lines of a 32-bit image. The only difference is in the definition of the sources of flooding when computing the watershed transformation of slice i. To get them, the watershed lines obtained at the end of step i-1 are extracted (note that they are not necessarily closed) from the status image. The corresponding pixels are "WTS LABELLED". Then, these lines are subtracted from the label image. Therefore, this label image contains only the interior of the current catchment basins and the source pixels which have not started to pour water yet (since their altitude is higher than the current level). This modified label image is used as the new label image for the watershed transform at step i. At the end of the operation, watershed lines which have been removed at the beginning of this step are now regenerated again together with new watershed lines produced during the flooding of the pixels belonging to slice i. When all the slices have been processed, the status image contains the final watershed transformation of the initial 32-bit image.

Geodesic reconstructions for 32-bit images

This slicing technique can also be used to realise geodesic reconstructions by HQ on 32-bit images, as illustrated in Figure 14. The approach is even simpler than for the watershed transformation. We just have to cut the mask image and the marker image into slices, to replace the pixel values under the floor or above the ceiling of each slice by 0 and 255 respectively and to perform the reconstruction of the new mask function by the new marker one using the HQ algorithm. Then, the successive reconstructions just have to be added to get the final result.

Figure 14: Reconstruction operators by slicing on 32-bit images Note that this approach can be used equally for the geodesic reconstruction and for the dual reconstruction. Note also that it is not necessary to process the different slices in a specific order contrary to what is done with the watershed transformation. This is due to the fact that, in geodesic reconstructions, propagations are horizontal, whereas they are upstream in the watershed transform.

Remind that there exists already general reconstruction operators in MAMBA (build and dualBuild) which are based on a recursive algorithm and which can handle 32-bit images without any problem. Therefore, it is not sure that the algorithm described above will be faster than the recursive one for 32-bit images. It will depend on the complexity of the images, as already discussed in section 4.3, and also on the range of pixel values: the larger this range, the greater the number of slices.

Hierarchical queue for 3D images

Contrary to what was possible for 32-bit images, you cannot use the watershed and build operators present in MAMBA on 3D images. Indeed these types of images can be seen as sequences of 2D images but because the coded operators completely ignore the neighbors in the previous and the next image the result cannot be correct.

However the algorithms descriptions are still correct, only the definition of neighbor pixels is extended (and as for 2D images, depends on the grid used). Thus the same algorithms could be applied using specifically coded low-level C functions.

Conclusion

We tried in this paper to describe in a fairly exhaustive way both the principle of HQ operations and their implementations in the MAMBA image library.

The WTS and reconstructions implementations made in the MAMBA library brought original solutions for assessing the efficiency, accuracy and speed of the algorithms.

Regarding the watershed transformation (with apparent watershed lines), the initialization phase is simple and does not require pre-processing steps. In addition, the management of the neighbors of a token is made via a buffer memory, which reduces bias and increases processing speed by reducing breaks in the data stream. The final post-treatment could certainly be avoided. However, it is not very penalising and does not increase significantly the processing time.

The HQ reconstruction has also been greatly improved by simplifying the initialization step which now consists simply in filling the HQ with all pixels of the image to be reconstructed. The initialization as it was succinctly presented in previous documents describing this reconstruction was indeed rather strange since it required features (the extrema) that are generally obtained by means of algorithms that use or are very similar to reconstructions. The queue structure designed for watershed transforms is also used in the HQ reconstruction. Its duplication allows to assign more than one token to the same pixel.

The implementation of the hierarchical queue in MAMBA was performed using a threaded list built in the form of a two-dimensional array of tokens which allows a one-to-one correspondence between each token and the pixel that it represents. This structure opens the way for new uses of the HQ for increasing the performances of other morphological transformations.

2

  Reminder : Hierarchical Queue general design Let us start by recalling what hierarchical queues are and how they are used.

Figure 1 :

 1 Figure 1: Simple queue (a) and hierarchical queue (b)

{

  For e a c h p r i o r i t y queue , from p r i o r i t y 0 t o N, do : { While t h e queue i s n o t empty , do : { -E x t r a c t a t o k e n x from t h e HQ -Determine t h e n e i g h b o r i n g p i x e l s o f x which a r e CANDIDATE i n s -For e a c h n e i g h b o r y w i t h a CANDIDATE tag , do : { -A s s i g n t o y i n image g t h e same l a b e l a s x . -I n s e r t t h e t o k e n y i n t o t h e queue w i t h p r i o r i t y c o r r e s p o n d i n g t o t h e g r e y l e v e l o f y i n f ( i f i t e x i s t s ) o r i n t o t h e queue o f h i g h e s t p r i o r i t y s t i l l e x i s t i n g .

Figure 3 :

 3 Figure 3: Management rule for the lower priority tokens

Figure 4 :

 4 Figure 4: Original image used to illustrate the algorithm

Figure 5 :

 5 Figure 5: Button hole configuration

Figure 6 :

 6 Figure 6: Computation time of the watershed transformation with HQ for various image sizes

Figure 8 :

 8 Figure 8: Examples of more and more complex reconstructions

Figure 9 :

 9 Figure 9: Reconstruction speed versus complexity

Figure 10 :

 10 Figure 10: Computation time for the minima extraction versus their depth

Figure 11 :

 11 Figure 11: Hierarchical queue structure. Each queue has an entry in a hierarchical list indicating the first and last token stored in the queue. The token array is a threaded list which allows to scan the entire queue. Coordinates (-1, -1) for the next token indicates the end of the queue

Figure 12 :

 12 Figure 12: Double hierarchical queue implementation for the geodesic reconstruction. The grey tokens correspond to the same pixel (same position in the two token arrays) and have been stored twice in the HQ

  illustrates the principle of the algorithm. The original image is split into successive slices. The first slice contains all the pixels with a value in the range [0, 255]. Pixels with values greater than 255 are also given the value 255. The catchment basins of this function are built but the flooding is stopped at level 254. Then a second slice is defined, comprising pixel values in the range [254, 509]. These values are shifted in the range [0, 255]. Pixel values higher than 509 take the value 255. The catchment basins of this new function are built up to the level 254. The flooding sources are simply the catchment basins produced at the previous step. The procedure is repeated until the maximum value of the initial 32-bit image has been treated (that is, contained in a slice).

This patent and all its international versions are lapsed since November

26, 2008.

For any p i x e l x i n t h e image : { -A s s i g n p o i n t x a new v a l u e i n image g e q u a l t o : g ( y ) = max ( g ( x ) , f ( x ) ) The t o k e n c o r r e s p o n d i n g t o p o i n t x i s t h e n i n s e r t e d i n t o t h e p r i o r i t y queue g ( x ) and t h e p o i n t i s l a b e l l e d "CANDIDATE" i n t h e l a b e l image . }

HQ processing for dual reconstruction

Once initialised, the hierarchical queue is processed according to the following algorithm: While t h e HQ i s n o t empty , do : { -R e t r i e v e a t o k e n x from t h e HQ ( h i g h e s t p r i o r i t y queue ) . { -I f t h e l a b e l o f x i n d i c a t e s t h a t x d o e s n o t a l r e a d y have i t s f i n a l v a l u e ( l a b e l o t h e r than "FINAL " ) , g i v e

x t h e l a b e l "FINAL" , i n d i c a t i n g t h a t i t s v a l u e i s now f i n a l .

-For e v e r y n e i g h b o r y o f x w i t h a "CANDIDATE" l a b e l ( z e r o l a b e l ) , a s s i g n a new v a l u e a t t h i s p o i n t y i n g image e q u a l t o : g ( y ) = max ( g ( x ) , f ( y ) )

The t o k e n c o r r e s p o n d i n g t o p o i n t y i s t h e n i n s e r t e d i n t o t h e g ( y ) p r i o r i t y queue and t h e p o i n t i s l a b e l l e d "QUEUED" i n t h e l a b e l i m a g e l . } }

This algorithm consists simply in performing a geodesic erosion of each pixel of g by starting from the lowest ones and by propagating this erosion to neighboring pixels. This propagation is very fast as it avoids to take into account unnecessary pixels.

Example of dual reconstruction with a HQ

We will illustrate this approach by using again the one-dimensional example provided in [START_REF] Beucher Serge | The Morphological approach of segmentation: the watershed transformation[END_REF]. Although this example is very simple, there is absolutely no difference in the HQ operation when 2D images are processed. We represent image g at each step (image f, in grey, is unchanged), the status image and the HQ.

Initialisation: all points of g are stored on the HQ in the stack corresponding to their value. The status image only contains "CANDIDATE" labels (colored in white).