A Block-Based Regularized Approach for Image Interpolation

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.Theproposed 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 squaresminimization, 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.


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][2][3][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][7][8][9][10].The edge-directed interpolation is tuned based on the covariance.Other methods, including fuzzy-based [11], PDEbased [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 lowresolution 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 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 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.

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  (0, 0) = ℎ (−1, −1)  (−1, −1) + ℎ (−1, 0)  (−1, 0) + ℎ (−1, 1)  (−1, 1) + ℎ (0, −1)  (0, −1) + ℎ (0, 0)  (0, 0) + ℎ (0, 1)  (0, 1) where ℎ(, ) is the weight that is proportional to the area of HR pixels (, ) in the LR pixel (0, 0).In this case, ℎ(−1, −1) = (1 − Δ)(1 − Δ) and ℎ(−1, 0) = (1 − Δ).ℎ(, ) 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 h = sz  , where Let (, ) and (, ) represent the observed LR image of size  ×  and the original HR image of size  × , respectively, where  is the decimation factor.The 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 f and g are the lexicographically ordered HR and LR images, respectively, n is the noise vector, and H and D 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 f given the observation g.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, f opt 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 f opt = [(DH)  (DH) + I] −1 (DH)  g.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 . In view of this, we propose a computationally efficient method to address this issue in the paper.

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.
The structure of the Kronecker product for the decimation matrix D in (3) depends on the decimation factor, .For instance, the structure of D when  = 2 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 DH 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 DH can be expressed as Next, we apply SVD on D 1 H 1 and D 2 H 2 to obtain where , and Σ 2 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 (U 1 ⊗ U 2 )  and multiply the second row of the matrix by ] Finally, we denote y g, to rewrite the above equation to obtain the optimal solution y opt as To find out the closed-form solution of ( 13), we first get , since the Kronecker product for matrices A, B, and vector c yields the property that (A ⊗ B)c = vec(Bivec(c)A  ).Then, according to (13), we can get 2 + ‖y‖ 2 2 and use pseudoinverse as the optimal solution as where . Finally, we can rewrite y opt 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 ( 3  3 ) operations are needed in the new scheme as compared to ( 6  6 ) for direct implementation of pseudoinverse for (4).Therefore, this makes the algorithm computationally attractive and feasible to be implemented in real-life applications.

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 f opt = P()g, where f opt is the estimated HR image and P() = ((DH)  (DH) + I) −1 (DH)  .Substituting ( 9) and ( 10) into DH and P(), we obtain where respectively, which can be expressed as where k1, , k2, , ũ1, , and ũ2, are the column orthonormal basis of V 1 , V 2 , U 1 , and U 2 , respectively.The residual error vector between the original image f and the interpolation result f opt is equal to If the noise is zero-mean additive white Gaussian with variance  2  and independent of the image, then the expected error vector and square error between f opt and f 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., (r The mean square error (MSE) (‖r‖ 22 ) 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 (ff  ) =  2   I, where  2  = (f  f)/ 2 is the image power.This assumption is used to enable (‖r‖ 2  2 ) 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

Comparison with Other Image Interpolation Approaches.
The proposed approach is compared with other conventional methods: Lagrange (2nd-order polynomial) and bicubic and The image is divided into blocks of size 16 × 16 with 4 × 4 overlapping area.There is a 4 × 4 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 10 −4 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.0dB 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.

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 10 −2 to 10 −5 .In addition, the estimated  based on the GCV function in ( 22) is also adopted.Table 2 summarizes the results   is much faster than the edge-directed interpolation, faster than the bicubic interpolation, but slightly slower than the Lagrange interpolation.

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.

Figure 1 :
Figure 1: The relationship between HR and LR pixels.

Table 1 :
PSNR performance comparison of various image interpolation approaches.andCalTrain(400 × 512), as shown in Figure2.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  = 2 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.

Table 2 :
Comparison of different regularization parameters in PSNR performance.