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

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.


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] 2 International Journal of Computational Mathematics 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][30][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.

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 + = max{ , 0}. The same projection idea can be applied to impose a box constraint to ensure low ≤ ≤ up . See [7,32].
There have been several other related ideas for enforcing nonnegativity in image processing [32][33][34][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.

The Earlier Blind Deconvolution Work. You and Kaveh
where the fitting term is a common choice for (1) and there is freedom to choose the regularisation terms 1 and 2 ; You and Kaveh used the 1 seminorm for these two terms. Chan and Wong [3] proposed an improvement to this using the total variation (TV) seminorm for regularisation given by 1 ( ) = ∫ Ω |∇ | Ω and 2 (ℎ) = ∫ Ω |∇ℎ| Ω, hence solving min ,ℎ where |∇ | = √ 2 + 2 . Minimising (4) with respect to the image and the kernel ℎ, we obtain the coupled partial differential equations given by International Journal of Computational Mathematics
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, ℎ( , ) > 0 and ( , ) > 0, the kernel should be symmetric (ℎ( , ) = ℎ(− , − )), and the kernel should have a unit integral ∫ Ω ℎ( , ) Ω = 1. These constraints are imposed exactly but only after each alternate minimisation step. The complete algorithm is given in Algorithm 1.
Adding the above 4 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 (4) 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]. 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. [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 ℎ( , ) = ( ( , )) 2 , ( , ) = ( ( , )) 2 . Then after Fourier transforms, one getŝ

The Blind Deconvolution Model by Miura
where the summations imply formulations after discretization [25]. Further a conjugate gradient type solver is utilised to computê,̂which will be used to yield the nonnegative solutions ℎ, . [19].

A Shock Filter Based Model by Money and Kang
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 which is a decoupled system and can be solved directly in a noniterative way between ℎ and . [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 1 , . . . , of the same object are available which is readily possible for video images in some situations. The minimization proposed is

A Multichannel Blind Deconvolution Model byŠroubek and Milanfar
International Journal of Computational Mathematics 5 where in a discrete setting with R Δ 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).

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.

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 = exp(− ( , )) for some function . Unfortunately this particular transform does not work as it is not capable of dealing with dark regions (where ≈ 0) in a stable way. Another choice could be = ( ) = ( ( , )) 2 as in [25,26]. This is valid for ≈ 0 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 , , > 0, max ≤ , and is a small tuning parameter which controls the spread of the function. Clearly = 0 poses no problems to the transform. For the usual intensity range [0, 255] for , we may take = 255, = 1/10, = 1/100: Since ( ) for ∈ R and −1 ( ) for ∈ (− , ∞) are monotone functions, we can work out their lower and upper bounds: Note, for (13)

Reformulation of the Blind Deblurring Model.
In order to apply the same transform to both the image and the kernel ℎ, we introduce the 8 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, 1 6 International Journal of Computational Mathematics Similar to the previous section, we need not vary the parameters to take into account varying noise or blur levels with the exception of 1 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 = a ( ) and ℎ = b ( ), respectively. Here from solving (17), the nonnegativity constraints are exactly and implicitly enforced; that is , ℎ ≥ 0, but the remaining symmetry and unit integral constraints on the kernel are still required. The advantage of realizing positivity , ℎ > 0 (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.
Minimising the above functionals from (19) with respect to and , giveñand̃, yields the Euler Lagrange equations (see Appendix B for the derivation): After estimating (0) and ℎ (0) of the image and the kernel, respectively, we apply the inverse transforms obtaining (0) and (0) . Then on convergence, , ℎ are finally obtained from , .
, then exit or continue. (10) end for (11) Accept the restored image = and the restored kernel ℎ = .
Discretisation of the Fitting (Integral) Term. We wish to discretise the quantity b (̃) * 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].
where 1 , 2 are of BTTB form and and denotes the discretised TV operator for (22).

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 1 and 2 . 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 1 , 2 are positive constants.

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 = a ( ), ℎ = b ( ). 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.
Minimising with respect to each of the arguments, we have = ℎ (ℎ − ) + 1 + 1 ( − ( 1 +̃1)) = 0, where 1 , 2 ,̃1 = 1 (̃),̃2 = 2 (̃) 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]).

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 2 norm to regularise the blurring kernel ℎ = b ( ) and the TV to regularise the image = a ( ), as a hybrid model of (3) and (4): Then the Euler-Lagrange equations corresponding to Method I of Section 3.3 can be derived 2 2 (V 2 ( ,̃)) (V 2 ( ,̃) * − 2 (̃, ,̃)) Other computational details and use of Method II can be considered similarly.

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 6 test images in 4 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 10 −3 with 10 3 maximum iterations although this number is rarely reached.

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.
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.
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.

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 mean 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.

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][52][53], further work is required to develop a method for the selection of optimal parameters for both regularisation terms in the blind model.

A. Parameter Selection for Nonnegativity
While the transform necessitates the selection of additional parameters 1 , . . . , 4 , 1 , . . . , 4 , it should be noted that they are easily chosen. In this section, we consider each pair , in turn. The parameters 1 , 1 are easily chosen, assuming knowledge of the bits-per-sample (bps) value of the true image and for 1 and setting 1 = 1 for the blur function since we assume it to have a unit integral. Better results may be achieved by selecting a lower value for 1 but if it is set too low, then recovery of the true kernel will not be possible. The small parameters 4 and 4 should be chosen to be proportional to 1 and 1 , respectively. Typically, 4 = 1 /255 is an appropriate value for the image and 4 = 10 −4 1 is used for the kernel since typical kernels have a small maximum value.
Of the remaining parameters, 2 and 2 determine the value of at = 0.
The parameters 2 and 2 may be selected to control the value of at = a ( ) = 0 and at ℎ = b ( ) = 0. We consider two cases for the image; similar cases apply for the kernel. The first case is given by a ( ) = 1 /2 and the second given by a ( ) = 1 at = 0 where 1 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 Figure 7: Test results on Images 5-6 with Blur 2. Row 1, l-r for Image 5, received data , restored image 2 using Algorithm 2 and 3 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. 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 2 = 1. It is therefoire advised that 2 = 1 and 2 = 1 are the appropriate values for these parameters.
In summary, once 1 , 1 and 4 , 4 are defined, the remaining parameters in the transforms a ( ) and b ( ) 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 TV( a ( )) = ∫ Ω |∇ a ( )| Ω that is ( / )(TV( a ( + )))| =0 . Letting 1 = 1 ( ) = 2 −2 / 3 , we have minTV ( a ( )) = ∫ Ω 2 ( 1 + 2 4 ) 1 √ 2 + 2 (1 + 1 )