[知识点]计算几何I——基础知识与多边形面积

// 此博文为迁移而来，写于2015年4月9日，不代表本人现在的观点与看法。原始地址：http://blog.sina.com.cn/s/blog_6022c4720102vxaq.html

1、前言
数学在讲解析几何，现在了解一下信息学中的计算几何也是极好的，和解析几何几乎是相通的。前面的内容比较简单，甚至自己看看就行都不值得写篇文章了，但是为了后面的内容作铺垫还是写写好了。
本章节从头到尾没有去看过别人的博文，基本上都是看的刘汝佳的《算法竞赛入门经典训练指南》，所以很多讲解是原创，有误敬请指出。几何问题背景知识多，内容杂乱，有些地方容易混淆，大家细看。

2、建模

开始做几何问题之前，必须先要将一些定义写好。我们都知道，向量是有方向，大小的量，而点是坐标系上的一个个点（词穷了= =。）但是，在信息学里面的平面坐标系中，向量和点的定义是一致的——只有x,y两个变量。因为一个点就相当于把向量的起点平移到坐标原点之后的终点坐标。下面就是他们的定义：
----------------------------------------------------------------------------------------------------

struct point
{
double x,y;
point(double x=0,double y=0):x(x),y(y) {} // 构造函数，方便代码编写
};
typedef point vector; //从程序实现上看，vector只是point的别名
----------------------------------------------------------------------------------------------------

3、符号

有了向量的表示，如何计算？这里要用到运算符重载，重载+、-、*、/、<、==。如何写数学里面都学了吧，我就不详说了。大家只要注意一下dcmp的作用。
----------------------------------------------------------------------------------------------------

vector operator + (vector a,vector b) { return vector(a.x+b.x,a.y+b.y); }

vector operator - (point a,point b) { return vector(a.x-b.x,a.y-b.y); }

vector operator * (vector a,double num) { return vector(a.x*num,a.y*num); }

vector operator / (vector a,double num) { return vector(a.x/num,a.y/num); }

bool operator < (point a,point b) { if (a.x!=b.x) return (a.x

const double eps=1e-10;
int dcmp(double x) { if (abs(x)<0)?-1:1; } // abs函数为绝对值
bool operator == (vector a,vector b) { return ( dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0); }

----------------------------------------------------------------------------------------------------

dcmp函数的作用是什么？当两个值差距很小时，同样会返回相同，即返回0，用以避免精度问题带来的麻烦。

4、基本运算

1、点积（Dot）。向量a和b的点积等于两者长度的乘积乘上它们夹角的余弦。夹角为从a到b逆时针旋转的角。因此，当夹角大于90°时，点积为负。顺便写上如何求向量长度（Length）和夹角（Angle）：

----------------------------------------------------------------------------------------------------

double Dot(vector a,vector b) { return a.x*b.x+a.y*b.y; }

double Length(vector a) { return sqrt(Dot(a,a)); }

double Angle(vector a,vector b) { return acos(Dot(a,b)/Length(a)/Length(b)); } // acos为反余弦，全称arccos

----------------------------------------------------------------------------------------------------

2、叉积（Cross）。数学里暂时没学，但是作用非常大！后面的面积都需要用到它，我开始看的时候还半天没看懂。两个向量a和b的叉积等于a和b组成的三角形的有向面积的两倍。有向面积的概念是，朝着向量a看，如果向量b在你左边，则叉积大于0；如果在你右边，则小于0。


// 标准图已丢失

如上图所示，Cross(a,b)>0，且它的值为灰色阴影部分的面积。
----------------------------------------------------------------------------------------------------

double Cross(vector a,vector b) { return a.x*b.y-a.y*b.x; } // 叉积
double Area2(point a,point b,point c) { return Cross(b-a,c-a); } // 由a，b，c三个点构成的平行四边形面积

----------------------------------------------------------------------------------------------------

根据上述的点积，叉积的公式以及作用，我们可以比较方便的得到两个向量的位置关系：



括号里第一个数是点积，第二个数是叉积。设向量a水平向右，根据向量a，b的点积叉积的正负零情况，可以得到向量b的大致方向。这在以后将有很大的作用。

5、直线

在数学里面，我们已经学了直线的多种表达方式——两点式，一般式等等，信息学里面采用参数式方程。任何一条直线可以用该直线上一点P0和方向向量v表示，则直线上所有点P满足P=P0+tv，t是一个参数。则直线的参数式为：[b] P=P0+tv[/b]

1、直线交点（lineIntersection）。设两条直线分别为P+t1v和Q+t2w，设向量u=P-Q，则可以得到：



----------------------------------------------------------------------------------------------------

point lineIntersection(point p,vector v,point q,vector w)
{
vector u=p-q;
double t=Cross(w,u)/Cross(v,w);
return p+v*t;
}

----------------------------------------------------------------------------------------------------

2、点到直线的距离（Point to line）。很简单，直接用叉积即平行四边形的面积除以底边就可得到。

----------------------------------------------------------------------------------------------------

double Distance_PTL(point x,point a,point b)
{
vector v1=a-x,v2=b-x;
return (abs(Cross(v1,v2)/Length(v1)));
}

----------------------------------------------------------------------------------------------------

3、点到线段的距离（Point to segement）。注意是线段！也就是说，如果点到线段所在直线的垂线并不在线段上，那么点到线段的距离并不是直接和到直线的距离一样（很拗口，慢慢看）。设投影点为Q，如果在线段AB上，即为P到直线AB的距离；如果在射线BA上，则所求距离为PA；否则为PB。如何判断Q的位置？如下图：



前面我们讲了利用点积叉积来判断两个向量的方向，这个地方就可以用到了。大家可以自己去试试v1，v2，v3之间的点积关系能够说明什么。

----------------------------------------------------------------------------------------------------

double Distance_PTS(point x,point a,point b)
{
if (a==b) return Length(x-a);
vector v1=b-a,v2=x-a,v3=x-b;
if (dcmp(Dot(v1,v2))<0) return Length(v2); // q在射线ba上
else if (dcmp(Dot(v1,v3))>0) return Length(v3); // q在射线ab上
else return fabs(Cross(v1,v2)/Length(v1)); // q在线段ab上
}

----------------------------------------------------------------------------------------------------

4、点在直线上的投影（lineProjection）。设AB的参数式为A+tv，P的投影点Q为Q=A+t0v，由于PQ⊥AB，则Dot(v,P-(A+t0v))=0，即Dot(v,P-A)-t0×Dot(v,v)=0，解出t0=Dot(v,x-a)/Dot(v,v)，带入参数式，得Q点。

----------------------------------------------------------------------------------------------------

point lineProjection(point x,point a,point b)
{
vector v=b-a;
return (a+v*(Dot(v,x-a)/Dot(v,v)));
}

----------------------------------------------------------------------------------------------------

5、判断线段相交（Point to line）。判断条件：线段a的两个端点在线段b的两侧（即叉积符号不同）。

----------------------------------------------------------------------------------------------------

bool segmentIntersection(point a1,point a2,point b1,point b2)
{
double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
return (dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0);
}

----------------------------------------------------------------------------------------------------

6、多边形面积

将前面的知识全部铺垫好了，求多边形面积可以算是实战演练的第一步。求面积很简单，任取多边形的一个顶点（一般取读入的第一个顶点），从该点出发将多边形分成n-2个三角形，然后将面积加起来。



对于凸多边形来说，似乎比较好理解，如上图，S=S1+S2+S3+S4。凹多边形连接了每个顶点之后，可能会出现重叠面积和多边形以外的面积啊？可以确定的一点是，重叠面积=以外的面积（原因可以问你们数学老师）。所以无需顾虑，同样地，S=S1+S2+S3。
----------------------------------------------------------------------------------------------------

double polygonArea(point p[100],int n)
{
double area=0;
for (int i=1;i<=n-2;i++)
area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}

----------------------------------------------------------------------------------------------------

7、总结

比较早就开始看这一章，但是搞了很久，主要是内容太杂了，搞得有点心烦。当然今天只讲了一些最基础的东西，后面将持续更新。
计算几何II——二维几何之凸包

计算几何III——二维几何之半平面交
计算几何IV——圆与圆的计算
计算几何V——三维几何（这玩意儿前提是我能搞懂= =）