Cause Analysis and Removal of Boundary Artifacts in Image Deconvolution

Abstract

In this paper, we conducted a cause analysis on boundary artifacts in image deconvolution. Results of the cause analysis show that boundary artifacts are caused not only by a misuse of boundary conditions but also by no use of the normalized backprojection. Results also showed that the correct use of boundary conditions does not necessarily remove boundary artifacts. Based on these observations, we suggest not to use any specific boundary conditions and to use the normalized backprojector for boundary artifact-free image deconvolution.

1. INTRODUCTION

The image deconvolution is a basic problem in image processing and has many applications. Examples include restoration of astronomical images, microscopy, and medical imaging, etc [1]. Many methods have been developed for the image deconvolution. For example, we have the classical Wiener filtering [2], Bayesian based approaches [3-5], wavelet deblurring [6,7], variational methods [8,9,10], etc.

The image deconvolution, however, often produces undesirable artifacts in deconvolved images. The paper [11] classifies obstacles of the image deconvolution into following four categories; noise, insufficient deconvolution, boundary artifacts, and incorrect blurring model. In this paper our main concern is the removal of boundary artifacts.

The boundary artifacts in image deconvolution problem are caused by various factors, depending on deconvolution methods. For instance, classical frequency domain filtering methods, which utilize the decoupling property of the circular convolution via fast Fourier transform (FFT), severely suffers from boundary artifacts, due to the mismatch in blurring models; the frequency domain filtering methods assumes a circular convolution as the blurring model, instead of a truncated convolution, which fits more accurately blurring models in many applications.

To remove boundary artifacts, various boundary condition methods have been proposed along with iterative deconvolution methods. The boundary condition (BC) method imposes certain condition on pixels across boundaries. Among them, periodic, reflective, and anti-reflective BCs have been most popular [12-15]. BC methods, however, often do not remove boundary artifacts effectively in case when imposed conditions are greatly mismatched to characteristic of images to be recovered.

In this paper, to investigate what really causes boundary artifacts, we conducted various simulation studies. Through a careful analysis on source of boundary artifacts, we suggest not to impose any specific BCs to avoid imposing mismatched BCs to characteristics of images to be restored. This result was preluded in the work of Calvetti et. al. [16]. The main result of the cause analysis is that no use of BC is not sufficient, while the use of the normalized backprojector which is designed to compensate for the non-uniformity in contributions from image pixels to the observation along with no use of BC removes boundary artifacts effectively. Simulation results in this paper will support this claim.

 

2. BACKGROUNDS

2.1 Problem of Image Deblurring

The image deconvolution problem is to find the true image f from the observed image g which is blurred and noised as

where P is the linear transform that determines the blurring process and n represents a mean zero Gaussian noise. In this paper, we assume that the blurring by P is space-invariant so that there exists a point spread function (PSF) k such that Pf = k ∗ f the convolution of k and f. The exact definition of k ∗ f will be given in the next section.

In (1), to include every pixel of f that affects the observed image g through blurring, it is natural to assume that the domain Ω of the true image f is larger than the domain ∧ of the observed image g. The blurring process makes some near-boundary pixels of the observed image g to be influenced by ‘unseen pixels’ (pixels in Ω - ∧). It is often regarded that the existence of unseen pixels causes boundary artifacts [12].

In this paper, we will use following notations and terminologies.

where Sk = {v|kv ≠ 0}. The blurring transform P will be called projector throughout this paper.

Throughout this paper, we assume that the PSF k satisfies following conditions:

2.2 Iterative Methods

Landweber Method (LM)

For g = Pf, the Landweber iteration takes

Here, β the step length must have the value of

where σmax is the largest singular value of P. For details and applications of Landweber method, see, e.g., [17].

Richardson-Lucy (RL)

For g = Pf, starting from f0 = 1Ω, where 1Ω is the all-one image defined on Ω, the Richardson-Lucy method takes

fn＋1 = fn.∗Ptsn./Pt1∧, sn = g./Pfn,

where, .∗ and ./ are pixel-by-pixel multiplication and division between two images. and 1∧ is the all-one images defined on ∧. For details and applications of Richardson-Lucy method, see, e.g., [4,5].

Conjugate Gradient Least Square (CGLS)

The conjugate gradient method is designed for the linear system governed by positive definite matrix [18]. For (1), we consider the Tikhonov regularization method formulated by

minf(║g - Pf║2＋λ║f║22),

where λ>0 is a regularization parameter and ║・║2 is the L2-norm. This variational problem leads to the following normal equation:

Since the governing matrix PtP ＋ λ1Ω in (5) is positive definite, the conjugate gradient method is applicable to (5). We will denote this method by Conjugate Gradient Least Square (CGLS).

2.3 Boundary Conditions

In this section, we will explain various BCs such as periodic, reflective, and anti-reflective boundary conditions.

To simplify the presentation, we assume that

for some positive integers N1, N2 and

for some positive integers L1, L2, M1, M2. In this case,

where ν = 1,2

We divide the true image f into 9 parts as follows

where fc = (fi1,fi2), 0 ≤ iν< Nν (ν = 1,2) represents the image part defined on ∧; the others represent the image part defined on Ω-∧, the set of unseen image pixels across the boundary of ∧. Boundary condition method impose certain restriction on fnw, fn, fne, fw, fe, fs, and fse in (9).

Periodic BC (PBC)

In 2-Dimensional, the i-th row of the periodic boundary condition imposed images is

and similarly for the column.

Reflective BC (RBC)

The reflective boundary condition, in 2-dimensional, imposed image is as follows

and similarly for the column.

Anti-Reflective BC (ABC)

The anti-reflective boundary condition, in 2-dimensional, imposed image is as follows

and similarly for the column.

Any set of BCs introduce an extension operator E: L2(∧)→L2(Ω) such that

where and f∗nw, f∗n, f∗ne, f∗w, f∗e, f∗sw, f∗s, and f∗se represent parts of the image imposed by boundary condition. By using extension operator E associated with PBC, RBC, and ABC can be defined by (10), (11), and (12), respectively. BC methods use the extension operator E to modify (1) to

Iterative deblurring methods described in section 2.3 can be applied to (14) by replacing f with fc and P with PE.

One of the main advantage of BC methods is the computational efficiency. For example, PBC, RBC, and ABC based models can be diagonalizable by FFT, discrete cosine transform (DCT), and sine transform (DST), respectively. For details, see. e.g., [13-15].

Free BC (FBC)

Instead of imposing specific BCs, works in [15,18,19] suggest not to use them. In [19], Calvetti et. al. claim that “In an Aristotelian approach to knowledge, when it is not known a priori which BCs should be chosen, by admitting our lack of information it is possible to let the data itself determine them.” In [19], this no use of BC approach is called free BC (FBC). Thus, FBC does not modify the deconvolution problem (1).

We will use the notation, 'iterative method'-'BC', to describe the iterative deblurring method and the boundary condition. For example, CGLSRBC and CGLS-ABC are CGLS iteration with reflective BC and anti-reflective BC, respectively.

 

3. PROPOSED METHOD

3.1 Cause Analysi

To investigate the source of boundary artifacts, we conducted simulations with test images in Fig. 1. Fig. 1 shows images satisfying various BCs ((a) PBC, (b) RBC, (c) ABC, (d) no BC) on 10 pixel rows on the upper boundary and 10 pixel columns on the left boundary. In this paper, we used 11×11 diagonal PSF kD whose diagonal elements are {30,29,...,20} and the first diagonal element is k0,0. Notice that the location of BC-imposed image pixels in Fig. 1 is determined by the support of the PSF kD.

Fig. 1.Test images. (a) (upper-left)‘clock’satisfying PBC, (b) (upper-right)‘peppers’satisfying RBC, (c) (lower-left)‘girl’satisfying ABC, (d) (lowerright)‘ bridge’satisfying no BCs. All images are of size 250×250.

Fig. 2 shows blurred versions of images in Fig. 1 by the 11×11 diagonal PSF kD. In this simulation, we did not add noise to exclude noise-related artifacts in deconvolved images for a better visual comparison on boundary artifact removal.

Fig. 2.Blurred images by the 11×11 diagonal PSF kD from images in Fig. 1. All blurred images are of size 240×240.

Fig. 3 shows deconvolved images by LM from images in Fig. 2. In Fig. 3, we presented restored images on ∧, if necessary (Fig. 3(d)), by cropping out image pixels on Ω - ∧. Here we note that LM-PBC (Fig. 3)(a), LM-RBC (Fig. 3)(b), and LM-ABC (Fig. 3(c)) recover images fc defined on ∧, while LM-FBC recovers f defined on Ω. In Fig. 3(c)(lower-left) satisfying ABC, the boundary artifacts were very severe even though the image was obtained by only 10 iterations. Other images in Fig. 3 were obtained by 1000 iterations. In this simulation, we used a fixed parameter β = 0.25 which is the step length of LM iterations in (6) for all pixels. Smaller step lengths β can reduce the artifacts in Fig. 3(c), but, at the same time, it decreases the convergence rate.

Fig. 3.Deconvolved images by LM from blurred images in Fig. 2 by using (a) (upper-left) PBC, (b) (upper-right) RBC, (c) (lower-left) ABC, and (d) (lower-right) FBC.

In the simulation for Fig. 3, original images to be recovered are assumed to satisfy imposed BCs, as seen in Fig. 1. Thus, in the simulation for Fig. 3, there were no unseen pixels across upper and left boundaries in observed images in Fig. 2. Some of restored images ((c), (d)) in Fig. 3, however, exhibited boundary artifacts. These simulation results show that correctly imposed BC does not guarantee boundary artifact removal in LM.

Fig. 4 shows deconvolved images by RL from images in Fig. 2. Again, in Fig. 3, we presented restored images only on ∧.

Fig. 4.Deconvolved images by RL from blurred images in Fig. 2 by using (a) (upper-left) PBC, (b) (upper-right) RBC, (c) (lower-left) ABC, and (d) (lower-right) FBC.

In the simulation for Fig. 4, original images to be recovered are also assumed to satisfy imposed BCs, as in the simulation for Fig. 3. Thus, again, in the simulation for Fig. 4, there were no unseen pixels across upper and left boundaries in observed images in Fig. 2. Fig. 4(c), however, exhibited severe boundary artifacts.

Restored images by RL showed different implications for ABC (Fig. 4(c)) as compared with PBC (Fig. 4(a)), RBC (Fig. 4(b)), and FBC (Fig. 4(d)). Boundary artifacts in RL-ABC (Fig. 4(c)) can be explained by following argument. Notice that RL-ABC takes the iteration

where E is the extension transform associated with ABC. It is well known that RL-type iterations does not provide accurate estimation for pixels v where the denominator (EtPt1∧)v = 0. The characteristic of the extension transform E associated with ABC, however, makes (EtPt1∧)v = 0 for some pixel v. This result shows that RL-ABC is not a good approach for the image deconvolution.

As compared with LM-FBC (Fig. 3(d)), RLFBC (Fig. 4(d)) did not show boundary artifacts. Since FBC means no use of BCs, this result show that the boundary artifact removal can be archived without using correct BCs.

As seen in Fig. 3(d), FBC alone, however, does not remove boundary artifacts [20]. Despite the failure of FBC in removing boundary artifacts, we suggest FBC as the first step to avoid boundary artifacts caused by the use of mismatched BCs, by noting that boundary artifacts in FBC can be easily removable by using the normalized backprojector. The exact definition of the normalized backprojector and detailed explanation will be given in the next section.

3.2 Normalized Backprojector

Simulation results in Section 3.1 show that the use of correct BC does not guarantee boundary artifact removal in LM and RL. Moreover, the result in Fig. 4(d) by RL-FBC shows that boundary artifacts can be removed without using correct BCs. The comparison of Fig. 3(d) by LM-FBC with Fig. 4(d) by RL-FBC indicates that the boundary artifact can be removed not by imposing correct BC but by using special characteristic of RL that is not in LM.

To find out what special characteristic of RL leads to boundary artifact removal, we reformulate the Landweber method (4) in the following generalized form [21]:

where Bn = diag{β1(n),β2(n),…,βΩ(n)} with positive elements controlling step size of each iteration and positive definite matrices Mn. Notice that RL also can be written in an addition form (16) with

This shows that the main difference between LM and RL is the use of the normalized term Pt1∧ in the denominator. We call 1Ω./Pt1∧.∗Pt normalized backprojector (NBP) by noting that it backprojects the all-one image 1∧ on ∧ to the all-one image 1Ω on Ω.

We denote the modified LM that uses the normalized backprojector 1Ω./Pt1∧.∗Pt instead of the standard backprojector Pt with FBC by LMNFBC. We also use LM-NPBC, LM-NRBC, LMNABC to refer NBP-based LM with PBC, RBC, and ABC, respectively. With the same argument, we define CGLS-NFBC. Here we note that RLNPBC, RL-NRBC, RL-NABC, RL-NFBC are same as RL-PBC, RL-RBC, RL-ABC, RL-FBC, respectively.

To check whether the normalized backprojector 1Ω./Pt1∧.∗Pt can remove boundary artifacts in LM, we conducted image deconvolution simulation by NBP based LM.

Fig. 5 shows deconvolved images by NBP based LM from images in Fig. 2. As in simulations for Fig. 3 and 4, in the simulation for Fig. 5, original images to be recovered are assumed to satisfy imposed BCs, as seen in Fig. 1, and hence there were no unseen pixels across upper and left boundaries in observed images in Fig. 2.

Fig. 5.Deconvolved images from blurred images in Fig. 2 by using (a) (upper-left) LM-NPBC, (b) (upper-right) LM-NRBC, (c) (lower-left) LMNABC, and (d) (lower-right) LM-NFBC.

Fig. 5(a) by LM-NPBC was exactly identical to Fig. 3(a) by LM-PBC. This phenomenon can be explained by the following argument. Notice that the standard backprojector (PE)t of LM-PBC, where E is the extension transform associated with PBC, is already normalized, i.e., (PE)t1∧ = 1Ω because of conditions (3) imposed on the PSF k.

Fig. 5(b) by LM-NRBC showed a slight improvement over Fig. 3(b) by LM-RBC.

Fig. 5(c) by LM-NABC exhibited severe boundary artifacts as in Fig. 4(c) by RL-ABC. Boundary artifacts in LM-NABC (Fig. 5(c)) can be explained by the same argument used for Fig. 4(c) by RL-ABC, since both methods use EtPt1∧ (here E is the extension operator associated with ABC. Notice that this term can be 0 for some pixels) as the normalized term in their backprojectors. This result shows that ABC is not good for NBP.

Fig. 5(d) by LM-NFBC showed a clear improvement over Fig. 3(d) by LM-FBC. All boundary artifacts were removed. Based on these observations, we suggest not to use any specific BCs and to use the normalized backprojector for boundary artifact removal. In other words, we suggest to use iterative methods based on NFBC.

3.3 Implementation

To explain how FBC is implemented, we will use commands of Matlab [1]. Before we do this, let us recall that if Ω = RM×M, ∧ = RN×N, and k is a L×L sized PSF, then M = N＋L- 1.

As mentioned earlier, FBC does not use an extension operator E. Thus, the implementation of FBC depends on how the projection (the multiplication of the matrix P and images on L2(Ω)) and the backprojection (the multiplication of the transpose matrix Pt and images on L2(∧)) are computed. For these computations, we first compute optical transfer functions K and K(r) of the PSF k and its reversely ordered PSF k(r) (i.e., = k_i1,_i2) by using following commands:

˃˃ Z = Z(r) = zeros(M,M);

˃˃ Z(1 : L, 1: L) = k;

˃˃ Z(r)(1 : L, 1: L) = k;

˃˃ K = fft2(Z);

˃˃ K(r) = fft2(Z(r));

With these optical transfer functions K and K(r), we can compute y = Pf and x = Ptg respectively by

˃˃ Y = real(ifft2(K.∗fft2(f)));

˃˃ y = Y(L: M, L: M);

(taking-out)

and

˃˃ G = zeros(M,M);

˃˃ G(1: N, 1: N) = g;

(zero-padding)

˃˃ x = real(ifft2(K(r).∗fft2(G)));

In case when PBC is used, we uses the extension operator E to change the problem g = Pf＋n (1) to (14), where For computations of projection and backprojection associated with PBC, we first compute PBC-related optical transfer functions Kp and Kp(r) of k and k(r), respectively, by using following commands:

˃˃ Zp = Zp(r) = zeros(N,N);

˃˃ Zp(1 : L, 1: L) = k;

˃˃ Zp(r)(1 : L, 1: L) = k(r);

˃˃ Kp = fft2(Zp);

˃˃ Kp(r) = fft2(Zp(r));

Here we note that the main difference between K, K(r) (in FBC) and Kp, Kp(r) (in PBC) is the size of optical transfer function; K and K(r) are of size M×M, while Kp and Kp(r) are of size N×N.

With these optical transfer functions Kp and Kp(r), we can compute and x = Ptg respectively by

˃˃ y = real(ifft2(Kp.∗fft2(fc)));

and

˃˃ x = real(ifft2(Kp(r).∗fft2(g)));

In case when RBC or ABC are used, we can have similar computational efficiencies to PBC-related one by replacing 'fft2' and 'ifft2' with 'dct2' and 'idct2' for RBC and 'dst2' and 'idst2' for ABC, where 'dct'and 'dst' represent DCT and DST, respectively.

The computational efficiency comparison between FBC and PBC clearly shows that FBC can be implemented with a fast algorithm, which requires only few extra computations for 'bigger sized optical transfer functions', 'taking-out' and 'zero-padding' than PBC based algorithms. Simulation results in this paper will show that those extra computations in FBC pay off by removing boundary artifacts more efficiently than PBC, RBC, and ABC.

 

4. SIMULATION STUDIES

We conducted simulations to test the performance of the proposed method NFBC. Fig. 6 shows test images (a) ‘moon’ and (b) ‘boat’. Fig. 7 (a) and (b) show blurred and noised images by 11×11 Gaussian PSF k and 11×11 diagonal PSF kD, respectively. We added mean zero Gaussian noises to blurred images. Here standard deviations of added noises were set to be 0.5% of means of blurred images

Fig. 6.Test images. (a)‘moon’, (b)‘boat’. The size of test images is 250×250.

Fig. 7.Blurred, by (a) k and (b) kD, and noised images. The noise term was generated by mean zero Gaussian noise with the standard deviation = 0.5% of the mean of the blurred image. Here, we used 11×11 Gaussian PSF kG with σ = 2.5 pixel unit length and 11×11 diagonal PSF kD which was explained in section 3.1.

In simulations, we chose the image that had the smallest relative square error (RSE) in 500 iterations as the deconvolved image by the tested method. Here the RSE is defined by

where, and fi1,i2 are pixel values of the deconvolved image and the true image, respectively.

Fig. 8 shows deconvoloved images by CGLSPBC. Results showed that boundary artifacts in the Gaussian deconvolution (Fig. 8(a)) and more severe ones in the diagonal deconvolution (Fig. 8(b)). The relative differences in boundary artifacts are caused by facts that the Gaussian PSF kG decays more rapidly than the diagonal PSF kD and the ‘moon’ image has boundaries which are relatively more matched to PBC than the ‘boat’ image.

Fig. 8.Deconvolved images by CGLS-PBC from (a)(left) Fig. 7(a) and (b)(right) Fig. 7(b). The image in (a) was obtained at 7 iterations with RSE = 0.52% and the image in (b) was obtained at 5 iterations with RSE = 3.70%.

CGLS-PBC attains the smallest RSE at 7 iterations with RSE = 0.52% for Fig. 8(a) and at 5 iterations with RSE = 3.70% for Fig. 8(b). These results indicate that boundary artifacts caused by the mismatch of PBC to the image characteristic prevents CGLS-PBC from attaining smaller RSE with longer iterations.

Fig. 9 shows deconvoloved images by CGLSRBC. Results showed that severe boundary artifacts in the diagonal deconvolution in Fig. 9(b), while boundary artifacts in the Gaussian deconvolution in Fig. 9(a) were hardly noticeable. The comparison between Fig. 8(a) by CGLS-PBC and Fig. 9(a) by CGLS-RBC indicates that RBC is relatively more matched to the boundary characteristic of the ‘moon’ image than PBC. As seen in Fig. 9(b) by CGLS-RBC, the mismatch between the boundary characteristic of the ‘boat’ image and RBC caused severe boundary artifacts.

Fig. 9.Deconvolved images by CGLS-RBC from (a)(left) Fig. 7(a), and (b)(right) Fig. 7(b). The image in (a) was obtained at 29 iterations with RSE = 0.42% and the image in (b) was obtained at 17 iterations with RSE = 1.47%.

CGLS-RBC attains the smallest RSE at 29 iterations with RSE = 0.42% for Fig. 9(a) and at 17 iterations with RSE = 1.47% for Fig. 8(b). The comparison of these results with those in Fig. 8 also indicates that the mismatch of RBC is less severe than that of PBC.

Fig. 10 shows deconvolved images by CGLSABC. Results showed that both the Gaussian deconvolution (a) and the diagonal deconvolution (b) suffer from boundary artifacts.

Fig. 10.Deconvolved images by CGLS-ABC from (a) (left) Fig. 7(a) and (b)(right) Fig. 7(b). The image in (a) was obtained at 486 iterations with RSE = 0.81% and the image in (b) was obtained at 124 iterations with RSE = 1.10%.

CGLS-ABC attains the smallest RSE at 486 iterations with RSE = 0.81% for Fig. 10(a) and at 124 iterations with RSE = 1.10% for Fig. 10(b). The comparison of these results with those in Fig. 8 and 9 indicates that ABC decreases the convergence rate.

Fig. 11 shows deconvoloved images by CGLSFBC (the first row) and CGLS-NFBC (the second row). Results using FBC only showed that both the Gaussian deconvolution (Fig. 11(a)) and the diagonal deconvolution (Fig. 11(b)) suffer from boundary artifacts. CGLS-FBC attains the smallest RSE at 68 iterations with RSE = 0.56% for Fig. 11(a) and at 56 iterations with RSE = 0.94% for Fig. 11(b).

Fig. 11.Deconvolved images by CGLS-FBC and CGLSNFBC from (a and c) Fig. 7(a) and (b and d) Fig. 7(b). The image in (a)(upper-left) was obtained at 68 iterations with RSE = 0.56% by CGLS-FBC, the image in (b)(upper-right) was obtained at 56 iterations with RSE = 0.94% by CGLS-FBC, the image in (c)(lower-left) was obtained at 54 iterations with RSE = 0.41% by CGLS-NFBC, and the image in (d)(lower-right) was obtained at 61 iterations with RSE = 0.47% by CGLS-NFBC,.

Fig. 11(c) and 11(d) show deconvoloved images by CGLS-NFBC. Boundary artifacts were effectively removed both in the Gaussian deconvolution (Fig. 11(c)) and the diagonal deconvolution (Fig. 11(d)). CGLS-NFBC attains the smallest RSE at 54 iterations with RSE = 0.41% for Fig. 11(c) and at 61 iterations with RSE = 0.47% for Fig. 11(d). Thus, CGLS-NFBC achieved the smallest RSE in this simulation.

These results show that the use of NFBC removes the chance of using mismatched BCs, as seen in Fig. 8(b) and 9(b), and efficiently removes boundary artifacts by using NBP. The comparison between Fig. 11(a) and 11(b) by FBC only and Fig. 11(c) and 11(d) by FBC+NBP supports the latter claim.

 

5. CONCLUSION AND DISCUSSION

In this paper, based on the cause analysis on boundary artifacts in image deconvolution, we suggest not to use any specific BCs and to use the normalized backprojector for the boundary artifact removal. By not using any specific BCs, the proposed method is immune to impose mismatched BCs on unseen pixels. By using the normalized backprojector, the proposed method can remove boundary artifacts effectively. Simulation results supported described claims by showing that the proposed method removes boundary artifacts more effectively than PBC, RBC, and ABC, with only few extra computations than PBC based methods.