首页 > 代码库 > POJ2451 Uyuw's Concert (半平面交)
POJ2451 Uyuw's Concert (半平面交)
POJ2451 给定N个半平面 求他们的交的面积。 N<=20000
首先参考 POJ1279
多边形的核 其实就是这里要求的半平面交 但是POJ1279数据较小 O(n^2)的算法 看起来是要TLE的
但是试着提交了一下 一遍就A了。。。 看来暴力的半平面切割法实际表现远远好于O(n^2)
如果数据再大一点呢? N=100000?
有一个办法可以优化到O(nlogn)
step1. 将所有半平面按极角排序,对于极角相同的,选择性的保留一个。 O(nlogn)
step2. 使用一个双端队列(deque),加入最开始2个半平面。
step3. 每次考虑一个新的半平面:
a.while deque顶端的两个半平面的交点在当前半平面外:删除deque顶端的半平面
b.while deque底部的两个半平面的交点在当前半平面外:删除deque底部的半平面
c.将新半平面加入deque顶端
step4.删除两端多余的半平面。
具体方法是:
a.while deque顶端的两个半平面的交点在底部半平面外:删除deque顶端的半平面
b.while deque底部的两个半平面的交点在顶端半平面外:删除deque底部的半平面
重复a,b直到不能删除为止。
step5:计算出deque顶端和底部的交点即可。
先给出O(n^2)的代码 以后来补 nlogn的
#include<iostream> #include<stdio.h> #include<stdlib.h> #include<string.h> #include<math.h> #include<algorithm> #include<queue> #include<vector> using namespace std; const double eps=1e-9; int cmp(double x) { if(fabs(x)<eps)return 0; if(x>0)return 1; else return -1; } const double pi=acos(-1.0); inline double sqr(double x) { return x*x; } struct point { double x,y; point (){} point (double a,double b):x(a),y(b){} void input() { scanf("%lf%lf",&x,&y); } friend point operator +(const point &a,const point &b) { return point(a.x+b.x,a.y+b.y); } friend point operator -(const point &a,const point &b) { return point(a.x-b.x,a.y-b.y); } friend bool operator ==(const point &a,const point &b) { return cmp(a.x-b.x)==0&&cmp(a.y-b.y)==0; } friend point operator *(const point &a,const double &b) { return point(a.x*b,a.y*b); } friend point operator*(const double &a,const point &b) { return point(a*b.x,a*b.y); } friend point operator /(const point &a,const double &b) { return point(a.x/b,a.y/b); } double norm() { return sqr(x)+sqr(y); } }; struct line { point a,b; line(){}; line(point x,point y):a(x),b(y) { } }; double det(const point &a,const point &b) { return a.x*b.y-a.y*b.x; } double dot(const point &a,const point &b) { return a.x*b.x+a.y*b.y; } double dist(const point &a,const point &b) { return (a-b).norm(); } point rotate_point(const point &p,double A) { double tx=p.x,ty=p.y; return point(tx*cos(A)-ty*sin(A),tx*sin(A)+ty*cos(A)); } bool parallel(line a,line b) { return !cmp(det(a.a-a.b,b.a-b.b)); } bool line_joined(line a,line b,point &res) { if(parallel(a,b))return false; double s1=det(a.a-b.a,b.b-b.a); double s2=det(a.b-b.a,b.b-b.a); res=(s1*a.b-s2*a.a)/(s1-s2); return true; } bool pointonSegment(point p,point s,point t) { return cmp(det(p-s,t-s))==0&&cmp(dot(p-s,p-t))<=0; } void PointProjLine(const point p,const point s,const point t,point &cp) { double r=dot((t-s),(p-s))/dot(t-s,t-s); cp=s+r*(t-s); } struct polygon_convex { vector<point>P; polygon_convex(int Size=0) { P.resize(Size); } }; struct halfPlane { double a,b,c; halfPlane(point p,point q) { a=p.y-q.y; b=q.x-p.x; c=det(p,q); } halfPlane(double aa,double bb,double cc) { a=aa;b=bb;c=cc; } }; double calc(halfPlane &L,point &a) { return a.x*L.a +a.y*L.b+L.c; } point Intersect(point &a,point &b,halfPlane &L) { point res; double t1=calc(L,a),t2=calc(L,b); res.x=(t2*a.x-t1*b.x)/(t2-t1); res.y=(t2*a.y-t1*b.y)/(t2-t1); //cout<<res.x<<" "<<res.y<<endl; return res; } polygon_convex cut(polygon_convex &a,halfPlane &L) { int n=a.P.size(); polygon_convex res; for(int i=0;i<n;i++) { if(calc(L,a.P[i])>-eps)res.P.push_back(a.P[i]); else { int j; j=i-1; if(j<0)j=n-1; if(calc(L,a.P[j])>-eps) res.P.push_back(Intersect(a.P[j],a.P[i],L)); j=i+1; if(j==n)j=0; if(calc(L,a.P[j])>-eps) res.P.push_back(Intersect(a.P[i],a.P[j],L)); } } return res; } double INF=10000; polygon_convex core(vector<point> &a) { polygon_convex res; res.P.push_back(point(0,0)); res.P.push_back(point(INF,0)); res.P.push_back(point(INF,INF)); res.P.push_back(point(0,INF)); int n=a.size(); for(int i=0;i<n-1;i+=2) { halfPlane L(a[i],a[(i+1)]); res=cut(res,L); } return res; } bool comp_less(const point &a,const point &b) { return cmp(a.x-b.x)<0||cmp(a.x-b.x)==0&&cmp(a.y-b.y)<0; } polygon_convex convex_hull(vector<point> a) { polygon_convex res(2*a.size()+5); sort(a.begin(),a.end(),comp_less); a.erase(unique(a.begin(),a.end()),a.end());//删去重复点 int m=0; for(int i=0;i<a.size();i++) { while(m>1&&cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)--m; res.P[m++]=a[i]; } int k=m; for(int i=int(a.size())-2;i>=0;--i) { while(m>k&&cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)--m; res.P[m++]=a[i]; } res.P.resize(m); if(a.size()>1)res.P.resize(m-1); return res; } bool is_convex(vector<point> &a) { for(int i=0;i<a.size();i++) { int i1=(i+1)%int(a.size()); int i2=(i+2)%int(a.size()); int i3=(i+3)%int(a.size()); if((cmp(det(a[i1]-a[i],a[i2]-a[i1]))*cmp(det(a[i2]-a[i1],a[i3]-a[i2])))<0) return false; } return true; } int containO(const polygon_convex &a,const point &b) { int n=a.P.size(); point g=(a.P[0]+a.P[n/3]+a.P[2*n/3])/3.0; int l=0,r=n; while(l+1<r) { int mid=(l+r)/2; if(cmp(det(a.P[l]-g,a.P[mid]-g))>0) { if(cmp(det(a.P[l]-g,b-g))>=0&&cmp(det(a.P[mid]-g,b-g))<0)r=mid; else l=mid; }else { if(cmp(det(a.P[l]-g,b-g))<0&&cmp(det(a.P[mid]-g,b-g))>=0)l=mid; else r=mid; } } r%=n; int z=cmp(det(a.P[r]-b,a.P[l]-b))-1; if(z==-2)return 1; return z; } long long int distant(point a,point b) { return (int(b.x)-int(a.x))*(int(b.x)-int(a.x))+(int(b.y)-int(a.y))*(int(b.y)-int(a.y)); } double convex_diameter(polygon_convex &a,int &First,int &Second) { vector<point> &p=a.P; int n=p.size(); double maxd=0; if(n==1) { First=Second=0; return maxd; } #define next(i)((i+1)%n) for(int i=0,j=1;i<n;++i) { while(cmp(det(p[next(i)]-p[i],p[j]-p[i])-det(p[next(i)]-p[i],p[next(j)]-p[i]))<0) j=next(j); double d=dist(p[i],p[j]); if(d>maxd) { maxd=d; First=i,Second=j; } d=dist(p[next(i)],p[next(j)]); if(d>maxd) { maxd=d; First=next(i),Second=next(j); } } return maxd; } double area(vector<point>a) { double sum=0; for(int i=0;i<a.size();i++) sum+=det(a[(i+1)%(a.size())],a[i]); return sum/2.; } vector<point> pp; int main() {freopen("t.txt","r",stdin); int T=1; while(T--) { int n; scanf("%d",&n); n*=2; pp.resize(n); for(int i=0;i<n-1;i+=2) pp[i].input(),pp[i+1].input(); polygon_convex pc=core(pp); printf("%.1lf\n",-area(pc.P)+0.005); } return 0; }
POJ2451 Uyuw's Concert (半平面交)
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。