{'fraction_non_alphanumeric': 0.03457631032094204, 'fraction_numerical': 0.01800969752943893, 'mean_word_length': 4.608449336354807, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hyperspectral-Multispectral Image Fusion with Weighted LASSO

Spectral imaging enables spatially-resolved identification of materials in remote sensing, biomedicine, and astronomy. However, acquisition times require balancing spectral and spatial resolution with signal-to-noise. Hyperspectral imaging provides superior material specificity, while multispectral images are faster to collect at greater fidelity. We propose an approach for fusing hyperspectral and multispectral images to provide high-quality hyperspectral output. The proposed optimization leverages the least absolute shrinkage and selection operator (LASSO) to perform variable selection and regularization. Computational time is reduced by applying the alternating direction method of multipliers (ADMM), as well as initializing the fusion image by estimating it using maximum a posteriori (MAP) based on Hardie's method. We demonstrate that the proposed sparse fusion and reconstruction provides quantitatively superior results when compared to existing methods on publicly available images. Finally, we show how the proposed method can be practically applied in biomedical infrared spectroscopic microscopy.


I. INTRODUCTION
Hyperspectral imaging (HSI) combines spectroscopy with traditional digital imaging to acquire a localized material spectrum at each pixel. Spectral images are acquired using either hyperspectral techniques [1]- [3], which provide high spectral resolution, or multi-spectral techniques [4], which provide a limited subset of spectral bands. Hyperspectral images provide the best molecular specificity, however their high acquisition time imposes limitations in spatial resolution and image quality, making adjacent image features difficult to distinguish even though they may be chemically distinct. Multispectral imaging (MSI) employs tunable sources or filters to provide a subset of spectral samples in a much shorter time. Multi-spectral images generally provide greater spatial resolution at the expense of chemical specificity, particularly when the compounds within the image are unknown a priori. We propose an approach that relies on fusing multi-modal images acquired using both HSI and MSI techniques to generate images with high spatial resolution, molecular specificity, and signal-to-noise (SNR).
Infrared (IR) microscopy is used to identify the distribution of chemical constituents in samples, such as cancer biopsies [5]- [8], without the need for destructive labeling [9]- [11]. However, Fourier transform infrared (FTIR) spectroscopic imaging [12], [13] is impractical for applications that require large samples (eg. 2 cm 2 ) or high spatial resolution (≈1 µm) [14] due to limits in light source intensity, mercury cadmium telluride (MCT) detector resolution, and diffraction-limited optics. While the recent availability of tunable quantum cascade lasers (QCLs) allows discrete-frequency infrared (DFIR) [15], [16], this comes at the cost of spectral resolution, since individual bands must be selected for imaging a priori. In addition, coherent illumination introduces imaging artifacts such as fringing.
Remote sensing applications leverage visible and nearinfrared (NIR) imaging that have similar trade-offs. The Hyperion Imaging Spectrometer on an EO-01 satellite covers a 400 nm to 2500 nm spectral range with a spatial resolution of 30 m [17]. The GeoEye-1 satellite sensor provides a much higher spatial resolution (1.84 m) with a limited spectral range of 450 nm to 920 nm [18].
Several methods have been proposed for multi-modal image fusion [19] including pansharpening, dictionary learning, and convex optimization. Pansharpening fuses an multispectral (MS) tensor with a single high-resolution panchromatic image (PAN) [20]- [22]. This is typically easier than fusing two hyperspectral images, which is computationally prohibitive with existing pansharpening algorithms. Dictionary learning is based on linear spectral unmixing [23], [24] with the following assumptions: (1) multiple scattering between compounds is insignificant, (2) spectral components are discrete and readily separable, and (3) each constituent spectrum is known [25]. Dictionary learning relies on selecting a support tensor containing sufficient sparsity and density to remove noise while accounting for all molecular constituents. The support tensor is learned through Bayesian inference [23], which provides a good prediction if the MSI and HSI have signal-to-noise levels above 30 dB. Previous approaches fuse images by minimizing a convex objective function [26] containing two quadratic data-fitting terms and an edge-preserving regularizer. The data fitting terms account for sampling rate, additive noise, and blur induced by the diffraction limit. The regularizer is a form of vector total variation (TV) which promotes a piece-wise smooth solution with discontinuities aligned across bands. TV regularization has the advantage of preserving edges, but also removes textures and fine-scale detail. TV is a strong prediction/fusion method when the spatial information has limited high-frequency spatial content (sharp edges) [27]. Since high spatial resolution is a priority in many imaging applications, maintaining high-frequency is often critical.
We propose a fusion-based approach that leverages the advantages of HSI and MSI to obtain benefits from both modalities. Our acquisition model integrates modality-specific features of the data, including noise and spectral/spatial sam- pling. Our proposed fusion model relies on least absolute shrinkage and selection operator (LASSO) to regulate and optimize the fusion image. Since LASSO is costly to compute with large hyperspectral images, we employ the alternating direction method of multipliers (ADMM) method for efficiency. Our fusion model is illustrated in Figure 1.
Our contributions are summarized as follows: • Image fusion using proposed LASSO: We use LASSO's Frobenius norm to minimize the differences between fusion and input images (MSI and HSI). LASSO's nuclear norm is used to further reduce noise and enforce sparsity in the result. • Fused image initialization via Maximum a Posteriori (MAP): We initialize our fused image based on previous work by Hardie et al. in [28]. • Implementation of ADMM: We apply ADMM decomposition to increase performance. • Performance Evaluation: We evaluate the proposed method using a variety of fusion metrics and compare results to the state-of-the-art in noise reduction, image fusion, and pansharpening.

II. SYSTEM MODEL
We describe each image as a third-order tensor representing two spatial and one spectral dimension, with spectra sampled in units of wavelength λ or wavenumberν (cm −1 ). The inputs have spatial pixel counts N H = X H × Y H and N M = X M × Y M for hyperspectral and multispectral images, respectively. We expect that N H ≤ N M since multispectral data is most useful with additional spatial detail. The spectral bands are given by Z H and Z M such that Z M < Z H . The following notation is used to express the model: y] T is a position in space and λ is a wavelength • S ∈ R N M ×Z H is the optimal object image at the highest possible spatial and spectral resolution given the pointspread-function p of the imaging system: S = p s The proposed algorithm uses two degraded images H and M of the object s to reconstruct F such that F − S → 0.

A. Image Generation
The ideal image S is degraded by the imaging process. A hyperspectral image with low spatial resolution is generated by: where the operator L ∈ R N H ×N M performs blurring and resampling of S. A multispectral image collected at high resolution is generated by: where B is a sparse matrix that extracts individual bands from S. Both η 1 and η 2 are independent sources of additive noise.

B. Principal Component Analysis
Subspace projection methods like principal component analysis (PCA) convert data with possibly correlated variables into a set of linearly uncorrelated values (i.e. principal components). Only the most energized (i.e. highest variance) components are kept. In HSI, noise is typically distributed isotropically throughout the spectrum, while actual signals lie within a smaller subspace. Furthermore, HS images are often large, with a 1 mm 2 biopsy image (2000 × 2000 × 371 pixels) requiring ≈5.8 GB at 32-bit precision. Dimension reduction reduces computational costs and memory requirements while increasing SNR by eliminating low-variance coefficients associated with noise.
Our goal is to obtain a fusion image .., f i,Z H ] represents the spectrum at each spatial pixel location. We represent f i in the subspace Q as: where r i ∈ R 1×Z is the projection of f i onto the orthonormal subspace spanned by the columns of Q ∈ RZ ×Z H such thatZ Z H . The image R ∈ R N M ×Z projected onto this subspace is composed of the reduced spectra: R = [r 1 , ..., r N M ] T . Integrating the dimension reduction transformation Q into (1) and (2) yields:

III. PROPOSED FUSION ALGORITHM
From equation (4) and (5), we propose using the Least Absolute Shrinkage and Selection Operator (LASSO) to optimize R. The fusion image is initialized with Maximum a Posteriori (MAP) [28] and the optimization is broken down into subproblems using the Alternating Direction Method of Multipliers (ADMM). The LASSO Frobenius norm minimizes the differences between the fusion image and accessible components of the observed images, while the nuclear norm further reduces noise and enforces sparsity as needed.

A. Regularization
A general LASSO optimization proposed by Tibshirani et al. [29] is given by: where is the observed data • η is a parameter that controls the degree of shrinkage We optimize R based on its relationships with the available HSI and MSI in equations (4) and (5) using the Frobenius norm: The HSI has lower resolution and SNR, therefore we implement weights for each term in Equation 7 to prioritize appropriate information from M. These weights are based on SNR, resolution, and blurriness. The MSI and HSI weights are on a logarithmic scale ranging from 30 to 50 and 5 to 30, respectively. After applying the weights, the equation (7) becomes: Equation (8) generates a fusion image that maximizes the integration of information from both H and M. The resulting fusion image will also contain a portion of the noise from the two source images. The variable selection component of LASSO, represented by the nuclear norm, is employed to remove these noise terms: The nuclear norm term minimizes the sum of magnitudes of R (Section III-C). Pixels with absolute magnitudes smaller than a threshold T , defined in Section III-C as η µ , are set to zero.
Since P (H, M) is not a function ofR, the optimization for P (R|H, M) (Equation 10) reduces to: Continuing to follow Hardie's work, we obtain the initialization ofR as: • ΛR |M is the row covariance matrix of R given M as: Solving (12) directly gives usR:

C. Alternating Direction Method of Multipliers (ADMM)
Since hyperspectral images tend to be large, with our application data in the range of (1128 × 1152 × 371), ADMM [30] is used to partition the optimization. ADMM is widely popular due to its ability to handle large data effectively [31]- [33] by breaking an expensive optimizations into less costly sub-problems. To solve the LASSO equation (9), we introduce the splitting variables W 1 = LR and scaled Lagrange multipliers J 1 for the first Frobenius norm term, W 2 = R and J 2 for the second Frobenius norm term, W 3 = R and J 3 for the Nuclear norm term. The augmented Lagrangian associated with the optimization of R is: where µ and η are the scalar regularization parameters, with η ≈ 1.25·10 −3 · H ∞ and µ ≈ 5·10 −2 · Λ H Z H . Minimizing (14) with respect to one term while holding others fixed allows the following update functions: For W 3 , we implemented soft-thresholding to solve it fast and efficiently, thus we first define the following term for the Frobenius term: Then, the optimization for W 3 becomes: which has an optimality condition at: 3i,j = 0, we have: Equation (20) for those pixels is re-written as: In equation (22), we see that when w u+1 3i,j ) LS . We then rewrite (22) as: If the pixel at location [i, j] has w u+1 3i,j = 0, the sub-differential of the nuclear norm is the interval [−1, 1]. We rewrite equation (20) for those pixels as: So we update W u+1 3 element-wise as: The final algorithm is given in Algorithm 1 with ADMM variables.

D. Complexity Analysis
Each iteration of the first for loop has a computational com- The following terms are independent of iteration and precomputed:

A. Data Collection and Simulation
In this section, we compare our algorithm to Bayes Sparse [23], convex optimization (HySure) [26], coupled sparse tensor factorization (CSTF) [34], and convex optimizationbased coupled nonnegative matrix factorization (CO-CNMF) [35]. We first describe our image acquisition and simulation methods, including the infrared microscopy approach that motivates this work (Section IV-A). We then describe the fusion metrics used to evaluate the resulting images (Section IV-B). Finally, we compare our proposed fusion framework on a variety of publicly available data sets and summarize the results (Section IV-C).
Arial images of the University of Pavia and Pavia Center were acquired from the Deutsches Zentrum fur Luft-und Raumfahrt under sponsorship from European Union's Hysens project. The raw images have a spatial resolution of 1.3 m with 103 spectral bands. HSI data was simulated by applying a 39×39 pixel Gaussian filter and downsampling every 4 pixels to reduce the spatial resolution to ≈5.2 m. We determine the Gaussian filter size by applying it on MSI from DFIR imaging system demonstrated in Figure 2.a. We increase the filter size until MSI's spatial information in Collagen and Epithelium areas cannot be distinguished similar to HSI taken by FTIR imaging system in Figure 2.b. MSI data is simulated by randomly selecting from the first 70 bands. Gaussian noise is added at 5dB to 10dB SNR (HSI) and at 50dB SNR (MSI). In University of Pavia and Pavia Center datasets, the MSI is 512 × 256 × 4 pixels and the HSI is 128 × 64 × 93 pixels. The University of Pavia dataset was tested with Gaussian noise, then Pavia Center dataset with Gaussian (Gau.) noise and Poisson (Poi.) noise.
Infrared spectroscopic images were obtained using two commercially available instruments: (1) HSI images were collected using an Agilent Cary 620 infrared spectroscopic microscope with 5.5 µm spatial resolution and a wavelength range of 1000 to 4000cm −1 , and (2) MSI images were acquired using a Daylight SPERO microscope with a spatial resolution of 1.3 µm and discrete sampling at wavelengths within a range of 900 to 1700cm −1 . Images of an identical breast biopsy core were collected using both systems. The HSI data was collected at high-speed using a single pass (without averaging). While traditionally designed for discreteband imaging, the SPERO microscope was used to collect a complete spectrum from 1000 nm to 1700 nm to provide a ground truth for comparison. Individual bands were drawn from this full data cube to represent MSI input to our algorithm. In this data set, MSI has a resolution of 1128×1152×10 pixels and HSI has a resolution of 282 × 288 × 371 pixels. We did not calculate for SNR in HSI and MSI, nor do we calculate for the blur kernel in them. We gave them a guessed estimate value for HSI at 30dB SNR and MSI at 50db SNR. For blur kernel, we apply a Gaussian blur kernel on MSI image until they resemble HSI image. The size of the Gaussian blur kernel is 39x39.
The Indian Pines (IP) dataset is obtained by AVIRIS sensor over the Indian Pines test site in North-western Indiana, sponsored by NASA's Jet Propulsion Laboratory. The dataset has a spatial resolution of 20 m, and 224 spectral bands. Salinas (Sal.) datasets is also obtained by AVIRIS sensor over Salinas Valley, California. It has a high spatial resolution at 3.7 m, and 224 spectral bands. Last dataset is Kennedy Space Center (KSC), Florida which is also taken by AVIRIS sensor. It has a spatial resolution of 18 m and 224 spectral bands. These three datasets are applied a similar blur filter and downsampling factor to Pavia University and Pavia Center. We then added them with Gaussian and Poisson noise in each experiment at 10 dB SNR. In Indian Pines dataset, the MSI has a resolution of 128 × 128 × 10, and HSI has a resolution of 32 × 32 × 200. In Salinas dataset, the MSI has a resolution of 256 × 128 × 10, and HSI has a resolution of 64 × 32 × 204. In Kennedy Space Center dataset, the MSI has a resolution of 256 × 128 × 10, and HSI has a resolution of 64 × 32 × 128.
All algorithms are adjusted to be fed with the same blurred and noise added data set in simulation case. In breast tissue case, SNR level of HSI is divided into two regions, one that overlapped with MSI's spectra (guessed SNR set at 30dB) and one that extends beyond MSI's spectra (guessed SNR set at 50dB). MSI has guessed SNR at 50dB. All algorithms are again provided with the same guessed SNR level of HSI and MSI. HS images are assumed to be blurry, so we also added an estimated blurred kernel to all algorithms.

B. Fusion Metrics
Peak signal to noise ratio (PSNR) is the ratio of peak signal power to noise power expressed in term of logarithmic a decibel scale: where S is the reference image andF is the fused image. Root-mean-square error (RMSE) is used to measure the differences between the predicted model F to a ground truth S. The smaller the value of RMSE, the better the fusion image quality: where N M is the number of pixels in MS image, Z H is the number of spectral bands in HS image. The spectral angle mapper (SAM), proposed by Alparone [36], quantifies similarities by calculating the spectral difference between fused image and the ground truth at each spatial pixel and then averaging them. SAM of the fusion image and ground truth vector The final SAM value is the average across all pixels and ranges from (-90, 90] degrees. An optimal SAM value is at 0 when the fusion image is exactly the same as ground truth. The universal image quality index (UIQI) was proposed by Wang and Bovik [37] to focus on the luminance, contrast, and struture of the fusion image. UIQI is calculated on each band where σ s,f is the covariance of (s, f ) while others values are means and variances. UIQI lies in the range [-1,1], where 1 implies that f = s. The comprehensive UIQI for an image is the average across all bands.
The relative dimensionless global error in synthesis (ER-GAS) , proposed by Du [38], measures changes in mean and dynamic range, while also providing a comprehensive estimate of image quality. ERGAS calculates the amount of spectral distortion in the image using: where µ i is the mean of the fused image at band i. A small ERGAS indicates less geometric or radiometric distortion, when ERGAS is 0 it means F = S. The degree of distortion (DD), introduced by Zhu and Bamler [39], which calculates the differences between two images S, F, is defined as: Note that the tensor S and F are vectorized in this expression. Smaller DD indicates a higher quality fusion image.

C. Result Comparison
In this section, we compare the differences in algorithms and fusion metrics between our method, Bayes Sparse and HySure due to their performance results are closer to ours than CSTF and CO-CNMF. Our algorithm outperforms these methods by a significant margin, with performance converging after 10dB (Table I). Similarly, our metrics values are much closer to the ideal value (zero) using RMSE, DD, and ERGAS metrics. High performance in the presence of low SNR is critical for expanding FTIR to medical imaging by providing data for a complete biopsy in a couple of hours. While image quality is an important metric, we also evaluate the tensor's pixel-level spectral response using the SAM metric, where our algorithm is able to produce an average of 2.4 for all pixels. The second best algorithm is Bayes Sparse with SAM average value at 4.8. In Table I, our algorithms has a significantly reduced computational time compared to others at an average of 5.16 s in University of Pavia dataset. The second best time belongs to HySure at an average of 38.32 s. Algorithm time is measured after HSI and MSI is simulated and fed to algorithms until a fused image is produced. In the case of Bayes Sparse, the Algorithm Time is the slowest of all methods. It takes about 105 s to train dictionaries, and about 70 s to fuse images. Figure 3 shows the difference between our results and fusion images produced using competing methods (Bayes Sparse, HySure, CO-CNMF, and CSTF) for a 10dB SNR and 39 × 39 Gaussian blur kernel. The two fusion images which are closest to Ground Truth are ours and Bayes Sparse.
HySure fusion image has a small amount of artifacts, while CO-CNMF and CSTF fusion image have a significant amount of artifacts. When we ran all algorithms with the same SNR and a 5 × 5 blur kernel, the artifacts are not present in CO-CNMF, CSTF, and HySure anymore. We concluded that the three algorithms (HySure, CO-CNMF, and CSTF) are not not capable of with resolving HSI's blur effect which HSI suffers seriously from. When we look at the color fidelity in LASSO and Bayes Sparse fusion images, both of them have accurate color reproduction except the area of landmass in the middle of the image, highlighted by Figure 3's Best Cases, where ours has slight reddish color and Bayes Sparse has a reddish color instead of bright yellow in Ground Truth.
We suspect HySure's low performance is caused by its inability to adjust weights based on to MSI and HSI noise. We present equation (6) in [26] here: The sparse regularizer is a vector Total Variation represented by the function ϕ(). Compared to our minimization functions in Equations 9 and 14, HySure uses similar scalar multipliers (such as µ and η) and a similar optimization scheme. These scalar multipliers act as weights on how important the regularization on those terms as well as defining the threshold for nuclear norm. However, HySure lacks penalty weights inside the Frobenius norms for HS and MS terms, which means the data from HSI has the same contribution as data from MSI The spectral resolution is 371 bands spanning from 1000cm −1 to 3700cm −1 , except DFIR which has 220 bands spanning from 1000cm −1 to 1760cm −1 . In this graph, the solid line is the mean of an area in the epithelium region of the tissue. The shaded area is the standard deviation of the pixels in the area. The blue graph depicts the DFIR image's spectral response, the grey graph depicts the FTIR image's spectral response, the orange graph depict dictionary learning image's spectral response, the green graph depict LASSO image's spectral response.
leading to a worse fused image. By utilizing these logarithmic scale weights of HSI and MSI, we dictate which image has higher contribution to our fusion image at a spectral range of our choice. In dictionary learning method, although we have similar initialization scheme for R and regularization for Frobenius norm, the sparse regularizer in Bayes Sparse is learned from the Bayesian MAPR. A support will be built based on this and it will identify nonzero elements on each row of R, meaning finding nonzero values of a pixel in spectral dimension. Thus, only nonzero elements defined by dictionaries learnt fromR will be updated. Bayes Sparse method works well if the correct support is found, but if the support is incorrect it will filter out both signal and noise or leaving noise behind. Instead of assumingR has the correct support, we adopt LASSO's nuclear norm as our sparse regularizer/noise reduction mechanism. In our algorithm, LASSO's Frobenius norms will grasp all the information available from HSI and MSI, then nuclear norm will reduce noise or enforce sparsity. By doing so, we lose less information compared to Bayes Sparse and giving us a globally optimal solution.
The spectra fidelity in fused images is best demonstrated in Figure 4. In this figure, SAM metric values at each pixels are graphed. The range of SAM value here is from 0 • to 57 • , with 0 • means no spectral distortion between fused image and ground truth at the particular pixel, while 57 • means the angle between the fused image and ground truth at that pixel's spectra is 57 • out of phase. In short, in this figure the darker the image, the less spectral distortion in fused image. Both LASSO and Bayes Sparse fused image is outperformed by HySure at the building with painted metal sheets (area with a highlighted SAM value), where HySure's SAM is at 12.07 • ,. However, in the rest of the image, our fusion methods outperforms the other two in every areas: asphalt, meadows, gravel, trees, baresoil, bitumen, self-locking bricks, shadows. Especially in meadows, bitumen, asphalt, self-locking bricks, and trees areas, where our fused image has a mean of SAM at 1 • and a standard deviation of 0.13 • .
In Figure 5, we plot the performance in multiple metrics of the fusion methods at each individual band. The behavior of Bayes Sparse is very interesting when looking at these graphs. At earlier bands between 1 and 68, Bayes Sparse performs very well and on par with our method. Bayes Sparse and our methods produce significantly better result The spectral resolution is 371 bands spanning from 1000cm −1 to 3700cm −1 , except DFIR which has 220 bands spanning from 1000cm −1 to 1760cm −1 . In this graph, the solid line is the mean of an area in the collagen region of the tissue. The shaded area is the standard deviation of the pixels in the area. The blue graph depicts the DFIR image's spectral response, the grey graph depicts the FTIR image's spectral response, the orange graph depict dictionary learning image's spectral response, the green graph depict LASSO image's spectral response. than HySure, with PSNR level difference ranging from 3dB to 10dB and occasional equalization. In ERGAS, DD and RMSE metric, there is a much wider difference between those two methods and HySure. In UIQI metrics, the differences between HySure, Bayes Sparse, and our proposed approach are reduced. However, Bayes Sparse's performance in all metrics after band 68th are unsteady. This behavior persists in all of the data sets contaminated with Gaussian noise we tested here. We suspect this is due to the dictionaries learned from MS images at bands which lies at the beginning of the spectra. Hence, the data it represents no longer valid at the end of the spectra.
From here on, HySure, CO-CNMF, and CSTF algorithms will not be included because their performance are lacking compared to the our proposed method and Bayes Sparse. In Figure 6, we are showing a set of images from a breast tissue at two different bands. On the first row, images from DFIR imaging system(MSI), FTIR imaging system (HSI), our method (LASSO), Wei's method (Bayes Sparse) are at wavelength 1550 cm −1 . Inside the tissue core, there are regions of epithelium and collagen. In this tissue data set, the DFIR's image, Figure 6.a, which shows a good contrast between collagen and epithelium, as well as showing patterns in the collagen and epithelium regions. DFIR's image, MSI, has a spatial resolution of 1128×1152 pixels. The HSI, Figure  6.b, has a very low spatial resolution at 282 × 288 pixels. In HSI, the distinction between collagen and epithelium is clear, but all spatial patterns within the two regions are indistinguishable. The next image in row one is our fusion method's result, Figure 6.c. The contrast in our fusion image is closer to the original image than Bayes Sparse's fusion image. Our fused image has similar patterns to DFIR image in both collagen and epithelium regions. The fusion image from Bayes Sparse method in [23] is on the upper right corner, Figure  6.d. In this image, the contrast is not close to the MSI. Both fusion images have fringing effect, but ours is less than Bayes Sparse's. On the second row, Bayes Sparse fusion image has their contrast distorted compared to HSI. Fringing effect in Bayes Sparse image is significantly worse here while ours is lessen. We suspect again this is due to the dictionaries learned by Bayes Sparse method are from earlier bands in the spectra, which is heavily based on DFIR image which has high spatial resolution with fringing effect.
In Figure 7 and 8, we compare the mean and standard deviation of spectra from the original and fusion images using our algorithm and Bayes Sparse. Reconstructed spectra represent the mean values of two tissue classes: epithelium and collagen. The performance of our algorithm is similar to Bayes Sparse in the fingerprint region (1000 to 1700cm −1 ) where data is available for both FTIR and DFIR source images (HSI and MSI). However, the statistical distribution of spectra in at larger wavenumbers conforms much better to the source FTIR image using the proposed algorithm. Bayes Sparse's statistical deviation corresponds to the error found in the Pavia data set and also correlates with large SAM values.
Due to the fact that we do not have a ground truth with high resolution in both spatial and spectral domain for this data set, we have to compare the fused image to MSI in 1000 cm −1 to 1700 cm −1 region, and HSI in the entire spectrum. Table  II and III help us showcase our performance numerically. In Table II, the matched wavelengths from 1000 cm −1 to 1700 cm −1 in HSI, LASSO, and Bayes Sparse images are extracted and compared. In this region, both algorithms perform very well with ours leads on all metrics but UIQI. We leads in the following metrics: PSNR (3.8 dB difference), RMSE (5 · 10 −3 difference), DD (2.5 · 10 −3 difference). This confirms the behaviors of our graphs in Figure 7, where both fused images' spectra follows DFIR's closely. In Table III, the fused images are down-sampled by averaging and compared to FTIR image on the entire spectral range. Our algorithm leads by a wide margin here. We leads in the following metrics: PSNR (4.4 dB difference), RMSE (4.5 · 10 −2 difference), UIQI (0.83 · 10 −3 difference), DD (2.4 · 10 −3 difference). The important metrics in this comparison is Degree of Distortion, because this metric's purpose is to evaluate fused image with no ground truth. In both Tables II and III, the metrics SAM and ERGAS are not included. The tissue image has background regions (dark blue areas in Figure 6) in which pixels values are zero which lead to SAM values at those pixel cannot be calculated. Also, in the spectral range of 1800 cm −1 to 2800 cm −1 , there are no spectral information so mean values at these wavelengths are close to zero. Therefore, ERGAS values in this region approaches infinity.
The last thing we want to discuss in this data set's experiment is computational time. The fusion time measures the time to solve the fusion problem only, while algorithm time measures time in all calculations steps after data (HSI and MSI) is fed into algorithms. In Fusion Time, our method requires 354.2 s to complete and Bayes Sparse requires 2170.9 s to complete. In Algorithm (Alg.) time, our method requires 479.3 s to complete and Bayes Sparse requires 6976.9 s to complete.
In Figure 9 (Supplemental Materials), we represent Pavia Center data set's fused images when the HSI and MSI are contaminated with Poisson noise at 10 dB SNR and 50 dB SNR, respectively. Noises were removed from both fused images, but the color fidelity in them is not as good as in our previous Pavia University experiment. In Figure 10 (Supplemental Materials), the different in spectral quality at each pixel is shown. In all materials, our fused image has much darker color than Bayes Sparse's fused image with an overall average SAM value at 3.25 • , while Bayes Sparse's fuse image has an overall average SAM value at 6.46 • . In Table IV (Supplemental Materials), we show our method and Bayes Sprase's performance in different metric for Pavia Center data set with Poisson noise. At all HSI noise level, our proposed method outperforms Bayes Sparse similarly to Table  I, but there is a reduction in fused images' quality in both methods when SNR decreases in HSI. The noise distribution type is the cause of the drop in fusion images' quality. In order to prove this, we ran the experiment with Pavia Center again with Gaussian noise. In Table V (Supplemental Materials), our method's performance is again consistent at all SNR level in HSI, which is similar to our method's performance in Pavia University data set with Gaussian noise. In PSNR metric, when switching from Gaussian noise to Poisson noise, our method and Bayes Sparse lose an average of 6.66 dB and 1.86 dB, respectively. In Figure 11 (Supplemental Materials), we once again look at the two methods' performance at each band with Pavia Center data set. The noise distribution here is Poisson and it provides us an interesting change in behavior of Bayes Sparse's performance. Bayes Sparse's graphs fluctuates heavily through the entire spectral bands, not just in the later bands. In Figure 12 (Supplemental Materials), we show the performance of our method and Bayes Sparse in Pavia Center data set with Gaussian noise. Bayes Sparse graphs fluctuates at the later bands again similar to our experiment with Pavia University data set in Figure 5.
In Table VI (Supplemental Materials), we performed similar experiments on Indian Pines (IP), Kennedy Space Center (KSC), and Salinas (Sal.) data sets in order to prove our method's superiority. Both Bayes Sparse and our method are fed with the same sets of HSIs and MSIs under two noise distributions. In Indian Pines and Salinas data sets, our method outperforms again with an average of 4 dB difference in PSNR metric. In KSC data set, our sparsest data set, the performance differences between the two fusion methods is negligible. Noise environments still have the same effect on both fusion methods. With Gaussian noise we have an average of 6.2 dB difference in PSNR between our method and Bayes Sparse. With Poisson noise we have an average of 1.85 dB difference in PSNR between our method and Bayes Sparse.

V. CONCLUSION
In this paper, we proposed an approach in image fusion with LASSO. With the implementation of ADMM and Bayesian MAP for initialization, our algorithm obtains a significantly better result when our inputs are badly blurred and heavily damaged with noise. Our method outperforms other fusion methods in every performance metric as well as computing time, which is largely due to our robustness of sparse regularization.