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

TL;DR

This paper introduces a convex optimization method for hyperspectral image restoration that uses a novel hybrid spatio-spectral regularization and hard data-fidelity constraints, simplifying hyperparameter tuning and improving restoration quality.

Contribution

The paper proposes a new convex optimization framework with hybrid regularization and data-fidelity constraints, avoiding complex hyperparameter tuning and enhancing noise/artifact removal in hyperspectral images.

Findings

Outperforms state-of-the-art hyperspectral restoration methods

Effectively removes noise and artifacts without oversmoothing

Simplifies hyperparameter setting through constraint-based data-fidelity

Abstract

We propose a new constrained optimization approach to hyperspectral (HS) image restoration. Most existing methods restore a desirable HS image by solving some optimization problem, which consists of a regularization term(s) and a data-fidelity term(s). The methods have to handle a regularization term(s) and a data-fidelity term(s) simultaneously in one objective function, and so 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, where we utilize a novel hybrid regularization technique named Hybrid Spatio-Spectral Total Variation (HSSTV) and incorporate data-fidelity as hard constraints. HSSTV has a strong ability of noise and artifact removal while avoiding oversmoothing and spectral distortion, without combining other regularizations such as low-rank modeling-based ones. In addition, the constraint-type data-fidelity enables us to translate the hyperparameters that balance between regularization and data-fidelity to the upper bounds of the degree of data-fidelity that can be set in a much easier manner. We also develop an efficient algorithm based on the alternating direction method of multipliers (ADMM) to efficiently solve the optimization problem. Through comprehensive experiments, we illustrate the advantages of the proposed method over various HS image restoration methods including state-of-the-art ones.

Tables

Table 1. Table 1: The feature of existing methods for HS image restoration.

	spatial correlation	spectral correlation	convexity	hyperparameters
HTV Yuan et al. (2012)	$○$	$\times$	convex	interdependent
SSAHTV Yuan et al. (2012)	$○$	$△$	convex	interdependent
SSTV Aggarwal and Majumdar (2016)	$△$	$○$	convex	interdependent
ASSTV Chang et al. (2015)	$○$	$○$	convex	interdependent
LRM Zhang et al. (2014)	$\times$	$○$	nonconvex	independent
LNWTV + LRM Liu et al. (2018); Li et al. (2018)	$○$	$○$	convex	interdependent
HTV + LRM He et al. (2016)	$○$	$○$	nonconvex	interdependent
ASSTV + LRM He et al. (2018); Kong et al. (2020)	$○$	$○$	nonconvex	interdependent
SSTV + LRM Cao et al. (2018); Wang et al. (2017, 2018)	$○$	$○$	convex	interdependent
SSTV + LRM Sun et al. (2018); Ince (2019)	$○$	$○$	nonconvex	interdependent
proposed	$○$	$○$	convex	independent

Table 2. Table 2: Parameter settings for ASSTV, LRMR, LRTV, and the proposed method.

		(i) $\left(\sigma ,s_p,l_v=l_h\right)=(0.05,0.04,0.04)$	(ii) $(0.1,0.05,0.05)$
ASSTV	${\tau }_v={\tau }_h$	1
	${\tau }_b$	3	2
LRMR	$r$	3
	$k$	$s_p+l_v+l_h-l_v{l}_h$ (the rate of sparse noise)
LRTV	$r$	2
	$\tau$	0.005	0.008
proposed	$\omega$	0.04

Table 3. Table 3: PSNR (top) and SSIM (bottom) in mixed noise removal experiments.

	HS image	noise level	HTV	SSAHTV	SSTV	ASSTV	LRMR	LRTV	proposed ($p=1$)	proposed ($p=2$)
	Beltsville	(i)	29.43	29.47	33.66	27.16	30.91	35.32	34.25	34.16
	$256\times 256\times 32$	(ii)	26.40	26.43	28.42	24.60	27.13	31.22	29.79	29.62
	Suwannee	(i)	30.14	30.18	34.59	32.60	30.30	36.20	35.15	36.01
	$256\times 256\times 32$	(ii)	26.70	26.74	29.55	28.71	26.90	31.95	31.08	31.22
	DC	(i)	26.46	26.51	33.03	28.80	31.71	34.78	33.36	33.08
	$256\times 256\times 32$	(ii)	23.84	23.88	27.71	25.25	27.35	29.53	28.57	28.32
	Cuprite	(i)	31.67	31.68	34.42	29.14	30.16	28.39	34.96	36.20
	$256\times 256\times 32$	(ii)	28.20	28.21	29.86	26.57	27.32	27.94	31.63	31.73
	Reno	(i)	28.53	28.57	34.37	30.49	32.21	37.06	35.11	34.96
	$256\times 256\times 32$	(ii)	25.56	25.61	28.11	26.95	28.47	31.00	29.83	29.72
	Botswana	(i)	27.98	28.05	33.32	26.47	31.62	29.00	33.61	33.53
	$256\times 256\times 32$	(ii)	25.21	25.25	28.55	24.01	28.31	27.33	29.39	29.35
PSNR	IndianPines	(i)	31.05	31.06	31.45	29.07	28.96	26.16	31.90	31.80
	$145\times 145\times 32$	(ii)	28.57	28.57	27.82	26.72	25.14	29.82	29.26	29.18
	KSC	(i)	30.17	30.25	34.74	31.64	33.74	35.74	36.39	36.33
	$256\times 256\times 32$	(ii)	28.03	28.06	29.23	28.62	30.19	30.22	31.82	31.72
	PaviaLeft	(i)	27.62	27.70	35.57	30.91	33.01	36.49	35.98	35.81
	$216\times 216\times 32$	(ii)	24.74	24.78	29.93	26.71	29.46	29.02	30.47	30.24
	PaviaRight	(i)	26.93	27.35	34.54	31.13	33.33	35.82	35.68	35.23
	$256\times 256\times 32$	(ii)	24.90	25.16	30.70	27.23	29.82	29.08	31.59	31.39
	PaviaU	(i)	27.92	28.04	35.52	31.65	33.00	36.72	36.31	36.17
	$256\times 256\times 32$	(ii)	25.24	25.29	30.21	27.42	29.43	28.90	31.04	30.80
	Salinas	(i)	32.59	32.64	35.86	32.83	31.82	36.74	37.60	37.65
	$217\times 217\times 32$	(ii)	28.88	28.91	28.19	28.99	28.02	32.73	32.01	32.12
	SalinaA	(i)	32.54	32.65	35.29	28.12	31.18	28.49	36.27	36.23
	$83\times 86\times 32$	(ii)	28.69	28.80	29.67	25.19	27.67	26.10	31.68	31.64
		(i)	0.7902	0.7904	0.8856	0.8111	0.8583	0.9372	0.9132	0.9085
	Beltsville	(ii)	0.6954	0.6959	0.7057	0.7177	0.7083	0.8568	0.8186	0.8088
		(i)	0.8406	0.8410	0.9353	0.9052	0.8689	0.9502	0.9559	0.9555
	Suwannee	(ii)	0.7542	0.7552	0.8146	0.8226	0.7470	0.8930	0.9125	0.9158
		(i)	0.7622	0.7633	0.9274	0.8676	0.9248	0.9613	0.9442	0.9394
	DC	(ii)	0.6189	0.6201	0.8092	0.7211	0.8214	0.8810	0.8611	0.8533
		(i)	0.8550	0.8552	0.9179	0.8632	0.8495	0.9396	0.9459	0.9426
	Cuprite	(ii)	0.7849	0.7852	0.7717	0.7953	0.7098	0.8814	0.9031	0.9058
		(i)	0.7818	0.7819	0.9322	0.8832	0.9012	0.9589	0.9531	0.9515
	Reno	(ii)	0.6640	0.6645	0.8045	0.7539	0.7905	0.8816	0.8679	0.8635
		(i)	0.7896	0.7900	0.9202	0.8199	0.9068	0.9282	0.9343	0.9344
	Botswana	(ii)	0.6810	0.6820	0.8175	0.7095	0.8201	0.8564	0.8745	0.8765
		(i)	0.8118	0.8120	0.8015	0.7671	0.7593	0.8190	0.8335	0.8243
SSIM	IndianPines	(ii)	0.7713	0.7713	0.6229	0.7303	0.7893	0.7939	0.7785	0.7689
		(i)	0.8271	0.8278	0.9116	0.8922	0.8890	0.9385	0.9542	0.9532
	KSC	(ii)	0.7598	0.7602	0.7885	0.8064	0.7529	0.8427	0.8809	0.8747
		(i)	0.7752	0.7770	0.9593	0.8828	0.9359	0.9612	0.9661	0.9645
	PaviaLeft	(ii)	0.6102	0.6116	0.8755	0.7267	0.8565	0.8791	0.8898	0.8815
		(i)	0.7769	0.7772	0.9494	0.8862	0.9256	0.9540	0.9616	0.9598
	PaviaRight	(ii)	0.6474	0.6471	0.8635	0.7493	0.8261	0.8507	0.9086	0.9006
		(i)	0.7973	0.7986	0.9452	0.8891	0.9124	0.9540	0.9622	0.9610
	PaviaU	(ii)	0.6776	0.6785	0.8444	0.7678	0.8103	0.8627	0.8935	0.8855
		(i)	0.8997	0.9002	0.9015	0.9163	0.8270	0.9509	0.9561	0.9564
	Salinas	(ii)	0.8570	0.8575	0.7117	0.8732	0.6670	0.9225	0.9223	0.9240
		(i)	0.9129	0.9137	0.9134	0.8468	0.8632	0.9384	0.9448	0.9416
	SalinaA	(ii)	0.8793	0.8803	0.7789	0.8110	0.7266	0.8951	0.9197	0.9195

Table 4. Table 4: PSNR (left) and SSIM (right) in the CS reconstruction experiment.

		PSNR						SSIM
	$m$	HTV	SSAHTV	SSTV	ASSTV	proposed ($p=1$)	proposed ($p=2$)	HTV	SSAHTV	SSTV	ASSTV	proposed ($p=1$)	proposed ($p=2$)
	0.4	27.46	27.49	27.53	26.51	31.15	30.71	0.6829	0.6940	0.6013	0.6836	0.8105	0.7948
Beltsville	0.2	26.23	26.25	24.34	24.12	29.63	29.18	0.6363	0.6493	0.4348	0.6108	0.7604	0.7427
	0.4	27.97	28.02	28.49	27.68	32.98	33.04	0.7332	0.7497	0.7377	0.7367	0.8902	0.8909
Suwannee	0.2	26.47	26.50	25.69	25.39	31.37	31.44	0.6810	0.7007	0.5739	0.6633	0.8531	0.8534
	0.4	24.69	24.73	27.33	24.71	29.70	29.29	0.6096	0.6242	0.7522	0.6245	0.8577	0.8460
DC	0.2	23.31	23.33	24.16	22.69	27.98	27.59	0.5215	0.5384	0.6120	0.5037	0.7970	0.7846
	0.4	29.94	29.96	28.21	28.59	34.36	34.34	0.7665	0.7804	0.6826	0.7652	0.8882	0.8895
Cuprite	0.2	28.77	28.77	25.79	26.38	32.97	32.95	0.7368	0.7525	0.5057	0.7207	0.8568	0.8578
	0.4	26.99	27.05	27.82	26.49	31.80	31.61	0.6769	0.6868	0.7414	0.6730	0.8733	0.8705
Reno	0.2	25.57	25.61	25.57	24.52	30.22	30.04	0.6202	0.6326	0.6276	0.5940	0.8263	0.8228
	0.4	26.10	26.15	27.81	25.13	30.32	30.15	0.6683	0.6803	0.7551	0.6460	0.8563	0.8598
Botswana	0.2	24.66	24.69	24.79	22.86	28.79	28.63	0.6014	0.6162	0.6225	0.5519	0.8119	0.8163
	0.4	30.54	30.55	27.65	29.55	31.36	31.04	0.7497	0.7777	0.5066	0.7617	0.7806	0.7655
IndianPines	0.2	29.99	29.99	25.11	28.19	30.71	30.46	0.7366	0.7658	0.3488	0.7465	0.7589	0.7491
	0.4	29.30	29.33	28.34	28.63	34.10	34.03	0.7660	0.7742	0.6814	0.7544	0.9019	0.9002
KSC	0.2	28.11	28.12	27.00	26.59	32.67	32.60	0.7318	0.7410	0.6008	0.7032	0.8698	0.8679
	0.4	25.66	25.69	29.66	25.41	31.96	31.83	0.6082	0.6205	0.8386	0.5932	0.8900	0.8857
PaviaLeft	0.2	24.26	24.27	27.17	23.24	30.20	30.08	0.5103	0.5251	0.7319	0.4434	0.8418	0.8364
	0.4	25.83	25.85	29.86	25.61	32.45	32.21	0.6357	0.6423	0.7962	0.6275	0.8937	0.8877
PaviaRight	0.2	24.30	24.30	27.54	23.61	30.56	30.38	0.5502	0.5584	0.6917	0.5069	0.8475	0.8392
	0.4	26.49	26.53	30.02	26.38	32.88	32.70	0.6867	0.6956	0.7830	0.6818	0.8950	0.8901
PaviaU	0.2	24.95	24.97	27.35	24.08	31.13	30.96	0.6138	0.6242	0.6623	0.5680	0.8557	0.8508
	0.4	31.19	31.24	27.69	30.18	35.43	35.51	0.8577	0.8672	0.6153	0.8566	0.9222	0.9245
Salinas	0.2	29.94	29.98	25.28	28.09	34.05	34.10	0.8404	0.8516	0.4620	0.8302	0.9052	0.9080
	0.4	30.67	30.82	27.93	28.19	34.45	34.14	0.8647	0.8871	0.6595	0.8489	0.9178	0.9208
SalinasA	0.2	28.68	28.75	24.15	24.94	32.71	32.36	0.8387	0.8655	0.4810	0.8005	0.8966	0.9002

Equations

$${\mathrm{prox}}_{\gamma f}(\mathbf{x}):=\underset{\mathbf{y}}{\mathrm{argmin}}f(\mathbf{y})+\frac{1}{2\gamma }\Vert \mathbf{y}-\mathbf{x}{\Vert }^2.$$

$$\displaystyle\iota_{C}({\mathbf{x}}):=\left\{\begin{array}[]{l l}0,&\mbox{if}~{}{\mathbf{x}}\in C,\\ \infty,&\mbox{otherwise}.\\ \end{array}\right.$$

$${\mathrm{prox}}_{\gamma {\iota }_C}(\mathbf{x})=P_C(\mathbf{x}):=\underset{\mathbf{y}\in C}{\mathrm{argmin}}\Vert \mathbf{x}-\mathbf{y}\Vert ,$$

$$\underset{\mathbf{x},\mathbf{z}}{\min }f(\mathbf{x})+g(\mathbf{z})\text{mbox}s.t.\mathbf{z}=\mathbf{G}\mathbf{x},$$

$$\displaystyle\left\lfloor\begin{array}[]{l}{\mathbf{x}}^{(n+1)}=\mathop{\rm argmin}\limits_{{\mathbf{x}}}f({\mathbf{x}})+\frac{1}{2\gamma}\|{\mathbf{z}}^{(n)}-{\mathbf{G}}{\mathbf{x}}-{\mathbf{d}}^{(n)}\|^{2},\\ {\mathbf{z}}^{(n+1)}=\mbox{prox}_{\gamma g}({\mathbf{G}}{\mathbf{x}}^{(n+1)}+{\mathbf{d}}^{(n)}),\\ {\mathbf{d}}^{(n+1)}={\mathbf{d}}^{(n)}+{\mathbf{G}}{\mathbf{x}}^{(n+1)}-{\mathbf{z}}^{(n+1)},\end{array}\right.$$

$$\mathbf{v}=\bar{\mathbf{u}}+\mathbf{s}+\mathbf{n},$$

$$\underset{\mathbf{u},\mathbf{s}}{\min }\Vert \mathbf{v}-\mathbf{u}-\mathbf{s}{\Vert }^2+{\lambda }_1{\mathcal{R}}_{\mathrm{TV}}(\mathbf{u})+{\lambda }_2\Vert \mathbf{s}{\Vert }_1,$$

$$\mathrm{HTV}(\mathbf{u})$$

$$\mathrm{ASSTV}(\mathbf{u})$$

$$\mathrm{SSTV}(\mathbf{u})$$

$$\underset{{\mathbf{U}}_{i,j},{\mathbf{S}}_{i,j}}{\min }\Vert {\mathbf{V}}_{i,j}-{\mathbf{U}}_{i,j}-{\mathbf{S}}_{i,j}{\Vert }_F^{2}s.t.\mathrm{rank}({\mathbf{U}}_{i,j})\le r,\mathrm{card}({\mathbf{S}}_{i,j})\le k,$$

$$\operatorname{HSSTV}({\mathbf{u}}):=\|{\mathbf{A}}_{\omega}{\mathbf{u}}\|_{1,p}\mbox{ with }{\mathbf{A}}_{\omega}:=\left(\begin{array}[]{c}{\mathbf{D}}{\mathbf{D}}_{b}\\ \omega{\mathbf{D}}\\ \end{array}\right),$$

$$\mathbf{v}=\mathbf{\Phi }\bar{\mathbf{u}}+\mathbf{n}+\mathbf{s},$$

$$\underset{\mathbf{u},\mathbf{s}}{\min }\mathrm{HSSTV}(\mathbf{u})$$

$$\displaystyle\mbox{ s.t. }\left[\begin{array}[]{l}{\mathbf{\Phi}}{\mathbf{u}}+{\mathbf{s}}\in{\mathcal{B}}^{{\mathbf{v}}}_{2,\varepsilon}:=\{{\mathbf{x}}\in{\mathbb{R}}^{M}|\|{\mathbf{v}}-{\mathbf{x}}\|\leq\varepsilon\},\\ {\mathbf{s}}\in{\mathcal{B}}_{1,\eta}:=\{{\mathbf{x}}\in{\mathbb{R}}^{M}|\|{\mathbf{x}}\|_{1}\leq\eta\},\\ {\mathbf{u}}\in[\mu_{\min},\mu_{\max}]^{NB},\end{array}\right.$$

$$\underset{\mathbf{u},\mathbf{s}}{\min }\Vert {\mathbf{A}}_{\omega }\mathbf{u}{\Vert }_{1,p}+{\iota }_{{\mathcal{B}}_{2,\epsilon }^{\mathbf{v}}}(\mathbf{\Phi }\mathbf{u}+\mathbf{s})+{\iota }_{{\mathcal{B}}_{1,\eta }}(\mathbf{s})+{\iota }_{[{\mu }_{\min },{\mu }_{\max }{]}^{NB}}(\mathbf{u}).$$

$$f:{\mathbb{R}}^{NB}\to {\mathbb{R}}^2:\mathbf{u}\mapsto (0,0),$$

$$g:{\mathbb{R}}^{5NB+2M}\to \mathbb{R}\cup \{\infty \}:({\mathbf{z}}_1,{\mathbf{z}}_2,{\mathbf{z}}_3,{\mathbf{z}}_4)\mapsto$$

$$\Vert {\mathbf{z}}_1{\Vert }_{1,p}+{\iota }_{{\mathcal{B}}_{2,\epsilon }^{\mathbf{v}}}({\mathbf{z}}_2)+{\iota }_{{\mathcal{B}}_{1,\epsilon }}({\mathbf{z}}_3)+{\iota }_{[{\mu }_{\min },{\mu }_{\max }{]}^{NB}}({\mathbf{z}}_4),$$

$$\mathbf{G}:{\mathbb{R}}^{NB}\to {\mathbb{R}}^{5NB+2M}:\mathbf{u}\mapsto ({\mathbf{A}}_{\omega }\mathbf{u},\mathbf{\Phi }\mathbf{u}+\mathbf{s},\mathbf{s},\mathbf{u}).$$

$${\mathbf{u}}^{(n+1)}=$$

$$\text{mbox}RHS=$$

$$-\frac{1}{2}({\mathbf{z}}_3^{(n)}-{\mathbf{d}}_3^{(n)})+({\mathbf{z}}_4^{(n)}-{\mathbf{d}}_4^{(n)}),$$

$${\mathbf{s}}^{(n+1)}=$$

$${\mathrm{prox}}_{\gamma {\iota }_{{\mathcal{B}}_{2,\epsilon }^{\mathbf{v}}}\circ \mathbf{\Phi }}(\mathbf{x})=\mathbf{x}+{\alpha }^{-1}{\mathbf{\Phi }}^{\top }(P_{{\mathcal{B}}_{2,\epsilon }^{\mathbf{v}}}(\mathbf{\Phi }\mathbf{x})-\mathbf{\Phi }\mathbf{x}).$$

$$[{\mathrm{prox}}_{\gamma \Vert \cdot {\Vert }_1}(\mathbf{x}){]}_i=\mathrm{sgn}(x_i)\max \left\{|x_i|-\gamma ,0\right\},$$

$$[{\mathrm{prox}}_{\gamma \Vert \cdot {\Vert }_{1,2}}(\mathbf{x}){]}_i=\max \left\{1-\gamma {\left(\sum _{j=0}^3{x}_{\tilde{i}+jNB}^2\right)}^{-\frac{1}{2}},0\right\}{x}_i,$$

$$P_{{\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}}}({\mathbf{x}})=\left\{\begin{array}[]{l l}{\mathbf{x}},&\mbox{if }{\mathbf{x}}\in{\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}},\\ {\mathbf{v}}+\frac{\varepsilon({\mathbf{x}}-{\mathbf{v}})}{\|{\mathbf{x}}-{\mathbf{v}}\|},&\mbox{otherwise},\end{array}\right.$$

$$P_{{\mathcal{B}}_{1,\eta }}=\mathrm{sgn}(\mathbf{x})\max (|\mathbf{x}|-\eta ,0),$$

$$[P_{[{\mu }_{\min },{\mu }_{\max }{]}^{NB}}(\mathbf{x}){]}_i=\min \{\max \{x_i,{\mu }_{\min }\},{\mu }_{\max }\}.$$

$$\text{mbox}SSIM(\mathbf{u},\bar{\mathbf{u}})$$

$$\text{mbox}{SSIM}_i(\mathbf{u},\bar{\mathbf{u}})$$

$$\min_{{\mathbf{u}}}\operatorname{HSSTV}({\mathbf{u}})\mbox{ s.t. }\left[\begin{array}[]{l}{\mathbf{\Phi}}{\mathbf{u}}\in{\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}},\\ {\mathbf{u}}\in[\mu_{\min},\mu_{\max}]^{NB}.\end{array}\right.$$

$${\mathbf{u}}^{(n+1)}=({\mathbf{A}}_{\omega }^{\top }{\mathbf{A}}_{\omega }+{\mathbf{\Phi }}^{\top }\mathbf{\Phi }+\mathbf{I}{)}^{-1}\text{mbox}RHS,$$

$$\text{mbox}RHS={\mathbf{A}}_{\omega }({\mathbf{z}}_1^{(n)}-{\mathbf{d}}_1^{(n)})+{\mathbf{\Phi }}^{\top }({\mathbf{z}}_2^{(n)}-{\mathbf{d}}_2^{(n)})+({\mathbf{z}}_3^{(n)}-{\mathbf{d}}_3^{(n)}).$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date \TitleA Constrained Convex Optimization Approach

to Hyperspectral Image Restoration

with Hybrid Spatio-Spectral Regularization

\AuthorSaori Takeyama 1,∗\orcidA, Shunsuke Ono 1\orcidB, and Itsuo Kumazawa 1\orcidC \AuthorNamesSaori Takeyama, Shunsuke Ono and Itsuo Kumazawa

\corresCorrespondence: takeyama.s.aa@m.titech.ac.jp

1 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 a wide range of fields, e.g., earth observation, agriculture, and medical and biological imaging Chang (2003); Plaza et al. (2009); Rasti et al. (2018).

Observed HS images are often affected by noise because of the small amount of light in narrow wavelength and/or sensor failure. Also, in compressive HS imaging scenarios Willett et al. (2014); Arce et al. (2014), we have to estimate a full HS image from a very small 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 a-priori 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 Yuan et al. (2012); Chang et al. (2015); Aggarwal and Majumdar (2016); Liu et al. (2018) have been proposed for HS image restoration. Besides, LRM-based approaches 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) Zhang et al. (2014).

Many recent methods Liu et al. (2018); Li et al. (2018); He et al. (2016, 2018); Cao et al. (2018); Wang et al. (2017, 2018); Sun et al. (2018); Ince (2019); Kong et al. (2020) 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, and so the methods require to carefully control the hyperparameter(s) balancing these terms. Specifically, such hyperparameters are interdependent, which means that a suitable value of a hyperparameter varies depending both 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.

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 ability of noise and artifact removal 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) Gabay and Mercier (1976); Eckstein and Bertsekas (1992a, b); Boyd et al. (2011).

The remainder of the paper is organized as follows. Section II introduces notation and mathematical ingredients. Section III reviews existing methods related to our method. In Section IV, 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 V, where we illustrate the advantages of our method over several state-of-the-art methods. Section VI concludes the paper. The preliminary versions of this work, without mathematical details, deeper discussion, new applications, nor comprehensive experiments have appeared in conference proceedings Takeyama et al. (2017, 2019).

2 Preliminaries

2.1 Notation and Definitions

In this paper, let ${\mathbb{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 $(\cdot)^{\top}$, and the Euclidean norm (the $\ell_{2}$ norm) of a vector by $\|\cdot\|$.

For notational convenience, we treat an HS image ${\boldsymbol{\mathcal{U}}}\in{\mathbb{R}}^{N_{v}\times N_{h}\times B}$ as a vector ${\mathbf{u}}\in{\mathbb{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 $i$th pixel in $k$th band is $i+(k-1)N$ (for $i=1,\ldots,N$ and $k=1,\ldots,B$).

2.2 Proximal Tools

A function $f:{\mathbb{R}}^{L}\rightarrow(-\infty,\infty]$ is called proper lower semicontinuous convex if $\mbox{dom}(f):=\{{\mathbf{x}}\in{\mathbb{R}}^{L}|\;f({\mathbf{x}})<\infty\}\neq\emptyset$, $\mbox{lev}_{\leq\alpha}(f):=\{{\mathbf{x}}\in{\mathbb{R}}^{L}|\;f({\mathbf{x}})\leq\alpha\}$ is closed for every $\alpha\in{\mathbb{R}}$, and $f(\lambda{\mathbf{x}}+(1-\lambda){\mathbf{y}})\leq\lambda f({\mathbf{x}})+(1-\lambda)f({\mathbf{y}})$ for every ${\mathbf{x}},{\mathbf{y}}\in{\mathbb{R}}^{L}$ and $\lambda\in(0,1)$, respectively. Let $\Gamma_{0}({\mathbb{R}}^{L})$ be the set of all proper lower semicontinuous convex functions on ${\mathbb{R}}^{L}$.

The proximity operatorMoreau (1962) plays a central role in convex optimization based on proximal splitting. The proximity operator of $f\in\Gamma_{0}({\mathbb{R}}^{L})$ with an index $\gamma>0$ is defined by

[TABLE]

We introduce the indicator function of a nonempty closed convex set $C\subset{\mathbb{R}}^{L}$, which is defined as follows:

[TABLE]

Then, for any $\gamma>0$, its proximity operator is given by

[TABLE]

where $P_{C}({\mathbf{x}})$ is the metric projection onto $C$.

2.3 Alternating Direction Method of Multipliers (ADMM)

ADMM Gabay and Mercier (1976); Eckstein and Bertsekas (1992a, b); Boyd et al. (2011) is a popular proximal splitting method, and it can solve convex optimization problems of the form:

[TABLE]

where $f\in\Gamma_{0}({\mathbb{R}}^{L_{1}})$, $g\in\Gamma_{0}({\mathbb{R}}^{L_{2}})$, and ${\mathbf{G}}\in{\mathbb{R}}^{L_{2}\times 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 ${\mathbf{G}}$ is a full-column rank matrix. For arbitrarily chosen ${\mathbf{z}}^{(0)},{\mathbf{d}}^{(0)}$ and a step size $\gamma>0$, ADMM iterates the following steps:

[TABLE]

Convergence property of ADMM is given as follows.

Theorem 1** (Convergence of ADMM Eckstein and Bertsekas (1992b)).**

Consider Prob. (1), and assume that ${\mathbf{G}}^{\top}{\mathbf{G}}$ is invertible and that a saddle point of its unaugmented Lagrangian ${\mathcal{L}}_{0}({\mathbf{x}},{\mathbf{z}},{\mathbf{y}}):=f({\mathbf{x}})+g({\mathbf{z}})-\langle{\mathbf{d}},{\mathbf{G}}{\mathbf{x}}-{\mathbf{z}}\rangle$ exists.111 A triplet $(\hat{\mathbf{x}},\hat{\mathbf{z}},\hat{\mathbf{d}})$ is a saddle point of an unaugmented Lagrangian ${\mathcal{L}}_{0}$ if and only if ${\mathcal{L}}_{0}(\hat{\mathbf{x}},\hat{\mathbf{z}},{\mathbf{d}})\leq{\mathcal{L}}_{0}(\hat{\mathbf{x}},\hat{\mathbf{z}},\hat{\mathbf{d}})\leq{\mathcal{L}}_{0}({\mathbf{x}},{\mathbf{z}},\hat{\mathbf{d}})$, for any $({\mathbf{x}},{\mathbf{z}},{\mathbf{d}})\in{\mathbb{R}}^{L_{1}}\times{\mathbb{R}}^{L_{2}}\times{\mathbb{R}}^{L_{2}}$. Then the sequence $({\mathbf{x}}_{n})_{n>0}$ generated by (5) converges to an optimal solution to Prob. (1).

3 Related Works

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

3.1 TV-based Methods

The methods proposed in Yuan et al. (2012); Aggarwal and Majumdar (2016); Chang et al. (2015) restore a desirable HS image by solving a convex optimization problem involving TV-based regularization. Let $\bar{{\mathbf{u}}}\in{\mathbb{R}}^{NB}$ be the desirable HS image, and the authors assume that an observation ${\mathbf{v}}\in{\mathbb{R}}^{NB}$ is modeled as follows:

[TABLE]

where ${\mathbf{n}}$ and ${\mathbf{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:

[TABLE]

where $\mathcal{R}_{\mathrm{TV}}$ is a regularization function based on TV, and $\lambda_{1}$ and $\lambda_{2}$ are hyperparameters. Here, The first and third terms evaluate data-fidelity on Gaussian and sparse noise, respectively. The hyperparameters $\lambda_{1}$ and $\lambda_{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 a very important but troublesome task.

In the following, we explain each TV. Let ${\mathbf{D}}=({\mathbf{D}}_{v}^{\top}{\mathbf{D}}_{h}^{\top})^{\top}\in{\mathbb{R}}^{2NB\times NB}$ be spatial differences operator with ${\mathbf{D}}_{v}$ and ${\mathbf{D}}_{h}$ being vertical and horizontal differences operator, respectively, and spectral differences operator are ${\mathbf{D}}_{b}\in{\mathbb{R}}^{NB\times NB}$. In Yuan et al. (2012); Aggarwal and Majumdar (2016); Chang et al. (2015), HTV, ASSTV, and SSTV are defined as follows:

[TABLE]

where $\|\cdot\|_{\mathrm{TV}}$ is a TV norm, which takes the $\ell_{2}$ norm of spacial difference vectors for all band and then summing up for all spatial pixels, and $\tau_{v}$, $\tau_{h}$, and $\tau_{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 Bresson and Chan (2008). HTV does not consider spectral correlation, resulting in spatial oversmoothing. To consider spectral correlation, the authors of Yuan et al. (2012) 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 (Fig. 1, blue line). The weights $\tau_{v}$, $\tau_{h}$, and $\tau_{b}$ in (8) balance the smoothness related to 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 $\tau_{v}$, $\tau_{h}$, and $\tau_{b}$. SSTV evaluate a-prior knowledge on HS images using spatio-spectral piecewise-smoothness. It is derived by calculating spatial differences through spectral differences (Fig. 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.

3.2 LRM-based Method

LRMR Zhang et al. (2014) 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\times q\times B$, and LRMR handles it as a matrix of size $q^{2}\times B$ that is obtained by lexicographically arranging the spatial vectors in the patch cube in the row direction. The observation model is expressed like Sec. 3.1, and the restoration problem is formulated as follows:

[TABLE]

where ${\mathbf{U}}_{i,j}$, ${\mathbf{V}}_{i,j}$, and ${\mathbf{S}}_{i,j}$ represents the patches of a restored HS image, an observation, and a sparse noise, respectively, which are centered at (i,j) pixel. Then, The $\|\cdot\|_{F}$ is a Frobenius norm, $\operatorname{rank}(\cdot)$ represents a rank function, and $\operatorname{card}(\cdot)$ 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 ${\mathbf{U}}_{i,j}$ and the cardinality of ${\mathbf{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 non-convex, and so it is a troublesome task to seek the global optimal solution of Prob. (10).

3.3 Combined Method

The methods Liu et al. (2018); Li et al. (2018); He et al. (2016, 2018); Cao et al. (2018); Wang et al. (2017, 2018); Sun et al. (2018); Ince (2019); Kong et al. (2020) combine TV-based and LRM-based approaches. Since they can evaluate multiple types of a-priori knowledge, i.e., piecewise-smoothness and low rankness, they can restore a more desirable HS image than the approaches only using TV-based or LRM-based regularization. Besides, some methods Liu et al. (2018); Li et al. (2018); Cao et al. (2018); Wang et al. (2017, 2018); Kong et al. (2020) approximate the rank and cardinality functions by their convex surrogates. As a result, the restoration problems are convex and can be solved by optimization methods based on proximal splitting.

However, the methods have to handle multiple regularization terms and/or a data-fidelity term(s) simultaneously in one objective function, and so 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.

4 Proposed Method

4.1 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

[TABLE]

where $\|\cdot\|_{1,p}$ is the mixed $\ell_{1,p}$ norm, and $\omega\geq 0$. We assume $p=1$ or $2$, i.e., the $\ell_{1}$ norm ($\|\cdot\|_{1,1}=\|\cdot\|_{1}$) or the mixed $\ell_{1,2}$ norm, respectively. We would like to remark that we can also see $\ell_{1}$-HSSTV ($p=1$) as anisotropic HSSTV and $\ell_{1,2}$-HSSTV ($p=2$) as isotropic HSSTV.

In (11), ${\mathbf{D}}{\mathbf{D}}_{b}{\mathbf{u}}$ and ${\mathbf{D}}{\mathbf{u}}$ correspond to local spatio-spectral and direct local spatial differences, respectively, as shown in Fig. 1 (red lines). The weight $\omega$ adjusts the relative importance of direct spatial piecewise-smoothness to spatio-spectral piecewise-smoothness. HSSTV evaluates two kinds of smoothness by taking the $\ell_{p}$ norm ($p=1$ or $2$) of these differences associated with each pixel and then summing up for all pixels, i.e., calculating the $\ell_{1}$ norm. Thus, it can be defined via the mixed $\ell_{1,p}$ norm. When we set $\omega=0$ and $p=1$, HSSTV recovers SSTV as (9), meaning that HSSTV can be seen as a generalization of SSTV.

As reviewed in Sec. 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. Fig. 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 $\omega$ is large. Thus, the weight $\omega$ should be set to less than one, as will be demonstrated in Sec. 5.

4.2 HS Image Restoration by HSSTV

We consider to restore a desirable HS image $\bar{{\mathbf{u}}}\in{\mathbb{R}}^{NB}$ from an observation ${\mathbf{v}}\in{\mathbb{R}}^{M}~{}(M\leq NB)$ contaminated by a Gaussian-sparse mixed noise. The observation model is given by the following form:

[TABLE]

where ${\mathbf{\Phi}}\in{\mathbb{R}}^{M\times NB}$ is a matrix representing a linear observation process, e.g., random sampling, ${\mathbf{n}}\in{\mathbb{R}}^{M}$ is a Gaussian noise with the standard deviation $\sigma$, and ${\mathbf{s}}\in{\mathbb{R}}^{M}$ is a sparse noise.

Based on the above model, we formulate HS image restoration using HSSTV as the following optimization problem:

[TABLE]

where ${\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}}$ is a ${\mathbf{v}}$-centered $\ell_{2}$-norm ball with the radius $\varepsilon>0$, ${\mathcal{B}}_{1,\eta}$ is a $\mathbf{0}$-centered $\ell_{1}$-norm ball with the radius $\eta>0$, and $[\mu_{\min},\mu_{\max}]^{NB}$ is a dynamic range of an HS image ($\mu_{\min}<\mu_{\max}$). This method simultaneously estimates the desirable HS image ${\mathbf{u}}$ and the sparse noise ${\mathbf{s}}$ for noise-robust restoration. The first and second constraints measure data fidelities to the observation ${\mathbf{v}}$ and the sparse noise ${\mathbf{s}}$, respectively. As mentioned in Afonso et al. (2011); Chierchia et al. (2015); Ono and Yamada (2015); Xie et al. (2016); He et al. (2016, 2018); Ono (2017); Takeyama et al. (2017, 2018a, 2018b); Chan et al. (2011); Takeyama et al. (2019), such a constraint-type data-fidelity enables us to translate the hyperparameter(s) balancing between regularization and data-fidelity like $\lambda_{1}$ and $\lambda_{2}$ in (6) to the upper bound of the degree of data-fidelity $\varepsilon$ and $\eta$ that can be set in a much easier manner.

Since all constraints are closed convex sets, and HSSTV is a convex function, Prob. (4.2) is a constrained convex optimization problem. In this paper, we adopt ADMM (see Sec. 2.3) for solving the problem. In what follows, we reformulate Prob. (4.2) into Prob. (1).

By using the indicator functions of the constraints, Prob. (4.2) can be rewritten as

[TABLE]

Note that from the definition of the indicator function, Prob. (17) exactly equals to Prob. (4.2). By letting

[TABLE]

Prob. (17) is reduced to Prob. (1). The resulting algorithm based on ADMM is summarized in Alg. 1.

The update of ${\mathbf{u}}$ and ${\mathbf{s}}$ in Alg. 1 come down to the following forms:

[TABLE]

Since the update of ${\mathbf{u}}$ and ${\mathbf{s}}$ in Alg. 1 is strictly-convex quadratic minimization, one can obtain this update forms by differentiating it. Here, we should consider the structure of ${\mathbf{\Phi}}$ because it affects the matrix inversion in (21). If ${\mathbf{\Phi}}$ is a block-circulant-with-circulant-blocks (BCCB) matrix Hansen et al. (2006), we can leverage 3DFFT to efficiently solve the inversion in Step 2 with the difference operators having periodic boundary, i.e., ${\mathbf{A}}_{\omega}^{\top}{\mathbf{A}}_{\omega}+{\mathbf{\Phi}}^{\top}{\mathbf{\Phi}}+{\mathbf{I}}$ can be diagonalized by the 3D FFT matrix and its inverse. If ${\mathbf{\Phi}}$ is a semi-orthogonal matrix, i.e., ${\mathbf{\Phi}}{\mathbf{\Phi}}^{\top}=\alpha{\mathbf{I}}$ $(\alpha>0)$, we leave it to the update of ${\mathbf{z}}_{2}$, which means that we replace $\iota_{{\mathcal{B}}^{{\mathbf{v}}}_{2,\varepsilon}}$ by $\iota_{{\mathcal{B}}^{{\mathbf{v}}}_{2,\varepsilon}}\circ{\mathbf{\Phi}}$ in (19) and ${\mathbf{\Phi}}{\mathbf{u}}$ by ${\mathbf{u}}$ in (20). This is because the proximity operator of $\iota_{{\mathcal{B}}^{{\mathbf{v}}}_{2,\varepsilon}}\circ{\mathbf{\Phi}}$ in this case can be computed by using (Combettes and Pesquet, 2011, Table 1.1-x) as follows:

[TABLE]

If ${\mathbf{\Phi}}$ is a sparse matrix, we offer to use a preconditioned conjugate gradient method Golub and Loan (2012) for approximately solving the inversion, or to apply primal-dual splitting methods Chambolle and Pock (2010); Combettes and Pesquet (2012); Condat (2013) instead of ADMM.222Primal-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 Ono et al. (2016); Chambolle et al. (2018); Combettes and Pesquet (2016); Ono (2019) might be useful for reducing the computational cost.

For the update of ${\mathbf{z}}_{1}$, the proximity operators are reduced to simple soft-thresholding type operations: for $\gamma>0$ and for $i=1,\ldots,4NB$, (i) in the case of $p=1$,

[TABLE]

where $\operatorname{sgn}$ is the sign function, and (ii) in the case of $p=2$,

[TABLE]

where $\tilde{i}:=((i-1)\mod NB)+1$.

The update of ${\mathbf{z}}_{2}$, ${\mathbf{z}}_{3}$, and ${\mathbf{z}}_{4}$ require the proximity operators of the indicator functions of ${\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}}$, ${\mathcal{B}}_{1,\eta}$ and $[\mu_{\min},\mu_{\max}]^{NB}$, respectively, which equal to the metric projections onto them (see Sec 2.2). Specifically, the metric projection onto ${\mathcal{B}}_{2,\varepsilon}^{{\mathbf{v}}}$ is given by

[TABLE]

that onto ${\mathcal{B}}_{1,\eta}$ is given by

[TABLE]

and that onto $[\mu_{\min},\mu_{\max}]^{NB}$ is given, for $i=1,\ldots,NB$, by

[TABLE]

5 Experiments

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 Spe , MultiSpec Mul and GIC GIC , where their dynamic range were normalized into $[0,1]$.

The proposed method was compared with HTV Yuan et al. (2012), SSAHTV Yuan et al. (2012), SSTV Aggarwal and Majumdar (2016), and ASSTV Chan et al. (2011). For a fair comparison, we replaced HSSTV in Prob. (4.2) with HTV, SSAHTV, SSTV, or ASSTV and solved the problem by ADMM. In the denoising experiments, we also compared our proposed method with LRMR Zhang et al. (2014) and TV-regularized low-rank matrix factorization (LRTV) He et al. (2016). Since LRMR and LRTV are customized to the mixed noise removal problem, we cannot adopt them for CS reconstruction. We did not compare our proposed method with a recent CNN-based HS image denoising method Liu and Lee (2019). The CNN-based method cannot be represented as explicit regularization functions and is 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 quantitively evaluate restoration performance, we used the peak signal-to-noise ratio (PSNR) [dB] index and the structural similarity (SSIM) Wang et al. (2004) index between a true HS image $\bar{\mathbf{u}}$ and a restored HS image ${\mathbf{u}}$. PSNR is defined by $10\log_{10}(NB/\|{\mathbf{u}}-\bar{\mathbf{u}}\|^{2})$, and the higher the value 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:

[TABLE]

where ${\mathbf{u}}_{i}$ and $\bar{{\mathbf{u}}}_{i}$ are the $i$th pixel-centered local patches of a restored HS image and a true HS image, respectively, $P$ is the number of patches, $\mu_{{\mathbf{u}}_{i}}$ and $\mu_{\bar{{\mathbf{u}}}_{i}}$ is the average values of the local patches of the restored and true HS images, respectively, $\sigma_{{\mathbf{u}}_{i}}$ and $\sigma_{\bar{{\mathbf{u}}}_{i}}$ represent the variances of ${\mathbf{u}}_{i}$ and $\bar{{\mathbf{u}}}_{i}$, respectively, and $\sigma_{{\mathbf{u}}_{i}\bar{{\mathbf{u}}}_{i}}$ denotes the covariance between ${\mathbf{u}}_{i}$ and $\bar{{\mathbf{u}}}_{i}$. Moreover, $C_{1}$ and $C_{2}$ are two constants, which avoid the numerical instability when either $\mu_{{\mathbf{u}}_{i}}^{2}+\mu_{\bar{{\mathbf{u}}}_{i}}^{2}$ or $\sigma_{{\mathbf{u}}_{i}}^{2}+\sigma_{\bar{{\mathbf{u}}}_{i}}^{2}$ is very close to zero. SSIM gives a normalized score between zero and one, where the maximum value means that ${\mathbf{u}}$ equals to $\bar{{\mathbf{u}}}$.

We set the max iteration number, the stepsize $\gamma$ and the stopping criterion of ADMM to 10000, 0.05 and $\|{\mathbf{u}}^{(n)}-{\mathbf{u}}^{(n+1)}\|<0.01$, respectively.

5.1 Denoising

First, we experimented on Gaussian-sparse mixed noise removal of HS images, where observed HS images included an additive white Gaussian noise ${\mathbf{n}}$ with the standard deviation $\sigma$ and sparse noise ${\mathbf{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 is $s_{p}$, $l_{v}$, and $l_{h}$, respectively. We generated noisy HS images by adding two types of mixed noise to true HS images: (i) $(\sigma,~{}s_{p},~{}l_{v}=l_{h})=(0.05,~{}0.04,~{}0.04)$, (ii) $(0.1,~{}0.05,~{}0.05)$. In the denoising case, ${\mathbf{\Phi}}={\mathbf{I}}$ in (12), and the radiuses $\varepsilon$ and $\eta$ in Prob. (4.2) were set to $0.83\sqrt{NB(1-(s_{p}(1-l_{v}-l_{h})+l_{v}+l_{h}-l_{v}l_{h}))\sigma^{2}}$ and $NB(0.45s_{p}+(l_{v}+l_{h})v_{ave}-l_{v}l_{h}v_{ave})$, respectively, where $v_{ave}$ is the average of the observed image. Table 2 shows the parameters settings for ASSTV, LRMR, LRTV, and the proposed method. We set these parameters to achieve the best performance for each method.

In Tab. 3, we show 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. In the LRTV case, one can see that some results by LRTV outperform them by the proposed method. LRTV utilizes both TV-based and LRM-based regularization techniques, leading to higher-quality restoration than the proposed method. Meanwhile, even though the proposed method uses only TV-based regularization, it outperforms LRTV over half situations.

Fig. 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 spacial details, and noise remains in the results by SSTV and LRMR. 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.

Fig. 4 plots PSNR or SSIM of the results by the proposed method versus various $\omega$ in (11) changed from $0.01$ to $0.2$, where the values of PSNR and SSIM are averaged over the 13 HS images. One can see that $\omega\in[0.03,~{}0.07]$ is a good choice. ASSTV and LRTV require to adjust the weight $\tau_{b}$ and the hyperparameter $\tau$ newly for difference noise intensity, respectively, but the suitable parameter $\omega$ in HSSTV is noise-robust.

Fig. 5 plots bandwise PSNR and SSIM (left) and spatial and spectral responses (right) of the denoised Suwannee HS image in the case of the noise level (ii). The graphs regarding bandwise 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 in the SSIM cases. The graph (c) plots the spatial response of the 243rd row of the 30th band. In the same way, the 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 true one. On the other hand, there exist undesirable variations in the spatial response of the result by SSTV and LRMR. In contrast, ASSTV, LRTV, and the proposed method restore similar responses to the true one. In the 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 the that of the true one, but the mean value is larger than the true one, (iii) the spectral response of the results by ASSTV is too smooth and different from the true one, and (iv) LRTV and the proposed method can restore a spectral response very similar to the true one.

To verify the sensitivity of the parameter $\varepsilon$ and $\eta$, we conducted additional experiments, where we examined various values of $\varepsilon$ and $\eta$. Specifically, we set $\varepsilon=0.83\alpha\sqrt{NB(1-(s_{p}(1-l_{v}-l_{h})+l_{v}+l_{h}-l_{v}l_{h}))\sigma^{2}}$ and $\eta=\beta 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 $\alpha$ and $\beta$ from $0.9$ to $1.1$ at $0.02$ interval (the DC and the KSC images and the noise level (ii)). Fig. 6 plots PSNR or SSIM of the results by HTV, SSAHTV, SSTV, ASSTV, and the proposed method versus $\alpha$ or $\beta$. For HTV, SSAHTV, ASSTV, and the proposed method, the graphs show that the suitable values of $\alpha$ and $\beta$ do not vary significantly for both image, and so the parameters $\varepsilon$ and $\eta$ 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 $\alpha,\beta>1$.

5.2 Real Noise Removal

We also examined HTV, SSAHTV, SSTV, ASSTV, LRMR, LRTV, 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 ${\mathbf{v}}$. To maximize the performance of each method, we searched for suitable values of $\sigma$, $s_{p}$, $l_{v}$, and $l_{h}$, and we set the parameters as with Sec. V. A. Specifically, we set the parameter $\varepsilon=3.1893$ and $\eta=31893$ for all TVs. Besides, the parameters $\tau_{v}$, $\tau_{h}$, and $\tau_{b}$ in ASSTV are set as $1$, $1$, $3$, respectively, $r$ and $k$ in LRMR are set as $3$ and $0.1204$, respectively, and $r$ and $tau$ in LRTV are set as $2$ and $0.008$.

Fig 7 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, and LRMR 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.

5.3 Compressed Sensing Reconstruction

We conducted on an experiment on compressed sensing (CS) reconstruction Baraniuk (2007); Candès and Wakin (2008). 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 Willett et al. (2014); Arce et al. (2014).

In this experiment, we assume that ${\mathbf{\Phi}}\in{\mathbb{R}}^{M\times NB}$ in (12) is a random sampling matrix with the sampling rate $m$ ($0<m<1$ and $M=mNB$). Here, since ${\mathbf{\Phi}}$ is a semi-orthogonal matrix, we can efficiently solve the problem as explained in Sec. 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 ${\mathbf{n}}$ with noise intensity $\sigma=0.1$.

We set the CS reconstruction problem as follows:

[TABLE]

The problem is derived by removing the second constraint and ${\mathbf{s}}$ from Prob. (4.2). Therefore, we can solve the above problem by removing ${\mathbf{s}}$, ${\mathbf{z}}_{3}$, and ${\mathbf{d}}_{3}$ in Alg. (1) and replacing ${\mathbf{z}}_{4}$ and ${\mathbf{d}}_{4}$ with ${\mathbf{z}}_{3}$ and ${\mathbf{d}}_{3}$, respectively. As in Sec. 4.2, the update of ${\mathbf{u}}$ is strictly-convex quadratic minimization, and so it comes down to

[TABLE]

We set $m=0.2$ or $0.4$ and $\varepsilon=\sqrt{mNB\sigma^{2}}$. In the ASSTV case, we set the parameters $(\tau_{v},\tau_{h},\tau_{b})=(1,1,0.5)$, which experimentally achieves the best performance.

Tab. 4 shows PSNR and SSIM of the reconstructed HS images. For all $m$ and HS images, both PSNR and SSIM of the results by the proposed method are almost higher than that by HTV, SSAHTV, SSTV, and ASSTV.

Fig. 8 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, where it appears as the difference from the color of the true 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.

Fig. 9 plots PSNR or SSIM of the results by the proposed method versus $\omega$ averaged over the 13 HS images for each $m$. The graphs show that $\omega\in[0.05,0.1]$ is a good choice in most cases. In comparison with Fig. (4), the suitable range of $\omega$ in CS reconstruction is almost the same as that in denoising.

Fig. 10 plots bandwise PSNR or SSIM (left) and spatial and spectral responses (right) (Suwannee, $m=0.2$). According to bandwise PSNR and SSIM, one can see that the proposed method achieves higher-quality reconstruction for all bands than HTV, SSAHTV, SSTV, and ASSTV. The graphs (c) and (d) plot the spatial and spectral responses of the same position in Sec. 5.1. The 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 true one. In the 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 true one, but the mean values are larger than the true one. The proposed method achieves the most similar reconstruction of spectral response among all the TVs.

6 Conclusion

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 constraints. 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.

\authorcontributions

Conceptualization, S.T., S.O., and I.K.; methodology, S.T. and S.O.; software, S.T.; validation, S.T.; formal analysis, S.T.; investigation, S.T.; writing—original draft, S.T.; writing—review and editing, S.O. and I.K.; supervision, S.O., and I.K.; project administration, S.T., S.O., and I.K.; funding acquisition, S.T., S.O., and I.K. All authors read and agreed to the published version of the manuscript.

\funding

This work was supported in part by JST CREST under Grant JPMJCR1662 and JPMJCR1666, and in part by JSPS KAKENHI under Grant 18J20290, 18H05413, and 20H02145.

\conflictsofinterest

The authors declare no conflict of interest.

\reftitle

References

\externalbibliography

yes