Convex Hull

2D Convex Hull

Definition

Convex Polygon

A convex polygon is a simple polygon where all interior angles are in the range \([0,\pi]\).

Convex Hull

The smallest convex polygon that can contain all given points in a plane is called the convex hull.

Formally, for a given set \(X\), the intersection \(S\) of all convex sets containing \(X\) is called the convex hull of \(X\).

Intuitively, it can be thought of as the shape formed by a rubber band stretched around all given points.

The convex hull encloses all given points with the smallest possible perimeter. If a concave polygon encloses all points, its perimeter is certainly not minimal, as shown in the figure below. By the triangle inequality, a convex polygon is optimal in terms of perimeter.

Andrew's Algorithm for Convex Hull

Common methods for computing the convex hull include Graham's scan and Andrew's algorithm. Here we mainly introduce Andrew's algorithm.

Properties

The time complexity of this algorithm is \(O(n\log n)\), where \(n\) is the size of the point set for which the convex hull is to be computed. The bottleneck lies in sorting all points by two keys (coordinates).

Process

First, sort all points with the x-coordinate as the primary key and the y-coordinate as the secondary key.

Clearly, the smallest and largest elements after sorting must lie on the convex hull. Moreover, because it is a convex polygon, if we travel counterclockwise from a point, the trajectory always "turns left". Once a right turn occurs, it indicates that this segment is not part of the convex hull. Therefore, we can use a monotonic stack to maintain the upper and lower hulls.

Since the upper and lower hulls rotate in different directions when viewed from left to right, to make the monotonic stack work, we first enumerate in ascending order to compute the lower hull, and then enumerate in descending order to compute the upper hull.

When constructing the hull, whenever the point about to be pushed onto the stack (\(P\)) and the two points at the top of the stack (\(S_1,S_2\), where \(S_1\) is the top) form a right turn, i.e., the cross product is less than \(0\): \(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0\), pop the top of the stack, go back one step, and continue checking until \(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}\ge 0\) or only one element remains in the stack.

Typically, points lying on the edges of the convex hull do not need to be preserved. Therefore, the "\(<\)" in the condition \(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0\) above can be changed to \(\le\) as appropriate, and the subsequent condition should be changed to \(>\).

Implementation

Code Implementation

C++Python

// stk[] is integer array storing indices
// p[] stores vectors or points
tp = 0;                       // initialize stack
std::sort(p + 1, p + 1 + n);  // sort points
stk[++tp] = 1;
// add first element to stack, and do not update used, so that 1 also updates the monotonic stack when closing the hull
for (int i = 2; i <= n; ++i) {
  while (tp >= 2  // next line * operator overloaded as cross product
         && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
    used[stk[tp--]] = 0;
  used[i] = 1;  // used indicates on the hull
  stk[++tp] = i;
}
int tmp = tp;  // tmp indicates size of lower hull
for (int i = n - 1; i > 0; --i)
  if (!used[i]) {
    // ↓ computing upper hull does not affect lower hull
    while (tp > tmp && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
      used[stk[tp--]] = 0;
    used[i] = 1;
    stk[++tp] = i;
  }
for (int i = 1; i <= tp; ++i)  // copy to new array
  h[i] = p[stk[i]];
int ans = tp - 1;

stk = []  # integer list storing indices
p = []  # stores vectors or points
tp = 0  # initialize stack
p.sort()  # sort points
tp = tp + 1
stk[tp] = 1
# add first element to stack, and do not update used, so that 1 also updates the monotonic stack when closing the hull
for i in range(2, n + 1):
    while tp >= 2 and (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0:
        # next line * operator overloaded as cross product
        used[stk[tp]] = 0
        tp = tp - 1
    used[i] = 1  # used indicates on the hull
    tp = tp + 1
    stk[tp] = i
tmp = tp  # tmp indicates size of lower hull
for i in range(n - 1, 0, -1):
    if used[i] == False:
        #      ↓ computing upper hull does not affect lower hull
        while tp > tmp and (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0:
            used[stk[tp]] = 0
            tp = tp - 1
        used[i] = 1
        tp = tp + 1
        stk[tp] = i
for i in range(1, tp + 1):
    h[i] = p[stk[i]]
ans = tp - 1

According to the above code, there are \(\textit{ans}\) elements on the convex hull (with point \(1\) stored additionally, so the \(h\) array has \(\textit{ans}+1\) elements), sorted in counterclockwise order. The perimeter is

\[ \sum_{i=1}^{\textit{ans}}\left|\overrightarrow{h_ih_{i+1}}\right| \]

Graham's Scan

Properties

Same as Andrew's algorithm, the time complexity of Graham's scan is \(O(n\log n)\), with the bottleneck also being sorting all points.

Process

First, find the point \(P\) with the smallest y-coordinate among all points. By the definition of convex hull, this point must lie on the convex hull. Then sort all points by the polar angle relative to point \(P\).

Similar to Andrew's algorithm, consider starting from point \(P\) and traveling counterclockwise on the convex hull. All nodes we pass must be "left turns". Formally, for any three consecutive points \(P_1, P_2, P_3\) in counterclockwise order on the convex hull, it must satisfy \(\overrightarrow{P_1 P_2} \times \overrightarrow{P_2 P_3} \ge 0\).

Create a new stack to store convex hull information. First push \(P\) onto the stack, then attempt to add each point in order of polar angle. If the point \(P_0\) to be pushed and the two points \(P_1, P_2\) at the top of the stack (where \(P_1\) is the top) form a "right turn", then pop \(P_1\) from the stack. Repeat this process until the point to be pushed and the two points at the top satisfy the condition, or only one element remains in the stack, then push \(P_0\) onto the stack.

Code Implementation

struct Point {
  double x, y, ang;

  Point operator-(const Point& p) const { return {x - p.x, y - p.y, 0}; }
} p[MAXN];

double dis(Point p1, Point p2) {
  return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}

bool cmp(Point p1, Point p2) {
  if (p1.ang == p2.ang) {
    return dis(p1, p[1]) < dis(p2, p[1]);
  }
  return p1.ang < p2.ang;
}

double cross(Point p1, Point p2) { return p1.x * p2.y - p1.y * p2.x; }

int main() {
  for (int i = 2; i <= n; ++i) {
    if (p[i].y < p[1].y || (p[i].y == p[1].y && p[i].x < p[1].x)) {
      std::swap(p[1], p[i]);
    }
  }
  for (int i = 2; i <= n; ++i) {
    p[i].ang = atan2(p[i].y - p[1].y, p[i].x - p[1].x);
  }
  std::sort(p + 2, p + n + 1, cmp);
  sta[++top] = 1;
  for (int i = 2; i <= n; ++i) {
    while (top >= 2 &&
           cross(p[sta[top]] - p[sta[top - 1]], p[i] - p[sta[top]]) < 0) {
      top--;
    }
    sta[++top] = i;
  }
  return 0;
}

Minkowski Sum

Definition

The Minkowski sum \(P+Q\) of point sets \(P\) and \(Q\) is defined as \(P+Q=\{a+b|a\in P,b\in Q\}\). That is, treat each point in set \(Q\) as a vector, translate each point in set \(P\) along these vectors, and the resulting set is \(P+Q\). Here we only discuss the Minkowski sum of convex hulls.

For example: For point set \(P=\{(0,0),(-3,3),(2,1)\}\) and point set \(Q=\{(0,0),(-1,3),(1,4),(2,2)\}\),

Translate \(P\) along each vector of \(Q\):

It is not difficult to see that the new shape is also a convex hull:

Properties

1. If point sets \(P\) and \(Q\) are convex, then their Minkowski sum \(P+Q\) is also convex.

   Proof

   Let \(e,f\in P+Q\), with \(a,b \in P\), \(c,d\in Q\) such that \(e=a+c,f=b+d\). Then for any \(t\in[0,1]\):

   \[ \begin{aligned} te + (1-t)f &= t(a+c)+(1-t)(b+d)\\ &=(ta+(1-t)b)+(tc+(1-t)d)\\ &\in P+Q. \end{aligned} \]

   Q.E.D.

   1. If point sets \(P\) and \(Q\) are convex, then the edge set of their Minkowski sum \(P+Q\) is formed by connecting the edges of \(P\) and \(Q\) after sorting by polar angle.
   Proof

   Assume without loss of generality that the slope of any edge in convex set \(P\) is different from that of any edge in \(Q\). Rotate the coordinate system so that an edge \(XY\) of \(P\) is horizontal and at the bottom.

   Let \(U\) be the lowest point in \(Q\) at this moment, and \(A\) be the lowest and leftmost point in \(P+Q\).

   It follows that \(\vec{A} = \vec{X} + \vec{U}\), so \(A\) must lie on the boundary of \(P+Q\).

   Similarly, the lowest and rightmost point \(B\) in \(P+Q\) satisfies \(\vec{B} = \vec{Y} + \vec{U}\), and also lies on the boundary of \(P+Q\).

   Therefore, \(\vec{AB} = \vec{XY} + \vec{U}\).

   If we rotate sequentially, the results continuously form each edge of \(P+Q\).

   Q.E.D.

Implementation

Based on property 2, we can sort the convex sets \(P\) and \(Q\) by polar angle to obtain their order of appearance in \(P+Q\). Treat \(P_1+Q_1\) as the starting point of \(P+Q\), then use a merge-like approach to place edges sequentially.

Time complexity: \(O(n+m)\)

Implementation

template <class T>
struct Point {
  T x, y;

  Point(T x = 0, T y = 0) : x(x), y(y) {}

  friend Point operator+(const Point &a, const Point &b) {
    return {a.x + b.x, a.y + b.y};
  }

  friend Point operator-(const Point &a, const Point &b) {
    return {a.x - b.x, a.y - b.y};
  }

  // dot product
  friend T operator*(const Point &a, const Point &b) {
    return a.x * b.x + a.y * b.y;
  }

  // cross product
  friend T operator^(const Point &a, const Point &b) {
    return a.x * b.y - a.y * b.x;
  }
};

template <class T>
vector<Point<T>> minkowski_sum(vector<Point<T>> a, vector<Point<T>> b) {
  vector<Point<T>> c{a[0] + b[0]};
  for (usz i = 0; i + 1 < a.size(); ++i) a[i] = a[i + 1] - a[i];
  for (usz i = 0; i + 1 < b.size(); ++i) b[i] = b[i + 1] - b[i];
  a.pop_back(), b.pop_back();
  c.resize(a.size() + b.size() + 1);
  merge(a.begin(), a.end(), b.begin(), b.end(), c.begin() + 1,
        [](const Point<i64> &a, const Point<i64> &b) { return (a ^ b) < 0; });
  for (usz i = 1; i < c.size(); ++i) c[i] = c[i] + c[i - 1];
  return c;
}

Example

Example [JSOI2018] War

There are two convex hulls \(P\) and \(Q\). Translate \(Q\) \(q\) times, and determine if there is an intersection after each move. \(1\le n,m\le 10^5,1\le q\le 10^5\).

Implementation

```cpp
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>
using namespace std;
using i64 = int64_t;
using isz = ptrdiff_t;
using usz = size_t;

template <class T>
struct Point {
  T x, y;

  Point(T x = 0, T y = 0) : x(x), y(y) {}

  friend Point operator+(const Point &a, const Point &b) {
    return {a.x + b.x, a.y + b.y};
  }

  friend Point operator-(const Point &a, const Point &b) {
    return {a.x - b.x, a.y - b.y};
  }

  // Dot product
  friend T operator*(const Point &a, const Point &b) {
    return a.x * b.x + a.y * b.y;
  }

  // Cross product
  friend T operator^(const Point &a, const Point &b) {
    return a.x * b.y - a.y * b.x;
  }

  friend istream &operator>>(istream &is, Point &p) { return is >> p.x >> p.y; }
};

template <class T>
vector<Point<T>> convex_hull(vector<Point<T>> p) {
  assert(!p.empty());
  sort(p.begin(), p.end(),
       [](const Point<i64> &a, const Point<i64> &b) { return a.x < b.x; });
  vector<Point<T>> u{p[0]}, d{p.back()};
  for (usz i = 1; i < p.size(); ++i) {
    while (u.size() >= 2 &&
           ((u.back() - u[u.size() - 2]) ^ (p[i] - u.back())) > 0)
      u.pop_back();
    u.push_back(p[i]);
  }
  for (usz i = p.size() - 2; (isz)i >= 0; --i) {
    while (d.size() >= 2 &&
           ((d.back() - d[d.size() - 2]) ^ (p[i] - d.back())) > 0)
      d.pop_back();
    d.push_back(p[i]);
  }
  u.insert(u.end(), d.begin() + 1, d.end());
  return u;
}

template <class T>
vector<Point<T>> minkowski_sum(vector<Point<T>> a, vector<Point<T>> b) {
  vector<Point<T>> c{a[0] + b[0]};
  for (usz i = 0; i + 1 < a.size(); ++i) a[i] = a[i + 1] - a[i];
  for (usz i = 0; i + 1 < b.size(); ++i) b[i] = b[i + 1] - b[i];
  a.pop_back(), b.pop_back();
  c.resize(a.size() + b.size() + 1);
  merge(a.begin(), a.end(), b.begin(), b.end(), c.begin() + 1,
        [](const Point<i64> &a, const Point<i64> &b) { return (a ^ b) < 0; });
  for (usz i = 1; i < c.size(); ++i) c[i] = c[i] + c[i - 1];
  return c;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  uint32_t n, m, q;
  vector<Point<i64>> a, b;
  cin >> n >> m >> q;
  a.resize(n), b.resize(m);
  for (auto &p : a) cin >> p;
  for (auto &p : b) cin >> p, p = 0 - p;
  a = convex_hull(a), b = convex_hull(b);
  a = minkowski_sum(a, b);
  a.pop_back();
  for (usz i = 1; i < a.size(); ++i) a[i] = a[i] - a[0];
  while (q--) {
    Point<i64> v;
    cin >> v;
    v = v - a[0];
    if (v.x < 0) {
      cout << "0\n";
      continue;
    }
    auto it = upper_bound(
        a.begin() + 1, a.end(), v,
        [](const Point<i64> &a, const Point<i64> &b) { return (a ^ b) < 0; });
    if (it == a.begin() + 1 || it == a.end()) {
      cout << "0\n";
      continue;
    }
    i64 s0 = *it ^ *prev(it), s1 = v ^ *prev(it), s2 = *it ^ v;
    cout << (s1 >= 0 && s2 >= 0 && s1 + s2 <= s0) << '\n';
  }
  return 0;
}

```

3D Convex Hull

Basics

Circle inversion: Inversion center is \(O\), inversion radius is \(R\). If a line through \(O\) passes through \(P\) and \(P'\), and \(OP\times OP'=R^{2}\), then \(P\) and \(P'\) are said to be inverses of each other with respect to \(O\).

Process

The process of computing the convex hull is as follows:

• First, apply a slight perturbation to avoid four points being coplanar.
• For a known convex hull, add a new point \(P\). Treat \(P\) as a point light source and cast rays toward the convex hull. It can be seen that the visible and invisible faces of the light must be separated by several edges.
• Delete the visible faces of the light and add the planes formed by the separating edges and \(P\). Repeat this process. By Pick's theorem, Euler's formula (for convex polyhedra, the number of vertices \(V\), edges \(E\), and faces \(F\) satisfy \(V−E+F=2\)), and circle inversion, the complexity is \(O(n^2)\).<sup>1</sup>

Template Problem

P4724【Template】3D Convex Hull

Repeat the above process to obtain the answer.

Code Implementation

#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
using namespace std;
constexpr int N = 2010;
constexpr double eps = 1e-9;
int n, cnt, vis[N][N];
double ans;

double Rand() { return rand() / (double)RAND_MAX; }

double reps() { return (Rand() - 0.5) * eps; }

struct Node {
  double x, y, z;

  void shake() {
    x += reps();
    y += reps();
    z += reps();
  }

  double len() { return sqrt(x * x + y * y + z * z); }

  Node operator-(Node A) const { return {x - A.x, y - A.y, z - A.z}; }

  Node operator*(Node A) const {
    return {y * A.z - z * A.y, z * A.x - x * A.z, x * A.y - y * A.x};
  }

  double operator&(Node A) const { return x * A.x + y * A.y + z * A.z; }
} A[N];

struct Face {
  int v[3];

  Node Normal() { return (A[v[1]] - A[v[0]]) * (A[v[2]] - A[v[0]]); }

  double area() { return Normal().len() / 2.0; }
} f[N], C[N];

int see(Face a, Node b) { return ((b - A[a.v[0]]) & a.Normal()) > 0; }

void Convex_3D() {
  f[++cnt] = {1, 2, 3};
  f[++cnt] = {3, 2, 1};

  for (int i = 4, cc = 0; i <= n; i++) {
    for (int j = 1, v; j <= cnt; j++) {
      if (!(v = see(f[j], A[i]))) C[++cc] = f[j];

      for (int k = 0; k < 3; k++) vis[f[j].v[k]][f[j].v[(k + 1) % 3]] = v;
    }

    for (int j = 1; j <= cnt; j++)
      for (int k = 0; k < 3; k++) {
        int x = f[j].v[k], y = f[j].v[(k + 1) % 3];

        if (vis[x][y] && !vis[y][x]) C[++cc] = {x, y, i};
      }

    for (int j = 1; j <= cc; j++) f[j] = C[j];

    cnt = cc;
    cc = 0;
  }
}

int main() {
  cin >> n;

  for (int i = 1; i <= n; i++) cin >> A[i].x >> A[i].y >> A[i].z, A[i].shake();

  Convex_3D();

  for (int i = 1; i <= cnt; i++) ans += f[i].area();

  cout << fixed << setprecision(3) << ans << '\n';
  return 0;
}

Exercises

• UVa11626 Convex Hull

• 「USACO5.1」Fencing the Cows

• POJ1873 The Fortified Forest

• POJ1113 Wall

• USACO22JAN Multiple Choice Test P

• 「SHOI2012」Credit Card Convex Hull

References and Notes

---

1. Notes on Learning 3D Convex Hull ↩

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用