Fast Predictive Simple Geodesic Regression
Zhipeng Dinga , Greg Fleishmanc,d , Xiao Yanga , Paul Thompsonc , Roland Kwitte , Marc Niethammera,b , The
Alzheimer’s Disease Neuroimaging Initiative∗
a Department

of Computer Science, University of North Carolina at Chapel Hill, USA
Research Imaging Center, University of North Carolina at Chapel Hill, USA
c Imaging Genetics Center, University of Southern California, USA
d Department of Radiology, University of Pennsylvania, USA
e Department of Computer Science, University of Salzburg, Austria

arXiv:1711.05766v1 [cs.CV] 15 Nov 2017

b Biomedical

Abstract
Deformable image registration and regression are important tasks in medical image analysis. However, they are computationally expensive, especially when analyzing large-scale datasets that contain thousands of images. Hence, cluster
computing is typically used, making the approaches dependent on such computational infrastructure. Even larger computational resources are required as study sizes increase. This limits the use of deformable image registration and regression
for clinical applications and as component algorithms for other image analysis approaches. We therefore propose using
a fast predictive approach to perform image registrations. In particular, we employ these fast registration predictions to
approximate a simplified geodesic regression model to capture longitudinal brain changes. The resulting method is orders
of magnitude faster than the standard optimization-based regression model and hence facilitates large-scale analysis on
a single graphics processing unit (GPU). We evaluate our results on 3D brain magnetic resonance images (MRI) from
the ADNI datasets.
Keywords: Fast prediction, image regression, ADNI dataset, longitudinal data

1. Introduction
Longitudinal image data provides us with a wealth of information to study aging processes, brain development and
disease progression. Such studies, for example ADNI [1]
and the Rotterdam study [2], involve analyzing thousands
of images. In fact, even larger studies will be available in
the near future. For example, the UK Biobank [3] targets on the order of 100,000 images once completed. With
the number of images increasing, large-scale image analysis typically resorts to using compute clusters for parallel
processing. While this is, in principle, a viable solution,
increasingly larger compute clusters will become necessary
for such studies. Alternatively, more efficient algorithms
can reduce computational requirements, which then facilitates computations on individual computers or much
smaller compute clusters, interactive (e.g., clinical) applications, efficient algorithm development, and use of these
efficient algorithms as components in more sophisticated
analysis approaches (which may use them as part of iterative processes).
∗ Data used in preparation of this article were obtained from
the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database
(adni.loni.usc.edu). As such, the investigators within the ADNI
contributed to the design and implementation of ADNI and/or
provided data but did not participate in analysis or writing of
this report. A complete listing of ADNI investigators can be
found at: http://adni.loni.usc.edu/wp-content/uploads/how_
to_apply/ADNI_Acknowledgement_List.pdf

Image registration is a key task in medical image analysis to study deformations between images. Building on
image registration approaches, image regression models [4,
5, 6, 7, 8, 9, 10, 11, 12, 13] have been developed to analyze
deformation trends in longitudinal imaging studies. One
such approach is geodesic regression (GR) [4, 7, 8] which
(for images) build on the large displacement diffeomorphic metric mapping model (LDDMM) [14]. In general,
GR generalizes linear regression to Riemannian manifolds.
When applied to longitudinal image data, it can compactly
express spatial image transformations over time. However, the solution to the underlying optimization problem is computationally expensive. Hence, a simplified, approximate, GR approach has been proposed [15] (SGR) to
decouple the computation of the regression geodesic into
pairwise image registrations. However, even such a simplified GR approach would require months of computation
time on a single graphics processing unit (GPU) to process
thousands of 3D image registrations for large-scale imaging
studies such as ADNI [1]. The primary reason computational bottleneck for SGR are the optimization required to
compute pair-wise registrations.
Recently, efficient approaches have been proposed for
deformable image registration [16, 17, 18, 19, 20, 21]. In
particular, for LDDMM, which is the basis of GR approaches for images, registrations can be dramatically sped
up, by either working with finite-dimensional Lie algebras [22] and frequency diffeomorphisms [21], or by fast

Input 3D MRI images

Encoder-Decoder Network

SGR

m1

Decoder
Conv> w/ str. 2
+ PReLU

Conv + PReLU

Conv + PReLU

Conv + PReLU

Conv> w/ str. 2
+ PReLU

Conv + PReLU

Conv + PReLU

Conv

Conv + PReLU

Conv + PReLU

Conv + PReLU

Conv> w/ str. 2
+ PReLU

Conv + PReLU

Conv + PReLU

Conv

Conv + PReLU

Conv + PReLU

Conv + PReLU

Conv> w/ str. 2
+ PReLU

Conv + PReLU

Conv + PReLU

Conv

Conv + PReLU
Conv + PReLU

Output

Conv> w/ str. 2
+ PReLU

Conv + PReLU
Conv + PReLU

Conv w/ str. 2
+ PReLU

Conv + PReLU
Conv + PReLU

Conv w/ str. 2
+ PReLU

Conv w/ str. 2
+ PReLU

Conv + PReLU

Conv w/ str. 2
+ PReLU

Conv + PReLU

I (t3 )

Conv + PReLU

I (t0 )

Conv + PReLU

I (t2 )

Conv + PReLU

Patches

Conv + PReLU

Encoder

Conv> w/ str. 2
+ PReLU

I (t1 )

I (t0 ) + m0
x-dim

y -dim

z-dim

mn

I (tn )
Figure 1: Principle of fast predictive simple geodesic regression (FPSGR). In the encoder-decoder network (middle), the inputs are patches
from the moving image and the target image at the same spatial location; the outputs are the predicted initial momenta (i.e., m1 , . . . , mn )
of the corresponding patches. Conv: Convolutional layer; ConvT : transpose of convolutional layer. In the simple geodesic regression (SGR)
part, all the pairwise initial momenta are averaged according to Eq. (9) to produce the initial momentum of the regression geodesic (marked
red).

able approach for brain image registration for the following reasons: First, FPIR predicts the initial momentum of
LDDMM and therefore inherits the theoretical properties
of LDDMM. Consequently, FPIR results in diffeomorphic
transformations, even though predictions are computed in
a patch-by-patch manner; this can not be guaranteed by
most other prediction methods. Second, patch-wise prediction allows for training of the prediction models based on a
very small number of images, containing a large number of
patches. Third, by using a patch-wise approach, even highresolution image volumes can be processed without running into memory issues on a GPU. Fourth, none of the existing predictive methods address longitudinal data. However, as both FPIR and SGR are based on LDDMM, they
naturally integrate and hence result in our proposed fast
predictive simple geodesic regression (FPSGR) approach.

predictive image registration (FPIR) [19, 20]. FPIR predicts the initial conditions (specifically, the initial momentum) of LDDMM, which fully characterize the geodesic
and the spatial transformation using a learned a patchbased deep regression model. Because numerical optimization of standard LDDMM registration is replaced by
a single prediction step, followed by optional correction
steps [20], FPIR is dramatically faster than optimizationbased LDDMM without compromising registration accuracy, as measured on several registration benchmarks [23].
Besides FPIR, other predictive image registration approaches have been proposed. Dosovitskiy et al. [24] use
a convolutional neural network (CNN) to directly predict
optical flow. Liu et al. [25] use an encoder-decoder network to synthesize video frames. Schuster et al. [26] investigate strategies to improve optical flow prediction via
a CNN. Cao et al. [16] use a sampling strategy and CNN
regression to directly learn the mapping from moving and
target image pairs to the final deformation field. Miao
et al. [17] use CNN regression for 2D/3D rigid registration. Sokooti et al.[18] use CNNs to directly predict a 3D
displacement vector field from input image pairs. An unsupervised approach for image registration was proposed
by de Vos et al. [27]; here, the loss function is the image
similarity measure between images themselves and a deformation is parameterized via a spatial transformer (which
essentially amounts to a parameterized model of deformation in image registration) which generates the sought-for
displacement vector field. In [28], Hong et al. employ
a low-dimensional band-limited representation of velocity fields in Fourier space [22] to speed up SGR [15] for
population-based image analysis.

Our contributions can be summarized as follows:

Predictive geodesic regression We use a fast predictive registration approach for image geodesic regression. Different to [20], we specifically validate that
our approach can indeed capture the frequently subtle deformation trends of longitudinal image data.
Large-scale dataset capability Our predictive regression approach facilitates large-scale image regression
within a short amount of time on a single GPU, instead of requiring months of computation time for
standard optimization-based methods on a single computer, or on a compute cluster.
Accuracy We assess the accuracy of FPSGR by (1) studying linear models of atrophy scores (which are derived from the nonlinear SGR model) over time, as

In this work, we will build on FPIR, as it is a desir2

−1
be computed by solving Φ−1
v = 0, Φ−1 (0) = id,
t + DΦ
where id denotes the identity map.

well as (2) correlations between atrophy scores and
various diagnostic groups.
Validation We demonstrate the performance of FPSGR
by analyzing > 6000 images of the ADNI-1 / ADNI-2
datasets. For comparison, we also perform SGR
using numerical optimization for the registrations,
again on the complete ADNI-1 / ADNI-2 datasets.

2.2. FPIR
Fast predictive image registration [19, 20] aims at predicting the initial momentum, m0 , between a source and a
target image patch-by-patch. Specifically, we use a deep
encoder-decoder network to predict the patch-wise momentum. As shown in Fig. 1, in 3D the inputs are two
layers of 15 × 15 × 15 image patches (15 × 15 in 2D), where
the two layers are from the source and target images respectively. Two patches are taken at the same position by
two parallel encoders, which learn features independently.
The output is the predicted initial momentum in the x,
y and z directions (obtained by numerical optimization
on the training samples). Basically, the network is split
into an encoder and a decoder part. An encoder consists
of 2 blocks of three 3 × 3 × 3 convolutional layers with
PReLU activations, followed by another 2 × 2 × 2 convolution+PReLU with a stride of two, serving as a “pooling”
operation. The number of features in the first convolutional layer is 64 and increases to 128 in the second. In
the decoder, three parallel decoders share the same input
generated from the encoder. Each decoder is the inverse
of the encoder except for using 3D transposed convolution
layers with a stride of two to perform “unpooling”, and
no non-linearity at the end. To speed up computations,
we use patch pruning (i.e., for brain imaging, e.g., patches
outside the brain are not predicted as the momentum is
expected to be zero there) and a large pixel stride (e.g.,
14 for 15 × 15 × 15 patches) for the sliding window of the
predicted patches.

This work is an extension of a recent conference paper [29].
In particular, all our experiments are now in 3D. We also
added significantly more results to further explore the behavior of FPSGR in comparison to optimization-based SGR.
In particular, we added (a) an example to visualize the
performance of regression models and associated quantitative comparisons (Sec. 4.1); (b) an analysis of local atrophy score correlated with clinical variables (Sec. 4.3);
(c) correlations within diagnostic groups (Sec. 4.3); (d) a
comparison with pairwise registration (Sec. 4.4); (e) and
experiments on extrapolation on unseen data (Sec. 4.5,
Sec. 4.6).
Organization. The remainder of this article is organized
as follows: Sec. 2 describes FPSGR, Sec. 3 discusses the
experimental setup and the training of the prediction models. In Sec. 4, we present experimental results for 3D MR
brain images. The paper concludes with a summary and
an outlook on future work.
2. Fast predictive simple geodesic regression
Our fast predictive simple geodesic regression approach is
a combination of two methods: First, fast predictive image
registration (FPIR) and, second, integration of FPIR with
simple geodesic regression (SGR). Both FPIR and SGR are
based on the shooting formulation of LDDMM [7]; Fig. 1
illustrates our overall approach. The individual components are described in the following.

2.3. Correction network

We follow [20] and use a two-step approach to improve
overall prediction accuracy. An additional correction step,
i.e., a correction network, corrects the prediction of the initial prediction network. Fig. 2 illustrates this two-step ap2.1. LDDMM
proach graphically. The correction network has the same
Shooting-based LDDMM and geodesic regression minimize
structure as the prediction network. Only the inputs and
X
outputs differ. For the prediction network, the inputs are
1
1
d2 (I(ti ), Y i ), (1) the original moving image and the original target image;
E(I0 , m0 ) = hm0 , Km0 i + 2
2
σ i
output is the predicted initial momentum. For the cors.t. mt + ad∗v m = 0, It + ∇I T v = 0, m − Lv = 0,
rection network, the inputs are the original moving image
and the warped target image; the output is the momentum
where I0 is the initial image (known for image-to-image
difference.
registration and to be determined for geodesic regression),
m0 is the initial momentum, K is a smoothing operator
2.4. SGR
that connects velocity v and momentum m as v = Km
Determining the initial image, I0 , and the initial momenand m = Lv with K = L−1 , σ > 0 is a weight, Y i is
tum, m0 , of Eq. (1) is computationally costly. However,
the measured image at time ti (there will be only one
in simple geodesic regression, the initial image is fixed to
such image for image-to-image registration at t = 1), and
2
the first image of a subject’s longitudinal image set (leftd (I1 , I2 ) denotes the image similarity measure between I1
∗
most part of Fig. 1). Furthermore, the similarity measure
and I2 (for example L2 or geodesic distance); ad is the
d(·, ·) is chosen as the geodesic distance between images
dual of the negative Jacobi-Lie bracket of vector fields:
and approximated so that the geodesic regression problem
adv w = −[v, w] = Dvw − Dwv and D denotes the Jaco−1
can be solved by computing pair-wise image registrations
bian. The deformation of the source image I0 ◦ Φ can
3

Final
momentum

Moving image

+
Initial
momentum mp

Prediction network

Correction network

Backward
warp Φ

Correction
momentum mc

Target image

Figure 2: Architecture of the prediction + correction network. Here, we use 2D images and the momentum in the x-direction for illustration.
All images are 3D in our experiments. (1) Predict the initial momentum mp and the corresponding backward deformation, Φ; (2) Predict
a correction of the initial momentum, mc , based on the difference between the moving image and the warped-back target image. The final
momentum is m = mp + mc . The correction network is trained based on the moving images and the warped-back target images of the
training dataset.

3D Longitudinal Test Case Deformation Error [pixel]
Data Percentile
0.3%
5%
25%
50%
75%
95%
Longitudinal Training
0.0156 0.0407 0.0761 0.1098 0.1559 0.2681
Cross-sectional Training 0.0544
0.1424
0.2641
0.3723
0.5067
0.7502
3D Cross-sectional Test Case Deformation Error [pixel]
Data Percentile
0.3%
5%
25%
50%
75%
95%
Longitudinal Training
0.1694
0.4802
1.0765
1.7649
2.7630
4.8060
Cross-sectional Training 0.1123 0.3024 0.5863 0.8737 1.2743 2.2659

99.7%
0.3238
0.8425
99.7%
5.6826
2.7836

Table 1: Deformation error of longitudinal and cross-sectional models tested on longitudinal and cross-sectional data. 2-norm deformation
errors in pixels w.r.t. the ground truth deformation obtained by numerical optimization for LDDMM. A prediction model trained with
longitudinal registration performs better for longitudinal registrations. Conversely, a model trained based on cross-sectional registration is
preferred for cross-sectional registrations.

with respect to the first image. Specifically, we define the
quadratic distance d2 between two images A and B as
Z
1 1 ∗ 2
2
kv kL dt,
(2)
d (A, B) =
2 0
Z 1
1
1
where v ∗ = arg min
kvk2L dt + 2 kQ(1) − Bk22 ,
2
σ
v
0

Hence, the squared geodesic distance between the two images can be approximated as

where v A
comes

s.t. Qt + ∇QT v = 0, and Q(0) = A .

1
E(I, m) = hm, Kmi
2
1 X
(ti − t0 )2 hK(m − mi ), m − mi i,
+ 2
2σ i

Assume we have an image I(t0 ) at time t0 as well as two
images A(ti ) and B(ti ). Further, assume that the spatial
transformation ΦA maps A(ti ) to I(t0 ) and ΦB maps B(ti )
−1
to I0 . Then A(ti ) = I(t0 ) ◦ Φ−1
A and B(ti ) = I(t0 ) ◦ ΦB .
Furthermore, assume that Φ maps A(ti ) to B(ti ), i.e.,
B(ti ) = A(ti ) ◦ Φ−1 . Then Φ = ΦB ◦ Φ−1
A . Assuming
that the geodesic between I(t0 ) and A(ti ) is parameterized by the initial velocity v A and between I(t0 ) and B(ti )
by the initial velocity v B and that we travel between I(t0 )
and A(ti ) in time ti − t0 (and similarly for B(ti )) we can
rewrite the map between A(ti ) and B(ti ) based on the
exponential map as
Φ = ExpId ((ti − t0 )v B ) ◦ ExpId (−(ti − t0 )v A ),

A

Φ ≈ ExpId ((ti − t0 )(v − v )).

(6)

where m is the sought-for initial momentum of the regression geodesic and mi are the initial momenta corresponding to the geodesic connecting I (the starting image of the
geodesic) and the measurements Yi in time ti − t0 . Differentiating Eq. (6) w.r.t. m results in
∇m E = K[m +

(3)

1 X
!
(ti − t0 )2 (m − mi )] = 0.
2
σ i

(7)

Thus,
2
i − t0 ) mi
i (t
P

P

which can be approximated to first order as
B

1
(ti − t0 )2 hK(mB − mA ), mB − mA i,
2
(5)
= KmA and v B = KmB . Hence, Eq. (1) be-

d2 (A(ti ), B(ti )) ≈

m=
(4)
4

σ2 +

i (ti

− t0 )2

.

(8)

In practice, σ 2 is very small and can thus be omitted. Furthermore, mi is obtained by either registering I to Y i in
unit time or, as in our FPSGR approach, by predicting the
momenta mi via FPIR, denoted as m
e i . As Equation 8 was
derived assuming that images are transformed into each
other in time ti − t0 instead of unit time, the obtained
unit-time predicted momenta m
e i correspond in fact to the
approximation m
e i ≈ (ti −t0 )mi . Finally, we obtain the approximated optimal m of the energy functional in Eq. (1),
for a fixed I = I0 as
P
(ti − t0 )m
ei
.
(9)
m ≈ Pi
2
(t
−
t
)
i
0
i

Training for registration. We randomly selected 120
subjects from ADNI-1 and registered their baseline images
to their 24 month follow-up images. We used the first 100
subjects for training and the remaining 20 subjects for
testing. For longitudinal training, we registered the baseline image of a subject to the subject’s 24-month image.
For cross-sectional training, we registered a subject’s baseline image to another subject’s 24-month image. To assess
the performance of prediction models trained on these two
types of paired data, we (1) perform the same type of
registrations on the held-out 20 subjects and (2) compare
the 2-norm of the deformation error computed from the
output of the prediction models with respect to the result
obtained by numerical optimization of LDDMM2 (which
serves as the “ground-truth”). Table 1 shows the results
which confirm our hypothesis that training the prediction
model with longitudinal registration cases is preferred for
longitudinal registration over training with cross-sectional
data. The deformation error is very small for longitudinal
training / testing which provides strong evidence that the
predictive method exhibits performance comparable to the
(costly) optimization-based LDDMM. Another interpretation of these results is, that it is beneficial to train a prediction model with deformations that are to be expected, i.e.,
relatively small deformations for longitudinal registrations
and larger deformations for cross-sectional registrations.
As we are interested in longitudinal registrations for the
ADNI data, we only train our 3D models using longitudinal
registrations in the following.

3. Setup / Training
All our experiments use 3D images from the ADNI dataset1
which consists of 6471 3D MR brain images of size 220 ×
220 × 220 voxels. In particular, ADNI-1 contains 3479 images from 833 subjects and ADNI-2 contains 2992 images
from 823 subjects. Images belong to various types of diagnostic categories which we will discuss later.
We perform the following two types of studies:
Registration We assess our hypothesis that training FPIR
on longitudinal data for longitudinal registrations
is preferred over training using cross-sectional data.
Vice versa, training FPIR on cross-sectional data for
cross-sectional registrations is preferred over training using longitudinal data. Comparisons are with
respect to registration results obtained by numerical
optimization (i.e., LDDMM).

Training for regression. The ADNI-1 dataset contains
228 normal controls, 257 subjects with mild cognitive impairment (MCI), 149 with late mild cognitive impairment
(LMCI), as well as 199 subjects suffering from Alzheimer’s
disease (AD). We randomly picked roughly 1/6 of patients
from each diagnostic category to form a set of 139 subjects
for training in ADNI-1, i.e., 38 normal controls, 43 MCI,
25 LMCI, as well as 33 AD subjects; this results in 139
subjects overall. The baseline images of each subject were
registered to all the later time-points within the same subject. To maintain the diagnostic ratio, we picked (out of
all registrations) 45 registrations from the normal group,
50 registrations from the MCI group, 30 registrations from
the LMCI group, and 40 registrations from the AD group,
resulting in 165 longitudinal registration cases for training.
The same strategy was applied to ADNI-2. In detail,
ADNI-2 contains 200 normal controls, 111 subjects with
significant memory complaint (SMC), 182 subjects with
early mild cognitive impairment (EMCI), 175 with late
mild cognitive impairment (LMCI), and 155 subjects with
Alzheimer’s disease (AD). We picked 150 subjects and
140 longitudinal registrations, consisting of 35 registrations from the control group, 20 registrations from the
SMC group, 30 registrations from the EMCI group, 30

Regression As for regression, we compare linear models
fitted to atrophy scores over time, where scores are
either obtained from FPSGR or optimization-based
SGR. Additionally, we study correlations between
atrophy scores and diagnostic groups. Our hypothesis is that FPSGR is accurate enough to achieve comparable performance to optimization-based SGR, at
much lower computational cost, in both situations.
3.1. Training of the prediction models
We use a randomly selected set of 120 patients’ MRI images from ADNI for training the prediction models and to
test the performance of FPIR. We use all of the ADNI data
for our regression experiments.

1 Data used in the preparation of this article were obtained from
the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database
(adni.loni.usc.edu). The ADNI was launched in 2003 as a publicprivate partnership, led by Principal Investigator Michael W. Weiner,
MD. The primary goal of ADNI has been to test whether serial
magnetic resonance imaging (MRI), positron emission tomography
(PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild
cognitive impairment (MCI) and early Alzheimer’s disease (AD).

2 LDDMM results are generated using a vector momentum formulation: https://bitbucket.org/scicompanat/vectormomentum

5

Distribution
Pred-1 6mo
NC
182
MCI∗
274
AD
153
Total
609
Pred-2 6mo
NC
182
MCI∗
272
AD
152
Total
606

registrations from the LMCI group, and 25 registrations
from the AD group. Note that there are fewer registrations than subjects (140 vs. 150) in this setup, as our
priority is to maintain the overall diagnostic ratio.
For both, ADNI-1 and ADNI-2, the remaining 5/6 of
the data is used for testing. We trained four prediction
models and their four corresponding correction models,
leading to eight prediction models in total, listed in Table
2. We also note that the training sets within ADNI-1 and
ADNI-2, resp., were not overlapping.
ADNI-1 Pred-1
ADNI-1 Pred+Corr-1
ADNI-1 Pred-2
ADNI-1 Pred+Corr-2
ADNI-2 Pred-1
ADNI-2 Pred+Corr-1
ADNI-2 Pred-2
ADNI-2 Pred+Corr-2

Model
Model
Model
Model
Model
Model
Model
Model

v1
v1
v2
v2
v1
v1
v2
v2

(no corr.)
+1x corr.
(no corr.)
+1x corr.
(no corr.)
+1x corr.
(no corr.)
+1x corr.

step

of prediction cases
12mo 18mo 24mo
172
8
151
221
165
122
173
66
163
566
239
436
12mo 18mo 24mo
168
9
144
224
169
124
168
64
160
560
242
428

in ADNI-1
36mo 48mo
128
38
80
11
69
20
277
69
36mo 48mo
119
33
70
10
67
22
256
65

Table 3: Distribution of Pred/Corr-1 and Pred/Corr-2 cases in
ADNI-1. MCI∗ is the combination of the MCI and LMCI diagnostic
groups. 18 month only has one diagnostic group.

step

Distribution of prediction cases in ADNI-2
Pred-1 3mo
6mo
12mo 24mo 36mo
NC∗
173
141
153
119
3
MCI∗
256
232
207
142
4
AD
93
95
105
66
1
Total
522
468
465
327
8
Pred-2 3mo
6mo
12mo 24mo 36mo
NC∗
172
142
159
122
3
MCI∗
257
230
202
149
4
AD
94
98
101
52
1
Total
523
470
462
323
8

step
step

Table 2: Overview of the trained prediction models.

3.2. Parameter selection
We use the regularization kernel
K = L−1 = (−a∇2 − b∇(∇·) + c)−2
with [a, b, c] set to [1, 0, 0.1]. The parameter σ, from equation (1), is set to 0.1. We train our network (using ADAM [30])
over 10 epochs with a learning rate of 0.0001.

Table 4: Distribution of Pred/Corr-1 and Pred/Corr-2 cases in
ADNI-2. Normal∗ denotes the combination of the Normal and SMC
diagnostic groups; MCI∗ denotes the combination of the EMCI and
LMCI diagnostic groups. Only a small number of images is available
for the 36 months time point.

3.3. Efficiency
Once trained, the prediction models allow fast computations of registrations. We use a TITAN X (Pascal) GPU
and PyTorch3 for our implementation of FPIR. For the
3D ADNI-1 dataset (220 × 220 × 220 MR images), FPSGR
took about one day to predict 2646 pairwise registrations
(i.e., 25 [s]/prediction) and to compute the regression result. Optimization-based LDDMM4 would require ≈ 40
days of runtime. Runtime for FPIR on ADNI-2 is identical
to ADNI-1 as the images have the same spatial dimension.
Compared to the-state-of-art fast geodesic regression
model [28], FPSGR is also at least twice as fast. The
model in [28] achieves ≈ 16 times speed-up compared with
SGR [15] for the same setting (parallel computing with the
same number of cores). In our case, we achieve a more
than 40 times speed-up compared with SGR for the same
setting (a single GPU).

Figure 3: Region of Interest (ROI) significantly associated with
atrophy in AD used to compute atrophy scores.

4. Experimental results for 3D ADNI data
For our experiments, we created 10 different (dataset, registration approach) combinations, each combination specifically designed to assess certain properties of our proposed
strategy. These combinations are described next.

3 http://pytorch.org
4 Here, we used 300 fixed iterations for each registration. 300
iterations can guarantee almost all the results converge. Note that
the optimization-based LDDMM also uses a GPU implementation.

1) All subjects from the ADNI-1 dataset in combination
6

with optimization-based LDDMM.

SMC), mild cognitive impairment (including EMCI and
LMCI) (1), and AD (2).

2) Two subgroups of ADNI-1 (i.e., different training data
portions) in combination with FPSGR without a correction network.

Specifically, we investigate the following six questions:
Q1) Can the prediction models for regression qualitatively capture similar trends to the regression model
obtained by numerical optimization? (Sec. 4.1)

3) The same two subgroups as in 2), but in combination
with FPSGR with a correction network.

Q2) Are atrophy measurements derived from FPSGR biased to overestimate or underestimate volume changes?
(Sec. 4.2)

4) The same five groups of 1-3, but for ADNI-2.
Our general hypothesis is that the prediction models
(for ADNI-1/2) show similar performance to optimizationbased LDDMM and that using the correction network for
the predictions improves results. To assess differences, we
compare differences in deformations. Specifically, for every deformation produced by the different approaches, we
compute its Jacobian determinant (JD). The JDs are then
warped to a common coordinate system for the entire ADNI
dataset using existing non-linear deformations from [31,
32]. Each such spatially normalized JD is then averaged
within a region where the rate of atrophy is significantly
associated with Alzheimer’s disease (AD), i.e., within a
statistical region of interest (stat-ROI) (see Fig. 3). Specifically, we quantify atrophy as


Z
1
det(Dφ(x)) dx × 100
(10)
1−
|ω| ω

Q3) Are FPSGR atrophy measurements consistent with
those derived from deformations via numerical optimization (LDDMM) which produced the training
dataset? (Sec. 4.3)
Q4) Are regression results more stable and hence capture
trends better than pairwise registrations? (Sec. 4.4)
Q5) Is the predictive power of the regression models strong
enough to forecast deformations for unseen future
timepoints (Sec. 4.5)
Q6) Do the prediction results capture expected trends in
deformation? (Sec. 4.6)
If these experiments resolve favorably, then the substantially improved computational efficiency of FPSGR
justifies its use for large-scale imaging studies. Tables 3
and 4 show the distributions of the prediction cases per
time-point and the diagnostic groups in ADNI-1 and ADNI-2,
respectively.

where det(·) denotes the determinant and | · | the cardinality/size of a set; ω is an area in the temporal lobes
which was determined in prior studies [31, 32] to be significantly associated with accelerated atrophy in Alzheimer’s
disease. The resulting scalar value is an estimate of the relative volume change experienced by that region between
the baseline and a follow-up image. Hence, its sign is positive when the region has lost volume over time and is
negative if the region has gained volume over time.
We limited our experiments to the applications in [33,
34], wherein nonlinear registration/regression is used to
quantify atrophy within regions known to be associated to
varying degrees with AD (2), mild cognitive impairment
(MCI) (1) (including LMCI5 ), and normal ageing (NC:
normal control) (0) in an elderly population. These are
the diagnostic groups for ADNI-1. For ADNI-2, there are
also 3 diagnostic categories6 : normal ageing (0) (including

4.1. Regression results
Table 1 indicates that FPIR can predict deformation
fields similar to the ones obtained using optimization-based
LDDMM, even for the subtle changes seen in longitudinal imaging data. However, it remains to be seen how
a predictive model performs for image regression. Fig. 4
shows an exemplary regression result. In this specific case,
large changes can be observed around the ventricles. To
illustrate differences between the methods, Fig. 4 shows
regression results based on optimization-based LDDMM,
for FPSGR without a correction network, and for FPSGR
with a correction network. All three methods successfully
capture the expanding ventricles and generally capture the
image changes. Both FPSGR methods show results that
are highly similar to SGR using optimization-based LDDMM. Hence, FPSGR is useful for longitudinal image regression. To further quantify the regression accuracy, we
compute the overlay error between measured images and
the images on the geodesic as

5 We combine MCI and LMCI mainly because (a) the diagnostic
changes available on the IDA website (https://ida.loni.usc.edu/
login.jsp) only provide these three diagnostic groups; (b) to be
consistent with the experiments conducted by Hua et al. [33], where
only Normal, MCI and AD were used as labels to classify ADNI-1.
Hereafter, in all discussions of ADNI-1, MCI is a combination of MCI
and LMCI of ADNI-1
6 Similar to ADNI-1, a detailed diagnosis for ADNI-2 is only available for the baseline images; MR images at later time points are
only labeled as NC, MCI, and AD. Thus, we combine SMC and
NC, as well as EMCI and LMCI to be consistent with the diagnostic changes in the ADNI Diagnosis Summary available on the IDA
website. Hereafter, in all discussions of ADNI-2, NC includes NC and
SMC and MCI includes EMCI and LMCI.

Eoverlay (I0 ◦ Φ−1
ti , Yi ) =

1
kI0 ◦ Φ−1
ti − Yi kL1
|Ω|

(11)

where Ω is the brain area, I0 ◦ Φ−1
ti is the regressed image
at time ti and Yi is the measured image at time ti . Table 5 shows the overlay error for the population of 100 subjects which includes all diagnostic groups in ADNI-1. Both
7

Figure 4: Exemplary regression result: one subject with 6 follow-up images from the ADNI-1 dataset. Image intensity range is [0, 2.49]. Top
row: Axial slices extracted from the 3D MR images at the same axial location for different months. Original: intensity differences between
the baseline image and its 6-month, 12-month, etc. follow-up image. LDDMM: intensity differences between the acquired images in the
top row and optimization-based regression results at each follow-up month(s). Pred-1: intensity differences between the acquired images in
the top row and the Pred-1 regression results at each follow-up month(s). Pred+Corr-1: Same as for Pred-1, but using the Pred+Corr1 regression model. Rectangles mark areas of major structural changes. Intensity differences are dramatically reduced, e.g., around the
ventricles, demonstrating that these structural changes are captured by all three methods. The prediction models (Pred-1, Pred+Corr-1) give
very similar results to the regression results obtained by numerical optimization (LDDMM).

FPSGR methods obtain results comparable with optimization- (5) starting with MCI and changing to AD at later time
based LDDMM. This justifies the use of the proposed
points (MCI-AD), and (6) AD for all the time points (ADmethods. The correction network generally increases the
AD)7 . In particular, we follow [33] and fit a straight line
prediction accuracy over using the prediction network only.
(i.e., linear regression) through all atrophy measurements
over time, conditioned on each diagnostic change category.
4.2. Bias
The intercept term is an estimate of the atrophy one would
measure when registering two scans acquired on the same
Estimates of atrophy are susceptible to bias [35]. To
day; hence it should be near zero and its 95% confidence
quantitatively assess this potential bias, we separately coninterval should contain zero. Quantitatively, Table 6 lists
sidered different diagnostic groups. Specifically, we considthe slopes, intercepts, and 95% confidence intervals for all
ered six diagnostic change groups in our experiments: (1)
NC for all time points (NC-NC), (2) starting with NC and
changing to MCI or AD at a later time point (NC-MCI),
7 In ADNI-1/ADNI-2, there are two patients who show a reversion
(3) MCI for all time points (MCI-MCI), (4) starting with
from AD to MCI. We omitted these cases in our experiment because
MCI and reversing to NC at later time points (MCI-NC),
the number of such cases is too small.
8

Measured Images
Original
LDDMM
Pred-1
Pred+Corr-1

I6mo
0.0770 ± 0.0212
0.0750 ± 0.0194
0.0754 ± 0.0213
0.0754 ± 0.0211

I12mo
0.0764 ± 0.0207
0.0686 ± 0.0176
0.0694 ± 0.0182
0.0691 ± 0.0182

Eoverlay (I0 ◦ Φ−1
ti , Yi )
I18mo
I24mo
0.0890 ± 0.0220 0.0810 ± 0.0223
0.0734 ± 0.0190 0.0609 ± 0.0168
0.0742 ± 0.0195 0.0621 ± 0.0188
0.0734 ± 0.0192 0.0615 ± 0.0166

I36mo
0.0899 ± 0.0341
0.0628 ± 0.0177
0.0654 ± 0.0184
0.0642 ± 0.0188

I48mo
0.0940 ± 0.0415
0.0663 ± 0.0221
0.0698 ± 0.0238
0.0688 ± 0.0235

Table 5: Mean+standard deviation of the overlay errors, see Eq. (11), over 100 patients in ADNI-1 dataset. Both prediction models
exhibit performance comparable to optimization-based regression results (LDDMM). Including a correction network generally improves the
performance of the prediction network.

ten groups of ADNI-1 and ADNI-2, respectively. LDDMMgroups. Note that our LDDMM optimization results and
1 and LDDMM-2 denote the optimization-based results
the prediction results show the same trends. Further, they
split into the same testing groups used for Pred-1 and
are directly comparable as the results are based on the
Pred-2 to allow for a direct comparison. All of the results
same test images (also for the atrophy measurements).
show intercepts that are near zero relative to the range of
changes observed and all prediction intercept confidence
4.3. Atrophy
intervals contain zero. For all diagnostic change groups
Atrophy estimates have also been shown to correlate
the prediction and prediction+correction models exhibit
with clinical variables [31]. To quantify this effect, we commore stable results than the optimization-based LDDMM
puted the Spearman rank-order correlation8 between our
method as indicated by the tighter confidence intervals.
atrophy estimates and the diagnostic groups (NC = 0, MCI
Furthermore, all slopes are positive, indicating average
= 1, AD = 2), and also between our atrophy estimates
volume loss over time. This is consistent with expectaand the scores of the mini-mental state exam (MMSE).
tions for an aging and neuro-degenerative population. The
We applied the Benjamini-Hochberg procedure [36] for all
slopes capture increasing atrophy with disease severity. In
the correlation results in this paper to reduce the false
ADNI-1/ADNI-2, we expect SlopeNC-NC < SlopeMCI-NC <
discovery rate for multiple comparisons. The overall false
SlopeNC-MCI < SlopeMCI-AD < SlopeAD-AD and all six exdiscovery rate was set to be 0.01, which resulted in an efperimental groups (i.e. LDDMM-1, Pred-1, Pred+Corrfective significance level of α ≈ 0.0093. Detailed results
1, LDDMM-2, Pred-2, and Pred+Corr-2) are generally
can be found in Table 7 and Fig. 6, respectively. In deconsistent with this expectation. Exceptions happen in
tail, for ADNI-1/2, we randomly selected 2009 cases from
ADNI-2 for the NC-MCI and MCI-NC cases. As the numeach diagnostic category at each month and calculated the
ber of subjects involved is relatively small, i.e., fewer than
Spearman rank-order correlation. Fig. 6 shows the re20, compared with the other cases (roughly 100), one may
sults for 50 repetitions. We observe median correlations
speculate that this observation is caused by the limited
for all four prediction models in the range of −0.36 to
number of data points for NC-MCI and MCI-NC as shown
−0.75 for MMSE and 0.36 to 0.65 for diagnostic category.
in the #data column of Table 6. However, the behavior
The correlations for all four prediction+correction modwithin each starting diagnostic category, is consistent, i.e.,
els were in the range of −0.40 to −0.75 for MMSE and
for NC SlopeNC-NC < SlopeNC-MCI and for MCI SlopeMCI-NC <0.36 to 0.65 for diagnostic category. Previous studies reSlopeMCI-MCI < SlopeMCI-AD . Hence, all six groups’ slope
ported Pearson correlations between comparable atrophy
results in ADNI-1/ADNI-2 are generally consistent with
estimates and clinical variables as high as −0.7 for MMSE
our expectation (and also consistent with results in [33]).
and 0.5 for diagnostic category for 100 subjects[31, 32].Our
The slope estimated from the prediction+correction retwo optimization-based LDDMM results achieve median
sults is larger than the slope estimated from the prediction
correlations ranging from −0.40 to −0.76 for MMSE and
model results and closer to the slope obtained from the
0.40 to 0.66 for diagnostic category, which is very similar
optimization-based LDDMM results. This indicates that
to the predction+correlation models. In general, the corthe correction network can improve prediction accuracy.
rection+prediction FPSGR models outperform the models
Fig. 5 shows linear regression results for the estimated
using only the prediction network. Further, using the coratrophy scores in ADNI-1/2 for the Pred+Corr-1 model.
rection network, FPSGR achieved comparable and someBoth the data points themselves (i.e., the atrophy scores),
times even slightly better performance compared to the
as well as kernel density estimates for the linear trends
optimization-based LDDMM SGR method, see Table 7 for
for each subject are shown. These results are consistent
additional quantitative results. Specifically, FPSGR using
with the results of Table 6 discussed above. We conclude
the prediction+correction network performs best in 8 out
that (1) neither LDDMM optimization nor FPSGR produced deformations with significant bias to overestimate or
8 We used Spearman rank-order correlation instead of Pearson corunderestimate volume change; (2) a linear model of atrorelation,
because the diagnostic groups imply an ordering only.
phy scores generated by FPSGR can capture intrinsic vol9 In ADNI-1 48 month, the number was 60 because there was not
ume change (i.e., slope) among different diagnostic change
enough data; ADNI-2 36 month was omitted due to lack of data.
9

Figure 5: Linear regression of atrophy scores with respect to time for different diagnostic changes of ADNI-1 Pred+Corr-1 and ADNI-2
Pred+Corr-1. Red line is the estimated regression line. green curves are the lower and upper bounds of the 95% confidence interval. Blue dots
indicate actual data points. Bright white / purple images indicate kernel density estimations for all real data points illustrating dominant
longitudinal trends in the data.

10

Figure 6: Boxplot of FPSGR-derived correlations with clinical variables in ADNI-1 and ADNI-2. Prediction results are comparable with
optimization-based LDDMM. Adding the correction network generally improves prediction results.

Figure 7: Boxplot of Spearman rank-order correlations between atrophy measures and MMSE with respect to time in ADNI-1 and ADNI-2.
Top row: ADNI-1 NC-NC group (left), ADNI-1 MCI-MCI group (middle), ADNI-1 AD-AD group (right). Bottom row: ADNI-2 NCNC group (left), ADNI-2 MCI-MCI group (middle), ADNI-2 AD-AD group (right). ADNI-1 MCI-MCI and ADNI-1 AD-AD show stronger
correlations with time. In comparison, correlations remain relatively stable over time for the diagnostic groups in ADNI-2.

11

Figure 8: Kernel density estimates of highest 10% local correlations of atrophy with MMSE within the ROI depicted in Fig. 3. Top row:
results of NC group, MCI group and AD group from ADNI-1. Bottom row: results of NC group, MCI group and AD group from ADNI-2.
Results show a shifting pattern for the ADNI-1 MCI case, the ADNI-1 AD case and the ADNI-2 AD case.

of 18 comparisons for MMSE and in 12 out of 20 comparisons for diagnostic group. In the cases where FPSGR with
prediction+correction network did not perform best its difference to the best method was generally very small. In
general FPSGR using the correction network performs better than FPSGR without the correction network. To check
for statistical differences in the performance of FPSGR, we
use a paired t-test. Table 8 shows the resulting p-values
for the three methods: optimization-based SGR (i.e., LDDMM), FPSGR without correction network (i.e., Pred)
and FPSGR with correction network (i.e., Pred+Corr).
In both correlation with MMSE and DX, FPSGR with
correction network shows significantly better performance
than LDDMM and FPSGR without correction network,
which justifies the use of the FPSGR method. In summary, FPSGR captures correlations between atrophy and
clinical measures well.
To further explore the correlations of atrophy with
MMSE scores, we visualize them separated by diagnostic group where diagnosis did not change (i.e., NC-NC,
MCI-MCI, AD-AD) in Fig. 7. For the ADNI-1 dataset, we
observe (as expected) very low correlations for the normal
diagnostic group (with no clear trend), and much stronger
correlations for the MCI and AD groups. MCI and AD
also exhibit increasingly stronger correlations with time.
In case of ADNI-2, the MCI group shows modest correlations, which remain consistent across time. Correlations
are relatively low for the normal groups. The AD groups
show increasingly strong correlations over time. In con-

trast to ADNI-1, ADNI-2 focuses mainly on earlier stages
of the diagnostic groups [34]. Hence, the deformations in
ADNI-2 are generally smaller than in ADNI-1. This may explain why the NC and MCI diagnostic groups show consistent correlation values over time (instead of stronger correlations as for AD in ADNI-2 or the MCI and AD groups
in ADNI-1).
To address the question how stat-ROI specific measures
behave over time, we explore how atrophy locally (i.e.,
voxel-by-voxel) correlates with MMSE. The local atrophy
is defined as
(1 − det(Dφ(x))) × 100 .
I.e., each voxel in a stat-ROI has an associated atrophy
score. Fig. 8 shows kernel density estimates of the highest 10% local correlations in a violin plot. For the ADNI-1
MCI and AD groups, a clear shift toward stronger correlations can be observed over time, consistent with the
boxplots of Fig. 7. This indicates the progression of the
disease. Correlations for the normal groups in ADNI 1/2
are mostly centered around a modest correlation (as expected). In ADNI-2, only the AD diagnostic group shows
a shift towards stronger correlations over time. All the
other diagnostic groups show a relatively consistent distribution over time. This is also consistent with Fig. 7.
4.4. Justification of SGR
For simple geodesic regression to be a useful model
it should outperform pairwise image registration. The
12

ADNI-1

NC-NC

NC-MCI

MCI-MCI

MCI-NC

MCI-AD

AD-AD

LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2

ADNI-2

NC-NC

NC-MCI

MCI-MCI

MCI-NC

MCI-AD

AD-AD

LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2

Slope
0.70,
0.44,
0.68,
0.66,
0.50,
0.58,
0.94,
0.58,
0.90,
1.19,
0.99,
1.07,
1.17,
0.80,
1.09,
1.00,
0.82,
0.90,
0.72,
0.44,
0.68,
0.79,
0.61,
0.70,
2.10,
1.40,
1.84,
1.92,
1.56,
1.64,
2.33,
1.50,
2.05,
2.28,
1.85,
1.95,
Slope
[0.55, 0.65,
[0.41, 0.48,
[0.50, 0.57,
[0.51, 0.62,
[0.47, 0.55,
[0.35, 0.44,
[0.56, 0.79,
[0.53, 0.68,
[0.63, 0.80,
[0.62, 0.90,
[0.58, 0.77,
[0.46, 0.68,
[0.71, 0.83,
[0.53, 0.61,
[0.64, 0.73,
[0.71, 0.82,
[0.58, 0.66,
[0.50, 0.59,
[0.03, 0.39,
[0.05, 0.29,
[0.08, 0.36,
[0.14, 0.40,
[0.24, 0.42,
[0.05, 0.26,
[1.65, 1.95,
[1.09, 1.27,
[1.39, 1.62,
[1.59, 1.91,
[1.15, 1.35,
[1.20, 1.45,
[2.49, 2.76,
[1.74, 1.90,
[2.14, 2.34,
[2.72, 2.99,
[1.97, 2.14,
[2.16, 2.36,
[0.62,
[0.37,
[0.61,
[0.57,
[0.43,
[0.51,
[0.72,
[0.39,
[0.71,
[0.88,
[0.72,
[0.80,
[0.97,
[0.65,
[0.92,
[0.83,
[0.69,
[0.77,
[0.48,
[0.26,
[0.51,
[0.54,
[0.40,
[0.49,
[1.94,
[1.28,
[1.70,
[1.75,
[1.42,
[1.49,
[1.97,
[1.23,
[1.74,
[1.92,
[1.56,
[1.65,

0.78]
0.50]
0.75]
0.75]
0.57]
0.65]
1.16]
0.78]
1.10]
1.50]
1.26]
1.34]
1.38]
0.96]
1.26]
1.17]
0.96]
1.04]
0.96]
0.62]
0.86]
1.03]
0.83]
0.91]
2.27]
1.53]
1.98]
2.09]
1.70]
1.78]
2.69]
1.77]
2.35]
2.65]
2.15]
2.24]
0.75]
0.55]
0.65]
0.72]
0.62]
0.52]
1.02]
0.82]
0.97]
1.18]
0.97]
0.91]
0.94]
0.68]
0.82]
0.92]
0.73]
0.67]
0.74]
0.52]
0.64]
0.67]
0.61]
0.48]
2.25]
1.46]
1.85]
2.23]
1.56]
1.69]
3.04]
2.07]
2.54]
3.27]
2.31]
2.56]

Intercept
[-0.25,-0.08, 0.09]
[-0.21, -0.08, 0.05]
[-0.15, -0.01, 0.13]
[-0.21, -0.04, 0.14]
[-0.16, -0.02, 0.11]
[-0.12, 0.01, 0.15]
[-0.45, -0.03, 0.39]
[-0.43, -0.05, 0.33]
[-0.40, -0.01, 0.37]
[-0.65, -0.05, 0.55]
[-0.68, -0.16, 0.36]
[-0.66, -0.14, 0.38]
[-0.28, 0.05, 0.39]
[-0.29, -0.03, 0.22]
[-0.14, 0.14, 0.42]
[-0.21, 0.06, 0.33]
[-0.20, 0.02, 0.24]
[-0.15, 0.07, 0.29]
[-0.85, -0.42, 0.01]
[-0.61, -0.29, 0.03]
[-0.52, -0.20, 0.13]
[-0.79, -0.36, 0.07]
[-0.62, -0.24, 0.14]
[-0.59, -0.21, 0.17]
[-0.28, 0.02, 0.31]
[-0.24, -0.02, 0.20]
[-0.17, 0.08, 0.33]
[-0.16, 0.14, 0.44]
[-0.11, 0.14, 0.39]
[-0.08, 0.17, 0.43]
[-0.17, 0.27, 0.70]
[-0.13, 0.21, 0.54]
[-0.04, 0.33, 0.70]
[-0.20, 0.24, 0.68]
[-0.13, 0.22, 0.57]
[-0.10, 0.25, 0.60]
Intercept
[-0.08, 0.03, 0.13]
[-0.03, 0.04, 0.12]
[-0.04, 0.05, 0.13]
[-0.10, 0.01, 0.12]
[-0.03, 0.05, 0.13]
[-0.09, -0.00, 0.08]
[-0.22, 0.01, 0.25]
[-0.14, 0.01, 0.16]
[-0.16, 0.02, 0.19]
[-0.32, -0.02, 0.28]
[-0.19, 0.01, 0.22]
[-0.25, -0.02, 0.22]
[-0.13, -0.00, 0.12]
[-0.06, 0.02, 0.10]
[-0.08, 0.02, 0.11]
[-0.14, -0.02, 0.09]
[-0.05, 0.03, 0.12]
[-0.12, -0.02, 0.07]
[-0.38, 0.05, 0.47]
[-0.24, 0.05, 0.33]
[-0.28, 0.05, 0.38]
[-0.28, 0.04, 0.35]
[-0.17, 0.05, 0.28]
[-0.22, 0.03, 0.29]
[-0.21, 0.13, 0.47]
[-0.12, 0.09, 0.30]
[-0.15, 0.11, 0.37]
[-0.16, 0.19, 0.53]
[-0.09, 0.14, 0.36]
[-0.13, 0.14, 0.41]
[-0.15, 0.07, 0.30]
[-0.09, 0.04, 0.18]
[-0.09, 0.08, 0.24]
[-0.15, 0.07, 0.29]
[-0.07, 0.07, 0.21]
[-0.15, 0.02, 0.18]

#data

ADNI-1

154
6mo

156

24
12mo

22

146
18mo

148

16
24mo

17

148
36mo

147

143
48mo

140

LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2

ADNI-2

170
3mo

175

16
6mo

17

184
12mo

183

16
24mo

21

70
36mo

65

101

LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2

MMSE
-0.4957
-0.4642
-0.5104
-0.4667
-0.4711
-0.4734
-0.5749
-0.5328
-0.5799
-0.5301
-0.5351
-0.5374
-0.4939
-0.4659
-0.4924
-0.4385
-0.4389
-0.4384
-0.6064
-0.5664
-0.6001
-0.5822
-0.5911
-0.5898
-0.5142
-0.4731
-0.5069
-0.4334
-0.4425
-0.4393
-0.7456
-0.7294
-0.7443
-0.6889
-0.6995
-0.7005
MMSE
N/A
N/A
N/A
N/A
N/A
N/A
-0.4989
-0.4768
-0.5128
-0.5072
-0.4718
-0.5066
-0.4756
-0.4530
-0.4908
-0.4937
-0.4626
-0.4987
-0.4120
-0.3670
-0.4109
-0.4095
-0.3411
-0.3943
-0.2474
-0.2474
-0.2474
0.0935
0.0935
0.0935

p-value
5.17e-39
8.09e-34
1.22e-41
4.17e-34
8.48e-35
3.54e-35
5.23e-51
9.46e-43
4.39e-52
6.81e-42
9.79e-43
3.73e-43
4.86e-16
3.18e-14
6.16e-16
9.50e-13
9.06e-13
9.75e-13
5.01e-45
2.83e-38
6.55e-44
4.11e-40
1.41e-41
2.28e-41
4.29e-20
7.38e-17
1.71e-19
3.79e-13
1.07e-13
1.67e-13
2.01e-13
1.18e-12
2.30e-13
2.25e-10
9.08e-11
8.31e-11
p-value
N/A
N/A
N/A
N/A
N/A
N/A
8.01e-31
6.22e-28
9.64e-33
4.29e-32
2.02e-27
5.25e-32
1.43e-27
7.32e-25
1.67e-29
1.07e-29
7.94e-26
2.35e-30
9.53e-15
8.51e-12
1.15e-14
2.09e-14
3.46e-10
2.20e-13
0.55
0.55
0.55
0.83
0.83
0.83

DX
0.5140
0.4754
0.5259
0.4814
0.4849
0.4890
0.5313
0.4898
0.5406
0.5055
0.5120
0.5155
0.4776
0.4313
0.4643
0.4000
0.3818
0.3790
0.5978
0.5607
0.5943
0.5534
0.5714
0.5709
0.5300
0.4926
0.5296
0.4815
0.4894
0.4863
0.6635
0.6458
0.6575
0.5927
0.6048
0.6067
DX
0.4254
0.4142
0.4353
0.4409
0.4280
0.4445
0.4688
0.4625
0.4846
0.4883
0.4742
0.4913
0.4859
0.4771
0.5064
0.5026
0.4913
0.5149
0.4476
0.4331
0.4632
0.4375
0.3940
0.4336
0.2869
0.2869
0.2869
0.1695
0.1695
0.1695

p-value
2.66e-42
1.30e-35
1.53e-44
1.75e-36
4.58e-37
9.67e-38
1.81e-42
1.97e-35
3.44e-44
1.17e-37
1.11e-38
2.89e-39
5.76e-15
3.37e-12
3.98e-14
1.12e-10
8.80e-10
1.19e-9
1.69e-43
2.18e-37
6.82e-43
1.24e-35
2.26e-38
2.65e-38
1.81e-21
2.42e-18
1.99e-21
2.93e-16
7.99e-17
1.34e-16
5.20e-10
2.08e-9
8.43e-10
1.98e-7
9.49e-8
8.49e-8
p-value
2.34e-24
4.72e-23
1.52e-25
2.77e-26
1.05e-24
9.64e-27
6.09e-27
3.47e-26
6.19e-29
1.58e-29
9.96e-28
6.33e-30
7.22e-29
9.39e-28
1.37e-31
7.05e-31
2.21e-29
1.44e-32
2.06e-17
2.71e-16
1.09e-18
1.93e-16
2.29e-13
3.79e-16
0.49
0.49
0.49
0.69
0.69
0.69

#data
608

606

565

560

238

241

435

427

277

256

69

65
#data
522

523

468

470

464

461

325

321

8

8

Table 7: FPSGR-derived correlations with clinical variables,
compared to correlations with clinical variables for SGR using
optimization-based LDDMM. The #data column lists the number
of data points analyzed. Green indicates that FPSGR using the prediction+correction network shows the strongest correlations; Yellow
indicates that FPSGR using the prediction network alone shows the
strongest correlations; Red indicates that LDDMM SGR shows the
strongest correlations. The MMSE column lists correlations between
atrophy scores and the mini-mental state exam scores; the DX column lists correlations between atrophy score and diagnostic category. Finally, the p-value column(s) list the p-values for the nullhypothesis that there is no correlation. Benjamini-Hochberg procedure was employed to reduce the false discovery rate and Purple
highlight indicates statistically significant. FPSGR using the prediction+correction network generally improves performance over using
the prediction network alone and frequently even performs slightly
better than the SGR results obtained by optimization-based LDDMM.

103

Table 6: Slope and intercept values for simple linear regression
of volume change over time. Our notation for slope and intercept
indicate [lower bound of 95% C.I., point estimate, upper bound of
95% C.I.]. The interval of intercept estimates all contain zero. The
slope changes between the different diagnostic groups. The #data
column lists the number of data points analyzed.

main conceptual difference is that the regression model
will recover an average trend based on multiple image
time-points, i.e., the resulting regression geodesic will be
a compromise between all the measurements. In contrast,
for pairwise image registration (which can be seen as a
trivial case of geodesic regression with two images only)
13

MMSE
LDDMM
Pred
Pred+Corr
MMSE
LDDMM
Pred
Pred+Corr
DX
LDDMM
Pred
Pred+Corr
DX
LDDMM
Pred
Pred+Corr

Normality Test
LDDMM
Pred
N/A
0.1507
0.1507
N/A
0.5361
0.0183
Paired t-test
LDDMM
Pred
N/A
0.0005484
0.9994516
N/A
0.0530827 0.0000282
Normality Test
LDDMM
Pred
N/A
0.1963
0.1963
N/A
0.2356
0.3208
Paired T -test
LDDMM
Pred
N/A
0.0010944
0.9989056
N/A
0.0186418 0.0000131

ADNI-1

Pred+Corr
0.5361
0.0183
N/A

NC-NC

NC-MCI

Pred+Corr
0.09469173
0.9999718
N/A

MCI-MCI

MCI-NC

MCI-AD

Pred+Corr
0.2356
0.3208
N/A

AD-AD

SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise

Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2

SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise

Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2

ADNI-2
NC-NC

Pred+Corr
0.9813582
0.9999869
N/A

NC-MCI

MCI-MCI

Table 8: Results of a Shapiro-Wilk normality test and a paired t-test
on MMSE and DX correlations among optimization-based LDDMM,
FPSGR without prediction network and FPSGR with correction network. The null-hypothesis for the Shapiro-Wilk normality test is that
the difference between column-method and row-method is normally
distributed. The null-hypothesis for the paired t-test is that the
column-method is statistically better than row-method. Green highlighted p-values indicate no rejection of the normality hypothesis (at
5% significance) and thus facilitate the paired t-test. p-values highlighted in red indicate a rejection of the normality null-hypothesis
and consequently do not allow a paired t-test.

MCI-NC

MCI-AD

AD-AD

Slope
0.44, 0.50]
0.52, 0.60]
0.50, 0.57]
0.57, 0.65]
0.58, 0.78]
0.63, 0.87]
0.99, 1.26]
0.96, 1.27]
0.80, 0.96]
0.86, 1.03]
0.82, 0.96]
0.85, 1.01]
0.44, 0.62]
0.45, 0.68]
0.61, 0.83]
0.56, 0.83]
1.40, 1.53]
1.42, 1.56]
1.56, 1.70]
1.60, 1.75]
1.50, 1.77]
1.55, 1.85]
1.85, 2.15]
1.85, 2.16]
Slope
[0.41, 0.48, 0.55]
[0.25, 0.33, 0.41]
[0.47, 0.55, 0.62]
[0.26, 0.35, 0.44]
[0.53, 0.68, 0.82]
[0.37, 0.57, 0.77]
[0.58, 0.77, 0.97]
[0.42, 0.65, 0.88]
[0.53, 0.61, 0.68]
[0.43, 0.52, 0.61]
[0.58, 0.66, 0.73]
[0.45, 0.54, 0.63]
[0.05, 0.29, 0.52]
[-0.10, 0.17, 0.45]
[0.24, 0.42, 0.61]
[0.03, 0.26, 0.49]
[1.09, 1.27, 1.46]
[0.88, 1.10, 1.32]
[1.15, 1.35, 1.56]
[0.89, 1.13, 1.37]
[1.74, 1.90, 2.07]
[1.57, 1.77, 1.96]
[1.97, 2.14, 2.31]
[1.79, 1.99, 2.19]
[0.37,
[0.44,
[0.43,
[0.48,
[0.39,
[0.39,
[0.72,
[0.65,
[0.65,
[0.69,
[0.69,
[0.70,
[0.26,
[0.21,
[0.40,
[0.29,
[1.28,
[1.28,
[1.42,
[1.44,
[1.23,
[1.25,
[1.56,
[1.53,

Intercept
[-0.21, -0.08, 0.05]
[-0.46, -0.30, -0.14]
[-0.16, -0.02, 0.11]
[-0.34, -0.18, -0.01]
[-0.43, -0.05, 0.33]
[-0.63, -0.16, 0.30]
[-0.68, -0.16, 0.36]
[-0.69, -0.10, 0.50]
[-0.29, -0.03, 0.22]
[-0.43, -0.15, 0.12]
[-0.20, 0.02, 0.24]
[-0.29, -0.04, 0.21]
[-0.61, -0.29, 0.03]
[-0.74, -0.31, 0.12]
[-0.62, -0.24, 0.14]
[-0.61, -0.14, 0.34]
[-0.24, -0.02, 0.20]
[-0.31, -0.06, 0.19]
[-0.11, 0.14, 0.39]
[-0.22, 0.06, 0.33]
[-0.13, 0.21, 0.54]
[-0.23, 0.13, 0.49]
[-0.13, 0.22, 0.57]
[-0.15, 0.23, 0.60]
Intercept
[-0.03, 0.04, 0.12]
[0.15, 0.24, 0.33]
[-0.03, 0.05, 0.13]
[0.22, 0.32, 0.43]
[-0.14, 0.01, 0.16]
[-0.06, 0.14, 0.33]
[-0.19, 0.01, 0.22]
[-0.07, 0.18, 0.42]
[-0.06, 0.02, 0.10]
[0.04, 0.14, 0.23]
[-0.05, 0.03, 0.12]
[0.09, 0.19, 0.29]
[-0.24, 0.05, 0.33]
[-0.12, 0.21, 0.53]
[-0.17, 0.05, 0.28]
[0.02, 0.29, 0.57]
[-0.12, 0.09, 0.30]
[0.08, 0.33, 0.58]
[-0.09, 0.14, 0.36]
[0.18, 0.44, 0.70]
[-0.09, 0.04, 0.18]
[0.01, 0.17, 0.34]
[-0.07, 0.07, 0.21]
[0.05, 0.21, 0.37]

Table 9: SGR prediction model compared with a pairwise prediction model. Slope and intercept values for simple linear regression
of volume change over time. The notation for slope and intercept
columns indicates [Lower bound of 95% C.I., point estimate, Upper bound of 95% C.I.]. Green indicates that the intercept is closer
to zero (also, zero is within the 95% confidence interval) for SGR prediction model; Yellow indicates that the intercept is closer to zero for
pairwise prediction model; Red indicates that the point estimate is
either biased to overestimate or underestimate volume change. The
SGR prediction model performs better than the pairwise prediction
model.

the deformation will in general be able to match the target image well. However, just as in linear regression, this
may accentuate the effects of noise. In both setups, images
can be interpolated or extrapolated based on the estimated
geodesic.
Tables 9 and 10 justify the use of SGR. Specifically,
Table 9 shows linear regression results of atrophy measures over time as obtained via SGR (i.e., using an SGR
fit over all time-points followed by atrophy computations
based on the deformations of the regression geodesic) compared with atrophy measures obtained by pairwise registration. For both the ADNI-1 and the ADNI-2 datasets,
SGR outperforms the pairwise registration approach in
two aspects: (1) the estimated intercept of SGR is generally closer to zero than for the pairwise method and the
intercept 95% confidence interval is narrower; (2) 11 out of
24 of the 95% confidence intervals of the pairwise methods
show bias to either overestimate or underestimate volume
change, while none of the SGR results show such significant bias. Table 10 compares the correlations between
atrophy and clinical measures (MMSE and diagnostic category) of SGR and the pairwise approach. SGR performs
better than the pairwise approach in 13 out of 18 cases for
MMSE and in 15 out of 20 cases for the diagnostic cat-

egory. Furthermore, when the pairwise method is better
than SGR, the difference is much smaller compared to the
differences observed for the cases where SGR is better than
the pairwise method. Also note that the pairwise method
shows better performance in later months compared to
earlier months. This could, for example, be because the
deformations are larger for later time-points and hence the
registration result becomes more stable, or because SGR is
also heavily influenced by the last time-point. To address
the above observation, we used a Shapiro-Wilk normality test and a Wilcoxon signed-rank test. From Table 11
we see that we can reject the null-hypothesis of normality and hence, a paired t-test is not appropriate. As an
14

ADNI-1
6mo

12mo

18mo

24mo

36mo

48mo

SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise

Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2

ADNI-2
3mo

6mo

12mo

24mo

36mo

SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise
SGR
Pairwise

Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2
Pred-1
Pred-1
Pred-2
Pred-2

MMSE
-0.4642
-0.3138
-0.4711
-0.3431
-0.5328
-0.4393
-0.5351
-0.4465
-0.4659
-0.4164
-0.4389
-0.4078
-0.5664
-0.5805
-0.5911
-0.5927
-0.4731
-0.4470
-0.4425
-0.4538
-0.7294
-0.7100
-0.6995
-0.6709
MMSE
N/A
N/A
N/A
N/A
-0.4768
-0.3378
-0.4718
-0.3312
-0.4530
-0.4305
-0.4626
-0.4223
-0.3670
-0.3772
-0.3411
-0.3517
-0.2474
-0.1650
0.0935
0.0935

p-value
8.09e-34
2.31e-15
8.48e-35
3.51e-18
9.46e-43
4.67e-28
9.79e-43
9.61e-29
3.18e-14
2.12e-11
9.06e-13
4.52e-11
2.83e-38
1.51e-40
1.41e-41
7.34e-42
7.38e-17
5.20e-15
1.07e-13
2.08e-14
1.18e-12
8.43e-12
9.08e-11
9.65e-10
p-value
N/A
N/A
N/A
N/A
6.22e-28
5.93e-14
2.02e-27
1.70e-13
7.32e-25
2.34e-22
7.94e-26
2.30e-21
8.51e-12
3.06e-12
3.46e-10
8.89e-11
0.55
0.70
0.83
0.83

DX
0.4754
0.3369
0.4849
0.3680
0.4898
0.3996
0.5120
0.4154
0.4313
0.3882
0.3818
0.3356
0.5607
0.5791
0.5714
0.5811
0.4926
0.4798
0.4894
0.4990
0.6458
0.6168
0.6048
0.5924
DX
0.4142
0.1744
0.4280
0.1503
0.4625
0.2633
0.4742
0.2849
0.4771
0.4472
0.4913
0.4374
0.4331
0.4515
0.3940
0.4239
0.4536
0.2869
0.1695
0.2608

p-value
1.30e-35
1.32e-17
4.58e-37
7.24e-21
1.97e-35
4.51e-23
1.11e-38
1.00e-24
3.37e-12
5.56e-10
8.80e-10
9.38e-8
2.18e-37
2.55e-40
2.26e-38
6.26e-40
2.42e-18
2.36e-17
7.99e-17
1.59e-17
2.08e-9
1.67e-8
9.49e-8
2.01e-7
p-value
4.72e-23
6.17e-5
1.05e-24
5.64e-4
3.47e-26
7.29e-9
9.96e-28
3.14e-10
9.39e-28
3.40e-24
2.21e-29
5.72e-23
2.71e-16
9.99e-18
2.29e-13
1.99e-15
0.26
0.49
0.69
0.53

#data

ADNI-1

608
60mo

606

Forecast

565
560
72mo

Forecast

238
241
435

84mo

Forecast

427

p-value
1.34e-13
5.16e-11
2.48e-13
2.32e-10
1.77e-10
9.97e-11
1.60e-10
1.45e-8
1.47e-10
3.18e-7
1.09e-7
1.09e-7
1.40e-6
2.84e-5
1.83e-6
4.65e-5
7.81e-5
8.35e-5

DX
0.5157
0.4816
0.5240
0.4761
0.4620
0.4652
0.4507
0.4364
0.4667
0.4233
0.4273
0.4259
0.4108
0.4018
0.4211
0.3622
0.3629
0.3609

p-value
3.85e-13
1.98e-11
1.38e-13
1.43e-11
6.63e-11
4.73e-11
4.37e-10
1.75e-9
8.52e-11
2.15e-9
1.46e-9
1.67e-9
6.17e-7
1.13e-6
2.98e-7
6.53e-6
6.25e-6
7.12e-6

#data
173

180

174

184

137

147

Table 12: Correlations of forecasting results. The #data column lists the number of data points analyzed. Green indicates
that FPSGR using the prediction+correction network shows the
strongest correlations; Yellow indicates that FPSGR using the prediction network alone shows the strongest correlations; Red indicates
that LDDMM SGR shows the strongest correlations. The BenjaminiHochberg procedure was employed to reduce the false discovery rate
(FDR). The Purple highlight indicates statistically significant results
after correction for multiple comparisons.

256
69
65
#data
522
523
468
470

ADNI-1

464
461

Original

325

Forecast
321

Replace
8

36mo
Original

8

Forecast
Replace
Original
Forecast
Replace
48mo
Original
Forecast
Replace

Shapiro-Wilk normality test
0.01943
0.03286

MMSE
-0.5242
-0.4727
-0.5193
-0.4501
-0.4527
-0.4582
-0.4607
-0.4132
-0.4615
-0.3662
-0.3793
-0.3793
-0.3986
-0.3495
-0.3946
-0.3293
-0.3199
-0.3187

277

Table 10: SGR prediction model compared with pairwise prediction model. Results show correlations with clinical variables. The
#data column lists the number of data points analyzed. Green indicates a stronger correlation for the SGR prediction method; Yellow
indicates a stronger correlation for the pairwise model. The p-value
column lists p-values for the null-hypothesis that there is no correlation. The Benjamini-Hochberg procedure was employed to reduce
the false discovery rate (FDR). The Purple highlight indicates statistically significant results after correction for multiple comparisons.
In general, SGR prediction performs better than pairwise prediction demonstrating that regression stabilizes the correlation results.
ADNI-2 36mo only has 8 data points and the p-value is greater than
0.1, thus we ignore this month in our comparison.

MMSE
DX

LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2

LDDMM-1
Pred-1
Pred+Corr-1
Pred-1
Pred+Corr-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
Pred-2
Pred+Corr-2
Pred-2
Pred+Corr-2
LDDMM-1
Pred-1
Pred+Corr-1
Pred-1
Pred+Corr-1
Pred-1
Pred+Corr-1
LDDMM-2
Pred-2
Pred+Corr-2
Pred-2
Pred+Corr-2
Pred-2
Pred+Corr-2

MMSE
-0.5142
-0.4731
-0.5069
-0.4583
-0.4708
-0.4923
-0.5097
-0.4334
-0.4425
-0.4393
-0.4078
-0.4005
-0.4202
-0.4164
-0.7456
-0.7294
-0.7443
-0.6332
-0.6541
-0.6446
-0.6668
-0.6889
-0.6995
-0.7005
-0.6528
-0.6403
-0.6334
-0.6307

p-value
4.29e-20
7.38e-17
1.71e-19
1.09e-15
1.42e-16
3.43e-18
1.37e-19
3.79e-13
1.07e-13
1.67e-13
1.36e-11
3.34e-11
2.75e-12
4.51e-12
2.01e-13
1.18e-12
2.30e-13
5.29e-9
1.10e-9
2.27e-9
3.98e-10
2.25e-10
9.08e-11
8.31e-11
3.79e-9
9.25e-9
1.49e-8
1.79e-8

DX
0.5300
0.4926
0.5296
0.4825
0.4980
0.5104
0.5375
0.4815
0.4894
0.4863
0.4398
0.4301
0.4635
0.4582
0.6635
0.6458
0.6575
0.6165
0.6317
0.6478
0.6800
0.5927
0.6048
0.6067
0.5568
0.5460
0.5970
0.5973

p-value
1.81e-21
2.42e-18
1.99e-21
1.93e-17
1.21e-18
1.21e-19
5.47e-22
2.93e-16
7.99e-17
1.34e-16
1.95e-13
7.40e-13
6.27e-15
1.38e-14
5.20e-10
2.08e-9
8.43e-10
1.70e-8
5.86e-9
1.78e-9
1.31e-10
1.98e-7
9.49e-8
8.49e-8
1.46e-6
2.55e-6
1.53e-7
1.50e-7

#data

277

256

69

65

Table 13: Forecast results compared with real data results. The
#data column lists the number of data points analyzed. The
Benjamini-Hochberg procedure was employed to reduce the false discovery rate (FDR). Purple highlight indicates statistically significant
results after corrections for multiple comparisons. Forecast results
are calculated by using SGR excluding 36mo and 48mo data points
and then predicting 36mo and 48mo correlations. Results are compared based on the same dataset except for two invalid data points
for the 36mo data.

Wilcoxon signed-rank test
0.0005226
0.0005083

Table 11: p-values for a Shapiro-Wilk normality test and Wilcoxon
signed-rank test on MMSE and DX correlations between the SGR
prediction model and the pairwise prediction model. The nullhypothesis for the Shapiro-Wilk normality test is that the difference
of two methods is normally distributed (at a significance level of
5%). The null-hypothesis for the Wilcoxon signed-rank test is that
the pairwise prediction method is statistically better than the SGR
prediction method (at a significance level of 5%).

method and in general also results in stronger correlations.
alternative, we conducted a Wilcoxon signed-rank test to
compare the SGR prediction model and the pairwise prediction model. Table 11 shows that the SGR prediction
model is statistically significantly better than the pairwise
prediction model. Based on the above points, we conclude that SGR is more stable over time than the pairwise

4.5. Forecasting
Another interesting question for SGR and geodesic regression in general is the suitability of the model for the
data. To address this question, we evaluate if SGR can
forecast unseen future time-points. Specifically we consider this question in two different scenarios:
15

Figure 9: Comparison of correlations among prediction results, Forecast results (-F) and Replace results (-C) in MCI converter groups
(MCI-NC, MCI-MCI and MCI-AD). In Pred-1, Forecast results outperform Replace results; In Pred-2, Forecast results and Replace results
are comparable.

Q1) Extrapolate-clinical: Can we extrapolate the SGR
results into the future (to time-points that do not exist in the ADNI image dataset, but for the clinical
data) while still obtaining strong correlations.

the closest measured time-point. For example, if we have
images at 6, 12, and 18 month, but we want to forecast
at 24 month, we use the 18 month image as the imputed
24 month image and then perform SGR on the 6, 12, 18,
and the imputed 24 month images. We then obtain the
deformation at 24 months from the SGR result.

Q2) Extrapolate-image: How well can correlations between atrophy and clinical measures be predicted for
time-points when we do or do not use image data at
that very time-point. We artificially leave out image
measurements so that we can compare prediction results to results when we have the image measurement.

ad Q1. Table 12 shows correlations between atrophy
and the clinical measures for the Forecast results for 60
month, 72 month and 84 month.The resulting correlations
of atrophy with diagnostic category are all above 0.3 (or
below -0.3). Furthermore, the Forecast correlations show
a downward trend with respect to time, which means that
the prediction of “far-away“ points is not as accurate as
for the “near” future. On the other hand, SGR using the
6 month to 48 month time points results in correlation
around -0.5 for MMSE and 0.5 for DX on average. Hence
the correlation with the diagnostic category is consistent
for that of 60 months. In other words, using 6 month
to 48 month data, our prediction model can predict accu-

For both scenarios we use two different forecasting approaches. In the first approach (Forecast) we simply compute SGR results with the available image time-points and
then extrapolate using the resulting regression geodesic to
the desired time-point in the future. In the second approach (Replace), we artificially impute the missing image time-points by simply replacing them by the image at
16

rately up to 60 month. Our prediction+correction network
performs as well as and even slightly better than SGR using optimization-based LDDMM. Fig. 10 shows that these
forecasting results capture the trends of the changes in the
temporal lobes near the hippocampus and changes in the
ventricles.
ad Q2. Table 13 and Fig. 9 show Forecast and Replace results for correlations between atrophy and clinical
measures in comparison to using all images. Specifically,
for the Forecast and Replace results we did not use the
available images at 36 and 48 month so we could compare against the results obtained when using these images. If FPSGR is a good model, it should results in
correlation results as close to the correlation results using all images as possible. The Forecast correlations are
only slightly weaker (0.02 to 0.05 lower) than the original
correlations using all images illustrating that FPSGR can
approximately forecast future changes.
The overall correlations in Table 13 show that the Replace group performs better than the Forecast group. In
particular, we are also interested in the prediction of MCI
converters, namely, MCI to NC, MCI to MCI, and MCI to
AD. The boxplots in Fig. 9 show the correlations for such
predictions. The Replace group in Fig. 9 show relatively
worse correlation performance than the Forecast group
in ADNI-1 Pred-1 and consistent performance in ADNI1 Pred-2. Hence SGR on a longitudinal image data can
achieves good forecasting result for MCI converters.
Thus, both Extrapolate-clinical and Extrapolateimage experiments justify the use of FPSGR in predicting
near future longitudinal trends especially for MCI converters.

Figure 10:
Average Jacobian determinant over time and
diagnostic category for ADNI-1 Prediction-1 and ADNI-1
Prediction+Correction-1 (experiments in ADNI-2 show similar
results). A value < 1 means shrinkage and value > 1 means
expansion. The 60 month - 84 month results contained in the purple
rectangle are forecasts using the data from 6 month - 48 month.
Results show consistent volume loss over time near the temporal
lobes and expansion over time near the ventricles/cerebrospinal
fluid.

4.6. Jacobian Determinant (JD)
The average JD images qualitatively agree with prior
results [33, 34]: severity of volume change increases with
severity of diagnosis and time. Change is most substantial
in the temporal lobes near the hippocampus (see Fig. 10).
In Fig. 10, 6 month to 48 month are existing data points;
60 month to 84 month are forecast results. Blue indicates
volume loss. Red indicates expansion. Results are consistent with expectations: volume loss increases with time
and severity of diagnosis in temporal lobes; volume expansion increases with respect to time and severity of diagnosis
around the ventricles / cerebrospinal fluid. The forecast
results capture visually sensible volume loss or expansion
over time, qualitatively illustrating the performance of our
method.

5. Conclusion and future work
In this work, we proposed a fast approach for geodesic regression (FPSGR) to study longitudinal image data. FPSGR
incorporates the recently proposed FPIR [19, 20] into the
SGR [15] framework, thus leading to a computationally
efficient solution to geodesic regression. Since FPSGR replaces the computationally intensive intermediate step of
computing pairwise initial momenta via a deep-learning
prediction method, it is orders of magnitude faster than existing approaches [15, 28], without compromising accuracy.
Consequently, FPSGR facilitates the analysis of large-scale
imaging studies. Experiments on the ADNI-1/ADNI-2 datasets
demonstrate that FPSGR captures expected atrophy trends
of normal aging, MCI and AD. It further (1) exhibits negligible bias towards volume changes within stat-ROIs, (2)
17

shows high correlations with clinical variables (MMSE and
diagnosis) and (3) produces consistent forecasting results
on unseen data.
In future work, it will be be interesting to explore
FPSGR for the task of classifying stable Mild Cognitive
Impairment (sMCI) and progressive Mild Cognitive Impairment (pMCI). Currently, FPSGR only shows modest accuracy for distinguishing these types within MCI.
Extending our approach to time-warped geodesic regression models [9] might improve the accuracy in this context. Furthermore, end-to-end prediction of averaged initial momenta would be an interesting future direction, as
this would allow learning representations that characterize
the geodesic path among multiple time-series images, not
only based on averages of momenta for two images as in
FPIR [19, 20].

oratory for Neuro Imaging at the University of Southern
California.
References

[1] C. R. Jack, J. Barnes, M. A. Bernstein, B. J. Borowski,
J. Brewer, S. Clegg, A. M. Dale, O. Carmichael, C. Ching,
C. DeCarli, et al., Magnetic resonance imaging in ADNI 2,
Alzheimer’s & Dementia 11 (7) (2015) 740–756.
[2] M. A. Ikram, A. van der Lugt, W. J. Niessen, P. J. Koudstaal, G. P. Krestin, A. Hofman, D. Bos, M. W. Vernooij, The
Rotterdam scan study: design update 2016 and main findings,
European journal of epidemiology 30 (12) (2015) 1299–1315.
[3] Biobank website: www.ukbiobank.ac.uk.
[4] M. Niethammer, Y. Huang, F.-X. Vialard, Geodesic regression
for image time-series, in: MICCAI, 2011, pp. 655–662.
[5] Y. Hong, Y. Shi, M. Styner, M. Sanchez, M. Niethammer, Simple geodesic regression for image time-series., in: WBIR, Vol. 12,
Springer, 2012, pp. 11–20.
[6] Y. Hong, S. Joshi, M. Sanchez, M. Styner, M. Niethammer,
Metamorphic geodesic regression, Medical Image Computing
Acknowledgements
and Computer-Assisted Intervention–MICCAI (2012) 197–205.
[7] N. Singh, J. Hinkle, S. Joshi, P. T. Fletcher, A vector momenta
Research reported in this publication was supported
formulation of diffeomorphisms for improved geodesic regression
and atlas construction, in: ISBI, 2013, pp. 1219–1222.
by the National Institutes of Health (NIH) and the Na[8] P. T. Fletcher, Geodesic regression and the theory of least
tional Science Foundation (NSF) under award numbers
squares on Riemannian manifolds, IJCV 105 (2) (2013) 171–
NIH R01AR072013, NSF ECCS-1148870, and EECS-1711776.
185.
The content is solely the responsibility of the authors and
[9] Y. Hong, N. Singh, R. Kwitt, M. Niethammer, Timewarped geodesic regression, in: Medical image computing and
does not necessarily represent the official views of the NIH
computer-assisted intervention–MICCAI, Vol. 17, 2014, p. 105.
or the NSF. We also thank Nvidia for the donation of a
[10] N. Singh, M. Niethammer, Splines for diffeomorphic image reTitanX GPU.
gression, in: Medical image computing and computer-assisted
intervention – MICCAI, Vol. 17, 2014, p. 121.
Data collection and sharing for this project was funded
by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) [11] Y. Hong, R. Kwitt, N. Singh, B. Davis, N. Vasconcelos, M. Niethammer, Geodesic regression on the Grassmannian, in: Euro(National Institutes of Health Grant U01 AG024904) and
pean Conference on Computer Vision, Springer, 2014, pp. 632–
DOD ADNI (Department of Defense award number W81XWH646.
[12] N. Singh, F.-X. Vialard, M. Niethammer, Splines for diffeomor12-2-0012). ADNI is funded by the National Institute on
phisms, Medical image analysis 25 (1) (2015) 56–71.
Aging, the National Institute of Biomedical Imaging and
[13] Y. Hong, R. Kwitt, N. Singh, N. Vasconcelos, M. NiethamBioengineering, and through generous contributions from
mer, Parametric regression on the Grassmannian, IEEE transthe following: AbbVie, Alzheimer’s Association; Alzheimer’s
actions on pattern analysis and machine intelligence 38 (11)
(2016) 2284–2297.
Drug Discovery Foundation; Araclon Biotech; BioClinica,
[14]
M.
F. Beg, M. I. Miller, A. Trouvé, L. Younes, Computing large
Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir,
deformation metric mappings via geodesic flows of diffeomorInc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli
phisms, International journal of computer vision 61 (2) (2005)
Lilly and Company; EuroImmun; F. Hoffmann-La Roche
139–157.
[15] Y. Hong, Y. Shi, M. Styner, M. Sanchez, M. Niethammer, SimLtd and its affiliated company Genentech, Inc.; Fujireple geodesic regression for image time-series, in: WBIR, 2012,
bio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Impp. 11–20.
munotherapy Research & Development, LLC.; Johnson &
[16] X. Cao, J. Yang, J. Zhang, D. Nie, M. Kim, Q. Wang, D. Shen,
Johnson Pharmaceutical Research & Development LLC.;
Deformable image registration based on similarity-steered CNN
regression, in: MICCAI, Springer, 2017.
Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diag[17] S. Miao, Z. J. Wang, Y. Zheng, R. Liao, Real-time 2D/3D regnostics, LLC.; NeuroRx Research; Neurotrack Technoloistration via CNN regression, in: Biomedical Imaging (ISBI),
gies; Novartis Pharmaceuticals Corporation; Pfizer Inc.;
2016 IEEE 13th International Symposium on, IEEE, 2016, pp.
1430–1434.
Piramal Imaging; Servier; Takeda Pharmaceutical Com[18]
H.
Sokooti, B. d. Vos, F. Berendsen, B. P. Lelieveldt, I. Isgum,
pany; and Transition Therapeutics. The Canadian InstiM. Staring, Nonrigid image registration using multi-scale 3D
tutes of Health Research is providing funds to support
convolutional neural networks, in: MICCAI, Springer, 2017.
ADNI clinical sites in Canada. Private sector contribu[19] X. Yang, R. Kwitt, M. Niethammer, Fast predictive image registration, in: International Workshop on Large-Scale Annotation
tions are facilitated by the Foundation for the National Inof Biomedical Data and Expert Label Synthesis, Springer, 2016,
stitutes of Health (www.fnih.org). The grantee organizapp. 48–57.
tion is the Northern California Institute for Research and
[20] X. Yang, R. Kwitt, M. Styner, M. Niethammer, Quicksilver:
Education, and the study is coordinated by the Alzheimer’s
Fast predictive image registration–a deep learning approach,
NeuroImage 158 (2017) 378–396.
Therapeutic Research Institute at the University of South[21] M. Zhang, R. Liao, A. V. Dalca, E. A. Turk, J. Luo, P. E. Grant,
ern California. ADNI data are disseminated by the LabP. Golland, Frequency diffeomorphisms for efficient image regis-

18

[22]

[23]

[24]

[25]
[26]

[27]

[28]

[29]

[30]
[31]

[32]

[33]

[34]

[35]

[36]

tration, in: International Conference on Information Processing
in Medical Imaging, Springer, 2017, pp. 559–570.
M. Zhang, P. T. Fletcher, Finite-dimensional Lie algebras for
fast diffeomorphic image registration, in: IPMI, 2015, pp. 249–
260.
A. Klein, J. Andersson, B. A. Ardekani, J. Ashburner,
B. Avants, M.-C. Chiang, G. E. Christensen, D. L. Collins,
J. Gee, P. Hellier, et al., Evaluation of 14 nonlinear deformation
algorithms applied to human brain MRI registration, Neuroimage 46 (3) (2009) 786–802.
A. Dosovitskiy, P. Fischer, E. Ilg, P. Hausser, C. Hazirbas,
V. Golkov, P. van der Smagt, D. Cremers, T. Brox, Flownet:
Learning optical flow with convolutional networks, in: ICCV,
2015, pp. 2758–2766.
Z. Liu, R. Yeh, X. Tang, Y. Liu, A. Agarwala, Video frame synthesis using deep voxel flow, arXiv preprint arXiv:1702.02463.
T. Schuster, L. Wolf, D. Gadot, Optical flow requires
multiple strategies (but only one network), arXiv preprint
arXiv:1611.05607.
B. D. de Vos, F. F. Berendsen, M. A. Viergever, M. Staring, I. Išgum, End-to-end unsupervised deformable image registration with a convolutional neural network, arXiv preprint
arXiv:1704.06065.
Y. Hong, P. Golland, M. Zhang, Fast geodesic regression for
population-based image analysis, in: International Conference
on Medical Image Computing and Computer-Assisted Intervention, Springer, 2017, pp. 317–325.
Z. Ding, G. Fleishman, X. Yang, P. Thompson, R. Kwitt, M. Niethammer, A. D. N. Initiative, et al., Fast predictive simple
geodesic regression, in: Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support,
Springer, 2017, pp. 267–275.
D. Kingma, J. Ba, Adam: A method for stochastic optimization,
arXiv preprint arXiv:1412.6980.
G. Fleishman, P. M. Thompson, The impact of matching functional on atrophy measurement from geodesic shooting in diffeomorphisms, in: ISBI, 2017.
G. Fleishman, P. M. Thompson, Adaptive gradient descent optimization of initial momenta for geodesic shooting in diffeomorphisms, in: ISBI, 2017.
X. Hua, D. P. Hibar, C. R. K. Ching, C. P. Boyle, P. Rajagopalan, B. Gutman, A. D. Leow, A. W. Toga, C. R. J.
Jr., D. J. Harvey, M. W. Weiner, P. M. Thompson, Unbiased
tensor-based morphometry:improved robustness & sample size
estimates for Alzheimer’s disease clinical trials, NeuroImage 66
(2013) 648–661.
X. Hua, C. R. K. Ching, A. Mezher, B. Gutman, D. P. Hibar,
P. Bhatt, A. D. Leow, C. R. J. Jr., M. Bernstein, M. W. Weiner,
P. M. Thompson, MRI-based brain atrophy rates in ADNI phase
2: acceleration and enrichment considerations for clinical trials,
Neurobiology of Aging 37 (2016) 26–37.
P. A. Yushkevich, B. B. Avants, S. R. Das, J. Pluta, M. Altinay, C. Craige, Bias in estimation of hippocampal atrophy
using deformation-based morphometry arises from asymmetric
global normalization: An illustration in ADNI 3T MRI data,
NeuroImage 50 (2) (2010) 434 – 445.
Y. Benjamini, Y. Hochberg, Controlling the false discovery rate:
a practical and powerful approach to multiple testing, Journal
of the royal statistical society. Series B (Methodological) (1995)
289–300.

19

