计算几何总结

终于开坑了,慢慢来QAQ

其实就是对着蓝书刷

由于博主菜爆了,半平面交和平面区域可能会咕咕咕。

模板

//有一些公式可以画图证,很容易的。
#include<vector>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define gt getchar()
#define ll long long
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
inline int in()
{
	int k=0;char ch=gt;
	while(ch<'-')ch=gt;
	while(ch>'-')k=k*10+ch-'0',ch=gt;
	return k;
}
using std::vector;
const double eps=1e-10,PI=acos(-1);const int N=605;
inline int dcmp(double x){return fabs(x)<eps?0:(x>0?1:-1);}
inline double torad(double deg){return deg/180*PI;}
inline double tojd(double rad){return rad*180/PI;}

struct Point
{
	double x,y;
	Point(double _x=0,double _y=0):x(_x),y(_y){}
	inline void out(char c=' '){printf("%.5lf %.5lf%c",x,y,c);}
};
inline Point Point_in(){static double x,y;scanf("%lf%lf",&x,&y);return Point(x,y);}
typedef Point Vec;
inline Vec operator+(Vec a,Vec b){return Vec(a.x+b.x,a.y+b.y);}
inline Vec operator-(Point a,Point b){return Vec(a.x-b.x,a.y-b.y);}
inline Vec operator*(Vec a,double p){return Vec(a.x*p,a.y*p);}
inline Vec operator/(Vec a,double p){return Vec(a.x/p,a.y/p);}
inline Vec operator-(Vec a){return a*(-1);}
inline bool operator<(const Point &a,const Point &b){return dcmp(a.x-b.x)<0||(dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)<0);}
inline bool operator==(const Point &a,const Point &b){return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}
inline bool operator!=(const Point &a,const Point &b){return a==b?0:1;}
inline double ang(Vec a){return atan2(a.y,a.x);}//极角:x轴旋转到该向量的角度
inline double Dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
inline double Len(Vec a){return sqrt(Dot(a,a));}
inline double ang(Vec a,Vec b){return acos(Dot(a,b)/Len(a)/Len(b));}
inline double Cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
inline double Area2(Point a,Point b,Point c){return Cross(b-a,c-a);}//面积的两倍不是平方
inline Vec Rotate(Vec a,double rad){return Vec(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}
inline Vec Normal(Vec a){double L=Len(a);return Vec(-a.y/L,a.x/L);}//使用前确保不是零向量

struct Line
{
	Point p;Vec v;
	Line(){}
	Line(Point a,Point b){p=a;v=b-a;}
	inline Point get(double t){return p+v*t;}
};
inline Point Line_jd(Line a,Line b)
{
	Vec u=a.p-b.p;
	double t=Cross(b.v,u)/Cross(a.v,b.v);
	return a.p+a.v*t;	
}
inline double dis_Point_Line(Point a,Line b){return fabs(Cross(b.v,a-b.p))/Len(b.v);}
inline double dis_Point_Line(Point P,Point a,Point b){return dis_Point_Line(P,Line(a,b));}
inline double dis_Point_Seg(Point P,Point a,Point b)
{
	if(a==b)return Len(P-a);
	Vec v1=b-a,v2=P-a,v3=P-b;
	if(dcmp(Dot(v1,v2))<0)return Len(v2);
	if(dcmp(Dot(v1,v3))>0)return Len(v3);
	return fabs(Cross(v1,v2))/Len(v1);
}
inline Point touying_Point_Line(Point a,Line b){return b.p+b.v*(Dot(b.v,a-b.p)/Dot(b.v,b.v));}
inline Point touying_Point_Line(Point P,Point a,Point b){return touying_Point_Line(P,Line(a,b));}
inline bool OnSeg(Point P,Point a,Point b){return dcmp(Cross(P-a,P-b))==0&&dcmp(Dot(P-a,P-b))<0;}
//opt=1表示可以在端点处相交,但不能是端点和端点相交
inline bool Seg_jiao_Seg(Point a1,Point a2,Point b1,Point b2,int opt=0)
{
	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);
	int d1=dcmp(c1),d2=dcmp(c2),d3=dcmp(c3),d4=dcmp(c4);
	if(opt)
	{
		if(d1==0&&d2==0)return OnSeg(b1,a1,a2)^OnSeg(b2,a1,a2);
		if(d1==0)return OnSeg(b1,a1,a2);
		if(d2==0)return OnSeg(b2,a1,a2);
	}
	return d1*d2<0&&d3*d4<0;
}

inline double PolygonArea(Point *p,int n)
{
	double area=0;
	for(int i=2;i<n;++i)
		area+=Cross(p[i]-p[1],p[i+1]-p[1]);
	return area/2;	
}
inline double abs_PolygonArea(Point *p,int n){return fabs(PolygonArea(p,n));}
struct Polygon
{
	Point p[N];int n;
	Polygon(){}
	Polygon(Point *a,int n):n(n){for(int i=1;i<=n;++i)p[i]=a[i];}
	inline double Area(){return PolygonArea(p,n);}
	inline double abs_Area(){return fabs(Area());}
	inline int Point_in(Point P)
		{
			int wn=0;
			for(int i=1;i<=n;++i)
			{
				int j=i==n?1:i+1;
				if(OnSeg(P,p[i],p[j]))return -1;
				int k=dcmp(Cross(p[j]-p[i],P-p[i]));
				int d1=dcmp(p[i].y-P.y),d2=dcmp(p[j].y-P.y);
				if(k>0&&d1<=0&&d2>0)++wn;
				if(k<0&&d2<=0&&d1>0)--wn;
			}
			return wn!=0;
		}
};

struct Circle
{
	Point c;double r;
	Circle(){}
	Circle(Point c,double r):c(c),r(r){}
	inline Point get(double a){return Point(c.x+cos(a)*r,c.y+sin(a)*r);}//传进来极角,传出点。
};
inline Circle Circle_in(){double x,y,r;scanf("%lf%lf%lf",&x,&y,&r);return Circle(Point(x,y),r);}
//返回交点个数,交点存在sol里,t1和t2是解在直线上的参数
int Line_Circle_jd(Line L,Circle C,double &t1,double &t2,vector<Point>& sol)
{
	double a=L.v.x,b=L.p.x-C.c.x,c=L.v.y,d=L.p.y-C.c.y;
	double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-C.r*C.r,delta=f*f-4*e*g;
	int t=dcmp(delta);if(t<0)return 0;if(t==0){return t1=t2=-f/(e*2),sol.push_back(L.get(t1)),1;}
	delta=sqrt(delta);t1=(-f-delta)/(e*2),t2=(-f+delta)/(e*2);
	sol.push_back(L.get(t1)),sol.push_back(L.get(t2));return 2;
}
//-1无限多交点0没有1,2是个数,b是用余弦定理算的
int Circle_Circle_jd(Circle C1,Circle C2,vector<Point>& sol)
{
	double d=Len(C1.c-C2.c);if(dcmp(d)==0)return dcmp(C1.r-C2.r)?0:-1;
	if(dcmp(C1.r+C2.r-d)<0||dcmp(fabs(C1.r-C2.r)-d)>0)return 0;
	double a=ang(C2.c-C1.c),b=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
	Point p1=C1.get(a-b),p2=C1.get(a+b);sol.push_back(p1);
	if(p1!=p2)sol.push_back(p2);return 2-(p1==p2);
}
//返回的是第一个圆上的弧度
int Circle_Circle_jd(Circle C1,Circle C2,vector<double>& sol)
{
	double d=Len(C1.c-C2.c);if(dcmp(d)==0)return dcmp(C1.r-C2.r)?0:-1;
	if(dcmp(C1.r+C2.r-d)<0||dcmp(fabs(C1.r-C2.r)-d)>0)return 0;
	double a=ang(C2.c-C1.c),b=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
	double p1=a-b,p2=a+b;sol.push_back(p1);
	if(dcmp(p1-p2)!=0)sol.push_back(p2);return 2-(dcmp(p1-p2)==0);
}
//返回切线条数
int get_qie(Point p,Circle C,vector<Point>&sol)
{
	Vec u=C.c-p;long double dis=Len(u);if(dcmp(dis-C.r)<0)return 0;
	if(dcmp(dis-C.r)==0)return sol.push_back(p),1;
	double a=asin(C.r/dis);dis=sqrt(dis*dis-C.r*C.r);
	Vec t1=Rotate(u,-a);t1=t1/Len(t1)*dis;sol.push_back(p+t1);
	Vec t2=Rotate(u, a);t2=t2/Len(t2)*dis;sol.push_back(p+t2);
	return 2;
}
int get_qie(Circle A,Circle B,vector<Point> *a,vector<Point> *b)
{
	int cnt=0;if(A.r<B.r)std::swap(A,B),std::swap(a,b);
	double d=Len(A.c-B.c),rdel=A.r-B.r,rsum=A.r+B.r;
	if(dcmp(d-rdel)<0)return 0;double base=ang(B.c-A.c);
	if(dcmp(d)==0&&dcmp(A.r-B.r)==0)return -1;
	if(dcmp(d-rdel)==0)
		return a->push_back(A.get(base)),b->push_back(B.get(base)),1;
	double ta=acos((A.r-B.r)/d);
	a->push_back(A.get(base+ta)),b->push_back(B.get(base+ta));++cnt;
	a->push_back(A.get(base-ta)),b->push_back(B.get(base-ta));++cnt;
	if(dcmp(d-rsum)==0)
		return a->push_back(A.get(base)),b->push_back(B.get(PI+base)),3;
	if(dcmp(d-rsum)>0)
	{
		double ta2=acos(rsum/d);
		a->push_back(A.get(base+ta2)),b->push_back(B.get(PI+base+ta2));++cnt;
		a->push_back(A.get(base-ta2)),b->push_back(B.get(PI+base-ta2));++cnt;	
	}
	return cnt;
}

//球面的经纬度转化为关于圆心的坐标
inline void get_coord(double R,double lat,double lng,double &x,double &y,double &z)
{
	lat=torad(lat),lng=torad(lng);
	x=R*cos(lat)*cos(lng);
	y=R*cos(lat)*sin(lng);
	z=R*sin(lat);
}
inline double get_dis_on_ball(double R,double t1,double g1,double t2,double g2)
{
	double x1,y1,z1,x2,y2,z2;
	get_coord(R,t1,g1,x1,y1,z1);
	get_coord(R,t2,g2,x2,y2,z2);
	double d=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
	return 2*asin(d/(2*R))*R;
}

inline int get_tb(Point *p,int n,Point *ch)
{
	++p,++ch;std::sort(p,p+n);int m=0,k;
	for(int i=0; i<n;ch[m++]=p[i],++i)
		while(m>1&&dcmp(Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<=0)--m;
	k=m;
	for(int i=n-2;~i;ch[m++]=p[i],--i)
		while(m>k&&dcmp(Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<=0)--m;
	if(n>1)--m;return m;
}
int main()
{
	return 0;
}



一些题

UVA11178

直接模拟即可。

核心代码

inline void Getans(Point a,Point b,Point c)
{
	Vec v1=b-a,v2=c-a,v=c-b;double t=ang(v1,v2),t2=ang(-v1,v);
	Vec v3=Rotate(v1,t/3),v4=Rotate(v,t2*2/3);
	Line l1(a,a+v3),l2(b,b+v4);(Line_jd(l1,l2)).out();
}
inline void work()
{
	Point a=Point_in(),b=Point_in(),c=Point_in();
	Getans(b,c,a),Getans(c,a,b),Getans(a,b,c);puts("");
}

LA3263

欧拉定理:平面图顶点数(V)、边数(E)、面数(F)满足(V+F-E=2)
所以(F=E+2-V)
核心代码

int n,c,e,ca=0;
	while(n=in())
	{
		for(int i=1;i<=n;++i)v[i]=p[i]=Point_in();
		--n;c=n,e=n;
		for(int i=1;i<=n;++i)
			for(int j=i+1;j<=n;++j)
				if(Seg_jiao_Seg(p[i],p[i+1],p[j],p[j+1],1))
					v[++c]=Line_jd(Line(p[i],p[i+1]),Line(p[j],p[j+1]));
		std::sort(v+1,v+c+1);c=std::unique(v+1,v+c+1)-v-1;
		for(int i=1;i<=c;++i)
			for(int j=1;j<=n;++j)
				if(OnSeg(v[i],p[j],p[j+1]))++e;
		printf("Case %d: There are %d pieces.
",++ca,e+2-c);
	}

UVA11796

先考虑只在一条线段上走,那么可以把一个的位移向量减到另一个上,就可以视为这个点不动,另外那个是一条线段。
折线的化就把它当很多段,一段一段地走就ok辣。注意选用合适的速度,这里用(frac 1{Len})来表示。
核心代码

Point p[N],q[N];
double mx,mi;
inline void Getans(Point P,Point a,Point b)
{
	mi=std::min(mi,dis_Point_Seg(P,a,b));
	mx=std::max(mx,std::max(Len(P-a),Len(P-b)));
}
inline void work(int ca)
{
	int n=in(),m=in();double LenA=0,LenB=0;
	for(int i=1;i<=n;++i)p[i]=Point_in();
	for(int i=1;i<=m;++i)q[i]=Point_in();
	for(int i=1;i<n;++i)LenA+=Len(p[i+1]-p[i]);
	for(int i=1;i<m;++i)LenB+=Len(q[i+1]-q[i]);
	mx=-1e9,mi=1e9;Point A=p[1],B=q[1];int An=1,Bn=1;
	while(An<n&&Bn<m)
	{
		double La=Len(p[An+1]-A),Lb=Len(q[Bn+1]-B);
		double T=std::min(La/LenA,Lb/LenB);
		Vec Va=(p[An+1]-A)/La*T*LenA,Vb=(q[Bn+1]-B)/Lb*T*LenB;
		Getans(A,B,B+Vb-Va);A=A+Va,B=B+Vb;
		if(A==p[An+1])++An;if(B==q[Bn+1])++Bn;
	}
	printf("Case %d: %.0lf
",ca,mx-mi);
}

UVA10674

这道题精度卡的特别紧,我一开始比较没有用dcmp挂烂了。。。
核心代码

int id[N];vector<Point>A,B;
inline bool cmp(const int &x,const int &y){return A[x]<A[y];}
int main()
{
	double a,b,c,d,e,f;
	while(~scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f))
	{
		if(!dcmp(a)&&!dcmp(b)&&!dcmp(c)&&!dcmp(d)&&!dcmp(e)&&!dcmp(f))break;
		Circle C1(Point(a,b),c),C2(Point(d,e),f);A.clear(),B.clear();
		int ans=get_qie(C1,C2,&A,&B);printf("%d
",ans);
		if(ans>0)
		{
			for(int i=1;i<=ans;++i)id[i]=i-1;
			std::sort(id+1,id+ans+1,cmp);
			for(int i=1;i<=ans;++i)
				A[id[i]].out(),B[id[i]].out(),printf("%.5lf
",Len(B[id[i]]-A[id[i]]));
		}
	}
	return 0;
}

UVA12304

超恶心的计算几何六合一,但有了模板就不难了。
第四个子问题要特判圆到直线距离,要不然就嘿嘿嘿了。
前面两个太弱智了,第三个也很简单,注意输出格式。
四、五、六分别可以转化为圆与直线相交,直线与直线相交,圆与圆相交的问题。
核心代码

char str[155];
void work1()
{
	Point p1=Point_in(),p2=Point_in(),p3=Point_in();
	Point A=p1+(p2-p1)/2;Vec v1=Normal(p2-p1);Line a(A,A+v1);
	Point B=p1+(p3-p1)/2;Vec v2=Normal(p3-p1);Line b(B,B+v2);
	Point C=Line_jd(a,b);double dis=Len(C-p1);
	printf("(%lf,%lf,%lf)
",C.x,C.y,dis);
}
void work2()
{
	Point p1=Point_in(),p2=Point_in(),p3=Point_in();
	double a=Len(p2-p3),b=Len(p3-p1),c=Len(p1-p2);
	Point P=(p1*a+p2*b+p3*c)/(a+b+c);double dis=dis_Point_Line(P,p1,p2);
	printf("(%lf,%lf,%lf)
",P.x,P.y,dis);
}
void work3()
{
	Point c=Point_in();double r;scanf("%lf",&r);Point P=Point_in();
	Circle t(c,r);vector<Point>sol;int ans=get_qie(P,t,sol);
	if(!ans)return puts("[]"),void();
	if(ans==1)
	{
		double r=ang(sol[0]-P);if(dcmp(r)<0)r=PI+r;
		r*=180;r/=PI;
		printf("[%lf]
",r);
	}
	else
	{
		double r1=ang(sol[0]-P),r2=ang(sol[1]-P);
		if(dcmp(r1)<0)r1=PI+r1;if(dcmp(r2)<0)r2=PI+r2;
		r1*=180;r1/=PI;r2*=180;r2/=PI;if(dcmp(r1-r2)>0)std::swap(r1,r2);
		printf("[%lf,%lf]
",r1,r2);
	}
}
void work4()
{
	Point P=Point_in(),a=Point_in(),b=Point_in();
	double r;scanf("%lf",&r);
	Vec v=Normal(a-b)*r;
	if(dcmp(Cross(P-a,P-b))==0)
	{
		Point A=P+v,B=P-v;
		if(B<A)std::swap(A,B);
		printf("[(%lf,%lf),(%lf,%lf)]
",A.x,A.y,B.x,B.y);
		return;
	}
	int dd=dcmp(dis_Point_Line(P,a,b)-r*2);
	if(dd==0)
	{
		Point t=Line_jd(Line(a,b),Line(P,P+v));
		v=(t-P)/2;P=P+v;printf("[(%lf,%lf)]
",P.x,P.y);
		return;
	}
	if(dd>0){return printf("[]
"),void();}
	vector<Point>sol;sol.clear();
	Line a1(a+v,a+v+(a-b)),a2(a-v,a-v+(a-b));
	double t1,t2;Circle C(P,r);
	Line_Circle_jd(a1,C,t1,t2,sol);Line_Circle_jd(a2,C,t1,t2,sol);
	std::sort(sol.begin(),sol.end());
	printf("[(%lf,%lf),(%lf,%lf)]
",sol[0].x,sol[0].y,sol[1].x,sol[1].y);
}
void work5()
{
	vector<Point>sol;sol.clear();double r;
	Point A1=Point_in(),A2=Point_in(),B1=Point_in(),B2=Point_in();scanf("%lf",&r);
	Vec v1=A2-A1,v2=Normal(v1)*r,v3=B2-B1,v4=Normal(v3)*r;
	Line a1(A1+v2,A1+v2+v1),a2(A1-v2,A1-v2+v1),b1(B1+v4,B1+v4+v3),b2(B1-v4,B1-v4+v3);
	sol.push_back(Line_jd(a1,b1));
	sol.push_back(Line_jd(a1,b2));
	sol.push_back(Line_jd(a2,b1));
	sol.push_back(Line_jd(a2,b2));
	std::sort(sol.begin(),sol.end());
	printf("[(%lf,%lf),(%lf,%lf),(%lf,%lf),(%lf,%lf)]
",sol[0].x,sol[0].y,sol[1].x,sol[1].y,sol[2].x,sol[2].y,sol[3].x,sol[3].y);
}
void work6()
{
	Point a=Point_in();double r1;scanf("%lf",&r1);
	Point b=Point_in();double r2;scanf("%lf",&r2);
	double r;scanf("%lf",&r);
	Circle C1(a,r1+r),C2(b,r2+r);
	vector<Point>sol;sol.clear();
	int ans=Circle_Circle_jd(C1,C2,sol);
	if(!ans)return puts("[]"),void();
	if(ans==1)
		return printf("[(%lf,%lf)]
",sol[0].x,sol[0].y),void();
	if(sol[1]<sol[0])std::swap(sol[0],sol[1]);
	printf("[(%lf,%lf),(%lf,%lf)]
",sol[0].x,sol[0].y,sol[1].x,sol[1].y);
}
int main()
{
	while(~scanf("%s",str+1))
	{
		if(str[1]=='T')work3();
		else if(str[1]=='I')work2();
		else if(str[5]=='u')work1();
		else if(str[8]=='h')work4();
		else if(str[19]=='L')work5();
		else work6();
	}
	return 0;
}

LA2572

精度要掌握好,非常恶心。
保证数据在(5^{-13})范围内扰动不变,说明可见部分不会太小。
对于每个圆,找出所有的小圆弧,然后覆盖,具体实现看代码。
核心代码

Circle a[N];
int n,o[N];
inline int getans(Point P)
{
	for(int i=n;i;--i)
	{
		int d=dcmp(Len(P-a[i].c)-a[i].r);
		if(d==0)o[i]=1;
		if(d<0)return i;
	}
	return 0;
}
int main()
{
	vector<double>sol;
	while(n=in())
	{
		for(int i=1;i<=n;++i)a[i]=Circle_in();
		memset(o,0,sizeof o);
		for(int i=1;i<=n;++i)
		{
			sol.clear();sol.push_back(0),sol.push_back(PI*2);
			for(int j=1;j<=n;++j)
				if(i!=j)Circle_Circle_jd(a[i],a[j],sol);
			std::sort(sol.begin(),sol.end());
			for(int j=0,sz=sol.size()-1;j<sz;++j)
			{
				double r=(sol[j]+sol[j+1])/2.0;
				for(int t=-1;t<=1;t+=2)
				{
					int res=getans(a[i].get(r+eps*t));
					if(res>0)o[res]=1;
				}
			}
		}
		int ans=0;
		for(int i=1;i<=n;++i)ans+=o[i];
		printf("%d
",ans);
	}
	return 0;
}

UVA10652

凸包的模板题。

Point p[N*4],ch[N*4];int cnt;
inline void work()
{
	int n=in();cnt=0;double sum=0;
	for(int i=1;i<=n;++i)
	{
		Point o=Point_in();
		double w,h,j,Ang;scanf("%lf%lf%lf",&w,&h,&j);Ang=-torad(j);
		p[++cnt]=o+Rotate(Vec(-w/2,-h/2),Ang);
		p[++cnt]=o+Rotate(Vec(-w/2, h/2),Ang);
		p[++cnt]=o+Rotate(Vec( w/2,-h/2),Ang);
		p[++cnt]=o+Rotate(Vec( w/2, h/2),Ang);
		sum+=w*h;
	}
	int m=get_tb(p,cnt,ch);
	double res=abs_PolygonArea(ch,m);
	printf("%.1lf %%
",sum*100/res);
}
原文地址:https://www.cnblogs.com/cx233666/p/10090490.html