Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

shape-preserving piecewise cubic interpolation for 3D curve in python

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
like image 470
Nader Avatar asked Jun 10 '12 18:06

Nader


1 Answers

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.

like image 79
fraxel Avatar answered Nov 08 '22 14:11

fraxel