Digital image warping
Download 2.54 Mb.
|
|
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. 256 SCANLINE ALGORITHMS (a) (b) (c) (d)
Figure 7.39: (a)I; (b)lx; (c)S; (d)XLUTand YLUT. 7.7 $EPAP, ABL IMAGE WARPING 257 (a) (b) (c) (d)
Figure 7.40: Undersampled spatial lookup tables. (a) Ixy; (b) lr; (c) Undersampled LUTs; (d) Output. 258 SCANLINE ALGORITHMS (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. 260 SCANLINE ALGORITHMS 7.8. DISCUSSION 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. 8 EPILOGUE 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 261
262 EPILOGUE 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 tables.
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 263 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 FAST FOURIER TRANSFORMS 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- cess.
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. 265
266 ALl, DISCRETE FOURIER TRANSFORM 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 k-O
where i =x/. Equations (Al.l.la) and (Al.l.lb) 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. (Al.l.la), 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. (Al.l.la) as
?n = f,w n O k=0
where
For reasons described later, we assume that N =2 r
where r is a positive integer. That is, N is a power of 2. 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. (Al.l.la), 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. 267
A1.2. DANIELSON-LANCZOS LEMMA 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 (A1.2.4) 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- plexity. 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. .I 268
Fn+N = fk e-i2g(n+N)k/N k-O
= fk e-i2nnklN e-i2klN k=0
= fk e-i2nnklN k=0
= Fn (A1.2.5) 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 Fn. 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 Therefore, 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. 269
A1.2.1. Butterfly Flow Graph 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)
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,
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.
270
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 271
an N of 8 for the W factor. Therefore, the listed numbers are normalized accordingly. 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 transform. 0 ,i 0 l10 l10 Figure A1.3: Application of the Danielson-Lanczos Lemma. 272 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 Begin
2. Replace f0 by f0 +fl andfl by f0 -fl. 3. Return End 4. Else do: Begin 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. End End
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. 273
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) '/ |