Abstract
This paper presents a new efficient algorithm for image interpolation based on regularization theory. To render a high-resolution (HR) image from a low-resolution (LR) image, classical interpolation techniques estimate the missing pixels from the surrounding pixels based on a pixel-by-pixel 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 closed-form 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 high-fidelity 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 high-resolution (HR) image by estimating the missing pixel information from the surrounding pixels of the observed low-resolution (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 higher-order continuity of image intensity in the spatial domain. These techniques, however, often lead to oversmoothness in the edge and textured regions. Some edge-directed 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 edge-directed interpolation is tuned based on the covariance. Other methods, including fuzzy-based [11], PDE-based [12], and regression-based 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 pixel-by-pixel 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 low-resolution image.
In view of this, a block-based regularized image interpolation approach is proposed in this paper by imposing the regularization constraint on the reconstructed high-resolution 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 closed-form 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 high-fidelity 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 two-dimensional (2D) integration operator characterizing the averaging process from the HR to LR image, and is the additive noise. Rewriting (2) in the matrix-vector 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 closed-form solution to the least squares problem in (4) using pseudoinverse is given by . However, the closed-form solution is impractical in many real-life 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 block-circulant 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 closed-form 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 real-life 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 zero-mean 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 zero-mean centering and assumed to be white, then , where is the image power. This assumption is used to enable to be upper-bounded 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 cross-validation (GCV) function since it has been shown to be robust in other regularization schemes [16]: Cross-validation is a widely used technique in the field of statistical data analysis. It uses the leave-one-out 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 (2nd-order polynomial) and bicubic and edge-directed 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 low-resolution is artificially generated by filtering the original ground truth high-resolution 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 signal-to-noise 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 1-2 dB over the edge-directed 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 ill-conditioned 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 P4-3.4 GHz, and 1 G RAM. Each approach is run for 10 times and their average run-time is presented in Table 3, where one can see that the proposed approach is much faster than the edge-directed interpolation, faster than the bicubic interpolation, but slightly slower than the Lagrange interpolation.
5. Conclusion
In this paper, a new regularization-based 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).