# Bayesian Image Restoration Using a Large-Scale Total Patch Variation Prior

Authors: Yang Chen; Weimin Yu; Yinsheng Li; Zhou Yang; Limin Luo; Wufan Chen
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

Edge-preserving Bayesian restorations using nonquadratic priors are often inefficient in restoring continuous variations and tend to produce block artifacts around edges in ill-posed inverse image restorations. To overcome this, we have proposed a spatial adaptive (SA) prior with improved performance. However, this SA prior restoration suffers from high computational cost and the unguaranteed convergence problem. Concerning these issues, this paper proposes a Large-scale Total Patch Variation (LS-TPV) Prior model for Bayesian image restoration. In this model, the prior for each pixel is defined as a singleton conditional probability, which is in a mixture prior form of one patch similarity prior and one weight entropy prior. A joint MAP estimation is thus built to ensure the iteration monotonicity. The intensive calculation of patch distances is greatly alleviated by the parallelization of Compute Unified Device Architecture(CUDA). Experiments with both simulated and real data validate the good performance of the proposed restoration.


---

## Body

## 1. Introduction

Image restorations have wide applications in remote sensing, radar imaging, tomographic imaging, microscopic imaging, astronomic imaging, digital photography, and so forth [1–4]. For linear and shift invariant imaging systems, the transformation from f to g is well described by the following additive linear degradation model [3, 4]:(1.1)g=A*f+ɛ,
where g, f, and ɛ represent the degraded observed image, the original true image, and the corrupting white Gaussian noise with variance σ2, respectively. A is the PSF (Point Spread function) of the imaging system and * is the linear convolution operator. Throughout this paper, we assume that the PSF A is known or could be numerically determined or estimated.Bayesian or maximum a posteriori (MAP) approach, within Markov Random Fields (MRFs) framework, can provide a stable solution of the ill-posed inverse image restorations through incorporating a priori information about the geometrical properties of an image [3–8]. Via modelling the unknown parameters in the prior probability density functions, prior information can also be interpreted as the regularization metrics, which measures the extent to which the contextual constraint assumption is violated. With the incorporated prior information, Bayesian approach is able to distinguish solutions from less desirable ones by transforming the original ill-posed problem into a well-posed one. Generally, we can build the following posterior probability P(f∣g) for image restorations,(1.2)P(f∣g)∝P(g∣f)P(f),(1.3)P(g∣f)=exp⁡(L(g,f))=exp⁡(-12‖g-A*f‖2),(1.4)P(f)=Z-1×exp⁡(-βU(f))=Z-1×exp⁡(-β∑jUj(f)).
Here, P(g∣f) and P(f) denote the likelihood distribution and the prior distribution, respectively. L(g,f) and U(f) denote the corresponding likelihood energy and the prior energy function. Uj(f) is the notation of the energy function U evaluated at pixel index j. The partition function Z is a normalizing constant. β is the global parameter controlling the balance between the likelihood energy and the prior energy. We can build the posterior energy function in the logarithm form(1.5)ψβ(f)=log⁡⁡P(f∣g)=L(g,f)-βU(f)=-12‖g-A*f‖2-β∑jUj(f).
The restored or deconvolved image f can be obtained through the maximization of function ψβ(f).The widely used quadratic membrane (QM) prior or the Tikhonov L2 regularization tend to smooth both noise and edge details and often lead to unfavorable oversmoothing in restorations [5]. On the other side, some edge-preserving Bayesian restoration methods have been proposed, which can be classified into three main categories: wavelet relevant regularization restorations, Bayesian restorations with nonquadratic prior energies, and Bayesian restorations with total variation (TV) prior. Wavelet relevant regularization restorations, with multiscale stochastic prior models, have been proposed in the research of edge-preserving restorations [8–11]. In [10], priors based on wavelet decompositions and heavy-tailed pdfs were incorporated into the EM restoration algorithm. In [11], Neelamani et al. proposed a Fourier-Wavelet hybrid restoration, which exploits the potentials of the sparse representation of noise for Fourier transform and the sparse representation of the coefficients in wavelet transform. All these wavelet regularization methods are based on the distribution modeling and thresholding of the decomposed wavelet coefficients [11]. Edge-preserving priors with nonquadratic energies were also proposed to preserve the edge details by using nonlinear inverse-proportional functions between edge existences and intensity differences [5, 12, 13]. The weighting matrices in nonquadratic prior Bayesian restorations preserve the edges by turning off or suppressing smoothing at appropriate locations. Another impressive research direction in this area is the total variation (TV)-based image restorations [7, 14, 15]. This kind of approach can also be viewed as TV priors with the half-quadratic regularization prior energies. With no enforced global image continuity, TV prior restorations often demonstrate good property of edge preserving.The nonquadratic priors preserve the edges by determining the pixel-wise regularizations based on the intensity-difference information within local fixed neighborhoods. The regularizations on TV priors also come from the information of the local intensity gradients. For all these nonquadratic and TV priors, the local information of intensity difference can by no means provide effective regularizations to discriminate noise from image textures in different scales. So these edge-preserving approaches often fail to cope well with some complicated restorations with high noise contaminations, and tend to produce block artifacts around the continuous edges in those situations [16–21].Objects in the images are reasonably assumed to be a composition of edges and backgrounds with intensities changing coherently over different scales. Thus, for the traditional prior model, the intensity differences between individual pixels within local neighborhoods are often insufficient to characterize objects [19, 20]. This limit elicited our previous work in [21] which introduced a spatial adaptive (SA) prior restoration approach. This SA prior works by adaptively including the relevant neighboring pixels and excluding the negative irrelevant ones within a large neighborhood. Though showing expressive results, the iterative restoration in [21] updates the prior weight using the latest image estimate in a one step late (OSL) way, which leads to inconsistent posterior energy function and nonconvergent iterations. In addition, the approach in [21] is greatly limited by the high computational cost in calculating all the patch distances for all the pixel pairs in each prior neighborhood. In this paper, to overcome the first problem, we applied a large-scale total patch variation (LS-TPV) prior model, which is built through the specification of the conditional probability for each pixel [22–25]. The joint maximization estimation of this constrained entropy function will lead to convergent image/weight restorations. To overcome the second problem, we introduce the parallel computation structure of CUDA to calculate the patch distances [26–28]. Furthermore, another benefit from the CUDA-accelerated parallelization is that it allows further performance enhancement for the LS-TPV prior restoration by using larger search neighborhoods.In Section2, a review of some previous prior models is illustrated, and after that we introduce the proposed LS-TPV prior model. In Section 3, we give a joint estimation algorithm for the proposed restoration. In Section 4, we perform comparative experiments with both simulated and real data. Relevant visual/quantitative results are also presented. Conclusion and discussion are given in Section 5.


## 2. Prior Model

Conventionally, the value ofUj(f) is commonly computed through a weighted sum of potential functions v of the differences between the pixels in the neighborhood Nj:(2.1)Uj(f)=∑b∈Njwbjv(fb-fj).Generally, different choices of the potential functionv lead to different priors. For the simple space-invariant QM prior, the potential function has the form v(t)=t2/2. Some Edge-preserving nonquadratic priors could be chosen by adopting a nonquadratic potential function v, such as the Huber potential function(2.2)v(t)={t22,|t|≤γ,γ|t|-γ22,|t|>γ,
where γ is the threshold parameter [5, 13]. Such edge-preserving nonquadratic priors preserve structure information by choosing nonquadratic potential functions that increase less as the differences between the adjacent pixels become bigger. Weight wbj is a positive value that denotes the interaction degree between pixel b and pixel j. In traditional prior models, it is simply considered to be inversely proportional to the distance between pixels b and j. So on a square lattice of image f, in which the 2D positions of pixels b and j (b≠j) are, respectively, (bx,by) and (jx,jy), wbj is usually calculated by the geometric distance 1/(bx-jx)2+(by-jy)2.Restorations using total variation (TV) prior take the prior energy as the following expression:(2.3)UTV(f)=∑j(Δjhf)2+(Δjvf)2,
where Δjh and Δjv are linear operators corresponding to the first-order horizontal and vertical intensity differences at pixel j, respectively [7, 14, 15]. Since both the smooth and sharp edges tend to have similar prior energies, this TV prior energy given by (2.3) does not penalize weak discontinuities and often favors images with bounded variations.Traditional quadratic priors, nonquadratic priors and TV prior, only provide local intensity information for Bayesian restorations. And the edge-preserving nonquadratic and TV priors, though being able to preserve edge information, tend to produce unfavorable block artifacts or false edges in high-noise situations.To improve restoration quality, in [21] we proposed a spatial adaptive (SA) prior model, in which a spatial adaptive prior neighborhood is implemented by setting the weight wbj to an 1-0 binary function to classify the pixels in neighborhood N into the pixels relevant to the center pixel j and those not. The building of the SA prior can be formalized as follows:(2.4)USA(f)=∑jUj(f)=∑j∑b∈Nj(wbj(fb-fj)2|Nrj|),(2.5)wbj={1,|nb(f)-nj(f)|2≤δ0,|nb(f)-nj(f)|2>δ,(2.6)|nb(f)-nj(f)|2=∑l(fl∈nb-fl∈nj)2.
Here, USA is the energy function for the SA prior. wbj, which represents the classification of the neighboring pixels in the search neighborhood Nj, can be computed via (2.5) and (2.6). |Nrj|, the number of the neighboring pixels with nonzero wbj in the neighborhood Nj, is the normalization factor for the different sizes of non-1 neighboring pixels in each Nj. Parameter δ in (2.5) is the threshold parameter. The distance |nb(f)-nj(f)|2 is determined by the distance measure (2.6) between the two translated comparing patches nb and nj, whose centers are, respectively, located at pixel b and pixel j.Though proved to be effective in suppressing both noise and block artifacts, all the weightsw in SA prior need to be heuristically inferred using current image estimate in the way of one step late (OSL), which often leads to inconsistent posterior energy function with no guaranteed convergence [8]. To overcome this problem, in this paper, we proposed an LS-TPV prior model, in which regularization takes effect through penalizing the weighted distances/variations between the patches surrounding each neighboring pixel pair. The weights w in each N are considered variables with a nonlinear dependence on image f. For the central pixel j and its neighboring pixel b in each prior neighborhood Nj, wbj is considered the similarity probability between the two patches nb and nj surrounding pixels b and j, respectively [24]. An entropy prior of all the weights wbj for each Nj is introduced into our model in the form of -∑j∑b∈Njwbjlnwbj. Here as the “Frame” model in [25], the constraint is ∑b∈Njwbj=1. Then, based on joint MAP theorem [12], we can obtain(2.7)P(f,w,η∣g)∝P(g∣f)P(f∣w,η)P(w,η).
Here, η is the introduced pixel-wise Lagrange-multiplier parameter. For each pixel j, with f̂ denoting the current image estimate, the prior energy ULS-TPV(fj,w,η∣f̂) of the LS-TPV prior is in the combined form of a weighted patch similarity energy, an entropy energy −∑j∑b∈Njwbjlnwbj and a normalization constraint ∑b∈Njwbj=1. We can thus obtain the following posterior energy function:(2.8)ψβ(f,w,η∣g)=log⁡⁡P(g∣f)+log⁡⁡P(f∣w,η)+log⁡⁡P(w,η)=L(g,f)-β(∑jULS-TPV(fj,w,η∣f̂))=-12‖g-A*f‖2-β∑j((∑b∈Njwbjdbj+h∑b∈Njwbjln⁡⁡wbj)+ηj(∑b∈Njwbj-1)).
Here, as (1.4), β is the global parameter controlling the balance between the likelihood energy and the prior energy. The patch distance dbj is calculated as ∥Gα*nb-Gα*nj∥2, in which Gα denotes the Gaussian convolution kernel with standard deviation α. ηj is the value of Lagrange-multiplier parameter at pixel j. The parameter h controls the maximum entropy constraint on w, which is routinely modulated with respect to the noise levels (large h is routinely set to suppress more noise for high-noise situations).


## 3. Restoration Algorithm

To obtain tractable maximization of the posterior probability using the proposed LS-TPV prior energyULS-TPV in (2.8), we use the algorithm of iterative coordinate descent (ICD), in which each element is sequentially updated by maximizing the log-style local conditional probability function [29]. In this way, the objective image f and the weight w are jointly estimated through alternative maximization of the ψβ(f,w,η∣g) with respect to f and w.(1)argmaxw-ULS-TPV(w,η∣f̂).
f̂ denotes the fixed current image estimation, and each weight wbj for each pixel j is updated by solving the following singleton optimization problem:
(3.1)argmax⁡wbj(-∑b∈Njwbj‖Gα*nb-Gα*nj‖2-h∑b∈Njwbjln⁡⁡wbj-ηj(∑b∈Njwbj-1)).
With Lagrange-multiplierηj introduced, by joint solving (3.1) with respect to each wbj and ηj, we can then obtain wbj:
(3.2)wbj=exp⁡(-‖Gα*nb-Gα*nj‖2/h)∑b′∈Njexp⁡(-‖Gα*nb′-Gα*nj‖2/h),
where ∥Gα*nb-Gα*nj∥2 can be easily calculated using currently available f̂. Each wbj is determined not only by the similarity metrics between patches nb and nj, but also by the sum of all the similarity metrics between patches with the central pixels located in Nj. So wbj is different from wjb.(2)argmaxfψβ(f,ŵ,η̂∣g).
Hereŵ denotes the fixed current weight estimate. Omitting the terms with no f, the ψβ(f,ŵ∣g) becomes
(3.3)ψβ(f,ŵ∣g)=L(g,f)-β∑j(∑b∈Njŵbj‖Gα*nb-Gα*nj‖2).
To make a tractable maximization of (2.8), based on the theories of pseudolikelihood and iterated conditional mode (ICM) in [29], we can restore each fj by solving the following factorized problem:
(3.4)argmax⁡fjψβ(fj,ŵbj∣g,f̂)⟹argmax⁡fjL(f,g)-β(∑b∈Njŵbj‖Gα*nb-Gα*nj‖2).Considering the convexity for∥g-A*f∥2 and ∥Gα*nb-Gα*nj∥2, we conclude that, given the fixed weight ŵ, the Hessian matrix of each ψβ(fj,ŵbj∣g,f̂) in (3.4) is a negative definite. In this step of image updating, we applied the restoration algorithm in [4] to solve (3.4) to obtain the image estimate of each pixel fj. The ψβ(f,w,η∣g) in (2.8) is separately convex with respect to each f and w. Denoting f̂n,ŵn and f̂n+1,ŵn+1 as the nth and (n+1)th iterated estimates, we can obtain following relation: ψβ(f̂n+1,ŵn+1∣g)≥ψβ(f̂n,ŵn+1∣g)≥ψβ(f̂n,ŵn∣g), which implies convergence to one local maximum and can be obtained for this joint updating strategy [30].In each iteration in the proposed restoration, it is very expensive to perform the pixel-wise calculations of all thedbj and wbj for each pair of translated patches n in each N over the whole image region. In the experiments, to save computation cost, we do not extend the sizes of N over the whole image region, and set them to some appropriate sizes. We achieve significant acceleration by using parallel Compute Unified Device Architecture (CUDA) [26–28]. For each iteration in the restoration using the proposed LS-TPV prior model, all threads in this block grid structure execute simultaneously to perform all the involved pixel-wise operations. Furthermore, considering the symmetry property that distance dbj equals to distance djb, we can further save one half computation cost by only calculating one of the two distances dbj and djb.


## 4. Experiments

Both simulated and real data are used in the experiments. To provide comparative results, we also performed Wiener deconvolution, ForWaRD restoration, Bayesian restoration using TV prior, Huber prior, and the SA prior in [21]. The Wiener deconvolution was performed by using the degraded image to estimate the power spectrum in Fourier domain. The TWIST algorithm in [7] and Newton-Raphson algorithm in [12] are used in the restorations using TV prior and Huber prior, respectively. Restored images from above Wiener filter method are taken as the initial images for all the Bayesian restorations. In experiment, the parameters needed to be preset include the β for all Bayesian restorations, the threshold γ for Huber prior, the wavelet threshold, and decomposition level in the ForWaRD method, the δ for the SA prior, the h and the n and N settings for the proposed LS-TPV prior. In the simulation experiment, all the parameters are set based on SNR maximization criteria with SNR calculated as(4.1)SNR=10log⁡10⁡(∑jM(Fj-F¯)2∑jM(Fj-f̂j)2).
Here M, f̂, F, and F¯ denote the total pixel number in image, the restored image, the original true image, and the mean of the original true image, respectively. As to the real data experiment, which has no available true image for SNR calculations, all the parameters are set based on the tradeoff between edge preservation and noise suppression. N and n of sizes (11×11N and 7×7n) is used in our work. A full study of the effects of the sizes of N and n on the resulting restorations is performed in Section 4.5.Generally, the solution to the restoration is considered optimized when the iteration goes stable. So for the Bayesian restorations using TV prior, Huber prior, and SA prior in which the energy functions are maximized with respect only to image, we stop the iterations when the current image estimates satisfy the conditionψβ(f̂n+1∣g)-ψβ(f̂n∣g)<0.995(ψβ(f̂n∣g)-ψβ(f̂n-1∣g)). For the restorations using the proposed LS-TPV prior in which the energy functions are maximized with respect to both image and weights, we stop the iterations when the current image and weight estimates satisfy the condition ψβ(f̂n+1,ŵn+1∣g)-ψβ(f̂n,ŵn∣g)<0.995(ψβ(f̂n,ŵn∣g)-ψβ(f̂n-1,ŵn-1∣g)). In practical experiments, to reach the above iteration conditions, it is found that 411, 433, 475, and 486 iterations are required for the Bayesian restorations using TV prior, Huber prior, SA prior, and the proposed LS-TPV prior, respectively.


### 4.1. Simulation Experiment

In this section, simulated experiments with256×256 “Lena” image (intensity range [−98    117]) are performed. In the simulation, a uniform 9×9 PSF is assumed, and Gaussian noise with variance = 0.1663 is added. The PSF used in this experiment is normalized to 1. Figure 1(b) is the degraded image simulated by (1.1) using the PSF with the white Gaussian noise. The parameters for different restorations are listed in Table 1. All the restored images are illustrated in Figures 1(c)–1(h), respectively.Table 1
Parameter settings in simulated data experiment.
TV priorHuber priorForWaRDSA priorLS-TPV priorβ=0.0038β=0.026,γ=0.35Threshold = 3.0,Levels = 3β=0.015,δ=900,7×7n,11×11Nβ=0.075,α=0.5,h=150,7×7n,11×11N(a) Original “Lena” image. (b) Degraded “Lena” image (SNR = 6.65) with9×9 uniform blur and additive white Gaussian noise (variance = 0.1663). (c) Wiener deconvolution. (d) TV prior Bayesian restoration. (e) Nonquadratic Huber prior Bayesian restoration. (f) The ForWaRD restoration in [11]. (g) The SA prior Bayesian restoration proposed in [21]. (h) The proposed LS-TPV prior Bayesian restoration (The left arrows: block-suppressing illustration, the right arrows: edge-preserving illustration).
(a)
Original image(b)
Degraded image(c)
Wiener restoration(d)
TV prior restoration(e)
Huber prior restoration(f)
ForWaRD restoration [11](g)
SA prior restoration [21](h)
LS-TPV prior restorationThe restored images from TV prior restoration and Huber prior restoration are, respectively, shown in Figures1(d) and 1(e). Figure 1(f) shows the restored image using the ForWaRD method. The result of Bayesian restoration using the SA prior is illustrated in Figure 1(g). For the Bayesian restoration using the proposed LS-TPV prior, the algorithm in Section 3 is used, and the corresponding restored image is shown in Figure 1(h).Through the results, we find the ForWaRD method and all the Bayesian restorations (Figures1(d)–1(h)) can overcome the ring effects in Wiener deconvolution (Figure 1(c)). The block artifacts (the left arrows), which are observed in the ForWaRD approaches and the restorations using Huber prior and TV prior, are effectively suppressed using the restorations using the SA prior and the proposed LS-TPV prior. The LS-TPV prior (Figure 1(h)) presents a further visual enhancement over the SA prior (Figure 1(g)) in block suppressing (the left arrows) and edge preserving (the right arrows). Figure 2 plots the profiles along one specified track (the red line in the hair region in “Lena” image in Figure 2) for different restorations, and we can see that the profile from the proposed restoration (the blue profile) has the closest match with the profile of the true “Lena” image (the red profile) of those from other restorations. Table 2 shows that, compared to other methods, the proposed LS-TPV prior restoration can lead to restored images with a higher SNR.Table 2
Signal to noise ratio (SNR) of the observed degraded images and the restored images in the simulated experiment with “Lena” image.
DegradedWienerTV priorHuber PriorForWaRDSA priorLS-TPV prior6.6515.0615.5415.5915.5015.8416.02Figure 2
Profile illustrations of different restorations.


### 4.2. Real Data Experiment

In this section, a real blurred astronomical256×256 “moon” image with intensity range [−77    170] is obtained from “http://www.iceinspace.com.au/”. Here we use a 5×5 normalized Gaussian PSF, which was tested to be effective. The parameters for different restoration approaches are listed in Table 3.Table 3
Parameter settings in real data experiment.
TV priorHuber PriorForWaRDSA priorLS-TPV priorβ=0.059β=0.045,γ=0.45Threshold = 3.5,Levels = 3β=0.032,δ=700,7×7n,11×11Nβ=0.095,α=0.5,h=325,7×7n,11×11NFigure3(b) is the restored result from Wiener deconvolution. The restored images from TV prior restoration and Huber prior restoration are, respectively, illustrated in Figures 3(c) and 3(d). Figures 3(e) and 3(f) show the restored images from the ForWaRD method and the Bayesian restoration using SA prior. As to the Bayesian restoration using the LS-TPV prior, the restored image is shown in Figure 3(g). Also, we can note that the ForWaRD method and all the Bayesian restorations (Figures 3(c)–3(g)) are free of the ring effects in Wiener deconvolution (Figure 3(b)). Also, the block artifacts observed in the restorations using Huber prior and TV prior are effectively suppressed in the restorations using SA prior and the LS-TPV prior (the lower red arrows in Figures 3(f) and 3(g)). The oversmoothing in the ForWaRD restoration can not be observed in the restorations using the SA prior and the proposed LS-TPV prior (the upper arrows in Figures 3(e), 3(f), and 3(g)). Compared to the restoration using the SA prior, the restoration using the LS-TPV prior also shows a better visual performance in preserving edges (the upper arrows in Figures 3(f) and 3(g)).(a) Original degraded “moon” image. (b) Wiener deconvolution. (c) TV prior Bayesian restoration. (d) Nonquadratic Huber prior Bayesian restoration. (e) The ForWaRD restoration in [11]. (f) The SA prior Bayesian restoration proposed in [21]. (g) The proposed LS-TPV prior Bayesian restoration (the lower arrows: block-suppressing illustration, the upper arrows: edge-preserving illustration).
(a)
Degraded image(b)
Wiener restoration(c)
TV prior restoration(d)
Huber prior restoration(e)
ForWaRD restoration [11](f)
SA prior restoration [21](g)
LS-TPV prior restoration


### 4.3. Computation Cost Comparisons

Table4 lists the recorded CPU times (in seconds) of different restorations for the experiments with simulated data and real data, respectively. Please note that the notations simu and real in Table 4 correspond to experiments with simulated data and real data, respectively. We can see that the restorations using SA prior and the LS-TPV prior are much computationally intensive than other methods. We define LS-TPV prior res1, LS-TPV prior res2, and LS-TPV prior res3 as the original serial restorations, the restorations after optimization, and the restorations after both optimization and parallelization. Table 4 shows that, with respect to the original serial version, the proposed restorations after optimization are about 3 times faster, and the restorations after both optimization and parallelization are about 50 times faster. But we can not claim that the LS-TPV prior restoration needs less computation costs than other methods since that we did not design the parallelized versions for them.Table 4
CPU times (seconds) needed for different restorations for real data experiment.
Wiener filtering:simu: 58.16, real: 64.48TV prior res:simu: 214.96, real: 228.92Huber prior res:simu: 270.52, real: 292.36ForWaRD res:simu: 2128.27, real: 2411.15SA prior res:simu: 2128.27, real: 2382.43LS-TPV prior res1: simu: 2885.61, real: 3058.07LS-TPV prior res2: simu: 910.48, real: 954.48LS-TPV prior res3: simu: 58.61, real: 64.45


### 4.4. Algorithm Monotonicity

Here, to analyze the monotonicity of the joint image/weight algorithm proposed in Section3, we calculated the posterior function energy ψβ(f,w,η∣g) and SNR with respect to the iteration numbers for the proposed restorations. The total iteration number is set to 500. As illustrated in Figure 4, we can see that both the calculated posterior function energy and the SNR increase monotonically during the whole iteration.For the proposed restoration approach, the calculated SNR (a) and the posterior function energy (b) with respect to iteration number (1–500 iterations)
(a)(b)


### 4.5. Study on NeighborhoodN and n of Different Sizes

To study the roles of the sizes ofN and n in restorations, we perform restorations using the proposed prior with N and n of different sizes in above simulation experiment with “Lena” image. To characterize local image structures, we do not set n larger than 7×7. Considering the fact that the total pixel numbers are different for different comparing patches n (49 for 7×7n, 9 for 3×3n, and 1 for 1×1n), parameter h should also change with the restorations using different patches n. The hyperparameter β and parameter h for the proposed prior are chosen by hand to produce the best stable images in terms of SNR maximization. Combinations of N and n with different sizes and the corresponding computed SNR and the recorded CPU time cost (with CUDA parallelization) are all listed in Table 5.Table 5
SNR and CPU time cost (seconds) for the simulation restorations using the proposed method.
N and nSNRCPU timesN and nSNRCPU timesN and nSNRCPU times(a)N3×3n1×115.4817.98 sec(b)N3×3n3×315.6528.59 sec(c)N3×3n7×715.7953.69 sec(d)N7×7n1×115.5421.16 sec(e)N7×7n3×315.7031.37 sec(f)N7×7n7×715.8850.91 sec(g)N11×11n1×115.5524.59 sec(h)N11×11n3×315.7238.42 sec(i)N11×11n7×716.0258.61 sec(j)N13×13n1×115.5529.54 sec(k)N13×13n3×315.7577.67 sec(l)N13×13n7×716.0795.31 secWe can see that higher SNR can be obtained when enlargingN and n from 1×1 to 11×11 and 7×7. We also note in Table 5 that more CPU time is needed for the restorations using the proposed prior with larger N and n. We should set the sizes of N and n based on the tradeoff between performance and computation cost. We find in Table 5 that no significantly quantitative SNR differences are made between 11×11N and larger 13×13N. In fact, as the iteration proceeds, more global information beyond the 11×11N will be incorporated into the regularization of the pixels.


## 4.1. Simulation Experiment

In this section, simulated experiments with256×256 “Lena” image (intensity range [−98    117]) are performed. In the simulation, a uniform 9×9 PSF is assumed, and Gaussian noise with variance = 0.1663 is added. The PSF used in this experiment is normalized to 1. Figure 1(b) is the degraded image simulated by (1.1) using the PSF with the white Gaussian noise. The parameters for different restorations are listed in Table 1. All the restored images are illustrated in Figures 1(c)–1(h), respectively.Table 1
Parameter settings in simulated data experiment.
TV priorHuber priorForWaRDSA priorLS-TPV priorβ=0.0038β=0.026,γ=0.35Threshold = 3.0,Levels = 3β=0.015,δ=900,7×7n,11×11Nβ=0.075,α=0.5,h=150,7×7n,11×11N(a) Original “Lena” image. (b) Degraded “Lena” image (SNR = 6.65) with9×9 uniform blur and additive white Gaussian noise (variance = 0.1663). (c) Wiener deconvolution. (d) TV prior Bayesian restoration. (e) Nonquadratic Huber prior Bayesian restoration. (f) The ForWaRD restoration in [11]. (g) The SA prior Bayesian restoration proposed in [21]. (h) The proposed LS-TPV prior Bayesian restoration (The left arrows: block-suppressing illustration, the right arrows: edge-preserving illustration).
(a)
Original image(b)
Degraded image(c)
Wiener restoration(d)
TV prior restoration(e)
Huber prior restoration(f)
ForWaRD restoration [11](g)
SA prior restoration [21](h)
LS-TPV prior restorationThe restored images from TV prior restoration and Huber prior restoration are, respectively, shown in Figures1(d) and 1(e). Figure 1(f) shows the restored image using the ForWaRD method. The result of Bayesian restoration using the SA prior is illustrated in Figure 1(g). For the Bayesian restoration using the proposed LS-TPV prior, the algorithm in Section 3 is used, and the corresponding restored image is shown in Figure 1(h).Through the results, we find the ForWaRD method and all the Bayesian restorations (Figures1(d)–1(h)) can overcome the ring effects in Wiener deconvolution (Figure 1(c)). The block artifacts (the left arrows), which are observed in the ForWaRD approaches and the restorations using Huber prior and TV prior, are effectively suppressed using the restorations using the SA prior and the proposed LS-TPV prior. The LS-TPV prior (Figure 1(h)) presents a further visual enhancement over the SA prior (Figure 1(g)) in block suppressing (the left arrows) and edge preserving (the right arrows). Figure 2 plots the profiles along one specified track (the red line in the hair region in “Lena” image in Figure 2) for different restorations, and we can see that the profile from the proposed restoration (the blue profile) has the closest match with the profile of the true “Lena” image (the red profile) of those from other restorations. Table 2 shows that, compared to other methods, the proposed LS-TPV prior restoration can lead to restored images with a higher SNR.Table 2
Signal to noise ratio (SNR) of the observed degraded images and the restored images in the simulated experiment with “Lena” image.
DegradedWienerTV priorHuber PriorForWaRDSA priorLS-TPV prior6.6515.0615.5415.5915.5015.8416.02Figure 2
Profile illustrations of different restorations.


## 4.2. Real Data Experiment

In this section, a real blurred astronomical256×256 “moon” image with intensity range [−77    170] is obtained from “http://www.iceinspace.com.au/”. Here we use a 5×5 normalized Gaussian PSF, which was tested to be effective. The parameters for different restoration approaches are listed in Table 3.Table 3
Parameter settings in real data experiment.
TV priorHuber PriorForWaRDSA priorLS-TPV priorβ=0.059β=0.045,γ=0.45Threshold = 3.5,Levels = 3β=0.032,δ=700,7×7n,11×11Nβ=0.095,α=0.5,h=325,7×7n,11×11NFigure3(b) is the restored result from Wiener deconvolution. The restored images from TV prior restoration and Huber prior restoration are, respectively, illustrated in Figures 3(c) and 3(d). Figures 3(e) and 3(f) show the restored images from the ForWaRD method and the Bayesian restoration using SA prior. As to the Bayesian restoration using the LS-TPV prior, the restored image is shown in Figure 3(g). Also, we can note that the ForWaRD method and all the Bayesian restorations (Figures 3(c)–3(g)) are free of the ring effects in Wiener deconvolution (Figure 3(b)). Also, the block artifacts observed in the restorations using Huber prior and TV prior are effectively suppressed in the restorations using SA prior and the LS-TPV prior (the lower red arrows in Figures 3(f) and 3(g)). The oversmoothing in the ForWaRD restoration can not be observed in the restorations using the SA prior and the proposed LS-TPV prior (the upper arrows in Figures 3(e), 3(f), and 3(g)). Compared to the restoration using the SA prior, the restoration using the LS-TPV prior also shows a better visual performance in preserving edges (the upper arrows in Figures 3(f) and 3(g)).(a) Original degraded “moon” image. (b) Wiener deconvolution. (c) TV prior Bayesian restoration. (d) Nonquadratic Huber prior Bayesian restoration. (e) The ForWaRD restoration in [11]. (f) The SA prior Bayesian restoration proposed in [21]. (g) The proposed LS-TPV prior Bayesian restoration (the lower arrows: block-suppressing illustration, the upper arrows: edge-preserving illustration).
(a)
Degraded image(b)
Wiener restoration(c)
TV prior restoration(d)
Huber prior restoration(e)
ForWaRD restoration [11](f)
SA prior restoration [21](g)
LS-TPV prior restoration


## 4.3. Computation Cost Comparisons

Table4 lists the recorded CPU times (in seconds) of different restorations for the experiments with simulated data and real data, respectively. Please note that the notations simu and real in Table 4 correspond to experiments with simulated data and real data, respectively. We can see that the restorations using SA prior and the LS-TPV prior are much computationally intensive than other methods. We define LS-TPV prior res1, LS-TPV prior res2, and LS-TPV prior res3 as the original serial restorations, the restorations after optimization, and the restorations after both optimization and parallelization. Table 4 shows that, with respect to the original serial version, the proposed restorations after optimization are about 3 times faster, and the restorations after both optimization and parallelization are about 50 times faster. But we can not claim that the LS-TPV prior restoration needs less computation costs than other methods since that we did not design the parallelized versions for them.Table 4
CPU times (seconds) needed for different restorations for real data experiment.
Wiener filtering:simu: 58.16, real: 64.48TV prior res:simu: 214.96, real: 228.92Huber prior res:simu: 270.52, real: 292.36ForWaRD res:simu: 2128.27, real: 2411.15SA prior res:simu: 2128.27, real: 2382.43LS-TPV prior res1: simu: 2885.61, real: 3058.07LS-TPV prior res2: simu: 910.48, real: 954.48LS-TPV prior res3: simu: 58.61, real: 64.45


## 4.4. Algorithm Monotonicity

Here, to analyze the monotonicity of the joint image/weight algorithm proposed in Section3, we calculated the posterior function energy ψβ(f,w,η∣g) and SNR with respect to the iteration numbers for the proposed restorations. The total iteration number is set to 500. As illustrated in Figure 4, we can see that both the calculated posterior function energy and the SNR increase monotonically during the whole iteration.For the proposed restoration approach, the calculated SNR (a) and the posterior function energy (b) with respect to iteration number (1–500 iterations)
(a)(b)


## 4.5. Study on NeighborhoodN and n of Different Sizes

To study the roles of the sizes ofN and n in restorations, we perform restorations using the proposed prior with N and n of different sizes in above simulation experiment with “Lena” image. To characterize local image structures, we do not set n larger than 7×7. Considering the fact that the total pixel numbers are different for different comparing patches n (49 for 7×7n, 9 for 3×3n, and 1 for 1×1n), parameter h should also change with the restorations using different patches n. The hyperparameter β and parameter h for the proposed prior are chosen by hand to produce the best stable images in terms of SNR maximization. Combinations of N and n with different sizes and the corresponding computed SNR and the recorded CPU time cost (with CUDA parallelization) are all listed in Table 5.Table 5
SNR and CPU time cost (seconds) for the simulation restorations using the proposed method.
N and nSNRCPU timesN and nSNRCPU timesN and nSNRCPU times(a)N3×3n1×115.4817.98 sec(b)N3×3n3×315.6528.59 sec(c)N3×3n7×715.7953.69 sec(d)N7×7n1×115.5421.16 sec(e)N7×7n3×315.7031.37 sec(f)N7×7n7×715.8850.91 sec(g)N11×11n1×115.5524.59 sec(h)N11×11n3×315.7238.42 sec(i)N11×11n7×716.0258.61 sec(j)N13×13n1×115.5529.54 sec(k)N13×13n3×315.7577.67 sec(l)N13×13n7×716.0795.31 secWe can see that higher SNR can be obtained when enlargingN and n from 1×1 to 11×11 and 7×7. We also note in Table 5 that more CPU time is needed for the restorations using the proposed prior with larger N and n. We should set the sizes of N and n based on the tradeoff between performance and computation cost. We find in Table 5 that no significantly quantitative SNR differences are made between 11×11N and larger 13×13N. In fact, as the iteration proceeds, more global information beyond the 11×11N will be incorporated into the regularization of the pixels.


## 5. Conclusions and Future Work Plan

In this paper on image restoration, we proposed an LS-TPV prior which penalizes the total distances between neighboring patches through a constrained mixture prior mode. A convergent joint image/weight updating algorithm, which estimates image and weights sequentially, is proposed to overcome the heuristical OSL weight determination for the SA prior restoration in [21]. We can see that, in addition to providing effective regularization, the proposed approach can lead to stable iteration with convergence guaranteed.The application of the proposed prior model needs a pixel-wise computation of the distances between neighboring patches over a large region, which implies high computation cost of restorations. In this paper, remarkable advance in shortening computation time is achieved by optimizing patch distance computation and replacing the original pixel-wise calculation of patch distances by CUDA framework.Further work includes applying the proposed approach in blind image deconvolutions, further justifying the proposed approach by using more image-quality measures based on human visual system, and further accelerating the computation by parallelize the PSF convolution in restorations.



---

*Source: 408241-2011-07-02.xml*