python - Pythonでの3D曲線の形状保存区分的3次補間

Question

3D 空間に曲線があります。matlab の pchip と同様に、形状を維持する区分的 3 次補間を使用したいと考えています。interp2d などの scipy.interpolate で提供される関数を調査しましたが、関数は一部の曲線構造では機能し、私が持っているデータ ポイントでは機能しません。それを行う方法のアイデアはありますか？

データポイントは次のとおりです。

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

Answer 1

おそらく、次のように 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()

プロデュース:

これは素晴らしくスムーズです。これは、投稿された完全なデータセットに対してもうまく機能します。

Answer 2

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)

3D データの場合は、3 つの座標をそれぞれ個別に補間するだけです。

Answer 3

これは、私が望んでいたことを行う別のソリューションです（形状を維持します）。
問題は、点を結ぶ明確な公式や方程式がないことです。その答えは、さまざまなポイントに共通する新しいデータ セットを作成することにあります。この新しいデータ セットは、パス (ノルム) に沿った長さです。次に、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()

その結果を下の図に示します。青い線はカーブ フィット データで、赤い点は元のデータ セットです。ただし、このアプローチで気づいたことの 1 つは、データ セットが長い場合、interpolate.interp1d は効率的でなく、時間がかかることです。