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

Question

我有一个复杂的几何问题。 我有两个系列的点，用它们的x、y和z表示为numpy数组（实际上，我在每个系列中有数千个点，但为了更好地理解，这里将其简化）。 第一个数据集（arr_1）表示切割面。该曲面通常是垂直的，并置换水平曲面。第二组是由曲面（表示arr_1的曲面）置换的垂直曲面的点。目前，我在这两组数据中有一些冗余数据。首先，在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

由于不可能有超过四个点的单一解决方案，因此可以运行^{}以获得根据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。设置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值设置阈值，需要一个点与线之间距离的公式。这可以通过例如postshere和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

下面是我编写的一些代码，它涵盖了问题的第一部分，即从点集合中查找和定义“曲面”

任何平面（平面）都可以用公式定义，其中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的a列表匹配。这来自于重新排列方程，其中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调用中使用的函数的例子

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