1

我一直在寻找找到两条线的交点的解决方案。我知道这可以通过找到他们的矢量产品来完成。

我在这里偶然发现了这个例子:

Numpy 和线的交叉点

def get_intersect(a1, a2, b1, b2):
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

我已经浏览了这个例子并在几个场景中使用了它,它似乎工作得很好。但是,有三件事我不太明白:

  1. 为什么向量需要是同质的(我们用一列填充的部分)?
  2. 与非均质解决方案(如果有的话)相比,均质解决方案有何不同?
  3. 为什么我们只检查 Z 轴的平行度而不检查 X 和 Y 轴的平行度?

我觉得我错过了一些非常明显的东西,但我无法理解它是什么。

4

2 回答 2

3

作为参考,来自维基百科的方程式:

[]

a1 = (x1, y1), a2 = (x2, y2), b1 = (x3, y3), b2 = (x4, y4).


计算解释

观察链接答案中的前两个叉积:

l1 = np.cross(h[0], h[1]) = (x1, y1, 1) ^ (x2, y2, 1)
   = (y1 - y2, x2 - x1, x1*y2 - x2*y1)

l2 = np.cross(h[2], h[3]) = (x3, y3, 1) ^ (x4, y4, 1)
   = (y3 - y4, x4 - x3, x3*y4 - x4*y3)

这两行就是计算上述等式中的 6 个不同项所需的全部内容。最后一个交叉产品:

x, y, z = np.cross(l1, l2)
      x = (x2 - x1) * (x3*y4 - x4*y3) - (x4 - x3) * (x1*y2 - x2*y1)
-->   y = (y3 - y4) * (x1*y2 - x2*y1) - (y1 - x2) * (x3*y4 - x4*y3)
      z = (y1 - y2) * (x4 - y3) - (y3 - y4) * (x2 - x1)

这些数字完全等于维基百科等式中的分子和分母。

像这样一个相当复杂的表达式需要数十条 FPU 指令来逐项计算。

均质化向量允许使用这种叉积方法,该方法可以针对少量 SIMD 指令进行优化——效率更高。


几何解释

假设您将同质化向量视为 3D 空间中的点,并将每一对与原点连接起来形成两个三角形:

在此处输入图像描述

所有 4 个点都位于平面Z = 1(灰色)上。

线L(绿色)是两个三角形(蓝色 + 红色)平面之间的交点,并通过原点O和所需的交点P(黄色)。

三角形的法线由其边向量的叉积给出。在这种情况下,边向量仅由 4 个点给出,因为另一个点是原点。在代码中,法线由l1和给出l2

平面的一个定义是位于其中的所有线都必须垂直于其法线。由于线L位于两个三角形的平面内,它必须垂直于l1l2,即它的方向由 给出np.cross(l1, l2)

同质化允许一个巧妙的最后一步,它使用相似的三角形来计算P

在此处输入图像描述

if z == 0:                          # lines are parallel
    return (float('inf'), float('inf'))
return (x/z, y/z)                   # Px = x / z, Py = y / z
于 2019-05-09T14:35:30.587 回答
1

作为对上述答案的轻微修正,使用 3D 叉积来计算 2D 线交点几乎是效率最低的。跨产品根本不进行 SIMD 优化(除非您竭尽全力使用 SOA 数据结构,但即使它是一种极其浪费的方法)。最好的方法是使用良好的旧联立方程。

从定义线的输入点开始:

line1: [(x1, y1), (x2, y2)]
line2: [(x3, y3), (x4, y4)]

计算一些方向向量:

// 1st direction
u1 = x2 - x1
v1 = y2 - y1
D1 = [u1, v1]

// 2nd direction
u2 = x4 - x3
v2 = y4 - y3
D2 = [u2, v2]

现在让我们将线方程重新表述为射线,并为这些射线上的任何点得出一个方程:

// coords of a point 'd1' distance along the 1st ray
Px = x1 + d1*u1
Py = y1 + d1*v1

// coords of a point 'd2' distance along the 2nd ray
Px = x3 + d2*u2
Py = y3 + d2*v2

让我们假设线相交,这意味着两个 P 将是相同的,允许我们声明:

x1 + d1*u1 = x3 + d2*u2
y1 + d1*v1 = y3 + d2*v2

我不会遍历每一步,而是根据 d1 重新排列两个方程,我们得到:

 d1 = x3 + d2*u2 - x1
      ---------------
            u1

 d1 = y3 + d2*v2 - y1
      ---------------
            v1

现在我们有 d1 的两个方程,所以让我们做另一个联立方程来获得 d2 的值:

x3 + d2*u2 - x1     y3 + d2*v2 - y1
---------------  =  ---------------
      u1                  v1

重新排列以隔离 d2:

 d2 = u1(y3 - y1) - v1(x3 - x1)
      -------------------------
           v1*u2 - u1*v2

如果 (v1*u2 - u1*v2) 恰好为零,则该方程没有解(我们称其为行列式,因为它就是这样!)。如果行列式不为零,则只需使用上述等式计算 d2,然后返回到我们之前的等式之一以找到点值:

Px = x3 + d2*u2
Py = y3 + d2*v2 

一些未经测试的 C++ 代码:

bool computeIntersectPoint(
  float x1, float y1,
  float x2, float y2,
  float x3, float y3,
  float x4, float y4,
  float outPoint[2])
{
  // compute direction vectors
  // in some cases, it can be worth 
  // storing the lines as rays as an 
  // optimisation. (avoids 4 subs)
  const float u1 = x2 - x1;
  const float v1 = y2 - x1;
  const float u2 = x4 - x3;
  const float v2 = y4 - x3;

  // check to see if we have a solution
  // 1 mul, 1 fmsub
  float determinant = v1*u2 - u1*v1; 
  if(determinant == 0)
    return false;

  // 2 sub, 1 mul, 1 fmsub
  float temp = u1 * (y3 - y1) - v1 * (x3 - x1);

  // 1 div
  float intersectDistance = temp / determinant;

  // 2 fma
  outPoint[0] = intersectDistance * u2 + x3;
  outPoint[1] = intersectDistance * v2 + y3;

  return true;
}

快速desmos证明:https ://www.desmos.com/calculator/gtlmnmzn6l

此时值得比较两种方法之间的指令数。3d 叉积需要 3 个 mult 和 3 个 fmsub 指令(如果 FMA 不可用,则需要 6 mul + 3 sub)。由于我们有 3 个,我们最多:9 个 mul 和 9 个 fmsub 操作。加上 2 个部门,我们得到:

9 mul
9 fmsub
2 div

我发布的方法需要:

1 div
6 sub
4 fma
2 mul

尽管如果您可以将线存储为射线,则可以保存其中的 4 个。

于 2019-05-14T07:12:36.093 回答