# Extension of p-Laplace Operator for Image Denoising

George Baravdish, Yuanji Cheng, Olof Svensson, Freddie Åström

► To cite this version:

George Baravdish, Yuanji Cheng, Olof Svensson, Freddie Åström. Extension of p-Laplace Operator for Image Denoising. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.107-116, 10.1007/978-3-319-55795-3_9. hal-01626927

HAL Id: hal-01626927

https://inria.hal.science/hal-01626927

Submitted on 31 Oct 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
# Extension of p-Laplace Operator for Image Denoising

George Baravdish¹, Yuanji Cheng², Olof Svensson¹, and Freddie Åström³*

¹ Linköping University, Sweden
{george.baravdish, olof.svensson}@liu.se

² Malmö University, Sweden
yuanji.cheng@mah.se

³ Heidelberg Collaboratory for Image Processing, Heidelberg University, Germany
freddie.astroem@iwr.uni-heidelberg.de

build

**Abstract.** In this work we introduce a novel operator $\Delta_{(p,q)}$ as an extended family of operators that generalize the p-Laplace operator. The operator is derived with an emphasis on image processing applications, and particularly, with a focus on image denoising applications. We propose a non-linear transition function, coupling p and q, which yields a non-linear filtering scheme analogous to adaptive spatially dependent total variation and linear filtering. Well-posedness of the final parabolic PDE is established via perturbation theory and connection to classical results in functional analysis. Numerical results demonstrates the applicability of the novel operator $\Delta_{(p,q)}$.

**Keywords:** p-Laplace operator, parabolic equations, image denoising, anisotropic diffusion, inverse problems.

## 1 Introduction

A well known inverse problem in image processing is image denoising [?]. In the last decades the energy functional approach together with its corresponding Euler-Lagrange (E-L) equation has attracted great attention in solving inverse problems applied to image reconstruction. One important case of E-L equations is the one which involves the p-Laplace operator

$$ \Delta_p u = \operatorname{div}(|\nabla u|^{p-2}\nabla u), \quad p \ge 1, \qquad (1) $$

associated with the evolution equation of p-Laplacian

$$ \left\{ \begin{array}{ll} \partial_t u - \Delta_p u = 0, & \text{in } \Omega \times (0, T) \\ u(0) = u_0, & \text{in } \Omega \\ \partial_n u = 0, & \text{on } \partial\Omega \times (0, T) \end{array} \right. \qquad (2) $$

<footnote>* Acknowledgment Support by the German Science Foundation and the Research Training Group (GRK 1653) is gratefully acknowledged.</footnote>
where $\Omega$ is a bounded domain in $\mathbb{R}^2$ and $u_0 : \Omega \to \mathbb{R}$ is a given degraded image [?], [?], [?] and $\nabla u$ is the gradient. The degenerate parabolic equation in (??) has been studied by many authors and we limit ourselves here to refer the reader to [?], It is well known that the case $p = 2$ gives the linear Gaussian filter, which however, impose strong spatial regularity and therefore image details such as lines and edges are oversmoothed. $p = 1$ is often refereed to as the method of total variation [?] and $p = 0$ is an instance of the so called balanced forward backward evolution [?].

In this work we study a decoupled form of the p-Laplace operator, expressed as a non-linear combination of the $\Delta_1$ and $\Delta_\infty$ operators, introduced below. We call our new operator $\Delta_{(p,q)}$. Via established existence theory we show that the corresponding perturbed parabolic equation is well-posed and close to the original operator.

In Section 2 we review some of the properties of the p-Laplace operator applied in image denoising and compare it with Perona-Malik models. Our main contribution is in section 3 where we extend the p-Laplace operator to a new operator $\Delta_{(p,q)}$ with focus on image denoising. We consider the corresponding variable version $\Delta_{(p(x),q(x))}$ in section 4. Finally, in section 5 we demonstrates the applicability of $\Delta_{(p,q)}$ by numerical results.

## 2 p-Laplacian for Image Denosing

An important feature in any evolution process for image denoising is preservation of certain geometrical features of the underlying image. In the case of image restoration these features include edges and corners. It is straight-forward to express the p-Laplace operator (??) as

$$ \Delta_p u = |\nabla u|^{p-1} \Delta_1 u + (p-1) |\nabla u|^{p-2} \Delta_\infty u, \quad (3) $$

where $\Delta_1 u = \text{div}(\frac{\nabla u}{|\nabla u|})$, $\Delta_\infty u = \frac{\nabla u}{|\nabla u|} \cdot (D^2 u) \frac{\nabla u}{|\nabla u|}$ and $D^2 u$ is the Hessian of $u$. However, an intuitive way to represent $\Delta_p$, giving direct interpretation of the diffusivity directions is to express $\Delta_p$ by using Gauge coordinates $(x, y) \to (\eta, \xi)$:

$$ \Delta_p u = |\nabla u|^{p-2} (u_{\xi\xi} + (p-1)u_{\eta\eta}) \quad (4) $$

where

$$ \eta = \frac{\nabla u}{|\nabla u|}, \quad \xi = \frac{\nabla^\perp u}{|\nabla u|}, \quad (5) $$

and

$$ u_{\xi\xi} = |\nabla u| \Delta_1 u, \quad u_{\eta\eta} = \Delta_\infty u. \quad (6) $$

From (??) it is now clear that $\Delta_p$ imposes the same diffusivity strength in both directions $\xi$ and $\eta$ independent of the magnitude of the gradient. In an attempt to resolve this drawback, Perona and Malik (P-M) [?] proposed to replace $|\nabla u|^{p-2}$ with $g(|\nabla u|^2)$ in (??). The idea is that the weight should satisfy $g(s) \to 0, s \to \infty$
<figure><img src="image_2.png" /><figcaption>Fig. 1: Synthetic test image with 20 standard deviations of noise (left) and the obtained results for TV (green/thick) and the $\Delta_{(p(x),q(x))}$-operator (red/dashed). We see that both error measures improve with our operator and we get less staircasing artifacts while preserving corner points and edges as shown in the detailed images.</figcaption></figure>

and $g(s) \to 1$, $s \to 0^+$. P-M studied weights like $g = k/(k+s)$ and $g = e^{-ks}$ and demonstrated the advantages of these weight functions for edge preservations. Rewriting the operator given by the P-M method in Gauge coordinates, we obtain

$$ \mathrm{div}(g(|\nabla u|^2)\nabla u) = g(|\nabla u|^2)u_{\xi\xi} + \phi(|\nabla u|^2)u_{\eta\eta}, \quad (7) $$

where $\phi(s) = (sg(s))'$. Thus the diffusion in the direction $\eta$ differs from the diffusion in the direction $\xi$. Since $\phi(s)$ is negative for large $s$, the evolution will be of backward diffusion effect near edges. This backward evolution cause problem for the well-posedness of the model and could also lead to staircasing problem [?].

The Perona-Malik PDE is a forward-backward type equation, and the diffusion is forward in the region $\{|\nabla u| < k\}$ and is backward and hence ill-posed in the region $\{|\nabla u| > k\}$. From evolution point of view, the forward-backward type equation may have infinitely many solutions and from variational minimization perspective the energy functional may have infinitely many minima [?]. To overcome the ill-posedness in the P-M model the authors in [?] introduced a regularization and proposed the following model

$$ \partial_t u = \mathrm{div}(g(|\nabla G_{\sigma} * u|)\nabla u), $$

where $G_{\sigma}$ is a Gaussian function with standard deviation $\sigma$. A similar but time dependent variance $\sigma(t)$ was used in [?]. However, it's rather tricky to choose $\sigma(t)$ since it should neither decay too fast nor too slow during the evolution.
# 3 Extending the p-Laplace Operator

## 3.1 Anisotropic decomposition via constant coefficients

The p-Laplace operator in (??), is an isotropic operator and expanding the divergence we obtain the equivalent form

$$ \Delta_p u = \partial_x ([(\partial_x u)^2 + (\partial_y u)^2]^{\frac{p-2}{2}} \partial_x u) + \partial_y ([(\partial_x u)^2 + (\partial_y u)^2]^{\frac{p-2}{2}} \partial_y u). \quad (8) $$

An anisotropic behavior of the $\Delta_p$-operator is induced by suppressing mixed derivatives, i.e., we define

$$ L_{(p,p)} u = \partial_x (|\partial_x u|^{p-2} \partial_x u) + \partial_y (|\partial_y u|^{p-2} \partial_y u), \quad 1 \le p < \infty. \quad (9) $$

In the case $p = 1$, (??) is known as isotropic TV and (??) is anisotropic TV [?]. By decoupling the exponents in (??) one obtains the operator $L_{(p_1,p_2)}$

$$ L_{(p_1,p_2)} u = \partial_x (|\partial_x u|^{p_1-2} \partial_x u) + \partial_y (|\partial_y u|^{p_2-2} \partial_y u), \quad 1 \le p_1, p_2 < \infty. \quad (10) $$

which has previously appeared in fluid mechanics and we refer to [?], [?].

Next, to see how the diffusion appears in the operator (??) we reformulate it in Gauge coordinates (??) by making the following definition.

**Definition 1.** The $L_{(p_1,p_2)}$-operator is given by

$$ L_{(p_1,p_2)} u = \partial_{\xi} (|\partial_{\xi} u|^{p_1-2} \partial_{\xi} u) + \partial_{\eta} (|\partial_{\eta} u|^{p_2-2} \partial_{\eta} u), \quad (11) $$

where $1 \le p_1, p_2 < \infty$.

The above operator $L_{(p_1,p_2)}$ is in fact a generalization of several known operators. We have

1. The case $p_1 = p_2 = 2$. The operator $L_{(2,2)}$ in (??) is the Laplacian, by now well studied. Due to the Laplacian's rotation invariance property we get

$$ L_{(2,2)} u = \partial_{\xi} (\partial_{\xi} u) + \partial_{\eta} (\partial_{\eta} u) = \Delta u = u_{xx} + u_{yy}. \quad (12) $$

2. The case $p_1 = 2$ and $p_2 = 1$. The operator $L_{(p_1,p_2)}$ in (??) is then given by

$$ L_{(2,1)} u = \partial_{\xi} (\partial_{\xi} u) + \partial_{\eta} (|\partial_{\eta} u|^{-1} \partial_{\eta} u). \quad (13) $$

Since we have $|\partial_{\eta} u|^{-1} \partial_{\eta} u = 1$, it follows that $L_{(2,1)} u = u_{\xi\xi}$. In Cartesian coordinates this corresponds to

$$ L_{(2,1)} u = |\nabla u| \Delta_1 u, \quad (14) $$

i.e. the mean curvature equation (see e.g., [?]). The corresponding regularized (weighted) mean curvature equation is given by

$$ \partial_t u = g(|\nabla G_{\sigma} * u|) |\nabla u| \Delta_1 u \quad (15) $$

previously studied in the context of image analysis, see e.g., [?].
3. The case $p_1 = 2$ and $p_2 = p \in (1, 2)$. It follows from (\ref{eq:operator}) that

$$L_{(2,p)}u = \partial_{\xi\xi}u + \partial_{\eta}(|\partial_{\eta}u|^{p-2}\partial_{\eta}u) = \partial_{\xi\xi}u + (p-1)u_{\eta}^{p-2}u_{\eta\eta} \quad (16a)$$

$$= |\nabla u| \Delta_1 u + (p-1) |\nabla u|^{p-2} \Delta_{\infty} u \quad (16b)$$

i.e. a mean curvature operator (\ref{eq:operator}) with a second order term corresponding to (\ref{eq:operator}). Note that the second order term induce invariant smoothing of the image data. Since (\ref{eq:operator}) is merely a special case of (\ref{eq:operator}), we further relax the mean curvature term next to better reflect the trade-off between edge-preservation and obtained smoothness. Although this modification appears straight-forward, its implications are non-trivial, however.

The anisotropic $L_{(2,p)}$ operator is given by

**Definition 2.** The operator $\Delta_{(p,q)}$ is

$$\Delta_{(p,q)}u = |\nabla u|^q \Delta_1 u + (p-1)|\nabla u|^{p-2} \Delta_{\infty} u, \quad p \in [1, 2], \quad q \ge 0, \quad (17)$$

*Remark 1.* Definition \ref{eq:operator} is a straight-forward relaxation of the exponents in (\ref{eq:operator}), motivated by the discussion in above point 3. The introduction of $q$ in (\ref{eq:operator}), defines an additional degree of freedom, allowing us to control the trade-off between $\Delta_1$ and $\Delta_{\infty}$, i.e., a trade-off between edge preservation and smoothness.

The corresponding evolution problem of the operator $\Delta_{(p,q)}$ is

**Definition 3.** The evolution problem of $\Delta_{(p,q)}$ is given by

$$\left\{ \begin{array}{l} \partial_t u - \Delta_{(p,q)} u = 0, \quad p \in [1, 2], \quad q \ge 0 \\ u(0) = u_0 \\ \partial_n u = 0 \end{array} \right. \quad (18)$$

We point out that $q = 0, p = 1$ yields the familiar isotropic TV regularizer and $q = 1, p = 2$ results in the heat equation, i.e., isotropic filtering. In the next section, we propose to couple $p$ and $q$ via a smooth non-linear transition function such that $\Delta_{(p,q)}$ can be thought of a spatially variant TV and isotropic regularizer.

# 4 Variable coefficients

## 4.1 Coefficient coupling

In the previous section we motivated the $\Delta_{(p,q)}$ operator and advocated to introduce variable coefficients, depending on the image data. The behavior of the operator $\Delta_{(p(x),q(x))}$ that we seek, is the edge-preserving effect of TV leading to $\Delta_{(p(x),q(x))} \to \Delta_{(1,0)}$ (the case of TV) as $|\nabla u| \to \infty$. In regions with small gradients we define the $(p(x),q(x))$-coefficients such that the operator is a linear filter, leading to $\Delta_{(p(x),q(x))} \to \Delta_{(2,1)} = \Delta$ (the case of isotropic diffusion) as $|\nabla u| \to 0$. To see how the diffusion appears in $\Delta_{(p(x),q(x))}$ we rewrite the operator in Gauge coordinates
**Lemma 1.** The operator $\Delta_{(p,q)}$ in (??) can be written as

$$ \Delta_{(p,q)} u = |\nabla u|^{q-1} \left( \xi^{\top} (D^2 u) \xi + (p-1) |\nabla u|^{p-q-1} \eta^{\top} (D^2 u) \eta \right). \quad (19) $$

*Proof.* Using the relations in (??) and by rotation invariance of Laplacian, we obtain

$$ \Delta u = u_{\xi\xi} + u_{\eta\eta} = \Delta_{\infty} u + |\nabla u| \Delta_1 u. \quad (20) $$

Now we observe that the operator $\Delta_1 u$ can be expanded as

$$ \Delta_1 u = \frac{u_y^2 u_{xx} - 2u_x u_y u_{xy} + u_x^2 u_{yy}}{|\nabla u|^3} = \frac{(\nabla^{\perp} u)^{\top} (D^2 u) \nabla^{\perp} u}{|\nabla u|^3}, \quad (21) $$

and using this in (??) gives

$$ \Delta_{\infty} u = |\nabla u|^2 \Delta u - (\nabla^{\perp} u)^{\top} (D^2 u) \nabla^{\perp} u = \nabla^{\top} u (D^2 u) \nabla u. \quad (22) $$

Thus, the $\Delta_{(p,q)}$ operator in (??) reformulates to

$$ \begin{aligned} \Delta_{(p,q)} u &= |\nabla u|^{q-3} (\nabla^{\perp} u)^{\top} (D^2 u) \nabla^{\perp} u + (p-1) |\nabla u|^{p-4} \nabla^{\perp} u (D^2 u) \nabla u \\ &= |\nabla u|^{q-1} \xi^{\top} (D^2 u) \xi + (p-1) |\nabla u|^{p-2} \eta^{\top} (D^2 u) \eta, \end{aligned} $$

which shows the result. $\square$

Next, we couple $p$ and $q$ via the relation

$$ p(|\nabla u|) = 1 + q(|\nabla u|) \quad (23) $$

and from which we derive the following result

**Lemma 2.** If $p(x) = 1 + q(x)$, then

$$ \Delta_{(1+q(x),q(x))} u = |\nabla u|^{q(x)-1} \left( \xi^{\top} (D^2 u) \xi + q(x) \eta^{\top} (D^2 u) \eta \right). \quad (24) $$

*Proof.* The proof follows immediately from Lemma ??.
$\square$

In this study, we couple $p$ and $q$ via the negative exponential function (although other selections are possible), i.e., we set

$$ q(|\nabla u|) = k_2 \exp(-|\nabla u|/k_1), \quad (25a) $$

$$ p(|\nabla u|) = 1 + q(|\nabla u|), \quad (25b) $$

where $k_1 > 0$ and $0 < k_2 < 1$. Under this selection we form the following parabolic PDE

$$ \partial_t u - |\nabla u|^{q(|\nabla u|)-1} \left( \xi^{\top} (D^2 u) \xi + q(|\nabla u|) \eta^{\top} (D^2 u) \eta \right) = 0. \quad (26) $$
One easily checks that (??) describes a smooth transition between total variation and linear filtering for the selection of p, q in (??). By using $|\nabla u|\Delta_1 u = \Delta u - \Delta_\infty u$, the operator $\Delta_{(1+q,q)}$ becomes

$$ \Delta_{(1+q,q)} u = |\nabla u|^{q(|\nabla u|)-1} \left( \mathrm{Tr}(D^2 u) + \frac{q(|\nabla u|)-1}{|\nabla u|^2} \nabla u^\top (D^2 u) \nabla u \right). \quad (27) $$

Remark 2. The operator $\Delta_{(1+q,q)}$ is non-linear and have unbounded coefficients. In this first study, we perturb the $\Delta_{(1+q,q)}$-operator to obtain a regularized version of the evolution problem (??). This regularization enables us pose the necessary conditions for well-posedness, introduced next.

## 4.2 Regularity of solutions

In order to set up the framework for numerical calculation, we define

**Definition 4.** Let $0 < \varepsilon, \delta < 1$ and $q$ as in (??). Then

$$ \Delta_{(1+q,q)}^{\varepsilon, \delta} u = (|\nabla u|^2 + \varepsilon^2)^{(q(|\nabla u|)-1)/2} \left( (1+\delta)\mathrm{Tr}(D^2 u) + \frac{q(|\nabla u|)-1}{|\nabla u|^2 + \varepsilon^2} \nabla u^\top (D^2 u) \nabla u \right). \quad (28) $$

Thus we study the regularized evolution equation

$$ u_t = \Delta_{(1+q,q)}^{\varepsilon, \delta} u = \sum_{i,j=1} a_{ij}^{\varepsilon, \delta}(\nabla u) u_{ij} \quad (29) $$

where

$$ a_{ij}^{\varepsilon, \delta}(\zeta) = (|\zeta|^2 + \varepsilon^2)^{(q(|\zeta|)-1)/2} \left[ (1+\delta)\delta_{ij} + \frac{q(|\zeta|)-1}{|\zeta|^2 + \varepsilon^2} \zeta_i \zeta_j \right], \quad (30) $$

and $\zeta = \nabla u$. Given $u^{(k)}$, we obtain the next update $u^{(k+1)}$ by solving the following, equivalent (see [?]), initial value problem iteratively

$$ \left\{ \begin{array}{ll} u_t^{(k+1)} = \sum a_{ij}^{\varepsilon, \delta}(\nabla u^{(k)}) u_{ij}^{(k+1)}, & \text{in } \Omega \times (0, T) \\ u^{(k+1)}(0) = u_0, & \text{in } \Omega \\ \partial_n u^{(k+1)} = 0, & \text{on } \partial\Omega \times (0, T) \end{array} \right. \quad (31) $$

**Proposition 1.** If the initial data $u^{(0)} = 0$ then the solution $u^{(k)}$ to (??) exists and is in $C^\infty(\Omega \times (0, T))$.

*Proof.* If the initial guess $u^{(0)} = 0$, then $\sum a_{ij}^{\varepsilon, \delta}(0)u_{ij} = (1+\delta)\varepsilon^{q(0)-1}\Delta u$ and we deduce Thm 4.31 in [?] that $u^{(1)}$ exists and is $C^\infty(\Omega \times (0, T))$. Given $u^{(k)} \in C^\infty$, then $||u^{(k)}||_{C^1(\Omega \times (0, T))}$ is bounded, $(|\zeta|^2 + \varepsilon^2)^{(q(|\zeta|)-1)/2}$ is bounded from below and $a_{ij}^{\varepsilon, \delta}$ are also $C^\infty(\Omega \times (0, T))$. Hence, there are constants $c = c(\varepsilon, \delta, k)$ and $C = C(\varepsilon, \delta, k) > 0$, depending only on $\varepsilon, \delta, ||u^{(k)}||_{C^1}$, such that

$$ c|\zeta|^2 \le a_{ij}^{\varepsilon, \delta}(\nabla u^{(k)})\zeta_i \zeta_j \le C|\zeta|^2. $$

It follows once again from Thm 4.31 in [?] that $u^{(k+1)}$ exists and also is $C^\infty(\Omega \times (0, T))$. □
<figure><img src="image_3.png" /><figcaption>Fig. 2: Example results for the $\Delta_{(p,q)}$-operator where we stopped the filtering process at the maximum SSIM value. “TV”-was obtained using the Split Bregman method [?]. We set $k_1 = k_2 = 0.1$ in $\Delta_{(p(x),q(x))}$. See text for details.</figcaption></figure>

# 5 Evaluation

## 5.1 Implementation

The parabolic PDE is discretized by using finite differences and a simple forward Euler scheme. We used a state of the art Split Bregman (SB) implementation of total variation. For details on SB see [?]. The regularization parameter for the SB was optimized in the range [0.01, 1] in increasing steps of 0.11. The regularization parameter of the Bregman-variables of SB was set to 1 and the scheme was terminated as $10^{-3} > ||u^{(k)} - u^{(k-1)}||_2/||u^{(k-1)}||_2$ and we choose the regularization parameter that produces highest SSIM value [?]. We also report the peak signal-to-noise value (PSNR). For the $(p(x), q(x))$-operator we found that $k_1 = k_2 = 0.1$ and the update stepsize as $\alpha = 10^{-5}$ and $\tau = 0.5$ works well for the considered noise level of 20 standard deviations of additive Gaussian noise. These values are ad-hoc and future work include methods for parameter estimation.

## 5.2 Results

First we test our algorithm on a synthetic test image seen in figure ?? . We see that the result of the proposed operator appears smoother than the result of
<figure><img src="image_4.png" /><figcaption>Fig. 3: Example denoising a grayscale image with 20 standard deviations of noise. In the close-up images to the right it can be seen that TV (thick/green) produces the characteristic staircasing effect while the operator $\Delta_{(p,q)}$ (dashed/red) shows good visual similarity with the noise free patch. TV shows good result in the sky (up right), but oversmooths, e.g., the window tiles seen in the close-up down right.</figcaption></figure>

TV in the center region, but yet preserves the corner point well. In figure ?? (cropped 256 × 256 pixels of image 35049.jpg [?]) we compare the visual quality and the SSIM-values for a range of *p*, *q*-values. As expected for $\Delta_{(2,1)}$ (isotropic filtering) performs the worst whereas the operator produces improved result w.r.t. SSIM as well as perceptual appearance. In the case of non-adaptive parameters, TV performs the best. However, the operator with adaptive coefficients produces both SSIM, PSNR values and less oversmoothing (down-right figure). We also include the result from a grayscale image “Castle” (cropped 256 × 256 pixels of image 102061.jpg [?]) in figure ?? . In this image TV performs very well in the sky (detail up right with green/thick frame) whereas the result from the $\Delta_{(p,q)}$ operator appears less noisy and looks visually more crisp. In both examples the proposed operator shows an improvement in PSNR as well as SSIM values.
# 6 Conclusion

In this paper we introduced a new family of operators, $\Delta_{(p,q)}$. Preliminary numerical results indicate that there could be a relationship between $p(x)$ and $q(x)$ that further improves the restoration effect. In forthcoming works we will investigate the operator $\Delta_{(p,q)}$ further regarding both regularity and different areas of applications.

# References

1. Alvarez, L., Lions, P.L., Morel, J.M.: Image selective smoothing and edge detection by nonlinear diffusion. SIAM J. Numer. Anal. 29(3), 845–866 (1992)

2. Antontsev, S.N., Shmarev, S.I.: Existence and uniqueness of solutions of degenerate parabolic equations with variable exponents of nonlinearity. Fundam. Prikl. Mat. 12(4), 3–19 (2006)

3. Antontsev, S., Shmarev, S.: Elliptic equations with anisotropic nonlinearity and nonstandard growth conditions. Handbook of differential equations. stationary partial differential equations. Amsterdam: Elsevier/North Holland 3, 1–100 (2006)

4. Aubert, G., Kornprobst, P.: Mathematical problems in image processing, Applied Mathematical Sciences, vol. 147. Springer, New York, second edn. (2006)

5. Baravdish, G., Svensson, O., Åström, F.: On Backward $p(x)$-Parabolic Equations for Image Enhancement. Numer. Funct. Anal. Optim. 36(2), 147–168 (2015)

6. Buades, A., Coll, B., Morel, J.M.: The staircasing effect in neighborhood filters and its solution. Image Processing, IEEE Trans. on 15(6), 1499–1505 (2006)

7. Chen, Y.G., Giga, Y., Goto, S.: Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations. J. Diff. Geom. 33(3), 749–786 (1991)

8. Does, K.: An evolution equation involving the normalized $p$-laplacian. Commun. Pure Appl. Anal. 10(1), 361–396 (2011)

9. Goldstein, T., Osher, S.: The split bregman method for 11-regularized problems. SIAM Journal on Imaging Sciences 2(2), 323–343 (2009)

10. Keeling, S.L., Stollberger, R.: Nonlinear anisotropic diffusion filtering for multiscale edge enhancement. Inverse Problems 18(1), 175–190 (2002)

11. Lieberman, G.M.: Second order parabolic differential equations. World Scientific Publishing Co. Inc., River Edge, NJ (1996)

12. Martin, D., Fowlkes, C., Tal, D., Malik, J.: A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In: ICCV. vol. 2, pp. 416–423 (2001)

13. Perona, P., Malik, J.: Scale-space and edge detection using anisotropic diffusion. IEEE Transactions, PAMI 12, 629–639 (1990)

14. Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. Physica D 60(1-4), 259–268 (1992)

15. Wang, Z., Bovik, A., Sheikh, H., Simoncelli, E.: Image quality assessment: from error visibility to structural similarity. Image Processing, IEEE Transactions on 13(4), 600–612 (2004)

16. Vogel and Oman: Iterative methods for total variation denoising, SIAM J. Sci. Comput. 17, pp.227-238, (1996)

17. Whitaker, R.T., Pizer, S.M.: A multi-scale approach to nonuniform diffusion. CVGIP: Image Understanding 57(1), 99–110 (1993)

18. You, Y.L., Xu, W., Tannenbaum, A., Kaveh, M.: Behavioral analysis of anisotropic diffusion in image processing. Image Processing, IEEE Trans. on 5(11), 1539–1553 (1996)