Joint 2D SSA with GSVD and binary linear programming based image super-resolution

- This paper proposes a joint two dimensional (2D) singular spectrum analysis (SSA) with the generalized singular value decomposition (GSVD) and the binary linear programming based method for performing the super-resolution. For a given low resolution image, first both the upsampling operation and a lowpass filtering are applied on each column of the image to obtain an enlarged image. Second, apply the 2D Hankelization to both the low resolution image and the enlarged image to obtain their corresponding trajectory matrices. Third, both the GSVD and the 2D de-Hankelization are applied to these two trajectory matrices to obtain their corresponding sets of the de-Hankelized 2D SSA components. Here, it is proved that the exact perfect reconstruction is achieved. In order to enhance the high frequency contents of the enlarged image, the selection of the de-Hankelized 2D SSA components is formulated as a binary linear programming problem. Computer numerical simulation results show that the proposed method outperforms the state of art methods.


Introduction
Super-resolution images are the high resolution images constructed using the low resolution images [1]. The low resolution images are usually obtained due to the limitations of the acquisition hardware [1]. In this case, the acquisition process can be modelled by the downsampling operation [2]. Since the downsampling operation results to the lost of the information and the occurrence of the aliasing, in general the downsampling operation is irreversible [3]. Hence, it is very challenges to reconstruct the high resolution images using these low resolution images. However, if the reconstruction error is suppressed, then the details of the images can be seen more clearly. This result is very useful to the image forensics and medical diagnosis [4].
To address the above difficulty, the low resolution images are first upsampled. Then, the upsampled images are processed via the lowpass interpolation filters [5]. Since the upsampling operation and the lowpass filtering are the linear operators, the obtained performances are limited by these linear operators. On the other hand, the deep learning approach [6] is proposed. Although the deep learning approach is a kind of nonlinear methods, the required computational power is very high.
The outline of this paper is as follows. Section 2 presents our proposed method based on the joint 2D SSA with the GSVD and binary linear programming. Section 3 presents the computer numerical simulation results. Finally, a conclusion is drawn in Section 4.

Method
This paper aims to perform a nonlinear and adaptive post processing to the upsampled and filtered images to reduce the reconstruction error. The research work is designed as follows. First, the trajectory matrix is constructed via the 2D Hankelization [7]. The details are discussed in Section 2.1. Second, the GSVD [8] is employed to decompose the trajectory matrices to obtain the 2D SSA components of both the low resolution images as well as the corresponding upsampled and filtered images. The details are discussed in Section 2.2. Third, the 2D SSA components are selected via a binary linear programming approach [9]. The details are discussed in Section 2.4. Finally, the post processing images are obtained via the 2D de-Hankelization. The details are discussed in Section 2.3.

2D Hankelization
Let a low resolution image be Y and the size of Y be 1 2 L L ×   . That is, 1 2 L L × ∈ ℜ Y   . First, Y is divided into some blocks. Let the sizes of these blocks be 1 Here, the superscript " T " denotes the transposition operator. Obviously, 1 2 , L L m n ∈ ℜ y  for Let the trajectory matrix be Y  . Now, the trajectory matrix is constructed by putting m Y horizontally for That is, Let the upsampled and filtered image be Z and the upsampling ratio be S . It is assumed that 1 S > as well as both 1 SL and 1 SL  are integers, but S is not necessary to be an integer. Obviously, the size of Z is Similarly, Z is also divided into some blocks.
Denote the sizes of these blocks as Similarly, the vectorization on  . That is, Instead of performing the singular value decomposition (SVDs) of Y  and Z  individually, this paper proposes to perform the GSVD on both Y  and Z  simultaneously. In particular, denote the superscript " H " as the conjugate transposition operator. Let a I be the a a × identity matrix as well as U and V be the 1   . That is,  as well as the column of X be l X for 1, , X l r =  . Define the Hankelized 2D SSA components as It is worth noting that the right decomposition matrices obtained via applying the GSVD to both Y  and Z  are X .
Hence, the Hankelized 2D SSA components of both Y  and Z  would be similar to each others. As a result, the super-resolution image would be close to the original image.
. Obviously, 1 2 , , Obviously, 1 , , . Define the de-Hankelized 2D SSA components as . Obviously, 1 2 , , and for Obviously, 1 , , It is worth noting that the diagonal averaging method guarantees the exact perfect reconstruction of the image if the conventional SVD is employed in the SSA and the blocks taken from both Y and Z are shifted by one pixel both horizontally and vertically. However, the exact perfect reconstruction condition is unknown if the GSVD is employed in the SSA as well as the blocks taken from both Y and Z are not shifted by one pixel vertically. To address these issues, we have the following results: this completes the proof.
 Theorem 1 reveals that the sum of the Hankelized 2D SSA components is equal to the trajectory matrix.
. On the other hand, since Besides, since . Moreover, as Likewise, since . On the other hand, since , ,1,1 , , This completes the proof.
 Theorem 2 reveals that the sum of the de-Hankelized 2D SSA components is equal to the original image.

Selection of de-Hankelized 2D SSA components via binary linear programming
It is worth noting that the size of Z is larger than that of Y , so it is difficult to define the objective function in the pixel domain to select l Z  for 1, , Z l r =  . On the other hand, the selection of is performed in the frequency domain so that the high frequency contents of Z can be enhanced.
Since l Z  for 1, , Z l r =  are processed in the block based manner and some spatial regions in Z contain more information than other spatial regions, the selection of l Z  for 1, , Z l r =  is also performed in the block based manner. This can provide a higher flexibility for reconstructing Z .
For the conventional SSA, the de-Hankelized 2D SSA components are either selected or dropped. Hence, the selection problem is formulated as a binary programming involves the 0 L norm operator which is highly nonconvex, in general there is more than one global optimal solution. In this case, the global optimal solution is not unique. To avoid this difficulty, the convex relaxation is performed. We define ( ) , ω ω π π π π ∈ − × − , respectively. Define  , , , ω ω π π π π ∈ − × − . This is equivalent to

Experimental results
This work does not involve any human participant as well as any data from the human tissue and the animal. A set of color images with various types of contents is employed for the demonstration. In particular, the set of images include the pictures of a building, a toy, a natural landscape, an airplane, a tree and a baby. Here, these are the low resolution images. Therefore, their sizes are very small. In particular, the sizes of the image "building", the image "toy", the image "natural landscape", the image "airplane", the image "tree", the image "baby" are chosen as 1  , respectively. Figure 1 shows some of these low resolution images. baby. To perform the super-resolution, the total numbers of the rows of the high resolution images are chosen as the double of those of the low resolution images. That is, 2 S = . To generate the enlarged images, the triangular interpolation is employed. That is, the values of the elements in the inserted rows are the average values of the elements of the rows before and after the inserted rows.
In the algorithm, the block size is chosen as 1 2 2 2 L L × = × . Also, the frequency domain is sampled using 225 points. That is, 225 K = . Moreover, the specifications on the upper bounds of the real parts and the imaginary parts of the constraint functions are set as 70 ε = . In order to compare our proposed method, a similar nonlinear adaptive approach is compared. In particular, the empirical mode decomposition (EMD) based method [10] is compared. Here, the problem of finding the high resolution image is formulated as an optimization problem in a similar manner. Figure 2 and Figure 3 show the super-resolution images obtained by our proposed approach and the EMD approach [10], respectively. It can be seen that our proposed approach can achieve the better super-resolution performances qualitatively.

Discussion
Since the proposed method yields an image with a higher resolution, the proposed method can be applied to the medical images so that the medical officers have the medical images with the higher resolutions to improve the accuracy of the diagnosis. However, as the proposed method requires the upsampled and the filtered images, it depends on the employed interpolation filter. To address this issue, the nonlinear and adaptive approaches such as the two dimensional singular spectrum analysis will be employed to generate the images with the larger sizes.

Conclusions
This paper proposes a joint 2D SSA with the GSVD and the binary linear programming based method for performing the super-resolution. In particular, the exact perfect reconstruction condition is derived when the GSVD is employed in the SSA as well as the blocks taken from both the original image and the enlarged image are not shifted by one pixel vertically. Moreover, as the Hankelized 2D SSA components of both the original image and the enlarged image would be similar to each others, the super-resolution image would be close to the original image. List of abbreviations 2D Two dimensional SSA Singular spectrum analysis GSVD Generalized singular value decomposition SVD Singular value decomposition EMD Empirical mode decomposition Availability of data and material All the images can be downloaded from the public domain.

Competing interests
There is no conflict of interest.

Authors' contributions
Ziyin Huang is responsible for the conceptualization, deriving the theory, implementation of the Matlab program and writing the draft paper. Bingo Wing-Kuen Ling is responsible for finding a funding to support this research, verification of the developed theory and revising the draft paper. Yui-Lam Chan and Huan Ye are responsible for checking the results.