Digital image warping

y

hn-s 2(hn-s+hn_2) hn-2

6(r 1 - r0)

6(r 2 - r 1 )

5(r._2 - r._3)

The system of equations becomes determinable once the boundary conditions are

specified.

289

A trivial choice for the boundary condition is achieved by setting y' =y'_ = O.

This is known as thefree-end condition that results in natural spline interpolation. Since

y' = 0, we know from Eq. (A2.2.6) that A2 = 0. As a result, we derive the following

expression from Eq. (A2.3.1).

Yl 3Ay0

Similarly, since Y-t = 0, 6A 3 + 2A2 = 0, and we derive the following expression

from Eq. (A2.3.1).

2y_2 + 4y;-1 = 6 'Aye-2 (A2.4.6)

AXn -2

Another condition is called the cyclic condition, where the derivatives at the end-

y) = y_ (A2.4.7)

y(; = y'

The boundary condition that we shall consider is the not-a-knot condition. This

requires y'" to be continuous across Xl and xn-2. In effect, this extrapolates the curve

from the adjacent interior segments [de Boor 78]. As a result, we get

A3 = A] (A2.4.8)

x2 [ yo , ,] 1 +y]

-2'x0 +Y+Yl j 'x2 [-Ayl '

1 = _2_Xl +y 1

Replacing y with an expression in terms of y and y allows us to remain consistent

with the stracture of a tridiagonal matrix already derived earlier. From Eq. (A2.4.2), we

isolate y« and get

lay0 Ayll ' AXl [ AXl+AX0] (A2.4.9)

y : 3AXlL+J -yo c

Substituting this expression into Eq. (A2.4.8) yields

yAxt[Axo+AXt I +Y [AJ0+AXl] 2= AY [3AXoAXl +2Ax121 + AAYx [Ax 1

Similarly, the last row is derived to be

2 , A

+ -- [ 3AXn-3A Xn-2 + 2 Xn-31

The version of this boundary condition expressed in terms of second derivatives is left to

290

the reader as an exercise.

Many simplifications are possible if we assume that the points are equispaced, i.e.,

Axk= 1. This is certainly the case for image reconstracfion, where cubic splines can be

used to compute image values between regularly spaced samples. The not-a-knot boun-

dary condition used in conjunction with the system of equations given in Eq. (A2.4.2) is

shown below. To solve for the polynomial coefficients, the column vector containing the

first derivatives must be solved and then substituted into Eq. (A2.3.1).

41

14

y

-Sy0 + 4Yl +Y2

3(y2 -Y0)

3(y3 -Yl)

3 (Vn - - y -3 )

--Yn-3 --4Yn-2 + 5yn-I

(A2.4.10)

A2.5. SOURCE CODE

Below we include two C programs for interpolating cubic splines. The first pro-

gram, called ispline, assumes that the supplied data points are equispaced. The second

program, ispline_gen, addresses the more general case of irregularly spaced data.

A2.5.1. Ispline

The function ispline takes Y 1, a list of len 1 numbers in double-precision, and

passes an interpolating cubic spline through that data. The spline is then resampled at

len2 equal intervals and stored in list Y2. It begins by computing the ulknQwn first

derivatives at each interval endpoint. It invokes the function getYD, which returns the

first derivatives in the list YD. Along the way, function tridiag is called to solve the tridi-

agohal system of equations shown in Eq. (A2.4.10). Since each derivative is coupled

only to its adjacent neighbors on both sides, the equations can be solved in linear time,

i.e., O(n). Once YD is initialized, it is used together with Y1 to compute the spline

coefficients. In the interest of speed, the cubic polynomials are evaluated by using

Horoer's rule for factoring polynomials. This requires three additions and three multipli-

cations per evaluated point.

Interpolating cubic spline function for equispaced points

Input: Y1 is a list of equispaced data points with lenl entries

Output: Y2 <- cubic spline sampled at len2 equispaced points

ispline(Y1 ,lenl ,Y2,1en2)

double*Y1, *Y2;

int lenl, len2;

int i, ip, oip;

double *YD, A0, A1, A2, A3, x, p, fctr;

/* compute 1st derivatives at each point -> YD */

YD = (double *) calloc(len1, sizeof(double));

getYD(Y1,YD,lenl);

/*

* p is real-valued position into spline

* oip is left endpoint of last visited interval

*/

oip = -1; /* force coefficient initialization '/

for(i=p=ip=0; i < len2; i++) {

/* check if in new interval */

if(ip != oip) {

/* update inte.'val */

oip = ip;

/* compute spline coefficients */

A0 = Yl[ip];

A1 = YD{ip];

A2 = 3.0*(Yl[ip+l]-Yl[ip]) - 2.0*YD[ip] - YD[ip+l];

A3 = -2.0*(Yl[ip+l]-Yl[ip]) + YD[ip] + YD[ip+l];

}

/* use Homer's rule to calculate cubic polynomial */

Y2[i] = ((A3*x + A2)*x + A1)*x + A0;

/* increment pointer */

ip = (p += fctr);

}

cfree((char *) YD);

291

YD <- Computed 1st derivative of data in Y (len enidos)

The not-a-knot boundary condition is used

getYD(Y,YD,len)

double *Y, *YD;

int en;

int i;

YD[0] = -5.0*Y[0] + 4.0'Y[1] + Y[2];

YD[len-1] = -Y[len-3] - 4.0*Y[len-2] + 5.0*Y[len-1];

/* solve for the tridiagonal matdx: YD=YD*inv(tridlag matrix) */

tridiag(YD,len);

Linear time Gauss Elimination with backsubstitution for 141

tddiagonal matrix with column vector D. Result goes into D

tridiag(D,len)

double *D;

int len;

int i;

? init first two entries; C is dght band of tridiagonal */

C = (double *) calloc(len, sizeof(double));

D[0] = 0.5*D[0];

Dill = (Dill - D[0]) / 2.0;

C[1] = 2.0;

/* Gauss elimination; forward substitution */

for(i = 2; i < len-1; i++) {

C[i] = 1.0 / (4.0 - C[i-1]);

D[i] = (Dill - D[i-1]) / (4.0- C[i]);

C[i] = 1.0/(4.0 - C[i-1]);

D[i]= (Dill - 4*D[i-1]) / (2.0- 4*C[i]);

/* backsubstitution */

for(i = len-2; i >= 0; i--) D[i] -= (D[i+l] * C[i+l]);

cfree((char *) C);

293

The function ispline_gen takes the data points in (X1,Y1), two lists of len 1

numbers, and passes an interpolating cubic spline through that data. The spline is then

resamplod at len 2 positions and stored in Y2. The resampling locations are given by X 2.

The function assumes that X2 is monotonically increasing and lies withing the range of

numbers in X 1.

As beforo, we begin by computing the unknown first derivativos at each interval

endpoint. The function getYD_gen is then invoked to return the first derivatives in the

list YD. Along the way, function tridiagen is called to solve the tridiagonal system of

equations given in Eq. (A2.4.2). Once YD is initialized, it is used together with Y 1 to

compute the spline coefficients. Note that in this general case, additional consideration

must now be givon to determine the polynomial interval in which the resampling point

lios.

Interpolating cubic spline function for irregularly-spaced points

Their x-coordinates are specified in Xl

Output: Y2 <- cubic spline sampled according to X2 (len2 entries)

Assume that Xl ,X2 entries are monotonically increasing

isplinegen(X 1 ,Y1 ,lenl ,X2,Y2,1en2)

double *Xl, *Y1, *X2, *Y2;

inf lenl, len2;

{

int i, j;

/* compute 1st derivatives at each point -> YD '/

YD = (double *) calloc(lenl, sizeof(double));

getYDgen(Xl ,Y1,YD,lenl);

/* error checking '/

if(X2[0] < Xl[0] II X2[len2-1] > Xl[lenl-1]) (

fprintf(stderr,"ispline_gen: Out of range0);

exit();

P

* pl is left endpoint of interval

* p3 is dght endpoint of interval

* j is input index for current interval

*/

for(i=j=0; i < len2; i++) {

/* check if in new interval "/

p2 = X2[i];

if(p2 > p3) {

/* find the interval which contains p2 */

for(; jXl[j]; J++);

if(p2 < Xl[I]) j--;

pl = XI[J]; /* update left endpoint */

p3 = XI[j+I]; /* update right endpoint */

/* compute spline coefficients */

dx = 1.0/(Xl[j+I] - Xl[j]);

dy = (Y1 [j+l] - YI[j]) * dx;

A0: YI[j];

A1 = YD[j];

A2 = dx ' (3.0*dy - 2,0*YD[j] - YD[j+I]);

A3 = dx*dx * (-2,0*dy + YD[j] + YD[j+I]);

}

/* use Homer's rule to calculate cubic polynomial '1

Y2[i] = ((A3*x + A2)*x + A1)*x + A0;

}

cfree((char *) ¾D);

The not-a-knot boundary condition is used

getYD_gen(X,Y,YD, len)

double *X, *Y, *YD;

int len;

{

int i;

/* allocate memory for tridiagonal bands A,B,C */

A = (double *) calloc(len, sizeof(double));

B = (double *) calloc(len, slzeof(double));

C = (double ') calloc(len, sizeof(double));

/* init first row data */

h0 = X[1]- X[0]; hl = X[2]- X[1];

r0 = (Y[1] - Y[0]) / h0; rl = (Y[2] - Y[1]) / hl;

B[0] = hl * (h0+hl);

C[0] = (h0+hl) * (h0+hl);

YD[0] = r0*(3*h0*hl + 2*h1*h1) + rl*h0*h0;

P init tddiagonal bands A, B, C, and column vector YD '/

/* YD will later be used to return the derivatives */

for(i = 1; i < len-1; i++) {

h0 = X[i]- X[i-1]; hl = X[i+l] - X[i];

r0: (Y[i] - Y[i-1]) / h0; rl = (Y[i+l] - Y[i]) / hl;

A[i] = hl;

B[i]= 2 * (h0+hl);

c[i] = h0;

YD[i] = 3 * (r0*hl + rl*h0);

I

/* last row */

B[i]= h0 * (h0+hl);

YD[i] = r0*hl*hl + rl*(3*h0*hl + 2*h0*h0);

/* solve for the tridiagonal matrix: YD=YD*inv(tridiag matrix) */

tridiag_ge n(A,B,C,Y D,le n);

cfree((char *) A);

cfree((char *) B);

dree((char *) C);

Gauss Elimination with backsubstitution for general

tridiagonal matrix with bands A,B,C and column vector D.

295

II

t ridiag_g en(A,B,C,D,len)

double *A, *B, *C, *D;

int len;

int i;

double b, *F;

/* Gauss elimination; forward substitution */

b = B[0];

D[O] = D[O] / b;

for(i = 1;i < len; i++) {

F[i] = C[i-1] / b;

b = B[i] - Alii*Eli];

if(b == 0) {

fpd ntf(stderr,"getY D_gen: divide-by-zero0);

exit();

D[i] = (D[i] - D[i-1]*A[i]) / b;

)

f' backsubstitution */

cfme((char *) F);

Appendix 3

FORWARD DIFFERENCE METHOD

The method of forward differences is used to simplify the computation of polynomi-

als. It basically extends the incremental evaluation, as used in scanline algorithms, to

higher-order functions, e.g., quadratic and cubic polynomials. We find use for this in

Chapter 7 where, for example, perspective mappings are approximated without costly

division operations. In this appendix, we derive the forward differences for quadratic and

cubic polynomials. The method is then generalized to polynomials of arbitrary degree.

The (first) forward difference of a function f (x) is defined as

Af(x) = f(x+)-f(x), >0 (A3.1)

It is used together with f (x), the value at the current position, to determine f (x + ), the

value at the next position. In our application, = 1, denoting a single pixel step size. We

shall assume this value for in the discussion that follows.

For simplicity, we begin by considering the forward difference method for linear

interpolation. In this case, the first forward difference is simply the slope of the line

passing through two supplied function values. That is, Af(x)=at for the function

f(x)=a]x+ao. We have already seen it used in Section 7.2 for Gouraud shading,

whereby the intensity value at position x+l along a scanline is computed by adding

Af (x) to f (x). Surprisingly, this approach readily lends itself to higher-order interpo-

lants. The only difference, however, is that Af (x) is itelf subject to update. That update

is driven by a second increment, known as the second forward difference. The extent to

which these increments are updated is based on the degree of the polynomial being

evaluated. In general, a polynomial of degree N requires N forward differences.

We now describe forward differencing for evaluating quadratic polyhomials of the

form

f(x) = a2x2 +alx +ao (A3.2)

297

Af(x) = f(x+l)--f(x) (A3.3)

= a2(2x+ 1)+at

Thus, Af (x) is a linear expression. If we apply forward differences to Af (x), we get

A2f (x) = A(Af (x)) (A3.4)

= Af(x+l)--Af(x)

= 2a2

Since A2f (x) is a constant, there is no need for further terms. The second forward differ-

added to the latest result to compute the new value. Each loop in the iteration can be

rewritten as

f(x+l) = f (x)+Af (x) (A3.5)

A/(x+l) = Af (x)+ A2f (x)

If computation begins at x = 0, then the basis for the iteration is given by f, A f, and A2f

evaluated at x=0. Given these three values, the second-degree polynomial can be

evaluated from 0 to lastx using the following C code.

for(x = 0; x < lastx; x++) {

f[x+l] = f[x] + Af; /* compute next point '/

Af += A2f; /* update 1st forward difference '/

}

Notice that Afis subject to update by A2f, but the latter term remains constant throughout

A similar derivation is given for cubic polynomials. However, an additional for-

ward difference constant must be incorporated due to the additional polynomial degree.

For a third-degree polynomial of the form

f(x) = a3x 3 +a2x 2 +alx +ao (A3.6)

the first forward difference is

299

Af (x) = f (x + 1) --f (x) (A3.7)

Since bf (x) is a second-degree polynomial, two more forward difference terms are

derived. They are

A2f(x) = z(Af(x)) = 6a3x+6a3 +2a2 (A3.8)

A3f(x) = A(A2f(x)) = 6a 3

The use of forward differences to evaluate a cubic polynomial between 0 and lastx is

demons'ated in the following C code.

tor(x = 0; x < lastx; x++} {

l[x+l] = f[2x] + zf; /* compute next point '/

Af += A f; /' update 1st forward difference*/

A2f += A3I; /* update 2nd forward difference */

}

In contrast to the earlier example, this case has an additional forward difference

a third-degree polynomial with only three additions per value. An alternate approach,

using Homer's role for factoring polynomials, requires three additions and three multipli-

cations. This makes forward differences the method of choice for evaluating polynomi-

als.

The forward difference approach for cubic polynomials is depicted in Fig. A3.1.

sion of this method in Chapter 7, texture coordinates are used as the function values.

Thus, we begin with u0, AU0, and A2U0 defined for position x0, where the subscripts

refer to the position along the scanline.

In order to compute our next texture coordinate at x = 1, we add AU 0 to u 0. This is

denoted by the arrows that are in contact with u0. Note that diagonal arrows denote

addition, while vertical arrows refer to the computed sum, Therefore, u 1 is the result of

adding AU 0 to u0. The following coordinate, u2, is the sum of ul and AU. The latter

term is derived from zu0 and A2U0. This regular structure collapses nicely into a com-

pact iteration, as demonstrated by the short programs given earlier.

Higher-order polynomials are handled by adding more forward difference terms.

This corresponds to augmenting Fig. A3.1 with additional columns to the right. The

order of computation is from the left to right. That is, the summations corresponding to

the diagonal arrows are executed beginning from the left coltann. This gives rise to the

300

adjacent elements directly below. Those elements are then combined in similar fashion.

points has been evaluated.

Position Value

XO UO

X2 U 2 AU! A2U0

x3 u3 au2 a2u a3uo

x4 u4 au3 2u2 3u

Figure A3.1: Forward difference method.

[Abdou 82]

[Abram 85]

[Akima 78]

[Akima 84]

[Anderson 90]

[Andrews 76]

[Antoniou 79]

[Atteia 66]

[Ballard 82]

[Barnhill 77]

REFERENCES

Abdou, Ikram E. and Kwan Y. Wong, "Analysis of Linear Interpolation

Schemes for Bi-Level Image Applications," IBM J. Res. Develop., vol.

26, no. 6, pp. 667-680, November 1982.

Abram, Greg, Lee Westover, and Turner Whirred, "Efficient Alias-free

Rendering Using Bit-Masks and Look-Up Tables," Computer Graphics,

(SIGGRAPH '85 Proceedings), vol. i9, no. 3, pp. 53-59, July 1985.

Aldma, H., "A Method of Bivariate Interpolation and Smooth Surface

Fitting for Irregularly Distributed Data Points," ACM Trans. Math.

Software, vol. 4, pp. 148-159, 1978.

Akima, H., "On Estimating Partial Derivatives for Bivariate Interpola-

tion of Scattered Data," Rocky Mountain J. Math., vol. 14, pp. 41-52,

1984.

Animation," ILM Technical Memo #1026, Computer Graphics Depart-

ment, Lucasfilm Ltd., 1990.

Andrews, Harry C. and Claude L. Patterson III, "Digital Interpolation of

Discrete Images," IEEE Trans. Computers, vol. C-25, pp. 196-202,

1977.

New York, 1979.

Atteia, M., "Existence et determination des foncfions splines a plusieurs

variables," C. R. Acad. Sci. Paris, vol. 262, pp. 575-578, 1966.

Ballard, Dana and Christopher Brown, Computer Vision, Prentice-Hall,

Englewood Cliffs, NJ, 1982.

Barnhill, Robert E., "Representation and Approximation of Surfaces,"

Mathematical Software III, Ed. by J.R. Rice, Academic Press, London,

pp. 69-120, 1977.

301

[Bennett 84a]

[Bennett 84b]

[Bergland 69]

[Bernstein 71]

[Bernstein 76]

[Bier 86]

[Bizais 83]

[Blake 87]

[Blinn 76]

[BooIt 86a]

[Boult 86b]

[Braccini 80]

[Bracewell 86]

[Briggs 74]

Bennett, Phillip P. and Steven A. Gabriel, "Spatial Transformation Sys-

tem Including Key Signal Generator," U.S. Patent 4,463,372, Ampex

Corp., July 3i, 1984.

Bennett, Phillip P. and Steven A. Gabriel, "System for Spatially

Transforming Images," U.S. Patent 4,472,732, Ampex Corp., Sep. 18,

1984.

Bergland, Glenn D., "A Guided Tour of the Fast Fourier Transform,"

Bernstein, Ralph and Harry Silverman, "Digital Techniques for Earth

Resource Image Data Processing," Proc. Amer. Inst. Aeronautics and

Astronautics 8th Annu. Meeting, vol. C21, AIAA paper no. 71-978,

October 1971.

Bemstein, Ralph, "Digital Image Processing of Earth Observation Sen-

sor Data," IBM J. Res. Develop., vol. 20, pp. 40-57, January 1976.

Bier, Eric A. and Ken R. Sloan, Jr., "Two-Part Texture Mappings,"

IEEE Computer Graphics and Applications, vol. 6, no. 9, pp. 40-53, Sep-

tember 1986.

Bizais, Y., I.G. Zubal, R.W. Rowe, G.W. Bennett, and A.B. Brill, "2-D

.Fitting and Interpolation Applied to Image Distortion Analysis," Pic-

torial Data Analysis, Ed. by R.M. Haralick, NATO ASI Series, vol. F4,

1983, pp. 321-333.

Blake, Andrew and Andrew Zisserman, Visual Reconstruction, MIT

Press, Cambridge, MA, 1987.

Blinn, James F. and Martin E. Newell, "Texture and Reflection in Com-

puter Generated Images," Comm. ACM, vol. 19, no. 10, pp. 542-547,

October 1976.

Boult, Terrarice E., "Visual Surface Reconstraction Using Sparse Depth

Data," Proc. IEEE Conference on Computer Vision and Pattern Recog-

nition, pp. 68-76, June 1986.

Boult, Terrance E., Information Based Complexity: Applications in Non-

linear Equations and Computer Vision, Ph.D. Thesis, Dept. of Computer

Sciense, Columbia University, NY, 1986.

Braccini, Carlo and Giuseppe Marino, "Fast Geometrical Manipulations

of Digital Images," Computer Graphics and Image Processing, vol. 13,

pp. 127-141, 1980.

Bracewell, Ron, The Fourier Transform and Its Applications, McGraw-

Hill, NY, 1986.

Briggs, I.C., "Machine Contouring Using Minimum Curvature," Geo-

physics, vol. 39, pp. 39-48, 1974.

[Brigham 88]

[Burt 88]

[Butler 87]

[Caelli 81]

[Casey 71]

[Catmull 74]

[Catmull 80]

[Chen 88]

[Clough 65]

[Cochran 67]

[Cook 84]

[Cook 86]

[Cooley 65]

[Cooley 67a]

[Cooley 67b]

303

Prentice-Hall, Englewood Cliffs, NJ, 1988.

Butt, Peter J., "Moment Images, Polynomial Fit Filters, and the Problem

of Surface Interpolation," Proc. 1EEE Conference on Computer Vision

and Pattern Recognition, pp. 144-152, June 1988.

Butler, David A. and Patricia K. Pierson, "Correcting Distortion in Digi-