rhbartels@watcgl.waterloo.edu (Richard Bartels) (03/12/90)
I am posting a pair of procedures to evaluate general B-splines
(the basis functions) and their derivatives. They make up in generality
what they lose in speed. They are useful for checking your own
spiffy optimized stuff, and for getting up and running on that next
unusual case you have not yet bothered to work out in special detail.
If your machine is fast enough, you won't ever need anything else.
The evaluation routine for the basis functions derives from the book by
Whatshisname, Beatty, and Barsky. The derivative routine is better
than the one in the book and was contributed by Al Vermeulen of
the University of Waterloo, one of those rare and gifted people who
understands Sablionnier's work.
If you are used to dealing with basis matrices, in the style of the
material in Foley and Van Dam, then you think of tracing out the i-th
segment of a spline curve thusly:
3 2
[ t t t 1 ] [ Matrix ] [ Vtx[i-3] ]
[ Vtx[i-2] ]
[ Vtx[i-1] ]
[ Vtx[i] ]
where Vtx[i] is the i-th control vertex, and this example assumes
cubics. The corresponding understanding you must have to use the code
successfully is that the evaluation procedure provides the product
of the t-posers vector with the basis matrix:
3 2
[ t t t 1 ] [ Matrix ] = [ bvals[3] bvals[2] bvals[1] bvals[0] ]
The desired point on the curve is given by the dot product:
[ bvals[3] bvals[2] bvals[1] bvals[0] ] [ Vtx[i-3] ]
[ Vtx[i-2] ]
[ Vtx[i-1] ]
[ Vtx[i] ]
For higher degrees the vectors are correspondingly longer; for lower
degrees they are shorter. The number of components in the vectors equals
the order of the spline. For the derivative routine, the dvals array stores
the derivative values of the basis functions according to the same scheme.
To figure out what knots to use with what vertices and what parameter
value, begin counting knot indices at zero, and start the first curve
segment after the d-th-indexed knot (where d is the degree of the spline).
Thus, cubics are degree 3, so the first segment of the curve would
begin "as soon as possible" after knots[3]. Segments of the curve
correspond to non-vacuous intervals between knots, so "as soon as possible"
means at the first index for which knots[index] < knots[index+1].
For cubics with single knots, index=3, so the first segment corresponds
to parameter values knots[3] <= parval < knots[4]. The control vertices
have corresponding indices: Vtx[3], Vtx[2], Vtx[1], Vtx[0], which is just
what corresponds to the indexing that most Foley and Van Dam readers
are used to.
The only extra consideration that is necessary is to have as many
additional knots on the right of the last segment as one has on the left
of the first; i.e. degree+1 = order. Thus, to create one segment of a
uniform quartic B-spline curve, for example, use knots:
0 1 2 3 4 5 6 7 8 9
Knots 0, 1, 2, 3, 4 are to the left of the parametric segment,
knots 5, 6, 7, 8, 9 are to the right of the parametric segment,
and the single segment of the curve corresponds to the interval
from parval=4 to parval=5.
As another example, two abutting cubic Bezier segments would have the knot
sequence:
0 0 0 0 1 1 1 2 2 2 2
The first Bezier segment corresponds to the parametric interval from
knots[3]=0 to knots[4]=1, the second segment runs from knots[6]=1 to
knots[7]=2. The control vertices for the first segment will be
Vtx[0], Vtx[1], Vtx[2], Vtx[3], and the control vertices for the
second segment will be Vtx[3], Vtx[4], Vtx[5], Vtx[6]. The reader can
check that these indices follow the rules outlined in the previous paragraphs.
In general, the value of the knots is not important; only the spacing
between the knots is important. The Bezier knots could just as
easily be:
5 5 5 5 6 6 6 7 7 7 7
If you want NURBS (non-uniform rational B-splines), then to each basis
function you will have an associated weight. You may also regard
this as a weight on the corresponding vertex, and this will tell you
which weight index to use with which bval. The dot product above
will be replaced by the quotient:
w[i-3]*bvals[3]*Vtx[i-3]
+ w[i-2]*bvals[2]*Vtx[i-2]
+ w[i-1]*bvals[1]*Vtx[i-1]
+ w[i]*bvals[0]*Vtx[i]
-----------------------------------
w[i-3]*bvals[3]
+ w[i-2]*bvals[2]
+ w[i-1]*bvals[1]
+ w[i]*bvals[0]
-Richard
-------------------------- cut here -----------------------------------------
/* ============================================================================
NUBSplineValues()
NUBSplineDerivs()
============================================================================ */
void NUBSplineValues( parval, order, knots, index, bvals )
int order, index;
double *knots, *bvals, parval;
{
/*-----------------------------------------------------------------------------
Purpose:
Carry out the Cox-deBoor recurrence to compute the nonzero B-spline
values at a given value of the parameter.
Input:
parval --- parameter value for evaluation
order ---- order of spline
knots ---- knot vector portion s.t. knots[index] <= parval < knots[index+1]
index ---- index for knots indicating position of parval
Output:
bvals[i] ---- for i=0,...,order-1 containing the nonzero B-spline values
B[index,order](parval),...,B[index-order+1,order](parval),
where bvals[j] = B[index-j,order](parval)
Caution:
No validity checking is done!
Make sure parval is not to the left of the first segment or
to the right of the last. Make sure that knots[index] < knots[index+1],
and make sure there are enough knots to the left of the first segment
and to the right of the last. If you don't, you'll be sooorrreeee.
-----------------------------------------------------------------------------*/
register int r,i,s;
register double omega;
if( order<=0 ) return;
bvals[0] = 1.0;
for ( r=2; r<=order; r++ ) {
i = index-r+1;
bvals[r-1] = 0.0;
for ( s=r-2; s>=0; s-- ) {
i++;
omega = ( parval - knots[i] ) / ( knots[i+r-1] - knots[i] );
bvals[s+1] += ( 1.0 - omega ) * bvals[s];
bvals[s] *= omega;
}
}
return;
} /* NUBSplineValues */
void NUBSplineDerivs( parval, order, deriv, knots, index, dvals )
int order, deriv, index;
double *knots, *dvals, parval;
{
/*-----------------------------------------------------------------------------
Purpose:
Carry out Sablionnier's tetrahedral recurrence to compute the
derivatives of the nonzero B-spline values at a given parameter.
Input:
parval --- parameter value for evaluation
order ---- order of spline
deriv ---- the desired derivative order
knots ---- knot vector portion s.t. knots[index] <= parval < knots[index+1]
index ---- index for knots indicating position of parval
Output:
dvals[i] ---- for i=0,...,order-1 containing the deriv-th
derivatives of the nonzero bspline values
B[index,order](parval),...,B[index-order+1,order](parval),
where dvals[j] = D^deriv B[index-j,order](parval)
Caution:
No validity checking is done!
Read the cautionary note for NUBSplineValues,
and consider well what thou doist.
-----------------------------------------------------------------------------*/
register int r,i,s;
register double omega;
for( r=0; r<order; r++ ) {
dvals[r] = 0.0;
}
if( (deriv>=order) || (deriv<0) ) return;
NUBSplineValues( parval, order-deriv, knots+deriv, index-deriv, dvals );
for( r=order-deriv+1; r<=order; r++ ) {
i = index-r+1;
for( s=r-2; s>=0; s-- ) {
i++;
omega = ((double)(r-1)) / (knots[i+r-1] - knots[i]);
dvals[s+1] += -omega * dvals[s];
dvals[s] *= omega;
}
}
return;
} /* NUBSplineDerivs */