7

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

以下是数据点： 

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
4

3 回答 3

11

您可能想像这样使用splprep() 和 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

这既漂亮又流畅，这也适用于您完整发布的数据集。

于 2012-06-10T22:31:55.507 回答
4

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-11T12:03:08.320 回答
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

于 2012-06-24T19:24:37.870 回答