# A consensus-constrained parsimonious Gaussian mixture model for clustering
hyperspectral images

Ganesh Babu<EMAIL_ADDRESS>School of Mathematics and Statistics,
Aoife Gowen<EMAIL_ADDRESS>School of Biosystems and Food Engineering
University College Dublin, Ireland.  Michael Fop<EMAIL_ADDRESS>School of
Mathematics and Statistics,  Isobel Claire Gormley<EMAIL_ADDRESS>
Isobel Claire Gormley and Michael Fop contributed equally to this work. School
of Mathematics and Statistics,

###### Abstract

The use of hyperspectral imaging to investigate food samples has grown due to
the improved performance and lower cost of spectroscopy instrumentation. Food
engineers use hyperspectral images to classify the type and quality of a food
sample, typically using classification methods. In order to train these
methods, every pixel in each training image needs to be labelled. Typically,
computationally cheap threshold-based approaches are used to label the pixels,
and classification methods are trained based on those labels. However,
threshold-based approaches are subjective and cannot be generalized across
hyperspectral images taken in different conditions and of different foods.
Here a consensus-constrained parsimonious Gaussian mixture model (ccPGMM) is
proposed to label pixels in hyperspectral images using a model-based
clustering approach. The ccPGMM utilizes available information on the labels
of a small number of pixels and the relationship between those pixels and
neighbouring pixels as constraints when clustering the rest of the pixels in
the image. A latent variable model is used to represent the high-dimensional
data in terms of a small number of underlying latent factors. To ensure
computational feasibility, a consensus clustering approach is employed, where
the data are divided into multiple randomly selected subsets of variables and
constrained clustering is applied to each data subset; the clustering results
are then consolidated across all data subsets to provide a consensus
clustering solution. The ccPGMM approach is applied to simulated datasets and
real hyperspectral images of three types of puffed cereal, corn, rice, and
wheat. Improved clustering performance and computational efficiency are
demonstrated when compared to other current state-of-the-art approaches.

##  1 Introduction

Hyperspectral imaging is a spectroscopic method combining digital imaging with
spectroscopy (Gowen et al., , 2019). Hyperspectral imaging collects an image
as a function of the light, in which each pixel is a vector of reflectance or
fluorescence information of the sample under study in the continuous spectrum.
The idea behind producing hyperspectral images lies in the interaction between
the photons emitted by a light source and the physical and chemical properties
of the sample (Amigo and Santos, , 2020). The interaction allows hyperspectral
images to capture critical information and traits of a sample not visible to
the naked eye. The critical information captured on the continuous spectrum
can be affected by spatial interference and spectral and redundant noise
(Amigo, , 2020). These issues are effectively handled with computationally
intensive spectral preprocessing techniques (Amigo, , 2010) like denoising and
scatter correction (Rinnan et al., , 2009; Bocklitz et al., , 2011; Geladi et
al., , 1985). The preprocessed information is then used to detect chemical
contamination, adulteration, and presence of foreign components in the food
samples under study (Feng and Sun, , 2012). Thanks to the increasing
availability of statistical methods for handling complex data, food engineers
are routinely employing classification methods for segmentation, food sample
detection, and analysis of large volumes of hyperspectral food images.

To train a classification method, each pixel of a hyperspectral image needs to
be appropriately labelled as e.g., either the background or the food sample.
Usually, the images are captured on appropriate backgrounds that facilitate
segmentation. In such cases, to select the regions of interest or the food
sample under study, a process called masking is performed to remove the
background pixels from the image (Amigo et al., , 2015). To create the mask,
multiple techniques are employed (Pal and Pal, , 1993), with the most widely
used being $k$-means clustering (Amigo et al., , 2015) and threshold-based
approaches (Pal and Pal, , 1993; Gowen et al., , 2019; Xu et al., , 2020). In
Amigo et al.,  (2015), $k$-means clustering with 2 clusters is applied to a
set of hyperspectral images to separate the background from the food sample.
Abdi and Williams,  (2010), Gowen et al.,  (2019), and Xu et al.,  (2020)
highlight some of the prominent threshold-based approaches used to label the
different food samples and background pixels in a given hyperspectral image.
In threshold-based approaches, the threshold for the segmentation of the
pixels is selected by producing a histogram of some useful pixel-related
summary of the entire spectrum, like the mean or the principal component
scores of the wavelength intensities for each pixel. The selection of which
summary to employ in the histogram depends on the type of background and the
food sample under study. In the histogram, typically a bimodal distribution is
observed and a threshold is subjectively chosen as the point that separates
the two modes. The quality of the threshold inferred labels in each
hyperspectral image depends on several decisions like the summary metric used
and which variables are selected for computing the summary. Furthermore, with
the threshold approach, only hard classification of pixels as background or
food type are possible, ignoring the inherent uncertainty in the labelling
process, particularly in certain areas of the image, like the pixels around
the edges of the food sample (Amigo and Santos, , 2020).

A major challenge in labelling the pixels of the hyperspectral images is that
each image has a large number of pixels and for each pixel, a large amount of
reflectance information is captured in a continuous $p$-dimensional spectrum,
whose wavelengths are highly correlated. Further, often threshold and
$k$-means clustering approaches are applied to each hyperspectral image
individually, which is impractical as the number of images increases and does
not allow for the borrowing of information across images. Also, the threshold
approach and $k$-means clustering with 2 clusters work well only when there is
exactly one type of food sample in a hyperspectral image. As the number of
food samples in a given image increases, multiple thresholds need to be
selected from the histogram, which is even more challenging. Overall,
threshold based approaches are subjective and are not easily generalisable
across hyperspectral images.

Approaches to clustering pixels in hyperspectral images are receiving
increasing attention (Zhai et al., , 2021). The density based DBSCAN
clustering approach (Khan et al., , 2014; Hahsler et al., , 2019) has been
widely used in imaging. Further, model-based clustering approaches have proved
effective, for example via Gaussian mixture models (GMM, Scrucca et al., ,
2016; McNicholas, , 2016; Bouveyron et al., , 2019) and parsimonious versions
(McNicholas and Murphy, , 2008). Bouveyron and Brunet-Saumard,  (2014) propose
a dimension reduction-based approach to segmenting hyperspectral images of
Mars. Further, Jacques and Ruckebusch,  (2016) consider co-clustering of
pixels and wavelengths in the context of hyperspectral imaging of an oil-in-
water emulsion using a latent block model.

Here, we propose a novel model-based clustering approach, a consensus-
constrained parsimonious Gaussian mixture model (ccPGMM), to label the pixels
of multiple hyperspectral images simultaneously in a computationally feasible
manner while accounting for uncertainty. A key aspect of ccPGMM is that it
involves a constrained parsimonious Gaussian mixture model (Melnykov et al., ,
2016) in which information regarding the types of food samples present in the
hyperspectral images along with information about some pixels is incorporated
as constraints to perform informed clustering. Further, as hyperspectral
images have a large number of pixels, each of which has an associated high-
dimensional spectrum of $p$ strongly correlated variables, ccPGMM employs a
parsimonious Gaussian mixture model (PGMM) which describes the high-
dimensional data using a small number of latent variables. To enable
computationally efficient implementation, ccPGMM adopts a consensus approach,
whereby the high-dimensional spectrum is divided into randomly selected
subsets of $d<p$ variables, and a constrained PGMM is fitted to each subset of
variables on all pixels. The clustering solutions of the constrained PGMMs
fitted to each subset are consolidated to provide a final consensus clustering
allocation for each pixel. The performance of the proposed ccPGMM approach is
demonstrated, and compared to other state-of-the-art approaches, through
thorough simulation studies and application to real hyperspectral images of
puffed cereals (Gowen et al., , 2019).

In what follows, Section 2 gives details of the motivating dataset containing
hyperspectral images of puffed cereals. Section 3 outlines the ccPGMM approach
and discusses its inference and implementation details. Section 4 delineates
thorough simulation studies that explore the performance of ccPGMM under
different settings. Section 5 discusses the application of ccPGMM to the
motivating puffed cereal images and Section 6 summarises the contributions of
ccPGMM and gives an overview of future research directions. The R (R Core
Team, , 2023) code to implement ccPGMM is available at this repository GitHub.

##  2 Hyperspectral images of puffed cereals

The ccPGMM approach is motivated by the need to simultaneously label pixels in
nine hyperspectral images of three types of puffed cereal, wheat, corn, and
rice (Gowen et al., , 2019). Each pixel in an image records the near-infrared
(NIR) spectrum, with $p=101$ equally spaced variables spanning the wavelengths
in the 880-1720nm interval. Each image $l=1,\ldots,L=9$ contains only one
puffed cereal and is a three-dimensional tensor of dimension $S_{l}\times
T_{l}\times p$ where $S_{l}$ and $T_{l}$ are the numbers of row and column
pixels respectively, which vary across the images as they are of different
dimensions (details are given in Appendix 8). Figure 1 shows a greyscale image
of one hyperspectral image of each puffed cereal type, where, for each pixel,
the intensity of the color is given by the average of the NIR spectrum.

For the purpose of ccPGMM, the pixels are assumed to be independent and so the
three-dimensional tensor of each image is represented by a rectangular dataset
of dimension $(S_{l}T_{l})\times p$. These $L=9$ rectangular datasets are then
collated to form a single dataset of dimension $N\times p$, where
$N=\sum_{l=1}^{L}S_{l}T_{l}$. Here, the final dataset of nine images results
in $N=28039$ pixels, each measured over $p=101$ wavelengths which are highly
correlated and contain spectral noise.

The objective is to simultaneously label all pixels in all nine images as
either background or a cereal type, taking into account available information
on the images (i.e., it is known which cereal type is in each image) and on
some pixels (e.g., corner pixels are background), in a computationally
feasible manner.

(a) Wheat

(b) Corn

(c) Rice

Figure 1: Greyscale images of one hyperspectral image of each puffed cereal
type. For each pixel, the intensity of the color corresponds to the average of
its NIR spectrum.

##  3 The ccPGMM and its inference

###  3.1 Parsimonious Gaussian mixture models (PGMM)

Factor analysis is a latent variable model that represents a large number $p$
of correlated variables using a smaller number $q\ll p$ of underlying latent
factors. In a single factor analysis model (McLachlan and Peel, , 2000), each
observation $\mathbf{x}_{n}\in\mathbb{R}^{p}$ for $n=1,\ldots,N$ in data
$\mathbf{X}=\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\\}$ of dimension $N\times
p$ is modelled as
$\mathbf{x}_{n}=\bm{\mu}+\bm{\Lambda}\mathbf{u}_{n}+\bm{\varepsilon}_{n}$
where $\bm{\mu}$ is the mean, $\bm{\Lambda}$ is the $p\times q$ loadings
matrix and $\mathbf{u}_{n}\sim\mathcal{N}_{q}(\mathbf{0},\mathbf{I})$ and
$\bm{\varepsilon}_{n}\sim\mathcal{N}_{p}(\mathbf{0},\mathbf{\Psi})$ are the
latent factor score and specific factor, respectively, for observation $n$.
Therefore, $\mathbf{x}_{n}\sim\mathcal{N}_{p}(\bm{\mu},\bm{\Sigma})$, where
${\bm{\Sigma}}=\bm{\Lambda\Lambda^{\top}+\Psi}$ provides a parsimonious
covariance structure.

The parsimonious Gaussian mixture model (PGMM) (McLachlan and Peel, , 2004;
McNicholas and Murphy, , 2008; McNicholas et al., , 2010) is a finite mixture
of $G$ factor analysis models such that each mixture component $g$ has a
multivariate Gaussian distribution with mean $\bm{\mu}_{g}$ and covariance
$\bm{\Lambda}_{g}\bm{\Lambda}_{g}^{\top}+\bm{\Psi}_{g}$. Denoting the
probability of membership of mixture component $g$ as $\tau_{g}$, under PGMM,
we have
${f}(\mathbf{x}_{n})=\sum_{g=1}^{G}\tau_{g}\>\mathcal{N}_{p}(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g}).$
By constraining $\bm{\Lambda}_{g}$ and $\bm{\Psi}_{g}$ to be the same or
different across components, a family of PGMMs have been developed (McNicholas
and Murphy, , 2008) with inference via an alternating expectation-conditional
maximization(AECM, Ghahramani et al., , 1996; McLachlan et al., , 2003). These
PGMMs facilitate model-based clustering of high-dimensional data.

###  3.2 Constrained parsimonious Gaussian mixture models

When labelling pixels in a set of hyperspectral images, there is often
information about the type of food sample(s) present in some of the
hyperspectral images, along with indirect information on the labels of some
pixels in the images. For example, in Figure 2, four blocks of pixels are
highlighted in three different colours in the two images. The pixels in the
blue blocks located at the corner of the images are assumed to belong to the
background. On the other hand, the pixels in the yellow and green blocks are
located in the center of their images, where the corresponding cereal is
observed to be situated. Pixels in the blue blocks should be clustered
together, as must the pixels in the yellow block and the pixels in the green
block; this information imposes what are known as positive constraints. Also,
the pixels in the blue blocks should not be placed in the same cluster as
pixels in the yellow or green blocks, and vice versa; this information imposes
what are known as negative constraints.

Melnykov et al.,  (2016) proposed a constrained Gaussian mixture model that
uses such additional information about some observations as constraints to
perform informed model-based clustering. However, the high-dimensionality of
the hyperspectral imaging data poses computational challenges for this
approach. Here, we build on Melnykov et al.,  (2016) and allow for constraints
to be incorporated when clustering using a PGMM. The resulting constrained
parsimonious Gaussian mixture model (constrained-PGMM) facilitates clustering
and therefore labelling of the pixels of the hyperspectral images.

(a) Wheat

(b) Corn

Figure 2: An illustration of positive ($+$) and negative ($-$) constraints.
Four blocks of three types of pixels are highlighted. Pixels in the blue
blocks (assumed to be the background) must be clustered together, as must the
pixels in the yellow block (wheat) and the pixels in the green block (corn).
However, pixels in the blue blocks must not be clustered with those in the
yellow block or with those in the green block, and vice versa.

Under the constrained-PGMM, the observed data log-likelihood is

$\ell_{o}(\mathbf{X})=\sum_{n=1}^{N}\log\sum_{g=1}^{G}\tau_{g}\>\mathcal{N}_{p}(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g}).$
(1)

As (1) is difficult to maximize directly, an AECM algorithm is used to fit the
constrained-PGMM. This requires the introduction of a latent component
indicator $\mathbf{z}_{n}=(z_{n1},\ldots,z_{nG})^{\prime}$, for
$n=1,\ldots,N$, where $z_{ng}=1$ if observation $n$ belongs to cluster $g$ or
0 otherwise. The AECM algorithm works in two cycles and allows for a different
definition of the complete data in each cycle. In the first cycle, the
component indicators $\mathbf{Z}$ = {$\mathbf{z}_{1}$,…, $\mathbf{z}_{N}$} are
assumed to be the missing data. The complete data log-likelihood is then

$\ell_{c}(\mathbf{X},\mathbf{Z})=\sum_{n=1}^{N}\sum_{g=1}^{G}{z}_{ng}[\log\tau_{g}+\log\>\mathcal{N}_{p}(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g})].$
(2)

In the E-step of the first cycle, the available constraints are incorporated
into PGMM. Let $B_{1}$ be the set of $J$ pixels (indexed by $j$) in the blue
blocks shown in Figure 2 which are known to be background and must be
clustered together. Let $B_{2}$ be the set of $O$ pixels (indexed by $o$) in
the yellow block of Figure 2 which are known to be wheat and should be
clustered together. These are positive constraints (denoted $+$). However, the
pixels in $B_{1}$ should not be clustered together with the pixels in $B_{2}$.
This is a negative constraint (denoted $-$). Thus the posterior probability
$\hat{{z}}_{{B_{1}},g}^{+\atop-}$ of the pixels in $B_{1}$ belonging to
cluster $g$, given the positive and negative constraints $+\atop-$ is

$\centering\indent\indent\hat{{z}}_{{B_{1}},g}^{+\atop-}=\frac{\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq g\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}{\displaystyle\sum_{g^{{}^{\prime}}=1}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g^{{}^{\prime}}}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq
g^{{}^{\prime}}\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})},\@add@centering$
(3)

where $\hat{\cdot}$ denotes an initial or current parameter value, as
relevant. In (3), the posterior probabilities of pixels in $B_{1}$ are
considered with only one set of pixels $B_{2}$ as a negative constraint,
however (3) can easily be extended to accommodate multiple such negative
constraints. The remaining pixels in Figure 2 which are not in the highlighted
blocks are considered as blocks with only one pixel and no negative
constraints. The estimated model parameters $\hat{\bm{\mu}}_{g}$ and
$\hat{\tau}_{g}$ are updated in the M-step of the first cycle as in PGMM(see
McNicholas and Murphy, , 2008). Further details are in Appendix 9.

In the E-step of the second cycle, the expected value of
$\hat{{z}}_{B_{1},g}^{+\atop-}$ is computed as in (3), with
$\hat{\bm{\mu}}_{g}$ and $\hat{\tau}_{g}$ updated to those from the M-step of
the first cycle. The expected value of $\mathbf{u}_{n}\>,n=1,\ldots,N$, and
the estimates of $\hat{\bm{\Lambda}}_{g}$ and
$\hat{\bm{\Psi}}_{g}\>,g=1,\ldots,G$ are computed as in PGMM; details are
reported in Appendix 9.

The AECM algorithm moves between these cycles, iteratively updating the model
parameters, until convergence. The values of $\hat{\mathbf{Z}}$ at convergence
are the posterior probabilities of cluster membership for each pixel; labels
can be obtained based, for example, on the _maximum a posteriori_ values of
$\hat{\mathbf{Z}}$. However, constrained-PGMM exhibits very long run times
when fitted to the hyperspectral image data and could not be fitted to the
entire $p=101$ variables of the motivating dataset in a computationally
efficient manner.

###  3.3 Consensus-constrained parsimonious Gaussian mixture models

To improve the computational efficiency of constrained-PGMM, we propose a
consensus-constrained parsimonious Gaussian mixture model (ccPGMM), inspired
by Russell et al.,  (2015). Russell et al.,  (2015) proposed a consensus
approach to avoid the selection of a single best model based on the Bayesian
information criterion (BIC, Neath and Cavanaugh, , 2012; Schwarz, , 1978). In
Russell et al.,  (2015), the clustering solutions of multiple models fitted on
the same data are consolidated to provide the final consensus clustering
solution. Previously, Strehl and Ghosh,  (2002) proposed three ensemble
clustering methods namely a cluster-based similarity partitioning algorithm, a
hyper-graph partitioning algorithm, and a meta-clustering algorithm to combine
multiple clustering solutions into one, but the uncertainty in cluster
membership was not accounted for. Fern and Brodley,  (2003) and Punera and
Ghosh,  (2008) extended these approaches, incorporating the posterior
probabilities into the ensemble.

Unlike Russell et al.,  (2015), the proposed ccPGMM takes a divide-and-conquer
approach. In ccPGMM, a constrained-PGMM is fitted to $M$ subsets (indexed by
$m$) of $d<p$ randomly selected variables from the high dimensional NIR
spectrum. The $M$ posterior probabilities which account for the uncertainty of
the cluster memberships of the $N$ pixels are then consolidated to provide the
final consensus clustering solution as follows. An $N\times N$ similarity
matrix $\mathbf{S}^{m}$ for $m=1,\ldots,M$ is computed based on the estimated
posterior probabilities $\hat{\mathbf{Z}}^{m}$ available on convergence of the
AECM after fitting the constrained-PGMM to the subset $m$. Each entry of
$\mathbf{S}^{m}$ is a similarity score $S_{ij}^{m}$ computed as

$S_{ij}^{m}=\begin{cases}\hat{\mathbf{z}}_{i}^{m}(\hat{\mathbf{z}}_{j}^{m})^{\top}&\text{if
}i\neq j,\\\ 1&\text{if }i=j,\end{cases}$

where $\hat{\mathbf{z}}_{i}^{m}$ and $\hat{\mathbf{z}}_{j}^{m}$ are the
posterior probabilities of cluster membership for observations $i$ and $j$
respectively. The $M$ similarity matrices are then averaged to compute the
final similarity matrix $\mathbf{S}$, whose entries $S_{ij}$ are given by

${S}_{ij}=\frac{1}{M}\sum_{m=1}^{M}\displaystyle\sum_{g=1}^{G}\hat{z}_{ig}^{m}\hat{z}_{jg}^{m},$

where $\hat{z}_{ig}^{m}$ and $\hat{z}_{jg}^{m}$ are the posterior
probabilities of observation $i$ and $j$ belonging to group $g$, inferred from
data subset $m$. A dissimilarity matrix $\mathbf{D}$ is then computed as
$\mathbf{D}=1-\mathbf{S}$ and hierarchical clustering with complete linkage
(Sokal, , 1963) is performed employing $\mathbf{D}$. As the number of clusters
$G$ is known based on the application at hand, the resulting dendrogram is cut
to give $G$ clusters and the final consensus clustering solution.

###  3.4 Technical details

In the proposed ccPGMM, a single constrained-PGMM with a fixed number of
clusters $G$ and a fixed number of factors $q$ is fitted to $M$ data subsets
with $d$ variables per subset; the values of $G$, $q$, $M$, and $d$ require
specification. The number of clusters $G$ is known and fixed in advance
according to the particular application. For example, in the motivating
application here (see Section 2) where there are three cereal types and the
background, $G=4$. The selection of the number of factors $q$ is based on
monitoring the proportion of variance explained by the first $q$ principal
components obtained from the application of principal component analysis (PCA)
to the entire dataset. There is a trade-off between the values of $M$, the
number of subsets of variables, and $d$, the number of variables per subset.
Typically, $d$ is smaller than $p$ but too small values may provide little
clustering information. Further, very large values of $d$ hamper the
computational feasibility. Also, the choice of $M$ must be cognisant of the
choice of $d$ to ensure good coverage of all the variables in the $M$ data
subsets. These choices are explored further in sections 4 and 5.

In terms of specifying the constraints, there are technically no limits on the
number of pixels that can be used. The size of the blocks of pixels
(constraints) highlighted in Figure 2 can vary and the number of pixels in one
block need not be the same as the number of pixels in the other blocks.
Usually, the hyperspectral images of food samples are plotted as greyscale
images, and the maximum region of pixels that are distinguishable as the
background and as the relevant sample under study are selected as constraints.
Increasing the number of positive constraints and the number of pixels within
the positive constraints does not significantly increase the computational
cost of fitting the ccPGMM, as in Melnykov et al.,  (2016). However, an
increase in the number of negative constraints will result in more complex
posterior probability calculations in (3), negatively impacting the
computational cost. Investigations on the choice of constrained pixels are in
sections 4 and 5.

The AECM algorithm, used to fit the constrained-PGMM to each data subset, is
initialized with the cluster labels obtained from applying $k$-means
clustering on the respective data subset. The initial values of
$\bm{\Lambda}_{g}$ and $\bm{\Psi}_{g}$ are computed based on the eigenvalue
decomposition of $\bm{\Sigma}_{g}$ (McNicholas and Murphy, , 2008). The
relative change in the log-likelihood is used to assess the convergence of the
AECM with a tolerance of $1e-6$.

##  4 Simulation studies

Simulation studies are conducted to assess the performance of ccPGMM and
compare it to the performance of the threshold approach and existing state-of-
the-art methods such as DBSCAN (Khan et al., , 2014; Hahsler et al., , 2019),
Gaussian mixture models (GMM, Scrucca et al., , 2016; McNicholas, , 2016;
Bouveyron et al., , 2019), PGMM (McNicholas and Murphy, , 2008) and consensus
PGMM (cPGMM) i.e., ccPGMM with no constraints. The convergence of each model-
based clustering approach is assessed with the default parameter settings
detailed in the respective R packages. Four scenarios are considered: well-
separated clusters (scenario 1), heavily overlapping clusters (scenario 2),
mildly overlapping clusters (scenario 3), and synthetic puffed cereal data
(scenario 4). For each scenario, five datasets are simulated.

In scenarios 1–3, each simulated dataset contains three simulated
hyperspectral images (one for each cereal type), giving a total of $N=7825$
pixels, where each pixel has an associated $p=101$ spectrum of strongly
correlated variables. The three images are simulated according to labels
resulting from applying the threshold approach to three of the hyperspectral
images in the motivating dataset detailed in Section 2. The $p$ variables for
each of the three types of cereal and the background pixels are generated from
different factor analysis models with $q=3$ and parameters $\bm{\mu}_{g}$,
$\bm{\Lambda}_{g}$, and $\bm{\Psi}_{g}$ for $g=1,\ldots,G=4$. To emulate the
real data, a strong correlation is induced between the $p$ variables at each
pixel by, for every set of $5$ consecutive variables, generating the
associated values in the loadings matrix from a
$\mathcal{N}_{q}(\mathcal{U}(0.3,0.9),0.03)$ and the values of
$\mbox{diag}(\bm{\Psi}_{g})$ from a $\Gamma(\mathcal{U}(0,0.1),1)$. The values
of the means $\bm{\mu}_{g}$ for $g=1,\ldots,G$ are used to control the degree
of overlap of the clusters in each scenario; details are given in the
respective subsections of Section 4. To determine the constraints, the
greyscale images of the simulated hyperspectral images are examined to decide
on the sets of pixels to be used as constraints, as shown in Figure 3. Here,
this resulted in $2,956$ pixels of the $N=7825$ pixels $(37.7\%)$ being used
to inform the constraints. The $1,356$ pixels in the blue blocks (background)
must be clustered together, as must the $750$ pixels in the yellow block
(wheat), the $450$ pixels in the green block (corn), and the $400$ pixels in
the red block (rice). These are the positive constraints. The pixels in one
coloured block must not be clustered together with the pixels in another
coloured block; these are the negative constraints.

(a) Wheat

(b) Corn

(c) Rice

Figure 3: Greyscale image of three simulated hyperspectral images for one of
the five well-separated cluster datasets. The pixels in the blue blocks are
background, the pixels in the yellow block are wheat, the pixels in the green
block are corn and the pixels in the red block are rice.

In scenario 4, a synthetic cereal dataset is considered which mirrors the
properties of the motivating hyperspectral images. Each dataset contains $L=9$
simulated hyperspectral images (three for each cereal type) with $N=28,039$
and where each pixel has an associated $p=101$ length spectrum. The $p$
variables for each cereal type and background pixels are generated from the
relevant factor analyzer model, fitted to the real data; details are discussed
in Section 4.4. In this scenario, $12,215$ pixels of the $N=28,039$ pixels
$(43.5\%)$ are selected as constraints by plotting the greyscale images of the
synthetic puffed cereals. Of the $12,215$ pixels selected as constraints,
$5,900$ are background, $3,130$ are wheat, $1,535$ are corn, and $1,650$ are
rice. Further, to assess the impact of the proportion of pixels used as
constraints on the clustering performance, a smaller set of $6,951$ pixels
($24.7\%$) are selected as constraints. Of the $6,951$ pixels, $3,801$ are
background, $1,540$ are wheat, $783$ are corn and $827$ are rice. The set of
the $43.5\%$ of pixels and the $24.7\%$ of pixels used as constraints are
illustrated in Figures 13 and 14 respectively in Appendix 10.

The simulation study is conducted on an Intel(R) core(TM) i7-10850H CPU @ 2.7
GHz processor with 16 GB RAM and 64-bit operating system. The code used to
generate the datasets for each scenario and to conduct the simulation studies
is available on GitHub.

###  4.1 Scenario 1: well-separated clusters

To generate the well-separated clusters, for $g=1,\ldots,4$, the cluster mean
parameters are generated as $\bm{\mu}_{g}\sim\mathcal{N}_{p}(a,0.5)$, where
$a\in\\{0,5,10,-5\\}$; the $p$ variables of the background, wheat, corn, and
rice pixels in each dataset are generated from the four well-separated factor
analyzer models, respectively. For DBSCAN, the minimum points required in the
neighborhood for core points are set to $2p$, and for the GMM, $G=4$ is
determined by the setting. For the PGMM approaches, $G=4$ and $q=1$, where $q$
is informed by fitting PCA to the entire dataset and monitoring the proportion
of variance explained by the components. For consensus based approaches,
$M=\\{10,25,50,100\\}$ randomly selected subsets of variables are considered
with the number of variables per subset $d=\\{10,20\\}$. The positive and
negative constraints for ccPGMM are selected as shown in Figure 3.

Figure 4: For each of the five simulated datasets with well-separated
clusters, the ARI between the known labels and the clustering solutions of
DBSCAN, GMM, and PGMM fitted on $p=101$ variables (first three panels from the
left) and cPGMM and ccPGMM (last two panels on the right) fitted with
different settings of $M$ and $d$. Figure 5: For each of the five simulated
datasets with well-separated clusters, the time taken to fit DBSCAN, GMM, and
PGMM on $p=101$ variables (first three panels from the left) and cPGMM and
ccPGMM (last two panels on the right) fitted with different settings of $M$
and $d$.

Figure 4 shows the adjusted Rand index(ARI, Hubert and Arabie, , 1985) between
the known labels and the clustering solutions for the five simulated datasets
in scenario 1 under the DBSCAN, GMM, PGMM, cPGMM, and ccPGMM approaches.
Figure 5 shows the time taken in minutes to fit the methods. All approaches
demonstrate strong clustering performance. In terms of computational cost,
DBSCAN clusters each dataset in less than a minute while GMM and PGMM take an
average of $27$ minutes and $15$ minutes respectively. The cPGMM and ccPGMM
fitted with different settings of $M$ and $d$ have lower computational costs
than GMM and PGMM, but similar to DBSCAN, with no loss of accuracy.

###  4.2 Scenario 2: heavily overlapping clusters

To generate data with heavily overlapping clusters, the cluster mean
parameters are generated as $\bm{\mu}_{g}\sim\mathcal{N}_{p}(a,1.5)$, where
$a\in\\{1,1.5,3,-1.5\\}$ for $g=1,\ldots,4$ respectively. Figures 6 and 7 show
the ARI between the known pixel labels and the clustering solutions of DBSCAN,
GMM, PGMM, cPGMM, and ccPGMM and the time taken in minutes by each method,
respectively. DBSCAN clusters the datasets with an average ARI of $0.34$ in
less than a minute on average. The GMM struggles to uncover the clustering
structure. The PGMM approach performs competitively, achieving an average ARI
of $0.40$ in $5$ minutes on average. The ccPGMM approach performs well in this
scenario: with $M\geq 25$ and $d=20$ an average ARI of $0.66$ is achieved,
which outperforms PGMM fitted to all $p$ variables. In terms of clustering
performance, ccPGMM performs better than cPGMM across all settings of $M$ and
$d$, in comparable time.

Figure 6: For each of the five datasets with heavily overlapping clusters, the
ARI of DBSCAN, GMM, and PGMM fitted on $p=101$ variables (first three panels
from the left) and cPGMM and ccPGMM (last two panels on the right) with
different settings of $M$ and $d$. Figure 7: For each of the five datasets
with heavily overlapping clusters, the time taken to fit DBSCAN, GMM, and PGMM
on $p=101$ variables (first three panels from the left) and cPGMM and ccPGMM
(last two panels on the right) with different settings of $M$ and $d$.

###  4.3 Scenario 3: mildly overlapping clusters

The $p$ variables of the background, wheat, corn, and rice pixels in each
dataset are generated from four mildly overlapping factor analyzer models,
with the cluster mean parameters generated as
$\bm{\mu}_{g}\sim\mathcal{N}(a,1.5)$ where $a=\\{0,2,5,-2\\}$ for
$g=1,\ldots,4$, respectively. Figures 8 and 9 show the ARI of the cluster
solutions of DBSCAN, GMM, PGMM, cPGMM, and ccPGMM fitted on the five datasets
with mildly overlapping clusters and the time taken in minutes to fit them.
While GMM again struggles to uncover the clustering structure, DBSCAN and PGMM
fitted to all $p$ variables clustered the pixels with an average ARI of $0.42$
and $0.64$ respectively. In terms of cPGMM and ccPGMM, Figure 8 shows that
increasing $d$ from $10$ to $20$ significantly improves the clustering
performance of both methods, with the performance of ccPGMM being the best
among the approaches considered. Run times are comparable overall, with cPGMM
with the highest values of $M$ and $d$ intuitively taking the longest. The
ccPGMM approach with $M\geq 10$ and $d=20$ is competitively computationally
efficient with an average ARI $>0.80$, which has again a much better
performance than PGMM fitted to all $p$ variables.

Figure 8: For each of the five datasets with mildly overlapping clusters, the
ARI of DBSCAN, GMM, and PGMM fitted on $p=101$ variables (first three panels
from the left) and cPGMM and ccPGMM (last two panels on the right) fitted with
different settings of $M$ and $d$. Figure 9: For each of the five datasets
with mildly overlapping clusters, the time taken to fit DBSCAN, GMM, and PGMM
on $p=101$ variables (first three panels from the left) and cPGMM and ccPGMM
(last two panels on the right) fitted with different settings of $M$ and $d$.

###  4.4 Scenario 4: synthetic puffed cereal data

Here a synthetic cereal dataset is considered which mirrors the properties of
the motivating hyperspectral images. To generate the data, PCA is fitted to
the real puffed cereal images and, based on the cumulative proportion of
variance explained (Table 2 in Appendix 11), $q=2$ is selected. Threshold
labels, generated as proposed by Xu et al.,  (2020), are used to identify the
pixels belonging to the 3 cereal types and the background in the real puffed
cereal images. Four individual factor analysis models with $q=2$ are then
fitted to pixels in the motivating dataset from each label type. Finally, to
generate the synthetic cereal images, for each pixel, $p$ variables are
generated from the resulting factor analysis models fitted to the background,
wheat, corn, and rice pixels in real data.

Figures 10 and 11 show, respectively, the ARI of the cluster solutions of the
different clustering methods fitted to the five synthetic cereal datasets, and
the time taken in minutes for estimation. Intuitively, since the synthetic
cereal datasets are generated from a mixture of factor analyzers, PGMM fitted
to all $p$ variables shows strong clustering performance (average ARI of
$0.94$), however, it takes an average of $50$ hours to cluster each synthetic
cereal dataset. The GMM approach has a much lower computational cost than PGMM
but the quality of the clustering solutions is inconsistent. For DBSCAN (not
illustrated), while computationally cheap, clustering performance is poor with
an average ARI of $0.32$ with a variance of $0.11$. The ccPGMM with $43.5\%$
of pixels as constraints and ccPGMM with $24.7\%$ of pixels as constraints
fitted with $M\geq 10$ and $d=20$ shows strong clustering performances.
Intuitively, ccPGMM with larger set of constraints ($43.5\%$) performs
slightly better than ccPGMM with smaller set of constraints ($24.7\%$) and
cPGMM across different settings of $M$ and $d$. For the ccPGMM approach with
$43.5\%$ pixels as constraints, fitted with $M\geq 25$ and $d=20$ has an
average ARI of $0.92$. Overall, ccPGMM shows comparable clustering performance
to PGMM fitted to all $p$ variables, in less than one-tenth of the
computational time.

Figure 10: For each of the five synthetic cereal datasets, the ARI of GMM, and
PGMM fitted on $p=101$ variables (first two panels from the left) and cPGMM
and ccPGMM (last two panels on the right) fitted with different settings of
$M$ and $d$. For ccPGMM, two settings with $24.7\%$ and $43.5\%$ of pixels
used as constraints are shown. Figure 11: For each of the five synthetic
cereal datasets, the time taken to fit GMM, and PGMM on $p=101$ variables
(first two panels from the left) and cPGMM and ccPGMM (last two panels on the
right) fitted with different settings of $M$ and $d$. For ccPGMM, two settings
with $24.7\%$ and $43.5\%$ of pixels used as constraints are shown.

##  5 Application to hyperspectral images of puffed cereals

To simultaneously cluster the pixels of the $9$ hyperspectral images in the
motivating dataset, we apply ccPGMM with $G=4$, as there are $3$ puffed cereal
types (wheat, corn, and rice) and the background, and $q=2$ based on Table 2
in Appendix 11. On the basis of the performance observed in simulation study
scenario 4 (Section 4.4) which mirrored the real data, $M=25$ subsets and
$d=20$ variables per subset were employed. As detailed in Section 4, the
greyscale version of the hyperspectral images were considered and, based on
pixels for which the labelling was clear, $12,215$ or $(43.5\%)$ of the
$N=28,039$ pixels were selected; these pixels informed the constraints. Also,
a smaller set of $6,951$ pixels out of the $N=28,039$ pixels ($24.7\%$) are
selected as constraints to assess the impact of the proportion of pixels used
as constraints on the clustering performance. For comparison purposes, DBSCAN,
GMM with $G=4$, and PGMM with $G=4$ and $q=2$ are fitted to the motivating
dataset with all $p$ variables. The threshold approach of Xu et al.,  (2020)
is also applied, individually, to each hyperspectral image to label the
pixels.

Figure 12 shows the greyscale images for three of the nine hyperspectral
images, and the pixel labels as generated using the threshold approach,
DBSCAN, GMM, PGMM, ccPGMM with $43.5\%$ pixels as constraints and ccPGMM with
$24.7\%$ pixels as constraints; Figures 15 and 16 in Appendix 12 illustrate
the same for the remaining six hyperspectral images. Figures 12(d), 12(e) and
12(f) show the labels generated for each hyperspectral image by the threshold
approach. Each type of pixel is coloured in a different shade of blue
according to its label; the resulting labels are very close to the respective
greyscale images (12(a), 12(b), and 12(c)). However, the threshold approach is
ad-hoc involving subjective choice, and does not provide any clustering
uncertainties. Figures 12(g), 12(h) and 12(i) show the cluster solutions of
DBSCAN where each pixel is coloured in a different shade of green depending on
its label. While DBSCAN correctly distinguishes the background pixels from the
cereal pixels in each hyperspectral image, the pixels of the three different
cereals are incorrectly allocated to the same cluster. Similarly, GMM clusters
the pixels of the three types of cereals into a single cluster as shown in
Figures 12(j), 12(k), and 12(l). In addition, GMM allocates the uncertain edge
pixels of the cereal grains to a separate cluster and allocates some pixels
that likely correspond to spectral noise into another separate cluster.
Figures 12(m), 12(n), and 12(o) illustrate the pixel labels based on the
application of PGMM to all $p$ variables. Each label obtained from the fitted
PGMM is coloured in a different shade of purple. Though PGMM has correctly
allocated several background and cereal pixels to their respective separate
clusters, it also appears to have incorrectly clustered some together, giving
non-contiguous images. In Figure 12(n), which shows a corn cereal, many
background pixels are allocated to the corn cluster. In Figure 12(m) and
12(o), corresponding to wheat and rice cereals respectively, several edge
pixels of the cereals are incorrectly allocated to the corn cluster. In
addition, PGMM takes nearly $72$ hours to cluster the motivating dataset.

Finally, Figures 12(p), 12(q) and 12(r) show the labels of the pixels, in
different shades of brown, obtained by applying ccPGMM with $43.5\%$ of pixels
as constraints. Applying ccPGMM to all nine hyperspectral images
simultaneously took only $6$ hours and the cluster labels are close to labels
obtained by applying the threshold approach to each image individually.
Compared to the other clustering solutions considered, ccPGMM distinguishes
the cereals well from the background while, overall, pixels are allocated to
different clusters corresponding to corn, wheat, and rice. However, some
pixels are still incorrectly mixed across clusters, particularly around the
edges of the wheat and rice cereals. Figures 12(s), 12(t) and 12(u) show the
cluster labels obtained by applying ccPGMM with $24.7\%$ pixels as constraints
in different shades of brown. The ccPGMM fitted with $24.7\%$ of pixels as
constraints cluster the pixels as well as the ccPGMM with $43.5\%$ of pixels
as constraints. However, some pixels around the edges of the wheat and corn
cereal are incorrectly clustered as rice. Importantly, the uncertainties
associated with the cluster membership for each pixel are available given the
model-based approach to clustering and the uncertainty plots are given in
Figures 17, 18, and 19 in Appendix 13; intuitively the pixels around the edges
of the cereals have the highest clustering uncertainties.

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

Figure 12: Greyscale images and pixel labels for three of the nine
hyperspectral images ($l=\\{1,4,7\\}$) of the cereals using the threshold
approach, DBSCAN, MBC, PGMM, ccPGMM with $43.5\%$ of pixels as constraints,
and ccPGMM with $24.7\%$ of pixels as constraints respectively for wheat
(12(a), 12(d), 12(g), 12(j), 12(m), 12(p), 12(s)), corn (12(b), 12(e), 12(h),
12(k), 12(n), 12(q), 12(t)) and rice (12(c), 12(f), 12(i), 12(l), 12(o),
12(r), 12(u)).

##  6 Discussion

A novel model-based approach to simultaneously labelling pixels of
hyperspectral images while availing of additional information about the
relationship between groups of pixels within the images is proposed. The
proposed ccPGMM demonstrates strong clustering performance when compared to
existing approaches both in terms of quality of the clustering solution, and
the computational time required to cluster such a large dataset.

While the ccPGMM approach is less subjective than the currently used threshold
approach, user decisions must still be made. For example, the selection of the
number of factors $q$, while well established, is subjective. Examining
alternative approaches to choosing $q$ either using model selection criteria
or using shrinkage priors (Bhattacharya and Dunson, , 2011) in a Bayesian
context would address this, but increase the computational cost. Indeed, under
the ccPGMM approach, allowing for different values of $q$ when fitting to the
$M$ different subsets of variables could bring improved performance but would
require further exploration. Further, the choice of $M$ and $d$ are currently
data-driven and therefore subjective. Additionally, similar to Jacques and
Ruckebusch,  (2016), considering co-clustering in the context of ccPGMM would
facilitate both clustering of pixels and wavelengths.

Overall, the proposed ccPGMM model works well when clustering a large number
of pixels in a set of hyperspectral images, and provides associated clustering
uncertainty, in a computationally efficient manner. As such, the approach
should have a broader application to other hyperspectral images where
labelling of the pixels is required.

##  7 Acknowledgement

This publication has emanated from research conducted with the financial
support of Science Foundation Ireland under Grant number 18/CRT/6049. For the
purpose of Open Access, the author has applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission.

## References

  * Abdi and Williams, (2010) Abdi, H. and Williams, L. J. (2010).  Principal component analysis.  Wiley interdisciplinary reviews: computational statistics, 2(4):433–459.
  * Amigo, (2010) Amigo, J. M. (2010).  Practical issues of hyperspectral imaging analysis of solid dosage forms.  Analytical and bioanalytical chemistry, 398(1):93–109.
  * Amigo, (2020) Amigo, J. M. (2020).  Hyperspectral and multispectral imaging: Setting the scene.  In Data Handling in Science and Technology, volume 32, pages 3–16. Elsevier.
  * Amigo et al., (2015) Amigo, J. M., Babamoradi, H., and Elcoroaristizabal, S. (2015).  Hyperspectral image analysis. a tutorial.  Analytica chimica acta, 896:34–51.
  * Amigo and Santos, (2020) Amigo, J. M. and Santos, C. (2020).  Preprocessing of hyperspectral and multispectral images.  In Data handling in science and technology, volume 32, pages 37–53. Elsevier.
  * Bhattacharya and Dunson, (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse Bayesian infinite factor models.  Biometrika, 98(2):291–306.
  * Bocklitz et al., (2011) Bocklitz, T., Walter, A., Hartmann, K., Rösch, P., and Popp, J. (2011).  How to pre-process Raman spectra for reliable and stable models?  Analytica chimica acta, 704(1-2):47–56.
  * Bouveyron and Brunet-Saumard, (2014) Bouveyron, C. and Brunet-Saumard, C. (2014).  Discriminative variable selection for clustering with the sparse fisher-em algorithm.  Computational Statistics, 29(3):489–513.
  * Bouveyron et al., (2019) Bouveyron, C., Celeux, G., Murphy, T. B., and Raftery, A. E. (2019).  Model-based clustering and classification for data science: with applications in R, volume 50.  Cambridge University Press.
  * Feng and Sun, (2012) Feng, Y.-Z. and Sun, D.-W. (2012).  Application of hyperspectral imaging in food safety inspection and control: a review.  Critical reviews in food science and nutrition, 52(11):1039–1058.
  * Fern and Brodley, (2003) Fern, X. Z. and Brodley, C. E. (2003).  Random projection for high dimensional data clustering: A cluster ensemble approach.  In Proceedings of the 20th international conference on machine learning (ICML-03), pages 186–193.
  * Geladi et al., (1985) Geladi, P., MacDougall, D., and Martens, H. (1985).  Linearization and scatter-correction for near-infrared reflectance spectra of meat.  Applied spectroscopy, 39(3):491–500.
  * Ghahramani et al., (1996) Ghahramani, Z., Hinton, G. E., et al. (1996).  The EM algorithm for mixtures of factor analyzers.  Technical report, Technical Report CRG-TR-96-1, University of Toronto.
  * Gowen et al., (2019) Gowen, A. A., Xu, J.-L., and Herrero-Langreo, A. (2019).  Comparison of spectral selection methods in the development of classification models from visible near infrared hyperspectral imaging data.  Journal of Spectral Imaging, 8.
  * Hahsler et al., (2019) Hahsler, M., Piekenbrock, M., and Doran, D. (2019).  DBscan: Fast density-based clustering with R.  Journal of Statistical Software, 91:1–30.
  * Hubert and Arabie, (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of classification, 2:193–218.
  * Jacques and Ruckebusch, (2016) Jacques, J. and Ruckebusch, C. (2016).  Model-based co-clustering for hyperspectral images.  Journal of Spectral Imaging.
  * Khan et al., (2014) Khan, K., Rehman, S. U., Aziz, K., Fong, S., and Sarasvady, S. (2014).  DBscan: Past, present and future.  In The fifth international conference on the applications of digital information and web technologies (ICADIWT 2014), pages 232–238. IEEE.
  * McLachlan and Peel, (2000) McLachlan, G. and Peel, D. (2000).  Mixtures of factor analyzers.  In In Proceedings of the Seventeenth International Conference on Machine Learning. Citeseer.
  * McLachlan and Peel, (2004) McLachlan, G. J. and Peel, D. (2004).  Finite mixture models.  John Wiley & Sons.
  * McLachlan et al., (2003) McLachlan, G. J., Peel, D., and Bean, R. W. (2003).  Modelling high-dimensional data by mixtures of factor analyzers.  Computational Statistics & Data Analysis, 41(3-4):379–388.
  * McNicholas, (2016) McNicholas, P. D. (2016).  Model-based clustering.  Journal of Classification, 33(3):331–373.
  * McNicholas and Murphy, (2008) McNicholas, P. D. and Murphy, T. B. (2008).  Parsimonious Gaussian mixture models.  Statistics and Computing, 18(3):285–296.
  * McNicholas et al., (2010) McNicholas, P. D., Murphy, T. B., McDaid, A. F., and Frost, D. (2010).  Serial and parallel implementations of model-based clustering via parsimonious Gaussian mixture models.  Computational Statistics & Data Analysis, 54(3):711–723.
  * Melnykov et al., (2016) Melnykov, V., Melnykov, I., and Michael, S. (2016).  Semi-supervised model-based clustering with positive and negative constraints.  Advances in data analysis and classification, 10(3):327–349.
  * Neath and Cavanaugh, (2012) Neath, A. A. and Cavanaugh, J. E. (2012).  The Bayesian information criterion: background, derivation, and applications.  Wiley Interdisciplinary Reviews: Computational Statistics, 4(2):199–203.
  * Pal and Pal, (1993) Pal, N. R. and Pal, S. K. (1993).  A review on image segmentation techniques.  Pattern recognition, 26(9):1277–1294.
  * Punera and Ghosh, (2008) Punera, K. and Ghosh, J. (2008).  Consensus-based ensembles of soft clusterings.  Applied Artificial Intelligence, 22(7-8):780–810.
  * R Core Team, (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rinnan et al., (2009) Rinnan, Å., Van Den Berg, F., and Engelsen, S. B. (2009).  Review of the most common pre-processing techniques for near-infrared spectra.  TrAC Trends in Analytical Chemistry, 28(10):1201–1222.
  * Russell et al., (2015) Russell, N., Murphy, T. B., and Raftery, A. E. (2015).  Bayesian model averaging in model-based clustering and density estimation.  arXiv preprint arXiv:1506.09035.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * Scrucca et al., (2016) Scrucca, L., Fop, M., Murphy, T. B., and Raftery, A. E. (2016).  mclust 5: clustering, classification and density estimation using Gaussian finite mixture models.  The R Journal, 8(1):289–317.
  * Sokal, (1963) Sokal, R. R. (1963).  The principles and practice of numerical taxonomy.  Taxon, pages 190–199.
  * Strehl and Ghosh, (2002) Strehl, A. and Ghosh, J. (2002).  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  Journal of machine learning research, 3(Dec):583–617.
  * Xu et al., (2020) Xu, J.-L., Riccioli, C., Herrero-Langreo, A., and Gowen, A. A. (2020).  Deep learning classifiers for near infrared spectral imaging: a tutorial.  Journal of Spectral Imaging, 9.
  * Zhai et al., (2021) Zhai, H., Zhang, H., Li, P., and Zhang, L. (2021).  Hyperspectral image clustering: Current achievements and future lines.  IEEE Geoscience and Remote Sensing Magazine, 9(4):35–67.

##  8 Appendix A

As discussed in Section 2, Table 1 details the dimensionality of the nine
hyperspectral images of the three puffed cereals.

Table 1: The dimensionality of motivating hyperspectral images of puffed cereals. Puffed cereal type |  $l$ |  $S_{l}$ |  $T_{l}$ |  $S_{l}\times T_{l}$
---|---|---|---|---
Wheat |  1 |  67 |  53 |  3551
Wheat |  2 |  62 |  64 |  3968
Wheat |  3 |  65 |  55 |  3575
Corn |  4 |  61 |  34 |  2074
Corn |  5 |  39 |  69 |  2691
Corn |  6 |  56 |  54 |  3024
Rice |  7 |  56 |  49 |  2744
Rice |  8 |  78 |  54 |  4212
Rice |  9 |  50 |  44 |  2200

##  9 Appendix B: AECM for constrained-PGMM

In constrained-PGMM, the observed data log-likelihood is

$\ell_{o}(\mathbf{X})=\sum_{n=1}^{N}\log\sum_{g=1}^{G}\tau_{g}\>f(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g}).$
(4)

As (4) is difficult to maximize directly, an AECM algorithm is used to fit the
constrained-PGMM. AECM works in two cycles and allows for a different
definition of the complete data in each cycle to estimate the parameters.

###  9.1 First cycle

In the first cycle, the component indicators $\mathbf{Z}$ =
{$\mathbf{z}_{1}$,…, $\mathbf{z}_{N}$} are assumed to be the missing data,
where $\mathbf{z}_{N}$ is a binary vector of length $G$. The complete data
log-likelihood is

$\ell_{c}(\mathbf{X},\mathbf{Z})=\sum_{n=1}^{N}\sum_{g=1}^{G}{z}_{ng}[\log\tau_{g}+\log\>f(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g})],$
(5)

where
$f(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g})\sim\mathcal{N}_{p}(\bm{\mu}_{g},\bm{\Lambda}_{g}\>\bm{\Lambda}_{g}^{\top}\>+\>\bm{\Psi}_{g})$.

####  9.1.1 E-step:

The expected complete data log-likelihood is,

$\begin{split}Q_{1}(\bm{\mu}_{g},\tau_{g})=\sum_{n=1}^{N}\sum_{g=1}^{G}E[{z}_{ng};\mathbf{x}_{n}]\Bigg{[}\log\tau_{g}-\frac{p}{2}\log(2\pi)\\\
-\frac{1}{2}\log(|\bm{\Lambda}_{g}\>\bm{\Lambda}_{g}^{\top}\>+\>\bm{\Psi}_{g}|)\\\
-\frac{1}{2}[(\mathbf{x}_{n}-\bm{\mu}_{g})^{\top}(\bm{\Lambda}_{g}\>\bm{\Lambda}_{g}^{\top}\>+\>\bm{\Psi}_{g})^{-1}(\mathbf{x}_{n}-\bm{\mu}_{g})]\Bigg{]}.\end{split}$

The posterior probability ${\hat{z}_{ng}}$ of observation $n$ belonging to
cluster $g$ is

$\begin{split}\indent
E[{z}_{ng}\mid\mathbf{x}_{n}]=\hat{{z}}_{ng}=\Pr({z}_{ng}=1;\mathbf{x}_{n},\hat{\tau}_{g},\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g},\hat{\bm{\Psi}}_{g}),\\\
\hat{{z}}_{ng}=\frac{\hat{\tau}_{g}\>\mathcal{N}_{p}(\mathbf{x}_{n};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})}{\sum_{g^{{}^{\prime}}}^{G}\hat{\tau}_{g^{{}^{\prime}}}\>\mathcal{N}_{p}(\mathbf{x}_{n};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})}.\end{split}$

where $\hat{\cdot}$ denotes an initial or current parameter estimate, as
relevant.

Let $B_{1}$ be the set of $J$ pixels (indexed by $j$) in the blue blocks shown
in Figure 2 which are known to be background and must be clustered together.
This is a positive constraint (denoted by $+$). The posterior probability
$\hat{z}_{{B_{1}}g}^{+}$ of the pixels in $B_{1}$ belonging to cluster $g$
given the positive constraints ($+$) is

$\centering\indent\indent\hat{{z}}_{{B_{1}},g}^{+}=\frac{\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})}{\displaystyle\sum_{g^{{}^{\prime}}=1}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g^{{}^{\prime}}}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})}.\@add@centering$

Let $B_{2}$ be the set of $O$ pixels (indexed by $o$) in the yellow block of
Figure 2 which are known to be wheat and should be clustered together. These
are positive constraints (denoted $+$). However, the pixels in $B_{1}$ should
not be clustered together with the pixels in $B_{2}$. This is a negative
constraint (denoted $-$). Let $A$ be the event that pixels in $B_{1}$ belong
to cluster $g$ and $B$ be the event that pixels in $B_{2}$ belong to any
cluster $f$ in $G$ but $f\neq g$ and vice versa. Then the probability of event
$A$ and $B$ occurring together is

$\Pr(A\cap B)=\Pr(A)\Pr(B).$ (6)

Based on equation (6), the posterior probability
$\hat{{z}}_{{B_{1}},g}^{+\atop-}$ of the pixels in $B_{1}$ belonging to
cluster $g$, given the positive and negative constraints $+\atop-$ is

$\centering\indent\indent\hat{{z}}_{{B_{1}},g}^{+\atop-}=\frac{\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq g\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}{\displaystyle\sum_{g^{{}^{\prime}}=1}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g^{{}^{\prime}}}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq
g^{{}^{\prime}}\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}.\@add@centering$
(7)

Similarly, the posterior probability $\hat{{z}}_{{B_{2}},g}^{+\atop-}$ of the
pixels in $B_{2}$ belonging to cluster $g$, given the positive and negative
constraints $+\atop-$ is

$\centering\indent\indent\hat{{z}}_{{B_{2}},g}^{+\atop-}=\frac{\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{g}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq g\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}{\displaystyle\sum_{g^{{}^{\prime}}=1}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{J}\hat{\tau}_{g^{{}^{\prime}}}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq
g^{{}^{\prime}}\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}.\@add@centering$

The remaining pixels in Figure 2 which are not in the highlighted blocks are
considered as blocks with only one pixel and no negative constraints.

####  9.1.2 CM-step:

The expected complete data log-likelihood can be rewritten as

$\begin{split}Q_{1}(\bm{\mu}_{g},\tau_{g})=\sum_{n=1}^{N}\sum_{g=1}^{G}\hat{z}_{ng}^{+\atop-}\Bigg{[}\log\tau_{g}-\frac{p}{2}\log(2\pi)\\\
-\frac{1}{2}\log(|\bm{\Lambda}_{g}\>\bm{\Lambda}_{g}^{\top}\>+\>\bm{\Psi}_{g}|)\\\
-\frac{1}{2}[(\mathbf{x}_{n}-\bm{\mu}_{g})^{\top}(\bm{\Lambda}_{g}\>\bm{\Lambda}_{g}^{\top}\>+\>\bm{\Psi}_{g})^{-1}(\mathbf{x}_{n}-\bm{\mu}_{g})]\Bigg{]}.\end{split}$
(8)

Differentiating (8) with respect to $\bm{\mu}_{g}$,

$\indent\indent\hat{\bm{\mu}}_{g}=\frac{\sum_{n=1}^{N}\hat{z}_{ng}^{+\atop-}\mathbf{x}_{n}}{\sum_{n=1}^{N}\hat{z}_{ng}^{+\atop-}}.$

Differentiating (8) with respect to $\tau_{g}$,

$\indent\indent\hat{\tau}_{g}=\frac{\sum_{n=1}^{N}\hat{z}_{ng}^{+\atop-}}{N}.$

###  9.2 Second cycle

In the second cycle, both component indicators
$\mathbf{Z}=\\{\mathbf{z}_{1},\ldots,\mathbf{z}_{n}\\}$ and latent factors
$\mathbf{U}=\\{\mathbf{u}_{1},\ldots,\mathbf{u}_{n}\\}$ are assumed to be the
missing data. Therefore, the complete data log-likelihood is

$\ell_{c}(\mathbf{X},\mathbf{Z},\mathbf{U})=\sum_{n=1}^{N}\sum_{g=1}^{G}z_{ng}[\log\tau_{g}+\log
f(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g},\mathbf{u}_{n})+\log
f(\mathbf{u}_{n})],$

where
$f(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g},\mathbf{u}_{n})\sim\mathcal{N}_{p}(\bm{\mu}_{g}+\bm{\Lambda}_{g}\mathbf{u}_{n},\bm{\Psi}_{g})$
and $f(\mathbf{u}_{n})\sim\mathcal{N}_{q}(0,I)$.
The complete data log-likelihood can be rewritten as,

$\begin{split}\ell_{c}(\mathbf{X},\mathbf{Z},\mathbf{U})=C+\sum_{g=1}^{G}\Bigg{[}\mathbf{n}_{g}\log\tau_{g}-\frac{\mathbf{n}_{g}}{2}\log(|\bm{\Psi}_{g}|)-\frac{1}{2}\>\Bigg{[}\sum_{n=1}^{N}z_{ng}[(\mathbf{x}_{n}-\bm{\mu}_{g})^{\top}\bm{\Psi}_{g}^{-1}(\mathbf{x}_{n}-\bm{\mu}_{g})]\Bigg{]}\\\
+\sum_{n=1}^{N}z_{ng}[(\mathbf{x}_{n}-\bm{\mu}_{g})^{\top}\bm{\Psi}_{g}^{-1}(\bm{\Lambda}_{g}\mathbf{u}_{n})]\\\
-\frac{1}{2}(\bm{\Lambda}_{g}^{\top}\bm{\Psi}_{g}^{-1}\bm{\Lambda}_{g}\sum_{n=1}^{N}z_{ng}\mathbf{u}_{n}\mathbf{u}_{n}^{\top})^{\top}\Bigg{]},\end{split}$

where $C$ is the constant and $\mathbf{n}_{g}=\sum_{n=1}^{N}z_{ng}$.

####  9.2.1 E-step:

The expected complete data log-likelihood is,

$\begin{split}Q_{2}(\bm{\Lambda}_{g},\bm{\Psi}_{g})=C+\sum_{g=1}^{G}\Bigg{[}\mathbf{n}_{g}\log\tau_{g}-\frac{\mathbf{n}_{g}}{2}\log(|\bm{\Psi}_{g}|)-\frac{\mathbf{n}_{g}}{2}(\bm{\Psi}_{g}^{-1}\mathbf{S}_{g})^{\top}\\\
+\sum_{n=1}^{N}E[z_{ng};\mathbf{x}_{n}][(\mathbf{x}_{n}-\bm{\mu}_{g})^{\top}\bm{\Psi}_{g}^{-1}(\bm{\Lambda}_{g}E[\mathbf{u}_{n};\mathbf{x}_{n}])]\\\
-\frac{1}{2}(\bm{\Lambda}_{g}^{\top}\bm{\Psi}_{g}^{-1}\bm{\Lambda}_{g}\sum_{n=1}^{N}E[z_{ng};\mathbf{x}_{n}]E[\mathbf{u}_{n}\mathbf{u}_{n}^{\top};\mathbf{x}_{n}])^{\top}\Bigg{]},\end{split}$
(9)

where
$\mathbf{S}_{g}=\frac{1}{\mathbf{n}_{g}}\sum_{n=1}^{N}{z}_{ng}(\mathbf{x}_{n}-{\bm{\mu}}_{g})(\mathbf{x}_{n}-{\bm{\mu}}_{g})^{\top}$.

The expected value $E[z_{ng}\mid\mathbf{x}_{n}]=\hat{z}_{ng}$ is computed with
$\bm{\mu}_{g}=\hat{\bm{\mu}}_{g}$ and $\tau_{g}=\hat{\tau}_{g}$ updated in the
CM-step of the first cycle with the known positive and negative constraints.
The posterior probability $\hat{{z}}_{{B_{1}},g}^{+\atop-}$ of the pixels in
$B_{1}$ belonging to cluster $g$, given the positive and negative constraints
$+\atop-$ is

$\centering\indent\indent\hat{{z}}_{{B_{1}},g}^{+\atop-}=\frac{\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g},\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq g\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}{\displaystyle\sum_{g^{{}^{\prime}}=1}^{G}\prod\limits_{\begin{subarray}{c}j\in{B_{1}}\\\
j=1\end{subarray}}^{J}\hat{\tau}_{g^{{}^{\prime}}}\mathcal{N}_{p}(\mathbf{x}_{j};\hat{\bm{\mu}}_{g^{{}^{\prime}}},\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}\>\hat{\bm{\Lambda}}_{g^{{}^{\prime}}}^{\top}\>+\>\hat{\bm{\Psi}}_{g^{{}^{\prime}}})\sum\limits_{\begin{subarray}{c}f=1\\\
f\neq
g^{{}^{\prime}}\end{subarray}}^{G}\prod\limits_{\begin{subarray}{c}o\in{B_{2}}\\\
o=1\end{subarray}}^{O}\hat{\tau}_{f}\mathcal{N}_{p}(\mathbf{x}_{o};\hat{\bm{\mu}}_{f},\hat{\bm{\Lambda}}_{f}\>\hat{\bm{\Lambda}}_{f}^{\top}\>+\>\hat{\bm{\Psi}}_{f})}.\@add@centering$
(10)

Similarly, the computation is repeated for estimating the posterior
probabilities of pixels in $B_{2}$ belonging to cluster $g$ with pixels in
$B_{1}$ as negative constraints. The expected value
$E[\mathbf{u}_{n}\mid\mathbf{x}_{n}]$ is computed by considering a joint
Gaussian distribution between $\mathbf{x}_{n}$ and $\mathbf{u}_{n}$. By
Gaussian conditioning formulas,

$E[\mathbf{u}_{n}\mid\mathbf{x}_{n}]=\hat{\bm{\Lambda}}_{g}^{\top}(\hat{\bm{\Lambda}}_{g}\>\hat{\bm{\Lambda}}_{g}^{\top}\>+\>\hat{\bm{\Psi}}_{g})^{-1}(\mathbf{x}_{n}-\hat{\bm{\mu}}_{g}).$

Also,

$E[\mathbf{u}_{n}\mathbf{u}_{n}^{\top}\mid\mathbf{x}_{n}]=\mathbf{I}_{q}-\hat{\bm{\beta}}_{g}\hat{\bm{\Lambda}}_{g}+\hat{\bm{\beta}}_{g}(\mathbf{x}_{n}-\hat{\bm{\mu}}_{g})(\mathbf{x}_{n}-\hat{\bm{\mu}}_{g})^{\top}\hat{\bm{\beta}}_{g}^{\top},$

where
$\hat{\bm{\beta}}_{g}=\hat{\bm{\Lambda}}_{g}^{\top}(\hat{\bm{\Lambda}}_{g}\hat{\bm{\Lambda}}_{g}^{\top}+\hat{\bm{\Psi}}_{g})^{-1}$.

####  9.2.2 CM-step:

The expected complete log-likelihood can be rewritten as

$\begin{split}Q_{2}({\bm{\Lambda}_{g}},{\bm{\Psi}_{g}})=C+\sum_{g=1}^{G}\mathbf{n}_{g}\Bigg{[}-\frac{1}{2}\log|\bm{\Psi}_{g}|-\frac{1}{2}[\bm{\Psi}_{g}^{-1}\mathbf{S}_{g}]^{\top}+[\bm{\Psi}_{g}^{-1}\bm{\Lambda}_{g}\hat{\bm{\beta}}_{g}\bm{S}_{g}]^{\top}\\\
-\frac{1}{2}[\bm{\Lambda}_{g}^{\top}\bm{\Psi}_{g}^{-1}\bm{\Lambda}_{g}\bm{\Theta}_{g}]\Bigg{]},\end{split}$
(11)

where,

$\centering\begin{split}\bm{\Theta}_{g}=\mathbf{I}_{q}-\hat{\bm{\beta}}_{g}{\bm{\Lambda}}_{g}+\hat{\bm{\beta}}_{g}\mathbf{S}_{g}\hat{\bm{\beta}}_{g}^{\top},\\\
\mathbf{S}_{g}=\frac{1}{\mathbf{n}_{g}}\sum_{n=1}^{N}\hat{z}_{ng}^{+\atop-}(\mathbf{x}_{n}-\hat{\bm{\mu}}_{g})(\mathbf{x}_{n}-\hat{\bm{\mu}}_{g})^{\top},\\\
\mathbf{n}_{g}=\sum_{n=1}^{N}\hat{z}_{ng}^{+\atop-}.\end{split}\@add@centering$

Differentiating 11 with respect to $\bm{\Lambda}_{g}$

$\hat{\bm{\Lambda}}_{g}=\mathbf{S}_{g}\hat{\bm{\beta}}_{g}^{\top}\bm{\Theta}_{g}^{-1}.$

Differentiating 11 with respect to $\bm{\Psi}_{g}^{-1}$

$\hat{\bm{\Psi}}_{g}=diag[\mathbf{S}_{g}-\hat{\bm{\Lambda}}_{g}\hat{\bm{\beta}}_{g}\mathbf{S}_{g}].$

The algorithm iteratively updates the parameters until convergence and the
posterior probabilities ${z}_{ng}$ at convergence are used to compute the
similarity matrix in the consensus approach.

##  10 Appendix C

As discussed in Section 4, Figures 13 and 14 highlight the $43.5\%$ of pixels
and $24.7\%$ of pixels used as constraints respectively.

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 13: Blocks of pixels ($12,215$ pixels of $N=28,039$ pixels ($43.5\%$))
selected as constraints in three greyscale images for each of the three types
of synthetic puffed cereal, wheat (13(a), 13(b), 13(c)), corn (13(d), 13(e),
13(f)), and rice (13(g), 13(h), 13(i)). The pixels in the blue blocks are
background, the pixels in the yellow blocks are wheat, the pixels in the green
blocks are corn and the pixels in the red blocks are rice.

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 14: Blocks of pixels ($6,951$ pixels of $N=28,039$ pixels ($24.7\%$))
selected as constraints in three greyscale images for each of the three types
of synthetic puffed cereal, wheat (14(a), 14(b), 14(c)), corn (14(d), 14(e),
14(f)), and rice (14(g), 14(h), 14(i)). The pixels in the blue blocks are
background, the pixels in the yellow blocks are wheat, the pixels in the green
blocks are corn and the pixels in the red blocks are rice.

##  11 Appendix D

As discussed in Section 4.4, Table 2 details the cumulative proportion of the
variance explained by the first $5$ principal components of the motivating
hyperspectral image dataset.

Table 2: Cumulative proportion of variance explained by the principal components of the motivating hyperspectral image dataset on puffed cereals. Number of principal components |  1 |  2 |  3 |  4 |  5
---|---|---|---|---|---
Cumulative proportion of variance explained |  94.2 |  99.7 |  99.92 |  99.94 |  99.95

##  12 Appendix E

As discussed in Section 5, Figures 15, and 16 shows the greyscale images for
the remaining six hyperspectral images, and the respective pixel labels as
generated using different clustering approaches.

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

Figure 15: Greyscale images and pixel labels for three of the nine
hyperspectral images ($l=\\{2,5,8\\}$) of the cereals using the threshold
approach, DBSCAN, MBC, PGMM, ccPGMM with $43.5\%$ of pixels as constraints,
and ccPGMM with $24.7\%$ of pixels as constraints respectively for wheat
(15(a), 15(d), 15(g), 15(j), 15(m), 15(p),15(s)), corn (15(b), 15(e), 15(h),
15(k), 15(n), 15(q), 15(t)) and rice (15(c), 15(f), 15(i), 15(l), 15(o),
15(r), 15(u)).

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

Figure 16: Greyscale images and pixel labels for three of the nine
hyperspectral images ($l=\\{3,6,9\\}$) of the cereals using the threshold
approach, DBSCAN, MBC, PGMM, ccPGMM with $43.5\%$ pixels as constraints, and
ccPGMM with $24.7\%$ pixels as constraints respectively for wheat (16(a),
16(d), 16(g), 16(j), 16(m), 16(p), 16(s)), corn (16(b), 16(e), 16(h), 16(k),
16(n), 16(q), 16(t)) and rice (16(c), 16(f), 16(i), 16(l), 16(o), 16(r),
16(u)).

##  13 Appendix F

As discussed in Section 5, Figures 17, 18, and 19 illustrates the uncertainty
associated with the pixel labels for the nine hyperspectral images of the
puffed cereals based on the cluster solutions of MBC, PGMM, ccPGMM with
$43.5\%$ pixels as constraints and ccPGMM with $24.7\%$ pixels as constraints.

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

Figure 17: Uncertainty associated with pixel labels for three of the nine
hyperspectral images ($l=\\{1,4,7\\}$) of the puffed cereals based on the
cluster solutions of MBC, PGMM, ccPGMM with $43.5\%$ of pixels as constraints,
and ccPGMM with $24.7\%$ of pixels as constraints respectively for wheat
(17(a), 17(e), 17(i),17(m)), corn (17(b), 17(f), 17(j), 17(n)) and rice
(17(c), 17(g), 17(k), 17(o)).

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

Figure 18: Uncertainty associated with pixel labels for three of the nine
hyperspectral images ($l=\\{2,5,8\\}$) of the cereals based on the cluster
solutions of MBC, PGMM, ccPGMM with $43.5\%$ of pixels as constraints, and
ccPGMM with $24.7\%$ of pixels as constraints respectively for wheat (18(a),
18(e), 18(i),18(m)), corn (18(b), 18(f), 18(j), 18(n)) and rice (18(c), 18(g),
18(k), 18(o)).

---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

Figure 19: Uncertainty associated with pixel labels for three of the nine
hyperspectral images ($l=\\{3,6,9\\}$) of the cereals based on the cluster
solutions of MBC, PGMM, ccPGMM with $43.5\%$ of pixels as constraints, and
ccPGMM with $24.7\%$ of pixels as constraints respectively for wheat (19(a),
19(e), 19(i),19(m)), corn (19(b), 19(f), 19(j), 19(n)) and rice (19(c), 19(g),
19(k), 19(o)).