AIOJ

2D Convex Hull

Definitions

Convex Polygon

A convex polygon is a simple polygon whose every interior angle lies in the range $[0,\pi]$ .

Convex Hull

The smallest convex polygon that contains all given points in the plane is called the convex hull.

Its definition is as follows: for a given set $X$ , the intersection of all convex sets containing $X$ , denoted by $S$ , is called the convex hull of $X$ .

Indeed, it can be understood as the shape formed by enclosing all given points with a rubber band.

The convex hull encloses all the given points with the minimum perimeter. If a concave polygon encloses all the points, its perimeter cannot be the minimum, as shown below. By the triangle inequality, a convex polygon is always optimal in terms of perimeter.

Andrew's Algorithm for Finding Convex Hulls

Common methods include the Graham 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 whose convex hull is to be found. The bottleneck is sorting all point coordinates by two keys.

Process

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

Clearly, after sorting, the minimum and maximum elements must both lie on the convex hull. Moreover, because it is a convex polygon, if we start from a point and walk counterclockwise, the path always turns 「left」; once it turns right, that segment is not on the convex hull. Therefore, we can use a monotone stack to maintain the upper and lower convex hulls.

Because the upper and lower convex hulls turn in different directions when viewed from left to right, to make the monotone stack work, we first enumerate in ascending order to find the lower convex hull, and then in descending order to find the upper convex hull.

When constructing a convex hull, suppose 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 element) are found to turn right, that is, their cross product is less than $0$ : $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0$ . Then pop the top element and return to the previous step to continue checking, until $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}\ge 0$ or only one element remains in the stack.

Usually, points lying on the edges of the convex hull do not need to be retained. Therefore, in the preceding paragraph, the condition $\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0$ containing the 「 $<$ 」 may, as appropriate, be changed to $\le$ , and the later condition should be changed to $>$ .

Andrew

Implementation

Code Implementation

C++

// stk[] is an integer array storing indices
// p[] stores vectors or points
tp = 0;                       // initialize the stack
std::sort(p + 1, p + 1 + n);  // sort the points
stk[++tp] = 1;
// Add the first element to the stack without updating used, so that 1 also updates the monotone stack when the convex hull is closed at the end
for (int i = 2; i <= n; ++i) {
  while (tp >= 2  // The * operator on the next line is overloaded as the cross product
         && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
    used[stk[tp--]] = 0;
  used[i] = 1;  // used indicates that the point is on the convex hull
  stk[++tp] = i;
}
int tmp = tp;  // tmp denotes the size of the lower convex hull
for (int i = n - 1; i > 0; --i)
  if (!used[i]) {
    // ↓ Finding the upper convex hull does not affect the lower convex 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 into a new array
  h[i] = p[stk[i]];
int ans = tp - 1;

Python

stk = []  # is an integer array storing indices
p = []  # stores vectors or points
tp = 0  # initialize the stack
p.sort()  # sort the points
tp = tp + 1
stk[tp] = 1
# Add the first element to the stack without updating used, so that 1 also updates the monotone stack when the convex hull is closed at the end
for i in range(2, n + 1):
    while tp >= 2 and (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0:
        # The * operator on the next line is overloaded as the cross product
        used[stk[tp]] = 0
        tp = tp - 1
    used[i] = 1  # used indicates that the point is on the convex hull
    tp = tp + 1
    stk[tp] = i
tmp = tp  # tmp denotes the size of the lower convex hull
for i in range(n - 1, 0, -1):
    if used[i] == False:
        #      ↓ Finding the upper convex hull does not affect the lower convex 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 code above, the final convex hull has $\textit{ans}$ elements (point $1$ is stored one extra time, so the $h$ array has $\textit{ans}+1$ elements), and they are sorted in counterclockwise order. Its perimeter is

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

Graham Scan

Properties

As with Andrew's algorithm, the Graham scan has time complexity $O(n\log n)$ , and its complexity bottleneck is also sorting all points.

Process

First find a point $P$ with the smallest y-coordinate among all the points. By the definition of a convex hull, we know that this point must lie on the convex hull. Then sort all points by their polar angle with respect to point P.

Similarly to Andrew's algorithm, consider starting from point $P$ and walking counterclockwise along the convex hull. Then all nodes we pass must turn 「left」. Formally, for any three consecutive points $P_1, P_2, P_3$ encountered in counterclockwise order on the convex hull, we must have $\overrightarrow{P_1 P_2} \times \overrightarrow{P_2 P_3} \ge 0$ .

Create a stack to store information about the convex hull. First push $P$ onto the stack, and then try to add each point in polar-angle order. If the point $P_0$ being pushed and the two points $P_1, P_2$ at the top of the stack (where $P_1$ is the top element) turn 「right」, pop the top element $P_1$ . Repeat this process until the point being pushed and the two points at the top of the stack satisfy the condition, or only one element remains in the stack, and 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 of point set $P$ and point set $Q$ , denoted by $P+Q$ , is defined as $P+Q=\{a+b|a\in P,b\in Q\}$ ; that is, regard each point in point set $Q$ as a vector, translate every point in point set $P$ along each of these vectors, and the set of results is point set $P+Q$ . Here we discuss only the Minkowski sum of convex hulls.

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

Translate $P$ along every vector of $Q$ :

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

Properties

If point sets $P$ and $Q$ are convex sets, then their Minkowski sum $P+Q$ is also a convex set.

??? note "Proof" Let $e,f\in P+Q$ , where $a,b \in P$ , $c,d\in Q$ , and $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}
$$

This proves the result.
```

If point sets $P$ and $Q$ are convex sets, then the edge set of their Minkowski sum $P+Q$ is obtained by connecting the edges of convex sets $P$ and $Q$ after sorting them by polar angle.

??? note "Proof" Without loss of generality, assume that the slope of every edge in convex set $P$ differs from the slope of every edge in $Q$ . Rotate the coordinate system so that one of the edges of $P$ , namely $XY$ , is parallel to the $x$ -axis and lies at the very bottom.

The lowest point of $Q$ is $U$, and the **lowest** and **leftmost** point of $P+Q$ is $A$.

We have $\vec{A} = \vec{X} + \vec{U}$, so $A$ must lie on the boundary of $P+Q$.

Similarly, in $P+Q$, the **lowest** and **rightmost** point is $B$, with $\vec{B} = \vec{Y} + \vec{U}$, and it must also lie on the boundary of $P+Q$.

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

If the rotation is carried out in order, the result successively forms every edge of $P+Q$.

This proves the result.

Implementation

Based on Property 2, we can sort convex sets $P,Q$ by polar angle to obtain their order of appearance on $P+Q$ , take $P_1+Q_1$ as the starting point of $P+Q$ , and then add the edges one by one using a method similar to merge.

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<T> &a, const Point<T> &b) { return (a ^ b) < 0; });
  for (usz i = 1; i < c.size(); ++i) c[i] = c[i] + c[i - 1];
  return c;
}

Example Problems

Example [JSOI2018] War

Given two convex hulls $P,Q$ , translate $q$ times the convex hull $Q$ and determine whether they have an intersection after each move. $1\le n,m\le 10^5,1\le q\le 10^5$ .

Implementation

#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<T> &a, const Point<T> &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

Basic Knowledge

Circle inversion: let the inversion center be $O$ and the inversion radius be $R$ . If a line passing 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 finding a convex hull is as follows:

First apply a small perturbation to the points to avoid the case of four coplanar points.
Given a known convex hull, add a new point $P$ . Treat $P$ as a point light source and cast rays toward the convex hull. The visible and invisible faces must be separated by several edges.
Delete the visible faces of the light and add the planes formed by their separating edges and $P$ . Repeating this process suffices. By Pick's theorem, Euler's formula (in a convex polyhedron, its vertices $V$ , number of edges $E$ , and number of faces $F$ satisfy $V−E+F=2$ ), and circle inversion, the complexity is $O(n^2)$ . ¹

Template Problem

P4724 Template 3D Convex Hull

Repeating the process above yields 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

References and Notes

Footnotes

A Short Note on Learning 3D Convex Hulls ↩