Graham Scan

At around the same time of the Jarvis March, R. L. Graham was also developing an algorithm to find the convex hull of a random set of points [1]. Unlike the Jarvis March, which is an operation, the Graham Scan is , where is the number of points and is the size for the hull. This means that the complexity of the Graham Scan is not output-sensitive; moreover, there are some cases where the Jarvis March is more optimal, depending on the size of the hull and the number of points to wrap.

Rather than starting at the leftmost point like the Jarvis March, the Graham scan starts at the bottom. We then sort the distribution of points based on the angle between the bottom-most point, the origin, and each other point. After sorting, we go through point-by-point, searching for points that are on the convex hull and throwing out any other points. We do this by looking for counter-clockwise rotations. If an angle between three points turns inward, the shape is obviously not convex, so we can throw that result out. We can find whether a rotation is counter-clockwise with trigonometric functions or by using a cross-product, like so:

function ccw(a::Point, b::Point, c::Point)
    return ((b.x - a.x)*(c.y - a.y) - (b.y - a.y)*(c.x - a.x))
end
ccw :: Point -> Point -> Point -> Double
ccw (xa, ya) (xb, yb) (xc, yc) = (xb - xa) * (yc - ya) - (yb - ya) * (xc - xa)
double ccw(struct point a, struct point b, struct point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
function ccw(a, b, c) {
  return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
}
#Is the turn counter clockwise?
def counter_clockwise(p1, p2, p3):
    return (p3[1]-p1[1])*(p2[0]-p1[0]) >= (p2[1]-p1[1])*(p3[0]-p1[0])
func counterClockwise(p1, p2, p3 point) bool {
    return (p3.y-p1.y)*(p2.x-p1.x) >= (p2.y-p1.y)*(p3.x-p1.x)
}
static double ccw(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

If the output of this function is 0, the points are collinear. If the output is positive, then the points form a counter-clockwise "left" turn. If the output is negative, then the points form a clockwise "right" turn. We basically do not want clockwise rotations, because this means we are at an interior angle.

To save memory and expensive append() operations, we ultimately look for points that should be on the hull and swap them with the first elements in the array. If there are elements on the hull, then the first elements in our output random distribution of points will be the hull. In the end, the code should look something like this:

function graham_scan!(points::Vector{Point})
    N = length(points)

    # Place the lowest point at the start of the array
    sort!(points, by = item -> item.y)

    # Sort all other points according to angle with that point
    other_points = sort(points[2:end], by = item -> atan(item.y - points[1].y,
                                                         item.x - points[1].x))

    # Place points sorted by angle back into points vector
    for i in 1:length(other_points)
        points[i+1] = other_points[i]
    end

    # M will be the point on the hull
    M = 2
    for i = 1:N
        while (ccw(points[M-1], points[M], points[i]) <= 0)
            if (M > 2)
                M -= 1
            # All points are collinear
            elseif (i == N)
                break
            else
                i += 1
            end
        end

        # ccw point found, updating hull and swapping points
        M += 1
        points[i], points[M] = points[M], points[i]
    end

    return points[1:M]
end
grahamScan :: [Point] -> [Point]
grahamScan [] = []
grahamScan pts = wrap sortedPts [p0]
  where [email protected](x, y)= minimumBy (compare `on` snd) pts
        sortedPts = sortOn (\(px, py) -> atan2 (py-y) (px-x) ) $ filter (/=p0) pts
        wrap [] ps = ps
        wrap (s:ss) [p] = wrap ss [s, p]
        wrap (s:ss) (p1:p2:ps)
          | ccw s p1 p2 > 0 = wrap (s:ss) (p2:ps)
          | otherwise       = wrap ss (s:p1:p2:ps)
size_t graham_scan(struct point *points, size_t size) {
    qsort(points, size, sizeof(struct point), cmp_points);
    polar_angles_sort(points, points[0], size);

    struct point tmp_points[size + 1];
    memcpy(tmp_points + 1, points, size * sizeof(struct point));
    tmp_points[0] = tmp_points[size];

    size_t m = 1;
    for (size_t i = 2; i <= size; ++i) {
        while (ccw(tmp_points[m - 1], tmp_points[m], tmp_points[i]) <= 0) {
            if (m > 1) {
                m--;
                continue;
            } else if (i == size) {
                break;
            } else {
                i++;
            }
        }

        m++;
        struct point tmp = tmp_points[i];
        tmp_points[i] = tmp_points[m];
        tmp_points[m] = tmp;
    }

    memcpy(points, tmp_points + 1, size * sizeof(struct point));

    return m;
}
function grahamScan(points) {
  // First, sort the points so the one with the lowest y-coordinate comes first (the pivot)
  points = [...points].sort((a, b) => (a.y - b.y));
  const pivot = points[0];

  // Then sort all remaining points based on the angle between the pivot and itself
  const hull = points.slice(1).sort((a, b) => polarAngle(a, pivot) - polarAngle(b, pivot));

  // The pivot is always on the hull
  hull.unshift(pivot);

  let n = hull.length;
  let m = 1;
  for (let i = 2; i < n; i++) {
    while (ccw(hull[m - 1], hull[m], hull[i]) <= 0) {
      if (m > 1) {
        m -= 1;
      } else if (m === i) {
        break;
      } else {
        i += 1;
      }
    }

    m += 1;
    [hull[i], hull[m]] = [hull[m], hull[i]];
  }

  return hull.slice(0, m + 1);
}
def graham_scan(gift):
    gift = list(set(gift)) #Remove duplicate points
    start = min(gift, key=lambda p: (p[1], p[0])) #Must be in hull
    gift.remove(start)

    s = sorted(gift,key=lambda point: polar_angle(start, point))
    hull = [start,s[0],s[1]]

    #Remove points from hull that make the hull concave
    for pt in s[2:]:
        while not counter_clockwise(hull[-2], hull[-1], pt):
            del hull[-1]
        hull.append(pt)
func grahamScan(points []point) []point {
    sort.Slice(points, func(a, b int) bool {
        return points[a].y < points[b].y || (points[a].y == points[b].y && points[a].x < points[b].x)
    })

    start := points[0]
    points = points[1:]

    sort.Slice(points, func(a, b int) bool {
        return polarAngle(start, points[a]) < polarAngle(start, points[b])
    })

    hull := []point{start, points[0], points[1]}
    for _, p := range points[2:] {
        for !counterClockwise(hull[len(hull)-2], hull[len(hull)-1], p) {
            hull = hull[:len(hull)-1]
        }
        hull = append(hull, p)
    }

    return hull
}
static List<Point> grahamScan(List<Point> gift) {
    gift = gift.stream()
               .distinct()
               .sorted(Comparator.comparingDouble(point -> -point.y))
               .collect(Collectors.toList());

    Point pivot = gift.get(0);

    // Sort the remaining Points based on the angle between the pivot and itself
    List<Point> hull = gift.subList(1, gift.size());
    hull.sort(Comparator.comparingDouble(point -> polarAngle(point, pivot)));

    // The pivot is always on the hull
    hull.add(0, pivot);

    int n = hull.size();
    int m = 1;

    for (int i = 2; i < n; i++) {
        while (ccw(hull.get(m - 1), hull.get(m), hull.get(i)) <= 0) {
            if (m > 1) {
                m--;
            } else if (m == 1) {
                break;
            } else {
                i++;
            }
        }
        m++;

        Point temp = hull.get(i);
        hull.set(i, hull.get(m));
        hull.set(m, temp);
    }
    return hull.subList(0, m + 1);
}

Bibliography

1Graham, Ronald L, An efficient algorithm for determining the convex hull of a finite planar set, Elsevier, 1972.

Example Code

struct Point
    x::Float64
    y::Float64
end

function ccw(a::Point, b::Point, c::Point)
    return ((b.x - a.x)*(c.y - a.y) - (b.y - a.y)*(c.x - a.x))
end

function graham_scan!(points::Vector{Point})
    N = length(points)

    # Place the lowest point at the start of the array
    sort!(points, by = item -> item.y)

    # Sort all other points according to angle with that point
    other_points = sort(points[2:end], by = item -> atan(item.y - points[1].y,
                                                         item.x - points[1].x))

    # Place points sorted by angle back into points vector
    for i in 1:length(other_points)
        points[i+1] = other_points[i]
    end

    # M will be the point on the hull
    M = 2
    for i = 1:N
        while (ccw(points[M-1], points[M], points[i]) <= 0)
            if (M > 2)
                M -= 1
            # All points are collinear
            elseif (i == N)
                break
            else
                i += 1
            end
        end

        # ccw point found, updating hull and swapping points
        M += 1
        points[i], points[M] = points[M], points[i]
    end

    return points[1:M]
end

function main()
    # This hull is just a simple test so we know what the output should be
    points = [
        Point(-5.,2), Point(5,7), Point(-6,-12), Point(-14,-14), Point(9,9),
        Point(-1,-1), Point(-10,11), Point(-6,15), Point(-6,-8), Point(15,-9),
        Point(7,-7), Point(-2,-9), Point(6,-5), Point(0,14), Point(2,8)
    ]
    hull = graham_scan!(points)
    println(hull)
end

main()
import Data.List (sortOn, minimumBy)
import Data.Function (on)

type Point = (Double, Double)

ccw :: Point -> Point -> Point -> Double
ccw (xa, ya) (xb, yb) (xc, yc) = (xb - xa) * (yc - ya) - (yb - ya) * (xc - xa)

grahamScan :: [Point] -> [Point]
grahamScan [] = []
grahamScan pts = wrap sortedPts [p0]
  where [email protected](x, y)= minimumBy (compare `on` snd) pts
        sortedPts = sortOn (\(px, py) -> atan2 (py-y) (px-x) ) $ filter (/=p0) pts
        wrap [] ps = ps
        wrap (s:ss) [p] = wrap ss [s, p]
        wrap (s:ss) (p1:p2:ps)
          | ccw s p1 p2 > 0 = wrap (s:ss) (p2:ps)
          | otherwise       = wrap ss (s:p1:p2:ps)

main = do
  -- We build the set of points of integer coordinates within a circle of radius 5
  let pts = [(x,y) | x<-[-5..5], y<-[-5..5], x^2+y^2<=5^2]
  -- And extract the convex hull
  print $ grahamScan pts
#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

struct point {
    double x, y;
};

int cmp_points(const void *a, const void *b) {
    struct point* pa = (struct point*) a;
    struct point* pb = (struct point*) b;

    if (pa->y > pb->y) {
        return 1;
    } else if (pa->y < pb->y) {
        return -1;
    } else {
        return 0;
    }
}

double ccw(struct point a, struct point b, struct point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

double polar_angle(struct point origin, struct point p) {
    return atan2(p.y - origin.y, p.x - origin.x);
}

void polar_angles_sort(struct point *points, struct point origin, size_t size) {
    if (size < 2) {
        return;
    }

    double pivot_angle = polar_angle(origin, points[size / 2]);

    int i = 0;
    int j = size - 1;
    while (1) {
        while (polar_angle(origin, points[i]) < pivot_angle) {
            i++;
        }
        while (polar_angle(origin, points[j]) > pivot_angle) {
            j--;
        }

        if (i >= j) {
            break;
        }

        struct point tmp = points[i];
        points[i] = points[j];
        points[j] = tmp;

        i++;
        j--;
    }

    polar_angles_sort(points, origin, i);
    polar_angles_sort(points + i, origin, size - i);
}

size_t graham_scan(struct point *points, size_t size) {
    qsort(points, size, sizeof(struct point), cmp_points);
    polar_angles_sort(points, points[0], size);

    struct point tmp_points[size + 1];
    memcpy(tmp_points + 1, points, size * sizeof(struct point));
    tmp_points[0] = tmp_points[size];

    size_t m = 1;
    for (size_t i = 2; i <= size; ++i) {
        while (ccw(tmp_points[m - 1], tmp_points[m], tmp_points[i]) <= 0) {
            if (m > 1) {
                m--;
                continue;
            } else if (i == size) {
                break;
            } else {
                i++;
            }
        }

        m++;
        struct point tmp = tmp_points[i];
        tmp_points[i] = tmp_points[m];
        tmp_points[m] = tmp;
    }

    memcpy(points, tmp_points + 1, size * sizeof(struct point));

    return m;
}

int main() {
    struct point points[] = {{2.0, 1.9}, {1.0, 1.0}, {2.0, 4.0}, {3.0, 1.0},
                                {2.0, 0.0}};

    printf("Points:\n");
    for (size_t i = 0; i < 5; ++i) {
        printf("(%f,%f)\n", points[i].x, points[i].y);
    }

    size_t hull_size = graham_scan(points, 5);

    printf("\nHull:\n");
    for (size_t i = 0; i < hull_size; ++i) {
        printf("(%f,%f)\n", points[i].x, points[i].y);
    }

    return 0;
}
function grahamScan(points) {
  // First, sort the points so the one with the lowest y-coordinate comes first (the pivot)
  points = [...points].sort((a, b) => (a.y - b.y));
  const pivot = points[0];

  // Then sort all remaining points based on the angle between the pivot and itself
  const hull = points.slice(1).sort((a, b) => polarAngle(a, pivot) - polarAngle(b, pivot));

  // The pivot is always on the hull
  hull.unshift(pivot);

  let n = hull.length;
  let m = 1;
  for (let i = 2; i < n; i++) {
    while (ccw(hull[m - 1], hull[m], hull[i]) <= 0) {
      if (m > 1) {
        m -= 1;
      } else if (m === i) {
        break;
      } else {
        i += 1;
      }
    }

    m += 1;
    [hull[i], hull[m]] = [hull[m], hull[i]];
  }

  return hull.slice(0, m + 1);
}

function polarAngle(a, b) {
  return Math.atan2(a.y - b.y, a.x - b.x);
}

function ccw(a, b, c) {
  return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
}

const points = [
  { x: 1, y: 3 },
  { x: 2, y: 4 },
  { x: 4, y: 0 },
  { x: 1, y: 0 },
  { x: 0, y: 2 },
  { x: 2, y: 2 },
  { x: 3, y: 4 },
  { x: 3, y: 1 },
];

const convexHull = grahamScan(points);
convexHull.forEach(p => console.log(`(${p.x}, ${p.y})`));
from math import atan2


#Is the turn counter clockwise?
def counter_clockwise(p1, p2, p3):
    return (p3[1]-p1[1])*(p2[0]-p1[0]) >= (p2[1]-p1[1])*(p3[0]-p1[0])


#Find the polar angle of a point relative to a reference point
def polar_angle(ref, point):
    return atan2(point[1]-ref[1],point[0]-ref[0])


def graham_scan(gift):
    gift = list(set(gift)) #Remove duplicate points
    start = min(gift, key=lambda p: (p[1], p[0])) #Must be in hull
    gift.remove(start)

    s = sorted(gift,key=lambda point: polar_angle(start, point))
    hull = [start,s[0],s[1]]

    #Remove points from hull that make the hull concave
    for pt in s[2:]:
        while not counter_clockwise(hull[-2], hull[-1], pt):
            del hull[-1]
        hull.append(pt)

    return hull


def main():
    test_gift = [(-5, 2), (5, 7), (-6, -12), (-14, -14), (9, 9),
                (-1, -1), (-10, 11), (-6, 15), (-6, -8), (15, -9),
                (7, -7), (-2, -9), (6, -5), (0, 14), (2, 8)]
    hull = graham_scan(test_gift)

    print("The points in the hull are:")
    for point in hull:
        print(point)


main()
package main

import (
    "fmt"
    "math"
    "sort"
)

type point struct {
    x, y int
}

func counterClockwise(p1, p2, p3 point) bool {
    return (p3.y-p1.y)*(p2.x-p1.x) >= (p2.y-p1.y)*(p3.x-p1.x)
}

func polarAngle(ref, point point) float64 {
    return math.Atan2(float64(point.y-ref.y), float64(point.x-ref.x))
}

func grahamScan(points []point) []point {
    sort.Slice(points, func(a, b int) bool {
        return points[a].y < points[b].y || (points[a].y == points[b].y && points[a].x < points[b].x)
    })

    start := points[0]
    points = points[1:]

    sort.Slice(points, func(a, b int) bool {
        return polarAngle(start, points[a]) < polarAngle(start, points[b])
    })

    hull := []point{start, points[0], points[1]}
    for _, p := range points[2:] {
        for !counterClockwise(hull[len(hull)-2], hull[len(hull)-1], p) {
            hull = hull[:len(hull)-1]
        }
        hull = append(hull, p)
    }

    return hull
}

func main() {
    points := []point{{-5, 2}, {5, 7}, {-6, -12}, {-14, -14}, {9, 9},
        {-1, -1}, {-10, 11}, {-6, 15}, {-6, -8}, {15, -9},
        {7, -7}, {-2, -9}, {6, -5}, {0, 14}, {2, 8}}

    fmt.Println("The points in the hull are:")
    hull := grahamScan(points)
    for _, p := range hull {
        fmt.Println(p)
    }
}
import java.util.ArrayList;
import java.util.Comparator;
import java.util.List;
import java.util.stream.Collectors;

public class GrahamScan {

    static class Point {
        public double x;
        public double y;

        public Point(double x, double y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public boolean equals(Object o) {
            if (o == null) return false;
            if (o == this) return true;
            if (!(o instanceof Point)) return false;
            Point p = (Point)o;
            return p.x == this.x && p.y == this.y;
        }
    }

    static double ccw(Point a, Point b, Point c) {
        return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
    }

    static double polarAngle(Point origin, Point p) {
        return Math.atan2(p.y - origin.y, p.x - origin.x);
    }

    static List<Point> grahamScan(List<Point> gift) {
        gift = gift.stream()
                   .distinct()
                   .sorted(Comparator.comparingDouble(point -> -point.y))
                   .collect(Collectors.toList());

        Point pivot = gift.get(0);

        // Sort the remaining Points based on the angle between the pivot and itself
        List<Point> hull = gift.subList(1, gift.size());
        hull.sort(Comparator.comparingDouble(point -> polarAngle(point, pivot)));

        // The pivot is always on the hull
        hull.add(0, pivot);

        int n = hull.size();
        int m = 1;

        for (int i = 2; i < n; i++) {
            while (ccw(hull.get(m - 1), hull.get(m), hull.get(i)) <= 0) {
                if (m > 1) {
                    m--;
                } else if (m == 1) {
                    break;
                } else {
                    i++;
                }
            }
            m++;

            Point temp = hull.get(i);
            hull.set(i, hull.get(m));
            hull.set(m, temp);
        }
        return hull.subList(0, m + 1);
    }

    public static void main(String[] args) {
        ArrayList<Point> points = new ArrayList<>();

        points.add(new Point(-5, 2));
        points.add(new Point(5, 7));
        points.add(new Point(-6, -12));
        points.add(new Point(-14, -14));
        points.add(new Point(9, 9));
        points.add(new Point(-1, -1));
        points.add(new Point(-10, 11));
        points.add(new Point(-6, 15));
        points.add(new Point(-6, -8));
        points.add(new Point(15, -9));
        points.add(new Point(7, -7));
        points.add(new Point(-2, -9));
        points.add(new Point(6, -5));
        points.add(new Point(0, 14));
        points.add(new Point(2, 8));

        List<Point> convexHull = grahamScan(points);

        convexHull.forEach(p -> System.out.printf("% 1.0f, % 1.0f\n", p.x, p.y));
    }
}

results matching ""

    No results matching ""