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

Authors: Ke Chen; Simon P. Harding; Bryan M. Williams; Yalin Zheng
Journal: International Journal of Computational Mathematics
 (2015)

---

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


---

## Body

## 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 functionz=z(x,y) in domain Ω has been contaminated with additive noise η and convolution (blur) operator h(x,y):
(1)z=h*u+η,
where η is an unknown Gaussian white noise with zero mean and u is the image to be restored. When h 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 h is unknown, simultaneous restoration of u,h is of great interest.A model for blind deblurring by recovering both the kernel and the imagesimultaneously 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 u and the kernel h, 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 h 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. Section2 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 nonnegativityu+=max⁡{u,0}. The same projection idea can be applied to impose a box constraint to ensure ulow≤u≤uup. 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 h 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 u(x)=(E(x))2, with x∈[0,1], in restoring one dimensional signal u from the model z(x)=h(x)*u(x)+η(x), where h(x) is theknown blur function and η the noise, or z^(s)=h^(s)u^(s)+η^(s) after Fourier transform. The central idea here is that any E(x) or its Fourier transform E^ can lead to nonnegative restoration u. For s∈[0,s0] with some cut-off frequency s0, h^(s)≠0 so u^(s)=z^(s)/h^(s)-η^(s)/h^(s). Noting that u(x)=(E(x))2 leads to u^(s)=E^(s)*E^(s), the method of Biraud [26] is
(2)u^(s)=argminf1s0∫0s0‍f*f-z^(s)h^(s)2ds.
To solve (2), a parametric iterative approach is proposed by fn=fn-1+ω-n for n≥1. See [26]. Once a good approximation E^(s) is obtained, an inverse transform would yield E(x) and then a nonnegative restoration u(x).


### 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
(3)min⁡u,hfYKu,h=12h*u-zL2Ω2+α1J1u+α2J2h,
where the fitting term is a common choice for (1) and there is freedom to choose the regularisation terms J1 and J2; You and Kaveh used the H1 seminorm for these two terms.Chan and Wong [3] proposed an improvement to this using the total variation (TV) seminorm for regularisation given by J1(u)=∫Ω‍|∇u|dΩ and J2(h)=∫Ω‍|∇h|dΩ, hence solving
(4)min⁡u,hfCWu,h=12h*u-zL2(Ω)2+α1∫Ω‍∇udΩ+α2∫Ω‍∇hdΩ,
where |∇u|=ux2+uy2. Minimising (4) with respect to the image u and the kernel h, we obtain the coupled partial differential equations given by
(5)h(-x,-y)*((h*u)(x,y)-z(x,y))-α1∇·∇u(x,y)∇u(x,y)β=0,u(-x,-y)*((u*h)(x,y)-z(x,y))-α2∇·∇h(x,y)∇h(x,y)β=0,
where |∇u(x,y)|β=ux2(x,y)+uy2(x,y)+β2 (similarly for |∇h(x,y)|β also) with β a small positive parameter introduced to avoid division by zero, which is set to 10-3/2 for experimental results for both the image u and blur funcion h 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 L2 regularisation, given by J(u)=∫Ω‍|∇u|2dΩ, 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 kernelh and the image u, 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, h(x,y)>0 and u(x,y)>0, the kernel should be symmetric (h(x,y)=h(-x,-y)), and the kernel should have a unit integral ∫Ω‍h(x,y)dΩ=1. These constraints are imposed exactly but only after each alternate minimisation step. The complete algorithm is given in Algorithm 1.Algorithm 1:Chan-Wong [3] method for the solution of (4).
(1)function  Chan-Wong  (h(0),u(0),z,α1,α2,max⁡it)
(2)for  k←1tomax⁡it  do
(3)    h(k+(1/4))←SOLVEuk-x,-y*uk*hk-z-α2∇⋅∇h(k)|∇h(k)|=0
(4)    h(k+(1/2))←hk+(1/4)x,yifh(k+(1/4))(x,y)≥00otherwise
(5)    h(k+(3/4))←12h(k+(1/2))(x,y)+h(k+(1/2))(-x,-y)
(6)    h(k+1)←h(k+(3/4))/∑x‍∑y‍h(k+(3/4))(x,y)
(7)    u(k+(1/2))←SOLVEh(k+1)(-x,-y)*(h(k+1)*u(k)-z)-α1∇⋅∇u(k)|∇u(k)|=0
(8)    u(k+1)←u(k+(1/2))(x,y)ifu(k+(1/2))(x,y)≥00otherwise
(9) end for
(10)end functionAdding the above4 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:
(6)h(k+1/2)⟵h(k+1/4)(x,y),ifh(k+1/4)(x,y)≥κ,0,otherwise.
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].Good restoration results for Example 1 (a) from a corrupted image (b) using Algorithm1. This model is able to improve the edges of the restored image (c), though the restoration is not excellent.
(a)
True data
(b)
Blurred image
(c)
Algorithm1 resultIllustration of the failure of Algorithm1 for a retinal image (Example 2). (a) True image; (b) corrupted image by out of focus blur; (c) failed restoration u; (d) restored u with thresholding κ=10-2; (e) restored u with thresholding κ=1/3; (f) restored h with thresholding κ=1/3.
(a)
Ex2 true image
(b)
Ex2 received data
(c)
Failed restoration using Algorithm1
(d)
Restoredu with κ=10-2
(e)
Restoredu with κ=1/3
(f)
Restoredh with κ=1/3Our 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 u and h. Starting from (1), that is, z(x,y)=h(x,y)*u(x,y)+η, with both h and u as unknowns, he defined h(x,y)=(ϕ(x,y))2,u(x,y)=(ψ(x,y))2. Then after Fourier transforms, one gets
(7)z^(s,t)={ϕ^(s,t)*ϕ^(s,t)}{ψ^(s,t)*ψ^(s,t)}+η^(s,t).
Further similar to (2), it is proposed to solve
(8)min⁡ϕ^,ψ^∑ϕ^‍∑ψ^‍{ϕ^(s,t)*ϕ^(s,t)}{ψ^(s,t)*ψ^(s,t)}-z^(s,t)2,
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 h,u.


### 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 u in the second equation of (5) is replaced by a reference image uτ which is obtained by using a shock filter to capture image edges in the blurred z. Then (5) becomes
(9)uτ-x,-y*uτ*hx,y-zx,y-α2∇·∇hx,y∇h(x,y)β=0,h-x,-y*h*ux,y-zx,y-α1∇·∇ux,y∇u(x,y)β=0,
which is a decoupled system and can be solved directly in a noniterative way between h and u.


### 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 u if K blurred images z1,…,zK of the same object are available which is readily possible for video images in some situations. The minimization proposed is
(10)min⁡u,{hk}γ2∑k=1K‍hk*u-zk2+Qu+Rhk,
where in a discrete setting
(11)Rh=δ2hTRΔh+Ψ(h),Ψ(h)=∑‍ψ(hk(i)),ψ(t)=t,ift≥0,+∞,otherwise
with RΔ a Laplacian related operator, and γ is a parameter. Here Q denotes the total variation regularizer for u and the crucial choice of ψ(t) ensures positivity of hk. However the same treatment was not applied to u. The optimization problem was further solved by a splitting idea in an augmented Lagrangian method (ALM).


## 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
(3)min⁡u,hfYKu,h=12h*u-zL2Ω2+α1J1u+α2J2h,
where the fitting term is a common choice for (1) and there is freedom to choose the regularisation terms J1 and J2; You and Kaveh used the H1 seminorm for these two terms.Chan and Wong [3] proposed an improvement to this using the total variation (TV) seminorm for regularisation given by J1(u)=∫Ω‍|∇u|dΩ and J2(h)=∫Ω‍|∇h|dΩ, hence solving
(4)min⁡u,hfCWu,h=12h*u-zL2(Ω)2+α1∫Ω‍∇udΩ+α2∫Ω‍∇hdΩ,
where |∇u|=ux2+uy2. Minimising (4) with respect to the image u and the kernel h, we obtain the coupled partial differential equations given by
(5)h(-x,-y)*((h*u)(x,y)-z(x,y))-α1∇·∇u(x,y)∇u(x,y)β=0,u(-x,-y)*((u*h)(x,y)-z(x,y))-α2∇·∇h(x,y)∇h(x,y)β=0,
where |∇u(x,y)|β=ux2(x,y)+uy2(x,y)+β2 (similarly for |∇h(x,y)|β also) with β a small positive parameter introduced to avoid division by zero, which is set to 10-3/2 for experimental results for both the image u and blur funcion h 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 L2 regularisation, given by J(u)=∫Ω‍|∇u|2dΩ, 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 kernelh and the image u, 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, h(x,y)>0 and u(x,y)>0, the kernel should be symmetric (h(x,y)=h(-x,-y)), and the kernel should have a unit integral ∫Ω‍h(x,y)dΩ=1. These constraints are imposed exactly but only after each alternate minimisation step. The complete algorithm is given in Algorithm 1.Algorithm 1:Chan-Wong [3] method for the solution of (4).
(1)function  Chan-Wong  (h(0),u(0),z,α1,α2,max⁡it)
(2)for  k←1tomax⁡it  do
(3)    h(k+(1/4))←SOLVEuk-x,-y*uk*hk-z-α2∇⋅∇h(k)|∇h(k)|=0
(4)    h(k+(1/2))←hk+(1/4)x,yifh(k+(1/4))(x,y)≥00otherwise
(5)    h(k+(3/4))←12h(k+(1/2))(x,y)+h(k+(1/2))(-x,-y)
(6)    h(k+1)←h(k+(3/4))/∑x‍∑y‍h(k+(3/4))(x,y)
(7)    u(k+(1/2))←SOLVEh(k+1)(-x,-y)*(h(k+1)*u(k)-z)-α1∇⋅∇u(k)|∇u(k)|=0
(8)    u(k+1)←u(k+(1/2))(x,y)ifu(k+(1/2))(x,y)≥00otherwise
(9) end for
(10)end functionAdding the above4 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:
(6)h(k+1/2)⟵h(k+1/4)(x,y),ifh(k+1/4)(x,y)≥κ,0,otherwise.
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].Good restoration results for Example 1 (a) from a corrupted image (b) using Algorithm1. This model is able to improve the edges of the restored image (c), though the restoration is not excellent.
(a)
True data
(b)
Blurred image
(c)
Algorithm1 resultIllustration of the failure of Algorithm1 for a retinal image (Example 2). (a) True image; (b) corrupted image by out of focus blur; (c) failed restoration u; (d) restored u with thresholding κ=10-2; (e) restored u with thresholding κ=1/3; (f) restored h with thresholding κ=1/3.
(a)
Ex2 true image
(b)
Ex2 received data
(c)
Failed restoration using Algorithm1
(d)
Restoredu with κ=10-2
(e)
Restoredu with κ=1/3
(f)
Restoredh with κ=1/3Our 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 u and h. Starting from (1), that is, z(x,y)=h(x,y)*u(x,y)+η, with both h and u as unknowns, he defined h(x,y)=(ϕ(x,y))2,u(x,y)=(ψ(x,y))2. Then after Fourier transforms, one gets
(7)z^(s,t)={ϕ^(s,t)*ϕ^(s,t)}{ψ^(s,t)*ψ^(s,t)}+η^(s,t).
Further similar to (2), it is proposed to solve
(8)min⁡ϕ^,ψ^∑ϕ^‍∑ψ^‍{ϕ^(s,t)*ϕ^(s,t)}{ψ^(s,t)*ψ^(s,t)}-z^(s,t)2,
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 h,u.


## 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 u in the second equation of (5) is replaced by a reference image uτ which is obtained by using a shock filter to capture image edges in the blurred z. Then (5) becomes
(9)uτ-x,-y*uτ*hx,y-zx,y-α2∇·∇hx,y∇h(x,y)β=0,h-x,-y*h*ux,y-zx,y-α1∇·∇ux,y∇u(x,y)β=0,
which is a decoupled system and can be solved directly in a noniterative way between h and u.


## 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 u if K blurred images z1,…,zK of the same object are available which is readily possible for video images in some situations. The minimization proposed is
(10)min⁡u,{hk}γ2∑k=1K‍hk*u-zk2+Qu+Rhk,
where in a discrete setting
(11)Rh=δ2hTRΔh+Ψ(h),Ψ(h)=∑‍ψ(hk(i)),ψ(t)=t,ift≥0,+∞,otherwise
with RΔ a Laplacian related operator, and γ is a parameter. Here Q denotes the total variation regularizer for u and the crucial choice of ψ(t) ensures positivity of hk. However the same treatment was not applied to u. 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 u,h directly instead of solving for u^,h^ 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 isu=exp⁡(-ψ(x,y)) for some function ψ. Unfortunately this particular transform does not work as it is not capable of dealing with dark regions (where u≈0) in a stable way. Another choice could be u=T(ψ)=(ψ(x,y))2 as in [25, 26]. This is valid for u≈0 but this T(ψ) 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(12)u=Tψ=w+2a1+e-2ψ/ϵ-a,withiiiiiiψ=T-1u=ϵ2ln⁡a+ua+w-u,
where constants w,a,ϵ>0, max⁡u≤w, and ϵ is a small tuning parameter which controls the spread of the function. Clearly u=0 poses no problems to the transform. For the usual intensity range [0,255] for u, we may take w=255, a=1/10, ϵ=1/100:
(13)u=T(ψ)=255+1/51+e-200ψ-110,withψ=T-1(u)=1200ln⁡u+1/10255-u+1/10.
Since T(ψ) for ψ∈R and T-1(u) for u∈(-a,∞) are monotone functions, we can work out their lower and upper bounds:
(14)u=T(ψ)∈-a,w+aforψ∈R;u=T(ψ)∈[0,w+a),ifψ≥ϵ2ln⁡aa+w;ψ=T-1u∈ϵ2ln⁡aa+w,ϵ2ln⁡a+waifu∈0,w⊂-a,w+a;ψ=T-1u∈ϵ2ln⁡3aa+2w,ϵ2ln⁡a+2w3aifu∈-a2,w+a2⊂-a,w+a.
Note, for (13), (ϵ/2)ln⁡⁡((a+w)/a)=ln⁡(2551)/200≈0.039 and (ϵ/2)ln⁡((a+2w)/3a)=ln⁡(5103)/200≈0.043. 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 imageu and the kernel h, we introduce the 8 parameters with subscripts as follows:
(15)u=Ta(ψ)=a1+2a41+a2e-2ψ/a3-a4,h=Tb(ω)=b1+2b41+b2e-2ω/b3-b4,
for the image and kernel, respectively; here all constants can be fixed before proceeding (see Appendix A for more details). In particular, a1 and b1 are the expected upper limits of the image intensity values and kernel values; a2 and b2 are introduced to control the values of the image and kernel at ψ=0 and ω=0, and a3 and b3 control the spread of ψ and ω. To give one feasible set, for image u, we have a=(a1,a2,a3,a4)=(255,1,117.16,1) if u∈[0,255] and for kernel h we take b=(1,1,0.46,10-4) if h∈[0,1]. That is we use
(16)u=Ta(ψ)=2571+e-2ψ/117.16-1,h=Tb(ω)=1+2·10-41+e-2ω/0.46-10-4.
Similar to the previous section, we need not vary the parameters to take into account varying noise or blur levels with the exception of b1 which may be chosen depending on the perceived level of blur. We now reformulate our old problem (4) as the new variational model:
(17)min⁡ψ,ωfTψ,ω=12Tb(ω)*Ta(ψ)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ+α2∫Ω‍∇Tb(ω)dΩ,
letting the image and the kernel be represented by u=Ta(ψ) and h=Tb(ω), respectively. Here from solving (17), the nonnegativity constraints are exactly and implicitly enforced; that is u,h≥0, but the remaining symmetry and unit integral constraints on the kernel are still required.The advantage of realizing positivityu,h>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.


### 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 u=Ta(ψ)=A1+B1ψ+O(ψ2). First, defining R1(ψ)=Ta(ψ)-B1ψ, we split the transform by
(18)Ta(ψ)=B1ψ+R1(ψ),
where R1(ψ) is nonlinear. Second, treating the above two terms differently through lagging, we propose a fixed-point lagging technique by substituting Ta(ψ,ψ~)=B1ψ+R1(ψ~) into (17) and get it linearised, leaving the remaining nonlinearity in terms with known quantities. Similarly we also have Tb(ω,ω~)=B2ω+R2(ω~) for h. Clearly Ta(ψ,ψ)=Ta(ψ)=u and Tb(ω,ω)=Tb(ω)=h. In alternating minimization, this converts problem (17) to
(19)min⁡ψ12Tb(ω,ω)*Ta(ψ,ψ~)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ,min⁡ω12Ta(ψ,ψ)*Tb(ω,ω~)-zL2(Ω)2+α2∫Ω‍∇Tb(ω)dΩ.
Our idea is to repeat the iterations (ψ~,ω~)→(ψ,ω)→(ψ~,ω~)⋯ until ψ-ψ~22+ω-ω~22 is small. Here a key observation is that residual functions R1,R2 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):
(20)B1Tb(ω~)*Tb(ω~)*Ta(ψ,ψ~)-z+α1J1(ψ~)ψ=0,B2Ta(ψ~)*Taψ~*Tbω,ω~-z+α2J2(ω~)ω=0,
which are linearised equations due to linear functions Ta(ψ,ψ~), Tb(ω,ω~)(21)B1Tbω~*Tbω~*B1ψ+R1(ψ~)-z+α1J1(ψ~)ψ=0,B2Ta(ψ~)*Ta(ψ~)*B2ω+R2(ω~)-z+α2J2(ω~)ω=0,
or simplified as, after moving known terms to the right hand sides,
(22)B12Tbω~*Tb(ω~)*ψ+α1J1(ψ~)ψ=B1z-1,B22Ta(ψ~)*Ta(ψ~)*ω+α2J2ω~ω=B2z-2,
where B1=(∂/∂ψ)Ta(ψ)ψ=0, B2=∂/∂ωTb(ω)ω=0,
(23)z-1=Tbω~*z-Tbω~*R1ψ~,z-2=Taψ~*z-Taψ~*R2ω~,J1ψ~ψ=4E~1a1+2a4(E~1-1)∇ψ1+E~13a32-∇·2a1+2a4E~11+E~12a3∇ψ~∇ψ,J2ω~ω=4E~2(b1+2b4)(E~2-1)∇ω1+E~23b32-∇·2(b1+2b4)E~21+E~22b3∇ω~∇ω,E~1=E1ψ~=a2exp⁡-2ψ~a3,E~2=E2ω~=b2exp⁡-2ω~b3.
After estimating u(0) and h(0) of the image and the kernel, respectively, we apply the inverse transforms obtaining ψ(0) and ω(0). Then on convergence, u,h are finally obtained from ψ,ω.


#### 3.3.1. Constraints forh

Note that the previously mentioned constraintsh(x,y)=h(-x,-y) and ∫Ωh(x,y)dΩ=1 take the new forms: Tb(ω(x,y))=Tb(ω(-x,-y)) and ∫Ω‍Tb(ω(x,y))dΩ=1. We can satisfy the first condition by imposing
(24)ω=12ω^(x,y)+ω^(-x,-y),
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 S^=∑x‍∑y‍Tb(ω^), our constraint is satisfied by ω from solving
(25)Tb(ω)=S^-1Tb(ω^)
and then ∫Ω‍Tb(ω(x,y))dΩ=1. The solution method is given in Algorithm 2.Algorithm 2:The first transform method for model (17) via  (20).
(1)function  Transform_Method_1 (h(0),u(0),z,α1,α2,max⁡it,TOL)
(2)ω(0)=Tb-1h(0), ψ(0)=Ta-1u(0)
(3)for  k←1tomax⁡it  do
(4)ω(k+(1/3))←SOLVE(20)
(5)ω(k+(2/3))←12ω(k+(1/3))(x,y)+ω(k+(1/3))(-x,-y)
(6)S←∑x‍∑y‍Tbω(k+(2/3))
(7)    Obtainω(k+1) from solving (25) with ω^=ω(k+(2/3)).
(8)ψ(k+1)←SOLVE(20).
(9)If   ψ(k+1)-ψ(k)+ω(k+1)-ω(k)≤TOL, then exit or continue.
(10)end for
(11) Accept the restored imageu=ψ and the restored kernel h=ω.
(12)end function


#### 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 u|Γ=h|Γ=0 which become ψ|Γ=-a3/2ln⁡⁡((a1+a4)/a4a2) and ω|Γ=-b3/2ln⁡⁡((b1+b4)/b4b2). That is, we have nonzero Dirichlet boundary conditions for the transformed variables ψ,ω.Discretisation of the Fitting (Integral) Term. We wish to discretise the quantity Tb(ω~)*ψ which is similar to h*u 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 u,h, 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
(26)A1(ψ~,ω~)ψ→=(H1TH1+α1J1)ψ→=b→1(ψ~,ω~),A2(ψ~,ω~)ω→=(H2TH2+α2J2)ω→=b→2(ψ~,ω~),
where H1, H2 are of BTTB form and and Jj 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 P1 and P2. We implement the product preconditioner following the work of [46] and given by
(27)Pj=1γjvj~c*vj~c+γjI1/2γjI+αjJjψ~·vj~c*vj~c+γjI1/2,j=1,2
for (26), where vj~c is a circulant approximation [46, 47] to vj~ defined above and γ1,γ2 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 h*u in the original variables or to remove the nonlinearity in the fitting term while still imposing the two proposed transforms u=Ta(ψ), h=Tb(ω).This is achieved by treating the transforms as constraints foru,h,ψ,ω 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:
(28)f(u,h,ψ,ω,λ1,λ2)=12h*u-zL2(Ω)2+α1J1(Ta(ψ))+γ12u-Ta(ψ)L2Ω2+λ1,u-Taψ+α2J2Tbω+γ22h-Tb(ω)L2(Ω)2+〈λ2,h-Tb(ω)〉,
where J1 and J2 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 λ1 and λ2 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(29)∂f∂u=hT(hu-z)+λ1+γ1u-B1ψ+R~1=0,∂f∂ψ=-B1λ1-γ1B1u-B1ψ+R~1+J1(ψ~)ψ=0,∂f∂h=uT(uh-z)+λ2+γ2h-B2ω+R~2=0,∂f∂ω=-B2λ2-γ2B2h-B2ω+R~2+J2(ω~)ω=0,
where B1,B2,R~1=R1(ψ~), R~2=R2(ω~) are the same as in (22):
(30)J1ψ~ψ=4a2e-2ψ~/a3(a1+2a4)a2e-2ψ~/a3-1∇ψ1+a2e-2ψ~/a33a32-∇·2(a1+2a4)a2e-2ψ~/a31+a2e-2ψ~/a32a3∇ψ~∇ψ,(31)J2ω~ω=4b2e-2ω~/b3(b1+2b4)b2e-2ω~/b3-1∇ω1+b2e-2ω~/b33b32-∇·2(b1+2b4)b2e-2ω~/b31+b2e-2ω~/b32b3∇ω~∇ω.
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.Algorithm 3:The second transform based algorithm for positivity with ADM.
(1)function  Transform_Method_2(z,α1,α2,a,b,h0,γ1,γ2,tol,max⁡it)
(2)u(0)←z,  λ1(0)←0,   λ2(0)←0
(3)Calculate a={a1,a2,a3,a4},   and  b={b1,b2,b3,b4}
(4)ψ(0)←-a3/2log⁡(a1+a4-u0)/(a2(u(0)+a4))
(5)ω(0)←-b3/2log⁡(b1+b4-h0)/(b2(h(0)+b4))
(6)for  l←1  to  max⁡it  do
(7)   Solve∂f/∂h=0  (29) for h(l+1), i.e.
u(l)Tu(l)h(l+1)-z+λ2(l)+γ2h(l+1)-B2ω(l)+R~2(l)=0
(8)for  l2←1  to  max⁡it  do
(9)      Solve∂f/∂ω=0  (29) for ω(l2+1) given ω(l2), i.e.
-B2λ2(l)-γ2B2h(l+1)-B2ω(l2+1)+R2(l2)+J2ω(l2)ω(l2+1)=0
(10) end for
(11)λ2(l+1)←λ2(l)+γ2h(l+1)-Tbω(l+1)
(12) Solve ∂f/∂u=0  (29) for u(l+1), i.e.
h(l+1)Th(l+1)u(l+1)-z+λ1(l)+γ1u(l+1)-B1ψ(l)+R~1(l)=0
(13) for  l1←1  to  max⁡it  do
(14)   Solve ∂f/∂ψ=0  (29) for ψ(l1+1) given ψ(l1), i.e.
-B1λ1(l)-γ1B1u(l+1)-B1ψ(l1+1)+R1(l1)+J1ψ(l1)ψ(l1+1)=0
(15) end for
(16) λ1(l+1)←λ1(l)+γ1u(l+1)-Taψ(l+1)
(17)end for
(18)end functionThe 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 u but it is less ideal for smooth h; 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 Kernelh

In an attempt to improve the result of recovered smooth (such as Gaussian) kernels, we introduce the following functional which uses theL2 norm to regularise the blurring kernel h=Tb(ω) and the TV to regularise the image u=Ta(ψ), as a hybrid model of (3) and (4):
(32)fTMψ,ω=12Tb(ω)*Ta(ψ)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ+α24∫Ω‍∇Tb(ω)2dΩ.Then the Euler-Lagrange equations corresponding to Method I of Section3.3 can be derived
(33)B22v2ψ,ψ~Tv2ψ,ψ~*ω-z-2ψ~,ψ,ω~+J2ω~ω=0,B12v1ω,ω~Tv1ω,ω~*ψ-z-1ω~,ω,ψ~+J1ψ~ψ=0,
where v2(ψ,ψ~)=B1ψ+R1(ψ~), z-2(ψ,ψ~ω~)=B2-1z--B2-1B1R2(ω~)*ψ,
(34)J2ω~ω=2b1+2b4b1E1+b2E2b32·2b2E-1b31+b2E∇ω·∇ω+Δω,v1ω,ω~=B2ω+R2ω~,z-1ω,ω~ψ~=B1-1z--B1-1B2R1ψ~*ω,J1ψ~ψ=4E~1(a1+2a4)(E~1-1)∇ψ1+E~13a32-∇·2(a1+2a4)E~11+E~12a3∇ψ~∇ψ.
Other computational details and use of Method II can be considered similarly.


## 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 isu=exp⁡(-ψ(x,y)) for some function ψ. Unfortunately this particular transform does not work as it is not capable of dealing with dark regions (where u≈0) in a stable way. Another choice could be u=T(ψ)=(ψ(x,y))2 as in [25, 26]. This is valid for u≈0 but this T(ψ) 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(12)u=Tψ=w+2a1+e-2ψ/ϵ-a,withiiiiiiψ=T-1u=ϵ2ln⁡a+ua+w-u,
where constants w,a,ϵ>0, max⁡u≤w, and ϵ is a small tuning parameter which controls the spread of the function. Clearly u=0 poses no problems to the transform. For the usual intensity range [0,255] for u, we may take w=255, a=1/10, ϵ=1/100:
(13)u=T(ψ)=255+1/51+e-200ψ-110,withψ=T-1(u)=1200ln⁡u+1/10255-u+1/10.
Since T(ψ) for ψ∈R and T-1(u) for u∈(-a,∞) are monotone functions, we can work out their lower and upper bounds:
(14)u=T(ψ)∈-a,w+aforψ∈R;u=T(ψ)∈[0,w+a),ifψ≥ϵ2ln⁡aa+w;ψ=T-1u∈ϵ2ln⁡aa+w,ϵ2ln⁡a+waifu∈0,w⊂-a,w+a;ψ=T-1u∈ϵ2ln⁡3aa+2w,ϵ2ln⁡a+2w3aifu∈-a2,w+a2⊂-a,w+a.
Note, for (13), (ϵ/2)ln⁡⁡((a+w)/a)=ln⁡(2551)/200≈0.039 and (ϵ/2)ln⁡((a+2w)/3a)=ln⁡(5103)/200≈0.043. 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 imageu and the kernel h, we introduce the 8 parameters with subscripts as follows:
(15)u=Ta(ψ)=a1+2a41+a2e-2ψ/a3-a4,h=Tb(ω)=b1+2b41+b2e-2ω/b3-b4,
for the image and kernel, respectively; here all constants can be fixed before proceeding (see Appendix A for more details). In particular, a1 and b1 are the expected upper limits of the image intensity values and kernel values; a2 and b2 are introduced to control the values of the image and kernel at ψ=0 and ω=0, and a3 and b3 control the spread of ψ and ω. To give one feasible set, for image u, we have a=(a1,a2,a3,a4)=(255,1,117.16,1) if u∈[0,255] and for kernel h we take b=(1,1,0.46,10-4) if h∈[0,1]. That is we use
(16)u=Ta(ψ)=2571+e-2ψ/117.16-1,h=Tb(ω)=1+2·10-41+e-2ω/0.46-10-4.
Similar to the previous section, we need not vary the parameters to take into account varying noise or blur levels with the exception of b1 which may be chosen depending on the perceived level of blur. We now reformulate our old problem (4) as the new variational model:
(17)min⁡ψ,ωfTψ,ω=12Tb(ω)*Ta(ψ)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ+α2∫Ω‍∇Tb(ω)dΩ,
letting the image and the kernel be represented by u=Ta(ψ) and h=Tb(ω), respectively. Here from solving (17), the nonnegativity constraints are exactly and implicitly enforced; that is u,h≥0, but the remaining symmetry and unit integral constraints on the kernel are still required.The advantage of realizing positivityu,h>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.


## 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 u=Ta(ψ)=A1+B1ψ+O(ψ2). First, defining R1(ψ)=Ta(ψ)-B1ψ, we split the transform by
(18)Ta(ψ)=B1ψ+R1(ψ),
where R1(ψ) is nonlinear. Second, treating the above two terms differently through lagging, we propose a fixed-point lagging technique by substituting Ta(ψ,ψ~)=B1ψ+R1(ψ~) into (17) and get it linearised, leaving the remaining nonlinearity in terms with known quantities. Similarly we also have Tb(ω,ω~)=B2ω+R2(ω~) for h. Clearly Ta(ψ,ψ)=Ta(ψ)=u and Tb(ω,ω)=Tb(ω)=h. In alternating minimization, this converts problem (17) to
(19)min⁡ψ12Tb(ω,ω)*Ta(ψ,ψ~)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ,min⁡ω12Ta(ψ,ψ)*Tb(ω,ω~)-zL2(Ω)2+α2∫Ω‍∇Tb(ω)dΩ.
Our idea is to repeat the iterations (ψ~,ω~)→(ψ,ω)→(ψ~,ω~)⋯ until ψ-ψ~22+ω-ω~22 is small. Here a key observation is that residual functions R1,R2 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):
(20)B1Tb(ω~)*Tb(ω~)*Ta(ψ,ψ~)-z+α1J1(ψ~)ψ=0,B2Ta(ψ~)*Taψ~*Tbω,ω~-z+α2J2(ω~)ω=0,
which are linearised equations due to linear functions Ta(ψ,ψ~), Tb(ω,ω~)(21)B1Tbω~*Tbω~*B1ψ+R1(ψ~)-z+α1J1(ψ~)ψ=0,B2Ta(ψ~)*Ta(ψ~)*B2ω+R2(ω~)-z+α2J2(ω~)ω=0,
or simplified as, after moving known terms to the right hand sides,
(22)B12Tbω~*Tb(ω~)*ψ+α1J1(ψ~)ψ=B1z-1,B22Ta(ψ~)*Ta(ψ~)*ω+α2J2ω~ω=B2z-2,
where B1=(∂/∂ψ)Ta(ψ)ψ=0, B2=∂/∂ωTb(ω)ω=0,
(23)z-1=Tbω~*z-Tbω~*R1ψ~,z-2=Taψ~*z-Taψ~*R2ω~,J1ψ~ψ=4E~1a1+2a4(E~1-1)∇ψ1+E~13a32-∇·2a1+2a4E~11+E~12a3∇ψ~∇ψ,J2ω~ω=4E~2(b1+2b4)(E~2-1)∇ω1+E~23b32-∇·2(b1+2b4)E~21+E~22b3∇ω~∇ω,E~1=E1ψ~=a2exp⁡-2ψ~a3,E~2=E2ω~=b2exp⁡-2ω~b3.
After estimating u(0) and h(0) of the image and the kernel, respectively, we apply the inverse transforms obtaining ψ(0) and ω(0). Then on convergence, u,h are finally obtained from ψ,ω.


### 3.3.1. Constraints forh

Note that the previously mentioned constraintsh(x,y)=h(-x,-y) and ∫Ωh(x,y)dΩ=1 take the new forms: Tb(ω(x,y))=Tb(ω(-x,-y)) and ∫Ω‍Tb(ω(x,y))dΩ=1. We can satisfy the first condition by imposing
(24)ω=12ω^(x,y)+ω^(-x,-y),
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 S^=∑x‍∑y‍Tb(ω^), our constraint is satisfied by ω from solving
(25)Tb(ω)=S^-1Tb(ω^)
and then ∫Ω‍Tb(ω(x,y))dΩ=1. The solution method is given in Algorithm 2.Algorithm 2:The first transform method for model (17) via  (20).
(1)function  Transform_Method_1 (h(0),u(0),z,α1,α2,max⁡it,TOL)
(2)ω(0)=Tb-1h(0), ψ(0)=Ta-1u(0)
(3)for  k←1tomax⁡it  do
(4)ω(k+(1/3))←SOLVE(20)
(5)ω(k+(2/3))←12ω(k+(1/3))(x,y)+ω(k+(1/3))(-x,-y)
(6)S←∑x‍∑y‍Tbω(k+(2/3))
(7)    Obtainω(k+1) from solving (25) with ω^=ω(k+(2/3)).
(8)ψ(k+1)←SOLVE(20).
(9)If   ψ(k+1)-ψ(k)+ω(k+1)-ω(k)≤TOL, then exit or continue.
(10)end for
(11) Accept the restored imageu=ψ and the restored kernel h=ω.
(12)end function


### 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 u|Γ=h|Γ=0 which become ψ|Γ=-a3/2ln⁡⁡((a1+a4)/a4a2) and ω|Γ=-b3/2ln⁡⁡((b1+b4)/b4b2). That is, we have nonzero Dirichlet boundary conditions for the transformed variables ψ,ω.Discretisation of the Fitting (Integral) Term. We wish to discretise the quantity Tb(ω~)*ψ which is similar to h*u 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 u,h, 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
(26)A1(ψ~,ω~)ψ→=(H1TH1+α1J1)ψ→=b→1(ψ~,ω~),A2(ψ~,ω~)ω→=(H2TH2+α2J2)ω→=b→2(ψ~,ω~),
where H1, H2 are of BTTB form and and Jj 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 P1 and P2. We implement the product preconditioner following the work of [46] and given by
(27)Pj=1γjvj~c*vj~c+γjI1/2γjI+αjJjψ~·vj~c*vj~c+γjI1/2,j=1,2
for (26), where vj~c is a circulant approximation [46, 47] to vj~ defined above and γ1,γ2 are positive constants.


## 3.3.1. Constraints forh

Note that the previously mentioned constraintsh(x,y)=h(-x,-y) and ∫Ωh(x,y)dΩ=1 take the new forms: Tb(ω(x,y))=Tb(ω(-x,-y)) and ∫Ω‍Tb(ω(x,y))dΩ=1. We can satisfy the first condition by imposing
(24)ω=12ω^(x,y)+ω^(-x,-y),
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 S^=∑x‍∑y‍Tb(ω^), our constraint is satisfied by ω from solving
(25)Tb(ω)=S^-1Tb(ω^)
and then ∫Ω‍Tb(ω(x,y))dΩ=1. The solution method is given in Algorithm 2.Algorithm 2:The first transform method for model (17) via  (20).
(1)function  Transform_Method_1 (h(0),u(0),z,α1,α2,max⁡it,TOL)
(2)ω(0)=Tb-1h(0), ψ(0)=Ta-1u(0)
(3)for  k←1tomax⁡it  do
(4)ω(k+(1/3))←SOLVE(20)
(5)ω(k+(2/3))←12ω(k+(1/3))(x,y)+ω(k+(1/3))(-x,-y)
(6)S←∑x‍∑y‍Tbω(k+(2/3))
(7)    Obtainω(k+1) from solving (25) with ω^=ω(k+(2/3)).
(8)ψ(k+1)←SOLVE(20).
(9)If   ψ(k+1)-ψ(k)+ω(k+1)-ω(k)≤TOL, then exit or continue.
(10)end for
(11) Accept the restored imageu=ψ and the restored kernel h=ω.
(12)end function


## 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 u|Γ=h|Γ=0 which become ψ|Γ=-a3/2ln⁡⁡((a1+a4)/a4a2) and ω|Γ=-b3/2ln⁡⁡((b1+b4)/b4b2). That is, we have nonzero Dirichlet boundary conditions for the transformed variables ψ,ω.Discretisation of the Fitting (Integral) Term. We wish to discretise the quantity Tb(ω~)*ψ which is similar to h*u 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 u,h, 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
(26)A1(ψ~,ω~)ψ→=(H1TH1+α1J1)ψ→=b→1(ψ~,ω~),A2(ψ~,ω~)ω→=(H2TH2+α2J2)ω→=b→2(ψ~,ω~),
where H1, H2 are of BTTB form and and Jj 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 P1 and P2. We implement the product preconditioner following the work of [46] and given by
(27)Pj=1γjvj~c*vj~c+γjI1/2γjI+αjJjψ~·vj~c*vj~c+γjI1/2,j=1,2
for (26), where vj~c is a circulant approximation [46, 47] to vj~ defined above and γ1,γ2 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 h*u in the original variables or to remove the nonlinearity in the fitting term while still imposing the two proposed transforms u=Ta(ψ), h=Tb(ω).This is achieved by treating the transforms as constraints foru,h,ψ,ω 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:
(28)f(u,h,ψ,ω,λ1,λ2)=12h*u-zL2(Ω)2+α1J1(Ta(ψ))+γ12u-Ta(ψ)L2Ω2+λ1,u-Taψ+α2J2Tbω+γ22h-Tb(ω)L2(Ω)2+〈λ2,h-Tb(ω)〉,
where J1 and J2 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 λ1 and λ2 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(29)∂f∂u=hT(hu-z)+λ1+γ1u-B1ψ+R~1=0,∂f∂ψ=-B1λ1-γ1B1u-B1ψ+R~1+J1(ψ~)ψ=0,∂f∂h=uT(uh-z)+λ2+γ2h-B2ω+R~2=0,∂f∂ω=-B2λ2-γ2B2h-B2ω+R~2+J2(ω~)ω=0,
where B1,B2,R~1=R1(ψ~), R~2=R2(ω~) are the same as in (22):
(30)J1ψ~ψ=4a2e-2ψ~/a3(a1+2a4)a2e-2ψ~/a3-1∇ψ1+a2e-2ψ~/a33a32-∇·2(a1+2a4)a2e-2ψ~/a31+a2e-2ψ~/a32a3∇ψ~∇ψ,(31)J2ω~ω=4b2e-2ω~/b3(b1+2b4)b2e-2ω~/b3-1∇ω1+b2e-2ω~/b33b32-∇·2(b1+2b4)b2e-2ω~/b31+b2e-2ω~/b32b3∇ω~∇ω.
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.Algorithm 3:The second transform based algorithm for positivity with ADM.
(1)function  Transform_Method_2(z,α1,α2,a,b,h0,γ1,γ2,tol,max⁡it)
(2)u(0)←z,  λ1(0)←0,   λ2(0)←0
(3)Calculate a={a1,a2,a3,a4},   and  b={b1,b2,b3,b4}
(4)ψ(0)←-a3/2log⁡(a1+a4-u0)/(a2(u(0)+a4))
(5)ω(0)←-b3/2log⁡(b1+b4-h0)/(b2(h(0)+b4))
(6)for  l←1  to  max⁡it  do
(7)   Solve∂f/∂h=0  (29) for h(l+1), i.e.
u(l)Tu(l)h(l+1)-z+λ2(l)+γ2h(l+1)-B2ω(l)+R~2(l)=0
(8)for  l2←1  to  max⁡it  do
(9)      Solve∂f/∂ω=0  (29) for ω(l2+1) given ω(l2), i.e.
-B2λ2(l)-γ2B2h(l+1)-B2ω(l2+1)+R2(l2)+J2ω(l2)ω(l2+1)=0
(10) end for
(11)λ2(l+1)←λ2(l)+γ2h(l+1)-Tbω(l+1)
(12) Solve ∂f/∂u=0  (29) for u(l+1), i.e.
h(l+1)Th(l+1)u(l+1)-z+λ1(l)+γ1u(l+1)-B1ψ(l)+R~1(l)=0
(13) for  l1←1  to  max⁡it  do
(14)   Solve ∂f/∂ψ=0  (29) for ψ(l1+1) given ψ(l1), i.e.
-B1λ1(l)-γ1B1u(l+1)-B1ψ(l1+1)+R1(l1)+J1ψ(l1)ψ(l1+1)=0
(15) end for
(16) λ1(l+1)←λ1(l)+γ1u(l+1)-Taψ(l+1)
(17)end for
(18)end functionThe 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 u but it is less ideal for smooth h; 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 Kernelh

In an attempt to improve the result of recovered smooth (such as Gaussian) kernels, we introduce the following functional which uses theL2 norm to regularise the blurring kernel h=Tb(ω) and the TV to regularise the image u=Ta(ψ), as a hybrid model of (3) and (4):
(32)fTMψ,ω=12Tb(ω)*Ta(ψ)-zL2(Ω)2+α1∫Ω‍∇Ta(ψ)dΩ+α24∫Ω‍∇Tb(ω)2dΩ.Then the Euler-Lagrange equations corresponding to Method I of Section3.3 can be derived
(33)B22v2ψ,ψ~Tv2ψ,ψ~*ω-z-2ψ~,ψ,ω~+J2ω~ω=0,B12v1ω,ω~Tv1ω,ω~*ψ-z-1ω~,ω,ψ~+J1ψ~ψ=0,
where v2(ψ,ψ~)=B1ψ+R1(ψ~), z-2(ψ,ψ~ω~)=B2-1z--B2-1B1R2(ω~)*ψ,
(34)J2ω~ω=2b1+2b4b1E1+b2E2b32·2b2E-1b31+b2E∇ω·∇ω+Δω,v1ω,ω~=B2ω+R2ω~,z-1ω,ω~ψ~=B1-1z--B1-1B2R1ψ~*ω,J1ψ~ψ=4E~1(a1+2a4)(E~1-1)∇ψ1+E~13a32-∇·2(a1+2a4)E~11+E~12a3∇ψ~∇ψ.
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 u and the kernel h, given the received image z. 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 consider6 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 103 maximum iterations although this number is rarely reached.


### 4.1. Test of Robustness of Algorithm2

Set 1 (simple image with blurs).
Result Set1 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.Experimental results on Image 1 with Blur 1. Row 1, l-r: Image 1, received data corrupted by Blur 1, restored image using Algorithm2. 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.
(a)
True image
(b)
Received SNR = 13.5
(c)
Algorithm1 [3] SNR = 9.8
(d)
Algorithm2 SNR = 21.9
(e)
True image
(f)
Received PSNR = 7.9
(g)
Algorithm2 PSNR = 27.3Set 2 (detailed image containing many zero-points with blurs).
Result Set2 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.Experimental results on Images 2–4 with Blur 1. Row 1, l-r: Image 2, received data corrupted by Blur 1, restored image using Algorithm2. 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).
(a)
True image
(b)
Received data
(c)
Algorithm2 restored image
(d)
True image
(e)
Received data
(f)
Algorithm2 restored image
(g)
True image
(h)
Received data
(i)
Algorithm2 restored image
(j)
True image
(k)
Received data
(l)
Algorithm2 restored imageSet 3 (detailed photograph images with blurs).
Result Set3 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 Set4 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.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 Algorithm1 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)
Result by Algorithm1
(d)
Result by Algorithm2
(e)
True image
(f)
Received data
(g)
Restored imageExperimental results on Images 3-4 with Blur 2, restored using Algorithm2. 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)
Algorithm2 restored image
(d)
True image
(e)
Received data
(f)
Algorithm2 restored imageTest results on Images 5-6 with Blur 2. Row 1, l-r for Image 5, received dataz, restored image u2 using Algorithm 2 and u3 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.
(a)
True image
(b)
Received data
(c)
Algorithm2 PSNR = 23.4
(d)
Algorithm3 PSNR = 20.8
(e)
True image
(f)
Received data
(g)
Algorithm2 PSNR = 21.6


### 4.2. Comparison of Algorithms2 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, Algorithm3 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(35)ReRes=Nz(k)-g2Nz0-g2,SNR=10log⁡10u-umean22u-u~22,PSNR=10log⁡10mnmax⁡u~,u2u~-u22,
where u and u~ are the true image and the restored image, respectively, and umean is the mean value of the original image.Finally we are ready to show in Table1 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 Algorithms2 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
Receivedz
Algorithm2
Algorithm3
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


## 4.1. Test of Robustness of Algorithm2

Set 1 (simple image with blurs).
Result Set1 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.Experimental results on Image 1 with Blur 1. Row 1, l-r: Image 1, received data corrupted by Blur 1, restored image using Algorithm2. 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.
(a)
True image
(b)
Received SNR = 13.5
(c)
Algorithm1 [3] SNR = 9.8
(d)
Algorithm2 SNR = 21.9
(e)
True image
(f)
Received PSNR = 7.9
(g)
Algorithm2 PSNR = 27.3Set 2 (detailed image containing many zero-points with blurs).
Result Set2 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.Experimental results on Images 2–4 with Blur 1. Row 1, l-r: Image 2, received data corrupted by Blur 1, restored image using Algorithm2. 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).
(a)
True image
(b)
Received data
(c)
Algorithm2 restored image
(d)
True image
(e)
Received data
(f)
Algorithm2 restored image
(g)
True image
(h)
Received data
(i)
Algorithm2 restored image
(j)
True image
(k)
Received data
(l)
Algorithm2 restored imageSet 3 (detailed photograph images with blurs).
Result Set3 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 Set4 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.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 Algorithm1 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)
Result by Algorithm1
(d)
Result by Algorithm2
(e)
True image
(f)
Received data
(g)
Restored imageExperimental results on Images 3-4 with Blur 2, restored using Algorithm2. 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)
Algorithm2 restored image
(d)
True image
(e)
Received data
(f)
Algorithm2 restored imageTest results on Images 5-6 with Blur 2. Row 1, l-r for Image 5, received dataz, restored image u2 using Algorithm 2 and u3 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.
(a)
True image
(b)
Received data
(c)
Algorithm2 PSNR = 23.4
(d)
Algorithm3 PSNR = 20.8
(e)
True image
(f)
Received data
(g)
Algorithm2 PSNR = 21.6


## 4.2. Comparison of Algorithms2 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, Algorithm3 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(35)ReRes=Nz(k)-g2Nz0-g2,SNR=10log⁡10u-umean22u-u~22,PSNR=10log⁡10mnmax⁡u~,u2u~-u22,
where u and u~ are the true image and the restored image, respectively, and umean is the mean value of the original image.Finally we are ready to show in Table1 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 Algorithms2 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
Receivedz
Algorithm2
Algorithm3
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.



---

*Source: 860263-2015-02-15.xml*