A Constrained Convex Optimization Approach to Hyperspectral Image Restoration with Hybrid Spatio-Spectral Regularization

: We propose a new constrained optimization approach to hyperspectral (HS) image restoration. Most existing methods restore a desirable HS image by solving some optimization problems, consisting of a regularization term(s) and a data-ﬁdelity term(s). The methods have to handle a regularization term(s) and a data-ﬁdelity term(s) simultaneously in one objective function; therefore, we need to carefully control the hyperparameter(s) that balances these terms. However, the setting of such hyperparameters is often a troublesome task because their suitable values depend strongly on the regularization terms adopted and the noise intensities on a given observation. Our proposed method is formulated as a convex optimization problem, utilizing a novel hybrid regularization technique named Hybrid Spatio-Spectral Total Variation (HSSTV) and incorporating data-ﬁdelity as hard constraints. HSSTV has a strong noise and artifact removal ability while avoiding oversmoothing and spectral distortion, without combining other regularizations such as low-rank modeling-based ones. In addition, the constraint-type data-ﬁdelity enables us to translate the hyperparameters that balance between regularization and data-ﬁdelity to the upper bounds of the degree of data-ﬁdelity that can be set in a much easier manner. We also develop an efﬁcient algorithm based on the alternating direction method of multipliers (ADMM) to efﬁciently solve the optimization problem. We illustrate the advantages of the proposed method over various HS image restoration methods through comprehensive experiments, including state-of-the-art ones.


Introduction
Hyperspectral (HS) imagery has 1D spectral information, including invisible light and narrow wavelength interval, in addition to 2D spatial information and thus can visualize unseen intrinsic characteristics of scene objects and environmental lighting. This makes HS imaging a key technique in many applications in various fields, e.g., earth observation, agriculture, and medical and biological imaging [1][2][3].
Observed HS images are often affected by noise because of the small amount of light in narrow wavelength and/or sensor failure. In addition, in compressive HS imaging scenarios [4,5], we have to estimate a full HS image from a minimal number of measurements. Thus, we need some methods for restoring desirable HS images from such degraded observations in HS applications.
Most HS image restoration methods are established based on optimization: a desirable HS image is characterized as a solution to some optimization problem, which consists of a regularization term and a data-fidelity term. The regularization term evaluates apriori knowledge about underlying properties on HS images, and the data-fidelity term keeps the consistency with a given observation. Thanks to the design, these methods get a reasonable result under ill-posed or ill-conditioned scenarios typical in HS image restoration.
Regularization techniques for HS image restoration are roughly classified into two groups: total variation (TV)-based approach and low-rank modeling (LRM)-based one. TV models the total absolute magnitude of local differences to exploit the piecewise-smooth structures of an image. Many TV-based approaches [6][7][8][9] have been proposed for HS image restoration. Besides, LRM-based approaches [10,11] exploit the underlying low-rank structure in the spectral direction of an HS image. A popular example is the so-called Low-rank matrix recovery (LRMR) [10].
Many recent methods [9,[12][13][14][15][16][17][18][19][20][21][22] combine TV-based and LRM-based approaches, and in general, they perform better than approaches using either regularization. This is because TV-based approaches model the spatial structure of an HS image, whereas LRM-based approaches the spectral one. Naturally, the methods have to handle multiple regularization terms and a data-fidelity term(s) simultaneously in one objective function. So the methods must carefully control the hyperparameter(s) balancing these terms. Specifically, such hyperparameters are interdependent, which means that a suitable value of a hyperparameter varies depending on the multiple regularization terms used and the noise intensities on a given observation. Hence, the hyperparameter settings in such combined approaches are often troublesome tasks. Table 1 summarizes the features of the methods reviewed in this section. Table 1. The feature of existing methods for hyperspectral (HS) image restoration.

Feature Spatial Correlation Spectral Correlation Convexity Hyperparameters
HTV [6] × convex interdependent SSAHTV [6] convex interdependent SSTV [7] convex interdependent ASSTV [8] convex interdependent LRM [10,11] × nonconvex independent LNWTV + LRM [9,12] convex interdependent HTV + LRM [13] nonconvex interdependent HTV + LRM [14,15] convex interdependent ASSTV + LRM [16,17] nonconvex interdependent SSTV + LRM [18][19][20] convex interdependent SSTV + LRM [21,22] nonconvex interdependent proposed convex independent Based on the above discussion, we propose a new constrained convex optimization approach to HS image restoration. Our proposed method restores a desirable HS image by solving a convex optimization problem involving a new TV-based regularization and hard constraints on data-fidelity. The regularization, named Hybrid Spatio-Spectral Total Variation (HSSTV), is designed to evaluate two types of local differences: direct local spatial differences and local spatio-spectral differences in a unified manner to effectively exploit both the underlying spatial and spectral structures of an HS image. Thanks to this design, HSSTV has a strong noise and artifact removal ability while avoiding oversmoothing and spectral distortion without combining LRM. Moreover, the constrained-type data-fidelity in the proposed method enables us to translate interdependent hyperparameters to the upper bounds of the degree of data-fidelity that can be determined based only on the noise intensity. As a result, the proposed method has no interdependent hyperparameter. We also develop an efficient algorithm for solving the optimization problem based on the well-known alternating direction method of multipliers (ADMM) [23][24][25].
The remainder of the paper is organized as follows. Section 2 introduces notation and mathematical ingredients. Section 3 reviews existing methods related to our method. In Section 4, we define HSSTV, formulate HS image restoration as a convex optimization problem involving HSSTV and hard-constraints on data-fidelity, and present an ADMM-based algorithm. Extensive experiments on denoising and compressed sensing (CS) reconstruction of HS images are given in Section 5, where we illustrate the advantages of our method over several state-of-the-art methods. In Section 6, we discuss the impact of the parameters in our proposed method. Section 7 concludes the paper. We have published the preliminary versions of this work in conference proceedings [26,27], including only limited experiments and discussions. In this paper, we newly conduct real noise removal experiments (see Section 5.2) and give much deeper discussions on the performance, utility, and parameter settings of the proposed method.

Notation and Definitions
In this paper, let R be the set of real numbers. We shall use boldface lowercase and capital to represent vectors and matrices, respectively, and := to define something. We denote the transpose of a vector/matrix by (·) , and the Euclidean norm (the 2 norm) of a vector by · .
For notational convenience, we treat an HS image U ∈ R N v ×N h ×B as a vector u ∈ R NB (N := N v N h is the number of the pixels of each band, and B is the number of the bands) by stacking its columns on top of one another, i.e., the index of the component of the ith pixel in kth band is i + (k − 1)N (for i = 1, . . . , N and k = 1, . . . , B).

Proximal Tools
for every x, y ∈ R L and λ ∈ (0, 1), respectively. Let Γ 0 (R L ) be the set of all proper lower semicontinuous convex functions on R L .
The proximity operator [28] plays a central role in convex optimization based on proximal splitting. The proximity operator of f ∈ Γ 0 (R L ) with an index γ > 0 is defined by We introduce the indicator function of a nonempty closed convex set C ⊂ R L , which is defined as follows: Then, for any γ > 0, its proximity operator is given by where P C (x) is the metric projection onto C.

Alternating Direction Method of Multipliers (ADMM)
ADMM [23][24][25] is a popular proximal splitting method, and it can solve convex optimization problems of the form: min where f ∈ Γ 0 (R L 1 ), g ∈ Γ 0 (R L 2 ), and G ∈ R L 2 ×L 1 . Here, we assume that f is quadratic, g is proximable, i.e., the proximity operator of g is computable in an efficient manner, and G is a full-column rank matrix. For arbitrarily chosen z (0) , d (0) and a step size γ > 0, ADMM iterates the following steps: The convergence property of ADMM is given as follows.
Theorem 1 (Convergence of ADMM [24]). Consider Problem (1), and assume that G G is invertible and that a saddle point of its unaugmented Lagrangian L 0 (x, z, y) : Then the sequence (x n ) n>0 generated by (2) converges to an optimal solution to Problem (1).

Related Works
In this section, we elaborate on existing HS image restoration methods based on optimization.

TV-Based Methods
The methods proposed in [6][7][8] restore a desirable HS image by solving a convex optimization problem involving TV-based regularization. Letū ∈ R NB be the desirable HS image, and the authors assume that an observation v ∈ R NB is modeled as follows: where n and s are an additive white Gaussian noise and a sparse noise, respectively. Here, the sparse noise corrupts only a few pixels in the HS image but heavily, e.g., impulse noise, salt-and-pepper noise, and line noise. The observation and the restoration problem of the methods are given by the following forms: min where R TV is a regularization function based on TV, and λ 1 and λ 2 are hyperparameters. Here, The first and third terms evaluate data-fidelity on Gaussian and sparse noise, respectively. The hyperparameters λ 1 and λ 2 represent the priorities of each term. If we can choose suitable values of the hyperparameters, then this formulation yields high-quality restoration. However, the hyperparameters are interdependent, which means that suitable values of the hyperparameters vary depending on the used TV-based regularization term and the noise intensities on a given observation. Therefore, the settings of the hyperparameters are an essential but troublesome task. In the following, we explain each TV. Let D = (D v D h ) ∈ R 2NB×NB be spatial differences operator with D v and D h being vertical and horizontal differences operator, respectively, and spectral differences operator are D b ∈ R NB×NB . In [6][7][8], HTV, ASSTV, and SSTV are defined as follows: where · TV is a TV norm, which takes the 2 norm of spacial difference vectors for all band and then summing up for all spatial pixels, and τ v , τ h , and τ b are the weight of the vertical, horizontal, and spectral differences. HTV evaluates direct spatial piecewise-smoothness and can be seen as a generalization of the standard color TV [29]. HTV does not consider spectral correlation, resulting in spatial oversmoothing. To consider the spectral correlation, the authors of [6] proposed SSAHTV. SSAHTV is a weighted HTV, and the weight is determined by spectral information. However, since SSAHTV does not directly evaluate spectral correlation, it still causes spatial oversmoothing. ASSTV evaluates direct spatial and spectral piecewise-smoothness ( Figure 1, blue line). The weights τ v , τ h , and τ b in (5) balance the smoothness of vertical, horizontal, and spectral differences, respectively. Owing to the definition, ASSTV can evaluate spatial and spectral correlation, but it produces spectral oversmoothing even if we carefully adjust τ v , τ h , and τ b . SSTV evaluate a-prior knowledge on HS images using spatio-spectral piecewise-smoothness. It is derived by calculating spatial differences through spectral differences ( Figure 1, yellow line). SSTV can restore a desirable HS image without any weight, but it produces noise-like artifacts, especially when a given observation is contaminated by heavy noise and/or degradation. Figure 1. Calculation of local differences in SSTV, ASSTV, and our HSSTV. SSTV evaluates the 1 norm of spatio-spectral differences (yellow line). ASSTV evaluates the 1 norm of direct spatial and spectral differences (blue line). HSSTV evaluates the mixed 1,p norm of direct spatial and spatio-spectral differences (red line).

LRM-Based Method
LRMR [10] is one of the popular LRM-based methods for HS image restoration, which evaluates the low rankness of an HS image in the spectral direction. To preserve the local details, LRMR restores a desirable HS image through patch-wise processing. Each patch is a local cube of the size of q × q × B, and LRMR handles it as a matrix of size q 2 × B obtained by lexicographically arranging the spatial vectors in the patch cube in the row direction. The observation model is expressed like Section 3.1, and the restoration problem is formulated as follows: where U i,j , V i,j , and S i,j represent the patches of a restored HS image, an observation, and a sparse noise, respectively, which are centered at (i,j) pixel. Then, the · F is a Frobenius norm, rank(·) represents a rank function, and card(·) is a cardinality function. The method evaluates the low rankness of the estimated HS image and sparsity of the sparse noise by limiting the number of the rank of U i,j and the cardinality of S i,j using r and k, respectively. Thanks to the design, LRMR achieves high-quality restoration for especially spectral information. Meanwhile, since LRMR does not fully consider spatial correlation, the result by LRMR tends to have spatial artifacts when an observation is corrupted by heavy noise and/or degradation. Besides, the rank and cardinality functions are nonconvex, and so it is a troublesome task to seek the global optimal solution of Problem (7).
However, the methods have to handle multiple regularization terms and/or a data-fidelity term(s) simultaneously in one objective function; they require to carefully control the hyperparameters balancing these terms. Since the hyperparameters rely on both the regularizations and the noise intensity on an observation, i.e., the hyperparameters are interdependent, the hyperparameter settings are often troublesome tasks.

Hybrid Spatio-Spectral Total Variation
We propose a new regularization technique for HS image restoration, named HSSTV. HSSTV simultaneously handles both direct local spatial differences and local spatio-spectral differences of an HS image. Then, HSSTV is defined by where · 1,p is the mixed 1,p norm, and ω ≥ 0. We assume p = 1 or 2, i.e., the 1 norm ( · 1,1 = · 1 ) or the mixed 1,2 norm, respectively. We would like to mention that we can also see 1 -HSSTV (p = 1) as anisotropic HSSTV and 1,2 -HSSTV (p = 2) as isotropic HSSTV.
In (8), DD b u and Du correspond to local spatio-spectral and direct local spatial differences, respectively, as shown in Figure 1 (red lines). The weight ω adjusts the relative importance of direct spatial piecewise-smoothness to spatio-spectral piecewise-smoothness. HSSTV evaluates two kinds of smoothness by taking the p norm (p = 1 or 2) of these differences associated with each pixel and then summing up for all pixels, i.e., calculating the 1 norm. Thus, it can be defined via the mixed 1,p norm. When we set ω = 0 and p = 1, HSSTV recovers SSTV as (6), meaning that HSSTV can be seen as a generalization of SSTV.
As reviewed in Section 3, since SSTV only evaluates spatio-spectral piecewise-smoothness, it cannot remove similar noise in adjacent bands. The direct spatial differences in HSSTV help to remove such noise. Figure 2 is restored HS images from an observation contaminated by similar noise in adjacent bands (the upper half area) and random noise (the lower half area). One can see that large noise remains in the upper half area of the result by SSTV. In contrast, HSSTV effectively removes all noise. However, since minimizing the direct spatial differences strongly promotes spatial piecewise-smoothness, HSSTV produces spatial oversmoothing when the weight ω is large. Thus, the weight ω should be set to less than one, as demonstrated in Section 5.

HS Image Restoration by HSSTV
We consider restoring a desirable HS imageū ∈ R NB from an observation v ∈ R M (M ≤ NB) contaminated by a Gaussian-sparse mixed noise. The observation model is given by the following form: where Φ ∈ R M×NB is a matrix representing a linear observation process, e.g., random sampling, n ∈ R M is Gaussian noise with the standard deviation σ, and s ∈ R M is a sparse noise. Based on the above model, we formulate HS image restoration using HSSTV as the following optimization problem: where B v 2,ε is a v-centered 2 -norm ball with the radius ε > 0, B 1,η is a 0-centered 1 -norm ball with the radius η > 0, and [µ min , µ max ] NB is a dynamic range of an HS image (µ min < µ max ). This method simultaneously estimates the desirable HS image u and the sparse noise s for noise-robust restoration. The first and second constraints measure data fidelities to the observation v and the sparse noise s, respectively. As mentioned in [13,16,26,[31][32][33][34][35][36][37][38][39], such a constraint-type data-fidelity enables us to translate the hyperparameter(s) balancing between regularization and data-fidelity like λ 1 and λ 2 in (3) to the upper bound of the degree of data-fidelity ε and η that can be set in a much easier manner.
Since all constraints are closed convex sets and HSSTV is a convex function, Problem (10) is a constrained convex optimization problem. In this paper, we adopt ADMM (see Section 2.3) for solving the problem. In what follows, we reformulate Problem (10) into Problem (1).
By using the indicator functions of the constraints, Problem (10) can be rewritten as Note that from the definition of the indicator function, Problem (11) is exactly equal to Problem (10). By letting Problem (11) is reduced to Problem (1). The resulting algorithm based on ADMM is summarized in Algorithm 1.

end
The update of u and s in Algorithm 1 comes down to the following forms: Since the update of u and s in Algorithm 1 is strictly-convex quadratic minimization, one can obtain these updated forms by differentiating it. Here, we should consider the structure of Φ because it affects the matrix inversion in (15). If Φ is a block-circulant-with-circulant-blocks (BCCB) matrix [40], we can leverage 3DFFT to efficiently solve the inversion in Step 2 with the difference operators having a periodic boundary, i.e., A ω A ω + Φ Φ + I can be diagonalized by the 3D FFT matrix and its inverse. If Φ is a semiorthogonal matrix, i.e., ΦΦ = αI (α > 0), we leave it to the update of z 2 , which means that we replace ι B v 2,ε by ι B v 2,ε • Φ in (13) and Φu by u in (14). This is because the proximity operator of ι B v 2,ε • Φ, in this case, can be computed by using [41] (Table 1.1-x) as follows: If Φ is a sparse matrix, we offer to use a preconditioned conjugate gradient method [42] for approximately solving the inversion or to apply primal-dual splitting methods [43][44][45] instead of ADMM (Primal-dual splitting methods require no matrix inversion but in general their convergence speed is slower than ADMM. Otherwise, randomized image restoration methods using stochastic proximal splitting algorithms [46][47][48][49] might be useful for reducing the computational cost. For the update of z 1 , the proximity operators are reduced to simple soft-thresholding type operations: for γ > 0 and for i = 1, . . . , 4NB, (i) in the case of p = 1, where sgn is the sign function, and (ii) in the case of p = 2, The updates of z 2 , z 3 , and z 4 require the proximity operators of the indicator functions of B v 2,ε , B 1,η , and [µ min , µ max ] NB , respectively, which equate the metric projections onto them (see Section 2.2). Specifically, the metric projection to B v 2,ε is given by x−v , otherwise, that onto B 1,η is given by and that onto [µ min , µ max ] NB is given, for i = 1, . . . , NB, by

Results
We demonstrate the advantages of the proposed method by applying it to two specific HS image restoration problems: denoising and CS reconstruction. In these experiments, we used 13 HS images taken from the SpecTIR [50], MultiSpec [51], and GIC [52], with their dynamic ranges normalized into [0, 1]. All the experiments were conducted by MATLAB 2018a.
The proposed method was compared with HTV [6], SSAHTV [6], SSTV [7], and ASSTV [38]. For a fair comparison, we replaced HSSTV in Problem (10) with HTV, SSAHTV, SSTV, or ASSTV and solved the problem by ADMM. In the denoising experiments, we also compared our proposed method with LRMR [10], TV-regularized low-rank matrix factorization (LRTV) [13], and local low-rank matrix recovery and global spatial-spectral TV (LLRGTV) [16]. Since LRMR, LRTV, and LLRGTV are customized to the mixed noise removal problem, we cannot adopt them for CS reconstruction. We did not compare our proposed method with recent CNN-based HS image denoising methods [53,54]. The CNN-based methods cannot be represented as explicit regularization functions and are fully customized to denoising tasks. In contrast, our proposed method can be used as a building block in various HS image restoration methods based on optimization. Meanwhile, CNN-based methods strongly depend on what training data are used, which means that they cannot adapt to a wide range of noise intensity. Thus, the design concepts of these methods are different from TVs and LRM-based approaches.
To quantitatively evaluate restoration performance, we used the peak signal-to-noise ratio (PSNR) [dB] index and the structural similarity (SSIM) [55] index between a groundtruth HS imagē u and a restored HS image u. PSNR is defined by 10 log 10 (MAX 2 I /MSE), where MAX I is the max value of the dynamic range of HS images, and MSE is the mean square error between groundtruth HS images and restored HS images. The higher PSNR is, the more similar the two images are. SSIM is an image quality assessment index based on the human vision system, which is defined as follows: where u i andū i are the ith pixel-centered local patches of a restored HS image and a groundtruth HS image, respectively, P is the number of patches, µ u i and µū i is the average values of the local patches of the restored and groundtruth HS images, respectively, σ u i and σū i represent the variances of u i and u i , respectively, and σ u iūi denotes the covariance between u i andū i . Moreover, C 1 and C 2 are two constants, which avoid the numerical instability when either µ 2 u i + µ 2 u i or σ 2 u i + σ 2 u i is very close to zero.
SSIM gives a normalized score between zero and one, where the maximum value means that u equals toū.

Denoising
First, we experimented on the Gaussian-sparse mixed noise removal of HS images. The observed HS images included an additive white Gaussian noise n with the standard deviation σ and sparse noise s. In these experiments, we assumed that sparse noise consists of salt-and-pepper noise and vertical and horizontal line noise, with these noise ratio in all pixels being s p , l v , and l h , respectively. We generated noisy HS images by adding two types of mixed noise to groundtruth HS images: (i) (σ, s p , l v = l h ) = (0.05, 0.04, 0.04), (ii) (0.1, 0.05, 0.05). In the denoising case, Φ = I in (9), and the radiuses ε and η in Problem (10) were set to 0.83 where v ave is the average of the observed image. Table 2 shows the parameters settings for ASSTV, LRMR, LRTV, LLRGTV, and the proposed method. We set these parameters to achieve the best performance for each method.  In Table 3, we show the PSNR and SSIM of the denoised HS images by each method for two types of noise intensity and HS images. For HTV, SSAHTV, SSTV, ASSTV, and LRMR, the proposed method outperforms the existing methods. Besides, one can see that the PSNR and SSIM of the results by the proposed method are higher than those by LLRGTV for most situations. Meanwhile, LRTV outperforms the proposed method in some cases. This would be because LRTV combines TV-based and LRM-based regularization techniques. We would like to mention that the proposed method outperforms LRTV in over half of the situations, even though it uses only TV-based regularization. Figure 3 shows the resulting images on Salinas (the noise level (i), top) and PaviaU (the noise level (ii), bottom), with their PSNR (left) and SSIM (right). Here, we depicted these HS images as RGB images (R = 8th, G = 16th, and B = 32nd bands). One can see that the results by HTV, SSAHTV, and ASSTV lose spatial details, and noise remains in the results by SSTV, LRMR, and LLRGTV. Besides, since the restored images by SSTV and LRTV lose color with large noise intensity, SSTV and LRTV change spectral variation. In contrast, the proposed method can restore HS images preserving both details and spectral information without artifacts.   , and spatial and spectral responses (right) of the denoised Suwannee HS image in the case of the noise level (ii). The graphs regarding band-wise PSNR and SSIM show that the proposed method achieves higher-quality restoration than HTV, SSAHTV, and SSTV for all bands and ASSTV and LRMR for most bands. Besides, even though the proposed method only utilizes HSSTV, the results by the proposed method outperform those by LRTV for some bands of the SSIM case and those by LLRGTV for many bands. Graph (c) plots the spatial response of the 243rd row of the 30th band. Similarly, graph (d) plots the spectral response of the 243rd row and 107th col. We can see that the spatial response of the results by HTV and SSAHTV is too smooth compared with the groundtruth one. On the other hand, there exist undesirable variations in the spatial response of the results by SSTV, LRMR, and LLRGTV. In contrast, ASSTV, LRTV, and the proposed method restore similar responses to the groundtruth one. In graph (d), one can see that (i) HTV, SSAHTV, and LRMR produce spectral artifacts, (ii) the shape of the spectral responses of the results by SSTV is similar to that of the groundtruth one, but the mean value is larger than the groundtruth one, (iii) the spectral response of the results by ASSTV is too smooth and different from the groundtruth one, and (iv) LRTV, LLRGTV, and the proposed method can restore a spectral response very similar to the groundtruth one.

Real Noise Removal
We also examined HTV, SSAHTV, SSTV, ASSTV, LRMR, LRTV, LLRGTV, and the proposed method on an HS image with real noise. We selected noisy 16 bands from Suwannee and used it as a real observed HS image v. To maximize the performance of each method, we searched for suitable values of σ, s p , l v , and l h , and we set the parameters as with Section 5.1. Specifically, we set the parameter ε = 3.1893 and η = 31893 for all TVs. Besides, the parameters τ v , τ h , and τ b in ASSTV are set as 1, 1, 3, respectively, r and k in LRMR are set as 3 and 0.1204, respectively, r and tau in LRTV are set as 2 and 0.008, and r, λ, and τ in LLRGTV are set as 2, 0.1 and 0.01. Figure 5 shows the results, where the HS images are depicted as RGB images (R = 2nd, G = 6th, and B = 13rd bands). The results by HTV and SSAHTV have spatial oversmoothing, and SSTV, ASSTV, LRMR, and LLRGTV produce spatial artifacts. Besides, one can see that the results by LRMR and LRTV have spectral artifacts. On the other hand, the proposed method can restore a detail-preserved HS image without artifacts.

Compressed Sensing Reconstruction
We experimented on compressed sensing (CS) reconstruction [56,57]. The CS theory says that high-dimensional signal information can be reconstructed from incomplete random measurements by exploiting sparsity in some domains, e.g., the gradient domain (TV). In general, HS imaging captures an HS image by scanning 1D spatial and spectral information because it senses spectral information by dispersing the incident light. Therefore, capturing moving objects is very difficult in HS imaging. To overcome the drawback, one-shot HS imaging based on CS has been actively studied [4,5].
In this experiment, we assume that Φ ∈ R M×NB in (9) is a random sampling matrix with the sampling rate m (0 < m < 1 and M = mNB). Here, since Φ is a semiorthogonal matrix, we can efficiently solve the problem, as explained in Section 4.2. Moreover, since the main objective in the experiments is to verify CS reconstruction performance by HSSTV, we assume that the observations are contaminated by only an additive white Gaussian noise n with noise intensity σ = 0.1.
We set the CS reconstruction problem as follows: The problem is derived by removing the second constraint and s from Problem (10). Therefore, we can solve the above problem by removing s, z 3 , and d 3 in Algorithm (1) and replacing z 4 and d 4 with z 3 and d 3 , respectively. As in Section 4.2, the update of u is strictly-convex quadratic minimization, and so it comes down to We set m = 0.2 or 0.4 and ε = √ mNBσ 2 . In the ASSTV case, we set the parameters (τ v , τ h , τ b ) = (1, 1, 0.5), which experimentally achieves the best performance. Table 4 shows PSNR and SSIM of the reconstructed HS images. For all m and HS images, both PSNR and SSIM of the proposed method results are almost higher than that by HTV, SSAHTV, SSTV, and ASSTV.  Figure 6 is the reconstructed results on KSC and Reno with the random sampling ratio m = 0.4 and 0.2, respectively. Here, the HS images are depicted as RGB images (R = 8th, G = 16th, and B = 32nd bands). One can see that (i) HTV and SSAHTV cause spatial oversmoothing, (ii) SSTV produces artifacts and spectral distortion, appearing different from the color of the groundtruth HS images, and (iii) the results by ASSTV have spatial oversmoothing and spectral distortion. On the other hand, the proposed method reconstructs meaningful details without both artifacts and spectral distortion. Figure 7 plots band-wise PSNR or SSIM (left) and spatial and spectral responses (right) (Suwannee, m = 0.2). According to band-wise PSNR and SSIM, one can see that the proposed method achieves higher-quality reconstruction for all bands than HTV, SSAHTV, SSTV, and ASSTV. Graphs (c) and (d) plot the spatial and spectral responses of the same position in Section 5.1. Graph (c) shows that (i) the spatial response of the results by HTV, SSAHTV, and ASSTV are oversmoothing, (ii) SSTV produces undesirable variation, and (iii) the spatial response reconstructed by the proposed method is similar to the groundtruth one. In graph (d), HTV and SSAHTV generate undesirable variation, and ASSTV causes oversmoothing. Thanks to the evaluation of spatio-spectral piecewise-smoothness, SSTV reconstructs a similar spectral response to the groundtruth one. However, the mean value is larger than the groundtruth value. The proposed method achieves the most similar reconstruction of spectral response among all the TVs. KSC 18

Discussion
In this section, we discuss multiple parameters in our proposed method.

The Impact of The Weight ω
The weight ω in (8) balances spatio-spectral piecewise-smoothness and direct spatial piecewise smoothness. We verified the impact ω on the mixed noise removal and CS reconstruction experiments. The weight ω was varied from 0.01 to 0.2 at 0.01 interval. The other settings were the same as Sections 5.1 and 5.3. Figures 8 and 9 plot PSNR or SSIM of the results by the proposed method versus various ω on the mixed noise removal and CS reconstruction experiments, respectively, where the values of PSNR and SSIM are averaged over the 13 HS images. One can see that ω ∈ [0.03, 0.07] is a good choice in the mixed noise removal case, and ω ∈ [0.05, 0.1] is a good choice in most cases in the CS case. The results show that the suitable range of ω in CS reconstruction is almost the same as mixed noise removal.
ASSTV and LRTV require adjusting the weight τ b and the hyperparameter τ newly for difference noise intensity, respectively, but the suitable parameter ω in HSSTV is noise-robust.

The Sensitivity of The Parameters ε and η
To verify the sensitivity of the parameter ε and η in (10), we conducted additional mixed noise removal experiments, where we examined various values of ε and η. Specifically, we set ε = 0.83α NB(1 − (s p (1 − l v − l h ) + l v + l h − l v l h ))σ 2 and η = βNB(0.45s p + (l v + l h )v ave − l v l h v ave ), which are hand-optimized values of the parameters and changed α and β from 0.9 to 1.1 at 0.02 interval (the DC and the KSC images and the noise level (ii)). Figure 10 plots PSNR or SSIM of the results by HTV, SSAHTV, SSTV, ASSTV, and the proposed method versus α or β. For HTV, SSAHTV, ASSTV, and the proposed method, the graphs show that the suitable values of α and β do not vary significantly for both image, so the parameters ε and η are independent of both a regularization technique and an observed image. In the SSTV cases, the shapes of the plots are different between KSC and DC. This is because in the DC case, SSTV converges for all parameter settings, while in the case of KSC, it does not converge when α, β > 1.

Conclusions
We have proposed a new constrained optimization approach to HS image restoration. Our proposed method is formulated as a convex optimization problem, where we utilize a novel regularization technique named HSSTV and incorporate data-fidelity as hard constraint. HSSTV evaluates direct spatial piecewise-smoothness and spatio-spectral piecewise-smoothness, and so it has a strong ability of HS restoration. Thanks to the design of the constraint-type data-fidelity, we can independently set the hyperparameters that balance between regularization and data-fidelity. To solve the proposed problem, we develop an efficient algorithm based on ADMM. Experimental results on mixed noise removal, real noise removal, and CS reconstruction demonstrate the advantages of the proposed method over various HS image restoration methods.

Conflicts of Interest:
The authors declare that there is no conflict of interest.