测试两条线段相交时的算术准确性问题
Arithmetic accuracy issues while testing for intersection of two line segments
我写了一个代码测试平面中两条线段的交点。我不会用所有的细节来打扰你。
代码采用两条线段,每条线段由两个端点描述,然后通过在 y = a*x + b
中拟合 a
和 b
将每个线段拟合成一条直线。然后通过x = (b2 - b1) / (a2 - a1)
找到两条直线的交点。最后,它测试交点 x
是否包含在两条线段内。
相关部分如下所示:
# line parameterization by a = Delta y / Delta x, b = y - a*x
a1 = (line1.edge2.y - line1.edge1.y) / (line1.edge2.x - line1.edge1.x)
b1 = line1.edge1.y - a1 * line1.edge1.x
a2 = (line2.edge2.y - line2.edge1.y) / (line2.edge2.x - line2.edge1.x)
b2 = line2.edge1.y - a2 * line2.edge1.x
# The intersection's x
x = - (b2 - b1) / (a2 - a1)
# If the intersection x is within the interval of each segment
# then there is an intersection
if (isininterval(x, line1.edge1.x, line1.edge2.x) and
isininterval(x, line2.edge1.x, line2.edge2.x)):
return True
else:
return False
为简洁起见,我放弃了很多处理特定情况的测试,例如当边彼此平行时(a1==a2
),当它们在同一条线上时,当边的长度为 0 时,当边缘沿垂直轴(然后 a
变为无限大)等
函数isininterval
就是
def isininterval(x0, x1, x2):
"""Tests if x0 is in the interval x1 to x2"""
if x1 <= x0 <= x2 or x2 <= x0 <= x1:
return True
else:
return False
现在的问题是:我发现由于舍入误差,当交点与线段边缘重合时,测试会给出错误的结果。
例如,line1 在 (0,0) 和 (3,5) 之间,line2 在 (3,5) 和 (7,1) 之间,得到的交点 x 是 2.9999999999999996
,这给出了错误的答案。应该是 3.
你能提出一个解决方案吗?
这是一个problem/feature的浮点运算。有一些方法可以通过以某种方式对指令进行排序来最大限度地减少错误,但最终您会得到近似的答案,因为您可能用有限的位数来表示无限的数字。
您需要以能够容忍这些错误的方式定义您构建的任何函数。查看您的示例,"proper" 值与您得到的值之间的差异约为 1e-16
- 极低。
由于不平等,尤其是平等,放宽 exact/bitewise 匹配的限制是值得的。例如,如果您想测试 x == 3
,您可以将其写为 abs(x - 3) < EPSILON
,其中 EPSILON = 1e-6
或 EPSILON = 1e-9
。基本上,你想要拥有的和你拥有的之间的差异小于一个非常小的值。同上,对于不等式,您可以测试 3 - EPSILON <= x
或 x <= 3 + EPSILON
.
我写了一个代码测试平面中两条线段的交点。我不会用所有的细节来打扰你。
代码采用两条线段,每条线段由两个端点描述,然后通过在 y = a*x + b
中拟合 a
和 b
将每个线段拟合成一条直线。然后通过x = (b2 - b1) / (a2 - a1)
找到两条直线的交点。最后,它测试交点 x
是否包含在两条线段内。
相关部分如下所示:
# line parameterization by a = Delta y / Delta x, b = y - a*x
a1 = (line1.edge2.y - line1.edge1.y) / (line1.edge2.x - line1.edge1.x)
b1 = line1.edge1.y - a1 * line1.edge1.x
a2 = (line2.edge2.y - line2.edge1.y) / (line2.edge2.x - line2.edge1.x)
b2 = line2.edge1.y - a2 * line2.edge1.x
# The intersection's x
x = - (b2 - b1) / (a2 - a1)
# If the intersection x is within the interval of each segment
# then there is an intersection
if (isininterval(x, line1.edge1.x, line1.edge2.x) and
isininterval(x, line2.edge1.x, line2.edge2.x)):
return True
else:
return False
为简洁起见,我放弃了很多处理特定情况的测试,例如当边彼此平行时(a1==a2
),当它们在同一条线上时,当边的长度为 0 时,当边缘沿垂直轴(然后 a
变为无限大)等
函数isininterval
就是
def isininterval(x0, x1, x2):
"""Tests if x0 is in the interval x1 to x2"""
if x1 <= x0 <= x2 or x2 <= x0 <= x1:
return True
else:
return False
现在的问题是:我发现由于舍入误差,当交点与线段边缘重合时,测试会给出错误的结果。
例如,line1 在 (0,0) 和 (3,5) 之间,line2 在 (3,5) 和 (7,1) 之间,得到的交点 x 是 2.9999999999999996
,这给出了错误的答案。应该是 3.
你能提出一个解决方案吗?
这是一个problem/feature的浮点运算。有一些方法可以通过以某种方式对指令进行排序来最大限度地减少错误,但最终您会得到近似的答案,因为您可能用有限的位数来表示无限的数字。
您需要以能够容忍这些错误的方式定义您构建的任何函数。查看您的示例,"proper" 值与您得到的值之间的差异约为 1e-16
- 极低。
由于不平等,尤其是平等,放宽 exact/bitewise 匹配的限制是值得的。例如,如果您想测试 x == 3
,您可以将其写为 abs(x - 3) < EPSILON
,其中 EPSILON = 1e-6
或 EPSILON = 1e-9
。基本上,你想要拥有的和你拥有的之间的差异小于一个非常小的值。同上,对于不等式,您可以测试 3 - EPSILON <= x
或 x <= 3 + EPSILON
.