3
Diffeomorphic Metric Mapping of
High Angular Resolution Diffusion
Imaging based on Riemannian
Structure of Orientation Distribution
Functions
Due to the inter-subject anatomical variation, it is necessary to align ODF images of
different subjects into a common space so that group-level statistical inference can be
performed (see Figure 3.1). In this chapter, we propose a novel registration algorithm
to align HARDI data characterized by ODFs across subjects under the Riemannian
manifold of ODFs and the LDDMM framework introduced in Chapter 2. Our proposed
algorithm seeks an optimal diffeomorphism of large deformation between two ODF
fields in a spatial volume domain and at the same time, locally reorients an ODF in a
23
3. DIFFEOMORPHIC METRIC MAPPING OF HIGH ANGULAR
RESOLUTION DIFFUSION IMAGING BASED ON RIEMANNIAN
STRUCTURE OF ORIENTATION DISTRIBUTION FUNCTIONS
manner such that it remains consistent with the surrounding anatomical structure. To
this end, we first define the reorientation of an ODF when an affine transformation is
applied and subsequently, define the diffeomorphic group action to be applied on the
ODF based on this reorientation. We incorporate the Riemannian metric of ODFs for
quantifying the similarity of two HARDI images into a variational problem defined
under the LDDMM framework. We finally derive the gradient of the cost function in
both Riemannian spaces of diffeomorphisms and the ODFs, and present its numerical
implementation. Both synthetic and real brain HARDI data are used to illustrate the
performance of our registration algorithm.
HARDI
data
ODF
images
ODF

, when a non-singular
affine transformation
A
is applied. As illustrated in Figure 3.2, we denote the trans-
formed
√
ODF
as

ψ(

s)=Aψ(s)
, reflecting the fact that an affine transformation
induces changes in both the magnitude of ψ and the sampling directions of s. We will
now show how to derive the analytical form when a non-singular affine transformation
acts on an ODF.
Figure 3.2:
Illustration of affine transformation on square-root ODFs. (Similar to the shape
of ODF, the colors of ODF also indices the relative values of ODF in each direction, where
blue stands for low ODF value and red for high value.)
First of all, we denote
s =(r sin θ cos ϕ, r sin θ sin ϕ, r cos θ)
in Cartesian coor-
dinates and

s =(r sin

θ cos ϕ, r sin

θ sin ϕ, r cos



s


s

d

Ω, (3.2)
and in polar coordinates,
p(θ, ϕ)sinθdθdφ =

p(

θ, ϕ)sin

θd

θdϕ. (3.3)
When

s = As in Cartesian coordinates, for the small volume fraction, we have
dxdydz =det(A)dxdydz,
and in polar coordinates,
r
2
sin

θd

det(A)sinθdθdϕ. (3.4)
From Eqns. (3.3) and (3.4), we get
p(θ, ϕ)sinθdθdφ =

p(

θ, ϕ)sin

θd

θdϕ
=

p(

θ, ϕ)
s
3
As
3
det(A)sinθdθdϕ. (3.5)
26
3.1 Affine Transformation on Square-Root ODFs
Removing sin θdθdϕ from both sides yields
p

s
s

=

A
−1
s
A
−1
s

.
For an ODF, s ∈ S
2
, s =1, and therefore, we have

p (s)=
det A
−1
A
−1
s
3
p

A
−1
s
A
−1
s

⇒ Aψ(s)=



s
are normalized back into the unit sphere
S
2
. This is analogous to a pullback deformation. Notice that for
s ∈ S
2
, Eq.
(3.7)
defines an invertible function of
s
and therefore, we can find the ODF
Aψ(s)
using the
change-of-variable technique of PDF. Recall the fundamental theorem for PDF: let
X
be a continuous random variable having probability density function
f
X
(x)
. Suppose
g(x)
is one-to-one and differentiable function of
x
. Then the random variable
Y
defined
by
Y = g(X)

. In either case, the following theorem is obtained.
Theorem 3.1.1. Reorientation of ψ based on affine transformation A.
Let
Aψ(s)
be the result of an affine transformation
A
acting on a
√
ODF ψ(s)
. The following
analytical equation holds true
Aψ(s)=

det A
−1
A
−1
s
3
ψ

A
−1
s
A
−1
s

, (3.8)
where · is the norm of a vector.


=

det A
−1
det B
−1
A
−1
s
3
B
−1
s
3
ψ

A
−1
B
−1
s
A
−1
B
−1
s

=


Aψ(s)
varies when
A
is a rotation, shearing,
or scaling and
ψ(s)
is an isotropic ODF, an ODF with a single fiber, or an ODF with
crossing fibers. From Figure 3.3, one immediately observes that a shearing or scaling
introduces anisotropy under the reorientation scheme used here. The phenomena is in
line with what is observed in [
92
]. By construction,
Aψ(s)
fulfills the definition of the
√
ODF
, i.e.,
Aψ(s)
is positive and the integration of
(Aψ(s))
2
is equal to
1
. Hence,
the similarity of
Aψ(s)
to the square-root ODFs can be quantified in the Riemannian
structure given in §2.1 for the HARDI registration.
3.2 Diffeomorphic Group Action on Square-Root ODFs
We have shown in

x ∈ Ω
. We define the action of diffeomorphisms on
ψ(s,x)
in the
29
3. DIFFEOMORPHIC METRIC MAPPING OF HIGH ANGULAR
RESOLUTION DIFFUSION IMAGING BASED ON RIEMANNIAN
STRUCTURE OF ORIENTATION DISTRIBUTION FUNCTIONS
Figure 3.3:
Examples of local affine transformations on an isotropic ODF in the first
row, an ODF with a single orientation fiber in the middle row, and an ODF with crossing
fibers in the bottom row. From left to right, three types of affine transformations,
A
,on
the ODFs are demonstrated: in panel (a), a rotation with angle
θ
z
, where
A =[cosθ
z
−
sin θ
z
0; sin θ
z
cos θ
z
0; 001]
; in panel (b), a vertical shearing with factor
ρ

, i.e.,
A
x
= D
x
φ
. According to Eq.
(3.8)
, the action
of diffeomorphisms on ψ(s,x) can be computed as
φ · ψ(s,x)=





det

D
φ
−1
φ

−1




D
φ

Property 2. (The law of composition for diffeomorphic group action on ODF)
As-
sume
φ
and
ϕ
to be two diffeomorphisms and
ψ(s,x)
as the
√
ODF
with the orientation
30
3.2 Diffeomorphic Group Action on Square-Root ODFs
direction s ∈ S
2
located at x ∈ Ω. The following property holds true
ϕ · (φ · ψ)(s,x)=(ϕ ◦φ) · ψ(s,x) , (3.11)
where ◦ stands for composition between φ and ϕ.
Proof. Based on the equation (3.10), it yields
ϕ · (φ · ψ)(s,x)=ϕ ·

(A
x
ψ(s,x)) ◦ φ
−1
(x)

=


x
ϕ, and B
φ(x)
= D
φ(x)
ϕ.
Using the property from the equation (3.9), we have
ϕ · (φ · ψ)(s,x)=

B
φ(x)
A
x

ψ(s,x)

◦ φ
−1
◦ ϕ
−1
(x)
=[D
x
(ϕ ◦ φ) ψ(s,x)] ◦ (ϕ ◦φ)
−1
(x)
=(ϕ ◦ φ) · ψ(s,x) .
For the sake of simplicity, we denote φ · ψ(s,x) as
φ · ψ(s,x)=Aψ ◦ φ
−1

t
, which are solutions of ordinary differential equations
˙
φ
t
= v
t
(φ
t
),t∈ [0, 1],
starting
from the identity map
φ
0
= Id
. They are therefore characterized by time-dependent
velocity vector fields
v
t
,t∈ [0, 1]
. We define a metric distance between a target volume
ψ
targ
and a template volume
ψ
temp
as the minimal length of curves
φ
t
·ψ

·
V
.
To ensure solutions are diffeomorphisms,
V
must be a space of smooth vector fields
[
88
]. Using the duality isometry in Hilbert spaces, one can equivalently express the
lengths in terms of m
t
, interpreted as momentum such that for each u ∈ V ,
m
t
,u◦ φ
t

2
= k
−1
V
v
t
,u
2
, (3.13)
where we let
m, u
2
denote the

m
t
. Using the identity
v
t

2
V
= k
−1
V
v
t
,v
t

2
= m
t
,k
V
m
t

2
and the standard fact that energy-minimizing
curves coincide with constant-speed length-minimizing curves, one can obtain the
metric distance between the template and target
√
ODF

m
t
:
˙
φ
t
=k
V
m
t
(φ
t
),
φ
0
=Id
ρ(ψ
temp
, ψ
targ
)
2
+ λ

x∈Ω
E
x
(φ
1
· ψ

J(m
t
) = inf
m
t
:
˙
φ
t
=k
V
m
t
(φ
t
),
φ
0
=Id

1
0
m
t
,k
V
m
t

2

targ
(s,x)
as
ψ
targ
(x)
. Note that since we are dealing with vector fields in
R
3
, the kernel of
V
is
a matrix kernel operator in order to get a proper definition. We define this kernel as
k
V
Id
3×3
, where
Id
3×3
is an identity matrix, such that
k
V
can be a scalar kernel. In the
rest of the chapter, we shall refer to this LDDMM mapping problem as LDDMM-ODF.
3.3.1 Gradient of J with respect to m
t
The gradient of
J
with respect to

t
)=2m
t
+ λη
t
, (3.16)
where
η
t
= ∇
φ
1
E +

1
t

∂
φ
s
(k
V
m
s
)


(η
s
+ m

E
, where
E =

x∈Ω
E
x
dx
, which will be
discussed in the following.
34
3.3 Large Deformation Diffeomorphic Metric Mapping for ODFs
Gradient of E with respect to φ
1
:
The computation of
∇
φ
1
E
is not straightforward
and the Riemannian structure of ODFs has to be incorporated. Let’s first compute
∇
φ
1
E
x
at a fixed location,
x
. We consider a variation

directly give the expression of
∂

E
x
|
=0
and the reader is referred to
§
3.3.1.1 for the full
derivation of terms (A) and (B) in the following equation.
∂

E
x
|
=0
(3.18)
=2

log
Aψ
temp
◦φ
−1
1
(x)
ψ
targ
(x),

−1
1
(x)
ψ
targ
(x),
∂ log
Aψ
temp
◦φ
−1
1
(x)
A

ψ
temp
◦ (φ

1
)
−1
(x)
∂
|
=0

Aψ
temp
◦φ

∂
|
=0

Aψ
temp
◦φ
−1
1
(x)
  
term (A)
−2

log
Aψ
temp
◦φ
−1
1
(x)
ψ
targ
(x),
∂ log
Aψ
temp
◦φ
−1
1

(x)
ψ
targ
(x),


(D
x
φ
1
)
−
∇
x
(Aψ
temp
),h


◦ φ
−1
1
(x)

Aψ
temp
◦φ
−1
1
(x)

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
2
x

Aψ
temp
(x)

div

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
3
x


is a
3 × 1
vector with the
i
th element as
one and the rest as zero.
·, ·
Aψ
temp
(x)
is the Fisher-Rao metric defined in Eq.
(2.2)
.
∇
x
(Aψ
temp
)
in term (A) is the first derivative of the
√
ODF
,
Aψ
temp
, with respect to
x
. Since
Aψ
temp
also lies in the Riemannian manifold of

(x)
Aψ
temp
(x + e
1
)
1
|e
2
|
log
Aψ
temp
(x)
Aψ
temp
(x + e
2
)
1
|e
3
|
log
Aψ
temp
(x)
Aψ
temp
(x + e

L
i
x
=(D
x
φ
1
)
−1
su
i
−
1
2
Aψ
temp
(s,x)w
i
,
where
w
i
is the
i
th column of
(D
x
φ
1
)

φ
1
)
−
∇
s

ψ

s
s
,x


s
3

.
In sum,
∇
φ
1
E
can be computed by integrating
∇
φ
1
E
x
over the image space. With a

x
φ
1
)
−
(3.19)
⎡
⎢
⎢
⎢
⎢
⎣
 log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)) ,
1
|e
1
|
log
Aψ
temp
(x)
Aψ

temp
(x)
 log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)) ,
1
|e
3
|
log
Aψ
temp
(x)
Aψ
temp
(x + e
3
) 
Aψ
temp
(x)
⎤
⎥
⎥

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
2
x

Aψ
temp
(x)

div

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
3
x


φ

1
and thus, do not consider the variation in
A
, i.e.,
A

is ignored. Therefore, in [
2
], the gradient of
E
with respect to
φ
1
only incorporates
term (A) of Eq.
(3.18)
. This term is similar to the scalar image matching case and only
takes into account image shape difference in the volume space. We illustrate this in
Figure 3.5, where we have one template image and two target images. Figure 3.5 (a)
shows the template image, where its overall image shape is circular and the ODFs at
each voxel inside the circle are oriented horizontally. Figure 3.5 (b) shows the first target
image, where its overall image shape is an ellipsoid and the ODFs inside its voxels are
oriented horizontally. Figure 3.5 (c) shows the second target image, where its overall
image shape is circular as the template image but the ODFs at each voxel inside the
circle are oriented at
45
◦
. The results obtained using only term (A) as proposed in [

(3.18)
. The reader
can skip this subsection without any discontinuation by assuming that the derivation of
terms (A) and (B) in Eq. (3.18) holds true.
Term (A):
For the sake of simplicity, we denote term (A) of Eq.
(3.18)
as
E
A
and
rewrite
E
A
= −2

log
Aψ
temp
◦φ
−1
1
(x)
ψ
targ
(x),
∂ log
Aψ
temp
◦φ

=0
= −[(Dφ
1
)
−1
h] ◦ φ
−1
1
(x),wehave
E
A
=2

log
Aψ
temp
◦φ
−1
1
(x)
ψ
targ
(x),


(D
x
φ
1
)

Aψ
temp
◦φ
−1
1
(x)
ψ
targ
(x),
∂ log
Aψ
temp
◦φ
−1
1
(x)
A

ψ
temp
◦ (φ
1
)
−1
(x)
∂
|
=0

Aψ

temp
◦φ
−1
1
(x)
.
According to Theorem 3.1.1, we have
A

ψ
temp
(s,x) ◦φ
−1
1
(x)=
⎡
⎢
⎣





det

D
x
φ

1

φ

1

−1
s
,x

⎤
⎥
⎦
◦ φ
−1
1
(x) .
Denote s =

D
x
φ
1

−1
s.
Given
∂(D
x
φ

1

1
)
−1
,we
can now compute
∂A

ψ
temp
(s,x) ◦φ
−1
1
(x)
∂
|
=0
=

−

det

D
x
φ
1

−1

∇

2
Aψ
temp
(s,x)trace

D
x
h(D
x
φ
1
)
−1


◦ φ
−1
1
(x)
=


D
x
h(D
x
φ
1
)
−1

φ
1

−1
(D
x
φ
1
)
−
∇
s

ψ

s
s
,x


s
3

.
We now derive the above equation in order to express it in an explicit form of
h
. Before
doing so, we first define a
3 × 3
identity matrix as

3
]
, where
w
i
is the
i
th column of
(D
x
φ
1
)
−1
. Thus, the trace of
D
x
h(D
x
φ
1
)
−1
39
3. DIFFEOMORPHIC METRIC MAPPING OF HIGH ANGULAR
RESOLUTION DIFFUSION IMAGING BASED ON RIEMANNIAN
STRUCTURE OF ORIENTATION DISTRIBUTION FUNCTIONS
can be written as
trace


(x)
∂
|
=0
=


D
x
h(D
x
φ
1
)
−1
s,u

−
3

i=1
1
2
Aψ
temp
(s,x)

D
x
hw

⎢
⎢
⎢
⎢
⎣
div(u
1
w)
div(u
2
w)
div(u
3
w)
⎤
⎥
⎥
⎥
⎥
⎦
,h

,
where u
i
is the ith element of u.
As a consequence, when defining
L
i
x

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
1
x

Aψ
temp
(x)

div

log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)),L
2
x


◦ φ
−1
1
(x) .
40
3.3 Large Deformation Diffeomorphic Metric Mapping for ODFs
3.3.2 Euler-Lagrange Equation for LDDMM-ODF
To summarize, LDDMM-ODF can be written in the form of the variational problem as
Eq. (3.15) with its Euler-Lagrange optimality conditions as the following theorem.
Theorem 3.3.1.
Given the energy
J
in the following form of Eq.
(3.15)
, the gradient
of
J
with respect to
m
t
given by the following equation has fixed points satisfying the
necessary minimizer of the objective function J.
∇J(m
t
)=2m
t
+ λη
t
, (3.20)
where

V
m
s
)
is the partial derivative of
k
V
m
s
with respect to
φ
s
.
η
t
in Eq.
(3.17)
can be solved backward given
η
1
= ∇
φ
1
E
, where
∇
φ
1
E
is the gradient given in Eq.

m
t
to be the sum of Dirac measures such that
m
t
=

N
i=1
α
i
(t) ⊗δ
φ
t
(x
i
)
and
41
3. DIFFEOMORPHIC METRIC MAPPING OF HIGH ANGULAR
RESOLUTION DIFFUSION IMAGING BASED ON RIEMANNIAN
STRUCTURE OF ORIENTATION DISTRIBUTION FUNCTIONS
therefore:
ρ(ψ
temp
, ψ
targ
)
2
=


,
where
α
i
(t)
is the momentum vector at
x
i
and time
t
. In discretization of the spherical
domain
S
2
, we discretize it into
N
S
equally distributed diffusion sampling directions on
the sphere. For each diffusion sampling direction
k
, it can be represented as
3
D vector
with unit length
s
k
in Cartesian coordinates and
(r
k

(x
i
),φ
t
(x
j
))α
j
(t) . (3.22)
2. Compute ∇
φ
1
(x
i
)
E in Eq. (3.19), which is described in details below.
3.
Solve
η
t
=[η
i
(t)]
N
i=1
in Eq.
(3.17)
using the backward Euler integration, where
i
indices x

3.3 Large Deformation Diffeomorphic Metric Mapping for ODFs
We now discuss how to compute
∇
φ
1
(x
i
)
E
in Eq.
(3.19)
, which involves the
√
ODF
interpolation in the spherical coordinate for
Aψ
temp
(x
i
)
atafixed
x
i
and the
√
ODF
interpolation in the image spatial domain for
ψ
targ
(φ

1
(x
i
))
. For the
√
ODF
interpolation in the spherical coordinate for
Aψ
temp
(x
i
)
at a fixed
x
i
, we compute
Aψ
temp
(s
k
,x
i
)
according to Eq.
(3.8)
using angular interpolation on
S
2
based on

(x
i
)

=exp
ψ
targ

s
k
,φ
1
(x
i
)


j∈N
i
w
j
log
ψ
targ

s
k
,φ
1
(x

)
and
x
j
. The exponential maps and logarithm maps can be computed
via Eq. (2.4) and Eq. (2.5) respectively. Finally, the inner product in Eq. (3.19),
 log
Aψ
temp
(x)
ψ
targ
(φ
1
(x)) ,
1
|e
i
|
log
Aψ
temp
(x)
Aψ
temp
(x + e
i
) 
Aψ
temp

differ from the template image (Figure 3.6 (a)) to the target image (Figure 3.6 (b)).
We will compare the performance of LDDMM-ODF to the LDDMM algorithm based
on DTI (LDDMM-DTI) [
96
]. We refer the reader to [
96
] for detailed mathematical
derivation for LDDMM-DTI.
In the HARDI model, such orientation differences are encoded by the ODFs, while
in the DTI model, the diffusion tensors of both the template and target data look
like disks, where the first two eigenvalues being equal and the third eigenvalue being
almost zero. Although the overall image shapes are the same in both the template and
target HARDI data, the LDDMM-ODF algorithm is able to characterize the orientation
difference of the ODFs between them by generating the deformation shown in Figure
3.6 (d) with the help of term (B) of Eq.
(3.18)
. Note that a mask is used here and only
the ODFs in the circle is transformed. LDDMM-DTI fails to find any deformation
(Figure 3.6 (e)) even though the reorientation of the tensor is taken into account in the
tensor mapping.
44
3.5 HARDI Data of Children Brains
Figure 3.5:
The first and second rows respectively illustrate the original HARDI and their
enlarged images. Compared to the image on panel (a), the image on panel (b) has the
same ODFs but a different ellipsoidal image shape, while the image on panel (c) shows
different ODFs but the same circular image shape. Panels (d) and (e) show the deformations
(grid) and the corresponding momenta (arrows), calculated using
∇
φ

young children (
6
years old). All
three algorithms are developed under the LDDMM framework as given in
§
3.3 with
the exception that the matching functional,
E
, is the least square difference between
two image intensities for the image mapping, LDDMM-FA, and the Frobenius norm
between two tensors for the DTI mapping, LDDMM-DTI. More precisely, LDDMM-FA
is based on the method developed by [
97
] and LDDMM-DTI is based on the method
developed by [96]. In our implementation however, we optimize the deformation with
respect to the momentum rather than the velocity (see [
93
]). It is important to note that
all three mapping algorithms used in the following evaluation have the same numerical
scheme, such that any potential errors due to numerical related issues are avoided and
we can make a fair comparison.
Our image data are acquired using a
3T
Siemens Magnetom Trio Tim scanner with a
32
-channel head coil at the National University of Singapore. Diffusion weighted imag-
ing protocol is a single-shot echo-planar sequence with
55
slices of
2.3mm