Construction of a 4D Brain Atlas and Growth Model Using Diffeomorphic Registration - Longitudinal Registration and Shape Modeling - page 27

where

T ij =exp( ˄ v ij )

(4)

The non-negative constants ʻ 1 , 2 , 3 weigh the contribution of each term. Bend-
ing energy, BE, and a Jacobian-based penalty term, JAC, are used to prevent
folding and tearing, i.e.,

∂ 2 v ( x )
∂x m ∂x n

2

3 3

1
|ʩ|

BE ( v )= (5)

x ∈ʩ m =1 n =1

1
|ʩ|

log 2 (det ( J v ( x ))) JAC ( v )= (6)

x ∈ʩ

where ʩ denotes the finite set of positions on the transformed image lattice at
which the energy is evaluated, and J v ( x ) denotes the Jacobian matrix of the
velocity field v evaluated at x . This is similar to the approach of Modat et
al. [ 13 ], as the exponential map is only guaranteed to generate a diffeomorphism
when the velocity field is suciently smooth.
Our formulation differs from others by the use of a single parametric trans-
formation and only one similarity evaluation as opposed to separate forward and
backward transformations [ 3 , 12 ] or similarity evaluated twice [ 20 ]. The method
is similar to that of [ 11 ] in that we transform both images half-way and use
a single similarity term. The image similarity is therefore evaluated for images
which are equally affected by the deformation and interpolation. Additionally,
inverse consistency reduces the number of required pairwise registrations.
We use an approximate but fast scaling-and-squaring on the control point
lattice as presented in [ 13 ] for the computation of the exponential map ( 4 ). Given
the derivative of ( 3 ) with respect to ν ij , we perform a conjugate gradient descent

to find the set of parameters which minimize E .
The NMI gradient is first computed separately for each half transformation as
in [ 14 ]. The resulting gradient fields are then added up with their corresponding
weights ˄ =

0 . 5 to obtain the gradient field, ʴ u ij . Note that the scaling factor
˄ accounts for both the averaging of the separate NMI gradient fields as well as
the inversion of the gradient corresponding to the half backward transformation.
The obtained gradient field is then composed with the current velocity field.
This composition is approximated in the log-domain using the Baker-Campbell-
Hausdorff (BCH) formula [ 20 ], i.e.,

±

ʴ v ij = ʴ u ij + 1

2 [ v ij ,ʴ u ij ]+ 1 12 [ v ij , [ v ij ,ʴ u ij ]]

(7)

where the first term of the BCH formula is omitted in order to obtain the dif-
ference between the two velocity field estimates. This computation is similar to
the update step of the symmetric LogDemons [ 20 ]. By interpolating all vector
fields (incl. the Lie bracket [

]) by cubic B-spline functions with control points
defined on the same lattice as the local velocity field ( 2 ), the NMI gradient with
respect to

·, ·

ν ij

is approximated by the spline coecients of ( 7 ).