Digital image warping

the image proportions are altered resulting in a differentially scaled image.

Ix, y, 1] = [u, v, 1] Sv (3.3.5)

0

3.3.4. Shear

a22. We now consider the case where a 11 =a22 = 1, and a 12 =0. By allowing a21 to be

nonzero, x is made linearly dependent on both u and v, while y remains identical to v. A

similar operation can be applied along the v-axis to compute new values for y while x

remains unaffected. This effect, called shear, is therefore produced by using the off-

diagonal terms. The shear transform along the u-axis is

[x, y, 1] = [u, v, 1] 1 (3.3.6a)

0

where Hv is used to make x linearly dependent on v as well as u. Similarly, the shear

Ix, y, 1] = [u, v, 1] 1 (3.3.6b)

0

0 SPATIAL TRANSFOR MATIONS

Multiple wansforms can be collapsed into a single composite transformation. The

transforms are combined by taking the product of the 3 x 3 matrices. This is generally

not a commutative operation. An example of a composite transformation representing a

translation followed by a rotation and a scale change is given below.

[x,y, 1] = [u, v, 1] Mcon (3.3.7)

where

= COS0 Sv

TuTv 0 0

SucoS0 Svsin0

= -Susin0 Svcos0

Su(TucosO- TvsinO) Sv(TusinO+ TvcosO)

3.3.6. Inverse

The inverse of an affine transformation is itself affine. It can be readily computed

from the adjoint adj(T1) and determinant det(Tl) of transformation matrix T 1. From

linear algebra, we know that T/q = adj (T 1 ) / det (T 1 ) where the adjoint of a matrix is

simply the transpose of the matrix of cofactors [Strang 80]. This yields

[all a12 i]

[u, v, 1] = Ix, y, 1] lA21 A22 (3.3.8)

LA31 A32

[x, y, 1] 1 ' a22 -a12 0

= --a21 all 0

3.3.7. Inferring Affine Transformations

An affine transformation has six degrees of freedom, relating directly to coefficients

all, a21 , a31, a12 , a22 , and a32. In computer graphics, these coefficients are known by

virtue of the applied coordinate transformation. In areas such as remote sensing, how-

ever, it is usually of interest to infer the mapping given only a reference image and an

observed image. If an affine mapping is deemed adequate to describe the transformation,

the six coefficients may be derived by specifying the coordinate correspondence of three

3.3 AFFINE TRANSFORMATIONS 51

noncollinear points in both images. Let (uk,vk) and (xt,yk) for k =0,1,2 be these three

points in the reference and observed images, respectively. Equation (3.3.9) expresses

their relationship in the form of a matrix equation. The six unknown coefficients of the

affine mapping are determined by solving the system of six linear equations contained in

Eq. (3.3.9).

x! y = ul v a21 a22 (3.3.9)

X2 Y2 u2 v2 a31 a32

Let the system of equations given above be denoted as X = UA. In order to deter-

mine the coefficients, we isolate A by multiplying both sides with U -t , the inverse of the

matrix containing points (uk,v). As before, U -1 = adj (U) / der (U) where adj (U) is the

adjoint of U and der(U) is the determinant. Although the adjoint is always computable,

an inverse will not exist unless the determinant is nonzero. Fortunately, the constraint on

U to consist of noncollinear points serves to ensure that U is nonsingular, i.e., det (U) ; O.

Consequently, the inverse U -1 is guaranteed to exist. Solving for the coefficients in

terms of the known (u,v 0 and (x,yD pairs, we have

A = U-iX (3.3.10)'

or equivalently,

[all a12 i] v1--¾2 12--v 0 VO--Vl][XoYo!l

I

a21 a22 = de' tt2-ttl tt-tt2 ttl-tto Xl Yl

where

When more than three correspondence points are available, and when these points

are known to contain errors, it is common practice to approximate the coefficients by

solving an overdetermined system of equations. In that case, matrix U is no longer a

square 3 x 3 matrix and it must be inverted using any technique that solves a least-

squares linear system problem.

Since only three points are needed to infer an affine mapping, it is clear that affine

transformations realize a limited set of planar mappings. Essentially, they can map an

input triangle into an arbitrary triangle at the output. An input rectangle can be mapped

into a parallelogram at the output. More general distortions, however, cannot be handled

by affine transformations: For example, to map a rectangle into an arbitrary quadrilateral

requires a perspective, bilinear, or more complex transformation. Fast incremental

methods for computing aftinc mappings are discussed in Chapter 7.

SPATIAL TRANSFORMATIONS

3.4. PERSPECTIVE TRANSFORMATIONS

The general representation of aperspective transformation is

I alla12a13]

[x',y',w'] = [U,V,W] a21 a22 a23 (3.4.1)

a31 a32 a33

where x = x' /w' and y = y' /w'.

Aperspective transformation, orprojective mapping, is produced when [a13 a23 ]T

is nonzero. It is used in conjunction with a projection onto a viewing plane in what is

known as a perspective or central projection. Perspective transformations preserve

parallel lines only when they are parallel to the projection plane. Otherwise, lines con-

verge to a vanishing point. This has the property of fomshortening distant lines, a useful

technique for rendering realistic images. For perspective transformations, the forward

mapping functions are

X' alltt +a21v +a31

X = -- = (3.4.2a)

w' a13tt + a23v + a33

y = Y' = al2u+a22v+a32 (3.4.2b)

w' a13tt + a23v + a33

They take advantage of the fact that w' is allowed to vary at each point and division by

w' is equivalent to a projection using rays passing through the origin. Note that affine

xansformafions are a special case of perspective transformations where w' is constant

over the entire image, i.e., a 13 = a23 = 0.

Perspective Wansformations share several important properties with affine transfor-

mations. They are planar mappings, and thus their forward and inverse wansforms are

single-valued. They preserve lines in all orientations. That is, lines map onto lines

(although not of the same orientation). As we shall see, this desirable property is lacking

in more general mappings. Further/note, the eight degrees of freedom in Eq. (3.4.1) is

sufficient to permit planar quadrilateral-to-quadrilateral mappings. In contrast, affine

transformations offer only six degrees of freedom (Eq. 3.3.1) and thereby facilitate only

triangle-to-triangle mappings.

Examples of projective warps are shown in Fig. 3.6. Note that the intersections

along the edges are no longer equispaced. Also, in the rightmost image the horizontal

lines remain parallel because they lie parallel to the projection plane.

3.4.1, Inverse

The inverse of a projective mapping can be easily computed in terms of the adjoint

of the Ixansformation matrix T 1. Thus, Ti -1 =adj(Tl)/det(T1) where ad(Tl) is the

adjoint of T t and det(Tl) is the determinant. Since two matrices which are (nonzero)

scalar multiples of each other are equivalent in the homogeneous coordinate system,

3A PERSPECTIVE TRANSFORMATIONS 53

Figure 3.6: Perspective warps.

there is no need to divide by the determinant (a scalar). Consequently, the adjoint matrix

can be used in place of the inverse matrix. This proves to be a very useful result, espe-

cially since the adjoint will be well-behaved even if the determinant is very small when

the matrix is nearly singular. Note that if the matrix is singular, the inverse is undefined

and therefore the adjoint cannot be a scalar multiple of it. Due to these results from

linear algebra, the inverse is expressed below in terms of the elements in T 1 that are used

to realize the forward mapping.

[All A12 A13]

[u,v,w] = [x',y',w'l IA21 A22 A23 [ (3.4.3)

[A3 A32

= [x',y',w'] a23a31-a21a33 alla33-a13a31 a13a21_alla23

a21a32-a22a31 a 12a31 -alia32 alla22-a12a21

3.4.2. Inferring Perspective Transformations

A perspective transformation is expressed in trms of the nine coefficients in the

general 3 x 3 matrix T 1. Without loss of generality, T1 can be normalized so that

a33 = 1. This leaves eight degrees of freedom for a projective mapping. The eight

coefficients can be determined by establishing correspondence between four points in the

reference and observed images. Let (uk,vk) and (xk,yk) for k=0,1,2,3 be these four

points in the reference and observed images, respectively. Assuming a33 = 1, Eqs.

(3.4.2a) and (3.4.2b) can be rewritten as

4 SPATIALTRANSFORMATIONS

X = allu +a21v +a31 -a 13/.x -a23vx (3.4.4a)

y = a 12u + a 22v + a 32 - a 13/ly -- a23vy (3.4.4b)

Applying Eqs. (3.4.4a) and (3.4.4b) to the four pairs of correspondence points yields the

8 x 8 system of equations shown in Eq. (3.4.5).

-Uo v 0 1 0 0 0 -UoX 0 -VoX O-

ui vt 1 0 00,-utxt -vix

u2 !' 2 1 0 0 0 -u2x 2 -P2x2

/3 V3 1 0 0 0--/3X3 --V3X3

0 0 0 UO VO 1 --UoYo --VoYo A = X (3.4.5)

0 0 0 Ut V 1 -uy t -vy

0 0 0 u2 v2 1 -u2y2 -v2y2

0 0 0 u3 v3 1 -u3y 3 -v3y 3

where A = [all a21 a31 a12 a22 a32 a13 a23 iT are the unknown coefficients, and

X = [Xo x x2 x3 Yo Yt Y2 Y3 ]r are the known coordinates in the observed image.

The coefficients are determined by solving the linear system. This yields a solution

to the general (planar) quadrilateral-to-quadrilateral problem. Speedups are possible

when considering several special eases: square-to-quadrilateral, quadrilateral-to-square,

and quadrilateral-to-quadrilateral using the results of the last two cases. We now con-

sider each case individually. A detailed exposition is found in [Heckbert 89].

3.4.2.1. Case 1: Square-to. Quadrilateral

Consider mapping a unit/_stare onto an arbitrary quadrilateral. The following

four-point correspondences are'eslished from the uv-plane onto the .xy-plane.

(0,0) --> (Xo,Yo)

(1,0) -> (x,yt)

(1,1) --> (x2,Y2)

(0,1) . (x3,Y3)

In this case, the eight equations become

3.4 PERSPECTIVE TRANSFORMATIONS 55

a31 = x0

all +6/31 --a13x 1 = X 1

all +a21 +a31 -a13x2-a23x2 = X2

a21 +a31 --a23x 3 = X 3

a32 = Y0

a12+a32-a13Yl = Yl

a12 +a22 +a32 --aDY2 -- a23Y2 = Y2

a22 +a32-a23Y3 = Y3

The solution can take two forms, depending on whether the mapping is affine or perspec-

tive. We define the following terms for our discussion.

Xl = Xl -x2 x2 = x3 -x2 x3 = xo -Xl +x2 -x3

Ayl = yl -y2 Ay2 = y3-y2 Ay3 = yo-yl +y2-y3

If hx3 = 0 and Ay 3 = 0, then the .xy quadrilateral is a parallelogram. This implies that the

mapping is affine. As a result, we obtain the following coefficients.

ß all = Xl-x0

a21 = X 2 -x 1

a31 = x0

a12 = Yl -Y0

a22 = Y2-Yl

a32 = YO

a13 = 0

a23=0

perspective transformation are

a13 = Ay 2 / Ay 2

a23 = Ay 3 Ay2

56 S PATIA L TRANSFORMAT IONS

all = Xl--xoq-a13xl

a21 = x3--xo+a23x3

a31 = x0

at2 = Yl -yo+a13Y!

a22 = Y3-YO+a23Y3

a32 = Y0

This proves to be faster than the direct solution with a linear system solver. The compu-

tation may be generalized to map arbitrary rectangles onto quadrilaterals by pre-

multiplying with a scale and Ixanslation matrix.

3.4.2.2. Case 2: Quadrilateral-to-Square

This case is the inverse of the mapping already considered. As discussed earlier, the

adjoint of a projective mapping can be used in place of the inverse. Thus, the simplest

solution is to compute the square-to-quadrilateral mapping coefficients described above

to find the inverse of the desired mapping, and then take its adjoint to compute the

quadrilateral-to-square mapping.

3.4.2.3. Case 3: Quadrilateral-to-Quadrilateral

The results of the last two cases may be cascaded to yield a fast solution to the gen-

eral quadrilaterai-to-quadrilateral mapping problem. Figure 3.7 depicts this formulation.

case 3 

Figure 3.7: Quadrilateral-to-quadrilaterai mapping [Heckbert 89].

The general quadrilateral-to-quadrilateral problem is also known as four-corner

,tapping. Perspective Ixansformations offer a planar solution to this problem. When the

quadrilaterals become nonplanar, however, more general solutions are necessary. Bil-

inear transformations are an example of the simplest mapping functions that address

four-comer mappings for nonplanar quadrilaterals.

BILINEAR TRANSFORMATIONS 57

3.5. BILINEAR TRANSFORMATIONS

The general representation of a bilinear transformation is

a2 b2

a0 b0

problem for nonplanar quadrilaterals. It is most commonly used in the forward mapping

formulation where rectangles are mapped onto nonplanar quadrilaterals. It is pervaaive

in remote sensing and medical imaging where a grid of markings on the sensor are

imaged and registered with their known positions for calibration purposes. It is also

cormnon in computer graphics where it plays a central role in forward mapping algo-

rithms for texture mapping.

Bilinear mappings preserve lines that are horizontal or vertical in the source image.

This follows from the bilinear interpolation used to realize the transformation. Thus,

points along horizontal and vertical lines in the source image (including borders) remain

equispaced. This is a property shared with affine transformations. However, lines not

oriented along these two directions (e.g., diagonals) are not preserved as lines. Instead,

diagonal lines map onto quadratic curves at the output. Examples of bilinear warps are

shown in Fig. 3.8.

Figure 3.8: Bilinear warps.

Bilinear mappings are defined through piecewise functions that must interpolate the

coordinate assignments specified at the vertices. This scheme is based on bilinear inter-

polation to evaluate the X and Y mapping functions. We illustrate this method below for

computing X (u,v). An identical procedure is performed to compute Y (u,v).

58 S PATIA L TRA NSFOR MATIONS

3.5.1. Bilinear Interpolation

Bilinear interpolation utilizes a linear combination of the four "closest" pixel

values to produce a new, interpolated value. Given four points, (Uo,Vo), (Ul,VO,

(u2,v2), and (u3,v3), and their respective function values x0, Xl, x2. and x3, any inter-

mediate coordinate X (u,v) may be computed by the expression

X (u,v) = ao + a u + a2v + a3uv (3.5.2)

where the ai coefficients are obtained by solving

Xl = ul v UlVl at (3.5.3)

X2 /12 V2 /2V2 a2

X3 /13 V 3 /3V3 a3

Since the four points are assumed to lie on a rectangular grid, we rewrite them in the

above matrix in terms of u0, u, v 0, and v2. Namely, the points are (u0,v0), (Ul,V0),

(u0,v2), and (u 1 ,v2), respectively. Solving for ai and substituting into FA t. (3.5.2) yields

X(u',v') = Xo+(xl-xo)u' +(x2-xo)v'+(xs-x:-xl+xo)u' v' (3.5.4)

where/1' and v'  (0,1) are normalized coordinates that span the rectangle, and

u = Uo + (u t-uo) u'

V = Vo+(Vl-VO)V 

Therefore, given coordinates (u,v) and function values (xo,xl,x2,x3), the normalized

coordinates (u',v') are computed and the point correspondence (x,y) in the arbitrary qua-

drilateral is determined by Eq. (3.5.4). Figure 3.9 depicts this bilinear interpolation for

the X mapping function.

Figure 3.9: Bilinear interpolation.

3.5 BILINEAR TRANSFORMATIONS 59

3.5.2, Separability

The bilinear mapping is a separable transformation: it is the product of two 1-D

mappings, each operating along orthogonal axes. This property enables us to easily

extend 1-D linear interpolation into two dimensions, resulting in a computationally

efficient algorithm. The algorithm requires two passes, with the first pass applying 1-D

linear interpolation along the horizontal direction, and the second pass interpolating

along the vertical direction. For example, consider the rectangle shown in Fig. 3.10.

Points xol and x23 are interpolated in the first pass. These results are then used in the

second pass to compute the final value x.

0 u' 1 u'

1

Figure 3.10: Separable bilinear interpolation.

with the solution given in Eq. (3.5.4). In the first (horizontal) pass, we compute

Xo = Xo + (xt-xo) u' (3.5.5)

X23 = X2 + (X3--X2)t'

These two intermediate results are then combined in the second (vertical) pass to yield

the final value

X = X01 + (X23--X01) ;' (3.5.6)

= X 0 -I- (X l--X0) U' + [ (X2--X0) -I- (X3--X2--X l+X0)/1' ] ¾t

= X0 + (X l--X0) U' + (X2--X0) ¾' + (X3--X2--X l+X0) U' 1'

Notice that this result is identical with the classic solution derived in Eq. (3.5.4).

60 S PATIA L TRANSFOR MATIONS

3.5.3. Inverse

In remote sensing, the opposite problem is posed: given a normalized coordinate

(x',y') in an arbitrary (distorted) quadrilateral, find its position in the rectangle. Two

solutions are presented below.

By inverting Eq. (3.5.2), we can determine the normalized coordinate (u',v')

corresponding to the given coordinate (x,y). The derivation is given below. First, we

rewrite the expressions for x and y in terms of u and v, as given in Eq. (3.5.2).

x = ao +alu +a2 v +asuv (3.5.7a)

y = bo + blU + b2v + bsuv (3.5.7b)

Isolating u in Eq. (3.5.7a) gives us

x - a0 - a2v

u - -- (3.5.8)

In order to solve this, we must determine v. This can be done by substituting FA t. (3.5.8)

into Eq. (3.5.7b). Multiplying both sides by (al + asv) yields

y(al+a3v) = bo(al+asv) + bl(x-ao-a2v) + b2v(al+asv) + b3v(x-ao-a2v) (3.5.9)

This can be rewritten as

C2V2+ClV+Co = 0 (3.5.10)

where

Co = al (bo - y) + bl (x - ao)

C 2 = a3b2-a2b3

The inverse mapping for v,thasrequires the solution of a quadratic equation. Once v is

determined, it is plugged intq. (3.5.8) to compute u. Clearly, the inverse transform is

multi-valued and is more difficult to compute than the forward transform.

3.5.4. Interpolation Grid

Mapping from an arbitrary grid to a rectangular grid is an important step in per-

forming any 2-D interpolation within an arbitrary quadrilateral. The procedure is given

as follows.

1. To any point (x,y) inside an interpolation region defined by four arbfixary points, a

normalized coordinate (u',v') is associated in a rectangular region. This makes use

of the results derived above. A geometric interpretation is given in Fig. 3.11, where

the normalized coordinates can be found by determining the grid lines that intersect

at (x,y) (point P). Given the positions labeled at the vertices, the normalized coor-

dinates (u',v') are given as

s.s BILINEAR TRANSFORMATIONS 61

PolP0 PP2

P1Po P3P2

v' = P02P0 P13P1

P2Po PsP1

(3.5.11)

2. The function values at the four quadrilateral vertices are assigned to the rectangle

vertices.

3. A rectangular grid interpolation is then performed, using the normalized coordinates

to index the interpolation function.

4. The result is then assigned to point (x,y) in the distorted plane.

P3

P02 Pl3

P

It is important to note that the primary benefit of this procedure is that higher-order

interpolation methods (e.g., spline interpolation) that are commonly defined to operate on

rectangular lattices can now be extended into the domain of non-rectangular grids. This

thereby allows the generation of a continuous interpolation function for any arbitrary grid

[Bizais 83]. More will be said about this in Chapter 7, when we discuss separable mesh

warping,

3.6. POLYNOMIAL TRANSFORMATIONS

Geometric correction requires a spatial transformation to invert an unknown distor-

tion function. The mapping functions, U and V, have been almost universally chosen to

be global bivariate polynomial transformations of the form

U =   aijxly j (3.6.1)

i-0 j=0

i---o j=o

62 s PATIA L TRANSFORMATION S

where aij and blj are the constant polynomial coefficients. Since this formulation for

geometric correction originated in remote sensing [Markarian 71 ], the discussion below

will center on its use in that field. All the examples, though, have direct analogs in other

related areas such as medical imaging [Singh 79] and computer vision [Rosenfeld 82].

The polynomial Ixansformations given above are low-order global mapping func-

tions operating on the entire image. They are intended to account for sensor-related spa-

tial distortions such as centering, scale, skew, and pincushion effects, as well as errors

due to earth curvature, viewing geometry, and camera attitude and altitude deviations.

Due to dynamic operating conditions, these errors are comprised of internal and external

components. The internal errors are sensor-related distortions. External errors are due to

platform perturbations and scene characteristics. The effects of these errors have been

categorized in [Bernstein 71] and are shown in Fig. 3.12.

......  Scan Radially Tangentially

Centering Size Skew Nonlinearity Symmetric

Symmetric

TYPICAL SENSOR INTERNAL DISTORTIONS

Aspect Angle Distortion Scale Distortion Terrain Relief

(Attitude Effects) (Altitude Effect)

TYPICXTERNAL IMAGE DISTORTIONS

I 

I I

Figure 3.12: Common geometric image distortions.

These errors are characterized as low-frequency (smoothly varying) distortions.