Computing an interpolating cubic spline
From Wasteland
These page contains some pointers to and notes on practical issues about computing an interpolating cubic spline using C++ code.
Contents |
Description of the problem
The problem we are dealing with is that we have a set of ordered control points (typically in 1D, 2D or 3D), and we want to compute a univariate cubic spline f(t) that goes through them, i.e. we want to interpolates between them with a piecewise polynomial function of degree 3 (= order 4).
Reducing all problems to the 1D case
- If you have 1D points with an x-coordinate, then you solve
- an interpolation problem for (t, x)
- If you have 2D points with (x, y)-coordinates, then you solve
- an interpolation problem for (t, x)
- an interpolation problem for (t, y)
- If you have 3D points with (x, y, z)-coordinates, then you solve
- an interpolation problem for (t, x)
- an interpolation problem for (t, y)
- an interpolation problem for (t, z)
Note also that the matrix that describes the linear system you have to solve only depends on t, and doesn't need to be recomputed for each dimension.
Computing the knot vector
We can compute the corresponding knot value t(i) for each control point using Lee's (1989) centripetal scheme.
There are other schemes available (e.g. uniform), but are considered inferior in terms of the curve produced.
Also from a practical point of view, the centripetal scheme is used by Matlab's Spline Toolbox function cscvn, and it can be convenient to use the same method in your application if then you want to prototype or test things in Matlab.
Computing the spline coefficients
The polynomial for each segment of the spline looks like
f(t) = a + b * t + c * t² + d * t³
where t in [0, t(i+1) - t(i)].
Thus, the only thing missing is the value of the 4 coefficients a, b, c and d for each segment of the spline.
Theory to solve the problem
Computing the knot vector
The formula to compute the knot vector can be found in Lee (1989), as mentioned above, or in the help page for function cscvn.
Computing the spline coefficients
A good online description on how to compute the spline coefficients can be found in chapter 9, Global Cubic Space Curve Splines, of Knott (2000).
The bit of information missing is the relationship between the tangent vectors m(i) and the coefficients:
- a(i) = f(i)
- b(i) = m(i)
- c(i) = 3 (f(i+1)-f(i))/t²(i) - (2m(i) + m(i+1))/t(i)
- d(i) = (m(i) + m(i+1))/t²(i) - 2 (f(i+1)-f(i))/t³(i)
Software implementations to solve the problem
C++
Function QwtSpline from the Qwt C++ library computes the coefficients for the 1D case.
You need to compute the knot vector t externally, and then run QwtSpline for each dimension, e.g. (t, x), (t, y), etc.
Note that because a(i) is equal to the control points f(i), it only stores the other coefficients internally.
Note also that the coefficient notation is inverted with respect to Knott (2000):
- coefficientsA() in QwtSpline is d in Knott (2000)
- coefficientsB() in QwtSpline is c in Knott (2000)
- coefficientsC() in QwtSpline is b in Knott (2000)
Code snippet
This is the bit of code I use in Seg3D to implement the 3D spline interpolation
// declare arrays to compute the spline. The first dimension of
// the array is for the knot vector, while the second dimension
// will be for either the x-, y- or z-coordinates of the seed
// points, respectively
QwtArray<QwtDoublePoint> px( seeds_.size() + 1 );
/**
* Compute the knot vector for the spline using Lee's
* centripetal scheme
*
* E.T.Y. Lee. Choosing nodes in parametric curve interpolation,
* Computer-Aided Design, 21:363–370, 1989.
**/
// first knot is always 0.0
px[0].setX(0.0);
// compute inside knots
for ( size_t i = 1; i < seeds_.size(); ++i ) {
px[i].setX( sqrt( ( seeds_[i] - seeds_[i-1] ).length() ) );
}
// compute knot that closes the periodic spline
px[seeds_.size()].setX(
sqrt( ( seeds_[0] - seeds_[seeds_.size()-1] ).length() ) );
// cumulative summation of the knot values
for ( size_t i = 1; i < px.size(); ++i ) {
px[i].setX( px[i].x() + px[i-1].x() );
}
// the knot vector is the same for all coordinates note that we
// need to make a deep copy, the operator "=" would not work (it
// only makes a shallow copy)
QwtArray<QwtDoublePoint> py = px.copy(), pz = px.copy();
/**
* Compute the periodic cubic spline
**/
// assign seeds' coordinates to the point vectors
for ( size_t i = 0; i < seeds_.size(); ++i ) {
px[i].setY( seeds_[i][0] );
py[i].setY( seeds_[i][1] );
pz[i].setY( seeds_[i][2] );
}
px[seeds_.size()].setY( seeds_[0][0] );
py[seeds_.size()].setY( seeds_[0][1] );
pz[seeds_.size()].setY( seeds_[0][2] );
// create spline objects to interpolate the seeds
QwtSpline spx, spy, spz;
spx.setSplineType( QwtSpline::Periodic );
spy.setSplineType( QwtSpline::Periodic );
spz.setSplineType( QwtSpline::Periodic );
// compute spline coefficients for each coordinate
if ( !spx.setPoints( px ) ||
!spy.setPoints( py ) ||
!spz.setPoints( pz ) ) {
painter_->set_status("Internal error when creating spline.");
return;
}
/**
* Sample the periodic cubic spline so that we can draw an
* approximation with straight lines
**/
// compute sampling step on the spline parameter t
const double dt =
(px[px.size() - 1].x() - px[0].x())
/ (numPoints - 1);
// the points_ vector has the size of seeds_ when we are working
// with a linear spline, or until we have enough seed points in
// the cubic. But now points_ is going to contain the sampling
// of the cubic spline, so we need to make sure that it has the
// right length for it
points_.resize( this->numPoints );
// sample the spline
double t = px[0].x();
for( size_t i = 0; i < numPoints; i++ ) {
points_[i][0] = spx.value( t );
points_[i][1] = spy.value( t );
points_[i][2] = spz.value( t );
t += dt;
}
Matlab
Matlab's Spline Toolbox function cscvn takes as input the control points coordinates, (x, y, z), and does everything internally, to compute the corresponding spline.
References
- E.T.Y. Lee. Choosing nodes in parametric curve interpolation, Computer-Aided Design, 21:363–370, 1989.
- G.D. Knott, "Interpolating cubic splines", Brikaeuser, Boston, 2000.
