Spline to XYZ coordinates

[onslauth]

Hi @msteinbeck,

This is not an issue, just a quick question regarding the conversion of a BSpline using InterpolateCubicNatural to XYZ coordinates.

Looking at the examples, would I use the sample( ) function to get XYZ coordinates to plot the necessary lines to represent the spline?

[msteinbeck]

Hi @onslauth,

yes, sample is the easiest method to plot a spline. It creates a polyline consisting of n points. The default value for n is 100. The points are structured as follows: [x_0, y_0, z_0, x_1, y_1, z_1, ...., x_n-1, y_n-1, z_n-1], given that your spline has dimensionality 3.

[onslauth]

Thanks @msteinbeck for the quick response.

Is there a way to specify a distance between the points when using InterpolateCubicNatural?

At the moment I am trying to find a replacement for SciPy's splprep, and one of the options is to specify a distance between the input points using the u parameter. My use case is that I am receiving coordinates from a gyroscope, where it records the position every X metres apart ( not always uniform ). I then use those values to generate a smooth curve for the path of the gyro.

[msteinbeck]

Is there a way to specify a distance between the points when using InterpolateCubicNatural?

I don't understand this question. The points passed to InterpolateCubicNatural are interpolated. Do you mean the distance between the sampled points? If one of the components of the points to be interpolated is strictly monotonically increasing (I guess X always increases?), then Bisect may be what you are looking for. It let's you extract a point on the spline with respect to a certain component. From the docs:

Tries to find a point P on \p spline such that:

      ts_distance(P[index], value, 1) <= fabs(epsilon)

This function is using the bisection method to determine P. Accordingly, it
is expected that the control points of \p spline are sorted at component
\p index either in ascending order (if \p ascending != 0) or in descending
order (if \p ascending == 0). If the control points of \p spline are not
sorted at component \p index, the behaviour of this function is undefined.
For the sake of fail-safeness, the distance of P[index] and \p value is
compared with the absolute value of \p epsilon (using fabs).

Here you can find an example: https://github.com/msteinbeck/tinyspline/blob/master/examples/cxx/time_series.cpp

Using Bisect, you can sample points with certain distance.

[msteinbeck]

If I understand correctly, you want to evaluate the spline at the knots generated by splprep if parameter u is not supplied. You can create the list of knots as specified by the docs of splprep:

    An array of parameter values. If not given, these values are
    calculated automatically as M = len(x[0]), where

        v[0] = 0
        v[i] = v[i-1] + distance(x[i], x[i-1])
        u[i] = v[i] / v[M-1]

and pass this list to EvalAll.

[onslauth]

Its a bit hard to explain, but I'll try:

The gyro returns 3 values for each point ( distance, inclination, direction ), so a dataset would look like the following:

[ [ 0.0, -69.8, 333.2 ],
  [ 8.0, -69.9, 332,4 ],
  [ 32, -70.0, 331.7 ] ]

I use the 2nd and 3rd column to calculate a new array of shape[ 3, 3 ] which is passed into splprep and I use the first column to pass into splprep as the u parameter.

Also Bisect seems to be exactly what I am looking for, as I would iterate over the start and end values in column 0 in increments of 1 to get my values.

[msteinbeck]

Also Bisect seems to be exactly what I am looking for, as I would iterate over the start and end values in column 0 in increments of 1 to get my values.

Great! But keep in mind that distance must be strictly monotonically increasing. If Bisect solves your issue, please close this issue.

[onslauth]

Bisect looks like what I would use to get the values out, but I am still not sure how to pass in the u parameter as in the case of splprep. I will keep reading and see what I can find.

Thank you for the help.

[msteinbeck]

The DeBoorNet instance returned by Bisect has the property Knot. This is the u corresponding to value passed to Bisect. If you would pass Knot to BSpline::Eval, you would get the same point.

[onslauth]

Sorry to bother you again @msteinbeck,

It seems that I do not quite understand how Bisect works.

I am using the following input data:

0, 0.30820827556315195, -0.15568720276894915, -0.9384930227595559
8, 0.304552451379279, -0.15921617393201867, -0.9390942520947091
32, 0.30114099064220906, -0.16214771720730453, -0.9396926207859083
56, 0.30874984760754304, -0.1546104066159541, -0.9384930227595559
80, 0.3070157384642725, -0.15441282983585378, -0.9390942520947091
104, 0.31170775470673323, -0.1734943106106251, -0.9342044743210294
128, 0.31801356666820046, -0.17555292038329312, -0.9316912275855489
134, 0.31801356666820046, -0.17555292038329312, -0.9316912275855489
140, 0.3142673180673903, -0.18217439236858343, -0.9316912275855489
146, 0.32481798552268165, -0.1662177608727445, -0.9310558158625283
152, 0.3258590935752626, -0.17472411295760046, -0.9291325715340562
158, 0.3225818054076336, -0.18399741330034458, -0.9284858268809135
164, 0.3217070609189397, -0.18874441158804287, -0.92783625389892
170, 0.3237001699665234, -0.17869209405610512, -0.9291325715340562
176, 0.3285306482989197, -0.18663138351276823, -0.9258705848099947
182, 0.3315677418553308, -0.1777850633088033, -0.9265286308718373
188, 0.32112651975692796, -0.19601646482773286, -0.9265286308718373
194, 0.3248303548248393, -0.1898155335032606, -0.9265286308718373
200, 0.3311200562249719, -0.19504461895621233, -0.9232102171129809
206, 0.3235179605155646, -0.20137219509055682, -0.924546033612313
212, 0.33145996880649414, -0.19446640865969766, -0.9232102171129809
218, 0.330270109963688, -0.20878809648798174, -0.9205048534524404
224, 0.3376990686922145, -0.19654555139944607, -0.9205048534524404
230, 0.33716731495056373, -0.20986816996553617, -0.917754625683981
236, 0.33335095174710566, -0.20992229777101828, -0.9191353392552344
242, 0.34322008185265157, -0.20297979057485002, -0.917060074385124
248, 0.33900802171552613, -0.22184063322806827, -0.9142539552342637
254, 0.3493704650143137, -0.20826659596837205, -0.9135454576426009
260, 0.3507394901628449, -0.20908269874760388, -0.9128341772330428
266, 0.35335476498589724, -0.2139992523407996, -0.910683660806177
272, 0.35005506878257614, -0.22819802844612946, -0.9085081775267217
278, 0.35311034472227737, -0.22931253909923735, -0.9070440142914649
284, 0.3544374975981803, -0.23017440225418265, -0.9063077870366499
290, 0.3578958529496517, -0.23064813800916842, -0.9048270524660196
296, 0.36130207380018153, -0.22227451883922802, -0.9055687990111396
302, 0.3632083253832531, -0.22519870286295848, -0.9040825496607784
308, 0.35303960208351315, -0.24082728823698696, -0.9040825496607784
314, 0.3662836238948309, -0.22622137684478769, -0.9025852843498605
320, 0.36961868826904826, -0.2327613103664163, -0.8995577789551803
326, 0.3612731424051445, -0.24829615162922491, -0.898794046299167
332, 0.36818475526080374, -0.24927981980568478, -0.8957117602394129
338, 0.374233897845814, -0.2514747430683212, -0.8925858184521255
344, 0.34024115403206945, -0.22863256718325314, -0.912120116172273
350, 0.29338846280894926, -0.19714886526524067, -0.9354440308298673

I would now like to iterate over 0 to 350 and get the values per increment, i.e. 0, 1, 2, 3, 4.

The first 10 values from SciPy from the above is:

[0]: [0. 0. 0.]
[1]: [ 0.30795933 -0.15593875 -0.93853239]
[2]: [ 0.61542705 -0.31236974 -1.87714317]
[3]: [ 0.92241318 -0.46927646 -2.81583163]
[4]: [ 1.22892857 -0.62664214 -3.75459689]
[5]: [ 1.53498487 -0.78444969 -4.69343783]
[6]: [ 1.84059456 -0.94268174 -5.63235315]
[7]: [ 2.14577093 -1.10132063 -6.57134133]
[8]: [ 2.45052809 -1.26034842 -7.51040065]
[9]: [ 2.75488095 -1.41974685 -8.44952916]
[10]: [ 3.05884524 -1.57949741 -9.38872473]

I am trying to use bisect as follows:

for( var i = 0; i < 350; i++ )
{
    DeBoorNet net = spline.Bisect( i );
    var result    = net.Result;
    Console.WriteLine( "{0,-20} {1,-20} {2,-20} {3,-20}", result[ 0 ], result[ 1 ], result[ 2 ], result[ 3 ] );
}

But the values it provides do not match, and I expect I am using it incorrectly.

Is it possible to use TinySpline to get the same sort of results as from SciPy?

[msteinbeck]

There are a lot of different types of interpolation. TinySpline, currently, supports cubic natural interpolation and Catmull-Rom interpolation (with arbitrary alpha). If I understand correctly, splprep uses some sort of periodic interpolation. Periodic interpolation is primarily used to created closed splines (also known as splinegons). If your spline needs to be closed, neither cubic natural nor Catmull-Rom interpolation is suitable for you. Actually, I know the maths of how to implement periodic interpolation. But, alas, I don't have the time yet to implement it.

You should plot both splines and compare them with each other. You will notice that both methods, splprep and cubic natural, perfectly interpolate your input points. You can also try InterpolateCatmullRom. Maybe this type of interpolation suits you better.

[onslauth]

If I understand what you mean by closed splines, then that is not what I am looking for.

The input values represent a line with a slight curve. But the spline is required to pass through the input points.

[msteinbeck]

The input values represent a line with a slight curve. But the spline is required to pass through the input points.

Cubic natural interpolation creates such a spline. It is just a bit different from the spline interpolated by splprep.

[onslauth]

But how do I get the same values out of TinySpline?

Is the above code correct to call Bisect to get the values at every increment?

[msteinbeck]

But how do I get the same values out of TinySpline?

Well, by implementing the interpolation of splprep in TinySpline :).

Is the above code correct to call Bisect to get the values at every increment?

I need to see how you are calling InterpolateCubicNatural.

[onslauth]

Variable holding the points:

List<double> points = new List<double>();

Data being saved into points as:

points.Add( d );
points.Add( x );
points.Add( y );
points.Add( z );

The data in ( d, x, y, z ) formatted for readability:

0    , 0.30820827556315195 , -0.15568720276894915, -0.9384930227595559
8    , 0.304552451379279   , -0.15921617393201867, -0.9390942520947091
32   , 0.30114099064220906 , -0.16214771720730453, -0.9396926207859083
56   , 0.30874984760754304 , -0.1546104066159541 , -0.9384930227595559
80   , 0.3070157384642725  , -0.15441282983585378, -0.9390942520947091
104  , 0.31170775470673323 , -0.1734943106106251 , -0.9342044743210294
128  , 0.31801356666820046 , -0.17555292038329312, -0.9316912275855489
134  , 0.31801356666820046 , -0.17555292038329312, -0.9316912275855489
140  , 0.3142673180673903  , -0.18217439236858343, -0.9316912275855489
146  , 0.32481798552268165 , -0.1662177608727445 , -0.9310558158625283
...

And I create the spline with the following:

BSpline spline = BSpline.InterpolateCubicNatural( points, 4 );

Currently I am iterating over 0->350 as follows ( which I expect is incorrect ):

for( var i = 0; i < 350; i++ )
{
    DeBoorNet net = spline.Bisect( i );
    var result    = net.Result;
    Console.WriteLine( "{0,-20} {1,-20} {2,-20} {3,-20}", result[ 0 ], result[ 1 ], result[ 2 ], result[ 3 ] );
}

Output from the above:

0                    0.30820827556315195  -0.15568720276894915 -0.9384930227595559
1.000000000000051    0.307509344708701    -0.1564302768067063  -0.9385985913393957
2.000000000000058    0.3069295671226543   -0.15702464896882626 -0.9386891801810344
2.999999999999957    0.30644195959296894  -0.1575039958762714  -0.9387681849846752
4.000000000000136    0.3060104571863283   -0.15791332587451729 -0.9388401392878798
5.000000000000182    0.30561491456397566  -0.15827777971109028 -0.9389075750072239
6                    0.30524413016524626  -0.15861133651220172 -0.9389718998593849

The data seems to be slightly out, just looking at the first line, and comparing that with the SciPY output:

[0]: [0. 0. 0.]
[1]: [ 0.30795933 -0.15593875 -0.93853239]
[2]: [ 0.61542705 -0.31236974 -1.87714317]

It also seems to be that SciPy returns the summed up columns to the point. For example, if we take the data from TinySpline, and sum up the columns, we get a similar output:

0                    0.30820827556315195  -0.15568720276894915 -0.9384930227595559
1.000000000000051    0.307509344708701    -0.1564302768067063  -0.9385985913393957

If we sum up the value of the preceeding lines, we get a similar output to SciPy:

1.000000000000051  0.61571762027  -0.31211747957  -1.8770916141

Is my usage of bisect correct?

Also I tried to insert 4 0s to try to get the same starting point as SciPy, but the result is still out:

0                    0.30820828474847434  -0.15568720740878936 -0.9384930507288276

[msteinbeck]

Is my usage of bisect correct?

From what I see, your usage of Bisect looks good.

The data seems to be slightly out, just looking at the first line, and comparing that with the SciPY output:

[0]: [0. 0. 0.] [1]: [ 0.30795933 -0.15593875 -0.93853239] [2]: [ 0.61542705 -0.31236974 -1.87714317]

I'm wondering why SciPY interpolates 2 as [ 0.61542705 -0.31236974 -1.87714317]. You input data says:

0    , 0.30820827556315195 , -0.15568720276894915, -0.9384930227595559
8    , 0.304552451379279   , -0.15921617393201867, -0.9390942520947091

Based on 0 -> 0.30820827556315195 and 8 -> 0.304552451379279 there is no reason for the interpolation to map 2 -> 0.61542705. The same is true for the remaining components.

[onslauth]

SciPy's interpolation of 2 includes the previous entries of 0 and 1 added to the value at 2.

SciPy seems to return a rolling total for every interpolation, which makes it very easy to take the data and draw 3D constructs from it, as you don't have to sum up all the points between two interpolated values.

Can you suggest any further changes that I could make to the code to get the results closer to the returned values from SciPy? I would like to if possible replicate the data exactly.

[msteinbeck]

as you don't have to sum up all the points between two interpolated values.

This is not what I would call interpolation. Anyhow, SciPy is a mature library, so they know what they are doing :).

would like to if possible replicate the data exactly.

I doubt that. SciPy is using a special interpolation which is not implemented in TinySpline (yet). If you could provide the implementation details, I could have a look. As always, PRs are welcome. For now, you can only sum up the entries by your own. That said, the results will slightly differ as cubic natural spline interpolation has other objects than periodic interpolation.

[msteinbeck]

I'm closing this issue. Don't hesitate to reopen it if necessary.