Segment Segment 交集代码的精度问题

Question

我正在调试我编写的一种算法，该算法通过识别所有交点并绕着外围走动，将复杂的自相交多边形转换为简单的多边形。

我写了一系列随机数据生成压力测试，在这一次中，我遇到了一个有趣的情况，导致我的算法在正确运行数千次后失败。

double fRand(double fMin, double fMax)
{
    double f = (double)rand() / RAND_MAX; // On my machine RAND_MAX = 2147483647
    return fMin + f * (fMax - fMin);
}
// ... testing code below:  
srand(1);
for(int j=3;j<5000;j+=5) {
    std::vector<E_Point> geometry;
    for(int i=0;i<j;i++) {
        double radius = fRand(0.6,1.0);
        double angle = fRand(0,2*3.1415926535);
        E_Point pt(cos(angle),sin(angle));
        pt *= radius;
        geometry.push_back(pt); // sending in a pile of shit
    }
    // run algorithm on this geometry

这是对它的相关性以及我如何达到现在的位置的概述。还有很多细节我要省略。

我能够做的是将问题缩小到我正在使用的段-段交叉代码：

bool intersect(const E_Point& a0, const E_Point& a1,
               const E_Point& b0, const E_Point& b1,
               E_Point& intersectionPoint) {

    if (a0 == b0 || a0 == b1 || a1 == b0 || a1 == b1) return false;
    double x1 = a0.x; double y1 = a0.y;
    double x2 = a1.x; double y2 = a1.y;
    double x3 = b0.x; double y3 = b0.y;
    double x4 = b1.x; double y4 = b1.y;

    //AABB early exit
    if (b2Max(x1,x2) < b2Min(x3,x4) || b2Max(x3,x4) < b2Min(x1,x2) ) return false;
    if (b2Max(y1,y2) < b2Min(y3,y4) || b2Max(y3,y4) < b2Min(y1,y2) ) return false;

    float ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3));
    float ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3));
    float denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1);
    // check against epsilon (lowest normalized double value)
    if (fabs(denom) < DBL_EPSILON) {
        //Lines are too close to parallel to call
        return false;
    }
    ua /= denom;
    ub /= denom;

    if ((0 < ua) && (ua < 1) && (0 < ub) && (ub < 1)) {
        intersectionPoint.x = (x1 + ua * (x2 - x1));
        intersectionPoint.y = (y1 + ua * (y2 - y1));
        return true;
    }
    return false;
}

发生的情况是我有两个交点，此函数为其返回完全相同的交点值。这是相关几何的放大视图：

垂直线由点定义 (0.3871953044519425, -0.91857980824611341), (0.36139704793723609, 0.91605957361605106)

由(0.8208980020500205, 0.52853407296583088), (0.36178501611208552, 0.88880385168617226) 点组成的几乎水平的绿色线

和(0.36178501611208552, 0.88880385168617226), (-0.43211245441046209, 0.68034202227710472)的白线

如您所见，最后两行确实有一个共同点。

我的函数给了我(0.36178033094571277, 0.88880245640159794)的解决方案 对于这两个交叉口（您在图片中看到一个红点）。

这是一个大问题的原因是我的周界算法依赖于对每条边上的交点进行排序。由于这两个交点都被计算为具有相同的精确值，因此排序将它们置于错误的方向。周界的路径来自上方，沿着左边的白线，而不是左边的绿线，这意味着我不再沿着多边形的外周长。

要解决这个问题，我可能可以做很多事情，但我不想搜索所有交集列表来检查其他点以查看位置是否相等。更好的解决方案是尝试提高我的交集函数的准确性。

所以我要问的是，为什么解点如此不准确？是因为其中一条线几乎是垂直的吗？我应该先进行某种转换吗？在这两种情况下，几乎垂直的线都以a0 和a1 传递。

更新：嘿，看看这个：

TEST(intersection_precision_test) {
    E_Point problem[] = {

    {0.3871953044519425, -0.91857980824611341}, // 1559
    {0.36139704793723609, 0.91605957361605106}, // 1560
    {-0.8208980020500205, 0.52853407296583088}, // 1798
    {0.36178501611208552, 0.88880385168617226}, // 1799
    {-0.43211245441046209, 0.6803420222771047} // 1800
    };
    std::cout.precision(16);
    E_Point result;
    intersect(problem[0],problem[1],problem[2],problem[3],result);
    std::cout << "1: " << result << std::endl;
    intersect(problem[0],problem[1],problem[3],problem[2],result);
    std::cout << "2: " << result << std::endl;
    intersect(problem[1],problem[0],problem[2],problem[3],result);
    std::cout << "3: " << result << std::endl;
    intersect(problem[1],problem[0],problem[3],problem[2],result);
    std::cout << "4: " << result << std::endl;

    intersect(problem[2],problem[3],problem[0],problem[1],result);
    std::cout << "rev: " << result << std::endl;
    intersect(problem[3],problem[2],problem[0],problem[1],result);
    std::cout << "revf1: " << result << std::endl;
    intersect(problem[2],problem[3],problem[1],problem[0],result);
    std::cout << "revf2: " << result << std::endl;
    intersect(problem[3],problem[2],problem[1],problem[0],result);
    std::cout << "revfboth: " << result << std::endl;
}

输出：

Starting Test intersection_precision_test, at Polygon.cpp:1830
1: <0.3617803309457128,0.8888024564015979>
2: <0.3617803309457128,0.8888024564015979>
3: <0.3617803314022162,0.8888024239374175>
4: <0.3617803314022162,0.8888024239374175>
rev: <0.3617803635476076,0.8888024344185281>
revf1: <0.3617803313928456,0.8888024246235207>
revf2: <0.3617803635476076,0.8888024344185281>
revfboth: <0.3617803313928456,0.8888024246235207>

我是否真的用完了尾数位，或者我可以使用更智能的算法做得更好吗？

这里的问题是我没有简单的方法来确定何时将顶点设置为非常接近另一条线。我不介意移动它，甚至完全摧毁它，因为这些都不会搞砸我的排序！

Answer 1

如果您将浮点中间值（ua、ub、denom）更改为双精度并打印 ua 值（除法后），您将得到这些：

0x1.f864ab6b36458p-1 in the first case
0x1.f864af01f2037p-1 in the second case

我已将它们以十六进制打印，以便于查看这些位。这两个值在前 22 位中一致（1.f864a 加上 b 和 f 的高位）。浮点数只有 23 位有效位！毫不奇怪，如果您以浮点数计算中间值，它们会四舍五入得到相同的答案。

在这种情况下，您也许可以通过使用双精度而不是浮点数计算中间值来解决该问题。 （我的观点是 struct 对 x 和 y 使用双精度。如果您使用浮点数，我不知道计算双精度中的中间值是否会有所帮助。）

但是，如果垂直线段更靠近两个水平线段的交点，则可能仍然需要比双精度更高的精度。

如果垂直线段正好通过水平线段的共享端点，你会怎么做？我假设你正确地处理了这个案例。

如果您查看使用双精度计算的 ub 值（除法后），您会得到这些：

0.9999960389052315
5.904388076838819e-06

这意味着交叉点非常非常接近水平线段的共享端点。

所以这就是我认为你可以做的事情。每次计算交集时，请查看 ub。如果它足够接近 1，移动端点成为交点。然后将交叉点视为一开始就是一个精确的通过案例。您实际上不必更改数据，但您必须将共享端点的两条线段视为移动端点，这意味着当前相交测试和下一条线段上的测试。