IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 12, NO. 11, NOVEMBER 2003 1427
Fast Parametric Elastic Image Registration
Jan Kybic, Member, IEEE, and Michael Unser, Fellow, IEEE
Abstract—We present an algorithm for fast elastic multidimen-
sional intensity-based image registration with a parametric model
of the deformation. It is fully automatic in its default mode of op-
eration. In the case of hard real-world problems, it is capable of
accepting expert hints in the form of soft landmark constraints.
Much fewer landmarks are needed and the results are far superior
compared to pure landmark registration. Particular attention has
been paidtothe factorsinfluencing the speedof this algorithm.The
B-spline deformation model is shown to be computationally more
efficient than other alternatives.
The algorithm has been successfully used for several two-di-
mensional (2-D) and three-dimensional (3-D) registration tasks in
the medical domain, involving MRI, SPECT, CT, and ultrasound
image modalities. We also present experiments in a controlled
environment, permitting an exact evaluation of the registration
accuracy. Test deformations are generated automatically using a
random hierarchical fractional wavelet-based generator.
Index Terms—Elastic registration, image registration, land-
marks, splines.
I. INTRODUCTION
I
MAGE registration is the task of finding a correspondence
function mapping coordinates from a reference image to
coordinates of homologous points in a test image [1]. We call
the registration elastic [2] if the family of correspondence func-
tions is sufficiently general, capable of expressing essentially
arbitrary nonlinear relations.
1
Image registration is applied
in the areas of motion analysis [4]–[6], video compression
and coding [7], object tracking [8], or image stabilization. It
leads to algorithms for segmentation [9], depth reconstruction
from stereo images [10], [11], and for general 3-D recon-
struction. In the biomedical domain, there is a frequent need
for comparing images for analysis and diagnostic purposes.
This is accomplished by registering the images and aligning
them by warping using the correspondence function identified.
Applications include intra-subject [12], inter-subject [13], [14],
and inter-modality analysis [15], [16], [17], registration with
annotated atlases [18], [19], quantification and qualification of
feature shapes and sizes [20], distortion compensation [21],
[22] and motion detection [23], [24] and compensation [25].
Manuscript received September 6, 2001; revised February 12, 2003. The as-
sociate editor coordinating the review of this manuscript and approving it for
publication was Dr. Patrick Perez.
J. Kybic was with the Biomedical Imaging Group, LIB, Swiss Federal Insti-
tute of TechnologyLausanne,CH-1015LausanneEPFL,Switzerland.Heisnow
with the Center for Machine Perception, Department of Cybernetics, Faculty of
Electrical Engineering, Czech Technical University, Prague 6, Czech Republic
(e-mail: kybic@ieee.org).
M. Unser was with the Biomedical Imaging Group, LIB, Swiss Federal Insti-
tute of Technology Lausanne, CH-1015 Lausanne EPFL, Switzerland.
Digital Object Identifier 10.1109/TIP.2003.813139
1
Note that elasticity is used here in a wider sense than just the mechanical
linear elasticity [3].
Various nonlinear registration algorithms for brain warping
applications are presented by Warfield et al. [9]. Bayesian inter-
pretation of elastic matching was reviewed by Gee [19], also in
the context of human neuroanatomy. Articles by Van den Elsen
et al. [26] and Maintz and Viergever [27] contain a very com-
prehensive and detailed classification of available methods for
medical imaging applications. Lester and Arridge [28] treat the
hierarchical aspects of the algorithms.
The deformation models of elastic registration algorithms
fall into two basic categories. The first type are nonparametric,
local methods—the deformation function is basically uncon-
strained and belongs to a very large and unrestrictive functional
space. These methods can be formulated as variational, defining
a scalar criterion that completely determines the final solution
[2]. More generally, they can be also expressed using partial
differential equations (PDE) [29]–[32].
The presented algorithm belongs to a second group of
methods that use parametric models, representing the deforma-
tion by a moderate number of parameters, often in the multi-
scale setting. Specific examples include hierarchical basis func-
tions by Moulin et al. [7], quadtree-splines [5], multiresolution
subspaces [33], [34], and wavelets [35], [36]. Splines are well
suited for this kind of problems; they have appeared in various
incarnations. In this paper we use a multiresolution B-spline
representation, as was initially suggested in the pioneering
work of Szeliski et al. [10], [5].
A. Proposed Algorithm
The algorithm described in this article is a synthesis of sev-
eral ideas. First, it is a generalization to multiple dimensions of
the unidirectional registration algorithm we described in [22].
Its main features are the use of B-splines to describe both the
image and the deformation, a double multiresolution strategy
(for both the image and the deformation), a scalar pixel-based
difference measure, and an iterative multidimensional opti-
mization algorithm [37], [38]. The deformation model has been
generalized and the whole algorithm re-engineered for faster
execution.
Second, we present the idea of semi-automatic registration,
targeted to more difficult registration problems. We ask an
expert to identify a small number of corresponding points in
both images. The points are also called landmarks [3], [12],
[39], [40]. We add a term to the data part of the criterion, to
steer the algorithm toward the correct solution indicated by the
landmarks.
B. Organization of This Article
In Section II, we describe the concept of registration by min-
imization, the difference measure, the B-spline image model,
1057-7149/03$17.00 © 2003 IEEE
1428 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 12, NO. 11, NOVEMBER 2003
and the structure of the deformation model. In Section III we
justify our choice of B-splines as basis functions for the de-
formation model. We present the optimization method in Sec-
tion IV, where we also describe the multiresolution strategy.
Section V is devoted to the semi-automatic mode incorporating
landmark information into the global criterion. We deal with
implementation issues in Section VI and present experiments
and applications in Section VII. For more details, we refer the
reader to the first author’s thesis report [38] and its associated
web page.
II. P
ROBLEM FORMULATION
The input images are given as two -dimensional discrete
signals
and , where , and is an -
dimensional discrete interval representing the set of all pixel
coordinates in the image. We call
and reference and test
images, respectively. We suppose that the test image is a geo-
metrically deformed version of the reference image, and vice
versa.
2
This is to say that the points with the same coordinate
in the reference image and in the warped test image
should correspond. Here, is a continuous
version of the test image and
is a deformation (correspon-
dence) function to be identified
A. Cost Function
The two images
, will not be identical because of noise
and also because the assumption that there is a geometrical map-
ping between the two images is not necessarily correct. There-
fore, we define the solution to our registration problem as the
result of the minimization
, where is
the space of all admissible deformation functions
.Wehave
chosen the SSD (sum of squared differences) criterion
(1)
because it is fast to evaluate and yields a smooth criterion sur-
face which lends itself well to optimization. Minimization of (1)
yields the optimal solution
in the ML (maximum likelihood)
sense under the assumption that
is a deformed (warped) ver-
sion of
with i.i.d. (independent and identically distributed)
Gaussian noise added to each pixel. The SSD criterion proved
to be robust enough, especially if preprocessing was used to
equalize the image values—we mostly applied high-pass fil-
tering and histogram normalization [22]. In principle, there is no
difficulty in extending our method for more sophisticated pixel-
based similarity measures, such as information-based measures
[41], especially mutual information [17], or weighted
norms.
Only the evaluation of the criterion and its derivatives (gradient)
needs to be changed.
2
In the multimodal case, which we are not considering here, there can be also
an intensity mapping between the two images.
B. Image Interpolation
In accordance with [22], we choose to interpolate the image
using uniform B-splines:
3
(2)
where
is a tensor product of B-splines of degree , that is
, with .
C. Deformation Model Structure
So far, we have considered the deformation function
to be
an arbitrary admissible function
. We will restrict
it now to a family of functions described by a finite number of
parameters
:
(3)
where
is a set of parameter indexes and are the corre-
sponding basis functions. This transforms a variational problem
into a much easier finite-dimensional minimization problem, for
which numerous algorithms exist [43]. Moreover, the restriction
of the family
of all possible functions can already guarantee
some useful properties, such as the regularity (smoothness) of
the solution. Note that the addition of
in the above equation
makes the set of zero parameters correspond to identity.
D. Existence, Unicity, and Regularization
Note that the criterion
is nonnegative and continuous and
is periodic due to boundary conditions. Consequently, has
a minimum; i.e., the proposed problem has a solution. However,
depending on the images at hand, the solution does not have
to be unique and there can be local minima. Fortunately, this
does not pose problems in practice thanks to a multiresolution
approach (Section IV-B) which smoothes out images at coarse
levels and brings us sufficientlyclose to the solution at fine reso-
lution levels. The algorithm will find a solution if started within
the attraction basin of that solution. The virtual springs (Sec-
tion V) play a role of an a priori information and a regulariza-
tion term; extra regularization can be applied [44] if desired.
III. D
EFORMATION BASIS
The purposeofthis section is to motivate our choice of (cubic)
B-splines [42] as the most adequate basis functions
to rep-
resent the deformation in model (3). The alternative possibil-
ities that come to mind are polynomials [45], harmonic func-
tions [18], [46], radial basis functions [3], [47], and wavelets
[35], [48], [49].
It is highly desirable to have as few basis functions as pos-
sible to contribute to each particular point, while keeping the
3
Uniform symmetric B-splines [42] of degree
n
are piecewise polynomials
of degree
n
. The polynomial pieces are delimited by uniformly placed knots.
B-splines of degree
n
have continuous derivatives up to order
n
0
1
everywhere.
Their integer shifts form a basis. The first (degree zero) symmetric B-spline is
defined as
(
x
)=1
for
x
2
(
0
1
=
2
;
1
=
2)
and 0 otherwise. Higher order
B-splines are defined by recursively as
=
3
; and their support is
(
0
(
n
+1)
=
2
;
+(
n
+1)
=
2)
.
KYBIC AND UNSER: FAST PARAMETRIC ELASTIC IMAGE REGISTRATION 1429
Fig. 1. Basis functions involved in evaluatingthevalue of a 1-D function atone
point (denoted by a vertical line): (a) radial basis functions
j
x
j
, (b) harmonic
functions, (c) cubic B-splines, and (d) cubic B-spline wavelets.
approximation quality. First, short basis functions have small
overlap. This reduces the interdependency between the coeffi-
cients (parameters) and consequently makes the minimization
problem easier to solve. Small overlap also makes the Hessian
(the matrix of second partial derivatives, needed for some opti-
mizers) more sparse and therefore potentially faster to invert.
Second, the size of the support of the basis functions directly
influences the speed of the calculation. The evaluation of the
deformation function (3) at
points costs op-
erations, where
is the number of functions contributing
to a single point.
4
The cost of evaluating the gradient of the
criterion
with respect to the coefficients is higher but asymp-
totically equivalent, because each of the
pixels contributes
to exactly
components of the gradient. Note that this cost
is independent of the total number
of the basis functions
(unless
). The cost of evaluating the Hessian is
operations. (See also Section VI-A.)
Fig. 1 shows the generating functions needed to calculate a
valueat one point (denoted by the verticalbar) for variousbases;
only functions that are nonzero at that point are considered.
Except for the Fourier basis, we choose basis functions of the
4
We assume that the cost of evaluating the basis function itself is constant or
that their values can be precalculated.
same degree (cubic), generating the same space. We see clearly
that the least number of contributing functions (four) is in the
B-spline case. This effect turns out to be even more dramatic in
higher dimensions.
The reasoning above rules out the polynomials because no
fast algorithm is known for their evaluation and the brute-force
evaluation is slow due to their long support. As for the radial
basis functions, although there are algorithms with reduced
asymptotical complexity for evaluation of radial basis functions
[50]–[53], their overhead is still nonnegligible. We decided
against the harmonic (Fourier) basis functions because of their
lack of localization (the fact that any two of them overlap).
Another argument against the Fourier basis is that it cannot
express linear functions (affine deformations). The only two
remaining candidate basis are therefore B-splines and B-spline
wavelets.
A. Splines Versus Wavelets
To make a fair comparison between B-spline and wavelet
bases, we consider compactly supported cubic B-spline
wavelets [54] spanning the same cubic spline space. First, let
us analyze the task of evaluating the deformation at a single
point. For simplicity, we will work in 1-D. There are only
four participating B-splines altogether while there are four
participating B-spline wavelets at each level, plus four scaling
functions (cubic B-splines) at the coarsest level. Second, to
evaluate the deformation at a set of equally spaced points
(this corresponds to a regular grid in multiple dimensions),
the direct B-spline representation is also the most efficient,
the interpolation requiring only four multiplications per pixel.
This is better than all alternatives available when using the
B-spline wavelets, including iterative filterbank and FFT-based
algorithms.
Note that the complexity of evaluation of the gradient of the
criterion corresponds to the complexity of the evaluation of the
deformation because the same type of formula is involved (see
Section VI-A).
B. B-Spline Deformation Model
The B-spline deformation model is obtained by substituting a
scaled version of the B-spline (or tensor product thereof) in (3)
(4)
where
is the degree of splines used, is the knot spacing,
and the division is taken elementwise. This corresponds to
placing the knots on a regular grid over the image. We require
the node spacing
to be integer, which together with the
separability of
implies that the values of the B-spline
are only needed at a very small number of points
and can be precalculated. We can evaluate on the
whole grid with the cost of only
multiplications
per pixel.
The B-spline model has good approximation properties and
is fast to evaluate. It is physically plausible, for example cubic
splines minimize the ‘strain energy’
[55], [56]. It can
encode all affine transformations, including rigid body motion.
1430 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 12, NO. 11, NOVEMBER 2003
Moreover, B-splines are scalable in the sense that any coarse
level deformation can be represented at a finer scale without
any loss of information given an integer ratio between scales.
The expansion operator (Section VI-C) is therefore exact.
IV. O
PTIMIZATION STRATEGY
A. Optimization Algorithm
Recall from (1) and (4) that we need to minimize a crite-
rion
with respect to a finite number of parameters . To de-
termine which of the many available algorithms performs best
in our context, we tested four local iterative algorithms which
can be cast into a common framework: At each step
we take
the actual estimate
and calculate a proposed update .
If the step is successful, then the proposed point is accepted,
. Otherwise, a more conservative update
is calculated, and the test is repeated.
1) Gradient descent with feedback step size adjustment with
update rule:
. After a successful
step,
is multiplied by , otherwise it is divided by
.
5
2) Gradient descent with quadratic step size estimation.We
choose a step size
minimizing the following approxi-
mation of the criterion around
:
, where is identified from
the two last calculated criterion values
. As a fallback
strategy, the previous step size is divided by
, as above.
3) Conjugated gradient. This algorithm [43] chooses its de-
scent directions to be mutually conjugate so that moving
along one does not spoil the result of previous optimiza-
tions. To work well, the step size
has to be chosen opti-
mally. Therefore, at each step, we need to run another in-
ternal one-dimensional minimization routine which finds
the optimal
; this makes it the slowest algorithm in our
setting.
4) Marquardt–Levenberg. The most effective algorithm in
the sense of the number of iterations was a regularized
Newton method inspired by the Marquardt–Levenberg
algorithm (ML), as in [22]. Various approximations of
the Hessian matrix
were examined (see also Sec-
tion VI-A).
As the behavior of all optimizers is comparable at the beginning
of the optimization process (see Fig. 2), the main factor deter-
mining the speed is the cost of a single iteration. The evaluation
costs are presented in Table I; for the ML algorithm, the cost of
the Hessian matrix inversion (which grows with the cube of the
number of parameters) must be added. It follows that the gra-
dient descent (GD) iterations are the least costly, the difference
between the two variant being minimal. We therefore recom-
mend to use the GD algorithm with the quadratic step size esti-
mation (which works better than the feedback adjustment) and
we use it for experiments in the remainder of the paper. One ad-
ditional pleasant property of the GD algorithm is its tendency
to leave uninfluential coefficients intact, unlike the ML algo-
rithm. Consequently, less regularization is needed for the GD
algorithm.
5
We used
=
10
and
=15
.
Fig. 2. Evolution of the SSD criterion during first 18 iterations when regis-
tering the Lena image, artificially deformed with 2
2
4
2
4 cubic B-spline
coefficients and a maximum displacement of about 30 pixels, without multi-
resolution. The optimizers used were: Marquardt–Levenberg with full Hessian
(MLH), Marquardt–Levenberg with only the diagonal of the Hessian taken into
account (MLdH), and gradient descent (GD). The deformation was recovered
in all cases with an accuracy between 0.1 and 0.01 pixels (see also Section VII).
TABLE I
R
ELATIVE TIMES TO EVALUATE THE CRITERION
E
,ITS GRADIENT
r
E
, AND
HESSIAN
r
E
, FOR A VOLUME OF 64
2
64
2
17 VOXELS APPROXIMATED BY
CUBIC SPLINES, AS A
FUNCTION OF THE SPLINE DEGREE
n
USED TO MODEL
THE
DEFORMATION AND THE SIZE OF THE PARAMETER GRID
n
.(THE
ABSOLUTE TIME TO EVALUATE
E
WAS ABOUT 1s)
Fig. 3. Comparison of gradient descent (GD), conjugated gradient (CG),
and Marquardt–Levenberg (ML) optimization algorithm performances when
registering SPECT images with control grid of 6
2
6
2
6 knots. The graphs
give the value of the finest-level SSD criterion of all successful (i.e., criterion-
decreasing) iterations as a function of the execution time. The abrupt changes
are caused by transitions between resolution levels.
Under different constraints, when a small number of param-
eters is sought, the criterion is smooth, and high precision is
needed, the ML algorithm performs the best. This is because its
higher cost per iteration is compensated for by a smaller number
of iterations due to the quadratic convergence. An example of
such a situation is shown in Fig. 3. (See also [57].) Among Mar-
quardt–Levenberg (ML) algorithms, we found the performance
to be superior when using the full Hessian.
KYBIC AND UNSER: FAST PARAMETRIC ELASTIC IMAGE REGISTRATION 1431
B. Multiresolution
As in [22], we use the multiresolution approach for both
the image and deformation models. We start with the coarsest
resolution versions of both, and alternatively refine the image
and the deformation model every time convergence is reached,
until the finest level. The coarse versions of images are gener-
ated using a reduction operator (see Section VI-C). Conversely,
coarse level solutions are extrapolated to finer levels using an
expansion operator (cubic spline interpolation).
V. S
EMI-AUTOMATIC REGISTRATION
We realize that although the multiresolution approach leads
to a very robust registration algorithm, there are cases when it
is mislead by an apparent similarity of features which do not
correspond physically. Therefore, we developed an extension of
the algorithm which can use expert hints. The hints come in the
form of a set of landmarks and are used to gear the algorithm
toward the correct solution.
The idea of a hybrid registration algorithm combining both
local features (points or lines) with global ones (intensities) has
appeared for example in [58], [59]
6
and others. However, as
both [58], [59] use essentially local, nonparametric deformation
models, the landmark constraints need to be first interpolated
everywhere to serve as an a priori deformation field. This is in
contrast with our method which only imposes the landmark in-
formation at landmark points where it is really known. Another
difference is that thanks to our parametric deformation model
the additional overhead is negligible.
The landmark information is incorporated in the automatic
process using the concept of virtual springs, tying each pair
of corresponding points together. We augment the data part of
the criterion
with a term , corresponding to the poten-
tial energy of the springs, and minimize the sum of the two:
. The spring term is
(5)
where
is the number of springs, are weighting factors cor-
responding to their stiffnesses, and
, resp. , are the land-
mark positions in the reference, resp. test images. The spring
factors
control the influence of the particular landmark pairs.
We propose to start with all
and adjust them experi-
mentally to get the most satisfactory results. We should aim for
a compromise between
too small that does not succeed in
making the algorithm to converge to the right solution, and
too high that forces the solution to a landmark position that is
perhaps not sufficiently precise.
As an example, we tried to register an MRI slice from an
atlas
7
with a sample MRI test image.
8
The atlas is a labeled and
annotated collection of images. To identify the same structures
in the test image, we register it with the unlabeled version of
the atlas. Once the geometric correspondence is established, the
6
We thank the reviewers for bringing this to our attention.
7
Courtesy of Harvard Medical School, http://www.med.harvard.edu/
AANLIB/home.html.
8
We use a proton density MR image from the Visible Human project http://
www.meddean.luc.edu/lumen/meded/grossanatomy/cross_section/index.html.
structures and their labels from the atlas can be projected onto
the test image. Prior to registration, the histogram of the test
image was matched to that of the reference. The unsupervised
registration correctly registers some of the structures but misses
others; in particular the skull boundary (see Fig. 4). We then
identified several landmarks in both images (Fig. 5). Using this
minute hint, the semi-automatic algorithm could recover a plau-
sible deformation, even though the landmark information alone
(using e.g., thin-plate splines) would not have been enough [38].
We gave the weight 1.0 to all landmarks except the landmark
at the bottom left part of the skull which had a weight of 0.2.
This made the final positions of the landmarks coincide with
the target ones to within about 2 pixel for the least weighted
landmark and about 1 pixel for all the others.
Adding the spring term privileges likely solutions based on
our a priori knowledge and makes the problem better-posed.
The points need not to be image-dependent landmarks. For ex-
ample anchoring the four corners of the image prevents the so-
lution from degenerating. In this way, the springs play in part
the role of a regularization factor.
The landmarks are added when the automatic algorithm
cannot solve the problem by itself and an input from a human
expert is needed. For this reason, we decided to accept the land-
mark data as trustworthy and definitive. This is unlike in [58],
[59], where the landmarks come from an automatic process,
such as iterative closest-point algorithm (ICRP), and therefore
cannot be regarded as definitive. However, it is possible to give
a certain feedback to the expert, for example the value of the
criterion in landmark neighborhoods. This could be also used
to reject misplaced landmarks.
VI. I
MPLEMENTATION ISSUES
The purpose of this section is to describe some specific as-
pects of our implementation. These are mostly independent of
the main philosophy of the algorithm but can have a major im-
pact on its performance.
A. Explicit Derivatives
For the optimization algorithm, we need to calculate the par-
tial derivatives of
, as they form the gradient vector
and the Hessian matrix . Starting from equation (1),
we obtain the first partial derivatives
(6)
as well as the second partial derivatives
(7)
From (1) defining the SSD criterion, we get
and . The derivative
of the deformation function (4) is simply
. The deformation model is linear and all its
second derivatives are therefore zero; that is the reason for the
