5.7 与凸包的联系
(Connection to Convex Hulls)

1986年，Edelsbrunner 和 Seidel 发现了 Delaunay 三角分割与高出一维的凸包之间存在着美妙的联系。我将首先解释二维凸包与一维 Delaunay 三角分割之间的联系(诚然，这是微不足道的)， 然后将其推广到二维 Delaunay 三角分割与三维凸包。这种联系将为我们提供一种通过三维凸包来计算 Delaunay 三角分割，进而得到维诺图的简便方法。

5.7.1. 一维 Delaunay 三角分割 (One-Dimensional Delaunay Triangulations)

我们从一维开始，这里的数学原理是显而易见的。

设 \(P = \{x_1, \dots, x_n\}\) 是 x 轴上的一组点。一维 Delaunay 三角分割仅仅是连接 \(x_1\) 到 \(x_2\) 到 \(\dots\) 直到 \(x_n\) 的路径。 我们将把它视为一组坐标为 \((x_i, x_i^2)\) 的二维点集在 \(x\) 轴上的投影。这些点可以看作是 \(x_i\) 向上投影到抛物线 \(z = x^2\) 上的结果。 一个显而易见的事实是，这些二维点的凸包(convex hull)向下投影就是一维 Delaunay 三角分割，只要丢弃凸包的"顶部"边缘即可。 除了这一简单的现象之外还有很多其它内容，通过考虑抛物线的切线可以阐明更多深层含义。

抛物线 \(z = x^2\) 在点 \(x = a\) 处的斜率为 \(2a\)，因为 \(dz/dx = 2x\)。因此，抛物线在点 \((a, a^2)\) 处的切线方程为

为了在三维空间中采用相同的分析过程，我们在垂直方向上，将这条切线平移 \(r^2\) 的距离，然后研究它与抛物线的交点。移动后的切线方程变为

这条直线在哪里与抛物线相交呢?

因此，移动后的切线在原切点 \(a\) 左右 \(\pm r\) 的位置与抛物线相交。注意，\(x = a \pm r\) 可以被看作是中心在 \(a\) 半径为 \(r\) 的一维圆(即线段)的方程。 如图 5.23 所示，其中 \(a = 5, r = 3\)，因此该"圆盘"即为线段 \([2, 8]\)。

5.7.2. 二维 Delaunay 三角分割 (Two-Dimensional Delaunay Triangulations)

现在我们在二维空间中重复同样的分析。

如图 5.24 所示，抛物面方程为 \(z = x^2 + y^2\)。取平面中给定的站点，并将它们向上投影直到碰到抛物面，按如下方式映射每个点:

取这组三维点的凸包(convex hull)，如图 5.25 所示。现在丢弃该凸包的"顶部"面，即，所有外法向量指向上方的面，所谓的指向上方是说与 \(z\) 轴向量的点积为正。 结果是一个底部的"壳"。将其投影到 \(xy\) 平面上，所得到的就是 Delaunay 三角剖分! 如图 5.26 所示。我们现在形式化的描述这一惊人的联系。

在点 \((a, b)\) 上方的切平面方程为

这直接类比于方程 \(z = 2ax - a^2\)，因为 \(\partial z/\partial x = 2x\) 且 \(\partial z/\partial y = 2y\)。他想说啥。 现在将这个平面向上平移 \(r^2\)，就像我们在上一小节中将切线向上平移一样：

还是那个问题，这个平移后的平面在哪里与抛物面相交?

平移后的平面与抛物面相交于一条曲线，该曲线投影是一个圆! 如图 5.27 和 5.28 所示。

现在我们反转视角以引出 Delaunay 三角分割。假设平面 \(\pi\) 通过抛物面上三个点 \(\Delta = (p_i, p_j, p_k)\)，这三个点构成了三维凸包的一个面。 该平面穿过抛物面。如果我们把 \(\pi\) 垂直向下平移，在某一点它将不再与抛物面相交。假设它接触的最后一点是 \((a, b, a^2 + b^2)\)。 那么我们可以把 \(\pi\) 看作是该切平面 \(\tau\) 向上平移的结果，设平移量为 \(r^2\)。

由于 \(\Delta\) 是凸包的一个下表面，抛物面上的所有其他点都在 \(\pi\) 之上。由于它们在 \(\pi\) 之上，它们距离 \(\tau\) 的高度超过 \(r^2\)， 因为 \(\tau\) 在 \(\pi\) 下方 \(r^2\) 处。因此，这些点将被投影到 \(xy\) 平面上半径为 \(r\) 的圆之外，即，由 \(\Delta\) 在 \(xy\) 平面上确定的圆不包含所有其他站点。 所以它形成了一个 Delaunay 三角形。凸包的每个下三角面都对应一个 Delaunay 三角形，所以凸包"底部"的投影即为 Delaunay 三角分割! 请再次参考图 5.26。

让我再用另一种方式解释这一重要结论。从平面 \(\tau\) 开始，该平面与抛物面在点 \(p = (a, b)\) 上方相切，解除点向下投影到 \(p\)。 现在将 \(\tau\) 向上移动。它与抛物面交线的投影是一个以 \(p\) 为中心扩张的圆。当 \(\tau\) 碰到抛物面上位于某个站点正上方的点 \(q\) 时， 这个扩张的圆就会碰到那个站点。因此，在 \(\tau\) 到达 \(\pi\) 之前，该圆是空的。当它到达 \(\pi\) 时，它穿过了三个站点，这三个站点的投影构成了由 \(\pi\) 支撑的三角凸包面 \(\Delta\)。

根据上述讨论，我们可以得到一个有用推论:

这些点的共面性可以通过检查它们所确定的四面体体积(方程(1.12)和 式(4.1))是否为零来验证。

5.7.3 内涵 (Implications)

由于三维凸包可以在 \(O(n \log n)\) 时间内计算出来，参见4.2.2节， 这意味着 Delaunay 三角分割也可以在同样的时间内计算出来。一旦得出了 Delaunay 三角分割，计算维诺图就相对容易了(练习 5.7.5[2])。这正是另一种构建维诺图的 \(O(n \log n)\) 算法。

正如人们所希望的那样，维诺图与高一维空间中的凸包之间的这种关系，在任意维度都成立。因此，三维空间中的维诺图图和 Delaunay 三角分割都可以通过四维空间中的凸包来构建。 事实上，四维凸包代码最常见的用途可能就是用于构建 Delaunay 四面体的实体网格。一般而言，\(d\) 维点集的对偶维诺图是 \(d+1\) 维点集的"下"凸包的投影。

5.7.4 Delaunay 三角分割的实现 \(O(n^4)\) (Implementation of Delaunay Triangulation: \(O(n^4)\) Code)

如果不考虑时间复杂度的话，定理 5.7.2 使得计算 Delaunay 三角分割的代码极其简洁。特别是，如果 \(O(n^4)\) 的复杂度是可以接受的(虽然很少见)， 那么 Delaunay 三角分割可以用不到三十行 C 语言代码计算出来! 这在 Code5.1 中给出，部分是出于猎奇，但也为了强调对几何学的深刻理解，能带来如何干净的代码。

        main()
        {
            int x[NMAX], y[NMAX], z[NMAX];   /* 输入点 xy, z=x^2+y^2 */
            int n;                           /* 输入点的数量 */
            int i, j, k, m;                  /* 四个点的索引 */
            int xn, yn, zn;                  /* (i,j,k) 的外向法向量 */
            int flag;                        /* 如果 m 在 (i,j,k) 上方则为真 (true) */
        
            /* 输入点并计算 z = x^2 + y^2。 */
            scanf("%d", &n);
            for ( i = 0; i < n; i++ ) {
                scanf("%d %d", &x[i], &y[i]);
                z[i] = x[i] * x[i] + y[i] * y[i];
            }
        
            /* 对于每个三元组 (i,j,k) */
            for ( i = 0; i < n - 2; i++ ) {
                for ( j = i + 1; j < n; j++ ) {
                    for ( k = i + 1; k < n; k++ ) {
                        if ( j != k ) {
                            /* 计算三角形 (i,j,k) 的法向量。 */
                            xn = (y[j]-y[i])*(z[k]-z[i]) - (y[k]-y[i])*(z[j]-z[i]);
                            yn = (x[k]-x[i])*(z[j]-z[i]) - (x[j]-x[i])*(z[k]-z[i]);
                            zn = (x[j]-x[i])*(y[k]-y[i]) - (x[k]-x[i])*(y[j]-y[i]);
        
                            /* 只检查旋转抛物面底部的面：zn < 0。 */
                            if ( flag = (zn < 0) )
                                for ( m = 0; m < n; m++ )
                                    /* 对于每个其他点 m, 检查 m 是否在 (i,j,k) 上方。 */
                                    flag = flag && ((x[m]-x[i])*xn +
                                                    (y[m]-y[i])*yn +
                                                    (z[m]-z[i])*zn <= 0);
        
                            if (flag)
                                printf("%d\t%d\t%d\n", i, j, k);
            }   }   }   }
        }

代码的 \(O(n^4)\) 结构在四个嵌套的 for 循环中显而易见。对于每组三点 \((i, j, k)\)，程序检查所有其他点 \(m\) 是否都在包含 \(i, j, k\) 的平面上或平面上方。 如果是，\((i, j, k)\) 就被认为是一个 Delaunay 三角形输出。这与算法 3.2 中的二维凸包算法的很相似。 通过将指向三角形外部的法向量 \((x_n, y_n, z_n)\) 与从点 \(i\) 指向点 \(m\) 的向量进行点积，来判定 \(m\) 是否在平面上方。

尽管如此计算 Delaunay 三角分割的代码很简洁有趣，但由于以下几个原因，它不实用:

1. 代码会输出 Delaunay 面内部的所有三角形，而该面的边界可能由四个或更多共圆点组成。比如，一个正方形面会输出四个三角形，代表该正方形的两种三角分割形式。还要后处理一下才能用。
2. 代码中存在明显的低效之处，例如，当发现一个点位于 \(ijk\) 平面下方时，\(m\) 循环可以中断。这些无伤大雅，很容易处理，但代价是代码会稍微变长。
3. \(O(n^4)\) 的时间复杂度是不可接受的。我们曾取 \(n = 100\) 和 \(n = 200\) 进行过两次测试，计算时间分别为 12 秒和 239 秒，显示出点数加倍时， 时间增加了超过 \(2^4 = 16\) 倍。这表明 \(n = 1,000\) 时将需要几天时间。

3D 凸包到 Delaunay 三角分割: \(O(n^2)\) 代码

将第 4 章中开发的用于构建 3D 凸包的平方级代码转换为用于构建 Delaunay 三角分割的平方级代码是一项简单的任务。所有的复杂性都在于凸包代码(约 1,000 行代码)以及定理 5.7.2 中的推理。 只需要很小的修改。首先，ReadVertices(Code4.10)应读取 \(x\) 和 \(y\) 坐标， 并计算 \(z = x^2 + y^2\)。其次，在整个凸包构建完成后，调用 LowerFaces 来遍历所有面并识别哪些面是凸包的底部(Code5.2)。 同 Code5.1 一样，通过计算每个面 \(f\) 的法向量的 \(z\) 坐标来完成。Normz 只是对 Collinear(Code4.12)的微小改动。 如果 Normz(f) < 0，则该面位于凸包的底部，其在 \(xy\) 平面上的投影即为一个 Delaunay 三角形。

        void LowerFaces( void )
        {
            tFace f = faces;
            int Flower = 0;    /* Total number of lower faces. */
        
            do {
                if ( Normz( f ) < 0 ) {
                    Flower++;
                    printf("lower face indices: %d, %d, %d\n",
                           f->vertex[0]->vnum,
                           f->vertex[1]->vnum,
                           f->vertex[2]->vnum );
                }
                f = f->next;
            } while ( f != faces );
            printf("%d lower faces identified. \n", Flower);
        }
        
        int Normz( tFace f )
        {
            tVertex a, b, c;
        
            a = f->vertex[0];
            b = f->vertex[1];
            c = f->vertex[2];
            return
                ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) -
                ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] );
        }

我们以 \(n = 10\) 个点的示例进行说明，其坐标如表 5.1 所示。注意，\(z\) 坐标因平方运算而大幅增加了量级。 这个小示例完全处于第 4.3.5 节讨论的体积计算的安全范围内， 但与 3D 凸包代码相比，平方运算确实缩小了此代码的适用范围。代码的 Postscript 输出如图 5.29 所示。

正如预期的那样，这种 \(O(n^2)\) 的实现比 \(O(n^4)\) 的代码快得多: 在相同的 \(n = 100\) 和 \(n = 200\) 的示例中，\(O(n^4)\) 耗时 12 秒和 239 秒，\(O(n^2)\) 耗时 0.2 秒和 1.4 秒。 此外，按照4.5 节讨论的随机化加速技术可以将其转化为一种期望时间为 \(O(n \log n)\) 的算法。

5.7.5 练习

1. [Programming] dt2.c 的范围 Range of dt2.c。找到一个坐标尽可能小的点集，并对于该点集执行 dt2.c (Code5.2)，使其由于体积计算溢出而返回错误的结果(参见练习 4.3.6[12])。
2. [Programming] \(D(P) \Rightarrow V(P\))。修改 dt2.c 代码，使其根据 Delaunay 三角分割计算维诺图(参见练习 5.5.6[1])。这将需要反复通过三个给定点 \(a, b, c\) 构造圆。 该圆圆心坐标 \(p = (p_0, p_1)\) 可以按如下方式计算: $$ \begin{aligned} A &= b_0 - a_0, \\ B &= b_1 - a_1, \\ C &= c_0 - a_0, \\ D &= c_1 - a_1, \\ E &= A(a_0 + b_0) + B(a_1 + b_1), \\ F &= C(a_0 + c_0) + D(a_1 + c_1), \\ G &= 2(A(c_1 - b_1) - B(c_0 - b_0)), \\ p_0 &= (DE - BF)/G, \\ p_1 &= (AF - CE)/G. \end{aligned} $$ 这些方程略显繁琐的形式是为了减少计算最终坐标所需的乘法次数。如果 \(G = 0\)，那么这三个点共线，就不存在有限半径的圆通过它们。否则，圆的半径为 $$ r^2 = (a_0 - p_0)^2 + (a_1 - p_1)^2 $$ 输出所有维诺顶点的坐标。对于每条有限长度的维诺边，输出其两个端点(可以是它们的坐标，也可以是维诺顶点列表中的索引)。 对于每条无界的维诺边射线，输出其端点以及一个向量，表示沿该射线指向无穷远处的。
3. 最远点维诺图 Furthest-point Voronoi diagram。论证变换后点的凸包的"顶部"是最远点维诺图的对偶图。关于该图的定义，参见练习 5.5.6[11]。 "顶部"面是指那些外法向量具有正 \(z\) 分量的面。因此，从 \(z = +\infty\) 处观察抛物面凸包，显示的是 \(\mathcal{F}(P)\) 的对偶图! 参见图 5.30。
4. 圆形可分性 Circular separability。给定两个平面点集 \(A\) 和 \(B\)，设计一个寻找闭圆盘(如果存在的话)的算法，该圆盘包含 \(A\) 中的每一个点，但不包含 \(B\) 中的任何一个点。