Jesús Angulo
email: jesus.angulo@mines-paristech.fr

Generalised Morphological Image Diffusion

Keywords: counter-harmonic mean, mathematical morphology, quadratic structuring function, image diffusion, nonlinear scale-space theory, nonlinear partial differential equation

Relationships between linear and morphological scale-spaces have been considered by various previous works. The aim of this paper is to study how to generalise the discrete and continuous diffusion-based approaches in order to introduce nonlinear filters whose limit effects mimic the asymmetric behaviour of morphological dilation and erosion, as well as other multiscale filters, hybrid between the standard linear and morphological filters. A methodology based on the counter-harmonic mean is adopted here. Partial differential equations are formulated and details of numerical implementation are discussed to illustrate the various studied cases: isotropic, edge-preserving and coherence-enhancing diffusion. We also found a new way to derive the classical link between Gaussian scalespace and dilation/erosion scale-spaces based on quadratic structuring functions in the discrete and continuous setting. We have included some preliminary applications of the generalised morphological diffusion to solve image processing problems such as denoising and image enhancement in the case of asymmetric bright/dark image properties.

Introduction

Two fundamental paradigms of image filtering appear distinguished nowadays in the state-ofthe-art. On the one hand, differential methods inspired from the (parabolic) heat equation, including isotropic diffusion [START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF], nonlinear diffusion [START_REF] Perona | Scale-Space and Edge Detection Using Anisotropic Diffusion[END_REF][START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], anisotropic diffusion [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF], etc.

The main properties of these techniques are the appropriateness to deal with the notion of scale-space of image structures and the ability to process symmetrically the bright/dark image structures. Practical algorithms involve (local-adaptive) kernel convolution as well as PDE-formulation and subsequent numerical solutions. The interested reader should refer to basic references [START_REF]Geometry-Driven Diffusion in Computer Vision[END_REF], [START_REF] Florack | Image Structure[END_REF] and [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. On the other hand, mathematical morphology operators [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF] which are formulated in terms of geometric notions as well as in terms of complete lattice theory. Morphological filters entail mainly the computation of supremum and infimum values in neighbourhoods (or structuring elements) which correspond respectively to the dilation and the erosion, the two basic operators. Morphological operators present also good scale-space properties [START_REF] Jackway | Scale-Space Properties of the Multiscale Morphological Dilation-Erosion[END_REF][START_REF] Van Den Boomgaard | The morphological equivalent of Gaussian scale-space[END_REF][START_REF] Schmidt | The morphological equivalents of relativistic and alpha-scalespaces[END_REF], but by the natural duality of complete lattices, most operators appear by pairs and one acts on bright structures and the other one on dark structures. This latter property of asymmetry is in fact an advantage which allows defining evolved operators by composition product of a pair of dual ones. For instance, the opening (resp. closing) is obtained by the product of an erosion (resp. dilation) followed by a dilation (resp. erosion), then the product of openings and closings leads to the alternate filters and other families of morphological filters [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF]. Diffusion involves blurring image structures whereas morphological dilation and erosion involve enhancement of image structure transitions. In fact, morphological operators are related to geometric optics models and in particular to the (hyperbolic) eikonal equation. Hence, there exists also a well motivated formulation of morphological operators using PDEs [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Arehart | Mathematical morphology: The Hamilton-Jacobi connection[END_REF][START_REF] Brockett | Evolution equations for continuous-scale morphology[END_REF][START_REF] Maragos | Differential morphology and image processing[END_REF]. This differential or continuous-scale morphology can be solved using numerical algorithms for curve evolution [START_REF] Sapiro | Implementing continuous-scale morphology via curve evolution[END_REF]. Thus multiscale flat dilation/erosion by disks as structuring elements (resp. unflat dilation/erosion by parabolic structuring functions) can be modelled in a continuous framework. Morphological operators using geometric structuring elements are today one of the most successful areas in image processing. However, it is obvious that the soundness and maturity of numerical methods to implement the different versions of image diffusion constitute an advantage against continuous morphology implementation, which requires more specific numerical schemes to achieve robust results [START_REF] Van Den Boomgaard | Numerical solution schemes for continuous-scale morphology[END_REF][START_REF] Breuß | Highly accurate PDE-based morphology for general structuring elements[END_REF].

Motivation and Aim: Asymmetrisation of diffusion image filtering.

Let us consider the two images given in Fig. 1. The first image is an example of impulse-noise corrupted image, where the noise presents a mean value which is significatively smaller than the mean image value. In other words, the noise is asymmetrically shifted towards dark values. A successful denoising approach, based for instance on nonlinear adaptive diffusion, should be able to deal with this asymmetry. The second image is a low SNR retinal image and one of its typical applications involves to extract the vessels. The latter structures are dark with respect to their background. Vessel detection can be solved using morphological operators, but the result would be poor without a prior enhancement of the image using typically a coherenceenhancing diffusion-based filtering step. Obviously, the relation of intensities between the vessels and the background could be of help to adapt the nonlinearity of the appropriate In this context, the aim of this paper is to study how to generalise the diffusion-based approaches in order to introduce a family of nonlinear filters whose effects mimic the asymmetric behaviour of morphological dilation and erosion, as well as other evolved adaptive filters for edge-preserving and coherence-enhancement. A parameter P of nonlinearity is used to tune the asymmetrization (according to the sign of P ) and limit effects (according to the modulus of P ). Our goal is therefore to propose a new approach of generalised morphological image diffusion that expands and includes the Gaussian scale-space and the morphological erosion/dilation scale-spaces and which is implemented by means of standard (well possed) numerical algorithms.

Related work.

Theoretical investigation of the relation between the scale-space concepts of linear and morphological scale spaces has been carried out by various studies. A parallelism between the role which is played by the Fourier transform in the convolution with Gaussian kernels and the dilation using quadratic structuring function has been established using the notion of slope transform. More generally, the slope transform and its applications to morphological systems was developed independently and simultaneously by [START_REF] Dorst | Morphological Signal Processing and the Slope Transform[END_REF] and [START_REF] Maragos | Slope Transforms: Theory and Application to Nonlinear Signal Processing[END_REF]. Early work by van den Boomgaard [START_REF] Van De Boomgaard | The morphological equivalent of Gauss convolution[END_REF][START_REF] Van Den Boomgaard | The morphological equivalent of Gaussian scale-space[END_REF] on the morphological equivalent of Gaussian convolution and Gaussian scale-space should be also mentioned. In [START_REF] Burgeth | An Explanation for the Logarithmic Connection between Linear and Morphological System Theory[END_REF] it was shown that the slope transform in the (max, +)-algebra corresponds to the logarithmic multivariate Laplace transform in the (+, •)-algebra; and that the Cramer transform is the Legendre-Fenchel transform of the logarithmic Laplace transform. Following a similar line, but introducing a modified Cramer-Fourier transform, the recent work [START_REF] Schmidt | The morphological equivalents of relativistic and alpha-scalespaces[END_REF] has explored a set of connections between morphological evolution processes and relativistic scale-spaces and alpha-scale-spaces. The structural similarities between linear and morphological processes were studied in [START_REF] Florack | Pseudo-Linear Scale-Space Theory[END_REF] in order to construct a one-parameter semi-linear process that incorporates Gaussian scale-space, and both types of morphological scale-spaces by quadratic structuring elements as limiting processes of a one-parameter transformation of grey-level values. A morphological scale-space by deforming the algebraic operations related to Minkowski (or L p ) norms and generalised means was proposed in [START_REF] Welk | Families of generalised morphological scale spaces[END_REF], where the Gaussian scale-space is a limit case. We adopted here a different methodology in order to link diffusion-based image filtering and morphological image filtering. The starting point of our approach is the notion of counter-harmonic mean [START_REF] Bullen | Handbook of Means and Their Inequalities[END_REF]. In fact, the idea of using the counter-harmonic mean for constructing robust morphologicallike operators, without the notions of supremum and imfimum, was proposed in [START_REF] Van Vliet | Robust Local Max-Min Filters by Normalized Power-Weighted Filtering[END_REF]. Our purpose in this paper is to go further and to exploit the counter-harmonic mean to propose a more general framework which can be exploited for the discrete/continuous paradigms and algorithms of image diffusion.

Reader interested on a generalised axiomatic formulation of Gaussian scale-spaces beyond the rotationally symmetric Gaussian scale-space is referred to the modern approach developed in [START_REF] Lindeberg | Generalized Gaussian scale-space axiomatics comprising linear scale-space, affine scale-space and spatio-temporal scale-space[END_REF].

Paper organisation. This paper is an extended and improved version of a conference contribution [START_REF] Angulo | Pseudo-Morphological Image Diffusion using the Counter-Harmonic Paradigm[END_REF]. In particular, in the present manuscript, a complete proof of the different results and properties is given, as well as the formulation of the partial differential equations associated to the new filters. More illustrative examples included and potential applications in low level vision are explored.

The outline of the paper is as follows. In the next section we review the notion of counter-harmonic mean. The appropriateness of counter-harmonic mean to approximate flat dilation/erosion is considered in Section 3. Section 4 introduces the novel counter-harmonic Gaussian scale-space. The limit relationships with parabolic dilation/erosion as well as other theoretical properties are also considered in Section 4. Then, Section 5 proposes a continuous formulation of the counter-harmonic Gaussian scale-space by means of a PDE named counter-harmonic isotropic diffusion. Numeral schemes and approximation are discussed. Section 6 extends these investigations to the counter-harmonic locally adaptive diffusion, in particular giving the counterparts of the Perona and Malik model and to the Weickert model of coherence-enhanced diffusion. Some preliminary applications of the generalised morphological diffusion to solve image processing problems such as denoising and image enhancement in the case of asymmetric bright/dark image properties are discussed in Section 7. The paper is concluded with a summary and perspectives in Section 8.

Definition 1 Let a = (a 1 , a 2 , • • • , a n ) and w = (w 1 , w 2 , • • • , w n ) be non-negative real n-tuples, i.e., a, w ∈ R n + .
If r ∈ R then the r-th counter-harmonic mean (CHM) of a with weight w is given by [START_REF] Bullen | Handbook of Means and Their Inequalities[END_REF]

K [r] (a; w) =            ∑ n i=1 w i a r i ∑ n i=1 w i a r-1 i if r ∈ R max(a i ) if r = +∞ min(a i ) if r = -∞ (1)
The equal weight case will be denoted by K [r] (a). We notice that K [1] (a; w) is the weighted arithmetic mean A(a; w) and K [0] (a; w) is the weighted harmonic mean H(a; w). Remark also that

K [1/2] (a 1 , a 2 ) = G(a 1 , a 2 )
, where G(a; w) denotes the weighted geometric mean; however this result is only valid for a 2-tuple. It is well known that these classical means are ordered between them; i.e., H(a; w) ≤ G(a; w) ≤ A(a; w).

The following properties are useful for the rest of the paper.

Proposition 2 If a and w are non-negative real n-tuples, we have the following property:

K [r] (a; w) = ( K [-r+1] (a -1 ; w) ) -1
Proof. We have

K [r] (a -1 ; w) = ∑ n i=1 w i a -r i ∑ n i=1 w i a -r+1 i = ( ∑ n i=1 w i a -r+1 i ∑ n i=1 w i a -r i ) -1 . Hence K [r] (a -1 ; w) = ( K [-r+1] (a; w)
) -1 . By taking now s = -r + 1, so r = -s + 1, we obtain the desired result r] (a; w); and if -∞ ≤ r ≤ 1 then the following stronger results holds:

K [s] (a; w) = ( K [-s+1] (a -1 ; w) ) -1 . Proposition 3 If 1 ≤ r ≤ +∞ then K [r] (a; w) ≥ M [
K [r] (a; w) ≤ M [r-1] (a; w); where M [r] (a; w) = ( 1 W ∑ n i=1 w i a r i ) 1/r
is the r-th power-mean, or Minkowski weighted mean of order r, defined for r ∈ R * . Inequalities are strict unless r = 1, +∞, -∞ or a is constant.

Proof. Assume that r ∈ R, we can rewrite K [r] (a; w) = (( ∑ n i=1 w i a r i ) 1/r ) r (( ∑ n i=1 w i a r-1 i ) 1/(r-1) ) (r-1) = ( ∑ n i=1 w i a r i ) 1/r ( ( ∑ n i=1 w i a r i ) 1/r ( ∑ n i=1 w i a r-1 i ) 1/(r-1)
) r-1

. Consequently, we have

K [r] (a; w) = M [r] (a; w) ( M [r] (a; w) M [r-1] (a; w) ) r-1
Considering r ≥ 1 and taking into account that M [r] (a; w) ≥ M [r-1] (a; w), we have K [r] (a; w) ≥ M [r] (a; w). Now, assume that r ≤ 1, then using Proposition 2, we have

K [r] (a; w) = ( K [-r+1] (a -1 ; w) ) -1 ≤ ( M [-r+1] (a -1 ; w) ) -1 = M [r-1] (a; w)
The inequalities justify the cases r = ±∞.

Proposition 4 If a and w are n-tuples and if -∞ ≤ r < s ≤ +∞ then K [r] (a; w) ≤ K [s] (a; w), with equality if and only if a is constant.

Proof. Using the following result (see [START_REF] Bullen | Handbook of Means and Their Inequalities[END_REF], Chapter 3): the function m(r) = ∑ n i=1 w i a r i is strictly log-convex in r if a is not constant. Now if we write log ( K [r] (a; w) ) = log (m(r))log (m(r -1)), and if -∞ < r < s < +∞, we easily obtain by strict convexity log (m(r))log (m(r -1)) < log (m(s)) -log (m(s -1)). If either r = -∞ or s = +∞ the result is immediate.

Robust Pseudo-Morphological Operators using CHM

The CHM has been considered in the image processing literature as a suitable filter to deal with salt and pepper noise [START_REF] Gonzalez | Digital Image Processing[END_REF]. More precisely, let v = f (x, y) be a grey-level image: f : Ω → V. Typically, for digital 2D images, (x, y) ∈ Ω where Ω ⊂ Z 2 is the discrete support of the image. The pixel values are v ∈ V ⊂ Z or R, but for the sake of simplicity of our study, we consider a non-negative normalized valued image; i.e., V = [0, 1].

Definition 5

The CHM filter is obtained as

κ P B (f )(x, y) =

∑

(s,t)∈B(x,y) f (s, t) P +1 ∑ (s,t)∈B(x,y) f (s, t) P = K [P +1] ({f (s, t)} (s,t)∈B(x,y) )

(2)

where B(x, y) is the window of the filter, centered at point (x, y), i.e., the structuring element in the case of morphological operators.

This filter is well suited for reducing the effect of pepper noise for P > 0 and of salt noise for P < 0. In the pioneering paper [START_REF] Van Vliet | Robust Local Max-Min Filters by Normalized Power-Weighted Filtering[END_REF], starting from the natural observation that morphological dilation and erosion are the limit cases of the CHM, i.e., lim

P →+∞ κ P B (f )(x, y) = max (s,t)∈B(x,y) (f (s, t)) = δ B (f )(x, y) (3)
and lim

P →-∞ κ P B (f )(x, y) = min (s,t)∈B(x,y) (f (s, t)) = ε B (f )(x, y), (4)
the CHM was used to calculate robust nonlinear operators which approach the morphological ones but without using max and min operators. In addition, these operators are more robust to outliers (i.e., to noise) and consequently they can be considered as an alternative to rankbased filters in the implementation of pseudo-morphological operators.

Comparison with Minkowski Power Means

It is easy to see that for P ≫ 0 (P ≪ 0) the pixels with largest (smallest) values in the local neighbourhood B will dominate the result of the weighted sum. Of course, in practice, the range of P is limited due to the precision in the computation of the floating point operations. Proposition 3 of the previous Section justifies theoretically the suitability of CHM with respect to the alternative approach by high-order Minkowski mean, as considered by Welk [START_REF] Welk | Families of generalised morphological scale spaces[END_REF]. Let us illustrate empirically how both means converge to the supremum (resp. infimum) when positive P increases (negative P decreases). Fig. 2 depicts convergence with respect to the value of P for the erosion (in blue) and dilation (in red), using Minkowski mean in (a) and using CHM in (b). The convergence is measured as the average difference value between the CHM for each P and the exact dilation/erosion obtained by max and min. The practical advantage of CHM to approach morphological operators is obvious: even for P = 100 (or P = -100) the dilation (resp. erosion) is not reached for Minkowski mean whereas the error in the results for CHM is already negligible for P = 20 (resp. P = -20). It should be noted that powers P = ±100 can only be computed as double-precision numbers (encoded as 64-bit floating-point number) We notice also in the empirical curves, as expected from Proposition 3, that the convergence to the erosion with P ≪ 0 is faster than to the dilation with equivalent P ≫ 0, i.e., for P > 0

|κ P B (f )(x, y) -δ B (f )(x, y)| ≥ |κ -P B (f )(x, y) -ε B (f )(x, y)|, ∀(x, y) ∈ Ω; ∀B

Symmetry, (Self )-Duality and other Properties

The asymmetry of P vs. -P involves that κ P B (f ) and κ -P B (f ) are not dual operators with respect to the involution associated to the usual complement (or negative), i.e.,

κ P B (f ) ̸ = n κ -P B ( n f ),
with n f = 1 -f . However, using the Proposition 2, we have the following duality

κ P B (f ) = i κ -P -1 B ( i f ), (5)
where the involution is now associated to the reciprocal (or inversion), i.e., i f = 1/f (equivalent to the complement of the logarithms). In addition, the case P = -1/2 is particularly relevant since we have the following result

κ -1/2 B (f ) = ( κ -1/2 B (f -1 ) ) -1 ,
which becomes a self-dual operator [START_REF] Heijmans | Inf-Semilattice Approach to Self-Dual Morphology[END_REF]. We note that κ

-1/2 B
(f ) is a kind of middle (or "neutral") operator between the arithmetic mean and the harmonic mean.

Symmetry associated to the duality of dilation and erosion is a basic property of classical morphology since their combination leads to the opening and closing [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. Fig. 2(c) gives rates of convergence of CHM to the opening and closing approximated by the product of κ P B (f ) and κ -P B (f ) and Fig. 2(d) using the CHM operators κ P B (f ) and κ -P -1 B (f ). As expected, we observe that the second case leads to more asymmetric results of pseudo-opening/pseudoclosing. This is due to the facts that convergence of κ -P B (f ) is already faster than convergence of κ P B (f ) and that κ -P -1 B (f ) ≤ κ -P B (f ). We can conclude that, even if the pair κ P B (f ) and κ -P -1 B (f ) is dual with respect to the reciprocal, that does not imply that their products are a better choice than κ P B (f ) and κ -P B (f ) for evolved operators such as the opening/closing. As it was already pointed out in [START_REF] Van Vliet | Robust Local Max-Min Filters by Normalized Power-Weighted Filtering[END_REF], another drawback of κ P B (f ) (resp. κ -P B (f )) is the fact that f (x, y) κ P B (f )(x, y) with P > 0 (resp. f (x, y) κ -P B (f )(x, y) with P < 0). Or in other words, the extensitivity (resp. anti-extensitivity) for P > 0 (resp. P < 0) is not guaranteed. However, according to proposition 4, the following ordering relationship holds for P > 0: using counter-harmonic mean; (c) pseudo-opening and pseudo-closing using counter-harmonic mean defined respectively as the combinations κ P B (κ -P B ) and κ -P B (κ P B ); (d) pseudo-opening and pseudo-closing using counterharmonic mean defined respectively as the combinations κ P B (κ -P -1

κ -P B (f )(x, y) ≤ κ P B (f )(x, y).

B

) and

κ -P -1 B (κ P B ).

Counter-Harmonic Gaussian Scale-Space

Canonic multiscale image analysis involves obtaining the multiscale linear convolutions of the original image:

ψ(f )(x, y; t) = (f * K σ )(x, y) = ∫ Ω f (u, v)K σ (x -u, y -v)dudv, ( 6
)
where

K σ is the two-dimensional Gaussian function K σ (x, y) = 1 2πσ 2 exp ( -(x 2 +y 2 ) 2σ 2
) , whose variance (or width) σ 2 is proportional to the scale t; i.e., σ 2 = 2t. Larger values of t lead to image representations at coarser levels of scale.

We can define according to our paradigm the following generalised morphological scalespace.

Definition 6

The counter-harmonic Gaussian scale-space of order P is defined as the following transform parameterized by scale t η(f )(x, y; t; P ) = (f

P +1 * K √ 2t )(x, y) (f P * K √ 2t )(x, y) = ∫ E f (u, v) P +1 K √ 2t (x -u, y -v)dudv ∫ E f (u, v) P K √ 2t (x -u, y -v)dudv . ( 7
)
By choosing P > 0 (resp. P < 0), η(f )(x, y; t; P ) leads to a scale-space of pseudo-dilations (resp. pseudo-erosions), whose filtering effects for a given scale t depend on the "nonlinearity order" of P , which skew the Gaussian weighted values towards the supremum or infimum value.

Since the image is stored as a collection of discrete pixels we need to produce a discrete approximation to the Gaussian kernel. An infinitely large Gaussian kernel cannot be applied in practice, hence one has to truncate it to a finite window as small as possible to obtain a fast transformation. It is well known that the Gaussian kernel is effectively zero more than about three standard deviations from the center, and so one can truncate the kernel at this point. We have adopted for the examples given in this paper a support window of size 2t × 2t, or equivalently σ 2 × σ 2 . The rationale behind this choice is to make easier the comparison with the limit cases P = +∞ and P = -∞, which are respectively a dilation and an erosion using a flat square of size 2t × 2t as structuring element.

It should be noted that there is a parallel genuine scale-space theory for discrete signals introduced by Lindeberg [START_REF] Lindeberg | Scale-Space for Discrete Signals[END_REF][START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF], corresponding to semi-discrete diffusion equations and convolution with a discrete analogue of the Gaussian kernel. By replacing (6) by its discrete analogue, a counter-harmonic Gaussian scale-space can be also naturally considered.

Fig. 3 depicts a comparative example of image filtering using the CHM Gaussian scalespace for a fixed t and variable P . The terminology of "pseudo-dilation", or more generally "pseudo-morphological operators" is probably inappropriate: so many operators have been published as pseudo-morphological ones, without relation with the present approach. However, due to the fact that the nonlinear filter η(f )(x, y; t; P ) for P > 0 (resp. P < 0) is not extensive (resp. anti-extensive) and does not commute with the supremum (resp. with the infimum), it cannot be considered stricto sensu as a dilation (resp. erosion). We notice also that, for a fixed value of the nonlinear parameter P , the associated CHM Gaussian scale-space for a variable scale t leads to families of multi-scale pseudo-dilations/erosions.

Limit statements

We know that η(f )(x, y; t;

P = +∞) = δ B (f )(x, y) and η(f )(x, y; t; P = -∞) = ε B (f )(x, y),
i.e., flat dilation and flat erosion, where B is the square support of the Gaussian kernel. Let us consider in detail the limit cases for P and P .

Proposition 7 For a given scale parameter t, the limits of η(f )(x, y; t; P ) with respect to P exist and are given by

lim P →+∞ η(f )(x, y; t; P ) ≈ sup (u,v)∈Ω ( f (x -u, y -v) - (u 2 + v 2 ) 2P (2t) ) , ( 8
) lim P →-∞ η(f )(x, y; t; P ) ≈ inf (u,v)∈Ω ( f (x -u, y -v) + (u 2 + v 2 ) 2P (2t) ) , ( 9
)
which can be interpreted respectively as the dilation and the erosion of f (x, y) with a quadratic structuring function

b √ 2t (x, y; P ) = (x 2 + y 2 ) 2P (2t) ,
i.e., the numerical dilation and erosion [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF] are defined by

(f ⊕ b √ 2t )(x, y) = sup (u,v)∈Ω ( f (x -u, y -v) -b √ 2t (u, v) ) , (f ⊖ b √ 2t )(x, y) = inf (u,v)∈Ω ( f (x -u, y -v) + b √ 2t (u, v)
) .

Note that in these limiting cases the CHM framework involves a "normalization" by P of the original Gaussian kernel scale parameter during unlinearization, i.e., the nonlinear asymptotic scale parameter is σ = √ P σ = √ P (2t). This result is perfectly coherent with those obtained from totally different paradigms [START_REF] Welk | Families of generalised morphological scale spaces[END_REF][START_REF] Florack | Pseudo-Linear Scale-Space Theory[END_REF]. Also, we notice again that for P = +∞ the structuring function becomes flat and hence we obtain the flat dilation.

Proof. We only give the proof for the case P → +∞ since the case P → -∞ is solved by a similar technique. By rewriting f P = exp(P log(f )), taking first order Taylor expansion log(f ) ≈ f -1 and first order Taylor expansion of exponential function such that we have: lim

N D = exp ( log ( N D )) ≈ 1 + log(N ) -log(D),
P →+∞ η(f )(x, y; t; P ) = 1 + log ∫ Ω exp((P + 1)[f (x -u, y -v) - (u 2 + v 2 ) 2(P + 1)(2t) -1])dudv -log ∫ Ω exp(P [f (x -u, y -v) - (u 2 + v 2 ) 2P (2t) -1])dudv,
which can be rewritten as 1 + (P + 1) log

( ∫ Ω ( exp(f (x -u, y -v) - (u 2 + v 2 ) 2(P + 1)(2t) -1
)

) (P +1) dudv ) 1 (P +1) -P log ( ∫ Ω ( exp(f (x -u, y -v) - (u 2 + v 2 ) 2P (2t) -1)
) P dudv

) 1 P ,
Using now the standard result lim

P →+∞ {∫ Ω g P (x)dx } 1/P = sup x∈Ω g(x),
which holds for positive and bounded function g with support space Ω, and and considering continuity and monotonicity of the logarithm, we obtain: lim P →+∞ η(f )(x, y; t; P ) = 1 + (P + 1) sup

(u,v)∈Ω (f (x -u, y -v) - (u 2 + v 2 ) 2(P + 1)(2t) -1) - P sup (u,v)∈Ω (f (x -u, y -v) - (u 2 + v 2 ) 2P (2t) -1)
.

By considering that both sup-convolutions give approximately the same result, i.e., sup

(u,v)∈Ω (f (x -u, y -v) - (u 2 + v 2 ) 2P (2t) -1) ≈ sup (u,v)∈Ω (f (x -u, y -v) - (u 2 + v 2 ) 2(P + 1)(2t) -1),
we finally obtain the corresponding result.

Properties of CHM Gaussian Scale-Space

We discuss now in more detail several features of the t × P scale-space η(f )(x, y; t; P ).

Continuity

The generic kernel function

K √ 2t (x, y) is a continuous function in R 2 .

Proposition 8

The scale-space image η(f )(x, y; t; P ) is continuous for all (x, y) ∈ Ω, for any t ≥ 0 and for any P ∈ R.

Proof. Consider K √ 2(t+∆t) (x + ∆x, y + ∆y) = K √ 2t (x, y) + E(x, y)
where assuming the continuity of K(x, y), which implies the continuity of K √ 2t (x, y), we have

lim |∆x| → 0 |∆y| → 0 ∆t → 0 E(x, y) = 0 for all (x, y) ∈ Ω ( 10
)
For t > 0, we have

η(f )(x + ∆x, y + ∆y; t + ∆t; P ) = ∫ E f (x -u + ∆x, y -v + ∆y) P +1 K √ 2(t+∆t) (u, v)dudv ∫ E f (x -u + ∆x, y -v + ∆y) P K √ 2(t+∆t) (u, v)dudv
Setting w = u -∆x and z = v -∆y we have

η(f )(x + ∆x, y + ∆y; t + ∆t; P ) = ∫ E f (x -w, y -z) P +1 K √ 2(t+∆t) (w + ∆x, z + ∆y)dwdz ∫ E f (x -w, y -z) P K √ 2(t+∆t) (w + ∆x, z + ∆y)dwdz = ∫ E f (x -w, y -z) P +1 ( K √ 2t (w, z) + E(w, z) ) dwdz ∫ E f (x -w, y -z) P ( K √ 2t (w, z) + E(w, z)
) dwdz

So, in the limit as |∆x| → 0, |∆y| → 0, ∆t → 0 and assuming that the expression (10) holds, hence η(f )(x + ∆x, y + ∆y; t + ∆t; P ) → η(f )(x, y; t; P ) establishing the continuity.

Ordering

We have the following property of ordering w.r.t. scale parameter t.

Proposition 9 There exists a value of P such that for any Q ≥ P > 0, we have that for any pair of ordered scales

0 < t 1 < t 2 involves η(f )(x, y; t 1 ; Q) ≤ η(f )(x, y; t 2 ; Q). Similarly, there is a P such that for R ≤ P < 0, if 0 < t 1 < t 2 then η(f )(x, y; t 1 ; R) ≥ η(f )(x, y; t 2 ; R).
Proof. If we consider that for a large enough P > 0, with Q ≥ P we have

η(f )(x, y; t 1 ; Q) = sup (u,v)∈Ω ( f (x -u, y -v) - (u 2 + v 2 ) 4Qt 1
) and we have also for 0 < t 1 < t 2 :

(u 2 +v 2 ) 4Qt 2 < (u 2 +v 2 ) 4Qt 1 , ∀(u, v) ∈ D where D ⊂ Ω. Then we obtain sup (u,v)∈Ω ( f (x -u, y -v) - (u 2 + v 2 ) 4Qt 1 ) ≤ sup (u,v)∈Ω ( f (x -u, y -v) - (u 2 + v 2 ) 4Qt 2 )
thus proving the proposition. Similarly, the proof can be obtained for P < 0 Hence, the ordering properties associated to the spatial scale t appears only in the pseudomorphological behaviour of the CHM for large enough |P |. The fundamental ordering relationship w.r.t. P is directly inherited from Proposition 3 of the CHM.

Proposition 10 For any pair of nonlinearity orders -∞ ≤ R < S ≤ +∞, we have

η(f )(x, y; t; R) ≤ η(f )(x, y; t; S)
for any t ≥ 0, with equality for t = 0.

Homogeneity and dimensional functionals

Homogenous dimensionality is a fundamental concept in measurements. This principle was extended in [START_REF] Rivest | Dimensionality in Image Analysis[END_REF] to image processing, in particular, in the framework of mathematical morphology, as a way to take into account the homogeneity of physical properties computed from images. Later, the case of nonflat scale-space mathematical morphology using elliptic poweroid structuring functions was studied in [START_REF] Jackway | On Dimensionality in Multiscale Morphological Scale-Space with Elliptic Poweroid Structuring Functions[END_REF]. Making a measurement on an image consists in applying a functional on the image, where a "functional" is a global parameter associated with a function. More precisely, a functional W on a function f is said to be dimensional if and only if there exist constants k 1 , k 2 such that for all λ 1 , λ 2 > 0

W (f ′ ) = λ k 1 1 λ k 2 2 W (f )
where f ′ (x, y) = λ 1 f (λ 2 x, λ 2 y) and where k 1 is the intensity dimension and k 2 the space dimension. Hence the parameter λ 1 accounts for the affinity along the gray tone axis and λ 2 for the homothety of the support space of image f . The relation W (f ′ ) restricts the way in which affinities and homotheties of an image affect dimensional measurements on it and results in a decoupling between affinity and homothety measurements. We notice that if k 1 = 0 the functional W is invariant under intensity affinities and if k 2 = 0 involves invariance under homotheties. Let us consider under which conditions a dimensional functional of the CHM Gaussian scale-space is also a dimensional functional of the underlying image.

Let us define the transformed scale-space:

λ 1 η(f )(λ 2 x,
W ( η(f ′ )(x, y; t; P ) ) = λ k 1 1 λ k 2 2 W (η(f )(x, y; t; P ))
Proof. We start by writing

λ 1 η(f )(λ 2 x, λ 2 y; λ 3 t; λ 4 P ) = λ 1 ∫ E f (λ 2 (x -u), λ 2 (y -v)) λ 4 P +1 K √ 2λ 3 t (λ 2 u, λ 2 v)dudv ∫ E f (λ 2 (x -u), λ 2 (y -v)) λ 4 P K √ 2λ 3 t (λ 2 u, λ 2 v)dudv .
On the one hand, we have

K √ 2λ 3 t (λ 2 x, λ 2 y) = 1 4πλ 3 t exp ( -(x 2 +y 2 ) 4t λ 2 2 λ 3
) , on the other hand, we

define f ′ (x, y) = λ 1 f (λ 2 (x, y)). In addition, λ 2 1 f λ 4 P +1 (λ 2 (x, y)) = λ 1 f λ 4 P (λ 2 (x, y))λ 1 f (λ 2 (x, y)) = (f ′ ) λ 4 P (x, y)f ′ (x, y) = (f ′ ) λ 4 P +1 (x, y). Thus λ 1 η(f )(λ 2 x, λ 2 y; λ 3 t; λ 4 P ) = ∫ E f ′ (x -u, y -v) λ 4 P +1 e -(u 2 +v 2 ) 4t λ 2 2 λ 3 dudv ∫ E f ′ (x -u, y -v) λ 4 P e -(u 2 +v 2 ) 4t λ 2 2 λ 3 dudv
Consequently if we fix λ 2 2 = λ 3 , for the case P = 0 we have

λ 1 η(f )(λ 2 x, λ 2 y; λ 3 t; 0) = η(f ′ )(x, y; t; 0)
and for the case P → +∞ and fixing λ 2 2 = λ 3 λ 4 we have

λ 1 η(f )(λ 2 x, λ 2 y; λ 3 t; λ 4 P ) = η(f ′ )(x, y; t; P )
In conclusion, a functional W will be dimensional only when a modification of homothety is compensated by an affinity on the spatial scale, which involves also to increase the nonlinearity order. This point can also be studied under the viewpoint of the Heijmans theory of [START_REF] Heijmans | Algebraic Framework for linear and morphological scale-spaces[END_REF][START_REF] Heijmans | Morphological scale-spaces, scale invariance, and Lie groups[END_REF].

Invariance

By its linearity, Gaussian scale-space is invariant under linear transformations

f → f = αf +β ⇒ ( f * K √ 2t ) (x, y) = α ( f * K √ 2t ) (x, y) + β. Morphological dilation/erosion are invariant only under grey-value shifts f → f = f + β ⇒ ( f ⊕ b √ 2t ) (x, y) = ( f ⊕ b √ 2t ) (x, y) + β.
We notice that the flat dilation/erosion are invariant under anamorphosis [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Matheron | Les treillis compacts[END_REF], i.e., any strictly increasing mapping is an anamorphosis. Let us look now at invariance properties of the CHM Gaussian scale-space.

Proposition 12

The CHM Gaussian scale-space η(f )(x, y; t; P ) is invariant under scalar multiplication of grey-values for any P , -∞ ≤ P ≤ +∞, i.e., η(αf )(x, y; t; P ) = αη(f )(x, y; t; P ). In the limit cases P = 0 (equivalent to Gaussian scale-space) and P → ±∞, the CHM Gaussian scale-space is invariant under linear transformations: η(αf + β)(x, y; t; P = 0) = αη(f )(x, y; t; P = 0) + β and lim P →±∞ η(αf + β)(x, y; t; P ) = α lim P →±∞ η(f )(x, y; t; P ) +β.

Proof. The proofs are direct from the expressions of the limit cases.

Hence the CHM Gaussian scale-space is invariant to grey-level addition and multiplication asymptotically to the limit cases.

Asymptotic separability and semi-group law

In the classical state-of-the-art on linear scale-spaces, the fundamental semi-group property [START_REF] Lindeberg | Scale-Space for Discrete Signals[END_REF][START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF], or recursivity principle [START_REF] Pauwels | An Extended Class of Scale-Invariant and Recursive Scale Space Filters[END_REF], states as follows: ∀t 1 , t 2 ≥ 0,

( (f * K √ 2t 1 ) * K √ 2t 2 ) = ( f * K √ 2(t 1 +t 2 )
)

. This is due to the fact that the Gaussian convolution kernels are the only separable and rotationally invariant functions that preserve the shape under Fourier transform. That allows also a separability in the implementation of the 2D Gaussian kernel convolution at scale t, using two successive 1D Gaussian kernel convolutions (at directions x and y) at scale t. From the works by van den Boomgaard [START_REF] Van De Boomgaard | The morphological equivalent of Gauss convolution[END_REF][START_REF] Van Den Boomgaard | The morphological equivalent of Gaussian scale-space[END_REF], it is also well known that in mathematical morphology the parabolic structuring function (in fact, any quadratic structuring function) are the equivalent class of functions which are dimensionally separable and closed with respect to the dilation/erosion, i.e., ∀t 1 , t 2 ≥ 0,

( (f ⊕ b √ 2t 1 ) ⊕ b √ 2t 2 ) = ( f ⊕ b √ 2(t 1 +t 2 )
)

, and similarly for the erosion. This result is usually proved in the slope transform domain [START_REF] Dorst | Morphological Signal Processing and the Slope Transform[END_REF][START_REF] Maragos | Slope Transforms: Theory and Application to Nonlinear Signal Processing[END_REF].

The filter family η(f )(x, y; t; P ) for a given P ∈ R has no semi-group structure and is, therefore, not a scale-space in strict sense. However, using the limit expressions it is easy to prove that for the asymptotic cases P = 0 and P → ±∞, the following semi-group property holds:

η (η(f )(x, y; t 1 ; P )) (x, y; t 2 ; P ) = η(f )(x, y; t 1 + t 2 ; P ); ∀t 1 , t 2 ≥ 0 This result constitutes an advantage for an efficient numerical computation of an acceptable approximation to the limit cases of 2D CHM Gaussian scale-space using 1D separable Gaussian kernel for the convolutions of f P +1 and of f P . Nevertheless, since computational efficiency is often important, the result of the operator η(f )(x, y; t; P ) for any P can be computed using low-order recursive filters for the numerator and denominator 2D Gaussian convolutions of the CHM filter. For example, the approach introduced in [57] uses a thirdorder recursive filter with one real pole and a pair of complex poles, applied forward and backward to make a sixth-order symmetric approximation to the Gaussian with low computational complexity for any scale t.

Counter-harmonic Diffusion

A continuous formulation of the counter-harmonic Gaussian scale-space by means of a diffusion PDE is introduced in this section, including a study of the numerical scheme and approximation.

Linear (isotropic) diffusion and its numerical solution

A filtered image u(x, y; t) of f (x, y) is calculated by solving the diffusion equation with the original image as initial state, and reflecting boundary conditions:

∂ t u = div (c∇u) = c∆u = c ( ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ) (11) u(x, y; 0) = f (x, y) ∂ n u| ∂Ω = 0
where c is the conductivity and n denotes the normal to the image boundary ∂Ω. The popularity of the Gaussian scale-space is due to its linearity and the fact that the multiscale function ψ(f )(x, y; t) can be generated from the isotropic heat diffusion, i.e.,

u(x, y; t) = (f * K √ 2t )(x, y) = ψ(x, y; t), t > 0.
PDE model [START_REF] Brockett | Evolution equations for continuous-scale morphology[END_REF] can also be solved using finite differences in an explicit scheme. Pixel i represents the location (x i , y i ). Let h l denote the grid size in the l direction: working on square grid, we assume m = 2 with h l = 1 for both horizontal and vertical directionsand τ denote the time step size (to guarantee stability, the step size must satisfy [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF]: τ = 1/( ∑ m l=1 2/h 2 l )), considering discrete times t k = kτ (with k positive integer). By u k i we denote approximation to u(x i , y i ; t k ). The simplest discretization can be written in a compact way as [START_REF] Weickert | Efficient image segmentation using partial differential equations and morphology[END_REF]:

u k+1 i = u k i + τ   m ∑ l=1 ∑ j∈N l (i) c k j + c k i 2h 2 l (u k j -u k i )   , ( 12
)
where N l (i) consists of the two neighbours of pixel i along the l discretized direction. The conduction coefficients c k i are considered here constant in time and space, but they will be time/space variant in next Section.

Counter-harmonic diffusion equation

In order to introduce the counter-harmonic counterpart of problem [START_REF] Brockett | Evolution equations for continuous-scale morphology[END_REF], let us start from the expected time discretization of the iterative solution: v(x, y; t + δt) = v(x, y; t) P +1 + δt∆v P +1 v(x, y; t) P + δt∆v P .

We note that v(x, y; t + δt) P = v(x, y; t) P + δt∆v P . Thus we can rewrite it as v(x, y; t + δt)v(x, y; t + δt) P = v(x, y; t) P +1 + δt∆v P +1 , son finally we have v(x, y; t + δt) P +1 -v(x, y; t) P +1 δt = ∆v P +1 .

Therefore, we obtain the following continuous model associated to the counter-harmonic scale-space:

{ ∂ t (v P +1 ) = ∆(v P +1 ), (x, y) ∈ Ω ⊂ R 2 , t ∈ R + v(x, y; 0) = f (x, y), f (x, y) > 0 ( 13
)
which can be seen as a simple nonlinear change of variable of the linear diffusion equation, i.e., u(x, y; t) → u(x, y; t) P +1 . However this form of the model does not give us an easy interpretation. We notice by the way that the solution of ( 13) is not at all equal to the solution of ( 11) with initial condition u(x, y; 0) = f (x, y) P +1 . Actually, we can reformulate ( 13) by two equivalent PDE problems.

Proposition 13

The solution v(x, y; t), with v > 0 for any (x, y) 13) is equivalent to either the solution of the Cauchy problem:

∈ Ω ⊂ R 2 , t ∈ R + , of problem (
{ ∂ t v = P ∥∇v∥ 2 v + ∆v v(x, y; 0) = f (x, y) > 0 ( 14
)
or to the Cauchy problem with change of variable v(x, y; t) → log e v(x,y;z) :      ∂ t w = (P + 1)∥∇w∥ 2 + ∆w w(x, y; 0) = log v(x, y; 0) = log f (x, y) v(x, y; t) = e w(x,y;t) (15)

Proof. Without a loss of generality, we consider the proof for the one dimensional case v(x, t), (x, t) ∈ R × R + , since the notation is simplified. Our starting point is the linear heat equation after the change of variable u = v P +1 , i.e.,

(v P +1 ) t = (v P +1 ) xx . ( 16
)
Let us compute the first derivative with respect to time and the second one with respect to space using the chain rule:

(v P +1 ) t = (P + 1)v P v t , (v P +1 ) xx = (P + 1)P v P -1 (v x ) 2 + (P + 1)v P v xx .
Substituting these formulae into equation ( 16), we obtain (P + 1)v P v t = (P + 1)P v P -1 (v x ) 2 + (P + 1)v P v xx , and finally, by canceling common factor (P + 1)v P , the evolutionary form which corresponds to ( 14) is gotten

v t = P (v x ) 2 v + v xx . ( 17
)
Now, we introduce another change of variable w(x, t) = log v(x, y) ⇔ v(x, y) = e w(x,y) .

Thus, applying again the chain rule, we get

v t = e w w t , v x = e w w x , v xx = e w (w x ) 2 + e w w xx .
The substitution into equation ( 17) gives e w w t = P (e w w x ) 2 e w + e w (w x ) 2 + e w w xx , then, by removing common factor e w and grouping factors in (w x ) 2 , the equivalent model ( 15) is obtained

w t = (P + 1)w x + w xx (18)
Obviously, the change of variable should be applied to the initial condition. Model (15) involves a clear understanding of the counter-harmonic diffusion which in fact is a nonlinear process sum of a linear diffusion and a (P + 1) weighted Hamilton-Jacobi process. That is, for large positive (resp. negative) P , the Hamilton-Jacobi term with positive sign (reps. negative sign) dominates over the diffusion term. Hamilton-Jacobi equations are of great importance in physics and play a central role in continuous morphology [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Arehart | Mathematical morphology: The Hamilton-Jacobi connection[END_REF][START_REF] Brockett | Evolution equations for continuous-scale morphology[END_REF][START_REF] Maragos | Differential morphology and image processing[END_REF]. More precisely, consider the Hamilton-Jacobi equation

∂u t (x, t) ± H (x, ∇u(x, t)) = 0 in R d × (0, +∞)
with the initial condition u(x, 0) = f (x) in R d . Such equations usually do not admit classic (i.e., everywhere differentiable) solutions but can be studied in the framework of the theory of viscosity solutions [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF]. It is well known [START_REF] Bardi | On Hopf's formulas for solutions of Hamilton-Jacobi equations[END_REF] that if the Hamiltonian H(x, p) = H(p) is convex then the solution of the Cauchy problem are given for + and -respectively by the so-called Lax-Oleinik formulae

u(x, t) = inf y∈R d [ f (y) + tH * ( x -y t )] and u(x, t) = sup y∈R d [ f (y) -tH * ( x -y t )]
where H * is the Legendre-Fenchel transform of function H. In the special case of (15), H(p) = (P + 1)|p| 2 such that H * (q) = |p| 2 4(P +1) and therefore, the Lax-Oleinik formulae of ∂ t w = (P + 1)∥∇w∥ 2 , with w(x, y; t) ∈ Ω × R + , w(x, y; 0) = log f (x, y), corresponds just to the following quadratic dilation (for P ≥ 0) and for -sign the quadratic erosion:

w(x, y; t) = sup (r,s)∈Ω ( log f (x -r, y -s) - (r 2 + s 2 ) 4(P + 1)t ) , for P ≥ 0, w(x, y; t) = inf (r,s)∈Ω ( log f (x -r, y -s) + (r 2 + s 2 ) 4(P + 1)t
) , for P < 0.

The time scaling t ′ = (P + 1)t is consistent with the asymptotic scale parameter in the convolution-based case of counter-harmonic Gaussian scale-space considered in Section 4. It is important to note that the hybrid model ( 15) is done on the logarithmic domain of the image.

Obviously P = 0 corresponds to linear diffusion. We note by the way that this case illustrates how the diffusion u t = u xx of the original image is equal to compute the logarithm of the image, then to apply the nonlinear diffusion process u t = (u x ) 2 + u xx and finally to compute the exponential of the result. Other relevant particular cases, such as P = -1 or P = -1/2 are related respectively to harmonic diffusion and the self-dual diffusion with respect to the inversion.

The use of the logarithmic transformation of the image is closely related to the use of a logarithmic brightness scale in order to achieve invariance to local multiplicative intensity transformations [START_REF] Lindeberg | Image matching using generalized scale-space interest points[END_REF] as well as exposure control mechanisms for visual receptive fields [START_REF] Lindeberg | A computational theory of visual receptive fields[END_REF]. Hence, this model of nonlinear diffusion is able to cope with some of the well known nonlinear phenomena in computer vision and human vision.

Relations with previous models of nonlinear diffusion

An hybrid evolutionary PDE such as ∂ t u = α∆u ± β∥∇u∥ 2 , with variable coefficients α and β which determines the part of linear versus morphological smoothing, was already conjectured in [START_REF] Brockett | Evolution equations for continuous-scale morphology[END_REF] and considered in [START_REF] Florack | Pseudo-Linear Scale-Space Theory[END_REF] by a different nonlinearization approach. It is related also to Kuwahara-Nagao PDE introduced in [START_REF] Van Den Boomgaard | Decomposition of the Kuwahara-Nagao operator in terms of a linear smoothing and a morphological sharpening[END_REF].

We have obtained the final PDE model ( 15) from the heat equation ( 11) by a composed change of variable of type u = v (P +1) followed by v = e w . In fact, in the area of nonlinear partial differential equations, this kind of trick is known as the Hopf-Cole transformation [START_REF] Cole | On a quasi-linear paribolic equation occurring in aerodynamics[END_REF][START_REF] Hopf | The partial differential equation u t + uu x = u xx[END_REF], i.e., u(x, t) = e αϕ(x,t) . More precisely, from heat equation u t = γu xx , the equation ϕ t = γα(ϕ x ) 2 + γϕ xx is obtained by the Hopf-Cole transformation. Then, by a second change of variable v = ϕ x , a form of the well-known Burgers' equation v t = 2γαvv x +γv xx is achieved. Thus for us it is remarkable that a continuous setting of the counter-harmonic mean is related to a milestone nonlinearisation model of physics.

Numerical solution and its decoupled approximation

Based on the motivation of this section, numerical solution of the original PDE problem ( 13) is related to the following discretization:

v k+1 i = (v k i ) P +1 + τ ( ∑ m l=1 ∑ j∈N l (i) c k j +c k i 2h 2 l ((v k j ) P +1 -(v k i ) P +1 ) ) (v k i ) P + τ ( ∑ m l=1 ∑ j∈N l (i) c k j +c k i 2h 2 l ((v k j ) P -(v k i ) P ) ) . ( 19
)
By this coupled evolution of both power terms discretized using basic finite differences, we avoid the need of a numerical scheme for the hybrid PDE [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], which requires a more sophisticated discretization. The corresponding scale-space associated to v(x, y; t) with initial condition the image f (x, y) is denoted η(f )(x, y; t; P ) As we have mentioned above, the result of coupled CHM diffusion η(f )(x, y; t; P ) are different of those obtained of the CHM Gaussian scale-space, which can be obtained by

η(f )(x, y; t; P ) = [u(x, y; t)] P +1 [u(x, y; t)] P , ( 20
)
where [u(x, y; t)] Q is the solution of standard diffusion equation ( 12) with the initial condition

u(x, y; 0) = f (x, y) Q .
From a theoretical viewpoint, we can see that at short time scales (t → 0) both solutions η(f )(x, y; t; P ) and η(f )(x, y; t; P ) should be similar. We note however that nonlinearization involves a normalization t → (P + 1)t. Nevertheless, at this point, we study empirically if the decoupled solution ( 20) is a good approximation of the one obtained by [START_REF] Dorst | Morphological Signal Processing and the Slope Transform[END_REF]. Let start our analysis with the CHM diffusion of a 1D signal given in Fig. 4. It concerns a very simple ramp signal (black signal in both figures) where the decoupled iterative scheme η(f )(x, y; t; P ) and the coupled one η(f )(x, y; t; P ) are compared, at same spatial scale t = 5, with respect to various values of P ≤ 0. Obviously the filtered signal at order P = 0 is the same in both cases as well as the limit case P = -∞, which corresponds to a flat erosion (yellow signal in both figures). We observe that for small |P | (e.g., P = -1 or -3), their behaviour is similar; however, for high values of |P |, we notice that the decoupled scheme produces smoother results. Hence, the main qualitative property of the coupled scheme is to produce a nonlinear diffusion with a notably sharp effect on the regularization of the structures. In Fig. 5 a comparison of CHM isotropic diffusion using also both the decoupled iterative scheme vs. the coupled iterative scheme is given, at spatial scale t = 5 and order P = 10. We note from this experiment of pseudo-dilation that the same effect of sharpness on η(f )(x, y; t; P ) is observed. We can conclude that the iteration n times of a unitary (i.e., σ = √ 2) generalised morphological Gaussian convolution introduces a nonlinearity degree superior to the corresponding generalised morphological Gaussian convolution at scale σ = √ 2n.

Counter-Harmonic Locally-Adaptive Diffusion

Ideas introduced in previous section are extended now to two well-known cases of locally adaptive image diffusion.

Edge-preserving diffusion

The big disadvantage of isotropic diffusion is the fact that linear smoothers do not only smooth the noise but also blur and shift structure edges, which are important image features in some applications, i.e., segmentation and object recognition. The pioneering idea introduced by Perona and Malik [START_REF] Perona | Scale-Space and Edge Detection Using Anisotropic Diffusion[END_REF] to reduce these problems is to locally adapt the diffusivity, i.e., the value of the conductivity c, to the gradient of the image at each iteration. More precisely, this nonlinear diffusion involves replacing the diffusion equation ( 11) by the following model:

∂ t u = div ( g ( ∥∇u∥ 2 ) ∇u ) (21)
u(x, y; 0) = f (x, y)

∂ n u| ∂Ω = 0
In this model the diffusivity has to be such that g ( ∥∇u∥ 2 ) → 0 when ∥∇u∥ 2 → +∞ and g ( ∥∇u∥ 2 ) → 1 when ∥∇u∥ 2 → 0. One of the diffusivities Perona and Malik proposed is the function

g λ (s) = 1 1 + s/λ 2 , λ > 0, such that g ′ λ (s) = - 1 λ 2 (1 + s/λ 2 ) 2
where λ is a threshold parameter that separates forward and backward diffusion. In effect the diffusivity entails a forward diffusion at location where ∥∇u∥ < λ and backward diffusion where ∥∇u∥ > λ. This model accomplishes the aim of blurring small fluctuations (noise) while preserving edges (by preventing excessive diffusion). Backward diffusion is an ill-posed process. In fact, to see the complexity of the Perona and Malik, the right-side of model ( 21) can be rewritten as

∇ • ( g ( ∥∇u∥ 2 ) ∇u ) = g ( ∥∇u∥ 2 ) ∆u + 2g ′ ( ∥∇u∥ 2 ) D 2 u∇u • ∇u, ( 22
)
where

D 2 u is the Hessian of u, i.e., 2g ′ ( ∥∇u∥ 2 ) D 2 u∇u • ∇u = 2g ′ ( ∥∇u∥ 2 ) ( (∂ x u) 2 ∂ 2 x u + (∂ y u) 2 ∂ 2 y u + 2(∂ x u)(∂ y u)∂ x ∂ y u
) .

To avoid underlying numerical and theoretical drawbacks of this model, it was proposed in [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], a new version of Perona and Malik theory, based on replacing diffusivity g ( ∥∇u∥ 2 ) by a regularized version g

( ∥∇u σ ∥ 2 ) with ∇u σ = ∇ (K σ * u)
where K σ is a Gaussian kernel. This latter model is just used in our framework and the numerical solution can be obtained using the explicit scheme [START_REF] Bullen | Handbook of Means and Their Inequalities[END_REF], where the local conductivity is here the approximation to the regularized term g ( ∥∇u σ ∥ 2 ) , i.e.,

c k i = 1 1 + ∥∇uσ(x i ,y i ;t k )∥ 2 λ 2 , ( 23
)
with the gradient computed by central differences. A critical point here is the choice of parameter λ, which is an ill-posed problem. Sophisticated approaches in the state-of-the-art are based on noise estimation. We have determined empirically that the value

λ = 1 m ( m ∑ l=1 Mean(∇ l u(x, y; t)) ) , ( 24
)
where m is the number of discretized directions, produces stable results. -dilation/erosion by counter-harmonic nonlinear diffusion using coupled scheme, (e)/(f) pseudo-dilation/erosion by counter-harmonic nonlinear diffusion using decoupled scheme.

Counter-harmonic Perona and Malik diffusion model

Similarly to the isotropic case, we can now formulate the partial differential equation underlying Perona and Malik process associated to the numerical scheme [START_REF] Dorst | Morphological Signal Processing and the Slope Transform[END_REF]. The starting point will be now the model

{ ∂ t (v P +1 ) = ∇ • ( g ( ∥∇v P +1 ∥ 2 ) ∇v P +1 ) , (x, y) ∈ Ω ⊂ R 2 , t ∈ R + v(x, y; 0) = f (x, y), f (x, y) > 0 (25)
The solution v(x, y; t) of ( 25) can be obtained from the following problem, let us called generalised morphological Perona and Malik diffusion model:

             ∂ t w = (P + 1) [ g(ξ)∥∇w∥ 2 + 2g ′ (ξ)e 2(P +1)w ∥∇w∥ 4 ] + g(ξ)∆w + 2g ′ (ξ)e 2(P +1)w D 2 w∇w • ∇w with ξ = ∥e (P +1)w ∇w∥ 2 , w(x, y; 0) = log v(x, y; 0) = log f (x, y) v(x, y; t) = e w(x,y;t) (26)
The proof of the relationship between ( 25) and ( 26) is based first on developing the derivatives, using chain rule, to obtain

∂ t v = g ( ∥v P ∇v∥ 2 ) P ∥∇v∥ 2 v + g ( ∥v P ∇v∥ 2 ) ∆v + [ ∇ • g ( ∥v P ∇v∥ 2 )] • ∇v,
Note that we have also introduced a normalization by dividing v P +1 by (P + 1). This is useful to simplify the model. Then, using the change of variable v = e w , it is rewritten as the PDE:

∂ t w = g ( ∥e (P +1)w ∇w∥ 2 ) (P + 1)∥∇w∥ 2 + g ( ∥e (P +1)w ∇w∥ 2 ) ∆w+ [ ∇ • g ( ∥e (P +1)w ∇w∥ 2 )] • ∇w.
The last term of this partial differential equation can be developed as follows:

[ ∇ • g ( ∥e (P +1)w ∇w∥ 2 )] • ∇w = [ g ′ ( ∥e (P +1)w ∇w∥ 2 ) 2∥e (P +1)w ∇w∥∇∥e (P +1)w ∇w∥ 2 ] • ∇w = 2g ′ ( ∥e (P +1)w ∇w∥ 2
) [ (P + 1)e 2(P +1)w ∥∇w∥ 4 + e 2(P +1)w D w ∇w • ∇w ] .

We give now a preliminary qualitative analysis of the model [START_REF] Heijmans | Morphological scale-spaces, scale invariance, and Lie groups[END_REF]. First, if we consider a positive diffusivity function g λ (ξ) =

( 1 + ξ/λ 2 ) -1 , with ξ = ∥e (P +1)w ∇w∥ 2 , we observe that the selection of the forward ξ < λ or the backward ξ > λ modes for the morphological term (P + 1)∥∇w∥ 2 or the diffusion term ∆u, will depend on the product e (P +1)w ∇w. That means that ∥e (P +1)w ∇w∥ → 0 when either e (P +1)w → 0 or ∥∇w∥ → 0. In other words, for positive (resp. negative) P the pixels with low (resp. high) intensities will be dilated (resp. eroded) or diffused independently of the gradient term. As a consequence, the edges of bright (resp. dark) well contrasted areas are preserved, the dark (resp. bright) zones are regularized. In addition, we have two additional terms on the negative function g ′ λ (ξ). The diffusion related term D 2 w∇w • ∇w is well known in nonlinear PDE and it can seen as a second order approximation of the morphological parabolic laplacian [START_REF] Guichard | A Note on Two Classical Enhancement Filters and Their Associated PDE's[END_REF], which weighted by g ′ λ (ξ)e 2(P +1)w , being of negative sign, implies an enhancement of edges. The term (P +1)∥∇w∥ 2 , with the negative sign of g ′ λ (ξ)e 2(P +1)w , involves an erosion (resp. a dilation) for positive (resp. negative) P , where the structuring function is of shape -∥x 2 + y 2 ∥ 4/3 /(Ct 1/3 ). In fact, a PDE of type ∂ t u = ±∥∇u∥ q , q ≥ 1 involves multi-scale dilations (+) and erosions (-) with structuring functions of shape depending on q, see in [START_REF] Diop | Multiscale Image Analysis Based on Robust and Adaptive Morphological Scale-Spaces[END_REF]. Therefore, both terms in (P + 1) (which will be dominant over diffusion for large |P |), i.e., g(ξ)∥∇w∥ 2 + 2g ′ (ξ)e 2(P +1)w ∥∇w∥ 4 , involve a combination of dilation/erosion effect according to the nature of the pixel intensity and the edgeness. In addition, inhomogeneous Hamilton-Jacobi model of type ∂ t u = g(ξ)∥∇w∥ 2 leads to viscosity solutions with sup-convolution where the structuring function depends on a geodesic distance. This geodesic distance is driven by g(ξ), which can be seen as related to the local metric of a Riemannian manifold. This later point is related to the recently framework of Riemannian mathematical morphology [START_REF] Angulo | Riemannian Mathematical Morphology[END_REF]. Nevertheless, the complexity of this model requires a deeper analysis out of the scope of this paper.

The generalised morphological Perona and Malik model ( 25) can be solved using the coupled scheme [START_REF] Dorst | Morphological Signal Processing and the Slope Transform[END_REF], with diffusivities g λ P +1

( ∥∇v P +1 ∥ 2 ) and g λ P ( ∥∇v P ∥ 2 ) for the weights c k i , as in [START_REF] Gonzalez | Digital Image Processing[END_REF]. Note that parameter λ should includes a normalization respectively by the power as follows

λ Q = Q -1 1 m ( m ∑ l=1 Mean(∇ l u(x, y; t)) ) ,
and which of course is different for P and P + 1.

As for the case of the counter-harmonic isotropic diffusion, this model at order P can be also approximated by a decoupled numerical estimation as

ξ(f )(x, y; t; P ) = [u(x, y; t)] P +1 nl [u(x, y; t)] P nl , ( 27
)
where [u(x, y; t)] Q nl is the solution of regularized Perona and Malik diffusion equation ( 21) with the initial condition u(x, y; 0) = f (x, y) Q .

Fig. 6 provides a comparison of decoupled and coupled counter-harmonic Perona and Malik diffusion ξ(x, y; t; P ) at scale t = 5 (regularization parameter σ = 0.5). The original image is a noisy binary image which in all the cases is restored without blurring the edges of the objects. This example is interesting to illustrate how the generalised Perona and Malik diffusion can be used to emphasize the bright or the dark image structures during the diffusion process and simultaneously preserving the edges. The corresponding pseudo-morphological operators for the coupled approach produce sharper results than the decoupled ones. But the decoupled approach can be considered as good approximation. The results in both cases are consistent with the qualitative analysis of the model done just above.

Fig. 7 provides the results of pseudo-dilation and pseudo-erosion (P = 10 and P = -10) using the counter-harmonic Perona and Malik model on a noisy image. Comparing with respect to the standard filtering (P = 0), the good properties of denoising without blurring are conserved but in addition, an effect of dilation/erosion is obtained. This kind of pseudomorphology is useful for instance to compute morphological gradient, i.e., ξ(f )(x, y; t; P ) -ξ(f )(x, y; t; -P ), of noisy images or to construct filters as morphological openings, i.e., )(x, y; t; -P ))(x, y; t; P ). See in Fig. 9(c) and (g) the comparison of the standard Perona and Malik diffusion and the corresponding pseudo-opening morphological diffusion, which removes some small bright structures.

Coherence-enhanced diffusion

Many other models of anisotropic diffusion have been explored in the state-of-the-art, which deal with an edge-preserving diffusion effect similar to Perona and Malik model, but improving the combination of edge preserving and enhancement, see for instance [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. Instead of exploring these approaches, we prefer to consider how the counter-harmonic paradigm is extended to the case of coherence-enhancing model.

Coherence-enhanced diffusion, or tensor-driven diffusion, was introduced by Weickert [54] in order to achieve an anisotropic diffusion filtering based on directionality information. The idea is to adapt the diffusion process to local image structure using the following nonlinear diffusion equation:

∂ t u = div (D∇u) . ( 28
)
where the conductivity function becomes a symmetric positive definite diffusion tensor D, which is a function adapted to the structure tensor:

J ρ (∇u σ ) = K ρ * (∇u σ ⊗ ∇u σ ).
The eigenvectors of J ρ are orthonormal and the eigenvalues are no negative. The corresponding eigenvalues (let us call them µ 1 ≥ µ 2 ) describe the local structure. Flat areas give µ 1 ≈ µ 2 , straight edges give µ 1 ≫ µ 2 = 0 and corners give µ 1 ≥ µ 2 ≫ 0. In order to control the diffusion, J ρ is not used directly, but tensor D has the same eigenvectors as J ρ , but different eigenvalues, thus controlling the diffusion in both directions. The eigenvalues are

{ λ 1 = α λ 2 = α + (1 -α) exp(-C/κ)
where κ is the orientation coherence and C > 0 serves as a threshold parameter. Parameter α > 0 is quite important and serves as a regularization parameter which keeps D uniformly positive definite. For this diffusion, we have used in our tests a numerical implementation using the additive operator splitting (AOS) scheme [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF], which is particularly efficient and has the advantage of being rotationally invariant compared to their multiplicative counterparts.

Counter-harmonic Weickert diffusion model

In this case, the starting point for the counter-harmonic model will be the same as in the case of Perona and Malik but, instead of dealing with a scalar function g ( ∥∇v P +1 ∥ 2 ) , we have modified diffusion tensor D, which is related to the the structure tensor from v P +1 . Without doing a mathematical analysis of the corresponding model, which requires very technical issues out of the scope of the paper, let us at least write a first form of the counter-harmonic Weickert diffusion model as follows:

         ∂ t w = (P + 1) D∥∇w∥ 2 + D∆w + ∇ D • ∇w w(x, y; 0) = log v(x, y; 0) = log f (x, y) v(x, y; t) = e w(x,y;t) (29)
This generalised morphological Weickert diffusion can be solved using the coupled scheme (19) or the AOS scheme in its counter-harmonic formulation. As previously, its solution can also be approached by the decoupled schema

χ(f )(x, y; t; P ) = [u(x, y; t)] P +1 ce [u(x, y; t)] P ce , (30)
where [u(x, y; t)] Q ce is the solution of Weickert diffusion equation ( 28) with the initial condition u(x, y; 0) = f (x, y) Q . As in the previous case, we have to adapt the regularization parameter α to the dynamics of the power images f (x, y) Q in order to have numerical stable results. Empirically, we have observed that

α = { 0.01|Q| if Q ̸ = 0 0.005 if Q = 0 (31)
leads to satisfactory results. Moreover, we have also observed in experimental tests that numerical solution of the CHM coherence diffusion given in expression [START_REF] Guichard | A Note on Two Classical Enhancement Filters and Their Associated PDE's[END_REF] can introduce some instabilities and errors for P > 0, whereas for P < 0 the results are numerically stable.

To avoid this drawback, we propose to compute, by duality, the coherence-enhanced pseudodilation χ(f )(x, y; t; P ), with P > 0, by means of the operation χ(f -1 )(x, y; t; -P -1) -1 , based on property [START_REF] Bardi | On Hopf's formulas for solutions of Hamilton-Jacobi equations[END_REF]. In any case, a deeper study on the numerical implementation of counter-harmonic Weickert diffusion, in particular in the case of coupled diffusion, is required.

Examples of pseudo-morphological anisotropic diffusion are given in Fig. 8 and Fig. 9(d)-(h). Dark regions are anisotropically pronounced in pseudo-erosion schemes (P < 0) whereas bright regions are anisotropically emphasized in pseudo-dilation as well as in their products, respectively the pseudo-openings and pseudo-closings. Finally, let us consider the case study summarised in Fig. 13. The vessels are dark structures which can be extracted using the residue of a morphological closing (the dual top-hat transform [START_REF] Soille | Morphological Image Analysis[END_REF]). From this image, the edges can be obtained using for instance the Canny edge detector [START_REF] Canny | A computational approach to edge detection[END_REF] (we have also used the implementation available in [START_REF] Kovesi | MATLAB and Octave Functions for Computer Vision and Image Processing. Centre for Exploration Targeting[END_REF]). We note the vessels are anisotropic structures which can be enhanced using the Weickert diffusion model. In Fig. 13 are compared the edges obtained using the same detection parameters from the original image f , the image χ(x, y; t = 5; P = 0)(f ) (standard Weickert diffusion) and the image χ(x, y; t = 5; P = 5)(χ(x, y; t = 5; P = -5)(f ′′ )). We consider that this example shows well that a prior knowledge of the nature of the structures to be enhanced (e.g., anisotropic and dark ones) can be used to tailor the most suitable generalised morphological diffusion. Figure 11: Denoising assessment in terms of SNR with respect to order of the generalised morphological nonlinear diffusion ξ(x, y; t; P )(f ′ ) (fixed t = 5): (a) pseudo-erosion ξ(x, y; t; P < 0)(f ′ ) and pseudo-dilation ξ(x, y; t; P > 0)(f ′ ) for speckle noise at four values of σ; (b) pseudo-opening ξ(x, y; t; P > 0)(ξ(x, y; t; P < 0)(f ′ )) and pseudo-closing ξ(x, y; t; P < 0)(ξ(x, y; t; P > 0)(f ′ )) for speckle noise at four values of σ; (a) pseudo-erosion ξ(x, y; t; P < 0)(f ′′ ) and pseudo-dilation ξ(x, y; t; P > 0)(f ′′ ) for asymmetric impulse noise at four values of σ; (b) pseudo-opening ξ(x, y; t; P > 0)(ξ(x, y; t; P < 0)(f ′′ )) and pseudo-closing ξ(x, y; t; P < 0)(ξ(x, y; t; P > 0)(f ′′ )) for asymmetric impulse noise at four values of σ.

Conclusions and perspectives

We have introduced a mathematical framework of generalised morphological image diffusion using the notion of counter-harmonic mean. In particular, we have proposed the counterharmonic Gaussian scale-space, and we have studied in detail its properties. In the case of 2D images, the nonlinear smoothed signals η(f )(x, y; t; P ) taken as a whole can be considered as a continuous function on R 2 × R + × R ∪ {+∞, -∞}, that is, a double pseudo-scale-space with respect to t and with respect to P . This family of filters expands and includes the Gaussian scale-space and the morphological erosion/dilation scale-spaces.

We have introduced the nonlinear diffusion PDE which provides an hybrid diffusion/Hamilton-Jacobi model for nonlinear continous image filtering. We have also proposed the counterharmonic inspired PDE from other models of locally adaptive diffusion: on the one hand, Perona and Malik model of edge-preserving diffusion, on the other hand, Weickert model of coherence-enhanced diffusion.

As the dilation/erosion operators depend on the extreme values of the additively weighted signal in the neighbourhood of a point, impulse noise and other kind of high-frequency noise cause spurious results of the standard morphological processing of noisy images. Linear diffusion filtering is able to deal with the noise, however the averaging process blurs and limits also the effect of the filtering in the case of asymmetric perturbation. As we have shown in the applications, the fundamental property of the present generalised image diffusion is its ability to robustly filter out the image structures in a nonlinear way, which can selectively bias the result toward high (low) values such as the dilation (erosion) filter according to the specific nature of noisy or target structures.

A deeper analytical study on the numerical implementation of counter-harmonic Perona-Malik and Weickert diffusion models, and in particular in the case of coupled diffusion schemas, should be achieved. We have pointed out that, for data defined in the interval (0, 1], the behaviour of filters with P < 0 is generally more stable than the case P > 0. Consequently, the use of the duality by inversion using Proposition 2 seems an appropriate alternative to solve the cases P > 0 by means of the dual -P -1. However a better understating of these numerical effects is still to be investigated.

Concerning the discrete framework, we think that the genuine discrete theory of Gaussian scale-space by Lindeberg [START_REF] Lindeberg | Scale-Space for Discrete Signals[END_REF], can be used for defining a discrete analogue of the counterharmonic scale-space. This issue could be considered in the future since it would lead to better properties and more efficient discrete implementations.

Additionally, it should be considered in ongoing work how the present nonlinearization of image diffusion can be turned to the Laplace-Beltrami framework to formulate diffusion in Riemannian manifolds [START_REF] Sochen | A general framework for low level vision[END_REF][START_REF] Kimmel | Image Processing via the Beltrami Operator[END_REF]. In particular, we would like to explore connections with the paradigm introduced by Florack in [START_REF] Florack | Scale Space Representations Locally Adapted to the Geometry of Base and Target Manifold[END_REF]: i.e.,

∂ t U = α∥∇U ∥ 2 g + ∆ g U,
where U (x; t) is a function defined in M × R + , such that M is a Riemannian manifold of metric tensor g µν and where ∆ g U = ( √ g) -1 ∂µ(g µν √ g∂µU ) is the Laplace-Beltrami operator on M and ∥∇U ∥ 2 g = g µν ∂ µ U ∂ ν U . From a more practical viewpoint, the theory presented in this paper has been used in [START_REF] Masci | Learning Framework for Morphological Operators using Counter-Harmonic Mean[END_REF] for machine learning morphological operators. That means to automatically select the structuring function, or convolution kernel, and the parameter of nonlinearity P from a given set of training images. The framework is founded on the fact that the counter-harmonic scale-spaces are differentiable with respect to the kernel and to P . Therefore, optimization techniques, such as stochastic gradient descent algorithms, underlying convolutional neural networks training, can be applied to our operators. We believe that this paradigm deserves additional developments since such approach opens the door to more generic nonlinear convolutional neural networks, which are at the heart of the deep learning paradigm in computer vision.

Figure 1 :

 1 Figure 1: Two image examples presenting asymmetric bright/dark image structures. These images will be studied in Section 7.

Figure 2 :

 2 Figure 2: Convergence with respect to P > 0 of nonlinear power means-based operators to morphological operators: (a) pseudo-dilation M P B and pseudo-erosion M -P B using Minkowski mean; (b) pseudo-dilation κ P B and pseudo-erosion κ -P B

Figure 3 :

 3 Figure 3: Generalised morphological Gaussian convolution (isotropic counter-harmonic diffusion) η(f )(x, y; t; P ) at scale t = 5: First row, original image, standard Gaussian filtered image (P = 0) and cases P = 1 and P = -1; middle row, counter-harmonic Gaussian pseudo-dilations (P > 0); bottom row, counter-harmonic Gaussian pseudo-erosions (P < 0).

Figure 4 :

 4 Figure 4: Comparison of counter-harmonic diffusion of a 1D signal using in (a) the decoupled iterative scheme η(f )(x, y; t; P ) (a) and in (b) the coupled iterative scheme η(f )(x, y; t; P ), at spatial scale t = 5 and with respect to various values of P ≤ 0.

Figure 5 :

 5 Figure 5: Comparison of counter-harmonic isotropic diffusion using the decoupled iterative scheme (equivalent to the generalised Gaussian convolution) vs. the coupled iterative scheme, at spatial scale t = 5 and order P = 10: (a) original image f (x, y), (b) standard dilation, (c) decoupled η(f )(x, y; t; P ) and (d) coupled η(f )(x, y; t; P ). At the bottom are zoom-in frames of a square section cropped from the corresponding images.

Figure 6 :

 6 Figure 6: Comparison of decoupled and coupled generalised morphological Perona and Malik diffusion ξ(x, y; t; P ) at scale t = 5 (regularization parameter σ = 0.5 ): (a) original image, (b) standard nonlinear diffusion, (c)/(d) pseudo-dilation/erosion by counter-harmonic nonlinear diffusion using coupled scheme, (e)/(f) pseudo-dilation/erosion by counter-harmonic nonlinear diffusion using decoupled scheme.

  (a) Original (b) P = 0 (c) P = +∞ (d) P = -∞ (e) Difference (c') P = 10 (d') P = -10 (e') Difference

Figure 7 :

 7 Figure 7: Generalised morphological Perona and Malik diffusion ξ(x, y; t; P ) at scale t = 5 (regularization parameter σ = 0.5 ): (a) original image, (b) standard nonlinear diffusion, (c)/(d) standard dilation/(d) standard erosion, (e) image difference between dil./ero. (gradient), (c')/(d') pseudo-dilation/erosion by counter-harmonic nonlinear diffusion, (e') image difference between pseudo-dil./ero.

5 Figure 8 :Figure 9 :

 589 Figure 8: Generalised morphological Weickert diffusion χ(x, y; t; P ) at scale t = 10 (with regularization parameters: local scale σ = 1.5, integration scale ρ = 6) : (a) original image, (b) standard anisotropic diffusion, (c)-(d) opening and closing of size equivalent to t = 5, (e)/(f) pseudo-dilation/erosion by counter-harmonic anisotropic diffusion, (g)/(h) pseudoopening/closing by counter-harmonic anisotropic diffusion.

Figure 10 :

 10 Figure 10: Image denoising using generalised morphological nonlinear diffusion ξ(x, y; t; P ): (a) original unnoisy image, f ; (b) image corrupted with speckle noise σ = 0.01, f ′ ; (c) image corrupted with asymmetric impulse noise σ = 0.1, f ′′ ; (d) restoration of image (b) using ξ(x, y; t = 5; P = 0)(f ′ ); (e) restoration of image (b) using ξ(x, y; t = 5; P = -3)(f ′ ); (f) restoration of image (b) using ξ(x, y; t = 5; P = 3)(ξ(x, y; t = 5; P = -3)(f ′ )); (g) restoration of image (c) using ξ(x, y; t = 5; P = 0)(f ′′ ); (h) restoration of image (c) using ξ(x, y; t = 5; P = -5)(f ′′ ); (i) restoration of image (c) using ξ(x, y; t = 5; P = 5)(ξ(x, y; t = 5; P = -5)(f ′′ )).

  Pair of Counterharmonic Perona-Malik (P,-P) σ = 0.002 (-P,P) σ = 0.005 (-P,P) σ = 0.01 (-P,P) σ = 0.02 (-P,P) σ = 0.002 (P,-P) σ = 0.005 (P,-P) σ = 0.01 (P,-P) σ = 0.02 (P,-P) Pair of Counterharmonic Perona-Malik (P,-P) σ = 0.02 (-P,P) σ = 0.05 (-P,P) σ = 0.1 (-P,P) σ = 0.15 (-P,P) σ = 0.02 (P,-P) σ = 0.05 (P,-P) σ = 0.1 (P,-P) σ = 0.15 (P,-P)

Figure 12 :

 12 Figure 12: Example of Harris corner detector: (a) original asymmetric impulse noisy image f ′′ ; (b) denoised image ξ(x, y; t = 5; P = 0)(f ′′ ) (standard Perona and Malik diffusion); (c) denoised image ξ(x, y; t = 5; P = 5)(ξ(x, y; t = 5; P = -5)(f ′′ )) (pseudo-opening); (d) keypoints of image (a); (e) keypoints of image (b); (f) keypoints of image (c).

Figure 13 :

 13 Figure 13: Edge vessel detection on retinal image using Canny detector: (a) original image f ; (b) enhanced image χ(x, y; t = 5; P = 0)(f ) (standard Weickert diffusion); (c) enhanced image χ(x, y; t = 5; P = 5)(χ(x, y; t = 5; P = -5)(f ′′ )); (d) dual top-hat of image (a); (e) dual top-hat of image (b); (f) dual top-hat of image (c); (g) detected edges from image (d); (h) detected edges from image (e); (i) detected edges from image (f).

  λ 2 y; λ 3 t; λ 4 P ), for any λ 1 , λ 2 , λ 3 , λ 4 > 0. = λ 3 λ 4 and if |P | is large enough or if λ 2 2 = λ 3 and P = 0. In such a case, there exist constants k 1 and k 2 and we have

	Proposition 11 A functional W of the nonlinear scale-space η(f )(x, y; t; P ) is dimensional
	if and only if λ 2 2

Counter-Harmonic Mean (CHM)Let us start by presenting the basic notion of this paper.

Acknowledgements. The author would like to thank the anonymous Reviewer who pointed out the duality of the CHM by the reciprocal as well as the case of self-duality. Many thanks also to the Reviewer who stimulate us in the formulation of the counter-harmonic diffusion PDE.

Applications to asymmetric image processing

In the previous sections, we have illustrated with various examples the comparative effects of generalised morphological diffusion. Our aim now is to demonstrate the significance of the approach for practical image analysis. Let us start with the problem of images corrupted with non-Gaussian noise. Fig. 10(a) shows the original unnoisy f image which has been corrupted using two kinds of noise: i) standard speckle noise f ′ , i.e., multiplicative noise f ′ = f + nf , where noise n is uniformly distributed random noise with mean 0 and variance σ; ii) asymmetric (negative) impulse noise f ′′ , generated as f ′′ = f -n, where n depends on the absolute value of the logarithm of an uniform noise with variance σ. Figs. 10(b) and (c) show just an example of each noisy image. We have simulated for both types of noise 10 images at four different values of σ.

We consider image denoising using the generalised morphological nonlinear (Perona and Malik) diffusion ξ(x, y; t; P ), where the nonlinear order P has been discretised in the interval -20 ≤ P ≤ 0 (pseudo-erosions) and in the interval 0 ≤ P ≤ 20 (pseudo-dilations) as well as their products (pseudo-openings and pseudo-closings). Finally, we have computed the corresponding SNR. The averaged values of SNR with respect to P are given in the curves of Fig. 11. The interest of these experiments is just to study the behaviour of ξ(x, y; t; P ) with respect to P , or more precisely to find which P produces the highest SNR. Figs. 11(d As we can observe from the curves of speckle noise, for any σ, the maximum SNR corresponds to a value -3 ≤ P ≤ 0, that is a pseudo-erosion. This result is confirmed by the fact that pseudo-openings at these values perform better than the pseudo-closing. However, the results are not very impressive for such a noise. The case of asymmetric negative impulse noise is more conclusive, since for significant noise levels, it is observed that the maximum SN R is much better for a high positive value the of order: 5 ≤ P ≤ 10. As expected, the results of pseudo-closings are better than the pseudo-openings.

The interest of a pre-processing step such as the generalised morphological diffusion ξ(x, y; t; P ) is not only to improve the SNR, but also to produce a suitable image regularisation which improves other subsequent image analysis steps. Hence, we have performed a Harris corner detector [START_REF] Harris | A combined corner and edge detector[END_REF], computed using the implementation available in [START_REF] Kovesi | MATLAB and Octave Functions for Computer Vision and Image Processing. Centre for Exploration Targeting[END_REF]. Fig. 12 provides three examples of the obtained key-points; in particular, it is compared, using the same detection parameters, the original asymmetric impulse noisy image f ′′ , the image ξ(x, y; t = 5; P = 0)(f ′′ ) (standard Perona and Malik diffusion) and the image ξ(x, y; t = 5; P = 5)(ξ(x, y; t = 5; P = -5)(f ′′ )) (pseudo-opening). As we can observe the keypoints obtained from the generalised morphological diffusion are more significant than the ones obtained from the standard Perona and Malik diffusion.