测试两条线段相交时的算术精度问题答案

【问题标题】：Arithmetic accuracy issues while testing for intersection of two line segments测试两条线段相交时的算术精度问题
【发布时间】：2016-07-05 06:38:59
【问题描述】：

我写了一个代码测试平面中两条线段的交点。我不会用所有的细节来打扰你。

代码采用两条线段，每条线段由两个端点描述，然后通过在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

现在的问题：我发现由于舍入误差，当交点与段边缘重合时，测试会给出错误的结果。

例如，线 1 介于 (0,0) 和 (3,5) 之间，线 2 介于 (3,5) 和 (7,1) 之间，生成的交点 x 为2.9999999999999996，给出错误答案。应该是 3。

您能提出一个解决方案吗？

【问题讨论】：

使用十进制类：from decimal import Decimal。 stackoverflow.com/questions/2986150/python-floating-number 。您的 isininterval 函数也可以简单地返回：return x1 <= x0 <= x2 or x2 <= x0 <= x1
编写涉及浮点数的布尔条件的标准方法是在条件中包含一些容错。例如，将x1 <= x0 <= x1 替换为x1 - eps <= x0 <= x2+ eps，其中eps 类似于0.000001
首先，你有充分的理由担心 2.9999999999999996 不是 3 吗？
@YvesDaoust，好吧，我唯一担心的是它不符合条件......
是否有充分的理由让一个段的顶点正好落在另一段上？这些段是孤立的还是形成链或图？

标签： python floating-point precision

【解决方案1】：

这是浮点运算的问题/特征。有一些方法可以通过以某种方式对指令进行排序来最小化错误，但最终你会得到近似的答案，因为你用有限的位数表示可能无限的数字。

您需要以能够容忍这些错误的方式定义您构建的任何函数。查看您的示例，“正确”值与您得到的值之间的差异是 1e-16 的顺序 - 非常非常低。

对于不等式，尤其是等式，放宽精确/逐位匹配的约束是值得的。例如，如果你想测试x == 3，你可以写成abs(x - 3) < EPSILON，其中EPSILON = 1e-6 或EPSILON = 1e-9。基本上，你想要拥有的东西和你拥有的东西之间的差异小于一个非常小的值。同上，对于不等式，您可以测试 3 - EPSILON <= x 或 x <= 3 + EPSILON。

【讨论】：

是否有机器定义的 epsilon？是否适用于所有机器？
sys 模块中有 sys.float_info.epsilon，但这可能有点偏低，因为从技术上讲，它被定义为 1 和大于 1 的最小数之间的差表示为浮点数。所以 [1, 1+eps] 或多或少之间没有任何浮动。但是，这是特定于应用程序的，您绝对应该随意尝试。
实际上，您可以通过编写一些真正仔细的代码来缓解这个“问题/功能”。我相信这项技术归功于 Jonathan Shewchuk。见cs.cmu.edu/~quake/robust.html。
您可能想要使用“相对差异”，例如abs(x - target) / max(abs(x), abs(target)) < epsilon，这样如果数字非常小（1E-32），您就不会接受实际上相差很大的值，并且同样，对于非常大 (1E+32) 的值也不会遇到问题。