Comparison of our methods with various pansharpening methods on AF and AVIRIS.
We first briefly review recent papers for pansharpening of hyperspectral (HS) images. We then present a recent pansharpening approach called hybrid color mapping (HCM). A few variants of HCM are then summarized. Using two hyperspectral images, we illustrate the advantages of HCM by comparing HCM with 10 state-of-the-art algorithms.
- hyperspectral images
- hybrid color mapping
- image fusion
Hyperspectral (HS) images have found a wide range of applications in terrestrial and planetary missions. NASA is planning a HyspIRI mission [1–4] that will perform vegetation monitoring for the whole Earth. The spatial resolution of HyspIRI is 60 m. The Compact Reconnaissance Imaging Spectrometer for Mars (CRISM)  hyperspectral (HS) imager with 100-m resolution has been monitoring Mars surface since 2006. Although the above imagers have resolution good enough for their respective missions, many other applications such as drought monitoring, fire damage assessment, etc., require higher resolutions. Other notable applications of multispectral and hyperspectral images include target detection [6–13], anomaly and change detection [14–23], tunnel monitoring [24, 25], and Mars exploration [26, 27].
Pansharpening of hyperspectral images usually refers to the fusion of a high-resolution (HR) panchromatic (pan) band with a low-resolution (LR) hyperspectral image cube. A generalization of the above is the fusion of high-resolution multispectral bands with low-resolution hyperspectral bands. According to Loncan et al. , pansharpening techniques for HS images can be classified into the following categories. The first category is the component substitution (CS) approach. Well-known CS approaches include Principal Component Analysis (PCA) , Gram-Schmidt (GS) , GS Adaptive (GSA) , and others. The CS approach is based on the substitution of a component with the pan image. The second category is the multiresolution analysis (MRA) approach. The MRA approach relies on the injection of spatial details that are obtained through a multiresolution decomposition of the pan image into the resampled hyperspectral bands. Some well-known algorithms in this category are Modulation Transfer Function Generalized Laplacian Pyramid (MTF-GLP) , MTF-GLP with High-Pass Modulation (MTF-GLP-HPM) , Hysure [34, 35], and Smoothing Filter-based Intensity Modulation (SFIM) . The third category contains the hybrid approaches that use concepts from different classes of methods, namely from CS and MRA ones. For example, guided filter PCA (GFPCA)  belongs to this group. The fourth category involves Bayesian inference, which interprets the fusion process through the posterior distribution of the Bayesian fusion model. Due to the ill-posed fusion problem, the Bayesian methodology can easily regularize the problem by defining a prior distribution for the scene of interest. Exemplar algorithms include Bayesian naive  and Bayesian sparse . The fifth category is known as non-negative matrix Factorization (NMF). Representative methods in this category include the coupled non-negative matrix factorization (CNMF)  method.
In addition to the above five categories, we noticed that deep learning-based approaches have been investigated in recent years to address pansharpening for hyperspectral images. In Licciardi et al. , the authors proposed the first deep learning-based fusion method, which used autoencoder. In 2015, a modified sparse tied-weights denoising autoencoder was proposed by Huang et al. . The authors assumed that there exists a mapping function between LR and HR images for both pan and MS or HS. During the training process, an LR pan was generated by the interpolated MS, then the mapping function was learned by the given LR pan patches as the input and HR pan patches as the output. In 2016, a supervised three-layer convolutional neural network was proposed by  to learn the mapping function between the input HR pan with the interpolated LR MS, and the output HR MS. In Qu et al. , a new pansharpening algorithm based on deep learning autoencoder was proposed. Preliminary pansharpening results using both hyperspectral and multispectral images are encouraging.
From the application viewpoint, we categorize the various pansharpening algorithms into three groups. Group 1 methods include coupled non-negative matrix factorization (CNMF) , Bayesian naive , and Bayesian sparse . These methods require the point spread function (PSF) to be available and they perform better than other methods in most cases. Group 2 methods do not require PSF and contain Principal Component Analysis (PCA) , guided filter PCA (GFPCA) , Gram-Schmidt (GS) , GS Adaptive (GSA) , Modulation Transfer Function Generalized Laplacian Pyramid (MTF-GLP) , MTF-GLP with High-Pass Modulation (MTF-GLP-HPM) , Hysure [34, 35], and Smoothing Filter-based Intensity Modulation (SFIM) . Group 3 methods use only the LR HS images and contain the super-resolution (SR) , the bicubic method , and plug-and-play alternating direction multiplier method (PAP-ADMM) .
This chapter is organized as follows. In Section 2, we review the basic idea of color mapping and its variants. In Section 3, we include extensive experiments to illustrate the performance of various pansharpening algorithms. Finally, we conclude our chapter with some future research directions.
2. Proposed hybrid color mapping algorithm and its variants
2.1. Basic idea of color mapping for generating HR HS images
As shown in Figure 1 , the idea of color mapping is to map a color pixel c (i, j) at location (i, j) with R, G, B bands to a hyperspectral pixel X (i, j) at the same location. This mapping is represented by a transformation matrix T, that is
where X (i, j) ∈ RN is a single hyperspectral pixel with N spectral bands, T ∈ R N × M , c (i, j) ∈ RM is a color pixel with M spectral bands, and N>>M. Here, M can be just one band such as the pan band. Hence, color mapping is quite general as it encompasses pan, color, and MS images. Our goal is to generate a HR HS image given a HR color image and an LR HS image. To determine T in Eq. (1), we simulate an LR color image by down-sampling the HR color image. The LR color image and the LR HS image are then used to determine T, which is then used to generate the HR HS image pixel by pixel.
Let us denote H as the set of all hyperspectral pixels X (i, j) for all (i,j) in the image and C as the set of all color pixels c (i, j) for all (i,j) in the image, i =1, …, NR , j = 1, …, NC with NR the number of rows and NC the number of columns in the image. Since X (i, j) and c (i, j) are vectors, H and C can be expressed as
We call the mapping (1) the global version and all pixels in C and H are used in estimating T.
To estimate T, we use the least-square approach, which minimizes the error
To avoid instability, we can add a regularization term in Eq. (3). That is,
And the optimal T becomes 
where λ is a regularization parameter and I is an identity matrix with the same dimension as CCT .
2.2. Hybrid color mapping
For many hyperspectral images, the band wavelengths range from 0.4 to 2.5μm. For color images, the R, G, and B wavelengths are 0.65, 0.51, and 0.475μm, respectively. Because the three color bands may have little correlation with higher-number bands in the hyperspectral image, we found that it is beneficial to extract several higher-number bands from LR HS image and stack them with the LR color bands. This idea is illustrated in Figure 2 . Details can be found in . Moreover, we also noticed that by adding a white band, that is, all pixel values are 1, we can deal with atmospheric and other bias effects.
Using the same treatment earlier, T can be obtained by minimizing the mean square error 
where H is the set of hyperspectral pixels and Ch is the set of hybrid color pixels. All the pixels in H and Ch are used.
The optimal T can be determined as
2.3. Local HCM
We further enhance our method by applying color mapping patch by patch. As shown in Figure 3 , a patch of size p × p is a sub-image in the original image. The patches can be overlapped. This idea allows spatial correlation to be exploited. Our experiments showed that the mapping will become more accurate using this local patch idea. Another advantage of using patches is to split the task into many small jobs so that parallel processing is possible.
2.4. Incorporation of PSF into HCM
Based on observations in  and our own investigations , some pansharpening algorithms with the incorporation of PSF can yield better performance than those without PSF. This motivates to incorporate PSF into our HCM. Our idea can be illustrated by using Figure 4 . The first component is the incorporation of a single-image super-resolution algorithm to enhance the LR hyperspectral image cube. Single-image super-resolution algorithms are well known [45, 46]. The idea is to improve the resolution of an LR image by using internal image statistics. Our proposed method enhances the LR HS bands and then fuses the results using the HCM algorithm. The second component utilizes the HCM algorithm that fuses a high-resolution color image with an enhanced hyperspectral image coming out of the first component. Recently, HCM has been applied to several applications, including enhancing Worldview-3 images , fusion of Landsat and MODIS images , pansharpening of Mastcam images , and fusing of THEMIS and TES .
2.5. Sparsity-based variants using L1 and L0 norms
In , we propose two variants of HCM. From Eq. (5), we can treat the HCM method as an L2 regularization problem. In , we investigate some variants of the regularization. In particular, we would like to apply L1 and L0 norms for the regularization term in Eq. (5). We favor the following two approaches: orthogonal matching pursuit (OMP)  and l1-minimization via augmented Lagrangian multiplier (ALM-l1) , in which OMP is described as an l0-based minimization problem:
where K is the sparsity level of the matrix α (K<<M × N ), and ALM-l1 solves for the l1-minimization convex relaxation problem:
where the positive-weighting parameter λ provides the trade-off between the two terms.
2.6. Application of pansharpening algorithms to debayering
Debayering or demosaicing refers to the reconstruction of missing pixels in the Bayer pattern [57, 58] as shown in Figure 5 . A variant of the Bayer pattern, known as CFA2.0 ( Figure 6 ), was introduced in . Recently, a new approach  to debayering was introduced based on pansharpening ideas. A thorough comparative study was performed using benchmark datasets. It was found that the pansharpening-based debayering approach holds good promise.
Two hyperspectral image datasets were used in our experiments. One was from the Air Force (AF) and the other one was one of the NASA AVIRIS’ images. The AF image ( Figure 7 ) has 124 bands (461–901 nm). The AVIRIS image ( Figure 8 ) has 213 bands (380–2500 nm). Both of them are natural scenes. The AF image size is 267×342×124 and the AVIRIS image size is 300×300×213.
The down-sampled image was used as the low-resolution hyperspectral image that needs to be improved. We picked R, G, and B bands from the original high-resolution hyperspectral image for color mapping. The bicubic method in the following plots was implemented by up-sampling the low-resolution image using bicubic interpolation. The results of bicubic method were used as baseline for comparison study.
3.2. Performance metrics
Similar to , five performance metrics are included here.
Time: It is the computational time in seconds. This metric is machine dependent and varies with runs. However, it gives a relative measure of the computational complexity of different algorithms.
Root mean-squared error (RMSE) : Given two matrices X and , the RMSE is calculated by using
The ideal value of RMSE is 0 if there is perfect reconstruction. To show the performance for each band, we also used RMSE(λ), which is the RMSE value between X(λ) and for each band λ.
Cross-correlation (CC) : It is defined as
where mλ is the number of bands in the hyperspectral image and CCS is the cross-correlation for a single-band image, given by
The ideal value of CC is 1. We also used CC(λ) = , which is the CC value between X(λ) and for each band, to evaluate the performance of different algorithms.
Spectral Angle Mapper (SAM) : It is defined as
where, for two vectors a, b ∈ Rmλ ,
〈a, b〉 is the inner product between two vectors and ‖•‖ denotes the two norms of a vector. The ideal value of SAM is 0.
Erreur relative globale adimensionnelle de synthèse (ERGAS) : It is defined as
where d is the ratio between the linear resolutions of the PAN and HS images, μk is the sample mean of the kth band of X. The ideal value of ERGAS is 0.
3.3. Advantages of the local mapping approach as compared to the global mapping method
We first compared global color mapping (Section 2.1) with local color mapping (Section 2.2) and bicubic interpolation method. Figures 9 and 10 show the RMSE between the ground-truth hyperspectral images and the super-resolution images produced by different methods. We can see that local color mapping is always better than global color mapping. In the AF image, for lower-number bands where R, G, and B bands reside, both global and local color mapping produce better performance than bicubic method. See Figure 9 . However, for higher-number bands, bicubic method is better. The reason is that the spectral correlation between higher-number bands and R, G, B bands is weak. For the AVIRIS image, the local color mapping results shown in Figure 10 are better than both global color mapping and bicubic results across almost all spectral bands.
3.4. Advantages of hybrid color mapping
Figures 11 and 12 show the performance of hybrid color mapping as well as RGBW (R, G, B, and white bands) color mapping described in Section 2. W refers to the white band, that is, an image with all 1s. We can see that adding a white band improves the performance across all bands. Moreover, the hybrid method, which combines more bands from the bicubic interpolated higher-number bands, performed the best in all of the bands. Here, all methods used local patches. The patch size is 4 × 4 and there is no overlapping between the patches.
3.5. Comparison with state-of-the-art algorithms
Now, we include a comparison with a state-of-the-art pan-sharpening algorithm  and another single-image SR algorithm . Both algorithms are recent and compared favorably than their counterparts in their categories. Figures 13 and 14 show the performance of variational wavelet pan-sharpening (VWP) , bicubic interpolation, single-image SR , and our local hybrid color mapping. We used red color band from ground-truth image as the pan image. For NASA data, we observe that in some bands, VWP is better than bicubic and single-image SR methods. However, for bands far away from the reference band, the error is large. In addition, other methods are always worse than our hybrid color mapping method. The reason that VWP  did not perform well in this study is perhaps due to the lack of a pan image which has a spectrum that extends to the high-wavelength regions. If the pan image extends to higher-wavelength regions, it is believed that VWP will perform well. We are grateful to the authors of  for sharing their source codes with us. The reason that the single-image SR  method did not perform well is because it was designed for color images, not for hyperspectral images.
3.6. Pixel clustering enhancement using the SR images
The goal of our research is not only to improve the visual performance of the hyperspectral images by enhancing the spatial resolution but also to enhance the pixel clustering performance as well. Figures 15 – 18 show clustering results using end members extracted from ground-truth AVIRIS hyperspectral image. We used K-means end member extraction technique to determine the clusters. It should be emphasized that we are not performing land-cover classification, which usually requires land-cover spectral signatures as well as atmospheric compensation. The physical meaning of each cluster is not our concern. Comparing Figures 15 and 16 , one can see that hybrid color mapping is significantly better than the bicubic method in terms of the fine details of the images. Moreover, Figure 16 is a lot closer to the ground-truth image in Figure 17 . The images also show that the hybrid color mapping produces much better clustering than the bicubic method as shown in Figure 18 .
3.7. Additional studies and observations
It should be noted that in all of the above studies, we have the following observations:
In the AVIRIS image, our local hybrid color mapping algorithm performed consistently well in all bands. That is, our method is better than bicubic and other methods in all bands.
In the AF image, our local hybrid color mapping algorithm performed also better than others across all bands. However, the performance in lower-number bands is better than that in higher-number bands. This is because the correlation between lower-number bands and higher-number bands is not that good.
We are curious about what will happen if bands other than R, G, B are chosen as reference image in generating the transformation T. Figures 19 and 20 show the comparison. Figure 19 compares the RMSEs using only R, G, B and the case of using bands 20, 50, and 100, for the AF image. Figure 20 compares the RMSEs using only RGB and the case of using bands 20, 90, and 170, for the AVIRIS image. We can observe that if we pick bands from low-, middle-, and high-number bands assuming that they are available, then the overall performance can be much better than that of using R, G, B bands only. The observation shows that we should utilize not only R, G, B bands information but also other bands information such as infrared image and thermal image if they are available.
We investigated the impact of adding regularization to the solution of finding T. In our simulations, we used λ = 10−5 v and v is the maximum singular value of CCT in Eq. (5) or ChCh T in Eq. (8). From Figures 21 and 22 , it can be seen that regularization can also further enhance the performance of our algorithm.
We also investigated the impact of different scaling factors in the down-sampling process. As can be seen in the supplemental materials section, results indicate that bicubic method performs well in low-scaling factors of 1.5 where spatial correlation is high. However, in our applications, we are interested in airborne and satellite images where spatial correlation is weak.
3.8. Comparison of sparsity-based approach with other methods
Here, we focus on a comparison with Groups 1–3 algorithms in the literature. We also compare the performance of different variants of HCM.
3.8.1. Comparison with Group 1 methods
Group 1 methods (BN, BS, CNMF) require PSF. In general, Group 1 methods performed better than Groups 1 and 2. In Section 2.4, we also incorporated PSF into our HCM algorithm. Based on results in Table 1 , we can see that the HCM (L2-norm) yielded better results than those Group 1 methods for the AF data. For the AVIRIS data, BS performed the best.
|Bayes Naive ||0.58||0.4357||0.9881||1.2141||1.6588||0.86||67.2879||0.9474||0.8136||2.1078|
|Bayes Sparse ||208.82||0.4133||0.9900||1.2395||1.5529||235.50||51.7010||0.9619||0.7635||1.8657|
|2||SFIM ||0.99 +||0.7132||0.98489||1.4936||2.2087||1.56y||62.02333||0.94911||0.918||2.0404|
|MTF-GLP ||1.38 +||0.8177||0.98314||1.6095||2.4541||2.25y||55.604||0.95451||0.91031||1.9703|
|MTF-GLP-HTM ||1.40 +||0.8050||0.98356||1.546||2.4214||2.23y||55.641||0.95451||0.90544||1.9718|
|GS ||1.05 +||2.1783||0.85792||2.4421||7.0807||1.83y||54.8851||0.95597||0.93488||1.9524|
|GSA ||1.21 +||0.7435||0.98764||1.5156||2.1764||1.98y||32.2331||0.97016||0.85211||1.6525|
|PCA ||2.37 +||2.3816||0.83824||2.6355||7.7176||2.98y||48.9125||0.96087||0.9173||1.8605|
|GFPCA ||1.17 +||0.6482||0.98615||1.5382||2.0607||2.17y||62.5283||0.93858||1.1736||2.2566|
|Hysure [28, 29]||117.06 +||0.8717||0.98059||1.7882||2.6294||62.47y||38.3131||0.9598||1.0181||1.8586|
|Ours||L2 Norm||3.70 +||0.3889||0.9965||1.1928||1.0986||4.47||32.2094||0.9553||0.9453||1.9101|
|L0 Norm||33.36 +||0.4288||0.9958||1.3677||1.2095||81.11||32.1875||0.9552||0.9362||1.9109|
|L1 Norm||1293.94 +||0.3917||0.9965||1.2259||1.1064||6.22||32.1875||0.9552||0.9362||1.9109|
3.8.2. Comparison with Group 2 methods
Only the GSA method in Group 2 performed better than others for the AVIRIS. Other methods in Group 2 did not perform well as compared to HCM and Group 2 methods. However, Group 2 methods are generally more computationally efficient, which may be advantageous in some real-time applications.
3.8.3. Comparison with Group 3 methods
Group 3 methods do not require PSF or high-resolution color images. Consequently, the performance is generally poor as compared to others. This is understandable.
3.8.4. Comparison among the HCM variants
From Table 1 , we can see that the L2 version of HCM performed better than other variants. This can also be seen from Figures 23 – 26 . However, one key advantage of the L1 and L0 variants is that the models are simpler, as the coefficients of L0 are clustered and the coefficients of L1 are much fewer than that of L2. The above can be confirmed by inspecting Figures 27 and 28 . We would like to mention that another advantage of the sparsity formulation is that it can handle noisy measurements in the hyperspectral images .
In this chapter, we review the various pansharpening algorithms for hyperspectral images. We then focus on one recent algorithm known as hybrid color mapping. Several variants are described. The performance of HCM is thoroughly compared with other methods in the literature.
One future direction is to investigate the performance of different pansharpening algorithms in the presence of noise. Another direction is to apply the pansharpened images to different applications such as target detection [11–13], border monitoring [24, 25], and anomaly detection [14–22, 62].