三角剖分

在幾何中，三角剖分是指將平面對象細分為三角形，並且通過擴展將高維幾何對象細分為單純形。 對於一個給定的點集，有很多種三角剖分，如：

OI 中的三角剖分主要指二維幾何中的完美三角剖分（二維 Delaunay 三角剖分，簡稱 DT）。

Delaunay 三角剖分

定義

在數學和計算幾何中，對於給定的平面中的離散點集 \(P\)，其 Delaunay 三角剖分 DT(\(P\)) 滿足：

1. 空圓性：DT(\(P\)) 是 唯一 的（任意四點不能共圓），在 DT(\(P\)) 中，任意 三角形的外接圓範圍內不會有其它點存在。
2. 最大化最小角：在點集 \(P\) 可能形成的三角剖分中，DT(\(P\)) 所形成的三角形的最小角最大。從這個意義上講，DT(\(P\)) 是 最接近於規則化 的三角剖分。具體的説是在兩個相鄰的三角形構成凸四邊形的對角線，在相互交換後，兩個內角的最小角不再增大。

性質

1. 最接近：以最接近的三點形成三角形，且各線段（三角形的邊）皆不相交。
2. 唯一性：不論從區域何處開始構建，最終都將得到一致的結果（點集中任意四點不能共圓）。
3. 最優性：任意兩個相鄰三角形構成的凸四邊形的對角線如果可以互換的話，那麼兩個三角形六個內角中最小角度不會變化。
4. 最規則：如果將三角剖分中的每個三角形的最小角進行升序排列，則 Delaunay 三角剖分的排列得到的數值最大。
5. 區域性：新增、刪除、移動某一個頂點只會影響鄰近的三角形。
6. 具有凸邊形的外殼：三角剖分最外層的邊界形成一個凸多邊形的外殼。

構造 DT 的分治算法

DT 有很多種構造算法，在 \(O(n \log n)\) 的構造算法中，分治算法是最易於理解和實現的。

分治構造 DT 的第一步是將給定點集按照 \(x\) 座標 升序 排列，如下圖是排好序的大小為 \(10\) 的點集。

一旦點集有序，我們就可以不斷地將其分成兩個部分（分治），直到子點集大小不超過 \(3\)。然後這些子點集可以立刻剖分為一個三角形或線段。

然後在分治回溯的過程中，已經剖分好的左右子點集可以依次合併。合併後的剖分包含 LL-edge（左側子點集的邊）。RR-edge（右側子點集的邊），LR-edge（連接左右剖分產生的新的邊），如圖 LL-edge（灰色），RR-edge（紅色），LR-edge（藍色）。對於合併後的剖分，為了維持 DT 性質，我們 可能 需要刪除部分 LL-edge 和 RR-edge，但我們在合併時 不會 增加 LL-edge 和 RR-edge。

合併左右兩個剖分的第一步是插入 base LR-edge，base LR-edge 是 最底部 的不與 任何 LL-edge 及 RR-edge 相交的 LR-edge。

然後，我們需要確定下一條 緊接在 base LR-edge 之上的 LR-edge。比如對於右側點集，下一條 LR-edge 的可能端點（右端點）為與 base LR-edge 右端點相連的 RR-edge 的另一端點（\(6, 7, 9\) 號點），左端點即為 \(2\) 號點。

對於可能的端點，我們需要按以下兩個標準檢驗：

1. 其對應 RR-edge 與 base LR-edge 的夾角小於 \(180\) 度。
2. base LR-edge 兩端點和這個可能點三點構成的圓內不包含任何其它 可能點。

如上圖，\(6\) 號可能點所對應的綠色圓包含了 \(9\) 號可能點，而 \(7\) 號可能點對應的紫色圓則不包含任何其它可能點，故 \(7\) 號點為下一條 LR-edge 的右端點。

對於左側點集，我們做鏡像處理即可。

當左右點集都不再含有符合標準的可能點時，合併即完成。當一個可能點符合標準，一條 LR-edge 就需要被添加，對於與需要添加的 LR-edge 相交的 LL-edge 和 RR-edge，將其刪除。

當左右點集均存在可能點時，判斷左邊點所對應圓是否包含右邊點，若包含則不符合；對於右邊點也是同樣的判斷。一般只有一個可能點符合標準（除非四點共圓）。

當這條 LR-edge 添加好後，將其作為 base LR-edge 重複以上步驟，繼續添加下一條，直到合併完成。

代碼

實現

#include <algorithm>
#include <cmath>
#include <cstring>
#include <list>
#include <utility>
#include <vector>

const double EPS = 1e-8;
const int MAXV = 10000;

struct Point {
  double x, y;
  int id;

  Point(double a = 0, double b = 0, int c = -1) : x(a), y(b), id(c) {}

  bool operator<(const Point &a) const {
    return x < a.x || (fabs(x - a.x) < EPS && y < a.y);
  }

  bool operator==(const Point &a) const {
    return fabs(x - a.x) < EPS && fabs(y - a.y) < EPS;
  }

  double dist2(const Point &b) {
    return (x - b.x) * (x - b.x) + (y - b.y) * (y - b.y);
  }
};

struct Point3D {
  double x, y, z;

  Point3D(double a = 0, double b = 0, double c = 0) : x(a), y(b), z(c) {}

  Point3D(const Point &p) { x = p.x, y = p.y, z = p.x * p.x + p.y * p.y; }

  Point3D operator-(const Point3D &a) const {
    return Point3D(x - a.x, y - a.y, z - a.z);
  }

  double dot(const Point3D &a) { return x * a.x + y * a.y + z * a.z; }
};

struct Edge {
  int id;
  std::list<Edge>::iterator c;

  Edge(int id = 0) { this->id = id; }
};

int cmp(double v) { return fabs(v) > EPS ? (v > 0 ? 1 : -1) : 0; }

double cross(const Point &o, const Point &a, const Point &b) {
  return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

Point3D cross(const Point3D &a, const Point3D &b) {
  return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
                 a.x * b.y - a.y * b.x);
}

int inCircle(const Point &a, Point b, Point c, const Point &p) {
  if (cross(a, b, c) < 0) std::swap(b, c);
  Point3D a3(a), b3(b), c3(c), p3(p);
  b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
  Point3D f = cross(b3, c3);
  return cmp(p3.dot(f));  // check same direction, in: < 0, on: = 0, out: > 0
}

int intersection(const Point &a, const Point &b, const Point &c,
                 const Point &d) {  // seg(a, b) and seg(c, d)
  return cmp(cross(a, c, b)) * cmp(cross(a, b, d)) > 0 &&
         cmp(cross(c, a, d)) * cmp(cross(c, d, b)) > 0;
}

class Delaunay {
 public:
  std::list<Edge> head[MAXV];  // graph
  Point p[MAXV];
  int n, rename[MAXV];

  void init(int n, Point p[]) {
    memcpy(this->p, p, sizeof(Point) * n);
    std::sort(this->p, this->p + n);
    for (int i = 0; i < n; i++) rename[p[i].id] = i;
    this->n = n;
    divide(0, n - 1);
  }

  void addEdge(int u, int v) {
    head[u].push_front(Edge(v));
    head[v].push_front(Edge(u));
    head[u].begin()->c = head[v].begin();
    head[v].begin()->c = head[u].begin();
  }

  void divide(int l, int r) {
    if (r - l <= 2) {  // #point <= 3
      for (int i = l; i <= r; i++)
        for (int j = i + 1; j <= r; j++) addEdge(i, j);
      return;
    }
    int mid = (l + r) / 2;
    divide(l, mid);
    divide(mid + 1, r);

    std::list<Edge>::iterator it;
    int nowl = l, nowr = r;

    for (int update = 1; update;) {
      // find left and right convex, lower common tangent
      update = 0;
      Point ptL = p[nowl], ptR = p[nowr];
      for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
        Point t = p[it->id];
        double v = cross(ptR, ptL, t);
        if (cmp(v) > 0 || (cmp(v) == 0 && ptR.dist2(t) < ptR.dist2(ptL))) {
          nowl = it->id, update = 1;
          break;
        }
      }
      if (update) continue;
      for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
        Point t = p[it->id];
        double v = cross(ptL, ptR, t);
        if (cmp(v) < 0 || (cmp(v) == 0 && ptL.dist2(t) < ptL.dist2(ptR))) {
          nowr = it->id, update = 1;
          break;
        }
      }
    }

    addEdge(nowl, nowr);  // add tangent

    for (int update = 1; true;) {
      update = 0;
      Point ptL = p[nowl], ptR = p[nowr];
      int ch = -1, side = 0;
      for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
        if (cmp(cross(ptL, ptR, p[it->id])) > 0 &&
            (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
          ch = it->id, side = -1;
        }
      }
      for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
        if (cmp(cross(ptR, p[it->id], ptL)) > 0 &&
            (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
          ch = it->id, side = 1;
        }
      }
      if (ch == -1) break;  // upper common tangent
      if (side == -1) {
        for (it = head[nowl].begin(); it != head[nowl].end();) {
          if (intersection(ptL, p[it->id], ptR, p[ch])) {
            head[it->id].erase(it->c);
            head[nowl].erase(it++);
          } else {
            it++;
          }
        }
        nowl = ch;
        addEdge(nowl, nowr);
      } else {
        for (it = head[nowr].begin(); it != head[nowr].end();) {
          if (intersection(ptR, p[it->id], ptL, p[ch])) {
            head[it->id].erase(it->c);
            head[nowr].erase(it++);
          } else {
            it++;
          }
        }
        nowr = ch;
        addEdge(nowl, nowr);
      }
    }
  }

  std::vector<std::pair<int, int> > getEdge() {
    std::vector<std::pair<int, int> > ret;
    ret.reserve(n);
    std::list<Edge>::iterator it;
    for (int i = 0; i < n; i++) {
      for (it = head[i].begin(); it != head[i].end(); it++) {
        if (it->id < i) continue;
        ret.push_back(std::make_pair(p[i].id, p[it->id].id));
      }
    }
    return ret;
  }
};

Voronoi 圖

Voronoi 圖由一組由連接兩鄰點直線的垂直平分線組成的連續多邊形組成，根據 \(n\) 個在平面上不重合種子點，把平面分成 \(n\) 個區域，使得每個區域內的點到它所在區域的種子點的距離比到其它區域種子點的距離近。

Voronoi 圖是 Delaunay 三角剖分的對偶圖，可以使用構造 Delaunay 三角剖分的分治算法求出三角網，再使用最左轉線算法求出其對偶圖實現在 \(O(n \log n)\) 的時間複雜度下構造 Voronoi 圖。

題目

SGU 383 Caravans 三角剖分 + 倍增

ContestHunter. 無盡的毀滅 三角剖分求對偶圖建 Voronoi 圖

Codeforces Gym 103485M. Constellation collection 三角剖分之後建圖進行 Floodfill

參考資料與拓展閲讀

1. Wikipedia - Triangulation (geometry)
2. Wikipedia - Delaunay triangulation
3. Samuel Peterson -Computing Constrained Delaunay Triangulations in 2-D (1997-98)

---

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