On the Bijectivity of Thin-plate Splines.
Anders P Erikson and Kalle
˚
Astr¨om
Centre for Mathematical Sciences, Lund University
Lund, Sweden
{anderspe,kalle}@maths.lth.se
Abstract
The thin-plate spline (TPS) has been widely used
in a number of areas such as image warping, shape
analysis and scattered data interpolation. Intro-
duced by Bookstein [1], it is a natural interpolating
function in two dimensions, parameterized by a fi-
nite number of landmarks. However, even though
the thin-plate spline has a very elegant intuitive in-
terpretation as well as mathematical formulation it
has no inherent restriction to prevent folding, i.e.
a non-bijective interpolating function. In this pa-
per we discuss some of the properties of the set
of parameterizations that form bijective thin-plate
splines, such as convexity and boundness. Meth-
ods for finding sufficient as well as necessary con-
ditions for bijectivity are also presented.
1 Introduction
The thin-plate spline (TPS) is a natural choice of
interpolating function in two dimensions and has
been a commonly used tool in computer vision for
over a decade. Its attraction might include a the el-
egant mathematical formulation along with a very
natural interpretation, as thin-plate splines can be
viewed as modeling the bending of a thin metal
plate under point constraints. However there are
disadvantages, firstly it might be a very computa-
tionally expensive tool. Secondly, a TPS-mapping
has no inherent restriction to prevent folding from
occurring, i.e. there are no constraints on when the
mapping is non-bijective. Since, in the context of
computer vision, non-bijective deformations of im-
ages are quite uncommon in natural images (see
figs. 1,2), it is the latter issue that will be ad-
dressed in this paper.
..
..
Figure 1: Bijective mapping.
The motivation for this work is that we wanted
to perform registration between images using thin-
plate splines. That is, the optimization problem
of finding a bijective deformation that minimizes
some similarity function between image pairs. There-
fore, to obtain the optimization constraints we need
to learn more about this set of bijective thin-plate
spline deformations. Furthermore, as are working
within an optimization framework we are also in-
terested in other properties of the set defined by
the bijectivity constraints, such as if this set is con-
vex, bounded and/or star-shaped.
2 The Thin-Plate Spline and
bijectivity constraints.
The thin-plate spline mapping g : R
2
→ R
2
given
by a set of k control points (or landmarks) T and a
set of k destination points Y such that g(T ) = Y is,
with the metal plate analogy in mind, the transfor-
mation that minimizes the roughness penalty (or
bending energy)
Z
R
2
(
∂
2
g
∂x
2
)dx (1)
It has been shown by Kent and Mardia [2] that
..
..
Figure 2: Non-bijective mapping.
a bivariate function φ in the form (for details see
[1])
φ
T,Y
(x) = (φ
1
(x), φ
2
(x))
T
= c + Ax +
+W (σ(x − T
1
), ..., σ(x − T
k
)) =
=
W
T
c A
s(x)
1
x
(2)
with
s(h) = ||h||
2
log(||h||), (3)
W
T
c A
=
Y
T
0 0
Γ
−1
T
(4)
minimizes eq. (1)
Combining this, the transformation can be writ-
ten as
φ
T,Y
(x) =
Y
T
0 0
Γ
−1
T
s(x)
1
x
= Y
T
N
T
(x) (5)
This gives us a deformation φ that for a fixed set
of control points T is parameterized (linearly) by
the destination points Y. And we are interested
in knowing for which Y do we get a bijective de-
formation, i.e the set
Ω
T
= {Y ∈ R
2k
|φ
T,Y
(x) is bijective }
Such a mapping φ : R → R is bijective if and
only if its functional determinant |J(φ)| is non-
zero. Using (5) we get
|J(φ)| =
∂φ
1
∂x
1
∂φ
1
∂x
2
∂φ
2
∂x
1
∂φ
2
∂x
2
= ... =
ˆ
Y
T
B
T
(x)
ˆ
Y (6)
with
B
T
(x) =
0 D(x)
−D(x) 0
(7)
D(x) = b
x
1
(x)b
x
2
(x)
T
− b
x
2
(x)b
x
1
(x)
T
b
x
1
(x) = Γ
11
∂s
∂x
1
(x) + Γ
1
b
x
2
(x) = Γ
11
∂s
∂x
2
(x) + Γ
2
(8)
Where
ˆ
Y is the vectorized version of the k-by-2
matrix Y. B
T
(x) is a 2k-by-2k symmetric, indef-
inite, rank-four matrix with zeros in the diagonal
and non-zero eigenvalues λ
T
(x), λ
T
(x), −λ
T
(x),
−λ
T
(x). So for each point x ∈ R
2
we get a quadratic
constraint on Y, (
ˆ
Y
T
B
T
(x)
ˆ
Y 6= 0) for local bijec-
tivity. For φ to be globally bijective this constraint
must either be > 0, ∀x ∈ R
2
or < 0, ∀x ∈ R
2
.
Ω
T
can thus be written
Ω
T
= {Y ∈ R
2k
|
ˆ
Y
T
B
T
(x)
ˆ
Y > 0, ∀x ∈ R
2
or
ˆ
Y
T
B
T
(x)
ˆ
Y < 0, ∀x ∈ R
2
}
Seeing that if Y ∈ Ω
T
then −Y ∈ Ω
T
, it does,
without loss of generality, suffice to examine
Ω
+
T
= {Y ∈ R
2k
|
ˆ
Y
T
B
T
(x)
ˆ
Y > 0, ∀x ∈ R
2
}
(with Ω
−
T
defined similarly, we can write Ω
T
=
Ω
+
T
∪ Ω
−
T
.)
0 5 10 15 20 25 30 35 40
0
5
10
15
20
25
30
35
40
..
−2000 −1000 0 1000 2000 3000
−2000
0
2000
4000
6000
8000
10000
Figure 3: The resulting quadratic constraints on a
subset of Y imposed by three points in R
2
. (x):
The control points, (*): The three arbitrarily cho-
sen points in R
2
.
So the sought after set is the intersection of an infi-
nite number of indefinite (and non-convex) quadratic
forms each given by eq.(6).
3 Convexity of Ω
+
T
and other
Properties.
As previously mentioned, we wish to use these bi-
jectivity constraints within an optimization con-
text and that therefore the convexity of Ω
+
T
is of
great interest. In general, one would not expect
that the intersection of non-convex sets would re-
sult in a convex set. Empirical observations made
seem to agree with this suspicion. For certain sim-
ple control configurations T, non-convexity of Ω
+
T
can quite easily be shown. A natural continuation
is to ask if imposing some further restrictions on Y
can result in Ω
+
T
becoming convex? For instance,
if one allows all but three (linearly independent)
points in Y freedom, (one can view this as elimi-
nating affine transformations of Y from Ω
+
T
). Once
again experiments have indicated that our set still
is non-convex. Other similar constraints on the
destination points, such as “locking” points close
to the border of the convex hull of Y, have also
resulted in a non-convex set. However, instances
when convexity occur do exist. If only one control
point is let loose, then Ω
+
T
, owing to the special
form of B
T
(x) in eq.(6), becomes a polytope, see
fig. 4.
0 20 40 60 80 100 120
−30
−20
−10
0
10
20
30
40
50
Figure 4: The resulting bijectivity constraints
when only allowing one destination point to move.
Even though many of these observations remain to
be proven they do provide strong indications that
Ω
+
T
in general is a non-convex set.
Regarding the boundness of Ω
+
T
. It is known
that for affine transformations a(Y) the following
holds
φ
T,a(Y)
(x) = a(φ
T,Y
(x)) (9)
So if Y ∈ Ω
+
T
and a(Y) nonsingular then a(Y) ∈
Ω
+
T
. Hence Ω
+
T
is clearly unbounded. We are how-
ever convinced that it can be proven that under an
affine elimination, as previously described, the set
becomes bounded.
Experiments have also indicated that Ω
+
T
might
be star-shaped, that is that the intersection be-
tween any line passing through T and Ω
+
T
is a
convex set. Further study of this property is still
required though.
Although the properties discussed in this sec-
tion have mainly been observed through experi-
mentation on a very limited number of different
control- and destination configurations it does pro-
vide us with the opportunity to get a better un-
derstanding of what kind of object we are working
with. Which, in summation, is a high-dimensional,
non-convex, unbounded monstrosity defined by an
infinite number of indefinite quadratic constraints
(which in addition also are non-linear in x).
4 Necessary and Sufficient Con-
ditions for Bijectivity
Given the complexity of the set of bijective thin-
plate spline deformations, the task finding a defin-
ing expression for it analytically is a formidable
one. One approach could for instance be to try
and solve the envelope [3] equations
Y
T
B
T
(x)Y = 0 (10)
Y
T
B(x)
T
(x)
0
x
1
Y = 0 (11)
Y
T
B(x)
T
(x)
0
x
2
Y = 0 (12)
It can be shown that since the constrain tangents
Ω
+
T
at its boundary ∂Ω
+
T
, the points on ∂Ω
+
T
is
a subset of the envelope of the family of all the
quadratic functions.
Instead we have chosen to use numerical meth-
ods to derive conditions on Ω
+
T
. By finding the
minimum-volume ellipsoid E
1
covering Ω
+
T
and the
maximum-volume ellipsoid E
2
inscribed in Ω
+
T
we
obtain a necessary and a sufficient conditions on
Y. That is
Y /∈ E
1
⇒ Y /∈ Ω
+
T
(13)
Y ∈ E
2
⇒ Y ∈ Ω
+
T
(14)
Finding such extremal volume ellipsoids can be
formulated as optimization problems [4, 5]. But
since we have finitely many variables and infinite
number of constraints we have a semi-infinite pro-
gram on our hand [7]. In order to avoid this we
simply approximate Ω
+
T
by the intersection of a
finite subset of these constraints.
˜
Ω
+
T
= {Y ∈ R
2k
|
ˆ
Y
T
B
T
(x)
ˆ
Y > 0,
for a finite number of x ∈ R
2
}
With E
1
and E
2
defined by
E
i
= {p ∈ R
n
| p
T
A
i
p +2b
T
i
p +c
i
≥ 0} and the bi-
jectivity constraints in the form {Y ∈ R
n
|Y
T
F (x)Y +
2g
T
Y + h(x) ≥ 0} we proceed.
The minimum volume ellipsoid is always a con-
vex optimization problem regardless the set it cov-
ers. Using the introduced notation. E
1
will be the
solution to the following convex program
minimize − log det E
1
s.t.
A
1
−
P
τ
i
F (x
i
) A
1
b
1
−
P
τ
i
g(x
i
)
(A
1
b
1
−
P
τ
i
g(x
i
))
T
c
i
−
P
τ
i
h(x
i
)
0
The maximum volume inscribed ellipsoid is a
convex optimization problem if the covering set it-
self is convex. Since Ω
+
T
in general is non-convex
−40 −35 −30 −25 −20 −15 −10 −5 0
−20
−10
0
10
20
30
40
50
60
70
80
90
Figure 5: Example of a subset of Ω
+
T
(correspond-
ing to letting Y free in two dimensions).
−40 −35 −30 −25 −20 −15 −10 −5 0
−50
0
50
100
Figure 6: Example of necessary (larger ellipsoid)
and sufficient conditions (smaller ellipsoid) on the
same subset as above.
we must employ other optimization algorithms. Nev-
ertheless, E
2
will be the solution to the following
optimization problem
maximize − log det E
2
s.t.
F (x
i
) − A
2
g(x
i
) − A
2
b
2
(g(x
i
) − A
2
b
2
)
T
h(x
i
) − c
i
0
5 Conclusion and Future Work
Even though this paper very much describes a work
in progress, it does contain a formulation of how to
characterize the set of bijective thin-plate splines.
It also includes a discussion of some experimentally
derived indications of some other properties of this
set as well as a method for finding necessary and
sufficient conditions for bijectivity. Future work in-
cludes finding such conditions analytically as well
as attempting to prove its convexity and bound-
ness properties.
References
[1] F. L. Bookstein, ”Principal Warps: Thin-
Plate Splines and the Decomposition of Defor-
mations”, IEEE trans. Pattern Analysis and
Machine Intelligence, Vol. 11, No. 6, 1989.
[2] Kent J. T. and Mardia K. V., ”The link be-
tween kriging and thin-plate splines.”, Proba-
bility, statistics and Optimization: a Tribute
to Peter Whittle, p325-329.
[3] Bruce J.W. and Giblin P.J., “Curves and Sin-
gularities”, Cambridge University Press.
[4] Vandenberghe L. and Boyd S., “Convex Op-
timization”, Cambridge University Press.
[5] Vandenberghe L., Boyd S. and Wu S.W.,
“Determinant Maximization with Linear Ma-
trix Inequality Constraints”, SIAM Journal
on Matrix Analysis and Applications, Volume
19, Number 2, pp. 499-533.
[6] Vandenberghe L. and Boyd S., “Semidefinite
Programming”, SIAM Review, Volume 38,
Number 1, pp. 49-95.
[7] Hettich R. and Kortanek K.O., ”Semi-infinite
programming: theory, methods, and applica-
tions”, SIAM Review, Volume 35, Number 3,
pp. 380-429.
