Shape preserving piecewise cubic interpolation for a 3D curve in python

Question

I have a curve in three-dimensional space. I want to use piecewise cubic interpolation that preserves a shape on it, similar to pchip in matlab. I researched the functions provided in scipy.interpolate, for example. interp2d, but functions work for some curve structures, not for the data points that I have. Any ideas on how to do this?

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 

Answer 1

You probably want to use splprep () and splev () like this (basic explanation 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() 

It produces:

Which is nice and smooth, also great for your complete dataset.

Answer 2

Scipy has pchip in scipy.interpolate --- but, well, 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 3D data, simply interpolate each of the three coordinates separately.

Answer 3

Here is another solution that does what I wanted (saving the form).
The problem is that there is no clear formula or equation for connecting the dots. The answer is to create a new dataset that is common to different points. This new dataset is the path length (normal). Then we 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 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 in the figure below. the blue line is the curve fit data, and the red dots are the original dataset. One thing that I noticed with this approach is that if the dataset is long, then interpolating .interp1d is inefficient and takes a lot of time.

Answer 4

I found an email from a person who was also looking for a version of the Matlab pchip () function in the native version of Python. He added his code , which, unfortunately, he wants to download as "attachment-001.bin". If you save the file and rename it to the pychip.py file, you will find that this is exactly what you are asking for.