Digital image warping

Download 2.54 Mb.
Size2.54 Mb.
1   ...   22   23   24   25   26   27   28   29   30

ple in which the traditional 2-pass method would clearly fail since a wide range of rota-

tion angles are represented. A vortex warp is shown in Fig. 7.44.



(a) (b)

(c) (d)

Figure 7.39: (a)I; (b)lx; (c)S; (d)XLUTand YLUT.


(a) (b)

(c) (d)

Figure 7.40: Undersampled spatial lookup tables. (a) Ixy; (b) lr; (c) Undersampled

LUTs; (d) Output.


(a) (b)

Figure 7,41: (a) Foldover; (b)XLUT and YLUT.

(a) (b)

Figure 7.42: Magnified foldover. (a) No filtering. (b) Filtered result.

(a) (b)

Figure 7.43: Bending rows. (a) Checkerboard; (b) Madonna.

(a) (b)

Figure 7.44: Vortex warp. (a) Checkerboard; (b) Madonna.



Scanline algorithms all share a common theme: simple interpolation, antialiasing,

and data access are made possible when operating along a single dimension. Using a 2-

pass transform as an example, the first pass represent a forward mapping. Since the data

is assumed to be unidirectional, a single-element accumulator is sufficient for filtering

purposes. This is in contrast to a full 2-D accumulator array for standard forward map-

pings. The second pass is actually a hybrid mapping function, requiring an inverse map-

ping to allow a new forward mapping to proceed. Namely, auxiliary function H must be

solved before G, the second-pass forward mapping, can be evaluated.

A benefit of this approach is that clipping along one dimension is possible. For

instance, there is no need to compute H for a particular column that is known in advance

to be clipped. This results in some timesavings. Te principal difficulty, however, is the

bottleneck problem which exists as a form of aliasing. Tis is avoided in some applica-

tions, such as rotation, where it has been shown that no scaling is necessary in any of the

1-D passes. More generally, special attention must be provided to counteract this degra-

dation. This has been demonstrated for the case of arbitrary spatial lookup tables.



Digital image warping is a subject of widespread interest. It is of practical impor-

tance to the remote sensing, medical imaging, computer vision, and computer graphics

communities. pical applications can be grouped into two classes: geometric correction

and geometric distortion. Geometric correction refers to distortion compensation of

imaging sensors, decalibralion, and geometric normalization. This is applied to remote

sensing, medical imaging, and computer vision. Geometric distortion refers to texture

mapping, a powerful computer graphics tool for realistic image synthesis.

All geometric transformations have three principal components: spatial transforma-

tion, image resampling, and antiallasing. They have each received considerable atten-

tion. However, due to domain-dependent assumptions and constraints, they have rarely

received uniform treatment. For instance, in remote sensing work where there is usually

no severe scale change, image reconstrantion is more sophisticated than antialiasing.

However, in computer graphics where there is often more dramatic image compression,

antialiasing plays a more significant role. This has served to obscure the single underly-

ing set of principles that govern all geometric transformations for digital images. The

goal of this book has been to survey the numerous contributions to this field, with special

emphasis given to the presentation of a single coherent framework.

Various formulations of spatial transformations have been reviewed, including

affine and perspective mappings, polynomial transformations, piecewise polynomial

transformations, and four-comer mapping. The role of these mapping functions in

geometric correction and geometric distortion was discussed. For instance, polynomial

transformations were introduced to extend the class of mappings beyond affine transfor-

mations. Thus, in addition to performing the common translate, scale, rotate, and shear

operations, it is possible to invert pincushion and barrel distortions. For more local con-

trol, piecewise polynomial transformations are widespread. It was shown that by estab-

lishing several correspondence points, an entire mapping function can be generated

through the use of local interpolants. This is actually a surface reconstruction problem.

There continues to be a great deal of activity in this area as evidenced by recent papers

on multigd relaxation algorithms to iterafively propagate constraints throughout the



surface. Consequently, the tools of this field of mathematics can be applied direcdy to

spatial transformations.

Image resampling has been shown to primarily consist of image reconstruction, an

interpolation process. Various interpolation methods have been reviewed, including the

(truncated) sinc function, nearest neighbor, linear interpolation, cubic convolution, 2-

parameter cubic filters, and cubic splines. By analyzing the responses of their filter ker-

nels in the frequenay domain, a comparison of interpolation methods was presented. In

particular, the quality of interpolation is assessed by examining the performance of the

interpolation kernel in the passbands and stopbands. A review of sampling theory has

been included to provide the necessary background for a comprehensive understanding of

image resampling and anfiaiiasing.

Antiaiiasing has recently attracted much attention in the computer graphics com-

munity. The earliest antiaiiasing algorithms were restrictive in terms of the preimage

shape and filter kernel that they supported. For example, box filtering over rectangular

preimages were common. Later developments obtained major performance gains by

retaining these restrictions but permitting the number of computations to be independent

of the preimage area. Subsequent improvements offered fewer restrictions at lower cost.

In these instances the preimage areas were extended to ellipses and the filter kernels, now

stored in lookup tables, were allowed to be arbitrary. The design of efficient filters that

operate over an arbitrary input area and accommodam arbitrary filter kernels remains a

great challenge.

Development of superior filters used another line of attack: advanced sampling stra-

tegies. They include supersampling, adaptive sampling, and stochastic sampling. These

techniques draw upon recent results on perception and the human visual system. The

suggested sampling patterns that are derived from the blue noise criteria offer promising

results. Their critics, however, point to the excessive sampling densities required to

reduce noise levels to unobjecfionable limits. Determining minimum sampling densities

which satisfy some subjective criteria requires addifionai work.

The final section has discussed various separable aigorithrns introduced to obtain

large performance gains. These algorithms have been shown to apply over a wide range

of transformations, including perspective projection of rectangles, bivariate patches, and

superquadrics. Hardware products, such as the Ampex ADO and Quantel Mirage, are

based on these techniques to produce real-time video effects for the television industry.

Recent progress has been made in scanline algorithms that avoid the botfieoeok problem,

a degradation that is particular to the separable method. These modifications have been

demonstrated on the speciai case of rotation and the arbitrary case of spatial lookup


Despite the relatively short history of geometric transformation techniques for digi-

tai images, a great deal of progress has been made. This has been accelerated within the

last decade through the proliferation of fast and cost-effective digital hardware. Algo-

rithms which were too costly to consider in the early development of this area, are either

commonplace or am receiving increased attention. Future work in the areas of recon-

strantion and antiaiiasing will most likely integrate models of the human visual system to


achieve higher quality images. This has been demonstrated in a recent study of a family

of filters defined by piecewise cubic polynomials, as well as recent work in stochastic

sampling. Related problems that deserve attention include new adaptive filtering tech-

niques, irregular sampling algorithms, and reconstruction from irregular samples. In

addition, work remains to be done on efficient separable schemes to integrate sophisti-

cated reconstraction and antiaiiasing filters into a system supporting more general spatial

transformations. This is likely to have great impact on the various diverse communities

which have contributed to this broad area.

Appendix 1


The purpose of this appendix is to provide a detailed review of the Fast Fourier

Transform (FFT). Some familiarity with the basic concepts of the Fourier Transform is

assumed. The review begins with a definition of the discrete Fourier Transform (DFT) in

Section AI.1. Directly evaluating the DFr is demonstrated there to be an O(N 2) pro-


The efficient approach for evaluating the DFr is through the use of FFT algorithms.

Their existence became generally known in the mid-1960s, stemming from the work of

J.W. Cooley and J.W. Tukey. Although they pioneered new FFT algorithms, the original

work was actually discovered over 20 years earlier by Danielson and Lanczos. Their for-

mulation, known as the Danielson-Lanczos Lemma, is derived in Section A1.2. Their

recursire solution is shown to reduce the computational complexity to O (N log2 N).

A modification of that method, die Cooley-Tukey algorithm, is given in Section

A1.3. Yet another variation, the Cooley-Sande algorithm, is described in Section A1.4.

These last two techniques are also known in the literature as the decimation-in-time and

decimation-in-frequency algorithms, respectively. Finally, source code, written in C, is

provided to supplement the discussion.




Consider an input function f (x) sampled at N regularly spaced intervals. This

yields a list of numbers, fk, where 0 -

taken to be complex numbers, i.e., having real and imaginary components. The DFT off

is defined as

F n = flce -i2nnlN O


__ x-, F i2nnktN 0-< n fn = N ,_]=_O , e -

where i =x/. Equations ( and ( define the forward and inverse DFTs,

respectively. Since both DITFs share the same cost of computation, we shall confine our

discussion to the forward DFT and shall refer to it only as the DFT.

The DFT serves to map the N input samples of f into the N frequency terms in F.

From Eq. (, we see that each of the N frequency terms are computed by a linear

combination of the N input samples. Therefore, the total computation requires N 2 com-

plex multiplications and N(N-1) complex additions. The straightforward computation of

the DFT thereby give rise to an O (N 2) process. This can be seen more readily if we

rewrite Eq. ( as

?n = f,w n O



W = e -i2nlN = cos(-2/N) + isin(-2;/N) (Al.l.3)

For reasons described later, we assume that

N =2 r

where r is a positive integer. That is, N is a power of 2.

Equation (A 1.i.2) casts the DFT as a matrix multiplication between the input vector

f and the two-dimensional array composed of powers of W. The entries in the 2-D array,

indexed by n and k, represent the N equally spaced values along a sinusoid at e/ch of the

N frequencies. Since straightforward matrix multiplication is an O(N 2) process, the

computational complexity of the DFT is bounded from above by this limit.

In the next section, we show how the DFT may be computed in O (N log2 N) opera-

tions with the FFT, as originally derived over forty years ago. By properly decomposing

Eq. (, the reduction in proportionality from N 2 to N log2 N multiply/add opera-

tions represents a significant saving in computation effort, particularly when N is large.



In 1942, Danielson and Lanczos derived a recursive solution for the DFr. They

showed that a DFT of length N t can be rewritten as the sum of two DFTs, each of length

N/2, where N is an integer power of 2. The first DFr makes use of the even-numbered

points of the original N; the second uses the odd-numbered points. The following proof

is offered.

Fn =  f,e -12nnk/N (A1.2.1)

(N/2)-I (N/2)-I

 f2,te-12n(2k)"N +  f2,+le -i2nn(2'+l)"N (AI.2.2)

k=0 k=0

(N/2)-I (N/2)-I

=  f2ke -12zCn/(N"2) + W n  f2k+l e-12r/(NI2) (A1.2.3)

k=0 k=0


n o

= F, + W F n

Equation (A1.2.1) restates the original definition of the DFT. The summation is

expressed in Eq. (A1.2.2) as two smaller summations consisting of the even- and odd-

numbered terms, respectively. To properly access the data, the index is changed from k

to 2k and 2k+l, and the upper limit becomes (N/2)-l. These changes to the indexing

variable and its upper limit give rise to Eq. (A1.2.3), where both sums are expressed in a

form equivalent to a DFr of length N/2. The notation is simplified further in Eq.

(A1.2.4). There, Fe,, denotes the n tn component of the Fourier Transform of length N/2

formed from the even components of the original f, while Fn is the corresponding

transform derived from the odd components.

The expression given in Eq. (A1.2.4) is the central idea of the Danielson-Lanczos

Lemma and the decimation-in-time FFT algorithm described later. It presents a divide-

and-conquer solution to the problem. In this manner, solving a problem (Fn) is reduced

to solving two smaller subproblems (Fg and Fn% However, a closer look at the two

sums, Fen and Fn , illustrates a potentially troublesome deviation from the original

definition of the DFT: N/2 points of fare used to generate N points. (Recall that n in F,

and Fn is still made to vary from 0 to N-1 ). Since each of the subproblems appears to be

no smaller than the original problem, this would thereby seem to be a wasteful app:oach.

Fortunately, there exists symmetries which we exploit to reduce the computational com-


The first simplification is found by observing that F n is periodic in the length of'the

transform. That is, given a DFT of length N, Fn+N = Fn. The proof is given below.

This is also known as an N-point DFT.



Fn+N =  fk e-i2g(n+N)k/N


=  fk e-i2nnklN e-i2klN


=  fk e-i2nnklN


= Fn


In the second line of Eq. (A1.2.5), the last exponential term drops out because the

exponent -i2Nk/N is simply an integer multiple of 2 and e -i2nk = 1. Relating this

result to Eq. (AI.2.4), we note that Fe and F have period N/2. Thus,

Fne+Ni2 = F 0 -< n < N/2 (A1.2.6)

Fn+Ni2 = Fn 0 -< n < N/2

This permits the N/2 values of F and F to trivially generate the N numbers needed for


A similar simplification exists for the W n factor in Eq. (A1.2.4). Since W has

period N, the first N/2 values can be used to trivially generate the remaining N/2 values

by the following relation.

cos((2/N)(n-}-N/2)) = -cos(2n/N) 0 -< n < N/2 (A1.2.7)

sin((2/N)(n+N/2)) = -sin(2n/N) 0 < n < N/2


W nqv/2 = -W n 0 -< n < N/2 (AI.2.8)

Summarizing the above results, we have

Fn = Fne +W"F O_

Fn+N/2 = Fn e -WnF 0 -

where N is an integer power of 2.


A1.2.1. Butterfly Flow Graph

Equation (A1.2.9) can be represented by the butterfly flow graph of Fig.,

where the minus sign in _+W n arises in the computation of Fn+N/2. The terms along the

branches represent multiplicative factors applied to the input nodes. The intersecting

node denotes a summation. For convenience, this flow graph is represented by the

simplified diagram of Fig. A 1. lb. Note that a butterfly performs only one complex multi-

plication (WnFn). This product is used in Eq. (A1.2.9) to yield Fn and FnV/2.

F Fn+Ni2 Fn  Fn+N/2

(a) (b)

Figure AI.I: (a) Butterfly ttow graph; (b) Simplified diagram.

The expansion of a butterfly ttow graph in terms of the computed real and imaginary

terms is given below. For notational convenience, the real and imaginary components of

a complex number are denoted by the subscripts r and i, respectively. We define the fol-

lowing variables.

g = F,

h = Fn

Wr = cos(-2n/N)

w i = sin(-2n/N)

Expanding Fn, we have

F n = g +Wnh (A1.2.10)

= [gr + igi] + [Wr q- iWi] [hr + ihi]

= [gr + igi] + [wrhr - wihi + iwrhi + iwihr]

= [gr + Wrhr -- wihi] + i [gi + wrhi + wihr]

The real and imaginary components of Wnh are thus wrh r -wih i and wrh i + wihr,

respectively. These terms are isolated in the computation so that they may be subtracted

from gr and gi to yield Fn+N/2 without any additional transform evaluations.


A1.2.2. Putting It All Together

The recursive formulation of the Danielson-Lanczos Lemma is demonstrated in the

following example. Consider list fof 8 complex numbers labeled f0 through f7 in Fig.

A1.2. In order to reassign the list entries with the Fourier coefficients Fn, we must evalu-

ate F and F,. As a result, two new lists are created containing the even and odd com-

ponents of f. The e and o labels along the branches denote the path of even and odd

components, respectively. Applying the same procedure to the newly created lists, suc-

cessive halving is performed until the lowest level is reached, leaving only one element

per list. The result of this recursive subdivision is shown in Fig. A1.2.

fo fl f2 f3 f4 f5 f6 f7

eee eeo eoe eoo oee oeo ooe ooo

Figure A1.2: Recursive subdivision into even- and odd-indexed lists.

At this point, we may begin working our way back up the tree, building up the

coefficients by using the Danielson-Lanczos Lemma given in Eq. (A1.2.9). Figure A1.3

depicts this process by using butterfly flow graphs to specify the necessary complex addi-

tions and multiplications. Note that bold lines are used to delimit lists in the figure.

Beginning with the 1-element lists, the 1-point DFTs are evaluated first. Since a 1-point

DFT is simply an identity operation that copies its one input number into its one output

slot, the 1-element lists remain the same.

The 2-point transforms now make use of the 1-point transform results. Next, the 4-

point transforms build upon the 2-point results. In this case, N is 4 and the exponent of

W is made to vary from 0 to (N/2)-l, or 1. In Fig. A1.3, all butterfly flow graphs assume


an N of 8 for the W factor. Therefore, the listed numbers are normalized accordingly.

For the 4-point transform, the exponents of 0 and 1 (assuming an N of 4) become 0 and 2

to compensate for the implied N value of 8. Finally, the last step is the evaluation of an

8-point transform. In general, we combine adjacent pairs of 1-point transforms to get 2-

point transforms, then combine adjacent pairs of pairs to get 4-point transforms, and so

on, until the first and second halves of the whole data set are combined into the final


0 ,i 0 l10 l10

Figure A1.3: Application of the Danielson-Lanczos Lemma.


A1.2.3. Recursive FFT Algorithm

The Danielson-Lanczos Lemma provides an easily programable method for the

DFT computation. It is encapsulated in Eq. (A1.2.9) and presented in the FFT procedure

given below.

Procedure FFT(N,f)

1. If N equals 2, then do


2. Replace f0 by f0 +fl andfl by f0 -fl.

3. Return


4. Else do:


5. Define g as a list consisting of all points of f which have an even index

and h as a list containing the remaining odd points.

6. Call FFT(N/2, g)

7. Call FFT(N/2, h)

8. Replace fn by gn + Wnh for n=0 to N-1.



The above procedure is invoked with two arguments: N and f. N is the number of

points being passed in array f. As long as N is greater than 2, f is split into two halves g

and h. Array g stores those points of f having an even index, while h stores the odd-

indexed points. The Fourier Transforms of these two lists are then computed by invoking

the FFT procedure on g and h with length N/2. The FFT program will overwrite the con-

tents of the lists with their DFT results. They are then combined in line 8 according to

Eq. (A1.2.4).

The successive halving proceeds until N is equal to 2. At that point, as observed in

Fig. A1.3, the exponent of W is fixed at 0. Since W is 1, there is no need to perform the

multiplication and the results may be determined directly (line 2).

Returning to line 8, the timesavings there arises from using the N/2 available ele-

ments in g and h to generate the N numbers required. This is a realization of Eq.

(A1.2.9), with the real and imaginary terms given in Eq. (A1.2.10). The following seg-

ment of C code implements line 8 in the above algorithm. Note that all variables, except

N, N 2, and n, are of type double.


ang = 0;

inc = -6.2831853 / N;

N2 =N/2;

1or(n=0; n

Wr = cos(ang);

wl: $in(ang);

ang += inc;

a = wr*hr[n] - wi*hi[n];

fr[nl = grin] + a:

fr[n+N2] = gr[n] - a;

a = wi*hr[n] + wr*hi[n];

fi[nl = gi[n] + a:

fi[n+N2] = gi[n] - a;

/* initialize angle '/

/* angle increment: 2/N '/

/* real part of W n '1

/* imaginary part of W n '1

/* next angle in W n '1

/* real part of Wnh (Eq. A1.2.1 0) */

/* Danielson-Lanczos Lemma (Eq. A1.2.9) */

/* imaginary part of Wnh (Eq. A1.2.10) '/

/* Danielson-Lanczos Lemma (Eq. A1.2.9) '/

Directory: filedownload

Download 2.54 Mb.

Share with your friends:
1   ...   22   23   24   25   26   27   28   29   30

The database is protected by copyright © 2024
send message

    Main page