- Email: [email protected]
Image restoration from a downsampled image by using the DCT
ARTICLE IN PRESS
Signal Processing 87 (2007) 2370–2380 www.elsevier.com/locate/sigpro
Image restoration from a downsampled image by using the DCT Yoshinori Abe, Youji Iiguni Department of Systems Innovation, Graduate School of Engineering Science, Osaka University, Toyonaka 560-8531, Japan Received 5 September 2006; received in revised form 18 January 2007; accepted 20 March 2007 Available online 30 March 2007
Abstract The high-resolution (HR) image restoration from a downsampled low-resolution (LR) image by using the discrete cosine transform (DCT) is proposed. The downsampling process is modeled in matrix form and the DCT is applied to the downsampling matrix to obtain the sparse matrix representation. The restored HR image is then efficiently computed from the LR image by using the sparsity. The distortion caused by the observation noise is theoretically analyzed, and the power of distortion is shown to be closely equal to the variance of the observation noise. Computer simulations show that the proposed method is superior to the cubic spline interpolation in the HR image restoration performance as far as the amount of the additive noise is small. r 2007 Elsevier B.V. All rights reserved. Keywords: High resolution image restoration; Discrete cosine transform; Lagrange multiplier method; Similar transformation
1. Introduction High-resolution (HR) image restorations, which are acquiring a HR image from downsampled lowresolution (LR) images, attract much attention [1–3]. In the HR image restoration methods, the downsampling process, which expresses the relation between the pixel values of the LR and HR images, is mathematically modeled. When using a simple charge coupled device (CCD) image capturing model, we can express the pixel value of the LR image as the sum of corresponding pixel values of the HR image [4,5]. Then the HR image restoration results in an inverse problem of the downsampling process. Unfortunately, the inverse problem beCorresponding author. Tel./fax: +81 6 6850 6377.
E-mail addresses: [email protected] (Y. Abe), [email protected] (Y. Iiguni).
comes ill defined when only one LR image can be used to restore the HR image. Nevertheless, if the LR image is under the Nyquist condition, the sinc interpolation can reconstruct the HR image perfectly from the LR image. However, it is not suitable for practical use because of the large computational cost. Other image interpolation methods that approximate the sinc interpolation with less complexity, such as the cubic spline interpolation [6–8], and image resizing methods in the frequency domain [9,10], have been proposed. However, the actual LR image does not necessarily satisfy the Nyquist condition, and therefore the restoration performance is low. To achieve a good restoration performance, HR image restoration methods which make the illdefined restoration problem well defined have been derived [4,11–13] from the analogy with the blurred image restoration method [14,15]. In the blurred
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2007.03.010
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
image restoration, a smoothness constraint such as the square sum of the discrete Laplacian of the restored image is imposed on the restored image to uniquely restore the blurred image [14]. Similarly, in the HR image restoration, the smoothness constraint is imposed on the restored HR image to determine it uniquely from the downsampled LR image [11]. However, the direct solution of the HR image restoration problem requires a large computation time, because a large-scale linear equation needs to be solved. A block segmentation method [16] or the Kronecker product method [13] has been proposed for the fast computation, but they decrease the restoration performance. In this paper, we derive a fast computation method for the HR image restoration by using the discrete cosine transform (DCT). In the blurred image restoration, Hunt has developed a fast computation method with the fast Fourier transform (FFT) [17]. The main idea of this method is the diagonalization of the blurring matrix by the discrete Fourier transform (DFT) matrix. Unfortunately, the fast computation method cannot be straightly applied to the current problem because the downsampling matrix is not square. We therefore model the downsampling process in the absence of the observation noise, perform the DCT to obtain a sparse representation, and then derive the analytical solution to the restoration problem. Since we can express the analytical solution in a scalar form, we can quickly generate the HR image from the LR image by using the DCT. We show through computer simulations that the restoration capability of the proposed method is superior to that of the cubic spline interpolation method. The proposed method considers a small amount of observation noises, such as quantization errors, in the downsampling process. We thus neglect the observation noise, and derive the analytical solution in the absence of noise. However, there are observation noises in practice. We theoretically analyze the distortion caused by the observation noise, and show that the power of distortion is closely equal to the variance of the observation noise. We also show through computer simulations that the proposed method is still superior to the cubic spline interpolation in the restoration performance if the amount of the observation noise is small. Finally, we combine the noise reduction method [18] with the proposed method to decrease the noise effect.
2371
2. Mathematical formulation of the image restoration 2.1. Observation model Let the ð2M 2MÞ HR image be fxij g and the ðM MÞ LR image be fykl g. The HR and LR images are obtained from the same scene captured by the different resolution cameras at the same position. The array structures of the photodetectors of the cameras are described in Fig. 1. The density of the photodetectors of the HR image is 4 times higher than that of the LR image. Then we can express the relation between HR and LR images as y ¼ Hx þ v,
(1)
where x, y, and v are lexicographically ordered vectors consisting of the HR image fxij g, the LR image fykl g, and the additive noise with mean zero and variance 2 , respectively. Here we assume that the additive noise is very small. The ðM 2 4M 2 Þ matrix H represents the filtering and down-sampling process and has the form H ¼ H s H s,
(2)
where ‘‘’’ denotes the Kronecker product and H s is the ðM 2MÞ matrix defined by 1 01 1 0 2 2 C B 1 1 C B C B 2 2 C B ð3Þ Hs ¼ B C: .. C B . C B A @ 1 1 0 2 2
Fig. 1. Array structure of photodetectors: small squares represent the photodetectors for HR image and hatched squares represent the ones for LR image.
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2372
2.2. HR image restoration problem
where the ð4M 2 4M 2 Þ matrix Q is defined by
We shall consider the problem of restoring the HR image x from the LR image y. Since the problem is ill defined, we impose a smoothness constraint on x to uniquely determine x. The constraint we imposed is the square sum of the discrete Laplacian of x. The discrete Laplacian of each xij is defined by xi1;j þ xi;j1 4xij þ xi;jþ1 þ xiþ1;j . We here suppose that xij s at image boundaries i ¼ 0, i ¼ 2M 1, j ¼ 0, or j ¼ 2M 1 are under the Neumann boundary condition, that is,
Q ¼ PT P
x1;j ¼ x0;j ,
y ¼ Hx.
(7)
Under the condition (7) we shall estimate the HR image x from the LR image y by minimizing the square sum of the discrete Laplacian xT Qx. The estimation problem is formulated as x
subject to y ¼ Hx.
(8)
We solve the problem (8) by using the Lagrange multiplier method. Let the Lagrangian L be
xi;2M ¼ xi;2M1 . The boundary condition is illustrated in Fig. 2. The lexicographically ordered vector consisting of the discrete Laplacians of fxij g is then expressed as [19] Px ¼ ðI Ps þ Ps IÞx,
(4)
where I is the ð2M 2MÞ identity matrix, and Ps is the ð2M 2MÞ matrix defined by 0
0
ð6Þ
When the additive noise v is a quantization error, it is negligibly small. Then Eq. (1) is approximately rewritten as
arg min xT Qx
xi;1 ¼ xi;0 , x2M;j ¼ x2M1;j ,
1 1 B 1 2 B B .. Ps ¼ B . B B @
¼ I ðPTs Ps Þ þ 2Ps Ps þ ðPTs Ps Þ I.
0 1 .. 1
.
..
. 2 1
1
C C C C. C C 1 A
(5)
1
Then the square sum of the discrete Laplacians of x can be represented by the quadratic form xT Qx,
Lðx; kÞ ¼ xT Qx þ kT ðy HxÞ,
(9)
where k ¼ ðl00 ; l01 ; . . . ; l0;M1 ; l10 ; . . . ; lM1;M1 ÞT is the M 2 -dimensional vector consisting of the Lagrange multipliers. Then the following equations must hold: q Lðx; kÞ ¼ 2Qx H T k ¼ 0, ð10Þ qx q Lðx; kÞ ¼ y Hx ¼ 0. ð11Þ qk Since Q1 ¼ PT P1 ¼ 0, the matrix Q is found to be singular. Therefore, it is difficult to straightly solve the system of Eqs. (10) and (11). We thus perform the similar transformations of Q and H by the DCT matrix, and we show that they become sparse matrices.
Fig. 2. Neumann boundary condition for the image x (rectangular region).
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2.3. Similar transformations of Q and H
with
We put the one-dimensional (1-D) DCT matrix of size M as W M . The ij element of W M is given by pffiffiffi 2 piðj þ 0:5Þ ðW M Þij ¼ pffiffiffiffiffiffi cðiÞ cos (12) M M with
8 < p1ffiffiffi ; cðiÞ ¼ 2 : 1
i ¼ 0;
(13)
otherwise:
We should note that W M is a unitary matrix. We ~ y~ , and k~ by transform the variables x, y, and k into x, the two-dimensional (2-D) DCT matrix as follows: x~ ¼ ðW 2M W 2M Þx, y~ ¼ ðW M W M Þy, k~ ¼ ðW M W M Þk.
2373
ð14Þ ð15Þ
q~ ij ¼ ðp~ i þ p~ j Þ2 .
(24)
e as We can rewrite H e ¼ ðW M W M ÞHðW 2M W 2M ÞT H ¼ ðW M W M ÞðH s H s ÞðW 2M W 2M ÞT ¼ ðW M H s W T2M Þ ðW M H s W T2M Þ. e s, es H ¼H
ð25Þ
where we put e s ¼ W M H sW T H 2M 0~ h0 B B h~1 B ¼B .. B B . @
ð16Þ
0
h~M1
0
h~Mþ1
C h~2M1 C C C, C C c A
Then Eqs. (10) and (11) can be rewritten as T
ex~ H e k~ ¼ 0, 2Q e x~ ¼ 0, y~ H
ð26Þ ð17Þ ð18Þ
e and H e are the similar where the matrices Q transformations of Q and H defined by e ¼ ðW 2M W 2M ÞQðW 2M W 2M ÞT , Q e ¼ ðW M W M ÞHðW 2M W 2M ÞT . H
1
ð19Þ
with
8 1 pi > > p ffiffi ffi cos ; > > 4M > 2 < h~i ¼ 0; > > 1 pi > > > p ffiffi ffi cos ; : 4M 2
0pipM 1; i ¼ M; 2M; M þ 1pip2M 1:
ð20Þ
(27)
e and H e are sparse matrices. We shall show that Q We note that Ps can be diagonalized by the DCT matrix as [19]
The derivation of Eq. (26) is shown in Appendix. e and H e are sparse We see from (23) and (26) that Q matrices.
es ¼ W 2M Ps W T ¼ diagfp~ 0 ; p~ 0 ; . . . ; p~ 2M1 g P 2M
(21)
with
pi p~ i ¼ 2 1 cos . 2M
(22)
Using (6) and (21), we can straightforwardly show e is a diagonal matrix as follows: that Q e ¼ ðW 2M W 2M ÞQðW 2M W 2M ÞT Q ¼ ðW 2M W 2M ÞðI ðPTs Ps Þ þ 2Ps Ps þ
ðPTs Ps Þ
IÞðW 2M W 2M Þ
2.4. Fast computation by the DCT Here we show that the solution of the system of Eqs. (17) and (18) can be written in a scalar form by using (23) and (26). Using (26), we can rewrite Eq. (18) in a scalar form as 8 ~ ~ > < h0 h0 x~ 00 ; y~ kl ¼ h~k h~l x~ kl þ h~k h~2Ml x~ k;2Ml þ h~2Mk h~l x~ 2Mk;l > : þh~2Mk h~2Ml x~ 2Mk;2Ml
k ¼ l ¼ 0; otherwise:
(28)
T
From the top equation in (28), we can straightforwardly determine x~ 00 from y~ 00 as
¼ I ðW 2M PTs Ps W T2M Þ þ 2ðW 2M Ps W T2M Þ ðW 2M Ps W T2M Þ þ ðW 2M PTs Ps W T2M Þ I eT P e e e eT e ¼ I ðP s s Þ þ 2Ps Ps þ ðPs Ps Þ I ¼ diagfq~ 00 ; q~ 01 ; ; q~ 0;2M1 ; q~ 10 ; ; q~ 2M1;2M1 g ð23Þ
x~ 00 ¼
y~ 00 . 2 h~
(29)
0
Therefore, we shall consider the case of i40 or j40 to determine x~ ij from y~ kl . Substituting (23) and (25)
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2374
into (17) gives ( 0; 2q~ ij x~ ij ¼ ~ ~ ~ hi hj lkðiÞkðjÞ
3. Consideration and reduction of additive noise i ¼ M or j ¼ M; otherwise
with the definition ( i; 0pipM; kðiÞ ¼ 2M i; M þ 1pip2M 1: Since q~ ij a0, Eq. (30) is rewritten as 8 0; i ¼ M or j ¼ M; > < ~ ~ ~ x~ ij ¼ hi hj lkðiÞkðjÞ otherwise: > : 2q~ ij
(30)
3.1. Robustness against additive noise
(31)
In the previous section, we have developed the HR image restoration method in the absence of the additive noise v. In practice, there are noises in the observations as in (1). We here analyze the influence of the additive noise v on the image distortion. We define the image distortion caused by the additive noise v as
(32)
Substituting (32) into the bottom equation of (28), we have 2y~ l~ kl ¼ kl Dkl
(33)
e ¼ x 1 x0 ,
(36)
where x0 and x1 are the restored HR images in the absence or presence of the additive noise v, respectively. Using an appropriate ð4M 2 M 2 Þ e we can rewrite Eq. (35) in matrix form as matrix B, ey. x~ ¼ B~
(37)
Then the DCT of e is computed by
with the definition 2 2
2 2
2
2
~ ~ 1 ~ ~ ~ ~ 1 Dkl ¼ h~k h~l q~ 1 kl þ hk h2Ml q k;2Ml þ h2Mk hl q 2Mk;l 2 2 þ h~2Mk h~2Ml q~ 1 2Mk;2Ml .
ð34Þ
Here we used the property that DkðiÞ;kðjÞ ¼ Dij for 0pi; jp2M 1. Substituting (33) into (30), and summarizing it together with the result for the case of i ¼ j ¼ 0, we have 8 y~ 00 > > ; i ¼ j ¼ 0; > > ~ ~ > > < h0 h0 i ¼ M or j ¼ M; (35) x~ ij ¼ 0; > ~ ~ > h h > i j > > otherwise: y~ > : q~ ij Dij kðiÞkðjÞ This equation shows that the DCT of the HR image x~ ij can be computed from the DCT of the LR image y~ kl by scalar computation. The algorithm of the proposed method is summarized in Fig. 3. We can compute x~ from y~ by using (35) in OðM 2 Þ. The DCT of y and the inverse DCT of x~ can be computed in OðM 2 log MÞ. Therefore, the proposed method can compute the HR image x from the LR image y in OðM 2 log MÞ processing time. Step 1. Compute the DCT of y by Equation (15) to obtain ~y Step 2. Compute x~ from y~ by Equation (35) Step 3. Compute the inverse DCT of x~ by Equation (14) to obtain x
Fig. 3. Algorithm of the proposed method.
e~ ¼ ðW 2M W 2M Þe ¼ ðW 2M W 2M Þx1 ðW 2M W 2M Þx0 e y þ v~ Þ B~ ey ¼ Bð~ ev, ¼ B~
ð38Þ
where v~ is the DCT of v defined by v~ ¼ ðW M W M Þv.
(39)
We can formulate the mean square error (MSE) of e as follows: 1 E½eT e 4M 2 1 E½~eT e~ ¼ 4M 2 1 evÞT ðB~ evÞ. E½ðB~ ð40Þ ¼ 4M 2 Assuming that the additive noise v is spatially uncorrelated, v~ is also spatially uncorrelated as
MSE ¼
E½~vv~ T ¼ E½vvT ¼ 2 I.
(41)
Therefore, the MSE can be written as 1 evÞT ðB~ evÞ E½ðB~ 4M 2 1 eT Þ evv~ T B E½TrðB~ ¼ 2 4M ! X 2M1 X 2 1 2M1 b~ 2 , ¼ 4M 2 i¼0 j¼0 ij
MSE ¼
ð42Þ
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
where Tr() stands for the trace of the matrix, and b~ij is defined by 8 1 > > ; i ¼ j ¼ 0; > > ~ ~ > > < h0 h0 i ¼ M or j ¼ M; (43) b~ij ¼ 0; > ~ ~ > h h > i j > > otherwise: > : q~ ij Dij From Eq. (42), we see that the MSE is in proportion to the variance of the additive noise 2 . Table 1 summarizes the proportional constant P P 2M1 ~2 b ð1=4M 2 Þ 2M1 for M ¼ 128, 256, 512, ij i¼0 j¼0 and 1024. It is interesting to see that the proportional constants are almost the same for various values of M. We have analytically derived the proportional constant of the proposed method. Here we evaluate the proportional constant of the cubic spline interpolation by computer simulations since it is difficult to compute it analytically. We generated LR Gaussian noise image, interpolated the noise image to obtain HR image by the cubic spline interpolation, and then computed the ratio of the variances of HR and LR images. Then the proportional constant was found to be 0.758. Since the constant of the proposed method is larger than 0.758, the proposed method is less robust against additive noise than the cubic spline. Therefore, as the noise variance becomes larger, the restoration performance of the proposed method gets worse. We should apply a noise reduction method to the proposed method when the additive noise is not so small. 3.2. Noise reduction in DCT domain Soon et al. has proposed a noise reduction method in the DCT domain [18]. Since the DCT of the LR image is computed in the proposed algorithm, the noise reduction method is immediately applicable to the proposed method by inserting it after step 1 of Fig. 3. Here we summarize the basic concept of the noise reduction method in [18]. Table 1 Proportional constant for various values of M M
128
256
512
1024
Constant
1.176444
1.176218
1.176105
1.176049
2375
We put the DCT of noisy LR image as y~ kl ¼ zkl þ v~kl ,
(44)
where zkl is the DCT of clean LR image and v~kl is the DCT of additive noise. With the assumption that DCT coefficients are statistically independent, the minimum mean square error (MMSE) estimated coefficient z^kl is expressed as z^kl ¼ E½zkl jy~ kl .
(45)
Here we define the lexicographically ordered vector consisting of z^kl as z^ . From Bayes’ theorem and the Gaussian distribution assumptions, we can rewrite it as z^kl ¼
Zz ðk; lÞ y~ , Zz ðk; lÞ þ Zv ðk; lÞ kl
(46)
where Zz ðk; lÞ ¼ E½jzkl j2 and Zv ðk; lÞ ¼ E½j~vkl j2 . When combining the noise reduction method with the proposed method, we compute z^ from y~ by using Eq. (46) after step 1, and then replace y~ by z^ in step 2. The increase of computation time is trivial since the noise reduction method requires only M 2 multiplications. 4. Simulation results In this section, we compared the restoration performances of the conventional and proposed methods. All simulations were done on an IBM PC/ AT compatible computer with an Intel Pentium 4 2.4 GHz and 512 Mbyte DRAM’s. We used eight images of size ð256 256Þ with 8 bit grayscale: ‘‘barbara’’, ‘‘boat’’, ‘‘bridge’’, ‘‘camera’’, ‘‘girl’’, ‘‘lena’’, ‘‘mandrill’’, and ‘‘peppers’’. First, we generated the ð128 128Þ LR image y from the original ð256 256Þ HR image x according to Eq. (7). Then we restored the ð256 256Þ HR image from y by the linear interpolation, the cubic, and the 5th order spline interpolations [8], and the proposed methods. Fig. 4 shows the original ‘‘lena’’ and the restored images by the five methods. The restored image by the proposed method seems to be more ‘‘high-passed’’ than the other three images. We next quantitatively measured the restoration performance of each method by the peak signal-tonoise ratio (PSNR) defined by ! 2552 PSNR ¼ 10 log10 P P2M1 2 , ð1=4M 2 Þ 2M1 i¼0 j¼0 eij (47)
ARTICLE IN PRESS 2376
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
Fig. 4. Original and restored images: (a) original ‘‘lena’’, (b) linear interpolation, (c) cubic spline interpolation, (d) 5th order spline interpolation, and (e) proposed method.
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
proposed method is much improved by including the noise reduction method, and that it is better than that of the cubic spline interpolation method. The computation time of the proposed method with
3 Difference in PSNR [dB]
barbara
2.5
boat
2
bridge camera
1.5
girl
1
lena mandrill
0.5
peppers
0 -0.5 -1 -1.5 0
20
40
60
80
100
2 Fig. 5. Differences of the PSNRs between the proposed and cubic spline interpolation methods versus the noise variance 2 .
3 barbara
Difference in PSNR [dB]
where eij is the difference of the pixel value between the original and restored images. Table 2 summarizes the PSNRs of eight images restored by the four methods. This shows that the proposed method is better than the linear interpolation method by 1.4 dB, and is better than both the cubic and the 5th order spline interpolation methods by 0.3 dB. The computation times of the proposed, linear, cubic spline, and 5th order spline interpolation methods were 0.011 [s], 0.0014 [s], 0.0052 [s], and 0.010 [s], respectively. The computation time of the 5th order spline is the same as that of our method while the PSNR is almost the same as that of the cubic spline. The proposed method achieves better restoration performance at the expense of an increase of computation time. We also evaluated the robustness of the proposed method against the additive noise. Here we generated the LR image y by Eq. (1), where we used a white Gaussian noise of the variance 2 as the additive noise v. Then we restored the HR image x from the noisy LR image y, and measured the PSNR between the original and restored images. We increased the noise variance 2 from 0 to 100 by 0.5. Fig. 5 shows the differences between the PSNRs of the proposed and cubic spline interpolation methods for the eight images. The PSNRs of the proposed method are better than those of the cubic spline interpolation method in all images as far as 2 is small. However, the cubic spline interpolation method is better when 2 is large, as we have theoretically demonstrated in Section 3.1. We combined the proposed method with the noise reduction method described in Section 3.2, and performed the same simulation to see the effectiveness of the noise reduction method. In the simulations, we assume that Zz ðk; lÞ and Zv ðk; lÞ are known, although they should be properly estimated in practice. Fig. 6 shows the differences between the PSNRs of the proposed method with noise reduction and the cubic spline interpolation method for the eight images. We see that the PSNR of the
2377
2.5
boat bridge
2
camera
1.5
lena
1
peppers
girl
mandrill
0.5 0 0
20
40
60
80
100
2
Fig. 6. Differences of the PSNRs between the proposed method with noise reduction and the cubic spline interpolation method versus the noise variance 2 .
Table 2 PSNRs for the three interpolation methods and the proposed method
Linear Cubic spline Spline (5th) Proposed
Barbara
Boat
Bridge
Camera
Girl
Lena
Mandrill
Peppers
Average
27.03 27.86 27.90 28.06
29.30 30.48 30.44 30.69
22.97 23.51 23.50 23.65
25.66 26.64 26.65 26.89
31.81 33.59 33.67 34.03
28.45 29.65 29.71 29.95
23.84 24.73 24.73 25.00
26.44 28.14 28.27 28.59
26.94 28.08 28.11 28.36
ARTICLE IN PRESS 2378
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
noise reduction was almost the same as the one without the noise reduction. 5. Conclusion We derived the HR image restoration method from a downsampled LR image by using the discrete cosine transform (DCT). The restoration performance of the proposed method is superior to that of the cubic spline interpolation at the expense of an increase of computation time. We theoretically evaluated the influence of the observation noise on the restoration performance, and showed
that the proposed method is still better than the cubic spline interpolation as far as the amount of the additive noise is small. We also showed that the noise reduction method is straightforwardly applicable to the proposed method, and it can efficiently reduce the noise effect. In this paper, we only considered the downsampling matrix in the form of Eq. (3). There may be a further work to generalize the downsampling process in the proposed method. For this purpose, we should derive the condition under which a downsampling matrix can be transformed into a sparse matrix by the DCT.
Appendix A. Derivation of Eq. (26) e s are The ij elements of H e s Þij ¼ ðW M H s W T Þij ðH 2M pffiffiffi M1 2cðiÞcðjÞ X piðl þ 0:5Þ ¼ cos M M l¼0 1 pjð2l þ 0:5Þ cos 2 2M 1 pjð2l þ 1 þ 0:5Þ þ cos . 2 2M Using the property of the cosine function, we can establish the following two relations: piðl þ 0:5Þ pjð2l þ 0:5Þ cos cos M 2M 1 pðið2l þ 1Þ þ jð2l þ 1 0:5ÞÞ pðið2l þ 1Þ jð2l þ 1 0:5ÞÞ cos ¼ þ cos 2 2M 2M 1 pðð2l þ 1Þði þ jÞ 0:5jÞ pðð2l þ 1Þði jÞ þ 0:5jÞ cos ¼ þ cos , 2 2M 2M piðl þ 0:5Þ pjð2l þ 1 þ 0:5Þ cos cos M 2M 1 pðið2l þ 1Þ þ jð2l þ 1 þ 0:5ÞÞ pðið2l þ 1Þ jð2l þ 1 þ 0:5ÞÞ cos ¼ þ cos 2 2M 2M 1 pðð2l þ 1Þði þ jÞ þ 0:5jÞ pðð2l þ 1Þði jÞ 0:5jÞ cos ¼ þ cos . 2 2M 2M Substituting these equations into Eq. (A.1) yields pffiffiffi X 1 pðð2l þ 1Þði þ jÞ 0:5jÞ pðð2l þ 1Þði jÞ þ 0:5jÞ 2cðiÞcðjÞ M1 e cos ðH s Þij ¼ þ cos 4 2M 2M M l¼0 pðð2l þ 1Þði þ jÞ þ 0:5jÞ pðð2l þ 1Þði jÞ 0:5jÞ þ cos þ cos . 2M 2M
ðA:1Þ
ðA:2Þ
ðA:3Þ
ðA:4Þ
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
Here we also have the following two equations: pðði þ jÞð2l þ 1Þ 0:5jÞ pðði þ jÞð2l þ 1Þ þ 0:5jÞ cos þ cos 2M 2M pði þ jÞð2l þ 1Þ pj ¼ 2 cos cos , 2M 4M pðði jÞð2l þ 1Þ þ 0:5jÞ pðði jÞð2l þ 1Þ 0:5jÞ cos þ cos 2M 2M pði jÞð2l þ 1Þ pj ¼ 2 cos cos . 2M 4M Substituting these equations into (A.4) gives M1 X pði þ jÞð2l þ 1Þ e s Þij ¼ cðiÞcðjÞ pffiffiffi ðH cos 2M 2M l¼0 pði jÞð2l þ 1Þ pj þ cos cos . 2M 4M
2379
ðA:5Þ
ðA:6Þ
ðA:7Þ
Here we note that
8 k ¼ 4mM; > < M; pkð2l þ 1Þ cos ¼ M; k ¼ ð4m þ 2ÞM; > 2M : l¼0 0 otherwise: e s are Therefore, the ij elements of H 8 1 pj > > pffiffiffi cos ; i j ¼ 0; > > > 4M < 2 e s Þij ¼ ðH 1 pj > ; i þ j ¼ 2M; pffiffiffi cos > > 4M 2 > > : 0 otherwise: M1 X
(A.8)
(A.9)
This is expressed in the matrix form as (26).
References [1] S.C. Park, M.K. Park, M.G. Kang, Super-resolution image reconstruction—A technical overview, Signal Process. Mag. 20 (3) (2003) 21–36. [2] S. Rhee, M.G. Kang, Discrete cosine transform based regularized high-resolution image reconstruction algorithm, Opt. Eng. 38 (8) (1999) 1348–1356. [3] H. Sakane, K. Nishikawa, H. Kiya, A super-resolution method based on the discrete cosine transforms, Eusipco’96, 1996. [4] R.C. Hardie, K.J. Barnard, J.G. Bognar, E.E. Armstrong, E.A. Watson, High-resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system, Opt. Eng. 37 (1) (January 1998) 247–260. [5] J.H. Shin, J.H. Jung, J.K. Paik, Regularized iterative image interpolation and its application to spatially scalable coding, IEEE Trans. Consum. Electron. 44 (3) (August 1998) 1042–1047.
[6] M. Unser, A. Aldroubi, M. Eden, Enlargement or reduction of digital images with minimum loss of information, IEEE Trans. Image Process. 4 (3) (March 1995) 247–258. [7] M. Unser, A. Aldroubi, M. Eden, Polynomial spline signal approximations: filter design and asymptotic equivalence with Shannon’s sampling theorem, IEEE Trans. Inform. Theory 38 (1) (1992) 95–103. [8] P. The´venaz, T. Blu, M. Unser, Image interpolation and resampling, in: I.N. Bankman (Ed.), Handbook of Medical Imaging, Processing and Analysis, Academic Press, San Diego CA, USA, 2000, pp. 393–420. [9] S.A. Martucci, Image resizing in the discrete cosine transform domain, in: ICIP’95, vol. 2, 1995, pp. 224–227. [10] E. Shinbori, M. Takagi, High quality magnification applying the Gerchberg–Papoulis iterative algorithm with DCT, Trans. IEICE J76-D (9) (September 1993) 1932–1940. [11] S.E. El-Khamy, M.M. Hadhoud, M.I. Dessouky, B.M. Salam, F.E. Abd El-Samie, Efficient implementation of image interpolation as an inverse problem, Digit. Signal Process. 15 (2005) 137–152.
ARTICLE IN PRESS 2380
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
[12] M. Ng, N. Bose, Mathematical analysis of super-resolution methodology, IEEE Signal Process. Mag. 20 (3) (May 2003) 62–74. [13] L. Chen, K. Yap, Regularized interpolation using Kronecker product for still images, in: ICIP 2005, vol. 2, 2005 pp. II1014–1017. [14] A. Rosenfeld, A.C. Kak, Digital Picture Processing, Academic Press, New York, 1976. [15] K. Horiike, Y. Yoshida, K. Fujita, Image restoration based on a Wiener filtering technique using discrete cosine transform, Syst. Comput. Japan (SCJAEP) 27 (13) (1996) 23–35. [16] S.E. El-Khamy, M.M. Hadhoud, M.I. Dessouky, B.M. Salam, F.E. Abd El-Samie, Sectioned implementation
of regularized image interpolation, in: Proceedings of the 46th IEEE International Midwest Symposium on Circuits and Systems (MWSCAS ’03), vol. 2, 2003, pp. 656–659. [17] B.R. Hunt, The application of constrained least squares estimation to image restoration by digital computer, IEEE Trans. Comput. C-22 (9) (1973) 805–812. [18] I.Y. Soon, S.N. Koh, C.K. Yeo, Noisy speech enhancement using discrete cosine transform, Speech Comm. 24 (1998) 249–257. [19] F. O’Sullivan, Discretized Laplacian smoothing by Fourier methods, J. Amer. Statist. Assoc. 86 (September 1991) 634–642.
Signal Processing 87 (2007) 2370–2380 www.elsevier.com/locate/sigpro
Image restoration from a downsampled image by using the DCT Yoshinori Abe, Youji Iiguni Department of Systems Innovation, Graduate School of Engineering Science, Osaka University, Toyonaka 560-8531, Japan Received 5 September 2006; received in revised form 18 January 2007; accepted 20 March 2007 Available online 30 March 2007
Abstract The high-resolution (HR) image restoration from a downsampled low-resolution (LR) image by using the discrete cosine transform (DCT) is proposed. The downsampling process is modeled in matrix form and the DCT is applied to the downsampling matrix to obtain the sparse matrix representation. The restored HR image is then efficiently computed from the LR image by using the sparsity. The distortion caused by the observation noise is theoretically analyzed, and the power of distortion is shown to be closely equal to the variance of the observation noise. Computer simulations show that the proposed method is superior to the cubic spline interpolation in the HR image restoration performance as far as the amount of the additive noise is small. r 2007 Elsevier B.V. All rights reserved. Keywords: High resolution image restoration; Discrete cosine transform; Lagrange multiplier method; Similar transformation
1. Introduction High-resolution (HR) image restorations, which are acquiring a HR image from downsampled lowresolution (LR) images, attract much attention [1–3]. In the HR image restoration methods, the downsampling process, which expresses the relation between the pixel values of the LR and HR images, is mathematically modeled. When using a simple charge coupled device (CCD) image capturing model, we can express the pixel value of the LR image as the sum of corresponding pixel values of the HR image [4,5]. Then the HR image restoration results in an inverse problem of the downsampling process. Unfortunately, the inverse problem beCorresponding author. Tel./fax: +81 6 6850 6377.
E-mail addresses: [email protected] (Y. Abe), [email protected] (Y. Iiguni).
comes ill defined when only one LR image can be used to restore the HR image. Nevertheless, if the LR image is under the Nyquist condition, the sinc interpolation can reconstruct the HR image perfectly from the LR image. However, it is not suitable for practical use because of the large computational cost. Other image interpolation methods that approximate the sinc interpolation with less complexity, such as the cubic spline interpolation [6–8], and image resizing methods in the frequency domain [9,10], have been proposed. However, the actual LR image does not necessarily satisfy the Nyquist condition, and therefore the restoration performance is low. To achieve a good restoration performance, HR image restoration methods which make the illdefined restoration problem well defined have been derived [4,11–13] from the analogy with the blurred image restoration method [14,15]. In the blurred
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2007.03.010
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
image restoration, a smoothness constraint such as the square sum of the discrete Laplacian of the restored image is imposed on the restored image to uniquely restore the blurred image [14]. Similarly, in the HR image restoration, the smoothness constraint is imposed on the restored HR image to determine it uniquely from the downsampled LR image [11]. However, the direct solution of the HR image restoration problem requires a large computation time, because a large-scale linear equation needs to be solved. A block segmentation method [16] or the Kronecker product method [13] has been proposed for the fast computation, but they decrease the restoration performance. In this paper, we derive a fast computation method for the HR image restoration by using the discrete cosine transform (DCT). In the blurred image restoration, Hunt has developed a fast computation method with the fast Fourier transform (FFT) [17]. The main idea of this method is the diagonalization of the blurring matrix by the discrete Fourier transform (DFT) matrix. Unfortunately, the fast computation method cannot be straightly applied to the current problem because the downsampling matrix is not square. We therefore model the downsampling process in the absence of the observation noise, perform the DCT to obtain a sparse representation, and then derive the analytical solution to the restoration problem. Since we can express the analytical solution in a scalar form, we can quickly generate the HR image from the LR image by using the DCT. We show through computer simulations that the restoration capability of the proposed method is superior to that of the cubic spline interpolation method. The proposed method considers a small amount of observation noises, such as quantization errors, in the downsampling process. We thus neglect the observation noise, and derive the analytical solution in the absence of noise. However, there are observation noises in practice. We theoretically analyze the distortion caused by the observation noise, and show that the power of distortion is closely equal to the variance of the observation noise. We also show through computer simulations that the proposed method is still superior to the cubic spline interpolation in the restoration performance if the amount of the observation noise is small. Finally, we combine the noise reduction method [18] with the proposed method to decrease the noise effect.
2371
2. Mathematical formulation of the image restoration 2.1. Observation model Let the ð2M 2MÞ HR image be fxij g and the ðM MÞ LR image be fykl g. The HR and LR images are obtained from the same scene captured by the different resolution cameras at the same position. The array structures of the photodetectors of the cameras are described in Fig. 1. The density of the photodetectors of the HR image is 4 times higher than that of the LR image. Then we can express the relation between HR and LR images as y ¼ Hx þ v,
(1)
where x, y, and v are lexicographically ordered vectors consisting of the HR image fxij g, the LR image fykl g, and the additive noise with mean zero and variance 2 , respectively. Here we assume that the additive noise is very small. The ðM 2 4M 2 Þ matrix H represents the filtering and down-sampling process and has the form H ¼ H s H s,
(2)
where ‘‘’’ denotes the Kronecker product and H s is the ðM 2MÞ matrix defined by 1 01 1 0 2 2 C B 1 1 C B C B 2 2 C B ð3Þ Hs ¼ B C: .. C B . C B A @ 1 1 0 2 2
Fig. 1. Array structure of photodetectors: small squares represent the photodetectors for HR image and hatched squares represent the ones for LR image.
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2372
2.2. HR image restoration problem
where the ð4M 2 4M 2 Þ matrix Q is defined by
We shall consider the problem of restoring the HR image x from the LR image y. Since the problem is ill defined, we impose a smoothness constraint on x to uniquely determine x. The constraint we imposed is the square sum of the discrete Laplacian of x. The discrete Laplacian of each xij is defined by xi1;j þ xi;j1 4xij þ xi;jþ1 þ xiþ1;j . We here suppose that xij s at image boundaries i ¼ 0, i ¼ 2M 1, j ¼ 0, or j ¼ 2M 1 are under the Neumann boundary condition, that is,
Q ¼ PT P
x1;j ¼ x0;j ,
y ¼ Hx.
(7)
Under the condition (7) we shall estimate the HR image x from the LR image y by minimizing the square sum of the discrete Laplacian xT Qx. The estimation problem is formulated as x
subject to y ¼ Hx.
(8)
We solve the problem (8) by using the Lagrange multiplier method. Let the Lagrangian L be
xi;2M ¼ xi;2M1 . The boundary condition is illustrated in Fig. 2. The lexicographically ordered vector consisting of the discrete Laplacians of fxij g is then expressed as [19] Px ¼ ðI Ps þ Ps IÞx,
(4)
where I is the ð2M 2MÞ identity matrix, and Ps is the ð2M 2MÞ matrix defined by 0
0
ð6Þ
When the additive noise v is a quantization error, it is negligibly small. Then Eq. (1) is approximately rewritten as
arg min xT Qx
xi;1 ¼ xi;0 , x2M;j ¼ x2M1;j ,
1 1 B 1 2 B B .. Ps ¼ B . B B @
¼ I ðPTs Ps Þ þ 2Ps Ps þ ðPTs Ps Þ I.
0 1 .. 1
.
..
. 2 1
1
C C C C. C C 1 A
(5)
1
Then the square sum of the discrete Laplacians of x can be represented by the quadratic form xT Qx,
Lðx; kÞ ¼ xT Qx þ kT ðy HxÞ,
(9)
where k ¼ ðl00 ; l01 ; . . . ; l0;M1 ; l10 ; . . . ; lM1;M1 ÞT is the M 2 -dimensional vector consisting of the Lagrange multipliers. Then the following equations must hold: q Lðx; kÞ ¼ 2Qx H T k ¼ 0, ð10Þ qx q Lðx; kÞ ¼ y Hx ¼ 0. ð11Þ qk Since Q1 ¼ PT P1 ¼ 0, the matrix Q is found to be singular. Therefore, it is difficult to straightly solve the system of Eqs. (10) and (11). We thus perform the similar transformations of Q and H by the DCT matrix, and we show that they become sparse matrices.
Fig. 2. Neumann boundary condition for the image x (rectangular region).
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2.3. Similar transformations of Q and H
with
We put the one-dimensional (1-D) DCT matrix of size M as W M . The ij element of W M is given by pffiffiffi 2 piðj þ 0:5Þ ðW M Þij ¼ pffiffiffiffiffiffi cðiÞ cos (12) M M with
8 < p1ffiffiffi ; cðiÞ ¼ 2 : 1
i ¼ 0;
(13)
otherwise:
We should note that W M is a unitary matrix. We ~ y~ , and k~ by transform the variables x, y, and k into x, the two-dimensional (2-D) DCT matrix as follows: x~ ¼ ðW 2M W 2M Þx, y~ ¼ ðW M W M Þy, k~ ¼ ðW M W M Þk.
2373
ð14Þ ð15Þ
q~ ij ¼ ðp~ i þ p~ j Þ2 .
(24)
e as We can rewrite H e ¼ ðW M W M ÞHðW 2M W 2M ÞT H ¼ ðW M W M ÞðH s H s ÞðW 2M W 2M ÞT ¼ ðW M H s W T2M Þ ðW M H s W T2M Þ. e s, es H ¼H
ð25Þ
where we put e s ¼ W M H sW T H 2M 0~ h0 B B h~1 B ¼B .. B B . @
ð16Þ
0
h~M1
0
h~Mþ1
C h~2M1 C C C, C C c A
Then Eqs. (10) and (11) can be rewritten as T
ex~ H e k~ ¼ 0, 2Q e x~ ¼ 0, y~ H
ð26Þ ð17Þ ð18Þ
e and H e are the similar where the matrices Q transformations of Q and H defined by e ¼ ðW 2M W 2M ÞQðW 2M W 2M ÞT , Q e ¼ ðW M W M ÞHðW 2M W 2M ÞT . H
1
ð19Þ
with
8 1 pi > > p ffiffi ffi cos ; > > 4M > 2 < h~i ¼ 0; > > 1 pi > > > p ffiffi ffi cos ; : 4M 2
0pipM 1; i ¼ M; 2M; M þ 1pip2M 1:
ð20Þ
(27)
e and H e are sparse matrices. We shall show that Q We note that Ps can be diagonalized by the DCT matrix as [19]
The derivation of Eq. (26) is shown in Appendix. e and H e are sparse We see from (23) and (26) that Q matrices.
es ¼ W 2M Ps W T ¼ diagfp~ 0 ; p~ 0 ; . . . ; p~ 2M1 g P 2M
(21)
with
pi p~ i ¼ 2 1 cos . 2M
(22)
Using (6) and (21), we can straightforwardly show e is a diagonal matrix as follows: that Q e ¼ ðW 2M W 2M ÞQðW 2M W 2M ÞT Q ¼ ðW 2M W 2M ÞðI ðPTs Ps Þ þ 2Ps Ps þ
ðPTs Ps Þ
IÞðW 2M W 2M Þ
2.4. Fast computation by the DCT Here we show that the solution of the system of Eqs. (17) and (18) can be written in a scalar form by using (23) and (26). Using (26), we can rewrite Eq. (18) in a scalar form as 8 ~ ~ > < h0 h0 x~ 00 ; y~ kl ¼ h~k h~l x~ kl þ h~k h~2Ml x~ k;2Ml þ h~2Mk h~l x~ 2Mk;l > : þh~2Mk h~2Ml x~ 2Mk;2Ml
k ¼ l ¼ 0; otherwise:
(28)
T
From the top equation in (28), we can straightforwardly determine x~ 00 from y~ 00 as
¼ I ðW 2M PTs Ps W T2M Þ þ 2ðW 2M Ps W T2M Þ ðW 2M Ps W T2M Þ þ ðW 2M PTs Ps W T2M Þ I eT P e e e eT e ¼ I ðP s s Þ þ 2Ps Ps þ ðPs Ps Þ I ¼ diagfq~ 00 ; q~ 01 ; ; q~ 0;2M1 ; q~ 10 ; ; q~ 2M1;2M1 g ð23Þ
x~ 00 ¼
y~ 00 . 2 h~
(29)
0
Therefore, we shall consider the case of i40 or j40 to determine x~ ij from y~ kl . Substituting (23) and (25)
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
2374
into (17) gives ( 0; 2q~ ij x~ ij ¼ ~ ~ ~ hi hj lkðiÞkðjÞ
3. Consideration and reduction of additive noise i ¼ M or j ¼ M; otherwise
with the definition ( i; 0pipM; kðiÞ ¼ 2M i; M þ 1pip2M 1: Since q~ ij a0, Eq. (30) is rewritten as 8 0; i ¼ M or j ¼ M; > < ~ ~ ~ x~ ij ¼ hi hj lkðiÞkðjÞ otherwise: > : 2q~ ij
(30)
3.1. Robustness against additive noise
(31)
In the previous section, we have developed the HR image restoration method in the absence of the additive noise v. In practice, there are noises in the observations as in (1). We here analyze the influence of the additive noise v on the image distortion. We define the image distortion caused by the additive noise v as
(32)
Substituting (32) into the bottom equation of (28), we have 2y~ l~ kl ¼ kl Dkl
(33)
e ¼ x 1 x0 ,
(36)
where x0 and x1 are the restored HR images in the absence or presence of the additive noise v, respectively. Using an appropriate ð4M 2 M 2 Þ e we can rewrite Eq. (35) in matrix form as matrix B, ey. x~ ¼ B~
(37)
Then the DCT of e is computed by
with the definition 2 2
2 2
2
2
~ ~ 1 ~ ~ ~ ~ 1 Dkl ¼ h~k h~l q~ 1 kl þ hk h2Ml q k;2Ml þ h2Mk hl q 2Mk;l 2 2 þ h~2Mk h~2Ml q~ 1 2Mk;2Ml .
ð34Þ
Here we used the property that DkðiÞ;kðjÞ ¼ Dij for 0pi; jp2M 1. Substituting (33) into (30), and summarizing it together with the result for the case of i ¼ j ¼ 0, we have 8 y~ 00 > > ; i ¼ j ¼ 0; > > ~ ~ > > < h0 h0 i ¼ M or j ¼ M; (35) x~ ij ¼ 0; > ~ ~ > h h > i j > > otherwise: y~ > : q~ ij Dij kðiÞkðjÞ This equation shows that the DCT of the HR image x~ ij can be computed from the DCT of the LR image y~ kl by scalar computation. The algorithm of the proposed method is summarized in Fig. 3. We can compute x~ from y~ by using (35) in OðM 2 Þ. The DCT of y and the inverse DCT of x~ can be computed in OðM 2 log MÞ. Therefore, the proposed method can compute the HR image x from the LR image y in OðM 2 log MÞ processing time. Step 1. Compute the DCT of y by Equation (15) to obtain ~y Step 2. Compute x~ from y~ by Equation (35) Step 3. Compute the inverse DCT of x~ by Equation (14) to obtain x
Fig. 3. Algorithm of the proposed method.
e~ ¼ ðW 2M W 2M Þe ¼ ðW 2M W 2M Þx1 ðW 2M W 2M Þx0 e y þ v~ Þ B~ ey ¼ Bð~ ev, ¼ B~
ð38Þ
where v~ is the DCT of v defined by v~ ¼ ðW M W M Þv.
(39)
We can formulate the mean square error (MSE) of e as follows: 1 E½eT e 4M 2 1 E½~eT e~ ¼ 4M 2 1 evÞT ðB~ evÞ. E½ðB~ ð40Þ ¼ 4M 2 Assuming that the additive noise v is spatially uncorrelated, v~ is also spatially uncorrelated as
MSE ¼
E½~vv~ T ¼ E½vvT ¼ 2 I.
(41)
Therefore, the MSE can be written as 1 evÞT ðB~ evÞ E½ðB~ 4M 2 1 eT Þ evv~ T B E½TrðB~ ¼ 2 4M ! X 2M1 X 2 1 2M1 b~ 2 , ¼ 4M 2 i¼0 j¼0 ij
MSE ¼
ð42Þ
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
where Tr() stands for the trace of the matrix, and b~ij is defined by 8 1 > > ; i ¼ j ¼ 0; > > ~ ~ > > < h0 h0 i ¼ M or j ¼ M; (43) b~ij ¼ 0; > ~ ~ > h h > i j > > otherwise: > : q~ ij Dij From Eq. (42), we see that the MSE is in proportion to the variance of the additive noise 2 . Table 1 summarizes the proportional constant P P 2M1 ~2 b ð1=4M 2 Þ 2M1 for M ¼ 128, 256, 512, ij i¼0 j¼0 and 1024. It is interesting to see that the proportional constants are almost the same for various values of M. We have analytically derived the proportional constant of the proposed method. Here we evaluate the proportional constant of the cubic spline interpolation by computer simulations since it is difficult to compute it analytically. We generated LR Gaussian noise image, interpolated the noise image to obtain HR image by the cubic spline interpolation, and then computed the ratio of the variances of HR and LR images. Then the proportional constant was found to be 0.758. Since the constant of the proposed method is larger than 0.758, the proposed method is less robust against additive noise than the cubic spline. Therefore, as the noise variance becomes larger, the restoration performance of the proposed method gets worse. We should apply a noise reduction method to the proposed method when the additive noise is not so small. 3.2. Noise reduction in DCT domain Soon et al. has proposed a noise reduction method in the DCT domain [18]. Since the DCT of the LR image is computed in the proposed algorithm, the noise reduction method is immediately applicable to the proposed method by inserting it after step 1 of Fig. 3. Here we summarize the basic concept of the noise reduction method in [18]. Table 1 Proportional constant for various values of M M
128
256
512
1024
Constant
1.176444
1.176218
1.176105
1.176049
2375
We put the DCT of noisy LR image as y~ kl ¼ zkl þ v~kl ,
(44)
where zkl is the DCT of clean LR image and v~kl is the DCT of additive noise. With the assumption that DCT coefficients are statistically independent, the minimum mean square error (MMSE) estimated coefficient z^kl is expressed as z^kl ¼ E½zkl jy~ kl .
(45)
Here we define the lexicographically ordered vector consisting of z^kl as z^ . From Bayes’ theorem and the Gaussian distribution assumptions, we can rewrite it as z^kl ¼
Zz ðk; lÞ y~ , Zz ðk; lÞ þ Zv ðk; lÞ kl
(46)
where Zz ðk; lÞ ¼ E½jzkl j2 and Zv ðk; lÞ ¼ E½j~vkl j2 . When combining the noise reduction method with the proposed method, we compute z^ from y~ by using Eq. (46) after step 1, and then replace y~ by z^ in step 2. The increase of computation time is trivial since the noise reduction method requires only M 2 multiplications. 4. Simulation results In this section, we compared the restoration performances of the conventional and proposed methods. All simulations were done on an IBM PC/ AT compatible computer with an Intel Pentium 4 2.4 GHz and 512 Mbyte DRAM’s. We used eight images of size ð256 256Þ with 8 bit grayscale: ‘‘barbara’’, ‘‘boat’’, ‘‘bridge’’, ‘‘camera’’, ‘‘girl’’, ‘‘lena’’, ‘‘mandrill’’, and ‘‘peppers’’. First, we generated the ð128 128Þ LR image y from the original ð256 256Þ HR image x according to Eq. (7). Then we restored the ð256 256Þ HR image from y by the linear interpolation, the cubic, and the 5th order spline interpolations [8], and the proposed methods. Fig. 4 shows the original ‘‘lena’’ and the restored images by the five methods. The restored image by the proposed method seems to be more ‘‘high-passed’’ than the other three images. We next quantitatively measured the restoration performance of each method by the peak signal-tonoise ratio (PSNR) defined by ! 2552 PSNR ¼ 10 log10 P P2M1 2 , ð1=4M 2 Þ 2M1 i¼0 j¼0 eij (47)
ARTICLE IN PRESS 2376
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
Fig. 4. Original and restored images: (a) original ‘‘lena’’, (b) linear interpolation, (c) cubic spline interpolation, (d) 5th order spline interpolation, and (e) proposed method.
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
proposed method is much improved by including the noise reduction method, and that it is better than that of the cubic spline interpolation method. The computation time of the proposed method with
3 Difference in PSNR [dB]
barbara
2.5
boat
2
bridge camera
1.5
girl
1
lena mandrill
0.5
peppers
0 -0.5 -1 -1.5 0
20
40
60
80
100
2 Fig. 5. Differences of the PSNRs between the proposed and cubic spline interpolation methods versus the noise variance 2 .
3 barbara
Difference in PSNR [dB]
where eij is the difference of the pixel value between the original and restored images. Table 2 summarizes the PSNRs of eight images restored by the four methods. This shows that the proposed method is better than the linear interpolation method by 1.4 dB, and is better than both the cubic and the 5th order spline interpolation methods by 0.3 dB. The computation times of the proposed, linear, cubic spline, and 5th order spline interpolation methods were 0.011 [s], 0.0014 [s], 0.0052 [s], and 0.010 [s], respectively. The computation time of the 5th order spline is the same as that of our method while the PSNR is almost the same as that of the cubic spline. The proposed method achieves better restoration performance at the expense of an increase of computation time. We also evaluated the robustness of the proposed method against the additive noise. Here we generated the LR image y by Eq. (1), where we used a white Gaussian noise of the variance 2 as the additive noise v. Then we restored the HR image x from the noisy LR image y, and measured the PSNR between the original and restored images. We increased the noise variance 2 from 0 to 100 by 0.5. Fig. 5 shows the differences between the PSNRs of the proposed and cubic spline interpolation methods for the eight images. The PSNRs of the proposed method are better than those of the cubic spline interpolation method in all images as far as 2 is small. However, the cubic spline interpolation method is better when 2 is large, as we have theoretically demonstrated in Section 3.1. We combined the proposed method with the noise reduction method described in Section 3.2, and performed the same simulation to see the effectiveness of the noise reduction method. In the simulations, we assume that Zz ðk; lÞ and Zv ðk; lÞ are known, although they should be properly estimated in practice. Fig. 6 shows the differences between the PSNRs of the proposed method with noise reduction and the cubic spline interpolation method for the eight images. We see that the PSNR of the
2377
2.5
boat bridge
2
camera
1.5
lena
1
peppers
girl
mandrill
0.5 0 0
20
40
60
80
100
2
Fig. 6. Differences of the PSNRs between the proposed method with noise reduction and the cubic spline interpolation method versus the noise variance 2 .
Table 2 PSNRs for the three interpolation methods and the proposed method
Linear Cubic spline Spline (5th) Proposed
Barbara
Boat
Bridge
Camera
Girl
Lena
Mandrill
Peppers
Average
27.03 27.86 27.90 28.06
29.30 30.48 30.44 30.69
22.97 23.51 23.50 23.65
25.66 26.64 26.65 26.89
31.81 33.59 33.67 34.03
28.45 29.65 29.71 29.95
23.84 24.73 24.73 25.00
26.44 28.14 28.27 28.59
26.94 28.08 28.11 28.36
ARTICLE IN PRESS 2378
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
noise reduction was almost the same as the one without the noise reduction. 5. Conclusion We derived the HR image restoration method from a downsampled LR image by using the discrete cosine transform (DCT). The restoration performance of the proposed method is superior to that of the cubic spline interpolation at the expense of an increase of computation time. We theoretically evaluated the influence of the observation noise on the restoration performance, and showed
that the proposed method is still better than the cubic spline interpolation as far as the amount of the additive noise is small. We also showed that the noise reduction method is straightforwardly applicable to the proposed method, and it can efficiently reduce the noise effect. In this paper, we only considered the downsampling matrix in the form of Eq. (3). There may be a further work to generalize the downsampling process in the proposed method. For this purpose, we should derive the condition under which a downsampling matrix can be transformed into a sparse matrix by the DCT.
Appendix A. Derivation of Eq. (26) e s are The ij elements of H e s Þij ¼ ðW M H s W T Þij ðH 2M pffiffiffi M1 2cðiÞcðjÞ X piðl þ 0:5Þ ¼ cos M M l¼0 1 pjð2l þ 0:5Þ cos 2 2M 1 pjð2l þ 1 þ 0:5Þ þ cos . 2 2M Using the property of the cosine function, we can establish the following two relations: piðl þ 0:5Þ pjð2l þ 0:5Þ cos cos M 2M 1 pðið2l þ 1Þ þ jð2l þ 1 0:5ÞÞ pðið2l þ 1Þ jð2l þ 1 0:5ÞÞ cos ¼ þ cos 2 2M 2M 1 pðð2l þ 1Þði þ jÞ 0:5jÞ pðð2l þ 1Þði jÞ þ 0:5jÞ cos ¼ þ cos , 2 2M 2M piðl þ 0:5Þ pjð2l þ 1 þ 0:5Þ cos cos M 2M 1 pðið2l þ 1Þ þ jð2l þ 1 þ 0:5ÞÞ pðið2l þ 1Þ jð2l þ 1 þ 0:5ÞÞ cos ¼ þ cos 2 2M 2M 1 pðð2l þ 1Þði þ jÞ þ 0:5jÞ pðð2l þ 1Þði jÞ 0:5jÞ cos ¼ þ cos . 2 2M 2M Substituting these equations into Eq. (A.1) yields pffiffiffi X 1 pðð2l þ 1Þði þ jÞ 0:5jÞ pðð2l þ 1Þði jÞ þ 0:5jÞ 2cðiÞcðjÞ M1 e cos ðH s Þij ¼ þ cos 4 2M 2M M l¼0 pðð2l þ 1Þði þ jÞ þ 0:5jÞ pðð2l þ 1Þði jÞ 0:5jÞ þ cos þ cos . 2M 2M
ðA:1Þ
ðA:2Þ
ðA:3Þ
ðA:4Þ
ARTICLE IN PRESS Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
Here we also have the following two equations: pðði þ jÞð2l þ 1Þ 0:5jÞ pðði þ jÞð2l þ 1Þ þ 0:5jÞ cos þ cos 2M 2M pði þ jÞð2l þ 1Þ pj ¼ 2 cos cos , 2M 4M pðði jÞð2l þ 1Þ þ 0:5jÞ pðði jÞð2l þ 1Þ 0:5jÞ cos þ cos 2M 2M pði jÞð2l þ 1Þ pj ¼ 2 cos cos . 2M 4M Substituting these equations into (A.4) gives M1 X pði þ jÞð2l þ 1Þ e s Þij ¼ cðiÞcðjÞ pffiffiffi ðH cos 2M 2M l¼0 pði jÞð2l þ 1Þ pj þ cos cos . 2M 4M
2379
ðA:5Þ
ðA:6Þ
ðA:7Þ
Here we note that
8 k ¼ 4mM; > < M; pkð2l þ 1Þ cos ¼ M; k ¼ ð4m þ 2ÞM; > 2M : l¼0 0 otherwise: e s are Therefore, the ij elements of H 8 1 pj > > pffiffiffi cos ; i j ¼ 0; > > > 4M < 2 e s Þij ¼ ðH 1 pj > ; i þ j ¼ 2M; pffiffiffi cos > > 4M 2 > > : 0 otherwise: M1 X
(A.8)
(A.9)
This is expressed in the matrix form as (26).
References [1] S.C. Park, M.K. Park, M.G. Kang, Super-resolution image reconstruction—A technical overview, Signal Process. Mag. 20 (3) (2003) 21–36. [2] S. Rhee, M.G. Kang, Discrete cosine transform based regularized high-resolution image reconstruction algorithm, Opt. Eng. 38 (8) (1999) 1348–1356. [3] H. Sakane, K. Nishikawa, H. Kiya, A super-resolution method based on the discrete cosine transforms, Eusipco’96, 1996. [4] R.C. Hardie, K.J. Barnard, J.G. Bognar, E.E. Armstrong, E.A. Watson, High-resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system, Opt. Eng. 37 (1) (January 1998) 247–260. [5] J.H. Shin, J.H. Jung, J.K. Paik, Regularized iterative image interpolation and its application to spatially scalable coding, IEEE Trans. Consum. Electron. 44 (3) (August 1998) 1042–1047.
[6] M. Unser, A. Aldroubi, M. Eden, Enlargement or reduction of digital images with minimum loss of information, IEEE Trans. Image Process. 4 (3) (March 1995) 247–258. [7] M. Unser, A. Aldroubi, M. Eden, Polynomial spline signal approximations: filter design and asymptotic equivalence with Shannon’s sampling theorem, IEEE Trans. Inform. Theory 38 (1) (1992) 95–103. [8] P. The´venaz, T. Blu, M. Unser, Image interpolation and resampling, in: I.N. Bankman (Ed.), Handbook of Medical Imaging, Processing and Analysis, Academic Press, San Diego CA, USA, 2000, pp. 393–420. [9] S.A. Martucci, Image resizing in the discrete cosine transform domain, in: ICIP’95, vol. 2, 1995, pp. 224–227. [10] E. Shinbori, M. Takagi, High quality magnification applying the Gerchberg–Papoulis iterative algorithm with DCT, Trans. IEICE J76-D (9) (September 1993) 1932–1940. [11] S.E. El-Khamy, M.M. Hadhoud, M.I. Dessouky, B.M. Salam, F.E. Abd El-Samie, Efficient implementation of image interpolation as an inverse problem, Digit. Signal Process. 15 (2005) 137–152.
ARTICLE IN PRESS 2380
Y. Abe, Y. Iiguni / Signal Processing 87 (2007) 2370–2380
[12] M. Ng, N. Bose, Mathematical analysis of super-resolution methodology, IEEE Signal Process. Mag. 20 (3) (May 2003) 62–74. [13] L. Chen, K. Yap, Regularized interpolation using Kronecker product for still images, in: ICIP 2005, vol. 2, 2005 pp. II1014–1017. [14] A. Rosenfeld, A.C. Kak, Digital Picture Processing, Academic Press, New York, 1976. [15] K. Horiike, Y. Yoshida, K. Fujita, Image restoration based on a Wiener filtering technique using discrete cosine transform, Syst. Comput. Japan (SCJAEP) 27 (13) (1996) 23–35. [16] S.E. El-Khamy, M.M. Hadhoud, M.I. Dessouky, B.M. Salam, F.E. Abd El-Samie, Sectioned implementation
of regularized image interpolation, in: Proceedings of the 46th IEEE International Midwest Symposium on Circuits and Systems (MWSCAS ’03), vol. 2, 2003, pp. 656–659. [17] B.R. Hunt, The application of constrained least squares estimation to image restoration by digital computer, IEEE Trans. Comput. C-22 (9) (1973) 805–812. [18] I.Y. Soon, S.N. Koh, C.K. Yeo, Noisy speech enhancement using discrete cosine transform, Speech Comm. 24 (1998) 249–257. [19] F. O’Sullivan, Discretized Laplacian smoothing by Fourier methods, J. Amer. Statist. Assoc. 86 (September 1991) 634–642.