GENERALIZED LINEAR MIXING MODEL ACCOUNTING FOR ENDMEMBER
VARIABILITY
Tales Imbiriba, Ricardo Augusto Borsoi, José Carlos Moreira Bermudez
Federal University of Santa Catarina, Florianópolis, SC, Brazil

arXiv:1710.07723v1 [] 20 Oct 2017

ABSTRACT
Endmember variability is an important factor for accurately unveiling vital information relating the pure materials and their distribution
in hyperspectral images. Recently, the extended linear mixing model
(ELMM) has been proposed as a modification of the linear mixing
model (LMM) to consider endmember variability effects resulting
mainly from illumination changes. In this paper, we further generalize the ELMM leading to a new model (GLMM) to account for
more complex spectral distortions where different wavelength intervals can be affected unevenly. We also extend the existing methodology to jointly estimate the variability and the abundances for the
GLMM. Simulations with real and synthetic data show that the unmixing process can benefit from the extra flexibility introduced by
the GLMM.
Index Terms— Hyperspectral data, GLMM, endmember variability.
1. INTRODUCTION
Hyperspectral imaging has attracted formidable interest of the scientific community in the past two decades, where hyperspectral images
(HIs) have been explored in an increasing number of applications in
different fields [1]. The limited spatial resolution of hyperspectral
devices often mixes the spectral contribution of different pure materials, termed endmembers, in the scene [2]. This phenomenon is
more prominent in remote sense applications due to the distance between airborne and spaceborne sensors and the target scene. The
mixing process must be well understood to accurately unveil vital
information relating the presence of pure materials and their distribution in the scene. Hyperspectral unmixing (HU) aims to solve this
problem by decomposing the hyperspectral image (HI) into a collection of endmembers and their fractional abundances [3].
Different mixing models have been proposed to explain the interaction between light and the endmembers [4]. The simplest and
most widely used model is the Linear Mixing Model (LMM) [2].
The LMM assumes that the observed reflectance vector (i.e. a HI
pixel) can be modeled as a convex combination of the spectral signatures of the endmembers present in the scene. This assumption
imposes positivity and sum-to-one constraints on the linear combination coefficients, leading to their interpretation as fractional abundances of the corresponding endmembers in each pixel. The simplicity of the LMM and the convexity constraints on the fractional
abundances naturally lead to fast and reliable unmixing strategies.
This work has been supported by the National Council for Scientific and
Technological Development (CNPq).
T. Imbiriba, R. A. Borsoi and J.C.M. Bermudez are with the Department of Electrical Engineering, Federal University of Santa Catarina,
Florianópolis, SC, Brazil. e-mail: talesim@gmail.com; raborsoi@ucs.br;
j.bermudez@ieee.org.

Given the intrinsic limitations of this simple mixing model, more
refined unmixing approaches have been more recently proposed to
account nonideal effects such as nonlinearity [4–6] and endmember
variability, often found in a typical scenes [7, 8].
Endmember variability can be caused by a myriad of factors
including environmental, illumination, atmospheric and temporal
changes [7], and its occurrence may result in significant estimation
errors being propagated throughout the unmixing process [9]. The
most common approaches to deal with spectral variability can be
divided in three basic classes. 1) to group endmembers in variational sets, 2) to model endmembers as statistical distributions, and
3) to incorporate the variability in the mixing model, often using
physically motivated concepts [8]. This work follows the third approach. Recently, [9] and [10] introduced variations of the LMM
to cope with spectral variability. The model called Perturbed LMM
model (PLMM) introduced in [9] defined an additive perturbation
to the endmember matrix. Such perturbation matrix then needs to
be estimated jointly with the abundances. Though the perturbation
matrix can model arbitrary endmember variations, it lacks physical
motivation. The Extended Linear Mixing Model (ELMM) proposed
in [10] extended the LMM model by using one pixel-dependent
multiplicative term for each endmember, a generalization that can
efficiently model changes in the observed reflectances due to illumination, an important effect [10]. This model addresses a physically
motivated problem, with the advantage of estimating a variability
parameter vector of much lower dimension when compared with the
additive perturbation matrix in PLMM. Although the ELMM performs well in situations where spectral variability is mainly caused
by illumination variations, it lacks a necessary flexibility when the
endmembers are subject to more complex spectral distortions. For
instance, experimental measurements on vegetation spectra under
different conditions have shown a significant dependence of the
spectral variation on wavelength intervals [11–13]. This type of
variability is not supported by the ELMM model in [10], which
assumes a fixed scaling across all wavelengths.
In this work we introduce a generalization of the ELMM model
proposed in [10] to accounts for endmember variability in arbitrary
regions of the measured spectrum. We call the resulting model the
Generalized Linear Mixing Model (GLMM). The estimation of the
required parameters is realized by generalizing the methodology
used in [10] through the use of three-dimensional tensors to accommodate the new model without significantly effecting the simplicity
of the proposed solution or its computational complexity. Simulation results using synthetic and real data indicate that the extra
flexibility introduced by the GLMM model can improve the results
of existing methods for different types of endmember variability.
This paper is organized as follows. In Section 2 we briefly review the LMM and its extended version ELMM. Section 3 introduces the proposed GLMM. In Section 4 we define new tensor variables and extend the solution in [10] to the GLMM. The performance

of the proposed method is compared with competing algorithms in
Section 5. Finally, the conclusions are presented in Section 6.

functional:
J(A, M, Ψ) =

2. EXTENDED LINEAR MIXING MODEL

N
1X
kr n − M n αn k2
2 n=1

+λM kM n − M 0
The Linear Mixing Model (LMM) [2] assumes that a given a pixel
r n = [rn,1 , . . . , rn,L ]> , with L bands, is represented as
r n = M αn + en
subject to 1> αn = 1 and αn  0

(1)

where M = [m1 , . . . , mR ] is an L × R matrix whose columns are
the R endmember spectral signatures mi = [mi,1 , . . . , mi,L ]> ,
αn = [αn,1 , . . . , αn,R ]> is the abundance vector, e ∼ N (0, σn2 I L )
is an additive white Gaussian noise (WGN), I L is the L × L identity
matrix, and  is the entrywise ≥ operator. The LMM assumes that
the endmember spectra are fixed for all pixels r n , n = 1, . . . , N ,
in the HI. This assumption can jeopardize the accuracy of estimated
abundances in many circumstances due to the spectral variability
existing in a typical scene. The Extended Linear Mixing Model
(ELMM) [10] partially mitigates such limitation by including a multiplicative diagonal weight matrix ψ n = diag(ψ1,n , . . . , ψR,n ) in
the LMM such that
r n = M ψ n αn + en



(4)

+ R(A) + R(Ψ).
where M and Ψ are L × R × N tensors with entries [M]:,:,n = M n ,
and [Ψ]:,:,n = Ψn respectively, A = [α1 , . . . , αN ] is the abundance matrix, λM is a parameter that controls the strictness of the
regularization over M n , and R(A) and R(Ψ) are spatial regularizations over A and Ψ. Thus, the optimization problem becomes
(A∗ , M∗ , Ψ∗ ) = arg min J(A, M, Ψ).

(5)

A, M, Ψ

The problem defined in (5) is non-smooth and non-convex with respect to all variables A, M, and Ψ, but is convex with respect to each
one of them. Thus, we follow the same approach used in [10] and
find a local stationary point minimizing (5) iteratively with respect
to each variable, leading to the strategy presented in Algorithm 1.
Algorithm 1: Global algorithm for solving (4)

4

Input : R, λM , λA , λΨ , A(0) , Ψ(0) and M 0 .
Output: A∗ , M∗ and Ψ∗ .
Set i = 0 ;
while stopping criterion is not satisfied do
i=i+1;
M(i) = arg min J(A(i−1) , M, Ψ(i−1) ) ;

5

A(i) = arg min J(A, M(i) , Ψ(i−1) ) ;

6

Ψ(i) = arg min J(A(i) , M(i) , Ψ) ;

(2)

with ψk,n ∈ R+ , k = 1, . . . , R. Each coefficient ψk,n scales the
whole spectrum of endmember mi in pixel n, leading to a simple
strategy to model variability resulting from illumination effects.

Ψn k2F

1
2
3

M

3. GENERALIZED LINEAR MIXING MODEL (GLMM)

A
Ψ

As explained in Section 1, we propose a generalization of the
ELMM model to allow for spectral variabilities per wavelength intervals. To this end, we propose to employ a band-dependent scaling
factor, enabling the new model to adapt to arbitrary variations of
the endmember spectra. In the new GLMM model each pixel r n is
written as
r n = (M Ψn )αn + en
(3)
where Ψn ∈ RL×R is a scaling matrix with entries [Ψn ]`,k =
ψn`,k ≥ 0, and is the Hadamard product. This model is a generalization of the ELMM where the scaling matrix Ψn acts on each
wavelength of each endmember individually. Such feature leads to
a more flexible model that allows to consider variabilities that are
not uniform along each endmember spectrum. ELMM is clearly a
particular case of GLMM, and the new model can be employed for
any level of granularity of variability per wavelength ranges, to the
limit of an independent scaling of each wavelength component of
each endmember in each pixel.
4. THE UNIMIXING PROBLEM
Assuming the availability of a reference endmember matrix M 0
(which can be obtained using any endmember extraction method),
the HU problem reduces to estimating the free parameters minimizing a given risk functional defined for the whole HI R =
[r 1 , . . . , r N ]. For this purpose the methodology presented in [10]
can be extended for the GLMM by defining three-dimensional tensors. Thus, we propose to minimize the following regularized cost

7
8

end
return A∗ = A(i) , M∗ = M(i) , Ψ∗ = Ψ(i) ;

The regularization functionals R(A) and R(Ψ) in (4) are selected in order to provide spatial smoothness to the abundances and
scaling factors and to enforce physical constraints (i.e. positivity and
sum-to-one constraints on the abundances). They are selected as

R(A) = λA kHh (A)k2,1 + kHv (A)k2,1
+ ιR+ (A) + µ> (A> 1R×1 − 1N ×1 )
and
R(Ψ) =

L
R
λΨ X X
(kHh ([Ψ]`,k,: )k2 + kHv ([Ψ]`,k,: )k2 )
2
`=1 k=1

where [·]`,k,: is a slice of a tensor for band `, endmember k and all
N pixels, and the parameters λA and λΨ control the weights of the
regularization terms in the cost function. The linear operators Hh
and Hv compute the first-order horizontal and vertical gradients of a
bidimensional signal, acting separately for each material of A. The
spatial regularization in the abundances is promoted by a mixed L2,1
P
norm of their gradient, where kXk2,1 = N
n=1 kxn k2 . This norm
is used to promote sparsity of the gradient across different materials
(i.e. to force neighboring pixels to be homogeneous in all constituent
endmembers). The L1 norm can also beP
used, leading to the Total
Variation regularization, where kXk1 = N
n=1 kxn k1 [14]. ιR+ (·)

is the indicator function of R+ (i.e. ιR+ (a) = 0 if a ≥ 0 and
ιR+ (a) = ∞ if a < 0) acting component-wise on its input, and
enforces the abundances positivity constraint. µ ∈ RN is a vector
of Lagrange multipliers associated with the sum-to-one constraint.
Note that the cost function must also be optimized with respect to µ
for this constraint to be enforced.
4.1. Optimiztion with respect to M
Rewriting the problem (5) using only the terms in (4) that depend on
M, the problem becomes
N
1X
kr n − M n αn k2
2 n=1
M0

−λM kM n − M 0 Ψn k2F .

M∗ = arg min

(6)

The problem in (6) can be solved individually for each pixel r n .
Thus, relaxing the positivity constraint on the elements of M, the
solution can be found as
M ∗n = (r n α> + λM M 0

−1
Ψn )(αn α>
n + λs I R )

(7)

where I R is the R × R identity matrix. Then, an approximate solution to the constrained problem can be obtained by projecting M ∗n
onto the nonnegative orthant RL×R
by thresholding the negative en+
tries to zero [10].
4.2. Optimization with respect to A
Restating the problem (5) only considering the terms in (4) that depend on A leads to the following optimization problem for the abundance matrix
A∗ = arg min
A

N
1X
kr n − M n αn k2
2 n=1

+ λA (kHh (A)k22,1 + kHv (A)k22,1 )

(8)

+ ιR+ (A) + µ> (A> 1R×1 − 1N ×1 )
This problem is clearly not separable with respect to the pixels in
the image. However, problem (8) can be efficiently solved using the
Alternating Direction Method of the Multipliers (ADMM) [15]. The
procedure is well described in [10] and will be suppressed here for
conciseness.
4.3. Optimization with respect to Ψ

Ψ∗ = arg min
Ψ

λM
2

kM n − M 0

5. SIMULATIONS
In this section, the performance of the proposed methodology is illustrated through simulations with both synthetic and real data. We
compare the proposed method based on the GLMM with the the
fully constrained least squares (FCLS), the scaled constrained least
squares (SCLS), the PLMM [9], and the ELMM [10].
To measure the accuracy of the unmixing methods we consider
the Root Means Squared Error (RMSE)
r
1
kvec(X) − vec( X ∗ )k2
(12)
RMSEX =
NX
where vec(·) is the vectorization operator X → x, Ra×b×c 7→
Rabc , NX = abc. In this work we apply the RMSE to evaluate
the estimates of the abundances (RMSEA ), of the endmembers tensor (RMSEM ) and of the reconstructed images (RMSER ). We also
consider the Spectral Angle Mapper for the HI
SAMR



N
∗
1 X
r>
n rn
=
arccos
N n=1
kr n kkr ∗n k

(13)

N
R
1 XX
SAMM =
arccos
N n=1

∗
m>
k,n mk,n
kmk,n kkm∗k,n k

!
.

(14)

Ψn k2F

n=1

(9)

which can be rewritten as

Ψ

where F and F −1 represents the bi-dimensional discrete Fourier
transform and its inverse respectively, 1p×q is a p × q matrix of
ones, where p and q are the number of rows and columns of the HI
cube, and hh and hv are convolution masks represented as p × q
matrices.

k=1

λΨ
+
(kHh (Ψ)k2F + kHv (Ψ)k2F )
2

Ψ∗ = arg min

>
where P = λM (m0`,k )2 I N + λΨ (H >
h H h + H v H v ). The solution in (11) involves the inverse of the N × N matrix P which can
be computationally intensive or intractable. However, if we assume
periodic boundary conditions for the differential operators Hh and
Hv , the corresponding matrices H v , H h and consequently P will
have the structure of a block circulant matrix with circulant blocks
(BCCB). Since BCCB matrices can be diagonalized using the bidimensional Discrete Fourier Transform, problem (11) can be solved
efficiently as follows [16]
!
F(λM m0`,k [M]`,k,: )
∗
−1
[Ψ]`,k,: = F
λM (m0`,k )2 1p×q + λΨ (|F(hh )|2 + |F(hv )|2 )

and for the endmembers tensor

Rewriting the optimization problem (5) considering only the terms
in (4) that depend on Ψ leads to
N
X

where [·]`,k,: is a slice of a tensor for band `, endmember k and all N
pixels and m0`,k = [M 0 ]`,k is a scalar. The problem can be solved
for each endmember k and band ` individually, and its solution is
given by

[Ψ]∗`,k,: = P −1 λM m0`,k [M]`,k,:
(11)

L
R

λM X X
k[M]`,k,: − m0`,k [Ψ]`,k,: k2
2
`=1 k=1

λΨ
+
(kHh ([Ψ]`,k,: )k2 + kHv ([Ψ]`,k,: )k2 )
2

(10)

5.1. Synthetic data
For a comprehensive comparison among the different methods we
created two synthetic datasets, namely Data Cube 0 (DC0) and Data
Cube 1 (DC1), represented in Figure 1. These datasets were built
using endmembers extracted from the USGS Spectral Library [17],
and different strategies were used to generate the abundance maps,
which exhibit spatial correlation between neighboring pixels. For
DC0, we adopted the variability model used in [10] (a multiplicative factor acting in each endmember), while for DC1 we considered

Table 1. Simulations with synthetic data.
Data Cube 0 – DC0
FCLS
SCLS
PLMM
ELMM
GLMM

Data Cube 1 – DC1

Houston Data

RMSEA

RMSEM

SAMM

RMSER

SAMR

RMSEA

RMSEM

SAMM

RMSER

SAMR

RMSER

SAMR

0.0968
0.0642
0.0641
0.0540
0.0512

0.0673
0.0689
0.0608
0.0587

0.0625
0.0566
0.0568
0.0601

0.0420
0.0403
0.0250
0.0254
0.0011

0.0566
0.0548
0.0341
0.0346
0.0011

0.0243
0.0509
0.0476
0.0209
0.0193

0.0457
0.0401
0.0425
0.0391

0.0617
0.0578
0.0609
0.0564

0.0385
0.0383
0.0257
0.0367
0.0202

0.0492
0.0489
0.0329
0.0469
0.0257

0.006082
0.006082
0.002918
0.003217
0.000313

0.057695
0.025361
0.013617
0.013225
0.000472

Met. Roofs

Concrete

Asphalt

SCLSU

FCLS

Vegetation

5.2. Real data
For simulations with real data we considered the Houston dataset
discussed in [10]. This dataset is known to have four endmembers
which were extracted using the VCA algorithm [18]. Figure 2 shows
the reconstructed abundance maps for all tested methods while Table 1 presents the results in terms of RMSER and SAMR . Figure 2
shows that the proposed GLMM method provided smooth and accurate abundance estimation, comparable with the results obtained
using the ELMM. In fact, for the Concrete endmember, the GLMM
abundance map shows stronger components in the stadium stands

GLMM

the variability following the GLMM where correlation was imposed
over Ψn using a 3-D Gaussian filter. White Gaussian noise was
added to both datasets resulting in a SNR of 30dB.
To find the optimal parameters for the selected algorithms
we performed a grid search for each dataset. The parameter
ranges were chosen based on the ranges tested and discussed by
the authors in the original publication of each algorithm. For
the PLMM we used γ = 1, since the authors fixed this parameter in all simulations, and searched for α and β in the range
[0.35, 0.7, 1.4, 25] and [10−9 , 10−5 , 10−4 , 10−3 ], respectively.
For both the ELMM and the GLMM, the parameters were selected
among the following values: λS , λM ∈ [0.01, 0.1, 1, 5, 10, 15],
λA ∈ [0.001, 0.01, 0.05], and λψ , λΨ ∈ [10−6 , 10−3 , 10−1 ].
The results are presented in Table 1. In terms of RMSE for
the abundance vectors, RMSEA , the proposed strategy clearly outperformed the competing algorithms for both datasets. This behavior can be verified for almost all metrics considered. The only
exception is the SAMM for DC0 where PLMM and ELMM presented smaller spectral angles. Regarding the increase of computational complexity introduced by the GLMM when compared with
the ELMM, the simulations point out that the GLMM approach demanded 3.62 × TimeELMM for DC0 and 1.92 × TimeELMM for DC1,
where TimeELMM is the CPU time elapsed during the ELMM unmixing process. The results show that the extra flexibility of the GLMM
can be beneficial for the HU problem at the expense of a reasonable
increase in the computational complexity.

ELMM

PLMM

Fig. 1. Synthetic data cubes DC0, left, and DC1, right.

Fig. 2. Abundance maps of the Houston dataset for all tested algotithms where the abundance values are represented by colors ranging
from blue (αk = 0) to red (αk = 1).

when compared with the other methods considering spectral variability. Although the results presented in Table 1 indicate better fitting for the GLMM method in both RMSE and SAM, these results
should be taken with the proper care, since the connection of reconstruction error and abundance estimation is not straightforward. The
GLMM demanded a computational time of 0.92 × TimeELMM .
6. CONCLUSIONS
This paper proposed a new Generalized Linear Mixing Model
(GLMM) that accounts for endmember spectral variability. The new
model generalizes the Extended Linear Mixing Model (ELMM) to
allow for the consideration of band dependent scaling factors for the
endmember signatures. This way the GLMM model can represent
a larger variety of realistic spectral variations of the endmembers,
generalizing the representation capability of the ELMM. To solve
the resulting optimization problem, we extended the variable splitting methodology used in [10] by including new tensor variables.
Simulation results with both synthetic and real data suggest that the
extra flexibility introduced by the GLMM can be beneficial for the
unmixing process, resulting in improvements in both the abundance
estimation and the reconstruction error.

7. REFERENCES
[1] J. M. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders,
N. Nasrabadi, and J. Chanussot, “Hyperspectral remote sensing data analysis and future challenges,” IEEE Geoscience and
Remote Sensing Magazine, vol. 1, no. 2, pp. 6–36, 2013.
[2] N. Keshava and J. F. Mustard, “Spectral unmixing,” IEEE
Signal Processing Magazine, vol. 19, no. 1, pp. 44–57, 2002.
[3] 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 Journal of Selected Topics in Applied Earth
Observations and Remote Sensing, vol. 5, no. 2, pp. 354– 379,
2012.
[4] N. Dobigeon, J.-Y. Tourneret, C. Richard, J. C. M. Bermudez,
S. McLaughlin, and A. O. Hero, “Nonlinear unmixing of hyperspectral images: Models and algorithms,” IEEE Signal Processing Magazine, vol. 31, no. 1, pp. 82–94, Jan 2014.
[5] T. Imbiriba, J. C. M. Bermudez, C. Richard, and J.-Y.
Tourneret, “Nonparametric detection of nonlinearly mixed pixels and endmember estimation in hyperspectral images,” IEEE
Transactions on Image Processing, vol. 25, no. 3, pp. 1136–
1151, March 2016.
[6] T. Imbiriba, J. C. M. Bermudez, and C. Richard, “Band selection for nonlinear unmixing of hyperspectral images as a maximal clique problem,” IEEE Transactions on Image Processing,
vol. 26, no. 5, pp. 2179–2191, May 2017.
[7] A. Zare and K. C. Ho, “Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral
unmixing,” Signal Processing Magazine, IEEE, vol. 31, pp.
95, January 2014.
[8] Lucas Drumetz, Jocelyn Chanussot, and Christian Jutten,
“Variability of the endmembers in spectral unmixing: recent
advances,” in 8th IEEE Workshop on Hyperspectral Image and
Signal Processing: Evolution in Remote Sensing, Los Angeles,
USA, 2016.
[9] P.-A. Thouvenin, N. Dobigeon, and J.-Y. Tourneret, “Hyperspectral unmixing with spectral variability using a perturbed
linear mixing model,” IEEE Trans. Signal Processing, vol. 64,
no. 2, pp. 525–538, Feb. 2016.
[10] L. Drumetz, M.-A. Veganzones, S. Henrot, R. Phlypo,
J. Chanussot, and C. Jutten, “Blind hyperspectral unmixing
using an extended linear mixing model to address spectral variability,” IEEE Transactions on Image Processing, vol. 25, no.
8, pp. 3890–3905, 2016.
[11] ZG Gao and LQ Zhang, “Multi-seasonal spectral characteristics analysis of coastal salt marsh vegetation in shanghai,
china,” Estuarine, Coastal and Shelf Science, vol. 69, no. 1,
pp. 217–224, 2006.
[12] Heike Schmidt and Arnon Karnieli, “Remote sensing of the
seasonal variability of vegetation in a semi-arid environment,”
Journal of arid environments, vol. 45, no. 1, pp. 43–59, 2000.
[13] Petr Lukeš, Pauline Stenberg, Miina Rautiainen, Matti Mõttus,
and Kalle M Vanhatalo, “Optical properties of leaves and needles for boreal tree species in europe,” Remote Sensing Letters,
vol. 4, no. 7, pp. 667–676, 2013.

[14] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Total variation spatial regularization for sparse hyperspectral unmixing,”
IEEE Transactions on Geoscience and Remote Sensing, vol.
50, no. 11, pp. 4484–4502, 2012.
[15] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating
direction method of multipliers,” Found. Trends Mach. Learn.,
vol. 3, no. 1, pp. 1–122, Jan. 2011.
[16] Per Christian Hansen, James G Nagy, and Dianne P O’leary,
Deblurring images: matrices, spectra, and filtering, SIAM,
2006.
[17] R. N. Clark, G. A. Swayze, K. E. Livo, R. F. Kokaly, S. J.
Sutley, J. B. Dalton, R. R. McDougal, and C. A. Gent, “Imaging spectroscopy: Earth and planetary remote sensing with the
usgs tetracorder and expert systems,” Journal of Geophysical
Research: Planets, vol. 108, no. E12, 2003.
[18] J. M. P. Nascimento and J. M. Bioucas-Dias, “Vertex Component Analysis: A fast algorithm to unmix hyperspectral data,”
IEEE Transactions on Geoscience and Remote Sensing, vol.
43, no. 4, pp. 898–910, April 2005.

