Cubic spline extrapolation
Asked Answered
N

5

7

I have a nice cubic spline code but it is for interpolation only. I need to extrapolate just a little into the future. Does anyone know of a good source of code, not a library, for doing this?

This is the code I wrote in basic (now ASM) for interpolation.

Nanananak answered 19/5, 2009 at 5:52 Comment(0)
A
6

You don't need new code for that.

To extrapolate the spline you can extrapolate the parameters of the first and last spline.

Depending on your existing code/library that might not be possible without modifying the code. In that case just prepend/append two other points to the beginning/end of your list of points. You can get those two points by linearily interpolating between the first/last two points.

Be careful: Depending on the original meaning of the points that extrapolation might be completely inappropriate, especially when it comes to statistical data. In that case you should consider using regression analysis.

Antispasmodic answered 19/5, 2009 at 6:19 Comment(2)
>You can get those two points by linearily interpolating between the first/last two points. That defeats the purpose of the exercise. I want to extend based on the CURVE produced by the spline.Nanananak
As I said, that would be the less-than-optimal fallback if your existing code/library does not allow you to extrapolate (i.e. use a parameter < 0 or > 1)Antispasmodic
T
4

For simplicity, I'm going to represent a cubic Bezier curve as 4 points (A, B, C, D), where A and D are the endpoints of the curve, and B and C are the "control handle points". (The actual curve usually does not touch the control handle points).

See "Don Lancaster's Guru's Lair Cubic Spline Library" for ways to convert this representation of a cubic Bezier curve into other popular representations.

interpolation

Given one cubic Bezier curve (P0, P1, P2, P3), we use De Casteljau's algorithm to chop a Bezier curve into a left half and a right half. This is super-easy even on a microcontroller that doesn't have a "multiply" instruction, because it only requires calculating a few averages until we get a midpoint:

P0
   F0 := average(P0, P1)
P1                       S0 := average(F0, F1)
   F1 := average(P1, P2)         Midpoint := average(S0, S1)
P2                       S1 := average(F1, F2)
   F2 := average(P2, P3)
P3

The entire Bezier curve is (P0, P1, P2, P3).

The left half of that entire Bezier curve is the Bezier curve (P0, F0, S0, M).

The right half of that entire Bezier curve is the Bezier curve (M, S1, F2, P3).

Many microcontrollers continue to divide each curve up into smaller and smaller little curves until each piece is small enough to approximate with a straight line.

But we want to go the other way -- extrapolate out to a bigger curve.

extrapolation

Given either the left half or the right half, we can run this in reverse to recover the original curve.

Let's imagine that we forgot the original points P1, P2, P3.

Given the left half of a Bezier curve (P0, F0, S0, M), we can extrapolate to the right with:

S1 := M + (M - S0)
F1 := S0 + (S0 - F0)
P1 := F0 + (F0 - P0)

then use those values to calculate

F2 := S1 + (S1 - F1)
P2 := F1 + (F1 - P1)

and finally

P3 := F2 + (F2 - P2)

to extrapolate and recover the extrapolated Bazier curve (P0, P1, P2, P3).

details

The extrapolated curve (P0, P1, P2, P3) passes through every point in the original curve (P0, F0, S0, M) -- in particular, starting at P0 and passing through the midpoint M -- and keeps on going until it reaches P3.

We can always extrapolate from any 4 points (P0, F0, S0, M), whether or not those 4 points were originally calculated as the left half (or right half) of some larger Bezier spline.

I'm sure you already know this, but just for clarity:

Midpoint = average(F0, F1)

means "find the Midpoint exactly halfway between points F0 and F1", or in other words,

Midpoint.x = (F0.x + F1.x)/2
Midpoint.y = (F0.y + F1.y)/2
Midpoint.z = (F0.z + F1.z)/2

The expression

S1 := M + (M - S0)

means "Given a line segment, with one end at S0, and the midpoint at M, start at S0 and run in a straight line past M until you reach the other end at S1", or in other words (unless you have a decent vector library) 3 lines of code

S1.x := M.x + (M.x - S0.x)
S1.y := M.y + (M.y - S0.y)
S1.z := M.z + (M.z - S0.z)

. (If you're doing 2D, skip all the "z" stuff -- it's always zero).

Trishtrisha answered 22/6, 2011 at 3:32 Comment(0)
W
2

You'll really have to expand that question a little. Also, "cubic spline" is a very wide term.

If you're interested in splines, I can heartly reccomend Carl de Boors "A Practical Guide to Splines". It is however a little mathematically oriented, but it has code examples included (they can be downloaded from the author's home page). Googling and wikiing for "cubic spline" can bring up some examples, maybe even in particular languages - another thing to add to the question (if you're looking for code).

If you're interested in extrapolation and curve fitting, googling those could help. Matlab package has a very nice curve fitting toolbox. Wikipedia has some links to useful references

Really, it is too wide a question, to even start guessing an answer.

Also, could you explain what exactly are you trying to do ? What kind of data ? Anything ?


Edit1: Here, try this: you may find something useful in here - link

Wifehood answered 19/5, 2009 at 6:17 Comment(4)
MATLAB is a library, im looking for code as mentioned. I have a seriesof data points producing a curve. I need to extend the curve foreward just a little. Its no more complicated than that.Nanananak
Well, just take the last 3rd order curve you've got and calculate the desired values then. If you've fitted it through, you've already got the coefficients.Wifehood
This code was written many many years ago when I could do this level of math easily. Now I would have to tear the whole thing apart and brush up on the concepts... non trivial. I was hoping this post would prompt a link to some code...Nanananak
But you have to know something about the algorithm used, and the code you wrote. After all, this is trivial if you wrote the spline in the first place. I cannot give you any help regarding code, because I don't know what you did. As I said in my post, "cubic spline" is not an exact description, really.Wifehood
Q
2

You need to write better requirements for requested code. Splines are usually used for interpolation of some unknown or complex function by using of some fixed data set. If you want to have an estimate of function's value outside of boundaries of this data set then you shouldn't use splines.

If your spline is function defined in the place where you really want to evaluate your value (cubic, but not piecewise-cubic) then you already can evaluate that value.

If you want to have ability to evaluate your spline outside of interpolation range, but leave it as piecewise-cubic function with the same values inside of interpolation range then you should extend spline range by some nodes, and add some logic of evaluation values at the new nodes (for example you want to have your spline be not only a continuous function, but also have some number of first derivatives be also continuous functions)

Really I suggest you to use some algorithm more suitable for extrapolation, like usage of Lagrange polynomial if everything you really need is single value not very far from points of original data set.

Quietism answered 19/5, 2009 at 6:33 Comment(1)
>you shouldn't use splines. Why not? The 3rd order polynominal is already calculated. I just need to use the"curvature" and extend foreward slightly. I will study the Lagrange polynominalNanananak
C
1

Generally for spline interpolation you use a variable t to interpolate over the line. As long as 0 <= t <= 1 you're interpolating. However, when t < 0 or t > 1 you're simply extrapolating the spline.

Coppage answered 19/5, 2009 at 6:18 Comment(2)
yes that is correct. I need datapoints for t > 1 in this example (but only slightly > tNanananak
Just evaluate the spline formula with t > 1, for small t you should be fine.Coppage

© 2022 - 2024 — McMap. All rights reserved.