Copyright 1984 by ABComputing July 15, 1984 ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º What The !@#$%& is a Bicubic Spline?? º º º º by º º º º Donald M. Esterling º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Introduction ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ Piecewise-cubic splines are a marvelous way to "fair" or, more accurately, interpolate a "smooth" curve through data points in a two-dimensional x-y coordinate frame. Interpolation should be contrasted with fitting, where the curve only approximates the data points by, e.g., a straight line. Splines received their name from a draftsman's spline, which is an elastic strip that can be adjusted to pass through a fixed set of points. Usually the result is a "nice" curve resembling a free-hand curve through these data points. There are many other interpolation methods, but they can often lead to unpleasant results (polynomial "wiggles.") But we don't live in a two-dimensional world. In addition, not every function has only one independent parameter. Suppose a set of data is defined on (or defines) a surface. As an example, think of a checkerboard. The board makes up an x-y plane. Now assign some height z(x,y) to each corner on the checkerboard (pull the corner up or down). If the checkerboard were made of reasonably stiff rubber, then we would have a smooth surface (rather than a smooth line) through our set of data points z=z(x,y). A bicubic spline is a reasonable approximation to that smooth surface. Bicubic splines are a widely used tool in the arsenal of finite-element modelers and CAD/CAM programmers. When I sought a formal description of them in numerical methods texts, I was surprised to find either no mention or that they were dealt with summarily as an exercise for the interested reader (a common subterfuge used by authors). I finally found a text that gave a neat trick to generate bicubic splines and that's what I'd like to share with you. If I had infinite space (and your attention), before moving on to bicubic splines I would first digress to remind you about interpolation in general and splines in particular. Lacking that, for more detailed information you'll just have to browse through a numerical methods text at your local bookstore. Generally, splines do not appear in these texts unless they are fairly recent ( <= 5 years old or so). ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Piecewise-cubic Spline ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ So on to our neat trick...First I do have to define a simple piecewise-cubic spline. It is a collection of cubic polynomials that (a) passes through the data points, (b) is continuous, and (c) has continuous first and second derivatives. Formally, for the interval x(i) <= x <= x(i+1), the interpolating function y(x) is y(x) = c(1,i) + c(2,i)*d + c(3,i)*(d**2)/2.0 + c(4,i)*(d**3)/6.0 (1) where d = x-x(i) and the c's are constants determined by conditions (a), (b), and (c). ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Cardinal Spline ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ The next ingredient in the making of a bicubic spline is called a "cardinal spline." This is a cubic spline that fits special data values of the form y(k) = 1, i equal to k (2) y(i) = 0, i not equal k where "k" has a fixed value. We will call this cardinal spline P(k,x). A different cardinal spline is generated for each value of k. Using these cardinal splines, a general set of data points (x(i), y(i)), i=1,2,...,n can be interpolated by the function U given by: n U(x) = ä P(i,x)*y(i) i=1 The routine QINTERP, provided on Diskette A as LIST1.FOR, returns the coefficients for a general cubic spline (the constants "c"). (This routine is one among many contained in the GRAFMATIC 2-D, 3-D graphics package.) The input to QINTERP (n,x,y,c) is: a. the number of data points (n), and b. the data values (x(i), y(i)) i=1,2,...,n The variable "n" must be typed as integer and the arrays x(n), y(n), and c(4,n) must be defined (dimensioned) in the calling program. To obtain a cardinal spline, input the special data specified in (2). The array c(4,n), used in (1), is returned as output. One last technical point - cubic splines come in various guises, depending on how one defines the end-point boundary conditions. I have used the oddball but useful boundary condition that fixes the coefficients for the polynomial for the first and second interval to be the same and the coefficients for the penultimate and last polynomials to be the same (the so-called "not-a-knot" boundary condition). This is discussed further in deBoor's book cited in the source code. The aforementioned source code also handles more standard boundary conditions as indicated in the comment lines. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Bicubic Spline ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ Where's the bicubic spline? This article is almost over and I'll bet you think I haven't delivered. But I have. If the data points (x(i,j), y(i,j), z(i,j)), i=1,2...,n, j=1,2,...,m, are defined on a surface, then a smooth interpolant (bicubic spline) through the data points can be obtained from n m U(x,y) = ä ä P(i,x) * Q(j,y) * z(i,j) i=1 j=1 where: P(i,x) is a cardinal spline on x, i.e., P(i,x(i))=1, P(i,x(k))=0 for i not equal to k and Q(j,y) is a cardinal spline on y, i.e. Q(j,y(j))=1, Q(j,y(k))=0 for j not equal to k. By construction U(x(r),y(s))=z(r,s) for all r=1,2,...,n and s= 1,2,...,m. In summary, to create a bicubic-spline interpolant, fashion a set of cardinal splines for your x and y coordinates with QINTERP. These functions stay fixed. Then any set of data x(i,j) can be interpolated by the smooth 2D surface generated by U(x,y). ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ DEMO Program ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ Diskette B of PCFL/PCUG includes a DEMO program on cubic-spline interpolation called "SPLINE.EXE", demonstrating "polynomial wiggles" and, not coincidentally, the graphics power of GRAFMATIC. EDITOR'S NOTE: Due to insufficient disk space, the source code for the DEMO was not included in Diskette B. Any references in the DEMO to this source code should be disregarded. As the DEMO was produced by the Microsoft FORTRAN compiler: "Portions copyright by Microsoft Corporation, 1982, 1983." ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ ³ ³ Donald M. Esterling is a veteran of two decades of ³ ³ programming, obtaining funds for computer simulation ³ ³ research from NSF, NBS, NASA and the USAF. He is founder of ³ ³ Microcompatibles, author of GRAFMATIC and PLOTMATIC ³ ³ (scientific/engineering graphics/pen plot support), and in ³ ³ his spare time is a Professor of Engineering, laboring in ³ ³ materials science, solid state physics, and, recently, ³ ³ (thanks to a generous grant from IBM) computer-aided ³ ³ manufacturing. ³ ³ ³ ³ He may be reached at Microcompatibles, 11443 Oak Leaf Drive, ³ ³ Silver Spring, MD 20901 (301) 593-0683. ³ ³ ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ File Name: ÛÛ fort1.txt ÛÛ ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ