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
where M>6. A similar equation holds for v and bij. Both of these expressions may be
written in matrix notation as
U = WA (3.6.3)
V=WB
In order to solve for A and B, we must compute the inverse of W. However, since W
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)
Expanding the matrix notation, we obtain Eq. (3.6.5), the following system of linear
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
66 SPATIAL TRA NSFORMAT IONS
M
M
= [ U (xk,Y,0- uk 12 (3.6.7)
k=l
M
= [aoo+aoxk+aotYk+anxy+aeox+ao2Y-Uk] 2
k=l
This is achiev by deteining the ptial derivatives of E with respect to cfficients
alj, d uatg them to zero. For each cfficietu aij, we have
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)
K
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)
k=l
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
68 SPAT IAL TRANSFORMATIONS
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
O;ii = -- M k=l i = 2,...,K (3.6.13a)
Pj(x,,y,O hi(xl,Yl)
[ ?j(x,y) l 2
k=l
Equations (3.6.13a) and (3.6.13b) are the results of isolating the o:ij coefficients
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
Applying the orthogonal property of Eq. (3.6.12), we obtain the following'simplification
M M
ai [ Pi(x,Yk) 12 = uPi(x,yk) (3.6.16)
k=l k=l
The desired coefficients are thereby given as
M
k=l
ai M (3.6.17a)
k=l
A similar procedure is repeated for computing bi, yielding
M
vPi(xI½,YI½)
k=l
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)
This represents a dramatic incre omputation over the nonweighted method since
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
aij t(x,y)x,y;[ IY,t = W, t u,txty (3.6.21)
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)
k=l
3.6 POLYNOMIAL TRA NSFOR MATIONS 71
and the parameters of the orthogonal polynomials are
M
W(x,y) [ (x,y) 12
k=l
Wk(x,y) P t (x,y0 hi(x
k=l
M
Wt,(x,y) Pj(xk,YD hi(xby&)
k=l
Wk(x,y) [ i(x,yO 12
k=l
Finally, the desired coefficients for Eq.
squares method are
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
ai(x,y ) =
M (3.6.24a)
wk(x,y) [ e(x,y0 ]
k=l
M
W(x,y) v ?i(x,yO
k=l
bi(x,y) = M (3.6.24b)
Wk(x,y) [ Pi(xk,y D ]2
k=l
The computational complexity of weighted least-squares is O (NMK3). It is N times
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 *X, *Y, *Z, *W, AiM XTERM S][MXTERMS];
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(ABS(Z[i] -f) > .5) break;
if(i == N) break; ? found terms such that error < .5 */
}
/* perform sudace approximation */
for(y=0; y
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]);
**Share with your friends:** |