Discovery and visualization of structural biomarkers from MRI using transport-based morphometry

TL;DR

Transport-based morphometry (TBM) is a novel MRI analysis technique that automatically uncovers subtle brain tissue changes related to aging and fitness without relying on predefined features, enabling detailed visualization and analysis.

Contribution

The paper introduces TBM, a new image transformation method that fully automates the detection and visualization of complex brain tissue changes from MRI data without prior feature specification.

Findings

TBM successfully identifies aging-related brain changes.

It reveals potential links between aerobic fitness and brain health.

The method enables visualization of tissue distribution shifts.

Abstract

Disease in the brain is often associated with subtle, spatially diffuse, or complex tissue changes that may lie beneath the level of gross visual inspection, even on magnetic resonance imaging (MRI). Unfortunately, current computer-assisted approaches that examine pre-specified features, whether anatomically-defined (i.e. thalamic volume, cortical thickness) or based on pixelwise comparison (i.e. deformation-based methods), are prone to missing a vast array of physical changes that are not well-encapsulated by these metrics. In this paper, we have developed a technique for automated pattern analysis that can fully determine the relationship between brain structure and observable phenotype without requiring any a priori features. Our technique, called transport-based morphometry (TBM), is an image transformation that maps brain images losslessly to a domain where they become much more separable. The new approach is validated on structural brain images of healthy older adult subjects where even linear models for discrimination, regression, and blind source separation enable TBM to independently discover the characteristic changes of aging and highlight potential mechanisms by which aerobic fitness may mediate brain health later in life. TBM is a generative approach that can provide visualization of physically meaningful shifts in tissue distribution through inverse transformation. The proposed framework is a powerful technique that can potentially elucidate genotype-structural-behavioral associations in myriad diseases.

Figures

Tables

Table 1. Table 1: Comparing methods of solving OMT in 2D

2D	OT Mapping Statistics
Method	Relative MSE	Mean curl	Mass transported
Our method	$0.23\pm 0.056\%$	$\left(8.7\pm 5.3\right)\times {10}^{-4}$	$1.63\pm 0.57$
Chartrand et al.	$1.8\pm 2.9\%$	0	$1.56\pm 0.46$
Haker et al.	$0.45\pm 0.59\%$	$\left(7.0\pm 0.16\right)\times {10}^{-6}$	$2.37\pm 0.79$

Table 2. Table 2: Comparing methods of solving OMT in 3D

3D	OT Mapping Statistics
Method	Relative MSE	Mean curl	Mass transported
Our method	$0.55\pm 0.0011\%$	$0.37\pm 0.61$	$1.3\pm 0.50$
Chartrand et al.	$3.0\pm 1.8\%$	0	$0.07\pm 0.04$
Haker et al.	$9.9\pm 2.5\%$	$\left(4.5\pm 1.6\right)\times {10}^6$	$12.5\pm 4.8$

Equations

$${\int }_{A}d\sigma (x)={\int }_{f(A)}d\mu (x),\forall A\subset \Omega .$$

$$\underset{f\in MP}{\min }{\int }_{\Omega }c(x,f(x))d\sigma (x)$$

$${\int }_{\Omega }|x-f(x){|}^{2}d\sigma (x)\le {\int }_{\Omega }|x-g(x){|}^{2}d\sigma (x),$$

$$\forall g\in MP$$

$${\int }_{\Omega }|x{|}^{2}d\mu (x)<\infty \text{and}{\int }_{\Omega }|x{|}^{2}d\sigma (x)<\infty ,$$

$${\int }_{\Omega }{I}_m(x)dx=1.$$

$$\begin{array}{r}{f}_m^{\ast }(x)=arg\underset{f_m\in \mathrm{MP}}{\min }{\int }_{\Omega }(f_m(x)-x{)}^2{I}_0(x)dx,\\ \text{s.t. }\det (D{f}_m(x))I_m(f_m(x))=I_0(x)\text{for}\forall x\in \Omega \end{array}$$

$$\begin{array}{r}I(x)=det(D{f}^{-1}(x))I_0(f^{-1}(x))\\ \text{where }{f}^{-1}(x)\text{ is the inverse mapping of }f(x)\end{array}$$

$$det(Df(x))I_1(f(x))=I_0(x),\forall f\in MP$$

$${\mathrm{argmin}}_f$$

$$M(f)$$

$$\frac{dM}{d{f}^i}=\frac{\partial \mathcal{L}}{\partial f^i}-\sum _{k=1}^n\frac{d}{d{x}^k}(\frac{\partial \mathcal{L}}{\partial f_{x^k}^i}),i=1,...,n$$

$$\frac{dM}{df}$$

$$\displaystyle\left\{\begin{array}[]{l}g^{k}=f^{(k-1)}+\frac{k-2}{k+1}(f^{(k-1)}-f^{(k-2)})\\ f^{k}=g^{k}-\alpha_{k}\frac{dM(g^{(k)})}{df}\end{array}\right.$$

$$w_{corr}=arg\underset{w}{\max }\frac{w^{T}Xv}{\sqrt{w^{T}w}}=\frac{Xv}{\sqrt{v^T{X}^{T}Xv}}$$

$$w_{PLDA}=arg\underset{||w||=1}{\max }\frac{w^T{S}_{T}w}{w^T(S_W+\alpha I)w}$$

$$S_T=U\Sigma U^T$$

$$M(f)$$

$$\frac{d{M}_1}{d{f}^i}=\frac{\partial {\mathcal{L}}_1}{\partial f^i}-\sum _{k=1}^n\frac{d}{d{x}^k}(\frac{\partial {\mathcal{L}}_1}{\partial f_{x^k}^i}),i=1,...,n$$

$$\frac{\partial {\mathcal{L}}_1}{\partial f^i}=det(Df)\frac{\partial I_1(f)}{\partial f^i}{I}_{error}$$

$$det(Df)$$

$$\frac{\partial {\mathcal{L}}_1}{\partial f_{x^k}^i}=C_{i,k}{I}_1(f)I_{error}$$

$$\frac{d{M}_1}{d{f}^i}=I_{error}(det(Df)\frac{\partial I_1(f)}{\partial f^i}-\sum _{k=1}^n\frac{d}{d{x}^k}{C}_{i,k}{I}_1(f))$$

$$\frac{d{M}_1}{df}=I_{error}(det(Df)\nabla I_1(f)-\nabla \cdot (adj(Df)I_1(f))).$$

$$\frac{\partial {\mathcal{L}}_1}{\partial f^i}=0$$

$$|\nabla \times f{|}^2=\sum _{p=1}^n(\sum _{l=1}^n\sum _{m=1}^n{ϵ}^{plm}{f}_{x^l}^m{)}^2$$

$$\frac{\partial {\mathcal{L}}_2}{\partial f_{x^k}^i}=\lambda \sum _{p=1}^n{ϵ}^{pki}(\sum _{l=1}^n\sum _{m=1}^n{ϵ}^{plm}{f}_{x^l}^m).$$

$$\frac{d{M}_2}{d{f}^i}$$

$$\frac{dM}{df}=I_{error}(det(Df)\nabla I_1(f)-\nabla \cdot (adj(Df)I_1(f)))$$

Full text

1 Introduction

Recent advances in magnetic resonance imaging (MRI) technology have enabled high-resolution imaging across many new modalities. Tissue properties can be now be measured at an unprecedented level of precision and detail. These developments hold promise to illuminate structural changes underlying diseases commonly considered medical mysteries. Unfortunately, the changes can often be subtle, spatially diffuse, and complex, escaping detection by visual inspection. For example, Figure 1 demonstrates how the common morphologic pattern that differentiates individuals who are most aerobically fit versus those who are least fit defies identification by gross inspection alone. Thus, there is a growing role for computer-aided techniques to aid in vision and detection of morphologic patterns from MRI. Computer-aided techniques are needed to answer the following questions: are there morphologic differences that differentiate these groups? If so, what are they?

Unfortunately, traditional techniques for MRI analysis have difficulty with analysis in the image domain as well, as they require features to be pre-specified and are prone to missing a multitude of physical changes that are not adequately assessed by these finite feature sets. For example, popular biomedical image analysis softwares such as WND-CHRM [36] or FreeSurfer [17] extract a number of pre-specified numerical descriptors, such thickness, volume, texture statistics, etc., from the images and test whether these quantities are statistically different between image sets using a trial and error approach. In fact, WND-CHRM extracts nearly 3000 generic features from the images for testing. However, not only is testing descriptors a tedious process, numerical descriptors such as SIFT, Gabor features, or histogram statistics often do not have direct biological meaning. Another major limitation of these approaches is that the analysis does not incorporate anatomic prior information, missing an opportunity to compare variations in terms of known anatomy. Deformation-based methods, which include deformation-ased morphometry (DBM) [4], tensor-based morphometry [3], and voxel-based morphometry (VBM) [2], also have limitations. Deformation-based methods rely on a nonrigid registration to align images before comparing them pixelwise. In practice, perfect structural and functional alignment cannot be ensured, and subtle changes in pixel alignment can vastly change the results obtained and undermine accuracy [9]. Another limitation is these techniques is that the deformation fields are not unique. Hence, the results of Tensor-Based Morphometry and deformation-based morphometry, which compare pixelwise across the determinant of Jacobian or deformation fields respectively, will vary depending on the particular field generated by the algorithm. While these methods incorporate anatomical information, projecting images onto individual pixels assumes that the changes are localized into clusters and misses spatially diffuse changes. Furthermore, deformation fields model changes in local volume expansion/contraction in order to match gross shapes. However, as Figure 2 illustrates, deformation fields cannot fully match images with zero error because they cannot quantify differences in tissue topology or texture. Figure 3 illustrates several neurologic diseases for which the main variation is in tissue texture rather than brain volume contraction/expansion - multiple sclerosis lesions and brain tumors. The authors of deformation-based techniques state that these registration-based methods are a way to index into pixels based on amount of gray matter per unit volume [34], but cannot offer insight into the physical meaning of these changes [18] as these methods are not generative. If a technique is not generative, then an observable datapoint, or a brain MR image cannot be generated given a feature set like a density map. However, a generative method would afford the ability to visualize the shifts in morphologic profile as dynamic changes across a series of brain MRIs to illuminate structural mechanisms.

In this paper, we describe transport-based morphometry (TBM) [41, 40, 5, 32, 25], which has the potential to enable fully automated MRI analysis without loss of information. Rather than analyzing images in the image domain where they may not be easily separable, we first transform them to a domain that enhances separability. Prior work demonstrates that when we transform 1D and 2D signals in other applications using TBM [41, 40, 5, 32, 25, 24], complex and nonlinear morphology in the image domain can be described by linear classification and regression models in the transform domain. Furthermore, the key advance of TBM is that it is generative and enables direct visualization of the interface between signal classes through inverse TBM transformation. The TBM technique computes the distance needed to morph one image with respect to a common template using the mathematics of optimal mass transport (OMT). As Figure 2 illustrates, unlike deformation-based approaches, OMT can match both shape and texture variations simultaneously; thus, information is not missed. However, TBM has never been applied for MRI-based detection and pattern analysis as current formulations and solutions to TBM are designed for smaller signals [41, 40, 5, 32, 25, 24].

In this paper, we demonstrate a TBM framework that is suitable for analysis of radiology data, the majority of which comprises three-dimensional data. We hypothesize that transforming MRI data using the new TBM approach can facilitate both discovery as well as visualization of discriminating differences in a manner similar to 1D and 2D signal analysis previously reported [41, 40, 5, 32, 25, 24]. Ultimately, if TBM could assess structural changes underlying clinical phenotypes in a fully automated manner without losing information, and in addition, visualize the shifts in tissue distribution as a series of radiology images, it could represent a breakthrough for scientific understanding as well as identification of objective clinical markers.

The specific contributions of this work are:

Novel, robust formulation and solution of transport-based morphometry (TBM) enabling its first application to radiology data and computational validation

• 2.

A new image transform to facilitate pattern analysis on MRI data, with equations for analysis and synthesis as well as description of how the TBM pipeline can be used for discrimination, regression, and unsupervised learning

• 3.

Demonstration of TBM on real world neuroimage analysis problems showing the advantage of a generative technique in identification of morphologic changes as well as visualization compared to current morphometry techniques

The remainder of this paper is organized as follows. In Section 2, we summarize key theoretical results and present equations for forward and inverse TBM transformation. In Section 3, we present our TBM solver suitable for MRI data. Section 4 describes a framework for regression, discrimination, and blind signal separation tasks in the transform domain. In Section 5, experimental methodology is presented. Section 6 presents results showing robustness of the proposed solver for 3D data and the ability of TBM to accurately assess dependent brain morphology changes with age in Section 6, whereas the traditional approaches based on diffeomorphic anatomic registration through exponentiated lie algebra algorithm (DARTEL) fails to detect these changes [1]. Finally, in Section 7, we present discussion of the results and in Section 8, we conclude this paper. The Appendix material details the derivation and experimental validation for our solver in Section 3.

2 Optimal Mass Transport for Signal Transformation

This section summarizes key theorems related to optimal mass transport, equations for signal transformation using TBM, and OMT minimization.

2.1 Overview of optimal transport theory

Let $\Omega$ be a measurable space. Let $\mu$ and $\sigma$ be probability measures defined on $\Omega$, with corresponding positive probability densities $I_{1}$ and $I_{0}$, respectively. A mass preserving transform $f$ that pushes $\sigma$ to $\mu$, or $f_{\#}\sigma=\mu$, satisfies the following,

[TABLE]

Figure 4 illustrates $\mu$ and $\sigma$, as well as the map $f_{\#}\sigma=\mu$. Such a mass preserving (MP) mapping $f$ is in general not unique; in fact, infinitely many MP mappings may exist that satisfy Equation 1. However, we are interested in finding the MP mapping that is optimal in the sense of mass transport, which we will define further in Equation 2. Let $MP$ be the set of all such mass preserving mappings, $MP:=\{f:\Omega\rightarrow\Omega|f_{\#}\sigma=\mu\}$. The optimal MP mapping in the mass transport sense can be written according to Monge’s formulation, which minimizes the following cost function,

[TABLE]

Here, $c:\Omega\times\Omega\rightarrow\mathbb{R}^{+}$ is the cost functional. The functional $c$ measures mass transportation cost and is often chosen to be the $L^{p}$-norm for which Equation (2) becomes the $L^{p}$-Wasserstein distance. The $L^{2}$-Wasserstein distance, $c(x,y)=|x-y|^{2}$, in particular has attracted rich attention in the image analysis, computer vision, and machine learning communities. For $c(x,y)=|x-y|^{2}$, Brénier [10] showed that there exists a unique optimal transportation map $f\in MP$ for which,

[TABLE]

when (i) $\Omega=\mathbb{R}^{n}$ and the probability measures have finite second-order moments (i.e. their densities vanish in the limit),

[TABLE]

and (ii) when $\sigma$ is absolutely continuous with respect to Lebesgue measure. Note that for certain measures (e.g. when $\sigma$ is not absolutely continuous) the Monge formulation of the optimal transport problem is ill-posed in the sense that there is no transport map that rearranges $\sigma$ into $\mu$. In such scenarios Kantorovich formulation of the problem is preferred.

Moreover, Brénier showed through polar factorization [10] that the transport map $f$ must be the gradient of a convex function $\phi:\Omega\rightarrow\mathbb{R}$, $f=\nabla\phi$. The preceding property implies that when $\Omega$ is a convex and connected subset of $\mathbb{R}^{n}$, the optimal transport map is curl free.

2.2 Linear optimal transport analysis framework

By considering magnetic resonance images to be smooth density functions, the similarity in spatial distribution of two tissues can be quantified based on the $L^{2}$-Wasserstein distance. Any MRI modality that generates scalar intensity maps (i.e. T1-weighted, T2-weighted, FLAIR, fractional anisotropy, etc.) is amenable to analysis by the TBM framework. In this work, we analyze T1-weighted images, where treating images as densities enables comparison between images where the absolute intensity value may not be physically meaningful.

Consider a set of magnetic resonance images $I_{1},...,I_{K}:\Omega\rightarrow\mathbb{R}^{+}$, corresponding to experimental subjects $1,...,K$, where $\Omega=[0,1]^{3}$, the images are first intensity normalized to produce densities such that

[TABLE]

where $m\in{1,...,K}$. A common reference image $I_{0}$, is chosen and the optimal transport mappings are calculated from the reference image to each subject’s MRI, $I_{m}$. Let $f_{m}:\Omega\rightarrow\Omega$ be a mass preserving mapping from $I_{0}$ to $I_{m}$. Then, the analysis equation [25] that transforms images to their corresponding representation in transform domain can be written based on

[TABLE]

Here, $Df_{m}$ is the Jacobian of the mapping $f_{m}$ and ${\operatorname{MP}}$ is the family of all mass preserving mappings from $I_{0}$ to $I_{1}$. The existence of a unique solution $f^{*}_{m}$ to above optimization was shown by Brénier [10].

The transport maps $f^{*}_{m}(x)$ are vector fields that define the direction and amount of mass transport needed to morph $I_{m}(x)$ into $I_{0}(x)$. OMT defines a nonlinear distance metric, as Figure 5 shows, where the arcs, or geodesics, on the manifold between two images $I_{0}$ and $I_{m}$ correspond to nonlinear OMT distances and are represented by $f^{*}_{m}(x)$. The metric space defined by the OMT-based distance metric is a Riemannian manifold, which is equipped with an inner product. Thus, projecting the manifold locally at $I_{0}$ to the tangent space maps the geodesics $f^{*}_{m}$ to linearized versions in the tangent space, called the linearized optimal transport (LOT) metric.

Then, it can be shown that $\hat{I}_{m}(x)=(f^{*}_{m}(x)-x)\sqrt{I_{0}(x)}$ provides a natural isometric linear embedding for image $I_{m}$ with respect to the LOT [41, 25]. This linear embedding is generative, thus, any arbitrary point in the LOT space can be directly inverted and visualized in the image domain [41, 25, 24] as a new image according to the synthesis equation [41, 25]

[TABLE]

LOT is a powerful technique for image analysis because simple Euclidean operations on the LOT-transformed embeddings, such as linear classification and linear regression correspond to nonlinear operations on the OMT manifold. Thus complex, spatially diffuse, nonlinear morphologic shifts in the image domain can be captured by simple Euclidean operations in the LOT domain, facilitating discovery of trends as well as visualization of these trends through inverse transformation [25].

3 Proposed approach

We have developed a method for solving optimal transport that enables the transport-based morphometry technique to be extended to large 3D volumetric images such as MRI. The authors offer this OMT approach as a viable option for carrying out TBM transformation for 3D volumetric images. While numerical OMT is a vast field, a detailed review or evaluation of OMT algorithms in general is beyond the scope of this paper.

3.1 Variational formulation of the problem

In order to find the optimal transport map, we reformulate the minimization in (2) by relaxing the MP constraint. Assuming that $\Omega$ is a convex and connected subset of $\mathbb{R}^{n}$, as it is the case for most image analysis problems (i.e. $\Omega=[0,1]^{n}$), and assuming that the probability measures $\mu$ and $\sigma$ are atomless and absolutely continuous, we can write the differential counterpart of Equation (1) as,

[TABLE]

where $D$ is the Jacobian matrix, and $det(.)$ denotes the determinant operator. The minimization in (2) for $c(x,y)=|x-y|^{2}$ can first be relaxed into the following optimization problem,

[TABLE]

for some small $\epsilon>0$. Next, we use the result from Brénier’s theorem which states that the optimal transport map is a curl free mass preserving map. Therefore we propose to modify the optimization problem in Equation (8) by regularizing the objective function with the curl of the mapping, $f$,

[TABLE]

where $\gamma$ is the regularization coefficient and $\nabla\times(.)$ is the curl operator. We note that modifying (8) to penalize the objective with the curl does not change the optimal solution, but solving (9) in practice helps guide the solution toward the curl free map. We can relax the optimization problem above further and write it as a regularized (or penalized) unconstrained optimization problem,

[TABLE]

Hence the formulation above contains terms explicitly signifying properties of MP mapping $\|det(Df)I_{1}(f)-I_{0}\|^{2}$ and a curl-free mapping $\|\nabla\times f\|^{2}$. The last term implicitly penalizes mappings that are not diffeomorphic when $det(Df(x))$ crosses zero.

The optimization problem in Equation (LABEL:eq:proposed) is not a convex problem. We use a multiscale variational optimization technique to help guide the solution toward the global optimum. We will see in the results section that the multiscale scheme is able to achieve solutions comparable to that obtained using convex methods when they apply. Section 3.3 describes the multiscale variational solver we devise for the optimization in (LABEL:eq:proposed).

3.2 Euler-Lagrange equations

The objective function in (LABEL:eq:proposed) can be written as,

[TABLE]

The Euler-Lagrange equations for the transport field $f$ then are of the form,

[TABLE]

where the superscripts denote the coordinate index for the vectors, and the subscripts denote partial derivatives, $f_{x^{k}}^{i}=\frac{\partial f^{i}}{\partial x^{k}}$. Writing the Euler-Lagrange equations for the objective function in (LABEL:eq:proposed) leads to,

[TABLE]

where $id(x)=x$ is the identity function, $adj(.)$ denotes the adjugate operator, $\nabla\cdot(.)$ is the divergence operator, and $I_{error}=det(Df)I_{1}(f)-I_{0}$. The derivation for the equation above is presented in the Appendix. Equation (13) is a key result from our formulation. The complexity for computing each gradient descent update step here is $\mathcal{O}(NlogN)$, where N is the number of pixels or voxels in the image. The computational complexity is determined by the cost of computing gradients $\mathcal{O}(N)$ in (13), but is dominated by the cost of cubic interpolation in computing the $det(Df)I_{1}(f)$ term.

3.3 Multiscale accelerated gradient descent

We can guide the solution toward the globally optimal solution by a multiscale scheme as depicted in Figure 6. Nesterov’s accelerated gradient descent method [29] is used at each scale to find the corresponding optimal transport map from Equation (LABEL:eq:proposed). The optimal transport map is then interpolated and used as the initial point for the accelerated gradient descent method in the next scale (finer scale).

The accelerated gradient descent update for $k$’th iteration ($k>1$) at each scale is as follows,

[TABLE]

where $\alpha_{k}$ is the gradient descent step size, and is automatically chosen at each gradient descent update such that the maximum displacement is fixed. The update at k = 1 is the usual gradient descent update.

Here, we have presented a viable approach for computing optimal transport minimization for MRI datasets, enabling a transport-based morphometry approach with radiology images.

4 Modeling shape and appearance of the brain

As previous work has established that treating signals as densities increases their separability in the transform domain [41, 40, 5, 32, 25, 24], we describe a framework for regression, discrimination, and blind signal separation in the transport space.

The data matrix $X\in\mathbb{R}^{d\times K}$ stores the vectorized transport maps $x_{m}$ corresponding to each subject $m\in 1,...,K$ where $d$ is the number of elements in the vectorized transport map and $K$ is the number of subjects. Figure 7 is the system diagram that illustrates the LOT transformation pipeline. In practice, the analysis is performed on the dimensionality-reduced data matrix for ease of computation.

4.1 Regression and correlation analysis with a clinical variable

The influence of an independent clinical variable $v\in\mathbb{R}^{K\times 1}$ on brain tissue distribution can be investigated by computing the direction in the transport domain $w_{corr}$ such that the linear correlation with age is maximized according to (17) [5]. Here, $X$ represents the reduced-dimension data matrix.

[TABLE]

Here, the direction $w=\bar{x}+\nu w_{corr}$ is a vector field that represents the direction and magnitude by which tissue is re-distributed due to $v$ and $\nu$ represents the increment or decrement to sample along the maximally correlated direction. Pearson’s correlation coefficient is computed on centered $v$ and $X$.

The images corresponding to the computed direction $w$ can be visualized through inverse TBM transformation by Equation (6) and illustrate the morphology that is associated with outcome $v$.

4.2 Discriminant analysis to differentiate groups of subjects

Another class of problems facilitated by the TBM technique is that of discriminating classes based on MRI appearance, such as the one posed in Figure 1. For these problems, penalized linear discriminant analysis (PLDA) [39] performed in the transport domain can find the direction in transport space that maximally separates $C$ classes. The PLDA direction is given by (18)

[TABLE]

where $S_{T}=\frac{1}{M}\sum_{m}(x_{m}-\bar{x})(x_{m}-\bar{x})^{T}$. Here, $\bar{x}=\frac{1}{M}\sum_{m=1}^{M}x_{m}$.

The within-class scatter matrix is $S_{W}=\sum_{C}\sum_{n\in C}(x_{n}-\bar{x_{c}})(x_{n}-\bar{x_{c}})^{T}$. The parameter $\alpha$ controls the tradeoff between the traditional linear discriminant analysis (LDA) direction and one that lies in the principal component analysis (PCA) subspace. The parameter $\alpha$ can be chosen by plotting the stability of the subspace as a function of $\alpha$.

Sampling along and inverting the direction $w_{PLDA}$ yields images showing the typical morphology of a class and how it changes as one progresses from one class to another.

4.3 Visualizing principal phenotypic variations in the brain

Given the covariance matrix $S_{T}$ defined in Section 4.2, the principal components are given by the eigenvectors of $S_{T}$. The eigenvectors represent the directions in the transport space that capture the main modes of variability in the dataset [5].

The factorization in Equation (19) gives both the principal components and eigenvalues, where the diagonal components of $\Sigma$ represent the variance for each principal component.

[TABLE]

For high dimensional data, the covariance matrix can be implicitly represented using the approach in [14]. Each principal component can be inverted and visualized to yield the principal phenotypic variations that comprise the images in the dataset.

5 Computational experiments

Here we describe image acquisition, preprocessing, morphometry analysis, and statistical learning steps. The code was prototyped in MATLAB (MathWorks, Natick, MA) using built-in libraries.

5.1 Datasets

5.1.1 MRI pattern analysis using transport-based morphometry

The ability of transport-based morphometry to aid in regression, discrimination, and signal separation tasks was assessed on images of 135 healthy subjects, ranging in age from 58 to 81 years (mean age 66.6 years, standard deviation 5.9 years). Both male and female subjects are included. T1-weighted brain images were collected using a 3D Magnetization Prepared Rapid Gradient Echo Imaging (MPRAGE) protocol with 144 contiguous slices. Images were acquired on a 3 T Siemens Allegra scanner with repetition time = 1,800 ms, echo time = 3.87 ms, field of view (FOV) 256 mm, and acquisition matrix $192\times 192$ mm , flip angle = 8 [16]. These images provide an expanded dataset of older subjects on which age-related brain morphology can be investigated.

5.2 Multiscale variational optimal transport

5.2.1 Image preprocessing

Images were skull-stripped and affinely registered to the MNI template using Statistical Parametric Mapping (SPM) software version 12 [30]. Images were normalized so that the sum of intensities was equal in both images (equal mass). By normalizing to a large positive number, $10^{6}$, numerical precision errors resulting from computations with small numbers are avoided. We also add a small constant 0.1 to the normalized images and renormalize so that they are strictly positive [13] to ensure that the OMT problem is well-posed.

5.3 Experiment 1: Modeling the effects of aging on brain tissue distribution

Regression analysis was performed to assess the relationship between brain tissue distribution and age using the approach outlined Section 4.1. The common reference image $I_{0}$ was computed by the Euclidean average of all the subjects. Statistical significance of the computed direction is assessed using permutation testing with T = 1000 tests.

The reults of regression analysis in the transport space were compared to those obtained using deformation-based analysis. The DARTEL [1] toolbox in SPM12 [30] was used to compute deformation fields. DARTEL is commonly used to perform standard VBM and DBM analysis. Images were skull-stripped, segmented, and affine registered to the MNI template similar to the OMT procedure before non-rigid registration by performed by DARTEL. The most correlated direction was computed on the tissue density maps of DARTEL-registered images using Equation (17) for VBM and using the deformation fields for DBM analysis. Modulated versions are used to compensate for the effects of spatial normalization [28].

The TBM analysis is also performed on segmented gray matter and white matter tissue maps separately to enable comparison to VBM.

5.4 Experiment 2: Assessing the effects of aerobic fitness on brain health

Discriminant analysis between high aerobic fitness vs. low aerobic fitness groups is performed using the PLDA approach in transport space described in Section 4.2. Aerobic fitness is measured by $vO_{2}$ L/min. The individuals were grouped into low-fit and high-fit groups based on those who had a $vO_{2}$ L/min greater than one standard deviation above the mean (high-fit: n = 22) and lower than one standard deviation below the mean (low-fit: n = 16).

Discriminant analysis using TBM was compared with that performed using deformation fields instead of transport maps generated by DBM. In VBM analysis, voxelwise comparison on modulated tissue maps was performed to seek the voxel clusters that are significantly different between the two classes using two-sample t-test, uncorrected for multiple comparisons or cluster thresholding.

5.5 Experiment 3: Visualizing principal phenotypic variations

Unsupervised learning using PCA was performed to visualize the top three principal phenotypic variations in the transport space using the approach described in Section 4.3.

6 Results

6.1 Modeling normal variability in the brain

Figure 8 shows the fraction of variance captured by principal components of the image domain (raw voxel values after affine registration) compared to modulated DARTEL registration and optimal mass transport solved using the described approach. Fewer components are needed to represent more of the variance in the transport space than either for image domain or DARTEL pixelwise comparison. Therefore, the information about variability in the dataset appears to be better captured by examining tissue distribution using OMT rather than comparing tissue intensities individually, before or after nonrigid registration. The intuition for how OMT better capture variability in the dataset was previously illustrated by Figure 5.

6.2 Modeling the effects of aging on brain tissue distribution through TBM regression

Aging is clinically known to be associated with tissue atrophy and disproportionate loss of tissue from frontotemporal regions [11]. In this section, TBM is compared with DBM and VBM in the ability to independently discover and model these changes.

6.2.1 Assessing global changes

The direction maximally correlated with age computed in the transport space using TBM is found statistically significant (Pearson’s r = 0.4605, p $<$ 0.001). Figure 9a shows the data when it is projected onto the maximally correlated direction, with each datapoint representing a subject’s image.

The most correlated direction shown in Figure 9a can be inverted to visualize the dynamic changes in morphology underlying the aging process. Figure 9b shows images generated by TBM inverse transformation (images are colorized to aid visual interpretation). We see that the changes captured by the TBM regression framework are well-corroborated by known changes in the clinical literature [11]. Specifically, the changes shown here are enlarging ventricles, especially in slices 75 and 66. There is global tissue thinning, and enlargement of the occipital horns of the lateral ventricles in slice 58. Normal anatomic landmarks characteristic of the brain are also clearly visible in Figure 9b, such as the internal capsule in slice 66 and thalamus in slice 58.

We also compare the results of our TBM regression analysis with that using deformation fields generated by DARTEL, commonly used in DBM analysis. The relationship between the deformation fields and age is found to be statistically significant (Pearson’s r = 0.2918, p = 0.0240), as Figure 10a shows, suggesting that there are significant shape changes with age that are captured by a DBM approach. Compared to the visualizations generated by TBM, we see that the those yielded by the deformation-based approach depict global shape changes but not texture changes. In all slices generated by deformation-based analysis, normal tissue landmarks are distorted. Especially, at the gray-white interface, there are is a ring-like texture that does not represent normal brain anatomy. Examining the associated images generated using the deformation-based approach, we see that while volume expansion of the ventricles is correctly identified, expected changes in tissue distribution are not well-captured using deformations alone. For example, in slice 75, the expected tissue thinning is not well-depicted in the frontal areas. In slice 66, there is an area of bilateral focal hyperintensity near the ventricles. This represents the distorted internal capsule that is correctly represented in slice 66 of Figure 9b. Other normal landmarks such as thalamus and putamen are notably absent in slice 58 as well as the occipital horns of the lateral ventricles. Thus, while DBM can capture global shape changes, such as enlargement of ventricles, texture information is not well modeled and many normal landmarks and are distorted.

6.2.2 Assessing gray matter and white matter changes

Transport-based morphometry was also applied to explore the effect of age on gray matter and white matter maps separately. Figure 11b shows the effect of age on gray matter distribution. The relationship is statistically significant with Pearson’s r = 0.4271 and p$<$0.001. There is thinning of the gray matter tissue when progressing from a 53 year old brain to a 79 year old brain, most markedly in the temporal lobe as can be seen in slice 75. Atrophy can be seen in all the slices by enlargement of the spaces.

The relationship between white matter distribution and age similarly shows atrophy and enlargement of the ventricles in Figure 12. The relationship is statistically significant with Pearson’s r = 0.4026 and p = 0.005. In addition, there appears to be disproportionate loss of white matter tissue from the frontal and temporal regions, which is best illustrated in slices 75 and 66 in Figure 12b.

In contrast, regression analysis performed on the modulated density maps registered by DARTEL, used for VBM, was unable to find a significant relationship between age and either gray matter morphology (Pearson’s r = 0.6787, p = 0.1980) or white matter morphology (Pearson’s r = 0.4965, p = 0.3870).

Figures 13 and 14 show the images generated by attempting to fit a regression model on individual pixel values on a fixed grid. We see that in both cases, progressing from age 53 to age 79, the intensity at voxels in the cortical gray matter is shown to decrease. However, no gross differences in shape are depicted. Similarly, examining the white matter images generated by regression on VBM maps (Figure 14), the intensity in the frontal white matter appears to grossly decrease, especially in slices 75 and 58. However, neither of these relationships were statistically significant, nor do they adequately depict atrophy and loss of tissue from frontotemporal regions.

Overall, the known effects of aging on the brain are best assessed and depicted by the transport-based morphometry technique. A DBM approach does not adequately model tissue texture and VBM can identify intensity changes at fixed voxel locations, but these do not appear to be statistically significant when it comes to modeling the effect of aging on the brain tissue.

6.3 Assessing the effects of aerobic fitness on brain health through TBM discrimination

The effects of aerobic fitness on the brain are assessed by separating high-fit individuals from low-fit individuals using the PLDA approach for discriminant analysis, comparing the ability of TBM, DBM, and VBM to discover and visualize the interface between the two groups.

Clear separation in the training subspace is an expected result whether raw pixels, deformation fields or transport maps are used, as Figures 15a and 16a show, but when visualizing the interface between classes, TBM demonstrates clear advantages in physical interpretability. Visualizing the interface between the high-fitness and low-fitness groups using TBM in Figure 15b, we see that brains corresponding to low-fit individuals appear to demonstrate changes in tissue distribution that are similar to those due to advancing age in Figure 9b. Similarly, those individuals belonging to the high-fit group have brain morphology that appears to be resemble those of younger subjects in an older adult population as the ventricles appear smaller and tissue architecture in the frontotemporal regions are better preserved. Thus, it appears that fitness preserves areas of the brain that are affected in normal aging.

Comparing the results to that obtained utilizing the deformation fields generated by DARTEL that are used in DBM analysis, the DBM visualizations show distortion of tissue topology. Figure 16b appears to depict enlargement of ventricles with low fitness, but normal anatomic landmarks are not easily visualized, including the interface between gray and white matter. Additionally, texture variations are not well-assessed.

The analysis is performed on gray matter and white matter maps separately as well in order to compare the performance of TBM with that of VBM. Figure 18 and 17 show the results when transport-based morphometry is performed on white matter maps and gray matter maps individually.

The interface between the groups is visualized using TBM, which shows loss of temporal lobe gray matter with low fitness in Figure 17b. White matter changes visualized by TBM shows loss of frontotemporal white matter with low fitness and enlarging ventricles in Figure 18b. The pattern of changes in brain tissue distribution seen is similar to that seen in aging.

For VBM analysis, the voxels were compared individually to identify those which had significant differences in intensity across the tissue density maps. The heat maps illustrating voxelwise differences are identified in Figure 19. The clusters here are uncorrected for multiple comparisons, but significance was selected at level $p<0.01$. The clusters identified by VBM appear to be spatially distributed across the entire brain. There are changes identified both in occipital and frontal regions of gray matter, as well as in the periventricular white matter affecting frontal regions predominantly. Interestingly, these are some of the same regions identified to be affected by fitness in Figure 15b. However, the global shifts in tissue profile such as atrophy are not captured or well-indexed by a voxelwise analysis, which is better suited for localizing to clusters.

Therefore, while VBM is better suited to localize changes to specific clusters or anatomic regions, and DBM does not adequately assess tissue topology, transport-based morphometry is able to more fully assess the role of fitness on brain health on tissue distribution later in life to generate direct visualizations of the interface between the two classes.

6.4 Visualizing principal phenotypic variations using TBM for unsupervised learning

Finally, TBM can be used to visualize the top three PCA directions generated in transport space using TBM to gain a sense of the principal modes of variation in the dataset, which show variations in brain size, level of tissue atrophy, and prominence of midbrain structures shown by Figure 20.

7 Discussion

We demonstrate a fully automated technique for MRI analysis that facilitates discovery of structural shifts associated with observable phenotypes called transport-based morphometry (TBM). The results confirm our hypothesis that designing a TBM framework suitable for analysis of radiology images can facilitate tasks of regression, discrimination, and signal separation in the transform domain. Our approach is able to assess and visualize aging-related morphologic changes in a fully automated manner. The changes discovered independently by TBM match those that are well-accepted clinically. Additionally, our technique is able to investigate morphological differences between high-fitness and low-fitness groups and yield physically meaningful visualization of the interface between the two classes, suggesting a new mechanism by which fitness may mediate brain health later in life. Finally, transport-based morphometry is shown to enable signal separation by allowing visualization of biologically interpretable principal phenotypic variations.

Traditional methods for assessing structural correlates in neuroimages, such as those that utilize numerical descriptors or pixelwise comparison are able to test only a subset of the information available. Deformation-based morphometry computes local volume expansion/contraction in terms of deformation fields, but cannot quantify differences in texture. We see that DBM is able to identify volume expansions, but the deformation fields lose information about tissue topology, distorting the texture of normal landmarks in the image. Voxel-based morphometry identifies voxels on a fixed grid, but cannot assess nonlinear or spatially diffuse changes such as atrophy or tissue thinning. Thus, results obtained using deformation-based analysis are influenced by limitations of the method, which confound biological insights.

Transport-based morphometry analyzes tissue spatial distribution in the transform domain, where we see that even linear regression and discrimination techniques in the transform domain are sufficient to assess and visualize changes in tissue distribution that are nonlinear, spatially diffuse, and affect various regions of the brain in unequal ways. There are several reasons why transforming images to transport space enhances a range of pattern analysis tasks. First, optimal mass transport provides a metric by which to compare nonlinear signals whereby distances between images in the image domain can be modeled in terms of geodesics on a Riemannian manifold in transport space [25, 41]. By projecting these geodesics locally to the tangent space, linearized versions of these metrics are available and as we see in this paper, Euclidean models in the transport space can capture a range of nonlinear morphologic changes. Second, the optimal mass preserving mapping has a unique minimizer with a bijective relationship with the source image with respect to a template. Therefore, in addition to enabling complex relationships to be more easily modeled, transport-based morphometry is a generative technique. An observable datapoint can be generated from any arbitrary point in the manifold. Because TBM is generative, a linear model can be directly visualized as a series of physically interpretable images through inverse TBM transformation. In contrast, VBM and DBM are not generative methods. We see indeed that compared to DBM and VBM, transport-based morphometry produces more direct visualizations. Using the TBM technique, we are able to adequately assess both global atrophy and local frontotemporal thinning with aging. In contrast, DBM was able to depict only local volume contraction and VBM was able to localize to individual voxels undergoing density changes. Furthermore, compared to VBM and DBM, transport-based morphometry coupled with discriminant analysis revealed a possible mechanism by which aerobic fitness mediates brain health later in life. By analyzing tissue spatial distribution using OMT, TBM can capture important phenomena that is not considered by VBM or DBM.

There are several limitations of this work. First, our approach for optimal transport minimization is non-convex. Although the approach does not guarantee theoretically that global minima will be achieved, the experimental results demonstrate that the multiscale scheme guides the minimization to the global minima and the results are comparable to those using convex formulations in 2D. We pose it as a future problem to couple the TBM framework presented in this paper with solvers that can overcome limitations with large 3D images and at the same time are convex. Another limitation is that analyzing the spatial distribution of voxels requires a normalization of images. Thus, the TBM transform does not directly consider whether there are statistically significant differences in the sum of voxel intensities. However, the latter limitation is easily remedied, as the sum of voxel intensities can be included as a feature when statistical analyses are performed in the feature domain.

Finally, our formulation and TBM solver are fully general to any image modality and encompasses a wide range of problems in regression, discrimination, and unsupervised learning. Thus, our approach opens the door to numerous research and clinical advances.

8 Conclusion

In conclusion, we presented a novel image transformation framework for MRI data to losslessly facilitate discovery of trends as well as yield biologically interpretable visualization of the morphologic changes associated with a variety of clinical outcomes. We demonstrate that our fully automated approach facilitates regression, discrimination, and blind signal separation with significant advancement over currently used techniques. Our approach is able to to independently discover aging-related changes that are well-corroborated clinically and provide new insight into the effects of fitness on the brain, unlike traditional methods. The results validate that our approach can be used as a statistical learning tool in diseases for which gene-structure-behavior relationships are not well-known.

Appendix A: Derivation of Euler-Lagrange Equation

Here we present the derivation of the Euler-Lagrange equation in (13). Starting from the objective function in Eq. (LABEL:eq:proposed) we have

[TABLE]

where the first term, $M_{1}(f)$, enforces $f$ to be mass preserving while the second term, $M_{2}(f)$, enforces $f$ to be curl free. Starting with the first term we can write the Euler-Lagrange equation as,

[TABLE]

where we have,

[TABLE]

Let $C$ be the cofactor matrix of $Df$. Then $det(Df)$ can be written as the sum of the cofactors of any columns or rows of $Df$,

[TABLE]

Using the cofactor matrix, C, we can write,

[TABLE]

And from Equations (22) and (24) we have,

[TABLE]

and writing the vector form of the above equation for all $i$ and using $C^{T}=adj(Df)$ we can write,

[TABLE]

For the second term, we have

[TABLE]

Furthermore, assuming that $n=2,3$ and using the Levi-Civita symbol $\epsilon$ we can write the norm squared of the curl of $f$ as follows,

[TABLE]

which leads to,

[TABLE]

Therefore we have,

[TABLE]

Finally, combining Equations (26) and (30) will lead to,

[TABLE]

Appendix B: Validating OMT registration on MRI Data

8.1 Optimal transport minimization

The equations for analysis and synthesis can be solved in closed form only for 1D signals [38], but for higher-dimensional signals, they must be solved using optimization techniques. Many solvers have been described in the OMT literature, although special challenges including drift, artifact, and computational time/complexity arise in numerical OMT of large image sets that may exceed millions of voxels.

For example, Haker et al.[20] solve for an initial MP map $f_{0}$ (not unique or optimal) through the Knothe-Rosenblatt rearrangement [38, 8] and then progressively update the initial map using composition (with another MP map $s$ that satisfies $s_{\#}\sigma=\sigma$) so that it becomes curl free to signify optimality. As pointed out by Haber et al.[19] and Rehman et al.[33], however, there exist two main shortcomings to the preceding numerical approaches. First, a robust method is needed to obtain an initial MP mapping, and the obtained initial map is often far from the optimal transport map. Second, and much more importantly, such methods update the transformation in a space which is tangential to the linearized MP constraint. Hence, for any finite step update used in the optimization, $f_{0}(s^{k})$, the mapping deviates or drifts from the set of mass preserving mappings $MP$. While the level of drift may or may not be acceptable in practice for a 2D solution, the drifting is amplified for 3D images as demonstrated in Section 6. The drift phenomenon necessitates solution by alternative methods for 3D images.

Convergent methods have also been proposed based on a fluid dynamics formulation of the problem [6, 26] or based on the solution of the Monge-Ampère equation [7], but these formulations come at the cost of an additional virtual time dimension, which is computationally expensive. Computational cost also becomes a challenge in approaches utilizing a system of linear equations that arise from the finite-difference implementation of the linearized Monge-Ampere equation [35, 31].

Another family of solvers [37, 15, 13] are based on Kantorovich’s formulation of the problem. In short, Kantorovich’s formulation searches for the optimal transport plan $\pi$ defined on $\Omega\times\Omega$ with marginals $\mu$ and $\sigma$ that minimizes the following,

[TABLE]

Here $\Pi(\mu,\sigma)$ is the set of all transport plans with marginals $\mu$ and $\sigma$. Chartrand et al [13] solves the dual problem to the Kantorovich formulation. Chartrand et al obtain the optimal transport map through a gradient descent solution. The obtained transport map as pointed out in [13] and shown in this paper see Figure 22 comes with the trade off of undesired artifacts, especially when the images are not smooth. Thus, additional work is needed to overcome challenges related to quality of MP match with the Chartrand et al.[13] approach.

In summary, in order to test the hypothesis that transport-based morphometry can both extract discriminant information and produce visualization of differences on MRI data, an OMT solution is required that can overcome computational challenges for large 3D data.

We validated our OMT approach on healthy, adult brain images obtained from the IXI dataset, Biomedical Image Analysis Group at the Imperial College in London [22].

10 images were selected at random from Guy’s Hospital in UK. Subjects were male and ranged from 41 to 86 years of age at the time of imaging (mean age 57.8 years, standard deviation 15.7 years). The images were T1-weighted images, obtained using a Philips Medical Systems Intera 1.5 T scanner, with the following imaging parameters. Repetition time = 9.813 ms, echo time = 4.603 ms, number of phase encoding steps = 192, echo train length = 0, reconstruction diameter = 240, flip angle = 8∘. The images are 128 x 128 x 128 matrix, with 1 $mm^{3}$ resolution.

8.1.1 Optimization parameters

The step size for accelerated gradient descent is chosen such that the maximum displacement per update is 0.01 of a pixel, the same for our method and comparison methods. At every step of gradient descent, there is a check to maintain that the mapping is diffeomorphic and the step size is reduced as necessary in order to ensure a diffeomorphic mapping. The parameters were obtained experimentally: $\lambda=100$, $\gamma=6.5\times 10^{4}$ when the MSE reaches 25% of the initial MSE to steer the solution towards a curl-free MP mapping, number of scales = 3. The multiscale approach was also implemented for the two other OMT methods in this paper, although the results did not significantly change when the scheme was used.

The termination criteria for all methods implemented is when MSE of the morphed source image relative to the template image reaches 0.55%. When the drifting phenomena in the Haker method [21] produces MSE 100x the initial MSE, we terminate the code. We report the mean $L_{2}$ norm of the curl per voxel and MSE relative to the template.

We use a numerical discretization scheme in which values are placed at pixel or voxel centers. A consistent second-order finite difference approximation was used for all differential operators, utilizing the DGradient toolbox for MATLAB [27].

8.2 Experiment: Validating OMT registration on MRI

We compare the approach of Haker et al [21] and Chartrand et al [13] to the approach described in this paper. The computational complexity of the gradient descent update step for all three methods implemented is $\mathcal{O}(NlogN)$. All methods implemented utilized the same preprocessed images.

All unique pairs of images were registered to each other for the 10 images, resulting in 45 total registration problems. The statistics reported in this paper are based on the registrations performed in turn with our method, that of Chartrand et al, and that of Haker et al.

The following three experiments investigate TBM for MRI pattern analysis enabled by our OMT approach.

8.3 Comparing MP registration methods

We report results for both 2D and 3D MP registration using optimal transport. The 2D image dataset was derived from the 3D dataset by extracting the same axial slice from the middle of every 3D brain image. The solver and accelerated gradient descent update equations are the same whether working with 2D or 3D images. We compare our OMT approach to current methods in the literature based on both the Monge formulation ([21]) and Kantorovich formulation ([13]) to demonstrate that our method is robust to the challenges of other OMT approaches.

8.3.1 2D optimal mass transport

Table 2 displays the mapping statistics for the set of 2D images. The mass transported is lowest for the Chartrand et al.’s method and compares favorably to the mass transport achieved using our method. The MSE for Chartrand et al indicates that this method also produces the poorest MP mapping of the three methods. The MSE achieved is 3-4 times that achieved by the other two methods, demonstrating that this method is prone to artifacts.

Our method achieves the lowest MSE in addition to mass transport distance. As reviewed in Section 2, the optimal MP mapping is the MP mapping that achieves minimum mass transport. The curl, a measure of optimality of the MP mapping, is 0 by design for the method of Chartrand *et al.*as expected (see Section 2). The mean curl is the highest for our method compared to the other two methods, but is still small in an absolute sense ($10^{-4}$).

Our method produces the best results in terms of mass-preservation and mass transport. In Figure 21, we see the optimal transport fields and their corresponding morphings for several 2D brain images. Visually, the transport fields and quality of morphings are similar for all three methods.

The proposed method was prototyped in MATLAB using built-in functions. The average runtime for $256\times 256$ brain images with $3$ scales was $19.36\pm 7.91$ seconds.

8.3.2 3D optimal mass transport

While in 2D all three methods seem to produce visually similar OMT mappings, we see in 3D that several phenomena become evident. Examining Figure 23, we compare the plots of curl and relative MSE over gradient descent iterations for several brain images. We see that the magnitude of the curl (on the order of $10^{6}$) is large for the Haker et al.’s method. The curl for the Chartrand et al.’s method remains at 0 by design. Our method produces curls for all images that tend toward zero with iterations of gradient descent.

We see that relative MSE with the Haker et al.’s approach increases significantly until we terminate the code when the MSE reaches 100x its initial value. Hence, starting after around 100 iterations of gradient descent, the phenomenon of drift with the Haker et al.’s approach becomes evident. The relative MSE of the Chartrand et al.’s approach decreases, but remains large in magnitude (5-10%) at termination. The large MSE results in visual artifacts in the quality of the MP match, which we can see in Figure 22. In contrast, for our method, all images are able to achieve the 0.55% termination criterion.

Table 2 corroborates the plots in Figure 23. Our method produces the lowest relative MSE (best MP mapping), and all brain images are able to achieve the termination criterion of 0.55%. Furthermore, our curl at termination is 8 orders of magnitude lower than that obtained using the Haker method.

In terms of mass transported, the Chartrand et al.’s method produces the lowest transport distance, although the MSE of the MP mapping is about 5-10x higher than that achieved using our method. We can also see artifacts visually in the mappings produced by the Chartrand et al.’s method compared to our method (Figure 22).

In Figure 22, we compare axial, sagittal and coronal slices mapped using our method and that of Chartrand et al.’s (The method of Haker *et al.*failed to produce a viable solution, which is why it is not shown.) We see that mappings produced by our method result in visually similar images to the target image $I_{0}$, whereas those produced by the Chartrand et al.’s method contain several artifacts.

Overall, our method outperforms both comparison methods for 3D images. Our method achieves the lowest MP mapping, while at the same time achieving small curl and mass transported.

The median runtime was under 20 minutes per brain in MATLAB using built-in libraries on a general purpose computer. There is significant opportunity for improvement with an implementation in native C.

Thus, we see that our approach is able to overcome traditional limitations of drift, artifact, and impractical computational complexity. Our approach enables the goal of pattern analysis on MRI using a transport-based morphometry approach.

Acknowledgment

This work was supported in part by NSF award CCF 1421502, and NIH awards R01 GM090033 as well as the National Institute on Aging awards R01 AG25667 and R01 AG25302. This material is also based upon work supported by the Dowd-ICES graduate fellowship. The authors would like to thank Shlomo Ta’asan, Misha Lavrov for stimulating conversations.