二维计算几何全家桶

快网络赛了但是计算几何还一点不会,所以最近狠狠恶补了一下计算几何的知识,但是由于本人实力有限,暂时还没有学会三维的计算几何,所以本文只介绍二维计算几何中一些比较常用到的知识

符号函数

符号函数是计算几何中经常用到的一个函数,它虽然很简单,但是在帮助我们判断几何中的位置和方向时有着重要的作用

int sgn(double x){
	if(fabs(x) < eps)return 0;
	else return x<0?-1:1;
}

点和向量

定义二维平面的点

struct Point{
	double x,y;
	Point(){};
	Point(double x,double y):x(x),y(y){}
};

用库函数hypot()计算直角三角形的斜边长

double Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}

计算两点之间的距离

double Dist(Point A,Point B){return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));}

向量的基本运算

Point operator + (Point B){return Point(x+B.x,y+B.y);}//向量的加法  
Point operator - (Point B){return Point(x-B.x,y-B.y);}//向量的减法
Point operator * (double k){return Point(x*k,y*k);}//向量与实数相乘
Point operator / (double k){return Point(x/k,y/k);}//向量与实数相除
Point operator == (Point B){return sgn(x-B.x) == 0&&sgn(y-B.y) == 0;}//向量相等

求向量A与B的点积

double Dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}

求向量A的长度

double Len(Vector A){return sqrt(Dot(A,A));}

由于在开根号的时候可能会造成精度的丢失,所以我们更经常去求向量A长度的平方

double Len2(Vector A){return Dot(A,A);}

求向量A与B的夹角

double Angle(Vector A,Vector B){return acos(Dot(A,B)/Len(A)/Len(B));}

求向量A与B的叉积

叉积有一个很重要的作用就是帮助我们判断位置关系,若 \(A×B>0\),B在A的逆时针方向;若\(A×B<0\),B在A的顺时针方向;若\(A×B=0\),B与A共线

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);}

利用叉积,我们也可以检查两个向量是否平行或重合

bool Parallel(Vector A,Vector B){return sgn(Cross(A,B)) == 0;}//返回true表示平行或重合

向量A逆时针旋转角度rad

Vector Rotate(Vector A,double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

求单位法向量,即逆时针旋转90°,然后取单位值

Vector Normal(Vector A){return Vector(-A.y/Len(A),A.x/Len(A));}

直线

直线的表示

struct Line{
	Point p1,p2;
	Line(){}
	Line(Point p1,Point p2):p1(p1),p2(p2){}
	//根据一个点和倾斜角确定直线
	Line(Point p,double angle){
		p1=p;
		if(sgn(angle-pi/2) == 0){p2 = (p1 + Point(0,1));}
		else{p2 = (p1 + Point(1,tan(angle)));}
	}
	//ax+by+c=0
	Line(double a,double b,double c){
		if(sgn(a) == 0){
			p1 = Point(0,-c/b);
			p2 = Point(1,-c/b);
		}
		else if(sgn(b) == 0){
			p1 = Point(-c/a,0);
			p2 = Point(-c/a,1);
		}
		else{
			p1 = Point(0,-c/b);
			p2 = Point(1,(-c-a)/b);
		}
	}
}

线段的表示

可以用两个点表示线段,起点P1,终点P2

typedef Line Segment;

点和直线的位置关系

int Point_Line_relation(Point p,Line v){
	int c = sgn(Cross(p-v.p1,v-p2-v.p1));
	if(c<0)return 1;//1:p在v的左侧
	if(c>0)return 2;//2:p在v的右侧
	return 0;//0:p在v上
}

点和线段的位置关系

bool Point_on_seg(Point p,Line v){//0:点p不在线段v上;1:点p在线段v上
	return sgn(Cross(p-v.p1,v.p2-v.p1)) == 0&&sgn(Dot(p-v.p1,p-v.p2) <= 0);
}

点到直线的距离

double Dis_point_line(Point p,Line v){
	return fabs(Cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2);
}

点在直线上的投影

Point Point_line_proj(Point p,Line v){
	double k = Dot(v.p2-v.p1,p-v.p1)/Len2(v.p2-v.p1);
	return v.p1+(v.p2-v.p1)*k;
}

点关于直线的对称点

Point Point_line_symmetry(Point p,Line v){
	Point q = Point_line_proj(p,v);
		return Point(2*q.x-p.x,2*q.y-p.y);
}

点到线段的距离

double Dis_point_seg(Point p,Segment v){
	if(sgn(Dot(p-v.p1,v.p2-v.p1) < 0) || sgn(Dot(p-v.p2,v.p1-v.p2)) < 0)
		return min(Dsitance(p,v.p1),Distance(p,v.p2));
		return Dis_point_line(p,v);
}

两条直线的位置关系

int Line_relation(Point p,Segment v){
	if(sgn(Cross(v1.p2-v1.p1,v2.p2-v2.p1)) == 0){
		if(Point_line_relation(v1.p1,v2) == 0)return 1;//1:重合
		else return 0;//0:平行
	}
	return 2;//2:相交
}

两条直线的交点

Point Cross_point(Point a,Point b,Point c,Point d){//线段1:ab,线段2:cd
		double s1 = Cross(b-a,c-a);
		double s2 = Cross(b-a,d-a);
		return Point(c.x*s2-d.x*s1,c.y*s2-d.y*s1)/(s2-s1);
}

判断两条线段是否相交

bool Cross_segment(Point a,Point b,Point c,Point d){
	double c1 = Cross(b-a,c-a),c2 = Cross(b-a,d-a);
	double d1 = Cross(d-c,a-c),d2 = Cross(d-c,b-c);
	return sgn(c1)*sgn(c2) < 0 && sgn(d1)*sgn(d2) < 0;
}

多边形

判断点和多边形的位置关系

int Point_in_polygon(Point pt,Point *p,int n){//点pt,多边形Point *p
	for(int i = 0;i < n;i++){
		if(p[i] == pt)return 3;//3:点在多边形的顶点上
	}
	for(int i = 0;i < n;i++){
		Line v = Line(p[i],p[(i+1)%n]);
		if(Point_on_seg(pt,v))return 2;//2:在多边形的边上
	}
	int num = 0;
	for(int i = 0;i < n;i++){
		int j = (i+1)%n;
		int c = sgn(Cross(pt-p[j],p[i]=p[j]));
		int u = sgn(p[i].y-pt.y);
		int v = sgn(p[j].y-pt.y);
		if(c > 0 && u < 0 && v >= 0)num++;
		if(c < 0 && u >= 0 && v < 0)num--;
	}
	return num != 0;//1:点在内部;0:点在外部
}

多边形的面积

double Polygon_area(Point *p,int n){
	double area = 0;
	for(int i = 0;i < n;i++)
		area += Cross(p[i],p[i+1]%n);
	return area/2;
}

多边形的重心

Point Polygon_center(Point *p,int n){
	Point ans(0,0);
	if(Polygon_area(p,n) == 0)return ans;
	for(int i = 0;i < n;i++)
		ans = ans+(p[i] + p[(i+1)%n])*Cross(p[i],p[(i+1)%n]);
	return ans/Polygon_area(p,n)/6;
}

凸包

struct Point{
	double x,y;
	Point(){}
	Point(double x,double y):x(x),y(y){}
	Point operator + (Point B){return Point(x+B.x,y+B.y);}
	Point operator - (Point B){return Point(x-B.x,y-B.y);}
	bool operator == (Point B){return sgn(x-B.x) == 0 && sgn(y-B.y) == 0;}
	bool operator < (Point B){
		return sgn(x - B.x) < 0 ||(sgn(x - B.x) == 0 && sgn(y - B.y) < 0);}
};
typedef Point Vector;
double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
double Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}
//Convex_hull 函数求凸包,凸包顶点放在ch中,返回值是凸包的顶点数
int Convex_hull(Point *p,int n,Point *ch){
	n=unique(p,p+n)-p;//去掉重复点
	sort(p,p+n);//对点排序:按x从小到大排序,如果x相同按y排序
	int v=0;
	//求下凸包,如果p[i]是右拐弯的,这个点就不在凸包上,退回
	for(int i = 0;i < n;i++){
		while(v > 1 && sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-1])) <= 0)
			v--;
		ch[v++]=p[i];
	}
	int j=v;
	//求上凸包
	for(int i= n-2 ;i >= 0;i--){
		while(v > j && sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-1])) <= 0)
			v--;
		ch[v++]=p[i];
	}
	if(n>1)v--;
	return v;
	}

平面最近点对

struct Point{double x,y;};
double Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}
bool cmpxy(Point A,Point B){return sgn(A.x - B.x) < 0||(sgn(A.x-B.x) == 0 && sgn(A.y - B.y) < 0);}//排序:先对x坐标排序,再对y坐标排序
bool cmpy(Point A,Point B){return sgn(A.y - B.y) < 0;}
Point p[maxn],tmp_p[maxn];
double Closest_Pair(int left,int right){
	double dis=INF;
	if(left==right)return dis;//只剩一个点
	if(left+1==right)return Distance(p[left],p[right]);//只剩两个点
	int mid=(left+right)/2;
	double d1=Closest_Pair(left,mid);//求s1内的最近点对
	double d2=Closest_Pair(mid+1,right);//求s2内的最近点对
	dis=min(d1,d2);
	int k=0;
	for(int i=left;i<=right;i++)//在s1和s2中间附近找可能的最小点对
		if(fabs(p[mid].x-p[i].x)<=dis)//按x坐标查找
			tmp_p[k++]=p[i];
		sort(tmp_p,tmp_p+k,cmpy);//按y坐标排序,用于剪枝,这里不能按x坐标排序
		for(int i=0;i<k;i++)
			for(int j=i+1;j<k;j++){
				if(tmp_p[j].y-tmp_p[i].y>=dis)break;//剪枝
				dis=min(dis,Distance(tmp_p[i],tmp_p[j]));
			}
	return dis;//返回最小距离
	}

旋转卡壳

struct Point{
    double x,y;
    Point(double X=0,double Y=0){x=X,y=Y;}
};
struct Vector{
    double x,y;
    Vector(double X=0,double Y=0){x=X,y=Y;}
};
inline Vector operator-(Point x,Point y){// 点-点=向量
    return Vector(x.x-y.x,x.y-y.y);
}
inline double cross(Vector x,Vector y){ // 向量叉积
    return x.x*y.y-x.y*y.x;
}
inline double operator*(Vector x,Vector y){ // 向量叉积
    return cross(x,y);
}
inline double len(Vector x){ // 向量模长
    return sqrt(x.x*x.x+x.y*x.y);
}

int stk[50005];
bool used[50005];
vector<Point> ConvexHull(Point* poly, int n){ // Andrew算法求凸包
    int top=0;
    sort(poly+1,poly+n+1,[&](Point x,Point y){
        return (x.x==y.x)?(x.y<y.y):(x.x<y.x);
    });
    stk[++top]=1;
    for(int i=2;i<=n;i++){
        while(top>1&&sgn((poly[stk[top]]-poly[stk[top-1]])*(poly[i]-poly[stk[top]]))<=0){
            used[stk[top--]]=0;
        }
        used[i]=1;
        stk[++top]=i;
    }
    int tmp=top;
    for(int i=n-1;i;i--){
        if(used[i]) continue;
        while(top>tmp&&sgn((poly[stk[top]]-poly[stk[top-1]])*(poly[i]-poly[stk[top]]))<=0){
            used[stk[top--]]=0;
        }
        used[i]=1;
        stk[++top]=i;
    }
    vector<Point> a;
    for(int i=1;i<=top;i++){
        a.push_back(poly[stk[i]]);
    }
    return a;
}

struct Line{
    Point x;Vector y;
    Line(Point X,Vector Y){x=X,y=Y;}
    Line(Point X,Point Y){x=X,y=Y-X;}
};

inline double DistanceToLine(Point P,Line x){// 点到直线的距离
    Vector v1=x.y, v2=P-x.x;
    return fabs(cross(v1,v2))/len(v1);
}

double RoatingCalipers(vector<Point> poly){// 旋转卡壳
    if(poly.size()==3) return len(poly[1]-poly[0]);
    int cur=0;
    double ans=0;
    for(int i=0;i<poly.size()-1;i++){
        Line line(poly[i],poly[i+1]);
        while(DistanceToLine(poly[cur], line) <= DistanceToLine(poly[(cur+1)%poly.size()], line)){
            cur=(cur+1)%poly.size();
        }
        ans=max(ans, max(len(poly[i]-poly[cur]), len(poly[i+1]-poly[cur])));
    }
    return ans;
	}

半平面的表示

struct Line{
	Point p;//直线上的一个点
	Vector V;//方向向量,它的左边是半平面
	double ang;//极角,从x正半轴旋转到v的角度
	Line(){};
	Line(Point p,Vector v):p(p),v(v){ang = atan2(v.y,v.x);}
	bool operator < (Line &L){return ang < L.ang;}//用于排序
};  

半平面交

struct Point(){
	double x,y;
	Point(){};
	Point(double x,double y):x(x),y(y){};
	Point operator + (Point B){return Point(x + B.x,y + B.y);}
	Point operator - (Point B){return Point(x - B.x,y - B.y);}
	Point operator * (double k){return Point(x*k,y*k);}
};
typedef Point Vector;
double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
struct Line{
	Point p;
	Vector v;
	double ang;
	Line(){};
	Line(Point p,Vector v):p(p),v(v){ang = atan2(v.y,v.x);}
	bool operator < (Line &L){return ang < L.ang;}
};
bool OnLeft(Line L,Point p){return sgn(Cross(L.v,p-L.p) > 0);}
//点p在线L的左边,即点p在线L的外面
Point Cross_point(Line a,Line b){//两直线交点
	Vector u = a.p-b.p;
	double t = Cross(b.v,u)/Cross(a.v,b.v);
	return a.p+a.v*t;
}
vector<Point> HPI(vector<Line>L){//求半平面交,返回凸多边形
	int n = L.size();
	sort(L.begin(),L.end());
	int first,last;//指向双端队列的第一个元素和最后一个元素
	vector<Point>p(n)//两个相邻的半平面的交点
	vector<Line>q(n);//双端队列
	vector<Point>ans//半平面交形成的凸包
	q[first = last = 0] = L[0];
	for(int i = 1;i < n;i++){
		//删除尾部的半平面
		while(first < last && !OnLeft(L[i],p[last-1]))last--;
		//删除首部的半平面
		while(first < last && !OnLeft(L[i],p[first]))first++;
		q[++last] = L[i];//将当前的半平面加入双端队列尾部
		//极角相同的两个半平面,保留左边
		if(fabs(Cross(q[last].v,q[last-1].v)) < eps){
			last--;
			if(OnLeft(q[last],L[i].p))q[last]=L[i];
		}
		//计算队列尾部半平面交点
		if(first < last)p[last-1] = Cross_point(q[last-1],q[last]);
	}
	//删除队尾无用的半平面
	while(first < last && !OnLeft(q[first],p[last-1]))last--;
	if(last-first <= 1)return ans;//空集
	p[last] = Cross_point(q[last],q[first]);//计算队列首尾部的交点
	for(int i = first;i <= last ;i++)ans.push_back(p[i]);//复制
	return ans;//返回凸多边形
	}

圆的定义

struct Circle{
	Point c;//圆心
	double r//半径
	Circle(){}
	Circle(Point c,double r):c(c),r(r){}
	Circle(double x,double y,double _r){c = Point(x,y);r = _r;}
};

点和圆的位置关系

int Point_circle_relation(Point p,Circle C){
	double dst = Distance(p,C.c);
	if(sgn(dst - C.r) < 0)return 0;//0:点在圆内
	if(sgn(dst - C.r) == 0)return 1;//1:点在圆上
	return 2;//2:点在圆外
}

直线和圆的位置关系

int Line_circle_relation(Line v,Circle C){
	double dst = Dis_point_line(C.c,v);
	if(sgn(dst - C.r) < 0)return 0;//0:直线和圆相交
	if(sgn(dst - C.r) == 0)return 1;//1:直线和圆相切
	return 2;//2:直线在圆外
}

线段和圆的位置关系

int Seg_circle_relation(Segment v,Circle C){
	double dst = Dis_point_seg(C.c,v);
	if(sgn(dst - C.r) < 0)return 0;//0:线段在圆内
	if(sgn(dst - C.r) == 0)return 1;//1:线段与圆相切
	return 2;//2:线段在圆外
}

直线和圆的交点

//pa和pb是交点,返回值是交点个数
int Line_cross_circle(Line v,Circle C,Point &pa,Point &pb){
	if(Line_circle_relation(v,C) == 2)return 0;//无交点
	Point q = Point_line_proj(C.c,v);//圆心在直线上的投影点
	double d = Dis_point_line(C.c,v);//圆心到直线的距离
	double k = sqrt(C.r*C.r - d*d);
	if(sgn(k) == 0){
		pa = q,pb = q;return 1;//一个交点,直线和圆相切
	}
	Point n =(v.p2 - v.p1)/Len(v.p2 - v.p1);//单位向量
	pa = q + n*k;pb = q - n*k;
	return 2;//两个交点
}

最小圆覆盖

struct Point{double x,y;};
double Distance(Point A,Point B){return hyhot(A.x - B.x,A.y - B.y);}
Point circle_center(const Point a,const Point b,const Point c){
	Point center;
	double a1 = b.x - a.x,b1 = b.y - a.y,c1 = (a1*a1 + b1*b1)/2;
	double a2 = c.x - a.x,b2 = c.y - a.y,c2 = (a2*a2 + b2*b2)/2;
	double d = a1*b2 - a2*b1;
	center.x = a.x + (c1*b2 - c2*b1)/d;
	center.y = a.y + (a1*c2 - a2*c1)/d;
	return center;
}
void min_cover_circle(Point *p,int n,Point &c,double &r){
	random_shuffle(p,p + n);//随机函数,打乱所有点
	c = p[0];r = 0;//从第一个点p[0]开始,圆心为p[0],半径为0
	for(int i = 1;i < n;i++){
		if(sgn(Distance(p[i],c) - r) > 0){//点p[i]在圆外部
		c = p[i]; r = 0;//重新设置圆心为p[i],半径为0
		for(int j = 0;j < i;j++){//重新检查前面所有点
			if(sgn(Distance(p[j],c)-r) > 0){//两点定圆
			c.x = (p[i].x + p[j].x)/2;
			c.y = (p[i].y + p[j].y)/2;
			r = Distance(p[j],c);
			for(int k = 0;k < j;k++){
				if(sgn(Distance(p[k],c) - r) > 0){
					c = circle_center(p[i],p[j],p[k]);
					r = Distance(p[i],c);
					}
				}
			}
		}
	}
}

以上是二维计算几何中经常会用到的一些模板,当然计算几何问题一直是算法竞赛中的难点,只有理解背后求解的原理并且勤加练习才能更好地掌握这类题目

posted @ 2024-08-03 20:20  isletfall  阅读(57)  评论(0)    收藏  举报