为什么我们需要均质化这些载体?
Why do we need to homogenise these vectors?
我一直在寻找找到两条线交点的解决方案。我知道这可以通过找到他们的向量积来完成。
我在这里偶然发现了这个例子:
Numpy and line intersections
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 填充列的部分)?
- 与非齐次解决方案(如果有的话)相比,齐次解决方案有何不同?
- 为什么我们只检查 Z 轴平行度的结果,而不检查 X 轴和 Y 轴的平行度?
我觉得我遗漏了一些非常明显的东西,但我无法理解它是什么。
供参考,来自维基百科的方程式:
让a1 = (x1, y1), a2 = (x2, y2), b1 = (x3, y3), b2 = (x4, y4)
.
计算解释
观察链接答案中的前两个 cross-products:
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 项不同项的全部内容。最后一个 cross-product:
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 指令来计算 term-by-term。
对向量进行均质化允许使用这种 cross-product 方法,该方法可以针对少量 SIMD 指令进行优化——效率更高。
几何解释
假设您将均化向量视为 3D 中的点 space,并将每对与原点连接起来形成两个三角形:
所有 4 个点都在平面上 Z = 1(灰色)。
直线L(绿色)是两个三角形(蓝色+红色)平面的交点,并通过原点O和想要的交点P(黄色).
三角形的 法线 由其侧向量的 cross-product 给出。在这种情况下,侧向量仅由 4 个点给出,因为另一个点是原点。在代码中,法线由 l1
和 l2
.
给出
平面的一个定义是平面内的所有直线都必须垂直于平面的法线。由于线 L 位于两个三角形的平面内,因此它必须垂直于 l1
和 l2
,即它的方向由 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
作为对上述答案的轻微修正,使用 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
我不会逐一进行,但是 re-arrange 两个关于 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 个潜艇。
我一直在寻找找到两条线交点的解决方案。我知道这可以通过找到他们的向量积来完成。
我在这里偶然发现了这个例子:
Numpy and line intersections
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 填充列的部分)?
- 与非齐次解决方案(如果有的话)相比,齐次解决方案有何不同?
- 为什么我们只检查 Z 轴平行度的结果,而不检查 X 轴和 Y 轴的平行度?
我觉得我遗漏了一些非常明显的东西,但我无法理解它是什么。
供参考,来自维基百科的方程式:
让a1 = (x1, y1), a2 = (x2, y2), b1 = (x3, y3), b2 = (x4, y4)
.
计算解释
观察链接答案中的前两个 cross-products:
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 项不同项的全部内容。最后一个 cross-product:
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 指令来计算 term-by-term。
对向量进行均质化允许使用这种 cross-product 方法,该方法可以针对少量 SIMD 指令进行优化——效率更高。
几何解释
假设您将均化向量视为 3D 中的点 space,并将每对与原点连接起来形成两个三角形:
所有 4 个点都在平面上 Z = 1(灰色)。
直线L(绿色)是两个三角形(蓝色+红色)平面的交点,并通过原点O和想要的交点P(黄色).
三角形的 法线 由其侧向量的 cross-product 给出。在这种情况下,侧向量仅由 4 个点给出,因为另一个点是原点。在代码中,法线由 l1
和 l2
.
平面的一个定义是平面内的所有直线都必须垂直于平面的法线。由于线 L 位于两个三角形的平面内,因此它必须垂直于 l1
和 l2
,即它的方向由 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
作为对上述答案的轻微修正,使用 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
我不会逐一进行,但是 re-arrange 两个关于 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 个潜艇。