# Automated Characterization and Monitoring of Material Shape using Riemannian
Geometry

Alexander Smith${}^{\text{a}}$, Steven Schilling${}^{\text{b}}$, Prodromos
Daoutidis${}^{\text{a}}$
${}^{\text{a}}$Department of Chemical Engineering and Materials Science
University of Minnesota, Minneapolis, MN, USA
${}^{\text{b}}$Technology and Application Development
Covia Corp, Independence, OH, USA
Corresponding Author: 421 Washington Ave SE #250, Minneapolis, MN 55455
<EMAIL_ADDRESS>
###### Abstract

Shape affects the physical and chemical properties of a material.
Characterizing the roughness, convexity, and general geometry of a material
can yield information on its catalytic efficiency, solubility, porosity, and
overall effectiveness in the application of interest. However, material shape
can be defined in conflicting ways where certain aspects of a material’s
geometry are emphasized over others, leading to measures of shape that are not
generalizable. We explore the use of a Riemannian geometric framework for
shape analysis that is generalizable, computationally scalable, and can be
integrated into common data analysis methods. Here, material shapes are
abstracted as points on a Riemannian manifold. This information can be used to
construct statistical moments (e.g., means, variances) and perform tasks such
as dimensionality reduction and statistical process control. We provide an
introduction to the mathematics of shape analysis and illustrate its
application on a manufactured/mined granular material dataset provided by
Covia Corp.

Keywords: Data science, Computer Vision, Statistical Process Control,
Dimensionality reduction, Shape Space Analysis

##  1 Introduction

The geometric shape of a manufactured material impacts both its chemical and
physical properties. This is evident in many areas of manufacturing such as
3-dimensional printing, crystal and co-crystal growth, cell morphology in
biomanufacturing, and in the production of granular materials such as sand and
ceramic microspheres [20, 42, 1, 35]. For example, granular materials like
sand and ceramic microspheres can be added into liquids such as paints,
coatings, and cements to improve hardness, decrease viscosity, increase
insulating capability, and reduce the overall amount of volatile organic
compounds (VOC’s) [52, 3]. However, these improvements are directly dependent
on the morphology of the individual grains, with spherical grains providing
the largest improvement and irregularly shaped particles resulting in
decreased performance due to how the particles coalesce at a larger scale [19,
2]. This also holds true for granular materials used in metal casting and
hydraulic fracturing where the flow of organic materials through the granular
medium needs to be closely controlled to prevent metal casting defects and
optimal extraction of natural gas and oil from deep rock formations [56, 27].
Thus, there is a need for methods to automatically quantify and monitor the
morphological characteristics of manufactured materials. Unfortunately, many
of the geometric methods employed in manufacturing are application specific.
An example is the use of rigid geometric structures for the analysis of
crystal morphology in the production of pharmaceuticals. Crystal structure
(i.e., crystal form) of an active pharmaceutical ingredient (API) impacts its
density, solubility, reactivity, and stability, among other properties [6].
The structure of crystals allows their shape to be quantified with measures
such as aspect ratio, form factor, and roundedness, or complex transforms such
as wavelet and Fourier transforms [44, 11, 4]. However, these methods are
difficult to apply to amorphous structures such as those found in cells used
in the production of bio-pharmaceuticals, in crystal polymorphisms, or in the
mining and refinement of natural materials such as sand [51, 30]. In these
cases, the shapes of the materials are often subjectively compared making it
difficult to directly quantify differences in their morphology. In the
analysis of natural granular material (e.g., sand), a visual chart known as
the Krumbein-Sloss chart is used to quantify shape [25]. The Krumbein-Sloss
chart is used in accordance with the EN:ISO 13503-2:2006 standard, in which 20
grains of sand are randomly chosen and manually compared to illustrations of
sand grains to estimate the roundness and sphericity of the particles [14].
Such methods are time-consuming and subjective to the interpretation of the
individual performing the analysis.

Machine learning (ML) methods, such as convolutional neural networks, have
recently been proposed as generalizable measures of morphology for materials
[53, 17, 8, 24, 46]. However, these methods require large amounts of well-
sampled training data. In the context of manufacturing and statistical quality
control this means that there must be a large number of defective/outlier
samples in the training data [54]. These data are sparse by nature and may
require significant interventions into operating process systems to obtain a
diverse enough sample of the various faults that could occur, resulting in
significant process downtime and lost revenue. Furthermore, ML methods are
difficult to physically interpret and do not provide a direct avenue for
robust statistical analysis of their output.

To address these challenges, we propose the use of a mathematical framework to
automatically characterize the morphology of manufactured materials using
Riemannian geometry. This framework, often referred to as _statistical shape
analysis_ , was first introduced by Kendall and defines shape as the geometric
information remaining after rotation, translation, and scaling have been
removed from an object [21, 23, 5]. It is based on the realization that
geometric shapes can be represented as points on a Riemannian manifold [13].
The structure of this Riemannian manifold can be used to directly quantify
differences between shapes based on geodesic distances on the manifold. These
geodesic distances can be used to develop statistical measures (e.g., means,
variances) of a material’s intrinsic morphology which can be used in process
monitoring and quality control. Furthermore, the Riemannian manifold structure
allows us to project data from the manifold onto a tangent (vector) space
which can be directly integrated in data analysis tasks such as dimensionality
reduction and classification [45]. This framework has been successfully
applied in many research areas such as evolutionary biology [33, 36],
neuroscience [38, 29], chemistry and biochemistry [43, 12, 9], cell morphology
[31, 34, 16], and design of experiments [7], but has found little application
in the areas of process monitoring and statistical process control.

In this work, we present the mathematical foundations of shape analysis
through a Riemannian geometric framework and illustrate its application on a
manufactured/mined granular material dataset provided by Covia Corp. We focus
on samples of sand and manufactured ceramic microspheres. We analyze
microscope images of these samples and develop an automated method for
extracting particle shapes from the images. We leverage the presented
Riemannian framework to perform dimensionality reduction to visualize the
structure of the data and hypothesis testing to understand the morphological
differences in the samples. The automated, computationally efficient nature of
this framework, coupled with its statistical power, suggest a powerful
statistical process control technique for the continuous improvement and
quality control of processes in which shape is a key factor. All code and data
needed to reproduce the results presented here can be found at the following
link: https://github.com/AdSmithPhd/Shape_Analysis.

##  2 Geometric Shape

We first cover the basic mathematics of the analysis of geometric shape.
Geometric shape is often defined as the geometric information that remains
when rotation, translation, and scaling are removed from a data object [13,
21]. The removal of translation, rotation, and scaling allows us to directly
compare the differences in the intrinsic morphology of different objects. For
example, Figure 1 illustrates an analysis of two triangles in which we
translate, scale, and rotate the triangles to best match the other. The
remaining geometric differences shown in Figure 1(d) represent the difference
in _shape_ of the triangles.

(a) Initialization

(b) Translation

(c) Scaling

(d) Rotation

Figure 1: Two triangles are translated, scaled, and rotated to emphasize the
difference in shape between the two triangles. Once these transformations are
completed, the triangles can be analyzed via the Riemannian geometric shape
space analysis framework.

We will cover the operations of translation, scaling, and rotation
individually based on the triangles shown in Figure 1 and then extend the
methods to general shapes. A key aspect of the analysis of these triangles is
that they can be represented by _landmark points_. A landmark point is defined
formally as a point of correspondence on each shape that matches between and
within populations [13]. For the triangles in our example the landmark points
are given by the location of the corners of each triangle. Thus, each triangle
can be represented by a set of $k=3$ points, ${x_{1},x_{2},x_{3}}$, where
$x_{i}\in\mathbb{R}^{2}$; this is known as the _configuration_. From this
configuration we can construct a $3\times 2$ matrix
$\mathbf{X}=[x_{1},x_{2},x_{3}]$ that represents the 3 landmark points in 2
dimensions. With this definition we can now understand the mathematics of
translation and leverage a centering matrix $\mathbf{C}$ to center our
individual triangles around the origin $(0,0)$. The centering matrix is often
applied in multivariate statistics and mean centers components of a vector
[28]. It is defined as:

$\displaystyle\mathbf{C}:=I_{k}-\frac{1}{k}J_{k}$ (2.1)

where $I_{k}$ is the identify matrix of size $k$ and $J_{k}$ is a $k\times k$
matrix of ones. For example, in our triangular case with three points,
$\mathbf{C}_{3}$ is given as:

$\displaystyle\mathbf{C}_{3}:=\left[{\begin{array}[]{rrr}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}}\right]-{\frac{1}{3}}\left[{\begin{array}[]{rrr}1&1&1\\\
1&1&1\\\
1&1&1\end{array}}\right]=\left[{\begin{array}[]{rrr}{\frac{2}{3}}&-{\frac{1}{3}}&-{\frac{1}{3}}\\\
-{\frac{1}{3}}&{\frac{2}{3}}&-{\frac{1}{3}}\\\
-{\frac{1}{3}}&-{\frac{1}{3}}&{\frac{2}{3}}\end{array}}\right]$ (2.11)

We apply $\mathbf{C}$ to the configuration matrices
$\mathbf{X}_{C}:=\mathbf{C}_{3}\mathbf{X}$ of both of our triangles in Figure
5(c). The next step in the shape analysis is to scale our shapes. A direct way
to scale a shape is to normalize by the _size_ of the shape. Size is a measure
that is formally defined as any positive real-valued function of a
configuration matrix $f(\mathbf{X})\rightarrow\mathbb{R}$ such that
$f(a\mathbf{X})=af(\mathbf{X})$ for any $a>0$ and $a\in\mathbb{R}$ [13]. An
often used measure of size is the centroid size, denoted as
$S(\mathbf{X})\rightarrow\mathbb{R}$. The centroid size is defined as [13]:

$\displaystyle S(\mathbf{X}):=||\mathbf{C}\mathbf{X}||$ (2.12)

where $||\cdot||$ represents the Frobenius (i.e. Euclidean) norm. We can then
scale our centered configurations using the centroid size. We will denote
these scaled and translated configurations as $\mathbf{Z}$:

$\displaystyle\mathbf{Z}:=\frac{\mathbf{X}_{C}}{||\mathbf{X}_{C}||}$ (2.13)

The scaled and translated versions of our triangles are found in Figure 5(c).
The final operation of rotation is a slightly more difficult operation, and
requires an understanding of _pre-shape space_.

###  2.1 Pre-Shape Space and Riemannian Geometry

Pre-shape space is an idea first introduced by Kendall and represents the
space of all k-point sets (that are not perfectly overlapping) in
$\mathbb{R}^{m}$ that have been properly translated and scaled [22]. This pre-
shape space is often denoted as $S_{m}^{k}$. For our triangle example in
Figure 1, the pre-shape space represents the space of all possible triangular
configurations of points and their rotations [13]. A key notion in the
analysis of shapes is that pre-shape space is in fact a hypersphere of unit
radius in $(k-1)m-1$ dimensions [22]. Thus, shapes are represented as single
points on the surface of this hypersphere. An illustration of the
$3$-dimensional pre-shape space sphere for all planar triangles is found in
Figure 2.

The pre-shape hypersphere represents a _Riemannian manifold_ on which we can
measure distance and angles. In general, a manifold $\mathcal{M}$ is a
topological space that _locally_ resembles a Euclidean space; this means that
the neighborhood of a point $p\in\mathcal{U}$ in an $n$-dimensional manifold
(with $\mathcal{U}\subseteq\mathcal{M}$) can be mapped to $n$-dimensional
Euclidean space through a continuous, bijective function. These neighborhoods
and associated mappings are also known as _charts_. For example, the surface
of the Earth is a 2D manifold and we can map the curved surface of the Earth
to a flat Euclidean plane (i.e., a 2D Euclidean space) using a chart in order
to measure properties such as distances or areas. A _Riemannian_ manifold
represents a differentiable manifold $\mathcal{M}$ on which there exists a
defined tangent space $T_{p}M$ at each point $p\in\mathcal{M}$ along with a
defined inner product $g_{p}:T_{p}M\times T_{p}M\rightarrow\mathbb{R}$. This
inner product allows us to compute distances and angles between points on the
manifold with respect to the geometry of the Riemannian manifold. A rigorous
introduction to the mathematics of Riemannian manifolds can be found in our
previous work and other resources [45, 37, 18].

Pre-shape space represents a simple Riemannian manifold, the unit sphere [22,
13]. For example, the distance between two points on a unit sphere is exactly
equal to the angle between the points. Thus, for our analysis we can compute
the distance between two points on the pre-shape space sphere through an inner
product. For two points (e.g., shapes) $\mathbf{Z}_{i},\mathbf{Z}_{j}\in
S_{m}^{k}$ we can compute their Frobenius inner product as:

$\displaystyle\langle\mathbf{Z}_{i},\mathbf{Z}_{j}\rangle:=\text{tr}(\mathbf{Z}_{i}^{T}\mathbf{Z}_{j})=\sum_{i=1}^{m}\lambda_{i}$
(2.14)

where $\lambda_{i}$ represent the eigenvalues of the matrix
$\mathbf{Z}_{i}^{T}\mathbf{Z}_{j}$. We can convert this inner product into an
angle, and thus a distance on the pre-shape space. We denote this distance as
a _geodesic_ distance $\rho(\mathbf{Z}_{i},\mathbf{Z}_{j})$ defined as:

$\displaystyle\rho(\mathbf{Z}_{i},\mathbf{Z}_{j}):=\text{arccos}(\sum_{i=1}^{m}\lambda_{i})$
(2.15)

An illustration of the geodesic distance between two points on the pre-shape
space sphere is found in Figure 2 represented by the red, dashed line and
blue, dash-dotted lines. The Riemannian geometry of the pre-shape space
provides a concise mathematical framework for the analysis of shape. However,
the pre-shape space does not account for rotation of shapes. Thus, our
previously defined geodesic distances on the pre-shape space sphere may
compute a difference between two shapes that reflects a difference in rotation
and shape rather than shape alone. This is illustrated in Figure 2, where the
blue, dash-dotted lines represent geodesics on the pre-shape space between
triangles that are not optimally rotated. The red, dashed line represents the
smallest geodesic distance between the two shapes, capturing the true
difference in shape of the triangles with optimal rotation.

Figure 2: Illustration of the pre-shape Riemannian manifold for triangles. We
illustrate geodesics between the solid triangle and three different rotations
of the dashed triangle. The dash-dotted, blue geodesics represent geodesics
between triangles that have not been optimally rotated. The dashed, red
geodesic represents the shortest geodesic between the two optimally aligned
triangles. This geodesic quantifies the true difference in shape between the
triangles.

Thus, there is a need to optimize the rotation of two shapes such that there
is optimal alignment on the pre-shape space and we are able to compute a
geodesic between the shapes that reflects changes in shape rather than
rotation. We can frame this as an optimization problem and use this framework
to identify an optimal _rotation matrix_ for one of the shapes. A rotation
matrix is formalized as an $m\times m$ matrix denoted as $\Gamma$ where
$\Gamma^{T}\Gamma=\Gamma\Gamma^{T}=I_{m}$ and $\text{det}(\Gamma)=+1$. $I_{m}$
is the identity matrix of size $m$ and $\text{det}(\cdot)$ represents the
matrix determinant [13]. The set of all $m\times m$ rotation matrices is
formally known as the special orthogonal group denoted SO(m). A simple method
for computing an optimal rotation (alignment) between two shapes
$\mathbf{X}_{i},\mathbf{X}_{j}$ is through the following formulation, which is
known as the _partial Procrustes distance_ denoted as $d_{p}(\cdot)$ [41]:

$\displaystyle d_{p}(\mathbf{X_{i}},\mathbf{X_{j}}):=\inf_{\mathbf{\Gamma}\in
SO(m)}\ ||\mathbf{Z}_{i}-\mathbf{Z}_{j}\Gamma||$ (2.16)

where $||\cdot||$ represents the Frobenius norm and
$\mathbf{Z}_{i},\mathbf{Z}_{j}\in S_{m}^{k}$ represent the shapes
$\mathbf{X_{i}},\mathbf{X_{j}}$ mapped to the pre-shape space. This can be
rewritten analytically as:

$\displaystyle
d_{p}(\mathbf{X_{1}},\mathbf{X_{2}})=\sqrt{2}\left(1-\sum_{i=1}^{m}\lambda_{i}\right)^{1/2}$
(2.17)

where $\lambda_{i}$ represent the square roots of the eigenvalues of
$\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}\mathbf{Z}_{2}^{T}\mathbf{Z}_{1}$. The
optimal rotation is given as $\hat{\Gamma}$:

$\displaystyle\hat{\Gamma}:=\mathbf{U}\mathbf{V}^{T}$ (2.18)

where $\mathbf{U}$ and $\mathbf{V}$ are from the singular value decomposition
of $\mathbf{Z}_{2}^{T}\mathbf{Z}_{1}$:

$\displaystyle\mathbf{Z}_{2}^{T}\mathbf{Z}_{1}=\mathbf{V}\mathbf{\Lambda}\mathbf{U}^{T}$
(2.19)

Using this information, we can formulate the geodesic distance between two
shapes that have been optimally rotated. We denote this distance as
$\hat{\rho}(\mathbf{X}_{i},\mathbf{X}_{j})$ and define it as:

$\displaystyle\hat{\rho}(\mathbf{X}_{i},\mathbf{X}_{j}):=2\
\text{arcsin}(d_{p}(\mathbf{X_{1}},\mathbf{X_{2}})/2)$ (2.20)

The optimal geodesic distance now represents a measure of the true geometric
difference between two shapes. Formally, this distance now represents a proper
geodesic distance on the _shape space_. The shape space represents the space
of all k-point sets (that are not perfectly overlapping) in $\mathbb{R}^{m}$
that have been properly translated, scaled, _and rotated_ [13, 22]. The shape
space contains all information about a given configuration matrix $\mathbf{X}$
that is invariant under translation, rotation and scaling. This is often
denoted by the set:

$\displaystyle[\mathbf{X}]:=\\{\mathbf{Z}\Gamma:\Gamma\in SO(m)\\}$ (2.21)

This means that a given pre-shape $\mathbf{Z}\in S_{m}^{k}$ and all its
possible rotations $\\{\mathbf{Z}\Gamma:\Gamma\in SO(m)\\}$ are mapped to
single point in shape space. Thus, we are able to now leverage the proper
geodesics in the shape space to compare the geometric shape difference between
objects. However, one limiting factor in our analysis is that the computed
geodesic distances represent non-linear relationships between different
shapes. This makes it difficult to directly apply many data-centric techniques
that assume that data lies in a linear (vector) space. One way to address this
issue is to leverage the tangent space of the pre-shape space manifold, which
is a linear (vector) space, and project our data to this linear space which
retains much of the geometric information of the pre-shape manifold.

###  2.2 Tangent Space Analysis

For a manifold $\mathcal{M}$, such as the pre-shape space, the _tangent space_
$T_{p}M$ at a point $p\in\mathcal{M}$ is formally defined as the space of
derivates of all _curves_ passing through $p$ [45]. A curve is defined as a
continuous function $\gamma:[a,b]\rightarrow\mathcal{M}$ which maps the
interval $[a,b]\in\mathbb{R}$ to the manifold $\mathcal{M}$. Thus, the tangent
space $T_{p}M$ is built from all vectors $v=\dot{\gamma(0)}\in T_{p}M$ where
$\gamma(0)=p\in\mathcal{M}$. The tangent space at a point $p\in\mathcal{M}$
represents a linear vector space that is of the same dimension as the manifold
from which it is constructed. A visualization of the tangent space for the
triangle pre-shape sphere is found in Figure 3.

(a) Tangent space at point $\mathbf{Z}_{i}\in S_{m}^{k}$

(b) Tangent space projection

Figure 3: (a) Illustration of the tangent space to a point $\mathbf{Z}_{i}\in
S_{m}^{k}$ with the tangent vector $v_{i,j}$ representing the projection of
the geodesic between $\mathbf{Z}_{i},\mathbf{Z}_{i}\in S_{m}^{k}$ to the
tangent space. (b) Illustration of how a point $\mathbf{Z}_{j}\in S_{m}^{k}$
can be mapped to the horizontal tangent space centered at $\mathbf{Z}_{i}\in
S_{m}^{k}$. To project to the horizontal component of the tangent space,
rather than the full tangent space, we must apply optimal rotation matrix
$\hat{\Gamma}$ to $\mathbf{Z}_{j}$ prior to the tangent space projection.

Another key aspect in the analysis of tangent spaces is the mapping from the
tangent space to the manifold, known as the _exponential map_ and inversely,
the _logarithmic map_. For a tangent space defined at a point
$p\in\mathcal{M}$ we can construct a curve
$\gamma:[0,1]\rightarrow\mathcal{M}$ such that $\gamma(0)=p$,
$\dot{\gamma}(0)=v\in T_{p}M$ and $\gamma(1)=q\in\mathcal{M}$. We define the
exponential map as the mapping from the vector $v\in T_{p}M$ to the endpoint
of the curve $q\in\mathcal{M}$:

$\displaystyle\text{exp}_{p}(v)=q$ (2.22)

We define the inverse of the exponential map as the logarithmic map, which
maps the point $q\in\mathcal{M}$ to the vector $v\in T_{p}M$:

$\displaystyle\text{log}_{p}(q)=v$ (2.23)

These operations can be defined for any constructed tangent space on the
surface of a Riemannian manifold and are illustrated in Figure 3. For our pre-
shape space sphere we will focus on an analytic formulation for the
logarithmic mapping as we are interested in mapping our data from the non-
linear surface of the sphere to the linear tangent space. There are many
choices for construction of a tangent space for the pre-shape sphere, but we
will focus on the _Procrustes tangent space_ as it is interpretable and has
proven effective in many application areas [13]. For two shapes on the pre-
shape manifold $\mathbf{Z}_{i},\mathbf{Z}_{j}\in S_{m}^{k}$ we can construct a
tangent space $T_{Z_{i}}M$ at $\mathbf{Z}_{i}$ and project $\mathbf{Z}_{j}$ to
the Procustes tangent space as:

$\displaystyle\hat{v}_{i,j}:=\mathbf{Z}_{j}\hat{\Gamma}-\text{cos}(\rho_{i,j})\mathbf{Z}_{i}$
(2.24)

where $\rho_{i,j}$ represents the geodesic distance between $\mathbf{Z}_{i}$
and $\mathbf{Z}_{j}$. The Procrustes tangent space differs from the full
tangent space at the point $\mathbf{Z}_{i}\in S_{m}^{k}$ as it accounts for
the rotation of the pre-shape $\mathbf{Z}_{j}\in S_{m}^{k}$ by applying the
optimal rotation matrix $\hat{\Gamma}$ derived from Procrustes analysis. This
distinction is often denoted by a decomposition of the tangent space of
$S_{m}^{k}$ into a _horizontal_ component that captures changes in shape and a
_vertical_ component that accounts for changes in rotation [41, 13]. Thus, we
are focused on the horizontal component of the tangent space as we are
interested in differences with respect to shape rather than rotation, which is
illustrated in Figure 3(b). Through this projection, we are able to project
the non-linear structure of the pre-shape manifold into a linear (vector)
space. This linear space can then be used to construct hypothesis tests,
execute statistical process control methods, and perform common data centric
tasks such as dimensionality reduction and regression/classification.

##  3 Application: Granular Material Morphology

We now focus on an application of the geometry of shape space in the analysis
of granular materials such as sand and ceramic microspheres. Granular material
particle geometry influences the function and bulk physical properties of many
manufactured materials. This holds true in many areas of manufacturing such as
pharmaceuticals, paints, coatings, metal casting, 3-dimensional printing, and
hydraulic fracturing, among others [52, 19, 6, 51]. A significant challenge in
the analysis of particle geometry is the characterization of particles that
are amorphous in nature and difficult to quantify using simple geometric
descriptors such as aspect ratios or sphericity. This is especially
problematic in the manufacture/mining of sand and production of ceramic
microspheres where geometry can fluctuate in ways that are not easily modeled.
Figure 4 shows example images from Covia corp of sand mined in different
geographic locations that exhibit angular and round geometry and ceramic
microsphere samples.

(a) Angular Sand

(b) Round Sand

(c) Ceramic Microspheres

Figure 4: Examples of sands mined at different geographic locations (a-b) and
manufactured ceramic microspheres (c).

Figure 4 illustrates the multitude of ways in which the geometry of particles
can fluctuate, demonstrating the difficulty in producing a framework for
quantifying the quality and variability of particle morphology. Currently, the
main method for characterizing the morphology of particle samples is through
manual comparison with the Krumbien-Sloss chart, which quantifies the
sphericity of a particle with its roundness (e.g., smoothness) [14, 25]. The
Krumbien-Sloss chart contains multiple illustrations of particles that are
compared against physical particles to obtain an estimate of both roundess and
sphericity. Particles that are both highly smooth and round are considered
high-quality material as the bulk physical properties of the material will be
homogeneous and predictable. Multiple problems exist with the Krumbien-Sloss
method. First, the method is manual, making it time consuming and dependent on
the individual performing the analysis. Second, there are multiple Krumbien-
Sloss charts, making the resulting quantification of the sample dependent on
the chart used. We will leverage the tools of shape space analysis to
automatically quantify the morphology of the samples shown in Figure 4 without
the use of a Krumbien-Sloss chart. We will also demonstrate how shape space
analysis can be used to develop statistical tests for changes in particle
morphology and for the implementation of dimensionality reduction techniques
such as principal component analysis (PCA).

###  3.1 Data Pre-Processing

We first take each of our images and extract an outline for each individual
particle. This is done through a common method using watershed segmentation in
the OpenCV package implemented in Python. We note in our analysis that these
images contain particles that are slightly overlapping or touching. This makes
it difficult to extract a perfect outline for each individual particle, thus,
we focus on particles in which we can extract a full outline. To compare each
particle using shape space analysis, we must first interpolate/resample the
outline of each particle so that each outline contains the same number of
points ensuring they are on the same pre-shape manifold. Two
resampled/interpolated outlines from the image in Figure 4(a) are shown in
Figure 5(a). Once the data is pre-processed we can translate, rotate, and
scale each particle outline, shown in Figure 5, so that we can measure
differences in their shape.

One difficulty faced in this analysis is that we are dealing with shapes that
do not have defined landmark points. Landmark points represent points on a
shape that can be consistently identified even when a shape is deformed. For
example, the triangle has three obvious landmarks which are the vertices of
the three corners of the triangle. Thus, for our data we must consider
translation, rotation, scaling, and reparameterization of the shape. This has
been deeply explored in the current literature, and one class of methods that
are highly effective in the analysis of data with no-defined landmarks are the
so-called _elastic metrics_ initially proposed in the work of Srivastava and
co-workers [49]. We will not cover this method in detail, but the elastic
metric proposed by Srivastava and co-workers represents parameterized curves
as points on a spherical pre-shape manifold through the square root velocity
function (SRVF). Given a real-valued curve
$f(t):[0,1]\rightarrow\mathbb{R}^{m}$ the SRVF is given as
$q:[0,1]\rightarrow\mathbb{R}^{m}$, which is defined as:

$\displaystyle q(t):=\frac{\dot{f(t)}}{\sqrt{||\dot{f(t)}||}}$ (3.25)

where $||\cdot||$ represents the Euclidean norm and $\dot{f(t)}$ represents
the derivative of $f$ with respect to $t$. The SRVF $q(t)$ is invariant to
translation, and can be further transformed to be invariant to scaling,
rotation, and reparameterization [49, 26]. Similar to our pre-shape space for
configuration matrices, this representation and its associated pre-shape
manifold can be used to construct geodesic distances and projections to a
tangent space, allowing us to perform statistical analysis and other data-
centric tasks like dimensionality reduction which we demonstrate in the
following section. To accomplish this we leverage the powerful Geomstats
Python package and its computational implementation of the elastic metric
within a discrete framework developed by Brigant [26, 32].

(a) Original contours

(b) Translated and scaled

(c) Reparameterized and rotated

Figure 5: Illustration of the translation, scaling, and optimal rotation of
two angular sand outlines.

###  3.2 Data Analysis

Once our data has been properly pre-processed, we can leverage the framework
introduced in the previous sections to peform both dimensionality reduction
and statistical analysis for our particle shapes. In our analysis, we are
interested in comparing our particles to a perfect circle, which represents a
particle with perfect sphericity and roundness. The geodesic distance from the
sphere provides a direct representation of the quality of the particles taken
from the sample. As can be seen in Figure 5(c) we should expect the ceramic
microspheres and the round sand samples to be geodesically close to the
perfect sphere with minimal variance in the distance and the angular sand
sample to be much farther away. First, we visualize the morphological
differences in the samples by performing principal component analysis (PCA) on
the horizontal tangent space of the pre-shape space manifold centered at the
circle. This form of PCA, introduced by Fletcher and co-workers, is often
referred to as principal geodesic analysis (PGA) as the method identifies
geodesic basis with respect to the manifold (via the tangent space) rather
than basis vectors in Euclidean space [15]. Figure 6 shows the first two
principal geodesics of PGA performed on the optimally aligned particle shapes.
We can see that the analysis reflects the intuition that the manufactured
ceramic microspheres are highly concentrated near the circle, with the round
sand being slightly farther from the circle, and the angular sand being of a
greater distance with a much larger variation in the distance.

Figure 6: Principal geodesic analysis of the particles outlines from angular
sand, round sand, and ceramic microsphere samples. The base point of defined
tangent space is the perfect circle. We project to the first two principal
geodesics that capture the most variance in the data with respect to the shape
space of our data. The projection illustrates the relative geodesic distances
to the perfect sphere for the three samples. The round sand and ceramic
microspheres have the tightest distribution around the perfect circle. We also
illustrate on the edges of the image the changes in morphology along the two
principal geodesics.

We can leverage these results to perform statistical analysis and understand
how the distance in the tangent space can be used to distinguish between the
samples, and whether or not there is a statistically significant difference
between each of the particle samples. Figure 7 shows box plots of the
distances in the tangent space from each particle to the perfect circle for
our samples. From the analysis, we see that both the mean and variance of the
geodesic distance from the circle for the ceramic microspheres and round sand
is much lower than the angular sand. To quantify this difference we can
perform a 2-sample t-test between the three different samples. The results of
the 2-sample t-test are shown in Table 1. We find that there is a
statistically significant difference in shape between the ceramic, the round
sand and the angular sand samples, with an $\alpha=0.05$. Figure 6 also
demonstates that we are able to capture morphological outliers in the data for
the ceramic microspheres and round sand samples. This suggests that the
automatic, computationally efficient nature of this method, coupled with its
statistical power, would make it amenable to the development of statistical
process control techniques for the continuous improvement and quality control
of processes in which shape is a key process parameter.

Figure 7: Box plots of the tangent space distance between each particle
outline and the perfect circle. We can see that the synthetic ceramic
microspheres have the tightest morphological distribution around the sphere,
with the round sound being slighty farther away, and the angular sand have the
largest mean distance and largest variation in morphology.

These results reflect the intuition around the shape of these three samples
gained from the photographs in Figure 4. Ceramic microspheres are specifically
manufactured to have consistent morphology, which is reflected in both the
mean distance to the perfect circle and in the tight grouping of the material,
but at the expense of a higher production cost. The round sand represents a
high-quality pure quartz sand that has consistent morphology often obtained
from sandstone, representing a more affordable alternative to man-made
materials. The angular sand is often not used in applications where the sand
requires homogeneous physical properties. However this can be very useful in
applications where resistance to shear and cushioning of the sand can be
beneficial, as is needed in equestrian arenas to reduce potential injuries to
horses [50].

| Ceramic | Round | Angular |
---|---|---|---|---
Ceramic | - |  |  |
Round | 0.045 | - |  |
Angular | $6.84e^{-9}$ | $2.40e^{-7}$ | - |
Table 1: t-test p-values between the different samples.

##  4 Discussion and Conclusion

We have presented a Riemannian geometric framework for the automatic
quantification of shape for materials. We show that shapes can be represented
as points on a non-linear Riemannian manifold. We can exploit the geometric
and mathematical properties of the Riemannian manifold to construct a measure
of the difference in shape between materials, known as shape space analysis.
Riemannian manifolds are equipped with a tangent space defined at each point
of the manifold, allowing us to project data from the non-linear manifold to a
linear (vector) space. This allows us to directly perform data analysis such
as hypothesis testing and dimensionality reduction which we demonstrate in a
test on granular material samples from Covia Corp. In our case study we show
that the Riemannian geometric framework can be applied automatically in the
analysis of granular material images. We show that this method provides
generalizable quantification of the differences in the shape of the individual
particles, which is in contrast with the subjective, manual, and bespoke
methods used commonly in the industry, such as the Krumbien-Sloss chart
method.

In future work we look to integrate this framework into the statistical
quality control of the shape of materials. This can be applied in the
mining/manufacturing of granular materials to control the overall quality of
the material, in the production of pharmaceuticals where the materials
morphology can dramatically impact the solubility and effectiveness of the
material, in the production of 3-dimensional printed materials, and in the
manufacturing of biomaterials such as cells or fungi where geometry of the
material is a key process parameter [39, 55, 10, 56, 6]. We will also explore
the use of the geometric shape of a material and statistics around its
morphology to predict bulk physical properties, such as porosity, hardness,
and solubility. To accomplish this we will also integrate these geometric
methods with _topological_ methods derived from the field of Topological Data
Analysis to provide a holistic, complementary analysis of a material’s
morphology [47, 48]. Furthermore, the shape of a material lies on a
differentiable manifold on which we can compute gradients and perform
optimization [40]. We look to leverage this structure to develop methods in
which we can change a material’s shape to perform multi-scale optimization and
design of materials and processes.

## References

  * [1] K. A. Alshibli, A. M. Druckrey, R. I. Al-Raoush, T. Weiskittel, and N. V. Lavrik.  Quantifying morphology of sands using 3d imaging.  Journal of Materials in Civil Engineering, 27(10):04014275, 2015\.
  * [2] F. Altuhafi, C. O’Sullivan, and I. Cavarretta.  Analysis of an image-based method to quantify the size and shape of sand particles.  Journal of Geotechnical and Geoenvironmental Engineering, 139(8):1290–1307, 2013.
  * [3] M. Bédérina, M. Khenfer, R. Dheilly, and M. Quéneudec.  Reuse of local sand: effect of limestone filler proportion on the rheological and mechanical properties of different sand concretes.  Cement and concrete research, 35(6):1172–1179, 2005.
  * [4] B. Bernard-Michel, S. Rohani, M.-N. Pons, H. Vivier, and H. S. Hundal.  Classification of crystal shape using fourier descriptors and mathematical morphology.  Particle & particle systems characterization, 14(4):193–200, 1997\.
  * [5] F. L. Bookstein.  Landmark methods for forms without landmarks: morphometrics of group differences in outline shape.  Medical image analysis, 1(3):225–243, 1997.
  * [6] D. Braga, L. Casali, and F. Grepioni.  The relevance of crystal forms in the pharmaceutical field: Sword of damocles or innovation tools?  International Journal of Molecular Sciences, 23(16):9013, 2022.
  * [7] E. Castillo and B. M. Colosimo.  Statistical shape analysis of experiments for manufacturing processes.  Technometrics, 53(1):1–15, 2011.
  * [8] A. Chowdhury, E. Kautz, B. Yener, and D. Lewis.  Image driven machine learning methods for microstructure recognition.  Computational Materials Science, 123:176–187, 2016.
  * [9] I. Czogiel, I. L. Dryden, and C. J. Brignell.  Bayesian molecular alignment using random fields.
  * [10] J. C. De Anda, X. Wang, X. Lai, K. Roberts, K. Jennings, M. Wilkinson, D. Watson, and D. Roberts.  Real-time product morphology monitoring in crystallization using imaging technique.  AIChE Journal, 51(5):1406–1414, 2005.
  * [11] H. Drolon, F. Druaux, and A. Faure.  Particles shape analysis and classification using the wavelet transform.  Pattern Recognition Letters, 21(6-7):473–482, 2000.
  * [12] I. L. Dryden, J. D. Hirst, and J. L. Melville.  Statistical analysis of unlabeled point sets: Comparing molecules in chemoinformatics.  Biometrics, 63(1):237–251, 2007.
  * [13] I. L. Dryden and K. V. Mardia.  Statistical shape analysis: with applications in R, volume 995.  John Wiley & Sons, 2016.
  * [14] EN:ISO standard 13503-2:2006, Petroleum and natural gas industries - Completion fluids and materials –Part: 2 Measurement of properties of proppants used in hydraulic fracturing and gravel-packing operations.  Standard, European Committee for Standarization, Mar. 2000.
  * [15] P. T. Fletcher, C. Lu, S. M. Pizer, and S. Joshi.  Principal geodesic analysis for the study of nonlinear statistics of shape.  IEEE transactions on medical imaging, 23(8):995–1005, 2004.
  * [16] A. M. Gabbert, J. A. Mondo, J. P. Campanale, N. P. Mitchell, A. Myers, S. J. Streichan, N. Miolane, and D. J. Montell.  Septins regulate border cell shape and surface geometry downstream of rho.  bioRxiv, pages 2021–04, 2021.
  * [17] M. Ge, F. Su, Z. Zhao, and D. Su.  Deep learning analysis on microscopic imaging in materials science.  Materials Today Nano, 11:100087, 2020.
  * [18] N. Guigui, N. Miolane, X. Pennec, et al.  Introduction to riemannian geometry and geometric statistics: from basic theory to implementation with geomstats.  Foundations and Trends® in Machine Learning, 16(3):329–493, 2023.
  * [19] H. Hafid, G. Ovarlez, F. Toussaint, P. Jezequel, and N. Roussel.  Effect of particle morphological parameters on sand grains packing properties and rheology of model mortars.  Cement and Concrete Research, 80:44–51, 2016.
  * [20] K. M. Z. Hossain, U. Patel, and I. Ahmed.  Development of microspheres for biomedical applications: a review.  Progress in biomaterials, 4:1–19, 2015.
  * [21] D. G. Kendall.  The diffusion of shape.  Advances in applied probability, 9(3):428–430, 1977.
  * [22] D. G. Kendall.  A survey of the statistical theory of shape.  Statistical Science, 4(2):87–99, 1989.
  * [23] D. G. Kendall, D. Barden, T. K. Carne, and H. Le.  Shape and shape theory.  John Wiley & Sons, 2009.
  * [24] H. Kim, J. Han, and T. Y.-J. Han.  Machine vision-driven automatic recognition of particle size and morphology in sem images.  Nanoscale, 12(37):19461–19469, 2020.
  * [25] W. C. Krumbein and L. L. Sloss.  Stratigraphy and sedimentation, volume 71.  LWW, 1951.
  * [26] A. Le Brigant.  A discrete framework to find the optimal matching between manifold-valued curves.  Journal of Mathematical Imaging and Vision, 61(1):40–70, 2019.
  * [27] M. S. Lim, D. C. Wijeyesekera, A. Zainorabidin, and I. Bakar.  The effects of particle morphology (shape and sizes) characteristics on its engineering behaviour and sustainable engineering performance of sand.  International Journal of Integrated Engineering, 4(4), 2012.
  * [28] J. I. Marden.  Analyzing and modeling rank data.  CRC Press, 1996.
  * [29] K. Mardia, F. Bookstein, and J. Kent.  Alcohol, babies and the death penalty: Saving lives by analysing the shape of the brain.  Significance, 10(3):12–16, 2013.
  * [30] S. M. Mercier, B. Diepenbroek, R. H. Wijffels, and M. Streefland.  Multivariate pat solutions for biopharmaceutical cultivation: current progress and limitations.  Trends in Biotechnology, 32(6):329–336, 2014.
  * [31] N. Miolane, M. Caorsi, U. Lupo, M. Guerard, N. Guigui, J. Mathe, Y. Cabanes, W. Reise, T. Davies, A. Leitão, et al.  Iclr 2021 challenge for computational geometry & topology: Design and results.  arXiv preprint arXiv:2108.09810, 2021.
  * [32] N. Miolane, N. Guigui, A. Le Brigant, J. Mathe, B. Hou, Y. Thanwerdas, S. Heyder, O. Peltre, N. Koep, H. Zaatiti, et al.  Geomstats: a python package for riemannian geometry in machine learning.  The Journal of Machine Learning Research, 21(1):9203–9211, 2020\.
  * [33] L. R. Monteiro, B. Bordin, and S. F. dos Reis.  Shape distances, shape spaces and the comparison of morphometric methods.  Trends in ecology & evolution, 15(6):217–220, 2000.
  * [34] A. Myers and N. Miolane.  Regression-based elastic metric learning on shape spaces of elastic curves.  arXiv preprint arXiv:2210.01932, 2022.
  * [35] S. Nuchitprasitchai, M. Roggemann, and J. M. Pearce.  Factors effecting real-time optical monitoring of fused filament 3d printing.  Progress in Additive Manufacturing, 2:133–149, 2017.
  * [36] P. O’HIGGINS.  The study of morphological variation in the hominid fossil record: biology, landmarks and geometry.  The Journal of Anatomy, 197(1):103–120, 2000.
  * [37] X. Pennec.  Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements.  Journal of Mathematical Imaging and Vision, 25(1):127–154, 2006\.
  * [38] X. Pennec.  Statistical computing on manifolds: from riemannian geometry to computational anatomy.  In Emerging Trends in Visual Computing: LIX Fall Colloquium, ETVC 2008, Palaiseau, France, November 18-20, 2008. Revised Invited Papers, pages 347–386. Springer, 2009.
  * [39] D. Quintanilla, T. Hagemann, K. Hansen, and K. V. Gernaey.  Fungal morphology in industrial enzyme production—modelling and monitoring.  Filaments in bioprocesses, pages 29–54, 2015.
  * [40] W. Ring and B. Wirth.  Optimization methods on riemannian manifolds and their application to shape space.  SIAM Journal on Optimization, 22(2):596–627, 2012.
  * [41] F. J. Rohlf.  Shape statistics: Procrustes superimpositions and tangent spaces.  Journal of classification, 16:197–223, 1999.
  * [42] S. Ross, D. Lamprou, and D. Douroumis.  Engineering and manufacturing of pharmaceutical co-crystals: a review of solvent-free manufacturing technologies.  Chemical Communications, 52(57):8772–8786, 2016.
  * [43] W. H. Sauer and M. K. Schwarz.  Molecular shape diversity of combinatorial libraries: a prerequisite for broad bioactivity.  Journal of chemical information and computer sciences, 43(3):987–1003, 2003.
  * [44] M. R. Singh, J. Chakraborty, N. Nere, H.-H. Tung, S. Bordawekar, and D. Ramkrishna.  Image-analysis-based method for 3d crystal morphology measurement and polymorph identification using confocal microscopy.  Crystal growth & design, 12(7):3735–3748, 2012.
  * [45] A. Smith, B. Laubach, I. Castillo, and V. M. Zavala.  Data analysis using riemannian geometry and applications to chemical engineering.  arXiv preprint arXiv:2203.12471, 2022.
  * [46] A. D. Smith, N. Abbott, and V. M. Zavala.  Convolutional network analysis of optical micrographs for liquid crystal sensors.  The Journal of Physical Chemistry C, 124(28):15152–15161, 2020\.
  * [47] A. D. Smith, P. Dłotko, and V. M. Zavala.  Topological data analysis: concepts, computation, and applications in chemical engineering.  Computers & Chemical Engineering, 146:107202, 2021.
  * [48] A. D. Smith and V. M. Zavala.  The Euler characteristic: A general topological descriptor for complex data.  Computers & Chemical Engineering, 154:107463, 2021.
  * [49] A. Srivastava, E. Klassen, S. H. Joshi, and I. H. Jermyn.  Shape analysis of elastic curves in euclidean spaces.  IEEE transactions on pattern analysis and machine intelligence, 33(7):1415–1428, 2010.
  * [50] M. Stavermann and E. Hernlund.  Equestrian surfaces part i–sand specifications.  2014\.
  * [51] M. Streefland, D. E. Martens, E. C. Beuvery, and R. H. Wijffels.  Process analytical technology (pat) tools for the cultivation step in biopharmaceutical production.  Engineering in Life Sciences, 13(3):212–223, 2013.
  * [52] A. K. Suryavanshi and R. N. Swamy.  Development of lightweight mixes using ceramic microspheres as fillers.  Cement and Concrete Research, 32(11):1783–1789, 2002.
  * [53] Y. Wang, G. Ma, J. Mei, Y. Zou, D. Zhang, W. Zhou, and X. Cao.  Machine learning reveals the influences of grain morphology on grain crushing strength.  Acta Geotechnica, 16(11):3617–3630, 2021.
  * [54] D. Weichert, P. Link, A. Stoll, S. Rüping, S. Ihlenfeldt, and S. Wrobel.  A review of machine learning for the optimization of production processes.  The International Journal of Advanced Manufacturing Technology, 104(5-8):1889–1902, 2019.
  * [55] H.-C. Wu and T.-C. T. Chen.  Quality control issues in 3d-printing manufacturing: a review.  Rapid Prototyping Journal, 24(3):607–614, 2018.
  * [56] W. Zheng and D. Tannant.  Frac sand crushing characteristics and morphology changes under high compressive stress and implications for sand pack permeability.  Canadian Geotechnical Journal, 53(9):1412–1423, 2016.