python - 如何让一个曲面通过四个点，求相邻点在曲面上的投影

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 中比 arr_1 的任何点更接近阈值(比方说 10)的点。我在第一张上传的照片中用黑色矩形显示了这两个点。提出一个好的解决方案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 坐标(从图中看这似乎是安全的)。那样的话可以关注this answer在 Math Stack Exchange 上。设置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 与阈值的平方进行比较，并在平方根上保存一些循环，因为这两个函数都是单调的。

长话短说

输入: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 指标。

附录二:点面距离

如果 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, :]