Deep unfolding Network for Hyperspectral Image Super-Resolution with Automatic Exposure Correction

In this paper, a scalar, a vector, a matrix, and a tensor are denoted by lowercase $x$, boldface lowercase $\mathbf{x}$, boldface capital $\mathbf{X}$, and calligraphic $\mathcal{X}$, respectively. $\mathbf{C}=\mathbf{A}\circ\mathbf{B}$ denotes the matrix-formed element-wise multiplication of $\mathbf{A}$ and $\mathbf{B}$, where matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ have the same dimension. $\|\mathbf{A}\|_{\text{F}}$ and $\|\mathbf{A}\|_{1}$ denote Frobenius norm and $\ell_{1}$ norm of matrix $\mathbf{A}$, respectively. For a tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, the mode-$i$ unfolding martix is defined as $\mathbf{X}_{(i)}\in\mathbb{R}^{I_{i}\times\prod_{j\neq i}I_{j}}$.

HSI-SR aims to recover an HR-HSI $\mathcal{Z}\in\mathbb{R}^{C\times W\times H}$ from an LR-HSI $\mathcal{X}\in\mathbb{R}^{C\times W_{\text{HSI}}\times H_{\text{HSI}}}$ and an HR-MSI $\mathcal{Y}\in\mathbb{R}^{C_{\text{MSI}}\times W\times H}$, where $C$ and $C_{\text{MSI}}$ are the spectral numbers in HSI and MSI, respectively. The degradation model for LR-HSI and HR-MSI can be expressed as follows:

		$\displaystyle\mathbf{X}_{(1)}=\mathbf{Z}_{(1)}\mathbf{H}$		(1)
		$\displaystyle\mathbf{Y}_{(1)}=\mathbf{P}\mathbf{Z}_{(1)}$

where $\mathbf{X}_{(1)}\in\mathbb{R}^{C\times N_{\text{HSI}}}$, $\mathbf{Y}_{(1)}\in\mathbb{R}^{C_{\text{MSI}}\times N}$, $\mathbf{Z}_{(1)}\in\mathbb{R}^{C\times N}$ are obtained by a mode-1 unfolding of $\mathcal{X}$, $\mathcal{Y}$, $\mathcal{Z}$, respectively, $N=W\times H$, $N_{\text{HSI}}=W_{\text{HSI}}\times H_{\text{HSI}}$ represent the number of pixels in $\mathbf{X}_{(1)}$ and $\mathbf{Y}_{(1)}$, respectively, $\mathbf{H}\in\mathbb{R}^{N\times N_{\text{HSI}}}$ is a spatial degeneracy matrix, which accounts for blurring and spatial downsampling operations, and $\mathbf{P}\in\mathbb{R}^{C_{\text{MSI}}\times C}$ is the spectral response matrix associated with the imaging sensor. For notational simplicity, in the following we will denote $\mathbf{X}_{(1)},\mathbf{Y}_{(1)},\mathbf{Z}_{(1)}$ by $\mathbf{X},\mathbf{Y},\mathbf{Z}$.

Considering the different exposure levels of the observational data $\mathcal{X}$ and $\mathcal{Y}$, we design a new HSI fusion model, which can be expressed as:

	$\displaystyle\mathop{\min}\limits_{\mathbf{L}_{1},\mathbf{L}_{2},\mathbf{Z}}$	$\displaystyle\frac{1}{2}\|\mathbf{X}-(\mathbf{Z}\circ\mathbf{L}_{1})\mathbf{H}\|_{\text{F}}^{2}+\frac{1}{2}\|\mathbf{Y}-\mathbf{P}(\mathbf{Z}\circ\mathbf{L}_{2})\|_{\text{F}}^{2}$		(2)
		$\displaystyle+\beta_{1}\Phi_{1}(\mathbf{L}_{1})+\beta_{2}\Phi_{2}(\mathbf{L}_{2})+\beta_{3}\Phi_{3}(\mathbf{Z})$

where $\mathbf{L}_{1},\mathbf{L}_{2}\in\mathbb{R}^{C\times N}$ represent the different levels of exposure on $\mathbf{Z}$, respectively, and $\Phi_{1}(\cdot)$, $\Phi_{2}(\cdot)$, $\Phi_{3}(\cdot)$ represent implicit regularizers for $\mathbf{L}_{1}$, $\mathbf{L}_{2}$, $\mathbf{Z}$, respectively. Note that $\mathbf{Z}\circ\mathbf{L}_{1}$ and $\mathbf{Z}\circ\mathbf{L}_{2}$ represent the light degraded $\mathbf{Z}$ at two different exposure settings, respectively.

An iterative method for solving problem (2) is proposed by solving three subproblems as follows:

		$\displaystyle\mathbf{L}_{1}=\mathop{\arg\min}\limits_{\mathbf{L}_{1}}\frac{1}{2}\|\mathbf{X}-(\mathbf{Z}\circ\mathbf{L}_{1})\mathbf{H}\|_{\text{F}}^{2}+\beta_{1}\Phi_{1}(\mathbf{L}_{1})$		(3)
		$\displaystyle\mathbf{L}_{2}=\mathop{\arg\min}\limits_{\mathbf{L}_{2}}\frac{1}{2}\|\mathbf{Y}-\mathbf{P}(\mathbf{Z}\circ\mathbf{L}_{2})\|_{\text{F}}^{2}+\beta_{2}\Phi_{2}(\mathbf{L}_{2})$
		$\displaystyle\begin{aligned} \mathbf{Z}=\mathop{\arg\min}\limits_{\mathbf{Z}}&\frac{1}{2}\|\mathbf{X}-(\mathbf{Z}\circ\mathbf{L}_{1})\mathbf{H}\|_{\text{F}}^{2}\\ +&\frac{1}{2}\|\mathbf{Y}-\mathbf{P}(\mathbf{Z}\circ\mathbf{L}_{2})\|_{\text{F}}^{2}+\beta_{3}\Phi_{3}(\mathbf{Z})\end{aligned}$

At the iteration $t+1$, applying the PGD algorithm to the three subproblems in (3) yields:

		$\displaystyle\mathbf{L}_{1}^{t+1}=\text{Prox}_{\beta_{1}\Phi_{1}}(\mathbf{L}_{1}^{t}+\beta_{1}\mathbf{Z}^{t}\circ((\mathbf{X}-(\mathbf{Z}^{t}\circ\mathbf{L}_{1}^{t})\mathbf{H})\mathbf{H}^{T}))$		(4)
		$\displaystyle\mathbf{L}_{2}^{t+1}=\text{Prox}_{\beta_{2}\Phi_{2}}(\mathbf{L}_{2}^{t}+\beta_{2}\mathbf{Z}^{t}\circ(\mathbf{P}^{T}(\mathbf{Y}-\mathbf{P}(\mathbf{Z}^{t}\circ\mathbf{L}_{2}^{t}))))$
		$\displaystyle\begin{aligned} \mathbf{Z}^{t+1}=\text{Prox}_{\beta_{3}\Phi_{3}}(\mathbf{Z}^{t}&+\beta_{3}((\mathbf{X}-(\mathbf{Z}^{t}\circ\mathbf{L}_{1}^{t+1})\mathbf{H})\mathbf{H}^{T}\circ\mathbf{L}_{1}^{t+1}\\ +&\mathbf{P}^{T}(\mathbf{Y}-\mathbf{P}(\mathbf{Z}^{t}\circ\mathbf{L}_{2}^{t+1}))\circ\mathbf{L}_{2}^{t+1}))\end{aligned}$

where $t=1,2,\cdots,T-1$, and $T$ is the maximum number of iterations, and $\text{Prox}_{\beta\Phi}(\cdot)$ denotes the proximal operator defined as

$$\text{Prox}_{\beta\Phi}(\mathbf{A})=\mathop{\arg\min}\limits_{\mathbf{M}}\frac{1}{2}\|\mathbf{A}-\mathbf{M}\|_{F}^{2}+\beta\Phi(\mathbf{M})$$

(5)

The iterative updating stage in (4) can be implemented by the proposed UHSR-AEC learning network as shown in Figure 2(a), through $T$ identical unfolding stages for regularizers’ learning after the initialization. Next, let us present the structure of the proposed UHSR-AEC in more detail, followed by the initialization module.

Network Structure for Implicit Regularizer: Figure 2 (b) shows the designed residual recovery network (RRN) for implicit regularizers. Three RRNs with identical structure for three respective regularizers $\Phi_{1}(\cdot),\Phi_{2}(\cdot),\Phi_{3}(\cdot)$ are used at each learning stage. Each RRN consists of two identical ResConvBlocks, where the residual structure is considered to avoid the gradient vanishing and the performance improvement.

Network Structure for Sampling Matrix: To overcome the mismatch between the spatial and spectral responses used in the training and test data, we design learning sampling matrices to replace $\mathbf{H},\mathbf{H}^{T},\mathbf{P},\mathbf{P}^{T}$, where $\mathbf{P}$ and $\mathbf{P}^{T}$ ($\mathbf{H}$ and $\mathbf{H}^{T}$) each contains one Conv (three ConvBlocks).

Randomly initialized $\mathbf{L}_{1}^{0}$, $\mathbf{L}_{2}^{0}$, $\mathbf{Z}^{0}$ may destroy the texture and details of the image, resulting in degradation of the recovery quality. Therefore, the initialization module (IM) is designed for better initial $\mathbf{L}_{1}^{0}$, $\mathbf{L}_{2}^{0}$, $\mathbf{Z}^{0}$.

Figure 2(c) shows the network structure of IM, composed of a feature extraction (FE) layer, a cross-attention (CA) layer, and a feature fusion (FF) layer. The FE first projects $\mathbf{X}$ and $\mathbf{Y}$ onto the high-dimensional space, and upsamples $\mathbf{X}$ and $\mathbf{Y}$ in the spatial domain and the spectral domain, respectively; the CA extracts the Channel Attention of $\mathbf{X}$ and the Spatial Attention of $\mathbf{Y}$ for high-quality spectral and spatial information; the FF layer formed by a Unet, fuses features of $\mathbf{X}$ and $\mathbf{Y}$ and maps them to $\mathbf{L}_{1}^{0}$, $\mathbf{L}_{2}^{0}$, $\mathbf{Z}^{0}$, by solving the following problem:

	$\displaystyle\mathop{\min}\limits_{\mathbf{L}_{1}^{0},\mathbf{L}_{2}^{0},\mathbf{Z}^{0}}$	$\displaystyle\|\hat{\mathbf{Z}}-\mathbf{Z}^{0}\|_{1}+\lambda\|\hat{\mathbf{Z}}_{x}-(\mathbf{Z}^{0}\circ\mathbf{L}_{1}^{0})\|_{1}$		(6)
	$\displaystyle+$	$\displaystyle\lambda\|\hat{\mathbf{Z}}_{y}-(\mathbf{Z}^{0}\circ\mathbf{L}_{2}^{0})\|_{1}$

where $\hat{\mathbf{Z}}$, $\hat{\mathbf{Z}}_{x}$ and $\hat{\mathbf{Z}}_{y}$ represent the ground truth, the associated exposure degraded HSI and MSI to mitigate the degradation of exposure, respectively. Note that IM is important in the training of UHSR-AEC, but it is virtual as the trained UHSR-AEC operates for the reconstruction of HSI-SR.

The proposed UHSR-AEC is trained through an end-to-end approach, while IM is fine-tuned at first during the training, together with the weights shared among $\Phi_{1}(\cdot)$, $\Phi_{2}(\cdot)$, and $\Phi_{3}(\cdot)$ at different sequential unfolding stages. The optimization problem with a $\ell_{1}$-norm based loss function, considered for training UHSR-AEC by its greater robustness against outliers than the Frobenius norm, is as follows:

	$\displaystyle\mathop{\min}\limits_{\mathbf{L}_{1}^{T},\mathbf{L}_{2}^{T},\mathbf{Z}^{T}}$	$\displaystyle\|\hat{\mathbf{Z}}-\mathbf{Z}^{T}\|_{1}+\eta_{1}\|\hat{\mathbf{Z}}_{x}-(\mathbf{Z}^{T}\circ\mathbf{L}_{1}^{T})\|_{1}$		(7)
	$\displaystyle+$	$\displaystyle\eta_{1}\|\hat{\mathbf{Z}}_{y}-(\mathbf{Z}^{T}\circ\mathbf{L}_{2}^{T})\|_{1}+\eta_{2}\|\hat{\mathbf{X}}-(\mathbf{Z}^{T}\circ\mathbf{L}_{1}^{T})\mathbf{H}\|_{1}$
	$\displaystyle+$	$\displaystyle\eta_{2}\|\hat{\mathbf{Y}}-\mathbf{P}(\mathbf{Z}^{T}\circ\mathbf{L}_{2}^{T})\|_{1}$

Datasets: The datasets CAVE¹¹1https://www1.cs.columbia.edu/CAVE/databases/multispectral/ and Harvard²²2http://vision.seas.harvard.edu/hyperspec/explore.html are used to evaluate the effectiveness of the proposed method. The CAVE dataset contains 32 HSIs, each with 31 spectral bands and 512$\times$512 pixels, and we select 20 HSIs as the training dataset and 12 HSIs as the test dataset; the Harvard dataset contains 50 HSIs, each with 31 spectral bands and 1392$\times$1040 pixels, and we select 30 HSIs as the training dataset and 20 HSIs as the test dataset. Due to the limited computational resources, we crop and downsample each training HSI into 64$\times$64 pixels and each testing HSI into 256$\times$256 pixels.
Comparison methods and evaluation indicators: Peer methods compared with the proposed UHSR-AEC include three LLIE methods: LIME [19], RetinexNet [16], and EFINet [20]; and three HSI-SR methods: LTMR [12], LTTR [21], and MoG-DCN [22]. For consistency, the deep learning-based methods are all trained with the same HSI datasets.

In the experiments, all the HSI data are normalized to [0,1]. Four metrics are chosen to evaluate reconstructed HSI-SR image quality: peak signal-to-noise ratio (PSNR), structural similarity index (SSIM), spectral angle mapper (SAM), and the dimensionless global relative error of synthesis (ERGAS). All the experiments are performed on a server with a GPU NVIDIA 4090 graphics card.
Implementation Details: For each reference image, i.e., HR-HSI $\mathbf{Z}$ (GT), we adjust the image exposure by two different values of Gamma correction ($\alpha\in U(0.2,2)$, $\gamma\in U(0.5,3)$), for generating $\hat{\mathbf{Z}}_{x}$ and $\hat{\mathbf{Z}}_{y}$; the observation $\mathbf{X}$ is generated from $\hat{\mathbf{Z}}_{x}$ via Gaussian blurring (Gaussian kernel size of 8, mean of 0, standard deviation of $\sqrt{3}$) and downsampling (downsampling ratio of $K=|\mathbf{Z}|/|\mathbf{X}|=4$); the observation $\mathbf{Y}$ is generated through the convolutional integration of $\hat{\mathbf{Z}}_{y}$ with a real simulated spectral response (based on a Nikon camera).

We use the Pytorch framework to implement the proposed method. In the proposed UHSR-AEC network, we use a four-layer Unet structure and the Adam optimizer with the learning rate $l_{r}=10^{-4}$, as well as the maximum number of iterations $I_{max}=250,000$. Moreover, $\beta_{1},\beta_{2},\beta_{3},\lambda,\eta_{1},\eta_{2},T$ are set to 0.001, 0.001, 0.005, 0.5, 0.3, 0.1, 3, respectively. We would like to mention that the loss functions in (6) and (7) are for one HSI training sample $(\hat{\bf Z},\hat{\bf Z}_{x},\hat{\bf Z}_{y})$, while for the case of multiple HSI training samples in our experiment, they are simply replaced by the sum of the corresponding loss functions associated with each sample [22].

The simulation HSI and MSI data are generated using CAVE and Harvard datasets for the following two exposure degradation cases:

Case 1: ($\alpha_{1}$ = 0.5,$\gamma_{1}$ = 0.7) for HSI, and ($\alpha_{2}$ = 1.3,$\gamma_{2}$ = 1.5) for MSI.

Case 2: ($\alpha_{1}$ = 0.5,$\gamma_{1}$ = 2.0) for HSI, and ($\alpha_{2}$ = 0.8,$\gamma_{2}$ = 1.5) for MSI.

The obtained simulation results are shown in Table 1. As can be seen from this table, UHSR-AEC is overall superior to the other nine methods. For the visual quality assessment, Figs. 3 and 4 show some results for all tested methods on the CAVE and Harvard datasets, respectively. For each sub-image, the reconstructed HR-HSI for three spectral bands ([30,15,10]) is shown. From these sub-images, it can be observed that: i) LIME-based methods yield overly bright images, ii) RetinexNet and EFINet-based methods expose some detail loss and colour imbalance in Case 1 and Case 2, respectively, iii) UHSR-AEC shows the best visualization in both Case 1 and Case 2.

Some results for Ablation experiments on the proposed UHSR-AEC are shown in Table 2, which are obtained using the CAVE dataset with exposure degradation parameters used in Case 1. One can see from this table, that the proposed UHSR-AEC performs worst without IM, while it (with IM) performs best only in terms of PSNR for $T=3$, and best for $T=4$ otherwise, indicating that IM is essential and a suitable value of $T$ is also needed for the proposed UHSR-AEC to operate in good shape.