1Hyperspectral pansharpening: a review

arXiv:1504.04531v1 [cs.CV] 17 Apr 2015

Laetitia Loncan, Luis B. Almeida, José M. Bioucas-Dias, Xavier Briottet, Jocelyn Chanussot, Nicolas Dobigeon,
Sophie Fabre, Wenzhi Liao, Giorgio A. Licciardi, Miguel Simões, Jean-Yves Tourneret, Miguel A. Veganzones,
Gemine Vivone, Qi Wei and Naoto Yokoya

Abstract—Pansharpening aims at fusing a panchromatic image
with a multispectral one, to generate an image with the highspatial resolution of the former and the high spectral resolutionof the latter. In the last decade, many algorithms have beenpresented in the literature for pansharpening using multispectral
data. With the increasing availability of hyperspectral systems,
these methods are now being adapted to hyperspectral images. In
this work, we compare new pansharpening techniques designed
for hyperspectral data with some of the state of the art methods
for multispectral pansharpening, which have been adapted for
hyperspectral data. Eleven methods from different classes (component substitution, multiresolution analysis, hybrid, Bayesianand matrix factorization) are analyzed. These methods are applied to three datasets and their effectiveness and robustness areevaluated with widely used performance indicators. In addition,
all the pansharpening techniques considered in this paper havebeen implemented in a MATLAB toolbox that is made available
to the community.
Index Terms—Pansharpening, hyperspectral image, image fusion

I. I NTRODUCTIONI

N the design of optical remote sensing systems, owing to
the limited amount of incident energy, there are criticaltradeoffs between the spatial resolution, the spectral resolution,and signal-to-noise ratio (SNR). For this reason, optical systems can provide data with a high spatial resolution but witha small number of spectral bands (for example, panchromatic
data with decimetric spatial resolution or multispectral datawith three to four bands and metric spatial resolution, like
PLEIADES [1]) or with a high spectral resolution but with
reduced spatial resolution (for example, hyperspectral data,
L. Loncan is with Gipsa-lab, Grenoble, France and ONERA(email:laetitia.loncan@onera.fr).
Luis B. Almeida and J. M. Bioucas-Dias are with Instituto de
Telecomunicações, Instituto Superior Técnico, Universidade de Lisboa,(email:{luis.almeida,bioucas@lx.it.pt}).
J. Chanussot, G. Licciardi, and M. Veganzones are with Gipsa-lab, Grenoble, France (email:{Jocelyn.chanussot,Giorgio-Antonino.Licciardi, miguelangel.veganzones}@gipsa-lab.grenoble-inp.fr).
X. Briottet and S. Fabre are with ONERA, Toulouse, France (email:{Xavier.briottet, Sophie.fabre}@onera.fr).
N. Dobigeon, J.-Y. Tourneret, and Q. Wei are with Universityof Toulouse, IRIT/INP-ENSEEIHT (email:{nicolas.dobigeon, jeanyves.tourneret,qi.wei}@enseeiht.fr).
W.
Liao
is
with
Ghent
University,
Ghent,
Belgium(email:wliao@telin.ugent.be).
M. Simões is with Instituto de Telecomunicações, Instituto Superior Técnico, Universidade de Lisboa and Gipsa-lab, Grenoble, France(email:miguel.simoes@lx.it.pt).
G. Vivone is with North Atlantic Treaty Organization (NATO) Scienceand Technology Organization (STO) Centre for Maritime Research andExperimentation (CMRE) (email:gvivone@unisa.it).
N. Yokoya is with University of Tokyo (email:yokoya@sal.rcast.utokyo.ac.jp).

subsequently referred to as HS data, with more than one hundred of bands and decametric spatial resolution like HYPXIM[2]). To enhance the spatial resolution of multispectral data,
several methods have been proposed in the literature under thename of pansharpening, which is a form of superresolution.
Fundamentally, these methods solve an inverse problem which
consists of obtaining an enhanced image with both highspatial and high spectral resolutions from a panchromatic
image and a multispectral image. The huge interest of thecommunity on this topic is evidenced by the existence ofsessions dedicated to this topic in the most important remotesensing and Earth observation conferences as well as by thelaunch of public contests, of which the one sponsored by thedata fusion committee of the IEEE Geoscience and Remote
Sensing society [3] is an example.
A taxonomy of pansharpening methods can be found in theliterature [4], [5], [6]. They can be broadly divided into fourclasses: component substitution (CS), multiresolution analysis(MRA), Bayesian, and variational. The CS approach relies on
the substitution of a component (obtained, e.g., by a spectral
transformation of the data) of the multispectral (subsequently
denoted as MS) image by the panchromatic (subsequently
denoted as PAN) image. The CS class contains algorithms
such as intensity-hue-saturation (IHS) [7], [8], [9], principal
component analysis (PCA) [10], [11], [12] and Gram-Schmidt(GS) spectral sharpening [13]. The MRA approach is based on
the injection of spatial details, which are obtained through amultiscale decomposition of the PAN image into the MS data.
The spatial details can be extracted according to several modalities of MRA: decimated wavelet transform (DWT) [14], undecimated wavelet transform (UDWT) [15], ”à-trous” wavelet
transform (ATWT) [16], Laplacian pyramid [17], nonseparable
transforms, either based on wavelets (e.g., contourlets [18]) ornot (e.g., curvelets [19]). Hybrid methods have been also proposed, which use both component substitution and multiscaledecomposition, such as guided filter PCA (GFPCA), describedin Section II-C. The Bayesian approach relies on the used ofposterior distribution of the full resolution target image given
the observed MS and PAN images. This posterior, which is theBayesian inference engine, has two factors: a) the likelihood
function, which is the probability density of the observed MS
and PAN images given the target image, and b) the priorprobability density of the target image, which promotes target
images with desired properties, such as being segmentallysmooth. The selection of a suitable prior allows us to copewith the usual ill-posedness of the pansharpening inverse
problems. The variational class is interpretable as particularcase of the Bayesian one, where the target image is estimatedby maximizing the posterior probability density of the full

2

resolution image. The works [20], [21], [22] are representativeof the Bayesian and variational classes. As indicated in Table
I, the CS, MRA, and Hybrid classes of methods are detailedin Sections II-A, II-B, and II-C, respectively. Herein, theBayesian class is not addressed in the MS+PAN context. It
is addressed in detail, however, in Section II-D in the context
of HS+PAN fusion.
With the increasing availability of HS systems, the pansharpening methods are now extending to the fusion of HS
and panchromatic images [23], [24], [25], [26]. Pansharpening of HS images is still an open issue, and very few
methods are presented in the literature to address it. The
main advantage of HS image with respect to MS one is themore accurate spectral information they provide, which clearly
benefits many applications such as unmixing [27], change
detection [28], object recognition [29], scene interpretation[30] and classification [31]. Several of the methods designedfor HS pansharpening were originally designed for the fusionof MS and HS data[32]–[36], the MS data constituting the highspatial resolution image. In this case, HS pansharpening can beseen as a particular case, where the MS image is composed ofa single band, and thus reduces to a PAN image. In this paper,
we divide these methods into two classes: Bayesian methods
and matrix factorization based methods. In Section II-D, we
briefly present the algorithms of [33], [36], and [35] of theformer class and in Section II-E the algorithm of [32] of thelatter class.
As one may expect, performing pansharpening with HSdata is more complex than performing it with MS data.
Whereas PAN and MS data are usually acquired almost in
the same spectral range, the spectral range of an HS image
normally is much wider than the one of the corresponding
PAN image. Usually, the PAN spectral range is close to
the visible spectral range of 0.4 − 0.8µm (for example, theadvanced land imager–ALI–instrument acquires PAN data in
the range 0.48 − 0.69µm). The HS range often covers thevisible to the shortwave infrared (SWIR) range (for example,
Hyperion acquires HS data in the range 0.4−2.5µm, the range
0.8−2.5µm being not covered by the PAN data). The difficultythat arises, consists in defining a fusion model that yields goodresults in the part of the HS spectral range that is not covered
by PAN data,in which the high resolution spatial information
is missing. This difficulty already existed, to some extent, in
MS+PAN pansharpening, but it is much more severe in the
HS+PAN case.
To the best of the authors’ knowledge, there is currently
no study comparing different fusion methods for HS data,
particularly on datasets where the spectral domain of the HS
image is larger than the one of the PAN image. This work
aims at addressing this specific issue. The remainder of thepaper is organized as follows. Section II reviews the methods
under study, i.e., CS, MRA, hybrid, Bayesian, and matrix
decomposition approaches. Section III summarizes the quality
assessment measures that will be used to assess the image
fusion results. Experimental results are presented in Section
IV. Conclusions are drawn in Section V.

TABLE IS UMMARY OF THE DIFFERENT CLASSES OF METHODS CONSIDERED IN
THIS PAPER . W ITHIN PARENTHESES , WE INDICATE THE ACRONYM OF
EACH METHOD , FOLLOWED BY THE NUMBER OF THE SECTION IN WHICH
THAT METHOD IS DESCRIBED .

M ETHODS ORIGINALLY DESIGNED FOR MS PANSHARPENING
Component substitution (CS,II-A)
Principal Component Analysis(PCA, II-A1)
Gram Schmidt (GS, II-A2)

Multiresolution analysis (MRA,
II-B)
Smoothing filter-based intensity
modulation (SFIM, II-B1)
Laplacian pyramid (II-B2)

Hybrid methods (II-C)
Guided Filter PCA (GFPCA)

Bayesian methods
Not discussed in this paper

M ETHODS ORIGINALLY DESIGNED FOR HS PANSHARPENING
Bayesian Methods (II-D)
Naive Gaussian prior (II-D1)
Sparsity promoting prior (II-D2)
HySure (II-D3)

Matrix Factorization (II-E)
Coupled Non-negative Matrix Factorization (CNMF)

II. H YPERSPECTRAL PANSHARPENING T ECHNIQUESThis section presents some of the most relevant methodsfor HS pansharpening. First, we focus on the adaptation of
the popular CS and MRA MS pansharpening methods for HS
pansharpening. Later, we consider more recent methods basedon Bayesian and matrix factorization approaches. A toolbox
containing MATLAB implementations of these algorithms canbe found online1 .
Before presenting the different methods, we introduce notation used along the paper. Bold-face capital letters refer tomatrices and bold-face lower-case letters refer to vectors. The
notation Xk refers to the kth row of X. The operator ()Tdenotes the transposition operation. Images are represented bymatrices, in which each row corresponds to a spectral band,
containing all the pixels of that band arranged in lexicographic
order. We use the following specific matrices:m ×n• X = [x1 , . . . , xn ] ∈ R λrepresents the full resolution
b represents
target image with mλ bands and n pixels; X
an estimate of that image.
m ×m• YH ∈ R λ, YM ∈ Rnλ ×n , and P ∈ R1×n
represents, respectively, the observed HS, MS, and PAN
images, nλ denoting the number of bands of the MS
image and m the total number of pixel in the YH image.
e H ∈ Rmλ ×n represents the HS image YH interpolated• Y
at the scale of the PAN image.
p
We denote by d =m/n the down-sampling factor,
assumed to be the same in both spatial dimensions.
A. Component Substitution
CS approaches rely upon the projection of the higher
spectral resolution image into another space, in order toseparate spatial and spectral information [6]. Subsequently, thetransformed data are sharpened by substituting the component
that contains the spatial information with the PAN image (orpart of it). The greater the correlation between the PAN image
1 http://openremotesensing.net

3and the replaced component, the less spectral distortion will
be introduced by the fusion approach [6]. As a consequence,a histogram-matching procedure is often performed before
replacing the PAN image. Finally, the CS-based fusion process
is completed by applying the inverse spectral transformation
to obtain the fused image.
The main advantages of the CS-based fusion techniques are
the following: i) high fidelity in rendering the spatial details
in the final image [37], ii) fast and easy implementation [8],and iii) robustness to misregistration errors and aliasing [38].
On the negative side, the main shortcoming of this class oftechniques is the generation of a significant spectral distortion,
cause by the spectral mismatch between the PAN and the HS
spectral ranges [6].
Following [4], [39], a formulation of the CS fusion scheme
is given bybkX

e k + gk (P − OL ) ,
=Y
H

(1)

b k denotes the kth band of thefor k = 1, . . . , mλ , where X
estimated full resolution target image, g = [g1 , . . ., gmλ ]T isa vector containing the injection gains, and OL is defined as
OL =

mλ
X

ei ,
wi Y
H

(2)

i=1

where the weights w = [w1 , . . . , wi , . . . , wmλ ]T measure
the spectral overlap among the spectral bands and the PAN
image [6], [40].
The CS family includes many popular pansharpening approaches. In [26], three approaches based on principal component analysis (PCA) [9] and Gram-Schmidt [13], [37] transformations have been compared for sharpening HS data. A brief
description of these techniques follows.
1) Principal Component Analysis: PCA is a spectral
transformation widely employed for pansharpening applications [9]. It is achieved through a rotation of the original data
(i.e., a linear transformation) that yields the so-called principal
components (PCs). The hypothesis underlying its application
to pansharpening is that the spatial information (shared by all
the channels) is concentrated in the first PC, while the spectral
information (specific to each single band) is accounted forthe other PCs. The whole fusion process can be described by
the general formulation stated by Eqs. (1) and (2), where thevectors w and g of coefficient vectors are derived by the PCA
procedure applied to the HS image.
2) Gram-Schmidt: The Gram-Schmidt transformation, often exploited in pansharpening approaches, was initially proposed in a patent by Kodak [13]. The fusion process starts byusing, as the component, a synthetic low resolution PAN image
IL at the same spatial resolution as the HS image2 . A complete
orthogonal decomposition is then performed, starting with thatcomponent. The pansharpening procedure is completed bysubstituting that component with the PAN image, and inverting
the decomposition. This process is expressed by (1) using the
2 GS is a more general method than PCA. PCA can be obtained, in GS, byusing the first PC as the low resolution panchromatic image [41].

gains [37]gk =e k , OL )cov(YHvar(OL )(3)

for k = 1, . . . , mλ , where cov (·, ·) and var (·) denote thecovariance and variance operations. Different algorithms are
obtained by changing the definition of the weights in (2).
The simplest way to obtain this low-resolution PAN image
simply consists of averaging the HS bands (i.e., by setting
wi = 1/mλ , for i = 1, . . . , mλ ). In [37], the authors proposedan enhanced version, called GS Adaptive (GSA), in which IL
is generated by the linear model in (2) with weights estimatedby the minimization of the mean square error between theestimated component and a filtered and downsampled versionof the PAN image.
B. Multiresolution Analysis
Pansharpening methods based on MRA apply a spatial
filter to the PAN image for generating details to be injected
into the HS data. The main advantages of the MRA-based
fusion techniques are the following: i) temporal coherence [5](see Sect.27.4.4), ii) spectral consistency, and iii) robustness
to aliasing, under proper conditions [38]. On the negativeside, the main shortcomings are i) the implementation ismore complicated due to the design of spatial filters, ii) thecomputational burden is usually larger when compared to CS
approaches. The fusion step is summarized as [4], [39]bk = Ye k + Gk ⊗ (P − PL ) ,
X
H

(4)

for k = 1, . . . , mλ , where PL denotes a low-pass versionof P, and the symbol ⊗ denotes element-wise multiplication.
Furthermore, an equalization between the PAN image and the
HS spectral bands is often required. P − PL is often called
the details image, because it is a high-pass version of P, andEq. (4) can be seen as describing the way to inject details
into each of the bands of the HS image. According to (4), theapproaches belonging to this category can differ in i) the typeof PAN low pass image PL that is used, and ii) the definitionof the gain coefficients Gk . Two common options for defining
the gains are:1) Gk = 1 for k = 1, . . . , mλ , where 1 is an appropriately
sized matrix with all elements equal to 1. This choiceidentifies the so-called additive injection scheme;e k PL for k = 1, . . . , mλ , where the symbol2) Gk = YHdenotes element-wise division. In this case, the details
are weighted by the ratio between the upsampled HS
image and the low-pass filtered PAN one, in orderto reproduce the local intensity contrast of the PAN
image in the fused image [42]. This coefficient selection
is often referred to as high pass modulation (HPM)
method or multiplicative injection scheme. Some possible numerical issues could appear due to the divisione k and PL for low value of PL creating
between YH
fused pixel with very high value. In our toolbox thisproblem is addressed by clipping these values by using
the information given by the dynamic range.
In the case of HS pansharpening, some further considerations should be taken into account. Indeed, the PAN and4

HS images are rarely acquired with the same platform. Thus,
the ratio between the spatial resolutions of the PAN and HS
images may not always be an integer number, or a power of
two. This implies that some of the conventional approaches
initially developed for MS images cannot be extended in a
simple way to HS images (for example, dyadic wavelet-basedalgorithms cannot be applied in these conditions).
1) Smoothing Filter-based Intensity Modulation (SFIM):
The direct implementation of Eq. (4) consists of applying asingle linear time-invariant (LTI) low pass filter (LPF) hLP to
the PAN image P for obtaining PL . Therefore, we havebk = Ye k + gk (P − P ∗ hLP )
X
H

(5)

for k = 1, . . . , mλ , where the symbol ∗ denotes the convolution operator. The SFIM algorithm [43] sets hLP to a simple
box (i.e., an averaging) filter and exploits HPM as the details
injection scheme.
2) Laplacian Pyramid: The low-pass filtering needed toobtain the signal PL at the original HS scale can be performedin more than one step. This is commonly referred to as
pyramidal decomposition and dates back to the seminal workof Burt and Adelson [17]. If a Gaussian filter is used to lowpass filter the images in each step, one obtains a so-called
Gaussian pyramid. The differences between consecutive levelsof a Gaussian pyramid define the so-called Laplacian pyramid.
The suitability of the latter to the pansharpening problem
has been shown in [44]. Indeed, Gaussian filters can betuned to closely match the sensor modulation transfer function(MTF). In this case, the unique parameter that characterizes
the whole distribution is the Gaussian’s standard deviation,
which is determined from sensor-based information (usuallyfrom the value of the amplitude response at the Nyquist
frequency, provided by the manufacturer). Both additive andmultiplicative details injection schemes have been used in thisframework [42], [45]. They will be referred to as MTF Generalized Laplacian Pyramid (MTF-GLP) [45] and MTFGLP with High Pass Modulation (MTF-GLP-HPM) [42],respectively.
C. Hybrid Methods
Hybrid approaches use concepts from different classes ofmethods, namely from CS and MRA ones, as explained next.
1) Guided Filter PCA (GFPCA): One of the main challenges for fusing low-resolution HS and high-resolution
PAN/RGB data is to find an appropriate balance between spectral and spatial preservation. Recently, the guided filter [46] hasbeen used in many applications (e.g. edge-aware smoothing
and detail enhancement), because of its efficiency and strong
ability to transfer the structures of the guidance image to thefiltering output. Its application to HS data can be found in [47],where the guided filter was applied to transfer the structuresof the principal components of the HS image to the initialclassification maps.
Here, we briefly describe an image fusion framework whichuses a guided filter in the PCA domain (GFPCA) [48].
The approach won the “Best Paper Challenge” award at
the 2014 IEEE data fusion contest [48], by fusing a low

spatial resolution thermal infrared HS image and a high spatial
resolution, visible RGB image associated with the same scene.
Fig. 1 shows the framework of GFPCA. Instead of using
CS, which may cause spectral distortions, GFPCA uses ahigh resolution PAN/RGB image to guide the filtering process
aimed at obtaining superresolution. In this way, GFPCA does
not only preserve the spectral information from the original
HS image, but also transfers the spatial structures of the highresolution PAN/RGB image to the enhanced HS image. To
speed up the processing, GFPCA first uses PCA to decorrelate
the bands of the HS image, and to separate the information
content from the noise. The first p  mλ PCA channels
contain most of the energy (and most of the information) ofan HS image, and the remaining mλ −p PCA channels mainlycontain noise (recall that mλ is the number of spectral bands ofthe HS image). When applied to these noisy (and numerous)mλ − p channels, the guided filter amplifies the noise andcauses a high computational cost in processing the data, whichis undesirable. Therefore, guided filtering is used to enlarge
only the first k PCA channels, preserving the structures of
the PAN/RGB image, while cubic interpolation is used to
upsample the remaining channels.
Let PCi , with (i ≤ p), denote the ith PC channel obtained
from the HS image YH , with its resolution increased to thatof the guided image Y (Y may be a PAN or an RGB image)through bicubic interpolation. The output of the filtering, PC0i ,can be represented as an affine transformation of Y in a local
window ωj of size (2d + 1) × (2d + 1) as follows:
PC0i = aj Y + bj , ∀i ∈ ωj .

(6)

The above model ensures that the output PC0i has an edge
only if the guided image Y has an edge, since ∇(PC0i ) =a∇Y. The following cost function is used to determine thecoefficients aj and bj :X(aj Y + bj − PCi )2 + a2j ,(7)
E(aj , bj ) =i∈ωjwhere  is a regularization parameter determining the degreeof blurring for the guided filter. For more details about theguided filtering scheme, we invite the reader to consult [46].
The cost function E leads the term aj Y + bj to be as close
as possible to PCi , in order to ensure the preservation of theoriginal spectral information. Before applying inverse PCA,GFPCA also removes the noise from the remaining PCA
channels P CRest using a soft-thresholding scheme (similarly
to [49]), and increases their spatial resolution to the resolution

Fig. 1.

Fusion of HS and PAN/RGB images with the GFPCA framework.

5of the PAN/RGB image using cubic interpolation only (without
guided filtering).

The fusion of HS and high spatial resolution images, e.g.,
MS or PAN images, can be conveniently formulated within
the Bayesian inference framework. This formulation allows an
intuitive interpretation of the fusion process via the posterior
distribution of the Bayesian fusion model. Since the fusion
problem is usually ill-posed, the Bayesian methodology offersa convenient way to regularize the problem by defining anappropriate prior distribution for the scene of interest. Following this strategy, different Bayesian estimators for fusing
co-registered high spatial-resolution MS and high spectralresolution HS images have been designed [33]–[36], [50]–[54]. The observation models associated with the HS and MS
images can be written as follows [50], [55], [56](8)where X, YH , and YM were defined in Section II, andn×n• B ∈ Ris a cyclic convolution operator, corresponding to the spectral response of the HS sensor expressedin the resolution of the MS or PAN image,n×m• S ∈ Ris a down-sampling matrix with downsampling factor d,n ×mλ• R∈R λis the spectral response of the MS or PAN
sensor,• NH and NM are the HS and MS noises, assumed tohave zero mean Gaussian distributions with covariancematrices ΛH and ΛM , respectively.
For the sake of generality, the formulation in this section assumes that the observed data is the pair of matrices (YH , YM ).
Since a PAN image can be represented by YM with nλ = 1,
the observation model (8) covers the HS+PAN fusion problem
considered in this paper.
Using geometrical considerations well grounded in the HS
imaging literature devoted to the linear unmixing problem[27], the high spatial resolution HS image to be estimated isassumed to live in a low dimensional subspace. This hypothesis is very reliable when the observed scene is composed of afinite number of macroscopic materials (called endmembers).
Based on the model (8) and on the low dimensional subspaceassumptions, the distributions of YH and YM can be expressed
as follows
YH |U ∼ MN mλ ,m (HUBS, ΛH , Im ),YM |U ∼ MN nλ ,n (RHU, ΛM , In )p (U|YH , YM ) ∝ p (YH |U) p (YM |U) p (U)or equivalently

D. Bayesian Approaches

YH = XBS + NH
YM = RX + NM

of U can be written as(9)where MN represents the matrix normal distribution [57], theeλcontaining intarget image is X = HU, with H ∈ Rmλ ×mits columns a basis of the signal subspace of size me λ  mλe λ ×nand U ∈ Rm
contains the representation coefficients of X
with respect to H. The subspace transformation matrix H canbe obtained via different approaches, e.g., PCA [58] or vertex
component analysis [59].
According to Bayes’ theorem and using the fact that thenoises NH and NM are independent, the posterior distribution

(10)

3

2
. 1 − 21
ΛH (YH − HUBS) F +− log p (U|YH , YM ) ={z}|2
HS data term.=logp(YH |U)1 − 21
2
ΛM (YM − RHU) F + λφ(U)| {z }|2{z} regularizer.
. MS data term=log p(U)=log p(YM |U)

(11)
def p
where kXkF =Tr(XXT ) is the Frobenius norm of X.
An important quantity in the negative log-posterior (11) is thepenalization term φ(U) which allows the inverse problem (8)to be regularized. The next sections discuss different ways ofdefining this penalization term.
1) Naive Gaussian prior: Denote as ui (i = 1, · · · , n) thecolumns of the matrix U that are assumed to be mutuallyindependent and are assigned the following Gaussian priordistributionsp (ui |µi , Σi ) = N (µi , Σi )(12)

where µi is a fixed image defined by the interpolated HS
image projected into the subspace of interest, and Σi is an unknown hyperparameter matrix. Different interpolations can beinvestigated to build the mean vector µi . In this paper, we havefollowed the strategy proposed in [50]. To reduce the numberof parameters to be estimated, the matrices Σi are assumed tobe identical, i.e., Σ1 = · · · = Σn = Σ. The hyperparameterΣ is assigned a proper prior and is estimated jointly with theother parameters of interest. To infer the parameter of interest,
namely the projected highly resolved HS image U, fromthe posterior distribution p (U|YH , YM ), several algorithms
have been proposed. In [33], [34], a Markov chain Monte
Carlo (MCMC) method is exploited to generate a collectionof NMC samples that are asymptotically distributed according
to the target posterior. The corresponding Bayesian estimators
can then be approximated using these generated samples. For
instance, the minimum mean square error (MMSE) estimatorof U can be approximated by an empiricalof the
PNMC average
1(t)b MMSE ≈generated samples U
U
,
where
t=Nbi +1
NMC −Nbi
Nbi is the number of burn-in iterations required to reach thesampler convergence, and U(t) is the image generated in
the tth iteration. The highly-resolved HS image can finallyb MMSE = HUb MMSE . An extension of thebe computed as X
proposed algorithm has been proposed in [53] to handle thespecific scenario of an unknown sensor spectral response. In[60], a deterministic counterpart of this MCMC algorithm has
been developed, where the Gibbs sampling strategy of [33]has been replaced with a block coordinate descent method tocompute the maximum a posteriori (MAP) estimator. Finally,
very recently, a Sylvester equation-based explicit solution of.
3 We use the symbol =to denote equality apart from an additive constant.
The additive constants are irrelevant, since the functions under consideration
are to be optimized, and the additive constants don’t change the locations of
the optima.

6

the related optimization problem has been derived in [61],leading to a significant decrease of the computational complexity.
2) Sparsity promoted Gaussian prior: Instead of incorporating a simple Gaussian prior or smooth regularization for
the HS and MS fusion [34], [50], [51], a sparse representation can be used to regularize the fusion problem. Morespecifically, image patches of the target image (projected
onto the subspace defined by H) are represented as a sparselinear combination of elements from an appropriately chosen
over-complete dictionary with columns referred to as atoms.
Learning the dictionary from the observed images instead ofusing predefined bases [62]–[64] generally improves image
representation [65], which is preferred in most scenarios.
Therefore, an adaptive sparse image-dependent regularization
can be explored to solve the fusion problem (8). In [36], thefollowing regularization term was introduced:. 1
φ(U) ∝ − log p (U) =
2

m
eλ
X

sensors B and R from the data themselves. Therefore, thecomplete methodology can be classified as a blind superresolution method, which, in contrast to the classical blind linear
inverse problems, is tackled by solving two convex problems.
The VTV regularizer (see [68]) is given by
v
m
eλ n
n u
uX
X
2 o


2 
t
, (14)φ U =(UDh )kj + (UDv )kj
j=1

of matrix A, and the products by matrices Dh and Dv
compute the horizontal and vertical discrete differences of animage, respectively, with periodic boundary conditions.
The HS pansharpened image is the solution to the following
optimization problem
minimize

Uk − P D̄k Āk2F, (13)

k=1

where Akj denotes the element in the kth row and jth column

U

2
2
λm
1
YH − HUBS +YM − RHU
22FF+ λφ φ U),
(15)

k=1

where
n
m
e ×n
• Uk ∈ R is the kth band (or row) of U ∈ R λ, withk = 1, . . . , me λ,n ×npat7→ Rn×1 is a linear operator that averages• P(·) : R p
the overlapping patches4 of each band, npat being thenumber of patches associated with the ith band,n ×nat• D̄i ∈ R p
is the overcomplete dictionary, whose
columns are basis elements of size np (corresponding to
the size of a patch), nat being the number of dictionary
atoms, andn ×npat• Āi ∈ R atis the code of the ith band.
Inspired by hierarchical models frequently encountered in
Bayesian inference [67], a second level of hierarchy canbe considered in the Bayesian paradigm by including thecode A within the estimation, while fixing the support Ω̄ ,Ω̄1 , . . . , Ω̄me λ of the code A. Once D̄, Ω̄ and H have beenlearned from the HS and MS data, maximizing the posterior
distribution of U and A reduces to a standard constrained
quadratic optimization problem with respect to (w.r.t.) U and
A. The resulting optimization problem is difficult to solve
due to its large dimension and due to the fact that the linear
operators H(·)BD and P(·) cannot be easily diagonalized.
To cope with this difficulty, an optimization technique thatalternates minimization w.r.t. U and A has been introducedin [36] (where details on the learning of D̄, Ω̄ and H can befound). In [61], the authors show that the minimization w.r.t.
U can be achieved analytically, which greatly accelerates thefusion process.
3) HySure: The works [35], [54] introduce a convex regularization problem which can be seen under a Bayesianframework. The proposed method uses a form of vector total
variation (VTV) [68] for the regularizer φ(U), taking into
account both the spatial and the spectral characteristics ofthe data. In addition, another convex problem is formulated
to estimate the relative spatial and spectral responses of the
4 A decomposition into overlapping patches was adopted, to prevent theoccurrence of blocking artifacts [66].

where λm and λφ control the relative weights of the different
terms. The optimization problem (15) is hard to solve, essentially for three reasons: the downsampling operator BS is
not diagonalizable, the regularizer φ U) is nonquadratic andnonsmooth, and the target image has a very large size. Thesedifficulties were tackled by solving the problem via the split
augmented lagrangian shrinkage algorithm (SALSA) [69], an
instance of ADMM. As an alternative, the main step of the
ADMM scheme can be conducted using an explicit solutionof the corresponding minimization problem, following thestrategy in [61].
The relative spatial and spectral responses B and R were
estimated by solving the following optimization problem:
minimize
B,R

RYH − YM BS

2+ λB φB (B) + λR φR (R)(16)
where φB (·) and φr (·) are quadratic regularizers and λb , λR ≥0 are their respective regularization parameters.
E. Matrix factorization
The matrix factorization approach for HS+MS fusion essentially exploits two facts: 1) A basis or dictionary H for
the signal subspace can be learned from the HS observed
image YH , yielding the factorization X = HU; 2) using thisdecomposition in the second equation of (9) and for negligible
noise, i.e., NM ' 0, we have YH = RHU. Assuming thatthe columns of RH are full rank or that the columns of U
admit a sparse representation w.r.t. the columns of RH, thenb and use it to
we can recover the true solution, denoted by U,
b = HU.
b The works [32], [70]–compute the target image as X[74] are representative of this line of attack. In what follow,
we detail the application of the coupled nonnegative matrix
factorization (CNMF) [32] to the HS+PAN fusion problem.
The CNMF was proposed for the fusion of low spatial
resolution HS and high spatial resolution MS data to produce
fused data with high spatial and spectral resolutions [32]. It
is applicable to HS pansharpening as a special case, in which

7

the higher spatial resolution image has a single band [75].
CNMF alternately unmixes both sources of data to obtain theendmember spectra and the high spatial resolution abundancemaps.
To describe this method, it is convenient to first briefly
introduce linear mixture models for HS images. These models
are commonly used for spectral unmixing, owing to their
physical effectiveness and mathematical simplicity [27]. The
spectrum at each pixel is assumed to be a linear combinationof several endmember spectra. Therefore, X ∈ Rmλ ×n is
formulated as
X = HU

(17)

mλ ×p

is the signature matrix, containing thewhere H ∈ R
spectral representations of the endmembers, and U ∈ Rp×n
is the abundance matrix, containing the relative abundancesof the different endmembers at the various pixels, with p
representing the number of endmembers. By substituting (17)into (8), YH and YM can be approximated as
YH ≈ HUH
YM ≈ HM U(18)

where UH = UBS and HM = RH. CNMF alternately
unmixes YH and YM in the framework of nonnegative matrix
factorization (NMF) [76] to estimate H and U under theconstraints of the relative sensor characteristics. CNMF starts
with NMF unmixing of the low spatial resolution HS data.
The matrix H can be initialized using, for example, the vertex
component analysis (VCA) [59], and H and UH are thenalternately optimized by minimizing kYH − HUH k2F using
Lee and Seung’s multiplicative update rules [76]. Next, U
is estimated from the higher spatial resolution data. HM isset to RH and U is initialized by the spatially up-sampled
matrix of UH obtained by using bilinear interpolation. For
HS pansharpening (nλ =1), only U is optimized by minimizing
kYM − HM Uk2F with the multiplicative update rule, whereas
both HM and U are alternately optimized in the case of
HS+MS data fusion. Finally, the high spatial resolution HS
data can be obtained by the multiplication of H and U.
The abundance sum-to-one constraint is implemented using
a method given in [77], where the data and signature matrices
are augmented by a row of constants. The relative sensor
characteristics, such as BS and R, can be estimated fromthe observed data sources [78].
III. Q UALITY A SSESSMENT OF FUSION PRODUCTS
Quality assessment of a pansharpened real-life HS image
is not an easy task [79], [9], since a reference image is
generally not available. When such an image is not available,
two kinds of comparisons can be performed: i) Each bandof the fused image can be compared with the PAN image,with an appropriate criterion. The PAN image can also becompared with the PAN image reconstructed from the fused
image. ii) The fused image can be spatially degraded to theresolution of the original HS image. The two images can thenbe compared, to assess to what extent the spectral information
has been modified by the fusion method.

In order to be able to use a reference image for quality
assessment, one normally has to resort to the use of semisynthetic HS and PAN images. In this case, a real-life HS
image is used as reference. The HS and PAN images to beprocessed are obtained by degrading this reference image. A
common methodology for obtaining the degraded images isWald’s protocol, described in the next subsection. In orderto evaluate the quality of the fused image with respect tothe reference image, a number of statistical measures can becomputed. The most widely used ones are described ahead,and used in the experiments reported in Section IV-B.
A. Wald’s Protocol
A general paradigm for quality assessment of fused images
that is usually accepted in the research community was first
proposed by Wald et al. [79]. This paradigm is based on
two properties that the fused data have to have, as much as
possible, namely consistency and synthesis properties. The
first property requires the reversibility of the pansharpening
process: it states that the original HS image should be obtainedby properly degrading the pansharpened image. The second
property requires that the pansharpened image be as similar
as possible to the image of the same scene that would beobtained, by the same sensor, at the higher resolution. Thiscondition entails that both the features of each single bandand the mutual relations among bands have to be preserved.
However, the definition of an assessment method that fulfills
these constraints is still an open issue [80], [81], and closely
relates to the general discussion regarding image quality
assessment [82] and image fusion [83], [84].
Wald’s protocol for assessing the quality of pansharpening methods [79], depicted in Fig. 2, synthetically generates
simulated observed images from a reference HS image, andthen evaluates the pansharpening methods’ results against thatreference image. The protocol consists of the following steps:• Given a HS image, X, to be used as reference, a simulated observed low spatial resolution HS image, YH ,is obtained by applying a Gaussian blur to X, and thendownsampling the result by selecting one out of every
d pixels in both the horizontal and vertical directions,where d denotes the downsampling factor.
• A simulated PAN image, P, is obtained by multiplying
the reference HS image, on the left, by a suitably chosen
spectral response vector, P = rT X.
• The pansharpening method to be evaluated is applied
to the simulated observations YH and P, yielding theestimated superresolution HS image, X̂.
• Finally, the estimated superresolution HS image and thereference one are compared, to obtain quantitative quality
measures.
B. Quality Measures
Several quality measures have been defined in the literature,in order to determine the similarity between estimated andreference spectral images. These measures can be generallyclassified into three categories, depending on whether they8Hyperspectral
super-resolution image
X
Blurring anddownsampling

Pansharpening
spectral response

Observed
hyperspectral image
YH

Observed
pansharpened image
P

Hyperspectral
super-resolution
approach

Estimated hyperspectral
super-resolution image
X̂

Quality measures

Q

where, given the vectors a, b ∈ Rmλ ,ha, biSAM(a, b) = arccos,
kakkbk(21)

ha, bi = aT b is inner product between a and b, and k·k is the`2 norm. SAM is a measure of the spectral shape preservation.
The optimal value of SAM is 0. The values of SAM reported
in our experiments have been obtained by averaging the values
obtained for all the image pixels.
3) Root mean squared error: The RMSE measures the `2berror between the two matrices X and X
b − XkF
Xb X) = k√(22)
RMSE(X,n ∗ mλpwhere kXkF = trace(XT X) is the Frobenius norm of X.
The ideal value of RMSE is 0.
4) Erreur relative globale adimensionnelle de synthèse:
ERGAS offers a global indication of the quality of a fused
image. It is defined as
v
u
2
mλ 
u 1 X
RMSEk
t
b
,
(23)
ERGAS(X, X) = 100 d
mλ
µk
k=1

Fig. 2. Flow diagram of the experimental methodology, derived from Wald’s
protocol (simulated observations), for synthetic and semi-real datasets.

where d is the ratio between the linear resolutions of the PAN
and HS images, defined asd=attempt to measure the spatial, spectral or global quality ofthe estimated image. This review is limited to the most widelyused quality measures, namely the cross correlation (CC),which is a spatial measure, the spectral angle mapper (SAM),
which is a spectral measure, and the root mean squared
error (RMSE) and erreur relative globale adimensionnelle desynthèse (ERGAS) [85], which are global measures. Below
we provide the formal definitions of these measures operatingb ∈ Rmλ ×n and on the reference
on the estimated image X
mλ ×n
bj and xj denote. In the definitions, x
HS image X ∈ R
b and X, respectively, the matrices
the jth columns of X
A, B ∈ R1×n denote two generic single-band images, andAi denotes the ith element of A.
1) Cross correlation: CC, which characterizes the geometric distortion, is defined asmλXb X) = 1
b i , Xi ),
CC(X,
CCS(X
(19)
mλ i=1
where CCS is the cross correlation for a single-band image,defined asPn
j=1 (Aj − µA )(Bj − µB )CCS(A, B) = qP
,
Pn
n
2
2
j=1 (Aj − µA )
j=1 (Bj − µB )Pnwhere, µA = (1/n) j=1 Aj is the sample mean of A. The
ideal value of CC is 1.
2) Spectral angle mapper: SAM, which is a spectral measure, is defined asn
1Xb
SAM(X, X) =SAM(bxj , xj ),
(20)
n j=1

PAN linear spatial resolution,
HS linear spatial resolutionb k − Xk k FkX√, µk is the sample mean of the kthnband of X. The ideal value of ERGAS is 0.
RMSEk =

IV. E XPERIMENTAL R ESULTS
A. Datasets
The datasets that were used in the experimental tests were
all semi-synthetic, generated according to the Wald’s protocol.
In all cases, the spectral bands corresponding to the absorption
band of water vapor, and the bands that were too noisy, were
removed from the reference image before further processing.
Three real-life HS images have been used as reference images
for the Wald protocol. In the following, we describe thedatasets that were generated from these images. Table IIsummarizes their properties. These datasets are expressed inspectral luminance (nearest to the sensor output, without preprocessing) and are correctly registered.
TABLE IIC ARACTERISTIC OF THE THREE DATASETS

datasetMoffett
Camargue
Garonsdimensions
PAN 185 × 395
HS 37 × 79
PAN 500 × 500
HS 100 × 100
PAN 400 × 400
HS 80 × 80

spatial res
20m
100m
4m
20m
4m
20m

N

instrument

224

AVIRIS

125

HyMap

125

HyMap

9

1) Moffett field dataset: This dataset represents a mixed
urban/rural scene. The dimensions of the PAN are 185 × 395
with a spatial resolution of 20m whereas the size of the HS
image is 37 × 79 with a spatial resolution of 100m (whichmeans a spatial resolution ratio of 5 between the two images).
The HS image has been acquired by the airborne hyperspectral instrument airborne visible infrared image spectrometer(AVIRIS). This instrument is characterized by 224 bands
covering the spectral range 0.4 − 2.5µm.
2) Camargue dataset: This dataset represents a rural areawith different kind of crops. The dimensions of the PAN
image are 500 × 500 pixels with a spatial resolution of 4m
whereas the size of the HS image is 100 × 100 pixels with aspatial resolution of 20m, (which means a spatial resolution
ratio of 5 between the two images). The HS image has beenacquired by the airborne hyperspectral instrument HyMap(Hyperspectral Mapper) in 2007. The hyperspectral instrument
is characterized by 125 bands covering the spectral range
0.4 − 2.5µm.
3) Garons dataset: This dataset represents a rural area witha small village. The dimension of the PAN image are 400×400
with a spatial resolution of 4m whereas the size of the HS
image is 80 × 80 with a spatial resolution of 20m, (whichmeans a spatial resolution ratio of 5 between the two images).
This dataset has been acquired with the HyMap instrument in2009.
B. Results and Discussion
Methods presented in Section II have been applied on
the three datasets presented in Section IV-A and analyzed
following the Wald’s Protocol (Section III-B). Tables III, IV, V
report their quantitative evaluations with respect to the quality
measures detailed in Section III-B.
TABLE IIIQ UALITY MEASURES FOR THE M OFFETT FIELD DATASET

methodSFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureCC0.92955
0.939190.938170.905210.93857
0.89580
0.916140.954960.97785
0.98168
0.97059SAM
9.52719.45999.356714.1636
11.2758
14.6132
11.3363
9.41777.13086.6392
7.6351

RMSE
365.2577
352.1290
354.8167
443.4351
363.7090
463.2204404.2979
314.4632
220.0310
200.3365254.2005

ERGAS
6.54296.0491
6.1992
7.59526.23597.92837.06195.42003.78073.42814.3582

Figures 3, 4, 5 represent the RMSEs per pixel between
the image estimated by some methods and the reference
image for the three considered datasets. Note that, for sake ofconciseness, some methods have not been considered here but
only their improved versions are presented. More specifically,
GS has been removed since GSA is an improved version of
GS. Indeed, GSA is expected to give better results than GS

TABLE IVQ UALITY MEASURES FOR THE C AMARGUE DATASET

methodSFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureCC0.91886
0.92397
0.92599
0.912620.92826
0.903500.890420.93000.95195
0.95882
0.9465SAM
4.28954.33784.28214.4982
4.19505.16374.84724.41873.64283.33453.8767

RMSE
637.1451
622.4711
611.9161
665.0173
587.1322
710.3275
745.6006
591.3190
489.5634
448.1721
511.8525

ERGAS
3.41593.2666
3.24973.54903.19403.86804.0001
3.17622.62862.47122.8181TABLE VQ UALITY MEASURES FOR THE G ARONS DATASET

methodSFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureCC0.770520.80124
0.79989
0.803470.807170.81452
0.633900.82993
0.86857
0.87834
0.86080

SAM
6.73566.61556.69056.66276.77196.63437.44156.95225.87495.63776.0224

RMSE
1036.4695
956.3047
962.1076
1037.6446
928.6229
1021.8547
1312.0373
893.9194784.1298
750.3510
778.1051

ERGAS
5.17024.82454.82805.13734.70765.01666.34164.49273.91473.76294.0454

thanks to its adaptive estimation of the weight for generating
the equivalent PAN image from the HS image, which allows
the spectral distortion to be reduced. Bayesian naive approachhas been also removed since the sparsity-based approach relies
on a more complex prior and gives also better results. MTFGLP and MTF-GLP-HPM yield similar results so only thelatter has been considered.
Figures 6 and 7 show extracts of the final result obtainedby the considered methods on the Camargue dataset in thevisible (R= 704.39nm, G= 557.90nm, B= 454.5nm) andin the SWIR (R= 1216.7nm, G= 1703.2nm, B= 2159.8nm)
domains, respectively.
Figures 8, 9 and 10 show pixel spectra recovered by thefusion methods, which correspond to 10th, 50th and 90th
percentile of RMSE, respectively. Those spectra have been
selected by choosing GSA as the reference for RMSE value.
GSA have been chosen since it is a classical approach thathas been widely used in literature and also gives good results.
To ensure a reasonable number of figures, only visual results
and some spectra of the Camargue dataset has been reported
in this article. The results for the two other datasets can befound in the supporting document [86] available online5 . In
5 http://openremotesensing.net

10

(a)

(b)

(c)

(d)

(e)

(f )

(g)

(h)

Fig. 3. RMSE between the methods’ result and the reference image, perpixel for the Moffett field dataset

(i)

(j)

Fig. 6. Details of original and fused Camargue dataset HS image in thevisible domain. (a) reference image, (b) interpolated HS image, (c) SFIM,(d) MTF-GLP-HPM, (e) GSA, (f) PCA, (g) GFPCA, (h) CNMF, (i) BayesianSparse, (j) HySure
Fig. 4. RMSE between the methods’ result and the reference image, perpixel for the Camargue dataset

Fig. 5. RMSE between the methods’ result and the reference image, perpixel for the Garons dataset

particular, because of the nature of the Garons dataset (village
with lot of small buildings) and the chosen ratio of 5, worse
results have been obtained than for the two first datasets since

a lot of mixing is presented in the HS image.
A visual analysis of the result shows that most of the fusion
approaches considered in this paper give good results, excepted
two methods: PCA and GFPCA. PCA belongs to the class of
CS methods which are known to be characterized by their highfidelity in rendering the spatial details but their generation ofsignificant spectral distortion. This is clearly visible in Figure
6 (f), where significant differences of color can be observedwith respect to the reference image, in particular when examining the different fields. GFPCA here also performs poorly.
Compared with PCA, there is less spectral distortion but theincluded spatial information seems to be not sufficient, since
the fused image is significantly blurred. Spatial information
provided by PCA is better since the main information of HS
image (where the spatial information is contained) is replacedby the high spatial information contained in the PAN image.
When using GFPCA, the guided filter controls the amount of
spatial information added to the data, so not all the spatial
information may be added to avoid to modify the spectral
information too much. For the Moffett field dataset, GFPCA
performs a little bit better since, in this dataset, there is a lotof large areas. Thus blur is less present whereas, in the Garons

11

(a)

(c)

(b)

(d)

(e)

Fig. 9. Luminance of the pixel corresponding to the 50th percentile of the
RMSE (Camargue dataset)

(f )

(g)

(i)

(h)

(j)

Fig. 7. Details of original and fused Camargue dataset HS image in the
SWIR domain. (a) reference image, (b) interpolated HS image, (c) SFIM, (d)MTF-GLP-HPM, (e) GSA, (f) PCA, (g) GFPCA, (h) CNMF, (i) BayesianSparse, (j) HySure
Fig. 10. Luminance of the pixel corresponding to the 90th percentile of the
RMSE (Camargue dataset)

Fig. 8. Luminance of the pixel corresponding to the 10th percentile of the
RMSE (Camargue dataset)dataset, GFPCA performs worse since this image consists ofnumerous small features, leading to more blurring effects. As aconsequence, in this case, GFPCA performs worse than PCA.
To analyze the spectrum in detail, chosen thanks to RMSE

percentiles, some additional information about the corresponding pixels are needed. Fig. 9 corresponds to a pixel in thereference image which represents a red building. Since in the
HS image this building is mixed with its neighborhood, we do
not have the same information between the reference image(“pure” spectrum) and the HS image (“mixed” spectrum). Fig.
8 corresponds to a pixel in a vegetation area, no mixing ispresent and very good results have been obtained for all themethods. For Fig. 10, the pixel belongs to a small building
not visible in the HS image and spectral mixing is thenalso present. More generally, spectra in the HS and reference
images differ since some mixing processes occur in the HS
image. Thus, the HS pansharpening methods are expected toprovide spectra that are more similar to the HS spectra (whichcontains the available spectral information) than the reference(which has information missing in the HS which should not
be found in the result, unless successful unmixing has been
conducted). However, it is important to note that for Fig.9,
Bayesian methods and HySure successfully recover spectrathat are more similar to the reference spectrum than the HS
spectrum.

12

Table VI report the computational times required by eachHS pansharpening methods on the Camargue dataset thosevalues have been obtained with a Intel Core i5 3230M 2.6
GHz with 8 GB RAM. Based on this table, these methods can
be classified as follows:• Methods which do not work well for HS pansharpening:
PCA, GS, GFPCA• Methods which work well with a low time of computation(few seconds): GSA, MRA methods, Bayesian Naive• Methods which work well with an average time ofcomputation (around one minute): CNMF
• Methods which work well (slightly better) with an important time of computation (few minutes, depends greatly
on the size of the dataset): Bayesian Sparse and HySure.
TABLE VIC OMPUTATIONAL TIMES OF THE DIFFERENT METHODS ( IN SECONDS )methodSFIM
MTF-GLP
MTF-GLP-HPM
GS
GSA
PCA
GFPCA
CNMF
Bayesian Naive
Bayesian Sparse
HySureMoffett
1.26
1.86
1.71
4.77
5.52
3.46
2.58
10.98
1.31
133.61140.05

Camargue
3.47
4.26
4.25
8.29
8.73
8.92
8.51
47.54
7.35
485.13296.27

Garons
2.74
4.00
2.98
5.56
5.99
6.09
4.36
23.98
3.07
259.44177.60To summarize, the comparison of the different methods
performances for RMSE curves and quality measures confirms
than PCA and GFPCA does not provide good results for
HS pansharpening (GFPCA is know to perform much betteron HS+RGB data). The other methods perform well, with
Bayesian approaches having slightly better result but with ahigher computational cost. The favorable fusion performanceobtained by the Bayesian methods can be explained, in part,by the fact that they rely on a forward modeling of the PAN
and HS images and explicitly exploit the spatial and spectral
degradations applied to the target image. However, thesealgorithms may suffer from performance discrepancies when
the parameters of these degradations (i.e., spatial blurring
kernel, sensor spectral response) are not perfectly known. In
particular, when these parameters are fully unknown and need
to be fixed, they can be estimated jointly with the fused image,as in [53], or estimated from the MS and HS images in a preprocessing step, following the strategies in [78] or [35]. CS
methods are fast to compute and easy to implement. They
provide good spatial results but poor spectral results withsignificant spectral distortions, in particular when considering
PCA and GS. GSA provides better results than the two othermethods thanks to its adaptive weight estimation reducing thespectral distortion of the equivalent PAN image created fromthe HS image. MRA methods are fast, MTF-based methods
give better results than SFIM and perform as well as the mostcompetitive algorithms with higher computational complexity.
SFIM does not perform as well than the other MRA methods
since it used a box filter which should give less good result. In
our experimentations, results from SFIM are not so different
from those obtained with the MTF-based methods. This may
come from the fact that semi-synthetic datasets are used so
MTF may not be fully used to its potential. Table VII reportthese pro and cons associated with each HS pansharpening
method.
Finally, note that, in our experimentations, no registration
error and temporal misalignment have been considered, whichsuggests that the robustness of the different methods has
not been fully analyzed. When such problems may occur,
CS and MRA methods may perform better thanks to theirgreat robustness. In particular, CS methods are robust against
misregistration error and aliasing, whereas MRA approaches
are robust against aliasing and temporal misalignment. It isalso worthy to note that the quality of a fusion method should
also be related to a specific application (such as classification
or target detection). Indeed, a method providing images withgood performance metrics may or may not be the best for thisspecific application.
V. C ONCLUSION
In this paper a qualitative and quantitative comparisonof 11 different HS pansharpening methods was conducted,
considering classical MS pansharpening techniques adapted
to the HS context, and methods originally designed for HS
pansharpening. More precisely, five classes of methods were
presented: CS, MRA, Hybrid, Bayesian and matrix factorization. Those methods were evaluated on three different datasets
representative of various scenario: mixed urban/rural area,rural area and urban area.
A careful analysis of their performances suggested a classification of these methods into four groups: i) Methods with poor
fusion results (CS-based methods and GFPCA), ii) Methods
with good fusion performances and low computational costs(MRA methods, GSA and Bayesian naive) that may be suitable for fusing large scale images, which is often the casefor spaceborne hyperspectral imaging missions, iii) Methods
with good fusion performances and reasonable computationalcosts (CNMF), iv) Methods with slightly better fusion results
but with higher computational costs (HySure and BayesianSparse). Those results were obtained with semi-synthetic
datasets with no registration error or temporal misalignment.
Thus robustness of the methods against these issues were
not taken into account. When such problems may happen,
different results could be obtained and classical pansharpening
methods (CS and MRA) may give better results thanks to theirrobustness to these specific issues.
The experiments and the quality measures presented in thispaper were performed using MATLAB implementations of thealgorithms. A MATLAB toolbox is made available online6 tothe community to help improving and developing new HS
pansharpening methods and to facilitate comparison of thedifferent methods.
6 http://openremotesensing.net

13

TABLE VIIP ROS AND C ONS OF EACH METHODS

methodpros

cons1) Reduced performance
SFIM
1) Low computational when compared to MTF
II.B.1complexity
methods (since it uses abox filter)
1) Performs well
MTF-GLP
2) Low computational
II.B.2complexity
1) Performs well
MTF-GLP-HPM
2) Low computational
II.B.2complexity
1) Spatial information
1) Low performanceis well preservedGS
for HS images
2) Low computational
II.A.22) Significant spectral
complexity
distortion
3) Easy implementation
1) Spatial information
is well preserved2) Spectral distortion
GSA
is reduced (compared
II.A.2to GS)
3) Low computational
complexity
4) Easy implementation
1) Spatial information
1) Low performanceis well preserved
PCA
for HS images
2)Low computational
II.A.1
2) Significant spectral
complexity
distortion
3) Easy implementation
1) Low performance1) Spectral information for HS images (work
GFPCA
is well preserved
better with RGB images)
II.C.1
2) Low computational 2) Not enough spatialcomplexity
information added (lotof blur)
1) Sensor characteristics
CNMF
1) Good results
required
II.E.1(spatial and spectral) 2) Medium computational
cost1) Good results
Bayesian Naive (spatial and spectral) 1) Sensor characteristics
II.D.1
2) Low computational requiredcomplexity
1) high computational
Bayesian Sparse 1) Good results
cost
II.D.2(spatial and spectral) 2) Sensor characteristics
required
HySure1) Good results1) high computational
II.D.3(spatial and spectral) cost

VI. ACKNOWLEDGMENT
The Garons and Camargue datasets were acquired in
the frame of the PRF Enviro program (internal federative project lead at ONERA). This work was partially

supported by the Fundacc ao para a Ciência e Tecnologia, Portuguese Ministry of Science and Higher Education(UID/EEA/50008/2013), project PTDC/EEI-PRO/1470/2012,and grant SFRH/BD/87693/2012. Part of this work has been
also supported by the Hypanema ANR Project n◦ ANR-12BS03-003 and by ANR-11-LABX-0040-CIMI within the program ANR-11-IDEX-0002-02. This work was supported by
the SBO-IWT project Chameleon: Domain-specific Hyperspectral Imaging Systems for Relevant Industrial Applications,and FWO project G037115N “Data fusion for image analysis
in remote sensing”. This work is supported by China Scholarship Council. This work is supported by DGA (Direction
Generale de l’Armement). This work was finally supported by
the ERC CHESS (CHallenges in Extraction and Separationof Sources) and by ANR HYEP ANR 14-CE22-0016-01(Agence National de la Recherche, Hyperspectral Imagery for
Environmental Planning).
R EFERENCES[1] [Online]. Available: http://smsc.cnes.fr/PLEIADES/index.htm[2] S. Michel, P. Gamet, and M.-J. Lefevre-Fonollosa, “HYPXIM – ahyperspectral satellite defined for science, security and defence users,”
in Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal Process.:
Evolution in Remote Sens. (WHISPERS), Lisbon, Portugal, June 2011,
pp. 1–4.
[3] L. Alparone, L. Wald, J. Chanussot, C. Thomas, P. Gamba, and L. M.
Bruce, “Comparison of pansharpening algorithms: Outcome of the 2006
GRS-S data fusion contest,” IEEE Trans. Geosci. and Remote Sens.,
vol. 45, no. 10, pp. 3012–3021, Oct. 2007.
[4] G. Vivone, L. Alparone, J. Chanussot, M. Dalla Mura, Garzelli, and
G. Licciardi, “A critical comparison of pansharpening algorithms,” inProc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS), July 2014, pp.
191–194.
[5] B. Aiazzi, L. Alparone, S. Baronti, A. Garzelli, and M. Selva, “25 yearsof pansharpening: a critical review and new developments,” in Signal
and Image Processing for Remote Sensing, 2nd ed., C. H. Chen, Ed.
Boca Raton, FL: CRC Press, 2011, ch. 28, pp. 533–548.
[6] C. Thomas, T. Ranchin, L. Wald, and J. Chanussot, “Synthesis ofmultispectral images to high spatial resolution: a critical review of fusion
methods based on remote sensing physics,” IEEE Trans. Geosci. andRemote Sens., vol. 46, no. 5, pp. 1301–1312, May 2008.
[7] W. Carper, T. M. Lillesand, and P. W. Kiefer, “The use of intensityhue-saturation transformations for merging SPOT panchromatic andmultispectral image data,” Photogramm. Eng. Remote Sens., vol. 56,
no. 4, pp. 459–467, April 1990.
[8] T.-M. Tu, S.-C. Su, H.-C. Shyu, and P. S. Huang, “A new look at IHS-like
image fusion methods,” Information Fusion, vol. 2, no. 3, pp. 117–186,
Sept. 2001.
[9] P. S. Chavez Jr., S. C. Sides, and J. A. Anderson, “Comparison of
three different methods to merge multiresolution and multispectral data:
Landsat TM and SPOT panchromatic,” Photogramm. Eng. Remote Sens.,
vol. 57, no. 3, pp. 295–303, March 1991.
[10] P. S. Chavez and A. Y. Kwarteng, “Extracting spectral contrast inlandsat thematic mapper image data using selective principal component
analysis,” Photogramm. Eng. Remote Sens., vol. 55, no. 3, pp. 339–348,
1989.
[11] V. Shettigara, “A generalized component substitution technique forspatial enhancement of multispectral images using a higher resolution
data set,” Photogramm. Eng. Remote Sens., vol. 58, no. 5, pp. 561–567,
1992.
[12] V. P. Shah, N. Younan, and R. L. King, “An efficient pan-sharpening
method via a combined adaptative-PCA approach and contourlets,” IEEE
Trans. Geosci. and Remote Sens., vol. 56, no. 5, pp. 1323–1335, May
2008.
[13] C. Laben and B. Brower, “Process for enhacing the spatial resolution ofmultispectral imagery using pan-sharpening,” U.S. Patent US6 011 875,
2000.
[14] S. Mallat, “A theory for multiresolution signal decomposition: thewavelet representation,” IEEE Trans. Patt. Anal. Mach. Intell., vol. 11,
no. 7, pp. 674–693, July 1989.

14[15] G. P. Nason and B. W. Silverman, “The stationary wavelet transform andsome statistical applications,” in Wavelets and Statistics, A. Antoniadis
and G.Oppenheim, Eds. New York, NY: Springer-Verlag, 1995, vol.
103, pp. 281–299.
[16] M. J. Shensa, “The discrete wavelet transform: wedding the a trous andMallat algorithms,” IEEE Trans. Signal Process., vol. 40, no. 10, pp.
2464–2482, Oct. 1992.
[17] P. J.Burt and E. H. Adelson, “The Laplacian pyramid as a compact
image code,” IEEE Trans. Comm., vol. 31, no. 4, pp. 532–540, April
1983.
[18] M. N. Do and M. Vetterli, “The contourlet transform: an efficient
directional multiresolution image representation,” IEEE Trans. Image
Process., vol. 14, no. 12, pp. 2091–2106, Dec. 2005.
[19] J. F. J.-L. Starck and F. Murtagh, “The undecimated wavelet decomposition and its reconstruction,” IEEE Trans. Image Process., vol. 16, no. 2,
pp. 297–309, Feb. 2007.
[20] C. Ballester, V. Caselles, L. Igual, J. Verdera, and B. Rougé, “Avariational model for P+XS image fusion,” Int. J. Computer Vision, vol.
5969, no. 1, pp. 43–58, 2006.
[21] F. Palsson, J. Sveinsson, M. Ulfarsson, and J. A. Benediktsson, “A new
pansharpening algorithm based on total variation,” IEEE Geosci. andRemote Sensing Lett., vol. 11, no. 1, pp. 318–322, Jan. 214.
[22] X. He, L. Condat, J. Bioucas Dias, J. Chanussot, and J. Xia, “A new
pansharpening method based on spatial and spectral sparsity priors,”
IEEE Trans. Image Process., vol. 23, no. 9, pp. 4160–4174, Sept 2014.
[23] M. Moeller, T. Wittman, and A. L. Bertozzi, “A variational approachto hyperspectral image fusion,” in SPIE Defense, Security, and Sensing,
2009.
[24] A. Garzelli, B. Aiazzi, S. Baronti, M. Selva, , and L. Alparone,“Hyperspectral image fusion,” in Proc. Hyperspectral Workshop, 2010,
pp. 17–19.
[25] L. Alparone, B. Aiazzi, S. Baronti, and A. Garzelli, Remote Sensing
Image Fusion, ser. Signal and Image Processing of Earth Observations.
Boca Raton, FL: CRC Press, 2015.
[26] G. Vivone, L. Alparone, J. Chanussot, M. Dalla Mura, Garzelli, and
G. Licciardi, “Multi-resolution analysis and component substitution
techniques for hyperspectral pansharpening,” in Proc. IEEE Int. Conf.
Geosci. Remote Sens. (IGARSS), July 2014, pp. 2649–2652.
[27] J. M. Bioucas Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader,and J. Chanussot, “Hyperspectral unmixing overview: Geometrical,
statistical, and sparse regression-based approaches,” IEEE J. Sel. Topics
Appl. Earth Observ. Remote Sens., vol. 5, no. 2, pp. 354–379, Apr. 2012.
[28] C. Souza Jr, L. Firestone, L. M. Silva, and D. Roberts, “Mapping
forest degradation in the Eastern Amazon from SPOT 4 through spectral
mixture models,” Remote Sens. Environment, vol. 87, no. 4, pp. 494–506,
2003, large Scale Biosphere Atmosphere Experiment in Amazonia.
[29] A. Mohammadzadeh, A. Tavakoli, and M. J. V. Zoej, “Road extraction
based on fuzzy logic and mathematical morphology from pansharpened
IKONOS images,” The Photogrammetric Record, vol. 21, no. 113, pp.
44–60, Feb. 2006.
[30] F. Laporterie-Déjean, H. de Boissezon, G. Flouzat, and M.-J. LefèvreFonollosa, “Thematical and statistical evaluations of five panchromatic/multispectral fusion methods on simulated PLEIADES-HR image,” Information Fusion, vol. 6, pp. 193–212, 2005.
[31] G. A. Licciardi, A. Villa, M. M. Khan, and J. Chanussot, “Image fusionand spectral unmixing of hyperspectral images for spatial improvementof classification maps,” in Proc. IEEE Int. Conf. Geosci. Remote Sens.
(IGARSS), 2012, pp. 7290–729.
[32] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix
factorization unmixing for hyperspectral and multispectral data fusion,”
IEEE Trans. Geosci. and Remote Sens., vol. 50, no. 2, pp. 528–537,
Feb. 2012.
[33] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, “Bayesian fusion of multiband images,” IEEE J. Sel. Topics Signal Process., 2015, to appear.
[34] ——, “Bayesian fusion of hyperspectral and multispectral images,” inProc. IEEE Int. Conf. Acoust., Speech, and Signal Processing (ICASSP),
Florence, Italy, May 2014, pp. 3176–3180.
[35] M. Simões, J. Bioucas Dias, L. Almeida, and J. Chanussot, “A convex
formulation for hyperspectral image superresolution via subspace-based
regularization,” IEEE Trans. Geosci. and Remote Sens., 2015, to appear.
[36] Q. Wei, J. M. Bioucas Dias, N. Dobigeon, and J.-Y. Tourneret, “Hyperspectral and multispectral image fusion based on a sparse representation,” IEEE Trans. Geosci. and Remote Sens., vol. 53, no. 7, pp. 3658–
3668, Sept. 2015.
[37] B. Aiazzi, S. Baronti, and M. Selva, “Improving component substitution
Pansharpening through multivariate regression of MS+Pan data,” IEEE

Trans. Geosci. and Remote Sens., vol. 45, no. 10, pp. 3230–3239, Oct.
2007.
[38] S. Baronti, B. Aiazzi, M. Selva, A. Garzelli, and L. Alparone, “Atheoretical analysis of the effects of aliasing and misregistration onpansharpened imagery,” IEEE J. Sel. Topics Signal Process., vol. 5, no. 3,
pp. 446–453, June 2011.
[39] G. Vivone, L. Alparone, J. Chanussot, M. Dalla Mura, Garzelli, and
G. Licciardi, “A critical comparison among pansharpening algorithms,”
IEEE Trans. Geosci. and Remote Sens., vol. 53, no. 5, pp. 2565–2586,
may 2015.
[40] T.-M. Tu, P. S. Huang, C.-L. Hung, and C.-P. Chang, “A fast intensityhue-saturation fusion technique with spectral adjustment for IKONOS
imagery,” IEEE Geosci. and Remote Sensing Lett., vol. 1, no. 4, pp.
309–312, 2004.
[41] B. Aiazzi, S. Baronti, F. Lotti, and M. Selva, “A comparison between
global and context-adaptive pansharpening of multispectral images,”
IEEE Geosci. and Remote Sensing Lett., vol. 6, no. 2, pp. 302–306,
April 2009.
[42] G. Vivone, R. Restaino, M. Dalla Mura, G. Licciardi, and J. Chanussot,“Contrast and error-based fusion schemes for multispectral image pansharpening,” IEEE Geosci. and Remote Sensing Lett., vol. 11, no. 5, pp.
930–934, May 2014.
[43] J. G. Liu, “Smoothing filter based intensity modulation: a spectral
preserve image fusion technique for improving spatial details,” Int. J.
Remote Sens., vol. 21, no. 18, pp. 3461–3472, Dec. 2000.
[44] L. Alparone, B. Aiazzi, S. Baronti, and A. Garzelli, “Sharpening ofvery high resolution images with spectral distortion minimization,” inProc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS), July 2003, pp.
458–460.
[45] B. Aiazzi, L. Alparone, S. Baronti, A. Garzelli, and M. Selva, “MTFtailored multiscale fusion of high-resolution MS and Pan imagery,”Photogramm. Eng. Remote Sens., vol. 72, no. 5, pp. 591–596, May
2006.
[46] K. He, J. Sun, and X. Tang, “Guided image filtering,” IEEE Trans. Patt.
Anal. Mach. Intell., vol. 35, no. 6, pp. 1397–1409, 2013.
[47] X. Kang, J. Li, and J. A. Benediktsson, “Spectral-spatial hyperspectral
image classification with edge-preserving filtering,” IEEE Trans. Geosci.
and Remote Sens., vol. 52, no. 5, pp. 2666–2677, 2014.
[48] W. Liao, X. Huang, F. Coillie, S. Gautama, A. Pizurica, W. Philips,
H. Liu, T. Zhu, M. Shimoni, G. Moser, and D. Tuia, “Processing ofmultiresolution thermal hyperspectral and digital color data: Outcomeof the 2014 IEEE GRSS data fusion contest,” IEEE J. Sel. Topics Appl.
Earth Observ. Remote Sens., Submitted.
[49] W. Liao, B. Goossens, J. Aelterman, H. Luong, A. Pizurica, N. Wouters,
W. Saeys, and W. Philips, “Hyperspectral image deblurring with pca andtotal variation,” in Proc. IEEE GRSS Workshop Hyperspectral Image
SIgnal Process.: Evolution in Remote Sens. (WHISPERS), Florida, US,
June 2013, pp. 1–4.
[50] 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, Sept. 2004.
[51] Y. Zhang, S. De Backer, and P. Scheunders, “Noise-resistant waveletbased Bayesian fusion of multispectral and hyperspectral images,” IEEE
Trans. Geosci. and Remote Sens., vol. 47, no. 11, pp. 3834 –3843, Nov.
2009.
[52] M. Joshi and A. Jalobeanu, “MAP estimation for multiresolution fusion
in remotely sensed images using an IGMRF prior model,” IEEE Trans.
Geosci. and Remote Sens., vol. 48, no. 3, pp. 1245–1255, March 2010.
[53] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, “Bayesian fusion of multispectral and hyperspectral images with unknown sensor spectral response,” in Proc. IEEE Int. Conf. Image Processing (ICIP), Paris,
France, Oct. 2014, pp. 698–702.
[54] M. Simões, J. Bioucas Dias, L. B. Almeida, and J. Chanussot, “Hyperspectral image superresolution: An edge-preserving convex formulation,”in Proc. IEEE Int. Conf. Image Processing (ICIP), Paris, France, Oct.
2014, pp. 4166–4170.
[55] R. Molina, A. K. Katsaggelos, and J. Mateos, “Bayesian and regularization methods for hyperparameter estimation in image restoration,” IEEE
Trans. Image Process., vol. 8, no. 2, pp. 231–246, 1999.
[56] R. Molina, M. Vega, J. Mateos, and A. K. Katsaggelos, “Variational
posterior distribution approximation in Bayesian super resolution reconstruction of multispectral images,” Applied and Computational Harmonic Analysis, vol. 24, no. 2, pp. 251 – 267, 2008.
[57] A. K. Gupta and D. K. Nagar, Matrix Variate Distributions, ser.
Monographs and surveys in pure and applied mathematics. Boca Raton,
FL: Chapman & Hall/CRC, 2000, vol. 104.

15[58] M. D. Farrell Jr and R. M. Mersereau, “On the impact of PCA dimension
reduction for hyperspectral detection of difficult targets,” IEEE Geosci.
and Remote Sensing Lett., vol. 2, no. 2, pp. 192–195, 2005.
[59] J. M. Nascimento and J. M. Bioucas Dias, “Vertex component analysis:
a fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci. andRemote Sens., vol. 43, no. 4, pp. 898–910, April 2005.
[60] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, “Bayesian fusion of multispectral and hyperspectral images using a block coordinate descentmethod,” in Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal
Process.: Evolution in Remote Sens. (WHISPERS), Tokyo, Japan, June
2015.
[61] ——, “Fast fusion of multi-band images based on solving a
Sylvester equation,” 2015, submitted. [Online]. Available: http://arxiv.org/abs/1502.03121/[62] S. Mallat, A wavelet tour of signal processing. New York: Academic
Press, 1999.
[63] J.-L. Starck, E. Candes, and D. Donoho, “The curvelet transform forimage denoising,” IEEE Trans. Image Process., vol. 11, no. 6, pp. 670–
684, 2002.
[64] N. Ahmed, T. Natarajan, and K. Rao, “Discrete cosine transform,” IEEE
Trans. Computers, vol. C-23, no. 1, pp. 90–93, 1974.
[65] M. Elad and M. Aharon, “Image denoising via sparse and redundant
representations over learned dictionaries,” IEEE Trans. Image Process.,
vol. 15, no. 12, pp. 3736–3745, 2006.
[66] O. G. Guleryuz, “Nonlinear approximation based image recovery using
adaptive sparse reconstructions and iterated denoising – Part I: theory,”
IEEE Trans. Image Process., vol. 15, no. 3, pp. 539–554, 2006.
[67] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, andD. B. Rubin, Bayesian data analysis, 3rd ed. Boca Raton, FL: CRC
press, 2013.
[68] X. Bresson and T. Chan, “Fast dual minimization of the vectorial total
variation norm and applications to color image processing,” Inverse
Probl. and Imag., vol. 2, no. 4, pp. 455–484, 2008.
[69] M. Afonso, J. Bioucas Dias, and M. Figueiredo, “An augmented
Lagrangian approach to the constrained optimization formulation ofimaging inverse problems.” IEEE Trans. Image Process., vol. 20, no. 3,
pp. 681–95, 2011.
[70] O. Berné, A. Tielens, P. Pilleri, and C. Joblin, “Non-negative matrix
factorization pansharpening of hyperspectral data: An application to
mid-infrared astronomy,” in Proc. IEEE GRSS Workshop Hyperspectral
Image SIgnal Process.: Evolution in Remote Sens. (WHISPERS), 2010,
pp. 1–4.
[71] R. Kawakami, J. Wright, Y. Tai, Y. Matsushita, M. Ben-Ezra, andK. Ikeuchi, “High-resolution hyperspectral imaging via matrix factorization,” in Proc. IEEE Int. Conf. Comp. Vision and Pattern Recognition(CVPR), 2011, pp. 2329–2336.
[72] A. Charles, B. Olshausen, and C. Rozell, “Learning sparse codes for
hyperspectral imagery,” IEEE J. Sel. Topics Signal Process., vol. 5, no. 5,
pp. 963–978, 2011.
[73] B. Huang, H. Song, H. Cui, J. Peng, and Z. Xu, “Spatial and spectral
image fusion using sparse matrix factorization,” IEEE Trans. Geosci.
and Remote Sens., vol. 52, no. 3, pp. 1693–1704, 2014.
[74] M. Veganzones, M. Simes, G. Licciardi, J. M. Bioucas Dias, and
J. Chanussot, “Hyperspectral super-resolution of locally low rank images
from complementary multisource data,” in Proc. IEEE Int. Conf. Image
Processing (ICIP), Paris, France, Oct. 2014, pp. 703–707.
[75] N. Yokoya, T. Yairi, and A. Iwasaki, “Hyperspectral, multispectral,and panchromatic data fusion based on coupled non-negative matrix
factorization,” in Proc. IEEE GRSS Workshop Hyperspectral Image
SIgnal Process.: Evolution in Remote Sens. (WHISPERS), 2011, pp. 1–4.
[76] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol. 401, pp. 788–791, 1999.
[77] D. C. Heinz and C. -I Chang, “Fully constrained least-squares linear
spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Trans. Geosci. and Remote Sens., vol. 29, no. 3,
pp. 529–545, March 2001.
[78] N. Yokoya, N. Mayumi, and A. Iwasaki, “Cross-calibration for data
fusion of EO-1/Hyperion and Terra/ASTER,” IEEE J. Sel. Topics Appl.
Earth Observ. Remote Sens., vol. 6, pp. 419–2013, April 2013.
[79] L. Wald, T. Ranchin, and M. Mangolini, “Fusion of satellite images ofdifferent spatial resolutions: Assessing the quality of resulting image,”
IEEE Trans. Geosci. and Remote Sens., vol. 43, pp. 1391–1402, 2005.
[80] I. Amro, J. Mateos, M. Vega, R. Molina, and A. K. Katsaggelos,“A survey of classical methods and new trends in pansharpening ofmultispectral images,” EURASIP J. Adv. Signal Process., vol. 2011,
no. 1, pp. 1–22, 2011.

[81] Q. Du, N. H. Younan, R. L. King, and V. P. Shah, “On the performance
evaluation of pan-sharpening techniques,” IEEE Geosci. and Remote
Sensing Lett., vol. 4, pp. 518–22, Oct. 2007.
[82] Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE
Signal Process. Lett., vol. 9, no. 3, pp. 81–84, March 2002.
[83] L. Alparone, B. Aiazzi, S. Baronti, A. Garzelli, F. Nencini, and M. Selva,“Multispectral and panchromatic data fusion assessment without reference,” Photogramm. Eng. Remote Sens., vol. 74, no. 2, pp. 193–200,
Feb. 2008.
[84] G. Piella and H.Heijmans, “A new quality metric for image fusion,”in Proc. IEEE Int. Conf. Image Processing (ICIP), vol. 2, 2003, pp.
173–176.
[85] L. Wald, Data Fusion : Definitions and Architectures - Fusion of imagesof different spatial resolutions. Les Presses de l’Ecole des Mines, 2002.
[86] L. Loncan, L. B. Almeida, J. M. Bioucas-Dias, X. Briottet,
J. Chanussot, N. Dobigeon, S. Fabre, W. Liao, G. Licciardi, M. Simoes,
J.-Y. Tourneret, M. Veganzones, G. Vivone, Q. Wei, and N. Yokoya,“Hyperspectral pansharpening: a review – Complementary results
and supporting materials,” Tech. Rep., 2015. [Online]. Available:
http://openremotesensing.net