Digital image warping

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)

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)

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

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

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

0

0

a2

aM+

aM+

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.