About a problem of definition of the geodesic erosion

Serge Beucher

Centre de Morphologie Mathématique

Mines ParisTech September 2011

Introduction, context

This note can be seen as a complement of the course on geodesic transforms [START_REF] Beucher | Geodesy and geodesic transformations[END_REF]. It aims at clarifying an issue in the definition of the numerical geodesic erosion and how this definition is inferred from the definition of the geodesic dilation. We will see that the definition of the numerical geodesic erosion (on functions) is not an extension of the set geodesic erosion. We explain the reason for this difference. This induces some difficulties when it comes to define operators built from the elementary operators as the numerical geodesic opening and closing. In fact, with the classical definition of the numerical geodesic erosion, this is simply not possible! So we will explain the logical progression used to obtain a correct definition of the numerical geodesic operators from their definition on sets. This will allow both to understand where the problems lie but also to define a numerical geodesic erosion which is a true extension of the geodesic erosion given for sets. Most definitions and proofs will be provided in a digital frame. These definitions could be introduced in a continuous context, however this would complicate them unnecessarily.

Reminder: geodesic erosion and dilation for sets

Let us recall the definition of the elementary geodesic dilation of a set Y X ( Y ) included in the geodesic space X. We have:

(1) X ( Y ) = ( Y / B ) 3 X where B is the elementary euclidean ball (size-1 hexagon or square according to the digital grid in use). This operator is obviously extensive and increasing. From the geodesic dilation, we can define the elementary geodesic erosion (of size 1) X ( Y ) of the set Y included in the geodesic space X by using the duality between these two operators with respect to the set complementation in X. We get the following definition:

(

2) X ( Y ) = X # X ( X # Y ) = X 3 [(( X # Y ) / B ) 3 X ] c = (( Y 4 X c ) 0 B ) 3 X
The larger operations are obtained by iteration:

(3

) X n ( Y ) = n X . X ... X ( Y )
Note that in the case of geodesic operators, dualities by complementation or by adjunction [START_REF] Henk | Morphological Image Operators[END_REF] lead to identical definitions (because, in a geodesic context, the only structuring elements in use are geodesic balls).

The set geodesic erosion is increasing (the geodesic space X being fixed) and anti-extensive. Let us recall also the geometric interpretation of these two operators. The geodesic dilation of size n of Y corresponds to the space swept by the geodesic balls of size n whose center is included in Y. It is also the locus of the centers of the geodesic balls of size n in X that intersect Y. The geodesic eroded set corresponds to the locus of the centers of the geodesic balls of size n of X included in Y:

(4) X n ( Y ) = x c X : B X ( x, n ) _ Y (5) X n ( Y ) = x c X : B X ( x, n ) 3 Y ! -= 4 ycY B X ( y, n )

Geodesic opening and closing for sets

In the set framework, it is easy to define the geodesic opening and closing of a set Y included in a geodesic space X. We have:

(6 X n ( Y ) corresponds to the set of geodesic balls of size n included in Y. The geodesic closing corresponds to the complement in X of the geodesic opening of the complement of Y in X. This is an extensive operator. This operator is also idempotent. Indeed:

) X n ( Y ) = X n . X n ( Y ) (7) X n ( Y ) = X n . X n ( Y )
(8

) X n . X n . X n . X n ( Y ) = X n . X n . X n ( Y ) = X n . [ X n . X n ( Y )]
The expression between brackets is the geodesic closing of the geodesic erosion of Y, .

X n ( Y ) This closing, extensive, contains the geodesic eroded set. Therefore we have:

(the geodesic dilation is increasing)

X n . X n . X n . X n ( Y ) q X n . X n ( Y ) As
(increase and anti-extensivity of the geodesic opening), the

X n . X n . X n . X n ( Y ) _ X n . X n ( Y ) property is proven.
It is of course also true for the geodesic closing. Therefore, these operators can legitimately be called opening and closing.

Application example

The geodesic opening and closing are not much used. The use of the opening by reconstruction or of the closing by reconstrution is much more common. However, we can give a simple example of use of the geodesic closing (Figure 1). This example shows a set Y made of several connected components included in a set X. The geodesic closing allows to merge the close enough connected components of Y without deep modifications of the other components.

The geodesic dilation for functions

Extending the geodesic dilation for sets to the numerical domain is easy. Let f and g be two greyscale images (functions). Assume that . In order to define the elementary f [ g geodesic dilation of f in the geodesic space defined by g, we just apply the classical g ( f ) procedure which consists in performing a geodesic dilation on all the sections (thresholds) of the subgraphs of f and g followed by a piling of the results to obtain the subgraph of . g ( f ) We know that this approach works only if the set operator is increasing, otherwise the re-stacking does not generate the subgraph of a function (Figure 2):

X i ( f ) = x : f ( x ) m i i < j X j ( f ) _ X i ( f ) increasing: ( X j ) _ ( X i ) X j ( ) _ X i ( )
In the case of the numerical geodesic dilation, this stacking generates a subgraph because the set geodesic dilation is increasing not only with respect to Y, the set to be dilated, but also with respect to X, the geodesic space. Indeed, by setting: , , , (9) X i ( f ) = x : f ( x ) m i X j ( f ) = x : f ( x ) m j X i ( g ) = x : g ( x ) m i X j ( g ) = x : g ( x ) m j With i < j, we have:

(10) X i ( g ) q X j ( g ) q X j ( f ) and:

(11) X i ( g ) q X i ( f ) q X j ( f ) Therefore this leads to:

(12) X i ( g ( f )) = ( X i ( f ) / B ) 3 X i ( g ) q ( X i ( f ) / B ) 3 X j ( g ) q ( Xj ( f ) / B ) 3 X j ( g ) = X j ( g ( f )) We can then easily give a specific algebraic formulation of the geodesic dilation of a function f in the geodesic space defined by g (we say that the geodesic dilation of f is carried out under g). We have: : (13) g ( f ) = ( f / B ) . g = inf ( f / B, g ) (the transition from set operators to their numerical equivalents is straightforward). The extension of the geodesic dilation to the numerical domain is immediate thanks to the good increasing property of the set operator with respect to the geodesic space. We can give a geometric interpretation of this operator on the subgraphs of f and of g by introducing the concept of geodesic cylinder. A geodesic cylinder defined in the C g ( x, n ) subgraph G(g) of a function g, with origin x, point of the subgraph, and with a radius n, is made of points of G(g) at a geodesic distance less than or equal to n from every point y belonging to the vertical semi-axis (Figure 3). The subgraph of the geodesic dilated ] -∞, x ] function is therefore made of all the geodesic cylinders of G(g) whose origin is in G(f). g ( f )

The numerical geodesic dilation is obviously increasing and extensive. We have obviously B 1 _ B 2 _ B 3 .

The numerical geodesic erosion, analysis of a mix-up

The numerical geodesic erosion is generally defined from the numerical geodesic dilation by introducing a duality by inversion around a revolving value [START_REF] Beucher | Geodesy and geodesic transformations[END_REF]. We define the elementary numerical geodesic erosion (size 1) of the function f above the function g (we indeed have

) as the operation:

f m g (14) g ( f ) = p -p-g ( p -f ) = sup ( f 0 B, g )
where p is any revolving value (the result is independent of this value). This transformation, however, has nothing to do with a possible extension of the set geodesic erosion to greyscale images! Indeed, in this definition, we have surreptitiously changed the geodesic space. We switched from the space G(g), subgraph of g, to its ongraph. It is as if, in the set case, we had defined the geodesic erosion in the geodesic space X c (Figure 4). Y would then contain X and we would define the elemntary geodesic erosion as: X ( Y ) (15) X ( Y ) = ( Y 0 B ) 4 X Note also that the sections of the elementary numerical geodesic erosion are defined with this formula. This commonly used definition of the numerical geodesic erosion (formula ( 14)) is indeed interesting. In particular, it allows to introduce an operator such as the dual geodesic reconstruction which is widely used in Mathematical Morphology. However, as erosion and dilation do not apply to the same geodesic space, defining a geodesic opening or closing by concatenating them makes no sense. Indeed, consider the transformation:

(16

) g ( f ) = g . g ( f )
We see immediately that this operator is equal to (assuming that, when , the geodesic f [ g erosion is, by convention, equal to the identity):

g ( f ) (17) g ( f ) = g . g ( f ) = g ( f )
This operator is the geodesic dilation (the geodesic erosion has no effect), so it does not verify the idempotence property and it is therefore difficult to argue that it is a closing. A solution consisting in using the result of the dilation as a geodesic space for the erosion of the function g (to avoid the above mentioned problem by ensuring that the function to be eroded is above) does not bring any improvement. Indeed, if we write:

(18

) g ∏ ( f ) = g (f) ( g )
The operator thus defined is always greater than or equal to . g ( f ) Thus we see that, if we transpose in the numerical case the previous example given for sets (see Figure 1), by replacing the sets X and Y by two functions f and g defined by (valued indicator functions):

f ( x ) = m if x c Y f ( x ) = 0 otherwise g ( x ) = M if x c Y ( M ú m ) g ( x )
= 0 otherwise the operator does not produce the same result as in the set case since it stops after the g ( f ) dilation. Can we avoid this problem and define a geodesic erosion in the same numerical geodesic space as the dilation? The answer is yes, just by using the duality by adjunction [START_REF] Henk | Morphological Image Operators[END_REF] and the geodesic cylinders. Let us consider the subgraphs et of the functions f and g (with . We can G ( f ) G ( g ) f [ g) define the geodesic erosion of G(f) by geodesic cylinders of size n, , as the union of

C g ( x, n ) all sets Z included in G(g) such that their geodesic dilation is included in G(f): (19) g [ G ( f )] = 4 Z : G(g) n ( Z ) _ G ( f ) This formulation is equivalent to: (20) g [ G ( f )] = 4 x c G ( g ) : C g ( x, n ) _ G ( f )
The geodesic erosion of f under g, has the previous set as its subgraph.

g n ( f ) (21) G g n ( f ) = g [ G ( f )]
This definition, however, is not very operative, like all definitions based on the adjunction. It is therefore desirable to express this operator in a algebraic form. To do this, we will first write the relationship between each section of and the sections et of X i g n ( f ) X j ( f ) X j ( g ) the functions f and g.

If a point x of G(g) at altitude i belongs to

, the cylinder is included in X i g n ( f ) C g ( x, n ) G(f). Therefore, each section at altitude j of this cylinder corresponds to a geodesic ball ( j [ i ) of size n included in G(f), thus in . The intersection of the axis of the cylinder X j ( f ) ] -∞, x ] and of thus belongs to the geodesic erosion of in . We can therefore C g ( x, n ) X j ( f ) X j ( f ) X j ( g ) write:

(22) X i g n ( f ) = 3 j[i X j (g) n ( X j ( f ))
To belong to the section X i at altitude i of the geodesic erosion of f under g, a point x must also belong to all the set geodesic erosions of the sections located at an altitude in X j ( f ) j < i the corresponding section . X j ( g ) We recognize, in this formulation, an already known technique to force an operator to be increasing when it is extended to the numerical domain. This approach is particularly used to extend the concepts of thinning and thickening [START_REF] Beucher | Segmentation d'Images et Morphologie Mathématique[END_REF]. Note that there are two ways for obtaining this result (Figure 5). We can either intersect all sections less than or equal to i or perform the union of all sections greater than or equal to i. In this present case, the first solution is obvious while the latter one prevails for the thinnings and thickenings (the first approach leads to the definition of other operators called over-thickening under-thinning, cf. [1, pages 125-127]). Starting from the above formula, we can show that the elementary digital geodesic erosion of f under g can be written:

( f ) = [( f -m ) 0 B ] . g (23) g
where m is the indicator function defined by: m

( x ) = v max if f ( x ) = g ( x )
m ( x ) = 0 otherwise (assuming f and g positive and bounded:

)

. 0 [ f [ g [ v max
The proof of this formula is quite long, it was postponed in Appendix A. This operator is increasing and anti-extensive. It is increasing because, assuming , we f [ f ∏ have . Indeed, the locus of points where is included in the set of points where m

[ m ∏ f ∏ = g . It is also anti-extensive because: f = g (24) ( f -m ) V B [ f -m (25) ( f -m ) . g = ( f . g ) -( g . m ) = f -( g . m ) = f -( f . m ) = f Q.E.D.
Figure 6 illustrates the various stages of construction of this transformation. This operator produces a result identical to the set operator when applied to indicator functions.

Geodesic opening and closing

By concatenating the operators and , we can then define two new operators g ( f ) g ( f ) called geodesic opening and closing:

(26) g ( f ) = g . g ( f )

(27) g ( f ) = g . g ( f ) These operators are increasing.

is anti-extensive. Indeed, the subgraph is made of g ( f ) G ( g ) the union of all the geodesic cylinders whose centers are included in . All thes G g cylinders are included in G(f). Hence we have . g ( f ) [ f It can be shown, mutatis mutandis, that is extensive. These two operators are also g ( f ) idempotent. This can be easily proven in the same way as the proof for sets. This property comes from the increase and extensivity/anti-extensivity properties of the operators et .

This new geodesic erosion

, defined in the same geodesic space as , makes it g ( f ) g ( f ) possible to define true openings and closings exhibiting properties and behavior which are similar to the set operators. Thus, in the previous example (Figure 1), with indicator functions instead of sets, the reader will verify that he/she obtains the same result.

Conclusion

This extension to the numerical case of the set geodesic erosion had, to my knowledge, never been proposed before. The problem raised by the surreptitious change of geodesic space had certainly been seen before but no solution had been provided. It is true that the geodesic opening and closing operators are not widely used. Thus, in the example given above of a "driven" merging of connected components, we can simply use a geodesic dilation to solve the problem in a satisfactory way. However, a new tool leads sooner or later to new ways to use it. It is thus possible to define more complex numerical operators (gradients, top-hat transforms, thickenings, thinnings, etc.) in geodesic spaces defined by functions. In particular, these operators could demonstrate their usefulness for applications where multiple imaging modalities are available. Finally, this operator shows both the interest of adjunction to define a dual operator but also its limitations when, from this general definition, we need to design a more operative one.

Appendix A -algebraic formulation of the geodesic erosion of f under g

First recall some definitions and formulas linking the sections of a subgraph of a function and the algebraic addition and subtraction operators. These definitions are first given in the continuous case before transcribing them in the discrete case. Let us call : X ( f ) = x : f ( x ) m The (topological) closed section at level of the subgraph G(f). In the same way: Y ( f ) = x : f ( x ) > Is the (topological) open section at the same altitude. We have:

X ( f + g ) = x : ( f + g )( x ) m ( f + g )( x ) m e f ( x ) m et g ( x ) m -, ≤ X ( f + g ) = 4 X ( f ) 3 X -( g ) Or: X ( f + g ) = 4 X -( f ) 3 X ( g )
In the same way:

X ( f -g ) = 4 X -( f ) 3 X ( -g ) But: X ( -g ) = x : g ( x ) [ -= Y - c ( g ) X ( f -g ) = 4 X -( f ) 3 Y -c ( g ) X ( f -g ) = 4 X + ( f ) 3 Y c ( g )
We can also write:

Y ( f + g ) = X - c ( -f -g ) Y ( f + g ) = 4 X --( -f ) 3 X ( -g ) c Y ( f + g ) = 3 X -- c ( -f ) 4 X c ( -g ) Y ( f + g ) = 3 Y + ( f ) 4 Y -( g )
That is:

Y ( f + g ) = 3 Y -( f ) 4 Y ( g )
In the same way:

Y ( f -g ) = 3 Y -( f ) 4 Y ( -g ) Y ( f -g ) = 3 Y -( f ) 4 X -c ( g ) Y ( f -g ) = 3 Y + ( f ) 4 X c ( g )
In the discrete case (integer functions), we have:

Y i ( f ) = x : f ( x ) > i = x : f ( x ) m i + 1 = X i+1 ( f )
We can prove that:

X i ( f + g ) = 4 j X i-j ( f ) 3 X j ( g ) X i ( f -g ) = 4 j X i+j ( f ) 3 X j+1 c ( g )
Or:

X i ( f -g ) = 3 j X i+j ( f ) 4 X j c ( g )
Using formula (2), we can write:

Appendix B -Summary of the various geodesic operators

This appendix summarizes the various elementary geodesic erosion and dilation operators available. We show that there is another dilation operator which has not been highlighted yet in this document. We also attempt to define operations allowing to go from an operator (eg erosion) to another (the dilation). Only the numerical formulation of the operators will be considered. One can notice that despite a different formulation, the set geodesic erosion is identical to the numerical X ( Y ) operator . Indeed, the mask M (equivalent for sets to the function m) corresponds to g ( f )

points in space which belong to both Y et X or to their complements. M is therefore the complement of the symmetrical difference Δ between X and Y:

because M = ( X < Y ) c = ( X 3 Y ) 4 ( X c 3 Y c ) = Y 4 X c
Y _ X So, we have:

Y 4 M = Y 4 X c which is the initial definition of the elementary geodesic erosion for sets.

The following table summarizes the various elementary operations applicable to a function f in a geodesic space defined by g.

g ( f ) = ( f / B ) . g g ( f ) = [( f -m ) 0 B ] . g with m = x : f = g f [ g g ( f ) = ( f 0 B ) -g g [ f

Dilation Erosion

The analysis of this table brings two immediate observations: -The erosion is deduced, as we already knew, from the dilation by inversion. The g ( f ) g ( f ) two functions f and g are revolved, then the dilation is carried out. Finally the negation of the result provides the erosion. Writing the inversion of a function, we have:

f f / B . g = f 0 B . g If we invert this result: f 0 B . g = ( f 0 B ) -g = g ( f )
-The definition of an elementary geodesic dilation of f over g is missing. This definition would be the counterpart of the erosion when . However, this dilation can be built g ( f ) f m g by inversion starting from the erosion . It is denoted (beware of the confusion g ( f ) g ( f ) betweeen the highlighting of this operator and of the inversion). We have: 7 illustrates this operation. The table of elementary geodesic operators can then be completed as follows:

with g ( f ) = f -m 0 B . g m = x : f = g That is: f -m 0 B . g = f -m / B -g = [( f . m ) / B ] -g But: m = x : f < g = x : f > g We arrive to the following definition: with g ( f ) = [( f . m ) / B ] -g m = x : f > g Figure
g ( f ) = ( f / B ) . g g ( f ) = [( f -m ) 0 B ] . g with m = x : f = g f [ g g ( f ) = [( f . m ) / B ] -g with m = x : f > g g ( f ) = ( f 0 B ) -g g [ f Dilation Erosion
We changed the notation of these operators to make it more consistent. The geodesic erosion and dilation of f under g are respectively are denoted and while operators applied

g ( f ) g ( f )
to f over g are denoted and . g ( f ) g ( f )

Only the transformations et , on the one hand, together with et on the

g ( f ) g ( f ) g ( f ) g ( f )
other hand, are linked by an inversion operator. It is however not possible to concatenate these transformations to obtain openings and closings. Conversely, this is completely possible with and or with and . We obtain two openings, one of f under g g ( f ) g ( f ) g ( f ) g ( f ) and the other of f over g:

g ( f ) = g . g ( f )
Two equivalent definitions exist for the closing:

g ( f ) = g . g ( f ) g ( f ) = g . g ( f )
Nevertheless, it is not easy to exhibit a duality operator for building from or g ( f ) g ( f ) g ( f ) from . This comes in particular from the fact that a mask m is involved in some g ( f ) definitions whereas it does not appear in others. We can however modify the previous definitions of and . These new definitions (see the table below) use slightly g ( f ) g ( f ) modified definitions of the masks in order to enhance the symmetries between the various operators.

g ( f ) = [( f . m ) / B ] . g with m = x : f [ g g ( f ) = [( f -m ) 0 B ] . g with m = x : f m g f [ g g ( f ) = [( f . m ) / B ] -g with m = x : f > g g ( f ) = [( f -m ) 0 B ] -g with m = m : f < g g [ f Dilation Erosion
Of course, the masks introduced in et are "fictitious" (in the erosion, it is empty g ( f ) g ( f )

and it covers the whole space in the dilation). However, these new definitions better show the internal structure of these operators and their similarities. Nevertheless, achieving an erosion using a simple operator applied to the dilation in order to have a more operative bridge than the duality by adjunction between these two transforms is a difficult task, at least in the numerical case. Given that we have direct formulations to achieve these transformations, the interest of such a bridge is finally quite low..

Appendix C -Geodesic operators in the MAMBA Image library

This appendix is intended for the Mamba library users. The module "Geodesy" of the library has been modified. It now contains the newly defined elementary geodesic operators. Thus, there are two geodesic erosions and two geodesic dilations. These are:

-upperGeodesicErode corresponding to . g ( f ) -lowerGeodesicErode corresponding to .

Figure 1 :

 1 Figure 1 : Example of geodesic closing of a set Y in a geodesic space X. On the left, the initial sets, on the right, geodesic closing of set Y.

Figure 2 :

 2 Figure 2 : Principle of the extension to the numerical domain of set increasing operators

Figure 3 :

 3 Figure 3 : Geodesic cylinders defined in the subgraph of a function g. On the right, sections of a bidimensional geodesic cylinder. Each section corresponds to a geodesic ball.We have obviously B 1 _ B 2 _ B 3 .

Figure 4 :

 4 Figure 4 : Set equivalent of the geodesic erosion for functions. In this case, we have . Y q X

Figure 5 :

 5 Figure 5 : Two approaches used to generate the subgraph of a function by stacking sets obtained with a non increasing set operator ψ . The procedure on the right (intersection of sections at an altitude lower than or equal to i) is used for the geodesic erosion. The procedure using the union of upper sections (on the left) is not currently used (it is used, however, for the thinnings).

Figure 6 :

 6 Figure 6 : Example of elementary geodesic erosion of a function f under a function g. Top, left: the initial functions. Top, right: the m function (indicator of the equality of f and g). Bottom, left: euclidean erosion of the function . Bottom, right: geodesic erosion of f fm under g.

  g g

Figure 7 :

 7 Figure 7 : elementary geodesic dilation of a function f over g.

  These two operations are increasing. The geodesic opening is anti-extensive (any point of is included in a geodesic ball of size n itself included in Y). The open set

X n ( Y )

  The former geodesicErode and geodesicDilate operators were kept for compatibility reasons with the previous version of Mamba. The geodesicDilate operator is nothing but an alias of lowerGeodesicDilate, while geodesicErode corresponds either to lowerGeodesicErode or to upperGeodesicErode depending on the depth of the images, binary or greyscale.Here is the Python source code of these new operators. You can find it in the geodesy.py module. Performs a upper geodesic erosion of image 'imIn' above 'imMask'.The result is put inside 'imOut', 'n' controls the size of the erosion. 'se' specifies the type of structuring element used to perform the computation (DEFAULT_SE by default). Performs a lower geodesic erosion of image 'imIn' under 'imMask'. The result is put inside 'imOut', 'n' controls the size of the erosion. 'se' specifies the type of structuring element used to perform the computation (DEFAULT_SE by default).The binary lower geodesic erosion is realised using the fact that the dilation is the dual operation of the erosion. (imIn, imMask, imOut, "inf") for i in range(n): mamba.generateSupMask(imOut, imMask, imWrk2, False) mamba.convertByMask(imWrk2, imWrk1, 0, mamba.computeMaxRange(imWrk1)[START_REF] Beucher | Segmentation d'Images et Morphologie Mathématique[END_REF]) mamba.logic(imOut, imWrk1, imOut, "sup") mC.erode(imOut, imOut, se=se) mamba.logic(imOut, imMask, imOut, "inf") def geodesicErode(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE): """ This transformation is identical to the previous version and it has been kept for compatibilty purposes.

	for i in range(n):
	mamba.diff(imOut, imMask, imOut)
	mC.dilate(imOut, imOut, se=se)
	mamba.logic(imMask, imOut, imOut, "sup")
	else:
	imWrk1 = mamba.imageMb(imIn)
	imWrk2 = mamba.imageMb(imIn, 1)
	g mamba.generateSupMask(imOut, imMask, imWrk2, True) for i in range(n): ( f )
	-upperGeodesicDilate corresponding to mamba.convertByMask(imWrk2, imWrk1, 0, mamba.computeMaxRange(imWrk1)[1]) . Warning! 'imMask' and 'imOut' must be different. g ( f ) -lowerGeodesicDilate corresponding to . g mC.dilate(imOut, imOut, se=se) ( f ) mamba.logic(imOut, imWrk1, imOut, "inf") """
	mamba.logic(imOut, imMask, imOut, "sup") if imIn.getDepth() == 1:
	These operators work with binary, 8-bit and 32-bit images. They are incorporated in the latest mamba.diff(imMask, imIn, imOut)
	version of Mamba (version 1.1.1). def lowerGeodesicDilate(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE): lowerGeodesicDilate(imOut, imMask, imOut, n, se=se)
	""" mamba.diff(imMask, imOut, imOut)
	Performs a lower geodesic dilation of image 'imIn' below 'imMask'. else:
	The result is put inside 'imOut', 'n' controls the size of the dilation. imWrk1 = mamba.imageMb(imIn)
	'se' specifies the type of structuring element used to perform the imWrk2 =mamba.imageMb(imIn, 1)
	""" This module provides a set of functions to perform geodesic computations using Mamba based functions. It includes build and dualbuild operations, geodesic erosion and dilation, computation of maxima and minima... computation (DEFAULT_SE by default). Warning! 'imMask' and 'imOut' must be different. """ def geodesicDilate(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE): """ This operator is simply an alias of lowerGeodesicDilate. It is kept for compatibility reasons. mamba.logicNote that the binary and the greytone operators are different.
	it works with imageMb instances as defined in mamba. """ """
	""" lowerGeodesicDilate(imIn, imMask, imOut, n, se=se) if imIn.getDepth() == 1:
	lowerGeodesicErode(imIn, imMask, imOut, n, se=se) # Contributors: Serge BEUCHER, Nicolas BEUCHER else:
	import mamba """ from mambaCore import ERR_BAD_DEPTH def upperGeodesicErode(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE): upperGeodesicErode(imIn, imMask, imOut, n, se=se)
	import mambaComposed as mC
	def upperGeodesicDilate(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE):
	"""
	Performs a upper geodesic dilation of image 'imIn' above 'imMask'. Warning! 'imMask' and 'imOut' must be different. The result is put inside 'imOut', 'n' controls the size of the dilation. 'se' specifies the type of structuring element used to perform the """
	computation (DEFAULT_SE by default). mamba.logic(imIn, imMask, imOut, "sup")
	Warning! 'imMask' and 'imOut' must be different. for i in range(n):
	"""
	mamba.logic(imIn, imMask, imOut, "sup")
	if imIn.getDepth() == 1:

mamba.logic(imIn, imMask, imOut, "inf") for i in range(n): mC.dilate(imOut, imOut, se=se) mamba.logic(imMask, imOut, imOut, "inf") mC.erode(imOut, imOut, se=se) mamba.logic(imOut, imMask, imOut, "sup") def lowerGeodesicErode(imIn, imMask, imOut, n=1, se=mC.DEFAULT_SE):

"""

g ( f ) = g . g ( f )

X j (g) ( X j ( f )) = X j ( f ) 4 X j c ( g ) 0 B 3 X j ( g ) The section at level i of the elementary numerical geodesic erosion can be written:

We have . All the sections of the function m are identical and correspond to the same set f [ g containing all the points x where . This set corresponds to the section :

Let us determine the section at level i of the function . We have:

That is:

Term (a) can be written:

Term (b) is written:

Which gives:

We can write:

Which proves formula (23).


Serge Beucher

Towards a unification of waterfalls, standard and P algorithms

Introduction

This paper can be considered as an addendum to the previous paper entitled "P algorithm, a dramatic enhancement of the waterfall transformation" [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF]. In this former document, some suggestions have been made here and there, some ideas have emerged. These suggestions and ideas have unfortunately not be developped due to some lack of time and also because the authors did not wish to complicate the notions and concepts which were introduced.

The purpose of this present paper is to come back to these suggestions, to explore them in more details. We shall see that this exploration will lead to a generalisation of many segmentation transformations as watersheds, waterfalls, standard and P algorithms. This generalisation will give birth to a general segmentation procedure which unify the previously mentioned tools and which depends on two parameters only.

This paper aims also at describing an implementation of this general segmentation tool in the MAMBA library. It is a matter of fact that, since the publication of the P algorithm, a great number of requests have been forwarded to the authors regarding a usable implementation of this P algorithm. This MAMBA implementation brings an answer to these requests.

Notations

Let us recall the notations and operators which are used in this document.

Operators

-w is the valued watershed transform. The result is a function equal to 0 inside the catchment basins and to the altitude of the points belonging to the watershed lines.

-h is the hierarchical image generator. This operator can be applied to any image (positive function) f, but in this paper, it will be used mainly on segmentations (see below).

Images

-s i denotes the final valued segmentation for the hierarchical level i.

-s' i+1 denotes the initial valued segmentation for the hierarchical level i. Generally, s' i ≠ s i as the hierarchical process at level i aims at modifying this initial segmentation. s 0 corresponds to the initial valued watershed (most of the time, it is a gradient watershed, but not always).

-h' i and h i are the hierarchical images corresponding respectively to s' i and s i :

h' i = h(s' i ) ; h i = h(s i ) -m is a numerical mask. It is an image made of two kinds of pixels: pixels at 0 or pixels with a value equal to the maximal possible value which can be defined depending on the depth of the image (for instance, with a 8-bit image, this value is 255). We often denote m as m(<condition>). When, at point x, <condition> is true, m takes its maximum value. When <condition> is false, m is equal to 0 at point x.

-Binary images (sets) can also be used: S i , S' i and M. They correspond to the pixels of the above numerical images which are different from 0.

The various steps towards a unification

Remind that the algorithms described in [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF] share the same general principle: reintroducing some of the contours which have been removed at each step of a waterfall transform (Fig. 1).

Fig. 1: Reminder of the three types of contours which are at stake in the waterfalls, standard and P algorithms. Type-1 contours should be always preserved, which is not the case in the classical waterfalls transformation.

An enhanced waterfalls transformation

Warning! In [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF], standard and P algorithms were referenced as "enhanced waterfalls transformation". However, the enhanced waterfalls operator which is presented here is different. So, from now on, the expression "enhanced waterfalls" will designate only this operator.

Among the various types of contours which are analysed, it is a matter of fact that type-1 contours are always significant and should be kept. Indeed, standard and P algorithms both preserve these contours. However, it is not the case for the classical waterfalls transformation. Therefore, before entering more sophisticated procedures, it is easy to define an enhanced waterfalls transformation where type-1 contours are systematically preserved. In fact, this operator is defined as follows. If s i is the current level of segmentation in the hierarchy, we can define its corresponding hierarchical image h i . Then, the valued watershed w of h i is realised. Finally, the corresponding hierarchical image h i+1 = h(w(h i )) is computed. The contours of s i which are higher than or equal to h i+1 belong to the new segmentation level s i+1 .

This enhanced waterfalls transformation has some advantages. Firstly, it is not necessary to cope with possible biases occurring with the classical operator and which are described in [1, page 25]. Therefore, this algorithm is faster, easier to implement, especially when we look for the complete waterfalls hierarchy. It is, in particular, useless to calculate again the hierarchical image at each step. Indeed, the hierarchical image h i+1 , although it has been calculated from the initial segmentation at level i+1 can be used at the next step as its minima contain all the minima, with a one-to-one correspondence, of the hierarchical image that would be built from the final segmentation s i+1 . So, at the next step, their respective watershed transforms are identical. Moreover, the resulting segmentation produced by this enhanced waterfalls transform is much better, as it is illustrated on some examples at Fig. 2.

Note that, on both cases however, the number of hierarchies is the same. Note also that, in this case, the stopping criterion is the idempotence. When h i+1 is equal to 0, all the contours of s i are reinjected and therefore s i+1 = s i .

Fig. 2: Examples of classical waterfalls transforms (middle images) and of enhanced waterfalls transforms (right images) calculated from the valued watersheds of the gradient of the original images (on the left). Only the last levels of hierarchy are given.

The reintroduction of type-1 contours dramatically enhances the final result. But this does not avoid the removal of some significant contours simply because they are lower than the final hierarchical image. These contours correspond to the maxima-islands introduced in [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF] which appeared in the first levels of hierarchy and which were removed by the watershed transform (which, as a semihomotopic operator, is not controlled by the maxima). This can be observed in the second example: the road marking (white dotted line) has disappeared because it corresponds to a maximum-island with a rather low height in the hierarchisation process.

A first generalisation of standard and P algorithms

Standard and P algorithms deal with the reintroduction of type-2 contours. Remind that these contours are defined as contours (belonging either to the previous segmentation or to the initial one) which are closer to the current hierarchical image than to 0. Although it is claimed that sandard and P algorithms are non parametric transforms, the definition of type-2 contours leads, de facto, to the introduction of a parameter equal to 2, as it is discussed in [1, §8.2]. As suggested in this discussion, a first generalisation of these algorithms simply consists in replacing this hidden parameter by an explicit one, λ, called gain. Therefore, the definition of type-2 contours is modified accordingly. A type-2 contour is a contour C k of the current segmentation s i such that its initial altitude a 0 (C k ) fulfils the following inequality:

λ a 0 (C k )>h i+1
where h i+1 is the hierarchical image corresponding to the initial segmentation s' i+1 . λ can be any real positive value. In this extension, as it was already the case in the previous definition of P and standard algorithms, the relationship between s 0 and h i+1 is a strict inequality:

λ s 0 >h i+1
If the gain λ is strictly greater than 1, this definition includes type-1 and type-2 contours. Indeed type-1 contours correspond to points of s 0 such that λ s 0 ≥h i+1 (large inequality). Therefore, multiplying the altitude of these points by λ >1 insure that the strict inequality is fulfilled. However, when λ = 1, this definition does not correspond to the enhanced waterfalls transform described above. Indeed, some type-1 contours (those which are strictly equal to h i+1 ) are removed. This is very annoying since type-1 contours must always be preserved. To cope with this problem, one solution would be to replace the strict inequality by a large one:

λ s 0 ≥h i+1
But this solution arises a new difficulty when λ is greater than 1. It tends to reinforce the influence of low contrasted contours in the first levels of hierarchy. For instance, when λ =2, contours at altitude 1 are preserved when they are surrounded by at least one contour at altitude 2. This gives to these contours a too important weight which induces negative effects on the following hierarchies and on the quality of the final result. So, the function which h i+1 should be compared to is s 0 , if s 0 ≥h i+1 or λ s 0 if λ s 0 >h i+1 .

In the digital case, these two conditions can be mixed by defining a composite function. We have:

λ s 0 >h i+1 → λ s 0 ≥h i+1 +1 →(λ s 0 -1)≥h i+1
Therefore, the function used for the comparison is defined as (λ s 0 -1)∨s 0 and the condition which must be fulfilled to preserve contours is:

(λ s 0 -1)∨s 0 ≥h i +1
This extension can be implemented in the following way:

• Compute (λ s 0 -1)∨s 0 (s 0 is the initial valued watershed transform).

• From segmentation s i (step i of the hierarchisation), compute its corresponding hierarchical image h i : h i = h(s i ). • Compute the initial segmentation s' i+1 : s' i+1 = w(h i ).

• Compute the hierarchical image h' i+1 from s' i+1 : h' i+1 = h(s' i+1 ).

• Compute m ((λ s 0 -1)∨s 0 ≥h ' i+1 ) .

• Remove the open (broken) contours in m and compute its corresponding binary mask M.

• For standard algorithm, compute M '=M∩S i and the corresponding numerical mask m'.

• For P algorithm, compute M '=M∩S 0 and the corresponding numerical mask m'.

• The final segmentation s i+1 is given by: s i+1 =m '∧(s i ∨h i+1 ) As for the enhanced waterfalls transformation, the stopping criterion is the idempotence. This can be easily proved. If s N is the last non empty segmentation, we have s' N+1 = w(h N ) = 0 and h N+1 = 0 (as h N contains a unique minimum). Therefore, M covers the entire space and s N+1 = s N .

The trick described above shows that it is possible to use a large range of comparison functions (as it was already suggested in [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF], page 60). The general comparison inequality will then be, at each step:

ψ(s 0 )≥h i+1
where ψ is any function such that ψ≥I (extensive function). If ψ is not extensive, the result would be irrelevant. ψ can also be any spatial function. This means that the comparison of each point of the contours to the current hierarchical image can be controlled not only by the height of the point but also by its position in the image. This allows to introduce external information for a better control of this hierarchical segmentation. For instance, trends, illumination variations, regions of interest, textured regions could be taken into account. This extension offers a large scope of possibilities. Now, the real problem is to define correctly ψ in order that these features could be used.

Note that, although it may be indifferent to use h i+1 or h' i+1 in the standard segmentation (between h' i+1 and h i+1 , no new minima appear), it is not true in practice. This is due to parity problems in the implementation of the watershed transform similar to those described in [1, page 63]. However, in P algorithm, using h i+1 is compulsory as new minima may be introduced at each level.

Note also that M∩S 0 and M are identical for P algorithm. The intersection is not necessary.

Fig. 3: When some points of a removed contour (FOZ lower than the hierachical image in blue) are higher than the hierarchical image (left), they must be suppressed otherwise they will be used for building the next hierarchical image which is likely to be erroneous (right).

You may also notice that a step of the implementation consists in removing broken contours in the mask function. These broken contours appear after the comparison as some points of these contours may be higher than h' i+1 although their lowest altitude is lower than the hierarchical image (Fig. 3). Thus, these contours are not completely removed. However, they must be suppressed otherwise the next hierarchical image would be erroneous. The fastest way to remove these contours is through a watershed transform of the mask function m. This is faster than a clipping procedure and, moreover, it is not depending of the complexity of the broken structures. Although it seems that this step is necessary only when working on valued watershed functions, in fact, it is not true. This procedure must be applied also on mosaic-gradient image. Indeed, if the points of the initial contours have the same height, it is not the case when we deal with higher levels of hierarchy since a single region may be surrounded by contours at different altitudes.

The last level of hierarchy of the CAR image in the standard segmentation is presented at Fig. 4 for three different values of the gain λ: 1.5, 2.0 and 3.0. The number of remaining contours increases as λ increases. This is obvious as, whatever the value of the gain, the number of hierarchy levels remains the same (4 in this example). Therefore the unique minimum in the last hierarchical image has always the same value. As the standard segmentation is nothing but a threshold of the previous remaining contours at an altitude proportional to 1/λ of the value of this minimum, this explains the result. In the same way, Fig. 5 illustrates some P algorithm results obtained with different values of λ (1.5, 2.0 and 3.0). Contrary to the standard segmentation, there is not necessarily an increasing of the remaining contours when λ increases. This is due to the fact that, contrary to the standard algorithm, new minima may appear in the successive hierarchical images. For the same reason, the number of hierarchical levels may vary as λ varies.

A further extension

Introducing the gain λ allows to bridge the gap between the standard segmentation and the enhanced waterfalls transform. Another extension (which was also suggested in [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF], page 44) brings the missing link between the standard and P algorithms. Indeed, the only difference between standard and P algorithms lies on the fact that, in standard algorithm, only contours belonging to the previous hierarchy are taken into account whereas, in P algorithm, all the contours (that is contours belonging to the first level of hierarchy -the initial watershed) are considered. When computing the (i+1)th level of hierarchy in standard algorithm, contours of S i are used, while contours of S 0 are considered in P algorithm. So, it is possible to define a hierarchical segmentation where contours belonging to an intermediary level of hierarchy would be taken into account. This is achieved by introducing an offset k between the current level of hierarchy and the level which is used to select the contours to be compared to the current hierarchical image.

Implementing this extension consists in using simply S j = S i-k instead of S i or S 0 in the computation of the mask image M '=M∩S j . When k = 0, the standard segmentation is performed. When k = +∞ (sufficiently high in practice -255 for a greyscale 8-bit image), P algorithm is realised. Note that if (i -k) is negative, Sj (which does not exist) is replaced by S 0 . It is equivalent to define S j as: S j =S max(i-k ,0) Note also that, when k > 0, some contours may oscillate in the successive hierarchies. However, due to the way the successive hierarchies are stored and coded (see below), this oscillation is most of the time masked and consequently has no effect on the hierarchies. Fig. 6 shows the results (last levels of hierarchy) of this extended segmentation procedure applied to image TOOLS for different values of the gain λ and the offset k. When no image is given, this simply means that no change occurred since the previous segmentation obtained with a lower offset. Although it is difficult to bring into general from a single example, it seems nevertheless that the best results are obtained when the offset is maximum and for values of the gain around 2, which tends to validate the heuristic choice of these parameters made in [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF].

Mamba implementations and performance

All the segmentation procedures described above have been implemented with the Mamba Image library. This open source library (X11 license), written in C and Python, is entirely devoted to mathematical morphology. It is freely available (version 1.0) with its documentation [START_REF] Beucher | MAMBA documentation[END_REF] at http://www.mamba-image.org. These new operators are not implemented in the current version (they will be in the next one). Therefore, the Python source code of these operators is given in annex, together with some add-ons which are needful to use them with the current version of Mamba.

The purpose of this implementation is twofold. Firstly, it answers some requests made to the author for an efficient and simple implementation of P algorithm. Unfortunately, until now, the available implementations of waterfalls and P algorithms use non open mathematical morphology software libraries. This implementation fills a gap. Secondly, it provides a strong basis for further experimentations with these tools, especially for tests on public image databases.

Description of the operators

Here is a brief description of the operators available in this Mamba implementation (see the annex for details).

hierarchy generates the hierarchical image from the initial image (imIn) which can be a valued watershed or a mosaic-gradient or even a simple gradient. Image imMask (which is binary) contains the watershed lines of the initial image. The hierarchical image is put in imOut. The algorithm used to build the hierarchical image has been described in [START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF] and [1, page 20]. This operator is used by the other segmentation processes.

-hierarchicalLevel computes the next level of hierarchy from the current one, imIn (which is a valued watershed or a mosaic-gradient). From imIn, its hierarchical image is built, then the watershed transform is computed. Finally, a correction is applied (its description is given in [1, page 25]) in order to be sure that the next hierarchy is embedded in the previous one. This operator is only used by the classical waterfalls generator. It is useless for the other segmentation processes as the above mentioned correction is no more necessary.

waterfalls is the classical waterfalls algorithm. From the initial image imIn (valued watershed or mosaic-gradient), it iterates the hierarchicalLevel operator to obtain the successive hierarchies of the waterfalls transformation in imOut. Each point in imOut takes the index value of the hierarchy where it appears for the last time (index 0 denotes the first level of hierarchy).

-enhancedWaterfalls implements the enhanced waterfalls transform where type-1 contours are always reintroduced. This implementation is faster than the waterfalls algorithm as the correction used in the previous approach (in hierarchicalLevel) is no longer necessary. Contours which are partly below the hierarchical image at each step are totally removed by means of a watershed transformation which acts as a skeleton by influence zones (SKIZ).

-standardSegment performs the standard segmentation with a gain λ as parameter (it is equal to 2 by default). The operator starts by computing the image equal to (λ s 0 -1)∨s 0 , where s 0 is the initial valued watershed image in imIn. At each step, this image is compared to the current hierarchical image and contours which are above this hierarchical image and which were still present in the previous segmentation level are kept. The others are removed by a watershed transform.

-segmentByP computes P algorithm with a gain λ (2 by default). This algorithm is very similar to the previous one. The only difference (but it is of outmost importance) comes from the fact that all the initial contours are considered, not only those still present in the previous step of hierarchy. This small difference induces dramatic changes in the result of the segmentation.

-generalSegment mixes and extends the standard and P algorithms by introducing the offset k which allows to choose the hierarchy level used in the comparison with the current hierarchical image. More precisely, the offset indicates how far in the past we go to select the contours which will be compared. generalSegment could be the only implemented operator. Indeed, if the offset is equal to 1, this operator is identical to standardSegment. If offset is equal to 255, it is identical to segmentByP. Note also that the smallest value for the offset is 1 and not 0 in the software implementation (offset is defined from the hierarchical level to be built and not from the last already built one).

-extendedSegment is the last implemented operator. It can be considered as experimental. Contrary to the others, the function used in the comparison (namely (λ s 0 -1)∨s 0 ) is not computed at the beginning of the process. Instead, a function ψ (see above) must be determined and entered in the operator. It corresponds to the imTest image. This extended segmentation should allow to take into account variations of luminance in the background, to reduce the influence of textured regions or other features which may deteriorate the quality of the final segmentation.

All these operators store the successive hierarchies in the imOut image and they also return the number of hierarchy levels. Each level of hierarchy can be obtained by simply thresholding the imOut image.

Performances

These operators are built directly on the images. Although Mamba is fast, this could lead to slow operations compared to the performances which could be obtained by using approaches based on graphs [START_REF] Marcotegui Beatriz | build(imWrk1, imWrk3) se = mC.structuringElement(mamba.getDirections(grid), grid) mC.dilate(imWrk3, imWrk1, 1, se) mamba.diff(imWrk2, imWrk1, imWrk1) mamba.logic(imWrk1, imWrk3, imWrk1, "sup") mamba.convertByMask(imWrk1, imWrk0, 255, 0) mamba.logic([END_REF]. However, these operators widely use the watershed transform and the geodesic reconstruction which are implemented in Mamba by means of hierarchical queues [5,6]. These transformations are fast, so the computation speed of the hierarchical segmentation operators is not so bad as it is shown in the table below. These values have been obtained on 256x256 grey scale images with a dual core Intel processor running at 2.66 GHz on Windows XP. The average computation time per hierarchy level is less than 30 ms (Fig. 7). A better value could be obtained by using a more powerful processor and with a Linux operating system (mainly because it is a 64-bit OS). Preliminary tests indicate that a speed increase of about 50% can be expected. On the other side, it is difficult to be sure that an implementation with graphs would be more efficient because this implementation needs to generate the graph from the image at the beginning of the process and to do the reverse operation at the end. Anyway, at this time, this graph implementation does not exist for P algorithm yet.

Image

Conclusion and next steps

We have shown, through the introduction of two parameters, that various hierarchical segmentation algorithms belong to the same class of operators. By simply varying these two operators (gain λ and offset k), it is possible to sweep continuously a whole range of transformations, from the enhanced waterfalls transform (λ = 1, k = 0) to P algorithm (λ = 2, k = +∞) to standard segmentation (λ = 2, k = 0).

[5] MEYER Fernand. Procédé de traitement d'images par files d'attente hiérarchisées. French Patent n° 91/033841, Institut National de la Propriété Industrielle, 20 mars 1991 (original document and its english translation can be found in the esp@cenet database).

[6] BEUCHER Nicolas and BEUCHER Serge. Hierarchical Queues: general description and implementation in MAMBA Image library. Web publication available at http://cmm.ensmp.fr/ ~beucher/publi/HQ_algo_desc.pdf, April 2011.

This document is copyrighted under the Creative Commons "Attribution Non-Commercial No

Derivatives" license. For terms of use, see http://creativecommons.org/about/licenses/.

Annex

This annex contains the source scripts of the various operators described above. These scripts are somewhat redundant as the general structure of each operator is very similar. However, we preferred to keep this redundacy in order to clearly show the modifications which allow to extend these different operators.

To run these scripts, you will need to install Mamba, version 1.0. This software library is available at http://www.mamba-image.org. Read carefully the documentation before installation (Mamba requires other Python libraries, in particular the Python Imaging Library).

Important! These scripts require that you add in existing Mamba modules some operators which do not exist yet in version 1.0. You need in particular the multiplication of an image by a real constant value. This operator is available below. Just copy the Python script and add it at the end of the miscellaneous.py module:

def mulRealConst(imIn, v, imOut, nearest=False): """ Multiplies image imIn by a real positive constant value v and puts the result in image imOut. inIn and imOut can be 8-bit or 32-bit images.

If imOut is a greyscale image (8-bit), the result is saturated (results of the multiplication greater than 255 are limited to this value). The constant v is truncated so that only its two first decimal digits are taken into account.

If 'nearest' is true, the result is rounded to the nearest integer value.

If not (default), the result is simply truncated. """ if imIn.getDepth()==1 or imOut.getDepth()==1: mamba.raiseExceptionOnError(mambaCore.ERR_BAD_DEPTH) imWrk1 = mamba.imageMb(imIn, 32) imWrk2 = mamba.imageMb(imIn, 1) v1 = int(v * 100) if imIn.getDepth()==8: imWrk1.reset()

Fig. 4 :

 4 Fig.4: Results of the standard algorithm applied to the valued watershed of the gradient of image CAR (left) for different values of the gain λ (from left to right, 1.5, 2.0, 3.0).

Fig. 6 :

 6 Fig.6: Last levels of hierarchy for various gains (horizontally) and offsets (vertically).

Fig. 7 :

 7 Fig. 7: Computation time for the generation of the successive hierarchical levels in P algorithm. The average time for computing one level is about 29 ms with the configuration described above.

This whole procedure will not be necessary with the next version (1.1) of Mamba which will contain these new operators.

We have also provided a Mamba-based implementation of this class of transformations which is fast enough to allow the exploration of its efficiency, its properties and its characteristics. The first tests which have been realised reinforce the interest for P algorithm and for its parametric version. The introduction of the offset, although it seems promising, validates, at first, the choice of an infinite offset for potentially providing the best results for the segmentations. However, this path, together with the use of a general comparison image (extended segmentation) need to be further explored.

Introducing parameters arises the problem of the choice of their optimal values in particular for the gain λ. Here again, a gain equal to 2.0 produces very often good results. However, it would be interesting to verify this point and, moreover, to find ways to automatically determine its best value. To this regard, the next step will consist in using image databases to test these operators and to try to assess their robustness. Instead of using the Berkeley database as it was the case in [1], database which provides a two important semantic content so that it is difficult to compare segmentations provided by human beings and by a low level algorithm, we intend to use a simpler database, the Segmentation Evaluation Database (available at http://www.wisdom.weizmann.ac.il/~vision/ Seg_Evaluation_DB/index.html) which seems more appropriate to test the efficiency of these segmentation algorithms.

mamba.copyBytePlane(imIn, 0, imWrk1) else: mamba.copy(imIn, imWrk1) mamba.mulConst(imWrk1, v1, imWrk1) if nearest: mamba.addConst(imWrk1, 50, imWrk1) mamba.divConst(imWrk1, 100, imWrk1) if imOut.getDepth() ==8: mamba.threshold(imWrk1,imWrk2,255,mamba.computeMaxRange(imWrk1) [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF]) mamba.copyBytePlane(imWrk1, 0, imOut) imWrk2.convert(8) mamba.logic(imOut, imWrk2, imOut, "sup") else:

mamba.copy(imWrk1, imOut)

Note also that, if you want to apply these operators on mosaic-gradient images, you will have to patch this transformation in the module segment.py. Indeed, this operator is faulty in Mamba, version 1.0. To apply the patch, simply copy the script below and replace the corresponding function in segment.py.

def mosaicGradient(imIn, imOut, grid=mamba.DEFAULT_GRID): """ Builds the mosaic-gradient image of 'imIn' and puts the result in 'imOut'. The mosaic-gradient image is built by computing the differences of two mosaic images generated from 'imIn', the first one having its watershed lines valued by the suprema of the adjacent catchment basins values, the second one been valued by the infima. Here is below the main script. Just copy it entirely and store it in a file named hierarchies.py. Then, you simply need to run this script to have all the hierarchical operators available in the Image library.

""" This module provides a set of functions to perform hierarchical segmentations operations using mamba. it works with imageMb instance as defined in mamba. This module contains the waterfalls algorithm and various hierarchical operators (enhanced waterfalls, standard hierarchy and P algorithm).

" """ General standard segmentation. This algorithm keeps the contours of the watershed transform which are above or equal to the hierarchical image associated to the next level of hierarchy when the altitude of the contour is multiplied by a 'gain' factor (default is 2.0). This transform also ends by idempotence. All the hierarchical levels of image 'imIn'(which is a valued watershed) are computed. 'imOut' contains all these hierarchies which are embedded, so that hierarchy i is simply obtained by a threshold [i+1, 255] of image 'imOut'. 'imIn' and 'imOut' must be greyscale images. 'imIn' and 'imOut' must be different. This transformation returns the number of hierarchical levels.

""" imWrk0 = mamba.imageMb(imIn) imWrk1 = mamba.imageMb(imIn) imWrk2 = mamba.imageMb(imIn) imWrk3 = mamba.imageMb(imIn) imWrk4 = mamba.imageMb(imIn, 1) imWrk5 = mamba.imageMb(imIn, 1) imWrk6 = mamba.imageMb(imIn, 32) mamba.copy(imIn, imWrk1) mC.mulRealConst(imIn, gain, imWrk6) mC.floorSubConst(imWrk6, 1, imWrk6) mamba.threshold(imWrk6, imWrk4, 255, mamba.computeMaxRange(imWrk6) [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF]) mamba.copyBytePlane(imWrk6, 0, imWrk0) mamba.convert(imWrk4, imWrk2) mamba.logic(imWrk0, imWrk2, imWrk0, "sup") mamba.logic(imWrk0, imWrk1, imWrk0, "sup") imOut.reset() nbLevels = 0 mamba.threshold(imWrk1, imWrk4, 1, 255) flag = not(mamba.checkEmptiness(imWrk4)) while flag: hierarchy(imWrk1, imWrk4, imWrk2, grid=grid) mamba.add(imOut, imWrk4, imOut) mC.valuedWatershed(imWrk2, imWrk3, grid=grid) mamba.threshold(imWrk3, imWrk5, 1, 255) flag = not(mamba.checkEmptiness(imWrk5)) hierarchy(imWrk3, imWrk5, imWrk2, grid=grid) mamba.generateSupMask(imWrk0, imWrk2, imWrk5, strict=False) mamba.logic(imWrk4, imWrk5, imWrk4, "inf") mamba.convertByMask(imWrk4, imWrk3, 0, 255) mamba.logic(imWrk1, imWrk3, imWrk3, "inf") mamba.negate(imWrk4, imWrk4) mamba.label(imWrk4, imWrk6, grid=grid) mamba.watershedSegment(imWrk3, imWrk6, grid=grid) mamba.copyBytePlane(imWrk6, 3, imWrk3) mamba.logic(imWrk1, imWrk2, imWrk1, "sup") mamba.logic(imWrk1, imWrk3, imWrk1, "inf") mamba.threshold(imWrk1, imWrk4, 1, 255) nbLevels += 1 return nbLevels def segmentByP(imIn, imOut, gain=2.0, grid=mamba.DEFAULT_GRID):

""" General segmentation by P algorithm. This algorithm keeps or reintroduces the contours of the initial watershed transform which are above or equal to the hierarchical image associated to the next level of hierarchy when the altitude of the contour is multiplied by a 'gain' factor (default is 2.0). This transform also ends by idempotence. All the hierarchical levels of image 'imIn' (which is a valued watershed) are computed. 'imOut' contains all these hierarchies which are embedded, so that hierarchy i is simply obtained by a threshold [i+1, 255]


Serge Beucher

Basic Morphological Operators Applied On Partitions

This document describes morphological operators designed for partition images. Two approaches are addressed. In the first one, each cell of the partition is considered independently of the others and various anti-extensive operators can be defined on it. However, these operators can be applied directly on the whole partition. There is no need to extract each cell to transform it. The second approach, on the contrary, considers a partition as a representation of a Delaunay graph. Morphological operators can be defined on this graph. This document explains how these operators can be implemented directly on the partition image. These operators have been implemented with the MAMBA Image library.

Introduction

Partitions (also called partition images or mosaic images) are very common in image analysis and in Mathematical Morphology. They can be found as the result of various morphological operators: watershed transforms, hierarchical segmentations, labellings, residual transforms, especially with the associated functions of a residual operator. The classical way for handling partitions considers that each tile or cell of the partition is a vertex of the valued Delaunay graph associated to this partition. This representation allows to define basic morphological transforms (erosion and dilation) on this neighbourhood graph. In this case, the partition is not modified. Only the valuations of the tiles are changed. These operators have been seldom used in practice. It is surely due to the fact that they seem to require the Delaunay graph associated to the partition to be obtained. However, it is possible to define them directly on the partition image rather quickly. We shall discuss this point in the second part of this paper. The second approach is totally the opposite. Indeed, when a partition is available, it is very interesting to have the capability to apply basic operators to each cell of the partition as if it was alone and not adjacent to the other cells. However, as a partition is often represented as a numerical function where a single grey level labels all the pixels which belong to the same cell, it is not simple to realize efficiently, that is rapidly, such transforms. Moreover, some operators are not always meaningful. It is the case, in particular, when dealing with extensive transforms (as dilations for instance). We shall explain in the first part of this document which operators can be defined and how they can be obtained. These operators have been implemented with the MAMBA IMAGE library [START_REF] Beucher | MAMBA documentation (MambaApi Reference Manual)[END_REF]. The source code is given in the appendix.

Definitions

Let us start with some definitions.

Partition image and background cells

Let be a partition of the space E:

Y k , k c [ 1, K ] ≤Y k ,Y k ∏ , k ! k ∏ : Y k 3 Y k ∏ = - 4 kc[1,K] Y k = E
K is the number of connected components (tiles) of the partition. A 2D partition image is represented by a numerical image f defined as:

x c E : f ( x ) = v k if x c Y k v k is an integer value also called the label of the cell Y k . If Y k and Y k' are adjacent cells, we have . v k ! v k ∏ This definition simply means that all the pixels belonging to the same cell have the same value but different cells may have the same label provided that they are not adjacent. Note that, thanks to this definition (less restrictive than the definition where each cell has a unique label), any greyscale image can be considered as a partition image (Fig. 1). Therefore, as we may have different tiles with the same label i, we will note:

X i = 4 v k =i Y k
the union of all the cells Y k of the partition which share the same label i. We have:

E = 4 ic[0,n] X i
n is the highest label used in the partition. The set X 0 containing all the cells with a label equal to 0 will play a particular role in the following. This set is called the background of the partition. It is perfectly allowed that a partition does not contain any background. This background allows to define and to use extensive operators.

Fig.1:

Examples of partition images. On the left, example of a "classical" partition. However, any grey level image can be considered as a partition (on the right).

Indicator function

Another useful definition in the sequel is the following: Let M be a set included in E. We denote m the valued indicator function of M defined as:

m ( x ) = n if x c M m ( x ) = 0 if not
n is the highest label used in the partition.

Operators applied to the cells of a partition

This section will explain how to define and to realize basic operators applied on each cell of the partition. Due to the fact that all the cells are adjacent, the only operators which should be defined must be anti-extensive otherwise the cells would overlap after transformation. However, we shall see that it is possible to apply extensive operators on partitions provided that a background exists in this partition.

Basic cells erosion

It is obviously possible to compute this operation cell by cell. Each cell is eroded and a new function f' is then built from each eroded cell thus generating a new partition where each previous cell is separated from the others by a background region. However, this procedure would be extremely slow. Let us describe faster ways to obtain this result. We want to define an operator ε P applied on the function f representing the partition Y k which produces a new function f' representing the partition where each cell Y' k is the Y k ∏ erosion with a structuring element B of the initial cell Y k :

x c E :

f ∏ ( x ) = v k if x c Y k 0 - B Note that B must
contain its origin otherwise the erosion would not be anti-extensive (albeit we know that, if this condition is not fulfilled, this case can be considered as a mere translation of the transformed partition). Note also that the initial background cells (the set X 0 ) are also eroded. However, as every point x which does not belong to any eroded cell is given the value 0 in f', the background set X 0 ' of the final partition is larger than the initial set X 0 : X 0 _ X 0

∏

We may also remark that:

x c E : f ∏ ( x ) = i if x c X i 0 -B where X i is the union of all the cells Y k with a label equal to i. Indeed, every cell Y k belonging to X i does not touch the other cells of X i . Therefore, we have:

X i V - B= 4 v k =i Y k 0 - B Let us prove that f' is equal to: f ∏ = P ( f ) = f . m with m, valued indicator function of the set: M = x : f 0 - B ( x ) = f / - B ( x
) M is the mask of the points of E where the (classical) erosion of f by B is equal to the dilation. Denoting Z i (f) the threshold of f at level i, we have:

Z i ( f ) = x c E : f ( x ) m i That is: Z i ( f ) = 4 jmi ( X j )
and:

Z i ( f ∏ ) = 4 jmi X j 0 - B
We have also:

X j = Z j ( f ) # Z j+1 ( f ) = Z j ( f ) 3 Z j+1 c ( f )
Therefore, we can write:

Z i ( f ∏ ) = 4 jmi Z j ( f ) 3 Z j+1 c ( f ) 0 - B Z i ( f ∏ ) = 4 jmi Z j ( f ) 0 - B 3 Z j+1 c ( f ) 0 - B Z i ( f ∏ ) = 4 jmi Z j ( f ) 0 - B 3 Z j+1 ( f ) / - B c
and:

Z i ( f ∏ ) = 4 jmi Z j ( ( f )) 3 Z j+1 c ( ( f ))
where ε and δ are simply another notations for the erosion and the dilation.

The mask M can also be defined as the threshold at 0 of the difference between the erosion and the dilation (as the structuring element B contains its origin, we have always

). ( f ) [ ( f ) M = Z 0 ( ( f ) -( f ))
This threshold can be written (see [START_REF] Beucher | About a problem of definition of the geodesic erosion[END_REF], appendix A for more details):

Z 0 ( ( f ) -( f )) = 4 j Z j ( ( f )) 3 Z j+1 c ( ( f ))
Let us write the threshold of the function f'. We have indeed:

Z ( m ) = M ( i > 0 ) Therefore: f ∏ = f . m u Z i ( f ∏ ) = Z i ( f ) 3 Z i ( M ) = Z i ( f ) 3
M Note that when i is equal to 0, all the thresholds are equal to E. Therefore:

Z i ( f ∏ ) = Z i ( f ) 3 4 j Z j ( ( f )) 3 Z j+1 c ( ( f )) Z i ( f ∏ ) = 4 j Z j ( ( f )) 3 Z j+1 c ( ( f )) 3 Z i ( f )
We can split this expression into two parts. The first one corresponds to the sections higher than or equal to i. The second part corresponds to the sections lower than i.

Z i ( f ∏ ) = 4 jmi Z j ( ( f )) 3 Z j+1 c ( ( f )) 3 Z i ( f ) 4 4 j<i Z j ( ( f )) 3 Z j+1 c ( ( f )) 3 Z i ( f ) Z i ( f ∏ ) = 4 jmi Z j ( ( f )) 3 Z j+1 c ( ( f )) 3 Z i ( f ) 4 4 j<i Z j ( ( f )) 3 Z j+1 c ( ( f )) 3 Z i ( f ) (a) (b)
In term (a), we always have as the erosion is increasing and anti-extensive. Z j ( ( f )) _ Z i ( f ) This term can be reduced to:

( a ) = 4 jmi Z j ( ( f )) 3 Z j+1 c ( ( f ))
Regarding term (b), we always have as the dilation is increasing and

Z j+1 c ( ( f )) 3 Z i ( f ) = - extensive. So (b) is equal to the empty set: ( b ) = - Finally, we obtain: Q.E.D. Z i ( f ∏ ) = 4 jmi Z j ( ( f )) 3 Z j+1 c ( ( f ))
The erosion f' of the cells of the partition image f by a structuring element B containing its origin is therefore obtained by the following sequence of operations:

Compute the morphological gradient of f with B.

Extract the pixels with a gradient equal to 0 (by thresholding). This produces a mask M. Then compute the indicator function m of this set M. We obtain finally . f ∏ = f . m

Cells dilations and openings

Defining the dilation of the cells of a partition is not easy as there is generally no "room" between the different cells to insert the dilations. The only exception concerns partitions with a background. In this case, it is possible to dilate each cell and to keep only parts of the dilated cells which fall into the background. If we denote by b the indicator function of the background X 0 , we could define the dilation δ P as:

P = f / - B . b -f
However, this definition is quite dangerous because, if the structuring element B is large, the dilation of some cells, those with a high valuation, will invade the background region to the detriment of the dilations of the cells with a lower valuation. To cope with this problem, the only allowed structuring element B must be the elementary structuring element (square or hexagon). This is equivalent to use a geodesic dilation. At each step, the elementary operation is performed:

P = f / - B . b ( f ) -f Where b(f)
is the indicator function of the background of f. Then the next step of dilation is achieved by replacing f by δ P and by iterating the operation. This transformation however is not very useful. Therefore it has not be implemented.

There is another case where the dilation of the cells can be performed with any structuring element B. This is allowed in the definition of the opening g P of each cell by B. This operator is simply defined by: P = P / B Indeed, in this case, the dilation by , transposed structuring element, is always contained in -B the non background initial cells. The corresponding MAMBA operation (named cellsOpen) is given in the appendix. Fig. 3 shows an example of opening of a partition.

Fig. 3: Hexagonal opening of size 3 (right) of the cells of the partition (left).

Elementary operators applied on cells

Most of the morphological operators use combinations of basic transforms defined on simple structuring elements composed by doublets of points in a given direction. It is, for instance, the case in the MAMBA library. Considering a structuring element L made of two points, extremities of a vector a, two basic operators are defined:

L [ f ( x )] = inf ( f ( x ) , f ( x + a )) and L [ f ( x )] = sup ( f ( x ) , f ( x + a )
) Combining these two operators (respectively called infNeighbor / inFarNeighbor and supNeighbor / supFarNeighbor in MAMBA) leads to the definition of many morphological transformations. These two operators can also be written:

if L [ f ( x )] = f ( x ) f ( x ) [ f ( x + a ) if not L [ f ( x )] = f ( x + a ) if L [ f ( x )] = f ( x ) f ( x ) m f ( x + a ) if not L [ f ( x )] = f ( x + a )
However, when defining the erosion of the cells of a partition, the basic transformation P which is used can be written:

if P L [ f ( x )] = f ( x ) f ( x ) = f ( x + a ) if not P L [ f ( x )] = 0
The definition of this operator is similar to the two previous ones. Replacing the two inequalities by a strict equality when comparing the values of the two pixels amounts to consider any pixel with a value different from the value of the pixel at the origin of the structuring element as a pixel belonging to the outside of the cell. This is precisely what we want. A dual operator can also be defined. We can write:

P L if P L [ f ( x )] = f ( x ) f ( x ) ! f ( x + a ) if not P L [ f ( x )]
= 0 These two operators have been programmed with MAMBA. They are named equalNeighbor and nonEqualNeighbor and are given in the appendix. These two operators allow to define hit-or-miss and thinnings transforms applied to each non background cell of the partition. Let (T, T') be a doublet of structuring elements, with . We can define:

T 3 T ∏ = - if P T [ f ( x )] = f ( x ) f ( x ) = f ( x ∏ ) , ≤x ∏ c T x if not P T [ f ( x )] = 0 if P T ∏ [ f ( x )] = f ( x ) f ( x ) ! f ( x ∏ ) , ≤x ∏ c T x ∏ if not P T ∏ [ f ( x )] = 0
The Hit-or-Miss transform of the partition P by (T, T') is therefore:

P P T,T ∏ = P T . P T ∏
We see that, in order to get non zero and anti-extensive transforms, the common origin of both structuring elements T and T' must belong to T. A partition thinning is defined by: P P = I -P P ( f ) = f -P ( f ) These operators (cellsHMT, cellsThin and cellsFullThin) have been implemented and are given in the appendix. Fig. 4 gives the result of a complete thinning operator applied on a partition image when an homotopic structuring element is used (D element, see [START_REF] Beucher | Segmentation d'images et morphologie mathématique[END_REF]).

Fig. 4: Complete D thinning operator applied to partition cells. When a cell is not simply connected, this operator cannot reduce it to a single point (homotopy preservation).

Distance function A fast computation of the distance function

of each non background cell can be d P ( f ) achieved easily.

Fig. 5: Distance function of the cells (right) of the partition (left).

We simply need to perform a size one erosion of the partition image. Then, after a thresholding of the eroded partition, we compute the distance function of the binary set. The final distance function is obtained by adding the indicator function of the initial non background cells (Fig. 5). Here is the detailed algorithm:

Compute the size 1 cells erosion, P ( f ) Extract from the eroded partition the non zero pixels: X = x : P ( f ( x )) ! 0 Compute the distance function of X, d ( X ) Define the indicator function k: k ( x ) = 1 if f ( x ) > 0 Compute the cells distance function :

d P d P ( f ) = d ( X ) + k
This operator is named cellsComputeDistance in the appendix.

Cells reconstructions

Cells reconstruction is a very important tool when dealing with partitions. Two main operators are needed. The first one consists in reconstructing among the different cells of a partition those which are marked. The marker set (binary) M allows to select some cells in the initial partition f. The final partition R P (f, M) is composed of these marked cells (their label value is preserved) and of background cells. The second operator aims at replacing the label value of some cells by a new one. The cells to be modified are pointed by a marker, the new value being the grey level of the marker. These two operators are in fact very similar and they are built the same way. We can consider that the first operator is equivalent to the second one where each marker value is identical to the label value of the corresponding marked cell. If each cell of the partition had a unique value (as many labels as cells), these two transformations could be implemented easily through the definition of a lookup table where the initial labels would be replaced by new ones according to the values indicated by the markers. However, with the extended definition of a partition given above, the implementation is more complex. Indeed, using the geodesic reconstruction can be a solution but if, in the one hand, the size of the cells is very small (one pixel wide for instance) and/or, in the other hand, the markers are placed on the boundary of the cells, the geodesic reconstruction will not work because either this reconstruction will invade non marked adjacent cells or some cells will not be reconstructed despite the fact that they are marked. One solution to this issue consists in using directional reconstructions. This operator already exists in the MAMBA library (buildNeighbor). Fig. 6 to 14 illustrate how to use it. Note that, as the process is directional, the illustrations which are monodimensional describe nevertheless the exact algorithm as it uses buildNeighbor which is a true monodimensional operator.

Let us consider a profile in a given direction α of a partition f (Fig. 6). The first step of the algorithm will consist in detecting all the grey levels transitions in a given direction by means of a directional gradient for instance (Fig. 7):

g g = f / - L -f 0 - L
Where is a doublet of adjacent points in direction α. L Fig. 7: Directional gradient of the previous profile.

Fig. 8: Mask function h corresponding to the zeros of the previous gradient image.

Note that the transition is detected whatever its orientation (upward or downward). This can be observed in the previous figure when the thickness of a cell is equal to 1. In this case, two adjacent transitions appear.

Then, a mask function h is defined corresponding to the zero values of the previous one (Fig. 8): if h ( x ) = n g ( x ) = 0 if not. h ( x ) = 0 (Remind that n is the highest label value of the partition).

Fig. 9: Initial profile and marker function m. This function marks the cells which will be rebuilt and gives their new values.

Consider now the marker set M, that is the set which indicates the cells of the partition to be rebuilt and, more precisely, the associated valued function m where each marker is given a grey value corresponding to the new label of every marked cell (Fig. 9).

Fig. 10: Dilation of the previous markers to insure that the directional reconstruction will be initiated.

Fig. 11: Intersection of the previous markers with the mask (reconstruction will occur inside the mask).

It is necessary to perform an initial directional dilation of m in order to be sure that the following directional reconstruction will occur in any case, especially when a marker is placed at the boundary of a cell. The direction of the dilation is the same as the direction of the previous directional gradient. Note that this dilation will never mark unwanted cells for the reason explained below. We obtain a new marker function (Fig. 10). A directional geodesic reconstruction of the function h can then be achieved R ( h, m ∏ . h ) using the function . The direction used is obviously the direction of the previous m ∏ . h gradient (Fig. 12). The intersection with h prevents the possible marking of adjacent cells which have been addressed above.

Fig. 12: Directional reconstruction of the mask.

Note that some markers may have disappeared in . These markers correspond to thin m ∏ . h cells in the current direction α. However, it is not at all a problem as we will see in the sequel. Note also that the purpose of the directional gradient is to provide a gap between adjacent cells in order to prevent the propagation of the reconstruction into non marked cells. This gap also prevents the initial directional dilation of the markers to overlap on the adjacent non marked cells.

Fig. 13: Reintroduction of the markers which are too thin to be taken into account in the reconstruction.

At the end of this step, we just have to reintroduce the markers which have been suppressed during the process (Fig. 13). In fact, these markers have not contributed to the reconstruction in the current direction. This reintroduction is simply achieved by a supremum with the reconstruction: R ( h, m ∏ . h ) -m

The resulting image can be used as a new marker for iterating the entire process in another directions until idempotence (Fig. 14). Finally, by combining the cells erosion and the cells reconstruction, it is very simple to define a reconstruction opening operator which applies simultaneously on each cell of a partition.

This operator can be found in the appendix (it is named cellsOpenByBuild) and Fig. 17 illustrates it.

Operations on the partition graph

In the previous part, we considered the partition image as a patchwork of cells. The same operation is performed on every cell of the partition independently of the adjacent cells. Everything happens as if each cell was extracted, was processed and reinserted again in the partition.

Fig. 18: Representation of a partition as a valued graph.

In this section, on the contrary, a partition image will be considered as a representation of a valued Delaunay graph (Fig. 18). Each cell of the partition corresponds to a valued vertex of the graph, the grey value of the cell being equal to the value of the vertex. Vertices corresponding to adjacent cells are connected by an edge in the graph.

Erosion and dilation of the partition graph Elementary morphological operators can be defined on this graph ([8], [9]).

Let G(V,E) the Delaunay graph corresponding to a partition P. V is the set of vertices (that is the cells of the partition) and E the set of edges connecting the adjacent cells. The elementary dilation of G is a new valued graph G'(V, E) with the same structure (same vertices and ( G ) edges) but where the new value of each vertex i is given by: v i

∏ v i ∏ = sup kcN i ( v k )
N i being the neighborhood of i composed of i and of all its neighboring vertices (Fig. 19).

Similarly, we can define the eroded graph as the graph where the new value of

G ∏∏ = ( G ) v i ∏∏ each vertex i is given by: v i ∏∏ = inf kcN i ( v k )
Dilations and erosions of larger sizes can be defined by iteration: New partition images can be associated to the new graphs. If P is the partition corresponding to G and f, the corresponding image, we have:

n ( G ) = n ) .. ) ( G ) n ( G ) = n ) .. ) ( G )
f ( x ) = v i if x c Y i Y i
being a cell of the partition P. A new partition P' corresponds to the graph . This partition is represented by a G ∏ = ( G ) function f' defined by: f

∏ ( x ) = v i ∏ = sup kcN i ( v k ) if x c Y i

Fig. 20: Each of the two adjacent cells (a) and (b) is adjacent to a maximal cell with the same value v1 (upper left). The corresponding graph is shown in the upper right. After dilation of the graph, the cells (a) and (b) have merged (lower left) and the edge joining these two cells (in blue) has been contracted (lower right).

However, the cells of the partition P' are not identical to the cells Y i of the initial partition Y i ∏ P. Indeed, two or more adjacent cells in the partition P (with different values) may have the same value after dilation or erosion of the graph. So, these initial cells of the partition P will be merged into a single cell in P' as they were adjacent and they now have the same value. Therefore the partition P' is the representation of a new graph G * which is different from the graph . This new graph is obtained by merging some adjacent vertices and by G ∏ = ( G ) connecting the resulting vertex to the remaining ones (Fig. 20). This operation is called edge contraction.

Defining graph operators on partitions

Performing morphological transforms on graphs is an interesting approach (albeit not really new as any operator applied on a digital image is in fact realized on a neighborhood graph). However, these transformations are seldom used (except in hierarchical segmentations). This is surely due to the fact that the step consisting in generating a graph from a partition image followed by the reverse step to build the resulting partition image is a tedious and time-consuming task, especially if we want to go back and forth between the image and the graph representations. Indeed, using a graph representation allows to deeply reduce the amount of data to be processed, but this also prevents to see the result of the operations, which is a very annoying drawback.

Fig. 21: Graph dilation (upper right) and erosion (lower right) performed directly on the partition image (left).

Fortunately, it is possible to define graph operators by acting directly on the partition image without the need to explicitely build the corresponding graph. Let P be the partition image associated to the graph G. We can obtain directly the partition P' associated to the graph by the following operation: G ∏ = ( G ) Dilate the function f representing the partition P. Perform a reconstruction of the partition P with as marker: ( f )

P t f t R P ( f, ( f )) = f ∏ t P ∏
The dilation of the partition image f allows to mark each cell with the supremum value of all its neighbors. This partial marking is then propagated into the entire cell by means of the reconstruction. In the same way, the partition P'' associated to the graph is obtained by a simple G ∏∏ = ( G ) inversion of the function f at start followed by another inversion of the result (Fig. 21).

Iterating these operations, that is performing a size n dilation or erosion on the graph G, requires that the initial partition image f be kept in order to use it at each step of the operation. The algorithm is therefore the following (for dilation):

Let f' = f (initialisation) Perform the operation n times to get a size n dilation (in f'). R P ( f, ( f ∏ )) = f ∏ If, on the contrary, we simply iterate the transformation on the partition produced by the previous step, the edge contraction associated to the operation induces a faster propagation of the extremal cells values (maxima in case of dilation, minima for the erosion), as shown at Fig. 22.

Fig. 22: Difference between a classical size 2 graph dilation (left) and an iteration of two dilations performed on the partition (right).

These two basic operators can be found in the appendix. They are named partitionErode and partitionDilate.

Conclusion

The operators which have been presented here have been originally defined and designed in order to be able to apply transformations and to collect information (size, shape, topology, etc.) on the cells of a partition image as if each cell was processed individually. However, all the cells are transformed at the same time, which increases dramatically the processing speed by avoiding the excruciatingly slow individual analysis approach.

To achieve this, we showed that the definition of two new basic transformations applied on doublets of points ( and ) allows to define on the cells of a partition P a whole range of P L P L anti-extensive operators: erosion, opening, Hit-or-Miss transform, thinnings. Geodesic reconstruction algorithms applied to the cells of a partition have also been designed. These algorithms widely use the recursive directional reconstruction implementation which exists in the Mamba Image library. Although this implementation is slower than approaches based on hierachical queues, it is currently the only one which allows to cope with propagation issues which appear when the thickness and area of the cells are very small. However, this reconstruction remains very fast.

Remind also that these operators can be applied on any digital image, thanks to the extended definition of a partition image which was introduced and used for their design. So, these operators are unquestionably useful when applied on partition images produced by many morphological transformations: watershed transform (especially the one which does not generate boundaries between the different catchment basins) [START_REF] Beucher | Hierarchical Queues: general description and implementation in MAMBA Image library[END_REF], hierarchical segmentations (waterfalls, P algorithms, etc.) [START_REF] Beucher | Towards a unification of waterfalls, standard and P algorithms[END_REF]

Appendix

The various partition operators presented above have been implemented with the Mamba Image library. The Python source code is given below. These operators have been gathered in a Pyhon module named partitions.py. This module is available at http://cmm.ensmp.fr/~beucher/Mamba/partitions.py.

""" This module contains operators acting on partitions. Two types of operators are defined: operators acting on each cell of the partition independently, or operators considering each cell of the partition as a node of a weighted graph.

In the first case, each cell of the partition is considered as a binary set and the defined operation is applied on each cell to produce a new transformation (which is not necessarily a partition). In the second case, the partition is considered as a graph and morphological operators are defined on this graph. These operators provide morphological transformations on graphs, without the need to explicitly define this graph structure from the partition. """ This operator compares the value of each pixel of image 'imIn' with the value of its neighbor in direction given by 'nb'. If the values are equal, the pixel is unchanged. Otherwise, it takes value 0. This operator works on hexagonal or square 'grid' and 'edge' is set to FILLED by default. This operator works for 8-bit and 32-bit images. """ imWrk1 = mamba.imageMb(imIn) imWrk2 = mamba.imageMb(imIn, 1) mamba.copy(imIn, imWrk1) mamba.copy(imIn, imOut) mamba.supNeighbor(imIn, imWrk1, nb, 1, grid=grid, edge=edge) mamba.infNeighbor(imOut, imOut, nb, 1, grid=grid, edge=edge) mamba.generateSupMask(imOut, imWrk1, imWrk2, False) mamba.convertByMask(imWrk2, imWrk1, 0, mamba.computeMaxRange(imIn) [START_REF] Beucher | MAMBA documentation (MambaApi Reference Manual)[END_REF]) mamba.logic(imOut, imWrk1, imOut, "inf") def nonEqualNeighbor(imIn, imOut, nb, grid=mamba.DEFAULT_GRID, edge=mamba.FILLED): """

Fig. 2

 2 Fig.2 illustrates this operator. The corresponding MAMBA operation (named cellErode) is given in the appendix.

Fig. 2 :

 2 Fig.2: Cells erosion (size 1, hexagonal structuring element) of the previous partition images (see Fig. 1)..

Fig. 6 :

 6 Fig. 6: Profile of a partition image f in a given direction.

  m ∏ = m / -L Computing the function allows to keep only parts of m' which are less than h (Fig. ( m ∏ . h ) 11).

Fig. 14 :Fig. 15 :

 1415 Fig. 14: Final reconstruction when all the directions are used iteratively.

Fig. 15

 15 Fig. 15 illustrates this transformation applied to extract marked cells of a partition image.

Fig. 16 :

 16 Fig. 16: The upper right partition image is used as a marker to rebuild and to value the cells of the upper left partition. The result of the operation is the lower image.

Fig. 16

 16 Fig. 16 shows a cells reconstruction of a partition image which uses as markers the cells of another partition image.

Fig. 17 :

 17 Fig. 17: Opening by reconstruction of size 8 (right image) of the cells of the left partition image.

Fig. 19 :

 19 Fig. 19: Neighborhood of a vertex v0 in a graph which is considered for the definition of the dilation and the erosion.

  , [7], residual operators (in particular ultimate openings and granulometric functions)[START_REF] Beucher | Numerical residues[END_REF]. These operators bring efficient and fast tools for cells filtering, extraction and measurement. It would also be possible to use them to address some more complex issues: comparison and mixing of partitions, watermarking algorithms, markers extraction and tracking, features extraction... Finally, we have revisited the basic morphological transformations on graphs. As, in most cases in image analysis, graphs are built from partitions, it is interesting to be able to perform graph operations while staying in the partition image representation of the graph. Indeed graph dilation and erosion can be obtained directly on the partition associated to the graph by means of the cell reconstruction defined in the first part of this document.

(available at http://cmm.ensmp.fr/~beucher/publi/Unified_Segmentation.pdf).

[7] Beucher S., Marcotegui B.: P algorithm, a dramatic enhancement of the waterfall transformation..CMM/Mines Paristech publication, september 2009. (available at http://cmm.ensmp.fr/~beucher/publi/P-Algorithm_SB_BM.pdf).

[8] Burger F., Gandillot X., Treillard P.: Transformations morphologique sur un graphe, application à l'étude des lymphocytes. Internship report, Paris School of Mines / Ecole Polytechnique, 1982.

[9] Vincent L.: Graphs and Mathematical Morphology. Signal Processing, Vol. 16, n°4, april 1989, pp. 365-388.

This document is copyrighted under the Creative Commons "Attribution Non-Commercial

No Derivatives"license. For terms of use, see http://creativecommons.org/about/licenses/.

  .erode(imIn, imOut, n=n, se=se, edge=edge) mamba.generateSupMask(imOut, imWrk1, imWrk2, False) mamba.convertByMask(imWrk2, imWrk1, 0, mamba.computeMaxRange(imIn)[START_REF] Beucher | MAMBA documentation (MambaApi Reference Manual)[END_REF]) mamba.logic(imOut, imWrk1, imOut, "inf") def cellsOpen(imIn, imOut, n=1, se=mC.DEFAULT_SE, edge=mamba.FILLED):""" Simultaneous opening of size 'n' (default 1) of all the cells of the partition image 'imIn' with 'se' structuring element. The resulting partition is put in 'imOut'. 'edge' is set to FILLED by default. This operation works on 8-bit and 32-bit partitions. """ imWrk = mamba.imageMb(imIn) cellsErode(imIn, imWrk, n, se=se, edge=edge) mC.dilate(imWrk, imOut, n, se=se.transpose()) def cellsComputeDistance(imIn, imOut, grid=mamba.DEFAULT_GRID, edge=mamba.EMPTY): """ Computation of the distance function for each cell of the partition image 'imIn'. The result is put in the 32-bit image 'imOut'. This operator works on hexagonal or square 'grid' and 'edge' is set to EMPTY by default.

	"""
	imWrk1 = mamba.imageMb(imIn)
	imWrk2 = mamba.imageMb(imIn, 1)
	se = mC.structuringElement(mamba.getDirections(grid), grid)
	cellsErode(imIn, imWrk1, 1, se=se, edge=edge)
	mamba.threshold(imWrk1, imWrk2, 1, 255)
	mamba.computeDistance(imWrk2, imOut, grid=grid, edge=edge)
	mamba.addConst(imOut, 1, imOut)
	""
	# Contributor: Serge BEUCHER
	import mamba
	import mambaComposed as mC
	def cellsErode(imIn, imOut, n=1, se=mC.DEFAULT_SE, edge=mamba.FILLED):
	"""
	Simultaneous erosion of size 'n' (default 1) of all the cells of the partition
	image 'imIn' with 'se' structuring element. The resulting partition is put in
	'imOut'.
	'edge' is set to FILLED by default.
	This operation works on 8-bit and 32-bit partitions.
	"""
	imWrk1 = mamba.imageMb(imIn)
	imWrk2 = mamba.imageMb(imIn, 1)
	mC.dilate(imIn, imWrk1, n=n, se=se)

mCdef equalNeighbor(imIn, imOut, nb, grid=mamba.DEFAULT_GRID, edge=mamba.FILLED):

This operator compares the value of each pixel of image 'imIn' with the value of its neighbor in direction given by 'nb'. If the values are different, the pixel is unchanged. Otherwise, it takes value 0. This operator works on hexagonal or square 'grid' and 'edge' is set to FILLED by default. This operator works for 8-bit and 32-bit images. """ imWrk1 = mamba.imageMb(imIn) imWrk2 = mamba.imageMb(imIn) imWrk3 = mamba.imageMb(imIn, 1) mamba.copy(imIn, imWrk1) mamba.copy(imIn, imWrk2) mamba.supNeighbor(imWrk1, imWrk1, nb, 1, grid=grid, edge=edge) mamba.infNeighbor(imWrk2, imWrk2, nb, 1, grid=grid, edge=edge) mamba.generateSupMask(imWrk1, imWrk2, imWrk3, True) mamba.convertByMask(imWrk3, imWrk1, 0, mamba.computeMaxRange(imIn) [START_REF] Beucher | MAMBA documentation (MambaApi Reference Manual)[END_REF]) mamba.logic(imIn, imWrk1, imOut, "inf") def cellsHMT(imIn, imOut, dse, edge=mamba.EMPTY):

""" A Hit-Or-Miss transform is performed on each cell of the partition 'imIn'. '