计算几何模板

定义

const double eps = 1e-8;
int dcmp(double x){    //判断符号
	if(fabs(x) < eps)	return 0;
	else return x < 0 ? -1 : 1;
}
typedef struct Point{   //点和向量
	double x,y;
	Point(){};
	Point(double x,double y):x(x),y(y){};
	Point operator + (const Point& e)const{   //加
		return Point(x + e.x,y + e.y);
	}
	Point operator - (const Point& e)const{   //减
		return Point(x - e.x,y - e.y);
	}
	Point operator * (const double e)const{  //乘
		return Point(x * e,y * e);
	}
	Point operator / (const double e)const{  //除
		return Point(x / e,y / e);
	}
	bool operator == (const Point &e)const{  //相等
		return !dcmp(x - e.x) && !dcmp(y - e.y);
	}
	bool operator < (const Point& e)const{  //比较
		return x < e.x || x == e.x && y < e.y;
	}
}Vector;
double Dot(Vector a,Vector b){
	return a.x * b.x + a.y * b.y;
}
double Cross(Vector a,Vector b){
	return a.x * b.y - a.y * b.x;
}

常用

double Angle(Vector A,Vector B){   //余弦定理计算夹角
	return acos(Dot(A,B)) / Length(A) / Length(B));
}
Vector Rotate(Vector A,double rad){  //向量旋转
	return Vector(A.x * cos(rad) - A.y * sin(rad),A.x * cos(rad) + A.y * sin(rad));
}
double Length(Vector A){   //计算向量的长度
	return sqrt(Dot(A,A));
}
Vector Normal(Vector A){   //计算单位法向量
	double len = Length(A);
	return Vector(-A.y / len,A.x / len);   //左转90度。
} 
double DistanceToLine(Point p,Point a,Point b){   //点到直线的距离
	Vector v1 = b - a,	v2 = p - a; 
	return fabs(Cross(v1,v2)) / Length(v1);
}
double DistanceToSegment(Point p,Point a,Point b){  //点到线段的距离
	if(a == b)	return Length(p - a);
	Vector v1 = b - a,	v2 = p - a,		v3 = p - b;
	if(dcmp(Dot(v1,v2)) < 0)	return Length(v2);
	if(dcmp(Dot(v1,v3)) > 0)	return Length(v3);
	else return DistanceToLine(p,a,b);
}

直线与直线的关系判断

int Judge(Point a1,Point b1,Point a2,Point b2){
	if(dcmp(Cross(b1 - a1,b2 - a2)) == 0){    
		if(dcmp(Cross(a2 - a1,b2 - a1)) == 0)	return 1; 	//重合  
		else	return 2;   //平行
	}else 	return 3;   //相交
}
Point intersection(Point p,Vector v,Point q,Vector w){   //两直线(p+tv)和(q+tw)求交点
	Vector u = p - q;
	double t = Cross(w,u) / Cross(v,w);
	return p + Point(v.x * t,v.y * t);
}

线段相交

//非规范相交
bool onsegment(Point p,Point a1,Point a2){ //点p是否在线段a1a2上，端点不重合。如果否则取等号，则端点可以重合。
	return dcmp(Cross(a1 - p,a2 - p)) == 0 && dcmp(Dot(a1 - p,a2 - p)) <= 0;  
}
//规范相交
bool SegmentProperIntersection1(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);
	return dcmp(c1) * dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;
} 
bool SegmentProperIntersection2(Point a1,Point b1,Point a2,Point b2){ 
    if(SegmentProperIntersection1(a1,b1,a2,b2))	return true;
    if(onsegment(a1,a2,b2) || onsegment(b1,a2,b2))	return true;
    if(onsegment(a2,a1,b1) || onsegment(b2,a1,b1))	return true;
    return false;
}

多边形面积

double getarea(){   //点的坐标为0……n-1
	double area = 0;
	for(int i=1;i<n-1;i++)
	    area += ((p[i] - p[0]) ^ (p[i+1]-p[0]));
	return fabs(area / 2);
}

皮克定理

极角排序

bool cmp1(Point a,Point b){        //atan2(y,x)表示与x正方向的夹角，大小为(-π,π]
	double t1 = atan2(a.y-st[0].y,a.x-st[0].x);
    double t2 = atan2(b.y-st[0].y,b.x-st[0].x);
	return t1 != t2 ? t1 < t2 : a.x < b.x;
}
bool cmp2(Point a,Point b){    //利用叉乘的性质
    double t = cross(a,b);
    return t > 0 || t == 0 && a.x < b.x;
}

凸包

void convex(){  //Graham扫描法
	sort(p,p+n,cmp);   //找纵坐标最小的点
	st[0] = p[0];   
	sort(p+1,p+n,cmp2);  //极角排序
	st[top=1] = p[1];;
	for(int i=2;i<n;i++){
		while(i >= 1 && cross(st[top-1],st[top],p[i]) < 0)	top--; 
		st[++top] = p[i];              cross(a,b,c)计算向量ab✖ac
	}
}
int ConvexHull(Point *p,int n,Point* ch){ //Andrew算法求凸包
	sort(p,p+n);   //先找x小的，如果x相同，找y小的
	int m = 0;
	for(int i=0;i<n;i++){
		while(m > 1 && Cross(ch[m-1] - ch[m-2],p[i] - ch[m-2]) <= 0)	m--;
		ch[m++] = p[i];
	}
	int k = m;
	for(int i=n-2;i>=0;i--){
		while(m > k && Cross(ch[m-1] - ch[m-2],p[i] - ch[m-2]) <= 0)	m--;
		ch[m++] = p[i];
	}
	if(n > 1)	m--;
	return m;
}

旋转卡壳

int rotating_caliper(Point* p,int n){  //求最远两点的距离的平方
	int t = 1;   //必须从1开始,否则边01对应的最远点在0
	double ans = 0;
	for(int i=0;i<n;i++){
		int t1 = i,	t2 = (i+1) % n;  //枚举边p1p2
		while(cross(p[t1] - p[t],p[t2] - p[t]) < cross(p[t1] - p[(t+1)%n],p[t2] - p[(t+1)%n]))//取等号会进入死循环（比如只有两个点）	
			t = (t + 1) % n;
		ans = max(ans,max(dis(p[t],p[t1]),dis(p[t],p[t2])));
	}
	return ans;
}	

半平面交

struct Line
{
	Point a,b;
	double ang;
	void getang(){	ang = atan2(b.y - a.y,b.x - a.x);}
	bool operator == (const Line& e)const{
		return !dcmp(ang - e.ang);
	}
	bool operator < (const Line& e)const{
		return ang < e.ang || !dcmp(ang - e.ang) && dcmp(Cross(b - a,e.b - a)) < 0;
	}
}s[N],q[N];
bool onright(Point a,Point b,Point p){
	return dcmp(Cross(b - a,p - a)) < 0;
}
int halfplaneintersection(){
	for(int i=0;i<n;i++)	s[i].getang();
	sort(s,s+n);
	int l,r;
	q[l = r = 0] = s[0];
	for(int i=1;i<n;i++){
		if(s[i] == s[i-1])	continue;
		while(l < r && onright(s[i].a,s[i].b,t[r-1]))	r--;
		while(l < r && onright(s[i].a,s[i].b,t[l]))	l++;
		q[++r] = s[i];
		if(l < r)	t[r-1] = intersection(q[r].a,q[r].a - q[r].b,q[r-1].a,q[r-1].a - q[r-1].b);
	}
	while(l < r && onright(q[l].a,q[l].b,t[r-1]))	r--;
	t[r] = intersection(q[l].a,q[l].a-q[l].b,q[r].a,q[r].a-q[r].b);
	if(r - l <= 1)	return 0;    //不存在
	int m = 0;
	for(int i=l;i<=r;i++)	p[m++] = t[i];	
	return m;
}

点在多边形内部

bool isPointInPolygon(Point p,Point *poly,int n){   //判断点是否在凸多边形内部
    for(int i=0;i<n;i++){
        double c1 = (poly[i] - poly[(i+1)%n]) ^ (p - poly[(i+1)%n]);
        double c2 = (poly[(i+2)%n] - poly[(i+1)%n]) ^ (p - poly[(i+1)%n]);
        if(dcmp(c1) * dcmp(c2) >= 0)    return false;   //等号表示点在多边形的边上
    }
    return true;
}