International Journal of Computational Mathematics
Volume 2015 (2015), Article ID 860263, 17 pages
http://dx.doi.org/10.1155/2015/860263

Research Article

A New Study of Blind Deconvolution with Implicit Incorporation of Nonnegativity Constraints

Ke Chen,¹ Simon P. Harding,² Bryan M. Williams,¹ and Yalin Zheng²

¹Centre for Mathematical Imaging Techniques and Department of Mathematical Sciences, University of Liverpool, Liverpool L69 L7L, UK
²Department of Eye and Vision Science, University of Liverpool, Liverpool L69 3GA, UK

Received 1 July 2014; Revised 29 December 2014; Accepted 12 January 2015

Academic Editor: David Defour

Copyright © 2015 Ke Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The inverse problem of image restoration to remove noise and blur in an observed image was extensively studied in the last two decades. For the case of a known blurring kernel (or a known blurring type such as out of focus or Gaussian blur), many effective models and efficient solvers exist. However when the underlying blur is unknown, there have been fewer developments for modelling the so-called blind deblurring since the early works of You and Kaveh (1996) and Chan and Wong (1998). A major challenge is how to impose the extra constraints to ensure quality of restoration. This paper proposes a new transform based method to impose the positivity constraints automatically and then two numerical solution algorithms. Test results demonstrate the effectiveness and robustness of the proposed method in restoring blurred images.

1. Introduction

Among image preprocessing problems is the reconstruction of an image from a given degraded image, such as images corrupted by noise [1, 2] or blur [3, 4] or images with missing or damaged portions [5]. Such tasks have been widely studied in the last few decades; see [6] for decoupling noise and blur modeling, [7] for imposing box constraints, [8] for a fast iterative solver for noise and blur modeling, and [9, 10] for general surveys. However there are still many outstanding issues to be addressed, especially when both the noise type and the blur type are unknown. Image restoration is closely related to higher level tasks such as segmentation [11] and registration [12, 13]. It should be remarked that models for the latter tasks often become ineffective if the underlying image is blurred.

Adopting the usual notation from the literature, assume an observed image function in domain has been contaminated with additive noise and convolution (blur) operator :

where is an unknown Gaussian white noise with zero mean and is the image to be restored. When is known, there exist many effective studies to reconstruct the restoration; see, for example, [14, 15] for fixed point methods, [16] for a Krylov conjugate gradient method, and [17] for a multilevel method. However when is unknown, simultaneous restoration of is of great interest.

A model for blind deblurring by recovering both the kernel and the image simultaneously with no a priori information was given by You and Kaveh [18] and later improved by Chan and Wong [3]. In the latter paper, they proposed an energy minimising model, derived partial differential equations by minimising with respect to the image and the kernel , and presented an alternate minimisation scheme for solving the model. The model has extra constraints, including nonnegativity constraints, which were crucial but not exactly implemented. It is fair to say that the model by [3] is not yet reliable for general use (as remarked by [19, 20]) and there are no substantial improvements of it since 1998. There are several works trying to adapt for and extend to specific applications; see, for example, [20, 21] for using multichannel images to restore a single image, [22] for using two blurred images to restore an image, [23] for a nonvariational method, and [24] for implementing [3] by a splitting method. The particular ideas of leaving out the constraints altogether and trusting that the model will give a good result are unfortunately unreliable since they do not lead to good results.

However it is known [19, 20] that the general model [3] can only deal with very special images where the kernel can be accurately estimated. Recognizing the importance of nonnegativity constraints, Miura [25] generalized a one dimensional idea of using square functions from [26] to image deblurring and attempted to solve blind deconvolution in Fourier domain (but without any regularization). Šroubek and Milanfar [20] generalized the model [3] to the case of having multichannel blurred images of the same true image and incorporated nonnegativity into an minimization energy.

In this paper, we focus on image deblurring in the blind case where the blur operator is unknown. In this case, we are aiming to reconstruct the true image and the cause of the degradation with no prior information. Of course, if extra information of the blur operator is available, it should be used to derive the so-called semiblind models [27, 28]. Although there exist other approaches [23, 29–31] for deblurring, we shall focus on the variational framework to model the single channel image deconvolution through satisfying the nonnegativity constraints exactly and implicitly. The end product is a robust image deblurring model; we also present two methods of solving it.

The rest of this paper is organised as follows. Section 2 reviews four related variational models. Two test examples are shown to illustrate and highlight the problems and challenges faced by the first and earlier model by [3]. Section 3 first introduces our transformation approach where both the image and the blurring function are reconstructed with nonnegativity imposed implicitly and then describes the numerical solution of the model. Section 4 presents some experimental results. Section 5 concludes the paper.

2. The Inverse Problem of Deblurring and Some Current Models

Here we review some blind deconvolution models before we introduce our method in the next section. As seen shortly, imposing the constraint of nonnegativity is crucial for such models.

Before proceeding, we remark that for the traditional image restoration applying a projection is the simplest idea of imposing nonnegativity . The same projection idea can be applied to impose a box constraint to ensure . See [7, 32].

There have been several other related ideas for enforcing nonnegativity in image processing [32–35]. One such example was given by [32] where a model was proposed for image reconstruction using nonnegative constraints for astronomical imaging by minimising a regularised Poisson likelihood functional while the idea of backprojection is similarly used in [33, 35]. The case of a Tikhonov regularisation (a much simpler regulariser than what we use here) was considered in [36]. The method of [34] ensured a positive kernel by considering a parametric model and optimizing a scalar which is the standard deviation of the point spread function.

A more sophisticated idea by Biraud [26] is to use the transform , with , in restoring one dimensional signal from the model , where is the known blur function and the noise, or after Fourier transform. The central idea here is that any or its Fourier transform can lead to nonnegative restoration . For with some cut-off frequency , so . Noting that leads to , the method of Biraud [26] is

To solve (2), a parametric iterative approach is proposed by for . See [26]. Once a good approximation is obtained, an inverse transform would yield and then a nonnegative restoration .

2.1. The Earlier Blind Deconvolution Work

You and Kaveh [18] proposed a model for simultaneous recovery of both the degradation function and the image in a variational framework, by solving the problem

where the fitting term is a common choice for (1) and there is freedom to choose the regularisation terms and ; You and Kaveh used the seminorm for these two terms.

Chan and Wong [3] proposed an improvement to this using the total variation (TV) seminorm for regularisation given by and , hence solving

where . Minimising (4) with respect to the image and the kernel , we obtain the coupled partial differential equations given by

where (similarly for also) with a small positive parameter introduced to avoid division by zero, which is set to for experimental results for both the image and blur funcion in order to avoid overly smooth reconstructions resulting from too high and the staircasing effect arising from setting too low. It is worth noting that as alternatives to the total variation seminorm we may also consider other regularisation terms such as regularisation, given by , the nonlocal TV [37, 38], the total generalised variation (TGV) [39–41], or the mean curvature [42, 43], as well as others [44, 45].

In order to solve the system, an alternate minimisation scheme was proposed to recover the kernel and the image , including the following constraints which aim to deal with the lack of a unique solution since the system is not jointly convex. This leads to imposing the constraints that the image and kernel should both be positive, and , the kernel should be symmetric (), and the kernel should have a unit integral . These constraints are imposed exactly but only after each alternate minimisation step. The complete algorithm is given in Algorithm 1.

(1) function  CHAN-WONG  ()

(2)      for    do

(3)            

(4)            

(5)            

(6)            

(7)            

(8)            

(9)      end for

(10) end function

Algorithm 1: Chan-Wong [3] method for the solution of (4).

Adding the above constraints ensures a unique solution but not imposing them in the algorithm introduces inconsistency which is problematic. To test this remark, the algorithm yields a reasonable result for the example given in Figure 1(c), but the same algorithm gives rise to poor results such as Figures 2(c) and 2(d) due to this inconsistency. We may attempt to improve the results by implementing a better thresholding technique applied to the kernel. To do this, we adjust the filter in step of Algorithm 1 to be dependent on a small positive parameter as follows:

This adjustment with a problem dependent parameter may offer some improvement (Figure 2(e)) but does not always lead to a good solution. Such problems were reported in other subsequent studies [19, 20].

(a) True data

(b) Blurred image

(c) Algorithm 1 result

Figure 1: Good restoration results for Example 1 (a) from a corrupted image (b) using Algorithm 1. This model is able to improve the edges of the restored image (c), though the restoration is not excellent.

(a) Ex2 true image

(b) Ex2 received data

(c) Failed restoration using Algorithm 1

(d) Restored with

(e) Restored with

(f) Restored with

Figure 2: Illustration of the failure of Algorithm 1 for a retinal image (Example 2). (a) True image; (b) corrupted image by out of focus blur; (c) failed restoration ; (d) restored with thresholding ; (e) restored with thresholding ; (f) restored with thresholding .

Our aim is to satisfy exactly these constraints by achieving the positivity on the kernel and the image in the functional in an implicit manner.

2.2. The Blind Deconvolution Model by Miura [25]

Following the work of Biraud [26], Miura [25] considered generalizing it to the image case and more importantly to the blind deconvolution problem by imposing nonnegativity for both and . Starting from (1), that is, , with both and as unknowns, he defined . Then after Fourier transforms, one gets

Further similar to (2), it is proposed to solve

where the summations imply formulations after discretization [25]. Further a conjugate gradient type solver is utilised to compute which will be used to yield the nonnegative solutions .

2.3. A Shock Filter Based Model by Money and Kang [19]

To improve the method of [3], in particular the algorithm (5), an interesting idea was proposed in [19] to decouple the equations so that edge information of the restoration is ensured. Precisely in the second equation of (5) is replaced by a reference image which is obtained by using a shock filter to capture image edges in the blurred . Then (5) becomes

which is a decoupled system and can be solved directly in a noniterative way between and .

2.4. A Multichannel Blind Deconvolution Model by Šroubek and Milanfar [20]

To overcome the poor performance of [3], it is suggested in [20] that there may be a better chance of restoring the blurred image if blurred images of the same object are available which is readily possible for video images in some situations. The minimization proposed is

where in a discrete setting

with a Laplacian related operator, and is a parameter. Here denotes the total variation regularizer for and the crucial choice of ensures positivity of . However the same treatment was not applied to . The optimization problem was further solved by a splitting idea in an augmented Lagrangian method (ALM).

3. A Refined Blind Model and Its Numerical Solution Methods

We now consider the single image blind deconvolution problem (1) and propose a way to improve the Algorithm 1 by Chan and Wong [3], through use of a related and different idea from [25, 26]. The similarity to [25, 26] lies in that, instead of treating negative components directly as in a projection method, we seek a transform that converts the original model into a new one that can satisfy the nonnegativity constraints. There are three clear differences: (i) we use a different transform from previous choices; (ii) we apply regularization to the restored quantities while previous work use nonlinear least squares fitting without regularization; (iii) we solve for directly instead of solving for in the Fourier domain.

3.1. Choice of Positivity Transforms

We aim to impose nonnegativity in the functional by representing the kernel and the image as transformed quantities which do not permit negative values. One such idea might be to represent the image as the exponential function; that is for some function . Unfortunately this particular transform does not work as it is not capable of dealing with dark regions (where ) in a stable way. Another choice could be as in [25, 26]. This is valid for but this is unbounded.

Since our aim is to bound the function’s upper and lower values, we consider a generalisation of a differentiable approximation of the Heaviside step function. Thus, a suitable and bounded transform can be given by

where constants , , and is a small tuning parameter which controls the spread of the function. Clearly poses no problems to the transform. For the usual intensity range for , we may take , , :

Since for and for are monotone functions, we can work out their lower and upper bounds:

Note, for (13), and . It should be noted that since our model assumes that the image is represented by the function, the choice of parameters does not have to be altered to take into account varying levels of noise and blur with the exception of the scaling parameter of the blur function.

3.2. Reformulation of the Blind Deblurring Model

In order to apply the same transform to both the image and the kernel , we introduce the parameters with subscripts as follows:

for the image and kernel, respectively; here all constants can be fixed before proceeding (see Appendix A for more details). In particular, and are the expected upper limits of the image intensity values and kernel values; and are introduced to control the values of the image and kernel at and , and and control the spread of and . To give one feasible set, for image , we have if and for kernel we take if . That is we use

Similar to the previous section, we need not vary the parameters to take into account varying noise or blur levels with the exception of which may be chosen depending on the perceived level of blur. We now reformulate our old problem (4) as the new variational model:

letting the image and the kernel be represented by and , respectively. Here from solving (17), the nonnegativity constraints are exactly and implicitly enforced; that is , but the remaining symmetry and unit integral constraints on the kernel are still required.

The advantage of realizing positivity (for any ) is accompanied by a new challenge (or disadvantage) of having to deal with a nonlinear convolution kernel in (17). We next present two methods for solving the model and below show the solution method for to illustrate the idea as the solution for is similar.

3.3. Solution I by a Fixed Point Method

Our first method will deal with the nonlinearity in directly. To construct an iterative scheme, we consider a linear approximation of the transform by the Taylor expansion given by . First, defining , we split the transform by

where is nonlinear. Second, treating the above two terms differently through lagging, we propose a fixed-point lagging technique by substituting into (17) and get it linearised, leaving the remaining nonlinearity in terms with known quantities. Similarly we also have for . Clearly and . In alternating minimization, this converts problem (17) to

Our idea is to repeat the iterations until is small. Here a key observation is that residual functions are lagged in iterations, not approximated in any way.

Minimising the above functionals from (19) with respect to and , given and , yields the Euler Lagrange equations (see Appendix B for the derivation):

which are linearised equations due to linear functions ,

or simplified as, after moving known terms to the right hand sides,

where , ,

After estimating and of the image and the kernel, respectively, we apply the inverse transforms obtaining and . Then on convergence, are finally obtained from .

3.3.1. Constraints for

Note that the previously mentioned constraints and take the new forms: and . We can satisfy the first condition by imposing

where is the result of the previous step. For the second constraint in the discrete setting, we interpret the integral of a function over the domain as the sum of all values of the function in . Letting , our constraint is satisfied by from solving

and then . The solution method is given in Algorithm 2.

(1) function  TRANSFORM_METHOD_1 ()

(2)     ,

(3)    for    do

(4)          (20)

(5)          

(6)          

(7)          Obtain from solving (25) with .

(8)          (20).

(9)          If   , then exit or continue.

(10)    end for

(11)    Accept the restored image and the restored kernel .

(12) end function

Algorithm 2: The first transform method for model (17) via (20).

3.3.2. Numerical Discretisation

We shall briefly discuss the discretisation of the linearised operators in (20) by a finite difference method. We adopt the zero Dirichlet boundary conditions for the original variables which become and . That is, we have nonzero Dirichlet boundary conditions for the transformed variables .

Discretisation of the Fitting (Integral) Term. We wish to discretise the quantity which is similar to and as discussed in [46, 47] because all such operators are of the same type (spatially invariant, i.e., convolution). With Dirichlet boundary conditions for , we shall obtain Block-Toeplitz-with-Toeplitz-Blocks (BTTB) structure; however other boundary conditions can also be considered leading to different but similar structures [48, 49]. It should be noted that parameters may be chosen to suit the assumption of Dirichlet boundary conditions but the choice is not influential for boundary conditions such as Neumann or periodic.

Discretisation of the Total Variation Regularisation (Differential) Term. Using central differences, we shall derive a linearised matrix with 5-diagonal coefficient matrix structure [46, 47].

Thus after discretization, (22) leads to

where , are of BTTB form and and denotes the discretised TV operator for (22).

3.3.3. Iterative Solution of Linear Systems

We now consider a solution method for solving our system (26) of discrete versions of linearised PDEs (22). We use a preconditioned conjugate gradient algorithm; see [46, 47]. In order to improve the speed of convergence, we make use of preconditioners and . We implement the product preconditioner following the work of [46] and given by

for (26), where is a circulant approximation [46, 47] to defined above and are positive constants.

3.4. Solution II by an Augmented Lagrangian Method

The central idea of our second method for model (17) is to maintain the linear deblurring in the original variables or to remove the nonlinearity in the fitting term while still imposing the two proposed transforms , .

This is achieved by treating the transforms as constraints for and utilizing an augmented Lagrangian formulation, as was done in modeling other problems [6, 7, 43, 50]. More importantly for our transformed formulation, as remarked, the nonlinearity introduced by the transforms to the blurring term is removed by the method.

Starting with the unconstrained nonnegative problem given by (4) we propose the augmented Lagrangian minimising functional:

where and are regularisers representing either total variation (where we expect jumps in intensity) or L2 regularization where we expect smooth edges. Nonnegativity is imposed implicitly by the transform and by and for the image and kernel, respectively, which force them to be close to their respective nonnegative representations.

Minimising with respect to each of the arguments, we have

where , are the same as in (22):

We use alternate minimisation to solve the minimisation problem of the augmented Lagrangian functional. We solve the first equation in (29) efficiently using Fourier transforms and employ an iterative technique to solve the other equations of (29). We present our algorithm in Algorithm 3.

(1) function  TRANSFORM_METHOD_2

(2)      , ,

(3)     Calculate , and

(4)     

(5)     

(6)      for    to    do

(7)            Solve   (29) for , i.e.

(8)            for    to    do

(9)                  Solve   (29) for given , i.e.

(10)          end for

(11)            

(12)          Solve   (29) for , i.e.

(13)          for    to    do

(14)                Solve   (29) for given , i.e.

(15)          end for

(16)          

(17)      end for

(18) end function

Algorithm 3: The second transform based algorithm for positivity with ADM.

The above presented two methods realized Chan-Wong original model [3]. As demonstrated shortly, the realization is so effective that even problems beyond the intended tasks of motion blur or out of focus blur (which are piecewise constants) can also be modeled. In fact, the TV seminorm gives acceptable restoration for but it is less ideal for smooth ; this prompts us to consider a simple generalization model. Of course, a better model would be to use a high order regulariser which is capable of restoring both nonsmooth and smooth functions (such as a mean curvature [42]).

3.5. A Mixed Model Suitable for Smooth Blur Kernel

In an attempt to improve the result of recovered smooth (such as Gaussian) kernels, we introduce the following functional which uses the norm to regularise the blurring kernel and the TV to regularise the image , as a hybrid model of (3) and (4):

Then the Euler-Lagrange equations corresponding to Method I of Section 3.3 can be derived

where , ,

Other computational details and use of Method II can be considered similarly.

4. Experimental Results

The aim of our experimental tests is to demonstrate the effectiveness of our new transform model (17) for restoring both the image and the kernel , given the received image . The results will illustrate the capability of our new algorithm for potentially wide applications. Comparison with previous and competing methods have been shown earlier along with some comparisons here; here we primarily aim to show how the new algorithm can better restore these examples and a range of images.

For experimental testing, we consider test images in sets. Our aim is to show that our model is able to recover the edges of images as well as many of the fine details, in cases of both piecewise constant blurs and the more challenging Gaussian blur. We demonstrate this using a piecewise constant out of focus blur function of radius 4 and a Gaussian blur function of support equal to the size of the image and of standard deviation 5. All images have been corrupted by 3% random noise. For each algorithm, residual tolerances have been set at with maximum iterations although this number is rarely reached.

4.1. Test of Robustness of Algorithm 2

Set 1 (simple image with blurs). Result Set 1 consists of a synthetic (artificial) image corrupted by out of focus blur (Blur 1) or Gaussian blur (Blur 2). We see in Figure 3 that our model is able to reconstruct the edges and preserve the smoothness of the images in the case of out of focus blur and offers a significant improvement in the case of Gaussian blur.

(a) True image

(b) Received SNR = 13.5

(c) Algorithm 1 [3] SNR = 9.8

(d) Algorithm 2 SNR = 21.9

(e) True image

(f) Received PSNR = 7.9

(g) Algorithm 2 PSNR = 27.3

Figure 3: Experimental results on Image 1 with Blur 1. Row 1, l-r: Image 1, received data corrupted by Blur 1, restored image using Algorithm 2. The PSNR/SNR is lowered by Algorithm 1 by 2.25/3.7 from 20.15/13.54 to 17.90/9.84 but increased using Algorithm 2 by 8.21/8.4 to 28.36/21.94. Row 2, l-r: Image 1, received data corrupted by Blur 2, restored image using Algorithm 2. The PSNR/SNR is increased by 19.44/34.53 from 7.88/−13.4 to 27.32/21.13. Our model is capable of restoring edges and preserving the background in black.

Set 2 (detailed image containing many zero-points with blurs). Result Set 2 consists of Image 2 corrupted by out of focus blur (Blur 1) or Gaussian blur (Blur 2). We see in Figure 4(f) that our model is able to preserve the black space well and reconstruct details in the case of out of focus blur. It is also able to restore some detail in the more challenging case of Gaussian blur.

(a) True image

(b) Received data

(c) Algorithm 2 restored image

(d) True image

(e) Received data

(f) Algorithm 2 restored image

(g) True image

(h) Received data

(i) Algorithm 2 restored image

(j) True image

(k) Received data

(l) Algorithm 2 restored image

Figure 4: Experimental results on Images 2–4 with Blur 1. Row 1, l-r: Image 2, received data corrupted by Blur 1, restored image using Algorithm 2. The PSNR/SNR is increased by 4.53/5.01 from 20.76/6.58 to 25.29/11.59. Row 2, l-r: Image 2, received data corrupted by Blur 2, restored image using Algorithm 2. The PSNR/SNR is increased by 15.26/30.52 from 12.61/−13.85 to 27.87/16.67. Our model is capable of restoring details in both cases and of preserving the background in black. Row 3, l-r: Image 3, received data corrupted by Blur 1, restored image using Algorithm 2. The PSNR/SNR is increased by 11.5/11.66 from 24.26/18.5 to 35.76/30.16. Row 4, l-r: Image 4, received data corrupted by Blur 1, restored image using Algorithm 2. The PSNR/SNR is increased by 11.76/11.97 from 19.85/15.72 to 31.61/27.69. Our model is capable of restoring many detailed features and some fine details as well as sharpening edges. There are very few defects in the restored image, notably surrounding the rope in (l).

Set 3 (detailed photograph images with blurs). Result Set 3 consists of Image 3 and Image 4 corrupted by out of focus blur (Blur 1) or Gaussian blur (Blur 2). We see in Figure 4 that in the case of out of focus blur our model is able to sharpen the images and recover many detailed features including fine details but can introduce pattern defects. In the case of Gaussian blur, we see in Figure 6 that many features are recovered in the image and background objects can be distinguished. The intensity ranges are also preserved.

Set 4 (detailed retinal images with blurs). Result Set 4 consists of Image 5 and Image 6 corrupted by out of focus blur (Blur 1) or Gaussian blur (Blur 2). In Figure 5, we see that in the case of out of focus blur our model is able to sharpen the images and recover many fine details. In the case of Gaussian blur, we see in Figure 7 that many features are recovered in the image, including blood vessels which were unseen by the blur. The intensity ranges are also preserved.

(a) True image

(b) Received data

(c) Result by Algorithm 1

(d) Result by Algorithm 2

(e) True image

(f) Received data

(g) Restored image

Figure 5: Experimental results on Images 5-6 with Blur 1. Row 1, l-r for Image 5: the PSNR/SNR is lowered from 27.48/15.36 by 0.92/0.97 to 26.56/14.39 using Algorithm 1 but increased by 5.93/6.27 to 33.41/21.63 using Algorithm 2. Row 2, l-r for Image 6: the PSNR/SNR is increased by 11.6/11.96 from 20.19/12.11 to 31.79/24.07. Our model is capable of restoring many detailed features and sharpening edges. Several of the blood vessels become visible in (d) and very fine details can be distinguished in (g).

(a) True image

(b) Received data

(c) Algorithm 2 restored image

(d) True image

(e) Received data

(f) Algorithm 2 restored image

Figure 6: Experimental results on Images 3-4 with Blur 2, restored using Algorithm 2. Row 1, l-r for Image 3: the PSNR/SNR is increased by 16.77/5.71 from 10.07/15.5 to 26.84/21.21. Row 2, l-r for Image 4: the PSNR/SNR is increased by 4.36/4.67 from 17.79/13.5 to 22.15/18.17. In this more challenging case of Gaussian blur, our model is capable of restoring some detailed features, including the books in the background of (c) and the buildings in (f).

(a) True image

(b) Received data

(c) Algorithm 2 PSNR = 23.4

(d) Algorithm 3 PSNR = 20.8

(e) True image

(f) Received data

(g) Algorithm 2 PSNR = 21.6

Figure 7: Test results on Images 5-6 with Blur 2. Row 1, l-r for Image 5, received data , restored image using Algorithm 2 and by Algorithm 3. The PSNR/SNR is increased from 7.65/−12.97 by 15.7/30.42 to 23.35/17.45 using Algorithm 2 and by 13.12/27.75 to 20.77/14.78 using Algorithm 3. Row 2, l-r: Image 6, received data corrupted by Blur 2, restored image using Algorithm 2. The PSNR/SNR is increased by 3.37/3.77 from 18.18/9.91 to 21.55/13.68. Our model is capable of restoring detailed features in these challenging cases of Gaussian blur. Most of the details are restored in both cases. The slower Algorithm 2 performs better than the faster Algorithm 3.

4.2. Comparison of Algorithms 2 and 3

Since our primary aim of this paper is to illustrate what is achievable for a single channel blind deconvolution, we shall mainly focus on quality of restoration, rather than algorithms’ efficiency (which is also important). In fact, in our presented implementations, Algorithm 3 is a few times faster than Algorithm 2 but the quality comparison says the opposite. This observation is in line with findings of [43] where a faster 2-level augmented Lagrangian method (ALM) with more approximations produces less quality of restoration for denoising compared to the slower 1-level ALM.

To quantify the restorations, it is useful to use the well known measures of signal-to-noise ratio (SNR) and the peak signal-to-noise ratio (PSNR), defined, respectively, as

where and are the true image and the restored image, respectively, and is the mean value of the original image.

Finally we are ready to show in Table 1 some quantitative measures of the restorations using 3 test images with two kinds of blurs. Clearly Algorithm 2 is better than Algorithm 3, though the latter is faster.

Table 1: Comparison of the new Algorithms 2 and 3. Clearly the slower Algorithm 2 gives better quality than Algorithm 3 (other CPU tests show that the latter is up to 6 times faster).

Test image
Blur type
Received
Algorithm 2
Algorithm 3

PSNR
SNR
PSNR
SNR
PSNR
SNR

1
Out of focus
20.15
13.54
28.36
21.94
27.54
21.2

1
Gauss
7.88
−13.4
27.32
21.13
24.08
17.67

2
Out of focus
20.76
6.58
25.29
11.59
22.09
10.32

2
Gauss
12.61
−13.85
27.87
16.67
25.83
14.22

5
Out of focus
27.48
15.36
33.41
21.63
32.21
20.28

5
Gauss
7.65
−12.97
23.35
17.45
20.77
14.78

5. Conclusions and Future Work

We have presented a total variation based blind deconvolution model with solution positivity achieved by implicit transforms and two solution algorithms for reconstructing a deblurred image along with its blur kernel. We demonstrated that we can ensure positivity and keep the correct range of the image intensities in the case of several blur types, extending the original Chan-Wong model’s applicability. This model is particularly effective in reconstructing the kernel without significant defects which can significantly impair the results of previous blind deconvolution algorithms.

Further work involves integrating the remaining constraints into the functional and automatic selection of regularisation parameters. While there has been work in the parameter selection with nonblind imaging models [51–53], further work is required to develop a method for the selection of optimal parameters for both regularisation terms in the blind model.

Appendices

A. Parameter Selection for Nonnegativity

While the transform necessitates the selection of additional parameters , it should be noted that they are easily chosen. In this section, we consider each pair in turn.

The parameters are easily chosen, assuming knowledge of the bits-per-sample (bps) value of the true image and for and setting for the blur function since we assume it to have a unit integral. Better results may be achieved by selecting a lower value for but if it is set too low, then recovery of the true kernel will not be possible. The small parameters and should be chosen to be proportional to and , respectively. Typically, is an appropriate value for the image and is used for the kernel since typical kernels have a small maximum value.

Of the remaining parameters, and determine the value of at .

and attempt to keep close to and close to . To do this for the image, we define two points and where . Then we can calculate

For our model, we fix and . Substitute this into (A.1), we obtain an equation for selecting given in (A.2) which is dependent on parameters which have been discussed earlier and , which is typically set to :

The parameters and may be selected to control the value of at and at . We consider two cases for the image; similar cases apply for the kernel. The first case is given by and the second given by at where is the lower bound of . The first option will make the graph pass through zero at the midpoint of the intensity values and the second will make all values of naturally positive since the lower bound of will be equal to zero. We attempt to satisfy each of these cases by the selection of only. Letting , we have and so we have

to give the first and second cases, respectively.

In application, either of these will be sufficient to recover the image with similar results. In the case of the kernel, better results are obtained with . It is therefoire advised that and are the appropriate values for these parameters.

In summary, once and are defined, the remaining parameters in the transforms and can be determined automatically.

B. Derivation of Euler Lagrange Equations for Model (19)

Considering the minimisation, we do it in parts. Note that when we minimise with respect to one part, for example, when we minimise with respect to the image, the lagged component of the transform is assumed to be equal to the unlagged part; that is . We may therefore write for the minimisations with respect to and , respectively:

For the second term, we minimise with respect to as follows. We wish to calculate the minimisation of the TV seminorm that is .

Letting , we have

Now noting that and letting be defined as follows, we have

Using Green’s theorem we have

Therefore the minimisation of the TV term alone leads to

with boundary conditions where we used and denoted

Thus the above term is simplified to

Combining with the fitting term, we obtain the Euler Lagrange equations:

where

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1–4, pp. 259–268, 1992.

A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image denoising,” in IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '05), vol. 2020, pp. 60–65, June 2005.

T. F. Chan and C.-K. Wong, “Total variation blind deconvolution,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 370–375, 1998.

J. P. Oliveira, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Adaptive total variation image deblurring: a majorization-minimization approach,” Signal Processing, vol. 89, no. 9, pp. 1683–1693, 2009.

T. F. Chan and J. H. Shen, “Variational image inpainting,” Communications on Pure and Applied Mathematics, vol. 58, no. 5, pp. 579–619, 2005.

Y. M. Huang, M. K. Ng, and Y.-W. Wen, “A fast total variation minimization method for image restoration,” Multiscale Modeling and Simulation, vol. 7, no. 2, pp. 774–795, 2008.

R. H. Chan, M. Tao, and X. Yuan, “Constrained total variation deblurring models and fast algorithms based on alternating direction method of multipliers,” SIAM Journal on Imaging Sciences, vol. 6, no. 1, pp. 680–697, 2013.

Y. Shi, Q. Chang, and J. Xu, “Convergence of fixed point iteration for deblurring and denoising problem,” Applied Mathematics and Computation, vol. 189, no. 2, pp. 1178–1185, 2007.

P. C. Hansen, J. G. Nagy, and D. P. O'Leary, Deblurring Images: Matrices, Spectra, and Filtering, SIAM, Philadelphia, Pa, USA, 2007.

M. I. Sezan and A. M. Tekalp, “Survey of recent developments in digital image restoration,” Optical Engineering, vol. 29, no. 5, pp. 393–404, 1990.

L. Rada and K. Chen, “A new variational model with dual level set functions for selective segmentation,” Communications in Computational Physics, vol. 12, no. 1, pp. 261–283, 2012.

J. Modersitzkie, Numerical Methods for Image Registration, Oxford University Press, Oxford, UK, 2004.

P. Thévenaz, U. E. Ruttimann, and M. Unser, “A pyramid approach to subpixel registration based on intensity,” IEEE Transactions on Image Processing, vol. 7, no. 1, pp. 27–41, 1998.

C. R. Vogel and M. E. Oman, “Fast, robust total variation-based reconstruction of noisy, blurred images,” IEEE Transactions on Image Processing, vol. 7, no. 6, pp. 813–824, 1998.

C. R. Vogel, Computational Methods for Inverse Problems, SIAM, 2002.

J. O. Abad, S. Morigi, L. Reichel, and F. Sgallari, “Alternating Krylov subspace image restoration methods,” Journal of Computational and Applied Mathematics, vol. 236, no. 8, pp. 2049–2062, 2012.

T. F. Chan and K. Chen, “An optimization-based multilevel algorithm for total variation image denoising,” Multiscale Modeling and Simulation, vol. 5, no. 2, pp. 615–645, 2006.

Y.-L. You and M. Kaveh, “A regularization approach to joint blur identification and image restoration,” IEEE Transactions on Image Processing, vol. 5, no. 3, pp. 416–428, 1996.

J. H. Money and S. H. Kang, “Total variation semi-blind deconvolution using shock filters,” Image and Vision Computing, vol. 26, pp. 302–314, 2008.

F. Šroubek and P. Milanfar, “Robust multichannel blind deconvolution via fast alternating minimization,” IEEE Transactions on Image Processing, vol. 21, no. 4, pp. 1687–1700, 2012.

F. Sroubek and J. Flusser, “Multichannel blind iterative image restoration,” IEEE Transactions on Image Processing, vol. 12, no. 9, pp. 1094–1106, 2003.

A. Rav-Acha and S. Peleg, “Two motion-blurred images are better than one,” Pattern Recognition Letters, vol. 26, no. 3, pp. 311–317, 2005.

A. S. Carasso, “Direct blind deconvolution,” SIAM Journal on Applied Mathematics, vol. 61, no. 6, pp. 1980–2007, 2001.

W. Wang and M. K. Ng, “On algorithms for automatic deblurring from a single image,” Journal of Computational Mathematics, vol. 30, no. 1, pp. 80–100, 2012.

N. Miura, “Blind deconvolution under band limitation,” Optics Letters, vol. 28, no. 23, pp. 2312–2314, 2003.

Y. Biraud, “A new approach for increasing the resolving power by data processing,” Astronomy & Astrophysics, vol. 1, pp. 124–127, 1969.

M. S. C. Almeida and L. B. Almeida, “Blind and semi-blind deblurring of natural images,” IEEE Transactions on Image Processing, vol. 19, no. 1, pp. 36–52, 2010.

L. Bar, N. Sochen, and N. Kiryati, “Semi-blind image restoration via Mumford-Shah regularization,” IEEE Transactions on Image Processing, vol. 15, no. 2, pp. 483–493, 2006.

A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding and evaluating blind deconvolution algorithms,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '09), pp. 1964–1971, June 2009.

A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Efficient marginal likelihood optimization in blind deconvolution,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR '11), pp. 2657–2664, June 2011.

A. S. Carasso, “Linear and nonlinear image deblurring: a documented study,” SIAM Journal on Numerical Analysis, vol. 36, no. 6, pp. 1659–1689, 1999.

J. M. Bardsley and C. R. Vogel, “A nonnegatively constrained convex programming method for image reconstruction,” SIAM Journal on Scientific Computing, vol. 25, no. 4, pp. 1326–1343, 2003.

F. Benvenuto, R. Zanella, L. Zanni, and M. Bertero, “Nonnegative least-squares image deblurring: improved gradient projection approaches,” Inverse Problems, vol. 26, no. 2, Article ID 025004, 2010.

M. Jiang, G. Wang, M. W. Skinner, J. T. Rubinstein, and M. W. Vannier, “Blind deblurring of spiral CT images,” IEEE Transactions on Medical Imaging, vol. 22, no. 7, pp. 837–845, 2003.

D. L. Snyder, J. A. O'Sullivan, B. R. Whiting et al., “Deblurring subject to nonnegativity constraints when known functions are present with application to object-constrained computerized tomography,” IEEE Transactions on Medical Imaging, vol. 20, no. 10, pp. 1009–1017, 2001.

A. Cornelio, E. Loli Piccolomini, and J. G. Nagy, “Constrained numerical optimization methods for blind deconvolution,” Numerical Algorithms, vol. 65, no. 1, pp. 23–42, 2014.

H. Xu, Q. Sun, N. Luo, G. Cao, and D. Xia, “Iterative nonlocal totalvariation regularization method for image restoration,” PLoS ONE, vol. 8, no. 6, Article ID e65865, 2013.

X. Zhang, M. Burger, X. Bresson, and S. Osher, “Bregmanized nonlocal regularization for deconvolution and sparse reconstruction,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 253–276, 2010.

F. Knoll, K. Bredies, T. Pock, and R. Stollberger, “Second order total generalized variation (TGV) for MRI,” Magnetic Resonance in Medicine, vol. 65, no. 2, pp. 480–491, 2011.

K. Bredies and M. Holler, “A TGV regularized wavelet based zooming model,” in Scale Spaceand Variational Methods in Computer Vision, vol. 7893 of Lecture Notes in Computer Science, pp. 149–160, Springer, Berlin, Germany, 2013.

K. Bredies, Y. Dong, and M. Hintermüller, “Spatially dependent regularization parameter selection in total generalized variation models for image restoration,” International Journal of Computer Mathematics, vol. 90, no. 1, pp. 109–123, 2013.

C. Brito-Loeza and K. Chen, “Multigrid algorithm for high order denoising,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 363–389, 2010.

L. Sun and K. Chen, “A new iterative algorithm for mean curvature-based variational image denoising,” BIT Numerical Mathematics, vol. 54, pp. 523–553, 2014.

R. Kaftory, N. Sochen, and Y. Y. Zeevi, “Variational blind deconvolution of multi-channel images,” International Journal of Imaging Systems and Technology, vol. 15, no. 1, pp. 56–63, 2005.

M. Lysaker, A. Lundervold, and X.-C. Tai, “Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time,” IEEE Transactions on Image Processing, vol. 12, no. 12, pp. 1579–1590, 2003.

C. R. Vogel and M. E. Oman, “Fast total variation-based image reconstruction,” in Proceedings of the ASME Design Engineering Conferences, vol. 3, part C, pp. 1009–1015, 1995.

K. Chen, Matrix Preconditioning Techniques and Applications, Cambridge University Press, Cambridge, UK, 2005.

M. K. Ng, R. H. Chan, and W.-C. Tang, “A fast algorithm for deblurring models with Neumann boundary conditions,” SIAM Journal on Scientific Computing, vol. 21, no. 3, pp. 851–866, 1999.

M. K. Ng and N. K. Bose, “Mathematical analysis of super-resolution methodology,” IEEE Signal Processing Magazine, vol. 20, no. 3, pp. 62–74, 2003.

F. Wang, Alternating direction methods for image recovery [Ph.D. thesis], Hong Kong Baptist University, 2012.

P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Review, vol. 34, no. 4, pp. 561–580, 1992.

J. Zhang, K. Chen, and B. Yu, “An iterative lagrange multiplier method for constrained total-variation-based image denoising,” SIAM Journal on Numerical Analysis, vol. 50, no. 3, pp. 983–1003, 2012.

K. Chen, E. L. Piccolomini, and F. Zama, “An automatic regularization parameter selection algorithm in the total variation model for image deblurring,” Numerical Algorithms, vol. 67, no. 1, pp. 73–92, 2013.