计算几何中的浮点误差

大多数几何算法都是为精确的实数算法而设计的。用浮点运算替换它们可能会导致这些实现失败。不幸的是，没有关于几何算法可能出错的原因以及原因的参考文档。帖子的其余部分讨论了使用浮点运算执行几何算法时可能出现的问题。

为了使解释简单，考虑给定点集 S 的凸包算法，即可以包含 S 的所有点的最小多边形； S 的极值点将是凸包多边形的顶点；凸包内的所有点都被简单地丢弃。

通常，由浮点实现计算的凸包可能会遗漏一些极值点，计算非凸多边形，或者算法可能永远运行。

1、增量平面凸包算法

凸包增量算法维护目前所见点的当前凸包 (CH)。最初，通过选择 S 中的三个非共线点来形成凸包。然后逐个考虑其余点。在考虑点r时，首先判断r是否在当前凸包多边形之外。如果不是，则丢弃 r。否则，通过形成从 r 到 CH 的切线来更新凸包。该算法将当前船体维护为其极值点按逆时针顺序排列的循环列表 L = (v0, v1, . . .,vk−1)，其中线段 (vi, vi+1) 是该船体的边缘当前船体。

2、单步失败重现

违反算法正确性的实例。考虑一个序列：p1, p2, p3, . . . 点使得前三个点形成一个逆时针三角形，并且插入一些后面的点导致违反正确性属性（因为浮点）。导致违规的示例总是涉及几乎或真正共线的点；需要明确的是，足够多的非共线点不会造成任何问题。尽管这看起来示例不切实际，但它主要取决于可容忍舍入误差的要求，因为点集可能包含几乎共线的点或真正共线的点，这些点通过转换为浮点表示而变得几乎共线。

浮点方向

平面内的三个点 p = (px, py), q = (qx, qy), r = (rx, ry) 在同一条直线上或形成左转或右转。

三元组 (p, q, r) 的方向定义为：

orientation(p,q,r) = sign((qx − px)(ry − py)−(qy − py)(rx − px))

但是，对于浮点运算，由于可能存在舍入误差，float_orient 的结果可能会在以下三种情况下与正确的方向不同：

四舍五入为零：将 + 或 - 错误分类为 0
扰动零：将 0 错误分类为 + 或 -
符号反转：将 + 错误分类为 - 或反之亦然。

float_orient，换句话说，三重点被分类为左转、右转或共线。

失败一：当前包外的点看不到当前凸包的边

考虑点集：

p1 = (7.3000000000000194, 7.3000000000000167) 
p2 = (24.000000000000068, 24.000000000000071) 
p3 = (24.00000000000005, 24.000000000000053) 
p4 = (0.50000000000001621, 0.5000000000)
p5 = (8, 4) 
p6 = (4, 9) 
p7 = (15, 27) 
p8 = (26, 25) 
p9 = (19, 11)

方向：

float_orient(p1, p2, p3) > 0 
float_orient(p1, p2, p4) > 0 
float_orient(p2, p3, p4) > 0 
float_orient(p3, p1, p4) > 0 (??)

上图显示了计算出的凸包，其中一个明显极值的点被排除在凸包之外； p1 ≈ (17, 17), p2 ≈ (24, 24) ≈ p3

前四个点几乎位于直线上：y = x，float_orient 显示了上述结果。如（??）所示，最后一次评估是错误的。在几何上，这些评估意味着 p4 看不到三角形的边（p1、p2、p3）。点 p5, . . . , p9 然后被正确识别为极值点并添加到船体中。然而，该算法永远不会从考虑 p4 时所犯的错误中恢复过来，并且计算结果与正确的船体大不相同。

失败 2：当前船体内部的点看到当前船体的边缘

取一个逆时针三角形（初始）并选择三角形内的第四个点，但非常靠近其中一条边。有符号反转的机会。例如，考虑点集：

p1 = (27.643564356435643, −21.881188118811881 ) 
p2 = (83.366336633663366, 15.544554455445542 ) 
p3 = ( 4.0, 4.0 ) 
p4 = (73.415841584158414, 8.8613861386138595)

方向：

float orient(p1, p2, p3) > 0 
float orient(p1, p2, p4) < 0 (??) 
float orient(p2, p3, p4) > 0 
float orient(p3, p1, p4) > 0

凸包被正确初始化为 (p1, p2, p3)。点 p4 在当前凸包内，但算法错误地认为 p4 可以看到边 (p1, p2)，因此将包更改为 (p1, p4, p2, p3)，一个稍微非凸的多边形。

失败3：当前包外的点看到凸包的所有边

考虑点集：

p1 = ( 200.0, 49.200000000000003) 
p2 = ( 100.0, 49.600000000000001) 
p3 = (−233.33333333333334, 50.93333333333333 ) 
p4 = ( 166.66666666666669, 49.333333333333336)

方向：

float orient(p1, p2, p3) > 0 
float orient(p1, p2, p4) < 0 
float orient(p2, p3, p4) < 0 
float orient(p3, p1, p4) < 0 (??)

前三个点组成一个逆时针方向的三角形，根据float_orient，算法认为p4可以看到三角形的所有边。

失败 4：当前船体外的点看到一组不连续的边

考虑点集：

p1 = (0.50000000000001243, 0.50000000000000189) 
p2 = (0.50000000000001243, 0.50000000000000333) 
p3 = (24.00000000000005, 24.000000000000053) 
p4 = (24.000000000000068, 24.000000000000071) 
p5 = (17.300000000000001, 17.300000000000001)

方向：

float_orient(p1, p4, p5) < 0 (??) 
float_orient(p4, p3, p5) > 0 
float_orient(p3, p2, p5) < 0 
float_orient(p2, p1, p5) > 0

插入前三个点，然后插入第四个点 p4 得到凸四边形 (p1, p4, p3, p2)，这是正确的。最后一个点 p5 只看到边缘 (p3, p2) 而看不到其他三个。然而，float_orient 使 p5 也能看到边缘 (p1, p4)。可见边和不可见边的子序列不连续。由于错误分类的边 (p1, p4) 先出现，该算法在该边插入 p5，不删除任何其他顶点，并返回具有自相交的多边形。

3、结束语

如上所示，算法可能会因舍入误差而失败，以天真地使用浮点运算实现的凸包算法为例。

与往常一样，可以在此处找到对 C++ 代码的引用，以及由于舍入错误而失败的凸算法变体以及如何防止它。

作为扩展，考虑使用精确浮点数 (MP_Float) 来尝试 CGAL 库。

原文链接：Floating Point Round Off Errors in Geometric Algorithms

BimAnt翻译整理，转载请标明出处