MULTISCALE IMAGE ANALYSIS BASED ON ROBUST AND ADAPTIVE MORPHOLOGICAL SCALE-SPACES

Mathematical morphology is a powerful tool for image analysis; however, classical morphological operators suffer from lacks of robustness against noise, and also intrinsic image features are not accounted at all in the process. We propose in this work a new and different way to overcome such limits, by introducing both robustness and locally adaptability in morphological operators, which are now deﬁned in a manner such that intrinsic image features are accounted. Dealing with partial differential equations (PDEs) for generalized Cauchy problems, we show that proposed PDEs are equivalent to impose robustness and adaptability of corresponding sup-inf operators, to structuring functions. Accurate numerical schemes are also provided to solve proposed PDEs, and experiments conducted for both synthetic and real images, show the efﬁciency and robustness of our approach.


INTRODUCTION
In many image processing and computer vision tasks (e.g., data compression, feature detection, motion analysis/detection, multiband frequency analysis, . . .), it is important to perform a multiscale image analysis, i.e., analyze the image at multiple spatial scales or resolutions.Owing to that fact, mathematical morphology (Serra, 1982) appeared as a powerful tool in multiscale analysis (Soille, 1999), mainly due to its nonlinearity aspects, shape and geometry description properties.Despite its interesting properties and plentiful successful applications in various domains, both former discrete and continuous morphological operators suffer from a lack of robustness against noise, and do not take into account the spatial adaptability.The same global way in which all image pixels are treated is in fact main reasons for that limitations.This is avoided in Lerallut et al. (2007) by using morphological filters with non fixed shape kernels (or amoebas), and in Angulo (2011) by considering bilateral structuring functions.Using a PDE approach, in order to enhance coherence of flowlike structures, an adaptive method was obtained in Breuß et al. (2007) by multiplying the image gradient term that appeared in classical morphological PDEs for dilations/erosions, with a space-variant matrix.

Motivation and goal
The main motivation of this work is to propose alternatives for remedying such weaknesses of robustness and adaptability.This is done by considering Perona and Malik-like diffusivity term (Perona and Malik, 1990).Let us first recall some fundamentals of this nonlinear image diffusion model.The pioneering idea introduced by Perona and Malik (1990) to reduce blurring and edge shifting of uniform linear diffusion, involved a locally adaptive diffusion with respect to the gradient of the image at each iteration.More precisely, this nonlinear diffusion involves replacing the heat transfer diffusion equation ∂ u/∂t = ∆u, by the following model: In this model (Eq.1), the diffusivity has to be such that g ∇u 2 → 0 when ∇u 2 → ∞, and g ∇u 2 → 1 when ∇u 2 → 0. One of the diffusivities proposed by Perona and Malik is the following: where λ is a threshold parameter that separates forward and backward diffusion.This model accomplishes the aim of blurring small fluctuations (noise) while enhancing edges (by preventing excessive diffusion).To avoid theoretical as well as numerical drawbacks on this model, Catté et al. (1992) proposed a new version of Perona and Malik theory, by replacing the diffusivity term g ∇u 2 with its regularized version, i.e., g ∇u σ 2 , where ∇u σ = ∇ (G σ * u), and G σ is a Gaussian kernel with a standard deviation equal to σ .This latter model is nowadays extremely popular in PDE-based image processing, and many variants have then been proposed.
Concretely, we wish in this work to meet two major objectives.On a first hand, we aim at exploring how a gradient-based local adaptability such as g ∇u 2 , can be introduced in both morphological PDEs and sup-inf based operators.On a second hand, we would like to propose associated numerical schemes for solving efficiently proposed methods.
State-of-the art on discrete adaptive morphology.Original formulation of discrete dilation and erosion for gray-level images (Serra, 1982;Sternberg, 1986), as well as the other morphological operators was translation invariant in space and intensity.Later, mathematical morphology was formulated in the framework of complete lattices (Serra, 1988;Heijmans and Ronse, 1990), which leads to a general case of dilation and adjoint erosion compatible with spatially-variant operators.Nevertheless, most of current implementations and classical applications studied in the literature are based on morphological operators, which are translation invariant in space and intensity (Soille, 1999), i.e., a same structuring function (or structuring element) is considered for each pixel of the image.A current active field in mathematical morphology is the construction of appropriate adaptive operators; i.e., structuring functions become dependent on the position or on the input image itself.In previous works, adaptive operators are based on two main approaches.On a first hand, a variability is considered on the support space of pixels, meaning that spatially variable shape of structuring functions are adapted according to: positions in the image, as in Beucher et al. (1987); Cuisenaire (2006); local regularity including the morphological amoebas (Lerallut et al., 2007), adaptive geodesic neighborhoods (Grazzini and Soille, 2009), bilateral structuring function (Angulo, 2013), nonlocal structuring function (Velasco-Forero and Angulo, 2013); local orientations, including typically anisotropic gradient-driven operators (Verdú-Monedero et al., 2011).
On another hand, a variability on T is considered: variable size of structuring functions according to the local intensity or contrast (Vachier and Meyer, 2007;Maragos and Vachier, 2008).
For an overview on the state-the-art on adaptive morphology, the interested reader is invited to the paper Maragos and Vachier (2009) and the most recent one ( Ćurić et al., 2014).Another milestone study (Roerdink, 2009) is very interesting for understanding the theoretical limitations of inputadaptive morphological operators.
State-of-the art on PDEs-based morphological operators.In 1992, three independent milestones papers appeared on nonlinear PDEs for modeling continuous multiscale morphological operators: -In Alvarez et al. (1992;1993), authors obtained PDEs for multiscale flat dilations and erosions, by means of compact convex structuring sets as part of their general work on developing PDE-based models for multiscale image processing that satisfy certain axiomatic principles.
-Authors in Brockett and Maragos (1992;1994) developed nonlinear PDEs that model multiscale morphological dilation/erosion, opening/closing, by using compact support structuring elements that are either convex sets or concave functions and may have non smooth boundaries.Appropriate numerical schema were provided as well.
-PDEs for multiscale dilation and erosion were obtained (van den Boomgaard and Smeulders, 1992;1994) in by studying the propagation of 2D boundaries sets and signal graphs, under multiscale dilations and erosions.Their work applies to convex structuring elements whose boundaries contain no linear segments, are smooth, and possess a unique normal at each point.
Paper contributions.We provide herein a new approach to make classical multiscale dilation and erosion operators, robust and spatially adaptive as regards to intrinsic image features.To constraint both continuous and sup-inf based models to meet our target objectives, we choose to investigate PDEs for generalized Cauchy problems.Parameters of such models are then linked to image features in two ways; the last one being extremely robust in a noisy environment.Moreover, suitable numerical schemes are provided for solving efficiently proposed PDEbased models.
Manuscript organization.The paper is organized as follows.We recall in the next section some background on mathematical morphology.After that, we present our proposed robust and adaptive multiscale approach in another section.Afterwards, numerical schemes for solving proposed models are described, and obtained results are shown.The paper ends with some discussions and outline of some interesting perspectives.

BACKGROUND ON MATHEMATICAL MORPHOLOGY
Algebraic lattice theory.From a theoretical viewpoint, mathematical morphology is nowadays formulated in the algebraic framework complete lattice theory (Serra, 1988;Heijmans, 1994).More precisely, the notion of adjunction links two operators (ε, δ ), in such a way that for any given dilation δ , there is a unique erosion ε such that (ε, δ ) is an adjunction (Heijmans and Ronse, 1990).In this paper, we adopt the alternative paradigm based on functional analysis, which leads to the formulation of dilation and erosion as the solution of PDEs.
Basic morphological operators.Since its introduction, mathematical morphology (Matheron, 1975;Serra, 1982) appeared as a powerful tool in image analysis (Soille, 1999).This is mainly due to its nonlinearity, shape and geometry description properties.Let E ⊆ Z 2 .For any function f : E → R, elementary dilation and erosion operators (obtained by adjunction and duality) are respectively defined by: where (resp. ) represents the supremum (resp.infimum).The further convention for scalar addition in is considered: We notice that dilation and erosion are related respectively to the supremal and infimal convolution known in convex analysis (Rockafellar, 1970).The structuring function b : R 2 → R is typically a general concave function.A simple case of flat dilation and erosion results when b equals to 0 in a convex bounded set Both dilation and erosion are invariant under translations in E (spatial or "horizontal" direction) and R (intensity or "vertical" direction).
Definition 2.1 Let F be a family of real functions mapped on E ⊆ R 2 .An operator S : In addition to that, dilation and erosion are increasing operators, and also satisfy the following properties: -Duality: -Adjunction: Different morphological filters can be obtained by composition of above operators (Matheron, 1975;Serra, 1982;Catté et al., 1995;Cao, 1998).Using any structuring function b, two interesting operators, namely the opening and closing, can be obtained respectively in the following way: Opening and closing are also increasing operators with two main properties: -Ordering (anti-extensivity and extensivity): For any structuring function, one has -Idempotence: For any structuring function, one has . By duality, one also has: Semi-group property, multiscale morphological filters and morphological PDEs.For t ≥ 0, one can define then multiscale dilations (and erosions) by replacing b by the family of multiscale structuring functions (b t ) t≥0 defined for t > 0 by: which satisfies the semi-group property (b s ⊕ b t )(x) = b s+t (x, y).The canonical multiscale structuring function, which can be seen as the morphological counterpart of the Gaussian kernel in linear filtering, is the infinity support parabolic structuring function From the works by Van den Boomgaard (van den Boomgaard and Dorst, 1997), it is also well known that this structuring functions are the equivalent class of functions which are dimensionally separable and closed with respect to the dilation/erosion.This result is usually proved in the slope transform domain (Dorst and van den Boomgaard, 1994;Maragos, 1995).Multiscale dilation and erosion filters are then respectively given by: Let us mention that in the case of flat structuring functions, this is equivalent to consider sets B t = tB (i.e., homothetic compact convex planar sets B of size t) as multiscale structuring elements.
Note that multiscale openings and closings can be defined in a same way by considering a family (b t ) t≥0 of multiscale structuring functions, respectively by: Alternatives to perform multiscale continuous flat dilations and erosions (dual operators) by using partial differential equations (PDEs) were proposed in (Alvarez et al., 1993;Sapiro et al., 1993;Brockett and Maragos, 1994): with (+) (resp.(−)) sign in Eq. 13 stands for the multiscale dilation (resp.erosion).Useful cases of unit structuring sets B are obtained by unit balls B | p = x ∈ R n : x p ≤ 1 of metrics induced by the L p norms.There are three special cases of norms p = 1, 2 and p = ∞ which correspond to particular PDEs: for B | 1 (i.e., rhombus), one has ∂ u t = ± ∇u ∞ ; for B | 2 (i.e., disk), one has ∂ u t = ± ∇u 2 ; and for B | ∞ (i.e., square), one has ∂ u t = ± ∇u 1 .Similarly, the PDE for multiscale parabolic dilations and erosions is given by: It is well known that both dilation/erosion PDEs, Eqs. 13 and 14 are special cases of Hamilton-Jacobi equations, which are of great interests in physics.In fact, let us consider the general Hamilton-Jacobi family of problems: (15) Such equations usually do not admit classic (i.e., everywhere differentiable) solutions, but they can be studied in the framework of viscosity solutions theory (Crandall et al., 1992).It is well known (Lions, 1982;Bardi and Evans, 1984) that if the Hamiltonian H(x, p) = H(p) is convex, then, the solution of the Cauchy problem are respectively given for the (+) and (−) signs in Eq. 15 by: where H * is the Legendre-Fenchel transform of function H.
Finally, let us quote works in (Alvarez et al., 1993) in which authors proposed the following PDE for multiscale continuous openings, given for any scale s > 0 by: with t ∈ [0, 2s]; sign(• ) stands for the signum function and r + := max(r, 0).Indeed, for t ∈ [0, s], PDE in Eq. 18 acts as a multiscale erosion, while for t ∈ ]s, 2s] it is in fact a multiscale dilation.Multiscale closings can be obtained by switching (+) and (−) signs in the two terms of PDE in Eq. 18.

PROPOSED ROBUST AND SPATIALLY ADAPTIVE SCALE-SPACES
It is clear that neither the discrete, Eqs. 2 and 3, nor the continuous, Eq. 13, formulations are robust in a noisy environment.In fact, in the presence of noise, taking the supremum in Eq. 2 or the infimum in Eq. 3 will definitely lead to wrong values, while hyperbolic PDEs as the one in Eq. 13 will blow up.Main reason for that is the fact that all image pixels are treated in a same global way.As previously said, many works were proposed for avoiding this issue (Breuß et al., 2007;Lerallut et al., 2007;Angulo, 2011).Let us point out the fact that the proposed approach in (Breuß et al., 2007) is truly adaptive, but is not robust at all, and sophisticated strategies were used to implement the orientation information in the model.In a recent study (Diop and Angulo, 2012b), we overcome those drawbacks by proposing different robust and spatially adaptive PDEs.The first proposed model was a Gaussian regularization of Eq. 13, as follow: where u σ = u * G σ is the convolution with a Gaussian G σ of variance σ .This first model solves the robustness against noise; however, it has a major side effect, that is, it creates some blur, and the latter increases as one goes further through scales.This could also be expected; because, iterating the Gaussian convolution is asymptotically equivalent to solve the Heat equation, which is an isotropic diffusion.Indeed, in a geometrical view, level lines are smoothed both in the normal and tangent directions.We also proposed a second PDE model that behaved in better way than the previous one, as follow (Diop and Angulo, 2012b): where g is an edge detector function.In addition to noise robustness, such a model is also locally adaptive to intrinsic image features, e.g., edges.We looked at cases where g is either a decreasing or an increasing function, defined by: for g decreasing, (21) where • stands for the Euclidean norm in R 2 .The monotonicity effects of g will be illustrated later in the numerical results section.
We present in next sections a different and new approach in making robust and adaptive morphological multiscale image analysis methods.In fact, contrary to (Diop and Angulo, 2012b) where the morphological analysis was mainly focused on the processed image, we propose here to work on structuring functions.
Generalized Cauchy problems for continuous multiscale models.In order to force both the robustness and adaptability, we propose here1 a different morphological multiscale approach, by making structuring functions depending on intrinsic image features such as edges.To achieve this, we propose to use the generalized Cauchy problems: with p > 1.The following remark is the starting point of our motivation.
Remark 3.1 states also that different kinds of structuring functions can be obtained depending on the value of p.It is then interesting to see how this parameter p behaves when p → 1, p ∈]1, 2[ and p → ∞, and also how it is linked to the structuring functions in the case of sup-inf formulations.This is in fact done by looking at the viscosity solution of Eq. 22, for example with the + sign, which is given by the following Lax-Oleinik formula (Lions et al., 1987): where c p = p − 1 p p/(p−1) . ( Let us illustrate some dilations carried out by using Eq.23.In fact, we consider an image with ten pixels chosen at random and perform dilations at different scales t and for different p also; results are shown in For a better comprehension of the effects of p on the structuring function, let us have a closer look for at the family of concave functions (k t,p ) t>0 defined for t > 0 and p > 1 by First of all, notice that for a fixed t > 0, when p → 1, then k t,p → 0; when p → ∞, then for all x, k t,p (x) → x .What is particularly interesting is to see the behavior of the family of functions x → k t,p (x) in the unit neighborhood of x.This is illustrated in Fig. 2 where different plottings of k t,p (x) are obtained for t = 1 and x ∈ [−1, 1].-Firstly, as expected, one is dealing with flat structuring functions as p → 1, on the one hand.On the other hand, as p increases up to 2, supports of of structuring functions k t,p (x) diminish and their shapes tend to a parabola, see Fig. 2a.
-Secondly, for p > 2, as p increases and p → ∞, shapes of k t,p (x) evolve from a parabola to the limit case, i.e., |x|, see Fig. 2b.
Accounting these two facts and also Remark 3.1, the fundamental idea of our approach is to make p depending on most relevant image features, i.e., gradients, which is in fact a nice way to locally adapt structuring functions to intrinsic image features.Thus, we propose two ways in doing that.
Adaptability approach based on image features.The basic principle consists in making p = p(u) as a function depending on image gradients, i.e., p = f (∇u) and 1 < p < ∞, in a way such that p is close to 1 (around the neighborhood of the considered pixel) in homogeneous image regions, and p belongs to [2, ∞[ in inhomogeneous image areas, e.g., contours, noise, texture, • • • .Typical p = f (∇u) is given in Fig. 3. To this aim, we propose to set in the models of Eqs.22 and 23, where g is a decreasing function and given as in Eq. 21.Multiscale image analysis results obtained with this approach are shown and discussed previously.
Adaptive coupled model with image features and edge threshold.As shown in the previous section, potentially good results are obtained with the preceding approach.In this section, we wish to significantly increase the robustness against noise.For doing so, we propose to enhance effects of the edgebased parameter p = p(u) = f (∇u) by the means of an edge threshold α which determines whether or not one is on image contours.
Let h be a function that detects the image contours.The proposed method is formulated as follows: Results In this section, we show multiscale image analysis results obtained with the PDE of Eq. 22 and using either the approach based only on image features or coupled with an edge threshold.The proposed generalized Cauchy-based PDE is solved with suitable numerical schemes like in Eqs.32 and 33.Other considerations are also accounted; for instance, in the discretization of g, which intervenes in the edgedependent parameter p = f (∇u) in both the featurebased and coupled models.
Obtained results show the efficiency of our locally adaptive approaches for multiscale image analysis.In fact, we first apply such methods on the noisy binary Cat image, given in Fig. 5(a), and obtained by adding a Gaussian noise with an SNR = 33.29 dB.
Results of the proposed multiscale dilations and erosions are respectively shown in Figs. 6 and 7. Results are compared to ones carried out by applying the classical PDE-based approach, corresponding to model of Eq. 13.It is not new that pixels are treated in a same global way, which involves results in a melting between the noise and image features, as one is going through scales; first columns from the left in Figs. 6 and 7.This fact has been already discussed and illustrated in previous sections.On the contrary, one can see the improved behavior of both proposed approaches in the sense of adaptability and robustness against noise, and in an even better manner by using the edge threshold-based technique coupled with image features, see third columns from the left in Figs. 6 and 7).
Let us point out another interesting fact regarding the coupled method.In fact, as depicted in Fig. 7(i), one could notice the robustness against noise through scales, and also how the cat whiskers are preserved.By noise robustness, we mean that the noise does not affect at all the performed erosions through scales.To put more emphasis on that fact, we show more erosions in coarser scales in Fig. 8.In fact, due to the robustness against noise, thin image features are preserved through scales, which results in a skeletonization of the cat.Thin cat features have high gradients and are then kept similarly to the noise.
Note.Skeletal abstraction is a difficult problem that has been greatly studied over years; and obviously it is out of the paper scope.Briefly saying, existing methods for extracting skeletons concern broad research areas comprising topological thinning algorithms (Arcelli and Baja, 1985;Lee and Kashyap, 1994;Borgefors et al., 1999;Bertrand and Couprie, 2009) where Blum grassfire transform (Blum, 1973) were used, curve evolution, variational and wavefront propagation methods (Leymarie and Levine, 1992;Geiger et al., 2003;Tek and Kimia, 2003), Voronoi diagram (Schmitt, 1989;Ogniewicz, 1993;Sheehy et al., 1996), and methods using Euclidean distance function computed for example with the Eikonal equation or Hamilton-Jacobi systems (Siddiqi et al., 2002;Torsello and Hancock, 2003;2004).For more information on that subject, interested readers can have a look on those references.From top to bottom − First line, 5 iterations using (a) classical approach, Eq. 13, (b) generalized Cauchy method, Eq. 22, where p depends only on image features, (c) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.Second line, 15 iterations using (d) classical approach, Eq. 13, (e) generalized Cauchy method, Eq. 22, where p depends only on image features, (f) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.Third line, 50 iterations using (h) classical approach, Eq. 13, (i) generalized Cauchy method, Eq. 22, where p depends only on image features, (j) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.From top to bottom − First line, 5 iterations using (a) classical approach, Eq. 13, (b) generalized Cauchy method, Eq. 22, where p depends only on image features, (c) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.Second line, 15 iterations using (d) classical approach, Eq. 13, (e) generalized Cauchy method, Eq. 22, where p depends only on image features, (f) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.Third line, 50 iterations using (h) classical approach, Eq. 13, (i) generalized Cauchy method, Eq. 22, where p depends only on image features, (j) generalized Cauchy method, Eq. 22, based on coupled edge threshold-based method.

CONCLUSION
We have provided here some contributions concerning major issues on former morphological operators; for instance, lacks of robustness against noise and adaptability to image features.Our proposed approach is different from what we proposed in (Diop and Angulo, 2012b), in the sense that constraints of robustness and adaptability were not applied anymore on the image itself, but on structuring functions.Also, in this work, we have investigated PDEs for generalized Cauchy problems for modeling robust and adaptive morphological scale-spaces.In addition, the parameter p of the PDE model and the associated supinf operators given by Lax-Oleinik formula, have been linked to intrinsic edge image features in two different ways.As shown by the obtained results, the approach for choosing p thanks to an edge detector threshold has appeared, as expected, to be extremely robust in a noisy environment.Finally, numerical schemes have also been proposed as well for the resolution of all proposed PDEs.The efficiency of our approaches has been illustrated on binary and gray-level images.
As for future work, we plan to extend this study to color, and more generally, to multichannel images.This is not a straightforward extension of what we have proposed herein, because a channel-wise based approach will undoubtedly result to color artifacts, in the sense that new (artificial) colors that do not primary exist will appear in carried out result.We plan to investigate the channel-wise method, as well as a more challenging vectorial approach as discussed in Bresson and Chan (2008); Strekalovskiy et al. (2012), even if contexts were different.Finally, it would be interesting to use in edge functions, a time-delay regularization (Belahmidi and Chambolle, 2005) for avoiding the Gaussian regularization with fixed variances, on the one hand, and for locally adaptability too, on the another hand.

Fig. 2 .
Fig. 2. Plottings of structuring functions k t,p (x) for a fixed t = 1 and different values of p.

Fig. 3 .
Fig. 3. Typical p = f (∇u) functions for a PDE model based on image features.

Fig. 8 .Fig. 9 .
Fig. 8. Multiscale erosions of the noisy binary Cat image performed with the adaptive PDE of model Eq.22 based on the coupled edge threshold-based method: (a) 100 iterations, (b) 200 iterations and (c) 2000 iterations.