# Dynamic multi-object Gaussian process models: A framework for data-driven
functional modelling of human joints

Jean-Rassaire Fouefack, Bhushan Borotikar, Tania S. Douglas, Valérie Burdin,
and Tinashe E.M. Mutsvangwa, Jean-Rassaire Fouefack, Tania Douglas and Tinashe
Mutsvangwa are with the Division of Biomedical Engineering, University of Cape
Town, 7935, South Africa.
E-mail: ffcjea001myuct.ac.za<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Jean-Rassaire Fouefack, Valérie Burdin are with
the IMT-Atlantique and the Laboratory of Medical Information Processing (LaTIM
INSERM U1101), Brest, France.
E-mail<EMAIL_ADDRESS>Borotikar is with the
Laboratory of Medical Information Processing (LaTIM INSERM U1101), CHRU de
Brest, Brest, France and University of Brittany Occidental, Brest, France.
E-mail<EMAIL_ADDRESS>This work is based on research supported
by the National Research Foundation (NRF) of South Africa (grant no’s 105950
and 114393); the South African Research Chairs Initiative of the NRF and the
Department of Science and Technology (grant no 98788); the South African
Medical Research Council and the “Mission Enseignement Supérieur, Recherche et
Innovation”, Brest Métrople, France (grant no 17-178)

###### Abstract

Statistical shape models (SSMs) are state-of-the-art medical image analysis
tools for extracting and explaining features across a set of biological
structures. However, a principled and robust way to combine shape and pose
features has been illusive due to three main issues: 1) Non-homogeneity of the
data (data with linear and non-linear natural variation across features), 2)
non-optimal representation of the $3D$ motion (rigid transformation
representations that are not proportional to the kinetic energy that move an
object from one position to the other), and 3) artificial discretization of
the models. In this paper, we propose a new framework for dynamic multi-object
statistical modelling framework for the analysis of human joints in a
continuous domain. Specifically, we propose to normalise shape and dynamic
spatial features in the same linearized statistical space permitting the use
of linear statistics; we adopt an optimal 3D motion representation for more
accurate rigid transformation comparisons; and we provide a 3D shape and pose
prediction protocol using a Markov chain Monte Carlo sampling-based fitting.
The framework affords an efficient generative dynamic multi-object modelling
platform for biological joints. We validate the framework using a controlled
synthetic data. Finally, the framework is applied to an analysis of the human
shoulder joint to compare its performance with standard SSM approaches in
prediction of shape while adding the advantage of determining relative pose
between bones in a complex. Excellent validity is observed and the shoulder
joint shape-pose prediction results suggest that the novel framework may have
utility for a range of medical image analysis applications. Furthermore, the
framework is generic and can be extended to n$>$2 objects, making it suitable
for clinical and diagnostic methods for the management of joint disorders.

###### Index Terms:

Combined 3D shape and pose analysis, Generative models, Human joints, 3D
statistical motion representation, Shoulder, Pattern recognition.

##  1 Introduction

Finding a minimum subset $\mathbf{D}_{0}$ with $N_{s}$ patterns that can
generally and specifically represent a given set of structures $\mathbf{D}$
with $N$ feature classes is a frequently researched area [1]. A common
strategy is to estimate a subspace of homogeneous patterns (for example,
geometric shapes) through models. This inductive learning approach leads to
models that allow for automated observation and interpretation of the learned
patterns. In the medical image analysis domain, pattern analysis has to
contend with two existing problems. First, patterns may present across non-
homogeneous feature classes. An example is morphofunctional analysis where
determining the changes in relative spatial dynamics of different anatomical
structures adds an additional pattern class on top of the morphological
(geometric) pattern class [2]. Second, patterns in some feature classes may
belong to non-linear spaces. The current trend is to formalise the problem as
an inverse problem and use deep learning based generative models [3, 4].
However deep generative models tend to have high specificity but low
generality; in some cases explaining out-of-training data inputs with high
confidence while being wrong [5, 6, 7]. Furthermore deep generative models
require large numbers of training examples to be successful. Finally, the
interpretability of deep generative models remains complex. Thus, non-
parametric analysis of non-linear variations of non-homogeneous features for
different medical image analysis tasks (shape completion, morphometric
analysis, segmentation) with the same model, remains desirable.

With regard to analysis of shape, a well established and understood formalism
for analysing 3D geometric variation in a linearized statistical space exists
in the form of statistical shape modelling (SSM). SSMs typically model the
data in Euclidean vector space using principal component analysis (PCA) that
treats feature variation changes as a linear combination of local translations
only [9, 8, 10, 11]. While limited in capturing the non-linearity in shape
space, the validity of this linearization for rigid shapes has been codified
in the literature for single anatomical structures [9, 8, 11]. Efforts for
faithfully representing the non-linearity of shape space have been reported
but have not become mainstream due to computational inefficiency and a lack of
robustness [12]. Recently, Luthi et al (2018)[11] introduced a generalisation
of SSMs, referred to as Gaussian process morphable models (GPMMs). In this
framework, the parametric low-dimensional model was represented as a Gaussian
process over deformation fields obtained from training examples. In contrast
to discrete models that are dependent on artificial discretization, GPMMs are
inherently continuous, that is, permitting of the arbitrary discretization of
the domain on which the model is defined. However, GPMMs do not embed inter-
object shape correlation, nor the anatomo-physiological relationship between
articulating objects. This is because current SSMs and GPMMs are unable to
represent, in a statistically robust and intuitive way, an articulating
anatomical complex composed of several rigid substructures which can move
relative to each other.

Similar to the shape case, pose variation using both linear and non-linear
statistical descriptions have been reported [14, 17, 16, 18, 15, 13, 19].
Previous studies have reported the embedding of relative pose variation of
articulating bones into SSM frameworks, in the form of statistical shape and
pose models (SSPM) [14, 17]. Others have reported statistical pose models
(SPM) only; that is, the model of variation across spatial orientations
without inclusion of their shapes [20]. A straightforward approach to model
shape and pose together is to model joint flexibility implicitly by
incorporating joint motion in statistical space [18]. However, this assumes
that pose variations are a statistical property of anatomy. Another approach
is to concatenate pose and anatomical features into a joint vector and apply
linear PCA to obtain joint principal components [16]. However, this describes
the pose transformations in Euclidean vector (linear) space; contrary to
reality as they belong to a manifold and thus need to be described using non-
linear statistical methods.

Various efforts have been reported in the literature to model shape and pose
features together. Bossa et Olmos (2007) proposed parametric low-dimensional
models of several subjects in different poses with the pose variation
representation obtained through principal geodesic analysis (PGA); a non-
linear extension of PCA [19]. Shape and pose features of the brain cortex were
concatenated in a long vector and standard multivariate statistics were
extracted [17]. Fitzpatrick et al (2011) [15] reported on a characterisation
of the relationship between shape and contact mechanics of the patella and
femur using an SSM of the articulating surfaces between both bones. The
authors considered the patello-femoral contact space as a new object and
modelled it as a standard SSM. Chen et al (2014) proposed an articulated SSM
for segmentation of the wrist from computed tomographic (CT) images without
concatenating the shape and pose features; their model was a combination of an
SSM and an SPM developed using PGA [21]. All the above shape and pose PGA
based models used a standard representation (Rodrigues, quaternion and Euler
representation) for modelling the rigid transformations describing the pose.
Moreau et al (2017) reported that the standard representation is limited in
its ability to model rigid transformations of non-compact objects (with
respect to a shape compactness measure). They instead proposed a new norm for
statistical analysis of rigid transformations that is robust for analysing
non-compact objects. However, they only demonstrated a new norm for SPMs
without considering shape and pose variability analysis. A non-linear shape
modelling approach based on differential coordinates (thus avoiding global
rigid transformation through local rotation and stretching) has been reported
[22, 12]. While this approach captures the inherent non-linearity in shape
variation, it still suffers from artificial discretization which prevents
marginalization of the resultant probabilistic model. Additionally, it does
not allow for modelling multi-object structures.

We hypothesise that for a principled approach to multi-object shape and pose
analysis, the following limitations identified in the literature need to be
addressed : 1) features of different classes (morphology and pose) need to be
properly normalised to obtain an homogeneous analysis space; 2) $3D$ motion
requires an optimal representation to induce an efficient associated norm for
comparison of rigid transformations and 3) models need to be built in a
continuous domain. We propose a novel framework that addresses all the above
limitations by modelling shape and pose variation of complex articulated
objects composed of multiple rigid objects. We refer to the proposed modelling
framework as Dynamic Multi-object-Gaussian Process Modelling (DMO-GPM)111The
code for this project is available online for research use:
https://github.com/rassaire/DmoGpm.

The structure of the rest of the paper is as follows. In section 2, we
identify and describe previous and related work in the domain of SSPMs and
GPMMs. In section 3, we describe the DMO-GPM framework and explain 1) how the
framework is extended and generalised to a joint consisting of $N$ objects and
the spatial positions are projected in a linear space using PGA, 2) how the
new framework addresses the inconsistent representation of $3D$ motions
typical in modelling joints with non-compact structures such as long bones
using the standard representation, and 3) how we adopt an energy displacement
representation that is more suitable for reconstruction of the objects of
interest. In section 4, we validate the framework using synthetic data with
predefined and accurately known shape and spatial dynamics. We also propose an
automatic pipeline for dynamic multi-object prediction using a combination of
Markov Chain Monte Carlo sampling and DMO-GPM approach which allows for
registering a model to pathological or incomplete multi-objects while avoiding
local minima. In section 5, we apply the framework to complex shoulder joint
data to demonstrate the utility of the approach compared to the state-of-the-
art sample kernel-based GPMM.

Figure 1: Proposed multi-object generative model architecture. Multi-shape and
pose features are obtained after establishing anatomical correspondence. The
pose feature linearization is performed through tangent space projection using
$\log$ mapping. A Gaussian process is defined over deformation fields obtained
from combined shape and pose features. Multi-object samples are obtained by 1)
using the shape deformation field component to instantiate shape instances;
and 2) a retrieval of their spatial orientation using rigid transformations
associated to the pose deformation field component.

##  2 Previous and related work

Before describing our DMO-GPM framework, we summarise the main concepts behind
traditional statistical shape models, multi-object statistical shape and pose
models, and GPMMs on which we will capitalise to synthesize the dynamic multi-
object framework.

###  2.1 Statistical shape and pose models

Let’s consider a multi-object data-set $\\{\\{S_{i}^{j}\\}_{i=1\dots
n}\\}_{j=1\dots N}$ with $n$ examples, where each example is composed of $N$
objects. For example, the gleno-humeral joint with training data samples and
having $2$ rigid bony objects - scapula and humerus. The underlying assumption
in the development of multi-object statistical shape and pose models (SSPMs)
is that there exists a representational space where the relationship between
shape and pose can be learnt from the data-set. This concept of SSPMs is well-
adapted for modelling articulated objects in motion. The $j^{th}$ object of
the $i^{th}$ data sample can be represented by $S_{i}^{j}$, a vector of
discrete landmarks with known correspondence between the examples
$\\{S_{i}^{j}\\}_{i=1\dots n}$ of the $j^{th}$ object, the reference object
$\Gamma_{S}^{j}$ and $\vec{\Gamma}_{S}^{j}$ a coordinate vector of points of
$\Gamma_{S}^{j}$, typically, one selected object from the database that should
be closest to the mean. The rigid transformation that aligns the $j^{th}$
object of the $i^{th}$ example to its reference object $\Gamma_{S}^{j}$ is
defined as

$h^{j}_{i}=\underset{h\in SE(3)}{\arg
min}\|h(S^{j}_{i})-\vec{\Gamma}_{S}^{j}\|^{2}$

where $SE(3)$ is the space of rigid transformations defining a Lie group with
$\|.\|$ being the sum of distance metric. The Lie group (smooth manifold)
defined by the system of objects is the direct product of the $N$ rigid
transformation groups denoted as $SE(3)^{N}$. Shape variations are relatively
small and are traditionally assumed to belong to a vector space, while pose
(rigid transformations) variations belong to a Riemannian manifold which is
not a vector space. In contrast to shape models that are built using a
classical PCA, the pose models are built using PGA which is in essence a PCA
in a space tangential to the manifold[19, 23]. Several studies use PGA to
perform linear statistics on a Riemannian manifold [19, 21, 14, 13]. To
perform the PGA, the $i^{th}$ transformation for the $j^{th}$ object is
represented in this tangent space to $SE(3)$ with respect to the mean rigid
transformation $\bar{h}_{j}$ using the log mapping:

$\displaystyle log_{\bar{h}_{j}}h_{i}^{j}=[p_{1},\dots,p_{m}]^{T}$ (1)

where $(p_{i})_{i=1\dots m}$ are the parameters of the rigid transformation
$\bar{h}_{j}^{-1}\circ h_{i}^{j}$ representing the Euler’s, Rodrigues’s or
quaternion parameters. These representations describe the $3D$ rotations in
3D. Rodrigues and quaternion represent a $3$D rotation as an axis vector and a
rotation angle while Euler represents it with three angles. For the rest of
the paper these representations will be called standard representations (SR).
To map the $i^{th}$ transformation back to the manifold, an exponential
mapping is used:

$h_{i}^{j}=\bar{h}_{j}exp(log_{\bar{h}_{j}}h_{i}^{j})$

To build the SSPMs, shape and pose features of the $i^{th}$ example are
concatenated in joint feature vectors [17] by:

$X^{i}=[w^{s}X^{i}_{S},w^{p}X^{i}_{p}]^{T}$

where $X^{i}_{S}=[S_{i}^{1},\dots,S_{i}^{N}]^{T}$ and
$X^{i}_{p}=[log_{\bar{h}_{j}}h_{i}^{1},\dots,log_{\bar{h}_{j}}h_{i}^{N}]^{T}$,
$w^{s}$ and $w^{p}$ denote shape and pose weightings, respectively, that
balance the relative importance between shape and pose within the model. The
resulted norm of the SR is defined as:

$\displaystyle\|h_{i}^{j}\|_{SR}=t^{T}\;t+s\,r^{T}\,r$ (2)

where $r$ and $t$ are rotation and translation vectors respectively; and $s$
the scaling factor. To build the model, an additional assumption is that the
joint feature residual data-set $\\{X^{i}\\}_{1\dots n}$ can be modelled by a
normal distribution similar to point distribution models PDMs [9] as follows:

$(s,p)\sim\mathcal{N}(\mu,\Sigma)$

where $\mu$ is the joint mean and $\Sigma$, the covariance matrix, are
estimated from the data-set $\\{X^{i}\\}_{1\dots n}$ as:

$\mu=\bar{X}=n^{-1}\sum_{i=1}^{n}X^{i}$
$\Sigma=\bar{X}=n^{-1}\sum_{i=1}^{n}(X^{i}-\bar{X})(X^{i}-\bar{X})^{T}$

The model is obtained by performing PCA which in a probabilistic
interpretation leads to:

$X=\bar{X}+\sum_{k=1}^{n}\alpha_{k}\sqrt{\lambda_{k}}\,e^{k}$

where $e^{k}=[e^{k}_{S},e^{k}_{p}]^{T}$ is the combined $k^{th}$ eigenvector
of the covariance matrix (associated to the $k^{th}$ eigenvalue
$\lambda_{k}$), obtained from single value decomposition (SVD) and
$\alpha_{k}\sim\mathcal{N}(0,1)$. Then, the SSPM is split into two parts, the
shape component, typically a PDM, defined as follows:

$\displaystyle
X_{\alpha_{k}}=\bar{X_{s}}+\sum_{k=1}^{n}\alpha_{k}\sqrt{\lambda_{k}}\,e^{k}_{s}$
(3)

and the pose component, typically SPMs, defined as:

$\displaystyle
h_{\alpha_{k}}=\bar{h}\,exp\left(\sum_{i=1}^{n}\alpha_{k}\sqrt{\lambda_{k}}e^{k}_{p}\right)$
(4)

A new instance of a discrete SSPM is obtained as
$s_{p}=h_{\alpha_{k}}(X_{\alpha_{k}})$ which derives the shape at a relative
spatial position.

###  2.2 Gaussian process morphable models

Gaussian process morphable models[11] are an extension of PDMs in equation 3.
In contrast to PDMs, GPMMs model deformations defined on a continuous domain,
although they can always be discretized to obtain a model that is
mathematically equivalent to a PDM. The advantage of such model is that after
training, the statistical descriptions may be marginalized (or customized
kernels could be used) to a region of interest which will provide better
interpretability of the underlying information. Additionally, having freedom
in defining the covariance function may be useful because explicit bias
kernels can be added based on expert knowledge. This GPMM flexibility
overcomes a typical limitation of traditional SSM whereby insufficient
training examples result in loss of model generality [24]. In GPMMs, shape
correspondence, established across $n$ examples, leads to deformation fields
$u_{i},i=1\dots n$ defined from the reference shape $\Gamma_{S}$ to $i^{th}$
examples in the training data. These deformations are then modelled as
Gaussian Processes ($GP$) on the reference shape domain $\Omega_{S}$
($\mathbb{R}^{3}\supset\Omega_{S}\supset\Gamma_{S}$) thus making them
independent of discretization. The GPMMs are formulated as follows:

$\displaystyle\\{\Gamma_{R}+u(\Gamma_{R})\\}\text{ with }u\sim GP(\mu,k)$ (5)

where for all points $x,y$ of $\Omega$, the mean and covariance (kernel
function) are defined as:

$\displaystyle\mu(x)=n^{-1}\sum_{i=1}^{n}u_{i}(x)$ (6) $\displaystyle
k(x,y)=n^{-1}\sum_{i=1}^{n}(u_{i}(x)-\mu(x))(u_{i}(y)-\mu(y))^{T}$

The function $u_{i}$ is the deformation field associated with the dense
correspondence established across the shape parameters between the reference
and the $i^{th}$ example. The deformations can be represented in terms of an
orthogonal set of basis functions $\\{\Phi\\}_{m=1}^{M},M\leq n$:

$\displaystyle
u(x)=\sum_{m=1}^{M}\alpha_{m}\sqrt{\lambda_{m}}\Phi_{m}(x),\text{ with
}x\in\Omega_{S},\,\alpha\sim\mathcal{N}(0,1)$ (7)

where $(\lambda_{m},\Phi_{m})$ are the eigenvalue/eigenfunction couples. These
couples are obtained using Karhunen-Loève expansion and Nyström approximation
[25, 26], as explained in [11]. That leads to a probabilistic PCA-based model
of shape directly on the deformations.

##  3 Dynamic multi-object Gaussian process models

In this section we develop the dynamic multi-object Gaussian process model
(DMO-GPM) framework that incorporates:

  * •

an homogeneous analysis space (linear) for shape and spatial orientation
features,

  * •

a kinetic energy based $3D$ motion representation for more accurate comparison
of rigid transformations,

  * •

a continuous representation to describe feature variation across the training
dataset.

The DMO-GPM framework is aimed at modeling multi-objects (with $N$ objects)
and their relative spatial positions as a deformation field
$u^{MO}=\\{u^{j}_{S},u^{j}_{K}\\}_{j=1}^{N}$ from a reference joint
$\Gamma=\cup_{j=1}^{N}(\Gamma^{j}_{S}\cup\Gamma^{j}_{K})\subset\mathbb{R}^{3}$,
where a $j^{th}$ object shape of the joint, at its relative spatial position,
is represented as

$h^{j}(\\{x+u^{j}_{S}(x)|x\in\Gamma^{j}_{S}\\})$ (8)

with the rigid transformation $h^{j}$ defined as

$h^{j}=\underset{h\in SE(3)}{\arg
min}\|\vec{u}_{h}-\vec{\Gamma}_{K}^{j}\|^{2}$

for some deformation
$u^{MO}:\Omega=\cup_{j=1}^{N}(\Omega^{j}_{S}\cup\Omega^{j}_{K})\subset\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$,
with $\Gamma\subset\Omega$ and $\vec{u}_{h}$ being the coordinate vector of
the set of points $h(\\{x+u^{j}_{K}(x)|x\in\Gamma^{j}_{K}\\})$. We model
deformations as a Gaussian process $u^{MO}\sim GP(\mu,k)$ with mean function
$\mu:\Omega\longrightarrow\mathbb{R}^{d}$ and covariance (kernel) function
$K:\Omega\times\Omega\longrightarrow\mathbb{R}^{3\times 3}$.

The DMO-GPMs are built through 4 steps. First, the anatomical correspondence
needs to be established across the multi-objects to guarantee the homology of
shape-pose features. Second, the optimal representation of spatial dynamics
associated with shapes is obtained. Third, linear and non-linear
transformations are normalized to obtain an homogeneous space for linear
statistics. Finally, the dynamic multi-object is modelled in a continuous
domain.

###  3.1 Multi-object correspondence establishment

We define a dynamic multi-object (joint) as a set of several biological
structures (or objects) from different family of shapes (e.g humerus, scapula
etc.) and their relative spatial position. Establishing correspondence across
several structures while keeping their anatomo-physiological relationship is
not a straight-forward task. Here we extend the framework presented in our
previous work [27] for two biological structures to $N$ ($N\in\mathbb{N}$)
biological structures. For a single biological structure, amongst several
different approaches [28, 29, 30], correspondence can typically be achieved
through rigid-registration and subsequently non-rigid registration performed
using free-form deformation (FFD) models such as GPMMs [11]. For multi-
objects, not only shape correspondence should be established across objects,
but also the relative spatial position between them should be preserved. To
register $N$ objects while keeping their anatomo-physiological relationship,
we assume that one of the $N$ objects remains fixed while the others move
relative to the fixed one. Individual FFD models ($N$ models) are built for
each object family and used to register objects of their respective family.
After the registration, each inverse of the rigid registration transformation
is used to retrieve the original spatial position of moving objects. This
operation retrieves the anatomo-physiological relationship between moving
objects. Rigid transformation of the fixed object is finally applied to all
moving objects to obtain their anatomo-physiological relationship with the
fixed object. Dense correspondence (shape and spatial dynamics) is established
across the data-set by using the same FFD as reference for the registration.

Let $T_{j}$ be the rigid transformation that transforms the $j^{th}$ object
($S^{j},j=1\dots N-1$) to the reference $\Gamma_{S}^{j}$, $T$ the rigid
transformation that transforms the fixed object ($S^{j_{0}}$) to its reference
$\Gamma^{j_{0}}$. The $j^{th}$ registered object, $S_{r}^{j}$, with its
spatial position retrieved is obtained as:

$S_{r}^{j}=T\circ T_{j}^{-1}(\underset{u\in DF}{\arg
min}SM(u(\Gamma_{S}^{j}),S^{j})+\beta\|u\|^{2})$

where $DF$ is the space of deformation fields, $\|.\|$ is a norm in $DF$, $SM$
is a similarity measure and $\beta$ a regularisation coefficient.

###  3.2 Optimal representation of $3D$ motion

The literature on SSPMs usually emphasizes on the common shape-pose space
required in medical image processing. The standard representation of rigid
transformations may not always encapsulate the kinetic energy necessary to
move points on an object from one position to another, particularly for non-
compact shapes. This reduces the efficiency of this parameterization for
statistical analysis of rigid transformations [13]. Thus, we propose to extend
the kinetic energy based representation by Moreau et al (2017) in a continuous
space and use it to encode a rigid transformation in DMO-GPM. We refer to this
new transformation representation as the energy displacement representation
(EDR).

We represent a multi-object complex as a set of $N$ objects; each object
defined by its shape and its spatial position parameters. A multi-object is
represented as a concatenated vector space and denoted as
$X=[X_{S},X_{K}]^{T}$ where $X_{S}=[S^{1},\dots,S^{N}]^{T}$ and
$X_{K}=[K^{1},\dots,K^{N}]^{T}$ with $S^{i}$ and $K^{i}$ representing the
point domain of the shape and the spatial orientation features for the
$i^{th}$ object of the joint, respectively.

Now let us assume there are $n$ multi-objects examples in the training data-
set. Let us explicitly define the rigid transformation representation which is
the displacement field representation that extends [13] in a continuous domain
and will be used to define joint features representations. We define the EDR
of the $i^{th}$ example of the object $j$ on its reference spatial position
domain $\Omega_{K}^{j}\subset\mathbb{R}^{3}$ with their respective reference
shape $\Gamma^{j}_{S}$ and reference pose points $\Gamma_{K}^{j}$. The EDR
$\delta^{j}_{i}$ and its associated metric $d$ are defined as:

$\displaystyle\delta^{j}_{i}(x)=h^{j}_{i}(x)-Id(x),x\in\Omega_{K}^{j},i=1\dots
n$ (9) $\displaystyle h^{j}_{i}=\underset{h\in SE(3)}{\arg
min}\|h(S^{j}_{i})-\vec{\Gamma}_{S}^{j}\|^{2}$ $\displaystyle
d^{2}(h^{j}_{i},Id)=\delta^{j}_{i}(\Omega_{K}^{j})^{T}\delta^{j}_{i}(\Omega_{K}^{j})$

where $SE(3)$ is a lie group of rigid transformations.

###  3.3 Normalising dynamic multi-object features for homogeneous
representation

To normalise shape and pose features, we define unified representations of
shapes and their relative spatial dynamics. These representations are joint
deformation fields over which a Gaussian process (GP) will be defined (section
3.4). A reference multi-object
$\Gamma=\cup_{j=1}^{N}(\Gamma^{j}_{S}\cup\Gamma^{j}_{K})\subset\mathbb{R}^{3}$
is chosen from the training examples. This joint should be representative of
the data-set, that is, approximately similar to the mean of the data-set. A
joint deformation field of the $i^{th}$ example in the training data-set is a
pair of deformation fields ($u^{i}_{S}$, $u^{i}_{K}$), where $u^{i}_{S}$ maps
a point $x$ of the reference shape objects $\cup_{j=1}^{N}\Gamma^{j}_{S}$ to
the point $u^{i}_{S}(x)$ of its corresponding shape example; and $u^{i}_{K}$
maps a point $x$ of the reference pose objects $\cup_{j=1}^{N}\Gamma^{j}_{K}$
to the point $u^{i}_{K}(x)$ of its corresponding pose example. The $i^{th}$
joint deformation field is defined

$\begin{dcases}u^{i}_{S}(x)\text{ if
}x\in\Omega_{S}=\cup_{j=1}^{N}\Omega^{j}_{S},\\\ u^{i}_{K}(x)\text{ if
}x\in\Omega_{K}=\cup_{j=1}^{N}\Omega^{j}_{K}.\\\ \end{dcases}$ (10)

where

$\begin{dcases}u^{i}_{S}(x)=u^{j}_{i}(x),\text{ if }x\in\Omega_{S}^{j}\\\
u^{i}_{K}(x)=h^{j}_{i}(x),\text{ if }x\in\Omega_{K}^{j}.\\\ \end{dcases}$ (11)

Before obtaining residual features (around the mean) that will represent
normalised features, the joint function of the mean shape and spatial
transformations has to be defined. We estimate this mean function, $\mu_{MO}$,
using the mean pose deformation fields as:

$\begin{dcases}\frac{1}{n}\sum_{i=1}^{n}u^{j}_{i}(x),\text{ if
}x\in\Omega^{j}_{S},\\\ \mu^{j}_{K}(x)\approx\bar{h}^{j}(x),\text{ if
}x\in\Omega_{K}^{j}.\\\ \end{dcases}$ (12)

with $\bar{h}^{j}=\underset{h\in SE(3)}{\arg
min}\Big{\|}h\Big{[}\frac{1}{n}\sum_{i=1}^{n}\delta^{j}_{i}(\vec{\Gamma}_{K}^{j})\Big{]}-\vec{\Gamma}_{K}^{j}\Big{\|}^{2}$.
It should be noted that the mean $\bar{h}$ is obtained using the mean of the
vector field (thanks to EDR), in contrast to SR where the Fréchet mean is
estimated through numerical optimization.

Thus far, pose representations across the training data-set belong to a
manifold which is a non-linear space. Rigid transformations are projected
(that is linearization of relative spatial transformations) onto a linearized
space to the manifold at the mean rigid transformation (tangent space) using
the EDR (eq 9), in order to obtain the unified space of shape and pose
features. Relative spatial transformations are linearized through $\exp/\log$
bijective mapping.

Let us define the $\exp/\log$ functions associated to the EDR between two
rigid transformations $h^{j}_{i_{1}}$, $h^{j}_{i_{2}}$:

$\displaystyle\log[h^{j}_{i}]=\delta_{i}^{j},j=1\dots N\text{ and }i=1\dots n$
(13)
$\displaystyle\log_{h^{j}_{i_{1}}}[h^{j}_{i_{2}}]=\log[(h^{j}_{i_{2}})^{-1}\circ
h^{j}_{i_{1}}]$ $\displaystyle\exp[\delta^{j}_{i}]=\underset{h\in SE(3)}{\arg
min}\|h[\delta^{j}_{i}(\vec{\Gamma}_{K}^{j})]-\vec{\Gamma}_{K}^{j}\|^{2}$
$\displaystyle\exp_{h^{j}_{i_{1}}}[\delta^{j}_{i_{2}}]=h^{j}_{i_{1}}\circ\exp[\delta^{j}_{i_{2}}].$

The projected joint deformation fields onto the vector space over which the GP
will be defined are $\\{u_{i}^{MO}\\}_{i=1}^{n}$

$u_{i}^{MO}=\begin{dcases}u^{i}_{S}(x)-\mu_{MO}(x),x\in\Omega_{S}\\\
\log_{\mu_{MO}}[u^{i}_{K}](x),x\in\Omega_{K}.\end{dcases}$ (14)

###  3.4 Dynamic multi-object modeling in a continuous domain

From a training dataset of $n$ examples with $N$ objects each, the deformation
fields $\\{u_{i}^{MO}\\}_{i=1}^{n}$ defined on the continuous domain
$\left(\Omega=\cup_{j=1}^{N}(\Omega^{j}_{S}\cup\Omega^{j}_{K})\right)$ are
modelled by a Gaussian process defined as:

$u\sim GP(\mu_{MO},k_{MO})$ (15)

where $k_{MO}$ is the kernel function defined as:

$k_{MO}(x,y)=\frac{1}{n}\sum_{i=1}^{n}u^{MO}_{i}(x)u^{MO}_{i}(y)^{T},x,y\in\Omega\\\
$ (16)

The multi-object deformation fields $u^{MO}$ are represented by an orthogonal
set basis functions denoted by $\\{\phi_{m}^{MO}\\}_{m=1}^{M}$, $M\leq n$ and
defined as:

$\phi_{m}^{MO}(x)=\begin{dcases}\phi_{m}^{S}(x)\text{ if }x\in\Omega_{S},\\\
\phi_{m}^{K}(x)\text{ if }x\in\Omega_{K}.\\\ \end{dcases}$ (17)

where $\phi_{m}^{S}$ and $\phi_{m}^{K}$ represent the shape and the energy
displacement component of the $m^{th}$ basis function.

The $j^{th}$ object is modelled as a combination of shape deformation field
and a rigid transformation (obtained using exponential map) as defined below:

$\begin{dcases}\sum_{m=1}^{M}\alpha_{m}\sqrt{\lambda^{MO}_{m}}\phi_{m}^{S_{j}}(x),\text{
if }x\in\Omega^{j}_{S},\\\
h^{j}:\exp_{\mu_{MO}}\left(\sum_{m=1}^{n}\alpha_{m}\sqrt{\lambda^{MO}_{i}}\phi_{m}^{K_{j}}(x)\right),\text{
if }x\in\Omega_{K}^{j}.\\\ \end{dcases}$ (18)

where $\alpha\sim\mathcal{N}(0,1)$ and the pair
$(\lambda_{m}^{MO},\phi_{m}^{MO})_{m=1\dots M}$ representing the $m^{th}$
eigenvalue and eigenfunction couple. These couples are obtained using
Karhunen-Loève expansion and Nyström approximation as in equation 7. The
$m^{th}$ eigenfunction $\phi_{m}^{MO}$ is a continuous multi-object
deformation field defined on $\Omega$, in the tangent space to $SE(3)$.

For example, figure 2 shows how a $j^{th}$ object is sampled from the model. A
shape instance of the $j^{th}$ object and its associated spatial position can
be generated from our probabilistic model for
$\alpha=(\alpha_{m})_{m}\sim\mathcal{N}(0,1)$ as:

$\displaystyle\textrm{shape\&pose}^{j}=h^{j}(\textrm{shape}^{j})$ (19) with
$\displaystyle\textrm{shape}^{j}=\sum_{m=1}^{M}\alpha_{m}\sqrt{\lambda^{MO}_{m}}\phi_{m}^{S_{j}}(\Gamma_{S}^{j})$

where $\Gamma_{S}^{j}\subset\Omega_{S}^{j}$ is the reference shape of the
$j^{th}$ object.

Figure 2: Obtaining shape sample at its associated spatial orientation. Top:
the rigid transformation is obtained from pose deformation filed sample.
Bottom: the sample sample is obtain from shape deformation filed and its
spatial orientation is retrieved using the rigid transformation generated

####  3.4.1 Marginalising the model

As in GPMMs, our DMO-GPM framework allows to operate on a specific part of the
domain for a given model using marginal probability properties, DMO-GPMs still
have these properties.

Let us assume that the DMO-GPM is a joint probabilistic model denoted by
$p_{\textrm{MO}}(\Omega)$. The marginal probabilistic model is defined as

$\displaystyle
p_{\textrm{MO}}(\Omega_{i}):=(\mu_{MO},k_{MO})|_{\Omega_{i}\subset\Omega},$

obtained from the joint probability in a specific domain.

Another important aspect of the DMO-GPM framework is its potential of making
the model “class specific”, that is, making the variations across a specific
class constant. For example, this model includes static multi-object models,
that is, setting the pose parameters constant to obtain a classic SSM of the
joint. Conversely, it includes a motion model, that is, setting the shape
parameters constant which leads to a pose-only statistical model.

####  3.4.2 Modelling with kernels

Another advantage of the DMO-GPMs is the possibility of increasing the
variability around the mean without retraining the model as in GPMMs. The
number of training examples may be insufficient to provide a subspace of
shapes to represent all possible shape variations within the family, which
introduces bias. The extension of the variability using a bias kernel can
overcome this limitation. Our framework allows for the understanding of such
variability and defining an appropriate kernel that encodes the variability. A
new kernel can be added to the trained model kernel without a need of reusing
the training examples. An example of kernel that can increase smoothness in a
model is the Gaussian kernel defined in a specific domain
$\Omega^{j}_{S}\text{ or }\Omega^{j}_{K},j=1,\dots,N$ by

$k_{\Omega}^{(s,\sigma)}(x,y)=s\exp(\|x-y\|^{2}/\sigma^{2}),x,y\in\Omega$

where the parameter $s$ defines the scale of the average error and where
$\sigma$ the range over which deformations are correlated [11]. The new DMO-
GPM kernel can be a combination of the pre-trained kernel with different
Gaussian kernels in each domain defined by

$k_{new}=k_{\textrm{MO}}+\sum_{j=1}^{N}I_{3\times
3}k_{\Omega_{j}}^{(s_{j},\sigma_{j})},\Omega_{j}\in\\{\Omega^{j}_{S},\Omega^{j}_{K}\\}$

where the $I_{3\times 3}$ is a $3\times 3$ unit matrix.

###  3.5 Using the framework in analysis tasks

In order to demonstrate the utility of the framework in prediction from the
observation of target examples, we adopted the Monte Carlo Markov chain (MCMC)
and Metropolis Hastings algorithms as in [21, 24]. Markov Chain Monte Carlo
methods are computational tools to perform approximate inference with
intractable probabilistic models by making the use of another distribution
from which one can sample, in order to ultimately draw samples from the true
posterior distribution model.

Let us consider target objects with arbitrary relative positions to each
other. The goal of the MCMC fitting process is to find the best model
parameters $\theta$ (a multi-object deformation field in our case) that
optimally represents the target example. The uncertainties related to pose
estimation are embedded in the DMO-GPM and therefore need not be modelled
separately. The target multi-object example may consist of partial objects;
for instance, images of bones with slight cracks caused by osteomalacia or any
incomplete anatomical structure on images due to pathology or due to image
artifacts. The advantage of the DMO-GPM is that when confronted with a missing
part on one object, the other objects in the target example can be used as
additional observations since the correlation between the objects is embedded
in the model. The DMO-GPM with MCMC fitting could be used to obtain a
posterior model of one object while another is missing.

Let $T_{o}$ be the observation (target multi-object complex),
$T_{o}=\\{T_{o}^{j}\\}$ is a set of $N$ objects at various relative positions.
Note that $T_{o}$ could be partial i.e. we may have less than $N$ objects. The
posterior model is estimated using the Bayes rule:

$p(\theta|T_{o})=\dfrac{p(T_{o}|\theta)p(\theta)}{p(T_{o})}$

where $p(T_{o})$ is intractable, $p(\theta)$ is the model and
$p(T_{o}|\theta)$ the likelihood. To estimate the posterior DMO-GPM, a multi-
object likelihood is proposed.

####  3.5.1 Multi-object likelihood

Full shape and pose estimations require a likelihood that captures within-
object and between-object features. We propose the multi-feature class
comparison likelihood that has a global and local component. The global
component ensures the prediction of all feature classes while each local
component ensures the predominance of its feature class in the global
likelihood. The $j^{th}$ local likelihood $l_{j}$ compares the observation of
the target $T_{o}^{j}$ with its instance $\theta(\Omega^{j})$ generated from
the model using independent feature-wise comparison. It can be seen as a
probability distribution of possible shape-pose samples evaluated for the
target and it is defined as:

$l_{j}(\theta,T_{o}^{j})=p(T_{o}^{j}|\theta(\Omega^{j}=\\{X_{i}\\}_{1}^{P}))=\prod_{x_{i},x_{j}\in
T_{o}^{j}}^{P}p(x_{j}|\theta(x_{i}))$ (20)

The global likelihood $l$ compares the observation of target $T_{o}$ with the
model instance $\theta(\Omega)$ and it is defined as:

$l(\theta,T_{o})=p(T_{o}|\theta(\Omega=\\{X_{i}\\}_{1}^{PN}))=\prod_{x_{i},x_{j}\in
T_{o}}^{PN}p(x_{j}|\theta(x_{i}))$ (21)

####  3.5.2 Obtaining the optimal parameter set $\theta$

A sample $\theta^{{}^{\prime}}$ is accepted as a new state $\theta$ for the
$j^{th}$ object with the metropolis decision probability $\beta_{j}$ and the
target $T$ with the probability $\beta$ defined by:

$\beta_{j}=\min\left\\{\dfrac{l_{j}(\theta,T_{o}^{j})p(\theta)}{l_{j}(\theta^{{}^{\prime}},T_{o}^{j})p(\theta^{{}^{\prime}})}\dfrac{Q(\theta^{{}^{\prime}}|\theta)}{Q(\theta^{{}^{\prime}}|\theta)},1\right\\}$
(22)

$\beta=\min\left\\{\dfrac{l(\theta,T_{o})p(\theta)}{l(\theta^{{}^{\prime}},T_{o})p(\theta^{{}^{\prime}})}\dfrac{Q(\theta^{{}^{\prime}}|\theta)}{Q(\theta^{{}^{\prime}}|\theta)},1\right\\}$
(23)

where
$Q(\theta^{{}^{\prime}}|\theta)=\mathcal{N}(\theta^{{}^{\prime}}|\theta,\Sigma_{\theta})$
is the proposal generator. Finally, the proposal $Q$ is fed through a sequence
of Metropolis acceptance decisions:

$\\{\beta_{j},\beta\\},j=1\dots,N$ (24)

The predicted target is the instance of the posterior model with the highest
probability $\beta$.

##  4 Validating the framework on synthetic data

To validate the DMO-GPM framework, we created synthetic data set with
precisely defined modes of shape and pose variation. We avoided real
biological data because these typically exhibit complex shape and pose
variations making them difficult to use for bench-marking and validation. Our
synthetic data consisted of surface mesh data of a ”lollipop” as defined in
[31]. The lollipop has a shaft, a neck and head. We used this template to
generate $30$ surface data with meshes composed of $6000$ vertices; all with
identical surface topology. The only difference in shape between all generated
meshes was the length of the major axis defining the ellipsoid of the head of
the lollipop. This was varied uniformly in the range from $1~{}mm$ to
$−30~{}mm$, so the shape variation was entirely encoded using one degree of
freedom. The shaft and neck areas were identical for all synthetic data.

###  4.1 Synthetic experiments with static objects

To evaluate the potential of our framework in modelling shape correlation
between different objects, we created artificial joints using two lollipops
per scenario. Each joint was composed of two lollipops with major axes for
corresponding pairs of lollipops of $r_{1}$ and $r_{2}$, for object 1 and
object 2, respectively. The span of the data-set of joints was created by
varying $r_{1}$ and $r_{2}$ as
$(r_{1},r_{2})=(i/10,j/10),i=16,\dots,30,j=1,\dots,15$ creating a negative
correlation between objects 1 and 2 (figure 3 (left)).

Figure 3: Illustration of the lollipop data. Left) Shape variation of the
first and the second objects, the first object parameter $r_{1}$ increasing in
the direction of the arrow and the second object parameter $r_{2}$ decreasing
in the direction of the arrow. Right) An example motion of a lollipop joint,
the fixed object (grey) and the moved object (white) at different angle
$\theta$.

####  4.1.1 Empirical validation of correlation between objects

The correlation between samples from the DMO-GPM and the one between the
synthetically created examples were compared to know whether our model could
establish the correlation between objects in the training data if it exists.
Our model captured the correlation prescribed in the training data-set. The
DMO-GPM showed a non-zero slope at $0.0001$ significance level with a negative
correlation coefficient of $-0.99$ as in the training data-set. The regression
curve of the training examples and our model examples were the same,
illustrating the precision to which our model captured the correlation between
objects of the synthetic joints.

####  4.1.2 Prediction using static model and MCMC

Prediction using statistical models remains challenging for medical image
processing tasks because it requires observations that may not be available.
The models developed have the potential to make predictions of one object
without or with limited observations (which are the features available) of
that object through leveraging information from the other object in a complex.
The assumption is that there is a good correlation between the two objects.
For the rest of the study we define four types of observations of target
objects (see table I) that will be used for predictions.

TABLE I: Taxonomy of observations of target object during predictions Type I |  The full object to be predicted is available
---|---
Type II |  Partial object to be predicted is available
Type III |  Anatomical features (e.g. landmarks) of the object to be predicted are available
Type IV |  No knowledge of the object to be predicted is available

To test this hypothesis, we generated $5$ joints of lollipops with the
parameters $r_{1}$ and $r_{2}$ between the minimum and the maximum parameters
that were in the training data-set. For each joint, object 1 was predicted
from the observation of object 2 only and conversely. We compared the results
obtained for each object with Type I observation to the ones obtained with
Type IV observation. Table II shows the root mean square (RMS) distance and
the Hausdorff distance (HD) errors of the predictions. The average RMS error
was $0.4\pm 0.06$ and $0.22\pm 0.02$ mm for the prediction of object 1 and
object 2 with Type I observation, respectively. The average HD error was
$0.9\pm 0.22$ and $0.82\pm 0.2$ mm for the prediction of object 1 and object 2
with Type IV observation, respectively. These errors are relatively large
compared to those predicted with Type I observation. However, this result
obtained with type IV observation shows that the correlation embedded in the
model provides additional information that can be useful in getting accurate
prediction.

TABLE II: Prediction errors of the lollipop joint objects with Type I and IV observations using MCMC and DMO-GPM  |  RMS (mm) |  HD (mm)
---|---|---
Object 1 with Type I observation |  0.002$\pm$2.2e-6 |  0.014$\pm$1.7e-3
Object 1 with Type IV observation |  0.4$\pm$0.06 |  0.9$\pm$0.22
Object 2 with Type I observation |  0.001$\pm$2.9e-6 |  0.01$\pm$0.001
Object 2 with Type IV observation |  0.22$\pm$0.02 |  0.82$\pm$0.2

###  4.2 Synthetic experiments with dynamic joints

To evaluate the potential of our framework in modelling shape and spatial
dynamics together, we simulated joint motion. First, we created a data-set
where shape does not correlate with motion. For each joint generated above, we
rotated the second lollipop (object 2) relative to the first one (object 1)
using the Euler’s angle convention for describing rigid transformations
$(\varphi,\theta,\psi)$. The second lollipop was moved in the $yz-plane$
corresponding to Euler’s angles $(0,\theta,0)$; where only $\theta$ was varied
by four angles
($\theta=\frac{1}{5}\pi,\frac{2}{5}\pi,\frac{3}{5}\pi,\frac{4}{5}\pi$)
defining a motion (figure 3 (right)). Defining the same motion for each joint
created a data-set where shape was uncorrelated with motion because each pair
of shapes performed the same movement. Secondly, we created a data-set where
the motion is constrained by the shape. This was done by simulating the same
movement as above but the maximum angle was inversely proportional to $r_{2}$
( $\theta_{max}=-\frac{3}{4r_{2}}\pi$), thus creating correlation between
shape and motion.

####  4.2.1 Non-correlation and correlation of shape-pose variation

We built a model using the synthetic training data-set with no shape and pose
variation correlation. The model behaved as expected; the first two principal
components (PC) accounted for $98.2\%$ of the total variation. The first PC
explained shape-only variation without any pose variation and counted for
$87.7\%$ and the second explained only pose variation and counted for $10.5\%$
of the total variation. Figure 4 shows the shape variation of the first and
second PC from $-3\sigma$ to $+3\sigma$ standard deviations about the mean. It
can be observed in the first PC that the pose is fixed at $-12~{}rad$ and only
the shape parameters $r_{1}$ and $r_{2}$ vary from $41.5$ to $77.3~{}mm$ for
the first object and from $38.7$ to $4.7~{}mm$ for the second which
corresponds to the correlation induced in the training data-set. In the first
PC, the second lollipop head length ($r_{2}$) decreases up to its sharp where
it starts being unrealistic due to the large sampling range (up to $+3\sigma$)
as indicate by the star (*) in figure 4. In the second PC the shape parameters
are fixed at the mean shape which are $r_{1}=59.4~{}mm$ for the object 1 and
$r_{2}=21.6~{}mm$ for the second object. The pose parameter $\theta$ varies
from $-0.2$ to $-1.5~{}rad$ representing abduction motion in the training
data-set. The shape and pose variation explained by two different PCs mirrors
the shape and the pose variation across the training data-set confirming that
the model represents the training data-set.

Figure 4: First and second PCs of the DMO-GPM with no correlation between
shape and pose variation. The fist PC accounts only for shape variation while
the second PC account only for pose variation.

We also built a model using synthetic training data in which pose variation
was constrained by shape. The model also behaved as expected. The first two
PCs accounted for $99\%$ of the total variation. The first PC was shape
variation as well as pose variation and accounted for $92.3\%$ and the second
explained pose-only variation and accounted for $6.7\%$ of the total variation
as shown in Figure 5. A visual inspection of the first PC reveals that shape
parameters $r_{1}$ and $r_{2}$ vary from $44.2$ to $80.1~{}mm$ for the first
object and from $36.2$ to $1.6~{}mm$ for the second object and that the pose
variation is inversely proportional to the second shape parameter; mirroring
the training data-set. In the second PC the shape parameters are fixed at the
mean shape parameters which are $59.4~{}mm$ for the first object and $21.6$
for the second object and only the pose parameter $\theta$ varies from $-0.01$
to $-1.4~{}rad$. The shape and pose variation explained by the first PC show a
correlation between shape and pose variation similar to that in the training
data-set.

Figure 5: First and second PCs of the DMO-GPM with correlation between shape
pose variation. The first PC accounts for both shape and pose variation
explaining the correlation between the two set of features. The second PC
accounts for pose variation only.

####  4.2.2 Comparison of DMO-GPM with EDR and SR

A DMO-GPM using EDR was compared to the one using the SR to know whether the
model with EDR better explains the spatial dynamic of the objects. We sampled
the PC describing the motion simulated and compared this with motion
prescribed in the training data. Figure 6 shows first PC and the third PC
which account for motion in the DMO-GPM with EDR and DMO-GPM with the SR,
respectively. The motion described by the DMO-GPM is more similar to that
prescribed in the training data compared to the motion described by DMO with
the SR. For example, the position of the mean shape is at a lower angle in the
DMO with SR compared to that in the training motion. The sample poses from
$-2$ and $+2$ standard deviations about the mean are completely deviated from
the training data motion plane. It can be observed that the mean pose for the
model with SR is not different from that of the training data and the pose
sampls are not distributed around the mean; this is because the SR does not
use object features (is independent of the object features).

Figure 6: Comparison of the motion sampled from DMO-DPM with EDR with that
using DMO-GPR with SR. Left: the motion from the training data. middle: the
motion sampled from the model with EDR. Right: the motion sampled from the
model with the SR

.

We also used the DMO-GPM built with the EDR and the one built with the SR to
predict lollipop joints with the second object at various poses, that is,
predicting each with type I observation. The same synthetic training examples
used to test object correlation above (see table II) were used. The test data-
set includes two positions of object 2 relative to the first object: $0~{}rad$
and $\frac{4}{5}\pi~{}rad$ for the models with no shape-pose correlation, and
$0~{}rad$ and $-\frac{3}{4r_{2}}\pi~{}rad$ for the models with shape-pose
correlation. Figures 7 and 8 show the box-plots of the surface to surface
distance RMS and the HD errors in predicting the second object. The DMO-GPM
with EDR outperforms DMO-GPM with SR when there is shape-pose correlation. For
no shape-pose correlation, the DMO-GPM with SR performed slightly better than
the one with EDR. This slight outperformance occurs even when the motion
described by the model is not statistically derived from the training motion,
while some poses still belong to the training poses.

Figure 7: Comparison of the shape plus pose prediction using DMO-GPM with EDR
and DMO-GPM with SR. Red: the RMS errors of the prediction using DMO-GPM with
EDR. Blue: the RMS errors of the prediction using DMO-GPM with SR. Figure 8:
Comparison of the shape plus pose prediction using DMO-GPM with EDR and the
one with SR. Red: shows the HD errors of the prediction using DMO-GPM with
EDR. Blue: Shows the HD errors of the prediction using DMO-GPM with SR.

##  5 Experiments

To further demonstrate our framework in clinical and medical image analysis
tasks, we applied our method to the shoulder joint; the most complex joint of
the human body. The complexity of the shoulder joint can be seen from both the
shape and the dynamics points of view. The joint contains the scapula, one of
the most complex bones in shape terms; and the humerus, one of the longest
bones of the body. This complexity already confounds the separate morphometric
analyses of each bone, let alone when both bones need to be analysed together.
Furthermore, in motion terms, the shoulder joint has the most degrees of
freedom. This leads to difficulties in segmenting its constituent bones at
their relative spatial positions.

###  5.1 Prediction of shoulder joint mechanics using DMO-GPM and MCMC

We applied DMO-GPMs to shoulder joint data (composed of two articulating
bones; scapula and humerus) to demonstrate different scenarios in which the
framework could be used. The data consisted of $3D$ mesh surfaces segmented
from computed tomography (CT) images of the bilateral shoulders of fresh
cadavers collected from the Division of Clinical Anatomy and Biological
Anthropology, Faculty of Health Sciences, University of Cape Town, South
Africa. Institutional ethics approval was granted for this study (Approval No:
HREC 546/2017). A total of $36$ shoulder joints were used; of these data, $8$
were from female and $10$ were from male remains. None presented with shoulder
pathology at the time of death. The age of the imaged decendents ranged from
21 to 90 years. The CT volumes were imported into the medical modelling tool
(Amira v6.2.0, http://www.fei.com/) for segmentation and 3D surface rendering.

####  5.1.1 Correspondence establishment across shoulder joints

In order to establish correspondence across joints within the data-set, we
used the surface meshes to create a median virtual shape for the humerus and
scapula, separately, using iterative median closest point-Gaussian mixture
model (IMCP-GMM) [26]. The mean-virtual shapes were used as morphologically
unbiased templates to built separate GPMMs. These models were used to register
all the samples across the data-set using the joint registration method
described in section 3.1. Table III shows the average surface to surface RMS
and HD between the original meshes and their estimations after the
registration process for both humeri and scapulae. The registration errors
were low for both humerus and scapula indicating that adequate anatomical
correspondence was established across examples in the training data-sets.
However, the humerus registration errors were low relative to the scapula
registration errors which could be attributed to the greater shape complexity
of the scapula relative to the humerus [26].

Figure 9: Reconstruction of a partial target scapula and humerus. Observed part of the shoulder for prediction (white). The Hausdorff distance reconstruction errors show better reconstruction with DMO-GPM (middle) than with SSMs (right) TABLE III: Shoulder joint training data-set registration errors. The average RMS and the HD of scapula and humerus examples. |  RMS (mm) |  HD (mm)
---|---|---
Humeri |  0.73$\pm$0.18 |  2.62$\pm$0.84
Scapulae |  1.11$\pm$0.2 |  3.01$\pm$0.8

###  5.2 Simulation of shoulder joint motion

Since the CT scans obtained were from cadavers, the humeri were in random
positions relative to their corresponding scapulae. Some of these positions
were not realistic (physiologically), therefore using them for motion analysis
would have led to an unrealistic study. To overcome this issue, we simulated
artificial motion through a biomechanics protocol.

To simulate the motion, the humerus was considered to be the mobile object
while the scapula was fixed. A random shoulder joint from the in-
correspondence data-set was used as a template and artificial pose variations
were obtained in four steps [27]: 1) determining the humeral head centre of
rotation using a sphere fitting algorithm [32]; 2) defining the International
Society of Biomechanics (ISB) recommended coordinates for both scapula and
humerus bones [16]; 3) obtaining the neutral position of the humerus by
aligning the humerus and scapula coordinate systems and then shifting the
scapula ISB origin to the glenoid centre; and 4) generating ten poses by
moving the humerus around the glenoid centre in the $yz$-plane (simulating
abduction). Five poses were generated for each joint
$(0,\frac{\pi}{6},\frac{\pi}{6},\frac{\pi}{3},\frac{\pi}{2},\frac{2\pi}{3})$.

###  5.3 Predicting shoulder joint shape features

To demonstrate the performance of our model in predicting shoulder anatomy, we
compared the results obtained using DMO-GPM to those using single structure
SSMs built as GPMM of the humerus. The prediction was done using MCMC. For
this purpose, a marginalised shape-only DMO-GPM of the shoulder was extracted
from the global DMO-GPM since only shape features were being predicted. We
demonstrated the prediction of the humerus from the observation of the scapula
in different scenarios. We created a test data-set, consisting of $10$
shoulder joints ($5$ from the training data and $5$ outside the training set)
for shoulder joint shape feature prediction.

First, we predicted the humerus with Type IV observation, that is, only having
the scapula of the patient. Table IV shows the errors in predicting the humeri
without any information at all using DMO-GPM. Regardless of somewhat large
errors, it can be concluded that if there is a strong correlation between the
scapula and humerus, a meaningful prediction of one given the other is
possible. These predictions result from the correlation that may exist between
scapula and humerus shapes captured by our modelling framework. Such
predictions would not be possible with single structure models or SSMs.
Second, we predicted humeri with Type III observation; two anatomical
landmarks (tips of the lateral and medial epicondyle) were selected on the
distal humerus targeted for prediction. Our aim here was to predict the humeri
when presented with a scapula and the two humerus landmarks using the DMO-GPM.
Again we compared the results of the DMO-GPM predictions of the humeri with
predictions made using a humerus SSM. Figure 10 shows the box-plots of
prediction errors of the humerus using DMO-GPM and single structure humerus
SSM. The prediction errors with DMO-GPM were low compared to errors for
predictions with the SSM. The outperformance of the DMO-GPM relative to the
single structure SSM could be attributed to the additional humerus knowledge
afforded through its corresponding scapula. Finally, humeri were predicted
with Type II observation. We simulated humerus fracture scenarios by removing
some part of the humeral surfaces and obtained partial humeri. The same
testing was performed between DMO-GPM and single structure SSM of the humeri
for prediction of the full humeri from partial humeri. Figure 11 shows the
box-plots of the RMS and the HD prediction errors of the humerus using DMO-GPM
and SSM. The DMO-GPM outperformed the single structure SSM and figure 9 shows
a reconstructed scapula with DMO-GPM and individuals SSMs. Again this
performance could be due the shape correlation that exists between the humerus
and the scapula that can be considered to carry latent information about the
humerus.

We predicted scapulae from the observation of partial scapulae using DMO-GPM
and GPM; that is, scapula fractures were simulated on the distal part of the
scapula. The DMO-GPM and the SSM of the scapula were used to predict those
scapulae. Figure 13 shows the prediction errors, where DMO-GPM outperformed
the scapula SSM.

TABLE IV: Prediction errors of the humerus with Type IV observations; average
RMS and HD DMO-GPM
---
Humerus with Type IV observation | RMS (mm) | HD (mm)
| 2.9$\pm$1.18 | 14.0$\pm$5.23
Single structure GPMM
Humerus with Type IV observation | RMS (mm) | HD (mm)
| Not possible
Figure 10: Comparison of the prediction of the humerus with Type III
observation SSM with that using DMO-GPM. Blue: the HD (left) and the RMS
distance (right) of the prediction using SSM. Red: HD (left) and the RMS
distance (right) of the prediction using DMO-GPM. Figure 11: Comparison of the
prediction of the humerus with Type II observation SSM with that using DMO-
GPM. Blue: the HD (left) and the RMS distance (right) of the prediction using
SSM. Red: the HD (left) and the RMS distance (right) of the prediction using
DMO-GPM Figure 12: Reconstruction of target shoulders. Left: The DMO-GPM.
Middle: target shoulder at various poses (from $0$ to $\dfrac{2\pi}{3}~{}rad$)
to be predicted. Right: Hausdorff distance and angle errors for the predicted
shoulder using DMO-GPM. Figure 13: Comparison of the prediction of the scapula
with Type II observation SSM with that using DMO-GPM. Blue: the HD (left) and
the RMS distance (right) of the prediction using SSM. Red: the HD (left) and
the RMS distance (right) of the prediction using DMO-GPM

###  5.4 Predicting shoulder joint motion

Figure 14: Prediction of the shoulder joint motion. DMO-GPM with EDR fitting
compared to that with SR

To evaluate the performance of the DMO-GPM in predicting training motion, the
simulated abduction motion was predicted. The scapula and the humerus (type I
observation) were predicted with the humerus at different poses relative to
the scapula. Five poses were generated within the motion range
$(0,\frac{3}{10}\pi,\frac{1}{2}\pi,\frac{7}{10}\pi,\frac{4}{5}\pi)$, these
test poses simulated abduction as in the training data-set. The prediction was
done using DMO-GPM with EDR and with SR in order to compare the performance of
the two models in predicting shape at their relative positions. Figure 14
shows the box-plots of the prediction errors obtained using pose marginalised
DMO-GPM with EDR and with SR. The DMO-GPM with the EDR outperformed that with
SR. For both scapula and humerus, the median error of the prediction with the
EDR model was less than that with SR. Figure 12 shows the reconstructed
shoulder using the DMO-GPM.

##  6 Conclusion

Understanding patterns and explaining them in a statistical framework so that
we can induce meaningful observations is a vast area of research. This study
presents a framework that tackles the existing problems in a combined
statistical representation of multi-object shape and pose structures by
providing a homogeneous yet simple way of building data-driven generative
models. The common aim in statistical modelling of multi-objects is to
efficiently explain shape variability and spatial dynamics across a population
of multi-objects in a low dimensional space. This study focused bringing the
two sets of feature classes involved (shape and pose) into a common linearized
space, thus finding the homogeneous analysis space. While multiple methods
exist to statistically represent the shape variations in linear space, its
pose counterpart is often ignored in the literature, which may lead to
incorrect variability of such shape-pose models. The proposed framework not
only uses a homogeneous and linear statistical representation of both shape
and pose features, it also embeds the mode in a continuous domain thereby
making the model generative. This type of framework in especially needed in
the medical imaging domain where structures of interest can be human
articulating joints. Human joint structures represent non-linear and non-
homogeneous feature classes in terms of its bone shapes and poses. Optimal and
homogeneous statistical representation of the variation across these
structures using a data-driven modelling framework is of utmost importance to
maintain the clinical accuracy of the inferences drawn from the model. In this
study, we have illustrated that DMO-GPM modeling framework can have valid
representations of shape and pose combined for multi-object structures.

We have provided a systematic validation of the DMO-GPM framework using
closely-controlled synthetic data (lollipop) and illustrated its application
is experimental medical data (shoulder joint). This framework is able to
incorporate more than two objects in a compact representation using Gaussian
process which means the developed models provide a continuous analysis space
with variations of shape and spatial dynamics. Marginalizing DMO-GPMs allow us
to analyze and/or predict shape-only or pose-only or combined
extractions/generations. This important characteristic of DMO-GPMs allow us to
implement the framework in multiple applications ranging from pre-surgical
planning to population based implant designs.

Results of validation with a synthetic data-set provide the efficiency of the
DMO-GPM framework in explaining correlations across a given population of
dynamic multi-objects. The DMO-GPM framework can accurately explain shape
correlation between objects and correlation between shape and spatial
variation, an advantage of the multi-object models in prediction compared to
single structure SSMs. Furthermore, DMO-GPM also preserves the optimal
representation of 3D motion using the EDR approach as illustrated on the non-
compact lollipop data-set. We also applied our framework on human shoulder
data to address one of the challenges in medical image processing, that is,
prediction from observations. DMO-GPM illustrated the prediction ability of
one of the shoulder bones from the observation of the other. The prediction
would be more accurate if more observations - partial or full - are available.
For example, estimating the humerus from the observation of the clavicle and
the scapula would be better than from the observation of the scapula only.
This could be of interest in surgical settings such as shoulder arthroplasty,
especially when one of the joint bones is partially or completely missing (due
to trauma).

The challenge of automatic assessment of pose variation in medical image
analysis can be due to the inability to standardize pose during image
acquisition. This is true for the same individual with multiple images or
different individuals. We evaluated our model in this scenario by prediction
of a shoulder joint with the humerus at different relative positions to
scapula, around a predetermined joint center, simulating abduction. Results
suggest that DMO-GPM can be incorporated into such scenarios for automatic
segmentation pipelines. Future work would involve extension of DMO-GPM with an
intensity model [33]. This would allow the use of the model for automatic
segmentation of $3D$ images (MRI and CT scans), especially dynamic MRI, with
due consideration to pose variations. Furthermore, DMO-GPM could be used to
optimize surgery. For example, by exploiting anatomical correlation between
joint objects, accurate joint segmentation could be achieved for the design of
the better joint implants. Additionally, parametric optimization could be
investigated, which is computationally less expensive than MCMC and could
improve the clinical utility of the models. Moreover, the proposed framework
could be used to answer biological questions.

##  Appendix A Notations

$k_{\Omega}^{(s,\sigma)}$ | Gaussian kernel defined in $\Omega$ with
---|---
| the scale $s$ and variance $\sigma$

$\mathcal{N}(0,1)$ | normal distribution with
---|---
| zero and unit covariance
$S^{j}_{i}$ |  $j^{th}$ object of the $i^{th}$ example
| in the training data
$n$ | number of the training examples
$N$ | number of objects constituting a
| training example
$h^{j}$ | rigid transformation assisted to the
$j^{t}h$ | shape sample from the DMO-GPM
$SE(3)$ | space of similarity transformations
$X^{i}$ | concatenated vector of shape and pose
| features of the $i^{th}$ example
$e^{k}$ | the combined eigenvector of shape ($e^{k}_{S}$)
| and pose ($e^{k}_{p}$) for SSPMs
$\mu$ | population mean
$\Sigma$ | population covarinace
$\alpha$ | Gaussian random vector
$X_{\alpha_{k}}$ | shape sample from SSPMs
$h_{\alpha_{k}}$ | pose sample from SSPMs
$u_{i}$ | deformation field associated to
| the $i^{th}$ example
$\Omega_{S}$ | reference shape domain
$\Omega_{K}$ | reference pose domain
$\Gamma_{S}$ | reference shape
$\lambda_{m}$ |  $m^{th}$ eigenvalue of the GPMMs
$\phi_{m}$ |  $m^{th}$ eigenvector of the GPMMs
$S^{j}_{0}$ |  $j^{th}$ object of the fixed object family
$\Gamma^{j_{0}}$ | reference of the fixed object family
$u^{MO}$ | multi object deformation field
$S^{j}$ |  $j^{th}$ object of the training example
$T$ | rigid transformation between
| the reference $\Gamma_{S}^{j}$ and the object $S^{j}$
$T^{j}$ | rigid transformation between
| the reference $\Gamma^{j}_{0}$ and the object $S^{j_{0}}$
$S^{j}_{r}$ | registered object corresponding
| to the object $S^{j}$
$X$ | multi-object concatenated vector of shape
| and spatial position
$K^{i}$ | spatial position displacement field of the
|  $i^{th}$ example
$X_{K}$ | spatial position component of $X$
$X_{S}$ | shape component of $X$
$\sigma^{j}_{i}$ | energy displacement representation
| associated to $S_{i}^{j}$
$\Omega_{K}^{j}$ | reference associated to the $K^{j}$
$u_{i}^{MO}$ | multi-object deformation field
| associated to $i^{th}$ example
$\mu_{MO}$ | DMO-GPM mean
$\phi^{S}_{m}$ | shape component of $\Psi^{MO}_{m}$
$\phi^{K}_{m}$ | pose component of $\Psi^{MO}_{m}$
$T_{o}$ | target multi-object example to be predicted
$r_{1}$ | shape parameter of the $1^{st}$ synthetic object
$r_{2}$ | shape parameter of the $2^{nd}$ synthetic object
$\varphi$ | Euler angles in X-axis
$\theta$ | Euler angle in Y-axis
$\psi$ | Euler angle in Z-axis
$p_{\textrm{MO}}$ | DMO-GPM
$\Gamma_{S}^{j}$ | reference of the $j$’s objects
$\vec{\Gamma}_{S}^{j}$ | a vector coordinate of the points of $\Gamma_{S}^{j}$
$h^{j}_{i}$ | rigid transformation between
| the reference $\Gamma_{S}^{j}$ and the object $S^{j}_{i}$

## Acknowledgments

The authors would like to thank Marcel Lüthi (Department of Mathematics and
Computer Science, University of Basel, Basel, Switzerland) for invaluable
input and discussion in preparation of this study.

## References

  * [1] Li, Y. and Maguire, L., 2010. Selecting critical patterns based on local geometrical and statistical information. IEEE transactions on pattern analysis and machine intelligence, 33(6), pp.1189-1201.
  * [2] Gorczowski, K., Styner, M., Jeong, J.Y., Marron, J.S., Piven, J., Hazlett, H.C., Pizer, S.M. and Gerig, G., 2009. Multi-object analysis of volume, pose, and shape using statistical discrimination. IEEE transactions on pattern analysis and machine intelligence, 32(4), pp.652-661.
  * [3] Toshev, A. and Szegedy, C., 2014. Deeppose: Human pose estimation via deep neural networks. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 1653-1660).
  * [4] Ouyang, W., Chu, X. and Wang, X., 2014. Multi-source deep learning for human pose estimation. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (pp. 2329-2336).
  * [5] Huang, G., Liu, Z., Van Der Maaten, L. and Weinberger, K.Q., 2017. Densely connected convolutional networks. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 4700-4708).
  * [6] Webster, R., Rabin, J., Simon, L. and Jurie, F., 2019. Detecting overfitting of deep generative networks via latent recovery. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (pp. 11273-11282).
  * [7] Nalisnick, E., Matsukawa, A., Teh, Y.W., Gorur, D. and Lakshminarayanan, B., 2018. Do deep generative models know what they don’t know?. arXiv preprint arXiv:1810.09136.
  * [8] Blanc, R. and Székely, G., 2012. Confidence regions for statistical model based shape prediction from sparse observations. IEEE Transactions on medical imaging, 31(6), pp.1300-1310.
  * [9] Cootes, T.F., Taylor, C.J., Cooper, D.H. and Graham, J., 1992. Training models of shape from sets of examples. In BMVC92 (pp. 9-18). Springer, London.
  * [10] Mutsvangwa, T., Burdin, V., Schwartz, C. and Roux, C., 2015. An automated statistical shape model developmental pipeline: application to the human scapula and humerus. IEEE Transactions on Biomedical Engineering, 62(4), pp.1098-1107.
  * [11] Lüthi, M., Gerig, T., Jud, C. and Vetter, T., 2018. Gaussian process morphable models. IEEE transactions on pattern analysis and machine intelligence, 40(8), pp.1860-1873.
  * [12] Von Tycowicz, C., Ambellan, F., Mukhopadhyay, A. and Zachow, S., 2018. An efficient Riemannian statistical shape model using differential coordinates: With application to the classification of data from the Osteoarthritis Initiative. Medical image analysis, 43, pp.1-9.
  * [13] Moreau, B., Gilles, B., Jolivet, E., Petit, P. and Subsol, G., 2017. A New Metric for Statistical Analysis of Rigid Transformations: Application to the Rib Cage. In Graphs in Biomedical Image Analysis, Computational Anatomy and Imaging Genetics (pp. 114-124). Springer, Cham.
  * [14] Anas, E.M.A., Rasoulian, A., John, P.S., Pichora, D., Rohling, R. and Abolmaesumi, P., 2014, March. A statistical shape+ pose model for segmentation of wrist CT images. In Medical Imaging 2014: Image Processing (Vol. 9034, p. 90340T). International Society for Optics and Photonics.
  * [15] Fitzpatrick, C.K., Baldwin, M.A., Laz, P.J., FitzPatrick, D.P., Lerner, A.L. and Rullkoetter, P.J., 2011. Development of a statistical shape model of the patellofemoral joint for investigating relationships between shape and function. Journal of biomechanics, 44(13), pp.2446-2452.
  * [16] Smoger, L.M., Fitzpatrick, C.K., Clary, C.W., Cyr, A.J., Maletsky, L.P., Rullkoetter, P.J. and Laz, P.J., 2015. Statistical modeling to characterize relationships between knee anatomy and kinematics. Journal of Orthopaedic Research, 33(11), pp.1620-1630.
  * [17] Bossa, M.N. and Olmos, S., 2007, April. Multi-object statistical pose+ shape models. In 2007 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro (pp. 1204-1207). IEEE.
  * [18] Klinder, T., Wolz, R., Lorenz, C., Franz, A. and Ostermann, J., 2008, September. Spine segmentation using articulated shape models. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 227-234). Springer, Berlin, Heidelberg.
  * [19] Fletcher, P.T., Lu, C., Pizer, S.M. and Joshi, S., 2004. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE transactions on medical imaging, 23(8), pp.995-1005.
  * [20] Bossa, M.N. and Olmos, S., 2006, June. Statistical model of similarity transformations: Building a multi-object pose. In 2006 Conference on Computer Vision and Pattern Recognition Workshop (CVPRW’06) (pp. 59-59). IEEE.
  * [21] Chen, X., Graham, J., Hutchinson, C. and Muir, L., 2014. Automatic generation of statistical pose and shape models for articulated joints. IEEE transactions on medical imaging, 33(2), pp.372-383.
  * [22] Ambellan, F., Zachow, S. and von Tycowicz, C., 2019, October. A surface-theoretic approach for statistical shape modeling. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 21-29). Springer, Cham.
  * [23] Styner, M., Gorczowski, K., Fletcher, T., Jeong, J.Y., Pizer, S.M. and Gerig, G., 2006, August. Statistics of pose and shape in multi-object complexes using principal geodesic analysis. In International Workshop on Medical Imaging and Virtual Reality (pp. 1-8). Springer, Berlin, Heidelberg.
  * [24] Madsen, D., Morel-Forster, A., Kahr, P., Rahbani, D., Vetter, T. and Lüthi, M., 2019. A closest point proposal for MCMC-based probabilistic surface registration. arXiv preprint arXiv:1907.01414.
  * [25] Berlinet, A., Thomas-Agnan, Reproducing Kernel Hilbert Spaces in Probability and Statistics, vol. 3. Berlin, Germany: Springer, 2004.
  * [26] Rasmussen, C.E., 2003, February. Gaussian processes in machine learning. In Summer School on Machine Learning (pp. 63-71). Springer, Berlin, Heidelberg.
  * [27] Fouefack, J.R., Alemneh, T., Borotikar, B., B., Burdin, V.,Douglas, T and Mutsvangwa, T., 2019, July. Statistical shape-kinematics models of the skeletal joints: Application tothe shoulder complex. In 2019 41 th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE.
  * [28] Jacq, J.-J., Cresson, T., Burdin, V. and Roux, C. 2008, ‘Performing accurate joint kine-matics from 3-d in vivo image sequences through consensus-driven simultaneous regis-tration’,IEEE Transactions on Biomedical Engineering55(5), 1620–1633
  * [29] Besl, P. J. and McKay, N. D. 1992, Method for registration of 3-d shapes,in‘SensorFusion IV: Control Paradigms and Data Structures’, Vol. 1611, International Society for Optics and Photonics, pp. 586–607
  * [30] Van De Giessen, M., Streekstra, G. J., Strackee, S. D., Maas, M., Grimbergen, K. A.,Van Vliet, L. J. and Vos, F. M. 2009, ‘Constrained registration of the wrist joint’,IEEE Transactions on Medical Imaging28(12), 1861–1869.
  * [31] Gee, A.H. and Treece, G.M., 2014. Systematic misregistration and the statistical analysis of surface data. Medical image analysis, 18(2), pp.385-393.
  * [32] Eberly, D. H. 2000. 3D Game Engine Design: A Practical Approach to Real-time Computer Graphics. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA
  * [33] Cootes, T.F., Edwards, G.J. and Taylor, C.J., 2001. Active appearance models. IEEE Transactions on Pattern Analysis and Machine Intelligence, (6), pp.681-685.