Compact spatio-spectral algorithm for single image super-resolution in hyperspectral imaging

Hyperspectral imaging (HSI) is used in a wide range of applications such as remote sensing, space imagery, mineral detection, and exploration. Unfortunately, it is difficult to acquire hyperspectral images with high spatial and spectral resolution due to instrument limitations. The super-resolution techniques are used to reconstruct low-resolution hyperspectral images. However, traditional superresolution (SR) approaches do not allow direct use of both spatial and spectral information, which is a decisive for an optimal reconstruction. This paper proposes a single image SR algorithm for HSI. The algorithm uses the fact that the spatial and spectral information can be integrated to make an accurate estimate of the high-resolution HSI. To achieve this, two types of spatio-spectral downsampling, and a three-dimensional interpolation are proposed in order to increase coherence between the spatial and spectral information. The resulting reconstructions using the proposed method are up to 2 dB better than traditional SR approaches.


Introduction
Hyperspectral imaging (HSI) collects a concatenation of bidimensional images that entails different wavelengths in a certain spectral range. Pixels in hyperspectral images are therefore represented by vectors whose entries correspond to the intensity in the different spectral bands. HSI enables the detection, classification, and identification of objects and features based on the spectral characteristics (Chakrabarti, 2011). HSI is an area with a significant impact in civilian and military applications including remote sensing, aerial, space imagery, natural resource exploration, farming, and astronomy (Belluco, 2006), (Borengasser, 2007), (Castrodad, 2010), (Melgani, 2004), (Underwood, 2003), (Dicker, 2006), (Turk, 1991). In all of these applications, it is important to obtain the highest resolution in the spatial and spectral dimensions. Typically, the hyperspectral spectrometers are used to capture high-resolution hyperspectral images, because these provide hundreds of narrow contiguous bands over a wide range of the electromagnetic spectrum. Hyperspectral sensors measure the reflective properties of objects in the visible and short wave infrared regions of the spectrum. Unfortunately, atmospheric scattering, imperfect imaging optics, secondary illumination, changing viewing angles, and sensor noise degrade the quality of these images, making the spatial resolution one of the most expensive and hardest characteristic to be improved in imaging systems.
In practice, modifying the imaging optics or sensor array increases the production costs (Arguello H. , 2012)- (Rueda, 2013); in addition it is impossible to fabricate spectrometers with arbitrary resolution. Hence, researchers have investigated the use of resolution enhancement techniques by post processing as a better alter native to improve the image quality (Akgun, 2005), (Mianji, 2008), (Zhao, 2011). To improve the spatial resolution of hyperspectral images, traditional super-resolution (SR) techniques may be used. Super-resolution methods are currently a very active area of research, as it allows the implementation of low-cost imaging sensors in hyperspectral spectrometers. Many methodologies have been applied to the super-resolution problem, such as SR from single image methods (Zhang, 2012), (Takeda, 2007), (Li, 2001), multi-image-based SR methods (Protter, 2009), (Akgun, 2005), (Chan, 2010), and learning-based SR methods (Freeman, 2002). However, these methods ignore the spectral information that is a crucial parameter for an optimal reconstruction of high resolution (HR) images. A comprehensive background on super-resolution methods can be found in (Akgun, 2005).

Bidimensional image super-resolution
The single image super-resolution (SR) problem consists on recovering a high resolution (HR) image X ∈ ! s 1 N×s 2 M from a low resolution (LR) Y ∈ ! N×M version. Figure 1 shows an example of an HR image restoration from an LR version of it. Here, N × M represents the spatial resolution, and s 1 ,s 2 ∈ ! are the super-resolution factors in the dimensions N and M, respectively. Let X ∈ ! s 1 s 2 NM×1 and y ∈ ! NM×1 be the vector representations of a HR and its LR image version, respectively (Zhao, 2011), (Winter, 2002), (Tanaka, 2007). Then the sensed image can be expressed as where H∈ ! s 1 s 2 NM×s 1 s 2 NM denotes the blurring matrix, D ∈ ! NM×s 1 s 2 NM is the downsampling matrix, and ε ∈ ! NM×1 is the noise introduced by the sensing system. Examples of the structure of the matrices D and H are indicated in Figure 2. Since, NM ≪ s 1 s 2 NM , Equation (1) leads to an undetermined system of linear equations, which has infinite number of solutions on x, i.e. an ill-posed problem. To make the image recovery process less ill-posed (Akgun, 2005), (Dong, 2011), Equation (1), can be rewritten as the least squares formulation (2) Solution to Equation (2) is given by however, this method leads to a low quality spatial reconstruction image. To increase quality, several SR methods incorporate an effective prior (denoted as a regularization term) into the reconstruction process in order to improve the above solution. Accordingly, Equation (2) is reformulated as This paper focuses on the problem of recovering the superresolution version of a given low-resolution hyperspectral image. Specifically, this work develops a novel reconstruction approach by which the spatio-spectral information of lowresolution (LR) HSI input is exploited. We named this novel super-resolution algorithm the C2SR algorithm. We suggest a fast algorithm to integrate the spatial and spectral information of an HSI to exploit the spatial resolution of the image. Moreover, C2SR can work on HSIs captured using current hyperspectral spectrometers. This becomes an advantage over other techniques like (Rueda-Chacón, 2014), (Arguello H. &., 2012), which require the re-design of spectrometers based on compressive sensing approaches. (a)

MARQUEZ, VARGAS, AND ARGUELLO
Example of a downsampling matrix D H ∈ ! 54×216 , and a blurring matrix H H ∈ ! 216×216 for a hyperspectral image X ∈ ! 3s 1 ×3s 2 ×6s 3 , with a super-resolution factors of s 1 = s 2 = 2 , and s 3 = 1. (a) Downsampling matrix D H . (b) The sensing matrix H H whose effect is equivalent to a symmetric Gaussian lowpass filter of size 3 with standard deviation σ = 1. The dark and white pixels depicts values 0 and 1, respectively.
where F(x) is a regularization prior, and μ is a regularization parameter which represents the tradeoff between the reconstruction error and the regularization term. To evaluate the reconstruction performance of our method, we choose two traditional methods of SR. Non-local means (NLM), and steering kernel regression (SKR) are two important methods of super-resolution. These methods have received a substantial attention, being a family of methods based upon local smoothness assumption, i.e., the local structure is relatively stable (Takeda, 2007), (Protter, 2009).

Hyperspectral image super-Resolution
Hyperspectral imaging entails signals typically spanning hundreds of contiguous wavelength bands in a certain spectral range. Pixels in hyperspectral images are therefore represented as vectors whose entries correspond to the intensity in different spectral bands. Let X H ∈ ! s 1 N×s 2 M×s 3 L be a high-resolution hyperspectral image (HR-HSI), and Y H ∈ ! N×M×L its low-resolution (LR-HSI) version. Here, L represents the spectral resolution, and s 3 ∈ ! is the super-resolution factor in the dimension L. Also let x H ∈ ! s 1 s 2 s 3 NML×1 and y H ∈ ! NML×1 be the HR-HSI and LR-HSI vector representation of X H , and Y H respectively. The HSI acquisition process of y H from x H can be modeled as where D H ∈ ! NML×s 1 s 2 s 3 NML is the downsampling matrix in the spatial and spectral domain, H H ∈ ! s 1 s 2 s 3 NML×s 1 s 2 s 3 NML is a matrix that describes the blurring in each spectral band (no blurring across the spectrum), and ε H ∈ ! NML×1 is the noise introduced by the sensing system. Examples of the matrices D H and H H are indicated in Figure 2. The signal x H can be estimated by solving or approximated by using a gradient descent method such as Equation (5) is extremely ill-posed, considering that the product (D H H H ) T y H generates HR-HSI. The missing spectral bands are therefore set up with zeros; this effect is commonly known as zero-padding ( Figure 3). Thus, the recover signal x H t+1 is a hyperspectral image with zeropadding in its spectral field, and this is inconsistent with the requirements of high-resolution HSI. Figure 3 depicts an image with zero-padding in its spectral fields. ) by traditional methods using Equation (7).

The compact spatio-spectral super-resolution algorithm (C2SR)
This work revisits the HSI super-resolution problem, proposing the compact spatio-spectral super-resolution approach. This approach solves two main problems in current HSI super-resolution techniques: zero-padding, and waste of spectral information. In summary, three major steps can describe the proposed algorithm. First, an initial approximation of the super-resolved spectral bands is obtained by scaling the LR-HSI using a three-dimensional inter polation filter represented by Ψ ∈ ! s 1 s 2 s 3 NML×s 1 s 2 s 3 NML . This step avoids the zero padding effect by distributing information, con tained in the LR-HSI, among the super-resolved spectral bands. The second step consists on reducing the spatial resolution without changing the spectral dimension. This process is represented by the matrix D b ∈ ! s 3 NML×s 1 s 2 s 3 NML . This is carried on to keep the dimension of the HR-HSI search space, hence reducing the computational complexity. Also, keeping the dimension avoids image noise amplification. To make the Equation (5) well posed, it is rewritten as Result: Output of the final high-resolution HIS x H . Therefore, the SR reconstruction consists on solving the inverse problem in the Equation (8) to estimate the HR-HSI image x H . The super-resolution reconstruction is reduced to solving the least squared problem To obtain a solution to Equation (9), it is reformulated in a concise form using the gradient descent method, i.e., where t represents the iteration times, and τ stands for the step size for gradient descent. A detailed description of the C2SR algorithm is shown in Table I

Image acquisition
In this section, the performance of the proposed method is assessed in terms of the super-resolution factor, and the noise in the acquired data. For this purpose, the two traditional SR methods [NLM]-[SKR] are implemented and compared with the proposed C2SR algorithm. The SR methods were tested using the standard image databases used in (Arguello H. , 2012), (Rueda-Chacón, 2014) (Rueda, 2013). The first HSI is a Lego image with 512 × 512 pixels of spatial resolution and L = 24 spectral bands. The second is the Ribeira city image with 1000 × 1000 pixels of spatial resolution and L = 32 spectral bands. The Lego image was acquired using a wide-band Xenon lamp as the light source, and a visible monochromator which spans the spectral range between 450 nm and 650 nm. The image intensity was captured using a CCD camera AVT Marling F0033B, exhibiting 512 × 512 pixels, with pixel pitch of 9,9 μm, and 24 bits of pixel depth. The resulting test hyperspectral image x H has 512 × 512 pixels of spatial resolution and L = 24 spectral bands. The second, the Ribeira city image was acquired using a low-noise Peltier cooled digital camera providing an x − y spatial resolution of 1344 × 1024 pixels (Hamamatsu Photonics) with a fast tunable liquid-crystal filter mounted in front of the lens, together with an infrared blocking filter. The peak-transmission wavelength was varied in 10 nm steps over 400 −720 nm and the bandwidth was 10 nm at 510 nm, decreasing to 7 nm at 400 nm to 16 nm at 720 nm. The resulting test hyperspectral image x H has 1000 × 1000 pixels of spatial resolution and L = 32 spectral bands (Foster, 2004).

Experimental results
Simulations results were analyzed in terms of Peak-Signal-to-Noise-Ratio (PSNR) in the spatial, and spectral reconstructed images (Wang, 2004). PSNR is defined as 10log 10 max x H ( ) 2 / MSE ( ) wherein MSE is the mean squared error, and max (x H ) depicts the maximum intensity of the hyperspectral image. In order to illustrate the spectral reconstruction performance, two spatial points were randomly chosen, and the spectral signatures plotted in Figure 4, these points are indicated as P1 and P2. Again, it can be seen how the curves using the proposed method are closer to the original. It is important to remark the impact of the tridimensional filters in the reconstructed spectral signature curve for hyperspectral images.
The robustness of the reconstructions with respect to the effect of noise, and decimation of the measurements is studied in the Table 2-3, and illustrated in Figure 5 and Figure 6, respectively. Gaussian noise with zero mean was added to the measurements. The Signal to Noise Ratio (SNR) is calculated according to SNR = 10log 10 µ ν −y where μ v is the mean value of the LR-HSI y, and σ noise is the standard deviation of the signal noise. The factors of decimation used were 2,4 and 8, in the spatial and spectral dimension. The results using the proposed method indicate improvements up to 4dB in the PSNR of the reconstructed hyperspectral image compared with the results using Lanczos interpolation. The NLM and SKR methods achieve better results than Lanczos interpolation; however, these methods also have inferior performance than the proposed method. Therefore, in terms of the values of PSNR, the proposed method achieves higher performance.

Conclusions
In this paper, we present a novel SR reconstruction method for a single hyperspectral image. The proposed method incorporates the spectral reconstruction, the local spectral similarity, and the local structural regularity into a unified iterative framework for SR task. The thorough experimental results show the effectiveness of the proposed method. The proposed SR framework can be naturally extended by the following considerations: the design of three-dimensional downsampling matrix and three-dimensional interpolation matrix. The discrete mat ematical model of the LR image formation and HR reconstruction process for hyperspectral images has been proposed. This mathematical model outperforms current approaches up to 2 dB in terms of PSNR.