Unmixing urban hyperspectral imagery with a Gaussian mixture model on endmember variability
Yuan Zhou Erin B. Wetherley Paul D. Gader
cs.CV






Y. Zhou et al.: Unmixing urban hyperspectral
imagery with a Gaussian mixture model on endmember variability

0000–0000/00$00.00  2018 IEEE
Unmixing urban hyperspectral imagery with a Gaussian mixture model
on endmember variability
    Yuan Zhou, Student Member, IEEE, Erin B. Wetherley,
and Paul D. Gader, Fellow, IEEEY. Zhou and P. D. Gader are with the Department of Computer and Information
Science and Engineering, University of Florida, Gainesville, FL, USA.
E-mail: mailto:yuan,pgader@cise.ufl.eduyuan,pgader@cise.ufl.edu.E. B. Wetherley is with the Department of Geography, University of
California Santa Barbara, Santa Barbara, CA, USA. E-mail: mailto:wetherley@umail.ucsb.eduwetherley@umail.ucsb.edu.
    Received: date / Revised version: date
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


Spectral unmixing given a library of endmember spectra can be achieved
by multiple endmember spectral mixture analysis (MESMA), which tries
to find the optimal combination of endmember spectra for each pixel
by iteratively examining each endmember combination. However, as library
size grows, computational complexity increases which often necessitates
a laborious and heuristic library reduction method. In this paper,
we model a pixel as a linear combination of endmembers sampled from
probability distributions of Gaussian mixture models (GMM). The parameters
of the GMM distributions are estimated using spectral libraries. Abundances
are estimated based on the distribution parameters. The advantage
of this algorithm is that the model size grows very slowly as a function
of the library size.

To validate this method, we used data collected by the AVIRIS sensor
over the Santa Barbara region: two 16 m spatial resolution and two
4 m spatial resolution images. 64 validated regions of interest (ROI)
(180 m by 180 m) were used to assess estimate accuracy. Ground truth
was obtained using 1 m images leading to the following 6 classes:
turfgrass, non-photosynthetic vegetation (NPV), paved, roof, soil,
and tree. Spectral libraries were built by manually identifying and
extracting pure spectra from both resolution images, resulting in
3,287 spectra at 16 m and 15,426 spectra at 4 m. We then unmixed ROIs
of each resolution using the following unmixing algorithms: the set-based
algorithms MESMA and AAM, and the distribution-based algorithms GMM,
NCM, and BCM. The original libraries were used for the distribution-based
algorithms whereas set-based methods required a sophisticated reduction
method, resulting in reduced libraries of 61 spectra at 16 m and 95
spectra at 4 m. The results show that GMM performs best among the
distribution-based methods, producing comparable accuracy to MESMA,
and may be more robust across datasets.



spectral unmixing, endmember variability, Gaussian mixture model,
MESMA, hyperspectral image analysis





§ INTRODUCTION




Hyperspectral images have important applications
in astronomy, agriculture, geoscience, surveillance (such as object
identification), material identification, and detecting processes
<cit.>. Because limited photons enter the sensor
when collecting narrow bandwidth channels from a high altitude, the
spatial resolution of hyperspectral image is usually very coarse,
i.e. a pixel may correspond to a region with a diameter of several
meters. Hence, multiple materials may exist in this region and contribute
to the measured pixel spectrum, also known as a mixed pixel <cit.>.
One important problem in hyperspectral imagery is to decompose mixed
pixels to identify the constituting materials (endmember) and
their proportions (abundance) that form the pixel spectrum.

The most common model that relates endmembers and abundances to a
pixel is the linear mixing model (LMM), which assumes that
the reflectance measured within each pixel is a unique linear combination
of the reflectances of each sub-pixel endmember, weighted by its abundance,
plus some noise <cit.>. The intuition behind this
model is that the fractional area of a material determines its representation
in the measured signal. However, when unmixing a hyperspectral image
with LMM, we usually encounter an additional problem that spectral
reflectance for identical materials are often different. For example,
asphalt spectra can vary significantly based on age, shadowing, and
composite materials <cit.>. This is sometimes
called endmember variability <cit.>.



Several factors can contribute to endmember variability, including
both extrinsic factors and intrinsic factors. The most significant
extrinsic factor is illumination. When solar incidence and emergence
angles are different for a surface, the observed signal will be different
<cit.>. Material angle matters as well, for example
roofs can be present at a variety of angles relative to incoming solar
radiation, producing different spectral signatures for one material.
Atmospheric condition can be another extrinsic factor affecting reflectance,
however this is usually corrected during image processing. Measurement
scale represents an important intrinsic factor. Objects or materials
that may be considered “pure”
may in reality be composed of materials at smaller scales with varying
reflectances <cit.>. For example, a tree canopy
can be considered a single, pure endmember, however this ignores the
spectral variety of tree leaves, bark, branches, and substrate that
composes a single tree pixel <cit.>. Similarly,
soils are composed of particles with different shapes, sizes, and
chemical composition <cit.>. The larger scale
we use to define an endmember, the larger intrinsic variability we
may expect from its spectra. For example, trees and turfgrass can
be defined as individual endmembers, however if we we wish to define
a class of green vegetation comprised of both turfgrass and tree,
its variability will not be less than the component endmember.

Considering endmember variability, we can generalize the LMM to the
following equation:


    𝐲_n=∑_j=1^M𝐦_njα_nj+𝐧_n, n=1,…,N

where 𝐲_n∈ℝ^B is the spectrum of the nth
pixel in the image, B is the number of bands, N is the number
of pixels, M is the number of endmembers. 𝐦_nj∈ℝ^B
is the jth endmember for the nth pixel. α_nj∈ℝ
is the abundance that usually satisfies the positivity and sum-to-one
constraints, i.e. α_nj≥0, ∑_jα_nj=1. Finally,
we have some additive noise 𝐧_n.

When it comes to unmixing in terms of (<ref>), we are referring
to retrieving {𝐦_nj,α_nj} from
{𝐲_n}, or {α_nj}
from {𝐲_n} and a library of endmember
spectra. The former is sometimes called unsupervised unmixing,
and because it is undetermined this can be a difficult problem. Studies
that have worked to solve unsupervised unmixing usually require several
assumptions, such as spatial smoothness of the abundances and the
existence of contiguous pure pixels <cit.>.
The latter is called supervised unmixing and depends on a library
of known endmember spectra. If the library is small enough to easily
enumerate all possible spectral combinations, the task can be trivial.
However, applying this scheme on larger libraries becomes computationally
inefficient. This is the problem we are addressing in this study.

Previous studies that have worked to solve this problem have used
methods that can be categorized as set-based or distribution-based
<cit.>. Set-based methods treat the endmember library
as an unordered set and try to pick the best combination of endmembers
to model each pixel. A widely used set-based method is multiple endmember
spectral mixture analysis (MESMA) <cit.>. The
general idea of MESMA is to test every endmember combination and select
the one with the smallest error within set thresholds that limit pixel
complexity. There are many variations to MESMA. In multiple-endmember
linear spectral unmixing model (MELSUM), the solution for abundances
is obtained from directly solving the linear equations and discarding
the negative values <cit.>. In automatic Monte
Carlo unmixing (AutoMCU), pixels are unmixed using multiple sets of
random combinations, with the mean fractional values assigned as abundances
<cit.>. In alternate angle
minimization (AAM), projection is iteratively used to find the spectrum
index of one endmember given the other endmembers fixed. Besides MESMA
variants, there is sparse unmixing that used the full spectral library
with a sparsity constraint on the abundances forcing them having only
a few nonzero elements <cit.>.

Contrary to set-based methods, distribution-based methods assume that
the endmembers for each pixel are sampled from probability distributions,
hence the linear combinations of these endmembers (pixels) also follow
some distribution. It works by modeling the spectral library as statistical
distributions, extracting parameters to describe these distributions,
and unmixing the pixels based on the distribution parameters. The
most widely used distribution is Gaussian, and its application for
spectral unmixing is known as the normal compositional model (NCM)
<cit.>.
The popularity of NCM comes from the fact that a linear combination
of Gaussian random variables is also a Gaussian random variable whose
mean and covariance matrix are linear combinations from the endmember
means and covariance matrices. Hence, the resulting probability density
function of the pixels has a simple analytical form. Fitting the actual
pixel values to the pixel distribution, the abundances can be solved
by several techniques, such as expectation maximization <cit.>,
sampling methods <cit.>,
and particle swarm optimization <cit.>.

Following this philosophy, some have worked to extend the idea to
distributions beyond Gaussian. In <cit.>, the authors
propose Beta distributions to model the spectral library. The benefit
is that Beta distributions have a domain in the range 0 –
1, so are more suitable for the reflectance range, and the actual
library may have a skewed mode in the distribution. In <cit.>,
the idea is further extended to use Gaussian mixture models (GMM)
for distributions. The rationale comes from the observation that library
endmembers may have multiple modes, whose shape cannot be represented
by a simple Gaussian or Beta distribution. Since GMM is more flexible,
it can approximate any distribution found in the library.



 §.§ Our contribution


Many unmixing studies are not well evaluated in presence of ground
truth. Commonly used hyperspectral datasets include Pavia University,
Indian Pines, Cuprite, Mississippi Gulfport, etc., which are not validated
with ground truth endmembers and abundances. Hence, the primary method
for evaluating their results include:


  * Compare the estimated endmembers with spectra in the USGS spectra
library (e.g. in Cuprite dataset) <cit.>.

  * Compare the estimated abundances with assumed segmentation maps of
pure materials (e.g. in Indian Pines, Pavia University, Gulfport datasets)
<cit.>.

  * Calculate the reconstruction error of estimated endmembers and abundances
and assume that a lower reconstruction error implies a better result
<cit.>.

Each of these methods can be problematic. First, different conditions
(sensor, atmosphere, light source) during data collection will affect
measured reflectances, making library comparison less ideal. Second,
high spatial resolution hyperspectral images are primarily composed
of pure pixels, and segmentation like abundance maps do not necessarily
indicate good unmixing capability for mixed pixels. Third, reconstruction
error is more related to model complexity than unmixing accuracy since
small reconstruction error could be achieved by overfitting <cit.>.

Moreover, these datasets are not comprehensive with respect to spatial
scales, scene diversity and generalization. For example, the Pavia
University and Gulfport datasets have about 1 m spatial resolution
in which most are pure pixels. Also, they are focused on only a few
urban sites, which contain mostly man-made materials with segmentation
like abundance distribution. Developing unmixing algorithms on them
will have a bias on forcing smooth and sparse abundance maps. Hence,
it is unknown if the algorithms validated on these datasets can be
applied to datasets with generalized scenarios.

In this work, we introduce a supervised unmixing algorithm based on
modeling endmember variability by GMM distributions, and compare several
set-based and distribution-based algorithms with a highly validated,
comprehensive dataset of 128 images with different spatial scales.
The algorithm was first introduced in <cit.>
and we modified it for this application. The dataset was developed
in <cit.> but only used for evaluating MESMA.
It contains two types of images, one with about 16 m pixel size, the
other with about 4 m pixel size. It covers a wide range of landcover,
including various kinds of road, roof, vegetation, and soil. Validation
abundances were obtained by classifying high resolution images corresponding
to the hyperspectral images. Unlike MESMA, which requires a small
and well-curated spectral library, the GMM algorithm uses the original
source library without modification to unmix fractions using inferred
parameters from the library.



§ DATASET


We used two low-resolution images (16 m) and two high-resolution images
(4 m) in this study. The low-resolution images were collected by the
Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) <cit.>
over Santa Barbara, CA, on August 29, 2014. The spatial resolutions
are 15.6 m/pixel and 15.8 m/pixel. The spectral range measures wavelengths
from 380 – 2500 nm with 224 bands of approximately 10
nm bandwidth. High-resolution images were collected by AVIRIS-Next
Generation with 3.9 m/pixel and 3.6 m/pixel spatial resolutions. The
spectral resolution is also higher, recording 432 bands of about 5
- 6 nm bandwidth across a similar spectral range as the 16 m dataset.
We spectrally resampled the AVIRIS-Next Generation imagery to 224
bands to produce an image with identical spectral parameters to the
16 m AVIRIS image. We also removed certain bands from analysis due
to atmospheric interference, reducing the number of bands to 164.
Initial image processing was conducted by the Jet Propulsion Laboratory,
with additional processing in the lab to reduce the effects of elevation
change on pixel location. 

The study area includes the cities of Santa Barbara and Goleta as
well as the land between them, near the California coast. Urban composition
is typical of the southwestern United States, including man-made materials
such as asphalt, concrete, metal, gravel, and brick, as well as vegetation
in the forms of turfgrass, various tree species, and large areas of
undeveloped land covered in senesced vegetation <cit.>.



 §.§ Validation Polygons


We produced 64 polygons that represented the variety of landcover
within the study area. Each polygon was 180 m by 180 m in size, or
11-12 pixels wide in the 16 m images (46 or 50 pixels wide in the
4 m images). Validation polygons were randomly distributed across
the area with a minimum distance 400 m. If a polygon contained large
areas of open water or an undetermined material, it was discarded
and a new polygon randomly generated. Cover was determined within
each polygon using a 1 m NAIP high-resolution image. We used a combination
of image segmentation, using ECognition, and manual adjustments to
classify the cover within each polygon as turf, tree, paved, roof,
soil, or non-photosynthetic vegetation (NPV). Cover was further confirmed
by visually inspecting August 2014 Google Earth imagery. Fig <ref>
displays all polygons as they appear in the 16 m images. 



Fig. <ref> shows a scatter plot of the 64 ground
truth abundances when the 6 endmember classes are merged to 3 categories
of vegetation, impervious, and non-vegetated pervious. Most polygons
are dominated by a mixture of impervious and vegetation materials.
To improve the representation of less common mixtures in the scene,
we added 5 polygons with high proportions of soil. 






 §.§ Library Building


We produced 240 polygons across the 4 m scene to extract pure spectra
and build the full spectral libraries. The polygons were intended
to capture class material variability as much as possible, and so
included multiple roof types, asphalt, concrete, trees, turfgrass,
soil, and NPV, as well as less common materials like rubber, solar
panels, tennis courts, and plastic tarps. These materials were then
grouped into one of our 6 endmember classes: turfgrass, NPV, paved,
roof, soil, and tree.

The same polygons were used to extract spectra from the 16 m imagery,
with necessary modifications as described in <cit.>.
Together, we produced a library of 16 m spectra and a library of 4
m spectra. After removing duplicate spectra, the final 16 m library
was comprised of 3,287 spectra and the 4 m library contained 15,426
spectra. 

These full spectral libraries were used to train the parameters of
distribution-based algorithms. However, they were too large to be
used by MESMA, and required reduction. We performed reduction in two
steps. First, iterative endmember selection (IES) <cit.>
was used to automatically select a subset of spectra that represented
the larger library. This is achieved iteratively, by gradually selecting
the most representative spectra and evaluating their representativeness
using a kappa coefficient. IES reduced the 16 m and 4 m library sizes
to 226 and 187, respectively. Libraries were further reduced using
iterative classification reduction (ICR), which uses MESMA as a classifier
to quickly identify and remove spectra that tend to map materials
incorrectly <cit.>. This reduced the libraries
to a final size of 61 for 16 m images and 95 for 4 m images. The spectra
for each endmember class for all the cases are plotted in Fig <ref>,
and their numbers are shown in Table <ref>.








§ METHOD




 §.§ The Gaussian Mixture Model for Unmixing


Here we briefly introduce the GMM based unmixing <cit.>,
which is a generative model that models a distribution on the input
space <cit.>. Suppose we have M endmember classes,
each has numerous spectra in the library. A pixel can be assumed to
be generated by randomly picking one spectrum for each endmember,
and linearly mixing them based on some abundances. In this way, if
we use a probability density function to represent the spectral distribution,
the actual endmembers can be assumed to be sampled from this distribution.
Suppose the jth endmember for the nth pixel is sampled from
a distribution modeled by GMM

    p(𝐦_nj|Θ)=∑_k=1^K_jπ_jk𝒩(𝐦_nj|μ_jk,Σ_jk),

where Θ={π_jk,μ_jk,Σ_jk: j=1,…,M,k=1,…,K_j}
are the GMM parameters. Allowing GMM to represent the library, we
can get multiple Gaussian components for each endmember. Take the
dataset in Section <ref> as an example, which can be
viewed from two perspectives. Fig. <ref> shows the
pixels from all the validation ROIs, library endmembers, and corresponding
Gaussian components when projected to 2 dimensions. The method for
estimating GMM parameters will be discussed later, however we can
see that the ellipses formed by these parameters surround validation
pixels at multiple positions on the edge of the pixel cloud. The pixels
can be viewed as picking points within ellipses and combining these
points linearly. Fig. <ref> shows the Gaussian components
from the wavelength-reflectance perspective, where the centers of
Gaussian components and their variation patterns are shown as curves.
Compared to MESMA, which evaluates every spectrum in the library,
GMM tries to combine every center of Gaussian components, but allows
the center to move according to its corresponding variation pattern.





Following the distribution assumption, if {𝐦_nj: j=1,…,M}
are independent and the noise is also sampled from a Gaussian p(𝐧_n)=𝒩(𝐧_n|0,𝐃),
then 𝐲_n=∑_j𝐦_njα_nj+𝐧_n
implies that the pixel follows a distribution

    p(𝐲_n|α_n,Θ,𝐃)=∑_𝐤∈𝒦π_𝐤𝒩(𝐲_n|μ_n𝐤,Σ_n𝐤),

where 𝒦={ 1,…,K_1}×{ 1,…,K_2}×⋯×{ 1,…,K_M}
is the Cartesian product of the M index sets, 𝐤=(k_1,…,k_M)∈𝒦,
π_𝐤∈ℝ, μ_n𝐤∈ℝ^B,
Σ_n𝐤∈ℝ^B× B are
defined by 


    π_𝐤=∏_j=1^Mπ_jk_j, μ_n𝐤=∑_j=1^Mα_njμ_jk_j, Σ_n𝐤=∑_j=1^Mα_nj^2Σ_jk_j+𝐃.


If we assume each pixel is independently sampled, the probability
density function of all the pixels is the product as


    p(𝐘|𝐀,Θ,𝐃)=∏_n=1^Np(𝐲_n|α_n,Θ,𝐃),

where 𝐀=[α_1,…,α_N]^T∈ℝ^N× M.
Given 𝐘,Θ,𝐃, the abundances
𝐀 can be estimated by maximum likelihood estimation
(MLE). Specifically, we want to maximize p(𝐘|𝐀,Θ,𝐃),
or minimize -log p(𝐘|𝐀,Θ,𝐃),
which becomes the following optimization problem by combining the
above equations


    ℰ(𝐀)=-∑_n=1^Nlog∑_𝐤∈𝒦π_𝐤𝒩(𝐲_n|μ_n𝐤,Σ_n𝐤),



    s.t. α_nj≥0, ∑_j=1^Mα_nj=1, ∀ n.


The objective function can be minimized by a generalized expectation
maximization (EM) algorithm, which alternates between an E step and
an M step <cit.>. The E step calculates the posterior
probability of the latent variable given the observed data and old
parameters. The M step increases the expected value of the complete
data log-likelihood. In our case, the E step calculates


    γ_n𝐤=π_𝐤𝒩(𝐲_n|μ_n𝐤,Σ_n𝐤)/∑_𝐤∈𝒦π_𝐤𝒩(𝐲_n|μ_n𝐤,Σ_n𝐤).

The M step tries to minimize

    ℰ_M=-∑_n=1^N∑_𝐤∈𝒦γ_n𝐤{logπ_𝐤+log𝒩(𝐲_n|μ_n𝐤,Σ_n𝐤)} .

It does not have a closed form solution for 𝐀. But we
can use gradient descent to minimize ℰ_M, where the
derivative can be calculated by


    ∂ℰ_M/∂𝐀=-∑_𝐤∈𝒦Λ_𝐤𝐑_𝐤^T-2𝐀∘∑_𝐤∈𝒦Ψ_𝐤𝐒_𝐤^T,

where 𝐑_𝐤∈ℝ^M× B, 𝐒_𝐤∈ℝ^M× B^2
are defined by 

    𝐑_𝐤=[μ_1k_1,μ_2k_2,…,μ_Mk_M]^T,


    𝐒_𝐤=[vec(Σ_1k_1),vec(Σ_2k_2),…,vec(Σ_Mk_M)]^T,

and Λ_𝐤∈ℝ^N× B,
Ψ_𝐤∈ℝ^N× B^2 denote

    Λ_𝐤=[λ_1𝐤,λ_2𝐤,…,λ_N𝐤]^T,


    Ψ_𝐤=[vec(Ψ_1𝐤),vec(Ψ_2𝐤),…,vec(Ψ_N𝐤)]^T,

where λ_n𝐤∈ℝ^B×1
and Ψ_n𝐤∈ℝ^B× B are 


    λ_n𝐤=γ_n𝐤Σ_n𝐤^-1(𝐲_n-μ_n𝐤),


    Ψ_n𝐤=1/2γ_n𝐤Σ_n𝐤^-T(𝐲_n-μ_n𝐤)(𝐲_n-μ_n𝐤)^TΣ_n𝐤^-T-1/2γ_n𝐤Σ_n𝐤^-T.

Given an initial 𝐀, we can update γ_n𝐤
and 𝐀 alternately until convergence, which leads to a
local minimum of the objective function. This algorithm requires several
clarifications and we will explain them in the following subsections.



 §.§ Determining the GMM Parameters


Suppose we have a library of endmember spectra {𝐘_j∈ℝ^N_j× B: j=1,…,M},
with which we can estimate the GMM parameters Θ.
The difficulty comes from estimating the number of components K_j
for each endmember, as once we know K_j, {π_jk,μ_jk,Σ_jk}
can be estimated by the standard EM algorithm. Estimating this K_j
is sometimes called model selection and has several approaches
<cit.>. We will use cross-validation-based information
criterion (CVIC) <cit.> as our metric to select K_j.

Given a candidate K_j, we can evaluate CVIC in the following
way. Let 𝐘_j be the spectra for the jth endmember
in the library, we can divide them into V=5 subsets with equal
size. For each subset 𝐘_j^v, the remaining spectra
are input to a MLE with K_j Gaussian components, and the trained
parameters are used to evaluate the log-likelihood of 𝐘_j^v.
Then the sum of all these log-likelihood values is calculated as ℒ_K_j,
which is our CVIC. Finally, the optimal K_j is the one that maximizes
ℒ_K_j out of all the candidates. To avoid many components,
we tried K_j=1,2,3,4.

This approach can serve as an ideal model selection. However, the
number of combinations |𝒦|=∏_jK_j
can still be very large, especially in real datasets where the libraries
contain many spectra. Hence, we use a threshold T_CVIC to further
reduce K_j manually. Let ℒ_j^' be the maximum
CVIC for the jth endmember; we pick the smallest K_j such
that |ℒ_K_j-ℒ_j^'|≤ T_CVICℒ_j^'.
Hence when T_CVIC=0, we have the ideal CVIC-based model selection.
As T_CVIC increases, we can have a reduced number of components.



 §.§ Projection


Analyzing the computation in Section <ref>,
we see that the time complexity is O(|𝒦|NB^3)
<cit.>. In addition to the number of combinations,
the number of bands is also crucial to execution time. We can reduce
the time cost by reducing the dimensionality of the data.

We use PCA to reduce the dimensionality, which gives a center 𝐜∈ℝ^B
and a projection matrix 𝐄∈ℝ^B× d such
that all the spectra are processed as 𝐄^T(𝐲-𝐜).
Note that (<ref>) still holds if both the pixel spectra
and endmember spectra are projected in this way. Hence the estimated
abundances for the projected spectra are the original abundances for
the data. Also note that if an endmember follows a GMM distribution,
the projected endmember also follows a GMM distribution. So we can
directly estimate the GMM parameters from the projected library spectra.

As for the data input to find this projection, there are two possibilities.
One is to use all the pixel data. This works if the image is big enough,
such that all the endmembers have sufficient presence. However, if
the image contains fewer pixels (e.g. in the 16 m dataset) with limited
endmembers, the directions determined by PCA will present the variation
within the image, which may not distinguish different endmembers.
The other method is to use the spectral library, i.e. combine all
the spectra in the library and perform PCA on them. We adopted this
method for our dataset. Specifically, we selected an equal number
of spectra for each endmember class in the library and concatenated
them. This ensures that the relative sizes of endmember classes do
not affect the direction, and also ensures that the mean lies in the
center.



 §.§ The Algorithm


The implementation of the algorithm is described in Algorithm <ref>.
In step 1, the spectra in the library are concatenated to form an
input to PCA. We project the data to 10 dimensions in step 2. Step
3 is elaborated in Section <ref>. Step
4 involves initialization of 𝐀, which utilizes the information
of multiple means from the Gaussian components. To be specific, we
set α_n𝐤←(𝐑_𝐤𝐑_𝐤^T+ϵ𝐈_M)^-1𝐑_𝐤𝐲_n,
project α_n𝐤 onto the simplex space,
and initialize α_n←α_n𝐤̂
where 𝐤̂=min_𝐤‖𝐲_n-𝐑_𝐤^Tα_n𝐤‖^2.
Step 5 is the main body, in which the M step is the most complicated
part. Because of the constraints on α_n, we
use projected gradient descent here. The projection function can be
found in <cit.>. The step size
τ can be set adaptively by using a small initial value and gradually
increasing it by a multiplier of 10 as long as the objective function
keeps decreasing. Since the covariance matrix from endmember variability
is usually much larger than the noise covariance, the latter can be
negligible in the computation and we use 𝐃=0.001^2𝐈_B.






§ RESULTS




 §.§ Setup


We ran GMM on the 16 m and 4 m images after training it using the
same resolution spectral libraries. Since GMM takes spectral libraries
to estimate the Gaussian components, we used the same components on
all 64 images. For reproducibility, we ran model selection of GMM
15 times, selected the most frequent combination, and applied it to
the dataset.

For comparison, we ran 2 set-based methods and 4 distribution-based
methods:


  * MESMA <cit.>. It was implemented in IDL by the
original authors and provided as an extension Viper Tool to ENVI.
We used the same parameters as in <cit.>, i.e.
maximum RMSE 2.5%, threshold RMSE 0.7%, abundances constrained between
0 and 1, maximum shade threshold 80%, and a maximum of three endmembers
plus shade for each pixel. Also, it will not allow multiple spectra
from one endmember class in the mixture. The obtained fractions were
normalized to give the final abundances.

  * Alternate angle minimization (AAM) <cit.>.
Its code was implemented in Matlab and downloaded from Rob Heylen’s
website. It tries every subset of endmembers, iteratively updates
the spectrum index of each endmember such that the reconstruction
error is minimized given the rest selected spectra, and hence finds
the best combination and abundances. Since it uses projection to find
the combination, theoretically it should work faster. It is different
from MESMA in several ways. First, it may not find the global minimum
because of its alternate optimization. Second, it may find a pixel
mixed by many endmembers instead of maximum three. Finally, it does
not include a shade endmember to adjust for brightness differences
between endmember and measured spectra. 

  * Gaussian mixture model (GMM). We used T_CVIC=0.05 for the 16
m dataset, which produced 216 combinations from the library. Because
the 4 m dataset was about 20 times larger, and the number of spectra
in the library was 3 times larger, we used a larger T_CVIC=0.2
(18 combinations) such that the whole process could still run in a
few hours.

  * GMM-1. We set the number of combinations to be 1, i.e. one component
for each endmember, which makes it to be NCM theoretically. However,
it has the same implementation of GMM hence reflects the difference
introduced by bringing multiple components.

  * Normal compositional model with sampling optimization (NCM Sampling).
There are many variations of NCM, with different optimization approaches
<cit.>.
We chose the sampling strategy in <cit.> which does
not assume statistical independence between different bands.

  * Beta compositional model (BCM) <cit.>. It is available
from Alina Zare's website. Assuming the independence of bands, it
uses Beta distribution to model each band and unmixes the pixels.

Excluding MESMA, which was implemented in IDL, all methods were implemented
in Matlab. MESMA was run on a PC with Intel Core i7-2760QM CPU and
8 GB memory. The other methods were run on a PC with Intel Core i7-3820
CPU and 64 GB memory. For distribution-based methods, the original
libraries were input to train the parameters for unmixing while set-based
methods used the reduced libraries.

We used two metrics to measure the differences between the estimated
and ground truth fractions: mean absolute difference (MAD) and correlation
coefficient (R). They were calculated for each endmember class
based on the 64 pairs of values. To visualize the values, we used
a variation of the Bland-Altman plot where the x-axis is the ground
truth value and the y-axis is the estimated minus the ground truth
value <cit.>. When comparing different algorithms
for unmixing quality, we will mainly resort to MAD as correlation
coefficient itself is not sufficient (slope and intercept are needed
to accompany R).



 §.§ Accuracy and Efficiency


16 m Case. Table <ref> shows the MAD
and correlation coefficient for the 16 m images. Original errors for
6 classes implies that GMM and AAM have the best accuracy, followed
by MESMA. The difference comes from the paved, roof and tree classes,
where GMM outperformed MESMA. In general, MESMA, AAM, GMM and GMM-1
had similar accuracy. Comparing all the distribution-based methods
with same input, GMM has the best performance overall, with fewest
errors for NPV, paved and roof.

Merged errors show that GMM and AAM retain their higher accuracy,
with MESMA falling further behind due to poor impervious accuracy.
Since merged errors are differences between summed up quantities of
similar materials, such as paved and roof, if the errors are enlarged
compared to individual errors, it means that both of the similar materials
are overestimated or underestimated, i.e. the algorithm tends for
confuse even dissimilar materials, such as pervious and impervious.



Fig. <ref> compares the estimated total abundances
to validated abundances for each material. We can see that NCM Sampling
and BCM tend to ignore paved or roof when they have presence. The
difference statistics between the estimated and the ground truth are
plotted in the Bland-Altman plots in Fig. <ref>,
where the set-based methods and GMM appear to be better than the others.





4 m Case. The error statistics of 4 m data are shown in Table <ref>.
MESMA and GMM are the most accurate, with AAM is not as good as MESMA.
One possible reason for the high accuracy of MESMA is that MESMA inherently
takes shade into account while AAM only combines input library spectra.
Hence, the slightly better performance of GMM over AAM is more significant
since they both ignore shade.



Fig. <ref> and Fig. <ref>
show the scatter plots and Bland-Altman plots for this data. Similar
to the statistics in Table <ref>, GMM has an advantage
over the other distribution-based methods on paved and roof. Compared
to AAM, GMM presents stable results to some extremely bad outliers
for soil in AAM.





Efficiency. Since they were run on different machines with
different implementation (MESMA and NCM Sampling have multiple threads),
the time costs are not for comparison, but for reference. In general,
all the algorithms run in a few hours. The fastest algorithm is GMM-1,
which is our implemented GMM with only one combination. NCM Sampling
turns out to be the slowest algorithm. It is expected since sampling
algorithms are usually slower than deterministic algorithms. The times
costs on the 4 m dataset are usually more than 10 times slower than
those on the 16 m dataset. This is because the image size of the former
is 16-19 times larger than the latter and the library size is also
larger. The least gap comes from GMM because of a changed T_CVIC
leading to a significantly less number of combinations. In <cit.>
the authors show that with the same implementation AAM is much faster
than MESMA, but here the result is converse. One reason is that the
parameters of MESMA force it to pick at most 3 endmembers for a pixel
instead of all the combinations. Also, multi-threading and implementation
techniques impact the real world time costs significantly.



 §.§ Extend to Semi-realistic Images


We extended the experiments to semi-realistic images to check if the
algorithm implementation or library reduction was overfitted to this
particular dataset. The method was to test the algorithms on another
batch of synthetic images generated by the library spectra. Since
all the algorithms assumed that a pixel was a linear combination of
endmember spectra from the library, the creation of this synthetic
dataset would use this assumption.

We created this dataset following the literature that emphasizes realistic
simulation <cit.>. For each image
in the original dataset, we randomly sampled spectra from the full
library while the number of spectra for each endmember class is equal
to the number in the reduced library. Then we used AAM to unmix the
image with these sampled spectra. The obtained abundances were sorted
to keep the largest three while the other were set to 0, and rescaled
such that their summation was one. This is to conform with the assumption
of MESMA. The endmembers and abundances were combined according to
the LMM to generate pixels. In this way, we can have a dataset where
the endmembers are randomly picked from the full library, and the
spatial distribution of abundances looks similar to the original one.
Fig. <ref> shows all the 4 m synthetic images
generated in this way. Comparing it with Fig. <ref>,
we can see its similarity. But inherently, the synthetic images follow
exactly the LMM with at most 3 endmembers for each pixel and they
are randomly picked from the full library.



We ran all the algorithms on this simulated dataset. Table <ref>
shows the unmixing results on this dataset. We can see that GMM and
NCM Sampling turn out to be the best methods. The superior performance
of NCM contrasts sharply to its worst result in Table <ref>.
Since we evaluate the difference between total abundances for a material,
it is possible that the relatively large pixel abundance error is
mitigated by averaging them. This is more possible for sampling algorithms
because statistically they tend to sample values around the correct
ones. 



Compared to the negligible difference between GMM and set-based methods
in Section <ref>, the advantage of
GMM is more obvious. Since the endmembers were randomly sampled from
the big library, set-based methods were less capable to unmix the
pixels using a reduced library that was derived based on the another
dataset. It is possible that a different reduced library based on
this simulated dataset may lead to better results for set-based methods.
However, that means, set-based methods may not be as robust as GMM
across datasets.



§ DISCUSSION AND CONCLUSION


We have proposed an unmixing algorithm based on endmember variability
modeled by GMM distributions. We validated the algorithm on a dataset
consisting of 128 images across 2 scales, with ground truth abundances
obtained by inspecting high resolution images. The results show that
with large libraries, GMM achieves comparable accuracy to MESMA without
the need for guided manual library reduction. We will discuss several
issues regarding the dataset, algorithm and results in this Section.

The dataset was well developed with various scenes, but the ground
truth has intrinsic errors coming from UTM coordinates. It happens
when the universal coordinates are used in the hyperspectral images
and the other high-resolution images for region correspondence. Because
these airborne images are spatially calibrated from its unstable collection
process, the coordinates derived from the map information may not
accurately reflect the real coordinates. Therefore, the region for
calculating the total abundances may have a small shift compared to
the region in the hyperspectral image. This is more likely in the
16 m data, which may explain larger overall errors for all the algorithms,
and could be mitigated by registering the two images to find exact
correspondence in the future <cit.>.

In applying GMM on the dataset, there are several implementation details
that affected the results. First, the projection affected the results
of GMM. Hence, we used a carefully determined projection in Algorithm <ref>.
Second, the number of combinations impacted the performance. If some
training data were present, we could gradually decrease T_CVIC
until the error stopped improving during reasonable time. Otherwise,
we suggest gradually increasing T_CVIC from 0 (without running
the whole algorithm) until the number of combinations reduces to approximately
100, which is a reasonable time cost multiplier of GMM-1. Third, we
didn’t use a spatial prior in the original GMM paper
<cit.>. We found that the prior made it work better
on some images, while worse on some other images when applied to this
dataset. In total, it didn’t improve the results much.
This has two possible reasons: (i) the pixel size is big enough such
that smoothness and sparsity are not obvious on the abundance maps;
(ii) the dataset contains a variety of scenes in which some of them
violate this property.

We also have some remarks on the results. For the real dataset, GMM
is slightly better than MESMA for the 16 m data while slightly inferior
for the 4 m data, so they are close in accuracy. However, MESMA used
a library reduction method that relied on manual guidance and user
knowledge of the study area, while GMM used the original library without
refinement. Hence, GMM may be more applicable to datasets without
ground truth. Second, the results of MESMA were slightly different
from those reported in <cit.>, which used this
same imagery and validation dataset. In <cit.>,
the average R^2 for individual and merge cases are 0.642 and
0.867 for 16 m, 0.811 and 0.923 for 4 m. Comparing them to Table <ref>
and Table <ref>, they are very similar. Note that
we used Matlab to extract the polygon ROIs directly from the original
images while in <cit.> the images were resampled
to a uniform spatial resolution of 18 m and corrected for locational
errors using Delaunay triangulation. Additionally, that study used
178 spectral bands in contrast to 164 bands used here.

Future work could include developing a dataset with more images and
more scale variation. Also, validation abundances can be more accurately
obtained by registration of hyperspectral and high resolution images.
Further work can also be done to improve the efficiency of GMM using
covariance matrices with a simple form.



§ ACKNOWLEDGMENTS


The authors wish to thank Alina Zare and Sheng Zou for providing the
NCM sampling code. We wish to acknowledge funding from the NASA Earth
and Space Science Fellowship Program and the Belgian Science Policy
Office in the framework of the STEREO III Program— Project
UrbanEARS (SR/00/307). We also thank the Jet Propulsion Laboratory
for providing radiometrically calibrated, orthorectified AVIRIS imagery. 

11C__Users_Yuan_Zhou_Dropbox_YuanHyperspectral_document_GMM_SantaBarbara_IEEEbib1
10

vane1993airborne
G. Vane, R. O. Green, T. G. Chrien, H. T. Enmark, E. G. Hansen, and W. M.
  Porter,
“The airborne visible/infrared imaging spectrometer (AVIRIS),”
Remote Sensing of Environment, vol. 44, no. 2-3, pp. 127–143,
  1993.

small2001estimation
C. Small,
“Estimation of urban vegetation abundance by spectral mixture
  analysis,”
International journal of remote sensing, vol. 22, no. 7, pp.
  1305–1334, 2001.

settle1993linear
J. Settle and N. Drake,
“Linear mixing and the estimation of ground cover proportions,”
International Journal of Remote Sensing, vol. 14, no. 6, pp.
  1159–1177, 1993.

herold2005spectral
M. Herold and D. Roberts,
“Spectral characteristics of asphalt road aging and deterioration:
  implications for remote-sensing applications,”
Applied Optics, vol. 44, no. 20, pp. 4327–4334, 2005.

zare2014endmember
A. Zare and K. Ho,
“Endmember variability in hyperspectral analysis: Addressing
  spectral variability during spectral unmixing,”
IEEE Signal Processing Magazine, vol. 31, no. 1, pp. 95–104,
  2014.

somers2011endmember
B. Somers, G. P. Asner, L. Tits, and P. Coppin,
“Endmember variability in spectral mixture analysis: A review,”
Remote Sensing of Environment, vol. 115, no. 7, pp. 1603–1616,
  2011.

meister2000brdf
G. Meister, A. Rothkirch, H. Spitzer, and J. Bienlein,
“BRDF field studies for remote sensing of urban areas,”
Remote Sensing Reviews, vol. 19, no. 1-4, pp. 37–57, 2000.

wentz2012synthesizing
E. A. Wentz, D. A. Quattrochi, M. Netzband, and S. W. Myint,
“Synthesizing urban remote sensing through application, scale, data
  and case studies,”
Geocarto International, vol. 27, no. 5, pp. 425–442, 2012.

roberts2004spectral
D. A. Roberts, S. L. Ustin, S. Ogunjemiyo, J. Greenberg, S. Z. Dobrowski,
  J. Chen, and T. M. Hinckley,
“Spectral and structural measures of northwest forest vegetation at
  leaf to landscape scales,”
Ecosystems, vol. 7, no. 5, pp. 545–562, 2004.

halimi2015unsupervised
A. Halimi, N. Dobigeon, and J.-Y. Tourneret,
“Unsupervised unmixing of hyperspectral images accounting for
  endmember variability,”
IEEE Trans. on Image Processing, vol. 24, no. 12, pp.
  4904–4917, 2015.

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

zhou2018gmmJournal
Y. Zhou, A. Rangarajan, and P. D. Gader,
“A Gaussian mixture model representation of endmember variability
  in hyperspectral unmixing,”
IEEE Transactions on Image Processing, vol. To Appear, 2018.

roberts1998mapping
D. A. Roberts, M. Gardner, R. Church, S. Ustin, G. Scheer, and R. Green,
“Mapping chaparral in the Santa Monica Mountains using multiple
  endmember spectral mixture models,”
Remote Sensing of Environment, vol. 65, no. 3, pp. 267–279,
  1998.

combe2008analysis
J.-P. Combe, S. Le Mouelic, C. Sotin, A. Gendrin, J. Mustard, L. Le Deit,
  P. Launeau, J.-P. Bibring, B. Gondet, Y. Langevin, et al.,
“Analysis of OMEGA/Mars express data hyperspectral data using a
  multiple-endmember linear spectral unmixing model (MELSUM): Methodology and
  first results,”
Planetary and Space Science, vol. 56, no. 7, pp. 951–975,
  2008.

asner2000biogeophysical
G. P. Asner and D. B. Lobell,
“A biogeophysical approach for automated SWIR unmixing of soils
  and vegetation,”
Remote sensing of environment, vol. 74, no. 1, pp. 99–112,
  2000.

asner2002spectral
G. P. Asner and K. B. Heidebrecht,
“Spectral unmixing of vegetation, soil and dry carbon cover in arid
  regions: comparing multispectral and hyperspectral observations,”
International Journal of Remote Sensing, vol. 23, no. 19, pp.
  3939–3958, 2002.

castrodad2011learning
A. Castrodad, Z. Xing, J. B. Greer, E. Bosch, L. Carin, and G. Sapiro,
“Learning discriminative sparse representations for modeling, source
  separation, and mapping of hyperspectral imagery,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 49,
  no. 11, pp. 4263–4281, 2011.

eches2010bayesian
O. Eches, N. Dobigeon, C. Mailhes, and J.-Y. Tourneret,
“Bayesian estimation of linear mixtures using the normal
  compositional model: Application to hyperspectral imagery,”
IEEE Trans. on Image Processing, vol. 19, no. 6, pp.
  1403–1413, 2010.

eches2010estimating
O. Eches, N. Dobigeon, and J.-Y. Tourneret,
“Estimating the number of endmembers in hyperspectral images using
  the normal compositional model and a hierarchical Bayesian algorithm,”
IEEE Journal of Selected Topics in Signal Processing, vol. 4,
  no. 3, pp. 582–591, 2010.

stein2003application
D. Stein,
“Application of the normal compositional model to the analysis of
  hyperspectral imagery,”
in IEEE Workshop on Advances in Techniques for Analysis of
  Remotely Sensed Data, 2003, pp. 44–51.

zare2010pce
A. Zare and P. D. Gader,
“PCE: Piecewise convex endmember detection,”
IEEE Trans. on Geoscience and Remote Sensing, vol. 48, no. 6,
  pp. 2620–2632, 2010.

zhangpso
B. Zhang, L. Zhuang, L. Gao, W. Luo, Q. Ran, and Q. Du,
“PSO-EM: A hyperspectral unmixing algorithm based on normal
  compositional model,”
IEEE Trans. on Geoscience and Remote Sensing, vol. 52, no. 12,
  pp. 7782–7792, 2014.

du2014spatial
X. Du, A. Zare, P. D. Gader, and D. Dranishnikov,
“Spatial and spectral unmixing using the Beta compositional
  model,”
IEEE Journal of Selected Topics in Applied Earth Observations
  and Remote Sensing, vol. 7, no. 6, pp. 1994–2003, 2014.

zhougaussian
Y. Zhou, A. Rangarajan, and P. D. Gader,
“A Gaussian mixture model representation of endmember variability
  for spectral unmixing,”
in 2016 8th Workshop on Hyperspectral Image and Signal
  Processing: Evolution in Remote Sensing (WHISPERS), Aug 2016, pp. 1–5.

nascimento2005vertex
J. M. Nascimento and J. M. Bioucas Dias,
“Vertex component analysis: A fast algorithm to unmix hyperspectral
  data,”
IEEE Trans. on Geoscience and Remote Sensing, vol. 43, no. 4,
  pp. 898–910, 2005.

lu2013manifold
X. Lu, H. Wu, Y. Yuan, P. Yan, and X. Li,
“Manifold regularized sparse NMF for hyperspectral unmixing,”
IEEE Trans. on Geoscience and Remote Sensing, vol. 51, no. 5,
  pp. 2815–2826, 2013.

zare2013piecewise
A. Zare, P. D. Gader, O. Bchir, and H. Frigui,
“Piecewise convex multiple-model endmember detection and spectral
  unmixing,”
IEEE Trans. on Geoscience and Remote Sensing, vol. 51, no. 5,
  pp. 2853–2862, 2013.

murphy2012machine
K. P. Murphy,
Machine learning: a probabilistic perspective,
MIT press, 2012.

wetherley2017mapping
E. B. Wetherley, D. A. Roberts, and J. P. McFadden,
“Mapping spectrally similar urban materials at sub-pixel scales,”
Remote Sensing of Environment, vol. 195, pp. 170–183, 2017.

roberts2012synergies
D. A. Roberts, D. A. Quattrochi, G. C. Hulley, S. J. Hook, and R. O. Green,
“Synergies between VSWIR and TIR data for the urban environment:
  An evaluation of the potential for the hyperspectral infrared imager
  (HyspIRI) decadal survey mission,”
Remote Sensing of Environment, vol. 117, pp. 83–101, 2012.

schaaf2011mapping
A. N. Schaaf, P. E. Dennison, G. K. Fryer, K. L. Roth, and D. A. Roberts,
“Mapping plant functional types at multiple spatial resolutions
  using imaging spectrometer data,”
GIScience & Remote Sensing, vol. 48, no. 3, pp. 324–344,
  2011.

bishop2006pattern
C. M. Bishop,
Pattern recognition and machine learning,
springer New York, 2006.

meng1993maximum
X.-L. Meng and D. B. Rubin,
“Maximum likelihood estimation via the ECM algorithm: A general
  framework,”
Biometrika, vol. 80, no. 2, pp. 267–278, 1993.

mclachlan2014number
G. J. McLachlan and S. Rathnayake,
“On the number of components in a Gaussian mixture model,”
Wiley Interdisciplinary Reviews: Data Mining and Knowledge
  Discovery, vol. 4, no. 5, pp. 341–355, 2014.

smyth2000model
P. Smyth,
“Model selection for probabilistic clustering using cross-validated
  likelihood,”
Statistics and Computing, vol. 10, no. 1, pp. 63–72, 2000.

duchi2008efficient
J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra,
“Efficient projections onto the l1-ball for learning in high
  dimensions,”
in Proc. of the 25th International Conference on Machine
  learning (ICML). ACM, 2008, pp. 272–279.

zhou2016spatial
Y. Zhou, A. Rangarajan, and P. D. Gader,
“A spatial compositional model for linear unmixing and endmember
  uncertainty estimation,”
IEEE Trans. on Image Processing, vol. 25, no. 12, pp.
  5987–6002, 2016.

heylen2016hyperspectral
R. Heylen, A. Zare, P. Gader, and P. Scheunders,
“Hyperspectral unmixing with endmember variability via alternating
  angle minimization,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 54,
  no. 8, pp. 4983–4993, 2016.

zare2013sampling
A. Zare, P. Gader, and G. Casella,
“Sampling piecewise convex unmixing and endmember extraction,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 51,
  no. 3, pp. 1655–1665, 2013.

bland1986statistical
J. M. Bland and D. Altman,
“Statistical methods for assessing agreement between two methods of
  clinical measurement,”
The lancet, vol. 327, no. 8476, pp. 307–310, 1986.

gao2013comparative
L. Gao, Q. Du, B. Zhang, W. Yang, and Y. Wu,
“A comparative study on linear regression-based noise estimation for
  hyperspectral imagery,”
IEEE Journal of Selected Topics in Applied Earth Observations
  and Remote Sensing, vol. 6, no. 2, pp. 488–498, 2013.

hao2015semi
Z. Hao, M. Berman, Y. Guo, G. Stone, and I. Johnstone,
“Semi-realistic simulations of natural hyperspectral scenes,”
in IEEE International Geoscience and Remote Sensing Symposium
  (IGARSS). IEEE, 2015, pp. 1004–1007.

zhou2017nonrigidReg
Y. Zhou, A. Rangarajan, and P. D. Gader,
“Nonrigid registration of hyperspectral and color images with vastly
  different spatial and spectral resolutions for spectral unmixing and
  pansharpening,”
in 2017 IEEE Conference on Computer Vision and Pattern
  Recognition Workshops (CVPRW), July 2017, pp. 1571–1579.


[
    < g r a p h i c s >
]Yuan Zhou
 received the B.E degree in Software Engineering (2008), the M.E.
degree in Computer Application Technology (2011), both from Huazhong
University of Science and Technology, Wuhan, Hubei, China. Then he
worked in Shanghai UIH as a software engineer for two years. Since
2013, he has been a Ph.D. student in the Department of CISE, University
of Florida, Gainesville, FL, USA. His research interests include image
processing, computer vision and machine learning.


[
    < g r a p h i c s >
]Erin Wetherley
 is a doctoral candidate in the Geography Department at the University
of California, Santa Barbara. She specializes in characterizing urban
environments using hyperspectral imagery, thermal imagery, and sub-pixel
analyses. Prior to her doctoral work, Erin earned a bachelors degree
in Environmental Studies from Brown University, and worked for several
years as a GIS and database manager at a Washington, D.C. non-profit
organization.


[
    < g r a p h i c s >
]Paul Gader
 (M’86–SM’09–F’11)
received the Ph.D. degree in mathematics for image-processing-related
research from the University of Florida, Gainesville, FL, USA, in
1986. He was a Senior Research Scientist with Honeywell, a Research
Engineer and a Manager with the Environmental Research Institute of
Michigan, Ann Arbor, MI, USA, and a Faculty Member with the University
of Wisconsin, Oshkosh, WI, USA, the University of Missouri, Columbia,
MO, USA, and the University of Florida, FL, USA, where he is currently
a Professor of Computer and Information Science and Engineering. He
performed his first research in image processing in 1984 working on
algorithms for the detection of bridges in forward-looking infrared
imagery as a Summer Student Fellow at Eglin Air Force Base. He has
since worked on a wide variety of theoretical and applied research
problems including fast computing with linear algebra, mathematical
morphology, fuzzy sets, Bayesian methods, handwriting recognition,
automatic target recognition, biomedical image analysis, landmine
detection, human geography, and hyperspectral and light detection,
and ranging image analysis projects. He has authored/co-authored hundreds
of refereed journal and conference papers.



