跳转至

凸包

二維凸包

定義

凸多邊形

凸多邊形是指所有內角大小都在 \([0,\pi]\) 範圍內的 簡單多邊形。

凸包

在平面上能包含所有給定點的最小凸多邊形叫做凸包。

其定義為：對於給定集合 \(X\)，所有包含 \(X\) 的凸集的交集 \(S\) 被稱為 \(X\) 的 凸包。

實際上可以理解為用一個橡皮筋包含住所有給定點的形態。

凸包用最小的周長圍住了給定的所有點。如果一個凹多邊形圍住了所有的點，它的周長一定不是最小，如下圖。根據三角不等式，凸多邊形在周長上一定是最優的。

Andrew 算法求凸包

常用的求法有 Graham 掃描法和 Andrew 算法，這裏主要介紹 Andrew 算法。

性質

該算法的時間複雜度為 \(O(n\log n)\)，其中 \(n\) 為待求凸包點集的大小，複雜度的瓶頸在於對所有點座標的雙關鍵字排序。

過程

首先把所有點以橫座標為第一關鍵字，縱座標為第二關鍵字排序。

顯然排序後最小的元素和最大的元素一定在凸包上。而且因為是凸多邊形，我們如果從一個點出發逆時針走，軌跡總是「左拐」的，一旦出現右拐，就説明這一段不在凸包上。因此我們可以用一個單調棧來維護上下凸殼。

因為從左向右看，上下凸殼所旋轉的方向不同，為了讓單調棧起作用，我們首先 升序枚舉 求出下凸殼，然後 降序 求出上凸殼。

求凸殼時，一旦發現即將進棧的點（\(P\)）和棧頂的兩個點（\(S_1,S_2\)，其中 \(S_1\) 為棧頂）行進的方向向右旋轉，即叉積小於 \(0\)：\(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0\)，則彈出棧頂，回到上一步，繼續檢測，直到 \(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}\ge 0\) 或者棧內僅剩一個元素為止。

通常情況下不需要保留位於凸包邊上的點，因此上面一段中 \(\overrightarrow{S_2S_1}\times \overrightarrow{S_1P}<0\) 這個條件中的「\(<\)」可以視情況改為 \(\le\)，同時後面一個條件應改為 \(>\)。

實現

代碼實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// stk[] 是整型，存的是下標
// p[] 存儲向量或點
tp = 0;                       // 初始化棧
std::sort(p + 1, p + 1 + n);  // 對點進行排序
stk[++tp] = 1;
// 棧內添加第一個元素，且不更新 used，使得 1 在最後封閉凸包時也對單調棧更新
for (int i = 2; i <= n; ++i) {
  while (tp >= 2  // 下一行 * 操作符被重載為叉積
         && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
    used[stk[tp--]] = 0;
  used[i] = 1;  // used 表示在凸殼上
  stk[++tp] = i;
}
int tmp = tp;  // tmp 表示下凸殼大小
for (int i = n - 1; i > 0; --i)
  if (!used[i]) {
    // ↓求上凸殼時不影響下凸殼
    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)  // 複製到新數組中去
  h[i] = p[stk[i]];
int ans = tp - 1;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
stk = [] # 是整型，存的是下標
p = [] # 存儲向量或點
tp = 0 # 初始化棧
p.sort() # 對點進行排序
tp = tp + 1
stk[tp] = 1
# 棧內添加第一個元素，且不更新 used，使得 1 在最後封閉凸包時也對單調棧更新
for i in range(2, n + 1):
    while tp >= 2 and (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0:
        # 下一行 * 操作符被重載為叉積
        used[stk[tp]] = 0
        tp = tp - 1
    used[i] = 1 # used 表示在凸殼上
    tp = tp + 1
    stk[tp] = i
tmp = tp # tmp 表示下凸殼大小
for i in range(n - 1, 0, -1):
    if used[i] == False:
        #      ↓求上凸殼時不影響下凸殼
        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

根據上面的代碼，最後凸包上有 \(\textit{ans}\) 個元素（額外存儲了 \(1\) 號點，因此 \(h\) 數組中有 \(\textit{ans}+1\) 個元素），並且按逆時針方向排序。周長就是

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

Graham 掃描法

性質

與 Andrew 算法相同，Graham 掃描法的時間複雜度為 \(O(n\log n)\)，複雜度瓶頸也在於對所有點排序。

過程

首先找到所有點中，縱座標最小的一個點 \(P\)。根據凸包的定義我們知道，這個點一定在凸包上。然後將所有的點以相對於點 P 的極角大小為關鍵字進行排序。

和 Andrew 算法類似地，我們考慮從點 \(P\) 出發，在凸包上逆時針走，那麼我們經過的所有節點一定都是「左拐」的。形式化地説，對於凸包逆時針方向上任意連續經過的三個點 \(P_1, P_2, P_3\)，一定滿足 \(\overrightarrow{P_1 P_2} \times \overrightarrow{P_2 P_3} \ge 0\)。

新建一個棧用於存儲凸包的信息，先將 \(P\) 壓入棧中，然後按照極角序依次嘗試加入每一個點。如果進棧的點 \(P_0\) 和棧頂的兩個點 \(P_1, P_2\)（其中 \(P_1\) 為棧頂）行進的方向「右拐」了，那麼就彈出棧頂的 \(P_1\)，不斷重複上述過程直至進棧的點與棧頂的兩個點滿足條件，或者棧中僅剩下一個元素，再將 \(P_0\) 壓入棧中。

代碼實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
struct Point {
  double x, y, ang;

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

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

三維凸包

基礎知識

圓的反演：反演中心為 \(O\)，反演半徑為 \(R\)，若經過 \(O\) 的直線經過 \(P\),\(P'\)，且 \(OP\times OP'=R^{2}\)，則稱 \(P\)、\(P'\) 關於 \(O\) 互為反演。

過程

求凸包的過程如下：

  • 首先對其微小擾動，避免出現四點共面的情況。
  • 對於一個已知凸包，新增一個點 \(P\)，將 \(P\) 視作一個點光源，向凸包做射線，可以知道，光線的可見面和不可見面一定是由若干條稜隔開的。
  • 將光的可見面刪去，並新增由其分割稜與 \(P\) 構成的平面。 重複此過程即可，由 Pick 定理、歐拉公式（在凸多面體中，其頂點 \(V\)、邊數 \(E\) 及面數 \(F\) 滿足 \(V−E+F=2\)）和圓的反演，複雜度 \(O(n^2)\)。1

模板題

P4724【模板】三維凸包

重複上述過程即可得到答案。

代碼實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>
using namespace std;
const int N = 2010;
const 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) { return {x - A.x, y - A.y, z - A.z}; }

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

  double operator&(Node A) { 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();

  printf("%.3f\n", ans);
  return 0;
}

練習

參考資料與註釋

  1. 三維凸包學習小記 ↩

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