[−][src]Module resid::spline

Our objective is to construct a smooth interpolating single-valued function y = f(x).

Catmull-Rom splines are widely used for interpolation, however these are parametric curves [x(t) y(t) ...] and can not be used to directly calculate y = f(x). For a discussion of Catmull-Rom splines see Catmull, E., and R. Rom, "A Class of Local Interpolating Splines", Computer Aided Geometric Design.

Natural cubic splines are single-valued functions, and have been used in several applications e.g. to specify gamma curves for image display. These splines do not afford local control, and a set of linear equations including all interpolation points must be solved before any point on the curve can be calculated. The lack of local control makes the splines more difficult to handle than e.g. Catmull-Rom splines, and real-time interpolation of a stream of data points is not possible. For a discussion of natural cubic splines, see e.g. Kreyszig, E., "Advanced Engineering Mathematics".

Our approach is to approximate the properties of Catmull-Rom splines for piecewice cubic polynomials f(x) = ax^3 + bx^2 + cx + d as follows: Each curve segment is specified by four interpolation points, p0, p1, p2, p3. The curve between p1 and p2 must interpolate both p1 and p2, and in addition f'(p1.x) = k1 = (p2.y - p0.y)/(p2.x - p0.x) and f'(p2.x) = k2 = (p3.y - p1.y)/(p3.x - p1.x).

The constraints are expressed by the following system of linear equations

ⓘThis example is not tested
[ 1  xi    xi^2    xi^3 ]   [ d ]    [ yi ]
[     1  2*xi    3*xi^2 ] * [ c ] =  [ ki ]
[ 1  xj    xj^2    xj^3 ]   [ b ]    [ yj ]
[     1  2*xj    3*xj^2 ]   [ a ]    [ kj ]

Solving using Gaussian elimination and back substitution, setting dy = yj - yi, dx = xj - xi, we get

ⓘThis example is not tested
a = ((ki + kj) - 2*dy/dx)/(dx*dx);
b = ((kj - ki)/dx - 3*(xi + xj)*a)/2;
c = ki - (3*xi*a + 2*b)*xi;
d = yi - ((xi*a + b)*xi + c)*xi;

Having calculated the coefficients of the cubic polynomial we have the choice of evaluation by brute force

ⓘThis example is not tested
for (x = x1; x <= x2; x += res) {
  y = ((a*x + b)*x + c)*x + d;
  plot(x, y);
}

or by forward differencing

ⓘThis example is not tested
y = ((a*x1 + b)*x1 + c)*x1 + d;
dy = (3*a*(x1 + res) + 2*b)*x1*res + ((a*res + b)*res + c)*res;
d2y = (6*a*(x1 + res) + 2*b)*res*res;
d3y = 6*a*res*res*res;

for (x = x1; x <= x2; x += res) {
  plot(x, y);
  y += dy; dy += d2y; d2y += d3y;
}

See Foley, Van Dam, Feiner, Hughes, "Computer Graphics, Principles and Practice" for a discussion of forward differencing.

If we have a set of interpolation points p0, ..., pn, we may specify curve segments between p0 and p1, and between pn-1 and pn by using the following constraints: f''(p0.x) = 0 and f''(pn.x) = 0.

Substituting the results for a and b in

2b + 6a*xi = 0

we get

ki = (3*dy/dx - kj)/2;

or by substituting the results for a and b in

2b + 6a*xj = 0

we get

kj = (3*dy/dx - ki)/2;

Finally, if we have only two interpolation points, the cubic polynomial will degenerate to a straight line if we set

ki = kj = dy/dx;

Structs

Point
PointPlotter

Functions

interpolate

Evaluation of complete interpolating function. Note that since each curve segment is controlled by four points, the end points will not be interpolated. If extra control points are not desirable, the end points can simply be repeated to ensure interpolation. Note also that points of non-differentiability and discontinuity can be introduced by repeating points.