Spline Surface (in 3D)Written by Paul BourkeNovember 1996
Creating a spline surface involves taking the product of the same spline blending functions used for spline curves as follows
where the control points form a 2D array P_{ij}. Most of the properties of the spline curve also apply to spline surfaces. For example
As an example the following show the surface given a 4x5 matrix of control points, degree 3 in both directions, the surface is sampled on a 30x40 grid.
Zero order continuity between spline surfaces is ensured if the control points along the joining edge are the same. First order continuity can be obtained by matching the control point slope across the boundary. This also applies to attempts to form a closed spline surface. Forming a closed surface is easier if third order splines are being used because the surface is tangential to the end planes. In the following image the surface on the right starts and ends at the origin. The surface on the right has an extra control point so that the start and end planes have the same tangent.
C Source Example The following example show how to create a spline surface in the images above using Geomview. /* Create an example of a spline surfaces */ #define NI 3 #define NJ 4 XYZ inp[NI+1][NJ+1]; #define TI 3 #define TJ 3 int knotsI[NI+TI+1]; int knotsJ[NJ+TJ+1]; #define RESOLUTIONI 30 #define RESOLUTIONJ 40 XYZ outp[RESOLUTIONI][RESOLUTIONJ]; int main(argc,argv) int argc; char **argv; { int i,j,ki,kj; double intervalI,incrementI; double intervalJ,incrementJ; double bi,bj; /* Create a random surface */ srandom(1111); for (i=0;i<=NI;i++) { for (j=0;j<=NJ;j++) { inp[i][j].x = i; inp[i][j].y = j; inp[i][j].z = (random() % 10000) / 5000.0  1; } } /* Step size along the curve */ incrementI = (NI  TI + 2) / ((double)RESOLUTIONI  1); incrementJ = (NJ  TJ + 2) / ((double)RESOLUTIONJ  1); /* Calculate the knots */ SplineKnots(knotsI,NI,TI); SplineKnots(knotsJ,NJ,TJ); intervalI = 0; for (i=0;i<RESOLUTIONI1;i++) { intervalJ = 0; for (j=0;j<RESOLUTIONJ1;j++) { outp[i][j].x = 0; outp[i][j].y = 0; outp[i][j].z = 0; for (ki=0;ki<=NI;ki++) { for (kj=0;kj<=NJ;kj++) { bi = SplineBlend(ki,TI,knotsI,intervalI); bj = SplineBlend(kj,TJ,knotsJ,intervalJ); outp[i][j].x += (inp[ki][kj].x * bi * bj); outp[i][j].y += (inp[ki][kj].y * bi * bj); outp[i][j].z += (inp[ki][kj].z * bi * bj); } } intervalJ += incrementJ; } intervalI += incrementI; } intervalI = 0; for (i=0;i<RESOLUTIONI1;i++) { outp[i][RESOLUTIONJ1].x = 0; outp[i][RESOLUTIONJ1].y = 0; outp[i][RESOLUTIONJ1].z = 0; for (ki=0;ki<=NI;ki++) { bi = SplineBlend(ki,TI,knotsI,intervalI); outp[i][RESOLUTIONJ1].x += (inp[ki][NJ].x * bi); outp[i][RESOLUTIONJ1].y += (inp[ki][NJ].y * bi); outp[i][RESOLUTIONJ1].z += (inp[ki][NJ].z * bi); } intervalI += incrementI; } outp[i][RESOLUTIONJ1] = inp[NI][NJ]; intervalJ = 0; for (j=0;j<RESOLUTIONJ1;j++) { outp[RESOLUTIONI1][j].x = 0; outp[RESOLUTIONI1][j].y = 0; outp[RESOLUTIONI1][j].z = 0; for (kj=0;kj<=NJ;kj++) { bj = SplineBlend(kj,TJ,knotsJ,intervalJ); outp[RESOLUTIONI1][j].x += (inp[NI][kj].x * bj); outp[RESOLUTIONI1][j].y += (inp[NI][kj].y * bj); outp[RESOLUTIONI1][j].z += (inp[NI][kj].z * bj); } intervalJ += incrementJ; } outp[RESOLUTIONI1][j] = inp[NI][NJ]; printf("LIST\n"); /* Display the surface, in this case in OOGL format for GeomView */ printf("{ = CQUAD\n"); for (i=0;i<RESOLUTIONI1;i++) { for (j=0;j<RESOLUTIONJ1;j++) { printf("%g %g %g 1 1 1 1\n", outp[i][j].x, outp[i][j].y, outp[i][j].z); printf("%g %g %g 1 1 1 1\n", outp[i][j+1].x, outp[i][j+1].y, outp[i][j+1].z); printf("%g %g %g 1 1 1 1\n", outp[i+1][j+1].x,outp[i+1][j+1].y,outp[i+1][j+1].z); printf("%g %g %g 1 1 1 1\n", outp[i+1][j].x, outp[i+1][j].y, outp[i+1][j].z); } } printf("}\n"); /* Control point polygon */ for (i=0;i<NI;i++) { for (j=0;j<NJ;j++) { printf("{ = SKEL 4 1 \n"); printf("%g %g %g \n",inp[i][j].x,inp[i][j].y,inp[i][j].z); printf("%g %g %g \n",inp[i][j+1].x,inp[i][j+1].y,inp[i][j+1].z); printf("%g %g %g \n",inp[i+1][j+1].x,inp[i+1][j+1].y,inp[i+1][j+1].z); printf("%g %g %g \n",inp[i+1][j].x,inp[i+1][j].y,inp[i+1][j].z); printf("5 0 1 2 3 0\n"); printf("}\n"); } } }
Spline curves (in 3D)Written by Paul BourkeNovember 1996
Spline curves originate from flexible strips used to create smooth curves in traditional drafting applications. Much like Bezier curves they are formed mathematically from piecewise approximations of cubic polynomial functions with zero, first and second order continuity.
BSplines are one type of spline that are perhaps the most popular
in computer graphics applications, they are defined as follows:
where N(v) are called blending functions, u_{k} are known as break points, where they occur on the curve are known as knots. There are a number of possible options for the knot positions, for example a uniform spacing where u_{k} = k. More commonly the following function is chosen The blending functions are defined as
Such curves have many advantages
There is another curve in the family called BetaSplines which have additional parameters beta1 and beta2 which adjust the shape relative to the convex hull. BSplines are subsets of BetaSplines. The following example (created with the source below) shows the spline through 4 control points, with degree t = 2, 3, and 4. t = 2 is just linear interpolation, as the degree increases the smoother the curve becomes.
The pink line shows the polygon formed by the control points, the white line is the spline curve. Note that the curve lies within the convex hull of the points. Note also that the tangent of the curve at the endpoints is the same as the slope of the line between the first two or last two control points. This makes for convenient matching of the curve between sections. For example in the following the spline on the left has control points starting and ending at the origin. The curve on the right has control points on the same square but the start and end control point is along the bottom edge and the slope of the first two and last two control points is the same.
Unfortunately the curve always passes through the first and last point. The following has the same control points as for the curve on the right above but for increasing values of the degree, t.
Splines of degree 3 are by far the most commonly used in practice. Except for the first edge and last edge of the control point polygon these splines touch the midpoint of each other edge. If the control point polygon is concave the spline passes through the midpoint, if the polygon is concave the spline is tangential at the edge midpoint. C Source Code
In what follows it will assumed there are n + 1 control points, P_{k}
where k ranges from 0 to n. The parameter t determines the order
of the polynomials used for the blending function, typically 3
or 4.
Following are three functions which implement spline curves, at the end there is a small demonstration program which calls these routines to create the image shown earlier. /* This returns the point "output" on the spline curve. The parameter "v" indicates the position, it ranges from 0 to nt+2 */ void SplinePoint(int *u,int n,int t,double v,XYZ *control,XYZ *output) { int k; double b; output>x = 0; output>y = 0; output>z = 0; for (k=0;k<=n;k++) { b = SplineBlend(k,t,u,v); output>x += control[k].x * b; output>y += control[k].y * b; output>z += control[k].z * b; } } /* Calculate the blending value, this is done recursively. If the numerator and denominator are 0 the expression is 0. If the deonimator is 0 the expression is 0 */ double SplineBlend(int k,int t,int *u,double v) { double value; if (t == 1) { if ((u[k] <= v) && (v < u[k+1])) value = 1; else value = 0; } else { if ((u[k+t1] == u[k]) && (u[k+t] == u[k+1])) value = 0; else if (u[k+t1] == u[k]) value = (u[k+t]  v) / (u[k+t]  u[k+1]) * SplineBlend(k+1,t1,u,v); else if (u[k+t] == u[k+1]) value = (v  u[k]) / (u[k+t1]  u[k]) * SplineBlend(k,t1,u,v); else value = (v  u[k]) / (u[k+t1]  u[k]) * SplineBlend(k,t1,u,v) + (u[k+t]  v) / (u[k+t]  u[k+1]) * SplineBlend(k+1,t1,u,v); } return(value); } /* The positions of the subintervals of v and breakpoints, the position on the curve are called knots. Breakpoints can be uniformly defined by setting u[j] = j, a more useful series of breakpoints are defined by the function below. This set of breakpoints localises changes to the vicinity of the control point being modified. */ void SplineKnots(int *u,int n,int t) { int j; for (j=0;j<=n+t;j++) { if (j < t) u[j] = 0; else if (j <= n) u[j] = j  t + 1; else if (j > n) u[j] = n  t + 2; } } /* Create all the points along a spline curve Control points "inp", "n" of them. Knots "knots", degree "t". Ouput curve "outp", "res" of them. */ void SplineCurve(XYZ *inp,int n,int *knots,int t,XYZ *outp,int res) { int i; double interval,increment; interval = 0; increment = (n  t + 2) / (double)(res  1); for (i=0;i<res1;i++) { SplinePoint(knots,n,t,interval,inp,&(outp[i])); interval += increment; } outp[res1] = inp[n]; } /* Example of how to call the spline functions Basically one needs to create the control points, then compute the knot positions, then calculate points along the curve. */ #define N 3 XYZ inp[N+1] = {0.0,0.0,0.0, 1.0,0.0,3.0, 2.0,0.0,1.0, 4.0,0.0,4.0}; #define T 3 int knots[N+T+1]; #define RESOLUTION 200 XYZ outp[RESOLUTION]; int main(int argc,char **argv) { int i; SplineKnots(knots,N,T); SplineCurve(inp,N,knots,T,outp,RESOLUTION); /* Display the curve, in this case in OOGL format for GeomView */ printf("LIST\n"); printf("{ = SKEL\n"); printf("%d %d\n",RESOLUTION,RESOLUTION1); for (i=0;i<RESOLUTION;i++) printf("%g %g %g\n",outp[i].x,outp[i].y,outp[i].z); for (i=0;i<RESOLUTION1;i++) printf("2 %d %d 1 1 1 1\n",i,i+1); printf("}\n"); /* The axes */ printf("{ = SKEL 3 2 0 0 4 0 0 0 4 0 0 2 0 1 0 0 1 1 2 1 2 0 0 1 1 }\n"); /* Control point polygon */ printf("{ = SKEL\n"); printf("%d %d\n",N+1,N); for (i=0;i<=N;i++) printf("%g %g %g\n",inp[i].x,inp[i].y,inp[i].z); for(i=0;i<N;i++) printf("2 %d %d 0 1 0 1\n",i,i+1); printf("}\n"); }References Chankin, G. An algorithm for high speed curve generation. Computer Graphics and Image Processing, 3 (1974), 346349 Riesenfeld, R. On Chaikin's algorithm. IEEE Computer Graphics and Applications, 4, 3 (1975) 304310
