半平面交 n^2和nlogn的模板

n^2普通模板:

 

/*半平面相交(直线切割多边形)(点标号从1开始)*/
Point points[MAXN],p[MAXN],q[MAXN];
int n;
double r;
int cCnt,curCnt;
inline void getline(Point x,Point y,double &a,double &b,double &c){
    a = y.y - x.y;
    b = x.x - y.x;
    c = y.x * x.y - x.x * y.y;
}
inline void initial(){
    for(int i = 1; i <= n; ++i)p[i] = points[i];
    p[n+1] = p[1];
    p[0] = p[n];
    cCnt = n;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
    double u = fabs(a * x.x + b * x.y + c);
    double v = fabs(a * y.x + b * y.y + c);
    return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
    curCnt = 0;
    for(int i = 1; i <= cCnt; ++i){
        if(a*p[i].x + b*p[i].y + c >= EPS)q[++curCnt] = p[i];
        else {
            if(a*p[i-1].x + b*p[i-1].y + c > EPS){
                q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
            }
            if(a*p[i+1].x + b*p[i+1].y + c > EPS){
                q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
            }
        }
    }
    for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
    p[curCnt+1] = q[1];p[0] = p[curCnt];
    cCnt = curCnt;
}
inline void solve(){
    //注意:默认点是顺时针,如果题目不是顺时针,规整化方向
    initial();
    for(int i = 1; i <= n; ++i){
        double a,b,c;
        getline(points[i],points[i+1],a,b,c);
        cut(a,b,c);
    }
    /*
    如果要向内推进r,用该部分代替上个函数
    for(int i = 1; i <= n; ++i){
        Point ta, tb, tt;
        tt.x = points[i+1].y - points[i].y;
        tt.y = points[i].x - points[i+1].x;
        double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
        tt.x = tt.x * k;
        tt.y = tt.y * k;
        ta.x = points[i].x + tt.x;
        ta.y = points[i].y + tt.y;
        tb.x = points[i+1].x + tt.x;
        tb.y = points[i+1].y + tt.y;
        double a,b,c;
        getline(ta,tb,a,b,c);
        cut(a,b,c);
    }
    */
    //多边形核的面积
    double area = 0;
    for(int i = 1; i <= curCnt; ++i)
        area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;
    area = fabs(area / 2.0);
    //此时cCnt为最终切割得到的多边形的顶点数,p为存放顶点的数组
}
inline void GuiZhengHua(){
     //规整化方向,逆时针变顺时针,顺时针变逆时针
    for(int i = 1; i < (n+1)/2; i ++)
      swap(points[i], points[n-i]);//头文件加iostream
}
inline void init(){
     for(int i = 1; i <= n; ++i)points[i].input();
        points[n+1] = points[1];
}

 

接下来是大牛的一个nlogn模板,大意就是维护一个真正有用的点集,然后尽量把没用的删除以便把n降到最低。对于能否真正降到log还心存疑惑,如果数据不巧,nlogn时间复杂度不会比n^2的低,或者更高,总之 慎用!

 

#include<iostream>
#include<fstream>
#include<iomanip>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<string>
#include<vector>
#include<sstream>
#include<cassert>
#define LL long long
#define eps 1e-10
#define inf 10000
#define zero(a) fabs(a)<eps
#define N 20005
using namespace std;
struct Point{
    double x,y;
}p[N*2];
struct Segment{
    Point s,e;
    double angle;
    void get_angle(){angle=atan2(e.y-s.y,e.x-s.x);}
}seg[N];
int m;
//叉积为正说明,p2在p0-p1的左侧
double xmul(Point p0,Point p1,Point p2){
    return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
Point Get_Intersect(Segment s1,Segment s2){
    double u=xmul(s1.s,s1.e,s2.s),v=xmul(s1.e,s1.s,s2.e);
    Point t;
    t.x=(s2.s.x*v+s2.e.x*u)/(u+v);t.y=(s2.s.y*v+s2.e.y*u)/(u+v);
    return t;
}
bool cmp(Segment s1,Segment s2){
    //先按极角排序
    if(s1.angle>s2.angle) return true;
    //极角相等,内侧的在前
    else if(zero(s1.angle-s2.angle)&&xmul(s2.s,s2.e,s1.e)>-eps) return true;
    return false;
}
void HalfPlaneIntersect(Segment seg[],int n){
    sort(seg,seg+n,cmp);
    int tmp=1;
    for(int i=1;i<n;i++)
        if(!zero(seg[i].angle-seg[tmp-1].angle))
            seg[tmp++]=seg[i];
    n=tmp;
    Segment deq[N];
    deq[0]=seg[0];deq[1]=seg[1];
    int head=0,tail=1;
    for(int i=2;i<n;i++){
        while(head<tail&&xmul(seg[i].s,seg[i].e,Get_Intersect(deq[tail],deq[tail-1]))<-eps) tail--;
        while(head<tail&&xmul(seg[i].s,seg[i].e,Get_Intersect(deq[head],deq[head+1]))<-eps) head++;
        deq[++tail]=seg[i];
    }
    while(head<tail&&xmul(deq[head].s,deq[head].e,Get_Intersect(deq[tail],deq[tail-1]))<-eps) tail--;
    while(head<tail&&xmul(deq[tail].s,deq[tail].e,Get_Intersect(deq[head],deq[head+1]))<-eps) head++;
    if(head==tail) return;
    m=0;
    for(int i=head;i<tail;i++)
        p[m++]=Get_Intersect(deq[i],deq[i+1]);
    if(tail>head+1)
        p[m++]=Get_Intersect(deq[head],deq[tail]);
}
double Get_area(Point p[],int &n){
    double area=0;
    for(int i=1;i<n-1;i++)
        area+=xmul(p[0],p[i],p[i+1]);
    return fabs(area)/2.0;
}
int main(){
    int n;
    while(scanf("%d",&n)!=EOF){
        seg[0].s.x=0;seg[0].s.y=0;seg[0].e.x=10000;seg[0].e.y=0;seg[0].get_angle();//设置边界
        seg[1].s.x=10000;seg[1].s.y=0;seg[1].e.x=10000;seg[1].e.y=10000;seg[1].get_angle();
        seg[2].s.x=10000;seg[2].s.y=10000;seg[2].e.x=0;seg[2].e.y=10000;seg[2].get_angle();
        seg[3].s.x=0;seg[3].s.y=10000;seg[3].e.x=0;seg[3].e.y=0;seg[3].get_angle();
        for(int i=0;i<n;i++){
            scanf("%lf%lf%lf%lf",&seg[i+4].s.x,&seg[i+4].s.y,&seg[i+4].e.x,&seg[i+4].e.y);//读点
            seg[i+4].get_angle();
        }
        HalfPlaneIntersect(seg,n+4);
        printf("%.1f\n",Get_area(p,m));
    }
    return 0;
}

posted on 2016-03-24 10:52  very_czy  阅读(346)  评论(0编辑  收藏  举报

导航