Eigen library: spline interpolation vs spline smoothing

Question

The C++ library Eigen provides an "unsupported" splines module which is giving me troubles.

The task is typical: a data-acquisition device provides me with a time series that is sort-of regular and has some holes with missing data. I want to fill in the gaps, for maybe a dozen data points, using a high-order polynomial.

I'm familiar with a spline function defined by data points or knots which it passes through. The Eigen math library provides a Spline class, which accepts "knot vector" and a matrix of "control points." But the resulting function doesn't actually pass through the control points at the knot values. I'm not sure what the quantitative relationship is to the input data… would like to know.

It has another submodule which provides functions that do pass through the original, given data points. But for some reason that is implemented using QR decomposition, which is slow unless I break the data into chunks. Then the resulting function goes way out of bounds for any missing data point, and the problem gets worse if I raise the polynomial degree. The greatest deviation occurs at the first interpolated point and then the function exponentially returns back to the expected range. So there's something wrong with the coefficients. I would expect it to at least be symmetrical.

The library's source code cites the NURBS book (Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data).

Is there something fundamental I'm missing about splines? Is there a quick fix? I'm just trying to take the shortest route to interpolation. The function I'm interpolating isn't a polynomial, and harmonic analysis might be better. But I'm not aware of a suitable tool.

Answer 1

The function that accepts a "knot vector" and "control points" is probably using this data to construct a b-spline. B-spline curves don't pass through (interpolate) their control points.

To construct a spline curve that interpolates the input data points, you typically have to solve a linear system, so QR decomposition is not surprising.

If you use a low-degree spline (say cubic, degree = 3), then the linear system will be "banded" and therefore easier/faster to solve. Also, a spline with low degree is less likely to oscillate wildly.

This page has an example that shows creation of an interpolating spline with the Eigen package: http://forum.kde.org/viewtopic.php?f=74&t=98871 (I think).

Alternatively, there are many other packages available for creating splines. Just search for "spline interpolation code".

Answer 2

As you know, B-splines smooth, run between the data points. If you want to interpolate, go through the data points exactly, using a B-spline function, here's how.

In: data points $f_i$ at the integers $i = 0\ 1\ 2 \dots n$
$\ \ \ \ \ $ a cubic B-spline function $B( x; data_i )$ for $0 \leq x \leq n$
Out: a wrapper function $I( x; data_i )$ that interpolates, $I(i) = f_i$ .

Method: we know that $B(i) = (f_{i-1} + 4 f_i + f_{i+1}) \,/\, 6, i = 0\ 1\ 2 \dots$
or in matrix form, $\ b = A f$ with $A$ tridiagonal.
Solve $\ A y = f$
then $\ \ B(y) = B( A^{-1} f ) = f$
i.e. $\ \ \ \ \, I( x; f_i ) \equiv B( x; A^{-1} f )$ interpolates.

Notes:
Solving this tridiagonal system is fast and easy. (You need boundary conditions such as $f_{-1} = f_0$ or $f_{-1} = 0$ , and similarly for $f_{n+1}$.)

A function like $(B + 2 I)\, /\, 3$ is between $B()$ and $I()$: it smooths less than $B$, but is not as sensitive to noise as $I$ . This particular blend is called an M-N filter, after
Mitchell and Netravali, Reconstruction Filters in Computer Graphics, 1988, 8p.
M-N splines also have less overshoot than Interpolating splines; see what-is-the-maximum-overshoot-of-interpolating-splines-in-d-dimensions.

A fancier way of using B-splines to interpolate, using IIR filters, is given in
M. Unser, Splines: A perfect fit for signal and image processing , 1999, 17p .