1 Introduction

A well known inverse problem in image processing is image denoising [4]. 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

$$\begin{aligned} \varDelta _pu = \,\mathrm {div}(|\nabla u|^{p-2}\nabla u),\quad p\ge 1, \end{aligned}$$
(1)

associated with the evolution equation of p-Laplacian

$$\begin{aligned} \left\{ \begin{array}{rcll} \partial _t u - \varDelta _p u &{}=&{} 0,&{}\text{ in }\ \varOmega \times (0,T)\\ u(0) &{}=&{} u_0,&{}\text{ in }\ \varOmega \\ \partial _nu&{}=&{}0,&{}\text{ on }\ \partial \varOmega \times (0,T) \end{array}\right. \end{aligned}$$
(2)

where \(\displaystyle \varOmega \) is a bounded domain in \(\displaystyle \mathbf{R}^2\) and \(\displaystyle u_0:\varOmega \rightarrow \mathbf{R}\) is a given degraded image [5, 7, 13] and \(\displaystyle \nabla u\) is the gradient. The degenerate parabolic Eq. in (2) has been studied by many authors and we limit ourselves here to refer the reader to [10]. It is well known that the case \(\displaystyle p=2\) gives the linear Gaussian filter, which however, impose strong spatial regularity and therefore image details such as lines and edges are oversmoothed. \(\displaystyle p=1\) is often refereed to as the method of total variation [13] and \(\displaystyle p=0\) is an instance of the so called balanced forward backward evolution [9].

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

In Sect. 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 Sect. 3 where we extend the p-Laplace operator to a new operator \(\displaystyle \varDelta _{(p,q)}\) with focus on image denoising. We consider the corresponding variable version \(\displaystyle \varDelta _{(p(x),q(x))}\) in Sect. 4. Finally, in Sect. 5 we demonstrates the applicability of \(\displaystyle \varDelta _{(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 (1) as

$$\begin{aligned} \varDelta _pu= & {} |\nabla u|^{p-1}\varDelta _1 u +(p-1) |\nabla u|^{p-2}\varDelta _{\infty } u, \end{aligned}$$
(3)

where \(\displaystyle \varDelta _1 u= \,\mathrm {div}\Big (\frac{\nabla u}{|\nabla u|}\Big )\), \(\displaystyle \varDelta _{\infty }u = \frac{\nabla u }{|\nabla u|}\cdot (D^2u) \frac{\nabla u }{|\nabla u|}\) and \(\displaystyle D^2u\) is the Hessian of u. However, an intuitive way to represent \(\displaystyle \varDelta _{p}\), giving direct interpretation of the diffusivity directions is to express \(\displaystyle \varDelta _{p}\) by using Gauge coordinates \(\displaystyle (x, y) \rightarrow (\eta ,\xi )\):

$$\begin{aligned} \varDelta _pu= & {} |\nabla u|^{p-2} (u_{\xi \xi } +(p-1) u_{\eta \eta }) \end{aligned}$$
(4)

where

$$\begin{aligned} \eta = \frac{\nabla u}{|\nabla u|},\qquad \xi = \frac{\nabla ^{\perp } u}{|\nabla u|}, \end{aligned}$$
(5)

and

$$\begin{aligned} u_{\xi \xi } = |\nabla u| \varDelta _1 u,\qquad u_{\eta \eta } = \varDelta _{\infty } u. \end{aligned}$$
(6)

From (4) it is now clear that \(\displaystyle \varDelta _p\) imposes the same diffusivity strength in both directions \(\displaystyle \xi \) and \(\displaystyle \eta \) independent of the magnitude of the gradient. In an attempt to resolve this drawback, Perona and Malik (P-M) [12] proposed to replace \(\displaystyle |\nabla u|^{p-2}\) with \(\displaystyle g(|\nabla u|^2)\) in (1). The idea is that the weight should satisfy \(\displaystyle g(s)\rightarrow 0\), \(\displaystyle s\rightarrow \infty \) and \(\displaystyle g(s)\rightarrow 1\), \(\displaystyle s\rightarrow 0^+\). P-M studied weights like \(\displaystyle g=k/(k+s)\) and \(\displaystyle 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

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

where \(\displaystyle \phi (s)=(sg(s))'\). Thus the diffusion in the direction \(\displaystyle \eta \) differs from the diffusion in the direction \(\displaystyle \xi \). Since \(\displaystyle \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 [17].

The Perona-Malik PDE is a forward-backward type equation, and the diffusion is forward in the region \(\displaystyle \{ |\nabla u | < k \}\) and is backward and hence ill-posed in the region \(\displaystyle \{ |\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 [17]. To overcome the ill-posedness in the P-M model the authors in [17] introduced a regularization and proposed the following model

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

where \(\displaystyle G_{\sigma }\) is a Gaussian function with standard deviation \(\displaystyle \sigma \). A similar but time dependent variance \(\displaystyle \sigma (t)\) was used in [16]. However, it’s rather tricky to choose \(\displaystyle \sigma (t)\) since it should neither decay too fast nor too slow during the evolution.

Fig. 1.
figure 1

Synthetic test image with 20 standard deviations of noise (left) and the obtained results for TV (green/thick) and the \(\displaystyle \varDelta _{(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. (Color figure online)

Full size image

3 Extending the p-Laplace Operator

3.1 Anisotropic Decomposition via Constant Coefficients

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

$$\begin{aligned} \varDelta _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). \end{aligned}$$
(8)

An anisotropic behavior of the \(\displaystyle \varDelta _p\)-operator is induced by suppressing mixed derivatives, i.e., we define

$$\begin{aligned} L_{(p,p)} u =\partial _x(|\partial _x u|^{p-2} \partial _x u) + \partial _y(|\partial _y u|^{p-2} \partial _y u), \qquad 1\le p<\infty . \end{aligned}$$
(9)

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

$$\begin{aligned} 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 , \end{aligned}$$
(10)

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

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

Definition 1

The \(\displaystyle L_{(p_1,p_2)}\)-operator is given by

$$\begin{aligned} 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), \end{aligned}$$
(11)

where \(\displaystyle 1\le p_1,p_2<\infty \).

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

  1. 1.

    The case \(\displaystyle p_1=p_2=2\). The operator \(\displaystyle L_{(2,2)}\) in (11) is the Laplacian, by now well studied. Due to the Laplacian’s rotation invariance property we get

    $$\begin{aligned} L_{(2,2)} u =\partial _{\xi }(\partial _{\xi } u) + \partial _{\eta }(\partial _{\eta } u) = \varDelta u = u_{xx}+u_{yy}. \end{aligned}$$
    (12)
  2. 2.

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

    $$\begin{aligned} L_{(2,1)} u = \partial _{\xi }( \partial _{\xi } u) + \partial _{\eta }(|\partial _{\eta } u|^{-1} \partial _{\eta } u). \end{aligned}$$
    (13)

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

    $$\begin{aligned} L_{(2,1)} u = |\nabla u| \varDelta _1 u, \end{aligned}$$
    (14)

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

    $$\begin{aligned} \partial _t u = g(|\nabla G_{\sigma } * u|)|\nabla u| \varDelta _1 u \end{aligned}$$
    (15)

    previously studied in the context of image analysis, see e.g., [1].

  3. 3.

    The case \(\displaystyle p_1=2\) and \(\displaystyle p_2=p\in (1,2)\). It follows from (11) that

    $$\begin{aligned} 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 } \end{aligned}$$
    (16a)
    $$\begin{aligned}&= |\nabla u| \varDelta _1 u + (p-1) |\nabla u|^{p-2}\varDelta _{\infty }u \end{aligned}$$
    (16b)

    i.e. a mean curvature operator (14) with a second order term corresponding to (3). Note that the second order term induce invariant smoothing of the image data. Since (16) is merely a special case of (3), 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 \(\displaystyle L_{(2,p)}\) operator is given by

Definition 2

The operator \(\displaystyle \varDelta _{(p,q)}\) is

$$\begin{aligned} \varDelta _{(p,q)} u = |\nabla u|^q \varDelta _1 u + (p-1) |\nabla u|^{p-2}\varDelta _{\infty }u,\quad p \in [1, 2],\quad q\ge 0. \end{aligned}$$
(17)

Remark 1

Definition 2 is a straight-forward relaxation of the exponents in (3), motivated by the discussion in above point 3. The introduction of q in (17), defines an additional degree of freedom, allowing us to control the trade-off between \(\displaystyle \varDelta _1\) and \(\displaystyle \varDelta _\infty \), i.e., a trade-off between edge preservation and smoothness.

The corresponding evolution problem of the operator \(\displaystyle \varDelta _{(p,q)}\) is

Definition 3

The evolution problem of \(\displaystyle \varDelta _{(p,q)}\) is given by

$$\begin{aligned} \left\{ \begin{array}{rcl} \partial _t u - \varDelta _{(p,q)} u &{}=&{}0,\quad p\in [1, 2],\quad q\ge 0\\ u(0) &{}=&{} u_0 \\ \partial _n u &{}=&{}0 \end{array}\right. \end{aligned}$$
(18)

We point out that \(\displaystyle q=0, p=1\) yields the familiar isotropic TV regularizer and \(\displaystyle 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 \(\displaystyle \varDelta _{(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 \(\displaystyle \varDelta _{(p,q)}\) operator and advocated to introduce variable coefficients, depending on the image data. The behavior of the operator \(\displaystyle \varDelta _{(p(x),q(x))}\) that we seek, is the edge-preserving effect of TV leading to \(\displaystyle \varDelta _{(p(x),q(x))}\rightarrow \varDelta _{(1,0)}\) (the case of TV) as \(\displaystyle |\nabla u|\rightarrow \infty \). In regions with small gradients we define the (p(x), q(x))-coefficients such that the operator is a linear filter leading to \(\displaystyle \varDelta _{(p(x),q(x))}\rightarrow \varDelta _{(2,1)}=\varDelta \) (the case of isotropic diffusion) as \(\displaystyle |\nabla u|\rightarrow 0\). To see how the diffusion appears in \(\displaystyle \varDelta _{(p(x),q(x))}\) we rewrite the operator in Gauge coordinates.

Lemma 1

The operator \(\displaystyle \varDelta _{(p,q)}\) in (17) can be written as

$$\begin{aligned} \varDelta _{(p,q)} u = |\nabla u|^{q-1} \bigg (\xi ^\top (D^2u) \xi + (p-1)|\nabla u|^{p-q-1} \eta ^\top (D^2u) \eta \bigg ). \end{aligned}$$
(19)

Proof

Using the relations in (6) and by rotation invariance of Laplacian, we obtain

$$\begin{aligned} \varDelta u = u_{\xi \xi }+u_{\eta \eta } = \varDelta _{\infty } u + |\nabla u| \varDelta _1 u. \end{aligned}$$
(20)

Now we observe that the operator \(\displaystyle \varDelta _1 u\) can be expanded as

$$\begin{aligned} \varDelta _1 u = \frac{u_y^2 u_{xx} - 2u_xu_yu_{xy} + u_x^2u_{yy}}{|\nabla u|^3} = \frac{(\nabla ^\perp u)^\top (D^2u) \nabla ^\perp u}{|\nabla u|^3}, \end{aligned}$$
(21)

and using this in (20) gives

$$\begin{aligned} \varDelta _{\infty } u = |\nabla u|^{2}\varDelta u - (\nabla ^\perp u)^\top (D^2u) \nabla ^\perp u&= \nabla ^\top u (D^2u) \nabla u. \end{aligned}$$
(22)

Thus, the \(\displaystyle \varDelta _{(p,q)}\) operator in (17) reformulates to

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

which shows the result.    \(\displaystyle \square \)

Next, we couple p and q via the relation

$$\begin{aligned} p(|\nabla u|) = 1 + q(|\nabla u|) \end{aligned}$$
(23)

and from which we derive the following result.

Lemma 2

If \(\displaystyle p(x)=1+q(x)\), then

$$\begin{aligned} \varDelta _{(1+q(x),q(x))} u = |\nabla u|^{q(x)-1} \bigg (\xi ^\top (D^2u) \xi + q(x) \eta ^\top (D^2u) \eta \bigg ). \end{aligned}$$
(24)

Proof

The proof follows immediately from Lemma 1.    \(\displaystyle \square \)

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

$$\begin{aligned} q(|\nabla u|)&= k_2 \exp \left( - |\nabla u| / k_1 \right) , \end{aligned}$$
(25a)
$$\begin{aligned} p(|\nabla u|)&= 1 + q(|\nabla u|), \end{aligned}$$
(25b)

where \(\displaystyle k_1>0\) and \(\displaystyle 0<k_2 <1\). Under this selection we form the following parabolic PDE

$$\begin{aligned} \partial _t u - |\nabla u|^{q(|\nabla u|)-1} \bigg (\xi ^\top (D^2u) \xi + q(|\nabla u|) \eta ^\top (D^2u) \eta \bigg ) = 0. \end{aligned}$$
(26)

One easily checks that (26) describes a smooth transition between total variation and linear filtering for the selection of p, q in (). By using \(\displaystyle |\nabla u|\varDelta _{1}u=\varDelta u-\varDelta _{\infty }u\), the operator \(\displaystyle \varDelta _{(1+q,q)}\) becomes

$$\begin{aligned} \varDelta _{(1+q,q)} u = |\nabla u|^{q(|\nabla u|)-1} \bigg (Tr(D^2u) + \frac{q(|\nabla u|)-1}{|\nabla u|^{2}} \nabla u^\top (D^2u) \nabla u \bigg ). \end{aligned}$$
(27)

Remark 2

The operator \(\displaystyle \varDelta _{(1+q,q)}\) is non-linear and have unbounded coefficients. In this first study, we perturb the \(\displaystyle \varDelta _{(1+q,q)}\)-operator to obtain a regularized version of the evolution problem (18). This regularization enables us to 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 \(\displaystyle 0<\varepsilon ,\delta <1\) and q as in (25). Then

$$\begin{aligned}&\varDelta ^{\varepsilon ,\delta }_{(1+q,q)} u = \nonumber \\&(|\nabla u|^{2}+\varepsilon ^{2})^{(q(|\nabla u|)-1)/2} \bigg ((1+\delta )Tr(D^2u) + \frac{q(|\nabla u|)-1}{|\nabla u|^{2}+\varepsilon ^{2}} \nabla u^\top (D^2u) \nabla u \bigg ). \end{aligned}$$
(28)

Thus we study the regularized evolution equation

$$\begin{aligned} u_{t} = \varDelta ^{\varepsilon ,\delta }_{(1+q,q)} u = \sum _{i,j=1} a_{ij}^{\varepsilon ,\delta }(\nabla u)u_{ij} \end{aligned}$$
(29)

where

$$\begin{aligned} a_{ij}^{\varepsilon ,\delta }(\zeta ) = (|\zeta |^{2}+\varepsilon ^{2})^{(q(|\zeta |)-1)/2} \Big [(1+\delta )\delta _{ij}+ \frac{q(|\zeta |)-1}{|\zeta |^{2}+\varepsilon ^{2}}\zeta _{i} \zeta _{j}\Big ], \end{aligned}$$
(30)

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

$$\begin{aligned} \left\{ \begin{array}{rclllll} u^{(k+1)}_{t} &{}=&{} \sum a_{ij}^{\varepsilon ,\delta }(\nabla u^{(k)})u^{(k+1)}_{ij}, &{}\text{ in }\ \varOmega \times (0,T)\\ u^{(k+1)}(0) &{}=&{} u_0,&{}\text{ in }\ \varOmega \\ \partial _n u^{(k+1)}&{}=&{}0,&{}\text{ on }\ \partial \varOmega \times (0,T)\\ \end{array}\right. \end{aligned}$$
(31)

Proposition 1

If the initial data \(\displaystyle u^{(0)} =0\) then the solution \(\displaystyle u^{(k)}\) to (31) exists and is in \(\displaystyle C^{\infty }(\varOmega \times (0,T))\).

Proof

If the initial guess \(\displaystyle u^{(0)} =0, \) then \(\displaystyle \sum a_{ij}^{\varepsilon ,\delta }(0)u_{ij} =(1+\delta ) \varepsilon ^{q(0)-1}\varDelta u\) and we deduce from Theorem 4.31 in [10] that \(\displaystyle u^{(1)} \) exists and is \(\displaystyle C^{\infty }(\varOmega \times (0,T))\). Given \(\displaystyle u^{(k)} \in C^{\infty }\), then \(\displaystyle ||u^{(k)}||_{C^1(\varOmega \times (0,T))}\) is bounded, \(\displaystyle (|\zeta |^{2}+\varepsilon ^{2})^{(q(|\zeta |)-1)/2}\) is bounded from below and \(\displaystyle a_{ij}^{\varepsilon , \delta }\) are also \(\displaystyle C^{\infty }(\varOmega \times (0,T))\). Hence, there are constants \(\displaystyle c=c(\varepsilon , \delta , k)\) and \(\displaystyle C=C(\varepsilon , \delta , k) >0\), depending only on \(\displaystyle \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 Theorem 4.31 in [10] that \(\displaystyle u^{(k+1)}\) exists and also is \(\displaystyle C^{\infty }(\varOmega \times (0,T))\).    \(\displaystyle \square \)

Fig. 2.
figure 2

Example results for the \(\displaystyle \varDelta _{(p,q)}\)-operator where we stopped the filtering process at the maximum SSIM value. “TV”-was obtained using the Split Bregman method [8]. We set \(\displaystyle k_1 = k_2 = 0.1\) in \(\displaystyle \varDelta _{(p(x),q(x))}\). See text for details.

Full size image

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 [8]. 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 \(\displaystyle 10^{-3} > ||u^{(k)}-u^{(k-1)}||_2/||u^{(k-1)}||_2\) and we choose the regularization parameter that produces highest SSIM value [14]. We also report the peak signal-to-noise value (PSNR). For the (p(x), q(x))-operator we found that \(\displaystyle k_1=k_2=0.1\) and the update stepsize as \(\displaystyle \alpha =10^{-5}\) and \(\displaystyle \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 Fig. 1. We see that the result of the proposed operator appears smoother than the result of TV in the center region, but yet preserves the corner point well. In Fig. 2 (cropped \(\displaystyle 256\times 256\) pixels of image 35049.jpg [11]) we compare the visual quality and the SSIM-values for a range of p, q-values. As expected for \(\displaystyle \varDelta _{(2,1)}\) (isotropic filtering) performs the worst whereas the operator produce 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 improve both SSIM, PSNR values and produces less oversmoothing (down-right figure). We also include the result from a grayscale image “Castle” (cropped \(\displaystyle 256\times 256\) pixels of image 102061.jpg [11]) in Fig. 3. In this image TV performs very well in the sky (detail up right with green/thick frame) whereas the result from the \(\displaystyle \varDelta _{(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.

Fig. 3.
figure 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 \(\displaystyle \varDelta _{(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. (Color figure online)

Full size image

6 Conclusion

In this paper we introduced a new family of operators, \(\displaystyle \varDelta _{(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 \(\displaystyle \varDelta _{(p,q)}\) further regarding both regularity and different areas of applications.