三维凸包的Chan算法

Reference: A minimalist’s implementation of an approximate nearest neighbor algorithm in fixed dimensions, Timothy M. Chan

显然的，问题归约到求三维下凸壳。一个静态的三维问题可以被看成一个动态的二维问题，假设每个点 \((x,y,z)\) 在时间 \(t\) 时变为了二维点 \((x,z-ty)\) 。这组（动）二维点的凸包相当于三维凸包与某个动面的交（可能经过缩放变换）。

在 \(t\) 在 \((-\infty,\infty)\) 运动的过程中，我们考虑计算出这个二维凸包的动画，这族二维凸包就可以方便地求出原先的三维凸包。考虑怎么表示这个二维凸包动画：可以发现每个点都是以（不同的）竖直的速度在二维平面上运动，在运动过程中每个点可能进入（成为）下凸壳或离开下凸壳。把他们分别称作“加入事件”和“删除事件”。可以发现每个事件对应于三维凸包上的一个面。（显然地）每个点最多参与每类事件各一次。好了，现在我们知道如何 \(\mathrm O(n)\) 地存储一个下凸壳动画了（一个 \(\mathrm O(n)\) 的初始状态，\(\mathrm O(n)\) 个按时间排好序的事件）。

接下来考虑怎么求这个动画。我们对 \(x\) 轴排序、分治（注意 \(x\) 是不动的）。 \(n\leq 2\) 的情形是显然的。假设左半边和右半边做好了，考虑两边的合并。首先静态二维情况下，\(x\)轴分离的两个点集 \(L,R\) 的并的凸壳是它们的凸壳，加上两个凸壳的公切线（设 \(u,v\) 分别为其切点），再去掉公切线之上的部分。

我们用 \(x^-\) 表示凸壳上 \(x\) 左侧的邻居， \(x^+\) 表示其右侧的邻居。

动态情况下，我们只需要维护那条公切线和与它相关的事件就好了。我们模拟时间的进行，依次处理每个事件。这里有若干种待选的发生事件，我们分别计算出这几个事件哪个最先发生，然后将时间快进到这个事件发生的时候处理这一事件就好了。可能的待选事件有：

\(u\) 的左侧，\(L\) 上发生了一个插入或删除事件。
\(v\) 的右侧，\(R\) 上发生了一个插入或删除事件。
\(u^-\) 从上到下运动穿过了 \(uv\) 的延长线，此时 \(u\) 被删除，\(u^-v\) 成为新切线。
\(v^+\) 从上到下运动穿过了 \(uv\) 的延长线，同理。
\(u^+\) 从上到下运动穿过了 \(uv\) ，此时 \(u^+\) 被加入，\(u^+v\) 成为新切线。
\(v^-\) 从上到下运动穿过了 \(uv\) ，同理。

可以发现将两个大小相近的凸壳动画合并的复杂度是 \(\mathrm O(n)\) 的，因此分治的总复杂度是 \(\mathrm O(n\log{n})\) 的。

以下这一份是论文原作者提供的代码，不到90个非空行，这份实现将空间卡到了 \(6n\) 个指针，实际上（在稍微损失空间的情况下）代码应当可以更简洁。

// Timothy Chan    "ch3d.cc"    12/02    3-d lower hull (in C++)

// a simple implementation of the O(n log n) divide-and-conquer algorithm

//   input: coordinates of points
//     n x_0 y_0 z_0 ... x_{n-1} y_{n-1} z_{n-1}

//   output: indices of facets
//     i_1 j_1 k_1 i_2 j_2 k_2 ...

//   warning: ignores degeneracies and robustness
//   space: uses 6n pointers

#include <iostream>

struct Point {
  double x, y, z;
  Point *prev, *next;
  void act() {  
    if (prev->next != this) prev->next = next->prev = this;  // insert
    else { prev->next = next; next->prev = prev; }  // delete
  }
};

const double INF = 1e99;
static Point nil = {INF, INF, INF, 0, 0};
Point *NIL = &nil;

inline double turn(Point *p, Point *q, Point *r) {  // <0 iff cw
  if (p == NIL || q == NIL || r == NIL) return 1.0;
  return (q->x-p->x)*(r->y-p->y) - (r->x-p->x)*(q->y-p->y);
}

inline double time(Point *p, Point *q, Point *r) {  
  if (p == NIL || q == NIL || r == NIL) return INF;
  return ((q->x-p->x)*(r->z-p->z) - (r->x-p->x)*(q->z-p->z)) / turn(p,q,r);
}

Point *sort(Point P[], int n) {  // mergesort
  Point *a, *b, *c, head;

  if (n == 1) { P[0].next = NIL; return P; }
  a = sort(P, n/2);
  b = sort(P+n/2, n-n/2);
  c = &head;
  do
    if (a->x < b->x) { c = c->next = a; a = a->next; }
    else { c = c->next = b; b = b->next; }
  while (c != NIL);
  return head.next;
}

void hull(Point *list, int n, Point **A, Point **B) {  // main algorithm
  Point *u, *v, *mid;  
  double t[6], oldt, newt;  
  int i, j, k, l, minl;

  if (n == 1) { A[0] = list->prev = list->next = NIL; return; }

  for (u = list, i = 0; i < n/2-1; u = u->next, i++) ;
  mid = v = u->next;
  hull(list, n/2, B, A);  // recurse on left and right sides
  hull(mid, n-n/2, B+n/2*2, A+n/2*2);

  for ( ; ; )  // find initial bridge
    if (turn(u, v, v->next) < 0) v = v->next;
    else if (turn(u->prev, u, v) < 0) u = u->prev;  
    else break;

  // merge by tracking bridge uv over time
  for (i = k = 0, j = n/2*2, oldt = -INF; ; oldt = newt) {  
    t[0] = time(B[i]->prev, B[i], B[i]->next);  
    t[1] = time(B[j]->prev, B[j], B[j]->next);    
    t[2] = time(u, u->next, v);  
    t[3] = time(u->prev, u, v);
    t[4] = time(u, v->prev, v); 
    t[5] = time(u, v, v->next);
    for (newt = INF, l = 0; l < 6; l++) 
      if (t[l] > oldt && t[l] < newt) { minl = l; newt = t[l]; }
    if (newt == INF) break;
    switch (minl) {
    case 0:  if (B[i]->x < u->x) A[k++] = B[i];  B[i++]->act();  break;
    case 1:  if (B[j]->x > v->x) A[k++] = B[j];  B[j++]->act();  break;
    case 2:  A[k++] = u = u->next;  break;
    case 3:  A[k++] = u;  u = u->prev;  break;
    case 4:  A[k++] = v = v->prev;  break;
    case 5:  A[k++] = v;  v = v->next;  break;
    }
  }
  A[k] = NIL;

  u->next = v;  v->prev = u;  // now go back in time to update pointers
  for (k--; k >= 0; k--) 
    if (A[k]->x <= u->x || A[k]->x >= v->x) {
      A[k]->act();
      if (A[k] == u) u = u->prev; else if (A[k] == v) v = v->next;
    }
    else { 
      u->next = A[k]; A[k]->prev = u; v->prev = A[k]; A[k]->next = v;
      if (A[k]->x < mid->x) u = A[k]; else v = A[k];
    }
}

int main() {
  int n, i;  
  cin >> n;

  Point *P = new Point[n];  // input
  for (i = 0; i < n; i++) { cin >> P[i].x; cin >> P[i].y; cin >> P[i].z; }

  Point *list = sort(P, n);
  Point **A = new Point *[2*n], **B = new Point *[2*n];
  hull(list, n, A, B);

  for (i = 0; A[i] != NIL; A[i++]->act())  // output 
    cout << A[i]->prev-P << " " << A[i]-P << " " << A[i]->next-P << "\n";
  delete A;  delete B;  delete P;
}

EtaoinWu