A new image deconvolution method with fractional regularisation

Abstract

Image deconvolution is an important pre-processing step in image analysis which may be combined with denoising, also an important image restoration technique, and prepares the image to facilitate diagnosis in the case of medical images and further processing such as segmentation and registration. Considering the variational approach to this problem, regularisation is a vital component for reconstructing meaningful information and the problem of defining appropriate regularisation is an active research area. An important question in image deconvolution is how to obtain a restored image which has sharp edges where required but also allows smooth regions. Many of the existing regularisation methods allow for one or the other but struggle to obtain good results with both. Consequently, there has been much work in the area of variational image reconstruction in finding regularisation techniques which can provide good quality restoration for images which have both smooth regions and sharp edges. In this paper, we propose a new regularisation technique for image reconstruction in the blind and non-blind deconvolution problems where the precise cause of blur may or may not be known. We present experimental results which demonstrate that this method of regularisation is beneficial for restoring images and blur functions which contain both jumps in intensity and smooth regions.

Introduction

Blur in images is a significant problem in many areas of imaging including medical, astronomical and consumer photography. Typically occurring in the acquisition of medical images due to various factors such as refractive error, the skill of the photographer, media opacity, the age and physical condition of the patient and, in ophthalmic imaging, eye movement, severe blur imposes a significant resolution cost which hampers clinical assessment and prevents further processing. In an existing programme of diabetic retinopathy screening, approximately 5% of the colour fundus images collected are too blurred for assessment. This necessitates the development of effective techniques for obtaining as much information as possible of the underlying sharp and detailed image.

As a mathematical consideration, we view an image as a function over a domain Ω whose range is determined by the bits-per-sample rate of the image. The action of blurring such an image may be then viewed as a convolution $[\kappa *u](\mathbf{x})$ of the image $u(\mathbf{x})$ with a point spread function κ which describes how the image is degraded. We also take into account the presence of noise $\eta (\mathbf{x})$ in the blurred image which makes image deconvolution an ill-posed and challenging problem. We thus model the forward problem of a sharp image $u(\mathbf{x})$ being corrupted by blur $\kappa (\mathbf{x})$ and noise $\eta (\mathbf{x})$ to produce a blurred image $z(\mathbf{x})$ as

$z(\mathbf{x})=[\kappa *u](\mathbf{x})+\eta (\mathbf{x}),[\kappa *u](\mathbf{x})=\int \Omega \kappa (\mathbf{x}-\mathbf{y})u(\mathbf{y})\text{d}\mathbf{y}$

(1)

The idea of approximating the image $u(\mathbf{x})$ given the image $z(\mathbf{x})$ is known as image deconvolution, whose problems may be categorised as: (1) non-blind^1–3 in the case of known blur, so that $\kappa (\mathbf{x})$ is known and $u(\mathbf{x})$ is to be recovered; (2) semi-blind^4–6 in the case where some information of the blur function is known or may be assumed and the task is to complete this information in order to obtain the blur function while restoring the unknown sharp image; (3) blind^7–10 which makes no assumptions about the blur related to a particular problem. In the blind case, the challenge is to identify the blur function and the image using only the given blurred image data.

In this work, we have considered solutions to the issue of blur degradation in images using variational modelling, which is a popular and very successful approach to image processing problems. Such methods typically involve the minimisation of a variational functional consisting of data fitting and regularisation components, among others. In image reconstruction, data fitting terms arise from assumptions about the cause of degradation, such as the forward problem given by equation (1) while regularisation ideas arise from the need to control the ill-posedness of the problem and to improve the quality of the restored image. In forming such terms, assumptions are typically made about the underlying image which is to be restored and is not yet known. The question of how to effectively provide such regularisation has given rise to an active research area.^2,11–14

In this paper, we propose a new method for image deconvolution using fractional derivatives to build a regularisation technique which can outperform comparable regularisation methods. We build in ideas for implicitly constrained image restoration in order to form an effective image deconvolution technique. We next extend this idea to the blind case of restoring the blur function and image simultaneously using fractional regularisation.

The rest of this paper is organised as follows. In the next section, we review some relevant image deconvolution and regularisation techniques. In the section ‘A fractional regularised image deconvolution model’, we present the new method for image deconvolution which is proposed in this paper. In the section ‘A blind deconvolution model with fractional regularisation’, we extend this idea to the case of blind deconvolution. In the ‘Experimental results’ section, we present experimental results and in the final section, we conclude this work.

Existing image deconvolution methods

In this section, we review some existing methods for image reconstruction, focusing on regularisation techniques for data representing underlying smooth and sharp images which we are aiming to recover.

It may be tempting to try to invert the blur function κ of the forward problem (1) in order to recover the image $u(\mathbf{x})$ . After finding the singular value decomposition of the discrete counterpart of the blur function, the ill-posedness of the problem may be mitigated by attempting to filter out small singular values. This can be achieved by multiplying the singular values s_i of the decomposition by a function such as which yield the Truncated and Tikhonov Filter Singular Value Decomposition models, ¹⁵ respectively, where $\lambda \in ℝ>0$ is a small positive parameter. These models are quick to implement and may offer reasonable results but the quality of results can be significantly improved by building variational regularised models which may be optimised in order to provide a solution.

The Tikhonov Filter model may be reformulated as the following variational problem of minimising a functional including a fitting term which makes the assumption that the received data $z(\mathbf{x})$ is the result of a convolution of the blur function $\kappa (\mathbf{x})$ and the hidden true image $u(\mathbf{x})$ which we aim to recover. This assumption is included in the problem by the requirement that this energy be minimised. In order to make the problem well posed, the functional is regularised by a restriction on the square of the L²-norm of the image intensity values

$minu\{\mathcal{F}\mathrm{Tik}=\frac{1}{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+\frac{\lambda }{2}\int \Omega |u(\mathbf{x})|2\text{d}\mathbf{x}\}$

Assuming that the blur function is known, this functional is minimised by solving the below Euler–Lagrange equation resulting from deriving the first-order optimality condition with respect to the image $u(\mathbf{x})$ . We thus recover the image by solving the problem

$\kappa †(\mathbf{x})*([\kappa *u](\mathbf{x})-z(\mathbf{x}))+\lambda u(\mathbf{x})=0$

where $\kappa †(\mathbf{x})=\kappa (-\mathbf{x})$ . Discretising this equation, we may solve it efficiently using the Conjugate Gradient method by rewriting it as

$([\kappa †*\kappa ](\mathbf{x})+\lambda \delta (\mathbf{x}))*u(\mathbf{x})=[\kappa †*z](\mathbf{x})$

Since this equation is linear in terms of the image function $u(\mathbf{x})$ , it is very efficient to solve. However, the rather simple regularisation assumption that the intensity values should be minimised is not reliable and typically insufficient to provide a good quality restored image.

You and Kaveh ¹⁶ made use of the H¹ semi-norm in order to provide regularisation for the image using gradient information. In the non-blind case, this leads to the minimisation problem which is similar to the Tikhonov case except that here we are aiming to minimise the overall gradient of the image. A solution is found by solving the Euler–Lagrange equation

$\kappa †(\mathbf{x})*([\kappa *u](\mathbf{x})-z(\mathbf{x}))-\lambda \Delta u(\mathbf{x})=0$

where Δ denotes the Laplacian of the image. This equation is again linear in terms of $u(\mathbf{x})$ , which can be seen by rewriting it as

$([\kappa †*\kappa ](\mathbf{x})-\lambda D2(\mathbf{x}))*u(\mathbf{x})=[\kappa †*z](\mathbf{x})$

where D² represents the Laplacian as a convolution operation. This model typically improves on the Tikhonov model and can give good results for smooth images but struggles to accurately reconstruct edges.

To improve on H¹, mean curvature gives a popular regularisation technique in image registration and denoising ¹⁴ which may also be applied to image deconvolution. A regularisation term based on this may be given by

$\int \Omega \text{K}(\mathcal{K}(\kappa (\mathbf{x})))\text{d}\mathbf{x},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{K}(s)=\frac{1}{2}s2,\mathcal{K}(\kappa (\mathbf{x}))=\nabla ·\frac{\nabla \kappa (\mathbf{x})}{|\nabla \kappa (\mathbf{x})|\beta }$

(2)

where $|·|\beta =\sqrt{(·)2+\beta }$ for some small $\beta \in ℝ>0$ . The corresponding Euler–Lagrange equation for this term is given by

$\nabla ·(\frac{\nabla \text{K}\text{'}(\mathcal{K}(\kappa ))}{|\nabla \kappa |\beta }-\frac{\nabla \kappa ·\nabla \text{K}\text{'}(\mathcal{K}(\kappa ))}{|\nabla \kappa (\mathbf{x})|\beta 3}\nabla \kappa )=0$

Similar to H¹ regularisation, models using mean curvature for regularisation can struggle to reconstruct sharp edges in images.

In order to improve the quality of restoration at edges, the total variation (tv) semi-norm, popular in image denoising^17–19 as well as deconvolution,^7,11,12 may be used. This gives rise to the optimisation problem

In practice, we replace the absolute value $|·|$ with a smooth approximation $|·|\beta =\sqrt{(·)2+\beta }$ where $\beta \in ℝ>0$ is a small non-negative value used to avoid division by zero at zero-gradient points and which may be tweaked to improve the result. This regularisation term provides improved reconstruction of edges in images but the results are dependent on the value of β selected. Too high a choice of β can result in overly smooth results while too low a choice can result in a piecewise constant result, exhibiting the staircase effect. This model can be used to reconstruct sharp edges but there is typically a trade-off between sharpness and smoothness in the restored image.

Efforts to try to benefit from the edge-preserving properties of total variation regularisation and the smooth results from regularisation techniques such as H¹ within a single image include mixed models. In these cases, the aim is to use tv regularisation in regions with sharp jumps and smooth regularisation terms in more homogeneous regions. In Lysaker and Tai, ²⁰ the authors propose a noise removal method by denoising an image first using tv and then denoising the same image using the second-order LLT model ²¹ in order to obtain approximations $u\mathrm{tv}(\mathbf{x})$ and $u\mathrm{llt}(\mathbf{x})$ , respectively. The final denoised image $u(\mathbf{x})$ is created as a convex combination of the two images, allowing the best regions of each image to appear in the final restoration result. In Chang et al., ²² the authors propose a single model which combines tv and second-order regularisation, presenting the functional minimisation problem for deconvolution where θ is chosen depending on the absolute value of the gradient of intensity values at each point. At lower and larger gradient values, $\theta \to 1$ so that the regularisation is predominantly tv while at smoother regions, as assumed by mid-range gradients, $\theta \to 0$ so that the second-order regulariser has more influence. This is a useful idea but in the case of deconvolution the method relies on correct distinction between sharp and smooth regions in order to determine the regularisation yet such information is not reliable in blurred and noisy images.

A fractional regularised image deconvolution model

In this section, we develop a method for image deconvolution which makes use of fractional regularisation which can preserve both edges and smoothness. We then improve this model by enforcing intensity-range constraints implicitly in the problem functional. The total α-order variation of the function u is given by

$\mathrm{TV}\alpha (u)=sup\phi \in K\int \Omega (-u\text{div}\alpha \phi )\text{d}\mathbf{x},\text{div}\alpha \phi =\sum _{i=1}^d\frac{\partial \alpha \phi i}{\partial x_i^{\alpha }}$

where $\partial \alpha \phi i/\partial x_i^{\alpha }$ denotes a fractional α-order derivative $D_{[a,b]}^{\alpha }\phi i$ of $\phi i$ along the x_i direction. It has been proven that if the function $u(\mathbf{x})$ is contained in the Banach space $W_p^{\alpha }(\Omega )=\{u\in Lp(\Omega )|\left|\right|u\left|\right|W_p^{\alpha }(\Omega )<+\infty \}$ with norm

$\left|\right|u\left|\right|W_p^{\alpha }(\Omega )=(\int \Omega \left|u\right|p\text{d}\mathbf{x}+\int \Omega \left|\nabla \alpha u\right|p\text{d}\mathbf{x})1/p,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\nabla \alpha u=(\frac{\partial \alpha u}{\partial x1},\dots ,\frac{\partial \alpha u}{\partial xd})\top$

then the following equality holds

$\mathrm{TV}\alpha (u)=\int \Omega \left|\nabla \alpha u\right|\text{d}\mathbf{x}$

We build this idea into a deconvolution minimisation problem consisting of data fitting and α-order fractional regularisation as

$minu\{\mathcal{F}=\frac{\lambda }{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+\int \Omega \left|\nabla \alpha u\right|\text{d}\mathbf{x}\}$

(4)

where u is the image to be recovered, z is the available blurred and noisy image data and κ is the blur function which, in this case, we assume to be known. In order to solve this minimisation problem, we may derive the corresponding Euler–Lagrange equation by minimising the functional with respect to the image. First, we consider a method of improving the speed of finding a solution.

We transform equation (4) to the constrained optimisation problem, introducing the variable $\mathbf{d}(\mathbf{x})=(d1(\mathbf{x}),d2(\mathbf{x}))\top$ into the functional to form the new but equivalent problem

$minu,\mathbf{d}\{\mathcal{F}=\frac{\lambda }{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+\int \Omega |\mathbf{d}(\mathbf{x})|\text{d}\mathbf{x}\}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{s}.\text{t}.\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathbf{d}(\mathbf{x})=\nabla \alpha u(\mathbf{x})$

This constraint is enforced by incorporating it into the functional, transforming it to the Bregman formulation

$(uk+1,\mathbf{d}k+1)=\underset{u,\mathbf{d}}{min}\{\frac{\lambda }{2}\int \Omega (\kappa *u-z)2\text{d}\mathbf{x}+\int \Omega \left|\mathbf{d}\right|\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-\int \Omega <{\mathbf{p}}_{\mathbf{d}}^k,\mathbf{d}-\mathbf{d}k>\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-\int \Omega <{\mathbf{p}}_u^k,u-uk>\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\mu }{2}\int \Omega (\mathbf{d}-\nabla \alpha u)2\text{d}\mathbf{x}\},{\mathbf{p}}_u^{k+1}={\mathbf{p}}_u^k-\frac{\mu }{2}(\nabla \alpha )\top (\nabla \alpha uk+1-\mathbf{d}k+1),{\mathbf{p}}_d^{k+1}={\mathbf{p}}_d^k-\frac{\mu }{2}(\mathbf{d}k+1-\nabla \alpha uk+1)$

which can be simplified to

$minu,\mathbf{d}\{\frac{\lambda }{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+{\Re }_{\mu }^{\alpha }(\mathbf{d}(\mathbf{x}),u(\mathbf{x});\mathbf{p}(\mathbf{x}))\}$

(5)

${\Re }_{\mu }^{\alpha }(\mathbf{d}(\mathbf{x}),u(\mathbf{x});\mathbf{p}(\mathbf{x}))=\int \Omega |\mathbf{d}(\mathbf{x})|\text{d}\mathbf{x}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha u(\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|2\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\int \Omega |\mathbf{p}(\mathbf{x})|2|\mathbf{p}(\mathbf{x})|2\text{d}\mathbf{x}$

(6)

where the multiplier $\mathbf{p}(\mathbf{x})=(p1(\mathbf{x}),p2(\mathbf{x}))\top$ is updated by

$\mathbf{p}k+1(\mathbf{x})=\mathbf{p}k(\mathbf{x})-\gamma (\mathbf{d}(\mathbf{x})-\nabla \alpha u(\mathbf{x}))$

We may split equation (5) into two subproblems given by

$minu\{\frac{\lambda }{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha u(\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|)2\text{d}\mathbf{x}\}$

(7)

$min\mathbf{d}\{\int \Omega |\mathbf{d}(\mathbf{x})|\text{d}\mathbf{x}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha u(\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|2\text{d}\mathbf{x}\}$

(8)

It can be noted that equation (8) has a closed form solution. ²³ It remains to solve the minimisation problem (7) which can be achieved by finding the solution of the corresponding Euler–Lagrange equation

$(-1)r\mu C\text{div}\alpha (\nabla \alpha u(\mathbf{x})-\mathbf{d}(\mathbf{x})-\frac{\mathbf{p}(\mathbf{x})}{\mu })\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\lambda \kappa †(\mathbf{x})*([\kappa *u](\mathbf{x})-z(\mathbf{x}))=0$

Separation method for fitting and regularisation

Solving problems (5) and (6) by deriving optimality conditions and discretising the resulting Euler–Lagrange equations over an m × n grid $\Omega m,n$ requires forming a dense mn × mn matrix. This increases the computational time required to solve the non-linear equations considerably. We aim to improve the cpu time by introducing a dual function $v(\mathbf{x})$ which replaces the function $u(\mathbf{x})$ in the fractional regularisation functional $\Re$ . This allows us to separate the non-linear but sparse regularisation from the dense but linear data fitting problem. We solve the resulting minimisation problem subject to the constraint of equality of the image function and its dual. That is, we solve the problem subject to $u(\mathbf{x})-v(\mathbf{x})=0$ . We then transform this constrained problem to an unconstrained one by incorporating the equality requirement into the functional. This results in the unconstrained optimisation problem

$minu,v,\mathbf{d}\{\frac{\lambda 1}{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\lambda 2}{2}\int \Omega (u(\mathbf{x})-v(\mathbf{x}))2\text{d}\mathbf{x}+{\Re }_{\mu }^{\alpha }(\mathbf{d}(\mathbf{x}),v(\mathbf{x});\mathbf{p}(\mathbf{x}))\}$

(9)

where ${\Re }_{\mu }^{\alpha }$ is given above. We can now split equation (9) into three subproblems given by

$min\mathbf{d}\{\int \Omega |\mathbf{d}(\mathbf{x})|\text{d}\mathbf{x}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha v(\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|2\text{d}\mathbf{x}\}$

(12)

where equation (12) has a closed form solution. We find the functions $u(\mathbf{x})$ and v(x) which minimise equations (10) and (11), respectively, as the solutions of the respective Euler–Lagrange equations. The first is given by

$\lambda 1[\kappa †*\kappa *u](\mathbf{x})+\lambda 2u(\mathbf{x})=\lambda 1[\kappa †*z](\mathbf{x})+\lambda 2v(\mathbf{x})$

whose solution $u(\mathbf{x})$ is a minimiser of equation (10). Since the associativity and distributivity properties are satisfied by convolution and we can write $u(\mathbf{x})=[\delta *u](\mathbf{x})$ for the delta function $\delta (\mathbf{x})=1$ for $\mathbf{x}=0$ and $\delta (\mathbf{x})=0$ otherwise, this equation is linear in terms of the function $u(\mathbf{x})$ and can be solved efficiently using Fourier transforms. The Euler–Lagrange equation corresponding to equation (11), such that solution the $v(\mathbf{x})$ is a minimiser of the subproblem, is given by

Implicitly constrained deconvolution

The above models can give good performance but do not take into account the range of image intensity (or point spread function) values. It is typical in image restoration for values to fall outside of the expected range of the restored image which can have a significantly detrimental effect on the result, particularly in the case of blind deconvolution. Traditionally, this issue has been addressed by projection of the data by scaling or truncation at the end of or at points during the solution algorithm. This can result in the recovered image not being the minimiser of the optimisation problem and can lead to poor quality results. There have been several pieces of work recently^1,11,24–27 which consider methods of improving on this. We achieve constraints on both the upper and lower bounds of the image following the work of Williams et al. ²⁷ by introducing functions $\psi (\mathbf{x})$ and $\zeta (\mathbf{x})$ such that $u(\mathbf{x})=\tau \mathbf{a}(\zeta (\mathbf{x}))$ and $v(\mathbf{x})=\tau \mathbf{a}(\psi (\mathbf{x}))$ for a function $\tau \mathbf{a}$ which is naturally constrained. Such a function and its inverse $\xi \mathbf{a}$ are given as

$\tau \mathbf{a}(\psi )=\frac{a1+2a4}{1+a2e-\frac{2\psi }{a3}},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\xi \mathbf{a}(u)=-\frac{a3}{2}\ln \frac{a1-u+a4}{a2(u+a4)}$

Substituting this automatically constrained function for the image into the functional and including the condition that this must be equal to the image, we build the optimisation problem

$min\psi ,\zeta ,\mathbf{d}\{\frac{\lambda 1}{2}\int \Omega (\kappa *\tau \mathbf{a}(\zeta )-z)2\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\lambda 2}{2}\int \Omega (\tau \mathbf{a}(\zeta )-\tau \mathbf{a}(\psi ))2\text{d}\mathbf{x}+{\Re }_{\mu }^{\alpha }(\mathbf{d},\psi ;\mathbf{p})\}$

(13)

$\text{s}.\text{t}.u(\mathbf{x})-\tau \mathbf{a}(\zeta (\mathbf{x}))=0$

(14)

$v(\mathbf{x})-\tau \mathbf{a}(\psi (\mathbf{x}))=0$

(15)

Since both $\tau \mathbf{a}(\psi )$ and $\tau \mathbf{a}(\zeta )$ have the same range and we already have the requirement that these functions be equal, we may enforce the condition (14) explicitly in the functional (13) by replacing the range-constrained function $\tau \mathbf{a}(\zeta )$ with the unconstrained function $u(\mathbf{x})$ and rely on the subsequently implicitly imposed constraint of equality between u and $\tau \mathbf{a}(\psi )$ to enforce the restriction on the intensity range of u. Doing this, we can rewrite equations (13) to (15) as the unconstrained optimisation problem

$min\psi ,u,\mathbf{d}\{\frac{\lambda 1}{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+{\Re }_{\mu }^{\alpha }(\mathbf{d}(\mathbf{x}),\psi (\mathbf{x});\mathbf{p}(\mathbf{x}))\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\lambda 2}{2}\int \Omega (u(\mathbf{x})-\tau (\psi (\mathbf{x})))2\text{d}\mathbf{x}\}$

(16)

where ${\Re }_{\mu }^{\alpha }$ is given above. We can now split equation (16) into three subproblems given by

$min\psi \{\frac{\lambda 2}{2}\int \Omega (u(\mathbf{x})-\tau (\psi (\mathbf{x})))2\text{d}\mathbf{x}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha \psi (\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|2\text{d}\mathbf{x}\}$

(17)

$minu\{\lambda 12([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+\lambda 2\int \Omega (u(\mathbf{x})-\tau (\psi (\mathbf{x})))2\text{d}\mathbf{x}\}$

(18)

$min\mathbf{d}\{\int \Omega |\mathbf{d}(\mathbf{x})|\text{d}\mathbf{x}+\frac{\mu }{2}\int \Omega |\mathbf{d}(\mathbf{x})-\nabla \alpha \psi (\mathbf{x})+\frac{\mathbf{p}(\mathbf{x})}{\mu }|2\text{d}\mathbf{x}\}$

(19)

where equation (19) has a closed form solution and $\psi (\mathbf{x})$ and $u(\mathbf{x})$ which minimise equation (16) can be found by solving the corresponding Euler–Lagrange equations to equations (17) and (18), respectively, given by

$(-1)r\mu C\text{div}\alpha (\nabla \alpha \psi (\mathbf{x})-\mathbf{d}(\mathbf{x})-\frac{\mathbf{p}(\mathbf{x})}{\mu })\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\lambda 2\frac{\partial \tau \mathbf{a}(\psi (\mathbf{x}))}{\partial \psi (\mathbf{x})}(\tau \mathbf{a}(\psi (\mathbf{x}))-u(\mathbf{x}))=0,[(\lambda 1\kappa †*\kappa +\lambda 2\delta )*u](\mathbf{x})=\lambda 1[\kappa †*z](\mathbf{x})+\lambda 2\tau \mathbf{a}(\psi (\mathbf{x}))$

(20)

In order to solve the optimisation problem (16) for discrete data, we define a mesh $\Omega m,n=\{(x_1^i,x_2^j)|i\in [0,n+1],j\in [0,m+1]\}$ of the image domain Ω. At the point $(x_1^i,x_2^j)$ , the discretisation of the fractional derivative of order α is given along the x₁ direction as for $j\in [0,n+1]$ , where $u(x_1^i,x_2^j):=ui,j$ and where ${\theta }_k^{\alpha }$ is given by

${\theta }_k^{\alpha }=\{1\text{for}k=0(1-k-1(1+\alpha )){\theta }_{k-1}^{\alpha }\text{for}k>0$

Let Ψ denote the discrete form of the transformed image ψ to be restored, discretised over the mesh $\Omega m,n$ , and let $\Psi$ denote the corresponding vector given by ‘unpacking’ the matrix Ψ, that is $\Psi [s]=\Psi [i,j]$ where $s=i+n(j-1).$ From equation (21), we can write the discrete fractional α-order derivatives of Ψ in the x₁ and x₂ directions, respectively, as

$(Im\otimes B_n^{\alpha })\Psi ,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}(B_m^{\alpha }\otimes In)\Psi$

where $B_m^{\alpha }$ and $B_n^{\alpha }$ are symmetric and negative definite Toeplitz matrices derived from equation (21) as the derivative coefficient matrices. We can thus give the discrete equation to be solved for Ψ as

$\underset{\mathbf{w}\Psi }{\underbrace{((B_n^{\alpha })\top (B_n^{\alpha }\Psi )+\Psi (B_m^{\alpha })\top B_m^{\alpha })+\frac{(-1)r\lambda }{\mu }\Psi }}=\mathbf{f},\mathbf{f}=\frac{(-1)r\lambda }{\mu }\mathbf{u}+((B_n^{\alpha })\top D1+D2B_m^{\alpha })+\frac{1}{\mu }((B_n^{\alpha })\top P1+P2B_m^{\alpha })$

(22)

where $D1,D2,P1,P2$ are matrices given by the discretisation of the functions $d1,d2,p1,p2$ . The discrete restored transformed image Ψ can then be found by solving equation (22) using the conjugate gradient method.

In order to obtain the restored image $u(\mathbf{x})$ given the blurred and noisy image $z(\mathbf{x})$ and the blur function $\kappa (\mathbf{x})$ , we make an initial estimate of the restored image which, in the absence of a closer estimate, is typically taken as the blurred image. We give a value for the regularisation parameters and proceed with Algorithm 1 below. Results achieved using this method are presented in the ‘Experimental results’ section and compared with existing work.

Algorithm 1 Implicitly Constrained Fractional Deconvolution: $\mathbf{u}\to \text{A}\mathrm{ICFD}(\mathbf{z},\varepsilon ,\chi )$

1: Initialise $\mathbf{u}0\to \mathbf{z},\Psi 0\to \xi \mathbf{a}(\mathbf{u}0)$ and $P_1^0,P_2^0$ ;

2: for $\ell \to 1:\chi$ do

3: Solve (19) for d using the closed form solution:

$(D1D2)\ell +1\to \mathrm{shrink}((B_n^{\alpha }\Psi \ell +1\Psi \ell +1(B_m^{\alpha })\top )+(P1P2)\ell ,\frac{1}{\mu })$

4: Solve the subproblem (17) for the discrete solution Ψ by solving equation (22) using CG;

5: Solve the subproblem (18) by solving the discrete counterpart of the EL equation (20) for the discrete image u of u using discrete Fourier transforms;

6: Obtain the update for p:

$(P1P2)\ell +1\to (P1P2)\ell +\gamma ((B_n^{\alpha }\Psi \ell +1\Psi \ell +1(B_m^{\alpha })\top )-(D1D2)\ell +1),\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\gamma \in (0,1]$

7: Check stopping criterion:

8: if $\left|\right|\Psi \ell +1-\Psi \ell \left|\right|<\varepsilon$ then

9: break;

10: end if

11: end for

12: $\mathbf{u}\to \tau \mathbf{a}(\Psi \ell +1)$ ;

Given the excellent ability of this fractional regularised model to reconstruct both smooth regions and sharp edges in images, it would be interesting to consider its ability to reconstruct point spread functions. Since convolution is associative, we may do this by solving a similar optimisation problem to that given above by providing regularisation to the blur function. We present this problem as

$min\upsilon ,\kappa ,\mathbf{d}\{\int \Omega \lambda 12([u*\kappa ](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+{\Re }_{\mu }^{\alpha }(\mathbf{d}(\mathbf{x}),\upsilon (\mathbf{x});\mathbf{p}(\mathbf{x}))\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\lambda 22\int \Omega (\kappa (\mathbf{x})-\tau \mathbf{b}(\upsilon (\mathbf{x})))2\text{d}\mathbf{x}\}$

(23)

where $\upsilon (\mathbf{x})$ is the inverse transform of the point spread function such that we expect $\kappa (\mathbf{x})=\tau \mathbf{b}(\upsilon (\mathbf{x}))$ . By setting up a test problem where we have the blurred and noisy data $z(\mathbf{x})$ and the sharp true data $u(\mathbf{x})$ are known, we may attempt to identify the blur function which solves the optimisation problem (23) for the given z and u. This can be done with synthetic images for interest but we expect this to have a practical application to blind deconvolution, as considered in the following section.

A blind deconvolution model with fractional regularisation

Here we consider an application of the work given in the previous section to the case of blind deconvolution, where we do not know the cause of the blur degradation. That is, we are aiming to reconstruct the image $u(\mathbf{x})$ assuming that the blurred and noisy image data are given by $z(\mathbf{x})=[\kappa *u](\mathbf{x})+\eta (\mathbf{x})$ where κ is also unknown. Following the work of Chan and Wong ⁷ and You and Kaveh, ¹⁶ we form an optimisation model for blind deconvolution in the form where $\mathcal{L}1,\mathcal{L}2$ denote regularisation functions for the image and blur function, respectively, and $\alpha 1,\alpha 2\in ℝ>0$ are small, positive regularisation parameters which measure the trade-off between data fitting and regularisation. The blur and image functions are sought subject to the constraints

$\kappa (\mathbf{x})\ge 0,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}u(\mathbf{x})\ge 0,\int \Omega \kappa (\mathbf{x})\text{d}\mathbf{x}=1,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\kappa (\mathbf{x})=\kappa (-\mathbf{x})$

(24)

which are imposed in order to reduce the space of functions which may solve this problem since it is not jointly convex. Extending this idea for blind deconvolution to the work of the previous section, we incorporate the fractional regularisation and implicit constraints of equations (16) and (23). We present this new implicitly constrained, fractional regularised blind deconvolution problem as

$min\psi ,u,\mathbf{d},\upsilon ,\kappa ,\mathbf{e}\{\frac{\lambda 1}{2}\int \Omega ([\kappa *u](\mathbf{x})-z(\mathbf{x}))2\text{d}\mathbf{x}+{\Re }_{\mu 1}^{\alpha 1}(\mathbf{d}(\mathbf{x}),\psi (\mathbf{x});\mathbf{p}(\mathbf{x}))\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\lambda 21\int \Omega (u(\mathbf{x})-\tau \mathbf{a}(\psi (\mathbf{x})))2\text{d}\mathbf{x}+{\Re }_{\mu 2}^{\alpha 2}(\mathbf{e}(\mathbf{x}),\upsilon (\mathbf{x});\mathbf{q}(\mathbf{x}))\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\lambda 22\int \Omega (\kappa (\mathbf{x})-\tau \mathbf{b}(\upsilon (\mathbf{x})))2\text{d}\mathbf{x}\}$

(25)

Since the positivity constraints of equation (24) are automatically satisfied, we require the solution of this problem subject to only the symmetry and unit integral constraints on the blur function. Minimising the functional of equation (25) with respect to the arguments, there exist closed form solutions to the corresponding subproblems of minimising equation (25) with respect to $\mathbf{d},\mathbf{e}$ . We derive the Euler–Lagrange equations whose solutions are minimisers of equation (25) as

$\psi (\mathbf{x}):\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}(-1)r{\mu }_1^C\text{div}\alpha 1(\nabla \alpha 1\psi -\mathbf{d}-\frac{\mathbf{p}}{\mu 1})+\lambda 21\frac{\partial \tau \mathbf{a}(\psi )}{\partial \psi }(\tau \mathbf{a}(\psi )-u)=0$

(26)

$u(\mathbf{x}):[(\lambda 1\kappa †*\kappa +\lambda 21\delta )*u](\mathbf{x})=\lambda 1[\kappa †*z](\mathbf{x})+\lambda 21\tau \mathbf{a}(\psi (\mathbf{x}))$

(27)

$\upsilon (\mathbf{x}):\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}(-1)r{\mu }_2^C\text{div}\alpha 2(\nabla \alpha 2\upsilon -\mathbf{e}-\frac{\mathbf{q}}{\mu 2})+\lambda 22\frac{\partial \tau \mathbf{b}(\upsilon )}{\partial \upsilon }(\tau \mathbf{b}(\upsilon )-\kappa )=0$

(28)

$\kappa (\mathbf{x}):[(\lambda 1u†*u+\lambda 22\delta )*\kappa ](\mathbf{x})=\lambda 1[u†*z](\mathbf{x})+\lambda 22\tau \mathbf{b}(\upsilon (\mathbf{x}))$

(29)

Given initial estimates of the image and the blur function, we proceed to solve the minimisation problem following Algorithm 2 below.

Algorithm 2 Blind Fractional Deconvolution: $\mathbf{u},\mathbf{k}\to \text{A}\mathrm{BFD}(\mathbf{z},\varepsilon 1,\varepsilon 2,\chi ,\chi 1,\chi 2)$

1: Initialise $\mathbf{u}0,0\to \mathbf{z},\Psi 0,0\to \xi \mathbf{a}(\mathbf{u}0,0)$ and $P_1^{0,0},P_2^{0,0}$ ;

2: Initialise $\mathbf{k}0,0\to \delta ,\mathbf{Y}0,0\to \xi \mathbf{b}(\mathbf{k}0,0)$ and $Q_1^{0,0},Q_2^{0,0}$ ;

3: for $\ell \to 1:\chi$ do

4: for $\ell 1\to 0:\chi 1$ do

5: Solve the subproblem of equation (25) for $\mathbf{d}\ell ,\ell 1+1$ using the closed form solution;

6: Solve the subproblem of equation (25) for $\Psi \ell ,\ell 1+1$ by solving the discrete form of equation (26) using CG;

7: Solve the subproblem of equation (25) for $\mathbf{u}\ell ,\ell 1+1$ by solving the discrete counterpart of the EL equation (27) for the discrete image u of u using discrete Fourier transforms;

8: Obtain the update for $\mathbf{p}\ell ,\ell 1+1$

$(P1P2)\ell ,\ell 1+1\to (P1P2)\ell ,\ell 1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\gamma 1((B_n^{\alpha 1}\Psi \ell ,\ell 1+1\Psi \ell ,\ell 1+1(B_m^{\alpha 1})\top )-(D1D2)\ell ,\ell 1+1),\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\gamma 1\in (0,1]$

9: Check stopping criterion:

10: if $\left|\right|\Psi \ell ,\ell 1+1-\Psi \ell ,\ell 1\left|\right|<\varepsilon 1$ then

11: break;

12: end if

13: end for

14: $(\mathbf{u},\Psi ,P1,P2)\ell +1,0\to (\mathbf{u},\Psi ,P1,P2)\ell ,\ell 1+1$ ;

15: for $\ell 2\to 0:\chi 2$ do

16: Solve the subproblem of equation (25) for $\mathbf{e}\ell ,\ell 2+1$ using the closed form solution;

17: Solve the subproblem of equation (25) for $\mathbf{Y}\ell ,\ell 2+1$ by solving the discrete form of equation (28) using CG;

18: Solve the subproblem of equation (25) for $\mathbf{k}\ell ,\ell 2+1$ by solving the discrete counterpart of the EL equation (29) for the discrete image k of κ using discrete Fourier transforms;

19: Obtain the update for $\mathbf{q}\ell ,\ell 2+1$

$(Q1Q2)\ell ,\ell 2+1\to (Q1Q2)\ell ,\ell 2\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\gamma 2((B_n^{\alpha 2}\mathbf{Y}\ell ,\ell 2+1\mathbf{Y}\ell ,\ell 2+1(B_m^{\alpha 2})\top )-(E1E2)\ell ,\ell 2+1),\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\gamma 2\in (0,1]$

20: Check stopping criterion:

21: if $\left|\right|\mathbf{Y}\ell ,\ell 2+1-\mathbf{Y}\ell ,\ell 2\left|\right|<\varepsilon 2$ then

22: break;

23: end if

24: end for

25: Enforce the unit integral and symmetry constraints

$\mathbf{k}\ell +\frac{1}{2},\ell 2+1\to \frac{\mathbf{k}\ell ,\ell 2+1}{\sum _i\sum _j{\mathbf{k}}_{i,j}^{\ell ,\ell 2+1}},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{\mathbf{k}}_{i,j}^{\ell +1,\ell 2+1}\to \frac{{\mathbf{k}}_{i,j}^{\ell +\frac{1}{2},\ell 2+1}+{\mathbf{k}}_{-i,-j}^{\ell +\frac{1}{2},\ell 2+1}}{2}\forall i,j$

26: $(\mathbf{Y},Q1,Q2)\ell +1,0\to (\mathbf{Y},Q1,Q2)\ell ,\ell 1+1$ ;

27: end for

28: $\mathbf{u}\to \tau \mathbf{a}(\Psi \ell +1,\ell 1+1)$ ;

29: $\mathbf{k}\to \tau \mathbf{b}(\mathbf{Y}\ell +1,\ell 1+1)$ ;

Experimental results

In this section, we aim to show that the implicitly constrained fractional deconvolution method introduced in the section ‘A fractional regularised image deconvolution model’ is capable of providing good quality deconvolution results and that it can outperform competing methods for image deblurring. We also demonstrate the ability of the blind fractional model described in the section ‘A blind deconvolution model with fractional regularisation’ to restore image data from blur corruption without knowledge of the blur function. For comparison purposes, we attempt to restore a set of 10 test images consisting of images with sharp and smooth regions (examples 1–4), images which are predominantly smooth (examples 5–8) and images corresponding to Gaussian blur functions (examples 9–10). Examples were restored using the total variation model (3) with results denoted by $u_{\mathrm{tv}}^{\beta }$ for $\beta =10-3,10-6$ which are common values for this parameter, using the mean curvature model by employing the regulariser given by equation (2) with restoration results denoted by u_mc. For the fractional deconvolution model, the order α was considered in the interval $[1.1,1.9]$ and corresponding results are denoted $u_f^{\alpha }$ . Tests were carried out for regularisation parameters in the range $[10-9,0)$ and the parameter providing the best quality results was used in each case. All experiments were carried out on a Dell XPS 8700 with an Intel Core i7-4770 processor and 32GB RAM and a tolerance of $10-8$ was used in each case for the stopping criteria. Where possible, the error of the m × n discrete image u is measured comparative to the true image $\overline{\mathbf{u}}$ using the peak signal-to-noise ratio (PSNR) ²⁸ $\mathcal{P}(\mathbf{u})$ given as

$\mathcal{P}(\mathbf{u})=20\log 10(\frac{max(\mathbf{u}[i,j])}{\frac{1}{\mathrm{mn}}\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}(\mathbf{u}[i,j]-\overline{\mathbf{u}}[i,j])})$

It can be noted in Table 1 and Figures 1 to 3 that the fractional model consistently provides good deconvolution results and is capable of outperforming competing models. It has been observed that, in the case of zero-mean Gaussian noise, as the level (variance) of the noise increases, the fractional model continues to perform well and to outperform total variation and mean curvature for predominantly sharp and smooth images. It can also be noticed that the fractional model has been shown to outperform models which are typically suited to problems with sharp jumps in intensity as well as to those with smoother gradients.

Figure 1.

Figure showing the results of restoring Example 1 using the total variation model for different β $(u_{\mathrm{tv}}^{\beta })$ , mean curvature $(u\mathrm{mc})$ and the fractional model for different α $(u_f^{\alpha })$ . It can be seen that the fractional results, particularly for $\alpha =1.5$ give improved quality results and the fractional model is able to preserve both smooth and sharp regions. (a) True image. (b) Received image z. (c) $u_{\mathrm{tv}}^{10-3}$ . (d) $u_{\mathrm{tv}}^{10-6}$ . (e) $u\mathrm{mc}$ . (f) $u_f^{1.3}$ . (g) $u_f^{1.4}$ . (h) $u_f^{1.5}$ .

Figure 2.

Figure showing the results of restoring the Gaussian blur function Example 9 using the total variation model for different β $(u_{\mathrm{tv}}^{\beta })$ , mean curvature $(u\mathrm{mc})$ and the fractional model for different α $(u_f^{\alpha })$ . It is clear from the cross section images (g)–(l) that tv does not achieve very good results for this smooth function, as expected. Competition between mean curvature and the fractional model is close but the numerical results $\mathcal{P}(u\mathrm{mc})=26.32$ and $\mathcal{P}(u_f^{1.3})=26.69$ demonstrate that the fractional model outperforms mc. (a) True image. (b) Received image z. (c) $u_{\mathrm{tv}}^{10-3}$ . (d) $u_{\mathrm{tv}}^{10-6}$ . (e) $u\mathrm{mc}$ . (f) $u_f^{1.3}$ . (g) True Image. (h) Received Image z. (i) $u_{\mathrm{tv}}^{10-3}$ . (j) $u_{\mathrm{tv}}^{10-6}$ . (k) $u\mathrm{mc}$ . (l) $u_f^{1.3}$ .

Figure 3.

Example of the satellite image corrupted by out of focus blur (a) and restored using the implicitly constrained fractional blind deconvolution algorithm $\text{A}\mathrm{BFD}$ . Details in the reconstructed image have been restored and it is clear that smooth regions and sharp edges have been well reconstructed. (a) Received image z. (b) Restored image $u\mathrm{BFD}$ .

Table 1.

Table of error values for the restored images achieved using the total variation model for different β $(u_{\mathrm{tv}}^{\beta })$ , mean curvature $(u\mathrm{mc})$ and the fractional model for different α $(u_f^{\alpha })$ . Error was calculated using PSNR given above. Examples 1–4 involve images with sharp and smooth regions, examples 5–8 are predominantly smooth images and examples 9–10 involve images resembling Gaussian blur functions. All images were corrupted with Gaussian blur and 10% noise. In each case, the fractional model outperforms the competing models. Bold values indicate the best PSNR achieved for each example.

Example	z	$u_{\mathrm{tv}}^{10-3}$	$u_{\mathrm{tv}}^{10-6}$	$u\mathrm{mc}$	$u_f^{1.3}$	$u_f^{1.4}$	$u_f^{1.5}$	$u_f^{1.6}$	$u_f^{1.7}$
1	33.02	35.07	33.58	37.98	42.69	42.63	42.66	42.55	42.41
2	29.03	31.13	30.24	31.90	34.77	34.77	34.72	34.29	34.13
3	24.18	24.42	24.29	24.92	26.65	26.70	26.66	26.64	26.39
4	23.40	23.93	23.82	24.19	24.76	24.76	24.76	24.70	24.65
5	44.30	45.17	44.20	45.83	45.95	45.96	45.97	46.10	46.11
6	40.87	45.11	44.22	45.69	45.47	45.74	45.98	45.98	45.99
7	44.71	45.66	45.10	45.98	46.32	46.32	46.22	46.16	46.14
8	44.30	45.63	45.13	45.93	46.17	46.26	46.17	46.11	46.09
9	25.40	26.12	26.09	26.32	26.69	26.68	26.65	26.60	26.54
10	25.17	26.05	25.97	26.21	26.60	26.59	26.55	26.51	26.39

Conclusion

In this paper, we have considered the problem of image deconvolution where the image to be restored may contain both sharp and smooth regions. To achieve this, we have built a new variational model for image deconvolution using fractional derivatives in order to provide convex regularisation for the image. We have also extended this to include intensity range considerations by incorporating them implicitly in the problem functional and formed a solution algorithm using variable splitting. Given the good results achieved by attempting to restore blur functions using fractional regularisation, we extended this idea to the case of blind deconvolution where we defined a model for restoring the image without knowledge of the blur function, using fractional regularisation for the image and blur function. We have demonstrated that the work introduced in this paper is capable of achieving good results for image deconvolution and that it can outperform competing models.

Considering future work, the choice of fractional order used for experimental testing was chosen empirically. While this allowed for good results to be achieved in reasonable time, it would be useful to consider how the order might be chosen automatically to provide the optimal result in terms of result quality. It was observed during experimental testing that the computation time of the fractional model was comparable to that of the total variation model when formed in a similar manner. Given an effective order-selection technique, further work may be carried out to investigate the computation time and complexity more thoroughly.

This work may be further extended to the case of semi-blind deconvolution^4,29,30 by relaxing the space of blur functions considered to only those which may be given by a set of parametric functions. We may also extend this idea to colour images ³¹ and the case of deconvolution in the presence of Poisson noise.^3,12