# Digital image warping

Download 2.54 Mb.
 Page 27/30 Date 18.10.2016 Size 2.54 Mb.
 AI.2.4. Cost of Computation The Danielson-Lanczos Lemma, as give 0 in Eq. (A1.2.9), can be used to calculate the cost of the computation. Let C (N) be the cost for evaluating the tzansform of N points. Combining the transforms of N points in Eq. (A1.2.9) requires effort proportional to N because of the multiplication of the terms by W n and the subsequent addition. If c is a constant reflecting the cost of such operations, then we have the following result for C(N). C(N) = 2C()+cN (A1.2.11) This yields a recurrence relation that is known to result into an O(NlogN) process. Viewed another way, since there are 1og2N levels to the recursion and cost 0 (N) at each level, the total cost is 0 (N log2 N). 274 AL3. COOLEY-TUKEY ALGORITHM The Danielson-Lanczos Lemma presented a recursive solution to computing the Fourier Transform. The role of the recursion is to subdivide the original input into smaller lists that are eventually combined according to the lemma. The starting point of the computation thus begins with the adjacent pairing of 1-point DFTs. In the preceding discussion, their order was determined by the recurslye subdivision. An alternate method is available to determine their order directly, without the need for the recursive algorithm given above. This result is known as the Cooley-Tukey, or decimation-in-time algorithm. To describe the method, we define the following notation. Let F ee be the list of even-indexed terms taken from F e. Similarly, F e is the list of odd-indexed terms taken fromF e. In general, the suing of symbols in the superscript specifies the path traversed in the tree representing the recursive subdivision of the input data (Fig. A 1.2). Note that the height of the tree is log2 N and that all leaves denote 1 -point DFTs that are actually elements from the input numbers. Thus, for every pattern of e's and o's, numbering log2 N in all, Feeee'"ee = fn for some n (A1.3.1) The problem is now to directly find which value of n corresponds to which pattern ofe's and o's in Eq. (A1.3.1). The solution is surprisingly simple: reverse the pattern of e's and o's, then let e = 0 and o = 1, and the resulting binary string denotes the value of n. This works because the successive subdivisions of the data into even and odd are tests of successive low-order (least significant) bits of n. Examining Fig. A1.2, we observe that traversing successive levels of the tree along the e and o branches corresponds to successively scanning the binary value of index n from the least significant to the most significant bit. The strings appearing under the bottom row designates the traversed path. The procedure for N = 8 is summarized in Table AI.1. There we see the binary indices listed next to the corresponding array elements. The first subdivision of the data into even- and odd-indexed elements amounts to testing the let(st significant (rightmost) bit. If that bit is 0, an even index is implied; a 1 bit designates an odd index. Subsequent subdivisions apply the same bit tests to successive index bits of higher significance. Observe that in Fig. A1.2, even-indexed lists move down the left branches of the tree. Therefore, the order in which the leaves appear from left to fight indicate the sequence of ls and Os seen in the index while scanning in reverse order, from least to most significant bits. 275 Original Index ooo OOl OlO Oll lOO lol 11o 111 Original Array f0 fl f2 f3 f4 f5 f6 f7 Bit-reversed Index 000 100 010 110 001 101 011 ill Reordered Array Table AI.I: Bit-reversal and array reordering for input into FFT algorithm. The distinction between the Cooley-Tukey algorithm and the Danielson-Lanczos Lemma is subtle. In the latter, a recursire procedure is introduced in which to compute the DFT. This procedure is responsible for decimating the input signal into a sequence that is then combined, during the traversal back up the tree, to yield the 'ansform output. In the Cooley-Tukey algorithm, though, the recursion is unnecessary since a clever bit- reversal trick is introduced to achieve the same disordered input. Furthermore, directly reordering the input in this way simplifies the bookkeeping necessary in recombining terms. Source code for the Cooley-Tukey FFT algorithm, written in C, is provided in Section A1.5. A1.3.1. Computational Cost The computation effort for evaluating the FFT is easily determined frm this formu- lation. First, we observe that there are log2 N levels of recursion necessary in computing Fn. At each level, there are N/2 butterflies to compute the F, and Fn terms (see Fig. A1.3). Since each butterfly requires one complex multiplication and two complex addi- tions, the total number of multiplications and additions is (N/2) log2 N and N log2 N, respectively. This 0 (N log2 N) process represents a considerable savings in computation over the 0 (N 2) approach of direct evaluation. For example if N > 512, the number of multiplications is reduced to a fraction of 1 percent of that required by direct evaluation. 276 A1.4. COOLEY-SANDE ALGORITHM In the Cooley-Tukey algorithm, the given data sequence is reordered according to a bit-reversal scheme before it is recombined to yield the transform output. The reordering is a consequence of the Danielson-Lanczos Lemma that calls for a recursive subdivision into a sequence of even- and odd-indexed elements. The Cooley-Sande FFT algorithm, also known as the decimation-in-frequency algo- rithm, calls for recursively splitting the given sequence about its midpoint, N/2. Fn =  fk e-i2v'nklN (A1.4.1) k=0 (N/2)-I N-1  fke -i2v'nk/N +  fke -i2nnk/N k=0 k=NI2 (N/2)-t (N/2)-i  fke-i2nnk/N+  fk+Ni2 e-i2nn(k+N/2)lN k=0 k=0 (N/2)-I = k__O [fk + fk+N'2e-in] e-i2't'/ Noticing that the e -nin factor reduces to +1 and -1 for even and odd values of n, respec- tively, we isolate the even and odd terms by changing n to 2n and 2n+l. (N/2)-, [ F2n =  fk + fk+N! e-i2u(2n)klN 0 --< n < NI2 (A1.4.2) k=O (N/2)-I c k=O (N/2)-i r F2+l = Z [fk -f,t+V/2J e -12n(2n+l)'/v 0 < n < N/2 (Ai.4.3) k=0 (N/2)- 1 r 1 Thus, the even- and odd-indexed values of F are given by the DFTs of f and fff where f[ = f + fl*+N2 (A1.4.4) f,= [f,--f/+N/2] W  (A1.4.5) The same procedure can now be applied to f[ and fl. This sequence is depicted in Fig. A1.4. The top row represents input list fcontaining 8 elements. Again, note that lists am delimited by bold lines. Regarding the butterfly notation, the lower left branches denote Eq. (A1.4.4) and the lower right branches denote Eq. (A1.4.5). 277 Since all the even-indexed values of F need f,, a new list is created for that pur- pose. This is shown as the left list of the second row. Similarly, the ff list is generated, appearing as the second list on that row. Of course, the list sizes diminish by a factor of two with each level since generating them makes use of f, and fk+v/2 to yield one ele- ment in the new list. This process of computing Eels. (A1.4.4) and (A1.4.5) to generate new lists terminates when N = i, leaving us F, the transform output, in the last row. In contrast to the decimation-in-time  algorithm, in which the input is disor- dered but the output is ordered, the opposite is true of the decimation-in-frequency FFT algorithm. However, reordering can be easily accomplished by reversing the binary representation of the location index at the end of computation. The advantage of this algorithm is that the values of f are entered in the input array sequentially. f0 fl f2 f3 f4 f5 f6 f7 0 I I 0 I I 0 I I i I Figure A1.4: Decimation-in-frequency FFT algorithm. 278 A1.5. SOURCE CODE This section provides sottree code for the recursive FFT procedure given in Section A1.2, as well as code for the Cooley-Tukey algorithm described in Section A1.3. The programs e written in C and make use of some library routines described below. The data is passed to the functions in quads. A quad is an image contool block, con- raiding information about the image. Such data includes the image dimensions (height and width), pointers to the uninterleaved image channels (buf  ... buf ), and other necessary information. Since the complex numbers have real and imaginary com- ponents, they occupy 2 channels in the input and output quads (channels 0 and 1). A brief description of the library routines included in the listing is given below. 1) cpqd (q 1,q 2) simply copies quad q 1 into q 2. 2) cpqdinfo (q 1,q 2) copies the header information of q 1 into q 2. 3) NEWQD allocates a quad header. The image memory is allocated later when the dimensions are known. 4) getqd(h,w, type) returns a quad containing sufficient memory for an image with dimensions h xw and channel datatypes type. Note that FFT_TYPE is defined as 2 channels of type float. 5) freeqd (q) frees quad q, leaving it available for any subsequent getqd call. 6) divconst (q 1,num, q 2) divides the data in q 1 by num and puts the result in q 2. Note that hum is an array of numbers used to divide the corresponding channels in q 1. 7) Finally, PI2 is defined to be 2, where  = 3.141592653589793. AI.5.1. Recursive FFT Algorithm fftl D(ql ,dir,q2) /* Fast Fourlet Transform (1 D) */ int dir; /* dir=0: forward; dir=l:inverse */ qdP ql, q2; /* ql =input; q2=output */ ( int i, N, N2; float *rl, *il, *r2, *i2, *ra, *ia, *ca, *lb; double FCTR, fctr, a, b, c, s, num; qdP qa, qb; cpqdinfo(ql, q2); N = ql->width; rl = (float *) ql->buf; il = (float *) ql->buf; r2 = (float *) q2->buf; 12 = (float *) q2->buf; if(N == 2) { a = rl + r111]; b: i1101 + i111]; r211] = rl- r111]; i211] = i110]- i111]; r210] = a; i210] = b; } else { N2=N/2; qa = getqd(1, N2, FFT_TYPE); qb = getqd(1, N2, FFT_TYPE); ra = (float *) qa->buf; ia = (float *) qa->buf; ca = (float *) qb->buf; ib = (float *) qb->buf; /* split list into 2 halves: even and odd */ for(l=0; i ra[i]= *rl++; ia[i] = *i1++; Ca[i] = *rl++; ib[i] = *i1++; } /* compute fit on both lists */ ftl D(qa, dir, qa); fftlD(qb, dir, qb); /* build up coefficients */ if(!dir) /* forward */ FCTR = -PI2 / N; else FCTR = PI2/N; for(fctr=i=0; I c: cos(fctr); s = sin(fctr); a = c*rb[i] - s*ib[i]; /* F(0)=f(0)+f(1); F(1)=f(0)-f(1) */ /* a,b needed when rl=r2 */ 279 280 r2[il = ra[il + a; r2[i+N21: ra[il - a; a = s*rb[i] + c*ib[i]; i2[i] = ia[i] + a; i2[i+N2] = ia[i] - a; } freeqd(qa); freeqd(qb); if(dir) { /* inverse: divide by log N */ num = num = 2; divconst(q2, num, q2); A1.5.2. Cooley-Tukey FFT Algorithm fit1D(q 1, dir, q2) /* Fast Fourier Transform (1 D) */ int dir; /* dir=l: forward; dir= -1: inverse */ qdP ql, q2; /* Uses bit reversalto avoid recursion */ { /* and trig recurrence for sin and cos */ int i, j, IogN, N, N1, NN, NN2, itr, offst; unsigned int a, b, msb; float *rl, *r2, *il, *i2; double wr, wi, wpr, wpi, wtemp, theta, tempr, tempi, num; qdP qsrc; if(q1 == q2) { qsm = NEWQD; cpqd(ql, qsrc); } else qsm =ql; cpqdinfo(ql, q2); rl = (float *) qsrc->buf; il = (float *) qsrc->but; r2 = (float *) q2->buf; i2 = (float *) q2->buf; N = ql->width; N1 =N-l; for(IogN=0,i=N/2; i; IogN++,i/=2); /* # of bits sig digits in N */ msb = LSB << (IogN-1); for(i=1; i a= 1; b=O; for(l=O; a && j if(a & LSB) b I= (msb>>j); a>>= 1; } /* swap complex numbers: [i] <--> [b] */ r2[i] = rl[b]; i2[i] = il[b]; r2[b] = rl[i]; i2[b] = il[i]; } /* copy elements 0 and N1 since they don't swap */ r210] = rl[O]; i210] = i110]; r2[N1] = rl[N1]; i2[N1] = il[Nl]; /* NN denotes the number of points in the transform. It grows by a power of 2 with each iteration. NN2 denotes NN/2 which is used to trivially generate NN points from NN2 complex numbers. Computation of the sines and cosines of multiple angles is made through recurrence relations. wr is the cosine for the real terms; wi is sine for 281 282 the Imaginary termS. '/ NN=I; for(itr=0; ttr NN2 = NN; NN <<=1; /*NN*=2*/ theta = -PI2 / NN * dir; wtemp = sin(.5*theta); wpr = -2 * wtemp * wtemp; wpi = sin(theta); wi =04 for(offst=0; offstfor(i=offst; i j=i+NN2; tempr = wr*r2[j] - wi*i2[j]; tempi = wi*r2[I] + wr*i2[]; r2[j] = r2[i] - tempr; r2[i] = r2111 +ternpr; i2[j] = i2[i] - ternpi; i2[i] = i2[i] + tempi; /*trigonometric recurrence */ wr = (wtemp=wr)*wpr -wi'wpi + wr; wi = wi*wpr + wtemp*wpi + wi; } ) if(dir == -1) { /* inverse transform: divide by N '/ ' num = num = N; divconst(q2, hum, q2); If(qsrc I= ql) lreeqd(qsrc); Appendix 2 INTERPOLATING CUBIC SPLINES The purpose of this appendix is to review the fundamentals of interpolating cubic splines. We begin by defining a cubic spline in Section A2.1. Since we are dealing with interpolating splines, constraints are imposed to guarantee that the spline actually passes through the given data points. These constraints are described in Section A2.2. They establish a relationship between the known data points and the unknown coefficients used to completely specify the spline. Due to extra degrees of freedom, the coefficients may be solved in terms of the first or second derivatives. Both derivations are given in Sec- tion A2.3. Once the coefficients are expressed in terms of either the first or second derivatives, these unknown derivatives must be determined. Their solution, using one of several end conditions, is given in Section A2.4. Finally source code, written in C, is provided in Section A2.5 to implement cubic spline interpolation for uniformly and nonuniformly spaced data points. A2.1. DEFINITION A cubic spline f (x) interpolating on the partition x0 < x < .. ß < xn-1 is a func- tion for which f (xt,)=y,¾ It is a piecewise polynomial function that consists of n-i cubic polynomials ft* defined on the ranges [xk,xk+l]. Furthermore, ft, are joined at x (k=i,...,n-2) such that f and f[' are continuous. An example of a cubic spline pass- ing through n data points is illustrated in Fig. A2.1. The k t polynomial piece, f, is defined over the fixed interval [x,xk+ ] and has the cubic form f(x) = A3(x -x/) 3 +A2(x -x,0 2 +A l(X -x/) +A0 (A2.1.1) 283 284 f (x) f5 fo ft f X0 X 1 X2 X3 X4 X5 X6 Figure A2.1: Cubic spline. A2.2. CONSTRAINTS Given only the data points (x,,y,t), we must determine the polynomial coefficients, A, for each partition such that the resulting polynomials pass through the data points and are continuous in their first and second derivatives. These conditions require ft, to satisfy the following constraints y, = f,t(x,) = A 0 (A2.2.1) Yk+l = fJt(x,+O = A3A x + A2A x 2 + A A x,t + Ao y,[ = f(x&) = A i Y+I = J(xt+l) = 3A3zx& 2 + 2A2Ax. + A 1 (A2.2.2) y' = f'(xt,) = 2,42 (A2.2.3) Y\$/+t = '+t (xt,) = 6A3ax,+2A2 Note that these conditions apply at the data points (xk,Yk). If the xk's are defined on a regular grid, they are equally spaced and Axn =xe+i -x = 1. This eliminates all of the Ax,t terms in the above equations. Consequently, Eqs. (A2.2.1) through (A2.2.3) reduce to 285 yk = A o (A2.2.4) Yk+l = A3 +A2 +AI +A0 y = A  (A2.2.5) y+l = 3A3+2A2+A1 y' = 2A 2 (A2.2.6) Y+i = 6A3 + 2A2 In the remainder of this appendix, we will refrain from making any simplifying assump- tions about the spacing of the data points in order to treat the more general case. A2.3. SOLVING FOR THE SPLINE COEFFICIENTS The conditions given above are used to find A 3, A2, A 1, and A0 which are needed to define the cubic polynomial piece ft. Isolating the coefficients, we get A0 = y (A2.3.1) At = y 1 [3AYt-2y-y+,] A2 =  AXk i f 2AY ,+ , ] A3 = AX'- [-- x +Yk Y+iJ In the expressions for A 2 and A 3, k = 0,..., n-2 and Ay,t = y+l -Y,t. Y,t+l-A3Ax-yAx,t--y ' (A2.3.2a) A 2 = AXk 2 A2.3.L Derivation of A2 From (A2.2.1), From (A2.2.2), Y+I - 3A 3Ax - y 2A 2 - (A2.3.2b) Finally, A 2 is derived from (A2.3.2a) and (A2.3.2b) 286 A2.3.2. Derivation of A3 From (A2.2.1), From (A2.2.2), y+ - A A X 2 -- yA x -- ye A3 = AX (A2.3.2c) y+i - 2A 2/X x: - y 3A 3 = Xk2 (A2.3.2d) Finally, A 3 is derived from (A2.3.2c) and (A2.3.2d) Axk j The equations in (A2.3.1) express the coefficients of f/ in terms of xk, Yk, x/+, y+, (known) and y, y+ (unknown). Since the expressions in Eqs. (A2.2.1) through (A2.2.3) present six equations for the four A i coefficients, the A terms could alternately be expressed in terms of second derivatives, instead of the first derivatives given in Eq. (A2.3.1). This yields A0 = y (A2.3.3) a l A y Ax, [Y,+i + 2y') Ax 6 y' A2=- A3 = 6'-x A2.3.3. Derivation ofA 1 and A 3 From (A2.2.1), y+ - A 3Ax 3 - -Ax 2 - y A1 = Ax/c From (A2.2.3), y+ -y' A3 6Ax/ Plugging Eq. (A2.3.4b) into (A2.3.4a), a = Ay AXiy+  _y,,]-½Ax, Ax 6 (A2.3.4a) 6Axk (A2.3.4b) AYk s I y+l + 2y,J (A2.3.4c) Ax 287 A2.4, EVALUATING THE UNKNOWN DERIVATIVES Having expressed the cubic polynomial coefficients in terms of data points and derivatives, the unknown derivatives still remain to be determined. They are typically not given explicitly. Instead, we may evaluate them from the given constcaims. Although the spline coefficients require either the first derivatives or the second deriva- tives, we shall derive both forms for the sake of completenessß A2.4.1. First Derivatives We begin by deriving the expressions for the first derivatives using Eqs. (A2.2.1) through (A2.2.3). Recall that the A coefficients expressed in terms of y' made use of Eqs. (A2.2.1) and (A2.2.2). We therefore use the remaining constraint, given in Eq. (A2.2.3), to express the desired relation. Constcalnt Eq. (A2.2.3) defines the second derivative off at the endpoints of its interval. By establishing that f'[_ (x) = f'(x,0, we enfome the continuity of the second derivative across the intervals and give rise to a rela- tion for the first derivatives. 6A3-'Ax,_i +2A - = 2A2  (A2.4.1) Note that the superscripts refer to the interval of the coefficient. Plugging Eq. (A2.3.1) into Eq. (A2.4.1) yields A_i -12 Ayk-1 1 6 Ayk-1 --2y AXk_i Ak- 1 AXk_ 1  J 4 AXk_ 1 -- AXk_--  AXk After collecting the y' terms on one side, we have Eq. (A2.4.2): 1 +y,[2[ "'1 + L'I ] +y+!  [Ax_l AxkJ = 3/ Y- -1 L L x_l xj Equation (A2.4.2) yields a matrix of n-2 equations in n unknowns. We can reduee the need for division operations by multiplying both sides by AXe4 AXe. This gives us the following system of equations, with 1 h=Ax and rk=Ay/Ax. hi 2(ho+h ) ho h2 2(hi+h2) hl hn-2 2(hn-3+hn-2) hn-3 y Y [ 3(roh + rtho) Y ] 3(rh2+r2h)  _ [ 3 (r n -3 hn -2 - rn -2 hn -3 ) 288 When the two end tangent vectors y6 and y_ are specified, then the system of equations becomes determinable. One of several boundary conditions described later may be selected to yield the remaining two equations in the matrix. A2.4.2. Second Derivatives An alternate, but equivalent, course of action is to determine the spline coefficients by solving for the unknown second derivatives. This procedure is virtually identical to the approach given above. Note that while there is no particular benefit in using second derivatives rather than first derivatives, it is presented here for generality. As before, we note that the A coefficients expressed in terms of y" made use of Eqs. (A2.2.1) and (A2.2.3). We therefore use the remaining constraint, given in Eq. (A2.2.2), to express the desired relation. Constraint Eq. (A2.2.2) defines the first derivative offk at the endpoints of its interval. By establishing that -1 (xk) =f(xk) we enforce the con- tinuity of the first derivative across the intervals and give rise to a relation for the second derivatives. 3A-iAx2-! +2A-iAx,t-I +A -t = A1  (A2.4.3) Again, the superscripts refer to the interval of the coefficient. Plugging Eq. (A2.3.3) into Eq. (A2.4.3) yields AX_i 6 ' [ y'[' + 2y_ = AX 6 Y+i + 2y After collecting the y" terms on one side, we have Ax ax_l j Equation (A2.4.4) yields the following matrix of n-2 equations in n unknowns. Again, for notational convenience we let h,t =Ax,t and rk =Aye/Axe. Directory: filedownloadDownload 2.54 Mb.Share with your friends:

The database is protected by copyright ©ininet.org 2020
send message

Main page