Digital image warping

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)

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

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,

Feeee'"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

ooo

OlO

lOO

11o

Original Array

f0

fl

f3

f5

f7

000

010

001

011

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 frm 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

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

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

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

k=O

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

f[ = f + fl*+N2 (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

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

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 [0] ... buf [15]), 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;

double FCTR, fctr, a, b, c, s, num[2];

qdP qa, qb;

cpqdinfo(ql, q2);

N = ql->width;

rl = (float *) ql->buf[0];

il = (float *) ql->buf[1];

r2 = (float *) q2->buf[0];

12 = (float *) q2->buf[1];

if(N == 2) {

a = rl[0] + r111];

b: i1101 + i111];

r211] = rl[0]- r111];

i211] = i110]- i111];

r210] = a;

i210] = b;

} else {

N2=N/2;

qb = getqd(1, N2, FFT_TYPE);

ra = (float *) qa->buf[0]; ia = (float *) qa->buf[1];

ca = (float *) qb->buf[0]; ib = (float *) qb->buf[1];

/* 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 */

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[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);

if(dir) { /* inverse: divide by log N */

num[0] = num[1] = 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[2];

qdP qsrc;

if(q1 == q2) {

qsm = NEWQD;

cpqd(ql, qsrc);

} else qsm =ql;

cpqdinfo(ql, q2);

rl = (float *) qsrc->buf[0];

il = (float *) qsrc->but[1];

r2 = (float *) q2->buf[0];

i2 = (float *) q2->buf[1];

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;

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[0] = num[1] = 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(xt+l) = 3A3zx& 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-yAx,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 -- yA 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 AXiy+  _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 Ak- 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 xj

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(rh2+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-iAx2-! +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: filedownload

Download 2.54 Mb.

Share with your friends: