STRUCTURALLY ADAPTIVE MATHEMATICAL MORPHOLOGY BASED ON NONLINEAR SCALE-SPACE DECOMPOSITIONS

Standard formulation of morphological operators is translation invariant in the space and in the intensity: the same processing is considered for each point of the image. A current challenging topic in mathematical morphology is the construction of adaptive operators. In previous works, the adaptive operators are based either on spatially variable neighbourhoods according to the local regularity, or on size variable neighbourhoods according to the local intensity. This paper introduces a new framework: the structurally adaptive mathematical morphology. More precisely, the rationale behind the present approach is to work on a nonlinear multi-scale image decomposition, and then to adapt intrinsically the size of the operator to the local scale of the structures. The properties of the derived operators are investigated and their practical performances are compared with respect to standard morphological operators using natural image examples.


INTRODUCTION
Mathematical morphology is a well-known nonlinear image processing methodology based on the application of lattice theory to spatial structures (Serra, 1988;Heijmans and Ronse, 1990).Let E be a subset of the Euclidean R d or the discrete space Z d , considered as the support space of the image, and let T be a set of grey-levels, corresponding to the space of values of the image.It is assumed that T = R = R ∪ {−∞, +∞}.A grey-level image is represented by a function, i.e., f ∈ F (E, T ) maps each pixel x ∈ E into a grey-level value t ∈ T : t = f (x).The two basic morphological mappings F (E, T ) → F (E, T ) are the grey-level dilation and the grey-level erosion given respectively by where f ∈ F (E, T ) is the original grey-level image and b ∈ F (E, T ) is the fixed structuring function.
The further convention for ambiguous expression is considered: Particularly interesting in theory and in practical applications (Soille, 1999), the flat grey-level dilation and erosion is obtained when the structuring function is flat and becomes a structuring element.More precisely, a flat structuring function of support subspace B is defined as where B is a Boolean set, i.e., B ⊆ E or B ∈ P(E), which defines the "shape" of the structuring element.We notice that B c denotes the complement set of B (i.e., B ∩ B c = / 0 and B ∪ B c = E).The structuring element is defined at the origin o ∈ E, then to each point p of E corresponds the translation mapping o to p, and this translation maps B onto B p , i.e., B p = {b + p : b ∈ B}.Therefore, the flat dilation and erosion of a grey-level image f (x) with respect to the structuring element B are respectively where B is the reflection of B with respect to the origin, i.e., B = {−b | b ∈ B}.
Standard formulation of morphological operators is translation invariant in the space ("horizontal" direction invariance) and in the grey-level intensity ("vertical" direction invariance), see properties and proofs in (Heijmans and Ronse, 1990), i.e., for a given translated image g(z) = f (x + y) + r, it is obtained that δ b (g)(z) = δ b ( f )(x + y) + r.But of course, this result is true if and only if the same structuring function b(x) is considered for each point x of the image.A current challenging topic in mathematical morphology is the construction of adaptive operators; or in other words, operators whose structuring functions become dependent on position or on the input image itself.In previous works, the proposed adaptive operators have been based on two main approaches.
-On the one hand, a variability on the support space E: spatially variable shape of structuring functions according to i) the geometric position in the image: perspective-adapted Morphology by Beucher et al. (1987) and by Cuisenaire (2006), ii) the local regularity of image values: morphological amoebas by Lerallut et al. (2005), intrinsic structuring elements by Debayle and Pinoli (2005) morphological bilateral filtering by Angulo (2011), iii) the orientation: curvilinear morphohessian filter by Tankyevych et al. (2008), locally-variant anisotropic morphological filters by Verdu-Monedero et al. (2011).
-On the other hand, a variability on the value space T : variable size of structuring functions according to the local intensity or contrast: viscous watershed by Vachier and Meyer (2005) and viscous openings/closings by Vachier and Meyer (2007), together with a PDE formulation by Maragos and Vachier (2008).

Motivation of the new adaptive framework
The intensity-adaptive operators, or levelvarying operators, studied in Vachier and Meyer (2007) and Maragos and Vachier (2008) are typically implemented using a level-set decomposition, followed by the processing of each level set by a flat operator of a size which depends on the index of the particular level set and finally the image reconstitution from the processed level sets.Hence, the adaptavility for each point depends on the absolute intensity value.Two points of the image associated to structures of similar scale but of different intensity are processed differently.This work introduces a new framework of nonlinear adaptive operators: the structurally adaptive mathematical morphology.More precisely, the rationale behind the present approach is to work on a nonlinear multi-scale image decomposition based on a morphological family of operators called levellings (Meyer, 1998;Meyer and Maragos, 2000).Then, the size of the structuring function b(x) at point x is intrinsically adapted to the local scale of the structures.In other words, two regions of the image appearing in a similar scale will undergo a nonlinear processing of similar "size" independently of the absolute intensities.The value of the processed image will be obtained by the reconstitution of the processed scales.We notice also that our approach is not limited to flat morphological operators.The properties of the derived operators are investigated and their practical performances are compared with respect to standard morphological operators using natural image examples.The paper is an extended and improved version of the conference paper (Angulo and Velasco-Forero, 2010).

Paper organisation
The rest of the paper is organised as follows.In the first part of Methods Section is presented the structural image model, using a pyramid of levellings, which is the basic ingredient for the present adaptability framework.Then, the second part of Methods Section introduces the structurally adaptive pseudo-dilation and pseudo-erosion and other associated operators, including also the study of the algebraic properties of these new operators.Results Section discusses various examples which illustrate the behavior of the structurally adaptive operators for typical applications of mathematical morphology.Finally, Section on conclusions and perspectives closes the paper.

METHODS STRUCTURAL IMAGE MODEL USING A PYRAMID OF LEVELLINGS
The classical tool from mathematical morphology for multi-scale image decomposition is based on the granulometric analysis (Serra, 1982;Maragos, 1989).A granulometry is the study of the size distribution of the objects of an image (Serra, 1982).Formally, for the discrete case, a granulometry is a family of openings Γ = {γ n } n≥0 that depends on a positive parameter n (which expresses a size factor) such as: i) γ 0 ( f T ); and iv) γ n verifies the granulometric semi- group absorption law; i.e., Moreover, a granulometry by closings (or antigranulometry) can be defined as a family of increasing closings Φ = {ϕ n } n≥0 .In practice, the most useful granulometry and anti-granulometry are those associated to morphological flat openings and closings: respectively, where B 1 = B is a isotropic structuring element of unit size (unitary discrete disk), and B n is the homothetic structuring element of size n, with n = 1, 2, • • • .The granulometric analysis of an image f with respect to a finite Γ consists in computing the pyramid of opened images and then in calculating the image residue between two successive openings: The pyramid of residues of openings {ρ leads to a decomposition of bright structures into N scales.We recall that the subgraph of the opened image by B i is equivalent to the union of the translations of the structuring when it fits the subgraph of the original image (Serra, 1982), or in other terms, image structures in the foreground (i.e., bright ones) that cannot contain B i are removed by the opening.Therefore, the residue ρ + i extracts the structures removed between B i−1 and B i .By duality, the image decomposition into N scales of dark structures is obtained from the pyramid of residues of closings {ρ In summary, using morphological openings and closings, it is obtained a pair of nonlinear scale-space image decomposition, i.e., Furthermore, morphological openings/closings deal selectively with the grey-level structures according to their support size but they also "distort" the contours (level sets) of the image objects.In mathematical morphology, the powerful geodesic operators (Serra, 1988;Soille, 1999) (or connected filters) are perfectly able to deal with the objects, without distorting their contours.
For all these reasons, we have considered for the present structurally adaptive operators a multiscale framework which, on the one hand, decomposes simultaneously bright and dark structures (by means of alternate sequential filters); and on the other hand, uses geodesic operators (in particular, the levelling operator) to extract precisely the grey-level objects.Let us consider in detail this methodology.

Nonlinear scale-space based on viscous levellings
The levelling (Meyer, 1998) of the reference image f (x) according to the marker image m(x) (this latter is generally a rough simplification of the reference image), denoted λ ( f , m)(x), is a morphological geodesic filter which simplifies textures and eliminates small details of the reference image according to the marked structures, but preserving the contours of remaining objects.In practice, the levelling is obtained by iteration of geodesic dilations and erosions until idempotence.More precisely, a straightforward algorithm is as follows (Meyer, 1998): such that the idempotence has been reached, i.e., where the i-iteration is Levellings of different nature are obtained according to the pair of transformations (α, β ) used in Eq. 8 and 9.More precisely, two of the most interesting ones are (Meyer, 1998) -standard levelling: Fig. 1 shows a comparison of standard levelling vs. viscous levelling using three different marker images.We observe that even the structures which have been removed from the marker can be reconstructed (by strict geodesic propagation) by standard levelling.
On the contrary, viscous levelling produces a reconstruction more faithful to the simplification degree driven by the marker.Consequently, by its skilful properties on structure regularisation, we consider here the viscous levelling as the fundamental operator for our structural decomposition.
(a) f A pyramid of levellings is then obtained by associating to the image a pyramid of markers, i.e., which leads to a nonlinear scale-space (Meyer and Maragos, 2000).The marker at level i determines the degree of simplification of the corresponding levelled image; the simplification effect increases with i.
The levelling can be considered as the activity supremum between an opening by reconstruction and a closing by reconstruction (Meyer, 1998), leading to a brigh/dark symmetric operator.Let ∁ f stand for the image complement (or negative We have for the levelling the following property: Hence, to obtain a self-dual and symmetric processing of brigh/dark structures, a self-dual family of marker images is required.Classical examples of self-dual filters are the convolution with Gaussian kernels or the (weighted) median filters.In particular, the nonlinear scale spaces associated to levellings where the markers are Gaussian filters have been studied by Meyer and Maragos (2000).In addition, this kind of marker was also used in (Sofou et al., 2005) to extract the texture layer of an image as the residue of the corresponding levelling.
Nevertheless, we propose in this study to use as multi-scale markers the averaged alternate sequential filters (AASF): where the pair of alternate sequential filters (ASF) (Serra, 1988) are given by the following products of openings and closings The properties of ASF are well known in mathematical morphology (Serra, 1988).In particular, they are idempotent and they follow the semi-group absorption law, i.e., ∀n ≥ m ≥ 0, we have ,II  n=max(n,m) , and consequently this property is inherited by AASF B n .ASF are very useful for noise cancellation and smoothing since they produces progressive simplification of image structures.This behaviour is related with their property of compatibility with scale modifications (i.e., intuitively that means that the filter of order k works at a scale k as a filter of order 1 at a scale 1).Theoretically, ASF are not self-dual, but our symmetric averaged version AASF B n leads to an quasi self-dual operator: we have observed that in practice the difference between AASF B n (∁ f )(x) and ∁AASF B n ( f )(x) is always negligible.The main drawback of AASF B n is the computation time, which increases with the size n of the scale.However, using well-known optimized implementations of discrete dilation/erosion for the openings/closings, the final computation time is not nowadays a limiting point.The rationale behind this choice of AASF B n , instead of a Gaussian filter of equivalent scale n, is because we have observed empirically that the viscous levellings associated with alternate sequential filter are more effective for extracting the bright/dark structures of corresponding scale n.

Structural image model
We have now all the ingredients to introduce the following structural model image: where c is the continuous component, s(x) is the structure component, and {t i (x)} N i=1 are the texture components at N scales.
More precisely, given the levelling scale-space {λ ( f , m i )(x)} N i=1 , the components of the present image model can be calculated as follows.
-The continuous component is obtained as the mean value of the last levelled image, i.e., We can consider that c varies with the DC level of the image acquisition device.In the operators presented in this paper c does not give any useful information but is still necessary for ongoing research on adaptive self-dual operators (see perspectives section); -The structure component corresponds to the meancentred last levelled image, i.e., Therefore, the image s(x) represents the most simplified image according to the choice of scale N. In other words, the level-sets of the significant objects having a "support size" bigger than the structuring element B N .
-The texture components are obtained as the derivative of the levelling pyramid, i.e., where m 0 (x) = f (x).Each residue image t i represents the brigh/dark objects which were in scale image i − 1 but have been removed in scale image i; that is, level-sets of objects whose "support size" is comprised between both successive scales.The intensity at each point of t i is the relative intensity of the object (which is invariant to monotone increasing intensity transformation on original image).In fact, the residue image between two successive levellings can be considered as a "morphological laplacian pyramid".In natural images, we can suppose that the texture components are sparse images (i.e., only a limited number of pixels are different of 0) and that the complexity of the original image is split into a simpler image (structure component) and a scale-complexity bounded series of images (texture components).
Our additive image decomposition model can be considered as generalisation of the Cartoon + Texture decomposition by Meyer (2002): f (x) = u(x) + v(x) where u(x) is the cartoon component (homogeneous zones of the objects) and v(x) is the texture oscillation, which is usually solved using variational algorithms (Aujol et al., 2006).The terms can be identified as u(x) = c + s(x) and v(x) = ∑ N i=1 t i (x).The value N of the larger structuring element (scale of structure equivalent) in our model is equivalent to the "scale size" in other Cartoon + Texture decomposition algorithms (Buades et al., 2010).(Buades et al., 2009).
The original image is given on Fig. 1.
Fig. 2 gives an example of three levels of a viscous levelling pyramid and the corresponding residues images, used to define the texture components.We notice that the texture components, with respect to the structure component, are not defined in terms of periodicity or regularity as in the nowadays extended Cartoon+Texture decomposition.That can be observed in the comparative example of our structural decomposition (size of largest structure scale N = 8) and the Cartoon+Texture decomposition recently proposed by Buades et al. (2010) (scale size σ = 3) given in Fig. 3.The Cartoon + Texture decomposition has been obtained using the algorithm available online in (Buades et al., 2009).Therefore, the notion of texture matches here exclusively with the notion of scale.We remark also that, in comparison with a laplacian pyramid associated to a Gaussian scalespace, the main advantages of the derivative levelling pyramid are, on the one hand, the preservation of the contours of extracted structures and on the other hand, at each scale, the size of the structures is bounded by the size of the AASF filter used as marker.

STRUCTURALLY ADAPTIVE OPERATORS
The structural image model introduced above is the key element for the structurally adaptive operators studied in this section.As we will show below, the two basic adaptive operators present a lack of some algebraic properties; and consequently they are not dilations and erosions in stricto sensu.Therefore, we prefer to them structurally adaptive pseudo-dilation and pseudo-erosion.

Adaptive pseudo-dilation and pseudo-erosion
Given the image f ∈ F (E, T ), structurally decomposed into the scales {m, m + k, m + 2k, • • • , M} (denoted compactly by {m : k : M}) according to the model in Eq. 13, we define the corresponding structural adaptive dilation as and the structural adaptive erosion as where δ b n (g)(x) and ε b n (g)(x) are the standard dilation and erosion of image g(x) according to the fixed spatially invariant structuring function b n (x).
We have considered in this study examples using the two most useful families of isotropic multi-scale structuring functions: -parabolic shape function of width n, i.e., b n (x) = − x 2 /2n; -flat disk of radius n, denoted by B n .
For the examples given in the paper, the family of disks B n used in the flat operators are implemented digitally by hexagons of radius n pixels.This choice is based on the properties of isotropy of hexagonal grid with respect to the square one and the fact that the hexagons are an approximation to the disk which can be implemented using optimized algorithms.In any case, if the computation time is not critical, the flat disks can be obtained by thresholding the paraboloidal structuring functions.Fig. 4 depicts a comparison of three flat structural adaptive pseudo-dilations for M = 4, 8 and 12 (with m = 2 and k = 2), versus the corresponding standard flat dilations.We notice that the "size" of the operator is determined by the size of the biggest scale M; or in other words, the size of the scale of the structure component.The two other parameters from the structural decomposition: i) the size of the smallest scale m, and ii) the sampling parameter k (or step size between two successive scales) have a lower impact on the adaptive operator, see for instance the example given in Fig. 5.This effect allows to introduce a sampling parameter k > 1 in order to reduce the number of processed discrete scales (reduction of computation time).

Properties of structurally adaptive operators
Let us consider the algebraic properties of these two morphological operators.The proofs are not included but are easily obtained from the standard properties of erosion and dilation.
-The structural adaptive pseudo-dilation (pseudoerosion) is extensive (anti-extensive), i.e., ∀ f ∈ F (E, T ), -The structural adaptive pseudo-dilation and pseudo-erosion are dual operators with respect to the grey-level inversion, i.e., ∀ f ∈ F (E, T ), -However, the fundamental law which links the pair of dilation/erosion, the adjunction property, fails for the introduced structurally adaptive operators, i.e., ∀ f , g ∈ F (E, T ), As recently pointed out by Roerdink (2009), this drawback is common to any adaptive morphology dilation/erosion where the processing at each point depends on the input value.
-Moreover, as a consequence of the lack of adjunction property, other properties are also lost.Typically, the structural adaptive pseudo-dilation (pseudo-erosion) does not commute with the supremum (infimum), i.e., ∀ f , g ∈ F (E, T ), -Another fundamental property of standard morphology which fails is the increasiness.The structural adaptive pseudo-dilation/pseudo-erosion are generally not increasing operators, i.e., f ≤ g does not always involve that δ m:k:M ( f ) ≤ δ m:k:M (g).In fact, the property is verified if or in other terms, the ordering must be preserved in the derivative pyramid, which is not always the case.
In mathematical morphology, a pair of dual operators which are extensive and anti-extensive are respectively a thickening and thinning (Serra, 1982;Soille, 1999); therefore, the structurally adaptive operators are just algebraic grey-level thickenings/thinnings.In any case, despite the lack of some properties, the operators δ m:k:M ( f ) and ε m:k:M ( f ) can be used to construct structurally adaptive gradient, laplacian, toggle mapping, etc.

Adaptive pseudo-opening and pseudo-closing
Due to the fact that the introduced structurally adaptive operators are not formally a pair of adjunct dilation/erosion, their product δ m:k:M ε m:k:M ( f ) is not an algebraic opening.Hence, we prefer to formulate the structurally adaptive pseudo-opening as and dually, the structurally adaptive pseudo-closing is defined by The structural adaptive opening in Eq. 19 (closing in Eq. 20) is anti-extensive (extensive) but it is not increasing.The idempotency, γγ( f ) = γ( f ), is the property which guarantees the stability of standard opening/closing.Even if our adaptive operators are not strictly idempotent, their iteration is almost stable when AASF filters are used for the levelling decomposition.As we can observe in example given in Fig. 6, the adaptive opening γ b 2:2:8 removes bright objects smaller than the size of the structuring function but, with respect to the standard opening γ b 8 , the notion of "smaller than" depends on the scale of decomposition of the object.The product of operators γ m:k:M ( f ) and ϕ m:k:M ( f ) can be then used to define other more complex operators.

APPLICATIONS AND RESULTS
Let us illustrate the performance of the present operators by means of four typical applications of mathematical morphology.The first case study, depicted in Fig. 7, corresponds to a coronary network acquired by contrast-enhanced radiography.Besides the vessels, the image presents a strong and irregular background (which includes other structures such as the ribs).In order to extract the vessels, the classical top-hat (Serra, 1982) is particularly useful: i.e., residue between the original image f and the opening γ b m:k:M ( f ).The size of the isotropic opening should be bigger than the diameter of the thickest vessels.For the particular example, we consider a parabolic structuring function of width n = 10.As we can observe, the standard opening produces a rough regularisation by removing all the peaks smaller than b 10 whereas the structurally adaptive opening removes also the peaks at different scales but preserves better the secondary structures; hence, the comparison of the standard th + b 10 ( f ) against the adaptive top-hat th + b 2:2:10 ( f ) shows that the second one presents a better extraction of thin vessels as well as less texture and irregularities.
The second application, given in Fig. 8, deals with an example of image denoising and detail simplification.The mosaic image is quite textured and noisy and the aim is to filter it out in order to obtain a more regular image where the main contours are well preserved.This goal can be achieved by means of the morphological centre operator (Serra, 1988), defined in terms of two products of openings and closings of a chosen size: In order to compare quantitatively the behavior of standard vs. adaptive morphological centre, we have evaluated in an experiment the quantitative assessment of noisy images.Fig. 9 summarizes the results of this experiment: starting from an image corrupted with four different levels of Gaussian noise and with four different levels of impulsive noise ("salt and pepper" noise) , we have filtered out the eight images using standard parabolic centre ζ b 4 and structural adaptive parabolic centre ζ b 1:1:4 , and then, we have computed the PSNR value with respect to the original image.For all the noise levels, the value of PSNR is better for the structurally adaptive case than for the standard one.We note also that the "size" of the filter 4 is not particularly appropriate for low Gaussian noise levels or for high impulsive noise levels, but even in such a case the trade-off between denoising and image preservation is better for the adaptive version of the morphological center.
A simple but efficient image contrast enhancement is obtained by subtracting the laplacian image from the original one: using a full scale decomposition, with m = 1 and k = 1.
The last case study deals with an edge detector which can be implemented using only a pair of morphological erosion and dilation, and which is useful for complex images, like the one given in Fig. 11.In fact, the zero crossings of the laplacian l p b m:k:M ( f ), as in the classical Marr-Hildreth model, correspond to the edges of image f .In order to select only the most prominent edges, the zero crossing detector is "multiplied" by the binary image obtained with a threshold by hysteresis of the morphological gradient (Serra, 1982;Soille, 1999)  (26) Besides the size parameter of the laplacian and the gradient, the single parameter of this edge detection model is the threshold value for the gradient.We observe again in this example that the structurally adaptive gradient ξ b 1:1:4 ( f ), using a parabolic structuring function, is more appropriate than the standard counterpart.The comparison by simple visual assessment of edge detection is not conclusive, but it seems that in the adaptive case the edges of small objects are better detected.

CONCLUSIONS AND PERSPECTIVES
We have discussed an additive image model associated to a nonlinear multi-scale image decomposition using a family of viscous levellings.Then, working on this decomposition, we have introduced structurally adaptive morphological operators, where each component of the image is processed with a structuring function of size intrinsically adapted to the corresponding local scale.
We have shown that, in practical applications, the proposed operators perform better than standard ones for object extraction, image denoising, image enhancement, etc.This positive behaviour is justified by the fact that, using geodesic operators for image decomposition, the different objects are almost individually processed and the "interferences" between adjacent objects are notably reduced.
In addition, we have proved that another advantage of the adaptability is the fact that the choice of the "size" for the erosions/dilations and the openings/closings is much less critical.
However, we have also demonstrated that the underlying algebraic structure of the structurally adaptive operators is less rich than the standard ones.As we have discussed, in the present image model, the structure and texture terms are positive/negative signals and consequently, in ongoing work, we will formulate adaptive self-dual erosion and opening.Finally, sparsity properties of our image model suggest us to consider it as a starting point to explore the notion of sparse mathematical morphology in future research.
and that f (x + h) − b(h) = +∞ when f (x + h) = +∞ or b(h) = −∞.The other morphological operators, such as the opening and the closing, are obtained as products of dilation/erosion.For instance, starting from the adjunction pair {δ b , ε b }, the opening and closing of a grey-level image f according to the structuring function b are the mappings F (E, T ) → F (E, T ) given respectively by

Figure 2 .
Figure 2. Example of laplacian pyramid on viscous levelling scale-space (original image on Fig. 1): first row (a), pyramid of viscous levellings; second row (b), derivative of the viscous levelling pyramid.

Figure 3 .
Figure 3.Comparison of our structural decomposition (size of largest structure scale N = 8) and the Cartoon + Texture decomposition (Buades et al., 2010) (scale size σ = 3): (a) structure component, (b) sum of texture components, (c) cartoon image, (d) texture image.The Cartoon + Texture decomposition corresponds to the algorithm available online in(Buades et al., 2009).The original image is given on Fig.1.

Figure 4 .
Figure 4. Comparison of flat structural adaptive pseudo-dilations (second row, b) vs. the corresponding standard flat dilation (first row, a).

Figure 6 .
Figure 6.Comparison of parabolic structural adaptive pseudo-openings (second row, b) vs. the corresponding standard parabolic openings (first row, a).
ζ − b m:k:M ( f ) = γ b m:k:M ϕ b m:k:M γ b m:k:M ( f ) .(23) We compare in particular the denoising effect using a flat structuring element of size N = 4 and a parabolic structuring function of similar width.As we can observe in the images, the adaptive flat or parabolic versions of the operator, ζ B 1:1:4 ( f ) and ζ b 1:1:4 ( f ) outperforms clearly the effects of the standard ones ζ B 4 ( f ) and ζ b 4 ( f ); yielding a more efficient denoising/simplification image, especially the parabolic case, with an excellent contour preservation of remaining structures.