Digital image warping

The global effects of the polynomial mapping will not account for high-frequency defor-

mations that are local in nature. Since most sensor-related errors tend to be low-

frequency, modeling the spatial Ixansformation with low-order polynomials appears

justified. Common values of N that have been used in the polynomials of Eq. (3.6.1)

include N = 1 [Steiner 77], N =2 [Nack 77], N =3 [Van Wie 77], and N =4 [Leckie 80].

For many practical problems, a second-degree (N = 2) approximation has been shown to

be adequate [Lillestrand 72].

Note that a first-degree (N= 1) bivariate polynomial defines those mapping func-

tions that are exactly given by a general 3 x 3 affine transformation matrix. As discussed

3.6 POLYNOMIAL TRANSFORMATIONS 63

in the Orevious section, these polynomials characterize common physical distortions, i.e.,

affine transformations. When the viewing geometry is known in advance, the selection

of the polynomial coefficients is determined directly from the scale, translation, rotation,

and skew specifications. This is an example typical of computer graphics. For example,

given a mathematical model of the world, including objects and the viewing plane, it is

relatively straightforward to cascade transformation matrices such that a series of projec-

tions onto the viewing plane can be realized.

In the fields of remote sensing, medical imaging, and computer vision, however, the

task of computing the spatial Ixansformation is not so straightforward. In the vast major-

ity of applications, the polynomial coefficients are not given directly. Instead, spatial

information is supplied by means of tiepoints or control points, corresponding positions

in the input and output images whose coordinates can be defined precisely. In these

cases, the central task of the spatial transformation stage is to infer the coefficients of the

polynomial that models the unknown distortion. Once these coefficients are known, Eq.

(3.6.1) is fully specified and it is used to map the observed (x,y) points onto the reference

(u,v) coordinate system. The process of using tiepoints to infer the polynomial

coefficients necessary for defining the spatial Ixansformation is known as spatial interpo-

lation [Green 89].

Rather than apply the mapping functions over the entire set of points, an interpola-

tion grid is often introduced to reduce the computational complexity. This method

evaluates the mapping function at a relatively sparse set of grid, or mesh, points. The

spatial correspondence of points internal to the mesh is computed by bilinear interpola-

tion from the corner points [Bernstein 76] or by fitting cubic surface patches to the mesh

[Goshtasby 89].

3.6.1. Inferring Polynomial Coefficients

Auxiliary information is needed to determine the polynomial coefficients. This

information includes reseau marks, platform attitude and altitude data, and ground con-

trol points. Reseau marks are small cruciform markings inscribed on the faceplate of the

sensor. Since the locations of the reseau marks can be accurately calibrated, the meas-

ured differences between their Ixue locations and imaged (distorted) locations yields a

sparse sensor distortion mapping. This accounts for the internal errors.

Extemal errors can be directly characterized from platform attitude, altitude, and

ephemerides data. However, this data is not generally precisely known. Consequently,

ground conlxol points are used to determine the external error. A ground control point

(GCP) is an identifiable natural landmark detectable in a scene, whose location and

elevation are known precisely. This establishes a correspondence between image coordi-

nates (measured in rows and columns) and map coordinates (measured in

latitude/longitude angles, feet, or meters). Typical GCPs include airports, highway inter-

sections, land-water interfaces, and geological pattems [Bernstein 71, 76].

A number of these points are located and differences between their observed and

actual locations are used to characterize the external error component. Together with the

64 SPATIAL TRANSFORMATiONS $,6 P OLYNOMIA L TRANSFORMATIONS 65

internal distortion function, this serves to fully define the spatial transformation that

inverts the distortions present in the input image, yielding a corrected output image.

Since there are more ground control points than undetermined polynomial coefficients, a

least-squared-error fit is used. In the discussion that follows, we describe several tech-

niques to solve for the unknown polynomial coefficients. They include the pseudoin-

verse solution, least-squares with ordinary and orthogonal polynomials, and weighted

least-squares with orthogonal polynomials.

3.6.2. Pseudoinverse Solution

Let a correspondence be established between M points in the observed and refer-

ence images. The spatial transformation that approximates this correspondence is chosen

to be a polynomial of degree N. In two variables (e.g., x and y), such a polynomial has K

coefficients where

v v-i (N+ 1) (N+2)

i:=o j=o

For example, a second-degree approximation requires only six coefficients to be solved.

In this case, N = 2 and K =6. Wc thus have

u2 t x2 y2 x2y2 x2 2 y ] ao

u3 1 x3 Y3 X3Y3 x32 Y32 a01 (3.6.2)

. . . all

a20

1 XM YM XMYM X4 y21212M l a02

written in matrix notation as

 U = WA (3.6.3)

V=WB

has dimensions M xK, it is not a square matrix and thus it has no inverse. Instead, we

first multiply both sides by W :e before isolating the desired A and B vectors. This serves

to cast W into a K xK square matrix that may be readily inverted [Wong 77]. This gives

us

WrU = wTwA (3.6.4)

equations. For notafional convenience, we omit the limits of the summations and the

associated subscripts. Note that all summations run from k =1 to M, operating on data x,,

y,, and u,.

Zu M x y Zxy X 2 Zy 2 aoo

Zxu Zx Zx 2 Z Zx2Y Zx 3 Z 2 a0

[U  x2y 2 x2y2 x3y 3 al l

A similar pedure is peffoed for v and b 0. Solving for A and B, we have

A = (wrw)-wru (3.6.6)

B = (WrW)-iWrV

This technique is own as tbe pseoinverse soluffon to the line least-uares prob-

lem. It leaves us with K-element vectors A and B, the polynomifl cfficients for e U

d V mapping functions, spectively.

Unfonately, this meth suffers om several problems. The pfima difficulW

lies in multiplying . (3.6.3) with W T. is sques tbe condition number, theby

reducing the precision of the coefficients by one half. As a result, alternate solutions e

recommendS. For example, it is preferable to compute a decomposition of W rather

th solve for its inverse dirfly. e ader is refe to the linear algebra literathe

for a discussion of singul value decomposition, and LU and QR decomposition tech-

niques [Sang 80]. An exposition can  found in [ess 88], whe emphasis is given to

an intuitive understanding of the benefits and awbacks of these thniques. e text

also provides useful souse ce written in the C proming language. (Versions of

tbe book with Pascal and Fo pro,ams e available as well).

3.6.3. Least-Squares With Ordinary Polomials

The pseudoinverse solution pves to  identical to that of the classic least-sques

foulation with ordiny polynomifls. Although  approaches sha me of e

same problems, the least-sques meth is discuss he due to its prominence in the

solution of overdetein systems of line uafions. Funheore, it can  alte to

yield a stable clos-fo solution for the unown coefficients, as descfi in the next

section.

Refeng back to Eq. (3.6.1) with N = 2, cfficients aij e detein by miz-

ing

M

M

k=l

=  [aoo+aoxk+aotYk+anxy+aeox+ao2Y-Uk] 2

k=l

This is achiev by deteining the ptial derivatives of E with respect to cfficients

aij 2 5n  = 0 (3.6.8)

By considering e pial derivative of E with respect to all six efficients, we obta

Eq. (3.6.9), the following system of linear equations. For notafional convenience, we

have omitted the limits of summation and subscripts.

u =aM +aoX +amy +a11 +a:oX  +ao2y 

XU =ax +a10x 2 +a01 +allX2y +a20x 3 +a02 2

yu =ay +a10 +a01Y 2 +all 2 +a20x2y +a02Y 3

u =a +a10x2y +a01 2 +allx2y 2 +a20x3y +a02 3

X2U =ax 2 +a10x 3 +aolx2y +aux3y +ae0x 4 +ao2x2y 2

y2u =ay 2 +al0  +a01Y 3 +all 3 +a20x2y  +a02Y n

This is a symmec 6 x 6 system of linear uations, whose mefficients e all sum-

mations from k= 1 to M which e evaluated from the original dam. By inspection, this

result is equivalent to Eq. (3.6.5), the system of equations defiv elier in the pseudoin-

verse solution. Known as the nodal eqtiov, this system of line uations c be

compactly expressed in the following notation.

a x i 'x I m =   y (3.6.10)

i=0 j kk= J

for l = 0, ..., N and m = 0, ..., N - l. Note that (i,j) e running indices ong rows, where

each row is associated with a nstant (/,m) pair.

Tbe least-squs predure operates on  overdetein system of line ua-

tions, (i.e., M dam points are us to deteine K cfficients, where M >. As a

result, we have only an approximating mapping function. If we substitute the (x,y) con-

ol point crdinates back into the polynomial, the mputed results will lie ne, but

will generally not coincide exactly with their countet (u,v) crdinates. Stated intui-

tively, the polynomial of oer K must yield the st coromise among all M nol

points. As we have seen, the st fit is deteined by minimizing tbe sum of e

squar-eor [U (x,yn) - uk]  for k = l ..... M. Of course, if M = K then no compromise

3.6 POLYNOMIAL TRANSFORMATIONS 67

is necessary and the polynomial mapping function will interpolate the points, actually

passing through them.

The results of the least-squares approach and the pseudoinverse solution will yield

those coefficients that best approximate the true mapping function as given by the M con-

trol points. We can refine the approximation by considering the error at each data point

and throwing away that point which has maximum deviation from the model. The

coefficients are then recomputed with M- 1 points. This process can be repeated until

we have K points remaining (matching the number of polynomial coefficients), or until

some error threshold has been reached. In this manner, the polynomial coefficients are

made to more closely fit an increasingly consistent set of data. Although this method

requires more computation, it is recommended when noisy data are known to be present.

3.6.4. Least-Squares With Orthogonal Polynomials

The normal equations in Eq. (3.6.10) can be solved if M is greater than K. As K

increases, however, solving a large system of equations becomes unstable and inaccurate.

Numerical accuracy is further hampered by possible linear dependencies among the

equations. As a result, the direct solution of the normal equations is generally not the

best way to find the least-squares solution. In this section, we introduce orthogonal poly-

nomials for superior results.

The mapping functions U and V may be rewritten in terms of orthogonal polynomi-

als as

K

u = aiPi(x,y) (3.6.11)

v =  biPi(x,y )

where a i and bi are the unknown coefficients for orthogonal polynomials Pi. As we shall

see, introducing orthogonal polynomials allows us to determine the coefficients without

solving a linear system of equations. We begin by defining the orthogonality property

and then demonstrate how it is used to construct orthogonal polynomials from a set of

linearly independent basis functions. The orthogonal polynomials, together with the sup-

plied control points, are combined in a closed-form solution to yield the desired a i and b i

coefficients.

A set of polynomials P t (x,y), P 2 (x,y) ..... PK (x,y) is orthogonal over points (xt,,yt)

if

M

 Pi(x,Y)Pj(xk,yD = 0 i  j (3.6.12)

These polynomials may be constructed from a set of linearly independent basis functions

spanning the space of polynomials. We denote the basis functions here as h(x,y),

[l . I

h2(x,y), ..., hK(X,y). The polynomials can be constructed by using the Gram-Schmidt

orthogonalization process, a method which generates orthogonal functions by the follow-

ing incremental procedure.

Pl(x,y) = O:llhl(X,y )

P 2(x,y) = o;21P l(X,y) + q.22h2(x,Y)

P3 (x,y) = 31P l(x,Y) + ;32P2(x,Y) + ;33h3(x,Y)

PK(x,Y) = O;K1P I (x,y) + O;K2P 2(x,y) + ß ß ß + O;KKhK(X,y)

Basis functions hi(x,y) are free to be any linearly independent polynomials. For exam-

ple, the first six basis functions that we shall use are shown below.

hl(X,y) = 1

h2(x,y) = x

h3(x,y) = Y

hn(x,y) -- x 2

hs(x,y) = xy

h6(x,y) = y2

The o;ij parameters are determined by setting Oil = I and applying the orthogonaiity

property of Eq. (3.6.12) to the polynomials. That is, the o;ij's are selected such that the

product Pi Pj = 0 for i :j. The following results are obtained.

M

, [ P I(xl,YI:) ] 2

 Pj(x,,y,O hi(xl,Yl)

 [ ?j(x,y) l 2

k=l

once Pi has been multiplied with P 1 and P j, respectively. This multiplication serves to

eliminate all product terms except those associated with the desired coefficient.

Having computed the o:ij's, the orthogonal polynomials are determined. Note that

they are simply linear combinations of the basis functions. We must now solve for the al

3.6 POLY NOMIA L TRANSFORMATIONS 69

coefficients in Eq. (3.6.11). Using the least-squares approach, the error E is defined as

E =  aiPi(xl,,ylO-ul (3.6.14)

The coefficients are determined by taking the partial derivatives of E with respect to the

coefficients and setting them equal to zero. This results in the following system of linear

equations.

k=l

M M

k=l k=l

M

k=l

k=l

M

 vPi(xI½,YI½)

bi = M (3.6.17b)

 [Pi(xI,YI)] 2

Performing least-squares with orthogonal polynomials offers several advantages

worth noting. First, determining coefficients ai and bi in Eq. (3.6.11) does not require

solving a system of equations. Instead, a closed-form solution is available, as in Eq.

(3.6.17). This proves to be a faster and more accurate solution. The computational cost

of this method is O (MK). Second, additional orthogonal polynomial terms may be

added to the mapping function to increase the fitting accuracy of the approximation. For

instance, we may define the mean-square error to be

Em = '  aiP i(x,yk) - ul (3.6.18)

If Eros exceeds some threshold, then we may increase the number of orthogonal

polynomial terms in the mapping function to reduce the error. The orthogonality pro-

perty allows these terms to be added without recomputation of all the polynomial

coefficients. This facilitates a simple means of adaptively determining the degree of the

polynomial necessary to safisfy some error bound. Note that this is not true for ordinary

70 SPATIAL TRANSFORMATIONS

polynomials. In that case, as the degree of the polynomial increases the values of all

parameters change, requiring the recomputation of all the coefficients.

Numerical accuracy is generally enhanced by normalizing the data before perform-

ing the least-square computation. Thus, it is best to translate and scale the data to fit the

[-1,1] range. This serves to reduce the ill-conditioning of the problem. In this manner,

all the results of the basis function evaluations fall within a narrow range that exploits the

numerical properties of floating point computation.

3.6.5. Weighted Least-Squares

One flaw with the least-squares formulation is its global error measure. Nowhere in

Eqs. (3.6.7) or (3.6.14) is there any consideration given to the distance between control

points (xk,Y,O and approximation positions (x,y). Intuitively, it is desirable for the error

contributions of distant control points to have less influence than those which are closer

to the approximating position. This serves to confine local geometr/c differences, and

prevents them from averaging equally over the entire image.

The least-squares method may be localized by introducing a weighting function Wk

that represents the contribution of control point (xk,Yk) on point (x,y).

w = 1 (3.6.19)

/8 + (x-xD  + (y -yk)

The parameter 8 determines the influence of distant control points on approximating

points. As  approaches zero, the distant points have less influence and the approximat-

ing mapping function is made to follow the local control points more closely. As g

grows, it dominates the distance measurement and curtails the impact of the weighting

function to discriminate between local and distant control points. This serves to smooth

the approximated mapping function, making it approach the results of the standard least-

squares technique discussed earlier.

The weighted least-squares expression is given as

M

E (x,y ) =  [U (xt,yk) - u,t] 2 Wt(x,y ) (3.6.20)

the error term now becomes a function of position. Notice that each (x,y) point must

now recompute the squared-error summation. For ordinary polynomials, Eq. (3.6.10)

becomes

for l = 0, ..., N and rn = 0, ..., N - I. For orthogonal polynomials, the orthogonality pro-

perty of Eq. (3.6.12) becomes

M

 W(x,y) Pi(xl,Y,t) Pj(x,,yk) = 0 i j (3.6.22)

3.6 POLYNOMIAL TRA NSFOR MATIONS 71

and the parameters of the orthogonal polynomials are

M

 W(x,y) [  (x,y) 12

 Wk(x,y) P t (x,y0 hi(x

k=l

M

 Wt,(x,y) Pj(xk,YD hi(xby&)

 Wk(x,y) [ i(x,yO 12

k=l

Finally, the desired coefficients for Eq.

i= 1,...,K (3.6.23a)

i= I::j::KK;l(3.6.23b)

(3.6.11) as determined by the weighted least-

M

 W(x,y) u e(x,y0

M (3.6.24a)

 wk(x,y) [ e(x,y0 ] 

k=l

 W(x,y) v ?i(x,yO

k=l

bi(x,y) = M (3.6.24b)

k=l

greater than that of standard least-squares, where N is the number of pixels in the refer-

ence image. Although this technique is costlier, it is warranted when the mapping func-

tion is to be approximated using information highly in favor of local measurements.

Source code for the weighted least-squares method is provided below. The pro-

gram, written in C, is based on the excellent exposition of local approximation methods

found in [Goshtasby 88]. It expects the user to pass the list of correspondence points in

three global arrays: X, Y, and Z. That is, (X [i ],Y [i ]) --> Z [i ]. The full image of interpo-

lated or approximated correspondence values are stored into global array S.

........ 71 [111] 1l ' ii

72 SPATIAL TRANSFORMATIONS

#define MXTERMS 10

? global variables */

int N;

float init_alpha(), PO[Y0, coef0, basis();

Weighted least-squares with orthogonal polynomials

Input: X, Y, and Z are 3 float arrays for x, y, z=f(x,y) coords

delta is the smoothing factor (0 is not smooth)

Output: S <- fitted surface values of points (xsz by ysz)

(X, Y, Z, and N are passed as global arguments)

Based on algodthro described by Ardesh[r Goshtasby in

"Image registration by local approximation methods",

Image and Vision Computing, vol. 6, no. 4, Nov. 1988

wlsq(delta, xsz, ysz, S)

float delta, *S;

int xsz, ysz;

int i, j, k, x, y, t, terms;

float a, c, f, p, dx2, dy2, *s;

? N is already initialized with the nurober of control points '/

W = (float *) calloc(N, sizeof(float)); ? allocate memory for weights '/

? determine the number of terms necessary for error < .5 (optional) */

for(terms=3; terms < MXTERMS; terms++) {

for(i=0; i

/* init W: the weights of the N control points on x,y '/

for(j=0; j

dx2 = (X[i]-X[i]) ' (X[i]-X[j]);

dy2 = (Y[i]-Y[j]) * (YIi]-Y[j]);

} W[j] . 1.0 / sqrt(dx2 + dy2 + delta);

? init A; alak oeffs of the ortho polynomials */

for(j=0; j

for(j=0; j

for(k=0; k

}

for(t=f=0; t

a = coef(t);

p = poly(t, X[i], Y[i]);

f +=(a'p);

}

3,6 POLY NOMIA L TRANSFORMATIONS 73

if(i == N) break; ? found terms such that error < .5 */

}

/* perform sudace approximation */

for(x=0; x

/* init W: the weights of the N control points on x,y */

for(i=0; i

dx2 = (x-X[i]) * (x-X[i]);

Directory: filedownload

Download 2.54 Mb.

Share with your friends: