计算几何（具象的难题）

前言：每到周天，就比较怠惰，所以随便写点什么

二维几何基础

基础运算

既然说是 基础，那我们就要从计算几何最基本的表示开始：

一般的，在计算机中我们会广泛的 使用向量

个人而言比较偏爱用结构体

struct node{
    double x,y;
    node (double xx=0,double yy=0)
    {
    	x=xx;y=yy;
	}
};
node po[N];

//向量+向量=向量 
node operator + (const node &a,const node &b){return node (a.x+b.x,a.y+b.y);}
//向量-向量=向量 
node operator - (const node &a,const node &b){return node (a.x-b.x,a.y-b.y);}
//向量*实数=向量 
node operator * (const node &a,const double &b){return node (a.x*b,a.y*b);}
//向量/实数=向量 
node operator + (const node &a,const double &b){return node (a.x/b,a.y/b);}

在这里涉及到了向量的运算：

向量之间的加减就是x，y的加减

向量和实数之间的加减乘除，相应的变化x，y就可以了

然而我们漏掉了一种运算， 向量乘向量（我们基本上不会涉及除法）

向量的乘法有两种

• 点积

两个向量的点积等于二者长度的乘积再乘上ta们的夹角 余弦值

如果 两个向量垂直，点积等于0

因为余弦是偶函数，所以点积满足交换律，

我们如果用点的形式存储向量的话，点积就可以简单的写成 x1x2+y1y2

//求点积 
double Dot(node x,node y) {return x.x*y.x + x.y*y.y;}
//利用点积求向量的长度
double Length(node x) {return sqrt(Dot(x,x));}
//求两向量的夹角
double Angle(node x,node y) {return acos(Dot(x,y) / Length(x) / Length(y));} 

• 叉积

两个向量的叉积等于二者长度的乘积再乘上ta们的夹角 正弦值

如果 两个向量平行，叉积等于0

然而叉积是不满足交换率的，实际上 cross(u,v)=-cross(v,u)

形象一点说，叉积就是两个向量组成的三角形的有向面积的两倍

简单来说，如果两个向量逆时针转动，那么叉积为正，顺时针转动叉积为负

在坐标系中，我们可以简单的写成 x1*y2-x2*y1

//求叉积 
double Cross(node x,node y) {return x.x*y.y - x.y*y.x;}
//面积
double SS(node x,node y) {return Cross(x,y);} 

两个向量的位置关系

向量旋转

这是很重要的一种运算

直接上公式：

x'=xcosa-ysina

y'=xsina+ycosa

其中a是逆时针旋转的角度

node rotate(node x,double a)    //a是弧度 
{
	return node(x.x*cos(a)-x.y*sin(a),x.x*sin(x)+x.y*cos(a));
}

 点和直线

直线的参数表示：

一般的我们可以用Ax+By+C=0表示一条直线，计算机中我们也可以通过存储A,B,C来确定一条直线

不过在这里我们还要介绍一种更常用的方法：

用直线上一点P0和方向向量v表示

直线上所有点P都满足P=P0+tv，其中t成为参数

如果已知直线上两个不同的点A,B，则 方向向量为（B-A），所以参数方程为 A+（B-A）t

参数方程最方便的地方在于直线，射线和线段的方程形式是一样的，区别仅仅在于参数：

直线的t没有限制，射线的t>0，线段的t介于0~1之间

• 直线交点

设两直线的参数方程分别为P+tv和Q+tw

设向量u=P-Q，交点在第一条直线的参数为t1，在第二条直线的参数为t2

则x，y的坐标可以列出方程，解得：

node Lineintersection(node P,node v,node Q,node w)
{
	node u=P-Q;
	double t=Cross(w,u)/Cross(v,w);
	return P+v*t;
}

• 点到直线的距离

点到直线的距离可以通过叉积求出

double Dis1(node P,node A,node B)
{
	return fabs(Cross(P-A,B-A))/Len(B-A);   //如果不取绝对值，得到的是有向距离 
}

• 点到线段的距离

这个问题有一点复杂，主要分为两种情况

可见我们需要判断点和线段的关系，这时候我们就可以利用 点积的性质了： 夹角小于九十度，得到的点积为正；夹角大于九十度，得到的点积为负

注意向量的方向

int dcmp(double x)
{
	if (fabs(x)<eps) return 0;
	else if (x>0) return 1;
	else return -1;
}

double Dis2(node P,node A,node B)
{
	if (A==B) return Len(P-A);
	node v1=B-A,v2=P-A,v3=P-B;
	if (dcmp(Dot(v1,v2))<0) return Len(v2);
	else if (dcmp(Dot(v1,v3))>0) return Len(v3);
	else return fabs(Cross(v1,v2))/Len(v1);
}

• 点在直线上的投影

虽然前边介绍的运算避免了计算 P在AB上的投影点Q，但是我们有时候还是需要把ta求出来的

为此，我们把直线AB写成参数式A+tv（v为向量AB）

并且设Q的参数为t0，则 Q=A+t0*v

根据PQ垂直于AB，我们可以得到 两个向量的点积为0：

Dot(v，P-(A+t0*v))=0

根据分配率，有：Dot(v,P-A)-t0*Dot(v,v)

这样就可以解出Q点了

node lineprojection(node P,node A,node B)
{
	node v=B-A;
	return A+v*(Dot(v,P-A)/Dot(v,v));
}

• 线段相交判定

给定两条线段，判断是否相交

我们定义“规范相交”为两线段瞧好有一个公共点，且不在任一线段的端点

线段规范相交的充要条件是： 每条线段的两个端点都在另一条线段的两侧

bool XiangJiao(node a1,node a2,node b1,node b2)            //我看不出有什么起高大上函数名的必要 
{
    double c1=Cross(a2-a1,b1-a1);
	double c2=Cross(a2-a1,b2-a1);
	double c3=Cross(b2-b1,a1-b1);
	double c4=Cross(b2-b1,a2-b1);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

如果允许线段在端点处相交，那就有一点困难了：

如果c1和c2都等于0，那么两线段共线，这是可能会有部分重叠的情况

如果c1和c2都不为0，则只有一种相交方法，即 某一个点在另一条线段

还需要如下一段代码，判断是否一点在一条线段上

bool Onit(node p,node a,node b)
{
	return dcmp(Cross(a-p,b-p))==0 && dcmp(Dot(a-p,b-p))<0;    //叉积为0(平行)  点积为负(夹角大于90°) 
}

多边形

直接上代码：

node Area(node *P,int n)
{
	double ans=0;
	for (int i=1;i<n;i++)
	    ans+=Cross(P[i]-P[0],P[i+1]-P[0]);
	return fabs(ans)/2;
}

点在多边形内的判定

我们在竞赛中一般采用 转角法

基本思路就是看这个多边形相对于这个点旋转了多少度

具体来说，我们把多边形每条边的转角加起来，

如果是360°，说明在 多边形内

如果是0°，说明在 多边形外

如果是180°，说明在 多边形边上

这个方法是适用范围非常广泛，但是如果我们直接按照定义实现，

需要计算大量的反三角函数，不仅速度慢，而且容易产生精度误差，所以我们一般是这么实现的：

我们假想有一条向右的射线，统计多边形穿过这条射线正反多少次，我们把这个数记为绕数wn（Winding Number）

逆时针穿过时，wn++，顺时针穿过时，wn--

int Init(node P,node *po)
{
	int wn=0;
	po[0]=po[n]; po[n+1]=po[1];
	for (int i=1;i<=n;i++)
	{
		if (Onit(P,po[i],po[i+1])) return -1;     //在边界上
		
		int k=dcmp(Cross(po[i+1]-po[i],P-po[i]));
		int d1=dcmp(po[i].y-P.y); 
		int d2=dcmp(po[i+1].y-P.y);
		if (k>0&&d1<=0&&d2>0) wn++;               //逆时针 
		if (k<0&&d2<=0&&d1>0) wn--;               //顺时针 
	}
	if (wn!=0) return 1;                          //内部 
	return 0;                                     //外部 
}