凸多边形面积并 计算几何+扫描线

思路:

(本题现在可以在JDFZOJ1005找到)

我们首先来围观一篇神题解.

http://www.cnblogs.com/ch3656468/archive/2011/10/17/2215551.html

(果然\(O(n^3)\)不是最优的么...)

首先,将所有的凸多边形上的点按照逆时针排列,随后将所有的向量分为向左和向右两种,竖直方向的向量直接无视.

对于某些向量,它们可能在一条直线上,因此我们就将这些向量一起处理.

对于所有向左的向量(他们在上边界上),我们处理出所有在凸多边形边界上,但却不在凸多边形内部的部分与\(x\)轴形成的梯形面积的和.

然后对于向右的向量们也用同样的方法计算.

(然后看图不知道怎么样就发现)答案就是两者相减.时间复杂度仅为\(O(n^2logn)\).

我的代码:

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<climits>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef double f2;
 
static const f2 eps=1e-8;
static const f2 INF=1e50;
inline int dcmp(const f2&x){return fabs(x)<eps?0:(x<0?-1:1);}
 
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
    bool operator<(const Point&B)const{return x<B.x||(dcmp(x-B.x)==0&&y<B.y);}
    bool operator!=(const Point&B)const{return dcmp(x-B.x)!=0||dcmp(y-B.y)!=0;}
    bool operator==(const Point&B)const{return dcmp(x-B.x)==0&&dcmp(y-B.y)==0;}
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p)const{return Point(p*x,p*y);}
    Point operator/(const f2&p)const{return Point(x/p,y/p);}
};
f2 Cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
f2 Dot(const Point&a,const Point&b){return a.x*b.x+a.y*b.y;}
 
struct Line{
    Point p,v;
    Line(){}
    Line(const Point&p1,const Point&p2):p(p1),v(p1-p2){}
    f2 asky(f2 x){
        return dcmp(x-p.x)==0?p.y:p.y+(x-p.x)/v.x*v.y;
    }
    bool operator==(const Line&B)const{
        f2 b=dcmp(0-p.x)==0?p.y:p.y+(0-p.x)/v.x*v.y;
        f2 _b=dcmp(0-B.p.x)==0?B.p.y:B.p.y+(0-B.p.x)/B.v.x*B.v.y;
        return dcmp(Cross(v,B.v))==0&&dcmp(b-_b)==0;
    }
    bool operator<(const Line&B)const{
        f2 k1=dcmp(v.x)==0?INF:v.y/v.x,k2=dcmp(B.v.x)==0?INF:B.v.y/B.v.x;
        return k1<k2;
    }
 
};
Point GetLineIntersection(const Line&l1,const Line&l2){
    if(dcmp(Cross(l1.v,l2.v))==0)return Point(-INF,-INF);
    return l1.p+l1.v*(Cross(l1.p-l2.p,l2.v)/Cross(l1.v,l2.v));
}
 
#define N 110
int n;
struct Polygen{
    Point P[N];int siz;
    void Init(){scanf("%d",&siz);for(int i=1;i<=siz;++i)scanf("%lf%lf",&P[i].x,&P[i].y);P[siz+1]=P[1];}
}Poly[11];
 
struct Quest{
    f2 x;int type;
    Quest(){}
    Quest(f2 _x,int _type):x(_x),type(_type){}
    bool operator<(const Quest&B)const{return x<B.x;}
};
Line Lines[N*N],reallines[N*N],q1[N*N],q2[N*N],seq[N*N];int cnt,id,n1,n2;
 
inline f2 Solve(Line*L,int m){
    sort(L+1,L+m+1);
    register int i,j,k;
    for(id=0,i=1;i<=m;i=j+1){for(j=i;j!=m&&L[j]==L[j+1];++j);seq[++id]=L[i];}
     
    f2 res=0;
     
    vector<Quest>q;
    for(i=1;i<=id;++i){
        q.clear();
        for(j=1;j<=n;++j){
            bool touch=0;
            for(k=1;k<=Poly[j].siz&&!touch;++k)if(seq[i]==Line(Poly[j].P[k],Poly[j].P[k+1]))touch=1;
            if(touch)continue;
             
            vector<Point>get;
            for(k=1;k<=Poly[j].siz;++k){
                Point x=GetLineIntersection(seq[i],Line(Poly[j].P[k],Poly[j].P[k+1]));
                if(x!=Point(-INF,-INF)&&Dot(Poly[j].P[k]-x,Poly[j].P[k+1]-x)<=0)get.push_back(x);
            }
            sort(get.begin(),get.end());
            vector<Point>realget;
            for(k=0;k<get.size();++k)if(!k||get[k]!=get[k-1])realget.push_back(get[k]);
             
            if(realget.size()==2){
                q.push_back(Quest(realget[0].x,0)),q.push_back(Quest(realget[1].x,1));
            }
        }
         
        f2 Mx,Mn;
        for(j=1;j<=m;++j)if(dcmp(Cross(seq[i].v,L[j].v))==0&&dcmp(seq[i].asky(0)-L[j].asky(0))==0){
            Mx=max(L[j].p.x,(L[j].p+L[j].v).x),Mn=min(L[j].p.x,(L[j].p+L[j].v).x);
            q.push_back(Quest(Mn,2)),q.push_back(Quest(Mx,3));
        }
         
        sort(q.begin(),q.end());
         
        int in=0,cnt=0;
        for(j=0;j<q.size();++j){
            f2 y0=seq[i].asky(q[j==0?0:j-1].x),y1=seq[i].asky(q[j].x);
            if(!in&&cnt)res+=(y0+y1)*(j==0?0:q[j].x-q[j-1].x)*.5;
            if(q[j].type==0)++in;if(q[j].type==1)--in;
            if(q[j].type==2)++cnt;if(q[j].type==3)--cnt;
        }
    }
     
    return res;
}
int main(){
    //freopen("tt.in","r",stdin);
    scanf("%d",&n);register int i,j;
    for(i=1;i<=n;++i){
        Poly[i].Init();
        for(j=1;j<=Poly[i].siz;++j){
            if(dcmp(Poly[i].P[j].x-Poly[i].P[j+1].x)==0)continue;
            if(dcmp(Poly[i].P[j].x-Poly[i].P[j+1].x)<0)q2[++n2]=Line(Poly[i].P[j],Poly[i].P[j+1]);
            else q1[++n1]=Line(Poly[i].P[j],Poly[i].P[j+1]);
        }
    }
     
    printf("%.3lf",Solve(q1,n1)-Solve(q2,n2));
     
    return 0;
}

JDFZOJ2204:打地鼠 期望+凸包

思路:
首先\(O(n^3)\)的思路非常朴素,只需要考虑每条向量的贡献即可.向量的出现概率即为这两个点出现且向量右侧的点均不出现.于是我们需要暴力\(O(n)\)check右侧的点是否均不出现.

我们考虑以一个点为起点的所有向量的答案.

那么对于每一个向量右侧的点必定都是一段区间,且满足单调性.

 

我们枚举起点,然后把剩下的点按照关于起点的极角序排序,然后依次扫描即可.

当然有几个坑点:

概率有一些为0的不能直接乘除,而需要记录零的个数.

还有就是坑爹的三点共线问题.

 

我们考虑当角度相同时,按照距离从大到小排序.

然后扫描的时候进行一些特判.这里没太弄明白就先挖坑了.

(好像一开始先把点随机扰动一下应该也可以吧)

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
 
typedef double f2;
 
#define N 1010
struct Point{
    f2 x,y,p;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p)const{return Point(p*x,p*y);}
}P[N];
inline f2 sqr(const f2&x){return x*x;}
inline f2 cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
inline f2 dot(const Point&a,const Point&b){return a.x*b.x+a.y*b.y;}
inline f2 length(const Point&a){return sqrt(sqr(a.x)+sqr(a.y));}
 
static const f2 eps=1e-8;
inline int dcmp(const f2&x){return fabs(x)<eps?0:(x<0?-1:1);}
 
struct Seg{
    Point v;
    f2 Ang,area,p,len;
    inline bool operator<(const Seg&B)const{
        return dcmp(Ang-B.Ang)==0?len>B.len:Ang<B.Ang;
    }
}S[N];
 
int main(){
    int n;scanf("%d",&n);register int i,j,k;
    for(i=0;i<n;++i)scanf("%lf%lf%lf",&P[i].x,&P[i].y,&P[i].p),P[i].p/=1000.0;
    f2 ans=0;
    for(i=0;i<n;++i){
        int id=0;
        for(j=0;j<n;++j)if(i^j){
            S[id].v=P[i]-P[j];
            S[id].Ang=atan2(S[id].v.y,S[id].v.x);
            S[id].area=cross(P[i],P[j]);
            S[id].p=P[j].p;
            S[id].len=length(S[id].v);
            id++;
        }
        sort(S,S+id);
        long double nowp=P[i].p;
        int zero=0;
        for(j=0;j<id;++j)if(dcmp(S[j].p-1)==0)++zero;else nowp*=(1-S[j].p);
         
        k=0;
        for(j=0;j<id;++j){
            while(dcmp(cross(S[j].v,S[k].v))>0||(dcmp(cross(S[j].v,S[k].v))==0&&dcmp(dot(S[k].v,S[j].v-S[k].v))<=0)){
                if(dcmp(1-S[k].p))nowp/=1-S[k].p;else --zero;
                k=(k+1)%id;
                if(k==j)break;
            }
            if(!zero)ans+=S[j].area*nowp*S[j].p;
            if(dcmp(1-S[j].p)!=0)nowp*=(1-S[j].p);else ++zero;
        }
    }
     
    printf("%.6lf",ans*.5);
    return 0;
}