y
hi 2(hi+h2) h2 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
A2.4.3. Boundary Conditions: Free-end, Cyclic, and Not-A-Knot
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
y-t 2 - 2Ax (A2.4.5)
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-
points of the span are set equal to each other.
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
A Xn-3 Xn-2 ß ß
The version of this boundary condition expressed in terms of second derivatives is left to
290
the reader as an exercise.
Thus far we have placed no restrictions on the spacing between the data points.
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
141
14
4
y
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
* ip is interval's left endpoint (integer)
* oip is left endpoint of last visited interval
*/
oip = -1; /* force coefficient initialization '/
fctr = (double) (lenl-1) / (len2-1);
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 */
x=p-ip;
Y2[i] = ((A3*x + A2)*x + A1)*x + A0;
/* increment pointer */
ip = (p += fctr);
}
cfree((char *) YD);
}
291
292
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];
for(i = 1; i < len-1; i++) Y D[i]=3.0*(Y[i+ 1 ]-Y[i- 1]);
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;
double *C;
? 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
A2.5.2. Ispline_gen
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
Input: Y1 is a list of irregular data points (lenl entries)
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;
double *YD, A0, A1, A2, A3, x, dx, dy, pl, p2, p3;
/* 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
* p2 is resampling position
* p3 is dght endpoint of interval
* j is input index for current interval
*/
p3 = X210] - 1; /* force coefficient initialization */
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
x=p2-pl;
Y2[i] = ((A3*x + A2)*x + A1)*x + A0;
}
cfree((char *) ¾D);
YD <- Computed 1st derivative of data in X,Y (len entdes)
The not-a-knot boundary condition is used
getYD_gen(X,Y,YD, len)
double *X, *Y, *YD;
int len;
{
int i;
double h0, hl, r0, rl, *A, *B, *C;
/* 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 */
A[i] = (h0+hl) * (h0+hl);
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
296
t ridiag_g en(A,B,C,D,len)
double *A, *B, *C, *D;
int len;
int i;
double b, *F;
F = (double *) calloc(len, sizeof(double));
/* 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 */
for(i = len-2; i >= 0; i--) D[i] -= (D[i+l] * F[i+l]);
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)
The first forward difference forf (x) is expressed as
297
298
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-
ence is used at each iteration to update the first forward difference which, in turn, is
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
the iteration.
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)
= 3a3x 2 +(3a3 +2a2)x +a3 +a2+al
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
term that must be updated, i.e., A2f. Even so, this method offers the benefit of computing
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.
The basis of the entire iteration is shown in the top row. For consistency with our discus-
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.
This cycle continues until the last diagonal is reached, denoting that the entire span of
points has been evaluated.
Position Value
XO UO
Xl Ul Au 0
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.
Anderson, Scott E. and Mark A.Z. Dippe, ' 'A Hybrid Approach to Facial
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.
Antoniou, Andreas, Digital Filters: Analysis and Design, McGraw-Hill,
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
302 REFmESCgS
[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,"
IEEE Spectrum, vol. 6, pp. 41-52, July i969.
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
Brigham, E. Oran, The Fast Fourier Transform and Its Applications,
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-
**Share with your friends:** |