Ray-/segment-intersection 并行性和共线性测试失败 bcz python 中的浮点精度
Ray-/segment-intersection test for parallelism and collinearity fails bcz of float accuracy in python
我正在按照 Gareth Rees 的重要说明,尝试实现一个函数来查找 python 中的 ray-/segment 路口:
and
这是我的函数:
from math import radians, sin, cos
import numpy as np
def find_intersection(point0, theta, point1, point2):
# convert arguments to arrays:
p = np.array(point0, dtype=np.float) # ray origin
q = np.array(point1, dtype=np.float) # segment point 1
q2 = np.array(point2, dtype=np.float) # segment point 2
r = np.array((cos(theta),sin(theta))) # theta as vector (= ray as vector)
s = q2 - q # vector from point1 to point2
rxs = np.cross(r,s)
qpxs = np.cross(q-p,s)
qpxr = np.cross(q-p,r)
t = qpxs/rxs
u = qpxr/rxs
if rxs == 0 and qpxr == 0:
t0 = np.dot(q-p,r)/np.dot(r,r)
t1 = np.dot(t0+s,r)/np.dot(r,r)
return "collinear"
elif rxs == 0 and qpxr != 0:
return "parallel"
elif rxs != 0 and 0 <= t and 0 <= u and u <= 1: # removed t <= 1 since ray is inifinte
intersection = p+t*r
return "intersection is {0}".format(intersection)
else:
return None
该功能在有交叉路口时工作正常。但它不识别并行或共线性,因为条件 rxs == 0 和 qpxr == 0 没有(曾经?)满足。 运行 例如:
p0 = (0.0,0.0)
theta = radians(45.0)
p1 = (1.0,1.0)
p2 = (3.0,3.0)
c = find_intersection(p0,theta,p1,p2)
其中 returns None。在 if 块给出
之前为 rxs 和 qpxr 添加打印语句
rxs = 2.22044604925e-16 qpxr = -1.11022302463e-16
我的结论是,由于浮点问题,该函数未能捕捉到第一个 if 语句的条件。 2.22044604925e-16 和 -1.11022302463e-16 非常小,但不幸的是不完全是 0。我知道浮点数不能用二进制表示。
我的结论是正确的还是我漏掉了什么?有没有什么想法可以避免这个问题?
非常感谢!
是的,你的结论是对的,问题出在"parallel"谓词的数值稳定性上。
您可以将结果与较小的数字进行比较(例如,eps=1.0E-9
)。它的大小可能取决于坐标范围(请注意,叉积给出了双倍的三角形面积,因此将 eps
归一化为 MaxVecLen**2
看起来很合理)。
更复杂但更精确的选项 - 使用稳健的几何谓词,如 these ones。也许 Python/NumPy 计算几何库包含此类操作的一些实现。
有一种简单而安全的方法可以解决这个问题。
写出射线的隐式方程(S(X, Y) = a X + b Y + c = 0
)。当您在函数 S
中插入线段端点的坐标时,您会得到两个值,设 S0
和 S1
。如果符号相反,则射线的支撑线与线段有交点。
在这种情况下,沿线段的交点位置由参数值给出,等于
- S0 / (S1 - S0).
此表达式享有 属性 始终可计算(前提是符号发生变化)且在 [0, 1]
范围内的优点。它允许安全地计算交点。
为了 select 只有那些在所需半线(射线)上的交点,您只需计算射线原点处 S(Xo, Yo)
的符号。
此程序不会检测平行或共线光线,但没关系。无论如何,它都会产生良好的结果。
我正在按照 Gareth Rees 的重要说明,尝试实现一个函数来查找 python 中的 ray-/segment 路口: and
这是我的函数:
from math import radians, sin, cos
import numpy as np
def find_intersection(point0, theta, point1, point2):
# convert arguments to arrays:
p = np.array(point0, dtype=np.float) # ray origin
q = np.array(point1, dtype=np.float) # segment point 1
q2 = np.array(point2, dtype=np.float) # segment point 2
r = np.array((cos(theta),sin(theta))) # theta as vector (= ray as vector)
s = q2 - q # vector from point1 to point2
rxs = np.cross(r,s)
qpxs = np.cross(q-p,s)
qpxr = np.cross(q-p,r)
t = qpxs/rxs
u = qpxr/rxs
if rxs == 0 and qpxr == 0:
t0 = np.dot(q-p,r)/np.dot(r,r)
t1 = np.dot(t0+s,r)/np.dot(r,r)
return "collinear"
elif rxs == 0 and qpxr != 0:
return "parallel"
elif rxs != 0 and 0 <= t and 0 <= u and u <= 1: # removed t <= 1 since ray is inifinte
intersection = p+t*r
return "intersection is {0}".format(intersection)
else:
return None
该功能在有交叉路口时工作正常。但它不识别并行或共线性,因为条件 rxs == 0 和 qpxr == 0 没有(曾经?)满足。 运行 例如:
p0 = (0.0,0.0)
theta = radians(45.0)
p1 = (1.0,1.0)
p2 = (3.0,3.0)
c = find_intersection(p0,theta,p1,p2)
其中 returns None。在 if 块给出
之前为 rxs 和 qpxr 添加打印语句rxs = 2.22044604925e-16 qpxr = -1.11022302463e-16
我的结论是,由于浮点问题,该函数未能捕捉到第一个 if 语句的条件。 2.22044604925e-16 和 -1.11022302463e-16 非常小,但不幸的是不完全是 0。我知道浮点数不能用二进制表示。
我的结论是正确的还是我漏掉了什么?有没有什么想法可以避免这个问题? 非常感谢!
是的,你的结论是对的,问题出在"parallel"谓词的数值稳定性上。
您可以将结果与较小的数字进行比较(例如,eps=1.0E-9
)。它的大小可能取决于坐标范围(请注意,叉积给出了双倍的三角形面积,因此将 eps
归一化为 MaxVecLen**2
看起来很合理)。
更复杂但更精确的选项 - 使用稳健的几何谓词,如 these ones。也许 Python/NumPy 计算几何库包含此类操作的一些实现。
有一种简单而安全的方法可以解决这个问题。
写出射线的隐式方程(S(X, Y) = a X + b Y + c = 0
)。当您在函数 S
中插入线段端点的坐标时,您会得到两个值,设 S0
和 S1
。如果符号相反,则射线的支撑线与线段有交点。
在这种情况下,沿线段的交点位置由参数值给出,等于
- S0 / (S1 - S0).
此表达式享有 属性 始终可计算(前提是符号发生变化)且在 [0, 1]
范围内的优点。它允许安全地计算交点。
为了 select 只有那些在所需半线(射线)上的交点,您只需计算射线原点处 S(Xo, Yo)
的符号。
此程序不会检测平行或共线光线,但没关系。无论如何,它都会产生良好的结果。