justuswolff
justuswolff

Reputation: 61

Get Trajectory of Three Dimensional Cubic Spline Scipy

I am trying to approximate a given route (coordinates) with a three dimensional cubic spline. Example data:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import numpy as np

%matplotlib inline

x = np.array([1, 2, 2.3, 3, 4, 5, 5.5, 8, 9, 9.5])
y = np.arange(0, 10)
z = np.sin(x) * np.cos(y^2) + x
fig = plt.figure(figsize=(10,6))
ax = axes3d.Axes3D(fig)
ax.stem(x, y, z)

enter image description here

I now approximate the data points with RBF Interpolator from scipy. Is this the right approach?

from scipy.interpolate import RBFInterpolator
coord_data = np.stack([x, y], -1)
spline = RBFInterpolator(coord_data, z, kernel = 'cubic')

How do I get the resulting spline now (the points it is following through)? And how do I access it's derivatives?

Upvotes: 1

Views: 1138

Answers (1)

Iddo Hanniel
Iddo Hanniel

Reputation: 2066

You are trying to approximate a route, i.e., a curve in 3D, and not a surface. The approach you are trying results in a surface and is therefore not suitable for your case.

A suitable representation for a 3D curve is in parametric form as a tuple (x(u), y(u), z(u)) where u is some parameter and each coordinate is a function of u. The curve fitting problem is reduced to three 2D problems of fitting (ui, xi), (ui, yi), and (ui, zi) separately (see also my answer here).

So, in order to perform a curve fit you are required to provide the parameterization ui for each input point. A common parameterization in spline fitting is the chord-length parameterization. This parameterization is defined by the accumulated length of the distances between the ordered points (u0=0, u1=|p1-p0|, u2 = u1+|p2-p1|... etc.).

The following code implements a spline interpolation of your data with the chord-length parameterization.

from scipy import interpolate

xyz = np.vstack([x, y, z]).T
u = np.cumsum(np.r_[[0], np.linalg.norm(np.diff(xyz, axis=0), axis=1)])
# u is the chord-legth parameterization for each xyz point

sx = interpolate.InterpolatedUnivariateSpline(u, x)  # x(u) spline
sy = interpolate.InterpolatedUnivariateSpline(u, y)  # y(u) spline
sz = interpolate.InterpolatedUnivariateSpline(u, z)  # z(u) spline

The code below samples the resulting spline and plots the result (in blue) on your data (black polyline). The result looks like this:

enter image description here

uu = np.linspace(u[0], u[-1], 100)
xx = sx(uu)
yy = sy(uu)
zz = sz(uu)
plt.plot(xx, yy, zz, "b")

You can also use the RBF functions you suggested, for univariate interpolation. The following code is an example of how to do this:

from scipy.interpolate import Rbf

rbfi_x = Rbf(u, x, function='cubic')
rbfi_y = Rbf(u, y, function='cubic')
rbfi_z = Rbf(u, z, function='cubic')

Sampling the resulting functions in a similar way to the spline sample above, and plotting the results on the previous figure, we get the plot below. As can be seen the spline and RBF interpolations are similar but not the same (for example near the endpoints).

enter image description here

Upvotes: 3

Related Questions