计算几何学习笔记

一、常量&约定

为了方便地比较浮点数，我们需要设计一个三态比较函数，如下。

#define IL inline 
const double pi=acos(-1.0);
const double eps=1e-8;
IL int sign(double x){
	if(fabs(x)<eps)return 0;
	return (x<0)?-1:1;
}
IL int dcmp(double x,double y){return sign(x-y);}

另外，除非特殊说明，所有涉及旋转/角度的操作，都默认逆时针为正。

二、点&向量

2.1 加、减、数乘

在计算几何中，由于三角函数运算的值不够精确，我们一般采用坐标运算，我们首先需要重载最基本的坐标运算。

struct Point{
	double x,y;
	Point(double X=0,double Y=0):x(X),y(Y){}
	inline void in(){scanf("%lf%lf",&x,&y);}
	inline void out(){printf("%.2lf %.2lf
",x,y);}
	Point operator + (const Point &r) const{return Point(x+r.x,y+r.y);}
	Point operator - (const Point &r) const{return Point(x-r.x,y-r.y);}
	Point operator * (double p) const{return Point(x*p,y*p);}
	Point operator / (double p) const{return Point(x/p,y/p);}
	bool operator < (const Point &r) const{return x<r.x||(x==r.x&&y<r.y);}
	bool operator == (const Point &r) const{return sign(x-r.x)==0&&sign(y-r.y)==0;}
};
typedef Point Vector;

2.2 点积

(oldsymbol{a}cdot oldsymbol{b}=|oldsymbol{a}||oldsymbol{b}|cos heta( heta=leftlangle oldsymbol{a},oldsymbol{b} ight angle)=x_1 x_2 +y_1 y_2)

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

利用点积，我们可以比较方便地求出向量的模长，两点间的距离，向量的夹角，以及一个向量在另一个向量上的投影。

IL double length(Vector A){return sqrt(dot(A,A));}
IL double dis(Point A,Point B){return length(A-B);}
IL double angle(Vector A,Vector B){return acos(dot(A,B)/length(A)/length(B));}
IL double project(Point a,Point b,Point c){return dot(b-a,c-a)/length(b-a);}

2.3 叉积

(oldsymbol{a} imes oldsymbol{b}=|oldsymbol{a}||oldsymbol{b}|sin heta( heta=leftlangle oldsymbol{a},oldsymbol{b} ight angle)=x_1 y_2 -x_2 y_1)

IL double cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}

叉积的结果可以表示两个向量(oldsymbol a)，(oldsymbol b)共起点时，所构成平行四边形的有向面积，这个性质会在之后的许多操作中得以应用。此外，有向面积的正负可以帮助我们判断两个向量共起点以后的位置关系。

具体而言，如果(oldsymbol a)在(oldsymbol b)左侧，那么(oldsymbol{a} imes oldsymbol{b}<0)

IL double area(Point a,Point b,Point c){return cross(b-a,c-a);}

2.4 极角

利用C++函数库中的(mathtt{atan2(y,x)})可以帮助我们计算(arctan frac{y}{x})的值，其值域为((-pi,pi])，即当点位于三四象限时其返回一个负值。

IL double polar_angle(Vector A){return atan2(A.y,A.x);}

2.5 旋转

[egin{bmatrix}x yend{bmatrix} imesegin{bmatrix}cos heta & sin heta\ -sin heta &cos hetaend{bmatrix} ]

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

三、直线&线段

3.1 直线的表示&参数

我们常用点向式（直线的参数方程）来表示一个直线：

[l:a+oldsymbol{v}t ]

为了方便，我们提供两种构造直线的方式（传入两点，或者传入一点和一个向量），并根据需要记录(a,b,oldsymbol{v},r)(极角)。

struct Line{
	Point a,b;
	Vector v;
	double r;
	Line(){} 
	Line(const Point &A,const Point &B,int op=0){
		if(!op){a=A,b=B;v=b-a;r=polar_angle(v);}
		else a=A,v=B;b=a+v;r=polar_angle(v);
	}
	Point point(double p){return a+v*p;}
};

3.2 点与直线（线段）的位置关系

如果一个点(p)满足(vec{ap} imes oldsymbol{v}=0)，我们就可以判定(p)在直线上。

如果我们需要判定点在线段上，还需要额外满足(vec{pa}cdot vec{pb}<0)

bool line_inc(const Point &p) const{
   return sign(cross(p-a,v))==0;
}
bool seg_inc(const Point &p) const{
   return sign(cross(a-p,b-p))==0&&sign(dot(a-p,b-p))<=0;
}

3.3 点在直线上的投影点(垂足)

用2.2中提到的方法计算，记得考虑要使用(oldsymbol{v})对应的单位向量。

IL Point point_line_proj(Point p,Line l){
	return l.a+l.v*(dot(l.v,p-l.a)/dot(l.v,l.v));
}

3.4 点到直线（线段）的距离

我们采取面积/底边的方式计算点到直线的距离（高）。

如果要计算点到线段的距离，需要先判断高是否可以落在线段上，如果不在，应该取返回点与最近端点的距离。

IL double point_to_line(Point p,Line l){
	return fabs(cross(l.v,p-l.a))/length(l.v);
}
IL double point_to_seg(Point p,Seg s){
	if(s.a==s.b)return length(p-s.a);
	Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.b;
	if(sign(dot(v1,v2))<0)return length(v2);
	if(sign(dot(v1,v3))>0)return length(v3);
	return point_to_line(p,Line(s.a,s.b)); 
}

3.5 直线（线段）间的位置关系

我们利用面积比来得到直线与直线的交点。

值得注意的是，为了防止被除数为零，我们需要预先判断两直线是否平行。

IL Point line_line_inter(Line x,Line y){//记得先判平行 
	Vector U=x.a-y.a;
	double t=cross(y.v,U)/cross(x.v,y.v);//注意方向 
	return x.a+x.v*t;
}

如果两线段相交，则两线段必然相互跨立对方。我们通过跨立试验来判断，两线段是否相交（直线与线段是否相交类似）。我们同样需要预先判断两直线是否平行。

IL bool seg_seg_inter(Seg x,Seg y){
	Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
	double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),
		c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
	return sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0;
}
IL bool line_seg_inter(Line x,Seg y){
	Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
	double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
	return sign(c1)*sign(c2)<=0;
}

四、圆

4.1 圆的表示&参数

我们可以使用圆心和半径表示一个圆。

struct Circle{
	Point o;double r;
	Circle(){}
	Circle(const Point &O,double R):o(O),r(R){}
}

4.2 点与圆的位置关系

计算点与圆心的距离然后判断即可。

bool inc(const Point &p){
	return dcmp(r,dis(o,p))>=0;
}

4.3 求三角形的外接圆

以任意两边中垂线的交点为圆心，交点到一点处的距离为半径作圆。

IL Line get_line(Point a,Point b){return Line((a+b)/2,rotate(b-a,pi/2),1);};
IL Circle get_circle(Point a,Point b,Point c){
	Line u=get_line(a,b),v=get_line(a,c);
	Point o=line_line_inter(u,v);
	return Circle(o,dis(o,a));
}

4.4 最小圆覆盖

如果点(p)不在点集(S)的最小覆盖圆内，则(p)一定在({S}cup{p}) 的最小覆盖圆上。

利用这个定理，我们可以在(O(n))时间复杂度内求出一个点集的最小覆盖圆。

Circle min_circle(Point p[],int n){
	random_shuffle(p,p+n);
	Circle c=Circle(p[0],0);
	for(int i=1;i<n;i++)
		if(!c.inc(p[i])){
			c=Circle(p[i],0);
			for(int j=0;j<i;j++)
				if(!c.inc(p[j])){
					c=Circle((p[i]+p[j])/2,dis(p[i],p[j])/2);
					for(int k=0;k<j;k++)
						if(!c.inc(p[k]))
							c=get_circle(p[i],p[j],p[k]);
				}
		}
	return c;
}

4.5 伸缩变换

对于椭圆问题，我们可以对整个坐标系进行伸缩变换，再进行计算。

五、普通多边形

5.1 多边形的面积

我们用三角剖分求一个一般多边形的面积。

double poly_area(Point p[],int n){
	double res=0;
	for(int i=1;i<n-1;i++) 
		res+=cross(p[i]-p[0],p[i+1]-p[0]);
	return res/2;
}

5.2 点与多边形的位置关系

我们采用射线法：考虑以该点为端点引出一条射线，如果这条射线与多边形有奇数个交点，则该点在多边形内部，否则该点在多边形外部，这被称为奇偶规则。

六、凸多边形

6.1 Andrew算法求凸包

我们采用Andrew算法可以(O(n))地计算一个点集的凸包。

int top,st[N];
bool used[N];
void Andrew(Point p[],int cnt){
	sort(p+1,p+1+cnt);
	for(int i=1;i<=cnt;i++){
		while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)used[st[top--]]=0;
		st[++top]=i;used[i]=1;
	}
	used[1]=0;
	for(int i=cnt;i;i--){
		if(used[i])continue;
		while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)top--;
		st[++top]=i;
	}
}

6.2 点与凸多边形的位置关系

只要一个点在所有边的左侧，点就在凸包内，这可以遍历一遍实现或者用二分实现。

6.3 旋转卡壳

旋转卡壳可以用于求凸包的直径。

double rotating_caliper(Point p[]){
	double res=0;
	for(int i=1,a=3;i<top;i++){
		Point d=p[st[i]],e=p[st[i+1]];
		while(dcmp(area(d,e,p[st[a]]),area(d,e,p[st[a+1]]))<0)a=a%(top-1)+1;
		res=max(res,max(dis(d,p[st[a]]),dis(e,p[st[a]])));
	}
	return res;
}

至于求平面最近点对，这不是个计算几何问题，但也放在这里，这个问题除了放在下面的解法外，还有分治的解法，可以去OI-wiki上查看。

struct cmp_y{
  bool operator()(const Point &a, const Point &b)const{return a.y<b.y;}
};
multiset<Point,cmp_y>s;
double GetMin(){
	double res=1e16;
	for(int i=1,l=1;i<=n;i++){
		while(l<i&&dcmp(p[i].x-p[l].x,res)>=0)s.erase(s.find(p[l++]));
		for(auto it=s.lower_bound(Point(0,p[i].y-res));
			it!=s.end()&&dcmp(it->y-p[i].y,res)<0;it++)
				if(dcmp(dis(*it,p[i]),res)<0)res=dis(*it,p[i]);
		s.insert(p[i]);
	}
	return res;
}

另外，平面最近最远点对问题有一个随机化的贪心算法：将整个坐标系随机旋转一个角度，再贪心计算。

6.4 半平面交

计算每条直线左侧的平面的交集。

bool cmp(const Line& x,const Line& y){
	if(sign(x.r-y.r)==0)return area(x.a,x.b,y.b)<0; 
	return x.r<y.r;
}
IL bool on_right(Line a,Line b,Line c){
	Point o=line_line_inter(b,c);
	return sign(area(a.a,a.b,o))<=0;
}
double half_plane_inter(Line l[],int cnt){
	sort(l+1,l+1+cnt,cmp);
	int L=0,R=-1;
	vector<int>q(cnt);
	for(int i=1;i<=cnt;i++){
		if(i>1&&sign(l[i].r-l[i-1].r)==0)continue;
		while(L+1<=R&&on_right(l[i],l[q[R-1]],l[q[R]]))R--;
		while(L+1<=R&&on_right(l[i],l[q[L]],l[q[L+1]]))L++;
		q[++R]=i;
	}
	while(L+1<=R&&on_right(l[q[L]],l[q[R-1]],l[q[R]]))R--;
	while(L+1<=R&&on_right(l[q[R]],l[q[L]],l[q[L+1]]))L++;
	q[++R]=q[L];
	vector<Point>ans(R);
	int k=0;
	for(int i=L;i<R;i++)
		ans[++k]=line_line_inter(l[q[i]],l[q[i+1]]);
	double res=0;
	for(int i=2;i+1<=k;i++)res+=area(ans[1],ans[i],ans[i+1]);
	return res/2;
}

七、面积并的计算

7.1 扫描线

计算图形的面积并时，一种常用的思路是运用扫描线——将一维取出排序，然后挨个统计对应的另一维的答案。

7.2 自适应辛普森法

总体思路是根据辛普森公式不断地近似计算函数，直到得到精度合适的值。

辛普森公式：

[int_{a}^{b} f(x) dx approx frac{b-a}{6}[f(a)+4f(frac{a+b}{2})+f(b)] ]

代码模板如下：

IL double f(double x){
}
IL double simpson(double a,double b){
	return ((b-a)*(f(a)+f(b)+4*f((a+b)/2)))/6;
}
double divide(double L,double R,double ans){
	double mid=(L+R)/2,l=simpson(L,mid),r=simpson(mid,R);
	if(dcmp(l+r,ans)==0)return ans;
	return divide(L,mid,l)+divide(mid,R,r);
}

如果本题的(f (x))计算较慢，需要根据题目需要将已经计算过的函数值作为参数传下去，避免重复计算。

如果精度不好保证，还需要将精度不断衰减保证效率。

附 :【模板】计算几何全家桶

#include<bits/stdc++.h>
using namespace std;
 
#define IL inline 
const int N=1e5+5;
const double pi=acos(-1.0);
const double eps=1e-8;
IL int sign(double x){
	if(fabs(x)<eps)return 0;
	return (x<0)?-1:1;
}
IL int dcmp(double x,double y){return sign(x-y);}

struct Point{
	double x,y;
	Point(double X=0,double Y=0):x(X),y(Y){}
	inline void in(){scanf("%lf%lf",&x,&y);}
	inline void out(){printf("%.2lf %.2lf
",x,y);}
	Point operator + (const Point &r) const{return Point(x+r.x,y+r.y);}
	Point operator - (const Point &r) const{return Point(x-r.x,y-r.y);}
	Point operator * (double p) const{return Point(x*p,y*p);}
	Point operator / (double p) const{return Point(x/p,y/p);}
	bool operator < (const Point &r) const{return x<r.x||(x==r.x&&y<r.y);}
	bool operator == (const Point &r) const{return sign(x-r.x)==0&&sign(y-r.y)==0;}
};
typedef Point Vector;
IL double dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;} 
IL double cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
IL double length(Vector A){return sqrt(dot(A,A));}
IL double dis(Point A,Point B){return length(A-B);}
IL double angle(Vector A,Vector B){return acos(dot(A,B)/length(A)/length(B));}
IL double project(Point a,Point b,Point c){return dot(b-a,c-a)/length(b-a);}
IL double polar_angle(Vector A){return atan2(A.y,A.x);}
IL double area(Point a,Point b,Point c){return cross(b-a,c-a);}
IL Vector rotate(Vector A,double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

struct Line{
	Point a,b;
	Vector v;
	double r;
	Line(){} 
	Line(const Point &A,const Point &B,int op=0){
		if(!op){a=A,b=B;v=b-a;r=polar_angle(v);}
		else a=A,v=B;b=a+v;r=polar_angle(v);
	}
	Point point(double p){return a+v*p;}
	bool line_inc(const Point &p) const{return sign(cross(p-a,v))==0;}
	bool seg_inc(const Point &p) const{
		return sign(cross(a-p,b-p))==0&&sign(dot(a-p,b-p))<=0;
	}
};
typedef Line Seg;
IL Point point_line_proj(Point p,Line l){
	return l.a+l.v*(dot(l.v,p-l.a)/dot(l.v,l.v));
}
IL double point_to_line(Point p,Line l){
	return fabs(cross(l.v,p-l.a))/length(l.v);
}
IL double point_to_seg(Point p,Seg s){
	if(s.a==s.b)return length(p-s.a);
	Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.b;
	if(sign(dot(v1,v2))<0)return length(v2);
	if(sign(dot(v1,v3))>0)return length(v3);
	return point_to_line(p,Line(s.a,s.b)); 
}
IL Point line_line_inter(Line x,Line y){//记得先判平行 
	Vector U=x.a-y.a;
	double t=cross(y.v,U)/cross(x.v,y.v);//注意方向 
	return x.a+x.v*t;
}
IL bool seg_seg_inter(Seg x,Seg y){
	Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
    double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),
           c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
    return sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0;
}
IL bool line_seg_inter(Line x,Seg y){
	Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
    double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
    return sign(c1)*sign(c2)<=0;
}

struct Circle{
	Point o;double r;
	Circle(){}
	Circle(const Point &O,double R):o(O),r(R){}
	bool inc(const Point &p){
		return dcmp(r,dis(o,p))>=0;
	}
};
IL Line get_line(Point a,Point b){return Line((a+b)/2,rotate(b-a,pi/2),1);};
IL Circle get_circle(Point a,Point b,Point c){
	Line u=get_line(a,b),v=get_line(a,c);
	Point o=line_line_inter(u,v);
	return Circle(o,dis(o,a));
}
Circle min_circle(Point p[],int n){
	random_shuffle(p,p+n);
	Circle c=Circle(p[0],0);
	for(int i=1;i<n;i++)
		if(!c.inc(p[i])){
			c=Circle(p[i],0);
			for(int j=0;j<i;j++)
				if(!c.inc(p[j])){
					c=Circle((p[i]+p[j])/2,dis(p[i],p[j])/2);
					for(int k=0;k<j;k++)
						if(!c.inc(p[k]))
							c=get_circle(p[i],p[j],p[k]);
				}
		}
	return c;
}

double poly_area(Point p[],int n){
    double res=0;
    for(int i=1;i<n-1;i++) 
        res+=cross(p[i]-p[0],p[i+1]-p[0]);
    return res/2;
}

int top,st[N];
void Andrew(Point p[],int cnt){
	vector<bool>used(cnt);
	sort(p+1,p+1+cnt);
	for(int i=1;i<=cnt;i++){
		while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)used[st[top--]]=0;
		st[++top]=i;used[i]=1;
	}
	used[1]=0;
	for(int i=cnt;i;i--){
		if(used[i])continue;
		while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)top--;
		st[++top]=i;
	}
}

double rotating_caliper(Point p[]){
	double res=0;
	for(int i=1,a=3;i<top;i++){
		Point d=p[st[i]],e=p[st[i+1]];
		while(dcmp(area(d,e,p[st[a]]),area(d,e,p[st[a+1]]))<0)a=a%(top-1)+1;
		res=max(res,max(dis(d,p[st[a]]),dis(e,p[st[a]])));
	}
	return res;
}

bool cmp(const Line& x,const Line& y){
	if(sign(x.r-y.r)==0)return area(x.a,x.b,y.b)<0; 
	return x.r<y.r;
}
IL bool on_right(Line a,Line b,Line c){
	Point o=line_line_inter(b,c);
	return sign(area(a.a,a.b,o))<=0;
}
double half_plane_inter(Line l[],int cnt){
	sort(l+1,l+1+cnt,cmp);
	int L=0,R=-1;
	vector<int>q(cnt);
	for(int i=1;i<=cnt;i++){
		if(i>1&&sign(l[i].r-l[i-1].r)==0)continue;
		while(L+1<=R&&on_right(l[i],l[q[R-1]],l[q[R]]))R--;
		while(L+1<=R&&on_right(l[i],l[q[L]],l[q[L+1]]))L++;
		q[++R]=i;
	}
	while(L+1<=R&&on_right(l[q[L]],l[q[R-1]],l[q[R]]))R--;
	while(L+1<=R&&on_right(l[q[R]],l[q[L]],l[q[L+1]]))L++;
	q[++R]=q[L];
	vector<Point>ans(R);
	int k=0;
	for(int i=L;i<R;i++)
		ans[++k]=line_line_inter(l[q[i]],l[q[i+1]]);
	double res=0;
	for(int i=2;i+1<=k;i++)res+=area(ans[1],ans[i],ans[i+1]);
	return res/2;
}