很多计算几何问题都可以利用本章介绍的相关代码组合解决,不过,只解决固定模式的问题并不是几何问题的全部内容。很多情况下需要组合这些代码,以编写出更复杂的算法。此时,几何算法的设计范式将有助于编写出更高效的算法。本节将介绍两种比较著名的范式,并通过示例讲解如何利用这些范式提高解题效率。 平面扫描 平面扫描(plane sweeping)或线扫描(line sweeping)范式是设计计算几何算法的最著名的范式之一。利用了平面扫描的算法会沿着水平或垂直方向“扫过”给定平面以解决问题。虽然看似很简单,但能够利用这种方法解决的问题特别多。 矩形合集的面积 给出几个各边均与x轴和y轴相互平行的矩形时,求其合集面积。如果只是求出各矩形面积后相加,就会重复计算重叠部分的面积,所以此题并没有那么简单。 解决此题的一种方法就是,利用求合集大小时常使用的容斥原理(inclusion-exclusion principle)。首先,求出各矩形面积后相加,然后求出每一对矩形的交集。这些交集被相加了两次,所以从结果中减掉1次。此时,每3个矩形的交集会被减掉3次。而第一步相加过程中,此部分也只被相加了3次。因此,还需要给结果加上这些部分。反复执行这些操作,最终就能得到答案。此方法最大的问题在于,需要检索所有能够组合出的交集。因此,有n个矩形时,需要访问2n个交集。如果题目中给出10到20个四边形,还尚有可能完成运算。但超过这些数目后,运算速度会变得非常缓慢。 设计平面扫描算法前,先想象一条从左向右扫过的垂线。假设有函数cutLength(),给出垂线坐标x时,能够返回用垂线分割四边形后的断面 的长度。那么,期望得到的合集面积就是此函数对x进行定积分的结果。 当然,并不需要实际对cutLength()函数进行积分。请看图15-17(a)中的4条垂线。各位置上的断面长度全都相同,只有矩形左右边的断面长度才有变化。我们将想要的值开始变化的位置称为事件(event)。图15-17(b)表示相同图标中的事件位置。预先算出各事件的位置并排序后,各事件之间的cutLength()函数值不变。因此,在某个事件的位置计算出cutLength()函数后,乘以到下一个事件的距离,结果就是这两个事件之间合集的面积。图15-17(c)表示利用这种方法解题时需要算出面积的矩形。 代码15-15是实现这种方法的unionArea()函数。代码并没有为了求出某个事件位置上的cutLength而遍历所有矩形,而是保留了表示垂线哪个部分与四边形重叠的数组count。开始时把所有四边形的y1、y2集合并排序后,就能利用这些坐标值分割垂线。图15-17(d)表示分割垂线的横线和1个区间内的count值。用count值更新各事件后,就能轻松求出cutLength。 图15-17 利用平面扫描求矩形的合集面积 代码15-15 求矩形合集面积的unionArea()函数 // x1<x2, y1<y2 struct Rectangle { int x1, y1, x2, y2; }; // 计算矩形面积。 int unionArea(const vector<Rectangle>& rects) { if(rects.empty()) return 0; // 事件信息:(x坐标,左侧还是右侧,四边形的序号) typedef pair<int,pair<int,int> > Event; vector<Event> events; vector<int> ys; // 遍历各四边形并找出y坐标和事件的集合。 for(int i = 0; i < rects.size(); ++i) { ys.push_back(rects[i].y1); ys.push_back(rects[i].y2); events.push_back(Event(rects[i].x1, make_pair(1, i))); events.push_back(Event(rects[i].x2, make_pair(-1, i))); } // 对y坐标集合排序并去重。 sort(ys.begin(), ys.end()); ys.erase(unique(ys.begin(), ys.end()), ys.end()); // 对事件目录排序。 sort(events.begin(), events.end()); int ret = 0; // count[i]=ys[i]~ys[i+1]区间中重叠的四边形个数。 vector<int> count(ys.size()-1, 0); for(int i = 0; i < events.size(); ++i) { int x = events[i].first, delta = events[i].second.first; int rectangle = events[i].second.second; // 更新count[] int y1 = rects[rectangle].y1, y2 = rects[rectangle].y2; for(int j = 0; j < ys.size(); ++j) if(y1 <= ys[j] && ys[j] < y2) count[j] += delta; // 计算cutLength int cutLength = 0; for(int j = 0; j < ys.size()-1; ++j) if(count[j] > 0) cutLength += ys[j+1] - ys[j]; // 将cutLength乘以到下一个事件的距离的结果与ret相加。 if(i + 1 < events.size()) ret += cutLength * (events[i+1].first - x); } return ret; } event和ys的元素个数都是O(N),所以此算法的时间复杂度为O(N2)。如果不用数组,而使用灵活的数据结构(第24章树状数组),就能把时间复杂度减少到O(NlgN)。 求多边形交集的面积 用特定轴线分割输入的图形会更容易求解。这种想法可应用于很多问题的解题过程。本章的金银岛问题也可以利用平面扫描法解决。将两个多边形的顶点和各边交点的x坐标视为事件,那么就如图15-18所示,两个多边形的交集中,两个事件之间的形态呈梯形。移动垂线并求出所有梯形面积,相加后就能得到交集面积。此方法虽然没有直接求出交集,但也可以应用于给出两个凹多边形时求交集面积的问题。 图15-18 多边形的交集是梯形的合集 相交线段 下面,给出平面线段的集合时,确认这些线段中有无相交线段 。解决此问题的最通用的方法是,逐对确认每对线段是否相交。虽然这种算法的时间复杂度是O(N2),但很难想出更快的算法。不过令人惊讶的是,用平面扫描算法就能在O(NlgN)时间内解决此问题。 试想一条从最左侧开始扫过整个平面后到达最右侧的垂线,那么就能找出这条垂线与线段左端点和右端点相遇的事件。图15-19表示平面上的各条线段,以及将此平面用事件分割后的结果。星号表示两条线段的交点,由图可知,在星号前一事件中相交的两条线段之间,没有任何其他线段。沙莫斯-霍伊(Shamos-Hoey)算法正是利用了这种性质。 图15-19 判断平面线段是否相交的沙莫斯-霍伊算法的操作原理 从左端起始访问各个事件,将当前与垂线相交的线段添加到集合,并保持更新。更新的方法是,若垂线与线段的左端点相遇,就将线段添加到集合;若与线段的右端点相遇,就将线段从集合中删除。并且,集合始终以各线段与垂线相遇时的y坐标的升序排列。此时,按照下列方法判断是否相交。 每次向集合添加线段时,就确认集合上与之相邻的两条线段是否相交。 每次从集合删除线段时,就确认集合上与之相邻的两条线段是否相交。 算法确认出有1对线段相交时,就会马上终止运行。此算法的正确性并不难证明,不过,线段集合的排序不容忽视。线段与垂线接触点的y坐标会不断根据垂线的移动发生变化,即使这样,上述方法并没有对集合重新排序。这是为什么呢?因为此算法如果遇到1对相交线段就会终止运算,而还未出现相交时,线段的相对顺序总是保持不变。 若想让此算法在O(NlgN)时间内完成运算,就必须选择高效的数据结构。事件个数已经是2n,所以向集合添加和从集合中删除、检索插入位置等运算必须在O(lgN)时间内完成。第22章的二叉搜索树就能满足这种性质。垂线位置总在变化,所以每次比较两条线段时,需要重新计算线段位置。因此,虽然利用标准库中的二叉搜索比较麻烦,但并不难。 如果利用此算法的扩展版本Bentley-Ottmann算法,不仅能够确认相交线段的存在,还能够在O((N+K)lgN)时间内找出所有交点(K为交点个数)。 旋转卡尺 另一个比较著名的计算几何算法的设计范式就是旋转卡尺(Rotating Calipers)。卡尺是测量较小物体的直径、宽度等数据时使用的测量工具。测量时,先将物体放入两个平行的测量臂之间,夹紧后就能读出测量值。旋转卡尺范式模拟卡尺测量物体的过程,这种范式常用于解决几何问题。 凸多边形的直径 应用旋转卡尺的问题中,最著名的是求凸多边形直径的问题。先将凸多边形的直径定义为,被凸多边形完全包含的最长线段长度。解决此问题的最简单的方法是,遍历每对顶点并找出其中相隔最远的两个顶点。此方法对顶点个数N具有O(N2)的时间复杂度。实际上,利用卡尺测量多边形管道直径时,并不会采用这种方法。在卡尺之间夹入管道后,旋转卡尺就能测得最大测量值。 利用旋转卡尺法的算法也同样使用这种方法。将多边形夹入两条平行的直线之间,然后沿着多边形的边旋转直线并测得顶点之间的距离。图15-20表示两条直线对一个五边形的运动过程。每幅图中,两条直线都围绕以a和b表示的两个顶点为中心旋转。此过程中得到的ab最大间距就是此多边形的直径。 图15-20 利用旋转卡尺求多边形直径的过程 有多种方法能够实现此算法,接下来介绍其中最直观的一种。代码15-16是实现这种方法的算法源代码。首先找出最左端和最右端的顶点,然后分别向两个顶点添加垂线。这两个向量的方向总是相反,所以只需保存第一条直线的方向向量。循环语句中,首先将两条直线按照逆时针方向旋转,并求出哪条直线先遇到多边形的另一个顶点。为此,需要求出当前垂线到下一个顶点之间的角度,如图15-20的实线箭头所示。这两个角度中,选择角度更小的直线进行旋转。可更新直线的方向并将旋转点更换成下一个点,以此完成直线旋转。 代码15-16 求凸多边形直径的卡尺算法 // 以逆时针方向给出顶点信息的凸多边形中,返回顶点之间的最长距离。 double diameter(const polygon& p) { int n = p.size(); // 首先找出最左端和最右端的顶点。 int left = min_element(p.begin(), p.end()) - p.begin(); int right = max_element(p.begin(), p.end()) - p.begin(); // 分别向p[left]和p[right]添加垂线。两条垂线具有相反方向, // 所以只需标出A的方向。 vector2 calipersA(0, 1); double ret = (p[right] - p[left]).norm(); // toNext[i]=表示p[i]到下一个顶点方向的单位向量 vector<vector2> toNext(n); for(int i = 0; i < n; ++i) toNext[i] = (p[(i+1) % n] - p[i]).normalize(); // a和b分别表示两条直线正在以哪个顶点为中心旋转。 int a = left, b = right; // 一直执行到两条直线旋转半圈后方向互换为止。 while(a != right || b != left) { // 查看a和b哪个到下一个顶点的角度更小 double cosThetaA = calipersA.dot(toNext[a]); double cosThetaB = -calipersA.dot(toNext[b]); if(cosThetaA > cosThetaB) { // thetaA < thetaB calipersA = toNext[a]; a = (a + 1) % n; } else { calipersA = -toNext[b]; b = (b + 1) % n; } ret = max(ret, (p[a] - p[b]).norm()); } return ret; } 除以上讲解外,代码15-16还进行了简单的优化。利用“[0, π]的余弦值总是递减”的特性,就可以不必计算角度,而是直接使用内积。通过到下一个顶点的方向就能很容易地求出直线方向,所以并不需要实际的角度值,而只需要角度的大小关系。