计算几何模板

处理精度判断正负

1
2
3
4
5
6
7
const double eps = 1e-8;
int dcmp(double x)//处理精度，判断正负 
{
	if(fabs(x) < eps) return 0;
	if(x < 0) return -1;
	return 1;
}

定义

点或向量

1
2
3
4
5
6
7
8
9
10
11
struct point//点或向量
{
	double x,y;
	point(double X = 0, double Y = 0){x = X, y = Y;}
};
typedef point Vector
//pair<double,double> ponit;
double dis(point a,point b)//两点距离
{
	return (double) sqrt((a.x-b.x) * (a.x-b.x) + (a.y-b.y) * (a.y-b.y))
}

直线

1
2
3
4
5
struct Line
{
	point A,B;
	Line(point a,point b){A = a, B = b;}
};

线段

1
2
3
4
5
struct Segment
{
	point A,B;
	Segment(point a,point b) {A = a, B = b;} 
};

多边形

1
2
3
4
5
struct point
{
	double x,y;
	point(double X = 0, double Y = 0){x = X, y = Y;}
}p[N];//顺序存储每个顶点

圆

1
2
3
4
5
6
struct circle
{
	point p;//圆心 
	double r;
	circle(point A,double B){p = A, r = B;}
};

运算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
double len(point a)//向量的模长 
{
	return (double) sqrt(a.x * a.x + a.y * a.y);
}
point operator + (point a,point b) {return point(a.x+b.x,a.y+b.y);}
point operator - (point a,point b) {return point(a.x-b.x,a.y-b.y);}
point operator * (point a,double k) {return point(a.x*k,a.y*k);}//数乘
double Dot(point a,point b){return a.x*b.x+a.y*b.y;}//点积 
double Cro(point a,point b){return a.x*b.y-a.y*b.x;}//叉积的模长 
point retate(point a,double altha)//逆时针旋转altha弧度
{
	return point(a.x * cos(altha) - a.y * sin(altha),a.x * sin(altha) + a.y *cos(altha));
}

位置关系判断

p是否在线段AB上

1
2
3
4
5
bool pd_PS(point p,point A,point B)//判断 p 是否在 线段 AB 上 
{
//	return dis(p,A) + dis(p,B) == dis(A,B);//这种做法精度损失较大。
	return (dcmp(Cro(p-A,p-B)) == 0) &&(dcmp(Dot(p-A,p-B) <= 0));
}

两线段是否相交(不包括交点)

1
2
3
4
5
6
bool pd_SS(point A,point B,point C,point D)//线段 AB 和线段 CD 是否相交 
{
	double f1 = Cro(C-A,B-A), f2 = Cro(D-A,B-A);//第一条线段的两个端点是否在第二条线段的两侧 
	double g1 = Cro(A-C,D-C), g2 = Cro(B-C,D-C);//第二条线段的两个端点是否在第一条线段的两侧 
	return ((f1 < 0) ^ (f2 < 0)) && ((g1 < 0) ^ (g2 < 0));
}

点在直线上投影

1
2
3
4
5
point GetLineRoot(point p,point A,point B)//p在直线AB上的投影 
{
	point v = B - a;//方向向量v 
	return A + v * (Dot(P-A,v) / Dot(v,v));
}

点在直线的那一边

假设直线上一点 P 方向向量为 v ,则 Q 在直线的那一边可以用 PQ x v 来判断。叉积为负，Q 在直线上方，

叉积为零，在直线上，叉积为正，Q 在直线下方。

点关于直线对称

1
2
3
4
point Symmetry_PL(point p,point A,point B)//p关于直线AB的对称点
{
	return p + (GetLineRoot(p,A,B)-p) * 2;//把Q延长两倍	
}

直线相交平行

1
2
3
4
5
6
point jiao_LL(point A,point B,point C,point D)//直线AB和直线CD的交点 
{
	point v = B - A, w = D - C, PQ = A - C;//直接套用公式就ok了 
	double t = Cro(PQ,w) / Cro(w,v);
	return A + t * v;
}//用的时候记得要特判一下两直线是否平行(v,w是否共线)

点在多边形内

1
2
3
4
5
6
7
8
9
10
11
12
point p[1005];
int n;
bool InPloygon(point A)//精度损失会比较大一些 
{
	double alpha = 0;
	for(int i = 1; i <= n; i++)
	{
		if(i == n) alpha += fabs(Angle(p[i]-A,p[1]-A));
		else alpha += fabs(Angle(p[i]-A,p[i+1]-A));
	}
	return dcmp(alpha-2 * pai) == 0;
}

点在凸多边形内

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int InPolygn(point p,point *a)
{
	if(Cro(p-a[1],a[2]-a[1]) < 0 || Cro(p-a[1],a[n]-a[1]) > 0) return 0;//在边界外 
	if(pd_PS(p,a[1],a[2]) || pd_PS(p,a[n],a[1])) return 2;//在边界上 
	int l = 2, r = n;
	while(r - l > 1)
	{
		int mid = (l + r)>>1;
		if(Cro(p-a[1],a[mid]-a[1]) < 0) l = mid;
		else r = mid;
	}
	if(Cro(p-a[l],a[r]-a[l]) > 0) return 0;//在边界外 
	if(pd_PS(p,a[l],a[r])) return 2;//在边界上 
	return 1;
}

圆与圆相交

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
point get(Circle c,double a)//返回c圆上圆心角为a的点的坐标 
{
	return point(c.x+cos(a)*c.r,c.y+sin(a)*r);
}
double Angle(point a){return atan2(a.x,a.y);}
int getCC(Cirle c1,Cirle c2,vector<point>& sta)//返回交点个数，sta存交点的坐标 
{
	double d = len(c1.c-c2.c);
	if(cdmp(d) == 0)
	{
		if(dcmp(c1.r-c2.r) == 0) return -1//两个圆重合 
		return 0;//没有交点 
	}
	if(dcmp(c1.r+c2.r-d) < 0) return 0;//没有交点的情况 
	if(dcmp(fabs(c1.r-c2.r)-d) > 0) return 0;
	double alpha = Angle(c2.c-c1.c);
	double deta = acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//余弦定理算夹角 
	point p1 = get(c1,alpha-deta);
	point p2 = get(c1,alpha+deta);//旋转 
	sta.push_back(p1);
	if(p1 == p2) return 1;
	sta.push_back(p2);
	return 2;
}

过一点作圆切线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int get_PC(point p,Circle c,line *v)//返回交点个数，v存储切线 
{
	Vector u = c.p - p;//BA向量 
	double dis = dis(p,c.p);
	if(dis < r) return 0;
	else if(dcmp(dis,r) == 0)//点在圆上的情况 
	{
		v[0] = retate(u,PI/2);
		return 1;
	}
	else//点不在圆上的情况 
	{
		double ang = asin(c.r/dis);
		v[0] = retate(u,ang);
		v[1] = retate(u,-ang);
		return 2;
	}
}

距离

点到线段

1
2
3
4
5
6
double dis_PS(point p,point A,point B)//点到线段的距离
{
	if(Dot(p-A,B-A) < 0) return dis(p,A);//垂足离A更近
	if(Dot(p-B,A-B) < 0) return dis(p,B);//垂足离B更近
	return abs(Cro(A-p,B-p))/dis(A,B);//垂足的距离最小
}

点到直线

1
2
3
4
double dis_PL(point p,point A,point B)//点到直线 AB 的距离 
{
	return abs(Cross(p-A,p-B)) / len(A-B);//面积除以底 
}

面积

1
2
3
4
5
6
double PloyArea(point *a,int n)
{
	double res = 0;
	for(int i = 2; i < n; i++) res += Cro(a[i]-a[1],a[i+1]-a[1]);
	return abs(res / 2);//如果你提前按极角排一下序，就可以省去取绝对值这一步
}