Graham scan issue at high amount of points

Question

I have an issue in graham scan algorithm when my list have a lot of points, but works every time fine with low amount of points. I made some screenshots :

working : (300 points)

not working (5000 points)

angle calculation:

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

}

sorting points by polar angle depend on max Y point :

//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;
                }
            }   
        }

ccw method :

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));
}

Graham scan algorithm :

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;

    }

I really don't know why my program explode on 800+ points... It's hard to debug, because algorithm works very well with 300,400,500... points.

Ty for informations.

Answer 1

The algorithm on Wikipedia is broken. It doesn't handle the case where multiple points are collinear with each other and with the minimum point. For instance, the following test case will fail:

        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),

        };

The issue is that, when marching along the points, it may be necessary to discard the current point rather than either extending the hull or replacing the last point on the hull, if the point is between the last two points in the hull. (This could only happen if the points were also collinear with the minimum point, otherwise the preceding angle sorting would prevent such double-backs.) There's no logic for that in the psuedocode shown.

I also think the Wikipedia algorithm may be handling floating point roundoff errors badly. In particular checking ccw <= 0 looks problematic.

Here's an attempt to clean up the Wikipedia algorithm. I had to get rid of the (vaguely kludgy) "sentinal point" because, in case all the points are horizontally aligned, it would get chosen essentially randomly:

    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;
    }

By the way, your algorithm is order n-squared in a couple places:

1. The bubblesort is order O(N^2)
2. Removing points while finding the hull rather than swapping them to the front of the list can be O(N^2) because all subsequent points have to be shifted down.

Let me know if it solves your problems, I've tested it a bit but not completely.