3.8 9LOBAL SPLINES 83
There are obviously many heuristic definitions one could give for the basis func-
tions. While nothing in the definition requires the basis functions to be rotationally sym-
metric, most heuristic definitions are based on that assumption. In this case, each basis
function becomes a function of the radial distance from the respective data point. For
example, hi(x,y ) = g(r) for some function g and radial distance r = (x -- Xi) 2 q- (y -- yi) 2 .
One of the most popular radially symmetric functions uses the Gaussian basis function
g (r) = e -r/ (3.8.3)
While it is possible to allow to vary with distance, in practice it is held constant. If is
chosen correctly, this method can provide quite satisfactory results. Otherwise, poor
results are obtained. In [Franke 79], it is suggested that c = 1.008 d/x-, where d is the
diameter of the point set and n is the number of points.
A second heuristically defined set of radial basis functions suggested by even more
reseamhers uses the popular B-spline. This has the advantage of having basis functions
with finite support, i.e., only a small neighborhood of a point needs to be considered
when evaluating its interpolated value. However, this method is still global in the sense
that there is a chain of interdependencies between all points that must be addressed when
evaluating the interpolating B-spline coefficients (see Section 5.2). Nevertheless, the
design matrix can be considerably sparser and better conditioned. The basis function can
then be taken to be
g (r) = 2 (1 -r/)+ - (1 - (2r/))+ (3.8.4)
where may be taken as 2.4192d/x-, d is the diameter of the point set, and n is the
number of points. Note that the + in the subscripts denote that the respective terms am
forced to zero when they are negative.
These heuristically defined basis functions have the intuitively attractive property
that they fall off with the distance to the data point. This reflects the belief that distant
points should exert less influence on local measurements. Again, nothing in the method
requires this property. In fact, while it may seem counter-intuitive, fine results have been
obtained with basis functions that increase with distance. An example of this behavior is
Hardy's multiquadratic splines [Hardy 71, 75]. Here the radial basis functions are taken
as
g (,9 = Wg' 0.8.5)
for a constant & Hardy suggests the value b = 0.815 m, where rn is the mean squared dis-
tance between points. Franke suggests the value b = 2.5 d/ffi', where d and n take on the
same definitions as before. In general, the quality of interpolation is better when the scat-
tered points are approximately evenly distributed. If the points tend to cluster along con-
tours, the results may become unsatisfactory.
Franke reports that of all the global basis functions tested in his research, Hardy's
multiquadratic was, overall, the most effective [Franke 79]. One important point to note
is that the design matrix for a radial function that increases with distance is generally ill-
conditioned. Intuitively, this is true because with the growth of the basis function, the
84 SPATIAL TRANSFORMATIONS
resulting interpolation tends to require large coefficients that delicately cancel out to pro-
duce the correct interpolated results.
The methods described thus far sutter from the difficulty of establishing good basis
functions. As we shall see, algorithms based on regularization theory define the global
splines by their properties rather than by an explicit set of basis functions. In general,
this makes it easier to determine the conditions for the uniqueness and existence of a
solution. Also, it is often easier to justify the interpolation techniques in terms of their
properties mtber than their basis functions. In addition, for some classes of functions and
dense data, more efficient solutions exist in computing global splines.
3.8.2. Regularlzation
The term regularization commonly refers to the process of finding a unique solution
to a problem by minimizing a weighted sum of two terms with the following properties.
The first term measures the extent to which the solution conforms to some preconceived
conditions, e.g., smoothness. The second term measures the extent to which the solution
agrees with the available data. Related techniques have been used by numerical analysts
for over twenty years [Atteia 66]. Generalizations are presented in [Duchon 76, 77] and
[Meinguet 79a, 79b]. It has only recently been introduced in computer vision [Grimson
81, 83]. We now briefly discuss its use in that field, with a warning that some of this
material assumes a level of mathematical sophistication suitable for advanced readers.
Techniques for interpolating (reconstructing) continuous functions from discrete
data have been proposed for many applications in computer vision, including visual sur-
face interpolation, inverse visual problems involving discontinuities, edge detection, and
motion estimation. All these applications can be shown to be ill-posed because the sup-
plied data is sparse, contains errors, and, in the absence of additional constraints, lies on
an infinite number of piecewise smooth surfaces. This precludes any prior guarantee that
a solution exists, or that it will be unique, or that it will be stable with respect to measure-
ment errors. Consequently, algorithms based on regularization theory [Tikhonov 77]
have been devised to systematically reformulate ill-posed problems into well-posed prob-
lems of variational calculus. Unlike the original problems, the variational principle for-
mulations are well-posed in the sense that a solution exists, is unique, and depends con-
tinuously on the data.
In practice, these algorithms do not exploit the full power of regularization but
rather use one central idea from that theory: an interpolation problem can be made
unique by restricting the class of admissible solutions and then requiring the result to
minimize a norm or semi-norm. The space of solutions is restricted by imposing global
variational principles stated in terms of smoothness constraints. These constraints,
known as stabilizing functionals, are regularizing terms which stabilize the solution.
They are treated together with penalty functionals, which bias the solution towards the
supplied data points, to form a class of viable solutions.
3.8 GLOBAL SPLINES 85
3.8.2.1. Grimson, 1981
Grimson first applied regularization techniques to the visual surface reconstruction
problem [Grimson 8I]. Instead of defining a particular surface family (e.g., planes) and
then fitting it to the data z (x,y), Grimson proceeded to fit an implicit surface by selecting
from among all of the interpolating surfaces f (x,y) the one that minimizes
E(f) = S(f)+P(f) 0.8.6)
= [ f f (f2 q-2fx2 q- fy)dxdy ] 2 + . [ f (xi,Yi)-z(xi,Yi)] 2
where E is an energy functional defined in terms of a stabilizing functional S and a
penalty functional P. The integral S is a measure of the deviation of the solution ffrom a
desirable smoothness constraint. The form of S given above is based on smoothness pro-
perties that are claimed to be consistent with the human visual system, hence the name
visual surface reconstruction. The summation P is a measure of the discrepancy between
the solution land the supplied data.
The surfaces that are computable from Eq. (3.8.6) are known in the literature as
thin-plate splines. Thin-plate interpolating surfaces had been considered in previous
work for the interpolation of aircraft wing deflections [Harder 72] and digital terrain
maps [Briggs 74]. The stabilizing functional S which these surfaces minimize is known
as a second-order Sobolev semi-norm. Grimson referred to it as the quadratic va7ation.
This semi-norm has the nice physical analogy of measu'ing the bending energy in a thin
plate. For instance, as S approaches zero, the plate becomes increasingly planar, e.g., no
bending.
The penalty measure P is a summation carried over all the data points. It permits
the surface to approximate the data z (x,y) in the least-squares sense. The scale parame-
ter 13 determines the relative importance between a close fit and smoothness. As [
approaches zero, the penalty term has greater latitude (since it is suppressed by 13) and S
becomes more critical to the minimization of E. This results in an approximating surface
f that smoothly passes near the data points. Forcing f to become an approximating sur-
face is appropriate when the data is known to contain errors. However, if f is to become
an interpolating surface, large values of [3 should be chosen to fit the data more closely.
This approach is based on minimizing E. If a minimizing solution exists, it will
satisfy the necessary condition that the first variation must vanish:
E (f) = S(f) + OP (f) = 0 (3.8.7)
The resulting partial differential equation is known as the Euler-Lagrange equation. It
governs the form of the energy-minimizing surface subject to the boundary conditions
that correspond to the given sparse data.
The Euler-Lagrange equation does not have an analytic solution in most practical
situations. This suggests a numerical solution applied to a discrete version of the prob-
lem. Grimson used the conjugate gradient method for approximation and the gradient
86 SPATIAL TRANSFORMATIONS
projection method for interpolation. They are both classical minimization algorithms
sharing the following advantages: they are iterative, numerically stable, and have parallel
variants which arc considered to be biologically feasible, e.g., consistent with the human
visual system. In addition, the gradient projection method has tbe advantage of being a
local technique. However, these methods also include the following disadvantages: the
rate of convergence is slow, a good criterion for terminating the iteration is lacking, and
the use of a grid representation for the discrete approximation precludes a viewpoint-
invariant solution.
There are two additional drawbacks to this approach that are due to its formulation.
First, the smoothing functional applies to the entire image, regardless of whether there
are genuine discontinuities that should not be smoothed. Second, the failure to detect
discontinuities gives rise to undesirable overshoots near large gradients. This is a man-
ifestation of a Gibbs or Mach band phenomenon across a discontinuity in the surfaces.
These problems arc addressed in the work of subsequent researchers.
3.8.2.2. Terzopoulos, 1984
Terzopoulos extended Grimson's results in two important ways: he combined the
thin-plate model together with mcmbrune elements to accommodate discontinuities, and
he applied multigrid relaxation techniques to accelerate convergence. Terzopoulos finds
the unique surface f minimizing E where
and
P (f ) = i, ai(Li[f ]-Li[z ]-ei)2 (3.8.9)
The stabilizing functional S, now referred to as a controlled-continuity stabilizer, is an
integral measure that augments a thin-plate spline with a membrane model. The penalty
functional P is again defined as a weighted Euclidean norm. It is expressed in terms of
the measurement functionals Li, the associated measurement errors El, and nonnegative
real-valued weights o; i. L i can be used to obtain point values and derivative information.
In Eq. (3.8.8), p (x,y) and x (x,y) are real-valued weighting functions whose range is
[0,1]. They are referred to as continuity control functions, determining the local con-
tinuity of the surface at any point. An interpretation of x is surface tension, while that of
p is surface cohesion. Their correspondence with thin-plate and membrane splines is
given below.
lim x(x,y)-o0 S (f) --> membrane spline
lim(x,y)-! S(f) -- thin-platespline
limo(x,y)_40 S(f) -- discontinuous surface
3.8 GLOBAL SPLINES 87
A thin-plate spline is characterized as a C I surface which is continuous and has continu-
ous first derivatives. A membrane spline is a C O surface that need only be continuous.
Membrane splines are introduced to account for discontinuities in orientations, e.g.,
comers and creases. This reduces the Gibbs phenomena (oscillations) near large gra-
dients by preventing the smoothness condition to apply over discontinuities.
Terzopoulos formulates the discrete problem as a finite element system. Although
the finite element method can be solved by using iterative techniques such as relaxation,
the process is slow and convergence is not always guaranteed. Terzopoulos used the
Gauss-Seidel algorithm, which is a special case of the Successive Over-Relaxation
(SOR) algorithm. He greatly enhanced the SOR algorithm by adapting multigrid relaxa-
tion techniques developed for solving elliptic partial differential equations. This method
computes a coarse approximation to the solution surface, uses it to initiate the iterative
solution of a finer approximation, and then uses this finer approximation to refine the
coarse estimate. This pmcedure cycles through to completion at which point we reach a
smooth energy-minimizing surface that interpolates (or approximates) the sparse data.
The principle of multigrid operation is consistent with the use of pyramid and mul-
tiresolution data structures in other fields. At a single level, the SOR algorithm rapidly
smooths away high frequency error, leaving a residual low frequency error to decay
slowly. The rate of this decay is dependent on the frequency: high frequencies are
removed locally (fast), while low frequencies require long distance propagation taken
over many iterations (slow). Dramatic speed improvements are made. possible by pm-
jecting the low frequencies onto a coarse grid, where they become high frequencies with
respect to the new grid. This exploits the fact that neighbor-to-neighbor communication
on a coarse grid actually covers much more ground per iteration than on a fine grid. An
adaptive scheme switches the relaxation process between levels according to the fre-
quency content of the error signal, as measured by a local Fourier transformation.
We now consider some benefits and drawbacks of Terzopoulos' approach. The
advantages include: the methods are far more computationally efficient over those of
Grimson, discontinuities (given a priori) can be handled, error can be measured dif-
ferently at each point, a convenient pyramid structure is used for surface representation,
and local computations make this approach biologically feasible.
Some of the disadvantages include: a good criterion for terminating the iteration is
lacking, the use of a grid representation for the discrete approximation precludes a
viewpoint-invariant solution, there is slower convergence away from the supplied data
points and near the grid boundary, and the numerical stability and convergence rates for
the multigrid approach are not apparent.
3.8.2.3. Discontinuity Detection
Techniques for increasing the speed and accuracy of this approach have been inves-
tigated by Jou and Bovik [Jou 89]. They place emphasis on early localization of surface
discontinuities to accelerate the process of minimizing the surface energy with a finite
element approximation. Terzopoulos suggests that discontinuities are associated with
llll
88 SPATIAL TRANSFORMATIONS
places of high tension in the thin-plate spline. While this does detect discontinuities, it is
prone to error since there is no one-to-one correspondence between the two. For
instance, it is possible to have many locations of high tension for a single continuity
beacuse of the oscillatory behavior of Gibbs effect. On the other hand, it is possible to
miss areas of high tension if the data points around a discontinuity are sparse.
Grimson and Pavlidis describe a method for discontinuity detection based on a
hypothesis testing technique. At each point, they compute a planar approximation of the
data and use the statistics of the differences between the actual values and the approxi-
mations for detection of both steps and creases [Grimson 85]. If the distribution of the
residual error appears random, then the hypothesis that there is no discontinuity is
accepted. Otherwise, if systematic trends are found, then a discontinuity has been
detected.
Blake and Zisserman have developed a technique based on "weak" continuity con-
straints, in which continuity-like constraints are usually enforced but can be broken if a
suitable cost is incurred [Blake 87]. Their method is viewpoint-invariant and robustly
detects and localizes discontinuities in the presence of noise. Computing the global
minimum is difficult because invariant schemes incorporating weak continuity con-
straints have non-convex cost functions that are not amenable to naive descent algo-
rithms. Furthermore, these schemes do not give rise to linear equations. Consequently,
they introduce the Graduated Non-Convexity (GNC) Algorithm as an approximate
method for obtaining the global minimum. The GNC approach is a heuristic that minim-
izes the objective function by a series of convex curves that increasingly refine the
approximation of the function near the global minimum. The initial curve localizes the
area of the solution and the subsequent curves establish the value precisely.
3.8,2.4. Boult and Kender, 1986
Boult and Kender examine four formalizations of the visual surface reconstruction
problem and give alternate realizations for finding a surface that minimizes a functional
[Boult 86a]. These methods are:
1. Discretization of the problem using variational principles and then discrete minimi-
zation using classical minimization techniques, as in [Grimson 81].
2. Discretization of a partial differential equation formulation of the problem, again
using a variational approach, and then use of discrete finite element approximation
solved with a multigrid approach, as in [Terzopoulos 84].
3. Direct calculation using semi-reproducing kemel splines, as in [Duchon 76].
4. Direct calculation using quotient reproducing kernel splines, as in [Meinguet 79].
The authors conclude that reproducing kernel splines are a good method for interpo-
lating sparse data. The major computational component of the method is the solution of
a dense linear system of equations. They cite the following advantages: the solution of
the linear system is well-understood, the algorithm results in functional forms for the sur-
face allowing symbolic calculations (e.g., differentiation or integration), there is no prob-
lem with slower convergence away from information points or near the boundary, the
3.8 GLOBAL SPLINES 89
algorithm can efficiently allow updating the information (e.g., adding/deleting/modifying
data points), no iteration is needed since the computation depends only on the number of
information points (not their values), and the method is more efficient than those of
Grimson or Terzopoulos for sparse data.
The disadvantages include: the resulting linear system is dense and indefinite which
limits the approach with which it can be solved, the reproducing kemels may be difficult
to derive, and the method may not be biologically feasible due to the implicit global
communication demands. Full details about this method can be found in [Boult 86b]. In
what follows, we present a brief derivation of the semi-reproducing kernel spline
approach. The name of this method is due to formal mathematical definitions. These
details will be suppressed insofar as they lie outside the scope of this book.
The use of semi-reproducing kernel spline allows us to interpolate data by almost
the same four steps as used in the computation of global splines defined through heuristic
basis functions. This is in contrast to regularization, which employs a discrete minimiza-
tion method akin to that used by Terzopoulos. Unlike the basis functions suggested ear-
lier, though, the basis functions used to define the semi-reproducing kernel splines com-
pute a surface with minimal energy, as defined in Eq. (3.8.6). To interpolate M data
points for mapping function U, the expression is
M+3
U(x,y)= aihi(x,y)
where the basis functions h i are
hi(x,y) = 0 ß [(x-xi) 2 + (y _yl)2]. log[(x-xi) 2 + (y _yi)2], i = 1 ..... M
hM+i(x,y) = 1
hM+2(x,y) = x (3.8.10)
hst+3(x,y) = y
for a constant 0 (see below).
The above expression for U has more basis functions than data points. These extra
terms, hM+i, i = 1, 2, 3, are called the basis functions of the null space, which has dimen-
sion d = 3. They span those functions which can be added to any other function without
changing the energy term (Eq. 3.8.6) of that function. They are introduced here because
they determine the components of the low-order variations in the data which are not con-
strained by the smoothness functional (the norm). In our case, since the norm is twice-
differentiable in x and y, the low-order variation is a plane and our null space must have
dimension d = 3.
Since we have additional basis functions, the design matrix is insufficient to define
the coefficients of the interpolation spline. Instead, the M +d basis function coefficients
can be determined from the solution of (M + d) x (M + d) dense linear system:
90 s PATIA L TRANSFORMAT IONS
where
U2
0
0
0
a2
= A aM
aM+
aM+
aM+
(3.8.11)
Ai'j = hi(xj'YJ) for i <(M +d), j
Ai,j = -1 +hi(xj,yj) fori=j
Ai.j = hj(xi,Yi) for i
Ai,j = 0 fori>M, j>M
The system can be shown to be nonsingular if the data spans the null space. For this
case, the data must contain at least three non-collinear points. Due to the mathematical
properties of the basis functions, the above spline is referred to as an interpolating semi-
reproducing kernel spline. Note that the given above corresponds to that used in the
expression for the energy functional in Eq. (3.8.6). As it approaches infinity, the system
determines the coefficients of the interpolating spline of minimal norm.
**Share with your friends:** |