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 （半平面交）
时间： 2024-11-06 21:09:13