c# - 格雷厄姆扫描问题在高点

Question

当我的列表有很多点时，我在格雷厄姆扫描算法中遇到了问题，但每次都可以很好地处理少量点。我做了一些截图:

工作:(300 分)

不工作(5000 点)

角度计算:

public static double angle(MyVector3D vec1, MyVector3D vec2)
{
    return Math.Atan2(vec2.Y - vec1.Y, vec2.X - vec1.X) * 180 / Math.PI;

}

根据最大 Y 点按极角排序点:

//bubblesort
    private void sortList()
    {
        double temp = 0.0;
        MyVector3D tempVector = new MyVector3D();
        for (int i = 1; i < points.Count; i++)
        {
            for (int j = 1; j < points.Count - 1; j++)
            {
                if (angles[j] < angles[j + 1])
                {
                    temp = angles[j + 1];
                    tempVector = points[j + 1];
                    angles[j + 1] = angles[j];
                    points[j + 1] = points[j];
                    angles[j] = temp;
                    points[j] = tempVector;
                }
            }   
        }

逆时针方法:

private double ccw(MyVector3D vec1, MyVector3D vec2, MyVector3D vec3)
{
    // ccwTest = ((vec2.X - vec1.X) * (vec3.Y - vec1.Y)) - ((vec2.Y - vec1.Y) * (vec3.X - vec1.X));
    return ((vec2.X - vec1.X) * (vec3.Y - vec1.Y)) - ((vec2.Y - vec1.Y) * (vec3.X - vec1.X));
}

格雷厄姆扫描算法:

for (int i = 2; i < points.Count; i++)
    {
        while (ccw(points[M - 1], points[M], points[i]) > 0)
        {
            if (M > 1)
            {
                points.RemoveAt(M);
                M -= 1;
                i--;
            }
            else if (i == points.Count - 1)
            {
                break;
            }

            else
            {
                i += 1;
            }
        }
        //goodPoints.Add(points[M]);
        //listBoxInfos.Items.Add("g" + (int)points[M].X + "," + (int)points[M].Y + "," + 0);
        //listBoxInfos.Items.Add("ccw" + ccwTest);
        M += 1;

    }

我真的不知道为什么我的程序会在 800+ 点上爆炸...很难调试，因为算法在 300,400,500...点上工作得很好。

需要信息。

Answer 1

维基百科上的算法被破坏了。它不处理多个点彼此共线以及与最小点共线的情况。例如，以下测试用例将失败:

        Vector3[] points3 = new Vector3[] 
        {
            new Vector3( 1, 1, 0),
            new Vector3( 5, 5, 0),
            new Vector3( 3, 3, 0),
            new Vector3( 2, 2, 0),
            new Vector3( 1, 1, 0),
            new Vector3( 1, 10, 0),

        };

问题是，当沿着点行进时，可能需要丢弃当前点而不是扩展船体或替换船体上的最后一个点，如果该点在船体中的最后两个点之间. (这只有在点也与最小点共线时才会发生，否则前面的角度排序会阻止这种双回。)在显示的伪代码中没有逻辑。

我还认为 Wikipedia 算法可能无法很好地处理浮点舍入误差。特别是检查 ccw <= 0 看起来有问题。

这是清理维基百科算法的尝试。我不得不摆脱(含糊不清的)“哨兵点”，因为如果所有点都水平对齐，它基本上会随机选择:

    public static IList<Vector3> GrahamScanCompute(IList<Vector3> initialPoints)
    {
        if (initialPoints.Count < 2)
            return initialPoints.ToList();

        // Find point with minimum y; if more than one, minimize x also.
        int iMin = Enumerable.Range(0, initialPoints.Count).Aggregate((jMin, jCur) =>
        {
            if (initialPoints[jCur].Y < initialPoints[jMin].Y)
                return jCur;
            if (initialPoints[jCur].Y > initialPoints[jMin].Y)
                return jMin;
            if (initialPoints[jCur].X < initialPoints[jMin].X)
                return jCur;
            return jMin;
        });
        // Sort them by polar angle from iMin, 
        var sortQuery = Enumerable.Range(0, initialPoints.Count)
            .Where((i) => (i != iMin)) // Skip the min point
            .Select((i) => new KeyValuePair<double, Vector3>(Math.Atan2(initialPoints[i].Y - initialPoints[iMin].Y, initialPoints[i].X - initialPoints[iMin].X), initialPoints[i]))
            .OrderBy((pair) => pair.Key)
            .Select((pair) => pair.Value);
        List<Vector3> points = new List<Vector3>(initialPoints.Count);
        points.Add(initialPoints[iMin]);     // Add minimum point
        points.AddRange(sortQuery);          // Add the sorted points.

        int M = 0;
        for (int i = 1, N = points.Count; i < N; i++)
        {
            bool keepNewPoint = true;
            if (M == 0)
            {
                // Find at least one point not coincident with points[0]
                keepNewPoint = !NearlyEqual(points[0], points[i]);
            }
            else
            {
                while (true)
                {
                    var flag = WhichToRemoveFromBoundary(points[M - 1], points[M], points[i]);
                    if (flag == RemovalFlag.None)
                        break;
                    else if (flag == RemovalFlag.MidPoint)
                    {
                        if (M > 0)
                            M--;
                        if (M == 0)
                            break;
                    }
                    else if (flag == RemovalFlag.EndPoint)
                    {
                        keepNewPoint = false;
                        break;
                    }
                    else
                        throw new Exception("Unknown RemovalFlag");
                }
            }
            if (keepNewPoint)
            {
                M++;
                Swap(points, M, i);
            }
        }
        // points[M] is now the last point in the boundary.  Remove the remainder.
        points.RemoveRange(M + 1, points.Count - M - 1);
        return points;
    }

    static void Swap<T>(IList<T> list, int i, int j)
    {
        if (i != j)
        {
            T temp = list[i];
            list[i] = list[j];
            list[j] = temp;
        }
    }

    public static double RelativeTolerance { get { return 1e-10; } }

    public static bool NearlyEqual(Vector3 a, Vector3 b)
    {
        return NearlyEqual(a.X, b.X) && NearlyEqual(a.Y, b.Y);
    }

    public static bool NearlyEqual(double a, double b)
    {
        return NearlyEqual(a, b, RelativeTolerance);
    }

    public static bool NearlyEqual(double a, double b, double epsilon)
    {
        // See here: http://floating-point-gui.de/errors/comparison/
        if (a == b)
        { // shortcut, handles infinities
            return true;
        }

        double absA = Math.Abs(a);
        double absB = Math.Abs(b);
        double diff = Math.Abs(a - b);
        double sum = absA + absB;
        if (diff < 4*double.Epsilon || sum < 4*double.Epsilon)
            // a or b is zero or both are extremely close to it
            // relative error is less meaningful here
            return true;

        // use relative error
        return diff / (absA + absB) < epsilon;
    }

    static double CCW(Vector3 p1, Vector3 p2, Vector3 p3)
    {
        // Compute (p2 - p1) X (p3 - p1)
        double cross1 = (p2.X - p1.X) * (p3.Y - p1.Y);
        double cross2 = (p2.Y - p1.Y) * (p3.X - p1.X);
        if (NearlyEqual(cross1, cross2))
            return 0;
        return cross1 - cross2;
    }

    enum RemovalFlag
    {
        None,
        MidPoint,
        EndPoint
    };

    static RemovalFlag WhichToRemoveFromBoundary(Vector3 p1, Vector3 p2, Vector3 p3)
    {
        var cross = CCW(p1, p2, p3);
        if (cross < 0)
            // Remove p2
            return RemovalFlag.MidPoint;
        if (cross > 0)
            // Remove none.
            return RemovalFlag.None;
        // Check for being reversed using the dot product off the difference vectors.
        var dotp = (p3.X - p2.X) * (p2.X - p1.X) + (p3.Y - p2.Y) * (p2.Y - p1.Y);
        if (NearlyEqual(dotp, 0.0))
            // Remove p2
            return RemovalFlag.MidPoint;
        if (dotp < 0)
            // Remove p3
            return RemovalFlag.EndPoint;
        else
            // Remove p2
            return RemovalFlag.MidPoint;
    }

顺便说一句，您的算法在几个地方是 n 阶的:

1. 冒泡排序的顺序是 O(N^2)
2. 在找到船体时删除点而不是将它们交换到列表的前面可以是 O(N^2) 因为所有后续点都必须向下移动。

如果它能解决您的问题，请告诉我，我已经对其进行了一些测试，但并不完全。