计算几何remix

搬运的remix

一.误差

计算几何中由于计算和点坐标可能会是小数，一般实数的定义使用double。（读入用%lf 读出用%f)

二.判等

不能使用等号判等，一般方法是定义一个很小的浮点数，然后将运算结果与其进行比较。

const double eps = 1e-9;//趋向于0的数
int dcmp(double x,double y)
{
	if (fabs(x - y) < eps) return 0;
	else if (x > y) return 1;
	else return -1;
}

三.头文件&常用函数&变量

1. 头文件

#include <cstdio>
#include <cmath>

2. 常用函数

const double pi = acos(-1,0);
const double inf =  1e100;
const double eps = 1e-6;

3. 变量

int fx = floor(x);//向下取整
int cx = ceil(x);//向上取整
int rx = round(x);//四舍五入取整

四.定义

1. 点的定义
   （结构体）

struct Point{
	double x;
	double y;
};

2. 向量的定义
   （同样可以用结构体表示）

typedef Point Vector;

五.运算

1. 向量运算&重载运算符
   （向量之间的加减就是x，y的加减。向量和实数之间的加减乘除，相应的变化x，y就可以了）

//向量+向量=向量 
Vector operator + (Point A,Point 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 p){
return Vector(A.x * p,A.y * p);
}
//向量/实数=向量 
Vector operator / (Vector A,double p){
return Vector(A.x / p,A.y / p);
}
//小于运算（Left Then Low排序）
bool operator < (const Point& a,const Point& b){
if(a.x == b.x)
	return a.y < b.y;
return a.x < b.x;
}

2. 内积运算（数量积、点积）

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

3. 外积运算（向量积、叉积）
   两个向量的叉积等于二者长度的乘积再乘上ta们的夹角正弦值，如果两个向量平行，叉积等于0。
   然而叉积是不满足交换率的，实际上cross(u,v)= -cross(v,u)。
   形象一点说，叉积就是两个向量组成的三角形的有向面积的两倍。（平行四边形面积）
   如果两个向量逆时针转动，那么叉积为正，顺时针转动叉积为负。在坐标系中，我们可以简单的写成x1y2-x2y1（(x1,y1)和(x2,y2)的叉积）。

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

4. 向量旋转
   x’=xcosa-ysina
   y’=xsina+ycosa

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

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

六.点和直线的常用操作

1. 判断点是否在直线上
   利用三点共线的等价条件α x β == 0，判断直线上取两不同点与待测点构成向量求叉积是否为零，即可知该点是否在直线上。

2. 计算两直线交点
   设两直线的参数方程分别为P+tv和Q+tw，
   设向量u=P-Q，交点在第一条直线的参数为t1，在第二条直线的参数为t2。
   则x，y的坐标可以列出方程，解得：

t1 = Cross(w , u) / Cross(v , w)，t1 = Cross(v , u) / Cross(v , w)

判断代码如下
（必须保证直线相交，否则将会出现除以零的情况）

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

3. 判断点是否在线段上

bool Onsegment(Point p,Point a1,Point a2){
	return dcmp(Cross(a1-p,a2-p)) == 0 && dcmp(Dot(a1-p,a2-p)) < 0;
}

4. 判断两直线是否相交

int sgn(double d){ return d<-eps?-1:d>eps; } 
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
	double c1 = Cross(a2 - a1,b1 - a1),c2 = Cross(a2 - a1,b2 - a1);
	double c3 = Cross(b2 - b1,a1 - b1),c4 = Cross(b2 - b1,a2 - b1);
	//if判断控制是否允许线段在端点处相交，根据题意添加
	if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
		bool f1 = Onsegment(b1,a1,a2);
		bool f2 = Onsegment(b2,a1,a2);
		bool f3 = Onsegment(a1,b1,b2);
		bool f4 = Onsegment(a2,b1,b2);
		bool f = (f1|f2|f3|f4);
		return f;
	}
	return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

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

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

6. 点到线段的距离

判断点和线段的关系，可以利用点积的性质：夹角小于九十度，得到的点积为正；夹角大于九十度，得到的点积为负。
要注意向量的方向

double Dis2(Point P,Point A,Point 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);
}

7. 点在直线上的投影
   把直线AB写成参数式A+tv（v为向量AB），并且设Q的参数为t0，则Q=A+t0v。
   根据PQ垂直于AB，我们可以得到两个向量的点积为0：
   Dot(v，P-(A+t0v))=0
   根据分配率，有：Dot(v,P-A)-t0Dot(v,v)
   这样就可以解出Q点了。

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

8. 多边形面积

Vector Area(Point *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;
}

9. 判断点是否在多边形内

在竞赛中一般采用转角法
基本思路就是看这个多边形相对于这个点旋转了多少度
具体来说，我们把多边形每条边的转角加起来，
如果是360°，说明在多边形内
如果是0°，说明在多边形外
如果是180°，说明在多边形边上
这个方法是适用范围非常广泛，但是如果我们直接按照定义实现，
需要计算大量的反三角函数，不仅速度慢，而且容易产生精度误差，所以我们一般是这么实现的：
我们假想有一条向右的射线，统计多边形穿过这条射线正反多少次，我们把这个数记为绕数wn（Winding Number）
逆时针穿过时，wn++，顺时针穿过时，wn–

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

int Init(Point P,Point *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;                                     //外部 
}