Вопрос по python, numpy, curve-fitting, scipy – сохраняющая форму кусочно-кубическая интерполяция для трехмерной кривой в питоне

6

У меня есть кривая в 3D-пространстве. Я хочу использовать сохраняющую форму кусочно-кубическую интерполяцию, похожую на pchip в matlab. Я исследовал функции, предоставляемые в 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
Что именно "не работает" для этой кривой? Junuxx

Ваш Ответ

4   ответа
3

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 координат отдельно.

pchip это просто ярлык дляPchipInterpolator normanius
@ ФВ. Не могли бы вы объяснить, что вы подразумеваете под Просто интерполируйте каждую из 3 координат отдельно? Разве вы не имеете в виду, что вы должны интерполировать только дважды? Первая интерполяция междуx а такжеy даетyi для данногоxi (как показано в вашем примере) и вторая интерполяция междуx а такжеz дастzi (для того жеxi). normanius
Похоже, теперь это называетсяpchip_interpolate? endolith
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 неэффективен и занимает много времени.

10

Вы, наверное, хотите использовать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()

Produces:

Что хорошо и гладко, это также хорошо работает для вашего полного опубликованного набора данных.

Да, но эти сплайны не обязательно являются монотонными и могут не соответствовать данным. pv.
Эй, Фраксель, можешь посмотреть на этот пост, чтобы понять, можешь ли ты помочь: найти касательный вектор в точке Nader
0

электронное письм от парня, который также искал версию Pychip () Matlab для языка Python. Он прилагается его код который, к сожалению, хочет скачать как 'attachment-001.bin'. Если вы сохраните файл и переименуете его в pychip.py, вы обнаружите, что это именно то, что вы просите.

У меня был pychip.py, о котором вы упомянули, однако, когда я попробовал, у него были некоторые проблемы и мне нужно больше работать. Код, кажется, работает для простых кривых, а не для данных примера, которые я предоставил выше. Ответы выше работали хорошо, хотя. Nader

Похожие вопросы