计算几何相关

orz hhw posted @ 2015年8月14日 09:25 in 算法学习 with tags 计算几何 模板 凸包 bzoj 半平面交 旋转卡壳 算法学习 随机增量 , 1334 阅读

差不多是最后一块算法了。。以后要多刷题了

一、凸包

凸包很久以前就学了,用极角排序扫描法就可以了,只有30行非常简单

#include<cstdio>
#include<algorithm>
#define N 222222
using namespace std;
int n,u,i, top;
struct P{double x,y;}p[N],q[N];
double dis(P a,P b){return (b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y);}
double multi(P a,P b,P c){return(b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);}
P operator-(P a,P b){P t;t.x=a.x-b.x;t.y=a.y-b.y;return t;}
double operator*(P a,P b){return a.x*b.y-a.y*b.x;}
bool operator<(P a,P b){
	double t=(a-p[1])*(b-p[1]);//t<0,a在b左边不换,否则a在b右边换
	return t<0||(t==0&&dis(a,p[1])<dis(b,p[1]));
}
int main(){
	scanf("%d",&n);u=1;
	for(i=1;i<=n;i++){
		scanf("%lf%lf",&p[i].x,&p[i].y);
		if(p[i].x<p[u].x||p[i].x==p[u].x&&p[i].y<p[u].y)u=i;
	}
	swap(p[1],p[u]);//将左下角顶点换到1
	sort(p+2,p+n+1);//按极角排序
	q[1]=p[1];q[2]=p[2];top=2;
	for(i=3;i<=n;i++){
		while(top>1&&(p[i]-q[top-1])*(q[top]-q[top-1])<=0)top--;//不符合条件退栈
		q[++top]=p[i];//符合条件入栈
	}
	for(i=1;i<=top;i++)printf("%.2lf %.2lf\n",q[i].x,q[i].y);
}

二、旋转卡壳

以前学的时候看到100+行的代码就望而却步,不过自己写写30多行足矣,就是写完凸包以后再把凸包绕一圈维护一些信息就可以了

BZOJ1069 做完凸包以后维护两个点为对角直径,剩下两个点可以单调找出,这样时间复杂度是n^2的可以达到题目要求

#include<cstdio>
#include<algorithm>
#include<cmath>
#define N 2222
using namespace std;
int n,u,i,j,x,y,top;double ans;
struct P{double x,y;}p[N],q[N];
double dis(P a){return a.x*a.x+a.y*a.y;}
P operator-(P a,P b){P t;t.x=a.x-b.x;t.y=a.y-b.y;return t;}
double operator*(P a,P b){return a.x*b.y-a.y*b.x;}
bool operator<(P a,P b){
	double t=(a-p[1])*(b-p[1]);
	return t<0||(t==0&&dis(a-p[1])<dis(b-p[1]));
}
int main(){
	scanf("%d",&n);u=1;
	for(i=1;i<=n;i++){
		scanf("%lf%lf",&p[i].x,&p[i].y);
		if(p[i].y<p[u].y||p[i].y==p[u].y&&p[i].x<p[u].x)u=i;
	}
	swap(p[1],p[u]);sort(p+2,p+n+1);
	q[1]=p[1];q[2]=p[2];top=2;
	for(i=3;i<=n;i++){
		while(top>1&&(p[i]-q[top-1])*(q[top]-q[top-1])<=0)top--;
		q[++top]=p[i];
	}
	q[top+1]=p[1];
	for(int x=1;x<=top;x++){
		int a=x%top+1,b=(x+2)%top+1;
		for(int y=x+2;y<=top;y++){
	for(;a%top+1!=y&&(q[y]-q[x])*(q[a+1]-q[x])>(q[y]-q[x])*(q[a]-q[x]);a=a%top+1);
	for(;b%top+1!=x&&(q[b+1]-q[x])*(q[y]-q[x])>(q[b]-q[x])*(q[y]-q[x]);b=b%top+1);
			ans=max((q[y]-q[x])*(q[a]-q[x])+(q[b]-q[x])*(q[y]-q[x]),ans);
		}
	}
	printf("%.3lf\n",ans/2);
}

BZOJ1185 做完凸包以后对每条边维护左边、右边和上边最多能到的地方,然后统计答案就可以了

#include<cstdio>
#include<algorithm>
#include<cmath>
#define N 55555
#define eps 1e-8
using namespace std;
int n,u,i,l,r,w,top,fir;double ans,D,L,R,H;
struct P{double x,y;}p[N],q[N],t[5];
double dis(P a){return sqrt(a.x*a.x+a.y*a.y);}
double operator*(P a,P b){return a.x*b.y-a.y*b.x;}
P operator*(P a,double b){P t;t.x=a.x*b;t.y=a.y*b;return t;}
P operator-(P a,P b){P t;t.x=a.x-b.x;t.y=a.y-b.y;return t;}
P operator+(P a,P b){P t;t.x=a.x+b.x;t.y=a.y+b.y;return t;}
double operator/(P a,P b){return a.x*b.x+a.y*b.y;}
bool operator<(P a,P b){
	double t=(a-p[1])*(b-p[1]);
	return t>0||(t==0&&dis(a-p[1])<dis(b-p[1]));
}
int main(){
	scanf("%d",&n);u=1;ans=1e60;
	for(i=1;i<=n;i++){
		scanf("%lf%lf",&p[i].x,&p[i].y);
		if(p[i].y<p[u].y||p[i].y==p[u].y&&p[i].x<p[u].x)u=i;
	}
	swap(p[1],p[u]);sort(p+2,p+n+1);
	q[1]=p[1];q[2]=p[2];top=2;
	for(i=3;i<=n;i++){
		for(;top>1&&(p[i]-q[top-1])*(q[top]-q[top-1])>=0;top--);
		q[++top]=p[i];
	}
	q[top+1]=q[1];
	l=1;r=1;w=1;
	for(i=1;i<=top;i++){
		D=dis(q[i]-q[i+1]);
		for(;(q[i+1]-q[i])*(q[w+1]-q[i])-(q[i+1]-q[i])*(q[w]-q[i])>-eps;w=w%top+1);
		for(;(q[i+1]-q[i])/(q[r+1]-q[i])-(q[i+1]-q[i])/(q[r]-q[i])>-eps;r=r%top+1);
		if(i==1)l=r;
		for(;(q[i+1]-q[i])/(q[l+1]-q[i])-(q[i+1]-q[i])/(q[l]-q[i])<eps;l=l%top+1);
		L=(q[i+1]-q[i])/(q[l]-q[i])/D,R=(q[i+1]-q[i])/(q[r]-q[i])/D;
		H=(q[i+1]-q[i])*(q[w]-q[i])/D;if(H<0)H=-H;
        double tmp=(R-L)*H;  
        if(tmp<ans){
            ans=tmp;
            t[0]=q[i]+(q[i+1]-q[i])*(R/D);
			t[1]=t[0]+(q[r]-t[0])*(H/dis(t[0]-q[r]));
			t[2]=t[1]-(t[0]-q[i])*((R-L)/dis(q[i]-t[0]));
			t[3]=t[2]-(t[1]-t[0]);
        }
	}
	printf("%.5lf\n",ans);
	for(int i=1;i<=3;i++)if(t[i].y<t[fir].y||t[i].y==t[fir].y&&t[i].x<t[fir].x)fir=i;
	for(int i=0;i<=3;i++)printf("%.5lf %.5lf\n",t[(i+fir)%4].x,t[(i+fir)%4].y);
}

三、半平面交

差不多就是把很多个直线交一下求出一个核,不过我只是把模板打了一遍

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define N 2222
using namespace std;
typedef double DB;int T,n,top,pot,i,j,pn,q[N];
struct P{DB x,y;}p[N];struct L{P a,b;DB ang;}l[N];
DB mul(P a,P b,P c){return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);}
bool cmp(L a,L b){DB d=a.ang-b.ang;return d<0||d==0&&mul(a.a,b.a,b.b)<0;}//取更靠近法向量的 
void addline(L &l,DB x1,DB y1,DB x2,DB y2){
	l.a.x=x1;l.a.y=y1;l.b.x=x2;l.b.y=y2;
	l.ang=atan2(y2-y1,x2-x1);
}
void insert(L l1,L l2,P &p){
	DB a1=l1.b.y-l1.a.y,b1=l1.a.x-l1.b.x,c1=(l1.b.x-l1.a.x)*l1.a.y-(l1.b.y-l1.a.y)*l1.a.x;
	DB a2=l2.b.y-l2.a.y,b2=l2.a.x-l2.b.x,c2=(l2.b.x-l2.a.x)*l2.a.y-(l2.b.y-l2.a.y)*l2.a.x;
	p.x=(c2*b1-c1*b2)/(a1*b2-a2*b1);p.y=(c1*a2-c2*a1)/(a1*b2-a2*b1);
}
bool judge(L a,L b,L c){P p;insert(b,c,p);return mul(p,a.a,a.b)>0;}
int main(){
	scanf("%d",&T);
	while(T--){
		scanf("%d",&n);for(i=0;i<n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
		for(i=0;i<n-1;i++)addline(l[i],p[i].x,p[i].y,p[i+1].x,p[i+1].y);//加直线 
		addline(l[i],p[i].x,p[i].y,p[0].x,p[0].y);
		sort(l,l+n,cmp);DB d=0;//将半平面按极角排序
		for(i=0,j=0;i<n;i++)if(l[i].ang-l[j].ang>0)l[++j]=l[i];
		n=j+1;q[0]=0;q[1]=1;top=1;pot=0;
		for(i=2;i<n;q[++top]=i++){
			for(;top>pot&&judge(l[i],l[q[top]],l[q[top-1]]);top--);//while最顶端的2个半平面的交点在当前半平面外,删除最顶端的半平面 
			for(;top>pot&&judge(l[i],l[q[pot]],l[q[pot+1]]);pot++);//while最底端的2个半平面的交点在当前半平面外,删除最底端的半平面 
		}
		for(;top>pot&&judge(l[q[pot]],l[q[top]],l[q[top-1]]);top--);//while最顶端的2个半平面的交点在最底端半平面外,删除最顶端的半平面
		for(;top>pot&&judge(l[q[top]],l[q[pot]],l[q[pot+1]]);pot++);//while最底端的2个半平面的交点在最顶端半平面外,删除最底端的半平面 
		for(pn=0,q[++top]=q[i=pot];i<top;i++,pn++)insert(l[q[i+1]],l[q[i]],p[pn]);//用来求交集的半平面 
		if(pn>=3)for(i=1;i<pn-1;i++)d+=mul(p[0],p[i],p[i+1]);d/=2;//面积 
		printf("%.2lf\n",d<0?-d:d);
	}
}

四、随机增量法

感觉就是随机一下求解顺序然后暴力搞就可以了,而且感觉只有一个问题需要这么解决,就是最小圆覆盖

BZOJ2823 裸的最小圆覆盖

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
struct P{double x,y;}p[500005],c,ret;
double r;int n,i,j,k;//c为圆心,r为半径 
double dis(P a, P b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
P get(P a,P b,P c){//返回三角形的外心   
	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,d=a1*b2-a2*b1;  
    ret.x=a.x+(c1*b2-c2*b1)/d;ret.y=a.y+(a1*c2-a2*c1)/d;  
    return ret;   
}
int main(){
	scanf("%d",&n);for(i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
	random_shuffle(p+1,p+n+1);c=p[1];r=0;
	for(i=2;i<=n;i++)if(dis(p[i],c)>r){//第一个点
		c=p[i];r=0;
		for(j=1;j<i;j++)if(dis(p[j],c)>r){//第二个点
			c.x=(p[i].x+p[j].x)/2;c.y=(p[i].y+p[j].y)/2;
			r=dis(c,p[j]);
			for(k=1;k<j;k++)if(dis(p[k],c)>r){//第三个点 
				c=get(p[i],p[j],p[k]);
				r=dis(c,p[k]);
			}
		}
	}
	printf("%.2lf %.2lf %.2lf",c.x,c.y,r);
}

BZOJ3564 坐标变换后的最小圆覆盖

#include<cstdio>
#include<cmath>
#include<algorithm>
#define eps 1e-10
using namespace std;
struct P{double x,y;}p[500005],c,ret;
double r,ang,bs;int n,i,j,k;
double dis(P a, P b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
P get(P a,P b,P c){
    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,d=a1*b2-a2*b1;  
    ret.x=a.x+(c1*b2-c2*b1)/d;ret.y=a.y+(a1*c2-a2*c1)/d;  
    return ret;   
}
void change(int i,double k){
	double x=p[i].x*cos(k)+p[i].y*sin(k);  
    double y=-1*p[i].x*sin(k)+p[i].y*cos(k);  
    p[i].x=x;p[i].y=y;  
}
int main(){
	scanf("%d",&n);for(i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
	scanf("%lf%lf",&ang,&bs);
	ang=ang/180*acos(-1.0);
	for(i=1;i<=n;i++){
		change(i,ang);
		p[i].x/=bs;
	}
	random_shuffle(p+1,p+n+1);
	c=p[1];r=0;
	for(i=2;i<=n;i++)if(dis(p[i],c)>r+eps){
		c=p[i];r=0;
		for(j=1;j<i;j++)if(dis(p[j],c)>r+eps){
			c.x=(p[i].x+p[j].x)/2;
			c.y=(p[i].y+p[j].y)/2;
			r=dis(c,p[j]);
			for(k=1;k<j;k++)if(dis(p[k],c)>r+eps){
				c=get(p[i],p[j],p[k]);
				r=dis(c,p[k]);
			}
		}
	}
	printf("%.3lf",r);
}

BZOJ2280 二分+随机增量+倍增+二分

#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#define eps 1e-9
#define N 100010
using namespace std;
int n,m,i,cnt,ans[N][2];
double l,r,mid,R;
struct P{double x,y;}a[N],b[N],O;
double dis(P x,P y){return sqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y));}
P get(P a,P b,P c){
	double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2,a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2,d=a1*b2-a2*b1;
	return (P){a.x+(c1*b2-c2*b1)/d,a.y+(a1*c2-a2*c1)/d};
}
void cal(int l,int r){
	int i,j,k,n=0;
	for(i=l;i<=r;i++)a[++n]=b[i];
	for(i=2;i<=n;i++)swap(a[1+(rand()%n)],a[i]);
	for(O=a[1],R=0,i=2;i<=n;i++)if(R+eps<dis(O,a[i]))
		for(O=a[i],R=0,j=1;j<i;j++)if(R+eps<dis(O,a[j]))
			for(O=(P){(a[i].x+a[j].x)/2,(a[i].y+a[j].y)/2},R=dis(O,a[i]),k=1;k<j;k++)
				if(R+eps<dis(O,a[k]))O=get(a[i],a[j],a[k]),R=dis(O,a[k]);
}
bool ok(double x){
	int i,j,l,r,mid,t,sum=0;
	for(i=1;i<=n;i=t+1){
		for(j=1;i+(1<<j)-1<=n;j++){
			cal(i,i+(1<<j)-1);
			if(R>x+eps)break;
		}
		t=i,l=i+(1<<(j-1))-1;r=i+(1<<j)-1;if(r>n)r=n;
		while(l<=r){
			cal(i,mid=(l+r)>>1);
			if(R<x+eps)l=(t=mid)+1;else r=mid-1;
		}
		if(++sum>m)return 0;
		ans[++cnt][0]=i;ans[cnt][1]=t;
	}
	return 1;
}
int main(){
	scanf("%d%d",&n,&m);for(i=1;i<=n;i++)scanf("%lf%lf",&b[i].x,&b[i].y);
	cal(1,n);r=R;
	if(m>1)for(;r-l>eps;cnt=0)if(ok(mid=(l+r)/2))r=mid;else l=mid;
	ok(r);printf("%.8lf\n%d\n",r,cnt);
	for(i=1;i<=cnt;i++)cal(ans[i][0],ans[i][1]),printf("%.8lf %.8lf\n",O.x,O.y);
}

五、其它

Pick定理:顶点在格点上多边形面积=多边形内部点数+边上点数/2 -1,感觉做数学题挺有用的

三角剖分:不会


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter