Spatial algorithms under the hood: Graham Scan

Last week at the Recurse Center, my self-proclaimed goal was to better understand how “gift wrapping” algorithms work. In computational geometry, gift wrapping algorithms are a family of algorithms that compute the convex hull of a given set of points. What is a convex hull? Basically, it is a shape that encloses a set of points, such that no corners of the shape are bent inwards, hence the name “gift wrapping”. But why would you want to find this shape? Although the technique seems to have a variety of applications outside of GIS, such as image processing, I am most familiar with GIS applications. For example, a convex hull can be used as a way to better describe patterns, such as animal movements that were collected as point features.

This blog post will focus on the Graham Scan algorithm as a way of finding the convex hull of a set of points. Let’s go through the steps needed:

1. 1) Given a set of points, select point with lowest y value (Fig. 1)
2. 2) Calculate angles between the lowest point and all the other points (Fig. 1)
3. 3) The size of the angle will determine in which order to iterate through points. Thus, sort points by angle relative to lowest point in ascending order
4. 4) Iterate through points. Add a point to an output array if there is a counterclockwise turn to previous point. To figure out if a turn is clockwise or counterclockwise, we can leverage the cross product of two vectors (Fig. 2 and 3).
5. 5) If a point is located clockwise relative to the previous point, we pop it off the output array.
6. 6) When we are done with scanning the points, return the result (Fig. 4)

Step 4 is arguably the most difficult part, if you, like me, don't remember Trigonometry class. To get the right crossproduct, we will have to recalculate the vectors each time, v1 going from P _n-1 to P _n and v2 going from P _n to P _n+1 . The result of the crossproduct is the vector z, that I've drawn poorly in a small sketch next to the main graphic in Fig. 2 and Fig. 3. As you can see, the z vector points either up (+) or down (-), indicating a counterclockwise move if positive and a clockwise move if negative. You can find a code snippet for the main function below. The full code including the helper functions for the crossproduct isCcw() and the angle sort sortAngles() can be found in the Github repo.

def GrahamScan(data):
    #Select point with lowest y and append it to output array
    gift = []
    minPoint = min(data, key = lambda x : x[1])
    gift.append(minPoint)

    #get array of sorted angles
    angles = sortAngles(data, minPoint)

    for i in range(len(angles)):
        p = angles[i].coords
        while (len(gift) >= 2) and (isCcw(gift[-2], gift[-1], p) < 0): 
            gift.pop()  #isCcw() returns negative results, therefore remove point from output array
        gift.append(p)

    return gift

It is probably also a good idea to look at the time complexity of the Graham Scan: Most sorting algorithms take n log(n) time, while the actual scan takes n time. We select the dominant function, therefore the overall time complexity is O(n log (n)). How could we optimize the Graham scan? We could for example wipe out points in the interior that we for sure know are not in the convex hull by finding the farthest points in the SW, NW, NE, and SE direction and eliminte the points that are enclosed by these 4 points. I will try to update the Repo with an example the above optimization process in the next weeks!

Next up: The Traveling Salesperson Problem