Image Anal Stereol 2015;34:39-50 doi:10.5566/ias.993 Original Research Paper MULTISCALE IMAGE ANALYSIS BASED ON ROBUST AND ADAPTIVE MORPHOLOGICAL SCALE-SPACES EL HADJI S. DIOP[,† AND JES ´US ANGULO Center for Mathematical Morphology, Department of Mathematics and Systems, MINES ParisTech, 35 rue Saint-Honor´ e, 77305 Fontainebleau Cedex -France e-mail: ehsdiop@hotmail.com, jesus.angulo@mines-paristech.fr (Received January 14, 2013; revised January 20, 2014; accepted March 25, 2014) ABSTRACT 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. Keywords: adaptability, morphological operators, partial differential equations, robustness. 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 .lters with non .xed shape kernels (or amoebas), and in Angulo (2011) by considering bilateral structuring functions. Using a PDE approach, in order to enhance coherence of .ow-like 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 .rst 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: .u nn hh = div g I.uI2.u , .t (1) u(x,0)= f (x) . In this model (Eq. 1), the diffusivity has to be such that nh nh g I.uI2› 0 when I.uI2 › ., and g I.uI2› 1 when I.uI2 › 0. One of the diffusivities proposed by Perona and Malik is the following: nh g(s2)= 1/ 1 +(s2/.2) , . > 0 , where . is a threshold parameter that separates forward and backward diffusion. This model accomplishes the aim of blurring small .uctuations (noise) while enhancing edges (by preventing excessive diffusion). To avoid theoretical as well as numerical drawbacks on this model, Catte´et al. (1992) proposed a new version of Perona and Malik nh theory, by replacing the diffusivity term g I.uI2 nh with its regularized version, i.e., g I.u. I2, where .u. = .(G. * u), and G. is a Gaussian kernel † Presently with the Departement of Mathematics, University of Thies, Cit´es, S´egal e Malick SY, BP 967, Thi`en´ 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 .rst hand, we aim at exploring how a gradient-based local adaptability such as g n I.uI2h , 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 ef.ciently 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 .eld 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 .rst 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), non-local structuring function (Velasco-Forero and Angulo, 2013); – local orientations, including typically anisotropic gradient-driven operators (Verd´u-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 Curi´et al., 2014). Another milestone (´ c study (Roerdink, 2009) is very interesting for understanding the theoretical limitations of input-adaptive 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 .at 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 ef.ciently proposed PDE-based 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. Image Anal Stereol 2015;34:39-50 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 . Z2. For any function f : E › R, elementary dilation and erosion operators (obtained by adjunction and duality) are respectively de.ned by: ( f . b)(x) :=[ f (y)+ b(x - y)] , (2) y.E ( f 8 b)(x) :=[ f (y) - b(y - x)] , (3) y.E where(resp.) represents the supremum (resp. in.mum). The further convention for scalar addition in is considered: f (y)+b(x-y)= -. when b(x - y)= -. or f (y)= -., and that f (y) - b(y - x)=+. when b(y+x)=+. or f (y)= -.. We notice that dilation and erosion are related respectively to the supremal and in.mal convolution known in convex analysis (Rockafellar, 1970). The structuring function b : R2 › R—is typically a general concave function. A simple case of .at dilation and erosion results when b equals to 0 in a convex bounded set B . E called structuring element and -. outside, i.e., b(x)= 0 if x . B and b(x)= -. if x . Bc: ( f . B)(x) = [ f (x - y)] , (4) y.B ( f 8 B)(x) = [ f (x + y)] . (5) y.B Both dilation and erosion are invariant under translations in E (spatial or “horizontal” direction) and — R (intensity or “vertical” direction). De.nition 2.1 Let F be a family of real functions mapped on E . R2 . An operator S : F › F is increasing or monotone if for all f , g . F s.t. f (x) . g(x), .x . E, then, (Sf )(x) . (Sg)(x), .x . E. In addition to that, dilation and erosion are increasing operators, and also satisfy the following properties: – Duality: ( f . b)(x)= -(- f 8 b)(x), .x . E. – Adjunction: ( f .b)(x) . g(x), .x . E .. f (x) . (g 8 b)(x), .x . E. Different morphological .lters can be obtained by composition of above operators (Matheron, 1975; Serra, 1982; Catte´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: ( f . b) := [( f 8 b) . b] , (6) ( f • b) := [( f . b) 8 b] . (7) Opening and closing are also increasing operators with two main properties: – Ordering (anti-extensivity and extensivity): For any structuring function, one has ( f . b) . f . ( f • b) – Idempotence: For any structuring function, one has [( f . b) . b]= ( f . b). By duality, one also has: [( f • b) • b]=( f • b). Semi-group property, multiscale morphological .lters and morphological PDEs. For t . 0, one can de.ne then multiscale dilations (and erosions) by replacing b by the family of multiscale structuring functions (bt )t.0 de.ned for t > 0 by: . . tb(x/t) for t > 0 bt (x)= 0 for t = 0, x = 0 . -. otherwise, which satis.es the semi-group property (bs . bt )(x)= bs+t (x,y). The canonical multiscale structuring function, which can be seen as the morphological counterpart of the Gaussian kernel in linear .ltering, is the in.nity support parabolic structuring function IxI2 bt (x)= - . (8) 2t 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 .lters are then respectively given by: ( f . bt )(x) :=[ f (y)+ bt (x - y)] , (9) y.E ( f 8 bt )(x) :=[ f (y) - bt (y - x)] . (10) y.E Let us mention that in the case of .at structuring functions, this is equivalent to consider sets Bt = tB (i.e., homothetic compact convex planar sets B of size t) as multiscale structuring elements. Note that multiscale openings and closings can be de.ned in a same way by considering a family (bt )t.0 of multiscale structuring functions, respectively by: ( f . bt ) := [( f 8 bt ) . bt ] , (11) ( f • bt ) := [( f . bt ) 8 bt ] . (12) Alternatives to perform multiscale continuous .at 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): .u = ±I.uI , .t (13) u(x,0)= f (x) , 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 = sa x . Rn : IxIp . 1 of metrics induced by the Lp 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 . ut = ±I.uI.; for B |2 (i.e., disk), one has .ut = ±I.uI2; and for B |. (i.e., square), one has .ut = ±I.uI1. Similarly, the PDE for multiscale parabolic dilations and erosions is given by: .u 1 = ±I.uI2 , .t 2(14) u(x,0)= f (x) . 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: .u(x,t) = ±H (x,.u(x,t)) = 0, in Rn × (0,+.) .t u(·,0)= f in Rn . (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: x - y u(x,t)= inf f (y)+tH*, (16) y.Rnt x - y u(x,t)= sup f (y) -tH*, (17) y.Rnt 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: .u = -sign+(s -t)I.uI + sign+(t - s)I.uI .t u(x,0)= f (x), (18) 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 in.mum in Eq. 3 will de.nitely 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 .rst proposed model was a Gaussian regularization of Eq. 13, as follow: . u = ±I.u. I in RN × (0,T ) , .t (19) u(x,0)= I(x) in RN ×{t = 0} , Image Anal Stereol 2015;34:39-50 where u. = u * G. is the convolution with a Gaussian G. of variance .. This .rst 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): . u = ±g(I.u. I)I.uI in RN × (0,T ) , (20) .t – .at structuring function, for p = 1, – parabolic structuring function, for p = 2. 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): Ix - yIp/(p-1) I(y) - cp (23) u(x,t)= sup t1/(p-1), u(x,0)= I(x) in RN ×{t = 0} , y.RN where g is an edge detector function. In addition to noise robustness, such a model is also locally adaptive where to intrinsic image features, e.g., edges. We looked at p - 1 = (24) cp . cases where g is either a decreasing or an increasing function, de.ned by: pp/(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 g(x)= . . . . . 1 + .IxI2p, p > 0 for g increasing, 1 + . exp(IxI2p), p > 0 for g increasing, scales t and for different p also; results are shown in Fig. 1. 1 for g decreasing, 1 + .IxI2 (21) where I·I stands for the Euclidean norm in R2. 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: .u = ±IDuIp in RN × (0,T ) , .t (22) u(x,0)= I(x) in RN ×{t = 0} , with p > 1. The following remark is the starting point of our motivation. Remark 3.1 PDE in problem Eq. 22 corresponds to multiscale dilations/erosions with: (a) (b) (c) (d) Fig. 1. (a) Original image with 10 pixels chosen at random. Multiscale dilations performed using model of Eq. 23 with: (b) t = 15 and p = 1.01, (c) t = 20 and p = 1.5, and (d) t = 10 andp = 2. 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 (kt,p)t>0 de.ned for t > 0 and p > 1 by IxIp/(p-1) kt,p(x) := -cp . (25) t1/(p-1) 1The idea was .rst presented in S4G conference (Diop and Angulo, 2012a). 43 First of all, notice that for a .xed t > 0, when p › 1, then kt,p › 0; when p › ., then for all x, kt,p(x) › IxI. What is particularly interesting is to see the behavior of the family of functions x › kt,p(x) in the unit neighborhood of x. This is illustrated in Fig. 2 where different plottings of kt,p(x) are obtained for t = 1 and x . [-1, 1]. (b) Fig. 2. Plottings of structuring functions kt,p(x) for a .xed t = 1 and different values of p. Indeed, basically saying, two facts are shown: – Firstly, as expected, one is dealing with .at structuring functions as p › 1, on the one hand. On the other hand, as p increases up to 2, supports of of structuring functions kt,p(x) diminish and their shapes tend to a parabola, see Fig. 2a. – Secondly, for p > 2, as p increases and p › ., shapes of kt,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. Fig. 3. Typical p = f (.u) functions for a PDE model based on image features. To this aim, we propose to set p(u)= 1 + g(I.u. I) , (26) 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 signi.cantly increase the robustness against noise. For doing so, we propose to enhance effects of the edge-based 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: Image Anal Stereol 2015;34:39-50 – In homogeneous image regions, i.e., h(I.uI) < ., we set p(u)= p1 = 1 + g1(I.u. I) (27) in the models of Eqs. 22 and 23. – In high frequency image areas, i.e., h(I.II) . ., we set p(u)= p2 = c + g2(I.u. I) (28) in the models Eqs. 22 and 23, with c . 2. For instance, g1 and g2 are decreasing edge-based functions taken as in Eq. 21, and g2 is decreasing more rapidly than g1. Fig. 4 illustrates what we just discussed about the edge value . and typical p(u)= f (.u) function used in this work. Fig. 4. Typical p(u)= f (.u) functions for coupled with image features and threshold edge value .. NUMERICAL RESULTS Discretization Let us set up the following discretization grid: xi = iDx, yj = jDy, tn = nDt, i, j = 0,±1,±2,··· , and n = 0,1,2,··· (29) An approximate solution u on the above grid is denoted by u(xi,yj), and satis.es: un(xi,yj) . u(iDx, jDy,nDt) . (30) Discretization. Proposed PDEs in Eqs. 19 and 20 are solved using an explicit scheme: un+1(xi, yj)= un(xi,yj) ±Dt·F(un(xi,yj))(Dun)±(xi,yj) , (31) where F(un(xi,yj)) . 1 and un = un for Eq. 13, un = un * G. for Eq. 19, and in the case of dealing with Eq. 20, F(un(xi,yj)) corresponds to the discretization of the edge function g(I.u. I). Hence, the (+) (resp. (-)) sign in the term (Dun)± of Eq. 31 stands for multiscale dilations (resp. erosions). It is known that discretization of the hyperbolic term could be obtained by using classical schemas introduced in (Osher and Sethian, 1988; Rouy and Tourin, 1992). New schemes are also proposed in this work. Before going through that, let us .rst consider the following forward and backward derivatives: un(xi + 1,yj) - (un(xi, yj) D+un(xi,yj)= , x Dx un(xi,yj) - (un(xi - 1,yj) D-un(xi,yj)= , x Dx un(xi,yj + 1) - (un(xi,yj) D+un(xi,yj)= , y Dy un(xi,yj) - (un(xi,yj - 1) D- y un(xi,yj)= . Dy Thus, our proposed schemes are the followings: n (Dun)+ = max(max(0,D+un)2 ,max(0, -D-un)2)+ x x h1/2 max(max(0,Dy +un)2 ,max(0,-Dy +un)2, (32) n (Dun)- = max(min(0,D+un)2 ,min(0, -D-un)2)+ x x h1/2 max(min(0,D+un)2 ,min(0,-D+un)2. y y (33) Note that the discretization of .u. in g(I.u. I) is not the same as in the hyperbolic term. Indeed, we notice that best schemes were obtained by adding all spatial derivatives all together, i.e., n nnnn I.u. I2 =(D+u. )2 +(D-u. )2 +(D+u. )2 +(D-u. )2 . xxyy (34) (a) (b) Fig. 5. (a) Noisy Cat image. (b) Owl image. 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 edge-dependent parameter p = f (.u) in both the feature-based and coupled models. Obtained results show the ef.ciency of our locally adaptive approaches for multiscale image analysis. In fact, we .rst 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; .rst 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). (a) (b) (c) (d) (e) (f) (h) (i) (j) Fig. 6. Multiscale dilations of the noisy binary Cat image given in Fig. 5(a). 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. (a) (b) (c) (d) (e) (f) (h) (i) (j) Fig. 7. Multiscale erosions of the noisy binary Cat image given in Fig. 5(a). 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. Image Anal Stereol 2015;34:39-50 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 dif.cult problem that has been greatly studied over years; and obviously it is out of the paper scope. Brie.y 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 grass.re 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. (a) (b) (c) 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. Proposed methods are also applied on the noisy grayscale image displayed in Fig. 5(b). Multiscale dilations and erosions are respectively depicted in Figs. 9 and 10. Other morphological operators are also applied to this image; namely, openings and closings, where the corresponding results are given in Fig. 11. All these examples con.rm the improved behavior in terms of noise robustness and locally adaptability of our proposed approaches; especially, the edge threshold-based method presented in the previous section, which correspond to examples in Figs. 11(c) and 11(f). (h) (i) (j) Fig. 9. Multiscale dilations of the noisy Owl image given in Fig. 5(b). 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. (a) (b) (b) (a) (b) (b) (h) (i) (j) Fig. 10. Multiscale erosions of the noisy Owl image given in Fig. 5(b). 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. (d) (e) (f) (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, multiscale closings at scale 15 using: (d) classical PDE, 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. 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 sup­inf 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 ef.ciency 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 (arti.cial) 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 Image Anal Stereol 2015;34:39-50 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 .xed variances, on the one hand, and for locally adaptability too, on the another hand. REFERENCES Alvarez L, Guichard F, Lions PL, Morel JM (1993). Axioms and fundamental equations of image processing. Arch Rational Mech Anal 123:199–257. Alvarez L, Lions PL, Morel JM (1992). Image selective smoothing and edge detection by nonlinear diffusion. SIAM J Numer Anal 29:845–66. Angulo J (2011). Morphological bilateral .ltering and spatially-variant adaptive structuring functions. In: Soille P, Pesaresi M, Ouzonis GK, eds. Math Morphol Appl Image Signal Proc. Lect Not Comput Sci 6671:212–23. Angulo J (2013). Morphological bilateral .ltering. SIAM J Imaging Sci 6:1790–822. Arcelli C, Baja GSD (1985). A width-independent fast thinning algorithm. IEEE T Pattern Anal 7:463–74. Bardi M, Evans LC (1984). On Hopf’s formulas for solutions of Hamilton-Jacobi equations. Nonlinear Anal Theor 8:1373–81. Belahmidi A, Chambolle A (2005). Time-delay regularization of anisotropic diffusion and image processing. ESAIM Math Model Num 39:231–51. Bertrand G, Couprie M (2009). On parallel thinning algorithms: minimal non-simple sets, P-simple points and critical kernels. J Math Imaging Vis 35:23–35. Beucher S, Blosseville JM, Lenoir F (1987). Traf.c spatial measurements using video image processing. In: Casasent DP, Hall EL, eds. Proc Conf Intel Robot Comput Vis VI. Proc SPIE 0848 Blum H (1973). Biological shape and visual science, part I. J Theor Biol 38:205–87. Borgefors G, Nystr¨om I, Baja GSD (1999). Computing skeletons in three dimensions. Pattern Recogn 32:1225– 36. Bresson X, Chan TF (2008). Fast dual minimization of the vectorial total variation norm and applications to color image processing. Inverse Probl Imaging 2:455–184 Breuß M, Benedi B, Weickert J (2007). Anisotropic continuous-scale morphology. In: Mart´i J, Bened´i JM, Mendonc¸a AM, Serrat J, eds. Pattern Recogn Image Anal. Lect Not Comput Sci 4478:515–22. Brockett RW, Maragos P (1992). Evolution equation for continuous-scale morphology. In: Proc IEEE Conf Acoust Speech Signal Proc (ICASSP’92). 3:125–8. Brockett RW, Maragos P (1994). Evolution equations for continuous-scale morphological .ltering. IEEE T Signal Proc 42:3377–386. Cao F (1998). Partial differential equations and mathematical morphology. J Math Pure Appl 77:909– 41. Catte´F, Dibos F, Koep.er G (1995). A morphological scheme for mean curvature motion and applications to anisotropic diffusion and motion of level sets. SIAM J Numer Anal 32:1895–909. Catt´ e F, Lions PL, Morel JM, Coll T (1992). Image selective smoothing and edge detection by nonlinear diffusion. SIAM J Numer Anal 29:182–93. Crandall MG, Ishii H, Lions PL (1992). User’s guide to viscosity solutions of second order partial differential equations. Bull Amer Math Soc 27:1–67. Cuisenaire O (2006). Locally adaptable mathematical morphology using distance transformations. Pattern Recogn 39:405–16. ´ Curi´c V, Landstr¨om A, Thurley MJ, Luengo-Hendriks CL (2014). Adaptive mathematical morphology – a survey of the .eld. Pattern Recogn Lett 47:18–28 Diop EHS, Angulo J (2012a). Morphological scale-spaces image analysis with robust and adaptive pdes. In: Bene¡s, V, ed. Proc Int Conf S4G Geometry. Prague, Czech Republic. Diop EHS, Angulo J (2012b). Robust nonlinear PDEs for multiscale morphological image analysis. In: Proc 84rd Ann Meet Int Assoc Appl Math (GAMM). Darmstadt, Germany. Dorst L, van den Boomgaard R (1994). Morphological signal processing and the slope transform. Signal Process 38:79–98. Geiger D, Liu TL, Kohn RV (2003). Representation and self-similarity of shapes. IEEE T Pattern Anal 25:86–101. Grazzini J, Soille P (2009). Edge-preserving smoothing using a similarity measure in adaptive geodesic neighbourhoods. Pattern Recogn 42:2306–16. Heijmans HJAM (1994). Morphological image operators. Boston: Academic Press. Heijmans HJAM, Ronse C (1990). The algebraic basis of mathematical morphology I: Dilations and erosions. Comput Vision Graph 50:245–95. Lee TC, Kashyap RL (1994). Building skeleton models via 3-D medial surface/axis thinning algorithms. CVGIP Graph Model Im 56:462–78. Lerallut R, Decenci`ere E, Meyer F (2007). Image .ltering using morphological amoebas. Image Vision Comput 25:395–404. Leymarie F, Levine MD (1992). Simulating the grass.re transform using an active contour model. IEEE T Pattern Anal 14:56–75. Lions PL (1982). Generalized solutions of Hamilton-Jacobi equations. London: Pitman. Lions PL, Souganidis PE, V´asquez JL (1987). The relation between the porous medium and the eikonal equations in several space dimensions. Rev Matem Ibero 3:275– 340. Maragos P (1995). Slope transforms: Theory and application to nonlinear signal processing. IEEE T Signal Proc 43:864–77. Maragos P, Vachier C (2008). A PDE formulation for viscous morphological operators with extensions to intensity-adaptive operators. San Diego, USA. In: Proc 15th IEEE Int Conf Image Proc (ICIP’08). 2200–3. Maragos P, Vachier C (2009). Overview of adaptive morphology: Trends and perspectives. In: Proc 16th IEEE Int Conf Image Proc (ICIP’09). Cairo, Egypt. 2241–4. Matheron G (1975). Random sets and integral geometry. New York: John Wiley & Sons. Ogniewicz RL (1993). Discrete Voronoi skeletons. Konstanz: Hartung-Gorre Verlag. Osher S, Sethian JA (1988). Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations. J Comp Phys 79:12–49. Perona P, Malik J (1990). Scale-space and edge detection using anisotropic diffusion. IEEE T Pattern Anal 12:629–39. Rockafellar RT (1970). Convex analysis. Princeton: Princeton Univ Press. Roerdink JBTM (2009). Adaptivity and group invariance in mathematical morphology. In: Proc 16th IEEE Int Conf Image Proc (ICIP’09). Cairo, Egypt. 2253–6. Rouy E, Tourin A (1992). A viscosity solutions approach to shape-from-shading. SIAM J Numer Anal 29:867–84. Sapiro G, Kimmel R, Shaked D, Kimia BB, Bruckstein AM (1993). Implementing continuous-scale morphology via curve evolution. Pattern Recogn 26:1363–72. Schmitt M (1989). Some examples of algorithm analysis in computational geometry by means of mathematical morphological techniques. In: Boissonnat JD, Laumond JP, eds. Geometry and Robotics. Lect Not Comput Sci 391:225–46. Serra J (1982). Image analysis and mathematical morphology, vol. I. London: Academic Press. Serra J (1988). Image analysis and mathematical morphology: Theoretical advances, vol. II. London: Academic Press. Sheehy DJ, Armstrong CG, Robinson DJ (1996). Shape description by medial surface construction. IEEE Trans Vis and Comp Graph 2:62–72. Siddiqi K, Bouix S, Tannenbaum A, Zucker SW (2002). Hamilton-Jacobi skeletons. Int J Comp Vis 48:215–31. Soille P (1999). Morphological image analysis. Springer-Verlag. Sternberg SR (1986). Grayscale morphology. Comput Vision Graph 35:333–55. Strekalovskiy E, Chambolle A, Cremers D (2012). A convex representation for the vectorial Mumford-Shah functional. In: Proc IEEE Conf Comput Vision Pattern Recogn (CVPR’12). Providence RI, USA. 1712-9. Tek H, Kimia BB (2003). Symmetry maps of free-form curve segments via wave propagation. Int J Comp Vis 54:35–81. Torsello A, Hancock ER (2003). Curvature correction of the Hamilton-Jacobi skeleton. In: Proc IEEE Conf Comput Vision Pattern Recogn (CVPR’03). I-828–34. Torsello A, Hancock ER (2004). A skeletal measure of 2D shape similarity. Comput Vis Image Und 95:1–29. Vachier C, Meyer F (2007). News form viscousland. In: Banon GJF, Barrera J, Braga-Neto UdM, Hirata NST, eds. Proc Int Symp Math Morph. Brazil: Rio de Janeiro. 189–200. van den Boomgaard R, Dorst L (1997). The morphological equivalent of Gaussian scale-space. In: Sporring J, Nielsen M, Florack L, Johansen P, eds. Gaussian scale-space theory. Comput Imaging Vision 8:203–20. van den Boomgaard R, Smeulders A (1992). Towards a morphological scale-space theory. In: O YL, Toet A, Foster D, Heijmans HJAM, Meer P, eds. Shape in Picure. Springer-Verlag. NATO ASI 126:631–40. van den Boomgaard R, Smeulders A (1994). The morphological structures of images: The differential equations of morphological scale-space. IEEE T Pattern Anal 16:1101–13. Velasco-Forero S, Angulo J (2013). On nonlocal mathematical morphology. In: Luengo Hendriks CL, Borgefors G, Strand R, eds. Math Morphol Appl Signal Image Process. Lect Not Comput Sci 7883:219–30. Verd´u-Monedero R, Angulo J, Serra J (2011). Anisotropic morphological .lters with spatially-variant structuring elements based on image-dependent gradient .elds. IEEE T Image Proc 20:200–12.