scipy 获取插值平面(CloughTocher2DInterpolator)与平面之间的直线交点

m0rkklqb 于 2023-04-30 发布在其他

关注(0)|答案(1)|浏览(237)

我有一些点，我想从其中插入一个平面，我做如下：

rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)
X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)

在这个例子中，我定义了x，y和z为随机的。x、y和z也可以是稀疏的。
现在，我希望能够从插值中得到一个与我定义为x*a + y*b + c = z的平面的交集（带有a，b和c变量）
我在scipy里找不到任何可以做这个的东西。我已经尝试得到所有的点，得到一个三角测量，并使用这个相交（见：How to visualize 3D delaunay triangulation in Python?）但是我得到的三角测量真的很糟糕。它根本没有连接到最近的点。
我也试过使用curve_fit，但没有任何运气，因为当你请求船体之外的值时，插值器可能会返回nan值，然后优化器似乎卡住了：

from scipy import optimize

def func(x, y):
    z = uslope[0] * x + uslope[1] * y + uslope[2]
    dist = interp(x, y) - z
    dist[np.isnan(dist)] = np.Infinity
    return dist

xdata = np.linspace(0, 1000, 50)
ydata = np.zeros(xdata.shape)
ydata.fill(500)

popt, pcov = optimize.curve_fit(func, xdata, ydata, p0=500, check_finite=True, method="dogbox")
print(popt)
print(pcov)

如果CloughTocher2DInterpolator复杂，我也可以使用另一个连续的2D插值器。

scipy

来源：https://stackoverflow.com/questions/76042340/get-line-intersection-between-interpolated-plane-cloughtocher2dinterpolator

1条答案

按热度按时间

jv4diomz1#

可以使用plt.contour函数查找折线的xy坐标，该折线近似于插值曲面与平面的相交曲线。
下面的代码从您的示例开始，然后查找函数Z_interp - Z_plane的零轮廓，它给出相交曲线的xy多段线（在一般情况下可以有多个曲线分支）。

np.random.seed(7)  # so you can reconstruct my results :)
x = np.random.random(10) - 0.5
y = np.random.random(10) - 0.5
z = np.hypot(x, y)
X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)
Z = interp(X, Y)

plt.pcolormesh(X, Y, Z, shading='auto')
plt.plot(x, y, "ok", label="input point")
plt.axis("equal")

(up这里只是你的代码。.）

# plane coefficients
a = 0
b = 0
c = 0.03

Z_plane = a*X + b*Y + c

cc = plt.contour(X, Y, Z-Z_pl, levels=[0]) 
polys = cc.allsegs[0]  # always one level [0] in our case

下图显示了此代码的结果：

当设置a=1; b=1; c=0.01时，我们得到下图。

如果要使用坐标值，可以通过polys变量（它只是坐标列表列表的列表）访问它们。以下是提取每条折线的x，y坐标的示例：

for poly in polys:
    xs, ys = zip(*poly)  
    # do something with the coordinates...

此方法实际上类似于您尝试为曲线拟合实现的函数，但该实现使用离散数组进行轮廓提取。
从这里开始，它取决于您的应用程序需要什么。如果需要z坐标，可以将x，y值插回平面（或曲面）以获得它们。如果需要平滑曲线（而不是折线），可以将样条曲线（或任何其他smooth fitting method）拟合到生成的多段线。

赞(0）回复(0）举报 2023-04-30

我来回答

scipy 获取插值平面(CloughTocher2DInterpolator)与平面之间的直线交点

1条答案

相关问题

热门标签

最新问答