<p>让我们把这个问题分成几个部分</p>
<ol>
<li>你有一堆数据点描述一个嘈杂的平面,<code>arr_1</code></李>
<li>您想知道它与另一个平面相交的位置,如<code>arr_2</code>所述</李>
<li>您想在<code>arr_2</code>中建立关于该交叉点的一些阈值</李>
</ol>
<p>我在这里展示的方法是假设数据是某个真值的度量,您希望对该值的最佳猜测执行这些操作,而不是原始数据。为此目的:</p>
<p><strong>第1部分:平面最小二乘拟合</p>
<p>有两种不同的方法来描述平面,例如法向量和点。最小二乘拟合最简单的方法可能是</p>
<pre><code>a * x + b * y + c * z = 1
</code></pre>
<p>假设您的数据表现良好,那么使用该方程进行简单拟合应该没有问题</p>
<pre><code>arr_1 @ [[a], [b], [c]] = 1 # almost python pseudo code
</code></pre>
<p>由于不可能有超过四个点的单一解决方案,因此可以运行<a href="https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html" rel="nofollow noreferrer">^{<cd4>}</a>以获得根据MSE优化的值:</p>
<pre><code>plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()
</code></pre>
<p>如果<code>arr_2</code>也是一个平面,则对它执行相同的操作以获得<code>plane_2</code></p>
<p><strong>第2部分:平面相交</p>
<p>许多平面相交的解依赖于平面的法向量。我将假设这两个平面都依赖于Y坐标(从图上看,这似乎是安全的)。在这种情况下,您可以在数学堆栈交换上遵循<a href="https://math.stackexchange.com/a/475983/295281">this answer</a>。设置<code>y = t</code>,您可以从系统中求解该行</p>
<pre><code>a1 * x + c1 * z = 1 - b1 * t
a2 * x + c2 * z = 1 - b2 * t
</code></pre>
<p>这里,向量<code>[a1, b1, c1]</code>是<code>plane_1</code>。在解决了细节之后,你会得到</p>
<pre><code>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)
</code></pre>
<p>这是<code>t</code>的任何值的参数化</p>
<p><strong>第3部分:点到线的距离</strong></p>
<p>要根据<code>m</code>和<code>b</code>对上面计算的线的<code>arr_2</code>值设置阈值,需要一个点与线之间距离的公式。这可以通过例如posts<a href="https://math.stackexchange.com/q/1815397/295281">here</a>和<a href="https://math.stackexchange.com/q/371649/295281">here</a>中的方法来实现</p>
<p>例如,单点<code>p</code>可以按如下方式处理:</p>
<pre><code>t = (p - b).dot(m) / m.dot(m)
dist = np.sqrt(np.sum((p - (m * t + b))**2))
</code></pre>
<p>如果您只对阈值处理感兴趣,可以将<code>dist**2</code>与阈值的平方进行比较,并在平方根上保存一些循环,因为这两个函数都是单调的</p>
<p><strong>TL;DR</strong></p>
<p>输入:<code>arr_1, arr_2, distance_threshold</code></p>
<pre><code># 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, :]
</code></pre>
<p><strong>附录1:RMSE</strong></p>
<p>如前所述,这种方法的真正优点(除了它将算法限制在~O(n)之外)是它消除了数据中的噪声。您可以使用最小二乘拟合的结果来计算有关拟合的平面数据的RMSE,以了解测量值的实际平面度。<code>lstsq</code>的第二个返回值是RMSE度量</p>
<p><strong>附录2:点到平面的距离</p>
<p>如果<code>arr_2</code>中的数据确实不是平面的,您可以稍微改变它的子集。您可以直接使用单个点与平面之间的距离公式,而不是查找一对平面的交点,如<a href="https://math.stackexchange.com/a/493775/295281">here</a>所示:</p>
<pre><code>np.abs(p * plane_1 - 1) / np.sqrt(plane1.dot(plane_1))
</code></pre>
<p>然后代码变为</p>
<pre><code># 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, :]
</code></pre>