平面上的凸包计算

    陈道蓄

    

    

    

    在郊野公园中有一片林地，生长着一些古老的树木。管理部门希望建围栏把这些树木围起来加以保护。为了便于外围修建步行道，方便游人观赏，将保护区设计成凸多边形。当然也希望围栏总长度尽可能小，以降低建设成本。

    为简化计算，我们假设可以用部分树木作为围栏的桩柱，换句话说，部分树木处于保护区域的边界上。图1是这个问题的示意图，左边标出树木的平面位置分布，右边则显示完成的围栏，位于围栏上的树木用空心点表示，围在内部的为灰色点。我们能否设计一个算法让计算机帮我们确定该如何建围栏呢？

    ● 问题模型和基本思想

    问题可以抽象为：在平面上给定一组点，如何生成一个包含全部点的最小凸包。所谓凸包是满足如下条件的多边形：连接多边形中任意两点（包括边界上的点）的直线段完全位于多边形内部（包括边界）。“最小”是指在保持凸多边形性质不变的前提下，若缩小该图形，则给定的点中至少有一个位于外部。

    图形处理中大量的问题需要用到这样的计算，凸包算法是计算几何领域最早的成果之一。

    模型基于平面笛卡尔坐标系。输入是一组点v1，v2，…，vn，其中vi表示为（xi，yi），即该点在笛卡尔坐标系中的横坐标值与纵坐标值。我们不妨假设所有坐标值均为整数，并且约定任何两点不会“过于靠近”。这个要求并不明确，由于计算过程中可能涉及小数，如果点过于接近可能导致后面讨论的算法产生错误结果。我们避免烦琐的详细定义，这并不影响读者理解算法的思想。

    算法的输出是输入点的一个子集构成的序列。出现在序列中的点即位于凸多边形边界上的点。边界的轨迹即为需要计算的凸包。我们约定其顺序是：从某个指定点出发，按照顺时针方向沿凸包排列。如果用画图工具将输出序列中每个相邻点对之间的连接直线段（包括首尾两个点之间的）全部画出，则可显示计算得到的凸包。

    显然，解题的关键是确定哪些点应该在凸包上。当点数非常多，且随机分布时，这很难确定，但作为构成凸包的每条线段，其特征倒不难看出来。从图1右边的图很容易看出边界上的每条线段所在的直线将平面切分为两个半平面，所有的点一定位于其中一个半平面上（包括分界线），如图2所示。这就意味着满足这一条件的线段的两个端点在边界上。反之，考虑任意两点连接的直线段，如果其对应的直线分割成的两个半平面中都有输入点，则该线段不可能在凸包上。

    ● 从思想到算法

    利用解析几何知识，很容易判定相对于特定线段，任给的两点是否在同一个半平面（是否位于线段的同一侧）。

    假设线段uv所在的直线将平面分割为两个半平面（如图3），点s和t处在同一个半平面中，折线s-u-v和t-u-v在u点总是向同一方向偏转（这里是右转），而处于另一半平面中的点w决定的折线w-u-v一定是向相反方向偏转（这里是左转）。

    以图3中的折线t-u-v为例，设三个点的坐标值分别为（x1，y1），（x2，y2），（x3，y3），则折转方向完全由下面的行列式值的符号所确定：

    当行列式值为正时，则反时针（左）转，当行列式值为负时，则顺时针（右）转。行列式值为0当且仅当三点共线。这里的例子该行列式值为负。两个点相对于线段uv处于两个不同的半平面当且仅当相应行列式的值符号相反。我们不详细讨论数学背景，读者可参考文献【1】。

    利用上述公式，我们很容易实现如下页图4所示的line_turn函数。

    利用函数line_turn，我们先给出一个非常直觀的算法：对输入的任意两点连接得到的直线段，判断其他n-2个点是否全部在其同一侧，如若不然，则这条线段不可能是凸包一部分（如下页图5）。简而言之，采用穷尽法找出凸包上的线段。（注意：这里输出的是线段集合，而不是模型描述中提到的边界点序列）

    输入n个点，可生成的直线段数量为O（n2），对每条线段，最坏情况下要检查除该线段端点以外的n-2个点是否都在同一个半平面中，因此计算复杂度是O（n3）。

    上述穷尽法没有利用输入中可能有帮助的隐含信息，一般效率不高。如果能找到一些“窍门”，直接可判定某些线段是否在凸包上，而不用每次都去检查所有点，就有可能改进算法。

    ● 一个效率较高的贪心算法

    我们再仔细审视一下图1中右边那个图，特别观察下面标注出的点（如下页图6）。

    尽管图6中的凸包在算法未执行完成前并不存在，但图中特别标示出的三个点，即v0（靠坐标值即可确定）、v1和vn-1是可以确定的。后两个点的确定只需要对除v0外任意点与v0及X轴正方向之间的夹角大小排序。最有启发意义的一点在于，我们可以断定输入的所有点都位于这三点构成的折线的同一侧。

    确定了v0后，首先对除了v0以外所有点相应的夹角从大到小排序（且按此给v0以外的点编号），显然在最终计算出的凸包上，按照顺时针方向，点的序号严格递增。而且沿凸包轨迹，任何正方向上相邻的三个点构成右转折线。

    这为设计效率更高的凸包算法提供了更清晰的思路：首先确定v0，这很容易。接着就按照上述相应的夹角从大到小给其他n-1个点排序，如果夹角相等就按与v0的距离从小到大排（这是为了处理共线的情况）。选择v0和v1作为开始的两个边界点（按照顺时针方向）。后面循点序号的增序逐步增加凸包上的点。

    按照角度大小排序并不能保证任意连续三个序号的点一定构成向右的折线。这是算法中最不直观的部分。每当加入一个新的边界点，必须检查当前序列中最后三个点构成的折线是否左转，如果是，就删除当中的点，再继续回溯并做同样的检查。下页图7是选初始点以及按照夹角大小排序的结果。

    下页图8显示从初识点开始，沿顺时针方向构造凸包的过程。

    下页图8（a）为起始状态。为了计算方便，选定v0后将其作为坐标系原点（0，0），并根据解析几何知识将其他所有点的坐标值做相应修改（这只需要简单算术运算，总代价是O（n））。

    图8（b）显示计算过程。每条边上数字表示在整个操作序列中相应线段加入以及被删除的次序（删除操作次序表示在括号中，边界上的线段没有删除操作）。需要指出的是，算法并不对边（线段）操作，只处理点。这里是为了让读者容易跟踪算法执行过程，画出线段，其实其加入还是删除是对点操作的反映。有两条线段上删除操作执行次序相同，那是因为算法其实删除的是点，这就同时把与该点关联的两条边删除了。

    其实并不需要计算出夹角的大小，我们只是用这个数据来排序。在0～π范围内正切函数是在两个不连续的区间（0～π/2，π/2～π）内分别严格递增的，所以只要区分横坐标的正负，直接比较y/x就可以了（区分正负区时必须保证不能让正负号不同的坐标值放在一起比）。如果想避免浮点计算带来的不精确，可以利用有关的库函数进行分数值计算与比较，当然也可以自己写一段程序实现精确的分数计算，这并不难。

    从上面的例子可以看出，算法执行中是通过回溯检查是否有不正确的点被当作边界点，需要不断纠正前面的误判。显然这个过程用堆栈实现最为方便。

    定义栈CH，按照点的序号，首先是v0，v1，v2依次进栈，从v3开始，每次按照序号向栈中加一个点，随即检查栈顶的连续三个点是否构成左偏转的折线，如果是，则从栈中删除中间那个点，并继续检查新栈顶位置的连续三个点，以此类推。折转方向的判定直接调用函数line_turn即可。（读者可以考虑为什么v2进栈时不需要检查）算法过程描述如下页图9所示。

    这个算法的最主要代价就是排序（想不到吧！），除此之外對每个点的处理代价都是常量，所以总代价为O（n）。整个算法的复杂度是O（nlogn）（也就是排序的代价）。

    ● 结束语

    虽然从渐进复杂度上来说，上述算法已经达到“最优”了，但仍然可以设法改进。

    一个非常简单的想法：先找出四个“极点”，即横/纵坐标值最大和最小的点。用这四个点构成一个四边形，位于该四边形中（包括内部与边界）的点除了“极点”外，都不可能是边界点。图10中标出了上述例子中的四个“极点”，并按照上面的方法进行“预处理”，算法需要考虑的点数从16个下降到10个。利用解析几何的基础知识很容易识别位于四边形中的点。这是个“启发式”算法，对于输入点集的任意分布，我们无法分析预处理“收益”可能有多大，但经验数据表明，当输入点数很大，且随机分布时，计算代价会大大下降。

    计算角度（哪怕是利用三角知识间接计算）仍然比较麻烦。文献【2】中介绍了只需要对点的横坐标排序的算法。文献【2】与文献【1】互补性很好。后者对相关数学背景以及算法的分析有很好的讨论，而前者对算法的实现和比较讨论得非常详细，包括给出了程序代码。

    这里的算法利用了与平面紧密相关的特性，如两条直线夹角等，所以不能将这样的思想自然推广到更高维度。其实凸包问题除了与几何空间相关，在大数据分析等应用中也有很大价值。那里的“维数”与物理空间无关，往往是允许定义的对象属性的个数，所以“高维”是很常见的，但三维或更高维的凸包计算问题超出了本文的范围。

    参考文献：

    [1]Thomas Cormen， etc. Introduction to Algorithms[M].殷建平，等，编译.北京：机械工业出版社，2013.

    [2]George Heineman， etc.Algorithms-in a Nutshell， 2nd ed. OReilly 2016（算法技术手册，影印版第2版）[M].南京：东南大学出版社，2017.

    注：作者系南京大学软件学院原院长，计算机系原系主任。