Mathematical Model

Consider n pins hammered into a wooden board. If we stretch a rubber band around all the pins and then release it, the band will snap tight and enclose all the pins. This geometric enclosure represents a convex hull.

A convex hull is defined as the convex polygon with the smallest perimeter that contains all points in a given set. Even if a concave polygon happens to have the same perimeter while containing all points, it still does not qualify as the convex hull. The key property is convexity, not minimal perimeter alone.

This article focuses on the Andrew algorithm, an efficient method for constructing convex hulls.

Prerequisites

Vector Cross Product

For two vectors v1(x1, y1) and v2(x2, y2), their cross product (also called the wedge product or outer product) produces a scalar value calculated as:

cross(v1, v2) = x1 × y2 - x2 × y1

The cross product encodes rotational relationships between vectors:

• Cross product > 0: Vector v1 lies clockwise from v2
• Cross product = 0: Vectors v1 and v2 are colinear
• Cross product < 0: Vector v1 lies counterclockwise from v2

Algorithm Overview

Convex Hull Construction

The algorithm proceeds as follows:

Step 1: Sorting

Sort all points by x-coordinate primarily, then by y-coordinate as a tiebreaker. The first point after sorting becomes the pole point (leftmost-bottommost point). This point is guaranteed to lie on the convex hull.

Step 2: Building the Lower Hull

Starting from the pole point, we traverse points counterclockwise, attempting to add each to the hull. Suppose the current hull contains m points. Let a and b be the last two points on the hull, with b being the most recent addition, and let c be the candidate point.

If vector bc is counterclockwise from vector ab (meaning cross(ab, bc) > 0), point c is valid and can be added.

Otherwise, if cross(ab, bc) ≤ 0, point b violates convexity and must be removed. We continue removing points untill the hull satisfies the convexity condition or contains only the pole point. Then c is added.

Step 3: Building the Upper Hull

After processing all points, we have construtced only the lower half of the hull. To complete the hull, we iterate backward from the last point, applying the same convexity check. This builds the upper hull.

If the lower hull contains k points, we must avoid popping points from position k onward during the upper hull construction to preserve correctness. The final hull consists of both halves connected seamlessly.

When the input contains collinear points and we choose to exclude intermediate collinear points, the condition becomes cross(ab, bc) < 0 for removal. If collinear points should remain on the hull boundary, we use cross(ab, bc) ≤ 0 for removal instead.

Edge Cases

When all input points are collinear, every point would be included twice in the worst case. Therefore, allocate the hull array with at least 2n capacity to prevent overflow. If more than one point exists in the final hull, exclude the duplicate pole point from the count.

The overall time complexity is O(n log n) due to the initial sorting step.

Reference Implementation (Retaining Collinear Points)

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int MAXP = 100000 + 5;

struct Point {
    double px, py;
    Point() : px(0), py(0) {}
    Point(double x, double y) : px(x), py(y) {}
    
    bool operator<(const Point& other) const {
        if (px != other.px) return px < other.px;
        return py < other.py;
    }
};

Point input[MAXP];
Point hull[MAXP * 2];
int pointCount;
int hullSize;

double wedge(Point a, Point b, Point c) {
    Point ab(b.px - a.px, b.py - a.py);
    Point ac(c.px - a.px, c.py - a.py);
    return (ab.px * ac.py) - (ac.px * ab.py);
}

double euclideanDist(Point a, Point b) {
    double dx = a.px - b.px;
    double dy = a.py - b.py;
    return sqrt(dx * dx + dy * dy);
}

int buildConvexHull() {
    sort(input, input + pointCount);
    hull[0] = input[0];
    hull[1] = input[1];
    hullSize = 1;
    
    for (int i = 1; i < pointCount; i++) {
        while (hullSize > 0 && wedge(hull[hullSize - 1], hull[hullSize], input[i]) < 0) {
            hullSize--;
        }
        hull[++hullSize] = input[i];
    }
    
    int lowerSize = hullSize;
    hull[++hullSize] = input[pointCount - 2];
    
    for (int i = pointCount - 2; i >= 0; i--) {
        while (hullSize > lowerSize && wedge(hull[hullSize - 1], hull[hullSize], input[i]) < 0) {
            hullSize--;
        }
        hull[++hullSize] = input[i];
    }
    
    return hullSize;
}

int main() {
    double perimeter = 0.0;
    scanf("%d", &pointCount);
    
    for (int i = 0; i < pointCount; i++) {
        scanf("%lf%lf", &input[i].px, &input[i].py);
    }
    
    if (pointCount >= 2) {
        int vertexCount = buildConvexHull();
        for (int i = 0; i < vertexCount; i++) {
            perimeter += euclideanDist(hull[i], hull[i + 1]);
        }
        perimeter += euclideanDist(hull[0], hull[vertexCount]);
        printf("%.2lf\n", perimeter);
    } else {
        printf("0\n");
    }
    
    return 0;
}

Reference Implementation (Removing Collinear Points)

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int MAXP = 100000 + 5;

struct Point {
    double px, py;
    Point() : px(0), py(0) {}
    Point(double x, double y) : px(x), py(y) {}
    
    bool operator<(const Point& other) const {
        if (px != other.px) return px < other.px;
        return py < other.py;
    }
};

Point input[MAXP];
Point stackArr[MAXP];
int n;
int stackTop;

double crossProduct(Point a, Point b) {
    return (a.px * b.py) - (b.px * a.py);
}

double areaSign(Point a, Point b, Point c) {
    Point ab(b.px - a.px, b.py - a.py);
    Point ac(c.px - a.px, c.py - a.py);
    return crossProduct(ab, ac);
}

double distance(Point a, Point b) {
    double dx = a.px - b.px;
    double dy = a.py - b.py;
    return sqrt(dx * dx + dy * dy);
}

int constructHull() {
    sort(input, input + n);
    stackTop = 0;
    
    for (int i = 0; i < n; i++) {
        while (stackTop > 1 && areaSign(stackArr[stackTop - 2], stackArr[stackTop - 1], input[i]) <= 0) {
            stackTop--;
        }
        stackArr[stackTop++] = input[i];
    }
    
    int lowerBound = stackTop;
    
    for (int i = n - 2; i >= 0; i--) {
        while (stackTop > lowerBound && areaSign(stackArr[stackTop - 2], stackArr[stackTop - 1], input[i]) <= 0) {
            stackTop--;
        }
        stackArr[stackTop++] = input[i];
    }
    
    if (n > 1) {
        stackTop--;
    }
    
    return stackTop;
}

int main() {
    double result = 0.0;
    scanf("%d", &n);
    
    for (int i = 0; i < n; i++) {
        scanf("%lf%lf", &input[i].px, &input[i].py);
    }
    
    if (n >= 2) {
        int m = constructHull();
        for (int i = 0; i < m - 1; i++) {
            result += distance(stackArr[i], stackArr[i + 1]);
        }
        result += distance(stackArr[0], stackArr[m - 1]);
        printf("%.2lf\n", result);
    } else {
        printf("0\n");
    }
    
    return 0;
}

Rotating Calipers: Finding the Diameter

Another classic problem involves finding the maximum distance between any two points in a set, known as the convex hull diameter. The rotating calipers technique solves this efficiently after the convex hull has been constructed.

Geometric Insight

Imagine two parallel lines squeezing the convex hull from opposite sides, rotating counterclockwise around the hull. At any orientation, one of two configurations occurs: either one line coincides with a hull edge while the other touches a hull vertex, or both lines touch hull vertices.

The second case is automatically covered when processing the first case, so we focus on configurations where a line coincides with hull edge ab and the opposite line touches vertex c. Here, c is called an antipodal point with respect to edge ab.

Key Properties

For a fixed edge, there is at most one antipodal point. The perpendicular distance from hull vertices to a given edge follows a unimodal distribution: it increases, reaches a maximum, then decreases.

The antipodal point maximizes the triangle area formed with the edge's endpoints. Since the area equals |cross(ab, bc)| / 2, the antipodal point corresponds to the maximum cross product magnitude for a fixed edge.

Algorithm

After constructing the hull with m vertices, initialize an antipodal pointer at vertex 1. For each edge i to i+1, advance the antipodal pointer while the cross product increases. The maximum distance encountered during this process yields the diameter.

The rotating calipers approach runs in O(n) time after hull construction, giving an overall complexity of O(n log n).

Reference Implementation

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int MAXP = 50000 + 5;

struct Point {
    double px, py;
    Point() : px(0), py(0) {}
    Point(double x, double y) : px(x), py(y) {}
    
    bool operator<(const Point& other) const {
        if (px != other.px) return px < other.px;
        return py < other.py;
    }
};

Point rawPoints[MAXP];
Point hullPoints[MAXP];
int totalPoints;
int hullCount;

double squaredDist(Point a, Point b) {
    double dx = a.px - b.px;
    double dy = a.py - b.py;
    return dx * dx + dy * dy;
}

double crossOp(Point a, Point b) {
    return (a.px * b.py) - (b.px * a.py);
}

double signedArea(Point a, Point b, Point c) {
    Point ab(b.px - a.px, b.py - a.py);
    Point ac(c.px - a.px, c.py - a.py);
    return crossOp(ab, ac);
}

int computeHull() {
    sort(rawPoints, rawPoints + totalPoints);
    hullCount = 0;
    
    for (int i = 0; i < totalPoints; i++) {
        while (hullCount > 1 && signedArea(hullPoints[hullCount - 2], hullPoints[hullCount - 1], rawPoints[i]) <= 0) {
            hullCount--;
        }
        hullPoints[hullCount++] = rawPoints[i];
    }
    
    int lowerSize = hullCount;
    
    for (int i = totalPoints - 2; i >= 0; i--) {
        while (hullCount > lowerSize && signedArea(hullPoints[hullCount - 2], hullPoints[hullCount - 1], rawPoints[i]) <= 0) {
            hullCount--;
        }
        hullPoints[hullCount++] = rawPoints[i];
    }
    
    if (totalPoints > 1) {
        hullCount--;
    }
    
    return hullCount;
}

double findDiameter() {
    int antipode = 1;
    double maxDist = 0.0;
    
    for (int i = 0; i < hullCount; i++) {
        int nextI = (i + 1) % hullCount;
        
        while (signedArea(hullPoints[i], hullPoints[nextI], hullPoints[antipode]) < 
               signedArea(hullPoints[i], hullPoints[nextI], hullPoints[(antipode + 1) % hullCount])) {
            antipode = (antipode + 1) % hullCount;
        }
        
        maxDist = max(maxDist, squaredDist(hullPoints[i], hullPoints[antipode]));
        maxDist = max(maxDist, squaredDist(hullPoints[nextI], hullPoints[antipode]));
    }
    
    return maxDist;
}

int main() {
    scanf("%d", &totalPoints);
    
    for (int i = 0; i < totalPoints; i++) {
        scanf("%lf%lf", &rawPoints[i].px, &rawPoints[i].py);
    }
    
    computeHull();
    printf("%.0lf\n", sqrt(findDiameter()));
    
    return 0;
}

Summary

The Andrew algorithm provides an elegant O(n log n) solution for convex hull construction by sorting points and building the hull in two passes. When combined with the rotating calipers technique, it efficiently solves the diameter problem. Both algorithms rely on the cross product's ability to encode directional information about vectors in the plane.