1.6 三角分割实现
Triangulation: Implementation

1.6.1 对角线,内部或者外部(Diagnoals, Internal or External)

完成了判定线段相交的实现之后，我们距离实现多边形的三角分割就不远了。我们的第一个目标是要找到多边形的一条对角线。

引理 1.5.1 刻画了对角线的两个条件，不与多边形的边相交、在多边形的内部。如果我们忽略内部对角线和外部对角线的差异，那么寻找对角线的过程就是简单地重复调用 Intersect 函数。对于多边形中每条边 $e$，如果它不恰好与对角线 $s$ 的两个端点之一重合，就检查 $e$ 是否与 $s$ 相交。一旦检测到交点，即可确定 $s$ 不是对角线。如果没有这样的边与 $s$ 相交，那么 $s$ 可能就是对角线。我们不能完全下结论的原因是，与 $s$ 的某个端点相连的边可能与 $s$ 共线，而这种情况是不会被检测到的。我们稍后会讨论这种可能性。Code 1.10 给出了用于检查相交的代码实现。

展开 Code 1.10 Diagonalie.>>>

        bool Diagnalie(tVertex a, tVertex b)
        {
            tVertex c, c1;

            c = vertices;
            do {
                c1 = c->next;
                if ((c != a) && (c1 != a) && (c != b) && (c1 != b) && Intersect(a->v, b->v, c->v, c1->v))
                    return FALSE;
                c = c->next;
            } while (c != vertices);
            retrun TRUE;
        }

1.6.2 InCone

现在我们来看引理 1.5.1 的第二个条件。我们需要区分出对角线实在多边形的内部还是外部。同时我们还需要考虑与对角线共用端点的边。我们使用一个布尔函数 InCone 来解决这个两个问题。该函数可以确定一个向量 $B$ 是否严格位于另外两个向量 $A$ 和 $C$ 之间的逆时针方向的开锥里(open cone counterclockwise)。如果向量 $A$ 和 $C$ 是多边形上相邻的两条边，那么 $B$ 就在对角线上。这个方法足以确定对角线了，我们会在后面介绍细节。现在我们来设计函数 InCone。

如果圆锥的顶点是凸角(convex angle, ≤ 180° 的角)，这将是一个简单的任务。但如果它是凹角，就需要分情况处理或者需要一些巧妙的方法。这里我们先分类处理，而将巧妙的方法留给练习 1.6.4[5]。

图 1.25(a) 展示了凸角的情况。当且仅当 $s$ 的一个端点位于顶点 $a$ 处并且在与 $a-, a+$ 构成的圆锥内部，$s$ 才在 $P$ 的内部。我们可以通过函数 Left 轻松搞定: $a-$ 必须在 $ab$ 的左侧，$a+$ 必须在 $ba$ 的左侧。这两个左侧的判定都必须严格的，以避免 $ab$ 与圆锥边界共线重叠的情况。

如图 1.25(b) 所示，当 a 为凹角时，刚才的两个条件就不充分了。$a-$ 和 $a+$ 可能都在内部对角线的左侧，也可能都在右侧，或者一个在左侧而另一个在右侧。值得注意的是，$a$ 邻域的外部是凸角锥体，如果将凹角顶点的内部和外部互换，就变成了凸顶点。因此在这种情况下，判定 $s$ 是否在内部的最简单方法就是，判定它不在椎体的外部。 $a+$ 不能在 $ab$ 的左侧或 $ab$ 上，$a-$ 不能同时在 $ba$ 的左侧或 $ba$ 上。

最后，区分一个角是凸角还是凹角很简单，只需要调用一次 Left 函数就可以确定。当且仅当 $a-$ 在 $aa+$ 的左侧或 $aa+$ 上时，$a$ 就是一个凸角。注意，如果 $a-, a, a+$ 共线，那么 $a$ 处的内角就是 $\pi$，我们仍将其定义为凸角(1.1.2 节)。 Code 1.11 就是上述思想的直接实现。

展开 Code 1.11 InCone.>>>

        bool InCone(tVertex a, tVertex b)
        {
            tVertex a0 = a->prev;
            tVertex a1 = a->next;
            if (LeftOn(a->v, a1->v, a0->v)) {
                // 若 a 是凸顶点
                return Left(a->v, b->v, a0->v) && Left(b->v, a->v, a1->v);
            } else {
                // 否则 a 是一个凹顶点
                return !(LeftOn(a->v, b->v, a1->v) && LeftOn(b->v, a->v, a0->v));
            }
        }

虽然这个 InCone 函数很简单，但很容易出错。整个函数由五个有符号面积计算组成，这恰恰说明了该计算的实用性。按：实在没看出实用性的因果关系。

1.6.3 对角(Diagnoal)

现在，我们已经有了判断 ab 是否为对角线的代码，当且仅当 Diagonalie(a, b)、InCone(a, b) 和 InCone(b, a) 都为真时，ab 才是对角线。 InCone 函数的调用既可以确保 ab 在多边形内部，也可以覆盖 Diagonalie 函数未检查的与端点相连的边。此外，函数的调用是有顺序的。InCone 函数的调用应该放在最前面，因为它们都是常数时间，只涉及顶点 $a,b$ 的邻域，无需考虑多边形的其余部分。而 Diagonalie 函数有一个遍历所有 n 条边的循环。如果任何一个 InCone 函数的调用返回 FALSE，就不会执行可能耗时较长的 Diagonalie 检查。参见Code 1.12。

展开 Code 1.12 Diagnal.>>>

        bool Diagnal(tVertex a, tVertex b)
        {
            return InCone(a, b) && InCone(b, a) && Diagonalie(a, b);
        }

1.6.4 练习

整数溢出integer overflow。如果一台机器上的 int 是 32 位的，$a,b,c$ 的坐标有多大不会导致 Area2 计算溢出？ Code 1.5
函数 IntersectProp。详细分析函数 IntersectProp(Code 1.7)，如果删除它的 if 语句的判定，该函数计算的是什么。探讨一下，删除之后，函数 Intersect(Code 1.9)还能不能正常工作？
函数 Intersect 的低效率。手动追踪一下函数 Intersect(Code 1.9)，看看它最多调用了几次函数 Area2(Code 1.5)，是否存在冗余的调用，设计一个新的函数避免这些冗余的调用。
保存交叉信息Saving intersection information。设计一个新算法，避免对同两条线段进行两次交点测试。分析新算法的时间复杂度和空间复杂度。
改进 InCone (Andy Mirzian)。证明当且仅当 Left(a, a+, b), Left(a,b,a-), Left(a, a-, a+) 这三个 Left 中有一个返回 false 时，$ab$ 在顶点 $a$ 的椎体内。
改进 Diagonal。证明 Diagonal 函数中对 InCome 的两次调用，移除任何一次都不会改变结果。

1.6.5 通过移除耳朵实现的三角分割(Triangulation by Ear Removal)

现在我们实现多边形的三角分割。一种方法是参考三角分割定理(定理 1.2.3)的证明，先找到一条对角线，将多边形分割成两部分，然后递归每一个部分。这种方法的效率很低，最终我们基于 Meister 的双耳定理(定理 1.2.7)实现。但首先，我们还是分析一下递归方法的运行速度。我们假设读者对 $O$ 表示法很熟悉。

基于对角线的算法(Diagonal-Based Algorithm)

基于定理 1.2.3 的算法的复杂度是 $O(n^4)$。其中有 $\binom{n}{2} = O(n^2)$ 个候选对角线，而每个对角线的判定也需要 $O(n)$。对 $n - 3$ 个对角线重复上述 $O(n^3)$ 的计算，就得到了一个 $O(n^4)$ 的算法。

我们可以利用双耳定理将速度提高 $n$ 倍。一定存在一条内部对角线将一个耳朵分隔开。而这种"耳朵对角线ear diagonal"只有 $O(n)$ 的候选对象，$(v_i, v_{i+2}), i = 0, \cdots, n-1$。这样递归过程也简单，因为一次只需要处理其中一个部分。另一个部分就是一个三角形耳朵，显然已经是三角分割。因此这种方法最坏情况下的复杂度是 $O(n^3)$。

移除耳朵(Ear Removal)

现在我们将上述算法改进到 $O(n²)$。因为每次调用 Diagonal 函数的复杂度为 $O(n)$，为了达到 $O(n²)$，Diagonal 函数最多只能被调用 $O(n)$。改进的关键在于，移除一只耳朵并不会显著改变多边形的形状，尤其是潜在的耳尖。所以我们首先要判断每个顶点 $v_i$ 是否是潜在的耳尖，即 $v_{i-1}$ 到 $v_{i+1}$ 是否构成对角线。虽然这一步骤已经是 $O(n²)$ 的复杂度了，但无需重复执行这个耗时的步骤。

设 $(v_0, v_1, v_2, v_3, v_4)$ 为 $P$ 的五个连续顶点，其中 $v_2$ 是耳尖，并且已经删除了 $E_2 = \Delta(v_1, v_2, v_3)$，如图 1.26 所示。哪些顶点的耳尖状态可能会改变？只有 $v_1$ 和 $v_3$。

比如说 $v_4$。它是不是耳尖取决于 $v_3 v_5$ 是不是对角线。移除 $E_2$ 后 $v_3 v_5$ 的端点保持不变。显然，如果这两个端点本来能够互相看到，那么移除 $E_2$ 也不会阻挡它们之间的视线。如果它们之前就无法互相看到，那么移除 $E_2$ 后它们仍然无法互相看到，这一点可能不太明显。

正如引理 1.2.2 的证明所述，只有两种情况需要考虑。如果 $v_3, v_5$ 是外部线段，那么显然删除 $E_2$ 不会使其变为内部线段。否则，$\Delta(v_3, v_4, v_5)$ 必须包含一个凹顶点(图 1.12 中的 x)。但是移除 $E_2$ 只会移除一个凸顶点。因此，移除 $E_2$ 不会改变 $v_4$ 的状态，除了 $v_1$ 和 $v_3$ 之外的所有顶点的状态也不会改变，因为 $v_1$ 和 $v_3$ 的耳朵对角线(ear diagnals)都与被移除的顶点 V2 相连。

这意味着，在耗时的初始化步骤之后，每次迭代只需调用两次 Diagonal 函数即可更新耳尖状态信息。由此可得到Algorithm 1.1 中所示的用于构建三角剖分的伪代码。

我们将任务"多边形的三角分割"换个说过——"以任意顺序输出构成三角分割的对角线"。这主要是为了便于表述。通常，特定应用需要更结构化的输出。例如，对偶图中需要三角形邻接信息。虽然获得更结构化的输出在渐近时间复杂度上并不更难，但通常会显著增加代码的复杂性。我们将不再深入探讨这些替代的三角剖分输出。

三角分割代码(Triangulation Code)

首先需要初始化布尔标志 v->ear，它是顶点数据结构的一部分(Code 1.2)。对每个顶点调用一次函数 Diagonal 即可。参见 Earlnit，Code 1.13。

展开 Code 1.13 EarInit.>>>

        void EarInit(void)
        {
            tVertex v0, v1, v2;
            v1 = vertices;
            do {
                v2 = v1->next;
                v0 = v1->prev;
                v1->ear = Diagnal(v0, v2);
                v1 = v1->next;
            } while (v1 != vertices);
        }

三角分割代码的主题是一个两层循环。外循环每次迭代移除一个耳朵，直到 $n = 3$ 时停止。内循环通过检查预先计算的 v2->ear 标志来寻找耳朵，其中 $v_2$ 是潜在的耳朵尖。一旦找到耳朵尖，就通过调用函数 Diagonal 更新 $v_1$ 和 $v_3$ 的耳朵状态，打印代表耳朵根的对角线，并将耳朵从多边形中移除。移除操作通过重新连接 $v_1$ 和 $v_3$ 的 next 和 prev 指针来实现。此时，如果 $v_2$ 所指的顶点对象未被其他应用程序使用，则可以释放它。必须注意，要避免 $v_2$ 成为顶点列表头，即访问循环列表的入口点。因此需要把表头的指针 vertices 指向 v_3。参见代码 1.14。在分析时间复杂度之前，我们先逐步推演一个示例。

展开 Code 1.14 Triangulate.>>>

        void Triangulate(void)
        {
            tVertex v0, v1, v2, v3, v4; /* 五个连续的顶点 */
            int n = nvertices;
            EarInit();
            while (n > 3) {  // 外层循环，每次移除一个耳朵
                v2 = vertices;
                do {         // 内层循环，寻找耳朵尖
                    if (v2->ear) {
                        v3 = v2->next; v4 = v3->next;
                        v1 = v2->prev; v0 = v1->prev;
                        PrintDiagonal(v1, v3) // (v1, v3) 是一个对角线

                        // 更新对角线端点的状态
                        v1->ear = Diagonal(v0, v3);
                        v3->ear = Diagonal(v1, v4);

                        // 移除耳朵 v2
                        v1->next = v3;
                        v3->prev = v1;
                        vertices = v3; // 防止 v2 成为链表头
                        n--;
                        break;
                    }
                    v2 = v2->next;
                } while (v2 != vertices);
        }   }

1.6.6 例子

图 1.27 中显示了一个多边形以及由简单的 main 函数 Code 1.15 生成的三角分割。读取和打印的代码很简单，这里不再赘述。

展开 Code 1.15 main.>>>

        main()
        {
            ReadVertices();
            PrintVertices();
            Triangulate();
        }

Table 1.1. Figure1.27 中多边形各个顶点坐标

i	(x,y)	i	(x,y)	i	(x,y)	i	(x,y)
0	(0, 0)	5	(10, 12)	10	(10, 15)	15	(5, 8)
1	(10, 7)	6	(12, 14)	11	(7, 17)	16	(-2, 9)
2	(12, 3)	7	(14, 9)	12	(0, 16)	17	(5, 5)
3	(20, 8)	8	(8, 10)	13	(1, 13)
4	(13, 17)	9	(6, 14)	14	(3, 15)

第 11 的点坐标原文是 (7, 10) 不构成一个多边形，这里改成了 (7,17)。

Table 1.2. 依次输出的对角线顶点对

order	对角线	order	对角线	order	对角线	order	对角线
1	(17, 1)	5	(9, 11)	9	(15, 3)	13	(15, 8)
2	(1, 3)	6	(12, 14)	10	(3, 7)	14	(15, 9)
3	(4, 6)	7	(15, 17)	11	(11, 14)	15	(9, 14)
4	(4, 7)	8	(15, 1)	12	(15, 7)

现在来看一下，本例中对角线的输出结果，如表 1.2 所示。$v_0$ 是耳尖，因此第一个对角线输出为 (17, 1)。 $v_1$ 不是耳尖，所以指针 v2 移动到 $v_2$。$v_2$ 是一个耳尖，输出对角线 (1, 3)。 $v_3$ 和 $v_4$ 都不是耳尖，直到 $v_5$ 时才输出下一个对角线 (4, 6)。线段 $v_3v_8$ 与 $v_7$ 共线，因此下一个耳朵要到 $v_{10}$ 才能找到。图 1.27 中深色阴影的子多边形显示了输出第 9 条对角线 (15, 3) 后剩余的多边形。另外，$v_9$ 与 $v_{11}v_{15}$ 共线，因此在移除对角线 (15, 9) 后，$v_9$ 不再是耳朵。

1.6.7 分析

现在我们分析算法的时间复杂度。如前所述，EarInit 的复杂度为 O(n²)。Triangulate 函数的外层循环迭代次数等于对角线的数量，即 $n - 3 = O(n)$。内层循环搜索耳朵的复杂度也为 $O(n)$，因为可能需要检查每个顶点。内层循环内部的工作量为 $O(n)$。在最坏情况下，可能需要需要对整个多边形调用两次函数 Diagonal 检查对角线是否被阻塞。因此，我们得到的算法时间复杂度为 $O(n³)$，不足预期的 $O(n²)$。更深入的分析表明，$O(n²)$ 才是 Triangulate 函数的正确时间复杂度界限。

如图 1.28 所示。删除 $v_0$ 后，内循环需要搜索 $v_1, \cdots, v_6$ 才能找到下一个耳尖 $v_7$。然后，它必须搜索 $v_8, \cdots, v_{12}$ 才能找到耳尖 $v_{13}$。这个例子表明，内循环确实可能需要迭代 $n$ 次才能找到一只耳朵。但请注意，循环内的两个 $O(n)$ 的函数 Diagnal 仅在找到一只耳朵后才需要调用——它们并非在每次迭代中都调用。因此，尽管表面上其复杂度为 $n × n × n = O(n³)$，但实际上它是 $n × (n + n) = O(n²)$ 的。

还可以更进一步(练习 1.6.8(4))，将渐近时间复杂度降低到二次方以下，需要采用完全不同的方法，这将在下一章中讨论。

1.6.8 练习

[Programming] 重复相交测试Repeated intersection tests。函数 Triangulate (Code 1.14) 经常需要检查线段与线段之间是否相交。修改代码，以便确定进行了多少次不必要的相交检查。在图 1.27 上进行测试。
[easy] 凸多边形Convex polygons。分析函数 Triangulate 在凸多边形上的运行性能。
螺旋Spiral。继续分析图 1.28，函数 Triangulate 是否会继续遍历边界寻找耳朵尖？更具体地说，如果多边形有 $n$ 个顶点，在完成之前，指针 v2 会沿着边界完整地循环多少次？
[Programming] 耳朵列表Ear list。通过将耳尖链接到一个单独的(循环)列表中，可以避免函数 Triangulate 的内循环搜索。该列表将满足 v->ear == TRUE 的顶点 v 与顶点结构中的指针 nextear 和 prevear 连接起来。迭代的过程中，可以通过该列表找到下一次迭代的目标耳尖。实现此改进，并验证其速度提升是否明显。
重心center of gravity。设计一个算法来计算多边形的重心，假设该多边形由均匀密度的材料切割而成。重心是一个点，可以用一个向量表示。三角形的重心位于其质心(centroid)处，其坐标恰好是三角形各顶点坐标的平均值。两个不相交的集合 $A,B$ 的并集 $S$，其重心 $\gamma(S)$ 是 $A$ 和 $B$ 两个部分重心的加权和。记 $S$ 的重量为 $\omega(S) = \omega(A) + \omega(B)$，那么 $$ \gamma(S) = \frac{\omega(A) \gamma(A) + \omega(B) \gamma(B)}{\omega(S)} $$ 这里，每个三角形的权重是其在均匀密度假设下的面积。

1.6 三角分割实现Triangulation: Implementation