A NEW NONLOCAL NONLINEAR DIFFUSION EQUATION: THE ONE-DIMENSIONAL CASE

Abstract

We prove a result on the existence and uniqueness of the solution of a new feature-preserving nonlinear nonlocal diffusion equation for signal denoising for the one-dimensional case. The partial differential equation is based on a novel diffusivity coefficient that uses a nonlocal automatically detected parameter related to the local bounded variation and the local oscillating pattern of the noisy input signal.

1 Introduction

Nonlinear partial differential equations (PDEs) can be used in the analysis and processing of digital images or image sequences, for example, to extract features and shapes or to filter out the noise to produce higher quality images (see, for example, [3, 4, 14, 15] and the references therein). Arguably, the main application of PDE-based methods in this field is the smoothing and restoration of images. From the mathematical point of view, the input (grey scale) image can be modelled by a real function $u_{0}(x)$ , $u_{0}:\Omega \rightarrow \mathbb {R}$ , where $\Omega \subset \mathbb {R}^{d}$ represents the spatial domain. Typically, this domain $\Omega $ is rectangular and $d=1, 2$ or $3$ . The function $u_{0}$ is considered as initial data for a suitable evolution equation with some kind of boundary conditions. The simplest (and oldest) PDE method for smoothing images is to apply a linear diffusion process: the starting point is the simple observation that the so-called Gauss function is related to the fundamental solution of the linear diffusion (heat) equation.

The flow produced by the linear diffusion equation spreads the information equally in all directions. Although this property is good for a local noise reduction in the case of additive noise, this filtering operation also destroys the image content such as the boundaries of the objects and the subregions present in the image. This means that the Gaussian smoothing not only smooths noise, but also blurs important features in the signal.

Recently, a new anisotropic diffusion model was introduced in [11] to analyse experimental signals in neuroscience: the diffusivity coefficient uses a nonlocal parameter related to the local bounded variation and the local oscillating pattern of the noisy input signal. In [2], the model was extended to the multidimensional case with an analysis for the existence of the solution in the two-dimensional case (images) and the introduction of a suitable numerical scheme. In this note, we focus on the one-dimensional case providing a complete analysis of the nonlocal diffusion equation, including the uniqueness that was an open problem.

2 A one-dimensional nonlocal nonlinear model

There is a vast literature concerning nonlinear anisotropic diffusions with applications to image processing, which dates back to the seminal paper by Perona and Malik [12], who considered a discrete version of the problem

(2.1) $$ \begin{align} \begin{cases} \dfrac{\partial u}{\partial t} - \nabla\cdot (g(|\nabla u|)\nabla u) = 0 \quad &\textrm{in } \Omega_{T}=(0,T)\times \Omega , \\[5pt] u(x,0)=u_{0}(x) \quad &\textrm{on } \Omega, \\[2pt] \dfrac{\partial u}{\partial \vec{n}} (x,t) =0 \quad &\textrm{on } \Gamma\times (0,T), \end{cases} \end{align} $$

where $ \Gamma = \partial \Omega $ , the image domain $\Omega \subset \mathbb {R}^{2}$ is an open regular set (typically a rectangle), $\vec {n}$ denotes the unit outer normal to its boundary $\Gamma $ , $\nabla \cdot $ is the divergence operator, and $u(x,t)$ denotes the (scalar) image analysed at time (scale) t and point x. The initial condition $u_{0}(x)$ is, as in the linear case, the original image. To reduce smoothing at the edges, the diffusivity g is chosen as a decreasing function of the ‘edge detector’ $|\nabla u|$ . Here, we introduce a nonlocal diffusive coefficient that considers the ‘monotonicity’ of the signal. In other words, a high modulus of the gradient may lead to a small diffusion if the function is also locally monotone. At the same time, we want to reduce the noise present, as in the case of linear diffusion. We focus on the one-dimensional case, more precisely, where $u:[a,b] \rightarrow \mathbb {R}$ is a real function defined on a bounded interval $[a,b]$ , and on a subinterval $[c,d] \subset [a,b ]$ . We define the local variation $LV_{[c,d]}(u)$ of u on the interval $[c,d]$ by

$$ \begin{align*} LV_{[c,d]}(u) = |u(d) - u(c)|. \end{align*} $$

We also define the total local variation $TV_{[c,d]}(u)$ of u on the interval $[c,d]$ by

$$ \begin{align*} TV_{[c,d]}(u) = \sup_{\mathcal{P}} \sum_{i=0}^{n_{P}-1} | u ( x_{i+1}) - u( x_{i}) | \end{align*} $$

where $\mathcal {P} = \{ P = \{ x_{0}, \ldots , x_{n_{P}} \} \mid P \textrm { is a partition of } [ c,d ] \} $ is the set of all possible finite partitions of the interval $[c,d]$ .

Let $\varepsilon \in \mathbb {R}^{+}$ , $\varepsilon \ll 1, \, \varepsilon>0$ and let $\delta \in \mathbb {R}^{+}$ . We define the ratio,

$$ \begin{align*} R_{\delta ,u} = \frac{\varepsilon + LV_{[ x-\delta ,x+\delta ]}(u)}{\varepsilon + TV_{[x-\delta ,x+\delta ]}(u)}. \end{align*} $$

If the parameter $\delta $ is chosen appropriately, we can distinguish between oscillations caused by noise contained in a range of amplitude $\delta $ . As in the Perona–Malik model given by (2.1), we adapt the diffusivity coefficient by using the above ratio $R_{\delta ,u}$ . For small values of the latter, we have to reduce the noise, while for values close to $1$ , the upper bound of $R_{\delta ,u}$ , we have to preserve the signal variation (as the edges in the image). The resulting diffusivity coefficient $g(R_{\delta ,u})$ becomes nonlocal. We assume that $g:[0,+\infty )\rightarrow \mathbb {R}$ is a positive, nonincreasing, Lipschitz continuous function such that $g(0)=1$ and $g(1)=\alpha>0$ . In the following, we assume that the parameter $\varepsilon $ ( $0<\varepsilon \ll 1$ ) is fixed. In Figure 1, we show an illustrative example of a denoised signal using our nonlocal and nonlinear diffusion filter. In particular, we have numerically simulated (2.2) by adopting a semi-implicit method based on central finite differences (see [2]) and with the following numeric values of the parameters (see also (2.3)):

$$ \begin{align*} g(s) = \begin{cases} 1 \quad &\text{ if }s=0\text{,} \\[2pt] 1 - s^{2} e^{-3.315/(s/\lambda)^{8})} \quad &\text{ if }s\neq 0\text{,} \end{cases}\quad \lambda= 4,\ \varepsilon =10^{-3},\ \delta = 0.075, \end{align*} $$

and the time domain $t\in [0,0.6]$ and space domain $x\in [0, 255]$ . The signal in Figure 1 was obtained from a simulation of a biophysical model of a neuron with an additive Gaussian noise (mean equal to $0$ and variance equal to $3$ ) (see [1] for more details). The MATLAB code and the details are available from the authors.

Figure 1 An illustrative example of signal denoising using the new nonlocal and nonlinear diffusion equation (2.2) with data from [1]. Solid line, original signal; grey line, signal with noise; dotted line, reconstructed signal.

In the following, $I=(a,b)\subset \mathbb {R}$ denotes a bounded open interval and $H^{k}(I)$ , $k\in \mathbb {N}$ , the Sobolev space of all functions u defined in I such that u and its distributional derivatives of order $1,\ldots ,k$ all belong to $L^{2}(I)$ . Let $D^{s}$ denote the distributional derivative. Then $H^{k}(I)$ is a Hilbert space for the norm

$$ \begin{align*} \| u \|_{k} = \| u \|_{H^{k}} = \bigg( \sum_{|s|\leq k} \int_{I} | D^{s} u(x) |^{2} \,dx \bigg)^{1/2} , \quad \| u \|_{0} = \| u \|_{L^{2}}. \end{align*} $$

Let $L^{p}(0, T; H^{k}(I))$ be the set of all functions u, such that, for almost every t in $(0, T)$ with $T>0$ , $u(t)$ belongs to $H^{k}(I)$ . Then $L^{p}(0, T; H^{k}(I))$ is a normed space for the norm

$$ \begin{align*} \|u \|_{L^{p}(0, T; H^{k}(I ))} = \bigg( \int_{0}^{T}\| u \|^{p}_{k} \,dt \bigg)^{1/p}, \end{align*} $$

where $p \geq 1$ and $k\in \mathbb {N}$ . Finally, we denote by $(\cdot \, , \cdot )$ the scalar product in $L^{2}(I)$ .

We now establish our existence result. As initial conditions, we take the original signal $u_{0}$ but with some regularisation obtained with a standard smoothing filter, for example, a Gaussian filter, and we assume homogeneous Neumann conditions at the boundary.

Theorem 2.1 (Existence).

Let $u_{0} \in H^{1}(I)$ and $T>0$ , $\delta>0$ . Then there exists $u\in L^{2}(0,T;H^{1}(I))\bigcap C^{0}([0,T]; L^{2}(I))$ , satisfying u(x, 0) = u ₀ (x) on I, ${\partial u}/{\partial x} =0$ at $x=a,\, b$ , and

(2.2) $$ \begin{align} \frac{\partial u}{\partial t} - \frac{\partial}{\partial x} \bigg( g( R_{\delta ,u}) \frac{\partial u}{\partial x} \bigg)=0, \end{align} $$

on $(0,T]\times I$ in the distributional sense.

Proof. We show the existence of a weak solution of (2.2) by a classical fixed point theorem of Schauder (see, for example, [7, Theorem 2.2]). We introduce the space

$$ \begin{align*} V(0,T)=\bigg\{ v\in L^{2}(0,T;H^{1}(I)),\,\frac{dv}{dt}\in L^{2}(0,T;(H^{1}(I))^{\prime} ) \bigg\}. \end{align*} $$

The space $V(0,T)$ is a Hilbert space with the graph norm. Let v be a function in $V(0,T)\bigcap L^{\infty } (0,T; L^{2}(I))$ such that

$$ \begin{align*} \| v \|_{L^{\infty} (0,T; L^{2}(I))} \leq \| u_{0}\|_{L^{2}(I)}. \end{align*} $$

We consider the following variational problem $(P_{v})$ :

$$ \begin{align*} \bigg\langle \frac{\partial u}{\partial t}(t),w \bigg\rangle + \int_{I} g(R_{\delta ,v}) \frac{\partial u(t)}{\partial x}\frac{\partial w}{\partial x} \,dx & =0\quad \mbox{for all } w\in H^{1}(I)\,\,\text{(a.e.) in}\,\,[0,T] \\[2pt] u(0) & \in H^{1}(I). \nonumber \end{align*} $$

Here $\langle \cdot \, , \cdot \rangle $ represents the duality product. A function $u\in H^{1}(I)$ has locally bounded variation (see, for example, [10, Theorem 5.1]) and, moreover, is equal almost everywhere (a.e.) to an absolutely continuous function and $u^{\prime }$ exists a.e. and belongs to $L^{2}(I)$ . The term $R_{\delta ,v}$ can be represented as

(2.3) $$ \begin{align} R_{\delta ,v} = \frac{\varepsilon + \Big|\int_{x-\delta}^{x+\delta} u^{\prime} (s) \,ds \Big|}{\varepsilon + \int_{x-\delta}^{x+\delta} | u^{\prime} (s)| \,ds} \end{align} $$

and $0<R_{\delta ,v}\leq 1$ . So, $g(R_{\delta ,v})\geq \alpha>0$ .

Using classical results about parabolic equations (see, for example, [8, Theorem 10.1] and [9, Theorem 7.3]), the problem $(P_{v})$ has a unique solution $U(v)$ in $V(0,T)$ . We can deduce the following estimates:

(2.4) $$ \begin{align} \| U(v) \|_{L^{\infty} (0,T;L^{2}(I))} & \leq \| u_{0}\|_{L^{2}(I)}, \nonumber\\[2pt] \| U(v) \|_{L^{2} (0,T;H^{1}(I))} & \leq C_{1}, \\[2pt] \| U(v) \|_{L^{2} (0,T;(H^{1}(I))^{\prime} )} & \leq C_{2}, \nonumber \end{align} $$

for suitable constants $C_{1}$ and $C_{2}$ depending only on $u_{0}$ , T and the Lipschitz constant of the function g. We introduce the subset $V_{0}$ of $V(0,T)$ defined by functions $v\in V(0,T)$ such that these estimates are satisfied and $v(0)=u_{0}$ . Then U is a mapping from $V_{0}$ to $V_{0}$ . Moreover, $V_{0}$ is a nonempty, convex and weakly compact subset of $V(0,T)$ .

To use the Schauder theorem, we have to prove that the mapping $v \rightarrow U(v)$ is weakly continuous from $V_{0}$ to $V_{0}$ . Then, since $V(0,T)$ is contained in $L^{2}(0,T; L^{2}(I))$ with compact inclusion, this yields the existence of $u\in V_{0}$ such that $u=U(u)$ .

Let $(v_{j})$ be a sequence in $V_{0}$ which converges weakly to $v\in V_{0}$ and $u_{j}=U(v_{j})$ . From the classical theorems of compact inclusion (see, for example, [8, Theorem 9.16]), up to sub-sequences,

$$ \begin{align*} u_{j} & \rightarrow u\quad\ \ \text{weakly in }\ L^{2}(0,T; H^{1}(I)), \\[2pt] \frac{du_{j}}{dt} & \rightarrow \frac{du}{dt} \quad \text{weakly in}\ L^{2}(0,T; (H^{1}(I))^{\prime} ),\\[2pt] \frac{\partial u_{j}}{\partial x} & \rightarrow \frac{\partial u}{\partial x} \quad \text{weakly in}\ L^{2}(0,T; L^{2}(I) ). \end{align*} $$

Moreover, $u_{j} \rightarrow u$ in $L^{2}(0,T; L^{2}(I))$ and a.e. on $I\times (0,T)$ and $u_{j}(0) \rightarrow u(0)$ in $ (H^{1}(I))^{\prime }$ . For the $(v_{j})$ , from (2.4), there is a subsequence such that $v_{j} \rightarrow v$ in $L^{2}(0,T; L^{2}(I))$ and, from the Rellich–Kodrachov theorem (see, for example, [9, Theorem 5.1], and (2.3)), $g(R_{\delta ,v_{j}}) \rightarrow g(R_{\delta ,v})$ in $L^{2}(0,T; L^{2}(I))$ . By the uniqueness of the solution of $(P_{v})$ , the whole sequence $u_{j}=U(v_{j})$ converges weakly in $V(0,T)$ . Thus, the mapping U is weakly continuous from $V_{0}$ into $V_{0}$ and we can apply the Schauder theorem.

Remark 2.2. A similar proof could be carried through in a more general case by considering a different measure of local variation, for example, using the absolute value of the difference between the maximum and minimum value in subintervals of length $2\delta $ .

Under the hypotheses of Theorem 2.1, we have the following uniqueness result.

Theorem 2.3 (uniqueness).

The solution $u\in L^{2}(0,T;H^{1}(I))\bigcap C^{0}([0,T]; L^{2}(I))$ of (2.2), with $u(0)\in H^{1}(I)$ and homogeneous Neumann conditions, is unique.

Proof. Let $\bar {u}$ and $\hat {u}$ be two solutions of (2.2) and let $u=\bar {u}-\hat {u}$ . Then for almost all t in $[0,T]$ ,

(2.5) $$ \begin{align} \frac{{d}\bar{u}}{dt}-\frac{\partial}{\partial x} \bigg( g ( R_{\delta ,\bar{u}}) \frac{\partial \bar{u}}{\partial x} \bigg) & = 0, \quad \bar{u}(0)=u_{0}, \end{align} $$

(2.6) $$ \begin{align} \frac{{d}\hat{u}}{dt}-\frac{\partial}{\partial x} \bigg( g ( R_{\delta ,\hat{u}}) \frac{\partial \hat{u}}{\partial x} \bigg) & = 0, \quad \hat{u}(0)=u_{0}. \end{align} $$

By subtracting (2.6) from (2.5),

$$ \begin{align*} \frac{d(\bar{u}-\hat{u})}{dt}-\frac{\partial}{\partial x} \bigg( g ( R_{\delta ,\bar{u}}) \frac{\partial \bar{u}}{\partial x}\bigg)+ \frac{\partial}{\partial x} \bigg( g ( R_{\delta ,\hat{u}}) \frac{\partial \hat{u}}{\partial x} \bigg)=0. \end{align*} $$

Adding and subtracting the quantity ${\partial _{x}}( g ( R_{\delta ,\bar {u}}) {\partial _{x}} \hat {u})$ , we can rewrite the equation as

(2.7) $$ \begin{align} \frac{du}{dt}-\frac{\partial}{\partial x} \bigg( g ( R_{\delta ,\bar{u}}) \frac{\partial {u}}{\partial x}\bigg)= \frac{\partial}{\partial x} \bigg( [ g ( R_{\delta ,\bar{u}})- g ( R_{\delta ,\hat{u}})] \frac{\partial \hat{u}}{\partial x} \bigg). \end{align} $$

Multiplying (2.7) by $u=(\bar {u}-\hat {u})$ , integrating on the interval I, using the properties of the function g and the lower bound $g(1)=\alpha>0$ and the estimates (2.4), we obtain

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\|u(t)\|_{L^{2}(I)}^{2} + \alpha \bigg\| \frac{\partial}{\partial x}u(t) \bigg\|_{L^{2}(I)}^{2} \leq C \| u(t) \|_{L^{2}(I)} \bigg\| \frac{\partial}{\partial x}\hat{u}(t) \bigg\|_{L^{2}(I)} \bigg\| \frac{\partial}{\partial x}u(t) \bigg\|_{L^{2}(I)}, \end{align*} $$

for a suitable constant C. The term on the right-hand side can be estimated, using Young’s inequality, by

$$ \begin{align*} \frac{2}{\alpha} C^{2} \| u(t) \|_{L^{2}(I)}^{2} \bigg\| \frac{\partial}{\partial x}\hat{u}(t) \bigg\|_{L^{2}(I)}^{2} + \frac{\alpha}{2} \bigg\| \frac{\partial}{\partial x}u(t) \bigg\|_{L^{2}(I)}^{2}. \end{align*} $$

Subtracting the term $({\alpha }/{2}) \| {(\partial }/{\partial x)}u(t) \|_{L^{2}(I)}^{2}$ on both sides and using the a priori estimates (2.4), we get the inequality

(2.8) $$ \begin{align} \frac{1}{2}\frac{d}{dt}\|u(t)\|_{L^{2}(I)}^{2} +\frac{\alpha}{2} \bigg\| \frac{\partial}{\partial x}u(t) \bigg\|_{L^{2}(I)}^{2} \leq C^{*} \| u(t) \|_{L^{2}(I)}^{2}, \end{align} $$

where $C^{*}=2C^{2}C_{1}/\alpha $ . Since $\bar {u}(0)=\hat {u}(0)=u_{0}$ , by the inequality (2.8) and Gronwall’s lemma (see, for example, [13, Theorem 1.8], we obtain the uniqueness of the solution.

Remark 2.4. Similar nonlocal equations could be obtained as diffusive limits from different kinetic microscale descriptions of the interactions of active particles (see, for example, [5, 6]).