蟒蛇3D形状保持分段立方插值

我在3D空间中有一条曲线。我想在matlab中使用类似于pchip的形状保持分段三次插值。我研究了在scipy.interpolate中提供的函数，例如interp2d，但这些函数适用于某些曲线结构，而不适用于我拥有的数据点。任何想法如何做到这一点？蟒蛇3D形状保持分段立方插值

这里有个数据点：

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

来源

2012-06-10 Nader

这条曲线究竟“不起作用”？ – Junuxx

你可能想使用splprep() and splev()，像这样（在评论基本交代）：

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()

产地：

enter image description here

这是很好，顺利，这也可以工作e用于完整发布的数据集。

来源

2012-06-10 22:31:55 fraxel

是的，但这些样条不保证是单调的，并且可能会超出数据。 –

嘿fraxel，你可以看看这个帖子，看看你是否可以协助： [在一点找到切向量]（http://stackoverflow.com/questions/13391449/find-tangent-vector-at-a-point -for-discrete-data-points） – Nader

我发现an email来自一位也在寻找Matlab的pchip（）函数的本地Python版本的人。他附加his code，不幸的是想下载为'attachment-001.bin'。如果保存文件并将其重命名为pychip.py，您会发现它正是您要求的。

来源

2012-06-10 23:36:48

我有你提到的pychip.py，但是，当我尝试它时有一些问题，需要更多的工作。该代码似乎适用于简单的曲线，而不适用于上面提供的示例数据。上面的答案虽然运作良好。 – Nader

SciPy的确实有pchip在scipy.interpolate ---但是，呃，有人忘了它在文档中：

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)

对于3-d的数据，只是插每个3的分别坐标。

来源

2012-06-11 12:03:08

看起来现在叫做['pchip_interpolate']（http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pchip_interpolate.html）？ – endolith

'pchip'只是['PchipInterpolator']的快捷键（https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PchipInterpolator.html） – normanius

@pv。你能解释一下你的意思吗？_只是分别插入3个坐标中的每一个？你不是说你必须插入两次吗？在'x'和'y'之间的第一个插值为给定的'xi'产生'yi'（如你的例子所示），'x'和'z'之间的第二个插值将产生'zi'（对于相同的' xi'）。 – normanius

这是另一种解决方案，可以做我想做的（形状保留）。
问题是，没有明确的公式或公式来连接点。答案在于创建一个通用于不同点的新数据集。这个新的数据集是沿着路径的长度（标准）。然后我们使用interp1插入每个集合。

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()

结果如下图所示。蓝线是曲线拟合数据，而红点是原始数据集。但是，我注意到使用这种方法的一件事是，如果数据集很长，那么interpolate.interp1d效率不高，需要很长时间。

curve fit interpolation

来源

2012-06-24 19:24:37 Nader

蟒蛇3D形状保持分段立方插值

回答

相关问题