Evolutionary Integral Models for Image Restoration

(1)

$\begin{aligned} &\frac{\partial u}{\partial n}(t,{\mathbf{x}}) = 0, (t,{\mathbf{x}})\in [0,T]\times\partial\Omega,\end{aligned}$

(2)

where $\nabla u$ , and $\nabla^{2} u$ are, respectively, the gradient and Hessian matrices of u with respect to the space variable ${\mathbf{x}}=(x,y)$ ; $\Omega\subset \mathbb{R}^2$ is a bounded domain (typically a square) with boundary $\partial\Omega$ ; (2) is an homogeneous Neumann boundary condition; u ₀ stands for the image to be restored; and F is a second-order differential operator, [1, 2,27]. The possibilities for F should satisfy basic properties concerning, at least, three aspects of the model: the well-posedness of the continuous problem and discretizations (a way of controlling the stability of the process); the satisfaction of scale-space properties (as a way to have architectural properties of the multiscale analysis, to ensure that the evolved image is a regularized version of the original one or the preservation of important features of the image); finally, the control of this smoothing and also the edge-enhancing in the multiscale process.

In this sense, although the classical Gaussian filtering, with $F(\mathbf{x},u,\nabla u,\nabla^{2} u)={\rm div}\left(\nabla u\right)$ in (1), is well-posed, provides stable discretizations and satisfies several scale-space properties, sometimes it is not efficient in the control of the diffusion, mainly because of the oversmoothing effect. In order to overcome this and other drawbacks, the literature stresses two main ideas: a nonlinear control of the diffusion, and the inclusion of anisotropy to make this control local and capable to discriminate discontinuities and edges. Several proposals in this sense can be seen in, e.g. [1–3, 14, 22, 32, 33], and references therein.

More recent is the use of evolutionary integral equations of the form, [25]

$\begin{aligned} u(t,{\mathbf{x}}) =& u_0(\mathbf{x}) + \int_0^t k(t-s)L u(s,\mathbf{x}){\rm d} s, (t,\mathbf{x})\in[0,T]\times\Omega,\end{aligned}$

(3)

$\begin{aligned} \frac{\partial u}{\partial n}(t,{\mathbf{x}}) =&0, (t,{\mathbf{x}})\in \partial\Omega\times[0,T], \nonumber\end{aligned}$

as models for the multiscale analysis. In (3), $L=\Delta$ stands for the Laplace operator, and k(t) is a convolution kernel. The case

leads to the heat equation, and

to the wave equation with zero initial velocity. (A more general context can be seen e.g. in [9, 18, 23]). If k(t) is differentiable, and

, then (3) is equivalent to the integro-differential problem

$\begin{aligned} u_{t}(t,{\mathbf{x}}) =& \int_0^t k^{\prime}(t-s)L u(s,{\mathbf{x}}){\rm d} s, (t,{\mathbf{x}})\in[0,T]\times\Omega,\end{aligned}$

(4)

$\begin{aligned} u(0,{\mathbf{x}}) =& u_0({\mathbf{x}}),\quad{\mathbf{x}}\in\Omega,\nonumber\\ \frac{\partial u}{\partial n}(t,{\mathbf{x}}) =&0, (t,{\mathbf{x}})\in [0,T]\times\partial\Omega.\nonumber\end{aligned}$

In [8, 11] a control of the diffusion based on (4) with

$\begin{aligned} k(t)=k_\alpha(t)=t^{\alpha-1}/\Gamma(\alpha),\end{aligned}$

(5)

has been proposed, where $\alpha\in (1,2)$ , and Γ is the Gamma function. Model (4), (5) interpolates the linear heat equation ( $\alpha=1$ ), and the linear wave equation ( $\alpha=2$ ), leading α to take a role of viscosity parameter, a term to control the diffusion of the image through the scales [23]. It is also natural to try to handle the diffusion through a selection of α as function of the image at the scale. In [11] a numerical technique, consisting of discretizing (4), (5) with a possibly different value of α for each pixel of the image is introduced. This procedure is modified in [10] to allow to consider a nonconstant viscosity parameter. This approach forms part of the growing interest in the use of fractional calculus for signal processing problems, see [21] for a review of fractional linear systems and also [12,31], along with references therein.

The purpose of this paper is going more deeply into the evolutionary integral modelling for image restoration, generalizing [8, 11] in several ways, representing the following novelties:

Under several non-restrictive hypotheses on the kernel k, the continuous model (3) is proved to satisfy scale-space properties (grey-level shift invariance, reverse-contrast invariance, translational invariance, and conservation of average value). Furthermore, the solution is shown to behave as the constant average value for long times. (Although the application of the evolutionary model (1) to image restoration does not usually require long times of computation, a good behaviour in this sense should always be taken into account).
The semi-discrete (in space) version of (3) is also studied. Under some hypotheses on the discrete spatial operator, it is proved that the corresponding semi-discrete model also possesses several scale-space properties (grey-level shift invariance, reverse-contrast invariance, and conservation of a semi-discrete average value) as well as the constant behaviour as limit for long times. When the semi-discrete model is considered as an approximation to the continuous one (3), these properties enforce the relation between them.
From the computational point of view, the freedom to choose the kernel k is strongly emphasized, since it can be used to control several features of the image: restoration, noise removal, or edge detection. Such properties are illustrated here by means of some examples with medical images.

According to these new results, the structure of the paper is as follows. Section 2 is devoted to the analysis of the above mentioned properties of the continuous model (3). These properties are proved for the Laplace operator, although the way how to extend them to more general spatial operators, [13], is described. The study of the semi-discrete (in space) version is carried out in Sect. 3. Finally, Sect. 4 illustrates the performance of the model with numerical examples. Some details about the implementation are explained and the corresponding codes are applied to several images by using different choices of the kernel. Sect. 5 contains some conclusions and future lines of research.

2 Continuous Evolutionary Integral Models

With the purpose of investigating the degree of adaptation of the evolutionary integral approach to the image processing rules, derived here are some properties of the continuous model (3). The following hypotheses on the kernel function k are assumed:

(H1) if $t\le0$ , and 0$ ” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_IEq17.gif”>.
(H2) k(t) is piecewise differentiable, of subexponential growth.
(H3) The integral

$\int_{0}^{+\infty} k(t){\rm d} t,$
is divergent but k(t) is locally integrable on $(0,+\infty)$ .
(H4) k(t) is 2-regular, [25]; this means that there is a constant denotes the Laplace transform of k(t), then

$\begin{aligned} \left|z^{n}\frac{d^{n}}{dz^{n}}\mathcal{L}\left(k\right)(z)\right|\leq c\left|\mathcal{L}\left(k\right)(z)\right|,\end{aligned}$
for all z with $${\rm Re}(z)>0$ ” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_IEq21.gif”></SPAN>, and <SPAN id=IEq22 class=InlineEquation><IMG alt=$ .

2.1 Well-Posedness

A first point of analysis concerns the well-posedness of the problem, which is usually a nontrivial question for some nonlinear models in image processing, [1,22,33]. In this case, under the hypotheses (H1)–(H4), results about existence, uniqueness, and regularity of solutions are obtained directly from the general theory of Volterra equations. Let S(t) be the resolvent of (3), that is, the transitional operator such that

$\begin{aligned} u(t,\mathbf{x})=S(t)\, u_{0}(\mathbf{x}), t\geq 0, \mathbf{x}\in \Omega,\end{aligned}$

(6)

is the solution of (3) at $\mathbf{x}$ and time t, with initial condition u ₀. It can be proved, see e.g. [25], Theorem 3.1, that u in (6) is ${\cal C}^{1}$ , and there is $M\geq 1$ such that

$$$\begin{aligned} &||u(t,\mathbf{x})||\leq M||u_{0}(\mathbf{x})||,\\ &||u_{t}(t,\mathbf{x})||\leq \frac{M}{t}||u_{0}(\mathbf{x})||, t>0.\end{aligned}$$ ” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_Eque.gif”></DIV></DIV></DIV>These inequalities introduce the diffusion character of (<SPAN class=InternalRef><A href=$ 3), [25].

2.2 Scale-Space Properties

A second group of theoretical properties of the model consists of scale-space properties. They are collected in the following theorem.

Theorem 1

Let S(t) be the transitional operator defined in ( 6 ). Under the hypotheses (H1)–(H4) the following properties hold:

(P1) Grey level shift invariance: $S(t)(u_{0}+C)=S(t)u_{0}+C,$ for any constant C , and if , then .
(P2) Reverse contrast invariance: for $t\geq 0$ , $S(t)(-u_{0})=-S(t)u_{0}$ .
(P3) Translational invariance: if $\tau_{\mathbf{h}}(u_{0})(\mathbf{x})=u_{0}(\mathbf{x}+\mathbf{h})$ , for $\mathbf{x}, \mathbf{x}+\mathbf{h}\in \Omega$ , then

$S(t)(\tau_{\mathbf{h}}u_{0})=\tau_{\mathbf{h}}\left(S(t)u_{0}\right), t\geq0.$
(P4) Conservation of average value: if

where $A(\Omega)$ stands for the area of Ω.

Proof

Hereafter, for the sake of the simplicity of the notation, $u_0(\mathbf{x})$ will be denoted by u ₀. Properties (P1)–(P3) are consequence of the uniqueness of solution. It is clear that $S(t)u_0=0$

. On the other hand, if C is a constant, the functions $u_{1}(t)= S(t)(u_{0}+C)$ and $u_{1}(t)= S(t)(u_{0})+C$ are solutions of (3) with initial condition $u_{0}+C$ ; thus, uniqueness proves (P1). The same argument proves (P2). As far as (P3) is concerned, note that

$\begin{aligned} &{\tau_{\mathbf{h}}u_{0}(\mathbf{x})+\int_{0}^{t}k(t-s)\Delta\tau_{\mathbf{h}} \left(S(t)u_{0}\right){\rm d} s}\\ &=\tau_{\mathbf{h}}u_{0}(\mathbf{x})+\int_{0}^{t}k(t-s)\tau_{\mathbf{h}}\Delta \left(S(t)u_{0}\right){\rm d} s =\tau_{\mathbf{h}}\left(S(t)u_{0}\right).\end{aligned}$

Thus, $\tau_{\mathbf{h}}\left(S(t)u_{0}\right)$ satisfies (3) with initial condition $\tau_{\mathbf{h}}u_{0}$ , and therefore coincides with $S(t)(\tau_{\mathbf{h}}u_{0})$ . Finally, observe that if

$I(t)=\int_{\Omega} u(t,\mathbf{x}){\rm d} x,$

then the regularity of the solution implies that I(t) is continuous, for $t\geq 0$ , differentiable, for $t>0$
” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_IEq45.gif”></SPAN>, and<br />
<DIV id=Equj class=Equation><br />
<DIV class=EquationContent><br />
<DIV class=MediaObject><IMG alt=

$t>0$
” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_IEq45.gif”></SPAN>, and<br />
<DIV id=Equj class=Equation><br />
<DIV class=EquationContent><br />
<DIV class=MediaObject><IMG alt=

Now, the divergence theorem, and the boundary conditions imply that $\frac{\rm d}{{\rm d} t}I(t)=0$ and therefore I is constant, for all $t\geq 0$ .□

2.3 Long Time Behaviour

A final relevant property is the behaviour of the solution as $t\rightarrow +\infty$ . Assuming, for simplicity, that the square domain is $\Omega=(0,1)\times (0,1)$ , and using separation of variables, the solution of (3) can be expressed in the form

$\begin{aligned} u(t,\mathbf{x})=\sum_{l,m} T_{l,m}(t)V_{l,m}(\mathbf{x}),\end{aligned}$

(7)

where $\{V_{l,m}\}_{l,m}$ is an orthogonal basis of eigenfunctions of the eigenvalue problem for the Laplace operator

$\begin{aligned} \Delta V(\mathbf{x})=&\lambda V(\mathbf{x}), \mathbf{x}\in\Omega,\\ \frac{\partial V}{\partial n}(\mathbf{x}) =&0, \quad\mathbf{x}\in \partial\Omega,\end{aligned}$

that has the eigenvalues $\lambda_{l,m}=-(l\pi)^{2}-(m\pi)^{2}$ , with a complete, orthogonal system of eigenfunctions, $V_{l,m}(\mathbf{x})=\cos(l\pi x)\cos(m\pi y)$ , for $l,m\in\mathbb{N}\cup\{0\}$ . In particular $V_{0,0}(\mathbf{x})=1$ , $\lambda_{0,0}=0$ . The expansion of the initial condition is, using the orthogonality, of the form

$\begin{aligned} u_{0}(\mathbf{x})=\sum_{l,m}\gamma_{l,m}V_{l,m}(\mathbf{x}), \gamma_{l,m}=\frac{{\displaystyle}\int_{\Omega}u_{0}(\mathbf{x})V_{l,m}(\mathbf{x}){{\rm d}{\mathbf{x}}}} {{\displaystyle}\int_{\Omega}V_{l,m}^{2}(\mathbf{x}){\rm d} x},\end{aligned}$

where in particular $\gamma_{0,0}=(1/A(\Omega)){{\displaystyle}\int_{\Omega}u_{0}(\mathbf{x}){{\rm d}{\mathbf{x}}}}$ is the average grey value. The time-dependent components of the expansion (7) satisfy the integro-differential problems

$$$\begin{aligned} T_{l,m}^{\prime}(t)=&\lambda_{l,m}\int_{0}^{t}k^{\prime}(t-s)T_{l,m}(s){\rm d} s, t>0, l,m\in\mathbb{N}\cup\{0\},\end{aligned}$$ ” src=”http://basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_15_Chapter_Equ8.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(8)</DIV></DIV><br /> <DIV id=Equm class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt=$

that can be written in terms of the Laplace transform

$\begin{aligned} z\, \mathcal{L}\left(T_{l,m}\right)(z)= \gamma_{l,m}+\lambda_{l,m}\mathcal{L}\left(k^{\prime}\right)(z)\mathcal{L}\left(T_{l,m}\right)(z).\end{aligned}$

This leads to

$\begin{aligned} \mathcal{L}\left(T_{l,m}\right)(z)= \frac{\gamma_{l,m}}{z\left(1-\lambda_{l,m}\mathcal{L}\left(k\right)(z)\right)}.\end{aligned}$

Therefore, for $\lambda_{l,m}<0$ , and under the hypotheses (H1)–(H4), as $t\rightarrow +\infty$ , $T_{l,m}(t)\rightarrow 0$ and the solution u behaves like

$\begin{aligned} T_{0,0}(t)V_{0,0}(\mathbf{x})=\gamma_{0,0}=(1/A(\Omega)){\int_{\Omega}u_{0}(\mathbf{x}){{\rm d}{\mathbf{x}}}}.\end{aligned}$

Remark 1

Extension to more general spatial operators. The previous analysis has been performed by using the Laplacian as spatial operator. The model (3) can be generalized by considering an unbounded, closed, densely defined linear operator L with domain D(L) on some Hilbert space $H(\Omega)$ of functions defined on Ω. Some of the previous properties also hold in this general case. More explicitly, we assume the following:

(h1) L is negative, in the sense that $\langle Lu,u\rangle\leq 0$ , where $\langle\cdot,\cdot\rangle$ stands for the inner product in $H(\Omega)$ , and $u\in D(L)$ .
(h2) L is self-adjoint under Neumann boundary conditions.
(h3) Under Neumann boundary conditions, $\lambda=0$ is a simple eigenvalue of L with Ker generated by the constant functions.

It can be seen that under hypotheses (H1)–(H3), properties (P1), (P2) and (P4) of Theorem 1 also hold in this case. On the other hand, the satisfaction of property (P3) requires the additional hypothesis of commutation between the operator L and $\tau_{\mathbf{h}}$ . Finally, concerning the long-term behaviour, the conclusion above is also satisfied by using the Spectral Theorem (e.g. [5]) for L as follows: there is a spectral family of projections $\{P_{\lambda}\}_{\lambda}$ in such a way that

$u(t,\mathbf{x})=\int_{-\infty}^{0}{\rm d}(P_{\lambda}u)(t,\mathbf{x}).$

Then (8) becomes

$\begin{aligned} (P_{\lambda}u)_{t}(t)=&\lambda\int_{0}^{t}k^{\prime}(t-s)(P_{\lambda}u){\rm d} s,\\ (P_{\lambda}u)(0)=&P_{\lambda}u_{0}.\end{aligned}$

Only gold members can continue reading. Log In or Register to continue

Stay updated, free articles. Join our Telegram channel

Tags: Developments in Medical Image Processing and Computational Vision

Jun 14, 2017 | Posted by admin in GENERAL SURGERY | Comments Off

Basicmedical Key

Fastest Basicmedical Insight Engine

Evolutionary Integral Models for Image Restoration