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

def counter_clockwise(p1, p2, p3):
"""Is the turn counter-clockwise?"""
return (p3 - p1) * (p2 - p1) >= (p2 - p1) * (p3 - p1)

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

double ccw(const point& a, const point& b, const 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.y,
item.x - points.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 p0@(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, size); struct point tmp_points[size + 1]; memcpy(tmp_points + 1, points, size * sizeof(struct point)); tmp_points = 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; // 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, p)) # Must be in hull gift.remove(start) s = sorted(gift, key=lambda point: polar_angle(start, point)) hull = [start, s, s] # 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  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 points = points[1:] sort.Slice(points, func(a, b int) bool { return polarAngle(start, points[a]) < polarAngle(start, points[b]) }) hull := []point{start, points, points} 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); }  std::vector<point> graham_scan(std::vector<point>& points) { // selecting lowest point as pivot size_t low_index = 0; for (size_t i = 1; i < points.size(); i++) { if (points[i].y < points[low_index].y) { low_index = i; } } std::swap(points, points[low_index]); point pivot = points; // sorting points by polar angle std::sort( points.begin() + 1, points.end(), [&pivot](const point& pa, const point& pb) { return polar_angle(pivot, pa) < polar_angle(pivot, pb); }); // creating convex hull size_t m = 1; for (size_t i = 2; i < points.size(); i++) { while (ccw(points[m - 1], points[m], points[i]) <= 0) { if (m > 1) { m--; continue; } else if (i == points.size()) { break; } else { i++; } } m++; std::swap(points[i], points[m]); } return std::vector<point>(points.begin(), points.begin() + m + 1); }  ### Bibliography  1 Graham, 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.y, item.x - points.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 p0@(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, size); struct point tmp_points[size + 1]; memcpy(tmp_points + 1, points, size * sizeof(struct point)); tmp_points = 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[] = {{-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}}; size_t num_initial_points = 15; printf("Points:\n"); for (size_t i = 0; i < num_initial_points; ++i) { printf("(%f,%f)\n", points[i].x, points[i].y); } size_t hull_size = graham_scan(points, num_initial_points); 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; // 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: -5, y: 2 }, { x: 5, y: 7 }, { x: -6, y: -12 }, { x: -14, y: -14 }, { x: 9, y: 9 }, { x: -1, y: -1 }, { x: -10, y: 11 }, { x: -6, y: 15 }, { x: -6, y: -8 }, { x: 15, y: -9 }, { x: 7, y: -7 }, { x: -2, y: -9 }, { x: 6, y: -5 }, { x: 0, y: 14 }, { x: 2, y: 8 }, ]; const convexHull = grahamScan(points); console.log("The points in the hull are:"); convexHull.forEach(p => console.log((${p.x}, \${p.y})));

from math import atan2

def counter_clockwise(p1, p2, p3):
"""Is the turn counter-clockwise?"""
return (p3 - p1) * (p2 - p1) >= (p2 - p1) * (p3 - p1)

def polar_angle(ref, point):
"""Find the polar angle of a point relative to a reference point"""
return atan2(point - ref, point - ref)

start = min(gift, key=lambda p: (p, p))  # Must be in hull

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

# 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():
(-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),
]

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
points = points[1:]

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

hull := []point{start, points, points}
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);
}

.distinct()
.sorted(Comparator.comparingDouble(point -> -point.y))
.collect(Collectors.toList());

// Sort the remaining Points based on the angle between the pivot and itself
hull.sort(Comparator.comparingDouble(point -> polarAngle(point, pivot)));

// The pivot is always on the hull

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

List<Point> convexHull = grahamScan(points);

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

#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>

struct point {
double x;
double y;
};

std::ostream& operator<<(std::ostream& os, const std::vector<point>& points) {
for (auto p : points) {
os << "(" << p.x << ", " << p.y << ")\n";
}
return os;
}

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

double polar_angle(const point& origin, const point& p) {
return std::atan2(p.y - origin.y, p.x - origin.x);
}

std::vector<point> graham_scan(std::vector<point>& points) {
// selecting lowest point as pivot
size_t low_index = 0;
for (size_t i = 1; i < points.size(); i++) {
if (points[i].y < points[low_index].y) {
low_index = i;
}
}
std::swap(points, points[low_index]);
point pivot = points;

// sorting points by polar angle
std::sort(
points.begin() + 1,
points.end(),
[&pivot](const point& pa, const point& pb) {
return polar_angle(pivot, pa) < polar_angle(pivot, pb);
});

// creating convex hull
size_t m = 1;
for (size_t i = 2; i < points.size(); i++) {
while (ccw(points[m - 1], points[m], points[i]) <= 0) {
if (m > 1) {
m--;
continue;
} else if (i == points.size()) {
break;
} else {
i++;
}
}
m++;
std::swap(points[i], points[m]);
}
return std::vector<point>(points.begin(), points.begin() + m + 1);
}

int main() {
std::vector<point> points = {{-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}};
std::cout << "original points are as follows:\n" << points;
const std::vector<point> hull = graham_scan(points);
std::cout << "points in hull are as follows:\n" << hull;
return 0;
}


##### Text 