2012-06-10 43 views
6

我在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 
+2

這條曲線究竟「不起作用」? – Junuxx

回答

10

你可能想使用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用於完整發布的數據集。

+3

是的,但這些樣條不保證是單調的,並且可能會超出數據。 –

+0

嘿fraxel,你可以看看這個帖子,看看你是否可以協助: [在一點找到切向量](http://stackoverflow.com/questions/13391449/find-tangent-vector-at-a-point -for-discrete-data-points) – Nader

0

我發現an email來自一位也在尋找Matlab的pchip()函數的本地Python版本的人。他附加his code,不幸的是想下載爲'attachment-001.bin'。如果保存文件並將其重命名爲pychip.py,您會發現它正是您要求的。

+0

我有你提到的pychip.py,但是,當我嘗試它時有一些問題,需要更多的工作。該代碼似乎適用於簡單的曲線,而不適用於上面提供的示例數據。上面的答案雖然運作良好。 – Nader

3

SciPy的確實有pchipscipy.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的分別座標。

+0

看起來現在叫做['pchip_interpolate'](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pchip_interpolate.html)? – endolith

+0

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

+0

@pv。你能解釋一下你的意思嗎?_只是分別插入3個座標中的每一個?你不是說你必須插入兩次嗎?在'x'和'y'之間的第一個插值爲給定的'xi'產生'yi'(如你的例子所示),'x'和'z'之間的第二個插值將產生'zi'(對於相同的' xi')。 – normanius

1

這是另一種解決方案,可以做我想做的(形狀保留)。
問題是,沒有明確的公式或公式來連接點。答案在於創建一個通用於不同點的新數據集。這個新的數據集是沿着路徑的長度(標準)。然後我們使用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