7

I have a curve in 3D space. I want to use a shape-preserving piecewise cubic interpolation on it similar to pchip in matlab. I researched functions provided in scipy.interpolate, e.g. interp2d, but the functions work for some curve structures and not the data points I have. Any ideas of how to do it?

Here are the data points: 

x,y,z
0,0,0
0,0,98.43
0,0,196.85
0,0,295.28
0,0,393.7
0,0,492.13
0,0,590.55
0,0,656.17
0,0,688.98
0,0,787.4
0,0,885.83
0,0,984.25
0,0,1082.68
0,0,1181.1
0,0,1227.3
0,0,1279.53
0,0,1377.95
0,0,1476.38
0,0,1574.8
0,0,1673.23
0,0,1771.65
0,0,1870.08
0,0,1968.5
0,0,2066.93
0,0,2158.79
0,0,2165.35
0,0,2263.78
0,0,2362.2
0,0,2460.63
0,0,2559.06
0,0,2647.64
-0.016,0.014,2657.48
-1.926,1.744,2755.868
-7.014,6.351,2854.041
-15.274,13.83,2951.83
-26.685,24.163,3049.031
-41.231,37.333,3145.477
-58.879,53.314,3240.966
-79.6,72.076,3335.335
-103.351,93.581,3428.386
-130.09,117.793,3519.96
-159.761,144.66,3609.864
-192.315,174.136,3697.945
-227.682,206.16,3784.018
-254.441,230.39,3843.779
-265.686,240.572,3868.036
-304.369,275.598,3951.483
-343.055,310.627,4034.938
-358.167,324.311,4067.538
-381.737,345.653,4118.384
-420.424,380.683,4201.84
-459.106,415.708,4285.286
-497.793,450.738,4368.741
-505.543,457.756,4385.461
-509.077,460.955,4393.084
-536.475,485.764,4452.188
-575.162,520.793,4535.643
-613.844,555.819,4619.09
-624.393,565.371,4641.847
-652.22,591.897,4702.235
-689.427,631.754,4784.174
-725.343,675.459,4864.702
-759.909,722.939,4943.682
-793.051,774.087,5020.95
-809.609,801.943,5060.188
-820.151,820.202,5085.314
-824.889,828.407,5096.606
-830.696,838.466,5110.448
-846.896,867.72,5150.399
-855.384,883.717,5172.081
-884.958,939.456,5247.626
-914.53,995.188,5323.163
-944.104,1050.927,5398.708
-973.675,1106.659,5474.246
-1003.249,1162.398,5549.791
-1032.821,1218.13,5625.328
-1062.395,1273.869,5700.873
-1091.966,1329.601,5776.411
-1121.54,1385.34,5851.956
-1151.112,1441.072,5927.493
-1180.686,1496.811,6003.038
-1210.257,1552.543,6078.576
-1239.831,1608.282,6154.121
-1269.403,1664.015,6229.658
-1281.875,1687.521,6261.517
-1298.67,1720.451,6304.797
-1317.209,1760.009,6353.528
-1326.229,1780.608,6377.639
-1351.851,1844.711,6447.786
-1375.462,1912.567,6515.035
-1379.125,1923.997,6525.735
-1397.002,1984.002,6579.217
-1416.406,2058.808,6640.141
-1433.629,2136.794,6697.655
-1448.619,2217.744,6751.587
-1453.008,2244.679,6768.334
-1461.337,2301.426,6801.784
-1471.745,2387.628,6848.122
-1479.813,2476.093,6890.468
-1485.519,2566.597,6928.713
-1488.852,2658.874,6962.744
-1489.801,2752.688,6992.481
-1488.358,2847.765,7017.828
-1484.534,2943.865,7038.72
-1478.344,3040.704,7055.099
-1469.806,3137.966,7066.915
-1469.799,3138.035,7066.922
-1458.925,3235.574,7074.155
-1445.742,3333.07,7076.775
-1444.753,3339.757,7076.785
-1438.72,3380.321,7076.785
-1431.268,3430.42,7076.785
-1416.787,3527.779,7076.785
-1402.308,3625.128,7076.785
-1401.554,3630.192,7076.785
-1387.827,3722.487,7076.785
-1373.347,3819.836,7076.785
-1358.866,3917.195,7076.785
-1357.872,3923.882,7076.785
-1353.32,3954.485,7076.785
-1344.387,4014.544,7076.785
-1329.906,4111.903,7076.785
-1315.427,4209.252,7076.785
-1300.946,4306.611,7076.785
-1286.466,4403.96,7076.785
-1271.985,4501.319,7076.785
-1257.504,4598.678,7076.785
-1243.025,4696.027,7076.785
-1228.544,4793.386,7076.785
-1214.065,4890.735,7076.785
-1199.584,4988.094,7076.785
-1185.104,5085.443,7076.785
-1170.623,5182.802,7076.785
-1156.144,5280.151,7076.785
-1141.663,5377.51,7076.785
-1127.183,5474.859,7076.785
-1112.703,5572.218,7076.785
-1098.223,5669.567,7076.785
-1083.742,5766.926,7076.785
-1069.263,5864.275,7076.785
-1054.782,5961.634,7076.785
-1040.302,6058.983,7076.785
-1025.821,6156.342,7076.785
-1011.342,6253.692,7076.785
-996.861,6351.05,7076.785
-982.382,6448.4,7076.785
-967.901,6545.759,7076.785
-953.421,6643.108,7076.785
-938.94,6740.467,7076.785
-924.461,6837.816,7076.785
-909.98,6935.175,7076.785
-895.499,7032.534,7076.785
-895.234,7034.314,7076.785
-883.075,7130.158,7076.785
-874.925,7228.243,7076.785
-871.062,7326.579,7076.785
-871.491,7425,7076.785
-876.213,7523.299,7076.785
-885.218,7621.308,7076.785
-898.489,7718.822,7076.785
-916.003,7815.673,7076.785
-937.722,7911.659,7076.785
-963.61,8006.615,7076.785
-993.613,8100.342,7076.785
-1027.678,8192.681,7076.785
-1065.735,8283.437,7076.785
-1083.912,8323.221,7076.785
-1107.12,8372.742,7076.785
-1148.885,8461.861,7076.785
-1190.655,8550.989,7076.785
-1232.42,8640.108,7076.785
-1274.19,8729.236,7076.785
-1315.955,8818.354,7076.785
-1357.724,8907.482,7076.785
-1399.49,8996.601,7076.785
-1441.259,9085.729,7076.785
-1483.024,9174.848,7076.785
-1486.296,9181.829,7076.785
-1510.499,9231.861,7076.785
-1526.189,9263.304,7076.785
-1570.139,9351.377,7076.785
-1614.085,9439.441,7076.785
-1658.035,9527.514,7076.785
-1701.98,9615.578,7076.785
-1745.93,9703.651,7076.785
-1789.876,9791.715,7076.785
-1833.826,9879.788,7076.785
-1877.771,9967.852,7076.785
-1921.721,10055.925,7076.785
-1965.667,10143.989,7076.785
-2009.625,10232.059,7076.785
-2053.585,10320.115,7076.785
-2097.551,10408.18,7076.785
-2141.512,10496.237,7076.785
-2185.477,10584.302,7076.785
-2229.438,10672.359,7076.785
-2273.403,10760.424,7076.785
-2317.364,10848.481,7076.785
-2352.213,10918.285,7076.785
Nader
  • 660
  • 1
  • 11
  • 17

3 Answers3

11

You probably want to use splprep() and splev(), like this (basic explaination in comments):

import scipy
from scipy import interpolate
import numpy as np

#This is your data, but we're 'zooming' into just 5 data points
#because it'll provide a better visually illustration
#also we need to transpose to get the data in the right format
data = data[100:105].transpose()

#now we get all the knots and info about the interpolated spline
tck, u= interpolate.splprep(data)
#here we generate the new interpolated dataset, 
#increase the resolution by increasing the spacing, 500 in this example
new = interpolate.splev(np.linspace(0,1,500), tck)

#now lets plot it!
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(data[0], data[1], data[2], label='originalpoints', lw =2, c='Dodgerblue')
ax.plot(new[0], new[1], new[2], label='fit', lw =2, c='red')
ax.legend()
plt.savefig('junk.png')
plt.show()

Produces:

enter image description here

Which is nice and smooth, also this works fine for your full posted dataset.

fraxel
  • 34,470
  • 11
  • 98
  • 102
  • 3
    Yes, but these splines are not guaranteed to be monotonic, and can overshoot the data. – pv. Jun 11 '12 at 12:04
  • Hey fraxel, can you look at this post to see if you can assist: [find tangent vector at a point](http://stackoverflow.com/questions/13391449/find-tangent-vector-at-a-point-for-discrete-data-points) – Nader Nov 15 '12 at 04:20
4

Scipy does have pchip in scipy.interpolate --- but, uhh, someone forgot to list it in the documentation:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import pchip
x = np.linspace(0, 1, 20)
y = np.random.rand(20)
xi = np.linspace(0, 1, 200)
yi = pchip(x, y)(xi)
plt.plot(x, y, '.', xi, yi)

For 3-D data, just interpolate each of the 3 coordinates separately.

pv.
  • 33,875
  • 8
  • 55
  • 49
  • Looks like this is now called [`pchip_interpolate`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pchip_interpolate.html)? – endolith Mar 05 '14 at 14:39
  • `pchip` is just a shortcut for [`PchipInterpolator`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PchipInterpolator.html) – normanius Dec 29 '17 at 02:04
  • @pv. Could you please explain what you mean by _just interpolate each of the 3 coordinates separately_? Din't you mean that you have to interpolate just twice? A first interpolation between `x` and `y` yields `yi` for given `xi` (as shown in your example), and a second interpolation between `x` and `z` will yield `zi` (for the same `xi`). – normanius Dec 29 '17 at 02:15
1

Here is another solution that does what I wanted (shape-preserving).
The problem is that there is no clear formula or equation to connect the points. The answer lies in creating a new data set that is common to the different points. This new data set is the length along the path (norm). We then use interp1 to interpolate each set. 

import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# read data from a file
# as x, y , z

# create a new array to hold the norm (distance along path)
npts = len(x)
s = np.zeros(npts, dtype=float)
for j in range(1, npts):
    dx = x[j] - x[j-1]
    dy = y[j] - y[j-1]
    dz = z[j] - z[j-1]
    vec = np.array([dx, dy, dz])
    s[j] = s[j-1] + np.linalg.norm(vec)

# create a new data with finer coords 
xvec = np.linspace(s[0], s[-1], 5000)

# Call the Scipy cubic spline interpolator
from scipy.interpolate import interpolate

# Create new interpolation function for each axis against the norm 
f1 = interpolate.interp1d(s, x, kind='cubic')
f2 = interpolate.interp1d(s, y, kind='cubic')
f3 = interpolate.interp1d(s, z, kind='cubic')

# create new fine data curve based on xvec
xs = f1(xvec)
ys = f2(xvec)
zs = f3(xvec)

# now let's plot to compare
#1st, reverse z-axis for plotting
z = z*-1 
zs = zs*-1

dpi = 75
fig = plt.figure(dpi=dpi, facecolor = '#617f8a')
threeDPlot = fig.add_subplot(111, projection='3d')
fig.subplots_adjust(left=0.03, bottom=0.02, right=0.97, top=0.98)
mpl.rcParams['legend.fontsize'] = 10

threeDPlot.scatter(x, y, z, color='r')  # Original data as a scatter plot
threeDPlot.plot3D(xs, ys, zs, label='Curve Fit', color='b', linewidth=1)
threeDPlot.legend()
plt.show()

The result is shown int he figure below. the blue line is the curve fit data while the red dots are the original data set. One thing I noticed with this approach though, is that if the data set is long, then interpolate.interp1d is not efficient and takes long time.

curve fit interpolation

Nader
  • 660
  • 1
  • 11
  • 17