C++

Convex Hulls: Graham Scan, Monotone Chain and QuickHull

Build the convex hull with Andrew's monotone chain, compare it with Graham scan, Jarvis march and QuickHull by counting their work, and test every degenerate case.

The rubber band around a set of points

Stretch a rubber band around a set of nails in a board and let it go. The shape it snaps to is the convex hull of the nails: a convex polygon whose corners are some of the nails, with every other nail inside it or on its boundary.

The hull is a building block for many other computations. A collision test can first check the hulls of two objects, which are much simpler than the objects. The two points of a set that are farthest apart are always hull vertices, so only the hull needs searching. A cluster of GPS fixes is often summarized by the polygon around it. And computing the hull correctly is the classic first exercise in robust geometry, because every way of getting degenerate input wrong shows up in it: duplicates, collinear points, all points on one line.

This lesson implements four hull algorithms on the exact integer predicates from lesson 2: Andrew's monotone chain as the one to use, then Graham scan, Jarvis march and QuickHull for comparison. It counts the work each one does on different kinds of input, and tests all of them against a brute-force definition.

What exactly the hull is

A set is convex if, for any two of its points, the segment between them lies entirely in the set. The convex hull of a set of points is the smallest convex set that contains them all. For finitely many points in the plane it is a convex polygon (or, in degenerate cases, a segment or a single point), and its corners are input points. Those corners are the extreme points: the ones that are not inside the hull of the others.

An algorithm must return the hull in a well-defined form. This course uses one convention for polygons throughout:

  • the vertices are listed counter-clockwise,
  • the first vertex is the lexicographically smallest point (smallest x, then smallest y),
  • the first vertex is not repeated at the end.

The number of hull vertices is called h. It can be anything from 1 to n, and it matters for the running time of some algorithms. Points in a square drawn uniformly at random have few hull vertices; points on a circle are all hull vertices.

Two policy decisions remain, and both must be made explicitly. Duplicates: the same point given twice appears once in the output. Collinear boundary points: a point lying on a hull edge, between two corners, is either dropped (the hull is its corners only) or kept (the hull lists every input point on its boundary). Most uses want the corners only; some, such as reconstructing a boundary from survey points, want them all. The implementation takes a parameter.

Andrew's monotone chain

A. M. Andrew's 1979 algorithm splits the hull into a lower hull, running from the leftmost to the rightmost point along the bottom, and an upper hull, running back along the top. Sorted from left to right, the lower hull turns left at every corner. So:

  1. Sort the points lexicographically (by x, then y) and remove duplicates.
  2. Walk through them from left to right, keeping a stack. Before pushing a point p, look at the last two points on the stack: if they and p do not make a left turn, the last point cannot be a corner of the lower hull, so pop it. Repeat, then push p.
  3. Do the same from right to left for the upper hull.
  4. Join the two chains. Each ends at the point where the other starts, so drop one copy of each.

The only geometry is the orientation test from lessons 1 and 2. "Not a left turn" means orientation(...) <= 0 when collinear boundary points are dropped, and < 0 when they are kept. That one comparison is the policy parameter.

The program runs the algorithm on a small example with a trace of every push and pop. The input is the square with corners (0,0), (4,0), (4,4) and (0,4), plus (2,0) and (2,4) on two of its sides, (2,2) inside, and (4,0) a second time.

main.cppc++23
#include <algorithm>
#include <print>
#include <vector>

struct Point { long long x, y; };

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

bool trace = false;  // print every push and pop

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                if (trace) std::println("  pop  ({},{})", hull.back().x, hull.back().y);
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            if (trace) std::println("  push ({},{})", p.x, p.y);
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
        if (trace && pass == 0) std::println("  -- upper hull --");
    }
    return hull;
}

void show(const char* label, const std::vector<Point>& h) {
    std::print("{}", label);
    for (Point p : h) std::print(" ({},{})", p.x, p.y);
    std::println("");
}

int main() {
    // A square with a point on two of its sides, one inside, and a duplicate.
    std::vector<Point> pts{{0, 0}, {4, 0}, {2, 0}, {4, 4},
                           {0, 4}, {2, 2}, {2, 4}, {4, 0}};
    trace = true;
    std::println("lower hull, left to right:");
    auto hull = convexHull(pts);
    trace = false;
    show("vertices only:      ", hull);
    show("keeping collinear:  ", convexHull(pts, true));
}
Output

Press Run to compile and run this example.

lower hull, left to right:
  push (0,0)
  push (0,4)
  pop  (0,4)
  push (2,0)
  push (2,2)
  pop  (2,2)
  push (2,4)
  pop  (2,4)
  pop  (2,0)
  push (4,0)
  push (4,4)
  -- upper hull --
  push (4,4)
  push (4,0)
  pop  (4,0)
  push (2,4)
  push (2,2)
  pop  (2,2)
  push (2,0)
  pop  (2,0)
  pop  (2,4)
  push (0,4)
  push (0,0)
vertices only:       (0,0) (4,0) (4,4) (0,4)
keeping collinear:   (0,0) (2,0) (4,0) (4,4) (2,4) (0,4)

Read the lower hull part of the trace against the rule. (0,4) is pushed and popped at once: going from (0,0) up to (0,4) and then down to (2,0) is a right turn. (2,2) is popped when (2,4) arrives, because (2,0), (2,2), (2,4) are collinear. When (4,0) arrives, both (2,4) (a right turn) and (2,0) (collinear, on the bottom edge) are popped. The lower hull ends as (0,0), (4,0), (4,4), and the upper hull is built the same way from the other end. The duplicate (4,0) never appears in the trace: it was removed after sorting.

Every point is pushed once per pass and popped at most once, so the scan takes O(n) orientation tests. The sort takes O(n log n) comparisons, so the whole algorithm runs in O(n log n) time in the worst case. On input that is already sorted, the sort can be skipped and the rest is O(n).

Seeing the hull

The same function on random points in a disk. Hull vertices are drawn in blue, the hull in light blue, and all other points in dark gray. The hull starts at the leftmost point (the lowest of them if there is a tie) and runs counter-clockwise, which on the screen is the same as in maths coordinates, because the Svg helper flips the y-axis.

main.cppc++23
#include <algorithm>
#include <format>
#include <print>
#include <string>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

class Svg {
public:
    Svg(double width, double height) : w_(width), h_(height) {}
    void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
        add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
            a.x, fy(a.y), b.x, fy(b.y), color, width);
    }
    void dot(Point p, const char* color = "#1f2937", double r = 3) {
        add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
    }
    void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
                 const char* fill = "none") {
        std::string list;
        for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
        add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
    }
    void text(Point p, const std::string& s, double size = 8) {
        add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
    }
    void print() const {
        std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
        std::print("{}", body_);
        std::println("</svg>");
    }

private:
    double fy(double y) const { return h_ - y; }  // SVG's y axis points down
    template <class... Args>
    void add(std::format_string<Args...> fmt, Args&&... args) {
        body_ += std::format(fmt, std::forward<Args>(args)...);
        body_ += '\n';
    }
    double w_, h_;
    std::string body_;
};

int main() {
    Lcg rng{11};
    std::vector<Point> pts;
    for (int i = 0; i < 60; ++i) {  // random points in a disk, so the hull is round-ish
        Point p{rng.range(-90, 90), rng.range(-90, 90)};
        if (p.x * p.x + p.y * p.y <= 90 * 90) pts.push_back({p.x + 100, p.y + 100});
    }
    auto hull = convexHull(pts);
    Svg svg(200, 200);
    svg.polygon(hull, "#2563eb", "rgba(37,99,235,0.15)");
    for (Point p : pts) svg.dot(p, "#1f2937", 2);
    for (Point p : hull) svg.dot(p, "#2563eb", 3.5);
    svg.text({hull[0].x + 6, hull[0].y + 2}, "start", 8);
    svg.print();
}
Output

Press Run to compile and run this example.

Collinear points and duplicates: decide, then test

A hull function is only as good as its behavior on the inputs nobody thinks of. Here is the catalogue for this implementation, with both policies.

main.cppc++23
#include <algorithm>
#include <format>
#include <print>
#include <string>
#include <utility>
#include <vector>

struct Point { long long x, y; };

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

std::string str(const std::vector<Point>& h) {
    std::string s;
    for (Point p : h) s += std::format("({},{}) ", p.x, p.y);
    return s.empty() ? "(empty)" : s;
}

int main() {
    const std::pair<const char*, std::vector<Point>> cases[] = {
        {"no points", {}},
        {"one point", {{3, 3}}},
        {"one point, three times", {{3, 3}, {3, 3}, {3, 3}}},
        {"two points", {{5, 1}, {1, 2}}},
        {"collinear, diagonal", {{0, 0}, {3, 3}, {1, 1}, {2, 2}, {1, 1}}},
        {"collinear, vertical", {{2, 5}, {2, 1}, {2, 3}}},
        {"triangle plus copies", {{0, 0}, {4, 0}, {0, 3}, {4, 0}, {0, 0}}},
    };
    for (const auto& [name, pts] : cases) {
        std::println("{}:", name);
        std::println("  vertices only:     {}", str(convexHull(pts)));
        std::println("  keeping collinear: {}", str(convexHull(pts, true)));
    }
}
Output

Press Run to compile and run this example.

no points:
  vertices only:     (empty)
  keeping collinear: (empty)
one point:
  vertices only:     (3,3) 
  keeping collinear: (3,3) 
one point, three times:
  vertices only:     (3,3) 
  keeping collinear: (3,3) 
two points:
  vertices only:     (1,2) (5,1) 
  keeping collinear: (1,2) (5,1) 
collinear, diagonal:
  vertices only:     (0,0) (3,3) 
  keeping collinear: (0,0) (1,1) (2,2) (3,3) 
collinear, vertical:
  vertices only:     (2,1) (2,5) 
  keeping collinear: (2,1) (2,3) (2,5) 
triangle plus copies:
  vertices only:     (0,0) (4,0) (0,3) 
  keeping collinear: (0,0) (4,0) (0,3)

Every row is a decision, and each one is visible in the code.

  • Fewer than three distinct points: returned as they are (sorted). The hull of one point is that point; the hull of two is the segment between them.
  • All points on one line: the monotone chain would otherwise walk to the far end and back, listing the inner points twice when they are kept. The function detects this case up front and returns the two end points, or all points in order along the line.
  • Duplicates: removed after sorting, with std::ranges::unique. Without that, two equal points make an orientation of 0 that the keep-collinear policy would keep, and the same point would appear twice.

The vertical line in the catalogue matters. Lexicographic order sorts points with equal x by y, so a vertical edge is handled by the same code as any other: no special case for "infinite slope", because no slope is ever computed.

Testing the hull against its definition

The property test compares the monotone chain with an independent, slow definition. A pair of input points (p, q) forms a hull edge if no input point lies to the right of the line p→q, and every input point on that line lies between p and q. The corners are the end points of hull edges, and the boundary points are all points lying on hull edges. Testing every pair against every point takes O(n³) time, which is fine for a dozen points.

On top of matching that set, each result must start at the smallest point, turn left (or, with the keep policy, go straight) at every vertex, and have no input point outside any edge. The generator uses tiny coordinate ranges so that duplicates and collinear triples are common, and makes a third of the sets entirely collinear.

main.cppc++23
#include <algorithm>
#include <print>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

bool onSegment(Point a, Point b, Point p) {  // p collinear with a, b: between them?
    return std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x) &&
           std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}

// Brute force, O(n^3): (p, q) is a hull edge if no point is right of the line p->q
// and every point on that line lies between p and q. Returns the points on hull
// edges: all of them (keepCollinear) or only the edges' endpoints.
std::vector<Point> bruteForce(const std::vector<Point>& pts, bool keepCollinear) {
    std::vector<Point> result;
    for (Point p : pts)
        for (Point q : pts) {
            if (p == q) continue;
            bool edge = std::ranges::all_of(pts, [&](Point r) {
                int t = orientation(p, q, r);
                return t > 0 || (t == 0 && onSegment(p, q, r));
            });
            if (!edge) continue;
            for (Point r : pts)
                if ((r == p || r == q || keepCollinear) && orientation(p, q, r) == 0)
                    result.push_back(r);
        }
    if (pts.size() > 0 && result.empty()) result.push_back(pts[0]);  // all points equal
    std::ranges::sort(result, lexLess);
    auto dups = std::ranges::unique(result);
    result.erase(dups.begin(), dups.end());
    return result;
}

bool check(const std::vector<Point>& hull, const std::vector<Point>& pts, bool keep) {
    auto sorted = hull;
    std::ranges::sort(sorted, lexLess);
    if (sorted != bruteForce(pts, keep)) return false;      // the right points
    if (!hull.empty() && !(hull[0] == sorted[0])) return false;  // starts at the smallest
    const size_t n = hull.size();
    if (n < 3) return true;
    for (size_t i = 0; i < n; ++i) {
        Point a = hull[i], b = hull[(i + 1) % n], c = hull[(i + 2) % n];
        int turn = orientation(a, b, c);
        if (turn < 0 || (turn == 0 && !keep)) return false;  // convex, counter-clockwise
        for (Point p : pts)
            if (orientation(a, b, p) < 0) return false;      // nothing outside
    }
    return true;
}

int main() {
    int cases = 0, passed[2] = {};
    for (int seed = 0; seed < 3000; ++seed) {
        Lcg rng{static_cast<unsigned long long>(seed)};
        const int n = rng.range(0, 12), range = rng.range(0, 6), shape = rng.range(0, 2);
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i) {
            Point p{rng.range(0, range), rng.range(0, range)};
            if (shape == 1) p.y = 2 * p.x;  // every point on one line
            pts.push_back(p);
        }
        ++cases;
        for (bool keep : {false, true})
            passed[keep] += check(convexHull(pts, keep), pts, keep);
    }
    std::println("random point sets:            {}", cases);
    std::println("passed, vertices only:        {}", passed[0]);
    std::println("passed, keeping collinear:    {}", passed[1]);
}
Output

Press Run to compile and run this example.

random point sets:            3000
passed, vertices only:        3000
passed, keeping collinear:    3000

A test that always passes might not test anything, so it is worth breaking the code on purpose. Removing the duplicate-removal line makes the keep-collinear policy fail on 1,890 of the 3,000 sets and the corners-only policy on 540; changing the turn test to turn >= 0 for both policies makes 545 corners-only hulls fail. (Try it yourself, exercise 1, has you reproduce one of these.)

Graham scan: sorting by angle without computing angles

Ronald Graham's 1972 algorithm takes a different route. It picks a pivot that is certainly a hull vertex, the lowest point (the leftmost of the lowest, if there is a tie), sorts the other points by the angle they make around the pivot, and then runs the same stack scan as the monotone chain, once, around the whole hull.

The natural way to sort by angle is to compute atan2 for each point. Do not. The orientation test compares two angles exactly: b comes after a if orientation(pivot, a, b) > 0. This works as a sort comparator because of the choice of pivot. Every other point lies in the half-plane above it (or on the ray to its right), so all angles are in [0°, 180°), and on that range "counter-clockwise from" is a strict weak ordering. With an arbitrary pivot, angles would wrap around and the comparator would not be transitive, which std::sort does not tolerate.

Points in the same direction from the pivot compare as equal under orientation alone, so the comparator breaks ties by distance, nearer first, using exact squared distances.

main.cppc++23
#include <algorithm>
#include <print>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

long long dist2(Point a, Point b) {  // squared distance, exact
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

std::vector<Point> grahamScan(std::vector<Point> pts) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    // Pivot: the lowest point, leftmost among the lowest. Every other point is then
    // less than a half-turn counter-clockwise from the positive x direction.
    auto lower = [](Point a, Point b) { return a.y < b.y || (a.y == b.y && a.x < b.x); };
    std::iter_swap(pts.begin(), std::ranges::min_element(pts, lower));
    const Point o = pts[0];
    std::sort(pts.begin() + 1, pts.end(), [o](Point a, Point b) {
        int turn = orientation(o, a, b);
        if (turn != 0) return turn > 0;     // b is counter-clockwise from a
        return dist2(o, a) < dist2(o, b);   // same direction: nearer first
    });
    std::vector<Point> hull;
    for (Point p : pts) {
        while (hull.size() >= 2 &&
               orientation(hull[hull.size() - 2], hull.back(), p) <= 0)
            hull.pop_back();
        hull.push_back(p);
    }
    return hull;  // counter-clockwise from the pivot
}

// Rotate a hull so that it starts at its lexicographically smallest point.
std::vector<Point> normalized(std::vector<Point> h) {
    if (!h.empty()) std::ranges::rotate(h, std::ranges::min_element(h, lexLess));
    return h;
}

int main() {
    int same = 0, cases = 0;
    for (int seed = 0; seed < 5000; ++seed) {
        Lcg rng{static_cast<unsigned long long>(seed)};
        const int n = rng.range(0, 30), range = rng.range(0, 8);
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, range), rng.range(0, range)});
        ++cases;
        same += normalized(grahamScan(pts)) == convexHull(pts);
    }
    std::println("Graham scan equals monotone chain on {} of {} point sets", same, cases);
    std::vector<Point> example{{3, 0}, {6, 3}, {3, 3}, {0, 3}, {3, 6}, {1, 1}, {5, 5}};
    std::print("example, from the pivot:");
    for (Point p : grahamScan(example)) std::print(" ({},{})", p.x, p.y);
    std::println("");
}
Output

Press Run to compile and run this example.

Graham scan equals monotone chain on 5000 of 5000 point sets
example, from the pivot: (3,0) (6,3) (5,5) (3,6) (0,3) (1,1)

Graham scan and the monotone chain produce identical hulls on all 5,000 sets, after rotating the Graham result to start at the lexicographically smallest point instead of the pivot. Both are O(n log n) in the worst case, dominated by the sort.

The monotone chain is usually the one to choose. It needs no pivot, its sort is a plain coordinate comparison rather than a geometric comparator, and its tie handling is uniform. Graham scan's ties are subtle: nearer-first ordering is right for collinear points on the first edge from the pivot, but keeping collinear boundary points on the last edge back to the pivot needs farthest-first ordering there instead. The version above drops collinear points, which avoids the problem; the keep variant is a well-known source of bugs.

Jarvis march and QuickHull

Two more algorithms are worth knowing, because their running time depends on the input in instructive ways.

Jarvis march (R. A. Jarvis, 1973), also called gift wrapping, walks around the hull directly. Start at the leftmost point, which is a vertex. From the current vertex p, the next vertex is the point q such that every other point is on the left of p→q; finding it takes one pass over all n points. Repeat until back at the start. That is O(n) work per hull vertex, O(nh) in total: excellent when h is small, terrible when every point is on the hull.

QuickHull, published independently by W. F. Eddy in 1977 and A. Bykat in 1978, works like quicksort. The leftmost and rightmost points are vertices; the chord between them splits the other points into those below and those above. On each side, the point farthest from the chord is a vertex too (comparing distances needs no square root: the cross product cross(a, b, p) is proportional to the distance from p to the line). The triangle formed by the chord and that point contains no other hull vertices, so all points inside it are discarded at once, and the algorithm recurses on the points beyond each of the triangle's two new edges.

main.cppc++23
#include <algorithm>
#include <print>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc, exact
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right turn, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

long long dist2(Point a, Point b) {  // squared distance, exact
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

// Jarvis march (gift wrapping): from each hull vertex, find the point that has
// every other point on its left. O(n) per vertex, O(nh) in total.
std::vector<Point> jarvisMarch(std::vector<Point> pts) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    std::vector<Point> hull;
    Point p = pts[0];  // the leftmost point is on the hull
    do {
        hull.push_back(p);
        Point next = p == pts[0] ? pts[1] : pts[0];
        for (Point r : pts) {
            int turn = orientation(p, next, r);
            // r is clockwise from next, or in the same direction but farther away
            if (turn < 0 || (turn == 0 && dist2(p, r) > dist2(p, next))) next = r;
        }
        p = next;
    } while (!(p == hull[0]));
    return hull;
}

// QuickHull: the point farthest from a chord is a hull vertex; recurse on the two
// sides. Appends the vertices strictly right of a->b, in order from a to b.
void quickHullSide(const std::vector<Point>& pts, Point a, Point b,
                   std::vector<Point>& out) {
    if (pts.empty()) return;
    Point far = pts[0];
    for (Point p : pts) {  // most negative cross = farthest right; ties: lexicographic
        long long c = cross(a, b, p), best = cross(a, b, far);
        if (c < best || (c == best && lexLess(p, far))) far = p;
    }
    std::vector<Point> left, right;  // the parts beyond the edges a-far and far-b
    for (Point p : pts) {
        if (cross(a, far, p) < 0) left.push_back(p);
        else if (cross(far, b, p) < 0) right.push_back(p);
    }  // everything else is inside triangle (a, far, b) and is discarded
    quickHullSide(left, a, far, out);
    out.push_back(far);
    quickHullSide(right, far, b, out);
}

std::vector<Point> quickHull(std::vector<Point> pts) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    const Point a = pts.front(), b = pts.back();  // both are hull vertices
    std::vector<Point> below, above;
    for (Point p : pts) {
        long long c = cross(a, b, p);
        if (c < 0) below.push_back(p);
        else if (c > 0) above.push_back(p);
    }
    std::vector<Point> hull{a};
    quickHullSide(below, a, b, hull);
    hull.push_back(b);
    quickHullSide(above, b, a, hull);
    return hull;
}

int main() {
    int cases = 0, jarvisSame = 0, quickSame = 0;
    for (int seed = 0; seed < 5000; ++seed) {
        Lcg rng{static_cast<unsigned long long>(seed)};
        const int n = rng.range(0, 30), range = rng.range(0, 8);
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, range), rng.range(0, range)});
        auto expected = convexHull(pts);
        ++cases;
        jarvisSame += jarvisMarch(pts) == expected;
        quickSame += quickHull(pts) == expected;
    }
    std::println("point sets:                  {}", cases);
    std::println("Jarvis march equals chain:   {}", jarvisSame);
    std::println("QuickHull equals chain:      {}", quickSame);
}
Output

Press Run to compile and run this example.

point sets:                  5000
Jarvis march equals chain:   5000
QuickHull equals chain:      5000

Two details are easy to get wrong. When several points are equally far from a chord, they lie on a line parallel to it, and only the end points of that run are vertices; the program picks the lexicographically smallest, which is an end point. And points exactly on a chord are discarded, which implements the corners-only policy.

Previous
Segment Intersection and the Sweep Line
Next
Polygons: Area, Winding and Point-in-Polygon