論文 / 著書情報
Article / Book Information
題目(和文)

制約付き凸最適化による ロバストな画像復元・合成に関する研究

Title(English)

A study of robust image restoration and fusion via constrained convex
optimization

著者(和文)

武山彩織

Author(English)

Saori Takeyama

出典(和文)

学位:博士（工学）,学位授与機関:東京工業大学,
報告番号:甲第11939号,
授与年月日:2021年3月26日,
学位の種別:課程博士,
審査員:熊澤 逸夫,小野 峻佑,山口 雅浩,山田 功,篠﨑 隆宏

Citation(English)

Degree:Doctor (Engineering),
Conferring organization: Tokyo Institute of Technology,
Report number:甲第11939号,
Conferred date:2021/3/26,
Degree Type:Course doctor,
Examiner:,,,,

学位種別(和文)

博士論文

Type(English)

Doctoral Thesis

Powered by T2R2 (Tokyo Institute Research Repository)

Tokyo Institute of Technology
School of Engineering
Department of Information and Communications Engineering
PhD thesis

A study of robust image restoration andfusion via constrained convex optimization
Supervisor Professor Itsuo Kumazawa
Co-supervisor Associate Professor Shunsuke Ono
January, 2021.
Department Information and Communications Engineering
Name
Saori Takeyama

Contents
1 Introduction
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Organization of this paper . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2

2 Preliminaries
2.1 Proximal Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 ADMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Primal-Dual Splitting Method . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
5
6

3 Image Restoration
3.1 Hyperspectral Image Denoising . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Related Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
TV-Based Methods . . . . . . . . . . . . . . . . . . . . . . . . . . .
LRM-Based Method . . . . . . . . . . . . . . . . . . . . . . . . . .
Combined Method . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Proposed method . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Hybrid Spatio-Spectral Total Variation . . . . . . . . . . . . . . .
HS Image Denoising by HSSTV . . . . . . . . . . . . . . . . . . . .
3.1.3 Experimets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gaussian noise removal . . . . . . . . . . . . . . . . . . . . . . . .
Mixed noise removal . . . . . . . . . . . . . . . . . . . . . . . . . .
Real Noise Removal . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Compressed Sensing Reconstruction . . . . . . . . . . . . . . . . . . . . .
3.2.1 CS Reconstruction under Gaussian Noise Situation . . . . . . . . .
Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . .
Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 CS Reconstruction under Mixed Noise Situation . . . . . . . . . .
Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . .
Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
8
8
9
11
11
11
12
14
14
23
28
28
29
29
30
31
35
35
36
37

i

4 Image Fusion
4.1 Blur and Noisy Image Fusion . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Basic Performance Evaluation . . . . . . . . . . . . . . . . . . . . .
Robustness to Inaccurate Blur Kernels . . . . . . . . . . . . . . . .
Facilitation of Parameter Setting . . . . . . . . . . . . . . . . . . .
4.2 Hyperspectral Image Fusion for Super-resolution . . . . . . . . . . . . . .
4.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Hyperspectral Pansharpening Experiments . . . . . . . . . . . . .
4.2.4 Hyperspectral and Multispectral Image Fusion Experiments . . . .
4.2.5 Compressed Hyperspecteral Pansharpening . . . . . . . . . . . . .
Proposed Framework . . . . . . . . . . . . . . . . . . . . . . . . . .
Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39
39
40
41
43
44
46
47
49
49
50
51
56
59
59
60

5 Conclusion

64

6 Acknowledgement

66

References

67

Publications related to this dissertation

78

ii

List of Figures
3.1

Calculation of local diﬀerences in SSTV, ASSTV, and our HSSTV. SSTV
evaluates the ℓ1 norm of spatio-spectral diﬀerences (yellow line). ASSTV
evaluates the ℓ1 norm of direct spatial and spectral diﬀerences (blue line).
HSSTV evaluates the mixed ℓ1,p norm of direct spatial and spatio-spectral
diﬀerences (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Restored HS images from an observation contaminated by similar noise inadjacent bands (the upper half area) and random noise (the lower half). .
3.3 Resulting HS images with their PSNR[dB] (left) and SSIM (right) in theGaussian denoising experiment (top: Suwannee, σ = 0.1, bottom: DC,
σ = 0.2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 PSNR[dB] or SSIM versus ω in (3.6) in the Gaussian denoising experiment.
3.5 Bandwise PSNR[dB] and SSIM and spatial and spectral responses in theGaussian denoising experiment (Suwannee). . . . . . . . . . . . . . . . . .
3.6 Resulting HS images with their PSNR[dB] (left) and SSIM (right) in themixed noise removal experiment (top: Salinas, the noise level (i), bottom:
PaviaU, the noise level (ii)). . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 Band-wise PSNR[dB] and SSIM and spatial and spectral responses in themixed noise removal experiment (Suwannee). . . . . . . . . . . . . . . . .
3.8 PSNR[dB] (top) or SSIM (bottom) versus ω in (3.6) in the mixed noise
removal experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 PSNR[dB] or SSIM versus α or β on the mixed noise removal experiment(the noise level (ii), left: DC, right: KSC). . . . . . . . . . . . . . . . . . .
3.10 The denoising results on the real noise removal experiments. . . . . . . . .
3.11 Resulting HS images with their PSNR[dB] (left) and SSIM (right) on the
CS reconstruction experiment under Gaussian noise situation (top: KSC,m = 0.4, bottom: Reno, m = 0.2). . . . . . . . . . . . . . . . . . . . . . .
3.12 Band-wise PSNR[dB] and SSIM and spatial and spectral responses on the
CS reconstruction experiment under Gaussian noise situation (Suwannee).
3.13 PSNR[dB] (top) or SSIM (bottom) versus ω in (3.6) on the CS reconstruction experiment under Gaussian noise situation. . . . . . . . . . . . . . . .
3.14 Resulting HS images on CS reconstruction under mixed noise situation(Reno, top: ℓ2 fidelity constraint only, bottom: proposed). . . . . . . . . .
4.1

Test images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iii

10
12

17
18
19

21
22
25
26
28

33
34
35
37
43

4.2
4.3
4.4

PSNR gain of the proposed method over the single image deblurring/denoising. 44
Resulting images with their PSNR[dB] in the first experiment (σb /σn = 8). 45
Evolution of NRMSEn versus iterations (left) and the evolution of PSNR[dB]
versus iterations (right) of Alg. 4. . . . . . . . . . . . . . . . . . . . . . . . 46
4.5 Restored results with their PSNR[dB] by using inaccurate blur kernels in
the motion blur case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.6 Restored results with their PSNR[dB] by using an inaccurate blur kernelin the spatially-varying blur case. . . . . . . . . . . . . . . . . . . . . . . . 48
4.7 Quality measures versus λ in (4.10) (top) / ω in (3.6) (bottom). . . . . . . 53
4.8 Resulting HS images on the HS pansharpening experiments (σH = 0.1,
σg = 0.05). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.9 Resulting HS images with four quality measure (PSNR[dB], SAM, ERGAS,
Q2n ) (DC, σg = 0.04). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.10 Resulting HS images (m = 0.0625 (CHPAN case) and r = 16 (HS pansharpening case)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.11 Spectral response of the results. . . . . . . . . . . . . . . . . . . . . . . . . 62

iv

List of Tables
2.1

Correspondence table between the abbreviations . . . . . . . . . . . . . .

4

3.1
3.2
3.3
3.4
3.5
3.6

The feature of existing methods for HS image restoration. . . . . . . . . .
PSNR[dB] of the results on Gaussian noise removal experiments. . . . . .
SSIM of the results on Gaussian noise removal experiments. . . . . . . . .
Parameter settings for ASSTV, LRMR, LRTV, and the proposed method.
PSNR[dB] (top) and SSIM (bottom) in mixed noise removal experiments.
PSNR[dB] (top) and SSIM (bottom) in the CS reconstruction experiment
under Gaussian noise situation. . . . . . . . . . . . . . . . . . . . . . . . .
The four quality measures of the results on CS reconstruction under mixed
noise situation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
15
16
23
24

Results of the experiment using inaccurate blur kernel. . . . . . . . . . . .
Best values of λ1 and λ2 in the existing formulation (Prob. (4.7)). . . . . .
Quality measures for σg = 0.025 (top), σg = 0.05 (middle), and σ = 0.075(bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Quality measures for σg = 0.04 (left) and σg = 0.06 (right) averaged all HS
images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Performance comparison with four quality measures (left: m = 0.25 (CHPAN case) and r = 4 (HS pansharpening case); right: m = 0.0625 (CHPAN
case) and r = 16 (HS pansharpening case)). . . . . . . . . . . . . . . . . .

46
48

3.7

4.1
4.2
4.3
4.4
4.5

v

32
37

52
56

60

Chapter 1

Introduction
1.1

Background

Images are aﬀected by some degradations and noises through an observation process,e.g., blur, down-sampling, and missing pixel. The noises are caused by an electric current, heat, imaging device, and so on, and capturing a noiseless image would be impossible. Therefore, image restoration is a key technique and actively studied such as denoising, deblurring, super-resolution, and compressed sensing (CS) reconstruction. Since
the restoration problems assume various degradations and noises, they require a robust
image restoration method.
Most image restoration methods are established based on optimization: a desirable
image is characterized as a solution to some optimization problem, which consists of aregularization term and a data-fidelity term. The regularization term evaluates apriori
knowledge about underlying properties on images, and the data-fidelity term keeps theconsistency with a given observation. Thanks to the design, these methods get a reasonable result under ill-posed or ill-conditioned scenarios typical in image restoration.
However, since the methods are not enough robust for degradations or noises, the results have over-smoothing or artifacts. This is caused by (i) not enough utilizing apriori
knowledge on an image, (ii) not full consideration of outlier and modeling error.
Recently, high-resolution images, especially spatial one, is required by many application. The resolution of images depends on a sensor in an imaging device, and they have a
trade-oﬀ between spatial and spectral resolution. Therefore, estimating higher-resolution
images from observation is an essential task. Besides, since blur removes edge and detail
information in an image, it is limited to restore an edge-preserving image from an only
blur image. For these degradations, one eﬀective approach is using a guide image thathas a diﬀerent characteristic. This technique is named image fusion and helps to robust
estimation.
Based on the above discussion, we propose novel robust image restoration techniques
via constrained convex optimization. Specifically, we propose a new regularization thateﬀectively utilizes apriori knowledge on an image. The proposed regularization function
evaluates two types of piecewise-smoothness on images, so it helps high-quality restoration. Then, for some image restoration and fusion problems, we formulate each con-1

strained convex optimization problem considering the sparsity of outlier. The problem
evaluates data-fidelity between an observation and a restored image by hard constraints,
so one can set the hyperparameters with comparative ease. Finally, we solve the restoration problem by iterative proximal splitting methods. In the experiments, we demonstrate
the performance of the proposed methods compared with existing methods for four image restoration and fusion problems: denoising, CS reconstruction, blur and noisy image
fusion, and super-resolution.

1.2

Organization of this paper

The remainder of the paper is organized as follows.
Chapter 1 research background and the aim
Chapter 2 notation and mathematical ingredients
Chapter 3 robust image restoration methods for HS image denoising and CS reconstruction
Chapter 4 robust image fusion methods for blur and noisy image fusion and superresolution
Chapter 5 conclusion of the paper

2Chapter 2

Preliminaries
In this paper, let R be the set of real numbers. We shall use boldface lowercase andcapital 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 image U ∈ RNv ×Nh ×B as
a vector u ∈ RN B (N := Nv Nh is the number of the pixels of each band, and B is thenumber of the bands) by stacking its columns on top of one another, i.e., the index of thecomponent of the ith pixel in kth band is i + (k − 1)N (for i = 1, . . . , N and k = 1, . . . , B).
In Tab. 2.1, we show the correspondence between the abbreviations used in this paper.

2.1

Proximal Tools

A function f : RL → (−∞, ∞] is called proper lower semicontinuous convex if dom(f ) :={x ∈ RL | f (x) < ∞} ̸= ∅, lev≤α (f ) := {x ∈ RL | f (x) ≤ α} is closed for every α ∈ R, andf (λx + (1 − λ)y) ≤ λf (x) + (1 − λ)f (y) for every x, y ∈ RL and λ ∈ (0, 1), respectively.
Let Γ0 (RL ) be the set of all proper lower semicontinuous convex functions on RL .
The proximity operator [16] plays a central role in convex optimization based on proximal splitting. The proximity operator of f ∈ Γ0 (RL ) with an index γ > 0 is definedby1
proxγf (x) := argmin f (y) +∥y − x∥2 .
2γ
yWe introduce the indicator function of a nonempty closed convex set C ⊂ RL , which
is defined as follows:{
0, if x ∈ C,
ιC (x) :=∞, otherwise.
Then, for any γ > 0, its proximity operator is given byproxγιC (x) = PC (x) := argmin ∥x − y∥,
y∈C

3

(2.1)

Table 2.1: Correspondence table between the abbreviations
abbreviations
ADMM
ASSTV
BCCB
CC
CHPAN
CNMF
CS
ERGAS
GFPCA
GS
GSA
HR-HS image
HS
HSSTV
HTV
HySureLLRGTV
LNWTV
LR-HS image
LRM
LRMR
LRTV
MS
MTF-GLP
MTF-GLP-HPM
NRMSE
PAN
PCA
PSNR
RMSE
SAM
SFIM
SSAHTV
SSIM
SSTV
TV

formal namealternating direction method of multipliers
anisotropic spatio-spectral total variation
block-circulant-with-circulant-blocks
cross correlation
compressed hyperspectral pansharpening
coupled non-negative matrix factorization
compressed sensing
erreur relative globale adimensionnelle de synthèseguided filter principal component analysis
gram schmidt
gram schmidt adoptive
hyperspectral image of high spatial and spectral resolution
hyperspectral
hybrid spatio-spectral total variation
hyperspectral total variation
hyperspectral Stein’s unbiased risk estimator
local low-rank matrix recovery and global spatial-spectral
total variation
local spatial neighborhood weighted spectral-spatial
hyperspectral total variation
hyperspectral image of high spectral resolution
but low spatial resolution
low rank modeling
low-rank matrix recovery
total variation-regularized low-rank matrix factorization
multispectral
modulation transfer function-generalized Laplacian pyramid
modulation transfer function-generalized Laplacian pyramidwith high pass modulation
normalized root mean squared error
panchromatic
principal component analysis
peak signal-to-noise ratio
the root mean squared error
the spectral angle mappersmoothing filter-based intensity modulation
spectral-spatial adaptive hyperspectral total variation
signal similarity
spatio-spectral total variation
total variation4

where PC (x) is the metric projection onto C.
In this paper, we use proximity operators of multiple functions, and we introduce them
here. The proximity operators of the ℓ1 norm and the mixed ℓ1,2 norm are reduced tosimple soft-thresholding type operations: for γ > 0 and for i = 1, . . . , 4N B, (i) in the caseof the ℓ1 norm,[proxγ∥·∥1 (x)]i = sgn(xi ) max {|xi | − γ, 0} ,(2.2)
where sgn is the sign function, and (ii) in the case of the mixed ℓ1,2 norm,− 1 
2


3


∑
2
[proxγ∥·∥1,2 (x)]i = max 1 − γ xĩ+jN B  , 0 xi ,




j=0

(2.3)

where ĩ := ((i − 1) mod N B) + 1.
v := {x ∈ RN B |∥x − v∥ ≤ ε},
The proximity operators of the indicator functions of B2,εv is the v-centered ℓ -ball with
B1,η := {x ∈ RN B |∥x∥1 ≤ η}, and [µ, µ]N B . Here, B2,ε2
the radius ε, and B1,η is the 0-centered ℓ1 -ball with the radius η. In (2.1), the proximity
operator of the indicator functions equate the metric projections onto them. Specifically,v is given by
the metric projection to B2,ε{v ,x,
if x ∈ B2,εv (x) =
PB2,ε(2.4)
ε(x−v)
v + ∥x−v∥ , otherwise,that onto B1,η is given by
PB1,η = sgn(x) max(|x| − η, 0),(2.5)

and that onto [µ, µ]N B is given, for i = 1, . . . , N B, by[P[µ,µ]N B (x)]i = min{max{xi , µ}, µ}.

2.2

(2.6)

ADMM

ADMM [1–4] is a popular proximal splitting method, and it can solve convex optimization problems of the form:
min f (x) + g(z) s.t. z = Gx,
x,z

(2.7)

where f ∈ Γ0 (RL1 ), g ∈ Γ0 (RL2 ), and G ∈ RL2 ×L1 . Here, we assume that f is quadratic,
g is proximable, i.e., the proximity operator of g is computable in an eﬃcient manner, andG is a full-column rank matrix. For arbitrarily chosen z(0) , d(0) and a step size γ > 0,
ADMM iterates the following steps:
 (n+1)1
 x= argmin f (x) + 2γ
∥z(n) − Gx − d(n) ∥2 ,x
 (n+1)(2.8)
 z= proxγg (Gx(n+1) + d(n) ),
(n+1)
(n)
(n+1)
(n+1)
d
= d + Gx−z,
The convergence property of ADMM is given as follows.

5

Theorem 1 (Convergence of ADMM [3]). Consider Prob. (2.7), and assume that G⊤ G
is invertible and that a saddle point of its unaugmented Lagrangian L0 (x, z, y) := f (x) +g(z)−⟨d, Gx−z⟩ exists. A triplet (x̂, ẑ, d̂) is a saddle point of an unaugmented Lagrangian
L0 if and only if L0 (x̂, ẑ, d) ≤ L0 (x̂, ẑ, d̂) ≤ L0 (x, z, d̂), for any (x, z, d) ∈ RL1 × RL2 ×
RL2 . Then the sequence (xn )n>0 generated by (2.8) converges to an optimal solution toProb. (2.7).

2.3

Primal-Dual Splitting Method

A primal-dual splitting method [17] is a popular splitting method like ADMM, and it
can solve convex optimization problems of the form:
min g(x) + h(Lx)
x

(2.9)

where, g ∈ Γ0 (RL1 ) and h ∈ Γ0 (RL2 are proximable functions, and L ∈ RL2 ×L1 is a linear
operator. For arbitrarily chosen y(0) and step sizes γ1 , γ2 > 0, the primal-dual splitting
method iterates the following steps:⌊ (n+1)x:= proxγ1 g (x(n) − γ1 (x(n) + L⊤ y(n) )),(2.10)(n+1)
y:= proxγ2 h∗ (y(n) + γ2 L(2x(n+1) − x(n) )).
Here, h∗ is the convex conjugate function of h, and the proximity operator can be computed, in [18, Theorem 14.3 (ii)] as follows:( )xproxγh∗ (x) = x − γ prox 1 h.
(2.11)γγ
For the convergence of (2.10), the step sizes γ1 and γ2 satisfy γ1 γ2 ∥L∥2op ≤ 1, where ∥ · ∥opis the operator norm.

6

Chapter 3

Image Restoration
3.1

Hyperspectral Image Denoising

A hyperspectral (HS) image has 1D spectral information, including invisible light andnarrow 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 [19–21]. Observed HS images are oftenaﬀected by noise because of the small amount of light in narrow wavelength and/or sensor
failure. Thus, we need some robust methods for restoring desirable HS images from suchdegraded observations.
Most HS image restoration methods are based on optimization, which consists of aregularization term and a data-fidelity term. The regularization term evaluates apriori
knowledge on HS images, and plays an essential role for robust 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 diﬀerences to exploit the piecewise-smooth structures of an image. Many TV-based approaches [6,7,9,22] have been proposed for HS image
restoration. Besides, LRM-based approaches [5,23] exploit the underlying low-rank structure in the spectral direction of an HS image. A popular example is the so-called Low-rankmatrix recovery (LRMR) [5].
Many recent methods [9, 12, 13, 24–32] combine TV-based and LRM-based approaches,
and in general, they perform better than approaches using either regularization. This isbecause TV-based approaches model the spatial structure of an HS image, whereas LRMbased 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. Tab. 3.1 summarizes the features of the methods
reviewed in this section.

7

Table 3.1: The feature of existing methods for HS image restoration.
FeatureMethods
HTV [6]SSAHTV [6]SSTV [7]ASSTV [22]
LRM [5, 23]
LNWTV + LRM [9, 24]
HTV + LRM [12]HTV + LRM [25, 26]
ASSTV + LRM [27, 32]
SSTV + LRM [13, 28, 29]
SSTV + LRM [30, 31]
proposed

SpatialCorrelation
⃝
⃝
△
⃝
×
⃝
⃝
⃝
⃝
⃝
⃝
⃝

Spectral
Correlation
×
△
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝

Convexity

Hyperparametersconvexconvexconvexconvex
nonconvex
convex
nonconvex
convex
nonconvex
convex
nonconvex
convexinterdependentinterdependentinterdependentinterdependent
independentinterdependentinterdependentinterdependentinterdependentinterdependentinterdependent
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 hardconstraints on data-fidelity. The regularization, named Hybrid Spatio-Spectral Total Variation (HSSTV), is designed to evaluate two types of local diﬀerences: direct local spatial
diﬀerences and local spatio-spectral diﬀerences in a unified manner to eﬀectively 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 andspectral distortion without combining LRM. Moreover, the constrained-type data-fidelity
in the proposed method enables us to translate interdependent hyperparameters to theupper bounds of the degree of data-fidelity that can be determined based only on thenoise intensity. As a result, the proposed method has no interdependent hyperparameter.
We also develop an eﬃcient algorithm for solving the optimization problem based on thewell-known alternating direction method of multipliers (ADMM) [1, 3, 4].

3.1.1

Related Works

We elaborate on existing HS image denoising methods based on optimization.
TV-Based Methods
The methods proposed in [6, 7, 22] restore a desirable HS image by solving a convex
optimization problem involving TV-based regularization. Let ū ∈ RN B be the desirable
HS image, and the authors assume that an observation v ∈ RN B is modeled as follows:v = ū + s + n,
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,8salt-and-pepper noise, and line noise. The observation and the restoration problem of themethods are given by the following forms:
min ∥v − u − s∥2 + λ1 RTV (u) + λ2 ∥s∥1 ,
u,s

(3.1)

where RTV 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 highquality restoration. However, the hyperparameters are interdependent, which means thatsuitable 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 thehyperparameters are an essential but troublesome task.
⊤ ⊤ ∈ R2N B×N B be spatial
In the following, we explain each TV. Let D = (D⊤v Dh )
diﬀerences operator with Dv and Dh being vertical and horizontal diﬀerences operator,
respectively, and spectral diﬀerences operator are Db ∈ RN B×N B . In [6, 7, 22], HTV,ASSTV, and SSTV are defined as follows:
HTV(u) := ∥Du∥TV ,(3.2)

ASSTV(u) := τv ∥Dv u∥1 + τh ∥Dh u∥1 + τb ∥Db u∥1 ,SSTV(u) := ∥DDb u∥1 ,

(3.3)
(3.4)

where ∥ · ∥TV is a TV norm, which takes the ℓ2 norm of spacial diﬀerence vectors for allband and then summing up for all spatial pixels, and τv , τh , and τb are the weight ofthe vertical, horizontal, and spectral diﬀerences. HTV evaluates direct spatial piecewisesmoothness and can be seen as a generalization of the standard color TV [33]. 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 (Fig. 3.1, blue line). The
weights τv , τh , and τb in (3.3) balance the smoothness of vertical, horizontal, and spectral
diﬀerences, 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 apriori knowledge on HS images using spatio-spectral piecewisesmoothness. It is derived by calculating spatial diﬀerences through spectral diﬀerences(Fig. 3.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 byheavy noise and/or degradation.
LRM-Based Method
LRMR [5] 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 thelocal details, LRMR restores a desirable HS image through patch-wise processing. Each9

direct vertical
diﬀerence

HSSTV
SSTV

direct horizontal
diﬀerencevertical-spectral
diﬀerence

ASSTV

direct spectral
diﬀerence

horizontal-spectral
diﬀerence

Figure 3.1: Calculation of local diﬀerences in SSTV, ASSTV, and our HSSTV. SSTV
evaluates the ℓ1 norm of spatio-spectral diﬀerences (yellow line). ASSTV evaluates theℓ1 norm of direct spatial and spectral diﬀerences (blue line). HSSTV evaluates the mixedℓ1,p norm of direct spatial and spatio-spectral diﬀerences (red line).
patch is a local cube of the size of q × q × B, and LRMR handles it as a matrix of sizeq 2 ×B obtained by lexicographically arranging the spatial vectors in the patch cube in therow direction. The observation model is expressed like Section 3.1.1, and the restoration
problem is formulated as follows:min ∥Vi,j − Ui,j − Si,j ∥2F s.t. rank(Ui,j ) ≤ r, card(Si,j ) ≤ k,Ui,j ,Si,j(3.5)

where Ui,j , Vi,j , and Si,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 thesparse noise by limiting the number of the rank of Ui,j and the cardinality of Si,j using

10

r and k, respectively. Thanks to the design, LRMR achieves high-quality restoration forespecially 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
Prob. (3.5).
Combined Method
The methods [9,12,13,24–32,34] combine TV-based and LRM-based approaches. Since
they can evaluate multiple types of apriori knowledge, i.e., piecewise-smoothness andlow rankness, they can restore a more desirable HS image than the approaches only
using TV-based or LRM-based regularization. Some methods [9, 13, 24–26, 28, 29, 32]
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 datafidelity term(s) simultaneously in one objective function; they require to carefully control
the hyperparameters balancing these terms. Since the hyperparameters rely on both theregularizations and the noise intensity on an observation, i.e., the hyperparameters areinterdependent, the hyperparameter settings are often troublesome tasks.

3.1.2

Proposed method

Hybrid Spatio-Spectral Total Variation
First, we propose a new regularization technique for HS image restoration, named
HSSTV. HSSTV simultaneously handles both direct local spatial diﬀerences and local
spatio-spectral diﬀerences of an HS image. Then, HSSTV is defined by()
DDb
HSSTV(u) := ∥Aω u∥1,p with Aω :=
,
(3.6)
ωD
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 (3.6), DDb u and Du correspond to local spatio-spectral and direct local spatial
diﬀerences, respectively, as shown in Fig. 3.1 (red lines). The weight ω adjusts therelative importance of direct spatial piecewise-smoothness to spatio-spectral piecewisesmoothness. HSSTV evaluates two kinds of smoothness by taking the ℓp norm (p = 1
or 2) of these diﬀerences 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 weset ω = 0 and p = 1, HSSTV recovers SSTV as (3.4), meaning that HSSTV can be seen
as a generalization of SSTV.

11

observationSSTVℓ1 -HSSTVℓ1,2 -HSSTV

Figure 3.2: Restored HS images from an observation contaminated by similar noise inadjacent bands (the upper half area) and random noise (the lower half).
As reviewed in Section 3.1.1, since SSTV only evaluates spatio-spectral piecewisesmoothness, it cannot remove similar noise in adjacent bands. The direct spatial differences in HSSTV help to remove such noise. Fig. 3.2 is restored HS images from an
observation contaminated by similar noise in adjacent bands (the upper half area) andrandom 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 eﬀectively removes all noise. However, since minimizing the direct spatial diﬀerences strongly promotes spatial piecewisesmoothness, HSSTV produces spatial oversmoothing when the weight ω is large. Thus,
the weight ω should be set to less than one, as demonstrated in Sec. 3.1.3.
HS Image Denoising by HSSTV
We consider restoring a desirable HS image ū ∈ RN B from an observation v ∈ RM (M ≤N B) contaminated by a Gaussian-sparse mixed noise. The observation model is given by
the following form:v = ū + n + s,
(3.7)
where n ∈ RM is Gaussian noise with the standard deviation σ, and s ∈ RM is a sparsenoise.
Based on the above model, we formulate HS image restoration using HSSTV as thefollowing optimization problem:v := {x ∈ RM |∥v − x∥ ≤ ε},u + s ∈ B2,εmin HSSTV(u) s.t.  s ∈ B1,η := {x ∈ RM |∥x∥1 ≤ η},
(3.8)
u,s
B
N
u ∈ [µu , µu ] ,v is a v-centered ℓ -norm ball with the radius ε > 0, B
where B2,ε2
1,η is a 0-centered ℓ1 B
N
is a dynamic range of an HS image.
norm ball with the radius η > 0, and [µu , µu ]Here, µu and µu are the minimum and maximum value of HS images, and µu < µu .
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 tothe observation v and the sparse noise s, respectively. As mentioned in [8, 10, 12, 27, 35–42], such a constraint-type data-fidelity enables us to translate the hyperparameter(s)12balancing between regularization and data-fidelity like λ1 and λ2 in (3.1) to the upperbound 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, Prob. (3.8)
is a constrained convex optimization problem. In this paper, we adopt ADMM (seeSec. 2.2) for solving the problem. In what follows, we reformulate Prob. (3.8) into
Prob. (2.7).
We remove the hard constraints from Prob. (3.8) by using the indicator functions of
the constraints. Prob. (3.8) can be rewritten asv (u + s) + ιBmin ∥Aω u∥1,p + ιB2,ε1,η (s) + ι[µu ,µu ]N B (u).

u,s

(3.9)

Note that from the definition of the indicator function, Prob. (3.9) is exactly equal to
Prob. (3.8). By letting
f :RN B → R2 : u 7→ (0, 0),(3.10)

g :R5N B+2M → R ∪ {∞} : (z1 , z2 , z3 , z4 ) 7→v (z2 ) + ιB
∥z1 ∥1,p + ιB2,ε1,ε (z3 ) + ι[µu ,µu ]N B (z4 ),G :RN B → R5N B+2M : u 7→ (Aω u, u + s, s, u).

(3.11)(3.12)

Prob. (3.9) is reduced to Prob. (2.7). The resulting algorithm based on ADMM is summarized in Alg. 1.
Algorithm 1: ADMM method for Prob. (3.8)
(0)

(0)

(0)

(0)

(0)

(0)

(0)

(0)

input : z1 , z2 , z3 , z4 , d1 , d2 , d3 , d4
1
2

3
4

while A stopping criterion is not satisfied do(u(n+1) , s(n+1) ) =
(n)
(n)
(n)
(n)
(n)
(n)
(n)
(n)
1
(∥z1 −Aω u−d1 ∥2 +∥z2 −(u+s)−d2 ∥2 +∥z3 −s−d3 ∥2 +∥z4 −u−d4 ∥2 );argmin 2γ
u,s
(n+1)

(n)
= proxγ∥·∥1,p (Aω u(n+1) + d1 );
(n)
(n+1)
= proxγιBv (u(n+1) + s(n+1) + d2 );
z2

z1

2,ε

5
6
7
8
9
10
11

(n+1)

z3

(n+1)

z4

(n+1)

= proxγιB

1,η

= proxγι(n)(n)(s(n+1) + d3 );(n)[µu ,µu ]N B

(u(n+1) + d4 );
(n+1)

d1
= d1 + Aω u(n+1) − z1
;
(n)
(n+1)
(n+1)
(n+1)
(n+1)
;
d2
= d2 + u
+s
− z2
(n+1)
(n)
(n+1)
(n+1)
d3
= d3 + s
− z3
;
(n+1)
(n)
(n+1)
d4
= d4 + u(n+1) − z4;n ← n + 1;

The update of u and s in Alg. 1 comes down to the following forms:()
3 −1
u(n+1) = A⊤
A
+
I
RHS,
ω ω
2
1 (n)
1 (n)
(n)
(n)
(n)
(n)
(n)
(n)
RHS = A⊤ω (z1 − d1 ) + (z2 − d2 ) − (z3 − d3 ) + (z4 − d4 ),2213

s(n+1) =

1 (n)
(n)
(n)
(n)
(z − u(n+1) − d2 + z3 − d3 ),2 2Since the update of u and s in Alg. 1 is strictly-convex quadratic minimization, one canobtain these updated forms by diﬀerentiating it.
For the update of z1 , the proximity operators are reduced to simple soft-thresholding
type operations described in Sec. 2.1. The updates of z2 , z3 , and z4 use the proximity
operators of the indicator functions of each set, i.e., (2.4), (2.5), and (2.6), respectively.

3.1.3

Experimets

We demonstrate the advantages of the proposed method over three noise type: Gaussian
noise, Gaussian-sparse mixed noise, and real noise. All the experiments were conductedby MATLAB 2018a.
To quantitatively evaluate restoration performance, we used the peak signal-to-noise
ratio (PSNR) [dB] index and the structural similarity (SSIM) [43] index between a groundtruth ū and a restored HS image u. PSNR is defined by 10 log10 (M AXI2 /M SE), where
M AXI is the max value of the dynamic range of HS images, and M SE is the mean squareerror between ground-truth 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 visionsystem, which is defined as follows:
1 ∑ (2µui µūi + C1 )(2σui ūi + C2 ),P(µ2ui + µ2ūi + C1 )(σu2 i + σū2 i + C2 )P

SSIM(u, ū) =i=1

where ui and ūi are the ith pixel-centered local patches of a restored HS image and aground-truth HS image, respectively, P is the number of patches, µui and µūi is the average values of the local patches of the restored and ground-truth HS images, respectively,σui and σūi represent the variances of ui and ūi , respectively, and σui ūi denotes the covariance between ui and ūi . Moreover, C1 and C2 are two constants, which avoid thenumerical instability when either µ2ui + µ2ūi or σu2 i + σū2 i is very close to zero. SSIM gives
a normalized score between zero and one, where the maximum value means that u equals
to ū.
Gaussian noise removalFirst, we conducted on Gaussian noise removal experiments. In these experiments, we
generate a noisy observation based on (3.7) with s = 0 and restore a desirable HS image
u. We use 13 HS images taken from the SpecTIR [44], MultiSpec [45] and GIC [46] astest images, where their dynamic range were normalized into [0, 1]. Here, the proposed
denoising method assumes that the observation is contaminated by Gaussian and sparsenoise, and then we improved the problem and the algorithm for Gaussian noise removal.
Specifically, we formulate a Gaussian noise removal problem as follows:[v ,u ∈ B2,εmin HSSTV(u) s.t.
uu ∈ [µu , µu ]N B ,14

Table 3.2: PSNR[dB] of the results on Gaussian noise removal experiments.
Beltsvillen = 256m = 256b = 32
Suwanneen = 256m = 256b = 32DCn = 256m = 256b = 32Cupriten = 256m = 256b = 32Renon = 256m = 256b = 32
Botswanan = 256m = 256b = 32
IndianPinesn = 145m = 145b = 32KSC
n = 256m = 256b = 32
PaviaLeftn = 216m = 216b = 32
PaviaRightn = 256m = 256b = 32
PaviaUn = 256m = 256b = 32
Salinasn = 217m = 217b = 32
SalinaAn = 83m = 86b = 32

σ
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3

HTV
32.19
29.33
26.92
25.79
33.03
30.04
27.32
25.88
30.06
26.87
24.29
23.05
34.32
31.62
29.31
28.13
32.01
28.98
26.50
25.19
31.35
28.27
25.59
24.26
33.13
31.46
30.22
29.67
33.89
31.21
28.96
27.82
31.23
27.97
25.23
23.93
31.90
28.58
25.73
24.33
32.16
28.90
26.07
24.65
35.93
33.05
30.56
29.16
36.17
33.15
29.93
28.15

SSTV
35.64
31.11
26.03
22.81
37.28
31.88
27.54
25.03
35.56
30.93
25.86
22.63
37.64
31.97
27.22
24.50
36.76
31.79
26.57
24.47
35.32
31.17
26.37
23.21
33.70
30.67
26.38
24.25
37.67
32.08
27.41
25.21
37.66
32.55
28.36
25.48
37.62
32.25
28.59
25.76
37.65
32.77
28.60
25.62
37.61
31.93
26.58
24.09
37.54
32.25
26.63
23.65

ASSTV
31.99
28.58
25.36
23.90
33.34
30.44
27.83
26.42
30.87
27.65
24.92
23.52
33.40
30.34
27.41
25.71
32.25
29.23
26.78
25.54
30.98
27.62
24.66
23.12
32.75
30.40
28.28
26.79
34.29
31.40
28.86
27.52
32.07
28.78
26.04
24.69
33.00
29.61
26.71
25.26
33.18
29.88
27.06
25.63
35.34
32.37
29.62
28.15
34.05
30.20
26.20
24.11

15

BM4D
37.01
33.39
30.03
28.04
38.05
34.56
31.53
29.79
36.16
31.88
28.27
26.43
39.10
35.65
32.39
30.60
37.46
33.44
29.91
28.12
36.01
32.33
28.94
27.07
34.59
32.17
30.19
29.04
39.36
35.30
31.86
30.10
38.09
33.87
29.92
27.78
38.77
34.60
30.53
28.34
38.96
34.89
30.86
28.74
40.27
36.89
33.61
31.73
39.60
35.73
31.83
29.61

Proposed (p = 1)
36.41
33.15
30.23
28.68
38.69
35.17
31.94
30.21
35.96
32.16
28.72
26.93
39.58
36.29
33.32
31.75
37.65
34.07
30.83
29.15
35.73
32.46
29.43
27.81
34.34
32.41
30.95
30.26
39.44
36.15
33.14
31.59
38.23
34.39
30.83
28.94
38.63
34.94
31.43
29.50
38.81
35.20
31.74
29.82
40.42
37.30
34.31
32.63
39.92
36.72
33.64
31.96

Proposed (p = 2)
36.19
32.79
29.80
28.30
38.73
35.19
32.01
30.30
35.57
31.74
28.33
26.57
39.56
36.26
33.32
31.77
37.43
33.83
30.65
29.03
35.65
32.30
29.26
27.67
34.03
32.05
30.67
30.07
39.41
36.09
33.10
31.54
37.98
34.22
30.71
28.83
38.38
34.69
31.21
29.31
38.60
35.01
31.61
29.71
40.44
37.36
34.43
32.80
39.81
36.46
33.30
31.63

Table 3.3: SSIM of the results on Gaussian noise removal experiments.
Beltsville

SuwanneeDC

CupriteReno

Botswana

IndianPines

KSC

PaviaLeft

PaviaRight

PaviaU

SalinasSalinaA

σ
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3
0.05
0.1
0.2
0.3

HTV
0.83700.74830.65830.61500.87230.79590.70910.65650.85840.73070.57720.49410.87420.80950.74890.72000.86040.75600.65030.59730.84950.75280.64050.57650.82920.77520.74140.73060.88790.81800.75080.71710.86120.73540.57140.47820.86250.74840.60650.52520.87190.77840.65850.58580.92020.88370.84800.82930.92310.89150.85090.8242SSTV
0.90760.7724
0.52300.37090.95210.8646
0.67720.5401
0.94860.86790.69270.5456
0.94120.83540.61150.46090.95270.86350.69190.55880.93710.86290.6996
0.55740.84140.69030.44070.31980.93610.81850.63370.4786
0.96990.91690.78770.6546
0.96170.89110.73550.59270.96020.88790.72200.57590.9246
0.79490.54150.40670.93640.82380.58800.4288ASSTV
0.86020.77720.68110.6216
0.89370.82630.75120.70580.89290.7866
0.64170.55410.88700.82840.76990.73910.88590.79730.69350.63680.85920.76680.65780.59200.85200.78380.74410.72800.90350.83960.77270.73520.88550.7776
0.63170.54150.88940.79040.66050.58040.89430.81330.70590.63860.92970.89510.85780.83820.92520.88210.82550.7932

16

BM4D
0.93550.87370.78330.71700.95630.91040.8406
0.7866
0.96350.91400.81840.73350.95380.90810.83420.77980.95990.90830.81770.74690.93570.87210.78020.71170.88290.81220.74450.69940.96510.91620.83460.77570.96830.92030.82240.73430.96530.91750.81630.7256
0.96490.92260.83780.76820.96150.93320.88820.8486
0.95800.91640.8416
0.7837

Proposed (p = 1)
0.93050.86680.77980.72390.96480.9276
0.86680.82090.96330.91560.82550.75010.96020.92260.86530.82440.96490.92200.8452
0.78480.94680.90200.83130.77650.87580.81640.76730.74640.96900.93670.87990.8384
0.97200.93460.85890.79000.97030.93510.86590.80380.96900.93360.86930.81540.96380.94050.90920.88680.96700.94160.90490.8782

Proposed (p = 2)
0.92620.85620.76300.70800.96570.92850.86720.8211
0.95860.90650.81310.73770.96180.92430.86650.82580.96380.91930.84190.7821
0.94770.90480.83490.7806
0.86350.79890.75540.7401
0.96910.93580.87860.83690.97050.93160.85370.78340.96920.93210.8592
0.79420.96750.93000.86470.81070.9647
0.94230.91210.89080.96880.94430.90680.8797

Suwannee

20.01, 0.2747

30.04, 0.7959

31.88, 0.8646

30.44, 0.8263

ground-truth

observation

HTV

SSTV

ASSTV

34.56, 0.9104

35.13, 0.9276

35.19, 0.9285

BM4D

proposed

proposed(p = 1)(p = 2)

DC

13.98, 0.1787

24.29, 0.5772

25.86, 0.6927

24.92, 0.6417

ground-truth

observation

HTV

SSTV

ASSTV

28.27, 0.8184

28.72, 0.8255

28.33, 0.8131

BM4D

proposed

proposed(p = 1)(p = 2)

Figure 3.3: Resulting HS images with their PSNR[dB] (left) and SSIM (right) in theGaussian denoising experiment (top: Suwannee, σ = 0.1, bottom: DC, σ = 0.2).

17σ = 0.05

σ = 0.1

σ = 0.2

σ = 0.3
PSNR[dB]

SSIM

Figure 3.4: PSNR[dB] or SSIM versus ω in (3.6) in the Gaussian denoising experiment.

18

(a) Bandwise PSNR(b) Bandwise SSIM(c) Spatial response(d) Spectral response

Figure 3.5: Bandwise PSNR[dB] and SSIM and spatial and spectral responses in theGaussian denoising experiment (Suwannee).

and solve it by ADMM.
We compared the proposed method with HTV [6], SSTV [7] and ASSTV [8]. For a faircomparison, we replaced HSSTV in Prob. (3.8) with HTV, SSTV or ASSTV, and solved
the problem by ADMM. We also compare HSSTV with BM4D [47], which is known to bea state-of-the-art denoising method for 3D signals. Note that BM4D and a recent CNNbased HS image denoising method [48] cannot be represented as explicit regularization
functions and are fully customized to denoising tasks. In contrast, our HSSTV can beused as a building block in various HS image restoration methods based on optimization.
In addition, BM4D is a nonlocal method, and thus, it has several limitations. Meanwhile,
CNN-based methods strongly depend on what training data are used, which means thatthey cannot adapt to a wide range of noise intensity. Thus, the design concepts of thesemethods are diﬀerent from TVs, and the comparison of HSSTV and BM4D is just areference. We also remark that we do not compare HSSTV with LRM-based techniques.
This is because our focus is to evaluate the potential of local TV regularization techniques

19including HTV, SSTV, ASSTV, and our√HSSTV.
The radius ε in Prob. (3.8) was set to N Bσ 2 . In the ASSTV case, we set the weights
τv = τh = τb = 1. In the BM4D case, we used the program code distributed by the authorsof [47]. We set the max iteration number, the stepsize γ and the stopping criterion of
ADMM to 5000, 0.1 and ∥u(n) − u(n+1) ∥ < 0.01, respectively.
In Tab. 3.2 and 3.3, we show PSNR and SSIM of the denoised HS images by eachmethod for various σ and HS images, respectively. The balancing weight ω in HSSTV
is set to 0.08 for the ℓ1 -norm case and 0.06 for the ℓ1,2 -norm case. For all HS images, σ
and quality measures, HSSTV outperforms HTV, SSTV, and ASSTV. In addition, eventhough HSSTV does not utilize nonlocal information, the denoising ability is better than
BM4D for most cases. We also found that SSTV does not work well when σ is large. On
the other hand, HSSTV is eﬀective for a wide range of noise intensity. This would be
because HSSTV evaluates direct spatial piecewise smoothness.
Fig. 3.3 shows the resulting images on Suwannee (σ = 0.1, top) and DC (σ = 0.2,
bottom) with their PSNR (left) and SSIM (right). Here, we depict these HS images as
RGB images, where R, G and B are 8th, 16th, and 32rd bands of HS images. In Fig. 3.3,
we show both the results and the magnified image of the small range, because one can more
see the advantage of the results by the proposed method. One can see that the results by
HTV and ASSTV lose spacial details, and noise remains in the results by SSTV. In the caseof σ = 0.1, BM4D can restore a high quality HS image, but it produces spectral artifacts
in the case of σ = 0.2. This means that BM4D cannot preserve spectral information whenan observed image is heavily contaminated by noise. In contrast, HSSTV can restore HS
images preserving both details and spectral information without artifacts.
Fig. 3.4 plots PSNR or SSIM of the results by HSSTV versus various ω changed from0.01 to 0.2, where the values of PSNR and SSIM are averaged over the 13 HS images.
One can see that ω ∈ [0.05, 0.1] is a good choice in most cases.
Fig. 3.5 plots bandwise PSNR and SSIM (left) and spatial and spectral responses (right)of the denoised Suwannee HS image, where σ = 0.1. The graphs regarding bandwise
PSNR and SSIM show that HSSTV achieves higher-quality restoration than HTV, SSTV,
and ASSTV for all bands and BM4D for most bands. 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 responseof the results by HTV and ASSTV are too smooth compared with the true one. On
the other hand, there exist undesirable variations in the spatial response of the result by
SSTV. In contrast, BM4D and HSSTV restore similar responses to the true one. In thegraph (d), one can see that (i) the shape of the spectral responses of the results by HTV
and SSTV are similar to the that of true one, but the mean values are larger than the true
one, (ii) the spectral response of the results by ASSTV and BM4D are too smooth anddiﬀerent from the true one, and (iii) HSSTV can restore a spectral response very similar
to the true one.

20

Salinasground-truth

PaviaUground-truth

17.06, 0.1908
observation

32.59, 0.8997
HTV

32.66, 0.9003
SSAHTV

35.86, 0.9015
SSTV

32.83, 0.9163
ASSTV

31.82, 0.8270
LRMR

36.74, 0.9509
LRTV

34.57, 0.9305
LLRGTV37.60, 0.9561
proposed

37.65, 0.9564
proposed(p = 1)(p = 2)15.39, 0.1877
observation

25.24, 0.6776
HTV

25.29, 0.6785
SSAHTV

30.21, 0.8444
SSTV

27.42, 0.7678
ASSTV

29.43, 0.8103
LRMR

28.90, 0.8627
LRTV

29.18, 0.8655
LLRGTV

31.04, 0.8935
proposed

30.80, 0.8855
proposed(p = 1)(p = 2)

Figure 3.6: Resulting HS images with their PSNR[dB] (left) and SSIM (right) in the mixed
noise removal experiment (top: Salinas, the noise level (i), bottom: PaviaU, the noise level(ii)).

21

(a) band-wise PSNR(b) band-wise SSIM(c) Spatial response(d) Spectral response

Figure 3.7: Band-wise PSNR[dB] and SSIM and spatial and spectral responses in themixed noise removal experiment (Suwannee).

22

Table 3.4: Parameter settings for ASSTV, LRMR, LRTV, and the proposed method.
(i) (σ, sp , lv = lh )Noise Level(ii) (0.1, 0.05, 0.05)
Parameters= (0.05, 0.04, 0.04)
τv = τh
1
ASSTV
τb
3
2
r
3
LRMR
k
sp + lv + lh − lv lh (the rate of sparse noise)
r
2
LRTV
τ
0.005
0.008
r
2
LLRGTVλ0.1
τ
0.01
proposed
ω
0.04

Mixed noise removal
We conducted on the simulation experiments for Gaussian-sparse mixed noise removal.
In the experiments, we used 13 HS images as like Gaussian noise removal experiment.
The observed HS images included an additive white Gaussian noise n with the standard
deviation σ and sparse noise s, which is equal to (3.7). 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 sp , lv , and lh , respectively. We generated noisy HS images byadding two types of mixed noise to ground-truth HS images: (i) (σ, sp , lv = lh ) =(0.05, 0.04, 0.04), (ii) (0.1, 0.05, 0.05).
The proposed method was compared with HTV [6], SSAHTV [6], SSTV [7], and ASSTV [8]. For a fair comparison, we replaced HSSTV in Prob. (3.8) with HTV, SSAHTV,
SSTV, or ASSTV and solved the problem by ADMM. We also compared our proposed
method with LRMR [5], TV-regularized low-rank matrix factorization (LRTV) [12], and local low-rank matrix recovery and global spatial-spectral TV (LLRGTV) [27]. We did notcompare our proposed method with recent CNN-based HS image denoising methods [48].
The CNN-based methods cannot be represented as explicit regularization functions andare fully customized to denoising tasks. In contrast, our proposed method can be used asa building block in various HS image restoration methods based on optimization. Meanwhile, CNN-based methods strongly depend on what training data are used, which meansthat they cannot adapt to a wide range of noise intensity. Thus, the design concepts ofthese methods are diﬀerent from TVs and√ LRM-based approaches.
We set ε and η in Prob. (3.8) as 0.83 N B(1 − (sp (1 − lv − lh ) + lv + lh − lv lh ))σ 2 andN B(0.45sp +(lv +lh )vave −lv lh vave ), respectively, where vave is the average of the observed
image. Tab. 3.4 shows the parameters settings for ASSTV, LRMR, LRTV, LLRGTV, andthe proposed method. We set these parameters to achieve the best performance for eachmethod. We set the max iteration number, the stepsize γ, and the stopping criterion of
ADMM to 10000, 0.05, and ∥u(n) − u(n+1) ∥ < 0.01, respectively.

23

Table 3.5: PSNR[dB] (top) and SSIM (bottom) in mixed noise removal experiments.
HS image
Beltsville
SuwanneeDC
CupriteReno
Botswana
PSNR
IndianPines
KSC
PaviaLeft
PaviaRight
PaviaU
SalinasSalinaA
Beltsville
SuwanneeDC
CupriteReno
Botswana
SSIM

IndianPines
KSC
PaviaLeft
PaviaRight
PaviaU
SalinasSalinaA

Noise
Level
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)
(i)
(ii)

HTV

SSAHTV

SSTV

ASSTV

LRMR

LRTV

LLRGTV

29.43
26.40
30.14
26.70
26.46
23.84
31.67
28.20
28.53
25.56
27.98
25.21
31.05
28.57
30.17
28.03
27.62
24.74
26.93
24.90
27.92
25.24
32.59
28.88
32.54
28.69
0.79020.69540.84060.75420.7622
0.61890.85500.78490.78180.66400.7896
0.68100.81180.77130.82710.75980.77520.61020.77690.64740.79730.6776
0.89970.85700.91290.8793

29.47
26.43
30.18
26.74
26.51
23.88
31.68
28.21
28.57
25.61
28.05
25.25
31.06
28.57
30.25
28.06
27.70
24.78
27.35
25.16
28.04
25.29
32.64
28.91
32.65
28.80
0.79040.69590.84100.75520.76330.62010.85520.78520.78190.66450.79000.68200.81200.77130.82780.76020.77700.61160.77720.6471
0.7986
0.67850.90020.85750.91370.8803

33.66
28.42
34.59
29.55
33.03
27.71
34.42
29.86
34.37
28.11
33.32
28.55
31.45
27.82
34.74
29.23
35.57
29.93
34.54
30.70
35.52
30.21
35.86
28.19
35.29
29.67
0.88560.70570.93530.81460.92740.80920.91790.77170.93220.80450.92020.81750.80150.62290.9116
0.78850.95930.87550.94940.86350.94520.84440.90150.71170.91340.7789

27.16
24.60
32.60
28.71
28.80
25.25
29.14
26.57
30.49
26.95
26.47
24.01
29.07
26.72
31.64
28.62
30.91
26.71
31.13
27.23
31.65
27.42
32.83
28.99
28.12
25.19
0.81110.71770.90520.82260.8676
0.72110.86320.79530.88320.75390.81990.70950.76710.73030.89220.80640.88280.72670.88620.74930.88910.76780.91630.87320.84680.8110

30.91
27.13
30.30
26.90
31.71
27.35
30.16
27.32
32.21
28.47
31.62
28.31
28.96
25.14
33.74
30.19
33.01
29.46
33.33
29.82
33.00
29.43
31.82
28.02
31.18
27.67
0.85830.70830.86890.74700.92480.82140.84950.70980.90120.79050.90680.82010.75930.78930.88900.75290.93590.85650.92560.8261
0.91240.81030.82700.66700.86320.7266

35.32
31.22
36.20
31.95
34.78
29.53
28.39
27.94
37.06
31.00
29.00
27.33
26.16
29.82
35.74
30.22
36.49
29.02
35.82
29.08
36.72
28.90
36.74
32.73
28.49
26.10
0.93720.85680.95020.89300.96130.88100.9396
0.88140.95890.8816
0.92820.85640.81900.79390.93850.84270.96120.8791
0.95400.85070.95400.86270.95090.92250.93840.8951

33.61
28.51
33.35
28.50
33.50
28.14
32.67
27.99
34.95
27.99
32.02
28.08
29.74
26.70
34.75
29.53
34.75
29.22
34.17
29.24
34.59
29.17
34.36
29.09
34.07
27.87
0.92780.82480.94310.8622
0.95480.8722
0.94110.85240.95230.86400.93840.87560.82240.74330.93220.82160.96010.8882
0.94700.8551
0.94930.86490.92850.83330.95130.854924

Proposed(p = 1)
34.25
29.79
35.15
31.08
33.36
28.57
34.96
31.63
35.11
29.83
33.61
29.39
31.90
29.26
36.39
31.82
35.98
30.47
35.68
31.59
36.31
31.04
37.60
32.01
36.27
31.68
0.91320.8186
0.95590.91250.94420.86110.94590.9031
0.95310.86790.93430.87450.83350.77850.95420.88090.96610.88980.96160.9086
0.96220.89350.9561
0.92230.94480.9197

Proposed(p = 2)
34.16
29.62
36.01
31.22
33.08
28.32
36.20
31.73
34.96
29.72
33.53
29.35
31.80
29.18
36.33
31.72
35.81
30.24
35.23
31.39
36.17
30.80
37.65
32.12
36.23
31.64
0.90850.80880.95550.91580.93940.85330.94260.90580.95150.86350.93440.87650.82430.76890.95320.87470.96450.88150.95980.9006
0.96100.88550.95640.92400.94160.9195

noise level (i)noise level (ii)Figure 3.8: PSNR[dB] (top) or SSIM (bottom) versus ω in (3.6) in the mixed noise
removal experiment.
In Tab. 3.5, 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 cansee that the PSNR and SSIM of the results by the proposed method are higher thanthose by LLRGTV for most situations. Meanwhile, LRTV outperforms the proposed
method in some cases. This would be because LRTV combines TV-based and LRMbased 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.
Fig. 3.6 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 seethat 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

25

PSNR vs α

SSIM vs α

PSNR vs β

SSIM vs β
DC

KSC

Figure 3.9: PSNR[dB] or SSIM versus α or β on the mixed noise removal experiment(the noise level (ii), left: DC, right: KSC).

26

and spectral information without artifacts.
Fig. 3.7 plots band-wise PSNR and SSIM (top), and spatial and spectral responses(bottom) 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 andLRMR 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 responseof the 243rd row of the 30th band. Similarly, graph (d) plots the spectral response of the243rd row and 107th col. We can see that the spatial response of the results by HTV
and SSAHTV is too smooth compared with the ground-truth 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 ground-truth 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 resultsby SSTV is similar to that of the ground-truth one, but the mean value is larger than
the ground-truth one, (iii) the spectral response of the results by ASSTV is too smoothand diﬀerent from the ground-truth one, and (iv) LRTV, LLRGTV, and the proposed
method can restore a spectral response very similar to the ground-truth one.
To search the suitable ω in (3.6), we changed ω between 0.01 and 0.2 at 0.01 intervals.
Fig. 3.8 plots PSNR or SSIM of the results by the proposed method versus various ω.
Here, the values of PSNR and SSIM are averaged over the 13 HS images. One can seethat ω ∈ [0.03, 0.07] is a good choice. The results show that the suitable ω in the mixed
noise removal experiments is a little smaller than that in the Gaussian noise removal
experiments, and the proposed method can achieve high quality denoising regardless ofnoise type if one set ω ∈ [0.05, 0.07]. ASSTV and LRTV require adjusting the weightτb and the hyperparameter τ newly for diﬀerence noise intensity, respectively, but thesuitable parameter ω in HSSTV is noise-robust.
To verify the sensitivity of the parameter ε and η in Prob. (3.8), we conducted additional mixed noise removal experiments,
where we examined various values of ε and√η. Specifically, we set ε = 0.83α N B(1 − (sp (1 − lv − lh ) + lv + lh − lv lh ))σ 2 and η =βN B(0.45sp + (lv + lh )vave − lv lh vave ), which are hand-optimized values of the parametersand changed α and β from 0.9 to 1.1 at 0.02 interval (the DC and the KSC images andthe noise level (ii)). Fig. 3.9 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 notvary significantly for both image, so the parameters ε and η are independent of both aregularization technique and an observed image. In the SSTV cases, the shapes of theplots are diﬀerent 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.

27

observation

HTV

SSAHTVLRMRLRTV

LLRGTV

SSTV

ASSTV

proposed (p = 1) proposed (p = 2)

Figure 3.10: The denoising results on the real noise removal experiments.
Real Noise Removal
We illustrated the performance of the proposed method for real noise removal. In theexperiments, we restore a desirable HS image from an observation with real noise by the
proposed mixed noise removal method. We selected noisy 16 bands from Suwannee andused it as a real observed HS image v. To maximize the performance of each method,
we searched for suitable values of σ, sp , lv , and lh . We also compared the proposed
method with HTV, SSAHTV, SSTV, ASSTV, LRMR, LRTV, and LLRGTV. We set theparameter ε = 3.1893 and η = 31893 for all TVs. Besides, the parameters τv , τh , and τbin 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.
Fig. 3.10 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.

3.2

Compressed Sensing Reconstruction

Multichannel images with the large number of bands have high-resolution 1D spectraland 2D spatial information, i.e., HS or multispectral (MS) images. The imaging requires
a 2D spectral sensor for capturing 1D spatial information, and thus it is impossible tocapture the images in real-time. To resolve the above problem, CS reconstruction plays
a very important role [49–51]. The CS theory said that a high-dimensional signal canbe reconstructed from a much fewer number of random observations by leveraging thesparsity of the signal, and so one-shot multichannel imaging based on CS is actively

28

studied [52, 53]. Therefore, the CS theory requires a robust reconstruction method tonoise and outliers.
In this section, we focus on CS methods for an HS image. This is because HS images
have the huge number of bands in multichannel image, and CS is essential techniquefor real-time HS imaging. CS reconstruction can be considered as the one of the image
restoration method, and most CS methods are based on optimization. The reconstruction
problems mainly consists of a regularization and a data-fidelity terms. As like with thedenoising methods, the regularization techniques play an important role for the improvement of the CS performance. Besides, many CS methods for HS images [54–57] assume
that observation is contaminated by only Gaussian noise, but optical systems often cause
the other type of noise and/or outliers.
Based on the above discussion, we propose a new robust CS method to noise andoutliers. Specifically, the method estimates a desirable HS image by solving optimization
problem, which evaluates noise and outliers by multiple constraints and utilizes HSSTV
as regularization. To eﬃciently solving the problem, we also provide an iterative solver
based on a primal-dual splitting method [17]. In the experiments, we illustrate the utilityof the proposed framework.

3.2.1

CS Reconstruction under Gaussian Noise Situation

Problem Formulation
Consider to restore an clean HS image ū ∈ RN B from an observation v ∈ RM , whichis cast as inverse problems of the form:v = Φū + n,

(3.13)

where Φ ∈ RM ×N B (M ≤ N B) is a sensing matrix, e.g., random sampling, n is anadditive white Gaussian noise with the standard deviation σ.
Based on the above model, we newly formulate compressed HS image reconstruction
as the following convex optimization problem:[
v := {x ∈ RM |∥x − v∥ ≤ ε},Φu ∈ B2,εmin HSSTV(u) s.t.
(3.14)uu ∈ [µu , µu ]N B .
The first term in (3.14) is a regularization function for HS image restoration. HSSTV simultaneously evaluates spatio-spectral piecewise-smoothness and direct spatial piecewisesmoothness, and we explain HSSTV in Sec. 3.1.2. The first constraint in (3.14) serve as
v is the v-centereddata-fidelity regarding Gaussian noise to the observation v, and B2,εℓ2 -norm ball with the radius ε. The second constraint in (3.14) is the dynamic range ofu (µu < µu ).

29

Algorithm 2: ADMM method for Prob. (3.14)
(0)

(0)

(0)

(0)

(0)

(0)

input : z1 , z2 , z3 , d1 , d2 , d3
12

while A stopping criterion is not satisfied do
(n)

1
(∥z1
u(n+1) = argmin 2γ
u

(n+1)

= proxγ∥·∥1,p (Aω u(n+1) + d1 );(n+1)= proxγιBv (Φu(n+1) + d2 );

3

z1

4

z2

5
6
7
8
9

(n)

(n)

− Aω u − d1 ∥2 + ∥z2

(n)

(n)

− Φu − d2 ∥2 + ∥z3(n)− u − d3 ∥2 );

(n)

(n)

2,ε

(n+1)

z3

= proxγι

(n+1)

(n)

[µu ,µu ]N B

(u(n+1) + d3 );

(n)

(n+1)

d1
= d1 + Aω u(n+1) − z1
;
(n+1)
(n)
(n+1)
d2
= d2 + Φu(n+1) − z2
;
(n+1)
(n)
(n+1)
d3
= d3 + u(n+1) − z3;n ← n + 1;

Optimization
We solve Sec. (3.14) by ADMM. To solve it ADMM, we reformulate Prob. (3.14) into
Prob. (2.7) like Sec. 3.1.2. By letting
f : RN B → R : u 7→ 0,(3.15)

g:R(3.16)

5N B+Mv (z2 ) + ι→ R ∪ {∞} : (z1 , z2 , z3 ) 7→ ∥z1 ∥1,p + ιB2,ε[µu ,µu ]N B (z3 ),G : RN B → R5N B+M : u 7→ (Aω u, Φu, u).

(3.17)

Prob. (3.14) is reduced to Prob. (2.7). The resulting algorithm based on ADMM issummarized in Alg. 2.
The update of u in Alg. 2 comes down to the following forms:()−1⊤
u(n+1) = A⊤
A
+
Φ
Φ
+
I
RHS,
ω ω
⊤
RHS =A⊤ω (z1 − d1 ) + Φ (z2 − d2 ) + (z3 − d3 )
(n)

(n)

(n)

(n)

(n)

(n)

(3.18)

Since the update of u in Alg. 2 is strictly-convex quadratic minimization, one can obtainthese updated forms by diﬀerentiating it.
Here, we should consider the structure of Φ because it aﬀects the matrix inversionin (3.18). If Φ is a block-circulant-with-circulant-blocks (BCCB) matrix [58], we canleverage 3DFFT to eﬃciently solve the inversion in Step 2 with the diﬀerence operators⊤
having a periodic boundary, i.e., A⊤ω Aω + Φ Φ + I can be diagonalized by the 3DFFT
matrix and its inverse. If Φ is a semiorthogonal matrix, i.e., ΦΦ⊤ = αI (α > 0), we leavevv ◦ Φ in (3.16) and Φuit to the update of z2 , which means that we replace ιB2,εby ιB2,εv ◦ Φ, in this case, can beby u in (3.17). This is because the proximity operator of ιB2,εcomputed by using [18, Table 1.1-x] as follows:v (Φx) − Φx).
proxγιBv ◦Φ (x) = x + α−1 Φ⊤ (PB2,ε2,ε30

If Φ is a sparse matrix, we oﬀer to use a preconditioned conjugate gradient method [59] forapproximately solving the inversion or to apply primal-dual splitting methods [17, 60, 61]instead of ADMM. Here, 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 [62–65] might be usefulfor reducing the computational cost.
For the update of z1 , the proximity operators are reduced to simple soft-thresholding
type operations described in Sec. 2.1. The updates of z2 and z3 use the proximity operators of the indicator functions of each set, i.e., (2.4) and (2.6), respectively.
Experiments
We experimented on CS reconstruction under Gaussian noise situation. In this experiment, we generate an observation based on (3.13) and reconstruct a desirable HS image
from it. For test HS images, we took 13 HS images as like the Gaussian noise removal
experiment. We assume that Φ ∈ RM ×N B in (3.13) is a random sampling matrix withthe sampling rate m (0 < m < 1 and M = mN B). Here, since Φ is a semiorthogonal
matrix, we can eﬃciently solve the problem, as explained in the optimization section. We
set the standard deviation of Gaussian noise σ as 0.1 and the sampling rate m of Φ as
0.2 or 0.4.
The proposed method was compared with HTV, SSAHTV, SSTV, and ASSTV. For faircomparison, we replaced HSSTV in Prob. (3.14) with the TVs and solved the problem
by ADMM. In the ASSTV case, we set the parameters (τv , τh , τb ) = (1, 1, 0.5),√ whichexperimentally achieves the best performance. We set the hyperparameter ε = mN Bσ 2
in (3.14), and ω in (3.6) is set to 0.08/0.06 in the ℓ1 /ℓ1,2 -HSSTV case, respectively.
Tab. 3.6 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 thanthat by HTV, SSAHTV, SSTV, and ASSTV.
Fig. 3.11 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
diﬀerent from the color of the ground-truth 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. 3.12 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 seethat 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 responsesof the same position in the denoising experiments. 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 ground-truth one. In graph (d), HTV and SSAHTV generate
undesirable variation, and ASSTV causes oversmoothing. Thanks to the evaluation of31

Table 3.6: PSNR[dB] (top) and SSIM (bottom) in the CS reconstruction experiment
under Gaussian noise situation.
HS image
Beltsville
SuwanneeDC
CupriteReno
Botswana
PSNR

IndianPines
KSC
PaviaLeft
PaviaRight
PaviaU
SalinasSalinasA
Beltsville
SuwanneeDC
CupriteReno
Botswana

SSIM

IndianPines
KSC
PaviaLeft
PaviaRight
PaviaU
SalinasSalinasA

m
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2
0.4
0.2

HTV
27.46
26.23
27.97
26.47
24.69
23.31
29.94
28.77
26.99
25.57
26.10
24.66
30.54
29.99
29.30
28.11
25.66
24.26
25.83
24.30
26.49
24.95
31.19
29.94
30.67
28.68
0.68290.63630.73320.68100.6096
0.52150.76650.73680.67690.6202
0.66830.60140.74970.7366
0.76600.73180.60820.51030.63570.55020.68670.61380.85770.84040.86470.8387SSAHTV
27.49
26.25
28.02
26.50
24.73
23.33
29.96
28.77
27.05
25.61
26.15
24.69
30.55
29.99
29.33
28.12
25.69
24.27
25.85
24.30
26.53
24.97
31.24
29.98
30.82
28.75
0.69400.64930.74970.70070.6242
0.53840.78040.75250.68680.63260.68030.61620.77770.76580.77420.74100.62050.52510.64230.55840.69560.62420.86720.85160.88710.8655

SSTV
27.53
24.34
28.49
25.69
27.33
24.16
28.21
25.79
27.82
25.57
27.81
24.79
27.65
25.11
28.34
27.00
29.66
27.17
29.86
27.54
30.02
27.35
27.69
25.28
27.93
24.15
0.60130.43480.73770.57390.75220.61200.68260.50570.74140.62760.75510.62250.5066
0.34880.68140.60080.8386
0.73190.79620.69170.78300.66230.61530.46200.65950.4810

32

ASSTV
26.51
24.12
27.68
25.39
24.71
22.69
28.59
26.38
26.49
24.52
25.13
22.86
29.55
28.19
28.63
26.59
25.41
23.24
25.61
23.61
26.38
24.08
30.18
28.09
28.19
24.94
0.68360.61080.7367
0.66330.62450.50370.76520.72070.67300.59400.64600.55190.76170.74650.75440.70320.59320.44340.62750.50690.68180.56800.85660.83020.84890.8005

Proposed (p = 1)
31.15
29.63
32.98
31.37
29.70
27.98
34.36
32.97
31.80
30.22
30.32
28.79
31.36
30.71
34.10
32.67
31.96
30.20
32.45
30.56
32.88
31.13
35.43
34.05
34.45
32.71
0.81050.76040.89020.85310.85770.79700.88820.85680.87330.82630.85630.81190.78060.75890.90190.86980.89000.84180.89370.84750.89500.85570.92220.90520.91780.8966

Proposed (p = 2)
30.71
29.18
33.04
31.44
29.29
27.59
34.34
32.95
31.61
30.04
30.15
28.63
31.04
30.46
34.03
32.60
31.83
30.08
32.21
30.38
32.70
30.96
35.51
34.10
34.14
32.36
0.79480.74270.89090.85340.84600.78460.88950.85780.87050.82280.85980.81630.76550.74910.90020.86790.88570.83640.88770.83920.89010.85080.92450.90800.92080.9002

KSC

18.78, 0.1267

29.30, 0.7660

29.33, 0.7742

28.34, 0.6814

ground-truth

observation

HTV

SSAHTVSSTV

28.63, 0.7544

34.10, 0.9019

34.03, 0.9002

ASSTV

proposed

proposed(p = 1)(p = 2)Reno

14.11, 0.07624

25.57, 0.6202

25.61, 0.6326

25.57, 0.6276

ground-truth

observation

HTV

SSAHTVSSTV

24.52, 0.5940

30.22, 0.8263

30.04, 0.8228

ASSTV

proposed

proposed(p = 1)(p = 2)

Figure 3.11: Resulting HS images with their PSNR[dB] (left) and SSIM (right) on the CS
reconstruction experiment under Gaussian noise situation (top: KSC, m = 0.4, bottom:
Reno, m = 0.2).

33

(a) band-wise PSNR(b) band-wise SSIM(c) Spatial response(d) Spectral response

Figure 3.12: Band-wise PSNR[dB] and SSIM and spatial and spectral responses on the
CS reconstruction experiment under Gaussian noise situation (Suwannee).
spatio-spectral piecewise-smoothness, SSTV reconstructs a similar spectral response to
the ground-truth one. However, the mean value is larger than the ground-truth value.
The proposed method achieves the most similar reconstruction of spectral response among
all the TVs.
To search the suitable ω in (3.6), we changed ω between 0.01 and 0.2 at 0.01 intervals.
Fig. 3.13 plots PSNR or SSIM of the results by the proposed method versus various ω,
where the values of PSNR and SSIM are averaged over the 13 HS images. One can seethat ω ∈ [0.05, 0.1] is a good choice in most cases. The results show that the suitable
range of ω in CS reconstruction is the same as Gaussian noise removal.

34

m = 0.4m = 0.2

Figure 3.13: PSNR[dB] (top) or SSIM (bottom) versus ω in (3.6) on the CS reconstruction
experiment under Gaussian noise situation.

3.2.2

CS Reconstruction under Mixed Noise Situation

Problem Formulation
To consider to CS reconstruction under mixed noise situation, we model the observationv ∈ RM as follows:v = Φū + n + s.
(3.19)
Based on the above model, we newly formulate compressed HS image reconstruction
considering mixed noise as the following convex optimization problem:v := {x ∈ RM |∥x − v∥ ≤ ε},Φu + s ∈ B2,εmin HSSTV(u) s.t.  s ∈ B1,η := {x ∈ RM |∥x∥1 ≤ η},(3.20)u,sNBu ∈ [µu , µu ] .
The first term in (3.20) is a regularization function for HS image restoration as definedin Sec. 3.1.2. The first and second constraints in (3.20) serve as data-fidelity regardingv
Gaussian and sparse noise to the observation v, respectively. The first constraint set B2,ε35

Algorithm 3: Primal-dual splitting method for Prob. (3.20)
(0)

(0)

(0)

(0)

input : u(0) , s(0) , y1 , y2 , y3 ,y4 , γ1 , γ2
1 while A stopping criterion is not satisfied do⊤2
u(n+1) = u(n) − γ1 (A⊤ω y1 + Φ y2 + y4 );
(n+1)
(n)
3
s
= s − γ1 (y2 + y3 );(n)(n)
y1 ← y1 + γ2 Aω (2u(n+1) − u(n) );
4
(n)
(n)
5
y2 ← y2 + γ2 (Φ(2u(n+1) − u(n) ) + (2s(n+1) − s(n) ));
(n)
(n)
6
y3 ← y3 + γ2 (2s(n+1) − s(n) );
(n)
(n)
7
y4 ← y4 + γ2 (2u(n+1) − u(n) );( (n) )
y1
(n+1)
(n)
8
y1
= y1 − γ2 prox 1 ∥·∥1,p;γ2
γ2( (n) )
y2
(n+1)
(n)
;
9
y2
= y2 − γ2 prox 1 ιBv
γ2
γ2
2,ε
( (n) )
y3
(n+1)
(n)
10
y3
= y3 − γ2 prox 1 ιB
;
γ2
γ2
1,η
( (n) )
y4
(n+1)
(n)
11
y4
= y4 − γ2 prox 1 [µu ,µu ]N B
;
γ2
γ2

12

n ← n + 1;

is the v-centered ℓ2 -norm ball with the radius ε, and the second one is the 0-centeredℓ1 -norm ball with the radius η. The third constraint in (3.20) is the dynamic range of u.
Optimization
In this part, we introduce an algorithm for solving Prob. (3.20). Since the problem isa highly nonsmooth convex optimization problem, we solve it by an iterative algorithm
based on a primal-dual splitting method [17]. In Sec. 2.3, the primal-dual splitting method
can solve convex optimization problems of the form of Prob. (2.9) by the algorithm (2.10),
and then we reformulate Prob. (3.20) into Prob. (2.9).
We put the constraints in Prob. (3.20) in the objective function by the indicator functions. Then, Prob. (3.20) can be rewritten as follows:v (Φu + s) + ιBmin ∥Aω u∥1,p + ιB2,ε1,η (s) + ι[µu ,µu ]N B (u).

u,s(3.21)

Note that Prob. (3.21) is equivalent to Prob. (3.20) (there is no approximation) becauseof the definition of the indicator function. By lettingg : RN B → R2 : (u, s) 7→ 0,v (z2 ) + ιB
h : R5N B+2M → R ∪ {∞} : (z1 , z2 , z3 , z4 ) 7→ ∥z1 ∥1,p + ιB2,ε1,η (z3 ) + ι[µu ,µu ]N B (z4 ),L : RN B → R5N B+2M : (u, s) 7→ (Aω u, Φu + s, s, u),Prob. (3.21) is reduced to Prob. (2.9). The resulting algorithm based on the primal-dual
splitting method is summarized in Alg. 3 with (2.11).
For the update of y1 , we use a simple soft-thresholding type operation (2.2) or (2.3).
The update of y2 , y3 , and y4 can be calculated by (2.4), (2.5), and (2.6), respectively.

36

Table 3.7: The four quality measures of the results on CS reconstruction under mixed
noise situation.
quality measures
methodsℓ2 fidelity
constraint only

Proposed(ℓ1 and ℓ2
fidelity constraints)

ground-truth

observation

HTV [6]SSTV [7]ASSTV [22]ℓ1 -HSSTVℓ1,2 -HSSTVHTV [6]SSTV [7]ASSTV [22]ℓ1 -HSSTVℓ1,2 -HSSTVHTVPSNR[dB]

SAM

ERGAS

Q2n

23.78
24.33
21.81
28.72
28.40
25.90
27.27
26.12
30.98
31.03

7.477
5.229
10.59
4.177
4.237
4.846
5.375
6.649
3.609
3.608

28.04
24.65
30.36
16.07
16.97
21.24
18.75
19.30
12.34
12.40

0.22940.39070.28290.56090.5656
0.36890.58780.38650.6726
0.6733

SSTV

ASSTV

ℓ1 -HSSTVℓ1,2 -HSSTV

Figure 3.14: Resulting HS images on CS reconstruction under mixed noise situation(Reno, top: ℓ2 fidelity constraint only, bottom: proposed).

Experiments
To verify the utility of the proposed framework, we conducted compressed HS reconstruction experiments under mixed Gaussian-sparse noise contamination. We examinedvarious regularizations: HTV, SSTV, ASSTV, and our HSSTV.
In the experiments, we artificially generated an observation based on (3.19) and reconstructed a desirable HS image from it. The sensing matrix Φ in (3.19) is assumed to bea random-sampling operator of ratio m = 0.2. We assumed that n is an additive white
Gaussian noise with the standard deviation σ = 0.05, and s is a salt-and-peeper noiseof ratio sp = 0.04. For test HS images, we took five HS images from the SpecTIR [44]and MultiSpec [45], cropped a region of size 256 × 256 × 32 for each HS image, and normalized
their dynamic range into [0, 1]. The radiuses ε and η in Prob. (3.20) were set to√mN B(1 − sp )σ 2 and 0.5mN Bsp , respectively.

37

We reconstruct a desirable HS image by solving Prob. 3.20 replacing HSSTV with HTV,
SSTV, and ASSTV and set the parameters in ASSTV and HSSTV as τb = 4 and ω = 0.04,
respectively. To verify the performance of the proposed framework, we compared ourframework with the method
considering only Gaussian noise, i.e., Prob. (3.14). In the√above case, we set ε = rN B((1 − sp )σ 2 + 0.25sp ). To evaluate the quality of restored
HS images, we used four quality measures: PSNR [dB], SAM [14], ERGAS [15], andQ2n [66]. We set the max iteration number and the stopping criterion of the primal-dual
splitting method to 10000 and ∥u(n) − u(n+1) ∥/∥u(n+1) ∥ < 1.0 × 10−4 , respectively.
Tab. 3.7 shows the values of the four quality measures of the results, which is the averageof five HS images. The left of the table represents the results by the method considering
only Gaussian noise, and the right is our framework. For most quality measures, one
can see that considering mixed Gaussian-sparse noise improves reconstruction accuracy.
Besides, HSSTV achieves the highest quality of all the regularizations.
In Fig 3.14, we show the reconstructed HS images (Reno) depicted as RGB images,
where we choose R = 8th, G = 16th, and B = 32nd bands. One can see that ourframeworks with ℓ1 /ℓ1,2 -HSSTV achieve the best performance.

38

Chapter 4

Image Fusion
4.1

Blur and Noisy Image Fusion

Image deblurring, removing blur from a given photograph, has been a fundamental andlongstanding problem in image processing and computer vision. Many image deblurring
methods, e.g., [67–75], are categolized as single image blind deblurring, that is, estimating
both the blur kernel and the latent image from a single blurred image (for more infomation
on single image blind deblurring, see a comprehensive survey [76]). Although single image
blind deblurring assumes the most realistic situation, it is a very challenging task due tothe highly under-costrained nature, so that it often leads to inaccurate estimation of blurkernels, high sensitivity to noise, and heavy dependence on prior information used.
To overcome the inherent diﬃculty, blind deblurring methods with a blurred/noisy
image pair have been studied [77–81]. These methods consider such a situation that both
images are captured in low light conditions with diﬀerent settings. Specifically, the blurred
image is taken with a slow shutter speed and a low ISO setting. With enough light, it has
the correct color and intensity, but it is blurry due to camera shake. On the other hands,
the noisy image is taken with a fast shutter speed and a high ISO setting. It is sharp butvery noisy because of insuﬃcient exposure and high camera gain. In addition, since it has
low contrast, the colors of this image are also partially lost. Under this situation, existing
methods take a two-step approach: first estimating the blur kernel from the image pairand then restoring a sharp image using the estimated kernel. Essentially, they can yielda better and much more stable result in their kernel estimation step than single image
deblurring methods, since the diﬀerence between the two images is extremely informative
for kernel estimation.
On the other hand, there exists a room for improvement in the image restoration stepof these methods. Specifically, the methods [77–79] estimate the latent image based only
on a given blurred image, i.e., do not exploit the information on a given noisy image intheir image restoration step, so that they are sensitive to the estimation error of the blurkernel. In addition, if noise in the blurred image is not negligible, restoring a sharp image
from it becomes diﬃcult even when using the true kernel. Meanwhile, the methods [80,81]
exploit the information on both images in their image restoration step. In these methods,
the image restoration problem is formulated as the minimization of a regularization term,

39

reflecting prior information on the latent image, plus two data-fidelity terms, keeping
the consistency to a blurred/noisy image pair, where the balance among these terms iscontroled by multiple weights. However, the tuning of such multiple weights is a tedioustask because they are interdependent and have no physical meaning. Indeed, suitablevalues of them vary depending on the latent image and/or the regularization terms used.
Based on the above discussion, we propose a new image restoration method of using
both a blurred image and a noisy image, which can be integrated into any blind deblurring methods with a blurred/noisy image pair. In our method, the image restoration
problem is formulated as a constrained convex optimization problem: minimizing a (possibly nonsmooth) regularization function subject to multiple hard constraints. Two of thehard constraints correspond to data-fidelity to a blurred image and that to a noisy image,respectively, where the degree of fidelity to each image can be controlled by independent
parameters that are explicitly related to the noise intensity of the image pair. We also
prove the existence of an optimal solution of the problem under reasonable assumptions.
Since our problem formulation properly incorporates the information on a blurred/noisy
image pair, it achieves (i) high quality restoration when the blurred image also contains
noise; and (ii) robustness to the estimation error of the blur kernel. At the same time, theindependence and clear meaning of the parameters thanks to the hard constraints oﬀer(iii) easy parameter setting.
Through several reformulations, we also provide an eﬃcient algorithm with guaranteed
convergence for solving the constrained convex optimization problem. Our algorithm isbased on ADMM a celebrated optimization method based on proximal splitting.

4.1.1

Problem Formulation

Consider to estimate an unknown latent color image ū ∈ R3N (3 is the number ofcolor channels) from an observed blurred image vn ∈ R3N and an observed noisy image
vb ∈ R3N . Specifically, following the prior work [80, 81] we model them as
vn = Bū + nn ,(4.1)

vb = ū + nb ,(4.2)

where B ∈ R3N ×3N is a blur operator estimated in advance, and nn and nb are additive
white Gaussian noises with their standard deviations σn and σb , respectively. The model
assumes that the blurred image vn also contains noise (usually σn < σb ), which is arealistic setting as addressed in [80, 82].
Based on the above observation models, we newly formulate a convex optimization
problem with multiple hard constraints for image restoration using a blurred/noisy image
pair as follows:u ∈ [0, 255]3N ,(4.3)
min R(Ψu) s.t.  Bu ∈ Bvn ,εn := {x ∈ R3N |∥x − vn ∥ ≤ εn },uu ∈ Bvb ,εb := {x ∈ R3N |∥x − vb ∥ ≤ εb }.
Here, R ◦ Ψ : RN → (−∞, ∞] is a regularization function (Ψ ∈ RL×3N , R ∈ Γ0 (RL )).
We assume that the proximity operator of R (NOT R ◦ Ψ) can be computed eﬃciently.

40

This assumption is essential in solving the problem by ADMM, as will be explained in
an optimization section. The first constraint set [0, 255]3N ⊂ R3N is the dynamic rangeof eight-bit color images. The second and third ones Bvn ,εn , Bvb ,εb ⊂ R3N are ℓ2 -normballs that represent data-fidelity to a blurred image vn and that to a noisy image vb ,respectively, where εn ≥ 0 and εb ≥ 0 are their radiuses determined based on the noise
intensity (noise standard deviations) of vn and vb .
Remark 1 (Design of regularization). TV [83] and its vectorial variants, e.g., [33, 84–87],
are well-known edge-preserving regularizers for images, and they have been used in many
deblurring methods. In this case, R is some norm, e.g., the ℓ1 norm, the mixed ℓ1,2 normor the nuclear norm, and Ψ is a discrete gradient operator. The proximity operators ofsuch norms are computable (see, e.g., [18, 88]). Another well-known example is frame
regularization relying on the sparsity of images in some transformed domain. In this case,
R is the ℓ1 norm, and Ψ is a frame analysis operator, e.g., wavelet [89] and curvelet [90].
More involved regularization, such as nonlocal regularization [91–93], regularization using
learned operators [94, 95] and plug-and-play regularization [96, 97], can also be handledin our formulation by setting Ψ to the corresponding nonlocal/learned analysis operator.
Remark 2 (Benefits of incorporating data-fidelity as hard constraints). Since εn and εb ,the radiuses of the ℓ2 -norm balls in (4.3), are directly related to the noise intensity (noise
standard deviations) of a blurred image and a noisy image, respectively, one can determine
their values with the help of existing noise level estimation methods. More importantly,
these parameters can be determined (almost) independent of the latent image ū and theregularization function R◦Ψ. This means that once finding suitable values of εn and εb forsome noise intensity, they can be used for various types of latent images and regularization
functions under the same noise intensity, which makes the setting of parameters on datafidelity much easier than existing methods that requires the tuning of interdependent
parameters (see the section of facilitation of parameter setting). Such benefits of hardconstraints have also been addressed, for example, in [10, 36, 37, 39, 98–100].
The following statement is on the existence of a solution of Prob. (4.3)
Proposition 1. Assume that the intersection of the constraint sets in (4.3) is nonempty,i.e.,
S := [0, 255]3N ∩ ΦBvn ,εn ∩ Bvb ,εb ̸= ∅,and that there exists some x ∈ S such that R(Ψx) < ∞. Then, Prob. (4.3) has at least
one optimal solution, i.e., the function R ◦ Ψ has a minimizer over S.
Proof : Since R ◦ Ψ ∈ Γ0 (R3N ) and S is a bounded closed convex subset of R3N , thestatement is a direct consequence of [101, Proposition 11.14].

4.1.2

Optimization

Since Prob. (4.3) is a highly nonsmooth constrained problem, we need suitable iterative
optimization methods to solve it. In this paper, we adopt ADMM, reviewed in Sec. 2.2. In
what follows, we reformulate Prob. (4.3) into the ADMM-applicable form, i.e., Prob. (2.7).

41

First, let us define the indicator functions (see Sec. 2.1) of the closed convex sets
[0, 255]3N , Bvn ,εn and Bvb ,εb . Then, Prob. (4.3) can be rewritten as
min R(Ψu) + ι[0,255]3N (u) + ιBvn ,εn (Bu) + ιBvb ,εb (u).

(4.4)

u

Second, we replace the input variables of all the terms in (4.4) with auxiliary variables
z1 , . . . , z4 , and express the relation between the input and the auxiliary variables by linear
equality constraints, yielding
min R(z1 ) + ι[0,255]3N (z2 ) + ιBvn ,εn (z3 ) + ιBvb ,εb (z4 )us.t. z1 = Ψu, z2 = u, z3 = Φu, z4 = u.

(4.5)

Third, we defineg(z1 , . . . , z4 ) := R(z1 ) + ι[0,255]3N (z2 ) + ιBvn ,εn (z3 ) + ιBvb ,εb (z4 ).
Then, the function g becomes proximable thanks to the variable splitting, as long as eachterm of g is proximable. Indeed, R is proximable from the assumption, and the otherterms, the three indicator functions, are also proximable because the metric projections
onto the corresponding closed convex sets are available (see Sec. 2.1).
Finally, by letting I be the identity matrix of size 3N × 3N and defining
 
Ψ
I

f (u) := 0 and G := 
(4.6)
B ,
I
Prob. (4.5) is reduced to Prob. (2.7).
The resulting algorithm based on ADMM is summarized in Alg. 4. Since G in (4.6) isa full column rank matrix due to I, G⊤ G is invertible, so that the convergence of Alg. 4
is guaranteed if a saddle point of g(z) − ⟨d, Gu − z⟩ exists.
Now we discuss the computation of each step of Alg 4. Since the update of u (Step 2) isstrictly-convex quadratic minimization because of the full-column-rankness of G, it boils
down to solving the following matrix inversion:u(n+1) =(Ψ⊤ Ψ + B⊤ B + 2I)−1 RHS
RHS :=(Ψ⊤ (z1 − d1 ) + (z2 − d2 ) + B⊤ (z3 − d3 ) + (z4 − d4 )).
(n)

(n)

(n)

(n)

(n)

(n)

(n)

(n)

Here, the above inverse matrix can be eﬃciently computed as with the CS reconstruction
cases.
For the updates of z1 , . . . , z4 (Step 3-6), we need to compute the proximity operatorsof each term of g. The proximity operator of R depends on the design of regularization.
As addressed in Remark 1, it is indeed computable for many types of regularization. The
vn , ι vb , and ι
proximity operator of ιB2,ε[0,255]3N equal to the metric projection onto each
B
n

2,εb

set given by (2.4) and (2.6).

42

Algorithm 4: ADMM method for Prob. (4.3)
(0)

(0)

(0)

(0)

(0)

(0)

(0)

(0)

input : z1 , z2 , z3 , z4 , d1 , d2 , d3 , d4 , and γ > 0
1 while A stopping criterion is not satisfied do
(n)
(n)
(n)
(n)
(n)
1
2
u(n+1) = argmin 2γ(∥z1 − Ψu(n) − d1 ∥2 + ∥z2 − u(n) − d2 ∥2 + ∥z3 − Φu(n) −
(n)

u
(n)

(n)

3

d3 ∥2 + ∥z4 − u(n) − d4 ∥2 );
(n)
(n+1)
z1
= proxγR (Ψu(n+1) + d1 );

4

z2

5

z3

6

z4

7
8
9
10
11

(n+1)

(n+1)

= proxγι

(n)[0,255]3N(n)= proxγιB(Φu(n+1) + d3 );= proxγιB(u(n+1) + d4 );

vn ,εn(n+1)

vb ,εb(n+1)(u(n+1) + d2 );

(n)

(n)

(n+1)

d1
= d1 + Ψu(n+1) − z1
;
(n+1)
(n)
(n+1)
d2
= d2 + u(n+1) − z2
;
(n+1)
(n)
(n+1)
d3
= d3 + Φu(n+1) − z3
;
(n+1)
(n)
(n+1)
(n+1)
d4
= d4 + u− z4;n ← n + 1;

Figure 4.1: Test images

4.1.3

Experiments

We demonstrate the three advantages of the proposed method: (i) high quality restoration when the blurred image also contains noise, (ii) robustness to the estimation errorof the blur kernel, and (iii) easy parameter setting. In the following experiments, weused 20 color images used as test images, which are taken from the Berkley Segmentation
Database [102] (Fig. 4.1). Then, we generated a pair of an observations based on (4.1)
and (4.2), and estimated a desirable image from them.
We utilized a popular color TV [33] as regularization in our method, which is defined
as follows:
TV(u) := ∥Du∥1,2 .
In this case, R := ∥ · ∥1,2 and Ψ := D in Prob. (4.3) The proximity operatorof the√2
mixed√ ℓ1,2 norm can be computed by (2.3). We set εn and εb in (4.3) to 0.95 3N σn and0.95 3N σb2 , respectively.
Experiments results of using both a blurred image and a noisy image are comparedwith that of using a blurred image only and that of using a noisy image only. At this

43

9
8
7
6
5
4
3
2
1
0
2

3

4

5

6

7

8

9

10

Figure 4.2: PSNR gain of the proposed method over the single image deblurring/denoising.

time, in the case of using a blurred image only, we solved that optimization problem:{u ∈ [0, 255]3N ,min TV(u) s.t.
uΦu ∈ Bvn ,εn := {x ∈ R3N |∥x − vn ∥ ≤ εn },and in the case of using a blurred image only, we solved that one:{u ∈ [0, 255]3N .
min TV(u) s.t.
uu ∈ Bvb ,εb := {x ∈ R3N |∥x − vb ∥ ≤ εb }.
We adopted PSNR[dB] to evaluate the objective quality of a restored image u, which
is given by
3N × 255
20 · log 10.
∥u − ū∥
Basic Performance Evaluation
To evaluate the eﬀectiveness of incorporating multiple hard constraints on data-fidelity
to a blurred/noisy image pair, we compare our method with two single image restoration

4429.77

24.06

21.74
24.06
latent image blurred image noisy image

33.12

32.09

35.09

27.06
single image
deblurring

27.58
single image
denoising

29.97
proposed

Figure 4.3: Resulting images with their PSNR[dB] in the first experiment (σb /σn = 8).

methods. One is single image deblurring, i.e., solving Prob. (4.3) without the constraintu ∈ Bvb ,εb . The other is single image denoising, i.e., solving Prob. (4.3) without theconstraint Φu ∈ Bvn ,εn . As in our method, ADMM was used for both methods.
In this experiment, we generated blurred images as follows: clean test images areblurred by a horizontal motion blur of 9 pixels and then contaminated by an additive
white Gaussian noise nn in (4.1) with the standard deviation σn = 2. Meanwhile, noisy
images were generated by adding a white Gaussian noise nb in (4.2) to clean test images,
where the standard deviation σb was increased from 4 to 20 by 2.
Fig. 4.2 plots the PSNR gain of our method over the single image deblurring (circle
marker) and the single image denoising (asterisk marker), where PSNR is averaged overthe 20 test images. One can see that for all the ratio of the noise standard deviations
σb /σn , the proposed method outperforms both single image deblurring and denoising
methods. This observation suggests that for a given blurred/noisy image pair, exploiting
information on the noisy image in the image restoration step is very eﬀective when theblurred image also contains noise.
Fig. 4.3 depicts some resulting images with their PSNR (σb /σn = 8). One can seethat 1. details are lost in the images obtained by the single image deblurring, 2. color
artifact remains in the images obtained by the single image denoising, and 3. Our method
achieves detail-preserving restoration with much less artifact.
We also check the convergence behavior of our algorithm (Alg. 4). For evaluationof convergence, we define the normalized root mean square error (NRMSE) between thecurrent estimate u(n) and the optimal solution u⋆ of Prob. (4.3), i.e., NRMSEn := ∥u(n) −u⋆ ∥/∥u⋆ ∥. Since the optimal solution u⋆ is analytically unavailable, it was pre-computed
by Alg. 4 with 100000 iterations. Fig. 4.4 plots the evolution of NRMSEn versus iterations(left) and the evolution of PSNR versus iterations (right), where the stepsize γ of ADMM
was set to 0.01. These plots suggest that Alg. 4 properly works, and that exploiting

45

35

10 0

30

10 -1

25
10

-2

20
10 -3
15
10

-4

10
10 -5

10 -6

5

0
50

100

150

200

250

300

50

100

150

200

250

300

Figure 4.4: Evolution of NRMSEn versus iterations (left) and the evolution of PSNR[dB]
versus iterations (right) of Alg. 4.

Table 4.1: Results of the experiment using inaccurate blur kernel.
kernel error
PSNR [dB]gain over single
image deblurring

0◦
32.97

1◦
32.91

2◦
32.62

3◦
32.29

4◦
32.01

5◦
31.79

6◦
31.61

7◦
31.47

8◦
31.36

9◦
31.26

10◦
31.17

spatially-varying
30.12

2.85

2.88

3.06

3.40

3.81

4.18

4.47

4.70

4.88

5.03

5.15

6.58

the information on a blurred/noisy image pair makes the convergence of ADMM fasterthan the single image deblurring/denoising, which is a positive side eﬀect of using a
blurred/noisy image pair.
Robustness to Inaccurate Blur Kernels
To illustrate the robustness of our method to the estimation error of blur kernels, we
conducted the following experiment. First, we generated images blurred by a certain blurkernel, which we refer to as the true blur kernel, and then in the image restoration step, weused an inaccurate blur kernel. Specifically, we consider the two cases: motion blur andspatially-varying blur. In the motion blur case, the true blur kernel was set to a horizontal
motion blur of 9 pixels, and the inaccurate blur kernel was set to a motion blur of 9 pixels
with its angle θ > 0, where we examined θ = 1◦ to 10◦ by 1◦ . In the spatially-varying
blur case, the true blur matrix Φ was made from spatially-varying per-pixel kernels, as
visualized in Fig. 4.6, and the inaccurate blur kernel was set to be spatially invariant withits kernel being the center kernel of the second image from left in Fig. 4.6. For both cases,
the blurred images contain an additive white Gaussian noise with the standard deviation
σn = 2, and the noisy images σb = 16.
Tab. 4.1 shows PSNR of restored images and the PSNR gain over the single image
deblurring, where these valuers are averaged over the 20 test images. One can observe thatfor the motion blur case, the PSNR gain over the single image deblurring is proportional46

blurred image

noisy image

18.12

24.06

single image
deblurring(θ = 1◦ )single image
deblurring(θ = 5◦ )single image
deblurring(θ = 10◦ )

26.96

22.92

proposed(θ = 1◦ )30.71

latent image

blurred image

noisy image

20.67

24.06

single image
deblurring(θ = 1◦ )single image
deblurring(θ = 5◦ )single image
deblurring(θ = 10◦ )

20.61

27.43

23.97

22.36

proposed(θ = 5◦ )proposed(θ = 10◦ )proposed(θ = 1◦ )proposed(θ = 5◦ )proposed(θ = 10◦ )

29.10

28.55

30.92

29.49

28.95

latent image

Figure 4.5: Restored results with their PSNR[dB] by using inaccurate blur kernels in the
motion blur case.

to the angle error, implying the robustness of our method to inaccurate blur kernels
compared with the single image deblurring. For the spartially-varying blur case, the
PSNR gain is also significant.
Fig. 4.5 and Fig. 4.6 depict several resulting images with their PSNR. One can see thatthe single image deblurring leads to oversmoothing when the kernel error is large. By
contrast, our method can restore sharp images in such a situation.
Facilitation of Parameter Setting
In the final experiment, we demonstrate that the setting of the parameters on datafidelity in our method are much easier than the existing methods using a blurred/noisy
image pair in the restoration step [80, 81]. In the existing methods, image restorationis performed by minimizing the sum of some specific regularization term and two datafidelity terms on a blurred/noisy image pair, which can be expressed as the following

47

21.17

latent image

blur kernel

24.06

21.58

29.08

21.11

24.06

20.73

28.44

blur image

noisy image

single image

proposeddeblurring

Figure 4.6: Restored results with their PSNR[dB] by using an inaccurate blur kernel inthe spatially-varying blur case.

Table 4.2: Best values of λ1 and λ2 in the existing formulation (Prob. (4.7)).
image
λ1
λ2
image
λ1
λ2

img1
733
8
img11
770
6

img2
1047
9
img12
895
8

img3
578
11
img13
988
8

img4
603
7
img14
932
8

img5
563
9
img15
905
7

img6
801
8
img16
812
6

img7
916
7
img17
917
8

img8
446
9
img18
678
8

img9
589
6
img19
920
8

img10
645
7
img20
650
8

optimization problem:
min R(Ψu) +u

λ2
λ1∥Φu − vn ∥2 + ∥u − vb ∥2 ,
2
2

(4.7)

where λ1 , λ2 > 0 control the balance among the three terms.
Tab. 4.2 shows the (hand-optimized) best values of λ1 and λ2 (in terms of PSNR) in
Prob. (4.7) for each test image, where the regularization term was set to the color TV.
For every test image, the noise standard deviations of a blurred/noisy image pair (σn , σb )were fixed at (2, 16). One can see that the best values of λ1 and λ2 are diﬀerent for eachtest image even though the noise standard deviations are the same for every test image.
This is because λ1 and λ2 depend on the latent image, so that the tuning of them isvery diﬃcult. By contrast, the proposed method achieves almost the√same restoration√
performance with common parameters εn = 0.95 3N σn2 and εb = 0.95 3N σb2 for all thetest images, i.e., the parameter setting is much easier.

48

4.2

Hyperspectral Image Fusion for Super-resolution

An HS image has rich spectral information, and HS imaging is a promising research
topic and oﬀers many applications in a wide range of fields [19, 20]. These applications
require an HS image of high spatial and spectral resolution (HR-HS image). However,since the amount of incident energy is limited, and there are critical tradeoﬀs between thespatial resolution and the spectral resolution of HS imaging systems, it is a very diﬃculttask to capture an HR-HS image.
HS image fusion techniques [103–105] try to resolve this dilemma, and have been actively studied [106–121]. They estimate an HR-HS image using a pair of an observed
HS image of high spectral resolution but low spatial resolution (LR-HS image) and an
observed guide image of high spatial resolution but low spectral resolution. In general,the guide image is assumed a panchromatic (PAN) image or an MS image, where a PAN
image has only 2D spatial information, i.e., a grayscale image.
Most recent HS image fusion methods [113–121] utilize apriori knowledge on an HSimage and observed information on an LR-HS image and a guide image, and estimate
HR-HS images by solving optimization problems. Thanks to the setting, they achieve
high-quality fusion. However, since these methods do not consider that a guide image
may also contain noise, they cannot achieve high-quality estimation, when the guide image
is noisy. As a result, the resulting HS images often have artifacts and spectral distortion.
To resolve the above problems, we propose a new robust HS image fusion method,which considers noisy guide images. The proposed method estimates not only an HR-HS
image but also a clean guide image, leading to high quality and robust estimation. The
method is built upon a convex optimization problem, where its objective function consistsof regularization terms for HS and guide images, respectively, and an edge similarity term
between HS and guide images. Data-fidelity to an LR-HS and a guide image and theirdynamic ranges are evaluated by hard constraints. This problem fully utilizes observed
information and apriori knowledge of an HS and a guide image, so that it can estimate an
HR-HS image without artifacts and spectral distortion even if both observed images arecontaminated by severe noise. To solve the optimization problem, we adopt a primal-dual
splitting method.

4.2.1

Problem Formulation
NB

In hyperspectral image fusion, an LR-HS image vH ∈ R r and an observed guide
image vg ∈ RM are assumed to be given with the observation model:
vH = SBū + nH ,(4.8)

vg = Rū + ng ,(4.9)

where S ∈ R r ×N B is a downsampling matrix with a downsampling rate of r (r is divisorof N ), B is a blur matrix, nH and ng are additive white Gaussian noises with standard
deviations σH and σg , respectively, and R ∈ RM ×N B is a matrix representing the spectral
response of the observed PAN image (R calculates weighted average along the spectral
direction). In general, since HS images contain more noise than guide images, we assumeNB49

σH > σg . This model says that both the LR-HS image and the observed guide image
contain considerable noise, which is a natural situation in hyperspectral imaging.
Based on the model, we formulate a hyperspectral image fusion problem as a convex
optimization problem. This problem estimates not only an HR-HS image uH ∈ RN B but
also a clean guide image ug ∈ RM , leading to high quality and robust estimation.
min HSSTV(uH ) + λ∥DuH − DMug ∥1,2 + ∥Dug ∥1,2
NB
vH
SBuH ∈ B2,ε:= {x ∈ R r |∥x − vH ∥ ≤ ε}, u ∈ B vg := {x ∈ RN |∥x − v ∥ ≤ η},
 g
g
2,η
s.t. 
NB
 uH ∈ [µuH , µuH ] ,ug ∈ [µug , µug ]M ,uH ,ug(4.10)

where the parameter λ > 0 is the parameter adjusting evaluation degree of the second
term, and M ∈ RN B×M is a linear operator that replicates the estimated guide image
along the spectral direction so that the number of bands in the guide image is the sameas that in the HS image.
The first and third terms in Prob. (4.10) are regularization functions for HS image andguide imege restoration, respectively. HSSTV has been shown to be very eﬀective in HS
image restoration, and we describe HSSTV in Sec. 3.1.2. By using HSSTV, the proposed
hyperspectral image fusion method can do robust estimation when the LR-HS image andthe observed guide image contain noise. The regularization for guide image restoration isvectorial TV.
The second term in Prob. (4.10) evaluates edge similarity between the HR-HS image uH
and the estimated guide image ug , which is originally proposed in [122]. Specifically, we
can assume that the non-zero diﬀerences of the HR-HS image are sparse and correspondto edges, and that their positions should be the same as those of the estimated guide
image. Hence, evaluating their errors by the mixed ℓ1,2 norm is a reasonable approachfor exploiting the spatial information on the estimated guide image.
The first constraint in (4.10) serves as data-fidelity to the LR-HS image vH and isdefined as the vH -centered ℓ2 -norm ball with the radius ε > 0. Likewise, the second
constraint in (4.10) plays data-fidelity role to the observed PAN image vg and is definedas the vg -centered ℓ2 -norm ball with the radius η > 0. The third and fourth constraintin (4.10) represent the dynamic range of a HS image and a guide image, respectively.

4.2.2

Optimization

Since Prob. (4.10) is a convex but highly constrained nonsmooth optimization problem, we require a suitable iterative algorithm, e.g., an alternating direction method ofmultipliers, to solve it. In this paper, we adopt a primal-dual splitting method [17] andwe explain it in Sec. 2.3. To solve Prob. (4.10) by the primal-dual splitting method, wereformulate Prob. (4.10) into Prob. (2.9).
We put the four constraints in Prob. (4.10) into the objective function. By the indicator

50

Algorithm 5: A primal-dual splitting method for Prob. (4.10).
(0)

(0)

(0)

(0)

(0)

(0)

input : u(0) , ug , y1 , y2 , y3 , y4 , y5
1 while A stopping criterion is not satisfied do do
(n)
(n)
(n)
2
u(n+1) = proxγ1 ι(u(n) − γ1 (A⊤+ D⊤ y2 + B⊤ S⊤ y4 ));
ω y1
NB
[µu

H

,µu

H

]

(n)

(n+1)
ug
= proxγ1 ι3[µug ,µug ]N

(n)

(ug

− γ1 (−M⊤ D⊤ y2(n)+ D⊤ y 3(n)(n)+ y5 ));(n)y1 ← y1 + γ2 Aω (2u(n+1) − u(n) );
(n)
(n)
(n+1)
(n)
y2 ← y2 + γ2 (D(2u(n+1) − u(n) ) − DM(2ug− ug ));
(n)
(n)
(n+1)
(n)
y3 ← y3 + γ2 D(2ug− ug );
(n)
(n)
(n+1)
y4 ← y4 + γ2 SB(2u− u(n) );
(n)
(n)
(n+1)
(n)
y5 ← y5 + γ2 (2ug
− ug );( (n) )
y1
(n+1)
(n)
y1
= y1 − γ2 prox 1 ,∥·∥1,p;γ2
γ2( (n) )
y2
(n+1)
(n)
y2
= y2 − γ2 prox λ ,∥·∥1,2;γ2γ2( (n) )
y3
(n+1)
(n)
y3
= y3 − γ2 prox 1 ,∥·∥1,2;γ2γ2( (n) )
y4
(n+1)
(n)
;
y4
= y4 − γ2 prox 1 ,Bv
γ2
2,ε
γ2
( (n) )
y5
(n)
(n+1)
= y5 − γ2 prox 1 ,Bvg
y5
;
γ2

4
5
6
7
8
9

10

11

12

13

γ2

2,η

n ← n + 1;

14functions, Prob. (4.10) can be rewritten asmin ∥Aω uH ∥1,p + λ∥DuH − DMug ∥1,2 + ∥Dug ∥1,2 + ιBvH (SBuH ) + ιBvg (ug )+uH ,ug2,ει[µu ,µu ]N B (uH ) + ι[µug ,µug ]M (ug ).
H2,η

(4.11)

H

Prob. (4.10) and Prob. (4.11) are equivalent from the definition of the indicator function.
Then, by lettingg :RN B+M → R2 : (uH , ug ) 7→ (ι[µu ,µu ]N B (uH ), ι[µug ,µug ]M (ug ))
H

5+ r1

h :R((

H)B+3)N +M → R ∪ {∞} : (y , y , y , y , y ) 7→
1

2

3

4

5

∥y1 ∥1,p + λ∥y2 ∥1,2 + ∥y3 ∥1,2 + ιBvH (y4 ) + ιBvg (y5 ),
2,ε

L :R

N (B+1)

2,η

1→ R((5+ r )B+3)N +M : (uH , ug ) 7→ (Aω uH , DuH − DMug , Dug , SBuH , ug ),Prob. (4.11) is reduced to Prob. (2.9). Using (2.11), the resulting algorithm for solving(4.10) is summarized in Alg. 5. Some steps in Alg. 5 require the proximity operators, andwe explain them in Sec 2.1.

4.2.3

Hyperspectral Pansharpening Experiments

In this section, we demonstrate the advantages of the proposed method over existing
HS pansharpening methods. For HS pansharpening, a guide image vg is a PAN image,51

Table 4.3: Quality measures for σg = 0.025 (top), σg = 0.05 (middle), and σ = 0.075(bottom).
σg

0.025

0.05

0.075

methodSFIM [112]
MTF-GLP [110]MTF-GLP-HPM [111]
GS [108]
GSA [109]
PCA [123]
GFPCA [107]
CNMF [116]
Bayesian Naive [113]
Bayesian Sparse [114]
HySure [115]proposed (p = 1)
proposed (p = 2)
SFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureproposed (p = 1)
proposed (p = 2)
SFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureproposed (p = 1)
proposed (p = 2)CC0.45280.69200.46050.59460.68410.59130.90190.88630.84980.85260.92730.95150.9516
0.41700.62840.42860.51080.62010.5086
0.88130.78390.7782
0.78300.87040.94090.94100.42590.58260.44290.43100.54590.42970.86940.6811
0.69200.70030.78680.93500.9351

52

SAM
38.87
34.68
38.89
39.77
41.71
39.93
11.18
15.10
27.20
26.68
15.93
9.777
9.763
39.45
35.47
39.54
41.04
44.77
41.21
11.54
16.23
30.88
30.34
20.46
9.891
9.878
39.59
36.03
39.57
42.00
48.48
42.17
11.71
17.69
35.00
34.35
25.34
9.919
9.907

RMSE
1571
974.4
1576
1101
1083
1111
462.1
512.1
602.2
594.1
402.9
322.4
322.2
1691
1112
1680
1213
1303
1221
500.6
729.0
800.9
785.4
557.1
344.4
344.2
1680
1210
1643
1311
1601
1317
520.2
951.5
1052
1023
774
356.7
356.5

ERGAS
23.86
16.05
23.80
20.54
20.11
20.72
8.045
8.338
11.07
10.95
7.017
5.672
5.666
25.32
17.76
25.11
22.39
23.81
22.53
8.762
11.90
14.26
14.03
9.714
6.135
6.130
25.63
18.92
25.25
24.02
28.89
24.13
9.138
15.48
18.38
17.93
13.52
6.387
6.382

CC

SAM

RMSE

ERGAS
λ

ω

Figure 4.7: Quality measures versus λ in (4.10) (top) / ω in (3.6) (bottom).

53

ground-truth

observed
HS imageobserved
PAN image

SFIM

MTF-GLP

MTF-GLP-HPM

GS

GSA

PCA

GFPCA

CNMF

Bayes Naive

Bayes Sparse

HySureℓ1 -HSSTVℓ1,2 -HSSTV

Figure 4.8: Resulting HS images on the HS pansharpening experiments (σH = 0.1, σg =0.05).
54

i.e., the number of pixel M = N . In this experiments, we generated a pair of an LR-HS
and an observed PAN image based on (4.8) and (4.9), estimated the HR-HS image fromthe pair. We used a Moﬀett field dataset as the true HR-HS image, which it is clippedin a region of size 256 × 128 × 176 and normalized its dynamic range into [0, 1], i.e.,µuH = µug = 0 and µuH = µug = 1 in Prob. (4.10). In (4.8) and (4.9), the downsampling
rate of S was set as r = 4, B was set to a 9 × 9 Gaussian blur matrix, and R was set toan weighted-average matrix with its weights wi (i = 1, . . . , B) were defined by{
1, if 1 ≤ i ≤ 41
wi =0, otherwise.
Then, we experimented with three pair of the standard deviations, (σH , σg ) = (0.1, 0.025),(0.1, 0.05), (0.1, 0.075). The above procedures follow Wald’s protocol [124], so that one
can see that it is a standard quality assessment methodology of HS pansharpening.
As compared methods, we utilize 11 existing methods: SFIM [112], MTF-GLP [110],
MTF-GLP-HPM [111], GS [108], GSA [109], PCA [123], GFPCA [107], CNMF [116],
Bayesian Naive [113], Bayesian Sparse [114] and HySure [115]. To set all parametersof these methods other than HySure, we used setting in a MATLAB toolbox of HS
pansharpening1 . For HySure, we set its hyperparameter as λvg hi = 0.1σH to enhance
its performance, and other parameters were set in the same way with other methods.
For our method, the parameters ε and η in (4.10) were set to oracle value, i.e., ε =∥v − SBū∥, η = ∥vg − Rū∥. Moreover, we varied λ ∈ [0.01, 0.1] and ω ∈ [0, 0.1] toinspect suitable them. We set the stepsizes, the max iteration number and the stopping
criterion of the primal-dual splitting method to γ1 = 0.005, γ2 = 1/1100γ1 , 5000 and∥u(n) − u(n+1) ∥/∥u(n) ∥ < 1.0 × 10−4 , respectively.
For the evaluation of the estimation performance, we adopt cross correlation (CC),the spectral angle mapper (SAM) [14], the root mean squared error (RMSE), and erreur
relative globale adimensionnelle de synthèse (ERGAS) [15] as quality measures, which aredefined as follows: for i = 1, . . . , N and j = 1, . . . , B,
∑N
1 ∑
(u
−αu,j )(ūi+(j−1)N −αū,j )√∑ i=1 i+(j−1)N
,
∑N
N
2
2
B
(u
−α
u,j )
i=1 i+(j−1)N
i=1 (ūi+(j−1)N −αū,j )B

CC(u, ū) =

j=1

1 ∑
arccosSAM(u, ū) =
N
N

i=1

(

u⊤
i ūi
∥ui ∥∥ūi ∥),∥u − ū∥RMSE(u, ū) = √
,
NB
v
u
B
∑
∥u∗j − ū∗j ∥2
100 u
u1
ERGAS(u, ū) =
(
)2′ ,
r tB
j=1 1 1⊤ u∗jprespectively, where ui = [ui , ui+N , . . . , ui+(B−1)N ] ∈ RB is the i-th spectral vector ofu, u∗j = [uN (j−1)+1 , uN (j−1)+2 , . . . , uN (j−1)+N ] ∈ RN is the j-th spatial one, αu,j =1

http://openremotesensing.net/

55

Table 4.4: Quality measures for σg = 0.04 (left) and σg = 0.06 (right) averaged all HS
images.
MTF-GLP [110]CNMF [116]
HySure [115]
MAPSMM [120]proposed(HTV, p = 1)
proposed(HTV, p = 2)
proposed(HSSTV, p = 1)
proposed(HSSTV, p = 2)∑Ni=1 ui+(j−1)N , αū,j =

PSNR
22.80
26.40
28.65
24.06

σm = 0.04
SAM ERGAS
10.45
7.557
9.139
6.104
8.354
7.511
10.39
14.00

Q2n
0.49370.62640.68360.4973

PSNR
22.28
24.25
27.93
23.49

σm = 0.08
SAM ERGAS
11.09
7.905
10.68
7.211
7.676
7.448
11.12
14.53

Q2n
0.46490.52730.65500.4658

29.79

5.103

4.459

0.7110

28.64

5.277

4.931

0.6790

29.89

5.020

4.386

0.7163

28.70

5.209

4.866

0.6832

29.86

5.774

5.062

0.7179

28.56

6.062

5.569

0.6829

29.95

5.711

5.006

0.7231

28.64

6.002

5.513

0.6880∑N

i=1 ūi+(j−1)N and 1 = [1, . . . , 1] ∈ R

N . Moreover, the closer CC

is 1 and the smaller SAM, RMSE and ERGAS are, the more alike the estimated HR-HS
image u and the true HR-HS image ū.
Tab. 4.3 shows CC, SAM, RMSE and ERGAS of the HR-HS images estimated by theexisting and proposed methods (p = 1 or 2 in (3.6)) for σg = 0.025, 0.05 and 0.075. For all
the quality measures and all standard deviations, one can see that the proposed method
outperforms all the existing methods.
Fig. 4.7 plots CC, SAM, RMSE and ERGAS of the HR-HS images estimated by the
proposed method versus λ in (4.10) and ω in (3.6), respectively, where we set ω = 0.01
in the λ graphs and λ = 0.03 in the ω graphs. In CC, RMSE and ERGAS case, we foundthat λ ∈ [0.02, 0.05] and ω ∈ [0, 0.02] are good choices, and HSSTV almost need not
to evaluate the direct spatial piecewise smoothness of an HS image in this experimental
setting. This is because the second term in Prob. (4.10) can evaluate it not just edge
similarity. For SAM case, λ ∈ [0.06, 0.08] and ω ∈ [0.05, 0.08] are good choices.
Fig. 4.8 is the estimated HR-HS images in the (σH , σg ) = (0.1, 0.05) case, which depictsas RGB images (R = 16th, G = 32nd and B = 64th bands). One can see that the results
estimated by most of the existing methods remain noise in the observed PAN image
and include artifacts. In addition, since the color in the results by GFPCA, CNMF andHySure is diﬀerent from that in the original HS image, it shows that these methods
produce spectral distortion. In contrast, the proposed method can estimate the HR-HS
image without noise, artifacts and spectral distortion, and it is most similar to the true
HR-HS image.

4.2.4

Hyperspectral and Multispectral Image Fusion Experiments

We demonstrate the advantages of the proposed method over existing HS and MS fusion
methods. In the experience, first, we generated both LR-HS and MS images using a test

56

ground-truth

Observed HS image

Observed MS image

MTF-GLPCNMF

HySure

MAPSMM(19.74, 12.63, 10.32, 0.3373) (22.97, 9.713, 8.825, 0.4808) (24.65, 8.670, 22.40, 0.4339) (21.64, 10.26, 29.24, 0.2890)

Proposed(HTV, p = 1)

Proposed(HTV, p = 2)

Proposed(HSSTV, p = 1)

Proposed(HSSTV, p = 2)

(25.34, 7.520, 8.349, 0.9072)(25.36, 7.500, 8.340, 0.9100)(24.55, 8.369, 9.480, 0.9020) (24.61, 8.316, 9.462, 0.9058)

Figure 4.9: Resulting HS images with four quality measure (PSNR[dB], SAM, ERGAS,
Q2n ) (DC, σg = 0.04).

57

HR-HS image based on (4.8) and (4.9), where nH and ng are additive white Gaussian
noises with standard deviation σH and σg , respectively. At this time, we adopted 15
HS images from the SpecTIR [44], MultiSpec [45], GIC [46], a Moﬀett field dataset andChikusei [125] as the test HS image, which are normalized its dynamic range into [0, 1].
In addition, the downsampling rate of S was set as r = 4, the kernel of B was a 9 × 9
Gaussian blur matrix, and R was a partial average matrix, which divides bands of an HSimage into eight parts, and averages it in each group.
Second, we estimated an HR-HS image from the pair by each method. In this experiments, we utilized MTF-GLP [110], CNMF [116], HySure [115] and MAPSMM [120] ascompared methods. Here, with respect to MTF-GLP and MAPSMM, all parameters were
set to the recommended values in a MATLAB toolbox of HS and MS fusion [105]. In thecase of HySure, we decided that its parameter regarding data fidelity λϕ = 0.5(σh + σm ),and for CNMF, the parameter was set to θh = σh and θm = σm for fair comparison. The
other parameters in HySure and CNMF were set to the values of the same toolbox. For
the proposed method, we decided that the parameters ε = ∥vH −SBūH ∥, η = ∥vg −Rūg ∥,ω = 0.02, and λ = 0.3. Besides, for verification of the availability of HSSTV regarding
an MS image denoising, we experimented by the method replacing the third term in
Prob. (4.10) with HSSTV, i.e., ∥Dug ∥1,2 is replaced HSSTV(ug ) with ω = 0.1.
Finally, we evaluated the estimated high HS images based on four standard quality
measures: PSNR[dB], SAM, ERGAS, and Q2n [66]. Note that the higher value of PSNR
and Q2n are, the more similar uh and ūh , and the lower value of SAM and ERGAS are,the more similar both images. The settings of max iteration number and the stopping
(n)
(n+1)
(n)
criterion of Alg. 5 were 5000 and ∥uh − uh∥/∥uh ∥ < 1.0 × 10−4 , respectively.
Tab. 4.4 shows PSNR[dB], SAM, ERGAS and Q2n of the estimated HS images by theexisting and proposed methods for (σH , σg ) = (0.1, 0.04) and (0.1, 0.06). At the time, itsvalues are the average of the results for all HS images. In Tab. 4.4, one can see that for
all the quality measures and for both noise intensity, the proposed method outperforms
all the existing methods. By comparison between the proposed by HTV and HSSTV,
in the almost PSNR, SAM and ERGAS cases, HTV is higher quality than HSSTV. We
think that since an MS image has weaker spectral correlation than an HS image, it hardly
necessary to evaluate the spatio-spectral piecewise smoothness of an MS image. The otherhands, in Q2n , the result of HSSTV and p = 2 is best. Since Q2n evaluates the structure
information and spectral/spatial distortion, the result of HSSTV have more correct detailthan that of HTV.
Fig. 4.9 is the estimated high HS images in the case of (σH , σg ) = (0.1, 0.04) and DC.
They are depicted as RGB images, where R, G, and B bands were set to the 30th, 60thand 90th bands of them in the case of the HS images, and these three bands are set as(R, G, B) = (2, 4, 6) bands in the MS image. One can see that the results by existing
methods have artifacts and spectral distortion. Specifically, the results using MTF-GLP,
MAPSMM, and CNMF have many artifacts, so that the spatial smoothness of them
are lost. In the HySure case, the result almost does not have artifacts. However, since
the color of the result changes compared with the original HS image, HySure produces
the spectral distortion. In contrast, the proposed methods avoid artifacts and spectral
distortion. In addition, HSSTV preserves more edge than HTV.

58

4.2.5Compressed Hyperspecteral Pansharpening

HS images are usually captured by dispersing the incident light, HS imaging requires
a 2D spectral sensor for 1D spatial information. This means that capturing an HS image
requires much time, and it is impossible to acquire it in one shot. To tackle the above
problem, the concept of CS plays a very essential role [49, 50], and one-shot HS imaging
based on CS is actively studied [52, 53].
HS imaging has another dilemma: capturing an HR-HS image itself is a challenging
task described in Sec. 4.2. HS image fusion resolves this dilemma, but existing HS image
fusion techniques [114–116,126–130] do not consider CS-based one-shot imaging scenarios.
Some MS image fusion methods involving the CS model have been proposed in [131–133].
However, since they assume that a given MS image only has 4 channel information (red,green, blue, and infrared radiation), utilizing the spectral information of MS images would
be less important than that of HS images.
Based on the above discussion, we bridge the gap between the said two approaches for
HS imaging. Specifically, we propose a novel HS imaging methodology, an HS pansharpening involving the CS model, named compressed HS pansharpening (CHPAN). The concept
enables to generate an HR-HS image from a compressed LR-HS image with the help of a
PAN image. For realistic setting, the concept assumes that both a CS observation and a
PAN image are contaminated by noise. Then, an HR-HS and a clean PAN images are simultaneously estimated from the noisy pair. This process is formulated as an optimization
problem involving an eﬀective spatio-spectral regularization [10, 11]. After reformulation,
the problem can be eﬃciently solved by a primal-dual splitting method [17], which is a
proximal splitting algorithm applied to a variety of imaging problems [85, 97, 134, 135].
Proposed Framework
In this methodology, a compressed LR-HS image vCHS ∈ RM is assumed to be given
by vCHS = ΦSBū + nCHS , where Φ ∈ M × rN B(M ≤ rN B) is a compressed sensing
matrix and nCHS is an additive white Gaussian noise with standard deviations σCHS .
We assume a PAN observation p ∈ RN like (4.9) with M = N .
To estimate an HR-HS image u and a clear PAN image q from the above two observations vCHS and p, we formulate a CHPAN as the following optimization problem:min HSSTV(u) + λ∥Du − DMq∥1,2 + ∥Dq∥1,2
u,q

vCHS
ΦSBu ∈ B2,ε:= {x ∈ RM |∥x − vCHS ∥ ≤ ε},p q ∈ B2,η := {x ∈ RN |∥x − p∥ ≤ η},
s.t. 
 u ∈ [µu , µu ]N B ,q ∈ [µq , µq ]N ,(4.12)

where M ∈ RN B×N is a linear operator that replicates a PAN image B times along thespectral direction. This problem is equals to Prob. (4.10) replacing SBu with ΦSBu in
the first constraint. Therefore, the problem can be solved by Alg. 5 replacing B⊤ S⊤ inthe step 2 of with B⊤ S⊤ Φ⊤ and SB in the step 7 with ΦSB.

59

Table 4.5: Performance comparison with four quality measures (left: m = 0.25 (CHPAN
case) and r = 4 (HS pansharpening case); right: m = 0.0625 (CHPAN case) and r = 16
(HS pansharpening case)).
quality measures
methods
HS pansharpening (p = 1)
HS pansharpening (p = 2)
CHPAN (p = 1)
CHPAN (p = 2)

PSNR
31.95
31.97
33.12
33.14

r = 4 and m = 0.25
SSIM
SAM
Q2n
0.85268.038
0.84170.8528 8.014 0.8419
0.88178.284
0.8566
0.8819 8.259 0.8568

r = 16 and m = 0.0625
PSNR
SSIM
SAM
Q2n
29.23
0.760610.22
0.761629.24
0.7608 10.20 0.7617
30.24
0.800810.77
0.790830.25 0.8011 10.75 0.7910

Experiments
We illustrate the utility of the proposed methodology and the estimation method. In
the experiment, we generate a pair of a CS observation and a noisy PAN image from atrue HR-HS image based on the observation model, and estimate an HR-HS image fromthe pair.
For the true HR-HS image, we used a Moﬀett field dataset, where we cropped a regionof size 256 × 128 × 176 and normalized its dynamic range into [0, 1]. The matrix R in
(4.9) was set to an weighted-average matrix, whose weights wi (i = 1, . . . , B) were defined
as wi = 1 for 1 ≤ i ≤ 41 and 0 otherwise. We set the standard deviations σCHS andσg of additive white Gaussian noises nCHS and ng as 0.1 and 0.04, respectively, and thesensing matrix Φ is a direct random-sampling matrix, where we examined two sampling
rates m = 0.25 and 0.0625.
We compared the results by CHPAN with that by HS pansharpening. We examined
two downsampling rates of S of HS pansharpening in (4.8), which are r = 4 and 16, andthen the HS pansharpening method with these rates were compared with the CHPAN
method with m = 0.25 and m = 0.0625, respectively. By this setting, the amount of
information in the observed HS image on the HS pansharpening experiment is equal tothat on the CHPAN experiment. Then, B was set to 5 × 5 and 9 × 9 Gaussian blur
matrixes in the case of r = 4 and 16.
The parameters ω in (3.6) and λ in (4.12) are set as 0.01 and 0.1, respectively. For our
CHPAN method, the parameters in (4.12) were set to ε = ∥vCHS − Φū∥ and η = ∥p − p̄∥,
where p̄ is a noise-less PAN image. We set the max iteration number and the stopping
criterion of the primal-dual splitting method to 5000 and ∥u(n) − u(n+1) ∥/∥u(n) ∥ < 1.0 ×
10−4 , respectively. The parameters in HS pansharpening were also set in the same way.
We use four quality measures: PSNR[dB], SSIM, SAM, and Q2n . The larger the valuesof PSNR, SSIM, and Q2n are, the more similar the estimated and true HS images are.
Likewise, the smaller the values of SAM is, the closer the two images are.
Tab. 4.5 shows PSNR[dB], SSIM, SAM, and Q2n of the estimated HS images by HS
pansharpening and CHPAN. For PSNR, SSIM, and Q2n , one can see that the CHPAN
method achieves to estimate better images than the HS pansharpening method. For SAM,
we will discuss later with Fig. 4.11.
Fig. 4.10 depicts the estimated HR-HS images as RGB images, where R, G, and B

60

ground-truth

HS pansharpening

HS pansharpening

HS pansharpening

observation(p = 1)(p = 2)CHPAN observation

CHPAN (p = 1)

CHPAN (p = 2)

Figure 4.10: Resulting HS images (m = 0.0625 (CHPAN case) and r = 16 (HS pansharpening case)).

61

HS pansharpening (r = 2)

HS pansharpening (r = 4)

CHPAN (m = 0.25)

CHPAN (m = 0.0625)

Figure 4.11: Spectral response of the results.
bands were set to the 8th, 16th and 32nd bands of them. The figure shows that theresults of the CHPAN method are clearer and closer to the ground-truth than those by
the HS pansharpening method.
In Fig. 4.11, we show the spectral responses of the results in Fig. 4.10. Specially, we
picked the spectral vectors of 58th row of 78th col from the results. One can see that
HS pansharpening changes the spectral variation between the 90 and 100 bands in the
ground-truth. On the other hand, CHPAN results in spectral variation similar to theoriginal one. As in Tab. 4.5, we found that the value of SAM of the result by CHPAN isslightly larger than that by HS pansharpening. Here, Tab. 4.5 shows that SAM value by
the result on CHPAN is slightly larger than that on HS pansharpening. SAM evaluates
the similarity of the spectral angles of two HS images. In other words, it ignores the error(distance) between the spectral values of the two images. This means that even if theerror between two spectral responses is large, SAM is small (smaller is better) as longas the shapes of the spectral responses of two HS images are similar. As can be seen in
Fig. 4.11, CHPAN estimates a more similar spectral response to the ground-truth thanHS pansharpening, despite the fact that CHPAN is slightly inferior to HS pansharpening62

in terms of SAM.

63

Chapter 5

Conclusion
In this paper, we have proposed new robust methods for image restoration and fusion
via constrained convex optimization. In Chap. 1, we have explained the main technique,image restoration. Images are aﬀected by many types of degradations and noises, so image
restoration methods require robustness. Existing image restoration methods based on
optimization achieve high-quality restoration, but the methods have several problems: notfully evaluating apriori knowledge on image and not considering outlier and/or modelingerror. To resolve the issue, our proposed methods oﬀer the following merits: (i) the fullevaluation of apriori knowledge on an image by a new regularization technique, (ii) theevaluation of the sparsity of outlier, and (iii) the eﬀective employment of multiple images
with diﬀerent characteristics.
In Chap. 2, we have introduced mathematical notation and abbreviation used in thispaper, and then we have explained proximal tools and algorithms solving a convex optimization problem, e.g., ADMM and the primal-dual splitting method.
We have introduced our proposed image restoration and fusion methods in Chap. 3
and Chap. 4, respectively. In Chap. 3, we focus on HS images, because they are easily
aﬀected by degradations and noises. For robust restoration, at first, we have proposeda new regularization technique, named HSSTV. Our regularization evaluates spatial andspectral piecewise smoothness on HS images as apriori knowledge. Next, we have formulated constrained optimization problems using HSSTV for HS image denoising and CS
reconstruction. The problems evaluate the sparsity of outliers, so the proposed methods
achieve robust restoration for outliers. Besides, data-fidelity and the sparsity of outlier
are evaluated by hard constraints, leading to intuitive parameter settings. Then, we have
solved them with ADMM after reformulation. In the simulation experiments, we haveillustrated the advantage of the proposed methods over existing methods for several image
restoration problem: Gaussian noise removal, Gaussian-sparse mixed noise removal, andCS reconstruction. Besides, we have confirmed the performance of the proposed methods
for real noise removal.
The approach (iii) is called image fusion, and we have proposed new robust image
fusion methods for two situations: blur and noisy image fusion and image fusion forsuper-resolution. The former estimates a desirable image from a pair of observed blurand noisy images, and the latter estimates an HR-HS image from an observed LR-HS

64

image and a guide image. Our proposed methods assume that all observations have noiseand estimate a desirable image by solving constrained optimization problems. Since our
proposed methods fully utilize information in both observations by evaluating the datafidelity and/or the edge similarity between them, the methods achieve high-quality robust
image fusion. We have demonstrated the performance of the proposed methods for image
fusion problems. Besides, the results show that the proposed methods are robust formodeling error because they can eﬀectively utilize information in both observations.
In the above results and discussion, for image restoration and fusion, one can see thatour proposed methods are robust for noises, degradations, and modeling error. In thiswork, we have mainly conducted simulation experiments, and the degradation and sensing
operators are ideally set. Therefore, the performance of the methods should be verified fora real situation. In the future, we intend to confirm the performance using real datasets,
real estimated operators, and a real imaging system.

65

Chapter 6

Acknowledgement
I would like to thank Professor Itsuo Kumazawa and Associate Professor Shunsuke
Ono for kind supervision. I would like to express my heartfelt gratitude to all the peoplein Kumazawa laboratory, Mathematical Data Informatics Laboratory, and Augmented
Vision Laboratory who spent time with me during my student life and gave me various
supports. I would like also to deeply thank Professor Masahiro Yamaguchi, Professor Isao
Yamada, and Associate Professor Takahiro Shinozaki, Tokyo Institute of Technology,
Japan, for serving as members of the examining committee of this dissertation. I am
grateful to my family and friends for mentally supporting me.
This work was supported by Japan Society for the Promotion of Science (JSPS), Research Fellow (DC1) grant.

66

References[1] D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear variational
problems via finite elements approximations,” Comput. Math. Appl., vol. 2, pp. 17–40, 1976.
[2] J. Eckstein and D.P. Bertsekas, “On the Douglas-Rachford splitting method andproximal point algorithm for maximal monotone operators,” Math. Program., vol.
55, pp. 293–318, 1992.
[3] J. Eckstein and D. P. Bertsekas, “On the Douglas―Rachford splitting method andthe proximal point algorithm for maximal monotone operators,” Math. Program.,vol. 55, no. 1-3, pp. 293–318, 1992.
[4] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,”
Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
[5] H. Zhang, W. He, L. Zhang, H. Shen, and Q. Yuan, “Hyperspectral image restoration using low-rank matrix recovery,” IEEE Trans. on Geosci. and Remote Sens.,
vol. 52, no. 8, pp. 4729–4743, 2014.
[6] Q. Yuan, L. Zhang, and H. Shen, “Hyperspectral image denoising employing aspectral–spatial adaptive total variation model,” IEEE Trans. on Geosci. and Remote Sens., vol. 50, no. 10, pp. 3660–3677, 2012.
[7] H. K. Aggarwal and A. Majumdar, “Hyperspectral image denoising using spatiospectral total variation,” IEEE Geosci. and Remote Sens. Lett., vol. 13, no. 3, pp.
442–446, 2016.
[8] S. H. Chan, R. Khoshabeh, K. B. Gibson, P. E. Gill, and T. Q. Nguyen, “Anaugmented lagrangian method for total variation video restoration,” IEEE Trans.
Image Process., vol. 20, no. 11, pp. 3097–3111, 2011.
[9] H. Liu, P. Sun, Q. Du, Z. Wu, and Z. Wei, “Hyperspectral image restoration
based on low-rank recovery with a local neighborhood weighted spectral-spatial
total variation model,” IEEE Trans. on Geosci. and Remote Sens., vol. 57, no. 3,
pp. 1–14, 2018.

67

[10] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral image restoration by
hybrid spatio-spectral total variation,” Proc. IEEE Int. Conf. Acoust., Speech,
Signal Process. (ICASSP), pp. 4586–4590, 2017.
[11] S. Takeyama, S. Ono, and I. Kumazawa, “A constrained convex optimization approach to hyperspectral image restoration with hybrid spatio-spectral regularization,” Remote Sensing, vol. 12, no. 21, pp. 3541, 2020.
[12] W. He, H. Zhang, L. Zhang, and H. Shen, “Total-variation-regularized low-rankmatrix factorization for hyperspectral image restoration,” IEEE Trans. on Geosci.
and Remote Sens., vol. 54, no. 1, pp. 178–188, 2016.
[13] Q. Wang, Z. Wu, J. Jin, T. Wang, and Y. Shen, “Low rank constraint and spatial
spectral total variation for hyperspectral image mixed denoising,” Signal Processing,vol. 142, pp. 11–26, 2018.
[14] F.A. Kruse, A.B. Lefkoﬀ, J.W. Boardman, K.B. Heidebrecht, A.T. Shapiro, P.J.
Barloon, and A.F.H. Goetz, “The spectral image processing system (SIPS)―
interactive visualization and analysis of imaging spectrometer data,” Remote Sensing of Environment, vol. 44, no. 2, pp. 145–163, 1993.
[15] L. Wald, “Quality of high resolution synthesised images: Is there a simple criterion?,” in Third conference” Fusion of Earth data: merging point measurements,
raster maps and remotely sensed images”, 2000, pp. 99–103.
[16] J. J. Moreau, “Fonctions convexes duales et points proximaux dans un espacehilbertien,” C. R. Acad. Sci. Paris Ser. A Math., vol. 255, pp. 2897–2899, 1962.
[17] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imaging and Vision, vol. 40, no. 1,
pp. 120–145, 2010.
[18] P. L. Combettes and J.-C. Pesquet, “Proximal splitting methods in signal processing,” in Fixed-Point Algorithms for Inverse Problems in Science and Engineering,pp. 185–212. Springer-Verlag, 2011.
[19] C. I. Chang, Hyperspectral imaging: techniques for spectral detection and classification, vol. 1, Springer Science & Business Media, 2003.
[20] A. Plaza et al., “Recent advances in techniques for hyperspectral image processing,”
Remote sensing of environment, vol. 113, pp. S110–S122, 2009.
[21] Behnood Rasti, Paul Scheunders, Pedram Ghamisi, Giorgio Licciardi, and Jocelyn
Chanussot, “Noise reduction in hyperspectral imagery: Overview and application,”
Remote Sens., vol. 10, no. 3, pp. 482, 2018.
[22] Y. Chang, L. Yan, H. Fang, and C. Luo, “Anisotropic spectral-spatial total variation
model for multispectral remote sensing image destriping,” IEEE Trans. Image
Process., vol. 24, no. 6, pp. 1852–1866, 2015.

68

[23] L. Sun, C. Ma, Y. Chen, Y. Zheng, H.J. Shim, Z. Wu, and B. Jeon, “Low rankcomponent induced spatial-spectral kernel method for hyperspectral image classification,” IEEE Trans. Circuits Syst. Video Technol., pp. 4133–4148, 2019.
[24] H. Li, P. Sun, H. Liu, Z. Wu, and Z. Wei, “Non-convex low-rank approximation
for hyperspectral image recovery with weighted total varaition regularization,” inProc. IEEE Int. Geosci. and Remote Sens. Symposium (IGARSS), 07 2018, pp.
2733–2736.
[25] H. K. Aggarwal and A. Majumdar, “Hyperspectral unmixing in the presence of
mixed noise using joint-sparsity and total variation,” IEEE Sel. Top. Appl. Earth
Obs. Remote Sens., vol. 9, no. 9, pp. 4257–4266, 2016.
[26] L. Sun, F. Wu, T. Zhan, W. Liu, J. Wang, and B. Jeon, “Weighted nonlocal lowrank tensor decomposition method for sparse unmixing of hyperspectral images.,”IEEE Sel. Top. Appl. Earth Obs. Remote Sens., pp. 1174–1188, 2020.
[27] W. He, H. Zhang, H. Shen, and L. Zhang, “Hyperspectral image denoising using
local low-rank matrix recovery and global spatial–spectral total variation,” IEEE
Sel. Top. Appl. Earth Obs. Remote Sens., vol. 11, no. 3, pp. 713–729, 2018.
[28] W. Cao, K. Wang, G. Han, J. Yao, and A. Cichocki, “A robust pca approach withnoise structure learning and spatial–spectral low-rank modeling for hyperspectral
image restoration,” IEEE Sel. Top. Appl. Earth Obs. Remote Sens., vol. 11, no. 10,
pp. 3863–3879, 2018.
[29] Y. Wang, J. Peng, Q. Zhao, Y. Leung, X. Zhao, and D. Meng, “Hyperspectral
image restoration via total variation regularized low-rank tensor decomposition,”
IEEE Sel. Top. Appl. Earth Obs. Remote Sens., vol. 11, no. 4, pp. 1227–1243, 2017.
[30] L. Sun, T. Zhan, Z. Wu, L. Xiao, and B. Jeon, “Hyperspectral mixed denoising via
spectral diﬀerence-induced total variation and low-rank approximation,” Remote
Sens., vol. 10, no. 12, pp. 1956, 2018.
[31] T. Ince, “Hyperspectral image denoising using group low-rank and spatial-spectral
total variation,” IEEE Access, vol. 7, pp. 52095–52109, 2019.
[32] X. Kong, Y. Zhao, J. Xue, J. CW. Chan, Z. Ren, HX. Huang, and J. Zang, “Hyperspectral image denoising based on nonlocal low-rank and tv regularization,” Remote
Sensing, vol. 12, no. 12, pp. 1956, 2020.
[33] X. Bresson and T. F. Chan, “Fast dual minimization of the vectorial total variation
norm and applications to color image processing,” Inverse Probl. Imag., vol. 2, no.
4, pp. 455–484, 2008.
[34] L. Sun, C. He, Y. Zheng, and S. Tang, “Slrl4d: Joint restoration of subspace
low-rank learning and non-local 4-d transform filtering for hyperspectral image.,”
Remote Sens., p. 2979, 2020.

69

[35] M. Afonso, J. Bioucas-Dias, and M. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,”
IEEE Trans. Image Process., vol. 20, no. 3, pp. 681–695, 2011.
[36] G. Chierchia, N. Pustelnik, J.-C. Pesquet, and B. Pesquet-Popescu, “Epigraphicalprojection and proximal tools for solving constrained convex optimization problems,” Signal, Image and Video Process., pp. 1–13, 2014.
[37] S. Ono and I. Yamada, “Signal recovery with certain involved convex data-fidelity
constraints,” IEEE Trans. Signal Process., vol. 63, no. 22, pp. 6149–6163, 2015.
[38] Y. Xie, Y. Qu, D. Tao, W. Wu, Q. Yuan, and W. Zhang, “Hyperspectral image restoration via iteratively regularized weighted schatten p-norm minimization,”
IEEE Trans. on Geosci. and Remote Sens., vol. 54, no. 8, pp. 4642–4659, 2016.
[39] S. Ono, “L0 gradient projection,” IEEE Trans. Image Process., vol. 26, no. 4, pp.
1554–1564, 2017.
[40] S. Takeyama, S. Ono, and I. Kumazawa, “Robust and eﬀective hyperspectral pansharpening using spatio-spectral total variation,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP), 2018, pp. 1603–1607.
[41] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral pansharpening using noisy
panchromatic image,” in Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), 2018, pp. 880–885.
[42] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral and multispectral data
fusion by a regularization considering,” in Proc. IEEE Int. Conf. Acoust., Speech,
Signal Process. (ICASSP), 2019, pp. 2152–2156.
[43] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process.,
vol. 13, no. 4, pp. 600–612, 2004.
[44] “Spectir,” http://www.spectir.com/free-data-samples/.
[45] “Multispec,” https://engineering.purdue.edu/ biehl/MultiSpec.
[46] “GIC,” http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral Remote Sensing Scenes.
[47] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, “Nonlocal transform-domain
filter for volumetric data denoising and reconstruction,” IEEE Trans. Image Process., vol. 22, no. 1, pp. 119–133, 2013.
[48] W. Liu and J. Lee, “A 3-d atrous convolution neural network for hyperspectral
image denoising,” IEEE Trans. on Geosci. and Remote Sens., 2019.
[49] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory, vol. 52, no. 4,
pp. 1289–1306, 2006.

70

[50] E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans.
Inform. Theory, vol. 52, no. 2, pp. 489–509, 2006.
[51] J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Magazine,
vol. 25, no. 2, pp. 14–20, 2008.
[52] R. M. Willett, M. F. Duarte, M. A. Davenport, and R. G. Baraniuk, “Sparsity andstructure in hyperspectral imaging: Sensing, reconstruction, and target detection,”
IEEE Signal Process. Magazine, vol. 31, no. 1, pp. 116–126, 2014.
[53] G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded
aperture spectral imaging: An introduction,” IEEE Signal Process. Magazine, vol.
31, no. 1, pp. 105–115, 2014.
[54] M. Golbabaee and P. Vandergheynst, “Hyperspectral image compressed sensing
via low-rank and joint-sparse matrix recovery,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP), 2012, pp. 2741–2744.
[55] S. Yang, M. Wang, P. Li, L. Jin, B. Wu, and L. Jiao, “Compressive hyperspectral
imaging via sparse tensor and nonlinear compressed sensing,” IEEE Trans. on
Geosci. and Remote Sens., vol. 53, no. 11, pp. 5943–5957, 2015.
[56] J. Wei, Y. Huang, K. Lu, and L. Wang, “Nonlocal low-rank-based compressedsensing for remote sensing image reconstruction,” IEEE Geosci. and Remote Sens.
Lett., vol. 13, no. 10, pp. 1557–1561, 2016.
[57] H. Garcia, C. V. Correa, and H. Arguello, “Multi-resolution compressive spectral imaging reconstruction from single pixel measurements,” IEEE Trans. Image
Process., vol. 27, no. 12, pp. 6174–6184, 2018.
[58] P. C. Hansen, J. G. Nagy, and D. P. O Leary, Deblurring Images: Matrices, Spectra,
and Filtering, SIAM, 2006.
[59] G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins University
Press, 4th edition, 2012.
[60] P. L. Combettes and J.-C. Pesquet, “Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone
operators,” Set-Valued and Variational Analysis, vol. 20, no. 2, pp. 307–330, 2012.
[61] L. Condat, “A primal-dual splitting method for convex optimization involving
Lipschitzian, proximable and linear composite terms,” J. Optimization Theory andApplications, 2013.
[62] S. Ono, M. Yamagishi, T. Miyata, and I. Kumazawa, “Image restoration using astochastic variant of the alternating direction method of multipliers,” in Proc. IEEE
Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2016.

71

[63] A. Chambolle, M. J. Ehrhardt, P. Richtárik, and C.-B. Schonlieb, “Stochastic
primal-dual hybrid gradient algorithm with arbitrary sampling and imaging applications,” SIAM J. Optim., vol. 28, no. 4, pp. 2783–2808, 2018.
[64] P. L. Combettes and J.-C. Pesquet, “Stochastic forward-backward and primal-dual
approximation algorithms with application to online image restoration,” in Proc.
Eur. Signal Process. Conf. (EUSIPCO), Aug 2016, pp. 1813–1817.
[65] S. Ono, “Eﬃcient constrained signal reconstruction by randomized epigraphicalprojection,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP),
2019, pp. 4993–4997.
[66] A. Garzelli and F. Nencini,“Hypercomplex quality assessment ofmulti/hyperspectral images,” IEEE Geosci. and Remote Sens. Lett., vol. 6,
no. 4, pp. 662–665, 2009.
[67] Qi Shan, Jiaya Jia, and Aseem Agarwala, “High-quality motion deblurring from asingle image,” ACM Transactions on Graphics, vol. 27, no. 3, pp. 73:1–73:10, 2008.
[68] Sunghyun Cho and Seungyong Lee, “Fast motion deblurring,” ACM Transactions
on Graphics, vol. 28, no. 5, pp. 145:1–145:8, 2009.
[69] Jian-Feng Cai, Hui Ji, Chaoqiang Liu, and Zuowei Shen, “Blind motion deblurringfrom a single image using sparse approximation,” in Proc. IEEE Conf. Comput.
Vis. Pattern Recognit. (CVPR), 2009, pp. 104–111.
[70] Li Xu and Jiaya Jia, “Two-phase kernel estimation for robust motion deblurring,”
in Proc. Eur. Conf. Comput. Vis. (ECCV), 2010, pp. 157–170.
[71] Dilip Krishnan, Terence Tay, and Rob Fergus, “Blind deconvolution using a normalized sparsity measure,” in Proc. IEEE Conf. Comput. Vis. Pattern Recognit.
(CVPR), 2011, pp. 233–240.
[72] Jian-Feng Cai, Hui Ji, Chaoqiang Liu, and Zuowei Shen, “Framelet-based blindmotion deblurring from a single image,” IEEE Trans. Image Process., vol. 21, no.
2, pp. 562–572, 2012.
[73] Lin Zhong, Sunghyun Cho, Dimitris Metaxas, Sylvain Paris, and Jue Wang, “Handling noise in single image deblurring using directional filters,” in Proc. IEEE Conf.
Comput. Vis. Pattern Recognit. (CVPR), 2013, pp. 612–619.
[74] Tomer Michaeli and Michal Irani, “Blind deblurring using internal patch recurrence,” in Proc. Eur. Conf. Comput. Vis. (ECCV), 2014, pp. 783–798.
[75] Daniele Perrone and Paolo Favaro, “A clearer picture of total variation blind deconvolution,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 38, no. 6, pp. 1041–1055,
2016.

72

[76] Wei-Sheng Lai, Jia-Bin Huang, Zhe Hu, Narendra Ahuja, and Ming-Hsuan Yang,“A comparative study for single image blind deblurring,” in Proc. IEEE Conf.
Comput. Vis. Pattern Recognit. (CVPR), 2016, pp. 1701–1709.
[77] L. Yuan, J. Sun, L. Quan, and H. Y. Shum, “Image deblurring with blurred/noisy
image pairs,” ACM Transactions on Graphics, vol. 26, no. 3, pp. 1:1–1:10, 2007.
[78] S. Lee, H.M. Park, and S.Y. Hwang, “Motion deblurring using edge map withblurred/noisy image pairs,” Optics Communications, vol. 285, no. 7, pp. 1777–
1786, 2012.
[79] C. Je, H. S. Jeon, C.H. Son, and H.M Park, “Disparity-based space-variant image
deblurring,” Signal Process.: Image Commun., vol. 28, no. 7, pp. 792–808, 2013.
[80] H. Li, Y. Zhang, J. Sun, and D. Gong, “Joint motion deblurring with blurred/noisy
image pair,” in Proc. Int. Conf. Pattern Recognit. (ICPR), 2014, pp. 1020–1024.
[81] C. H. Son and X. P. Zhang, “Layer-based approach for image pair fusion,” IEEE
Trans. Image Process., vol. 25, no. 6, pp. 2866–2881, 2016.
[82] L. Zhong, S. Cho, D. Metaxas, S. Paris, and J. Wang, “Handling noise in single image deblurring using directional filters,” in Proc. IEEE Conf. Comput. Vis. Pattern
Recognit. (CVPR), 2013, pp. 612–619.
[83] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal
algorithms,” Phys. D, vol. 60, no. 1-4, pp. 259–268, 1992.
[84] B. Goldluecke, E. Strekalovskiy, and D. Cremers, “The natural vectorial total
variation which arises from geometric measure theory,” SIAM J. Imag. Sci., vol. 5,
no. 2, pp. 537–563, 2012.
[85] S. Ono and I. Yamada, “Decorrelated vectorial total variation,” in Proc. IEEE
Conf. Comput. Vis. Pattern Recognit. (CVPR), 2014.
[86] S. Lefkimmiatis and A. Roussos.and P. Maragos.and M. Unser, “Structure tensor
total variation,” SIAM J. Imag. Sci., vol. 8, no. 2, pp. 1090–1122, 2015.
[87] S. Ono, K. Shirai, and M. Okuda, “Vectorial total variation based on arranged
structure tensor for multichannel image restoration,” in Proc. IEEE Int. Conf.
Acoust., Speech, Signal Process. (ICASSP), 2016.
[88] S. Ono, “Distributed convex optimization via proximal splitting: A survey on
admm-based approaches,” Journal of The Society of Instrument and Control Engineers, vol. 55, no. 11, pp. 954–959, 2016, (in Japanese).
[89] S. Mallat, A wavelet tour of signal processing, Academic Press, 2nd edition, 1999.
[90] E. Candès, L. Demanet, D. L. Donoho, and L. Ying, “Fast discrete curvelet transforms,” SIAM J. Multi. Model. Simul., vol. 5, no. 3, pp. 861–899, 2006.

73

[91] G. Gilboa and S. Osher, “Nonlocal linear image regularization and supervised
segmentation,” Multiscale Model. Simul., vol. 6, no. 2, pp. 595–630, 2007.
[92] A. Danielyan, V. Katkovnik, and K. Egiazarian, “BM3D frames and variational
image deblurring,” IEEE Trans. Image Process., vol. 21, no. 4, pp. 1715–1728,
2012.
[93] G. Chierchia, N. Pustelnik, B. Pesquet-Popescu, and J.-C. Pesquet, “A nonlocal structure tensor-based approach for multicomponent image recovery problems,”
IEEE Trans. Image Process., vol. 23, no. 12, pp. 5531–5544, 2014.
[94] J. Mairal, F. Bach, J. Ponce, G. Sapiro, and A. Zisserman, “Non-local sparse models
for image restoration,” in Proc. IEEE Int. Conf. Comput. Vis. (ICCV), 2009.
[95] D. Zoran and Y. Weiss, “From learning models of natural image patches to whole
image restoration,” in Proc. IEEE Int. Conf. Comput. Vis. (ICCV), 2011.
[96] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, G. T. Buzzard, L. F. Drummy,
J. P. Simmons, and C. A. Bouman, “Plug-and-play priors for bright field electron
tomography and sparse interpolation,” IEEE Trans. Comput. Imag., vol. 2, no. 4,
pp. 408–423, 2016.
[97] S. Ono, “Primal-dual plug-and-play image restoration,” IEEE Signal Process. Lett.,
vol. 24, no. 8, pp. 1108–1112, 2017.
[98] M. V. Afonso, J. M. Bioucas-Dias, and M. AT. Figueiredo, “An augmented lagrangian approach to the constrained optimization formulation of imaging inverse
problems,” IEEE Trans. Image Process., vol. 20, no. 3, pp. 681–695, 2011.
[99] M. Carlavan and L. Blanc-Fèraud, “Sparse Poisson noisy image deblurring,” IEEE
Trans. Image Process., vol. 21, no. 4, pp. 1834–1846, 2012.
[100] T. Teuber, G. Steidl, and R. H. Chan, “Minimization and parameter estimation for
seminorm regularization models with I-divergence constraints,” Inverse Problems,
vol. 29, no. 3, 2013.
[101] H. H. Bauschke and P. L. Combettes, Convex analysis and monotone operator
theory in Hilbert spaces, Springer, New York, 2011.
[102] D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuringecological statistics,” in Proc. IEEE Int. Conf. Comput. Vis. (ICCV), 2001.
[103] L. Loncan, L. B. de Almeida, J. M. Bioucas-Dias, X. Briottet, J. Chanussot, N. Dobigeon, S. Fabre, W. Liao, G. A. Licciardi, M. Simoes, et al., “Hyperspectral
pansharpening: A review,” IEEE Geoscience and remote sensing magazine, vol. 3,
no. 3, pp. 27–46, 2015.

74[104] N. Yokoya, C. Grohnfeldt, and J. Chanussot, “Hyperspectral and multispectral
data fusion: A comparative review of the recent literature,” IEEE Geoscience andRemote Sensing Magazine, vol. 5, no. 2, pp. 29–56, 2017.
[105] N. Yokoya, C. Grohnfeldt, and J. Chanussot, “Hyperspectral and multispectral data
fusion: A comparative review of the recent literature,” IEEE Geosci. and Remote
Sens. Magazine, vol. 5, no. 2, pp. 29–56, 2017.
[106] P. Chavez, S. C. Sides, J. A. Anderson, et al., “Comparison of three diﬀerent
methods to merge multiresolution and multispectral data- landsat tm and spot
panchromatic,” Photogrammetric Engineering and remote sensing, vol. 57, no. 3,
pp. 295–303, 1991.
[107] W. Liao, X. Huang, F. Van Coillie, S. Gautama, et al., “Processing of multiresolution thermal hyperspectral and digital color data: Outcome of the 2014 ieee grss
data fusion contest,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 6, pp. 2984–2996, 2015.
[108] C. A. Laben and B. V. Brower, “Process for enhancing the spatial resolution ofmultispectral imagery using pan-sharpening,” Jan. 4 2000, US Patent 6,011,875.
[109] B. Aiazzi, S. Baronti, and M. Selva, “Improving component substitution pansharpening through multivariate regression of ms + pan data,” IEEE Trans. on Geosci.
and Remote Sens., vol. 45, no. 10, pp. 3230–3239, 2007.
[110] B. Aiazzi, L. Alparone, S. Baronti, A. Garzelli, and M. Selva, “Mtf-tailored multiscale fusion of high-resolution ms and pan imagery,” Photogrammetric Engineering& Remote Sensing, vol. 72, no. 5, pp. 591–596, 2006.
[111] G. Vivone, R. Restaino, M. Dalla Mura, G. Licciardi, and J. Chanussot, “Contrast
and error-based fusion schemes for multispectral image pansharpening,” IEEE Geoscience and Remote Sensing Letters, vol. 11, no. 5, pp. 930–934, 2014.
[112] J. G. Liu, “Smoothing filter-based intensity modulation: A spectral preserve image
fusion technique for improving spatial details,” International Journal of Remote
Sensing, vol. 21, no. 18, pp. 3461–3472, 2000.
[113] R. C. Hardie, M. T. Eismann, and G. L. Wilson, “Map estimation for hyperspectral image resolution enhancement using an auxiliary sensor,” IEEE Trans. Image
Process., vol. 13, no. 9, pp. 1174 – 1184, 2004.
[114] Q. Wei, J. Bioucas-Dias, N. Dobigeon, and J. Y. Tourneret, “Hyperspectral andmultispectral image fusion based on a sparse representation,” IEEE Trans. on
Geosci. and Remote Sens., vol. 53, no. 7, pp. 3658 – 3668, 2015.
[115] M. Simões, J. Bioucas-Dias, L. B. Almeida, and J. Chanussot, “A convex formulation for hyperspectral image superresolution via subspace-based regularization,”
IEEE Trans. on Geosci. and Remote Sens., vol. 53, no. 6, pp. 3373–3388, 2015.

75

[116] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization
unmixing for hyperspectral and multispectral data fusion,” IEEE Trans. on Geosci.
and Remote Sens., vol. 50, no. 2, pp. 528–537, 2012.
[117] H. N. Gross and J. R. Schott, “Application of spectral mixture analysis and image
fusion techniques for image sharpening,” Remote Sensing of Environment, vol. 63,
no. 2, pp. 85–94, 1998.
[118] C. Lanaras, E. Baltsavias, and K. Schindler, “Hyperspectral super-resolution bycoupled spectral unmixing,” in Proceedings of the IEEE International Conference
on Computer Vision, 2015, pp. 3586–3594.
[119] Y. Zhang, Y. Gao, Y. Liu, and M. He, “Hyperspectral and multispectral image
fusion based on constrained cnmf unmixing,” in 2015 7th Workshop on Hyperspectral
Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2015, pp.
1–4.
[120] M. T. Eismann, Resolution enhancement of hyperspectral imagery using maximum
a posteriori estimation with a stochastic mixing model, Ph.D. dissertation, Dept.
Electrical and Computer Engineering, Univ. Dayton, OH, 2004.
[121] Q. Wei, N. Dobigeon, and J. Y. Tourneret, “Fast fusion of multi-band images basedon solving a sylvester equation,” IEEE Transactions on Image Processing, vol. 24,
no. 11, pp. 4109–4121, 2015.
[122] C. Chen, Y. Li, W. Liu, and J. Huang, “Image fusion with local spectral consistency and dynamic gradient sparsity,” in Proc. IEEE Conf. Comput. Vis. Pattern
Recognit. (CVPR), 2014, pp. 2760–2765.
[123] C.I. Chang and Q. Du, “Interference and noise-adjusted principal components
analysis,” IEEE Trans. on Geosci. and Remote Sens., vol. 37, no. 5, pp. 2387–2396,
1999.
[124] L. Wald, T. Ranchin, and M. Mangolini, “Fusion of satellite images of diﬀerent
spatial resolutions: Assessing the quality of resulting images,” PhotogrammetricEngineering & Remote Sensing, vol. 63, no. 6, pp. 691–699, 1997.
[125] N. Yokoya and A. Iwasaki, “Airborne hyperspectral data over chikusei,” Tech.
Rep., Space Application Laboratory, University of Tokyo, Japan, May 2016.
[126] G. Scarpa, S. Vitale, and D. Cozzolino, “Target-adaptive cnn-based pansharpening,” IEEE Trans. on Geosci. and Remote Sens., vol. 56, no. 9, pp. 5443–5457,
2018.
[127] C. Kwan, J. H. Choi, S. H. Chan, J. Zhou, and B. Budavari, “A super-resolution
and fusion approach to enhancing hyperspectral images,” Remote Sensing, vol. 10,
no. 9, pp. 1416, 2018.

76

[128] P. Liu and L. Xiao, “Multicomponent driven consistency priors for simultaneous
decomposition and pansharpening,” IEEE Sel. Top. Appl. Earth Obs. Remote Sens.,
vol. 12, no. 11, pp. 4589–4605, 2019.
[129] X. Luo, J. Yin, X. Luo, and X. Jia, “A novel adversarial based hyperspectral andmultispectral image fusion,” Remote Sensing, vol. 11, no. 5, pp. 492, 2019.
[130] W. Xie, J. Lei, Y. Cui, Y. Li, and Q. Du, “Hyperspectral pansharpening with deeppriors,” IEEE Trans. Neural Netw. Learn. Syst., 2019.
[131] S. Li and B. Yang, “A new pan-sharpening method using a compressed sensing
technique,” IEEE Trans. on Geosci. and Remote Sens., vol. 49, no. 2, pp. 738–746,
2010.
[132] C. Jiang, H. Zhang, H. Shen, and L. Zhang, “A practical compressed sensing-based
pan-sharpening method,” IEEE Geosci. Remote Sensing Lett., vol. 9, no. 4, pp.
629–633, 2011.
[133] W. Wang, L. Jiao, S. Yang, and K. Rong, “Distributed compressed sensing-based
pan-sharpening with hybrid dictionary,” Neurocomputing, vol. 155, pp. 320–333,
2015.
[134] E. Y. Sidky, J. H. Jørgensen, and X. Pan, “Convex optimization problem prototyping for image reconstruction in computed tomography with the chambolle–pockalgorithm,” Physics in medicine and biology, vol. 57, no. 10, pp. 3065, 2012.
[135] C. Sutour, C.-A. Deledalle, and J.-F. Aujol, “Adaptive regularization of the nlmeans: Application to image and video denoising,” IEEE Trans. Image Process.,
vol. 23, no. 8, pp. 3506–3521, 2014.

77

Publications related to thisdissertation
Article in Journal Papers[J1] S. Takeyama, S. Ono, and I. Kumazawa, “Image restoration with multiple hardconstraints on data-fidelity to blurred/noisy image pair,” IEICE Transactions on
Information and Systems, vol. E100-D, no.9, pp. 1953-1961, Sep. 2017.
[J2] S. Takeyama, S. Ono, and I. Kumazawa, “A constrained convex optimization approach to hyperspectral image restoration with hybrid spatio-spectral regularization,” Remote Sensing, vol. 12, pp. 3541, 2020.

Peer-Reviewed Articles in International Conference[C1] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral image restoration by hybrid spatio-spectral total variation,” IEEE International Conference on Acoustics,
Speech, and Signal Processing (ICASSP), New Orleans, LA, pp. 4586-4590, Mar.
2017.
[C2] S. Takeyama, S. Ono, and I. Kumazawa, “Robust and eﬀective hyperspectral pansharpening using spatio-spectral total variation,” IEEE International Conference on
Acoustics, Speech, and Signal Processing (ICASSP), Calgary, Canada, pp. 16031607, Apr. 2018.
[C3] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral image restoration and pansharpening,” APSIPA Annual Summit and Conference (APSIPA ASC), Honolulu,
Hawaii, 1 pages, Nov. 2018.
[C4] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral pansharpening using noisy
panchromatic image,” APSIPA Annual Summit and Conference (APSIPA ASC),
Honolulu, Hawaii, 6 pages, Nov. 2018.
[C5] S. Takeyama, S. Ono, and I. Kumazawa, “Hyperspectral and multispectral data
fusion by a regularization considering,” IEEE International Conference on Acoustics,
Speech, and Signal Processing (ICASSP), Brighton, UK, pp. 2152-2156, May 2019.
[C6] S. Takeyama, S. Ono, and I. Kumazawa, “Mixed noise removal for hyperspectral
images using hybrid spatio-spectral total variation,” IEEE International Conference78

on Image Processing (ICIP), virtual conference, pp. 3128-3132, Sep. 2019.
[C7] S. Takeyama and S. Ono, “Joint mixed-noise removal and compressed sensing reconstruction of hyperspectral image via convex optimization,” IEEE International
Geoscience and Remote Sensing Symposium (IGARSS), virtual symposium, Sep.
2020.
[C8] S. Takeyama and S. Ono, “Compressed hyperspectral pansharpening,” IEEE International Conference on Image Processing (ICIP), virtual conference, pp. 2855-2859,
Oct. 2020.

79