首页 关于
树枝想去撕裂天空 / 却只戳了几个微小的窟窿 / 它透出天外的光亮 / 人们把它叫做月亮和星星
目录

3.5 Graham 算法
(Graham's Algorithm)

计算几何领域发表的第一篇论文的荣誉或许应该授予 Graham 算法，该算法能够在 \(O(n \log n)\) 时间内找到二维点集的凸包(Graham 1972)。 作者之所以说"或许"，是因为 Toussaint (1985a) 发现了一篇更早的论文，其中包含了许多后来凸包算法的思想。参见 Bass & Schubert(1967) 在 1960 年代末，贝尔实验室有一个应用需要计算 \(n \approx 10,000\) 个点的凸包，他们发现当时使用的 \(O(n^2)\) 的算法速度太慢。 为了满足这一需求，Graham 开发了一个简洁的算法(personal comm. 1992)。

3.5.1 顶层描述(Top Level Description)

Graham 算法的基本思想很简单。我们先作一些假设，用一个例子来解释它，后面会将这些假设去掉。假设我们已知一个点 \(x\) 在凸包内部，并且包含点 \(x\) 的集合中任意三个点都不共线。 如果按照角度逆时针方向对这些点进行排序，如图 3.4 所示，排序后的点被记为 \(a, b, ..., j\)。我们将按序处理这些点，并围绕该集合逐步扩展凸包。 在任何一步，都有一个关于目前已检查过的点的凸包。当然，之后遇到的点会导致之前的判断被重新评估。

我们会在一个栈结构 \(S\) 中维护已经生成的凸包(hull-so-far)。起初栈中包含前两个点，即 \(S = (b, a)\)，其中 \(b\) 位于栈顶。我们将沿用栈从上到下、从左到右的排列方式， 我们将在 3.5.5 节中实现它。添加点 \(c\) 是因为 \((a, b, c)\) 在栈顶 \(b\) 处有一个左转的形态。\(S = (c, b, a)\) 是一个凸链，我们将一致维持这个条件。 再来看点 \(d\)，由于 \((b, c, d)\) 在栈顶 \(c\) 处右转，不构成凸链，所以做出撤销添加 \(c\) 的决定，将它从栈顶弹出，使 \(S\) 再次变为 \((b, a)\)。 现在添加点 \(d\)，因为 \((a, b, d)\) 在栈顶 \(b\) 处左转。

继续以这种方式添加 \(e\) 和 \(f\)之后，栈变为 \(S = (f, e, d, b, a)\)。由于 \((e, f, g)\) 和 \((d, e, g)\) 都是右转，因此点 \(g\) 会导致 \(f,e\)被删除。 然后可以添加 \(g\)，栈变为 \(S = (g, d, b, a)\)。以此类推。

如果我们像示例那样幸运，第一个点 \(a\) 就在凸包上，那么凸链将自然闭合，最终得到凸包 \(S = (j, i, g, d, b, a)\)。 注意，从栈顶到栈底代表顺时针遍历，因为我们是通过逆时针扫描构建的。如果 \(a\) 不在凸包上，那么链的头部就会开始“吞噬”尾部，算法分析也会变得更加困难。 我们将看到，这种情况是可以避免的。

3.5.2 版本 A 的伪代码(Pseudocode, Version A)

在更详细的阐述算法之前，我们先用算法 3.5 中用伪代码粗略的了解一下它的全貌。假设栈结构有原语 Push(p, S) 和 Pop(S)， 分别用于将元素 \(p\) 入栈 \(S\) 和出栈。我们用 \(t\) 表示栈顶，用 \(i\) 索引按角度排序的点。还有许多问题需进一步探讨，特别是程序的开始和结束， 在这个伪代码里，有一个 while 循环迭代 \(O(n)\) 次，每次出栈都会永久移除一个点，因此迭代次数不能超过 \(n\)。 加上 \(n\) 次前向操作，循环最多迭代 \(2n\) 次。因此，在排序步骤，算法的运行时间是线性的，而排序步骤的时间复杂度为 \(O(n \log n)\)。 我们将在 3.6 节中看到，这已经最佳结果了，它的时间复杂度是"最坏情况下的最优(worst-case optimal)"。

3.5.3 细节:边界条件(Details: Boundary Conditions)

我们还忽略了一些细节，将分两个阶段来处理他们。首先，本节将探讨各种边界条件。其次，后续章节将探讨实现方面的问题。

循环的开始和结束(Start and Stop of Loop)

虽然 Version A 的循环很简单，但也很难正确地启动和停止。该算法在起始和结束都会遇到一些麻烦。我们已经提过，如果栈底的元素 \(a\) 不在凸包上，这个循环就很难终止。 当第二个入栈的元素 \(b\) 不在凸包上，循环就很难启动。假设 \((a,b,c)\) 是一个右转结构，那么 \(b\) 将被弹出，此时栈变成 \(S = (a)\)。 但是至少需要两个元素才能确定第三个元素是否与栈顶元素形成左转结构。

显然，如果 \(a\) 和 \(b\) 都在凸包上，那么就可以避免无法启动或者终止的问题。我们将在下一小节中介绍如何实现这点。

排序的起点(Sorting Origin)

我们假设点 \(x\) 位于凸包的内部，并且其它点围绕着 \(x\) 排序。Graham 设计了一个精巧的线性算法来计算这样的内部点。 然而，这种计算并不必要，而且当输入坐标都是整数的情况下还会强制使用浮点数。为了保证整数输入下的正确结果，我们希望避免所有浮点运算。 译者按:目前我还不太能理解作者为啥要规避浮点运算，可能那个年代浮点运算太奢侈了吧。

简化方法是根据集合中的某个点进行排序，特别是凸包上的某个点。我们将使用最低点，它显然位于凸包上。 如果有多个最低点都有相同的 \(y\) 坐标，我们将使用最右侧的作为排序原点，对于图 3.4 该点是 \(j\)。 排序后的结果如图 3.5 所示，图中的所有点都已重新标记了数字，这是它们在实现中的索引。我们将这些点称为 \(p_0, p_1, \cdots, p_{n-1}\)， 其中 \(p_0\) 为排序原点 \(p_{n-1}\) 为逆时针排序的最后一点。

现在我们来解决前面讨论的启动和终止问题。如果我们从原点 \(p_0\) 出发做一水平射线，按照逆时针方向对各点进行排序，那么 \(p_1\) 必定位于凸包上， 因为它与 \(p_0\) 形成了一个极角(extreme angle)。但是它可能不是极点(如 3.1.1 节所定义)， 这个问题我们将在后面讨论。如果我们将栈初始化为 \(S = (p_0, p_1)\)，那么它将始终至少包含两个点，这避免了启动困难。 而且当循环回到 \(p_0\) 时，栈中的点永远不会被消耗，从而避免了终止困难。

共线(Collinearities)

现在我们来考虑最后一个边界条件，即三个或更多点共线的问题。此前我们一直假设这种情况不会发生。实际上它会在多个方面影响算法。首先，我们先明确一下期望的输出结果。

凸包共线(Hull Collinearities)。这里我们采用最有用的输出(4):仅包含按照凸包顺序排列的极点(the extreme vertices only, orderd around the hull)。因此，如果输入包含一个正方形的四个角以及散布在其边界上的点，那么输出应仅包含正方形的四个角。规避凸包上的非极点很容易实现， 只需要求严格左翻转 \(p_{t-1}, p_t, p_i\) 将 \(p_i\) 压入栈中，其中 \(p_t\) 和 \(p_{t-1}\) 是栈顶的两个点。此外如果 \(p_t\) 与 \(p_{t-1}\) 和 \(p_i\) 共线，那么将其删除。

共线排序(Sorting Collinearities)。共线引出了另一个问题，如果点 \(a\) 和 \(b\) 与点 \(p_0\) 构成的线段与水平线的夹角相同，我们应该如何对其排序呢？ 通常我们会假设(或希望)这无关紧要，但实际情况要复杂得多。至少有两项工作。首先，使用一致的排序规则，然后再妥善处理起始点和终止点以及凸包的共线问题。 一个合理的规则是，如果 \(\text{angle}(a) = \text{angle}(b)\)，那么距离点 \(p_0\) 更近的点应当排在前面，即若 \(|a - p_0| < |b - p_0|\) 则 \(a < b\)。 使用此规则，我们可以得到图 3.6 中所示的排序结果。

注意，图中点 \(p_1\) 点并非极点，因此从 \(S = (p_1, p_0)\) 开始是会有问题的。虽然可以通过从 \(S = (p_{n-1}, p_0)\) 开始来规避这个问题(注意 \(p_{n-1} = p_{18}\) 是极点)， 但在此我们要换一种方法。我们观察后发现，如果 \(\text{angle}(a) = \text{angle}(b)\) 且根据上述规则排序 \(a < b\)，若 \(a\) 不是凸包的极点，则可以删除它。 在图 3.6 中，点 \(p_1, p_5, p_9, p_{12}, p_{17}\) 都将因此而被删除。

重合点。通常情况下，如果点集中有一些相同点时，不考虑他们的话，程序可能会崩。我们可以将此问题视为排序共线性的一种特殊情况，即删掉那些重复点，只保留其中一个。

3.5.4 版本 B 的伪代码(Pseudocode, Version B)

在讨论算法实现细节之前，我们用算法 3.6 中的伪代码总结前面的讨论，其中包含了如下的更改。

我们尚未讨论循环终止的细节。考虑共线的问题，应该用条件 \(i < n\) 终止? 或者应该是 \(i < n\)，还是 \(i < n - 1\)? 从伪代码中可以看出，进入 while 循环时，共线的问题已经被处理了。因此，一旦 \(p_{n-1}\) 被压入堆栈，我们就可以认为循环结束了。 图 3.6 中的 \(p_{17}\) 已经删掉了。因此，循环终止条件应该是 \(i < n\)。

3.5.5 Graham 算法的实现(Implementation of Graham's Algorithm)

现在我们来讨论算法 3.6 的实现。假设输入点的坐标为整数，并且我们坚持避免所有浮点运算，以确保输出的正确性。 实际上，只有当坐标不太大时，才能保证这一点，但撇开这个限制不谈，该实现能够得到正确的凸包。

我们从数据结构开始，解决排序问题，最终给出实现代码。

数据的表示(Data Representation)

通常，我们可以选择将点存储在数组或链表中。这里我们选择使用数组，因为我们需要排序，所以要把数据存储在连续的内存空间中。 我们将使用一个结构体来表示每个点，该结构体与第一章中描述顶点的结构体(Code 1.2)类似。 这些点存储在一个全局数组 \(P\) 中，其中 \(P[0], \dots ,P[n-1]\) 分别对应 \(p_{0},\dots ,p_{n-1}\)。每个 \(P[i]\) 都是一个结构体， 包含坐标字段、唯一的标识编号以及一个用于标记删除的标志位。参见Code3.1。

        typedef struct tPointStructure tsPoint;
        typedef tsPoint *tPoint;
        struct tPointStructure {
            int vnum;
            tPointi v;
            bool delete;
        };
        #define PMAX 1000 /* Max # of points */
        typedef tsPoint tPointArray[PMAX];
        static tPointArray P;
        int n = 0;        /* Actual # of points */

栈最自然的表示方法是使用一个单链表，链表中的每个单元格"包含"一个点，即包含指向该点记录的指针。参见Code3.2。根据这些定义， 栈顶可以声明为 tStack top，栈顶下方的元素为 top->next。

        typedef struct tStackCell tsStack;  /* Used only in NEW() */
        typedef tsStack *tStack;
        struct tStackCell {
            tPoint p;
            tStack next;
        };

我们需要三个栈接口，参见Code3.3: Pop,Push,PrintStack。栈顶始终是列表的头部(最左侧)元素。Pop 函数用于释放栈顶并返回指向下一个单元的指针。 Push(p, t) 函数分配新的存储空间，将其填充为 p，并设为新的栈顶。注意在 PrintStack 函数中，t->p->v 指向点的坐标，t 的类型为 tStack， p 的类型为 tPoint，v 的类型为 tPointi，后者是Code1.1 中用于表示点坐标的数据结构。

        tStack Pop(tStack s) {
            tStack top;
            top = s->next;
            FREE(s);
            return top;
        }
        tStack Push(tPoint p, tStack top) {
            tStack s;
            /* Get new cell and fill it with point. */
            NEW(s, tsStack);
            s->p = p;
            s->next = top;
            return s;
        }
        void PrintStack(tStack t) {
            while (t) {
                printf ("vnum=%d\tx=%d\ty=%d\n", t->p->vnum,t->p->v[X] ,t->p->v[Y]);
                t = t->next;
            }
        }

排序(Sorting)

查找最低点(FindLowest)。我们先来处理排序中最简单的部分，找到集合中最右侧的最低点。函数 FindLowest(Code3.4) 可以完成此操作，并将该点交换到 P[0] 中。 函数 Swap 太简单了此处省略。

        void FindLowest(void) {
            int i;
            int m = 0; /* index of lowest so far. */

            for (i = 1; i < n; i++)
                if ((P[i].v[Y] < P[m].v[Y]) || ((P[i].v[Y] == P[m].v[Y]) && (P[i].v[X] > P[m].v[X])))
                    m = i;
            Swap(0, m); /* Swap P[0] and P[m] */
        }

规避浮点运算(Avoiding Floats)。排序步骤看似简单，但要保证排序的准确性，还是有些坑的。首先，我们引入一些符号。 令 \(r_i = p_i - p_0\)，表示从 \(p_0\) 到 \(p_i\) 的向量。我们要给出一个精确的计算方法，以确定何时 \(p_i < p_j\)，其中 "<" 表示排序关系。

Atan2。最直接的选择是定义 \(p_i < p_j\)，如果 \(\text{angle}(r_i) < \text{angle}(r_j)\)，其中 \(\text{angle}(r)\) 是 \(r\) 相对于正 \(x\) 轴逆时针转动的角度。 参见图 3.7。我们稍后会讨论打破平局的规则。由于 \(p_0\) 是最低点，所有这些角度都在 \((0, \pi]\) 范围内，这会很方便，因为角度有正有负的时候，排序会很麻烦 (我们将在 8.4.4 节中看到)。C语言恰好提供了标准函数, \(\text{angle}(r) = \text{atan2}(r[Y], r[X])。 但并不建议使用此函数，至少有两个原因。那个时候不推荐使用，现在还不推荐?

斜率Slopes。对于第一象限，可以用斜率 \(r[Y] / r[X]\) 代替反正切函数。在第二象限，有 \(-r[X] / r[Y]\)。虽然这种方法比 atan2 的计算更简单，但它也存在一些类似缺陷。 比如，如果 \(r_i = c r_j\)，其中 \(c\) 为某个正数，那么我们期望得到 \(\text{angle}(r_i) = \text{angle}(r_j)\)，然后执行我们的平局判定代码。 这相当于 \(a / b = (ca) / (cb)\)，在数学上当然是正确的。但是，这个等式在 C 语言的浮点除法中并不能保证成立，它与机器实现有关。 特别是，在 Cray Y-MP/4 上，\(1.0/1.0 == 3.0/3.0\)的结果为 FALSE。这是 Cray 执行除法时会先进行倒数运算，然后再进行乘法运算，而这两个运算有时会导致微小的误差。 例如，\(3.0/3.0 == 0.999999999999996\)。

教训就是，即使某些机器的 C 编译器和库完全符合规范，任何浮点运算都可能出错。因此，我们尽可能选择整数运算。

左Left。有些读者可能早就意识到，其实可以用 Left 函数来判断一个点是否位于由另外两个点确定的直线的左侧，这完全能满足我们对比 \(r_i, r_j\) 的需求。 我们在第一章中完全使用整数运算实现了这一点。 我们实现的函数 Left 本身是对 Area2 输出的简单测试，Area2 计算由三个点确定的三角形的有符号面积。这里我们将使用这个面积函数而不是 Left，因为这样更容易区分并列的情况。

在 1.4.3 节中实现的 Area2 是有缺陷的。当输入坐标相对于机器字长较大时，计算可能会溢出。 然而，标准 C 语言在整数溢出时不会报错，这一点我们将在第四章(Code4.23)中再次讨论。 注意，现在大多数执行几何计算的 C 代码都面临着一个两难困境，如果代码使用浮点运算，则可能依赖于机器相关的特性。而如果仅使用整数运算，则除非严格限制整数的大小，否则无法保证其正确运行。 幸运的是，随着 LEDA 等鲁棒的计算软件包的出现(参见第 9 章)，这种情况正在开始改变。

接下来，我们先讨论排序算法，再展示如何使用 Area2 进行比较。

快排Qsort。标准 C 库中包含一个名为 qsort（“快速排序”）的排序函数，它是一个通用函数，需要费点事配置一下才能正确使用。 该函数需要输入数据的形状和位置信息，并提供一个比较函数。它会使用提供的比较函数对数据从小到大进行原地排序。有四个参数:

  1. 指向数据起始地址的指针，在本例中为 &P[1]，即第一个点的地址。(在我们的排序中 \(P[0]\) 是固定的。)
  2. 参与排序的元素数量，\(n-1\)
  3. 每个元素占用字节数量，sizeof(tsPoint)
  4. 二元比较函数的指针，Compare

比较函数Compare。qsort 要求比较函数能够，根据第一个参数与第二个参数的关系，返回一个小于、等于或大于 0 的整数，来反映第一个参数小于、等于或大于第二个参数。 此外，调用该函数时，需要传入指向待比较元素的指针。这带来了一个小小的技术难题，因为我们需要三个输入，\(p_0, p_i, p_j\)。 一种解决方法是计算一个关于 \(r_i\) 的数组出来。这里我们选择另一种方案，将 \(P\) 设为全局变量，这样就可以在 Compare 函数体内部引用 \(p_0\) 而无需将其作为参数传递。

我们终于可以写出 Compare 函数的内核了，如 Code 3.5所示。如果 \(p_j\) 位于有向直线 \(p_0p_i\) 的左侧，则 \(p_i < p_j\)。 换句话说，如果 \(\text{Area2}(p_0, p_i, p_j) > 0\)，则 \(p_i < p_j\)，此时，qsort 将返回 -1 以表示小于。参见图 3.7。

        int Compare(const void *tpi, const void *tpj) {
            int a;    /* area */
            int x, y; /* projs, of ri & rj in 1st quadrant */
            tPoint pi, pj;
            pi = (tPoint)tpi;
            pj = (tPoint)tpj;
        
            a = Area2(P[0].v, pi->v, pj->v);
            if (a > 0)
                return -1;
            else if (a < 0)
                return 1;
            else {    /* Collinear with P[0] */
                x = abs(pi->v[X] - P[0].v[X]) - abs(pj->v[X] - P[0].v[X]);
                y = abs(pi->v[Y] - P[0].v[Y]) - abs(pj->v[Y] - P[0].v[Y]);
                if ((x < 0) || (y < 0)) {
                    pi->delete = TRUE;
                    return -1;
                } else if ((x > 0) || (y > 0)) {
                    pj->delete = TRUE;
                    return 1;
                } else { /* pints are coincident */
                    if (pi->vnum > pj->vnum)
                        pj->delete = TRUE;
                    else
                        pi->delete = TRUE;
                    return 0;
                }
                ndelete++; // 这里有 bug 吧，这句应该执行不到的吧
            }
        }

当面积为零时，我们会遇到共线性问题，这需要按照之前讨论的方式进行处理(3.5.3 节)。首先，我们需要确定哪个点更接近 \(p_0\)。 如果 \(p_i\) 更近，那么向量 \(r_i = p_i - p_0\) 的投影长度比 \(r_j = p_j - p_0\) 的投影长度更短。我们可以利用这点来避免计算距离。 Compare 会计算它们在 \(x,y\) 方向上的投影，并进一步比较。如果 \(p_i\) 更接近，则标记它稍后删除，如果 \(p_j\) 更近，则标记\(p_j\)。 如果 \(x = y = 0\)，那么有相同点 \(p_i = p_j\)，我们标记索引较小的那个。(选择删除哪个并不重要，但 qsort 要求保持一致。) 在所有情况下，该函数都会返回 \(\{-1, 0, +1\}\) 中的一个。

将指针 const void *tpi 强制转换成 pi = (tPoint)tpi 满足了 qsort 的要求， 同时可以在函数 Compare 自然的使用。

主函数(main)

解决了排序细节之后，我们从顶层设计上讨论一下主函数(Code3.6)。读取点集后，通过函数 FindLowest 扫描到最右侧的最低点将其交换到 P[0]， 然后使用 qsort 对 \(P[1], \dots, P[n-1]\) 排序。函数 Compare 会通过 delete 字段来标记一些待删除的点，我们在一个简单的函数 Squash 中完成删除操作(Code3.7)。 该函数维护了两个索引，\(i\) 是查询索引，依次检查各点的 delete 字段，\(j\) 是非删除点索引。函数将所有非删除的点 \(P[i]\) 复制到 \(P[j]\) 上。 现在最棘手的情况已经解决了，我们可以调用 top = Graham() 继续进行 Graham 扫描。

        main() {
            tStack top;
        
            n = ReadPoints();
            FindLowest();
            qsort(&P[1], n-1, sizeof(tsPoint), Compare);
            Squash();
        
            top = Graham();
            PrintStack(top);
        }

        void Squash(void) {
            int i,j;
            i = 0; j = 0;
            while (i < n) {
                if (!P[i].delete) {
                    Copy(i, j);
                    j++;
                }
                /* else do nothing: delete by skipping. */
                i++;
            }
            n = j;
        }

Graham 扫描的代码(Code for the Graham Scan)

Graham 扫描的代码几乎是前面伪代码(Code3.6)中 while 循环的直接翻译，参见Code3.8。这里我们遇到了一个非常常见的情况，大部分编码都是在做一些周边的准备工作， 而几何算法的核心部分只需要一小段代码就能实现。

        tStack Graham() {
            tStack top;
            int i;
            tPoint p1, p2;  /* Top tow pionts on stack */
            /* Initialize stack */
            top = NULL;
            top = Push(&P[0], top);
            top = Push(&P[1], top);
            /* Bottom two elements will never be removed. */
            i = 2;
            while (i < n) {
                p1 = top->next->p;
                p2 = top->p;
                if (Left(p1->v, p2->v, P[i].v)) {
                    top = Push(&P[i], top);
                    i++;
                } else {
                    top = Pop(top);
                }
            }
            return top;
        }

示例(Example)

如图 3.8 所示(与图 3.6 中的点相同)，点集中包含了所有类型的共线问题，旨在对代码进行严格测试。表 3.1 中列出了，角度排序后的点，以及需要删除的标记。 表中所示的 vnum 索引与图 3.8 中的标签一致。调用函数 Squash 删除五个点后，还剩 \(n = 14\) 个点。

右侧是随着 i 的增长，栈结构的变化过程。

在这个例子中，总迭代次数为 \(19 < 2 \times n = 2 \times 14 = 28\)。最终结果为 \((p_{11}, p_9, p_8, p_7, p_{16}, p_6, p_{18}, p_{15})\)，即顺时针方向的极点。

3.5.6 总结(Conclusion)

我们得到一个定理

定理3.5.1.平面上 \(n\) 个点的凸包可以通过 Graham 算法在 \(O(n \log n)\) 时间内找到。排序需要 \(O(n \log n)\) 次迭代， 扫描最多需要 \(2n\) 次迭代。他的算法可以完全用整数运算实现。

3.5.7 练习(Exercises)

  1. 所有点都共线All points collinear。如果所有输入点都共线，代码会输出什么？
  2. 最佳情况Best case。如果输入点全部都在凸包上，那么扫描的 while 循环(算法 3.6)会执行多少次迭代？
  3. 最坏情况Worst case。为每个 \(n\) 构造一组点，使得扫描的 while 循环迭代次数最大(在算法 3.6 中)。
  4. [Programming] 随机凸包Random hulls。编写代码生成正方形内的随机点。修改 graham.c，使其仅输出读取的点数和凸包上的点数。 进行足够多的随机试验，并使用足够大的 \(n\) 值，以验证定理，期望值是 \(O(\log n)\)。



Copyright @ 高乙超. All Rights Reserved. 京ICP备16033081号-1