From: horst.kraemer@berlin.snafu.de (Horst Kraemer) Newsgroups: sci.math Subject: Re: Cubic Spline? Date: Wed, 29 Oct 1997 12:29:40 GMT On Tue, 28 Oct 1997 15:12:06 +0000, B Skinner wrote: >I need to know about cubic spline; basically how does one go about >aproximating a curve when given a (large) number of points? >Is it a case of taking four points at a time? Or what? The conceptional approach is as follows: We want to interpolate a sequence of N points P_j in R^n by a curve in parameter form K(u). In our case the conditions are: for every coordinate i the curve segment between P_j and P_j+1 ist a polynomial of degree <= 3. For two adjacent polynomials the value and the first and second derivative (i.e. value, tangent and curvature) are identical at the common point. These constraints still leave two degrees of freedom for open curves. One may assume that the curvatures of the curve at the first and last point are 0 (natural spline). For closed curves the curve is uniquely defined by the constraints. These conditions may be condensed into a system of N+2 linear equations with N+2 unknowns which may be solved via direct substitution without pivoting because the system is tridiagonal with a dominant main diagonal. The main problem is the representation of the curve, i.e. the representation of the coefficients of the polynomial pieces. The B-Spline representation offers an elegant and efficient way without redundancy. In theory you need 4 coefficients per dimension to represent a polynomial piece of degree <=3, i.e. dim*4*(N-1) for the whole curve. With the B-Spline representation you will need only dim*(N+2) coefficients. An answer to your last question: It is not enough to consider points P1,P2,P3,P4 in order to calculate the curve segment with a polynomial degree <=3 between P2 and P3 if you want the curve to be "smooth" in the sense that tangents _and_ curvature coincide at the endpoints of the segments. The spline interpolation scheme is a global interpolation scheme, i.e. for every segment one needs information about all points of the curve. The theory about B-Splines is rather involved although the calculations needed for an interpolation are quite simple. You should seek general written information in order to understand what's going on. Upon request I can give you an example program in Turbo Pascal or Borland C. The program will not document the theoretical background. Regards Horst ============================================================================== Date: Sat, 29 May 99 15:52:17 -0400 From: Ken Shoemake To: "Dave Rusin" Subject: [What are splines] Let me offer a different answer. A cubic spline curve is a piecewise polynomial curve. These come in several different flavors, depending on what's important to you. You may choose any two of three options: o Pass through the data points o Pieces connect with second derivative continuity o Each data point affects only a few pieces The answer given by Kraemer chooses the first two options. If your data points are dense and perhaps noisy, you may prefer to choose the last two options. Either way, you have some additional choices. One is how you handle a little extra freedom at the loose ends. A more important choice, if the data points are unevenly spaced, is how you parameterize the curve. To implement the approximation choice (continuity and locality) you can proceed as follows. Suppose you have successive interior data points p, q, r, and s. Further suppose the parametric cuts are at a, b, c, d, e, and f. Parameter x, with c < x <= d, will give a point pqrs blended in three stages: p [(d-x)p+(x-a)q]/(d-a) = pq q [(d-x)pq+(x-b)qr]/(d-b) = pqr [(e-x)q+(x-b)r]/(e-b) = qr [(d-x)pqr+(x-c)qrs]/(d-c) = pqrs r [(e-x)qr+(x-c)rs]/(e-c) = qrs [(f-x)p+(x-c)q]/(f-c) = rs s As x varies from c to d, the result point is a weighted sum of p, q, r, and s. The weights are polynomials in x that depend on the parameter cut points (called "knots") a, b, c, d, e, and f. The weight polynomials are called B-splines. Often the curve itself is also called a B-spline -- just to confuse everyone. These approximating curves have many nice properties, and are often preferred to the global interpolating curves Kraemer describes. Generally, they are cheaper and better behaved. Variations of this blending generate many kinds of spline curves. You can find a generic routine, DialASpline, in "Graphics Gems V", A. Paeth, ed., AP Professional. All source code is freely available online at various sites, including Beneath the simplicity of DialASpline is the elegant theory of polar forms. Recall that Pythagoras claims a planar right triangle with side lengths a and b and diagonal c satisfies c^2 = a^2 + b^2. Transfering this to vectors in R^2 written in terms of an orthonormal basis {x,y}, the length of a vector V = (Vx,Vy) is Vx^2 + Vy^2. From this we deduce that the sum of vectors U and V differs from the Pythagorean value by (U+V)^2 - U^2 - V^2, which turns out to be 2*(Ux*Vx + Uy*Vy). This is twice the dot product (a.k.a. scalar product). The dot product of U and V is symmetric and linear in both U and V, and equals the Pythagorean quadratic value when U = V. More generally, any quadratic has a polar form with these same properties. Less commonly seen in mathematics books, any multivariate polynomial of any degree can be polarized. For example, a*X^3 + b*X^2 + c*X + d becomes a*(X1*X2*X3) + b*(X1*X2 + X2*X3 + X3*X1)/3 + b*(X1 + X2 + X3)/3 + c The B-spline curve data points (a.k.a. control points) are intimately related to the polar forms of the polynomial pieces. This enormously simplifies the theory of spline-based curves and surfaces. ---------- If you feel like including code, this is it: /****** lincrv.c ******/ /* Ken Shoemake, 1994 */ #define MAXDIM 2 typedef float Vect[MAXDIM]; typedef float Knot; typedef int Bool; int DialASpline(Knot t, Knot a[], Vect p[], int m, int n, Vect work[], unsigned int Cn, Bool interp, Vect val); /* Perform a generic vector unary operation. */ #define V_Op(vdst,gets,op,vsrc,n) {register int V_i;\ for(V_i=(n)-1;V_i>=0;V_i--) (vdst)[V_i] gets op ((vsrc)[V_i]);} static void lerp(Knot t, Knot a0, Knot a1, Vect p0, Vect p1, int m, Vect p) { register Knot t0=(a1-t)/(a1-a0), t1=1-t0; register int i; for (i=m-1; i>=0; i--) p[i] = t0*p0[i] + t1*p1[i]; } /* DialASpline(t,a,p,m,n,work,Cn,interp,val) computes a point val at parameter t on a spline with knot values a and control points p. The curve will have Cn continuity, and if interp is TRUE it will interpolate the control points. Possibilities include Langrange interpolants, Bezier curves, Catmull-Rom interpolating splines, and B-spline curves. Points have m coordinates, and n+1 of them are provided. The work array must have room for n+1 points. */ int DialASpline(Knot t, Knot a[], Vect p[], int m, int n, Vect work[], unsigned int Cn, Bool interp, Vect val) { register int i, j, k, h, lo, hi; if (Cn>n-1) Cn = n-1; /* Anything greater gives one polynomial */ for (k=0; t> a[k]; k++); /* Find enclosing knot interval */ for (h=k; t==a[k]; k++); /* May want to use fewer legs */ if (k>n) {k = n; if (h>k) h = k;} h = 1+Cn - (k-h); k--; lo = k-Cn; hi = k+1+Cn; if (interp) { /* Lagrange interpolation steps */ int drop=0; if (lo<0) {lo = 0; drop += Cn-k; if (hi-lon) {hi = n; drop += k+1+Cn-n; if (hi-lo=j; i--) { lerp(t,a[lo+i],a[lo+i+tmp],work[lo+i],work[lo+i+1],m,work[lo+i+1]); } } V_Op(val,=,,work[lo+h],m); return (k); } -- Ken