计算几何总结

前置知识

三角函数

如图：

$tan(θ)=AB/BO$
$sin(θ)=AB/AO$
$cos(θ)=BO/AO$

值的正负性：

向量

严格的解释见这里
形象的解释就是一个表示位移的没有固定位置的箭头

一些约定

如没有特别说明单位，角度默认为弧度
如没有特别说明，多边形点的存储方式为逆时针储存。

二维几何

基础

基本定义

struct Point
{
	double x,y;
	Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;

const double Pi=acos(-1);
inline double Angle(const Vector &A)//向量极角
{ return atan2(A.y,A.x); }
inline double R_to_D(double rad)
{ return 180/Pi*rad; }
inline double D_to_R(double D)
{ return Pi/180*D; }

inline Vector Rotate(const Vector &A,double rad)//将向量A逆时针旋转rad
{
	double csr=cos(rad),sir=sin(rad);
	return Vector(A.x*csr-A.y*sir,A.x*sir+A.y*csr);
}

inline Vector Normal(const Vector &A)//单位法向量（单位向量左转90°）
{
	double L=Length(A);
	return Vector(-A.y/L,A.x/L);
}
inline Vector Format(const Vector &A)//单位向量
{
	double L=Length(A);
	return Vector(A.x/L,A.y/L);
}

基本四则运算

inline Vector operator+(const Vector &A,const Vector &B)
{ return Vector(A.x+B.x,A.y+B.y); }

inline Vector operator-(const Point &a,const Point &b)
{ return Vector(a.x-b.x,a.y-b.y); }

inline Vector operator*(const Vector &A,double p)
{ return Vector(A.x*p,A.y*p); }

inline Vector operator/(const Vector &A,double p)
{ return Vector(A.x/p,A.y/p); }

inline bool operator<(const Point &a,const Point &b)
{ return a.x<b.x||(a.x==b.x&&a.y<b.y); }

const double ops=1e-10;
inline int dcmp(double x)
{ return (x>0?x:-x)<ops?0:(x>0?1:-1); }

inline bool operator==(const Point &a,const Point &b)
{ return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; }

点积与叉积

点积定义：两个向量v和w的点积等于二者长度的乘积乘上它们夹角的余弦。
点积性质：

1.满足交换律
2.对于向量加减法满足分配律（$u*(v+w)=u*v+v*w$)
3.见“正负性总结”

叉积定义：两个向量v和w的叉积等于v和w组成的三角形的有向面积的两倍
　　　　如图：
　　　　面积就是红色平行四边形的面积，有向则是指：若w在v的左边，则叉积为正，反之为负。
叉积性质：

1.$cross(v,w)=-cross(w,v)$
2.见“正负性总结”

正负性总结：
代码：

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

inline double Length(const Vector &A)
{ return sqrt(Dot(A,A)); }

inline double Angle(const Vector &A,const Vector &B)//A和B之间的夹角
{ return acos(Dot(A,B)/Length(A)/Length(B)); }

inline double Cross(const Vector &A,const Vector &B)
{ return A.x*B.y-A.y*B.x; }

inline double Area2(const Point &a,const Point &b,const Point &c)
{ return Cross(b-a,c-a); }

两直线交点

设直线分别为$P+tv$和$Q+tw$，设向量$u=P-Q$，则交点在第一条直线的参数为$t_1$，解得：
$\Large t_1=\frac{cross(w,u)}{cross(v,w)}$
即图中蓝色部分面积除以红色部分面积
代码：

inline Point GetLineIntersection
(const Point &A1,const Point &B1,const Point &A2,const Point &B2)
{
	Vector u=A1-A2,v=B1-A1,w=B2-A2;
	double t=Cross(w,u)/Cross(v,w);
	return A1+v*t;
}

点到直线距离

点到直线距离：用平行四边形的面积除以底即可。
点到线段距离：注意用点积判断与从端点引出的线段法向量之间的相对位置
代码：

inline double DistanceToLine
(const Point &P,const Point &A,const Point &B)
{
	Vector v1=B-A,v2=P-A;
	return fabs(Cross(v1,v2))/Length(v1);
}

inline double DistanceToSegment
(const Point &P,const Point &A,const 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);
	else if(dcmp(Dot(v1,v3))>0) return Length(v3);
	else return fabs(Cross(v1,v2))/Length(v1);
}

点在直线上的投影（垂足）

首先，把直线$AB$写成参数式$A+tv$(v=$\overrightarrow{A B}$)，并设投影点$Q$的参数为$t_0$，显然$Q=A+t_0 v$。根据$PQ$垂直于$AB$，两个向量点积为0，即$Dot(v,P-(A+t_0 v))=0$。根据分配律有$Dot(v,P-A)-t_0 *Dot(v,v)=0$，这样就可以解出$\Large t_0 = \frac{Dot(v,P-A)}{Dot(v,v)} $，从而得到$Q$点。
代码：

inline Point GetLineProjection
(const Point &P,const Point &A,const Point &B)
{
	Vector v=B-A;
	return A+v*(Dot(v,P-A)/Dot(v,v));
}

线段相交、点于线段的关系

用叉积判断端点之间的相对位置即可
代码：

inline bool SegmentProperIntersection
(const Point &a1,const Point &a2,const Point &b1,const Point &b2)//必须规范相交
{
	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 dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}

inline bool OnSegment
(const Point &p,const Point &a1,const Point &a2)//不包含端点
{ return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0; }

多边形有向面积

用叉积计算即可。
代码：

inline double PolygonArea(Point* p,int n)//按逆时针储存
{
	double area=0;
	for(int i=1;i<n-1;i++)
		area+=Cross(p[i]-p[0],p[i+1]-p[0]);
	return area/2;
}

圆的基本定义

struct Circle
{
	Point c;
	double r;
	Circle(Point c=Point(),double r=0):c(c),r(r){}
	inline Point point(double a)//通过极角获取圆上一点
	{ return Point(c.x+cos(a)*r,c.y+sin(a)*r); }
};

直线与圆的交点

解个方程组即可。
代码：

inline int GetLineCircleIntersection
(const Point &A,const Point &B,const Circle &C,vector<Point> &sol)
{
	double d=DistanceToLine(C.c,A,B);
	int mode=dcmp(d-C.r);
	if(mode>0) return 0;//相离
	Point P=GetLineProjection(C.c,A,B);
	if(mode==0)//相切 
	{
		sol.push_back(P);
		return 1;
	}
	double dist=sqrt(C.r*C.r-d*d);
	Vector v=Format(B-A);
	sol.push_back(P-v*dist);
	sol.push_back(P+v*dist);
	return 2;
}

圆与的交点

解方程即可。
代码：

inline GetCircleCircleIntersection
(Circle C1,Circle C2,vector<Point> &sol)
{
	if(C1.r<C2.r) swap(C1,C2);//to make sure C1 is bigger than C2
	double D=Length(C1.c-C2.c);
	if(dcmp(D)==0)
		return dcmp(C1.r-C2.r)==0?-1:0;
	if(dcmp(C1.r+C2.r-D)<0) return 0;
	if(dcmp(fabs(C1.r-C2.r)-D)>0) return 0;
	
	double d1=((C1.r*C1.r-C2.r*C2.r)/D+D)/2;
	double x=sqrt(C1.r*C1.r-d1*d1);
	Point O=C1.c+Format(C2.c-C1.c)*d1;
	Point P1=O+Normal(O-C2.c)*x,P2=O-Normal(O-C2.c)*x;
	sol.push_back(P1);
	if(P1==P2) return 1;
	sol.push_back(P2);
	return 2;
}

过圆外一点作圆的切线

细节见代码：

inline int GetTangents
(const Point P,const Circle C,vector<Point> &v)
{
	Vector u=C.c-P;
	double dist=Length(u);
	int mode=dcmp(dist-C.r);
	if(mode<0) return 0;//点在圆内
	if(mode==0)//点在圆上
	{
		v.push_back(P+Normal(u));
		return 1;
	}
	double x=sqrt(dist*dist-C.r*C.r);//解出点到切点的距离
	Circle C2(P,x);//作圆
	return GetCircleCircleIntersection(C,C2,v);//求交点
}

求圆的公切线

详见代码：

inline int _getTangents
(Circle A,Circle B,Point *a,Point *b)
{
	int cnt=0;
	if(A.r<B.r) { swap(A,B); swap(a,b); }
	double d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);
	double rdiff=A.r-B.r;
	double rsum=A.r+B.r;
	if(d2<rdiff*rdiff) return 0;//内含
	
	double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);
	if(d2==0&&A.r==B.r) return -1;//无限多条切线
	if(d2==rdiff*rdiff)//内切，一条公切线
	{
		a[cnt]=A.point(base); b[cnt]=B.point(base); cnt++;
		return 1;
	}
	
	double ang=acos((A.r-B.r)/sqrt(d2));
	a[cnt]=A.point(base+ang); b[cnt]=B.point(base+ang); cnt++;
	a[cnt]=A.point(base-ang); b[cnt]=B.point(base-ang); cnt++;//两条外公切线
	if(d2==rsum*rsum)//两圆相切，一条内公切线
	{
		a[cnt]=A.point(base);b[cnt]=B.point(Pi+base); cnt++;
	}
	else if(d2>rsum*rsum)//两条内公切线
	{
		double ang=acos((A.r+B.r)/sqrt(d2));
		a[cnt]=A.point(base+ang); b[cnt]=B.point(Pi+base+ang); cnt++;
		a[cnt]=A.point(base-ang); b[cnt]=B.point(Pi+base-ang); cnt++;
	}
	return cnt;
}
inline int GetTangents
(const Circle &A,const Circle &B,vector<Point> &va,vector<Point> &vb)//接口
{
	Point a[4],b[4];
	int cnt=_getTangents(A,B,a,b);
	for(int i=0;i<cnt;i++) va.push_back(a[i]);
	for(int i=0;i<cnt;i++) vb.push_back(b[i]);
	return cnt;
}

三角形内接圆

代码：

Circle NeiJieYuan(Point p1,Point p2,Point p3)
{
	double a=Length(p2-p3);
	double b=Length(p3-p1);
	double c=Length(p1-p2);
	Point p=(p1*a+p2*b+p3*c)/(a+b+c);
	return Circle(p,DistanceToLine(p,p1,p2));
}

三角形外接圆

代码：

Circle WaiJieYuan(Point p1,Point p2,Point p3)
{
	double Bx=p2.x-p1.x,By=p2.y-p1.y;
	double Cx=p3.x-p1.x,Cy=p3.y-p1.y;
	double D=2*(Bx*Cy-By*Cx);
	double ansx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;
	double ansy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;
	Point p(ansx,ansy);
	return Circle(p,Length(p1-p));
}

算法

点在多边形内判定

凸多边形：直接判断点是否在每一条边左边（假设逆时针储存每一条边）
普通多边形：假想有一条向右的射线，统计多边形穿过这条射线正反多少次，逆时针穿过时绕数加1，反之减1。
代码：

inline int PointInPolygon(const Point &p,Point *poly,int n)
{
	int wn=0;
	for(int i=0;i<n;i++)
	{
		if(PointOnSegment(p,poly[i],poly[(i+1)%n])||p==poly[i]) return -1;//在边界上
		int k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
		int d1=dcmp(poly[i].y-p.y);
		int d2=dcmp(poly[(i+1)%n].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;
}

凸包

顾名思义，就是求一个包含所有给定点的面积最小的凸多边形。
首先把所有点以x坐标为第一关键字、y坐标为第二关键字排序，删除重复的点后得到序列$p_1,p_2,...,p_n$，然后把$p_1$和$p_2$放入凸包中。从$p_3$开始，当新点在凸包“前进”方向的左边时将其放入凸包，否则依次删除最近放入凸包的点直至新点在左边。
代码：

int ConvexHull(Point *p,int n,Point *ch)
{
	sort(p,p+n);
	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;
}

凸包直径 - 旋转卡壳

凸包直径：凸包中两点间距离的最大值
先求出凸包然后不断地模拟两条刚好卡住凸包的直线的旋转过程，细节见代码：

double Diameter(Point *p,int n)
{
	if(n<=1) return 0;
	if(n==2) return Length(p[0]-p[1]);//特判
	double res=0;
	for(int u=0,v=1;u<n;u++)//u:当前处理地点 v:u对应的边的代表点
	{
		while(true)
		{
			int diff=dcmp(Cross(p[(u+1)%n]-p[u],p[(v+1)%n]-p[v]));
			if(diff<=0)//到达最远点，距离开始变近
			{
				res=max(res,Length(p[u]-p[v]));
				if(diff==0) res=max(res,Length(p[u]-p[(v+1)%n]));//处理一下细节
				break;//处理下一个点
			}
			v=(v+1)%n;//调整
		}
	}
	return res;
}

半平面交

半平面：一条有向直线的左边的部分
半平面交：所有半平面的公共部分
类似于凸包，用双端队列维护，细节详见代码：

struct Line
{
	Point P;
	Vector v;
	double ang;//极角
	Line(Point P=Point(),Vector v=Vector()):P(P),v(v) { ang=atan2(v.y,v.x); }
	inline bool operator<(const Line &L) const
	{ return ang<L.ang; }
};

inline bool OnLeft(const Line &L,const Point &p)//判断点是否在有向直线的左边
{ return dcmp(Cross(L.v,p-L.P))>0; }

Point GetLineIntersection(const Line &a,const 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;
}

int HalfplaneIntersection(Line *L,int n,Point *poly)
{
	sort(L,L+n);//按极角排序
	int first,last;//双端队列 [first,last]
	Point *p=new Point[n];//p[i]为q[i]和q[i+1]的交点
	Line *q=new Line[n];//双端队列
	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(dcmp(Cross(q[last].v,q[last-1].v))==0)
		{//两向量平行且相同，取内侧的一个
			last--;
			if(OnLeft(q[last],L[i].P)) q[last]=L[i];
		}
		if(first<last) p[last-1]=GetLineIntersection(q[last-1],q[last]);
	}
	while(first<last&&!OnLeft(q[first],p[last-1])) last--;//删除无用平面
	if(last-first<=1) return 0;//空集
	p[last]=GetLineIntersection(q[last],q[first]);
	
	int m=0;
	for(int i=first;i<=last;i++) poly[m++]=p[i];
	return m;
}

参考资料：

刘汝佳《算法竞赛入门经典训练指南》
网络（太多了，贴不下。。。）

posted @ 2019-08-19 21:03 happyZYM 阅读(357) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

happyZYM的博客

前进四