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.
| ||||||||||||||||||||
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
(a) Ex2 true image
(b) Ex2 received data
(c) Failed restoration using Algorithm 1
(d) Restored with
(e) Restored with
(f) Restored with
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.
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.
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](https://static.hindawi.com/articles/ijcm/volume-2015/860263/figures/860263.fig.003c.jpg)
(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
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
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
(a) True image
(b) Received data
(c) Algorithm 2 restored image
(d) True image
(e) Received data
(f) Algorithm 2 restored image
(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
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.
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.