Manifold Learning in Medical Imaging

Abstract

Manifold learning theory has seen a surge of interest in the modeling of large and extensive datasets in medical imaging since they capture the essence of data in a way that fundamentally outperforms linear methodologies, the purpose of which is to essentially describe things that are flat. This problematic is particularly relevant with medical imaging data, where linear techniques are frequently unsuitable for capturing variations in anatomical structures. In many cases, there is enough structure in the data (CT, MRI, ultrasound) so a lower dimensional object can describe the degrees of freedom, such as in a manifold structure. Still, complex, multivariate distributions tend to demonstrate highly variable structural topologies that are impossible to capture with a single manifold learning algorithm. This chapter will present recent techniques developed in manifold theory for medical imaging analysis, to allow for statistical organ shape modeling, image segmentation and registration from the concept of navigation of manifolds, classification, as well as disease prediction models based on discriminant manifolds. We will present the theoretical basis of these works, with illustrative results on their applications from various organs and pathologies, including neurodegenerative diseases and spinal deformities.

1. Introduction

Learning on large medical imaging datasets is an emerging discipline driven from the availability of vast amounts of raw data in many of today’s biomedical studies. However, challenges such as unbalanced data distributions, complex multivariate data and highly variable structural topologies demonstrated by real-world samples makes it much more difficult to efficiently learn the associated representation. An important goal of scientific data analysis in medicine, particularly in neurosciences or oncology, is to understand the behavior of biological process or physiological/morphological alterations. This introduces the need to synthesize large amounts of multivariate data in a robust manner and raises the fundamental question of data reduction: how to discover meaningful representations from unstructured high-dimensional medical images.

Several approaches have attempted to understand how dimension reduction and regression establishes the relationship in subspaces and finally determine statistics on manifolds that optimally describe the relationships between the samples [1]. However, certain assumptions based on the representation of shapes and images using smooth manifolds are made in most cases, which frequently will not be adequate in the presence of medical imaging data and often perturbed by nuisance articulations, clutter or varying contrast.

High-dimensional classification methods have shown promise to measure subtle and spatially complex imaging patterns that have diagnostic value [2, 3]. Defining statistics on a manifold is not a straightforward process when simple statistics cannot be directly applied to general manifolds [4]. But while Euclidean estimators have been used for vector spaces, none have been adapted for multimodal data lying in different spaces. Still, there has been interest in the characterization of data in a Riemann space [5, 6]. Unfortunately, manifold-valued metrics based on the centrality theory or the geometric median [7] often lacks robustness to outliers.

A related topic lies in dimensionally reduced growth trajectories of various anatomical sites which have been investigated in neurodevelopment of newborns for example, based on geodesic shape regression to compute the diffeomorphisms with image time series of a population [8]. These regression models were also used to estimate spatiotemporal evolution of the cerebral cortex [9]. The concept of parallel transport curves in the tangent space from low-dimensional manifolds proposed by Schiratti et al. [10] was used to analyze shape morphology [11] and adapted for radiotherapy response [12]. Regression models were proposed for both cortical and subcortical structures, with 4D varifold-based learning framework with local topography shape morphing being proposed by Rekik et al. [13].

This chapter presents several manifold learning methodologies designed to address challenges encountered in medical imaging. In Section 2, we present an articulated shape inference model from nonlinear embeddings, expressing the global and local shape variations of the spine and vertebrae composing it, introduced in [14]. We then present in Section 3 a probabilistic model from discriminant manifolds to classify the neurodegenerative stage of Alzheimer’s disease. Finally, a piecewise-geodesic transport curve in the tangent space from low-dimensional manifolds designed for the prediction of correction in spinal surgeries is shown in Section 4, introducing a time-warping function controlling the rate of shape evolution. We conclude this article in Section 5.

2. Shape inference through navigating manifolds

Statistical models of shape variability have been successful in addressing fundamental vision tasks such as segmentation and registration in medical imaging. However, the high dimensionality and complex nonlinear underlying structure unfortunately makes the commonly used linear statistics inapplicable for anatomical structures. Manifold learning approaches map high-dimensional observation data that are presumed to lie on a nonlinear manifold, onto a single global coordinate system of lower dimensionality.

Inferring a model from the underlying manifold is not a novel concept but far from being trivial. In this section, we model both global statistics of the articulated model and local shape variations of vertebrae based on local measures in manifold space. We describe a spine inference/segmentation method from CT and MR images, where the model representation is optimized through a Markov Random Field (MRF) graph, balancing prior distribution with image data.

2.1. Data representation

Our spine model $\mathbf{S}=\left\{{\mathbf{s}}_1,\dots ,{\mathbf{s}}_L\right\}$ consists of an interconnection of $L$ vertebrae. For each vertebra ${\mathbf{s}}_i$, we recover a triangular mesh with vertices $\left\{{\mathbf{v}}_j^i|j=1,\dots ,V\right\}$, where the $j^{\mathit{th}}$ vertex corresponds to approximately the same location from one shape to another and $V$ the number of vertices. Additionally, every ${\mathbf{s}}_i$ is annotated with landmarks on each model to rigidly register each object to its upper neighbor. Hence, an articulated deformable model (ADM) is represented by a vector of local intervertebral rigid transformations $A=\left[T_1,T_2,\dots ,T_L\right]$. To perform global shape modeling of $\mathbf{S}$, we convert $A$ to an absolute representation ${\mathbf{A}}_{\mathrm{abs}}=\left[T_1,T_1\circ T_2,\dots ,T_1\circ T_2\circ \dots \circ T_L\right]$ using recursive compositions. The transformations are expressed in the local coordinate system (LCS) of the lower vertebra. Center of transformation is the intersection of all three vertebral axes, following anteroposterior, cranial-caudal and left-right directions. Rigid transformations described here are the combination of a rotation matrix $R$, a translation $t$ and scaling $s$. We formulate the rigid transformation $T=\left\{s,R,t\right\}$ of a triangular mesh model as $y=\mathit{sRx}+t$ where $x,y,t\in {\mathcal{\Re }}^3$.

2.2. Manifold embedding

For nonlinear embeddings, we rely on the absolute vector representation ${\mathbf{A}}_{\mathrm{abs}}$ as given previously. Let us now consider $N$ articulated shape models expressed by the feature vectors ${\mathbf{A}}_{\mathrm{abs}}^i$, of dimensionality $D$. The aim is to create a low-dimensional manifold consisting of $N$ points ${\mathbf{Y}}_i$, ${\mathbf{Y}}_i\in {\Re }^d$, $i\in \left[1,N\right]$ where $d\ll D$ based on [15]. In such a framework, if an adequate number of data points is available, then the underlying manifold $\mathcal{M}$ is considered to be “well-sampled.” Therefore, it can represent the underlying population structure. In the sub-cluster corresponding to a pathological population, each point of the training set and its neighbors would lie within a locally linear patch as illustrated in Figure 1.

Figure 1.

Representation of intervertebral transformations in manifold space.

The main limitation of embedding algorithms is the assumption of Euclidean metrics in the ambient space to evaluate similarity between sample points. Thus, a metric in the space of articulated structures is defined so that it accommodates for anatomical spine variability and adopts the intrinsic nature of the Riemannian manifold geometry allowing us to discern between articulated shape deformations in a topological invariant framework. For each point, the $K$ closest neighbors are selected using a distortion metric which is particularly suited for geodesics. The metric $d_{\mathrm{M}}\left({\mathbf{A}}_{\mathrm{abs}}^i,{\mathbf{A}}_{\mathrm{abs}}^j\right)$ estimates the distance of articulated models $i,j$ where ${\mathbf{A}}_{\mathrm{abs}}^i$. The distance measure for absolute representations can therefore be expressed as a sum of articulation deviations

$$d_M\left({\mathbf{A}}_{\mathrm{abs}}^i,{\mathbf{A}}_{\mathrm{abs}}^j\right)=\sum _{k=1}^L{d}_M\left(T_k^i,T_k^j\right)=\sum _{k=1}^L\parallel {\mathbf{t}}_k^i-{\mathbf{t}}_k^j\parallel +\sum _{k=1}^L{d}_G\left({\mathbf{R}}_k^i,{\mathbf{R}}_k^j\right).$$E1

While for the translation, the $L_2$ norm is chosen, geodesical distances are used between rotation neighborhoods. This is expressed as $d_G\left({\mathbf{R}}_k^i,{\mathbf{R}}_k^j\right)=\parallel log\left({\left({\mathbf{R}}_k^i\right)}^{-1}{\mathbf{R}}_k^j\right){\parallel }_F$ where the log map is used to map a point in the manifold to the tangent plane.

Afterwards, the manifold reconstruction weights are estimated by assuming the local geometry of the patches can be described by linear coefficients that permit the reconstruction of every model point from its neighbors. In order to determine the value of the weights, the reconstruction errors are measured using the following objective function:

$$\epsilon \left(\mathbf{W}\right)=\sum _{i=1}^N{‖{\mathbf{A}}_{\mathrm{abs}}^i-\sum _{j=1}^K{W}_{\mathit{ij}}{\mathbf{A}}_{\mathrm{abs}}^j‖}^2$$E2

$$\text{subject}\text{to}\left\{\begin{array}{l}{W}_{\mathit{ij}}=0\text{if}{\mathbf{A}}_{\mathrm{abs}}^i\text{not}\text{neighbor}{\mathbf{A}}_{\mathrm{abs}}^j\\ \sum _j{\mathbf{W}}_{\mathit{ij}}=1\text{for}\text{every}i.\end{array}\right.$$E3

Thus, $\epsilon \left(\mathbf{W}\right)$ sums the squared distances between all data points and their corresponding reconstructed points. The weights ${\mathbf{W}}_{\mathit{ij}}$ represent the importance of the $j^{\mathit{th}}$ data point to the reconstruction of the $i^{\mathit{th}}$ element.

The algorithm maps each high-dimensional ${\mathbf{A}}_{\mathrm{abs}}^i$ to a low-dimensional ${\mathbf{Y}}_i$. These internal coordinates are found with a cost function minimizing the reconstruction error:

$$\begin{array}{c}\mathrm{\Phi }\left(\mathbf{Y}\right)=\sum _{i=1}^N{‖{\mathbf{Y}}_i-\sum _{j=1}^K{\mathbf{W}}_{\mathit{ij}}{\mathbf{Y}}_j‖}^2\\ =\sum _{i=1}^N\sum _{j=1}^N{\mathbf{M}}_{\mathit{ij}}{\mathbf{Y}}_i^T{\mathbf{Y}}_j\end{array}$$E4

with $\mathbf{M}$ as a sparse and symmetric $N\times N$ matrix enclosing the reconstruction weights ${\mathbf{W}}_{\mathit{ij}}$ such that $\mathbf{M}={\left(\mathbf{I}-\mathbf{W}\right)}^T\left(\mathbf{I}-\mathbf{W}\right)$, and $\mathbf{Y}$ spanning the ${\mathbf{Y}}_i$’s. The optimal embedding, up to a global rotation, is obtained from the bottom $d+1$ eigenvectors of $\mathbf{M}$ and helps to minimize the cost function $\Phi \left(\mathbf{Y}\right)$ as a simple eigenvalue problem. The $d$ eigenvectors form the $d$ embedding coordinates. The coordinates ${\mathbf{Y}}_i$ can be translated by a constant displacement without affecting the overall cost $\Phi \left(\mathbf{Y}\right)$. The eigenvector corresponding to the smallest eigenvalue corresponds to the mean value of the embedded data ${\mathbf{Y}}^0=\left\{{\mathbf{y}}_1,\dots ,{\mathbf{y}}_d\right\},{\mathbf{y}}_i=0,\forall i$. This can be discarded with $\sum {\mathbf{Y}}_i=0$ to obtain an embedding centered at the origin. Hence, a new ADM can be inferred in the embedded $d$-space as a low-dimensional point ${\mathbf{Y}}^{\mathrm{new}}$ by finding its optimal manifold coordinates ${\mathbf{y}}_i$.

To obtain the articulation vector for a new embedded point in the ambient space (image domain), one has to determine the representation in high-dimensional space based on its intrinsic coordinates. We first assume an explicit mapping $f:\mathcal{M}\to {\Re }^D$ from manifold space $\mathcal{M}$ to the ambient space ${\Re }^D$. The inverse mapping of ${\mathbf{Y}}_i$ is then performed by estimating the relationship between ${\Re }^D$ and $\mathcal{M}$ as a joint distribution, such there exists a smooth functional which belongs to a local neighborhood. Theoretically the manifold should follow the conditional expectation:

$$f\left(Y_i\right)\equiv E\left({\mathbf{A}}_{\mathrm{abs}}^i|\mathcal{M}\left(A_i\right)={\mathbf{Y}}_i\right)=\int A_i\frac{p\left({\mathbf{Y}}_i,A_i\right)}{p_{\mathcal{M}\left(A_i\right)}\left({\mathbf{Y}}_i\right)}\mathit{dD}$$E5

which captures the overall trend of the data in $D$-space. Here, both $p_{\mathcal{M}\left(A_i\right)}\left({\mathbf{Y}}_i\right)$ (marginal density of $\mathcal{M}\left(A_i\right)$) and $p\left({\mathbf{Y}}_i,A_i\right)$ (joint density) are unknown. Based on the Nadaraya-Watson kernel regression [16], we replace densities by kernel functions as $p_{\mathcal{M}\left(A_i\right)}\left({\mathbf{Y}}_i\right)=\frac{1}{K}{\sum }_{j\in \mathcal{N}\left(i\right)}{G}_h\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)$ and $p\left({\mathbf{Y}}_i,A_i\right)=\frac{1}{K}{\sum }_{j\in \mathcal{N}\left(i\right)}{G}_h\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)G_g\left(A_i,A_j\right)$ [17]. The Gaussian regression kernels $G$ require the neighbors ${\mathbf{A}}_{\mathrm{abs}}^j$ of $j\in \mathcal{N}\left(i\right)$ to determine the bandwidths $h,g$ so it includes all $K$ data points ($\mathcal{N}\left(i\right)$ representing the neighborhood of $i$). Plugging these estimates in Eq.(5), this gives:

$$f_{\mathrm{NW}}\left({\mathbf{Y}}_i\right)=\int A_i\frac{\frac{1}{K}{\sum }_{j\in \mathcal{N}\left(i\right)}{G}_h\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)G_g\left(A_i,A_j\right)}{\frac{1}{K}{\sum }_{j\in \mathcal{N}\left(i\right)}{G}_h\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)}\mathit{dD}.$$E6

By assuming $G$ is symmetric about the origin, we propose to integrate in the kernel regression estimator, the manifold-based distortion metric $d_M$ which is particularly suited for geodesic metrics and articulated diffeomorphisms. This generalizes the expectation such that the observations $\mathbf{Y}$ are defined in manifold space $\mathcal{M}$:

$$f_{\mathrm{NW}}\left({\mathbf{Y}}_i\right)=\underset{{\mathbf{A}}_{\mathrm{abs}}^i}{\text{argmin}}\frac{{\sum }_{j\in \mathcal{N}\left(i\right)}G\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)d_M\left({\mathbf{A}}_{\mathrm{abs}}^i,{\mathbf{A}}_{\mathrm{abs}}^j\right)}{{\sum }_{j\in \mathcal{N}\left(i\right)}G\left({\mathbf{Y}}_i,{\mathbf{Y}}_j\right)}$$E7

which integrates the distance metric $d_M\left({\mathbf{A}}_{\mathrm{abs}}^i,{\mathbf{A}}_{\mathrm{abs}}^j\right)$ defined in Eq. (1) and updates $f_{\mathrm{NW}}\left({\mathbf{Y}}_i\right)$ using the closest neighbors of point ${\mathbf{Y}}_i$ in the manifold space. This constrains the regression to be valid for similar data points in its vicinity since locality around ${\mathbf{Y}}_i$ preserves locality in ${\mathbf{A}}_{\mathrm{abs}}^i$.

2.3. Optimization on manifold

Once an appropriate modeling of spine shape variations is determined with a manifold, a successful inference between the image and manifold must be accomplished. We describe here how a new model is generated. We search the optimal embedded manifold point $\mathbf{Y}=\left({\mathbf{y}}_1,\dots ,{\mathbf{y}}_d\right)$ of the global spine model. Such a strategy offers an ideal compromise between the prior constraints, as well as the individual shape variations described by the weight vector $\mathbf{W}=\left({\mathbf{w}}_1,\dots ,{\mathbf{w}}_n\right)$ in a localized sub-patch. The energy $E$ of inferring the model $\mathbf{S}$ in the image $\mathcal{I}$ is a function of the set of displacement vectors $\mathrm{\Delta }$ in the manifold space for global shape representation. This involves: (a) a data-related term expressing the image cost and (b) a global prior term measuring deformation between low-dimensional vectors with shape models. The third term represents (c) a higher-order term which is expressed by the reconstruction weights $\Omega$ for local vertebra modeling. The energy $E$ can be expressed as the following combination of a global and local optimization:

$$E\left({\mathbf{S}}^0,\mathcal{I},\mathrm{\Delta },\mathrm{\Omega }\right)=V\left({\mathbf{Y}}^0+\mathrm{\Delta },\mathcal{I}\right)+\alpha V\left(\mathbf{N},\mathrm{\Delta }\right)+\beta V\left(\mathbf{H},\mathrm{\Delta },\mathrm{\Omega }\right).$$E8

The global alignment of the model with the target image primarily drives the deformation of the model. The purpose is to estimate the set of articulations describing the global spine model by determining its optimal representation ${\mathbf{Y}}^0$ in the embedded space. This is performed by obtaining the global representation using the mapping in (7) so that: $f_{\mathrm{NW}}\left(Y_i+\mathrm{\Delta }\right)=f_{\mathrm{NW}}\left(\left\{y_1+{\mathit{\delta }}_1,\dots ,y_d+{\mathit{\delta }}_d\right\}\right)$. This allows to optimize the model in manifold space coordinates while retrieving the articulations in $\mathcal{I}$. The global cost can be expressed as:

$$V\left({\mathbf{Y}}^0+\mathrm{\Delta },\mathcal{I}\right)=V\left(f_{\mathrm{NW}}\left(\left\{y_1+{\mathit{\delta }}_1,\dots ,y_d+{\delta }_d\right\}\right),\mathcal{I}\right)).$$E9

The inverse transform allows to obtain ${\mathbf{A}}_{\mathrm{abs}}^i+\mathit{D}$, with $\mathit{D}$ as deformations in the image space. Since the transformations $T_i$ are implicitly modeled in the absolute representation ${\mathbf{A}}_{\mathrm{abs}}^0$, we can formally consider the singleton image-related term as a summation of costs associated with each $L$ vertebra of the model:

$$V\left({\mathbf{A}}_{\mathrm{abs}}^0+\mathit{D},\mathcal{I}\right)=\sum _{i=1}^L{V}_i\left({{\mathbf{s}}_i}^{\ast }\left(T_i^0+d_i\right),\mathcal{I}\right)$$E10

where $V_i\left(\mathbf{s},\mathcal{I}\right)={\sum }_{{\mathbf{v}}_i\in \mathbf{s}}{\mathbf{n}}_i^T\left({\mathbf{v}}_i\right)\nabla \mathcal{I}\left({\mathbf{v}}_i\right)$ minimizes the distance between mesh vertices of the inferred shape and gradient image $\mathcal{I}$ by a rigid transformation. Here, ${\mathbf{n}}_i$ is the normal pointing outwards and $\nabla \mathcal{I}\left({\mathbf{v}}_i\right)$ the image gradient at ${\mathbf{v}}_i$.

The prior constraint for the rigid alignment are pairwise potentials between neighboring models $y_i$ such that the difference in manifold coordinates is minimal with regards to a prior distribution of neighboring distances $P$:

$$\alpha V\left(\mathbf{N},\mathrm{\Delta }\right)=\alpha \sum _{i\in G}\sum _{j\in \mathcal{N}\left(i\right)}{V}_{\mathit{ij}}\left(y_i^0+{\mathit{\delta }}_i,y_j^0+{\mathit{\delta }}_j,P\right).$$E11

This term represents the smoothness term of the global cost function to ensure that the deformation ${\mathit{\delta }}_i$ applied to point coordinates are regular, with $V_{\mathit{ij}}=\left(0,1\right)$ a distance assigning function based on the distances to $P$.

One can integrate the global data and prior terms along with local shape terms parameterized as the higher-order cliques, by combining (9), (11):

$$\begin{array}{cc}E\left({\mathbf{S}}^0,\mathcal{I},\mathit{\Delta },\mathrm{\Omega }\right)& =V\left(f_{\mathrm{NW}}\left(\left\{y_1+{\mathit{\delta }}_1,\dots ,y_d+{\mathit{\delta }}_{\mathbf{d}}\right\}\right),\mathcal{I})\right)\\ & +\alpha \sum _{i\in G}\sum _{j\in \mathcal{N}\left(i\right)}{V}_{\mathit{ij}}\left(y_i^0+{\mathit{\delta }}_i,y_j^0+{\mathit{\delta }}_j\right)+\beta \sum _{c\in \mathcal{C}}{V}_c\left({\mathbf{w}}_c^0+{\omega }_c\right).\end{array}$$E12

The optimization strategy of the resulting MRF (12) in the continuous domain is not a straightforward problem. The convexity of the solution domain is not guaranteed, while gradient-descent optimization approaches are prone to nonlinearity and local minimums. We seek to assign the optimal labels ${\mathcal{L}}^{\Delta }=\left\{l_1,\dots ,l_d\right\}$ and ${\mathcal{L}}^{\Omega }=\left\{l_1,\dots ,l_n\right\}$ which are associated to the quantized space $\Delta$ of displacements and local weight parameters $\Omega$ respectively. We consider that displacing the coordinates of point $y_i^0$ by ${\mathit{\delta }}^{l_i}$ is equivalent to assigning label $l_i$ to $y_i^0$. An incremental approach is adopted where in each iteration $t$ we look for the set of labels that improves the current solution s.t. $y_i^t=y_i^0+{\sum }_t{\mathit{\delta }}^{l_{i}t}$, which is a temporal minimization problem. Then (12) can be rewritten as:

$$\begin{array}{cc}{E}^t\left({\mathcal{L}}^{\mathit{\Delta }}{\mathcal{L}}^{\mathrm{\Omega }}\right)& =V\left(f_{\mathrm{NW}}\left(\left\{y_1^{t-1},l_1^{\Delta },\dots ,y_d^{t-1},l_d^{\Delta }\right\}\right),\mathcal{I})\right)\\ & +\alpha \sum _{i\in G}\sum _{j\in \mathcal{N}\left(i\right)}{V}_{\mathit{ij}}\left(y_i^{t-1},y_j^{t-1},l_i^{\Delta },l_j^{\Delta }\right)+\beta \sum _{c\in \mathcal{C}}{V}_c\left({\mathbf{w}}_c^{t-1},l_c^{\mathrm{\Omega }}\right).\end{array}$$E13

We solve the minimization of the higher-order cliques in (13) by transforming them into quadratic functions [18]. We apply the FastPD method [19] which solves the problem by formulating the duality theory in linear programming.

2.4. Results

Manifold learning. The manifold was built from a database containing $711$ scoliotic spines demonstrating several types of deformities. Each spine model in the database was obtained from biplanar radiographic stereo-reconstructions. It is modeled with 12 thoracic and 5 lumbar vertebrae (17 in total), represented by 6 landmarks on each vertebra (4 pedicle extremities and 2 endplate center points) which were manually identified by an expert on the radiographic images. The resulting manifold is shown in Figure 2.

Figure 2.

Low-dimensional manifold embedding of the spine dataset comprising 711 models exhibiting various types of deformities. The sub-domain was used to estimate both the global shape pose costs and individual shape instances based on local neighborhoods.

Adaptation of the articulated model was done on two different data sets. The first consisted of volumetric CT scans ($512\times 512\times 251$, resolution: $0.8\times 0.8$ mm, thickness: $1-2$ mm) of the lumbar and main thoracic regions obtained from 21 different patients acquired for operative planning purposes. The MR dataset comprised multi-parametric volumetric data ($256\times 256\times 160$, resolution: $1.3\times 0.9$ mm, thickness: $1$ mm) of 8 patients acquired for diagnostic purposes. For this study, only the T1 sequence was selected for the experiments. All patients on both datasets (29 in total) had 12 thoracic and 5 lumbar vertebrae. Both CT and MR data were manually annotated with 3D landmarks by an expert in radiology, corresponding to left and right pedicle tips as well as midpoints of the vertebral body. Segmentation of the vertebrae from the CT and MR slices were also made by the same operator.

CT imaging experiments. We first evaluated the model accuracy in CT images by computing the correspondence of the inferred vertebral mesh models to the segmented target structures. As a preprocessing step, a rough thresholding was performed on the whole volume to filter out noise artifacts. The overall surface-to-surface comparison results between the inferred 3D vertebral models issued from the articulated model and from known segmentations were first calculated. The mean errors are $2.2\pm 1.5$ mm (range: $0.6-5.4$ mm) for thoracic vertebra and $2.8\pm 1.9$ mm (range: $0.7-8.1$ mm) for lumbar vertebra.

MR imaging experiments. For the experiments involving the segmentation of 3D spine models from MR images, the surface-to-surface comparison showed encouraging results (thoracic: $2.9\pm 1.8$ mm, lumbar: $3.0\pm 1.9$ mm) based on differences to ground-truth. As in the previous experiments with CT imaging, ground-truth data was generated by manually segmenting the structures models which were validated by an expert in radiology. As difficult as the CT inference is, the MR problem represent an even greater challenge as the image resolution is more limited and interslice spacing is increased compared to CT. Modeling of the statistical properties of the shape variations and global pose becomes even more important in this case, as it relies heavily in the nonlinear distribution of the patient morphology.

3. Probabilistic modeling of discriminant nonlinear manifolds in the identification of Alzheimer’s

Neurodegenerative pathologies, such as Alzheimer’s disease (AD), are linked with morphological and metabolic alterations which can be assessed from medical imaging and biological data. Recent advances in machine learning have helped to improve classification and prognosis rates, but lack a probabilistic framework to measure uncertainty in the data. In this section, we present a method to identify progressive mild cognitive impairment (MCI) and predict their conversion to AD from MRI and positron emitting tomography (PET) images. We show a discriminative probabilistic manifold embedding where locally linear mappings transform data points in low-dimensional space to corresponding points in high-dimensional space. A discriminant adjacency matrix is constructed to maximize the separation between different clinical groups, including MCI converters and nonconverters, while minimizing the distance in latent variables belonging to the same class.

3.1. Probabilistic model for discriminant manifolds

Manifold learning algorithms are based on the premise that data are often of artificially high dimension and can be embedded in a lower dimensional space. However the presence of outliers and multiclass information can on the other hand affect the discrimination and/or generalization ability of the manifold. We propose to learn the optimal separation between four classes (1) normal controls, (2) nonconverter MCI patients, (3) converter MCI patients and (4) AD patients, by using a discriminant graph-embedding. Here, $n$ labeled points $\mathbb{Y}={\left\{\left(y_i,l_i\right)\right\}}_{i=1}^n$ defined in ${\mathbb{R}}^D$ are generated from the underlying manifold $\mathcal{M}$, where $l_i$ denotes the label (NC, cMCI, nMCI or AD). For the labeled data, there exists a low-dimensional (latent) representation of the high-dimensional samples such that $\mathbb{X}={\left\{\left({\mathbf{x}}_i,l_i\right)\right\}}_{i=1}^n$ defined in ${\mathbb{R}}^d$. We assume here that the mapping ${\mathbf{M}}_i\in {\mathbb{R}}^{D\times d}$ between high and low-dimensional spaces is locally linear, such that tangent spaces in local neighborhoods can be estimated with ${\mathbf{y}}_j-{\mathbf{y}}_i$ and ${\mathbf{x}}_j-{\mathbf{x}}_i$, representing the pairwise differences between connected neighbors $i,j$. Therefore the relationship can be established as ${\mathbf{y}}_j-{\mathbf{y}}_i\approx {\mathbf{M}}_i\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)$.

In order to effectively discover the low-dimensional embedding, it is necessary to maintain the local structure of the data in the new embedding. The graph $G=\left(\mathit{V},\mathit{W}\right)$ is an undirected similarity graph, with a collection of nodes $\mathit{V}$ connected by edges, and the symmetric matrix $\mathit{W}$ with elements describing the relationships between the nodes. The diagonal matrix $\mathit{D}$ and the Laplacian matrix $\mathit{L}$ are defined as $\mathit{L}=\mathit{D}-\mathit{W}$, with $\mathit{D}\left(i,i\right)={\sum }_{j\ne i}{\mathit{W}}_{\mathit{ij}}\forall i$.

Using the theoretical framework from [20], we can determine a distribution of linear maps associated with the low-dimensional representation to describe the data likelihood for a specific model:

$$logp\left(\mathbb{Y}|G\right)=log\iint p\left(\mathbb{Y},\mathbf{M},\mathbb{X}|G\right)d\mathbf{x}d\mathbf{M}$$E14

This joint distribution can be separated into three prior terms: the linear maps, latent variables and the likelihood of the high dimensional points $\mathbb{Y}$:

$$p\left(\mathbb{Y},\mathbf{M},\mathbb{X}|G\right)=p\left(\mathbb{Y}|\mathbf{M},\mathbb{X},G\right)p\left(\mathbf{M}|G\right)p\left(\mathbb{X}|G\right)$$E15

We now define the discriminant similarity graphs establishing neighborhood relationships, as well define each of the three prior terms included in the joint distribution.

Within and between similarity graphs: In our work, the geometrical structure of $\mathcal{M}$ can be modeled by building a within-class similarity graph ${\mathit{W}}_w$ for feature vectors of same group and a between-class similarity graph ${\mathit{W}}_b$, to separate features from all four classes. When constructing the discriminant locally linear latent variable embedding, elements are partitioned into ${\mathit{W}}_w$ and ${\mathit{W}}_b$ classes. The intrinsic graph $G$ is first created by assigning edges only to samples of the same class (ex: nMCI). Each sample is therefore reconstructed only from feature vectors of the same clinical group. Local reconstruction coefficients are incorporated in the within-class similarity graph, such that ${\mathit{W}}_w$ is defined as:

$$W_{w_{i,j}}=\left\{\begin{array}{ll}1& \text{if}{\mathbf{y}}_i\in {\mathcal{N}}_w\left({\mathbf{y}}_j\right)\text{or}{\mathbf{y}}_j\in {\mathcal{N}}_w\left({\mathbf{y}}_i\right)\\ 0,& \text{otherwise}.\end{array}\right)$$E16

with ${\mathcal{N}}_w$ containing neighbors of the same class. Conversely, ${\mathit{W}}_b$ depicts the statistical properties to be avoided in the inference process. Distances between samples from different clinical groups are computed as:

$$W_{b_{i,j}}=\left\{\begin{array}{ll}1& \text{if}{\mathbf{y}}_i\in {\mathcal{N}}_b\left({\mathbf{y}}_j\right)\text{or}{\mathbf{y}}_j\in {\mathcal{N}}_b\left({\mathbf{y}}_i\right)\\ 0,& \text{otherwise}\end{array}\right)$$E17

with ${\mathcal{N}}_b$ containing neighbors having different class labels from the $i$th sample. The objective is to transform points to a new manifold $\mathcal{M}$ of dimensionality $d$, i.e., ${\mathbf{y}}_i\to {\mathbf{x}}_i$, by mapping connected samples from the same group in ${\mathit{W}}_w$ as close as possible to the class cluster, while moving NC, nMCI, cMCI and AD samples of ${\mathit{W}}_b$ as far away from one another.

Model components: The prior added on the latent variables $\mathbb{X}$ are located at the origin of the low-dimensional domain, while minimizing the Euclidean distance of neighboring points that are associated with the neighborhood of high-dimensional points and maximizing the distance between coordinates of different classes. In order to set the variables with an expected scale $\alpha$ and $H$ representing the probability density function, the following log prior is defined:

$$logp\left(\mathbb{X}|\mathit{W},\alpha \right)=-\frac{1}{2}\sum _{i=1}^n\left(\alpha \parallel {\mathbf{x}}_i\parallel +\sum _{j=1}^n{W}_{w_{i,j}}\parallel {\mathbf{y}}_i-{\mathbf{y}}_j{\parallel }^2-\sum _{j=1}^n{W}_{b_{i,j}}\parallel {\mathbf{y}}_i-{\mathbf{y}}_j{\parallel }^2\right)-log{H}_{\mathbb{X}}$$E18

The prior added to the linear maps defines how the tangent planes described in low and high dimensional spaces are similar based on the Frobenius norm. This prior ensures smooth manifolds:

$$logp\left(\mathbf{M}|\mathit{W}\right)=-\frac{1}{2}\left({‖\sum _{i=1}^n{\mathbf{x}}_i‖}_F^2-\sum _{i=1}^n\sum _{j=1}^n\left(W_{w_{i,j}}-W_{b_{i,j}}\right)\parallel {\mathbf{M}}_i-{\mathbf{M}}_j{\parallel }_F^2\right)-log{H}_{\mathbf{M}}$$E19

Finally, approximation errors from the linear mapping ${\mathbf{M}}_i$ between low and high-dimensional domains are penalized by including the following log likelihood:

$$\begin{array}{ll}logp\left(\mathbb{Y}|\mathbb{X},\mathbf{W},\gamma \right)=\parallel \sum _{i=1}^n{\mathbf{y}}_i{\parallel }^2& -\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{W}_{w_{i,j}}\Delta {\left(i,j\right)}^{\mathrm{T}}\gamma \mathbf{I}\Delta \left(i,j\right)\\ & +\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{W}_{b_{i,j}}\Delta {\left(i,j\right)}^{\mathrm{T}}\gamma \mathbf{I}\Delta \left(i,j\right)-log{H}_{\mathbf{y}}\end{array}$$E20

with $\Delta \left(i,j\right)$ the difference in Euclidean distance between pairs of neighbors in high and low-dimensional space and $\gamma$ the update parameter for the EM inference. Samples of $\mathbf{y}$ are drawn from a multivariate normal distribution.

3.2. Variational inference

The objective is to infer the low-dimensional coordinates and linear mapping function for the described model, as well as the intrinsic parameters of the model $\mathrm{\Phi }=\left(\alpha ,\gamma \right)$. This is achieved by maximizing the marginal likelihood of:

$$logp\left(\mathbb{Y}|\mathbf{W},\mathrm{\Phi }\right)=\iint \rho \left(\mathbf{M},\mathbb{X}\right)log\frac{p\left(\mathbb{Y},\mathbf{M},\mathbb{X}|\mathbf{W},\mathrm{\Phi }\right)}{\rho \left(\mathbf{M},\mathbb{X}\right)}d\mathbf{x}d\mathbf{M}.$$E21

By assuming the posterior $\rho \left(\mathbf{M},\mathbb{X}\right)$ can be factored in separate terms $\rho \left(\mathbf{M}\right)$ and $\rho \left(\mathbb{X}\right)$, a variational expectation maximization algorithm can be used to determine the model’s parameters, which are initialized with $\Phi$. The E-step updates the independent posteriors $\rho \left(\mathbb{X}\right)$ and $\rho \left(\mathbf{M}\right)$, while the parameters of $\Phi$ are updated in the M-step by maximizing Eq. (21).

The discriminant latent variable model can then be used to perform the mapping of new image feature vectors to the manifold. The variational EM algorithm described in the previous section can be used to transform a set of new input points ${\mathbf{y}}_q$ without changing the overall neighborhood graph structure, by finding the distribution of the local linear map ${\mathbf{y}}_q$ and it is low-dimensional coordinate using the E-step explained above. Once the manifold representation ${\mathbf{x}}_q$ is obtained, a cluster analysis finds the corresponding class in the manifold, yielding a prediction of the input feature vector ${\mathbf{y}}_q$.

3.3. Experiments

We used the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database with 1.5 or 3.0 T structural MR images (adni.loni.usc.edu) and FDG-PET images. For this study, 187 subjects with both MRI and PET images during a 24 month period were used to train the probabilistic manifold model, including 46 AD patients, 94 MCI patients, and 47 normal controls. During the follow-up period, 43 MCI subjects converted to AD and 56 remained stable. All groups are matched approximately by age (mean of $76.7\pm 5.4$) and gender. Images were non-rigidly registered to a standard template, which was then segmented using FSL-FIRST automatic segmentation [21].

A 9-fold cross-validation was performed to assess the performance of the method. The optimal manifold dimensionality was set at $d=8$, when the trend of the nonlinear residual reconstruction error stabilized for the entire training set. We evaluated the classification performance of the proposed method for discriminating between cMCI and nMCI patients, by training the model with MRI, PET and with MRI + PET biomarkers from the ROIs illustrated in Figure 3. Figure 4 presents ROC curves obtained by the proposed and comparative methods such as SVM (nonlinear RBF kernel), LLE and LL-LVM [20]. The discriminative nature of the proposed framework clearly shows an improvement to standard learning approaches models which were trained using MRI only, PET only and combined multimodal features. It illustrates that increased accuracy (77.4%) can be achieved by combining MRI and PET features, showing the benefit of extracting complementary features from the dataset for prediction purposes. When comparing the performance of the proposed method to the other learning methods (SVM, LLE, LL-LVM), the probabilistic model integrating similarity graphs shows a statistically significant improvement $\left(p<0.01\right)$ to all three approaches based on paired t-test.

Figure 3.

Selected FSL segmented brain regions for feature selection on (left) MRI and (right) PET images.

Figure 4.

ROC curves comparing the SVM, LLE and LL-LVM with the proposed method for cMCI/nMCI prediction using MRI, PET and multimodality data.

4. Spatiotemporal manifold prediction model for surgery prediction

In this final section, we present a statistical framework for predicting the surgical outcomes following spine surgery of adolescents with idiopathic scoliosis. A discriminant manifold is first constructed to maximize the separation between responsive and nonresponsive groups of patients. The model then uses subject-specific correction trajectories based on articulated transformations in order to map spine correction profiles to a group-average piecewise-geodesic path. Spine correction trajectories are described in a piecewise-geodesic fashion to account for varying times at follow-up exams, regressing the curve via a quadratic optimization process. To predict the evolution of correction, a baseline reconstruction is projected onto the manifold, from which a spatiotemporal regression model is built from parallel transport curves inferred from neighboring exemplars (Figure 5).

Figure 5.

Proposed prediction framework for spine surgery outcomes. In the training phase, a dataset of spine models are embedded in a spatiotemporal manifold $\mathcal{M}$, into responsive (R) or nonresponsive (NR) groups. During testing, an unseen baseline 3D spine reconstruction ${\mathbf{y}}_q$ is projected on $\mathcal{M}$ using $f_{\mathrm{NW}}$ based on Nadaraya-Watson kernels. The closest samples to the projected point $\mathbf{x}$ are selected to regress the spatiotemporal curve $\gamma$ used for predicting the correction due with surgery.

4.1. Discriminant embedding of spine models

We propose to embed a collection of nonresponsive (NR) and (2) responsive (R) patients to surgery which will offer a maximal separation between the classes, by using a discriminant graph-embedding. Here, $n$ labeled points $\mathbb{Y}={\left\{\left({\mathbf{y}}_i,l_i,t_i\right)\right\}}_{i=1}^n$ defined in ${\mathbb{R}}^D$ are embedded in the low-dimensional manifold $\mathcal{M}$, where $l_i$ describes the label (NR or R) and $t_i$ defines the time of follow-up. We assume that for the sampled data, an underlying manifold of the high-dimensional data exists such that $\mathbb{X}={\left\{\left({\mathbf{x}}_i,l_i,t_i\right)\right\}}_{i=1}^n$ defined in ${\mathbb{R}}^d$. We rely on the assumption that a locally linear mapping ${\mathbf{M}}_i\in {\mathbb{R}}^{D\times d}$ exists, where local neighborhoods are defined as tangent planes estimated with ${\mathbf{y}}_j-{\mathbf{y}}_i$ and ${\mathbf{x}}_j-{\mathbf{x}}_i$, describing the paired distances between linked neighbors $i,j$. Hence, the relationship can be established as ${\mathbf{y}}_j-{\mathbf{y}}_i\approx {\mathbf{M}}_i\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)$.

Because the discriminant manifold structure in ${\mathbb{R}}^d$ requires to maintain the local structure of the underlying data, a undirected similarity graph $\mathcal{G}=\left(\mathit{V},\mathit{W}\right)$ is built, where each node $\mathit{V}$ are connected to each other with edges that are weighted with the graph $\mathit{W}$. The overall structure of $\mathcal{M}$ is therefore defined with ${\mathit{W}}_w$ for feature vectors belonging to the same class and ${\mathit{W}}_b$, which separate features from both classes. During the embedding of the discriminant locally linear latent manifold, data samples are divided between ${\mathit{W}}_w$ and ${\mathit{W}}_b$.

4.2. Piecewise-geodesic spatiotemporal manifold

Once sample points ${\mathbf{x}}_i$ are in manifold space, the objective is to regress a regular and smooth piecewise-geodesic curve $\gamma :\left[t_1,t_N\right]$ that accurately fits the embedded data describing the spatiotemporal correction following surgery within a 2 year period. For each sample data ${\mathbf{x}}_i$, the $K$ closest individuals demonstrating similar baseline features are identified from the embedded data, creating neighborhoods $\mathcal{N}\left({\mathbf{x}}_q\right)$ with measurements at different time points, thus creating a low-dimensional Riemannian manifold where data points ${\mathbf{x}}_{i,j}$, with $i$ denoting a particular individual, $j$ the time-point measurement and $j=0$ the preoperative model. By assuming the manifold domain is complete and piecewise-geodesic curves are defined for each time trajectories, time-labeled data can be regressed continuously in ${\mathbb{R}}^D$, thereby creating smooth curves in time intervals described by samples in ${\mathbb{R}}^d$.

However, due to the fact the representation of the continuous curve is a variational problem of infinite dimensional space, the implementation follows a discretization process which is derived from the procedure in [22], such that:

$$\begin{array}{l}E\left(\gamma \right)=\frac{1}{K_d}\sum _{i=1}^{K_d}\sum _{j=0}^{t_N}{w}_i\parallel \gamma \left(t_{i,j}\right)-\left({\mathbf{x}}_{i,j}-\left({\mathbf{x}}_{i,0}-{\mathbf{x}}_q\right)\right){\parallel }^2\\ +\frac{\lambda }{2}\sum _{i=1}^{K_d}{\alpha }_i\parallel v_i{\parallel }^2+\frac{\mu }{2}\sum _{i=1}^{K_d}{\beta }_i\parallel a_i{\parallel }^2.\end{array}$$E22

This minimization process simplifies the problem to a quadratic optimization, solved with LU decomposition. The piecewise nature is represented by the term $K_d\in \mathcal{N}\left({\mathbf{x}}_q\right)$, defined as samples along $\gamma$. The first component of Eq.(22) is a penalty term to minimize the geodesic distance between samples ${\mathbf{x}}_{i,j}$ and the regressed curve, where $w_i$ are weight variables based on sample distances. This helps regress a curve that will lie close to ${\mathbf{x}}_{i,j}$, shifted by ${\mathbf{x}}_q$ in order to have the initial reconstructions co-registered. The second term represents the velocity of the curve (defined by $v_i$, approximating $\dot{\gamma }\left(t_i\right)$), minimizing the $L_2$ distance of the 1${}^{\mathit{st}}$ derivative of $\gamma$. By minimizing the value of the curve’s first derivatives, this prohibits any discontinuities or rapid transitions of the curve’s direction, and is modulated by ${\alpha }_i$. Finally, an acceleration penalty term (defined by $a_i$) focuses on the 2${}^{\mathit{nd}}$ derivative of $\gamma$ with respect to $t_i$ by minimizing the $L_2$ norm. The acceleration is modulated by ${\beta }_i$. Estimates for $v_i$ and $a_i$ (weighted by $\left\{\lambda ,\mu \right\}$, respectively), are generated using geometric finite differences. These estimates dictates the forward and backward step-size on the regressed curve, leading to directional vectors in $\mathcal{M}$ as shown in [22]. In order to minimize $E\left(\gamma \right)$, a nonlinear conjugate gradient technique defined in the low-dimensional space ${\mathbb{R}}^d$ is used, thus avoiding convergence and speed issues. The regressed curve $\gamma$ is therefore defined for all time points, originating at $t_0$. The curve creates a group average of spatiotemporal transformations based on individual correction trajectories.

4.3. Prediction of spine correction

Finally, to predict the evolution of spine correction from an unseen preoperative spine model, we use the geodesic curve $\gamma :{\mathbb{R}}^D\to \mathcal{M}$ modeling the spatiotemporal changes of the spine, where each point $\mathbf{x}\in \mathcal{M}$ is associated to a speed vector $\mathbf{v}$ defined with a tangent plane on the manifold such that $\mathbf{v}\in {\mathrm{T}}_{\mathbf{x}}\mathcal{M}$.

Based on Riemannian theory, an exponential mapping function at $\mathbf{x}$ with velocity $\mathbf{v}$ can be defined from the geodesics such that ${\mathrm{e}}_{\mathbf{x}}^{\mathcal{M}}\left(\mathbf{v}\right)$. Using this concept, parallel transport curves defined in ${\mathrm{T}}_{\mathbf{x}}$ can help define a series of time-index vectors along $\gamma$ as proposed by [10]. The collection of parallel transport curves allows to generate an average trajectory in ambient space ${\mathbb{R}}^D$, describing the spine changes due to the corrective forces of tethering. The general goal is to begin the process at the preoperative sample, and navigate the piecewise-geodesic curve describing correction evolution in time, where one can extract the appearance at any point (time) in ${\mathbb{R}}^D$ using the exponential mapping. For implementation purposes, the parallel transport curve are constrained within a smooth tubular boundary perpendicular to the curve (from an ICA) to generate the spatiotemporal evolution in the coordinate system of the preoperative model.

Hence, given the manifold at time $t_0$ with $\mathbf{v}$ defined in the tangent plane and the regressed piecewise-geodesic curve $\gamma$, the parallel curve is obtained as:

$${\eta }^{\mathbf{v}}\left(\gamma ,s\right)={\mathrm{e}}_{\gamma \left(s\right)}^{\mathcal{M}}\left({\mathbf{x}}_{\gamma ,t_0,s}\left(\mathbf{v}\right)\right),s\in {\mathbb{R}}^d.$$E23

Therefore by repeating this mapping for manifold points seen as samples of individual progression trajectories along $\gamma \left(s\right)$, an evolution model can be generated. Whenever a new sample is embedded, new samples points along $\gamma \left(s\right)$, denoted as ${\eta }^{\mathbf{v}}\left(\gamma ,\cdot \right)$ can be generated parallel to the regressed piecewise curve in $\mathcal{M}$, capturing the spatiotemporal changes in correction.

A time warp function allowing $s$ to vary along the geodesic curve is described as ${\varphi }_i\left(t\right)={\theta }_i\left(t-t_0-{\tau }_i\right)+t_0$. Here, we propose to incorporate a personalized acceleration factor based on the spine maturity and flexibility derived from the spine bending radiographs and Risser grade. A coefficient ${\theta }_i=C_i\times R_i$ describing the change in Cobb angle $C_i$ between poses, and modulated by the Risser grade $R_i$. This coefficient regulates the rate of correction based on the $K$ neighboring samples. Finally, to take under account the relative differences between the group-wise samples and the query model once mapped onto the regressed curve, a time-shift parameter ${\tau }_i$ is incorporated in the warp function.

For spine correction evolution, displacement vectors ${\mathbf{v}}_i$ are obtained by a PCA of the hyperplane crossing ${\mathrm{T}}_{{\mathbf{x}}_i}\mathcal{M}$ in manifold $\mathcal{M}$ [10]. Hence, for any query sample ${\mathbf{x}}_q$ which represents the mapped preoperative 3D reconstruction (prior to surgery), the predicted model at time $t_k$ can be regressed from the piecewise-geodesic curve generated from embedded samples $\mathbf{x}$ in $\mathcal{N}\left({\mathbf{x}}_q\right)$ such that:

$${\mathbf{y}}_{q,t_k}={\eta }^{{\mathbf{v}}_q}\left(\gamma ,{\varphi }_i\left(t_k\right)\right)+{\epsilon }_{q,t_k}$$E24

which yields a predicted postoperative model ${\mathbf{y}}_{q,t_k}$ in high-dimensional space ${\mathbb{R}}^D$, and ${\epsilon }_{q,t_k}$ a zero-mean Gaussian distribution. The generated model offers a complete constellation of interconnected vertebral models composing the spine shape $\mathbf{S}$, at first-erect (FE), 1 or 2-year visits, including landmarks on vertebral endplates and pedicle extremities, which can be used to capture the local shape morphology with the correction process.

4.4. Experiments

The discriminant manifold was trained from a database of $438$ 3D spine reconstructions generated from biplanar images [23], originating from $131$ patients demonstrating several types of deformities with immediate follow-up (FE), 1 and 2 year visits. Patients were recruited from a single center prospective study. Patients were divided in two groups, with the first group composed of 94 responsive patients showing a reduction in Cobb angle over or equal to 10${}^{\circ }$ between the FE and follow-up visit. The second group was composed of 37 nonresponsive (NP) patients with a reduction of less than 10${}^{\circ }$. We evaluated the geometrical accuracy of the predictive manifold for 56 unseen surgical patients (mean age $12\pm 3$, average main Cobb angle on the frontal plane at the first visit was $47\pm {10}^{\circ }$), with predictions at $t=0$ (FE), $t=12$ and $t=24$ months. For the predicted models, we evaluated the 3D root-mean-square difference of the vertebral landmarks generated, the Dice coefficients of the vertebral shapes and in the main Cobb angle. The results are shown in Table 1. Results were confronted to other techniques such as biomechanical simulations performed on each subject using finite element modeling with ex-vivo parameters [25], a locally linear latent variable model [20] and a deep auto-encoder network [24]. Results from the predicted geometrical models show the regressed spatiotemporal geodesic curve yields anatomically coherent structures, with accurate local vertebral morphology.

	FE visit			1-year visit			2-year visit
	3D RMS	Dice	Cobb	3D RMS	Dice	Cobb	3D RMS	Dice	Cobb
Biomec. sim	3.3 $\pm$ 1.1	85 $\pm$ 3.4	2.8 $\pm$ 0.8	3.6 $\pm$ 1.2	84 $\pm$ 3.6	3.2 $\pm$ 0.9	4.1 $\pm$ 2.3	82 $\pm$ 3.9	3.6 $\pm$ 1.0
LL-LVM [20]	3.6 $\pm$ 1.4	83 $\pm$ 4.0	3.8 $\pm$ 1.5	4.7 $\pm$ 3.3	79 $\pm$ 4.4	5.5 $\pm$ 2.6	6.6 $\pm$ 4.4	71 $\pm$ 5.9	7.0 $\pm$ 3.9
Deep AE [24]	4.1 $\pm$ 1.5	80 $\pm$ 4.4	5.1 $\pm$ 2.7	5.0 $\pm$ 1.9	77 $\pm$ 4.9	5.8 $\pm$ 3.0	6.3 $\pm$ 4.6	72 $\pm$ 5.7	6.6 $\pm$ 4.2
Proposed	2.4 $\pm$ 0.8	92 $\pm$ 2.7	1.8 $\pm$ 0.5	2.9 $\pm$ 0.9	90 $\pm$ 2.8	2.0 $\pm$ 0.7	3.2 $\pm$ 1.3	87 $\pm$ 3.1	2.1 $\pm$ 0.6

Table 1.

3D RMS errors (mm), dice (%) and cobb angles (${}^o$) for the proposed method, and compared with biomechanical simulations, locally linear latent variable models (LL-LVM) and deep auto-encoders (AE).

Predictions are evaluated at FE, 1 and 2-years.

5. Discussion

Algorithms capable of extracting clinically relevant and meaningful descriptions from medical imaging datasets have become of widespread interest to theoreticians as well as practitioners in the medical field, accelerating the pace in recent years involving varied fields such as in machine learning, geometry, statistics and genomics to propose new insights for the analysis of imaging and biologic datasets. Towards this end, manifold learning has demonstrated a tremendous potential to learn the underlying representation of high-dimensional, complex imaging datasets.

We presented frameworks describing longitudinal, multimodal image features from neuroimaging data using a Bayesian model for discriminant nonlinear manifolds to predict the conversion of progressive MCI to Alzheimer’s disease. This probabilistic method introduces class-dependent latent variables which is based on the concept that local structure is transformed from manifold to the high-dimensional domain. This variational learning method can ultimately assess uncertainty within the manifold domain, which can lead to a better understanding of relationships between converters and nonconverters for patients with MCI.

Finally, a prediction method for the outcomes of spine surgery using geodesic parallel transport curves generated from probabilistic manifold models was presented. The mathematical models allow to describe patterns in a nonlinear and discriminant Riemannian framework by first distinguishing nonprogressive and progressive cases, followed by a prediction of structural evolution. The proposed model provides a way to analyze longitudinal samples from a geodesic curve in manifold space, thus simplifying the mixed effects when studying group-average trajectories.