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

El Hadji Samba Diop, Jesus Angulo

► To cite this version:

El Hadji Samba Diop, Jesus Angulo. MULTISCALE IMAGE ANALYSIS BASED ON ROBUST AND ADAPTIVE MORPHOLOGICAL SCALE-SPACES. Image Analysis & Stereology, 2015, 34 (1), pp.39-59. 10.5566/ias.993. hal-00975728v2

HAL Id: hal-00975728

https://minesparis-psl.hal.science/hal-00975728v2

Submitted on 15 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MULTISCALE IMAGE ANALYSIS BASED ON ROBUST AND ADAPTIVE MORPHOLOGICAL SCALE-SPACES

EL HADJI S. DIOP,† AND JESÚS ANGULO

Center for Mathematical Morphology, Department of Mathematics and Systems, MINES ParisTech, 35 rue Saint-Honoré, 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 defined 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 efficiency 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 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 flow-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 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 $\partial u/\partial t = \Delta u$, by the following model:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = \operatorname{div} (g (\|\nabla u\|^2) \nabla u) , \\ u(x, 0) = f(x) . \end{array} \right. \qquad (1) $$

In this model (Eq. 1), the diffusivity has to be such that $g (\|\nabla u\|^2) \to 0$ when $\|\nabla u\|^2 \to \infty$, and $g (\|\nabla u\|^2) \to 1$ when $\|\nabla u\|^2 \to 0$. One of the diffusivities proposed by Perona and Malik is the following:

$$ g(s^2) = 1 / (1 + (s^2 / \lambda^2)) , \quad \lambda > 0, $$

where $\lambda$ 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 (\|\nabla u\|^2)$ with its regularized version, i.e., $g (\|\nabla u_\sigma\|^2)$, where $\nabla u_\sigma = \nabla (G_\sigma * u)$, and $G_\sigma$ is a Gaussian kernel

<footnote>† Presently with the Departement of Mathematics, University of Thies, Cité Malick SY, BP 967, Thiès, Sénégal</footnote>
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||²),
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), non-local 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 $\mathcal{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 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 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 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.
## 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 $(\epsilon, \delta)$, in such a way that for any given dilation $\delta$, there is a unique erosion $\epsilon$ such that $(\epsilon, \delta)$ 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 $\mathbb{E} \subseteq \mathbb{Z}^2$. For any function $f: \mathbb{E} \to \mathbb{R}$, elementary *dilation* and *erosion* operators (obtained by adjunction and duality) are respectively defined by:

$$ (f \oplus b)(x) := \bigvee_{y \in \mathbb{E}} [f(y) + b(x-y)], \quad (2) $$

$$ (f \ominus b)(x) := \bigwedge_{y \in \mathbb{E}} [f(y) - b(y-x)], \quad (3) $$

where $\bigvee$ (*resp.* $\bigwedge$) represents the supremum (*resp.* infimum). The further convention for scalar addition in $\mathbb{1}$ is considered: $f(y)+b(x-y) = -\infty$ when $b(x-y) = -\infty$ or $f(y) = -\infty$, and that $f(y)-b(y-x) = +\infty$ when $b(y+x) = +\infty$ or $f(y) = -\infty$. 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: \mathbb{R}^2 \to \mathbb{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 $B \subseteq \mathbb{E}$ called structuring element and $-\infty$ outside, i.e., $b(x) = 0$ if $x \in B$ and $b(x) = -\infty$ if $x \in B^c$:

$$ (f \oplus B)(x) = \bigvee_{y \in B} [f(x-y)], \quad (4) $$

$$ (f \ominus B)(x) = \bigwedge_{y \in B} [f(x+y)]. \quad (5) $$

Both dilation and erosion are invariant under translations in $\mathbb{E}$ (spatial or “horizontal” direction) and $\mathbb{R}$ (intensity or “vertical” direction).

**Definition 2.1** Let $\mathcal{F}$ be a family of real functions mapped on $\mathbb{E} \subseteq \mathbb{R}^2$. An operator $S: \mathcal{F} \to \mathcal{F}$ is increasing or monotone if for all $f, g \in \mathcal{F}$ s.t. $f(x) \ge g(x), \forall x \in \mathbb{E}$, then, $(Sf)(x) \ge (Sg)(x), \forall x \in \mathbb{E}$.

In addition to that, *dilation* and *erosion* are increasing operators, and also satisfy the following properties:

- Duality: $(f \oplus b)(x) = -(f \ominus b)(x), \forall x \in \mathbb{E}$.

- Adjunction: $(f \oplus b)(x) \le g(x), \forall x \in \mathbb{E} \iff f(x) \le (g \ominus b)(x), \forall x \in \mathbb{E}$.

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:

$$ (f \circ b) := [(f \ominus b) \oplus b], \quad (6) $$

$$ (f \circled{b}) := [(f \oplus b) \ominus b]. \quad (7) $$

*Opening* and *closing* are also increasing operators with two main properties:

- Ordering (anti-extensivity and extensivity): For any structuring function, one has $(f \circ b) \le f \le (f \circled{b})$

- Idempotence: For any structuring function, one has $[(f \circ b) \circ b] = (f \circ b)$. By duality, one also has: $[(f \circ b) \circled{b}] = (f \circled{b})$.

**Semi-group property, multiscale morphological filters and morphological PDEs.** For $t \ge 0$, one can define then multiscale dilations (and erosions) by replacing $b$ by the family of multiscale structuring functions $(b_t)_{t \ge 0}$ defined for $t > 0$ by:

$$ b_t(x) = \begin{cases} tb(x/t) & \text{for } t > 0 \\ 0 & \text{for } t = 0, x = 0 \\ -\infty & \text{otherwise,} \end{cases} $$

which satisfies the semi-group property $(b_s \oplus 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

$$ b_t(x) = - \frac{\|x\|^2}{2t}. \quad (8) $$

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:

$$ (f \oplus b_t)(x) := \bigvee_{y \in \mathbb{E}} [f(y) + b_t(x-y)], \quad (9) $$

$$ (f \ominus b_t)(x) := \bigwedge_{y \in \mathbb{E}} [f(y) - b_t(y-x)]. \quad (10) $$
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 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 \ge 0}$ of multiscale structuring functions, respectively by:

$$ (f \circ b_t) := [(f \ominus b_t) \oplus b_t], \quad (11) $$

$$ (f \bullet b_t) := [(f \oplus b_t) \ominus b_t]. \quad (12) $$

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):

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = \pm \|\nabla u\|, \\ u(x, 0) = f(x), \end{array} \right. \quad (13) $$

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 \in \mathbb{R}^n : \|x\|_p \le 1\}$ of metrics induced by the $L_p$ norms. There are three special cases of norms $p = 1, 2$ and $p = \infty$ which correspond to particular PDEs: for $B|_1$ (i.e., rhombus), one has $\partial u_t = \pm \|\nabla u\|_\infty$; for $B|_2$ (i.e., disk), one has $\partial u_t = \pm \|\nabla u\|_2$; and for $B|_\infty$ (i.e., square), one has $\partial u_t = \pm \|\nabla u\|_1$. Similarly, the PDE for multiscale parabolic dilations and erosions is given by:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = \pm \frac{1}{2} \|\nabla u\|^2, \\ u(x, 0) = f(x). \end{array} \right. \quad (14) $$

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:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u(x,t)}{\partial t} = \pm H(x, \nabla u(x,t)) = 0, \text{ in } \mathbb{R}^n \times (0, +\infty) \\ u(\cdot, 0) = f \text{ in } \mathbb{R}^n. \end{array} \right. \quad (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:

$$ u(x,t) = \inf_{y \in \mathbb{R}^n} \left[ f(y) + tH^* \left( \frac{x-y}{t} \right) \right], \quad (16) $$

$$ u(x,t) = \sup_{y \in \mathbb{R}^n} \left[ f(y) - tH^* \left( \frac{x-y}{t} \right) \right], \quad (17) $$

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:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = -\operatorname{sign}^+(s-t) \|\nabla u\| + \operatorname{sign}^+(t-s) \|\nabla u\| \\ u(x, 0) = f(x), \end{array} \right. \quad (18) $$

with $t \in [0, 2s]$; $\operatorname{sign}(\cdot)$ stands for the signum function and $r^+ := \max(r, 0)$. Indeed, for $t \in [0, s]$, PDE in Eq. 18 acts as a multiscale erosion, while for $t \in ]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:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = \pm \|\nabla u_\sigma\| \text{ in } \mathbb{R}^N \times (0, T), \\ u(x, 0) = I(x) \text{ in } \mathbb{R}^N \times \{t = 0\}, \end{array} \right. \quad (19) $$
where $u_{\sigma} = u * G_{\sigma}$ is the convolution with a Gaussian $G_{\sigma}$ of variance $\sigma$. 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):

$$ \left\{ \begin{array}{l} \frac{\partial u}{\partial t} = \pm g(\|\nabla u_{\sigma}\|)\|\nabla u\| \text{ in } \mathbb{R}^N \times (0, T), \quad (20) \\ u(x, 0) = I(x) \text{ in } \mathbb{R}^N \times \{t = 0\}, \end{array} \right. $$

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:

$$ g(x) = \begin{cases} 1 + \alpha \|x\|^{2p}, & p > 0 \quad \text{for } g \text{ increasing,} \\ 1 + \alpha \exp(\|x\|^{2p}), & p > 0 \quad \text{for } g \text{ increasing,} \\ \frac{1}{1 + \alpha \|x\|^2} & \text{for } g \text{ decreasing,} \end{cases} \quad (21) $$

where $\|\cdot\|$ stands for the Euclidean norm in $\mathbb{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 here<fnref n="1" /> 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:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} = \pm \|Du\|^p \text{ in } \mathbb{R}^N \times (0, T), \quad (22) \\ u(x, 0) = I(x) \text{ in } \mathbb{R}^N \times \{t = 0\}, \end{array} \right. $$

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:

- flat 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 \to 1$, $p \in ]1, 2[$ and $p \to \infty$, 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):

$$ u(x,t) = \sup_{y \in \mathbb{R}^N} \left\{ I(y) - c_p \frac{\|x-y\|^{p/(p-1)}}{t^{1/(p-1)}} \right\}, \quad (23) $$

where

$$ c_p = \frac{p-1}{p^{p/(p-1)}}. \quad (24) $$

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 Fig. 1.

<figure><img src="image_2.png" /><figcaption>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$ and $p = 2$.</figcaption></figure>

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_{l,p})_{t>0}$ defined for $t > 0$ and $p > 1$ by

$$ k_{l,p}(x) := -c_p \frac{\|x\|^{p/(p-1)}}{t^{1/(p-1)}}. \quad (25) $$

<footnote n="1">¹The idea was first presented in S⁴G conference (Diop and Angulo, 2012a).</footnote>
First of all, notice that for a fixed $t > 0$, when $p \to 1$, then $k_{t,p} \to 0$; when $p \to \infty$, then for all $x$, $k_{t,p}(x) \to \|x\|$. What is particularly interesting is to see the behavior of the family of functions $x \mapsto 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 \in [-1, 1]$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Plottings of structuring functions $k_{t,p}(x)$ for a fixed $t = 1$ and different values of $p$.</figcaption></figure>

Indeed, basically saying, two facts are shown:

- Firstly, as expected, one is dealing with flat structuring functions as $p \to 1$, on the one hand. On the other hand, as $p$ increases up to 2, supports 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 \to \infty$, 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(\nabla u)$ and $1 < p < \infty$, 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, \infty[$ in inhomogeneous image areas, e.g., contours, noise, texture, ... Typical $p = f(\nabla u)$ is given in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Typical $p = f(\nabla u)$ functions for a PDE model based on image features.</figcaption></figure>

To this aim, we propose to set

$$p(u) = 1 + g(\|\nabla u_{\sigma}\|), \qquad (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 significantly increase the robustness against noise. For doing so, we propose to enhance effects of the edge-based parameter $p = p(u) = f(\nabla u)$ by the means of an edge threshold $\alpha$ 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:
- In homogeneous image regions, i.e., $h(\|\nabla u\|) < \alpha$, we set

$$p(u) = p_1 = 1 + g_1(\|\nabla u_{\sigma}\|) \quad (27)$$

in the models of Eqs. 22 and 23.

- In high frequency image areas, i.e., $h(\|\nabla I\|) \ge \alpha$, we set

$$p(u) = p_2 = c + g_2(\|\nabla u_{\sigma}\|) \quad (28)$$

in the models Eqs. 22 and 23, with $c \ge 2$.

For instance, $g_1$ and $g_2$ are decreasing edge-based functions taken as in Eq. 21, and $g_2$ is decreasing more rapidly than $g_1$. Fig. 4 illustrates what we just discussed about the edge value $\alpha$ and typical $p(u) = f(\nabla u)$ function used in this work.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Typical $p(u) = f(\nabla u)$ functions for coupled with image features and threshold edge value $\alpha$.</figcaption></figure>

## NUMERICAL RESULTS

**Discretization** Let us set up the following discretization grid:

$$\begin{array}{l} x_i = i \Delta x, \quad y_j = j \Delta y, \quad t^n = n \Delta t, \\ i, j = 0, \pm 1, \pm 2, \dots, \quad \text{and } n = 0, 1, 2, \dots \end{array} \quad (29)$$

An approximate solution $u$ on the above grid is denoted by $u(x_i, y_j)$, and satisfies:

$$u^n(x_i, y_j) \approx u(i \Delta x, j \Delta y, n \Delta t). \quad (30)$$

**Discretization.** Proposed PDEs in Eqs. 19 and 20 are solved using an explicit scheme:

$$u^{n+1}(x_i, y_j) = u^n(x_i, y_j) \pm \Delta t \cdot F(u^n(x_i, y_j))(Du^n)_{\pm}(x_i, y_j), \quad (31)$$

where $F(u^n(x_i, y_j)) \equiv 1$ and $u^n = u^n$ for Eq. 13, $u^n = u^n \star G_\sigma$ for Eq. 19, and in the case of dealing with Eq. 20, $F(u^n(x_i, y_j))$ corresponds to the discretization of the edge function $g(\|\nabla u_\sigma\|)$. Hence, the (+) (resp. (-)) sign in the term $(Du^n)_\pm$ 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 first consider the following forward and backward derivatives:

$$\begin{align*} D_x^+ u^n(x_i, y_j) &= \frac{u^n(x_i + 1, y_j) - (u^n(x_i, y_j))}{\Delta x}, \\ D_x^- u^n(x_i, y_j) &= \frac{u^n(x_i, y_j) - (u^n(x_i - 1, y_j))}{\Delta x}, \\ D_y^+ u^n(x_i, y_j) &= \frac{u^n(x_i, y_j + 1) - (u^n(x_i, y_j))}{\Delta y}, \\ D_y^- u^n(x_i, y_j) &= \frac{u^n(x_i, y_j) - (u^n(x_i, y_j - 1))}{\Delta y}. \end{align*}$$

Thus, our proposed schemes are the followings:

$$(Du^n)_+ = (\max(\max(0, D_x^+ u^n)^2, \max(0, -D_x^- u^n)^2) + \max(\max(0, D_y^+ u^n)^2, \max(0, -D_y^- u^n)^2))^{1/2}, \quad (32)$$

$$(Du^n)_- = (\max(\min(0, D_x^+ u^n)^2, \min(0, -D_x^- u^n)^2) + \max(\min(0, D_y^+ u^n)^2, \min(0, -D_y^- u^n)^2))^{1/2}. \quad (33)$$

Note that the discretization of $\nabla u_\sigma$ in $g(\|\nabla u_\sigma\|)$ 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.,

$$\|\nabla u_{\sigma}^{n}\|^{2} = (D_{x}^{+} u_{\sigma}^{n})^{2} + (D_{x}^{-} u_{\sigma}^{n})^{2} + (D_{y}^{+} u_{\sigma}^{n})^{2} + (D_{y}^{-} u_{\sigma}^{n})^{2}. \quad (34)$$

<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Noisy Cat image. (b) Owl image.</figcaption></figure>
**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(\nabla u)$ in both the feature-based 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).

<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

<caption>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.</caption>
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.

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 confirm 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).

<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>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.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11. Multiscale morphological operations performed on noisy Owl image given in Fig. 5(b). From top to bottom — First line, multiscale openings at scale 15 using: (a) classical PDE, 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, 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.</figcaption></figure>

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 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.

## 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 filtering 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 filtering. 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). Traffic 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öm 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í J, Benedí JM, Mendonç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 filtering. IEEE T Signal Proc 42:3377-386.

Cao F (1998). Partial differential equations and mathematical morphology. J Math Pure Appl 77:909-41.

Catté F, Dibos F, Koepfler 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é 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.

Ćurić V, Landström A, Thurley MJ, Luengo-Hendriks CL (2014). Adaptive mathematical morphology – a survey of the field. Pattern Recogn Lett 47:18-28

Diop EHS, Angulo J (2012a). Morphological scale-spaces image analysis with robust and adaptive pdes. In: Beneš, 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ère E, Meyer F (2007). Image filtering using morphological amoebas. Image Vision Comput 25:395-404.
Leymarie F, Levine MD (1992). Simulating the grassfire 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ásquez 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 Morph Appl Signal Image Process. Lect Not Comput Sci 7883:219-30.

Verdú-Monedero R, Angulo J, Serra J (2011). Anisotropic morphological filters with spatially-variant structuring elements based on image-dependent gradient fields. *IEEE T Image Proc* 20:200-12.