A comparison of interpolation schemes using radial

MASCOT10-IMACS/ISGG Workshop
ULPGC, Las Palmas de Gran Canaria, Spain
A comparison of interpolation schemes
using radial basis functions and splines
for landmark-based image registration
Giampietro Allasia, Roberto Cavoretto, Alessandra De Rossi
Department of Mathematics “Giuseppe Peano”, University of Torino,
Via Carlo Alberto 10, I–10123 Torino, Italy
[email protected]; [email protected]; [email protected]
Abstract
In this paper we investigate landmark-based image registration using different interpolation schemes. More precisely, we make a comparison by
considering a few of the most known basis functions used in approximation theory, both globally and compactly supported. Among them, we
analyze the behavior of globally supported radial basis functions (RBFs)
such Gaussians and thin plate splines, as well as Wendland’s functions
and splines which have compact support. Numerical results are given to
compare them in the image registration context.
Keywords: Gaussians, Thin plate splines, Wendland’s functions, Shepard’s
formula, Lobachevsky splines.
1. Introduction
The problem of image registration consists in finding a transformation between
two or more images (or image data), called source and target images. In other
words, the aim is to determine a reasonable transformation, such that a transformed version of the first image is similar to the second one. The image
registration process may be based on a finite set of landmarks, i.e. sparse
data points located on images, usually not uniformly distributed, where each
landmark of the source image has to be mapped onto the corresponding landmark of the target image (see, e.g., [9, 10, 12]). Now, in order to give an idea
slightly more formal, we consider the two sets SN = {xi , i = 1, 2, . . . , N } and
TN = {ti , i = 1, 2, . . . , N } each containing N point-landmarks in the source and
1
target images, respectively. Thus, the registration problem involves a transformation F : Rm → Rm , such that
F(xi ) = ti ,
i = 1, 2, . . . , N,
(1)
where each coordinate of the transformation function F = (F1 , F2 , . . . , Fm )T
is separately calculated. This means that the interpolation problem involving
Fk : Rm → R is solved for each coordinate k = 1, 2, . . . , m, with the corresponding conditions Fk (xi ) = ti,k , i = 1, 2, . . . , N . However, this problem can
be formulated in the context of multivariate scattered data interpolation, and
solved by using radial basis functions (RBFs) (see, e.g., [13]). Since RBFs are
usually globally supported and a single landmark pair change may influence the
whole registration result, in the last decade several methods have been employed
to circumvent this disadvantage, such as compactly supported RBFs [7], elastic body splines [8], the modified inverse distance weighted method [4, 5], and
Lobachevsky spline functions [2]. A more specific application which involves registration and includes imaging techniques, such as computer tomography (CT)
and magnetic resonance imaging (MRI), can be found in [11].
In this paper we focus on the use of different interpolation schemes for
landmark-based image registration. In particular, we compare some of the most
known basis functions used in approximation theory, both globally and compactly supported. Thus, we analyze the behavior of globally supported radial
basis functions such Gaussians and thin plate splines, as well as Wendland’s
functions and Lobachevsky splines which have compact support (see [6, 13]).
Numerical experiments point out differences in accuracy of the considered methods.
The paper is organized as follows. In Section 2 we consider radial basis functions, like Gaussians and thin plate splines (globally supported) and Wendland’s
functions (compactly supported). Section 3 is devoted to briefly describe local
Shepard’s formula which uses RBFs as local approximants, whereas in Section
4 we focus on Lobachevsky splines. Finally, Section 5 contains numerical results
and conclusions.
2. Radial basis functions
2.1. Gaussians and thin plate splines
Using radial basis functions (RBFs), the general coordinate of the transformation function Fk (x), k = 1, 2, . . . m, which for convenience we will denote as
F (x) instead of Fk (x), is assumed to have the form
F (x) = ϕ(x) + p(x),
(2)
where ϕ is a radial basis function spanning an N -dimensional space of functions
m
≡ Pv−1 (Rm ), i.e. p
depending only on the source landmarks xi , and p ∈ Pv−1
m
is a sum of polynomials up to degree v − 1. The space Pv−1 = span{πk }U
k=1 has
2
dimension U = (m + v − 1)!/(m!(v − 1)!), which must be lower than N . Then
we can rewrite (2) in the extended form
F (x) =
N
X
i=1
ai Φ(||x − xi ||2 ) +
U
X
bk πk (x),
(3)
k=1
where Φ(|| · −xi ||2 ) is the radial basis function depending only on the Euclidean
distance r = || · −xi ||2 , and ai and bk are the coefficients to be determined.
Therefore, in order to compute the coefficients a = (a1 , a2 , . . . , aN )T and b =
(b1 , b2 , . . . , bU )T in (3), it is required to solve the following system of linear
equations
Ma + Qb = t
(4)
QT a = 0,
where M = {Φ(||xj − xi ||)} is a N × N matrix, Q = {πk (xj )} is a N × U
matrix, and t denotes the column vector of the j-th coordinate of the target
point-landmarks tj corresponding to the image point-landmarks xj . This result
is obtained by requiring that F satisfies the interpolation conditions (1) and the
PN
side conditions i=1 ai πk (xi ) = 0, for k = 1, 2, . . . , U , i.e. QT a = 0.
The most popular choices for Φ in landmark-based registration are
Φ(r)
Φ(r)
=
=
r2v−m log r,
2 2
e−α r ,
2v − m ∈ 2N,
(generalized thin plate spline)
(Gaussian)
where α ∈ R+ and v ∈ Z. The Gaussian (G) is a strictly positive definite function, while the thin plate spline (TPS) is a strictly conditionally positive definite
function of order v. The addition of a polynomial term of order v − 1 along with
side conditions, in order to guarantee existence and uniqueness of a solution in
the linear equation system (4), is required only for strictly conditionally positive
definite functions. It allows us to have an invertible interpolation matrix.
Note that the interpolation matrices generated by Gaussians and thin plate
splines are dense, since they are globally supported. Moreover, the phenomenon
of ill-conditioning could happen also in image registration applications, where it
is usually required to interpolate a relatively small number of point-landmarks
that are very close to each other.
2.2. Wendland’s functions
Wendland’s functions are a family of compactly supported radial functions obtained by using the truncated power function ϕs (r) = (1 − r)s+ , which is strictly
positive definite and radial on Rm for s ≥ ⌊m/2⌋ + 1, and repeatedly applying
the operator I. Then, we define the functions ϕm,k as follows
ϕm,k = I k ϕ⌊m/2⌋+k+1 .
The functions ϕm,k are strictly positive definite and radial on Rm , are all supported on the interval [0, 1] and have a polynomial representation. In particular,
3
they have the form
ϕm,k (r) =
pm,k (r),
0,
r ∈ [0, 1],
r > 1,
where pm,k is a univariate polynomial of degree ⌊m/2⌋ + 3k + 1. Moreover, any
other compactly supported C 2k polynomial functions, which are strictly positive
and radial on Rs , do not have a smaller polynomial degree; then, ϕ ∈ C 2k are
unique up to a constant factor, and the polynomial degree is minimal for given
space dimension m and smoothness 2k (see, e.g., [13]).
Since Wendland’s functions are compactly supported, the interpolation matrices can be made sparse by scaling the support of the basic function appropriately. Then, by setting Φ(r) = ϕm,k (r/c), we can define Wendland’s functions
depending on a shape parameter c ∈ R+ . Below, we list some of the most commonly used functions for m = 3 (that are strictly positive definite and radial on
Rm , m ≤ 3) along with their degree of smoothness 2k
Φ(r)
Φ(r)
Φ(r)
Φ(r)
=
=
=
=
2
(1 − cr)+ ,
4
(1 − cr)+ (4cr + 1) ,
6
(1 − cr)+ 35(cr)2 + 18cr + 3 ,
8
(1 − cr)+ 32(cr)3 + 25(cr)2 + 8cr + 1 ,
C0
C2
C4
C6
where (·)+ is the truncated power function. Hence, we can consider Wendland’s
functions Φ in the interpolant (3), omitting the polynomial part.
3. Modified Shepard’s formula
Shepard’s method is commonly used for multivariate interpolation and approximation of scattered data. In literature, many versions of such kind of formula
can be found (see [6] for an overview). Here, we consider a modified version
of Shepard’s method, which is local and uses radial basis functions as local
approximants.
Thus, referring to the context of the landmark-based image registration, we
define the k-th coordinate of a modified Shepard’s transformation F : Rm → Rm
in the form
Fk (x) =
N
X
¯ j (x),
Lj (x)W
k = 1, 2, . . . , m.
(5)
j=1
The nodal functions Lj , j = 1, 2, . . . , N , are local approximants at xj , constructed on the NL source landmarks closest to xj and satisfying the interpo¯ j , j = 1, 2, . . . , N , are
lation conditions Lj (xj ) = tj . The weight functions W
given by
¯ j (x) = P Wj (x)
,
W
N
k=1 Wk (x)
4
j = 1, 2, . . . , N,
(6)
with
Wj (x) = τ (x, xj )/α(x, xj ),
(7)
where α(·, xj ) = k · −xj k22 and τ (·, xj ) is a non-negative localizing function. It
is defined as follows
1, if x ∈ C(xj ; ρ),
τ (x, xj ) =
0, otherwise,
C(xj ; ρ) being a hypercube of centre at xj and side ρ.
Note that each component Fk , k = 1, 2, . . . , m, is evaluated at x considering
only the NW landmark points closest to x.
In (5) we take the local approximants Lj , j = 1, 2, . . . , N , as local RBF
interpolants, of the form
Lj (x) =
NL
X
i=1
ai Φ(kx − xi k2 ) +
U
X
bk πk (x),
(8)
k=1
where the radial basis functions Φ(k· − xi k2 ) here depend on the NL source
landmark points of the considered neighborhood of xj .
4. Spline functions
Let us consider a class of spline functions, called Lobachevsky splines, arising
in probability theory and converging to the normal density function. These
functions are compactly supported, have simple analytic expressions and are
strictly positive definite (see [1, 6]).
Lobachevsky spline transformations Fn : Rm → Rm are such that each
component is expressed in the following form
(Fn )k (x) =
N
X
k=1
ck fn∗ (x1 − x1k )fn∗ (x2 − x2k ) · · · fn∗ (xm − xmk ),
where n ∈ N is even, x = (x1 , x2 , . . . , xm ), (x1k , x2k , . . . , xmk ) ∈ Rm and
r
r
n
n
fn∗ (xi − xik ) = a
i = 1, 2, . . . , m,
fn a
(xi − xik ) ,
3
3
+
(9)
denote Lobachevsky spline functions, which are strictly positive definite as well
as their products (see [13]). The kth component of the spline transformation Fn
is a linear combination of the product of univariate shifted functions
fn∗√
(xi −xik ),
√
which are spline functions with compact supports [xik − 3n, xik + 3n] and
continuous up to n − 2 derivatives.
In particular, the function fn∗ given in (9) converges to the normal density
function with expectation 0 and standard deviation 1, i.e.
−x2 1
lim fn∗ (x) = √ exp
,
(10)
n→∞
2
2π
5
and moreover the convergence is uniform for all x ∈ R. To get explicit expressions of the fn∗ , we refer to fn defined by the convolution product
Z +a
Z +∞
1
fn−1 (x − u)du,
f1 (u)fn−1 (x − u)du =
fn (x) =
2a −a
−∞
because by definition f1 (x) = 1/(2a) for −a ≤ x ≤ a and f1 (x) = 0 elsewhere.
Setting x − u = t, we get the recursive formula
Z x+a
1
fn (x) =
fn−1 (t)dt,
n = 2, 3, . . .
(11)
2a x−a
From (11) it follows for −na ≤ x ≤ na
na+x
2a
X
1
k n
fn (x) =
(−1)
[x + (n − 2k)a]n−1 ,
(2a)n (n − 1)!
k
(12)
k=0
where ⌊·⌋ means the greatest integer less than or equal to the argument. Note
that the classic B-splines with equally spaced knots are directly connected to
the Lobachevsky splines. They are compactly supported functions with support
[−na, na] and smoothness n − 2, i.e. fn ∈ C n−2 [−na, na].
Considering the connection between fn (x) and fn∗ (x), the limit (10) becomes
r
r
−x2 1
n
n
.
fn a
x = √ exp
lim fn∗ (x) = lim a
n→∞
n→∞
3
3
2
2π
Noteworthy convergence properties are satisfied also by integrals and derivatives of the Lobachevsky splines. From the central limit theorem for the convergence in distribution we have
r
Z x r
Z x
Z x
−t2 1
n
n
∗
√ exp
a
fn (t)dt = lim
lim
dt,
fn a
t dt =
n→∞ −∞
n→∞ −∞
3
3
2
2π
−∞
while, assuming k ≤ n − 2 ∈ N, the asymptotic behavior of derivatives of the
Lobachevsky splines (see [3]) is given by
r
r
−x2 1
n
n
k ∗
k
lim D fn (x) = lim D a
.
fn a
x
= Dk √ exp
n→∞
n→∞
3
3
2
2π
It is well-known that using a Gaussian it is convenient to introduce a shape
parameter α by rescaling x to αx, i.e. considering the function exp(−α2 x2 ).
This trick can be conveniently applied for Lobachevsky splines too, considering
fn∗ (αx) instead of fn∗ (x). More or less localized basis functions can also be
obtained by considering fn (x) instead of fn∗ (x), by varying the value of the
parameter a. Numerical computations show that the approximations obtained
by fn∗ (x) or fn (x) are practically equivalent choosing conveniently the values of
the parameters α and a.
6
5. Numerical results and conclusions
In this section we make a comparison by considering all the above methods.
In order to test the different interpolation schemes, we have obtained several
numerical results. Here, for brevity, we refer to two test examples given in [7, 8].
They simulate typical medical cases where image portions scale on grids here
formed by 21 × 21 points. Source and target image landmarks are marked by a
circle (◦) and a star (⋆), respectively.
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
−0.2
−0.2
1.2
0
0.2
0.4
0.6
0.8
1
1.2
Figure 1: Source and target landmarks: square scaling (left) and circle scaling
(right).
In particular, the former case consists in an expansion of a square (case I)
using 32 landmarks (see Figure 1, left), while the latter is a contraction of a
circle (case II) using 20 equidistant landmarks located on the inner circle and, to
prevent an overall scaling, also 40 quasi-landmarks, i.e. landmarks at invariant
positions, at the outer circle in the source and target images (see Figure 1,
right).
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
−0.2
−0.2
1.2
0
0.2
0.4
0.6
0.8
1
1.2
Figure 2: Circle scaling: RBF-TPS (left) and Shepard-TPS with NL = 3, NW =
6 (right).
7
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Figure 3: Square scaling: RBF-G (left) and Shepard-G with NL = 10, NW = 10
(right) for α = 1.
More precisely, in the second case we simulate a resection model of a tumor in
surrounding elastic brain tissue: indeed, the outer circle corresponds to the skull
bone, which is assumed to be rigid; the inner circle represents the boundary of
the tumor, whereas the space between the inner and the outer circle is assumed
to be filled with elastic material, which corresponds to brain tissue. However, the
aim is here to determine a transformation function, which connects the points of
two images, so that the image is affected by the slightest possible deformation.
In Figure 2 we show registration results obtained by applying the global TPS
method (left) and Shepard’s method using a TPS as local approximant (right).
Figure 3 contains results using Gaussian (left) and Shepard’s method with a
Gaussian as nodal function (right).
Finally, in Figure 4 we compare registration results obtained by applying compactly supported functions such as Wendland’s function C 2 (left) and
Lobachevsky spline C 4 (right).
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
−0.2
1.2
0
0.2
0.4
0.6
0.8
1
1.2
Figure 4: Square scaling: Wendland-C 2 for c = 0.2 (left) and Lobachevsky-C 4
for α = 1 (right).
8
In general, these tests point out both higher level of accuracy and greater
stability of local methods, as confirmed by Table 1 and in accordance with
Figure 5.
RBF-G
RBF-TPS
Shepard-G
Shepard-TPS
2.1221E − 1
2.1409E − 1
1.7020E − 1
2.1260E − 1
Wendland-C 2
Wendland-C 4
Lobachevsky-C 2
Lobachevsky-C 4
1.2125E − 1
1.7671E − 1
1.4013E − 1
1.8532E − 1
Table 1: RMSEs (square scaling).
Figure 5: Comparison of different methods.
From this study we deduce that when the number N of point-landmarks is
relatively large Shepard’s type transformations are particularly suited in this
context (case I). On the other hand, if deformations are localized in a small
part of the image, then compactly supported Wendland’s functions as well as
spline functions give a greater accuracy and stability, locally preserving at the
same time properties of the image in examination (case II).
References
[1] G. Allasia, Scattered multivariate interpolation by a class of spline functions, in: J. Vigo-Aguiar et al. (Eds), Proceedings of the 9th International
Conference CMMSE09, 2009, vol. 1, pp. 73–79.
[2] G. Allasia, R. Cavoretto, A. De Rossi, B. Quatember, W. Recheis, M. Mayr,
S. Demertzis, Radial basis functions and splines for landmark-based registration of medical images, in: T. E. Simos, G. Psihoyios, C. Tsitouras
9
(Eds.), Proceedings of the ICNAAM10, AIP Conference Proceedings, vol.
1281, Melville, New York, 2010, pp. 716–719.
[3] R. Brinks, On the convergence of derivatives of B-splines to derivatives of
the Gaussian functions, Comput. Appl. Math. 27 (2008), 79–92.
[4] R. Cavoretto, A. De Rossi, A local IDW transformation algorithm for medical image registration, in: T.E. Simos, G. Psihoyios, Ch. Tsitouras (Eds.),
Proceedings of the ICNAAM08, AIP Conference Proceedings, vol. 1048,
Melville, New York, 2008, pp. 970–973.
[5] R. Cavoretto, A. De Rossi, B. Quatember, Landmark-based registration
using a local radial basis function transformation, JNAIAM J. Numer. Anal.
Ind. Appl. Math. (2010), accepted for publication.
[6] R. Cavoretto, Meshfree Approximation Methods, Algorithms and Applications, PhD Thesis, University of Turin, 2010.
[7] M. Fornefett, K. Rohr, H.S. Stiehl, Radial basis functions with compact
support for elastic registration of medical images, Image Vision Comput. 19
(2001), 87–96.
[8] J. Kohlrausch, K. Rohr, H.S. Stiehl, A new class of elastic body splines for
nonrigid registration of medical images, J. Math. Imaging Vision 23 (2005),
253–280.
[9] J. Modersitzki, Numerical Methods for Image Registration, Oxford Univ.
Press, Oxford, 2004.
[10] J. Modersitzki, FAIR: Flexible Algorithms for Image Registration, Fundam.
Algorithms 6, SIAM, Philadelphia, PA, 2009.
[11] B. Quatember, M. Mayr, W. Recheis, S. Demertzis, G. Allasia, A. De Rossi,
R. Cavoretto, E. Venturino, Geometric modelling and motion analysis of the
epicardial surface of the heart, Math. Comput. Simulation 81 (2010), 608–
622.
[12] K. Rohr, Landmark-Based Image Analysis, Using Geometric and Intensity
Models, Kluwer Academic Publishers, Norwell, MA, 2001.
[13] H. Wendland, Scattered Data Approximation, Cambridge Monogr. Appl.
Comput. Math., vol. 17, Cambridge Univ. Press, Cambridge, 2005.
10

Download Report