Delaunay 三角剖分

终于写完了，贴个代码先。题目是 公路修建。

只写了 Bowyer-Watson 和 Guibas-Stolfi。

Bowyer-Watson 在网上看了一圈找不到有 OI 选手写，看着一堆工程码风不知道是什么意思，只好自己琢磨，写了好久。写的时候还用到了在学 Guibas-Stolfi 时学到的，但没在 Guibas-Stolfi 里用到的 “三角形链表” 技巧。

花了整整两天是写出来了，代码可读性几乎为零，而且常数特大。但网上说复杂度可以做到 (mathcal{O}(n sqrt n))，我就不知道怎么做。目前瓶颈在于：对于一张三角剖分好的图，求一个点在哪一个三角形内，要求复杂度不超过 (mathcal{O}(sqrt n))。这玩意我只会 (mathcal{O}(n)) 求。求各位大佬指教。

Bowyer-Watson：

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <utility>
#include <vector>

const int MaxN = 5005;
const int MaxF = 200000;
const double eps = 1e-9;

#ifdef Tweetuzki
#define debug(arg...) fprintf(stderr, arg)
#else
#define debug(arg...) void(0)
#endif

typedef struct vec_t {
  double x, y;
  vec_t(double _x = 0, double _y = 0) { x = _x, y = _y; }
  inline friend vec_t operator+(const vec_t &a, const vec_t &b) { return vec_t(a.x + b.x, a.y + b.y); }
  inline friend vec_t operator-(const vec_t &a, const vec_t &b) { return vec_t(a.x - b.x, a.y - b.y); }
  inline friend vec_t operator*(const vec_t &a, double k) { return vec_t(a.x * k, a.y * k); }
  inline friend double dot(const vec_t &a, const vec_t &b) { return a.x * b.x + a.y * b.y; }
  inline friend double cross(const vec_t &a, const vec_t &b) { return a.x * b.y - a.y * b.x; }
  inline friend double mod(const vec_t &a) { return sqrt(a.x * a.x + a.y * a.y); }

  inline void shake() {
    if (rand() % 2 == 0) x += 1.0 * rand() / RAND_MAX * eps;
    else x -= 1.0 * rand() / RAND_MAX * eps;
    if (rand() % 2 == 0) y += 1.0 * rand() / RAND_MAX * eps;
    else y -= 1.0 * rand() / RAND_MAX * eps;
  }
} node_t;

typedef struct vec3_t {
  double x, y, z;
  vec3_t(double _x = 0, double _y = 0, double _z = 0) { x = _x, y = _y, z = _z; }
  inline friend vec3_t operator+(const vec3_t &a, const vec3_t &b) { return vec3_t(a.x + b.x, a.y + b.y, a.z + b.z); }
  inline friend vec3_t operator-(const vec3_t &a, const vec3_t &b) { return vec3_t(a.x - b.x, a.y - b.y, a.z - b.z); }
  inline friend vec3_t operator*(const vec3_t &a, double k) { return vec3_t(a.x * k, a.y * k, a.z * k); }
  inline friend vec3_t cross(const vec3_t &a, const vec3_t &b) { return vec3_t(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); }
  inline friend double dot(const vec3_t &a, const vec3_t &b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
} node3_t;

struct triangle_t {
  int a, b, c;
  int la, lb, lc;

  triangle_t(int _a = 0, int _b = 0, int _c = 0, int _la = 0, int _lb = 0, int _lc = 0) {
    a = _a, b = _b, c = _c;
    la = _la, lb = _lb, lc = _lc;
  }
};

struct edge_t {
  int u, v;
  double w;
  edge_t(int _u = 0, int _v = 0, double _w = 0) { u = _u, v = _v, w = _w; }
  inline friend bool operator<(const edge_t &a, const edge_t &b) { return a.w < b.w; }
};

int N, CntF;
node_t A[MaxN + 5];
triangle_t F[MaxF + 5];
bool Del[MaxF + 5];

inline vec3_t mapping(const vec_t &a) { return vec3_t(a.x, a.y, a.x * a.x + a.y * a.y); }

inline bool inCircumcircle(const node_t &a, const node_t &b, const node_t &c, const node_t &p) {
  node3_t _a = mapping(a), _b = mapping(b), _c = mapping(c), _p = mapping(p);
  if (cross(b - a, c - a) < 0) std::swap(_b, _c);
  node3_t normal = cross(_b - _a, _c - _a);
  if (dot(normal, _p - _a) > eps) return false;
  else return true;
}

inline bool inCircumcircle(const triangle_t &t, const node_t &p) { return inCircumcircle(A[t.a], A[t.b], A[t.c], p); }

void init() {
  scanf("%d", &N);
  for (int i = 1; i <= N; ++i) {
    scanf("%lf %lf", &A[i].x, &A[i].y);
    A[i].shake();
  }
  A[N + 1] = node_t(-1e6, -1e6), A[N + 2] = node_t(-1e6, 1e6);
  A[N + 3] = node_t(1e6, -1e6), A[N + 4] = node_t(1e6, 1e6);
  F[++CntF] = triangle_t(N + 1, N + 3, N + 2, 2, 0, 0);
  F[++CntF] = triangle_t(N + 4, N + 2, N + 3, 1, 0, 0);
}

std::vector< std::pair<int, int> > Lk[MaxN + 5];
int To[MaxN + 5], Tof[MaxN + 5];
int NodeStk[MaxN + 5], NodeTp;

void dfs(int u, int f, int beg) {
  debug("dfs(%d, %d, %d)
", u, f, beg);
  if (u == beg) {
    if (f != 0) return;
    NodeStk[NodeTp++] = u;
    auto p = *(Lk[u].begin());
    To[u] = p.first;
    Tof[u] = p.second;
    dfs(p.first, u, beg);
  } else {
    NodeStk[NodeTp++] = u;
    for (auto p : Lk[u]) {
      if (p.first == f) continue;
      To[u] = p.first;
      Tof[u] = p.second;
      dfs(p.first, u, beg);
    }
  }
}

inline bool cmp(int a, int b, int c, int d) {
  if (a == c && b == d) return true;
  if (a == d && b == c) return true;
  return false;
}

inline void insert(int insertNode) {
  static int stk[MaxF + 5];
  int tp = 0;
  for (int i = 1; i <= CntF; ++i)
    if (Del[i] == false && inCircumcircle(F[i], A[insertNode]) == true) {
      stk[++tp] = i;
      Del[i] = true;
    }
  for (int i = 1; i <= tp; ++i) debug("stk[%d] = %d
", i, stk[i]);
  static int e[MaxF + 5][3]; int cnte = 0;
  for (int i = 1; i <= tp; ++i) {
    int x = stk[i];
    if (F[x].la == 0) {
      cnte++;
      e[cnte][0] = F[x].b, e[cnte][1] = F[x].c, e[cnte][2] = 0;
    } else if (inCircumcircle(F[F[x].la], A[insertNode]) == false) {
      int y = F[x].la;
      cnte++;
      if (F[y].la == x) {
        F[y].la = -1;
        e[cnte][0] = F[y].b, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else if (F[y].lb == x) {
        F[y].lb = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else {
        F[y].lc = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].b, e[cnte][2] = y;
      }
    }
    if (F[x].lb == 0) {
      cnte++;
      e[cnte][0] = F[x].a, e[cnte][1] = F[x].c, e[cnte][2] = 0;
    } else if (inCircumcircle(F[F[x].lb], A[insertNode]) == false) {
      int y = F[x].lb;
      cnte++;
      if (F[y].la == x) {
        F[y].la = -1;
        e[cnte][0] = F[y].b, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else if (F[y].lb == x) {
        F[y].lb = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else {
        F[y].lc = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].b, e[cnte][2] = y;
      }
    }
    if (F[x].lc == 0) {
      cnte++;
      e[cnte][0] = F[x].a, e[cnte][1] = F[x].b, e[cnte][2] = 0;
    } else if (inCircumcircle(F[F[x].lc], A[insertNode]) == false) {
      int y = F[x].lc;
      cnte++;
      if (F[y].la == x) {
        F[y].la = -1;
        e[cnte][0] = F[y].b, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else if (F[y].lb == x) {
        F[y].lb = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].c, e[cnte][2] = y;
      } else {
        F[y].lc = -1;
        e[cnte][0] = F[y].a, e[cnte][1] = F[y].b, e[cnte][2] = y;
      }
    }
  }
  for (int i = 1; i <= cnte; ++i)
    debug("e[%d] = (%d, %d, %d)
", i, e[i][0], e[i][1], e[i][2]);
  for (int i = 1; i <= cnte; ++i) {
    Lk[e[i][0]].emplace_back(e[i][1], e[i][2]);
    Lk[e[i][1]].emplace_back(e[i][0], e[i][2]);
  }
  NodeTp = 0;
  dfs(e[1][0], 0, e[1][0]);
  for (int i = 0; i < NodeTp; ++i)
    debug("node = %d, to = %d, tof = %d
", NodeStk[i], To[NodeStk[i]], Tof[NodeStk[i]]);
  for (int i = 0; i < NodeTp; ++i) {
    int id = CntF + (i % NodeTp) + 1;
    F[id] = triangle_t(insertNode, NodeStk[i], To[NodeStk[i]], Tof[NodeStk[i]], CntF + ((i + 1) % NodeTp) + 1, CntF + ((i - 1 + NodeTp) % NodeTp) + 1);
    if (Tof[NodeStk[i]] != 0) {
      int y = Tof[NodeStk[i]];
      if (F[y].la == -1 && cmp(NodeStk[i], To[NodeStk[i]], F[y].b, F[y].c)) F[y].la = id;
      if (F[y].lb == -1 && cmp(NodeStk[i], To[NodeStk[i]], F[y].a, F[y].c)) F[y].lb = id;
      if (F[y].lc == -1 && cmp(NodeStk[i], To[NodeStk[i]], F[y].a, F[y].b)) F[y].lc = id;
    }
  }
  CntF += NodeTp;
  for (int i = 0; i < NodeTp; ++i) Lk[NodeStk[i]].clear();
  for (int i = 1; i <= CntF; ++i)
    debug("i = %d, a = %d, b = %d, c = %d, la = %d, lb = %d, lc = %d, del = %d
", i, F[i].a, F[i].b, F[i].c, F[i].la, F[i].lb, F[i].lc, Del[i]);
}

int CntE, Par[MaxN + 5];
edge_t ME[MaxF + 5];

int find(int x) { return x == Par[x] ? x : Par[x] = find(Par[x]); }

void solve() {
  for (int i = 1; i <= N; ++i) insert(i);
  for (int i = 1; i <= CntF; ++i)
    if (Del[i] == false) debug("i = %d, a = %d, b = %d, c = %d, l = (%d, %d, %d)
", i, F[i].a, F[i].b, F[i].c, F[i].la, F[i].lb, F[i].lc);
  for (int i = 1; i <= N; ++i) Par[i] = i;
  for (int i = 1; i <= CntF; ++i) {
    if (F[i].a <= N && F[i].b <= N) ME[++CntE] = edge_t(F[i].a, F[i].b, mod(A[F[i].a] - A[F[i].b]));
    if (F[i].a <= N && F[i].c <= N) ME[++CntE] = edge_t(F[i].a, F[i].c, mod(A[F[i].a] - A[F[i].c]));
    if (F[i].b <= N && F[i].c <= N) ME[++CntE] = edge_t(F[i].b, F[i].c, mod(A[F[i].b] - A[F[i].c]));
  }
  std::sort(ME + 1, ME + 1 + CntE);
  double ans = 0;
  for (int i = 1; i <= CntE; ++i) {
    int p = find(ME[i].u), q = find(ME[i].v);
    if (p != q) {
      Par[p] = q;
      ans += ME[i].w;
    }
  }
  printf("%.2lf
", ans);
}

int main() {
#ifdef Tweetuzki
  freopen("input.txt", "r", stdin);
  freopen("output.txt", "w", stdout);
  freopen("errorfile.txt", "w", stderr);
#endif
  init();
  solve();
  return 0;
}

Guibas-Stolfi：

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>

const int MaxN = 5000, MaxM = 100000;
const double eps = 1e-9;

#ifdef Tweetuzki
#define debug(arg...) fprintf(stderr, arg)
#else
#define debug(arg...) void(0)
#endif

typedef struct vec_t {
  double x, y;
  vec_t(double _x = 0, double _y = 0) { x = _x, y = _y; }
  inline friend vec_t operator+(const vec_t &a, const vec_t &b) { return vec_t(a.x + b.x, a.y + b.y); }
  inline friend vec_t operator-(const vec_t &a, const vec_t &b) { return vec_t(a.x - b.x, a.y - b.y); }
  inline friend vec_t operator*(const vec_t &a, double k) { return vec_t(a.x * k, a.y * k); }
  inline friend double dot(const vec_t &a, const vec_t &b) { return a.x * b.x + a.y * b.y; }
  inline friend double cross(const vec_t &a, const vec_t &b) { return a.x * b.y - a.y * b.x; }
  inline friend double mod(const vec_t &a) { return sqrt(a.x * a.x + a.y * a.y); }

  inline void shake() {
    if (rand() % 2 == 0) x += 1.0 * rand() / RAND_MAX * (1e-9);
    else x -= 1.0 * rand() / RAND_MAX * (1e-9);
    if (rand() % 2 == 0) y += 1.0 * rand() / RAND_MAX * (1e-9);
    else y -= 1.0 * rand() / RAND_MAX * (1e-9);
  }
} node_t;

typedef struct vec3_t {
  double x, y, z;
  vec3_t(double _x = 0, double _y = 0, double _z = 0) { x = _x, y = _y, z = _z; }
  inline friend vec3_t operator+(const vec3_t &a, const vec3_t &b) { return vec3_t(a.x + b.x, a.y + b.y, a.z + b.z); }
  inline friend vec3_t operator-(const vec3_t &a, const vec3_t &b) { return vec3_t(a.x - b.x, a.y - b.y, a.z - b.z); }
  inline friend vec3_t cross(const vec3_t &a, const vec3_t &b) { return vec3_t(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); }
  inline friend double dot(const vec3_t &a, const vec3_t &b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
} node3_t;

struct Graph {
  int cnte;
  int head[MaxN + 5], to[MaxM * 2 + 5];
  int pre[MaxM * 2 + 5], nxt[MaxM * 2 + 5];

  Graph() {
    cnte = 1;
    memset(head, 0, sizeof head);
    memset(to, 0, sizeof to);
    memset(pre, 0, sizeof pre);
    memset(nxt, 0, sizeof nxt);
  }

  inline void addEdge(int u, int v) {
    cnte++;
    to[cnte] = v;
    nxt[cnte] = head[u];
    pre[head[u]] = cnte;
    head[u] = cnte;
  }

  inline void delEdge(int id) {
    if (pre[id] != 0) nxt[pre[id]] = nxt[id];
    if (nxt[id] != 0) pre[nxt[id]] = pre[id];
    if (head[to[id ^ 1]] == id) head[to[id ^ 1]] = nxt[id];
  }
};

struct edge_t {
  int u, v;
  double w;
  edge_t(int _u = 0, int _v = 0, double _w = 0) { u = _u, v = _v, w = _w; }
  inline friend bool operator<(const edge_t &a, const edge_t &b) { return a.w < b.w; }
};

int N, M;
node_t A[MaxN + 5], MemoryA[MaxN + 5];
int Indice[MaxN + 5];
int Par[MaxN + 5];
edge_t E[MaxM + 5];
Graph Gr;

inline node3_t mapping(const node_t &a) { return vec3_t(a.x, a.y, a.x * a.x + a.y * a.y); }

inline bool inCircumcircle(const node_t &a, const node_t &b, const node_t &c, const node_t &p) {
  node3_t _a = mapping(a), _b = mapping(b), _c = mapping(c), _p = mapping(p);
  if (cross(b - a, c - a) < 0) std::swap(_b, _c);
  node3_t normal = cross(_b - _a, _c - _a);
  if (dot(normal, _p - _a) > eps) return false;
  else return true;
}

inline bool intersect(const node_t &a, const node_t &b, const node_t &c, const node_t &d) {
  if (cross(c - a, b - a) * cross(b - a, d - a) < eps) return false;
  if (cross(a - d, c - d) * cross(c - d, b - d) < eps) return false;
  return true;
}

void init() {
  srand(20040313U);
  scanf("%d", &N);
  for (int i = 1; i <= N; ++i) {
    scanf("%lf %lf", &A[i].x, &A[i].y);
    MemoryA[i] = A[i];
    A[i].shake();
  }
}

inline bool levelCompare(int x, int y) { return A[x].x < A[y].x; }

inline std::pair<int, int> getInitialBaseEdge(int l, int r) {
  int m = (l + r) >> 1;
  static int stk[MaxN + 5];
  int tp = 0;
  stk[++tp] = l, stk[++tp] = l + 1;
  for (int i = l + 2; i <= r; ++i) {
    while (tp > 1 && cross(A[Indice[stk[tp]]] - A[Indice[stk[tp - 1]]], A[Indice[i]] - A[Indice[stk[tp]]]) < eps) tp--;
    stk[++tp] = i;
  }
  for (int i = 1; i < tp; ++i)
    if (stk[i] <= m && stk[i + 1] > m) return std::make_pair(Indice[stk[i]], Indice[stk[i + 1]]);
  return std::make_pair(-1, -1);
}

void fun(int l, int r) {
  debug("fun(%d, %d)
", l, r);
  if (l == r) return;
  if (l + 1 == r) {
    Gr.addEdge(Indice[l], Indice[r]);
    Gr.addEdge(Indice[r], Indice[l]);
    return;
  }
  int m = (l + r) >> 1;
  fun(l, m), fun(m + 1, r);
  std::pair<int, int> base = getInitialBaseEdge(l, r);
  for (;;) {
    Gr.addEdge(base.first, base.second);
    Gr.addEdge(base.second, base.first);
    int leftFinal = -1, rightFinal = -1;
    for (int i = Gr.head[base.first]; i; i = Gr.nxt[i]) {
      int candidate = Gr.to[i];
      if (cross(A[base.second] - A[base.first], A[candidate] - A[base.first]) < eps) continue;
      if (leftFinal == -1 || inCircumcircle(A[base.first], A[base.second], A[leftFinal], A[candidate]) == true)
        leftFinal = candidate;
    }
    for (int i = Gr.head[base.second]; i; i = Gr.nxt[i]) {
      int candidate = Gr.to[i];
      if (cross(A[candidate] - A[base.second], A[base.first] - A[base.second]) < eps) continue;
      if (rightFinal == -1 || inCircumcircle(A[base.first], A[base.second], A[rightFinal], A[candidate]) == true)
        rightFinal = candidate;
    }
    if (leftFinal == -1 && rightFinal == -1) break;
    if (leftFinal != -1 && rightFinal != -1) {
      if (inCircumcircle(A[base.first], A[base.second], A[leftFinal], A[rightFinal]) == true) leftFinal = -1;
      else rightFinal = -1;
    }
    if (leftFinal != -1) {
      for (int i = Gr.head[base.first]; i; i = Gr.nxt[i]) {
        int linknode = Gr.to[i];
        if (intersect(base.second, leftFinal, base.first, linknode) == true) {
          Gr.delEdge(i);
          Gr.delEdge(i ^ 1);
        }
      }
      base = std::make_pair(leftFinal, base.second);
    } else {
      for (int i = Gr.head[base.second]; i; i = Gr.nxt[i]) {
        int linknode = Gr.to[i];
        if (intersect(base.first, rightFinal, base.second, linknode) == true) {
          Gr.delEdge(i);
          Gr.delEdge(i ^ 1);
        }
      }
      base = std::make_pair(base.first, rightFinal);
    }
  }
}

int find(int x) { return x == Par[x] ? x : Par[x] = find(Par[x]); }

void solve() {
  for (int i = 1; i <= N; ++i) Indice[i] = i;
  std::sort(Indice + 1, Indice + 1 + N, levelCompare);
  fun(1, N);

  for (int u = 1; u <= N; ++u)
    for (int i = Gr.head[u]; i; i = Gr.nxt[i]) {
      int v = Gr.to[i];
      if (u < v) {
        debug("%d - %d
", u, v);
        E[++M] = edge_t(u, v, mod(MemoryA[v] - MemoryA[u]));
      }
    }
  double ans = 0;
  std::sort(E + 1, E + 1 + M);
  for (int i = 1; i <= N; ++i) Par[i] = i;
  for (int i = 1; i <= M; ++i) {
    int u = E[i].u, v = E[i].v;
    double w = E[i].w;
    int p = find(u), q = find(v);
    if (p != q) {
      ans += w;
      Par[p] = q;
    }
  }
  printf("%.2lf
", ans);
}

int main() {
#ifdef Tweetuzki
  freopen("input.txt", "r", stdin);
  freopen("output.txt", "w", stdout);
  freopen("errorfile.txt", "w", stderr);
#endif
  init();
  solve();
  return 0;
}