计算几何

计算几何是一门兴起于二十世纪七十年代末的计算机科学的一个分支，主要研究解决几何问题的算法。在现代工程和数学领域，计算几何在图形学、机器人技术、超大规模集成电路设计和统计等诸多领域有着十分重要的应用。

计算几何问题的输入一般是关于一组几何对象的描述，如一组点、一组线段，或者一个多边形的按逆时针顺序排列的一组顶点。输出常常是对有关这些对象的问题的回答，如是否直线相交，是否为一个新的几何对象，如顶点集合的凸包。

本文将介绍一些平面上的计算几何算法。在这些算法中，每个输入对象都是一组点{p1,p2,p3...}，其中每个pi=(xi,yi)。

点与向量

向量也称为矢量，是一种既有大小又有方向的量。在计算几何中，某个向量从点A指向点B，记为向量AB。

二维直角坐标系中的向量可以沿坐标轴分解，向量AB=xi+yj。所以向量AB可以用实数(x,y)来表示，同理n维向量也可以用n个实数来表示。因此一个点也可以表示一个向量。

向量的基本运算包括加减法，数乘，点积，叉积，混合积等。

如果有向量a=(x1,y1)，向量b=(x2,y2)。那么：

1. 向量a的长度|a|=sqrt(x²+y²)。
2. 加法的定义为：向量a+向量b=(x1+x2,y1+y2)。减法同理。
3. 数乘是向量和实数之间的运算，如果是一个正实数与向量相乘，向量只发生大小变化，方向不变。乘以负数则会将原向量反向，乘以0则变成零向量。
4. 点积，也叫数量积或内积，得到的结果是一个实数。定义向量a·向量b=x1y1+x2y2。点积满足交换律不满足结合律。
5. 叉积，也叫矢量积或外积。定义向量aX向量b=x1y2-x2y1，即两个向量形成的平行四边行的面积。

因为计算几何中经常涉及精度问题，需要对一个很小的数判断正负，为了修正误差，需要引入一个极小量eps。

1 const double eps=1e-8;2 int dcmp(double x){3     if (fabs(x)<eps) return 0;4     if (x>0) return 1;5     return -1;6 }

计算几何误差修正

我们设计一个二维点类，使它可以进行一些常见的向量运算。点类也是其它计算几何问题的基础。

 1 const double pi=acos(-1.0); 2 inline double sqr(double x){return x*x;}// 平方 3 struct Point{ 4     double x,y; 5     Point(){} 6     Point(double a,double b):x(a),y(b){}// 构造 7     void input(){scanf("%lf%lf",&x,&y);}// 读入 8     double norm(){return sqrt(sqr(x)+sqr(y));}// 到原点的距离 9     /** 运算 **/10     friend Point operator+(const Point &a,const Point &b){11         return Point(a.x+b.x,a.y+b.y);12     }13     friend Point operator-(const Point &a,const Point &b){14         return Point(a.x-b.x,a.y-b.y);15     }16     friend bool operator==(const Point &a,const Point &b){17         return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;18     }19     friend Point operator*(const Point &a,const double &b){20         return Point(a.x*b,a.y*b);21     }22     friend Point operator*(const double &a,const Point &b){23         return Point(a*b.x,a*b.y);24     }25     friend Point operator/(const Point &a,const double &b){26         return Point(a.x/b,a.y/b);27     }28 };29 double det(const Point &a,const Point &b){return a.x*b.y-a.y*b.x;}// 叉乘30 double dot(const Point &a,const Point &b){return a.x*b.x+a.y*b.y;}// 点乘31 double dist(const Point &a,const Point &b){return (a-b).norm();}// 距离

点类

线段的性质

线段

我们称p1和p2为线段p₁p₂的端点。如果要考虑p1和p2之间的顺序，称线段为有向线段。如果有向线段p₁p₂的起点p1是坐标原点(0,0)，则可以称其为向量p2。

实现一个线段类，用有向线段上的两个端点a->b来表示一个线段。

1 struct Line{2     Point a,b;3     Line(){}4     Line(Point x,Point y):a(x),b(y){}5 };

线段类

叉积

叉积是线段相关算法的中心。考虑图示两个向量p1和p2。

可以把叉积p1Xp2看作是由点(0,0),p1,p2和p1+p2=(x1+x2,y1+y2)所形成的平行四边形的面积。也可以把叉积定义为一个矩阵的行列式。

p1Xp2=x1y2-x2y1=-p2Xp1

如果p1Xp2为正数，则相对于原点(0,0)来说，p1在p2的顺时针方向上；如果p1Xp2为负数，则p1在p2的逆时针方向上。

如果叉积为0的话，说明这两个向量是共线的，它们指向相同的方向或相反的方向。

如果公共端点不是原点而是p0，判断有向线段p₀p₁在有向线段p₀p₂的顺时针方向还是逆时针方向，只需定义p1‘=p1-p0，p2‘=p2-p0。

然后计算叉积p1‘Xp2‘=(p1-p0)X(p2-p0)=(x1-x0)(y2-y0)-(x2-x0)(y1-y0)

如果该叉积为正，则有向线段p₀p₁在有向线段p₀p₂的顺时针方向；叉积为负，在逆时针方向；叉积为0，它们共线。

折线的拐角方向判断

对于两条连续的有向线段p₀p₁和p₁p₂，我们可以判断它们在点p1处是左转还是右转。

如图所示，只需检查有向线段p₀p₂是在有向线段p₀p₁的顺时针方向还是逆时针方向即可。

计算叉积(p2-p0)X(p1-p0)，如果为负，则在p1左转；如果为正，在p1右转；如果为0，不转向是直线。

线段相交 

如果一个线段的两个端点在一条直线的两边，我们称这条线段跨越了这条直线。

两个线段相交，当且仅当下面两个条件中至少一个成立：

1. 一个线段的某一个端点在另一条线段上。
2. 每个线段都跨越另一条线段所在的直线。

下面的代码实现了这一思想。

如果线段规范相交返回2，非规范相交返回1，不相交返回0。

过程Direction，利用叉积计算出线段的相对方位。过程OnSegment，用于确定一个与某一线段共线的点是否位于线段上。

 1 double Direction(Point pi,Point pj,Point pk){ 2     return det(pk-pi,pj-pi); 3 } 4 bool OnSegment(Point pi,Point pj,Point pk){ 5     if (min(pi.x,pj.x)<=pk.x&&pk.x<=max(pi.x,pj.x)&&min(pi.y,pj.y)<=pk.y&&pk.y<=max(pi.y,pj.y)) return true; 6     return false; 7 } 8 int SegCrossSeg(Line u,Line v){ 9     Point p1=u.a,p2=u.b,p3=v.a,p4=v.b;10     int d1=dcmp(Direction(p3,p4,p1));11     int d2=dcmp(Direction(p3,p4,p2));12     int d3=dcmp(Direction(p1,p2,p3));13     int d4=dcmp(Direction(p1,p2,p4));14     if ((d1^d2)==-2&&(d3^d4)==-2) return 2;15     if ((d1==0&&OnSegment(p3,p4,p1))||16         (d2==0&&OnSegment(p3,p4,p2))||17         (d3==0&&OnSegment(p1,p2,p3))||18         (d4==0&&OnSegment(p1,p2,p4))) return 1;19     return 0;20 }

线段相交

处理流程：计算出每个端点pi相对于另一线段的相对方位di。如果所有的方位都非0，则可以容易的确定线段是否规范相交，1^-1=-2，如果d1与d2或d3与d3符号不同的话，它们的异或值就是-2。否则我们判断端点与线段共线的情况中是否有边界情况，如果某点在另一条线段上，则它们不规范相交。其它情况则不相交。

线段求交问题

基本概念

给定平面上n条线段，判断其中是否有两条线段相交。

容易想到枚举每两条线段，判断它们是否相交，复杂度O(n²)。这里介绍一种利用扫描线的O(nlogn)的算法。

算法描述

扫描算法在许多计算几何算法中都用到了。

扫描线是一个假想的垂直的线，它从左到右移动。

为了简化问题，我们先假设输入中没有一条线段是垂直的，没有三条线段共线。最后对算法稍作改动，使假设不成立时仍能正常工作。

排序线段

因为假定不存在垂直线段，所以任何与给定垂直扫描线相交的输入线段与其有一个交点。根据交点的y坐标对交点进行排序。

对于两条线段s1和s2。如果一条横坐标为x的垂直扫描线与这两条线段都相交，则说这两条线段在x是可比的。如果s1和s2在x处是可比的，并且在x处，s1与扫描线的交点比s2与扫描线的交点高，则说x处s1位于s2之上，写作s1>r^s2。

如图a>r^c，a>t^b，b>t^c，a>t^c，b>u^c，线段d与其他任何线段都不可比。

当扫描线穿过两条线段的交点时，如图，它们在扫描线中的顺序颠倒了，扫描线v和w分别位于线段e和f的交点的左边和右边，有e>v^f，f>w^e。

注意我们假定没有三条线段相交于同一点，所以必有某条垂直扫描线x，使得相交线段e和f在扫描线中时连续的。任何穿过图中阴影区域的扫描线，e与f都在其中连续。

扫描线的移动

典型的扫描算法要维护下列两组数据。

扫描线状态：与扫描线相交的物体之间的关系

事件点调度：从左向右排列的x坐标的序列，它定义了扫描线的暂停位置。称每个这样的暂停位置为事件点。扫描线状态仅在事件点处才会发生变化。

对于某些算法，事件点调度是随算法执行而动态确定的。我们现在讨论的仅仅是基于输入数据的简单性质静态的确定事件点，这里每条线段的端点都是事件点。

我们通过x坐标，从左向右对线段的端点进行排序。如果有两个或更多端点位于同一条垂直线上，即有相同的x坐标，就将所有的左端点放置在右端点之前。对于x坐标相同的左端点们，将y坐标小的放置在前面，右端点也是如此。

在扫描过程中，当遇到左端点时，就把改线段插入到扫描线状态中，并且当遇到其右端点时，就把它从扫描线状态中删除。每当两条线段在扫描线中第一次变为连续时，就检查它们是否相交。

因此在扫描中有以下操作：

• 把线段s插入到扫描线中。
• 把线段s从扫描线中删除。
• 返回紧邻线段s的上下两条线段。

可以用各种平衡二叉树用O(logn)的复杂度来维护这些操作。在树上用叉积来比较两条线段的相对顺序。

具体流程

1. 将所有线段的端点加入事件队列Q。
2. 设S为所有和扫描线L相交的线段集合，初始S为空。
3. 从Q中取出下一个横坐标x最小的端点E。
4. 如果E是一条线段b的左端点，将此线段加入S中，接着找出S中与b上下相邻的两条线段a和c，检查是否(a,b)或(b,c)相交。
5. 如果E是一条线段b的右端点，从S中找出b上下相邻的两条线段a和c，将b从S中删去，同时判断(a,c)是否相交。
6. 如果Q中还有事件，返回第3步，否则算法结束。