Research Article  Open Access
A BlockBased Regularized Approach for Image Interpolation
Abstract
This paper presents a new efficient algorithm for image interpolation based on regularization theory. To render a highresolution (HR) image from a lowresolution (LR) image, classical interpolation techniques estimate the missing pixels from the surrounding pixels based on a pixelbypixel basis. In contrast, the proposed approach formulates the interpolation problem into the optimization of a cost function. The proposed cost function consists of a data fidelity term and regularization functional. The closedform solution to the optimization problem is derived using the framework of constrained least squares minimization, by incorporating Kronecker product and singular value decomposition (SVD) to reduce the computational cost of the algorithm. The effect of regularization on the interpolation results is analyzed, and an adaptive strategy is proposed for selecting the regularization parameter. Experimental results show that the proposed approach is able to reconstruct highfidelity HR images, while suppressing artifacts such as edge distortion and blurring, to produce superior interpolation results to that of conventional image interpolation techniques.
1. Introduction
Image interpolation aims to render a highresolution (HR) image by estimating the missing pixel information from the surrounding pixels of the observed lowresolution (LR) image. This image processing task is commonly referred to as image interpolation (when only a single LR image is available) or superresolution (when multiple LR images are available) [1–4]. It has a wide range of applications such as computational photography and video surveillance.
The challenge of image interpolation is to enhance the image quality while suppressing artifacts such as blurring and jagged edges. Various algorithms have been developed to address image interpolation over the years. Conventional piecewise polynomial techniques, such as bicubic interpolation, assume higherorder continuity of image intensity in the spatial domain. These techniques, however, often lead to oversmoothness in the edge and textured regions. Some edgedirected approaches adjust the algorithmic parameters according to the local features. The HR local covariance coefficients are estimated from the LR counterpart based on their geometric duality [5]. For example, inverse gradient has been employed to determine the weights of bicubic interpolation. Different edge types are identified and used to determine different interpolation strategy for local image area [6–10]. The edgedirected interpolation is tuned based on the covariance. Other methods, including fuzzybased [11], PDEbased [12], and regressionbased method [13], have also been reported in color image interpolation and superresolution. These adaptive methods typically employ heuristic reasoning to estimate parameters such as threshold values or filter weights on a pixelbypixel basis. Therefore, they require extra computation to determine these local parameters, and the quality of the interpolated images may vary significantly with respect to changes in these parameters. The aforementioned deterministic image interpolation approaches cannot suppress the noise or blurring incurred in the input lowresolution image.
In view of this, a blockbased regularized image interpolation approach is proposed in this paper by imposing the regularization constraint on the reconstructed highresolution image. The relationship between the HR and LR images is exploited to formulate the interpolation problem as the optimization of a cost function. The cost function consists of a data fidelity term and Tikhonov regularization functional [14]. The closedform solution to the optimization problem is estimated using a new framework of constrained least squares minimization, Kronecker product, and singular value decomposition (SVD). A key feature of the method is its computational efficiency in reconstructing highfidelity HR image, while alleviating common artifacts encountered by other interpolation techniques. This allows the proposed method to be employed readily in the areas of digital photography, computer vision, and medical imaging, among others.
The rest of this paper is organized as follows. Section 2 provides the mathematical formulation on the image interpolation. In Section 3, the proposed regularized image interpolation approach is presented, which is further evaluated in experimental results presented in Section 4. Finally, Section 5 concludes this paper.
2. Problem Formulation
The mathematical model between the HR and LR images plays a key role in the formulation of interpolation algorithms. Due to the finite sampling grid of the sensor array, the image acquisition processes can be modeled as integration (averaging), followed by decimation. Decimation and interpolation are inverse processes. An image decimation/interpolation model is illustrated in Figure 1. The dotted squares denote the HR pixels, while the shaded solid square represents the LR pixel. The decimation factor is 2 in this example.
From a physical viewpoint of image acquisition, the response of each sensor is proportional to the integral of the light projected onto the surface of the sensor. Therefore, the intensity of the LR pixel in Figure 1 is determined by the corresponding effective area in the HR image grids. Consider where is the weight that is proportional to the area of HR pixels in the LR pixel . In this case, and . is a 2D integration operator, which reflects the regions of the HR grid that contribute to the formation of a single pixel in the LR grid. This process is commonly expressed as , where (if is even) or (if is odd).
Let and represent the observed LR image of size and the original HR image of size , respectively, where is the decimation factor. The relationship between the LR and HR images can be modeled as where is the intermediary signal before decimation, is a twodimensional (2D) integration operator characterizing the averaging process from the HR to LR image, and is the additive noise. Rewriting (2) in the matrixvector formation, we have where and are the lexicographically ordered HR and LR images, respectively, is the noise vector, and and are the corresponding matrices constructed from the integration and decimation processes [15, 16]. The interpolation problem can, therefore, be formulated as solving the least squares problem for given the observation . The Tikhonov regularization framework is used in this paper to address this problem, as it is able to offer numerically stable and visually pleasing solution. Under this setting, is the solution to where is the regularization parameter which controls the relative contributions between the least square error (first term) and the regularization functional (second term). The regularization term is to ensure smoothness of the solution. The closedform solution to the least squares problem in (4) using pseudoinverse is given by . However, the closedform solution is impractical in many reallife applications due to the high computational cost associated with the inversion of the large matrix (with a dimension of ). In view of this, we propose a computationally efficient method to address this issue in the paper.
3. Proposed Regularized Image Interpolation Approach
The formulation of the proposed interpolation algorithm involves the manipulation of Kronecker product and SVD to reduce computational cost. The regularized least square solution can achieve a good tradeoff between edge preservation and noise suppression.
3.1. Proposed New Framework of Image Interpolation through SVD and Kronecker Product
The matrix in the linear image acquisition model (3) has a blockcirculant structure [15]: where is the circulant matrix constructed from the 2D integration operator . In the context of image interpolation problem, the integration operator is usually characterized as a 2D separable averaging process. Therefore, the matrix can be decomposed into the Kronecker product of two matrices as where denotes the Kronecker product. Here both and are circulant matrices, whose entries are constructed from the vector decomposition of integration operator: , where and .
The structure of the Kronecker product for the decimation matrix in (3) depends on the decimation factor, . For instance, the structure of when is given as Other decimation matrices with different decimation factors can be constructed in a similar fashion to (7) by adjusting the spacing between the consecutive 1s.
The minimization of (4) can be expressed as the constrained least squares problem via Since the dimension of is very large, it may not be computationally feasible to solve the problem directly using matrix pseudoinverse. In view of this, we develop an efficient scheme to solve (8) based on Kronecker product and SVD as follows.
First, based on (6) and (7), the Kronecker product of can be expressed as Next, we apply SVD on and to obtain where , , , , , and are the standard matrices arising from SVD. Then, substituting (9) and (10) into (8), we obtain We further multiply the first row of the matrix by and multiply the second row of the matrix by to get Finally, we denote , , to rewrite the above equation to obtain the optimal solution as
To find out the closedform solution of (13), we first get and , since the Kronecker product for matrices , , and vector yields the property that . Then, according to (13), we can get and use pseudoinverse as the optimal solution as where , , and . Finally, we can rewrite as By utilizing the newly derived equation (15), the computational cost of the estimation can be reduced. For an image of size with decimation factor , only operations are needed in the new scheme as compared to for direct implementation of pseudoinverse for (4). Therefore, this makes the algorithm computationally attractive and feasible to be implemented in reallife applications.
3.2. Regularization Analysis
After establishing the efficient algorithm to obtain the regularized solution, we study the effect of regularization on the interpolation results. The mean error and mean square error between the interpolated and original images are analyzed as follows.
The direct solution to the least square problem in (4) is given as , where is the estimated HR image and . Substituting (9) and (10) into and , we obtain where , . and are the th and th column vectors of the orthogonal matrices and , respectively, which can be expressed as where are the column orthonormal basis of , respectively.
The residual error vector between the original image and the interpolation result is equal to If the noise is zeromean additive white Gaussian with variance and independent of the image, then the expected error vector and square error between and are given by It can be observed from (19) that the regularization functional has biased the interpolated image, as a tradeoff for numerical stability (i.e., if since ). The mean square error (MSE) consists of two parts: the noise term and the image term. If the image is preprocessed using zeromean centering and assumed to be white, then , where is the image power. This assumption is used to enable to be upperbounded explicitly in terms of as below: It can be seen from (19) and (20) that regularization will reduce the impact of noise by introducing a small bias in the interpolated image.
The regularization technique is instrumental in providing satisfactory results so long as the regularization parameter is properly selected. A simple yet effective scheme to estimate the regularization parameter is based on the generalized crossvalidation (GCV) function since it has been shown to be robust in other regularization schemes [16]: Crossvalidation is a widely used technique in the field of statistical data analysis. It uses the leaveoneout principle to address approximation of noisy data. For a fixed value of model parameter, an interpolated image is redecimated to predict the observation. GCV determines the parameter by minimizing the weighted sum of prediction errors. An advantage of GCV is that it allows the selection of the regularization parameter even when the noise power is unknown. We derive the GCV function under the context of our formulation to give where and are the singular values of and , respectively. is the th entry of the column vector in (15). In fact, experimental results in Section 4 show that our algorithm is robust towards different values of so long as it falls within a reasonable range.
4. Experimental Results
4.1. Comparison with Other Image Interpolation Approaches
The proposed approach is compared with other conventional methods: Lagrange (2ndorder polynomial) and bicubic and edgedirected image interpolation [5]. The standard benchmark test images are used as ground truth in the experiments, including Lena , Boat , Board , and CalTrain , as shown in Figure 2. The observed lowresolution is artificially generated by filtering the original ground truth highresolution image with an integration operator and decimating it at a rate of in the horizontal and vertical directions. The downsampled image is further degraded under different noise levels to produce different SNR values: noiseless, 40 dB, and 30 dB.
(a)
(b)
(c)
(d)
The image is divided into blocks of size with overlapping area. There is a overlapping in each block to avoid boundary effect. For each block, the proposed algorithm is applied to perform interpolation using (10)–(15), and the regularization parameter is chosen as in this experiment. For performance evaluation, peak signaltonoise ratio (PSNR) is used as the objective performance metric. Table 1 summarizes the results obtained for the test images using different methods. It can be observed that the proposed method achieves higher PSNR than the conventional methods consistently. On average, the proposed method offers a PSNR improvement of 0.5–1.0 dB over the Lagrange and bicubic methods and 12 dB over the edgedirected method. In Figures 3 and 4, the interpolated results using the Boat and Board images are given. It can be observed that the proposed method preserves the overall sharpness of the interpolated images, in particular near the edge and textured regions. The subjective human evaluation is confirmed by the objective performance measure given in Table 1.

(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4.2. Evaluation on the Choice of Regularization Parameter
In this section, the impact of the regularization parameter on the interpolation results is investigated. The Lena image in Figure 2(a) is selected as the test image. The same algorithm as in previous experiment is run but different values of regularization parameter are used, which range from to . In addition, the estimated based on the GCV function in (22) is also adopted. Table 2 summarizes the results obtained using different regularization parameters. It can be observed that this method provides consistently good results for to and for various noisy environments from noiseless to 30 dB noise. The estimated using GCV also offers comparable performance. The results suggest that the proposed algorithm is robust towards different noisy environments and regularization parameters as long as the value of falls within a reasonable range .

Note that there is an obvious improvement from to . This gap can be explained by (20) together with the simulated mean square error (MSE) curve in Figure 5. The upper bound of MSE in (20) consists of three terms, where the first and third terms are positive and the second is negative. When the value of increases, the second term will decrease; however, the third term will increase. We simulate this situation in terms of ranging from to in Figure 5. The image energy is normalized to 1 (i.e., ) and 30 dB SNR noise level (i.e., ). Consider the illconditioned problem of matrix inversion; half of the singular values are assumed to be evenly distributed in the range of and the other half are assumed to be 0. It can be observed from Figure 5 that the magnitude of second term drops quickly and dramatically before the third term begins to increase. This leads to a basin area in the overall mean square error when is between a certain range.
4.3. Evaluation on Computational Complexity
These image interpolation approaches are run on a PC with Windows XP, MATLAB 7.1, CPU P43.4 GHz, and 1 G RAM. Each approach is run for 10 times and their average runtime is presented in Table 3, where one can see that the proposed approach is much faster than the edgedirected interpolation, faster than the bicubic interpolation, but slightly slower than the Lagrange interpolation.

5. Conclusion
In this paper, a new regularizationbased image interpolation algorithm using Kronecker product and singular value decomposition has been proposed. The proposed approach reduces the computational cost of interpolation while offering significant performance improvement over other conventional methods. This paper also analyzes the effect of regularization on the interpolation results and shows that the proposed approach is fairly robust towards different values of regularization parameter. It is worthwhile to further study how to extend the proposed approach for the superresolution image reconstruction.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported by National Natural Science Foundation of China (nos. 61105010 and 61375017), Program for Outstanding Young Science and Technology Innovation Teams in Higher Education Institutions of Hubei Province, China (no. T201202).
References
 J. Zhao, H. Feng, Z. Xu, and Q. Li, “An improved image deconvolution approach using local constraint,” Optics and Laser Technology, vol. 44, no. 2, pp. 421–427, 2012. View at: Publisher Site  Google Scholar
 J. Tian and K.K. Ma, “A survey on superresolution imaging,” Signal, Image and Video Processing, vol. 5, no. 3, pp. 329–342, 2011. View at: Publisher Site  Google Scholar
 P. Shanmugavadivu and K. Balasubramanian, “Particle swarm optimized multiobjective histogram equalization for image enhancement,” Optics and Laser Technology, vol. 57, pp. 243–251, 2014. View at: Publisher Site  Google Scholar
 H. Yang, Z. Zhang, M. Zhu, and H. Huang, “Edgepreserving image deconvolution with nonlocal domain transform,” Optics and Laser Technology, vol. 54, pp. 128–136, 2013. View at: Google Scholar
 X. Li and M. T. Orchard, “New edgedirected interpolation,” IEEE Transactions on Image Processing, vol. 10, no. 10, pp. 1521–1527, 2001. View at: Publisher Site  Google Scholar
 J. W. Hwang and H. S. Lee, “Adaptive image interpolation based on local gradient features,” IEEE Signal Processing Letters, vol. 11, no. 3, pp. 359–362, 2004. View at: Publisher Site  Google Scholar
 D. D. Muresan and T. W. Parks, “Adaptively quadratic (AQua) image interpolation,” IEEE Transactions on Image Processing, vol. 13, no. 5, pp. 690–698, 2004. View at: Publisher Site  Google Scholar
 L. Wang, S. Xiang, G. Meng, H. Wu, and C. Pan, “Edgedirected singleimage superresolution via adaptive gradient magnitude selfinterpolation,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 23, no. 8, pp. 1289–1299, 2013. View at: Google Scholar
 H. Xu, G. Zhai, and X. Yang, “Single image superresolution with detail enhancement based on local fractal analysis of gradient,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 23, no. 10, pp. 1740–1754, 2013. View at: Google Scholar
 F. Zhou, W. Yang, and Q. Liao, “Interpolationbased image superresolution using multisurface fitting,” IEEE Transactions on Image Processing, vol. 21, no. 7, pp. 3312–3318, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 V. Magudeeswaran and C. G. Ravichandran, “Fuzzy logicbased histogram equalization for image contrast enhancement,” Mathematical Problems in Engineering, vol. 2013, Article ID 891864, 10 pages, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 T. Liu and Z. Xiang, “Image restoration combining the secondorder and fourthorder PDEs,” Mathematical Problems in Engineering, vol. 2013, Article ID 743891, 7 pages, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 X. Liu, D. Zhao, R. Xiong, S. Ma, W. Gao, and H. Sun, “Image interpolation via regularized local linear regression,” IEEE Transactions on Image Processing, vol. 20, no. 12, pp. 3455–3469, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 A. N. Tikhonov, Equations of Mathematical Physics, vol. 39, Dover Publications, New York, NY, USA, 1963.
 R. L. Lagendijk and K. Biemond, Iterative Identification and Restoration of Images, Kluwer Academic Publishers, Boston, Mass, USA, 1990.
 N. P. Galatsanos and A. K. Katsaggelos, “Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation,” IEEE Transactions on Image Processing, vol. 1, no. 3, pp. 322–336, 1992. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2014 Li Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.