如何通过一个曲面通过四个点并找到相邻点在曲面上的投影

Question

我有一个复杂的几何问题。 我有两个系列的点，由它们的 x、y 和 z 表示为 numpy 数组（实际上我在每个系列中有数千个，但为了更好地理解，在这里简单化）。 第一个数据集 (arr_1) 表示切割表面。该表面通常是垂直的并取代水平表面。第二组是被一个表面（arr_1 所代表的表面）置换的垂直表面的点。目前我在这两组中有一些冗余数据。首先，在arr_1 中，我更喜欢只有四个角点并自己使用它们制作一个曲面。其次，在arr_2 中我也不需要太接近arr_1 的数据。最后，我想终止由arr_1 的四个角构成的曲面两侧的点。

import numpy as np
arr_1= np.array ([[31.95548952,  7.5       , 12.5       ],
       [31.95548952, 22.5       , 12.5       ],
       [37.5       ,  7.5       , 24.20636043],
       [37.5       , 22.5       , 24.20636043],
       [43.78154278,  7.5       , 37.5       ],
       [43.78154278, 22.5       , 37.5       ],
       [55.32209575,  7.5       , 62.5       ],
       [55.32209575, 22.5       , 62.5       ],
       [62.5       ,  7.5       , 78.50696445],
       [62.5       , 22.5       , 78.50696445],
       [66.52446985,  7.5       , 87.5       ],
       [66.52446985, 22.5       , 87.5       ]])
arr_2= np.array([[87.5       ,  7.5       , 49.99997914],
       [87.5       , 22.5       , 50.00001192],
       [62.5       ,  7.5       , 50.00004172],
       [62.5       , 22.5       , 50.00007749],
       [46.8747884 ,  7.5       , 62.5       ],
       [46.87483609, 22.5       , 62.5       ],
       [37.5       ,  7.5       , 69.99973059],
       [37.5       , 22.5       , 69.99977231],
       [12.5       ,  7.5       , 70.00012398],
       [12.5       , 22.5       , 70.00015974]])

然后，我使用以下代码找到这两个数组的每组之间的距离：

from scipy.spatial import distance
distances=distance.cdist(arr_1, arr_2)

当我运行它时，distances 是一个包含 12 行和 10 列的数组。 现在，我想删除arr_2 的点，它们比阈值（比如说10）更接近arr_1 的任何点。在第一张上传的照片中，我在黑色矩形中显示了这两个点。提出了一个很好的解决方案here，但是并没有解决我的问题，因为我想比较arr_1的每一行与arr_2的每一行的距离。我很感激任何解决方案。 事实上，arr_1 包含曲面的点，但我不需要所有这些点来制作我的曲面。我只挑选这组的角落来制作我的表面。我使用了一个非常耗时的 for 循环，所以我很欣赏任何更快的方法来找到我的观点的角落：

corner1 = np.array(arr_1.min (axis=0)) # this gives the row containing MIN values of x,y and z
corner2 = np.array([])
corner4 = np.array([])
corner3 = np.array(arr_1.max (axis=0)) # this gives the row containing MAX values of x,y and z
# the next block will find the corner in which x and y are minimum and z is maximum
for i in arr_1[:,0]:
    if i == max (arr_1[:,0]):
        for j in arr_1[:,1]: 
            if j == min (arr_1[:,1]):
                for h in arr_1[:,2]:
                    if h == max (arr_1[:,2]):
                        corner2 = np.append(corner2, np.array([i,j,h]))
                        corner2=corner2[0:3]
# the next block will find the corner in which x and z are minimum and y is maximum
for m in arr_1[:,0]:
    if m == min (arr_1[:,0]):
        for n in arr_1[:,1]: 
            if n == max (arr_1[:,1]):
                for o in arr_1[:,2]:
                    if o == min (arr_1[:,2]):
                        corner4 = np.append(corner4, np.array([m,n,o]))
                        corner4=corner4[0:3]

最后，在提取四个角后，我想用它们做一个曲面，并在其上找到 arr_2 相邻点的垂直（绿色箭头）或水平（红色箭头）投影。我不知道如何找到表面和投影。 感谢您阅读并关注我的详细问题！如果有人对它的任何部分提出任何解决方案，我将不胜感激。

Answer 1

让我们把这个问题分成几个部分。

1. 您有一堆描述噪声平面的数据点，arr_1。
2. 您想知道它与另一个平面相交的位置，由 arr_2 描述。
3. 您想在arr_2 中设置关于该交叉口的阈值。

我在这里展示的方法是假设数据是某个真值的测量值，并且您希望根据对该值的最佳猜测执行这些操作，而不是原始数据。为此：

第 1 部分：适合平面的最小二乘

有几种不同的方式来描述平面，例如法线向量和点。最简单的最小二乘拟合可能是

a * x + b * y + c * z = 1

假设您的数据表现得相当好，使用方程进行简单拟合应该没有问题

arr_1 @ [[a], [b], [c]] = 1   # almost python pseudo code

由于没有超过四个点的单一解决方案，您可以运行np.linalg.lstsq 以根据 MSE 优化值：

plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()

如果arr_2 也是一架飞机，你对它做同样的事情来得到plane_2。

第 2 部分：平面相交

许多平面相交的解决方案都依赖于平面的法向量。我将假设两个平面都依赖于 Y 坐标（从图中看起来很安全）。在这种情况下，您可以在 Math Stack Exchange 上关注 this answer。设置y = t，可以从系统中解线

a1 * x + c1 * z = 1 - b1 * t
a2 * x + c2 * z = 1 - b2 * t

这里，向量[a1, b1, c1] 是plane_1。搞清楚细节后，你得到

m = np.cross(plane_1, plane_2)
b = np.array([plane_1[2] - plane_2[2], 0, plane_2[0] - plane_1[0]]) / m[1]
m /= m[1]
line = (m * t + b)

这是t 的任何值的参数化。

第 3 部分：点到线的距离

要对上面根据m 和b 计算的线的arr_2 值进行阈值化，您需要一个点和线之间距离的公式。这可以通过 here 和 here 帖子中的方法来完成。

例如单点p可以这样处理：

t = (p - b).dot(m) / m.dot(m)
dist = np.sqrt(np.sum((p - (m * t + b))**2))

如果您只对阈值处理感兴趣，您可以将dist**2 与阈值的平方进行比较，并在平方根上节省一些循环，因为这两个函数都是单调的。

TL;DR

输入：arr_1, arr_2, distance_threshold。

# Find planes in data
plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()
plane_2 = np.linalg.lstsq(arr_2, np.ones((arr_2.shape[0], 1), dtype=arr_2.dtype))[0].ravel()

# Find intersection line, assume plane is not y=const
m = np.cross(plane_1, plane_2)
b = np.array([plane_1[2] - plane_2[2], 0, plane_2[0] - plane_1[0]]) / m[1]
m /= m[1]

# Find mask for arr_2
t = ((arr_2 - b).dot(m) / m.dot(m))[:, None]
dist2 = np.sum((arr_2 - (m * t + b))**2, axis=1)
mask = dist2 >= distance_threshold**2

# Apply mask
subset = arr_2[mask, :]

附录 1：RMSE

如前所述，这种方法的真正好处（除了它将您的算法限制在 ~O(n) 的事实之外）是它对您的数据进行去噪。您可以使用最小二乘拟合的结果来计算有关拟合的平面数据的 RMSE，以了解您的测量值的真实程度。 lstsq 的第二个返回值是 RMSE 指标。

附录2：点到平面的距离

如果arr_2 中的数据确实不是平面的，您可以对它进行一些不同的子集化。您可以直接使用单个点与平面之间的距离公式，而不是找到一对平面的交点，如给定here：

np.abs(p * plane_1 - 1) / np.sqrt(plane1.dot(plane_1))

代码就变成了

# Find planes in data
plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()
plane_2 = np.linalg.lstsq(arr_2, np.ones((arr_2.shape[0], 1), dtype=arr_2.dtype))[0].ravel()

# Find mask for arr_2
dist2 = (arr_2 * plane_1 - 1)**2 / plane_1.dot(plane_1)
mask = dist2 >= distance_threshold**2

# Apply mask
subset = arr_2[mask, :]

Answer 2

这是我整理的一些代码 - 它涵盖了您问题的第一部分 - 即从您的点集合中查找和定义“表面”。

任何平面（平面）都可以用公式ax + by + cz + d = 0定义，其中a,b,c是系数，d是一些常数。我们可以将其写入如下函数：

def linear_plane(data, a, b, c):
    x = data[0]
    y = data[1]
    z = data[2]
    return (a * x) + (b * y) + (c*z)

然后，在一些导入之后，我们可以使用 scipy 的 curve_fit 来搜索并找到最适合您的每个数组的函数的参数。 curve_fit 接受 linear_plane 函数的输入值并尝试调整 a、b、c 以使该函数的结果与包含大量 -1 的列表相匹配。这来自于重新排列 d = 1

的情况的方程

from scipy.optimize import curve_fit

v1,err1 = curve_fit(linear_plane,arr_1.T,np.repeat(-1,len(arr_1)))
v2,err2 = curve_fit(linear_plane,arr_2.T,np.repeat(-1,len(arr_2)))

v1、v2 中的每一个都是最接近定义您正在谈论的两个平面的系数 a、b、c。 err1 和 err2 显示剩余的“错误”，因此请注意它们是否太大。

现在找到了系数，您可以使用它们来可视化两个平面：

%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(xs=arr_1[:,0], ys=arr_1[:,1], zs=arr_1[:,2], c="blue")
ax.scatter3D(xs=arr_2[:,0], ys=arr_2[:,1], zs=arr_2[:,2], c="red")
xx1, yy1 = np.meshgrid([min(arr_1[:,0]),max(arr_1[:,0])],[min(arr_1[:,1]),max(arr_1[:,1])])
xx2, yy2 = np.meshgrid([min(arr_2[:,0]),max(arr_2[:,0])], [min(arr_2[:,1]),max(arr_2[:,1])])

# Use the coefficients in v1[0], v1[1], v1[2] to calculate z-values for all xx/yy values
z1 = (-v1[0] * xx1 - v1[1] * yy1 - 1) * 1. /v1[2]


# Use the coefficients in v2[0], v2[1], v2[2] to calculate z-values for all xx/yy values
z2 = (-v2[0] * xx2 - v2[1] * yy2 - 1) * 1. /v2[2]

# Using xx,yy and z values, plot the surfaces fit into the graph.
ax.plot_surface(xx1, yy1, z1, alpha=0.2, color=[0,0,1])
ax.plot_surface(xx2, yy2, z2, alpha=0.2, color=[1,0,0])

在这里，我们计算所有 x 和 y 的 z 值集合，这些值落在每个平面的最大和最小范围内，然后将它们绘制为显示相交的曲面：

我注意到Mad Physicist 刚刚发布了一个更完整的答案，所以暂时将其留在这里。我想做的一件事是扩展表面拟合，这样平面就不必是线性的——就像类似于一些 3d 多项式，应该只是一个例子替换 curve_fit 调用中使用的函数。

这里缺少的另一部分是如何计算描述这两个平面交叉位置的直线方程，我相信可以在另一个（更好的）答案中找到。