TEXT SIZE

CrossRef (0)
Immediate solution of EM algorithm for non-blind image deconvolution

Seung-Gu Kim1,a

aDepartment of Big Data Science, Sangji University, Korea
Correspondence to: 1 Department of Big Data Science, Sangji University, Sangjidae-gil 83 (Woo San-dong), Wonju-si, Gangwon-do 26339, Korea. E-mail: sgukim@sangji.ac.kr
Received January 7, 2022; Revised February 9, 2022; Accepted March 10, 2022.
Abstract
Due to the uniquely slow convergence speed of the EM algorithm, it suffers form a lot of processing time until the desired deconvolution image is obtained when the image is large. To cope with the problem, in this paper, an immediate solution of the EM algorithm is provided under the Gaussian image model. It is derived by finding the recurrent formular of the EM algorithm and then substituting the results repeatedly. In this paper, two types of immediate soultion of image deconboution by EM algorithm are provided, and both methods have been shown to work well. It is expected that it free the processing time of image deconvolution because it no longer requires an iterative process. Based on this, we can find the statistical properties of the restored image at specific iterates. We demonstrate the effectiveness of the proposed method through a simple experiment, and discuss future concerns.
Keywords : image deconvolution, EM algorithm, immediate EM solution
1. Introduction

Image deconvolution refers to the operation of reconstructing a given blur image into a sharp image. In general, the relationship is assumed as the following model

$y=h⊗x+e,$

where ⊗ denotes finite discrete convolution operator, h is a psf (point spread function) blur kernel, x is an unknown true sharp 2D image, and e is a noise in which elements are independently followed by Gaussian distribution $N(0,σe2)$.

When h is also unknown, the simultaneous estimation task of x and h is called “blind deconvolution”, whereas when h is given as a known, it is called “non-blind deconvoution”. So the latter is treated as an easier problem than the former.

However, many blind deconvolution techniques, like in Fergus et al. (2006), Levin et al. (2011), Jin et al. (2018) among others, use a non-blind deconvoution because it is necessary to estimate x under a given x in the iterative process. Therefore, it is still important to develop an efficient non-blind deconvoution technique.

For non-blind deconvoution, the most popular algorithm is RL (Richardson, 1972 ; Lucy, 1974). This technique is based on the maximum likelihood estimation by the EM algorithm under the Poisson model. However, this technique requires a lot of time to obtain the desired result. Moreover, as shown in the experiment in Section 5, it contains undesirable artifacts especially when the iteration is not sufficient.

In this paper, we provide non-blind deconvolution using the EM algorithm under the Gaussian model. However, due to the inherent slow convergence speed of the EM algorithm, the processing time is impractical when the image size is large. Therefore, in this study, an immediate restored image of the EM algorithm at any iterates is presented so that an iterative process is not required.

The following section presents terms and concepts necessary for this paper. Section 3 presents an iterative EM algorithm. In Section 4, an immediate solution of the iterative EM algorithm is derived. And based on this result, we examine the statistical properties of the restored image in the middle of iteration. In Section 5, we demonstrates the effectiveness of the proposed method through a simple experiment. Finally, in Section 6 conclusions and discussions are presented.

2. Preliminaries

Equation (1.1) can be expressed as the following matrix notation for convenience

$y=Hx+e,$

where yn×1 = {yi}, xm×1 = {xj}, and en×1 = {ei} are vectors reshaped by vectorizing 2D arrays y, x, and e in a column-stacking manner, respectively, and Hn×m = {hi j} is the corresponding convolution matrix (the blocked-Toelplitz matrix) for the filter h. Equivalently, it can be expressed as

$yi=hiT x+ei=∑j=1mhijxj+ei, i=1,…,n,$

where hi is the ith row of H, and thereafter the superscript T denotes the transpose of the matrix. Then yi independently follows $N(∑j=1mhijxj,σe2)$. Our goal is to estimate the sharp image x given a blur image y.

The unbiased least squared estimator which maximizes the log-likelihood, $L(x)∝-1/2σe2(y-Hx)T(y-Hx)$

$x^LS=(HT H)-1HT y=x+(HT H)-1HT e,$

is a good choice if the noise e is negligible. However, it gives very undesirable results even for small noises since it amplifies the noise by the term (HTH)−1HTe in systematic pattern.

To alleviate this problem, for a positive definite matrix W = {wjk}, we need to regularize the solution by using a penalized log-likelihood

$Lλ(x)∝-12(y-Hx)T(y-Hx)-λ2xT Wx,$

λ > 0 is a regularization constant. One can use W = KTK where K is the convolution matrix of a kernel k (e.g Laplacian or gradient etc.) in the way k ⊗ x ≡ Kx.

By maximizing the equation (2.4) we readily obtain the best solution as follows,

$x^λ=(HT H+λW)-1HT y.$

It is well-known that xλ is the ridge estimator when W = I, and in addition when $λ=σe2/σx2$ where $σx2$ is the variance of x, it is the Wiener deconvolution.

However, the optimal solution xλ gives an undesirable image for an inaccurately guessed λ. Therefore, the method of stopping when a desired image is obtained is often used in the iterative process of the EM algorithm.

3. Iterative EM deconvolution

A natural setting for understanding our problem as a missing/observed data structure is to look at it as follows,

$yi=∑j=1mzij,i=1,…,n,$

where zi j are missing data which follows $N(hijxj,σe2/m)$ independently, and therefore (y, z) is complete data denoting z = {zi j}. It is same as in RL deconvolution but their zi j ~ Poisson(hi jxj). The setting equation (3.1) allows us to estimate xj individually without taking into account the superimposed relation between xjs due to convolution. However, since the parameter space is extended from ℝm to ℝn×m, the convergence speed of the EM algorithm will be extremely slow.

Noting that the density p(y, z) = p(z) from relationship of equation (3.1), the (penalized) log-likelihood for complete data is

$Lλc(x)∝ -12mσe2∑i=1n∑j=1m(zij-hijxj)2-12σx2∑j=1m∑k=1mxjwjkxk∝ -12∑i=1n∑j=1m(zij2-2hijxjzij+hij2xj2)-λ2m∑j=1m∑k=1mxjwjkxk,$

where $λ=σe2/σx2$. EM algorithm at the ith iteration maximizes $Q(x;x(t))=E[Lλc(x)∣y,x(t)]$. In E-step, we need to compute two conditional expectation E[zi j|y, x(t)] and $E[zij2∣y,x(t)]$. The latter is just a necessary one for estimating λ. But since it will be used as a given tuning constant in this paper, we only compute the former in order to estimate solely x as follows.

$zij(t+1)=E(zij∣yi,xj(t))=E(zij∣xj(t))+cov(zij,yi∣xj(t))var(yi∣xj(t))[yi-E(yi∣xj(t))]=hijxj(t)+1m (yi-∑k=1mhikxk(t)),i=1,…,n, ;j=1,…,m.$

In M-step, substituting $zij(t+1)$ with the equation (3.3) in the normal equation 0 = ∂Q(x; x(t))/∂xj, we have

$0=∑i=1n(hijzij(t+1)-hij2xj)-λm∑k=1mwjkxk=xj(t)∑i=1nhij2+1m∑i=1nhij(yi-∑k=1mhikxk(t))-xj∑i=1nhij2-λm∑k=1mwjkxk,j=1,…,m,$

where $λ=σe2/σx2$ is used a user-specified tuning constant in this paper. Equivalently, equation (3.4) can be expressed in matrix form as

$0=x(t)+vHT(y-Hx(t))-x-λvWx,$

where $v-1=m∑i=1nhij2$. Note that $∑i=1nhij2$ is the sum of square of all elements in the jth column of the convolution matrix H, and they are same for all columns, thus v does not depend on j.

At hence, at (t + 1)th iteration, we obtain the EM deconvolution solution as

$x(t+1)=(I+λvW)-1(x(t)+vHT(y-Hx(t)))=Dλ-1(I-vHT H) x(t)+vDλ-1HT y,$

where Dλ = I + λvW.

As an alternative, by applying the Green (1990)’s One-Step-Late (OSL) method, one can also use λvWx(t) instead of λvWx in equation (3.5). Then we have

$x(t+1)=x(t)+vHT(y-Hx(t))-λvWx(t),=[I-v(HT H+λW)] x(t)+vHT y.$

By the essential property of EM algorithm, proceeding iteration monotonically increases the penalized log-likelihood Lλ(x) for the incomplete data y.

If the image y is very small, the matrix-typed algorithms (3.6) and (3.7) can be used as they are, but otherwise, the calculation is performed via the Fourier transform ℱ because the matrix operation time is quite burdensome. Let’s = ℱ (y), = ℱ (h), t = ℱ (xt), and = ℱ (k). Calculate

$A˜={(1-vH˜*H˜)D˜,for (3.6)1-v(H˜*H˜+λW˜),for (3.7),$$b˜={vH˜*Y˜D˜,for (3.6)vH*Y˜,for (3.7),$$X˜t+1=A˜X˜t+b˜,$

where = *, = 1 + λvW̃, and ‘*’ indicates the conjugation operator. Then, convert into the spatial domain as

$x(t+1)=|ℱ-1(X˜t+1)|.$

As mentioned previously, however, algorithm (3.6) and (3.7) have very slow convergence speed. Moreover, if the image y is large, the processing time is impractically expensive. If a very long iteration is required to obtain a desired image, the practicality of the EM algorithm will be reduced. In order to overcome this problem, in the next section, a method to obtain an immediate solution of the EM algorithm without performing iterations is suggested.

4. Immediate solution of EM algorithm

We can rewrite the equation (3.6) and (3.7) as follows,

$x(t+1)=Aλx(t)+bλ,$

saying that the current estimates x(t+1) is a linear transformation of the just preceding estimates x(t), where for the equation (3.6)

$Aλ=Dλ-1(I-vHT H),bλ=vDλ-1HT y,$

and for the equation (3.7)

$Aλ=I-v(HT H+λW),bλ=vHT y.$

For an initial estimates x(0), one can show that

$x(t+1)=Aλx(t)+bλ=Aλt+1x(0)+(Aλt+Aλt-1+⋯+I)bλ,$

where $Aλt=Aλ⋯Aλ$ that is multiplied t times, and denote $Aλ0=I$.

Note that if 0 < |Aλ| < 1, then the series $∑i=0tAλi=(I-Aλt+1)(I-Aλ)-1$, and $Aλt→0$ as t → ∞ (Searle, 1982). Thus $∑i=0∞Aλi=(I-Aλ)-1$.

Accordingly, our immediate EM solution after t iteration is

$xEM(t)=Aλtx(0)+(I-Aλt)(I-Aλ)-1 bλ$

In our case, it is straightforward to show that |Aλ| < 1 is satisfied for a sufficiently small v. For the equation (4.2), since Dλ = I + λvW is positive definite, |Aλ| = |IvHTH|/|Dλ| < |IvHTH|. Here, since vHTH positive definite, its eigenvalues are within [0, 1] for a small v, therefore |IvHTH| ≤ 1, and then |Aλ| < 1. In similar manner for the equation (4.3), since HTH + λW is also symmetric and positive definite, |Aλ| = |Iv(HTH + λW)| < 1 for a sufficiently small v. In this paper, it is recommended to choose 0 < v ≤ 1 because it worked well without divergence in various experiments. Actually, in our algorithm v can be used as a user-specified convergence step.

From (4.5), we find xEM(∞) = (IAλ)−1bλ because $Aλ∞=0$, and easily show (for both algorithm (4.2) and (4.3)) that

$(I-Aλ)-1bλ=(HT H+λW)-1HT y=x^λ$

which is just the solution (2.5). Consequently, the EM solution says that

$xEM(t)=Aλtx(0)+(I-Aλt)x^λ,$

which is a mixture of the initial estimates and the limit estimates where the proportion rate is the “weight matrix” $Aλt$ at the tth iteration. As t increases, it moves away from the initial estimates x(0) and approaches the limit estimates $x^$λ. Also, the convergency of this algorithm does not depend on initial estimates even when x(0) = 0.

It is worth investigating the bias and variance of xEM(t). For the true image x

$Bias(xEM(t))=x-AλtE(x(0))-(I-Aλt)E (x^λ)={x-E (x^λ)}+Aλt[{x-E (x(0))}-{x-E (x^λ)}]=bias of the limit estimates+Aλt×bias increments from initials to limits.$

Since the limit estimator $x^$λ is always determined irrespective of the initial estimator x(0), they are independent of each other. Therefore, the variance of xEM(t) is

$var(xEM(t))=Aλtvar(x(0))(Aλt)T+(I-Aλt)var (x^λ) (I-Aλt)T=a mixture of variance between initial and limit estimator.$

From equation (4.8) and (4.9), it can be seen that the bias and variance of the initial estimator x(0) is heavily reflected at the beginning of the algorithm, but eventually those of the limit estimator $x^$λ are reached along with At0 as t→∞.

Consider an initial estimator where the bias is greater than $x^$λ and its variance is zero, for an example, an totally flat image like x(0) = 0. Then, since the bias increments are positive, a restored image with a greater bias than $x^$λ is always provided during iteration (i.e. more smooth image). On the other hand, its variance becomes gradually larger (i.e. more sharpen image).

From these opposite properties, it can be expected to obtain a visually more satisfactory restored image than the x(0) or the $x^$λ while compromising bias and variance during iteration. This explains why many users of EM algorithms for image deconvolution adopt a strategy of stopping in the middle of iterations. However, the problem is that it takes too much processing time to get to that point using iterative EM. The proposed immediate EM solves this problem.

Meanwhile, the Fourier version of xEM(t) of equation (4.5) is

$A˜t=(A˜)tX˜t=A˜tX˜0+(1-A˜t)(1-A˜)b˜,$

where Ã, and 0 are denoted in Section 3. Then we obtain the immediate EM restoration image as xEM(t) = |ℱ −1(t)|.

It should be further noted that our xEM(t) allows the number of iterations t to be used a real numbers, but also an integer. For example, it can be seen that a restoration image xEM(2.5) provides about halfway between the 2nd and 3rd iteration. Another advantage of our method is that it allows us to examine the results of the restored image sequentially for the desired number of iterations. For example, you can get the restoration results: starting with 1000 iteration image, up to 2000 iteration image with increasing by 100 iterations.

5. Experiments

In this section, we use real images to show that immediate EM works properly. In Figure 1, a true image x, a blur kernel h, and the observed blur image y are provided. The x is a 512 × 512 RGB image where pixels are scaled between [0, 1]. The h a 49 × 49 uniform blur kernel with hight 0.0016 and its sum is 1. And the y is an observed image obtained by adding Gaussian noise with a standard deviation of σe = 0.001 to the result of x ⊗ h ≡ Hx.

Our goal is to reconstruct the distorted image y under a given H (that is, the matrix version of the blur filter h) to be close to the true image x. Reconstruction was performed separately for each R, G, and B channel, and then the results were merged.

In this experiment, the penalty term in Equation (2.5) will be used a simple “roughness penalty” for the purpose of removing noise

$xT Wx=∑j=1m(xj-x¯∂j)2=‖k⊗x‖2,$

where χ̄j denotes the average of neighbors at the pixel j of x, and j denotes the 1st order neighborhood,

$k=14 (0-10-14-10-10).$

And the initial estimates of immediate EM is x(0) = 0.

Figure 2 demonstrates the log-likelihood of the EM algorithm over iteration, where it monotonically increases in all cases. In particular, the original EM of equation (3.6) and OSL-EM of equation (3.7) showed very insignificant difference. And the larger the value of v, the faster the convergence result at the beginning of iteration, but a similar convergence speed near the convergence point.

The first row of Figure 3 provides the RL deconvolution results for 50, 1500, and 2500 iterations (using function richardson-lucy in python module skimage). And similar results were obtained for more iterations. Although it worked well for deblur and denoising, “ringing artifacts”, a problem unique to this technique, appear. First of all, the processing time was long, taking 342.71 seconds for 2500 iterations.

In the second row of Figure 3, Wiener deconvolution results are provided (using function wiener in python module skimage). The first and third ones are reconstructed using λ = 0.1 and λ = 0.17 in common for R, G, and B channel images, respectively. And the second one is with exact but unknown channel-wise different λ,

$(λR,λG,λB)=σe2(σx(R)2,σx(G)2,σx(B)2)=0.0012(0.0378,0.0440,0.0183)≈(2.65, 2.27, 5.46)×10-5$

which shows better result. If the correct λ is not specified, it gives undesirable results.

In the third row of Figure 3, our immediate EM deconvolution results are provided using λ = 0.17 that showed a bad result with Wiener deconvolution. The 1st and 3rd image are xEM(1) and xEM(∞) where they are represent a high biased and high variance, respectively. A good restored image (the 2nd) was obtained by choosing tentatively t = 50000, although it not reached at maximum likelihood. As expected, it confirms that subjectively and visually satisfactory images can exist during the EM algorithm. However, since the processing time for one iteration is about 0.34 seconds, it is expected that it will take about 4.7 hours to obtain the image x(50000) with the iterative EM algorithm. On the other hand, the proposed immediate EM xλ(50000) was obtained in only 0.26 seconds.

The PC used in the experiment in this chapter is Intel(R) i7-8700 CPU (3.20GHz), and the program is implemented in python.

6. Conclusions and discussion

In this paper, a non-blind deconvolution by iterative EM algorithm under Gaussian model was provided. Then an immediate EM without iterative processing that provides the same result of iterative EM was developed. This significantly reduces the processing time of image deconvolution with EM algorithm, thereby increasing the practicality. The result makes it possible to reveal the statistical properties of bias and variance for the restored image. This has hitherto been vaguely recognized but not explicitly known.

The Richardson-Lucy method is an EM algorithm based on Poisson model, and it is necessary to study whether an immediate solution can be obtained for this method as well.

In this study, very limited experiments were provided and an excessively simple penalty term was used. The use of various initial images and edge preserving penaty is expected to provide more interesting results.

Figures
Fig. 1. True image (Left). Blur filter (Middle), where all white area has a pixel value of 0.0016. Observed image (Right) where is convolved with a blur fiter and then added gaussian noise with = 0.001.
Fig. 2. Penalized log-likelihood () over the iteration. Original EM (Blue), OSL?EM (Red). ○ : = 0.2, △ : = 0.5, and + : = 1.0.
Fig. 3. The 1st row displays LR deconvolution: iter=50(left), iter=150(middle), and iter=2500(right). The 2nd row displays Wiener deconvolution: = 0.1(left) in common, accurate , and = 0.1(right) in common. The 3rd row displays the immediate EM deconvolution: EM (1)(left), EM (50, 000)(middle), and the limit estimates EM (∞) = (right).
References
1. Fergus R, Singh B, Hertzmann A, Roweis ST, and Freeman WT (2006). Removing camera shake from a single photograph. ACM Trans Graph, 25, 787-794.
2. Green PJ (1990). On the use of EM algorithm for penalized likelihood estimation. Journal of Royal Statistical Society B, 52, 443-452.
3. Levin A, Weiss Y, Durand F, and Freeman WT (2011). Efficient marginal likelihood optimization in blind deconvolution. CVPR.
4. Lucy L (1974). Bayesian-based iterative method of image restoration. Journal of Astronomy, 79, 745-754.
5. Jin M, Roth S, and Favaro P (2018). Normalized blind deconvolution, Computer Vision ? ECCV.
6. Richardson W (1972). Bayesian-based iterative method of image restoration. Journal of the Optical Society of America A, 62, 55-59.
7. Searle SR (1982). Matrix Algebra Usefull for Statistics, Wiley.