# Image Restoration Using Anisotropic Stochastic
Diffusion Collaborated with Non Local Means

Dariusz Borkowski, Katarzyna Jańczak-Borkowska

► To cite this version:

Dariusz Borkowski, Katarzyna Jańczak-Borkowska. Image Restoration Using Anisotropic Stochastic Diffusion Collaborated with Non Local Means. 12th International Conference on Information Systems and Industrial Management (CISIM), Sep 2013, Krakow, Poland. pp.177-189, 10.1007/978-3-642-40925-7_18 . hal-01496064

HAL Id: hal-01496064

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

Submitted on 27 Mar 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
# Image Restoration using Anisotropic Stochastic Diffusion Collaborated with Non Local Means*

Dariusz Borkowski<fnref n="1" /> and Katarzyna Jańczak-Borkowska<fnref n="2" />

**Abstract.** In this paper we explore the problem of the reconstruction of images with additive Gaussian noise. In order to solve this inverse problem we use stochastic differential equations with reflecting boundary and famous non local means algorithm. Expressing anisotropic diffusion in terms of stochastic equations allows us to adapt the concept of similarity patches used in non local means. This novel look on the reconstruction problem is fruitful, gives encouraging results and compares favourably with other image denoising filters.

**Keywords:** stochastic anisotropic diffusion, non local means

## 1 Introduction

Let $D$ be a bounded, convex domain in $\mathbb{R}^2$, $u: \bar{D} \to \mathbb{R}$ be an original image and $u_0: \bar{D} \to \mathbb{R}$ be the observed image of the form

$$u_0 = u + \eta,$$

where $\eta$ stands for a white Gaussian noise. We assume that $u$ and $u_0$ are appropriately regular. We are given $u_0$, the problem is to reconstruct $u$. This is a typical example of an inverse problem [2].

Problem of image denoising using fully automatic and reliable methods is one of the most important issues of digital image processing and computer vision. Efficient and effective reconstruction of images is an essential element of most image processing and recognising algorithms. Reconstruction algorithms allow us to make initial treatment of data for further analysis, which is very important, especially in astronomy, biology or medicine.

Various techniques were proposed to tackle this inverse problem. One may quote the linear filtering, DCT [31], wavelets theory [12, 15], variational methods [12, 23] and stochastic modelling which are generally based on the Markov field

<footnote n="1">¹ Faculty of Mathematics and Computer Science, Nicolaus Copernicus University, Chopina 12/18, 87-100 Toruń, Poland, dbor@mat.umk.pl</footnote>

<footnote n="2">² Institute of Mathematics and Physics, University of Technology and Life Sciences, al. prof. S. Kaliskiego 7, 85-789 Bydgoszcz, Poland, kaja@utp.edu.pl</footnote>

<footnote>* This research was supported by the National Science Centre (Poland) under decision number DEC-2012/07/D/ST6/02534</footnote>
theory and Bayesian approach [12, 16, 22]. The most important methods in the last decade in image processing have been methods driven by nonlinear diffusion equation [11, 21, 29, 30], where some papers [3–7, 17, 26, 27] involve advanced tools of stochastic analysis such as stochastic differential equations. In another class, one could include methods that take advantage of the non-local similarity of patches in the image. Among the most famous, we can name non local means (in short NL-means) [8, 9], BM3D [13, 14, 18], NL-Bayes [19] and K-SVD [1, 20].

In this paper we focus on two methods using to image denoising: anisotropic diffusion and non local means. We propose a method which combines these two approaches. Reconstructed pixel is expressed in terms of stochastic anisotropic diffusion which is itself driven by similarity of patches of the noisy image.

## 2 Mathematical Preliminaries

Let $D \subset \mathbb{R}^n$ be a domain with closure $\bar{D}$ and boundary $\partial D$. Let $T > 0$ and by $\mathbf{C}([0, T]; \mathbb{R}^n)$ we denote a set of continuous functions $f : [0, T] \to \mathbb{R}^n$.

**Definition 1.** Let $y \in \mathbf{C}([0, T]; \mathbb{R}^n)$, $y_0 \in \bar{D}$. A pair $(x, k) \in \mathbf{C}([0, T]; \mathbb{R}^{2n})$ is said to be a solution to the Skorokhod problem associated with $y$ and $D$ if

1. $x_t = y_t + k_t$, $t \in [0, T]$,

2. $x_t \in \bar{D}$, $t \in [0, T]$,

3. $k$ is a function with bounded variation $|k|$ on $[0, T]$, $k_0 = 0$ and

$$k_t = \int_0^t n_s d|k|_s, \quad |k|_t = \int_0^t 1_{\{x_s \in \partial D\}} d|k|_s, \quad t \in [0, T],$$

where $n_s = n(x_s)$ is an inward normal unit vector at $x_s \in \partial D$.

It is known that if $D$ is a convex set, then there exists a unique solution to the Skorokhod problem [25].

**Definition 2.** Let $(\Omega, \mathcal{F}, \mathcal{P})$ be a probability space.

1. An $n$-dimensional stochastic process $X = \{X_t; t \in [0, T]\}$ is a parametrised collection of random variables defined on a probability space $(\Omega, \mathcal{F}, \mathcal{P})$ with values in $\mathbb{R}^n$.

For each fixed $\omega \in \Omega$ the function $X_t(\omega)$, $t \in [0, T]$ is called a trajectory of $X$ and is denoted by $X(\omega)$.

2. A filtration $(\mathcal{F}_t) = \{\mathcal{F}_t; t \in [0, T]\}$ is a nondecreasing family of sub-$\sigma$-fields of $\mathcal{F}$, i.e. $\mathcal{F}_s \subseteq \mathcal{F}_t \subseteq \mathcal{F}$ for $0 \le s < t \le T$.

By $(\mathcal{F}_t^X)$ we denote a filtration generated by a process $X$, i.e.

$$\mathcal{F}_t^X = \sigma(X_s; 0 \le s \le t).$$

3. A stochastic process $X$ is adapted to the filtration $(\mathcal{F}_t)$ ($X$ is $(\mathcal{F}_t)$ adapted) if for each $t \in [0, T]$, $X_t$ is $\mathcal{F}_t$ measurable random variable.

**Definition 3.** Let $Y$ be $(\mathcal{F}_t)$ adapted process with continuous trajectories, $Y_0 \in \bar{D}$. We say that a pair $(X, K)$ of $(\mathcal{F}_t)$ adapted processes is a solution to the Skorokhod problem associated with $Y$ and $D$, if for almost every $\omega \in \Omega$, $(X(\omega), K(\omega))$ is a solution to the Skorokhod problem associated with $Y(\omega)$ and $D$.
In what follows, by $W = \{W_t; t \in [0, T]\}$ we shall denote a Wiener process starting from zero. We assume that we are given a point $x_0 \in \bar{D}$ and some function $\sigma : \mathbf{R}^n \to \mathbf{R}^n \times \mathbf{R}^m$.

**Definition 4.** Let $Y$ be an $(\mathcal{F}_t)$ adapted process. A pair $(X, K^{\bar{D}})$ of $(\mathcal{F}_t)$ adapted processes is called a solution to reflected stochastic differential equation (in short reflected SDE)

$$X_t = x_0 + \int_0^t \sigma(X_s) dW_s + K_t^{\bar{D}}, \quad t \in [0, T], \qquad (1)$$

if $(X, K^{\bar{D}})$ is a solution to the Skorokhod problem associated with

$$Y_t = x_0 + \int_0^t \sigma(X_s) dW_s, \quad t \in [0, T] \text{ and } D.$$

The process $X$ is called the process with reflection. The proof of existence and uniqueness of the solution to reflected SDEs can be found in [25].

## 3 Stochastic Anisotropic Diffusion

Following [21, 30] we propose the following stochastic model of an anisotropic diffusion:

$$X_t = x + \int_0^t \left[ \begin{array}{c} -\frac{(G_\gamma * u_0)_{x_2}(X_s)}{|\nabla(G_\gamma * u_0)(X_s)|}, 0 \\ \frac{(G_\gamma * u_0)_{x_1}(X_s)}{|\nabla(G_\gamma * u_0)(X_s)|}, 0 \end{array} \right] dW_s + K_t^{\bar{D}},$$

where $u_{x_i}(y) = \frac{\partial u}{\partial x_i}(y)$. To avoid false detections due to noise, $u_0$ is convolved

with a Gaussian kernel $G_\gamma(x) = \frac{1}{2\pi\gamma^2}e^{-\frac{|x|^2}{2\gamma^2}}$ (in practice a 3 × 3 Gaussian mask).

The reconstruction pixel is given by

$$u(x) = \mathbf{E}[u_0(X_T)] \approx \frac{1}{M} \sum_{i=1}^{M} u_0(X_T^m(\omega_i)), \qquad (2)$$

where $X^m(\omega_i)$ is the approximation of trajectory of stochastic process $X$ and $M$ is the number of Monte Carlo method iterations.

### 3.1 Euler's Approximation

Consider the following numerical scheme

$$X_0^m = X_0, \qquad (3)$$

$$X_{t_k}^m = \Pi_{\bar{D}}[X_{t_{k-1}}^m + \sigma(X_{t_{k-1}}^m)(W_{t_k} - W_{t_{k-1}})], k = 1, \dots, m,$$

where $t_k = kh, h = \frac{T}{m}, k = 0, 1, \dots, m$ and $\Pi_{\bar{D}}(x)$ denotes a projection of $x$ on the set $\bar{D}$. Since $D$ is convex, the projection is unique.
**Theorem 1.** Let $(X, K^{\overline{D}})$ be the solution to the reflected SDE (1). If there exists $C > 0$ such that

$$ \|\sigma(x) - \sigma(y)\|^2 \leq C|x - y|^2, $$

then

$$ \lim_{m \to +\infty} |X_T^m - X_T| = 0 \quad \text{almost surely.} $$

The proof of the above theorem can be found in [24].

## 3.2 Modified Diffusion

The numerical scheme (3) gives good results, but only with a small value of the time-step parameter $h = \frac{T}{m}$ (for example $h = 0.05$). Calculating the mean value using Monte Carlo method for small $h$ is not effective and takes a long time. To omit this problem, we improve the scheme (3) by adding a controlled parameter $p$ [4].

$$ \begin{gathered} X_0^m = X_0, \\ H_{t_k}^m = \Pi_{\overline{D}}[X_{t_{k-1}}^m + \sigma(X_{t_{k-1}}^m)(W_{t_k} - W_{t_{k-1}})], \end{gathered} \quad (4) $$

$$ X_{t_k}^m = \begin{cases} H_{t_k}^m, & \text{if } \Theta, \\ X_{t_{k-1}}^m, & \text{elsewhere,} \end{cases} \quad k = 1, 2, \dots, m, $$

where by $\Theta$ we mean the condition

$$ |(G_{\gamma} * u_0)(H_{t_k}^m) - (G_{\gamma} * u_0)(X_{t_{k-1}}^m)| < p. $$

Note that the parameter $p > 0$ guarantees that if the image exhibits a strong gradient then the process $X^m$ diffuses as a process with small value of the parameter $h$ and at locations where variations of the brightness are small, the process $X^m$ can diffuse with a large value of $h$ (for example $h = 4$).

The figure Fig. 1. illustrates a difference between the scheme (3) and the scheme (4). There are shown three examples of trajectories of the process $u_0(X_t^m)$ from the pixel A to the pixel B. Trajectories (I) and (III) were generated using the scheme (3) for large and small value of the parameter $h$, respectively. Trajectory (II) was generated using the scheme (4) for large $h$. It is easy to see, that at locations where the image is constant, trajectory (II) diffuses as trajectory (I). At locations where the image has strong gradient, the trajectory (II) is similar to the trajectory (III).

For small $h$ or $p = +\infty$ (in practice $p > 255$) the numerical scheme (4) is equivalent to the scheme (3).

## 3.3 Modified Diffusion with Random Terminal Time

At locations where gradient is large in all directions it is possible that condition $\Theta$ does not hold as many times as we would expect. To avoid this we propose
<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of trajectories of the process $(G_γ * u_0)(X_t^m)$ from pixel A to B: (I) – with using the scheme (3) and large $h$, (II) – with using the scheme (4) and large $h$, (III) – with using the scheme (3) and small $h$.</figcaption></figure>

the following modification [6]:

$$ 
\begin{aligned}
X_0^m &= X_0, \\
H_{t_k}^m &= \Pi_D[X_{t_{k-1}}^m + \sigma(X_{t_{k-1}}^m)(W_{t_k} - W_{t_{k-1}})], \\
X_{t_k}^m &= \begin{cases} H_{t_k}^m, & \text{if } \Theta, \\ X_{t_{k-1}}^m, & \text{elsewhere,} \end{cases} \quad k = 1, 2, \dots, \tau_m,
\end{aligned}
$$

where $\tau_m = \min\{k; k \ge m \text{ and } \Theta \text{ is true } m \text{ times}\}$.

Terminal time $\tau_m$ guarantees that the numerical simulation of the diffusion trajectory gives at least $m$ values of $X_{t_k}^m$ which differ from the value in the previous step.

## 3.4 Implementation

The above schemes are simple to implement. Observe that the schemes works well only if the model of the digital image $G_γ * u_0$ is continuous. In practice, we can use a linear interpolation to get the value of the image $G_γ * u_0$, for any point $x \in \bar{D}$. We note also that since $W_{t_k} - W_{t_{k-1}} \sim \mathcal{N}(0, t_k - t_{k-1})$, it can be approximated by a random number generator of the normal distribution.

# 4 Non Local Means Algorithm

In this section we cite results from [8–10].

Let $v = \{v(i)|i \in I\}$ be a discrete noisy image and $\{w(i,j)\}$ be the weights that depend on the similarity between the pixels $i$ and $j$ and satisfy the usual conditions $0 \le w(i,j) \le 1$ and $\sum_j w(i,j) = 1$. The reconstructed value $NL(v)(i)$
for a pixel i is defined as a weighted average of all pixels in the image

$$NL(v)(i) = \sum_{j \in I} w(i, j)v(j).$$

The weight $w(i, j)$ depends on the similarity of the intensity gray level vectors of neighbourhoods centred at pixels i and j.

**Definition 5.** A neighbourhood system on I is a family $\mathcal{N} = \{\mathcal{N}_i\}_{i \in I}$ of subsets of I such that for all $i \in I$,

1. $i \in \mathcal{N}_i$,

2. $j \in \mathcal{N}_i \Rightarrow i \in \mathcal{N}_j$.

The subset $\mathcal{N}_i$ is called the neighbourhood or the similarity window of i.

The weights can be defined by

$$w(i, j) = \frac{1}{Z(i)} e^{-\frac{d(\mathcal{N}_i, \mathcal{N}_j)}{s^2}}$$

where $Z(i)$ is the normalising factor $Z(i) = \sum_j e^{-\frac{d(\mathcal{N}_i, \mathcal{N}_j)}{s^2}}$ and $d(\mathcal{N}_i, \mathcal{N}_j)$ is some measure of distance between intensity gray level vectors of similarity windows. The number $s$ is a parameter that controls the decay of the exponential function.

In original approach [8] the authors propose to use square windows of fixed size as similarity windows (see Fig. 2) and the distance between neighbourhoods was measured as a decreasing function of the weighted Euclidean distance i.e. $d(\mathcal{N}_i, \mathcal{N}_j) = \|v(\mathcal{N}_i) - v(\mathcal{N}_j)\|_{2,a}^2$, where $v(\mathcal{N}_i) = (v(j), j \in \mathcal{N}_i)$ and $a > 0$ is the standard deviation of the Gaussian kernel.

## 4.1 Patchwise Implementation

By $\|B_{i,r} - B_{j,r}\|_2$ we denote the Euclidean distance between $B_{i,r}$ and $B_{j,r}$, where patch $B_{i,r}$ means a neighbourhood of a size $2r + 1 \times 2r + 1$ pixels centred at i. Patchwise implementation [10] is based on a simple observation. When computing the Euclidean distance $\|B_{i,r} - B_{j,r}\|_2$, all pixels in the patch $B_{i,r}$ have the same importance, and therefore the weight $F(\|B_{i,r} - B_{j,r}\|_2)$, where $F$ is a decreasing function, can be used to denoise all pixels in the patch $B_{i,r}$ and not only i. For computational purposes the searching of similar windows can be restricted from all pixels in the image to some square window $B_{i,f}$. The denoising of an image $v$ and a certain patch $B_{i,r}$ is equal to

$$\hat{B}_{i,r} = \frac{1}{Z} \sum_{j \in B_{i,f}} v(B_{j,r})w(B_{i,r}, B_{j,r}),$$

where $Z = \sum_{j \in B_{i,f}} w(B_{i,r}, B_{j,r})$ and the weight function is given by

$$w(B, Q) = e^{-\frac{\max(\|B-Q\|_2^2 - 2\rho^2, 0, 0)}{s^2}}.$$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Idea of NL-means from [8]: Pixels *j* and *k* have large weights *w*(*i*, *j*) and *w*(*i*, *k*) because their similarity windows are similar to that of *i*. The weight *w*(*i*, *l*) is much smaller because the intensity grey values in the similarity windows are very different.</figcaption></figure>

Here by $\rho$ we denoted the standard deviation of the noise and by $s$ the filtering parameter set depending on the value of $\rho$. The weight function is chosen in order to average similar patches up to noise. That is, patches with square distances smaller than $2\rho^2$ are set to 1, while larger distances decrease rapidly accordingly to the exponential kernel.

By applying the procedure for all patches in the image, we will get $(2r+1)^2$ possible estimates for each pixel. These estimates can be finally averaged at each pixel location in order to build the final denoised image

$$NL(v)(i) = \frac{1}{(2r+1)^2} \sum_{j \in B_{i,r}} \hat{B}_{j,r}(i).$$

# 5 Anisotropic Stochastic Diffusion Collaborated with Non Local Means

In this section we propose a new method of the image reconstruction based on modified diffusion with random terminal time and patchwise implementation of non local means.

In the case of numerical scheme with random terminal time the reconstructed formula of anisotropic diffusion (2) can be written as

$$u(x) \approx \sum_{i=1}^{M} \frac{1}{M} u_0(X_{\tau_m}^m(\omega_i)),$$

which means that each pixel $u_0(X_{\tau_m}^m(\omega_i))$ is weighted with the value $\frac{1}{M}$. But since pixels have different intensities we may consider them with different weights
depending on their neighbourhood. We follow NL-means algorithm and propose a new method of the image restoration based on modified diffusion but such that the weights depend on patches similarity:

$$u(x) = \frac{1}{Z} \sum_{i=1}^{M} u_0(X_{\tau_m}^m(\omega_i)) w(B_{x,r}, B_{X_{\tau_m}^m(\omega_i),r}).$$

In the above formula we used the following notations:

$$Z = \sum_{i=1}^{M} w(B_{x,r}, B_{X_{\tau_m}^m(\omega_i),r}), w(B, Q) = e^{-\frac{\max(\|B-Q\|_2^2 - 2\rho^2, 0.0)}{\rho^2}}$$

$$X_0^m(\omega_i) = x,$$

$$H_{t_k}^m(\omega_i) = \Pi_{\bar{D}}[X_{t_{k-1}}^m(\omega_i) + \sigma(X_{t_{k-1}}^m(\omega_i))(W_{t_k} - W_{t_{k-1}})],$$

$$X_{t_k}^m(\omega_i) = \begin{cases} H_{t_k}^m(\omega_i), & \text{if } \Theta, \\ X_{t_{k-1}}^m(\omega_i), & \text{elsewhere,} \end{cases} \quad k = 1, 2, \dots, \tau_m.$$

The meaning of the parameters in the new method is the same as in original approaches. This method of the image reconstruction we will call the stochastic diffusion with non local means (in short SDNLM).

# 6 Experimental Results

Some measures of quality for our evaluation experiments regarding new method, non local means algorithm, anisotropic Perona-Malik model [21] and anisotropic stochastic diffusion are presented in Table 1, Table 2, Fig. 3 and Fig. 4. The results refer to greyscale images *pirate* and *cameraman* corrupted with the Gaussian noise with standard deviation $\rho$. The maximum values of Peak Signal to Noise Ratio (in short PSNR) and Structural Similarity Index (in short SSIM) obtained using tested methods are given in tables. Parameters of SSIM were set to the default values as recommended by [28].

The analysis of the measures of image quality shows that in most cases the new method performs better. Moreover, when comparing the figures one can observe that the image created by the SDNLM is visually more pleasant. The reason for this is that the NL-means approach shows clear evidence of a halo of noise effect around the edges whereas anisotropic diffusions smooth details too much.

# 7 Conclusion

In this paper we proposed a new method of digital image denoising. Expressing anisotropic diffusion in terms of stochastic equations allows us to adapt the idea
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** a) Original image: 512 × 512 (top) and 128 × 128 (bottom) b) Noisy image: $\rho = 15$ c) NL-means: SSIM = 0.9245 d) Stochastic anisotropic diffusion: SSIM = 0.9246 e) Perona-Malik: SSIM = 0.9204 f) New method: SSIM = 0.9257</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** a) Original image: 512 × 512 (top) and 96 × 96 (bottom) b) Noisy image: $ρ = 20$ c) NL-means: SSIM = 0.9214 d) Stochastic anisotropic diffusion: SSIM = 0.9144 e) Perona-Malik: SSIM = 0.8838 f) New method: SSIM = 0.9196</figcaption></figure>
<table><caption>**Table 1.** Maximum values of PSNR</caption><thead><tr><th>Image</th><th>Noise ρ</th><th>NL-means algorithm [10]</th><th>Stoch. anisotropic diffusion [6]</th><th>Perona-Malik [21]</th><th>New method</th></tr></thead><tbody><tr><td rowspan="3">Pirate</td><td>10</td><td>33.0394</td><td>32.0450</td><td>32.8268</td><td><b>33.1379</b></td></tr><tr><td>15</td><td>30.7474</td><td>30.5158</td><td>30.8167</td><td><b>30.8944</b></td></tr><tr><td>20</td><td>29.6731</td><td>29.3619</td><td>29.5095</td><td><b>29.6914</b></td></tr><tr><td rowspan="3">Cameraman</td><td>10</td><td>35.3814</td><td>35.2304</td><td>34.9002</td><td><b>35.8238</b></td></tr><tr><td>15</td><td>32.9025</td><td>33.7815</td><td>32.8394</td><td><b>33.8727</b></td></tr><tr><td>20</td><td>32.1931</td><td><b>32.2683</b></td><td>31.2357</td><td>32.1571</td></tr></tbody></table>

<table><caption>**Table 2.** Maximum values of SSIM</caption><thead><tr><th>Image</th><th>Noise ρ</th><th>NL-means algorithm [10]</th><th>Stoch. anisotropic diffusion [6]</th><th>Perona-Malik [21]</th><th>New method</th></tr></thead><tbody><tr><td rowspan="3">Pirate</td><td>10</td><td><b>0.9562</b></td><td>0.9533</td><td>0.9519</td><td>0.9558</td></tr><tr><td>15</td><td>0.9245</td><td>0.9246</td><td>0.9204</td><td><b>0.9257</b></td></tr><tr><td>20</td><td>0.8952</td><td>0.8935</td><td>0.8913</td><td><b>0.8974</b></td></tr><tr><td rowspan="3">Cameraman</td><td>10</td><td>0.9600</td><td>0.9583</td><td>0.9436</td><td><b>0.9606</b></td></tr><tr><td>15</td><td>0.9334</td><td>0.9369</td><td>0.9173</td><td><b>0.9393</b></td></tr><tr><td>20</td><td><b>0.9214</b></td><td>0.9144</td><td>0.8838</td><td>0.9196</td></tr></tbody></table>

from non local means approach. The new method takes what is the best both from anisotropic diffusion and non local means method: reconstructed image is smooth and at the same time details are preserved.

As a future work, the algorithm can be extended to vector valued images, in particular, colour images.

## References

1. Aharon, M., Elad, M., Bruckstein, A.: K-SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation. IEEE Trans. Image Process. 54(11), 4311–4322 (2006)

2. Aubert, G., Kornprobst, P.: Mathematical problems in image processing. Springer-Verlag, New York (2002)

3. Borkowski, D.: Chromaticity Denoising using Solution to the Skorokhod Problem. Image Processing Based on Partial Differential Equations. Mathematics and Visualization, Part II, 149–61 (2007)

4. Borkowski, D.: Modified diffusion to Image Denoising. Adv. Soft Comp. 45, 92–99 (2007)

5. Borkowski, D.: Smoothing, Enhancing Filters in Terms of Backward Stochastic Differential Equations. LNCS 6374, 233–240 (2010)

6. Borkowski, D.: Euler’s Approximations to Image Reconstruction. LNCS 7594, 30–37 (2012)

7. Borkowski, D., Jańczak-Borkowska, K.: Application of Backward Stochastic Differential Equations to Reconstruction of Vector-valued Images. LNCS 7594, 38–47 (2012)

8. Buades, A., Coll, B., Morel, J.M.: A non local algorithm for image denoising. IEEE Computer Vision and Pattern Recognition 2, 60–65 (2005)
9. Buades, A., Coll, B., Morel, J.M.: A review of image denoising algorithms, with a new one. Multiscale Model. Simul. 4(2), 490–530 (2006)

10. Buades, A., Coll, B., Morel, J.M.: Non-Local Means Denoising. Image Processing On Line (2011)

11. Catte, F., Lions, P.L., Morel, J.M., Coll, T.: Image selective smoothing and edge detection by nonlinear diffusion. SIAM J. Numer. Anal. 29(1), 182–193 (1992)

12. Chan, T. F., Shen, J. J.: Image Processing and Analysis – Variational, PDE, wavelet, and stochastic methods. SIAM, Philadelphia (2005)

13. Dabov, K., Foi, A., Katkovnik, V., Egiazarian, K.: Image denoising by sparse 3D transform-domain collaborative filtering. IEEE Trans. Image Process. 16(8), 2080–2095 (2007)

14. Danielyan, A., Katkovnik, V., Egiazarian, K.: Bm3d frames and variational image deblurring. IEEE Trans. Image Process. 21(4), 1715–1728 (2012)

15. Donoho, D.L., Johnstone, I.M.: Ideal spatial adaptation via wavelet shrinkage. Biometrika 81(3), 425–455 (1994)

16. Geman, S., Geman, D.: Stochastic relaxation, gibbs distributions and the bayesian restoration of images. IEEE Pat. Anal. Mach. Intell. 6, 721–741 (1984)

17. Juan, O., Keriven, R., Postelnicu, G.: Stochastic Motion and the Level Set Method in Computer Vision: Stochastic Active Contours. Int. J. Comput. Vision 69(1), 7–25 (2006)

18. Katkovnik, V., Danielyan, A., Egiazarian, K.: Decoupled inverse and denoising for image deblurring: variational BM3D-frame technique. In Proceedings of IEEE International Conference on Image Processing (2011)

19. Lebrun, M., Buades, A., Morel, J.M.: Implementation of the Non-local Bayes image denoising. Image Processing on Line (2011)

20. Mairal, J., Elad, M., Sapiro, G.: Sparse representation for color image restoration. IEEE Trans. Image Process. 17(1), 53–69 (2008)

21. Perona, P., Malik, J.: Scale-space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Anal. Mach. Intell. 12(7), 629–639 (1990)

22. Richardson, W.H.: Bayesian-based iterative method of image restoration. JOSA 62(1), 55–59 (1972)

23. Rudin, L. I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. Phys. D 60, 259–268 (1992)

24. Slomiński, L.: Euler's approximations of solutions of SDEs with reflecting boundary. Stoch. Proc. Appl. 94, 317–337 (2001)

25. Tanaka, H.: Stochastic differential equations with reflecting boundary condition in convex regions. Hiroshima Math. J. 9(1), 163–177 (1979)

26. Unal, G., Krim, H., Yezzi, A.: Stochastic differential equations and geometric flows. IEEE Trans. Image Process. 11(12), 1405–1416 (2002)

27. Unal, G., Ben-Arous, G., Nain, D., Shimkin, N., Tannenbaum, A., Zeitouni, O.: Algorithms for stochastic approximations of curvature flows. Image Processing, Proceedings ICIP 2003 2–3, 651–4 (2003)

28. Wang, Z., Bovik, A. C., Sheikh, H. R., Simoncelli, E. P.: Image quality assessment: From error visibility to structural similarity. IEEE Trans. Image Process. 13(4), 600–612 (2004)

29. Weickert, J.: Theoretical Foundations Of Anisotropic Diffusion In Image Processing. Computing Supplement 11, 221–236 (1996)

30. Weickert, J.: Coherence-Enhancing Diffusion Filtering. Int. J. Comput. Vision 31(2/3), 111–127 (1999)

31. Yaroslavsky, L.P.: Local adaptive image restoration and enhancement with the use of DFT and DCT in a running window. In Proceedings of SPIE 2825, 2–13 (1996)