半平面交模板（O（n*n）&& O(n*log(n)）

摘自http://blog.csdn.net/accry/article/details/6070621

首先解决问题：什么是半平面？顾名思义，半平面就是指平面的一半，我们知道，一条直线可以将平面分为两个部分，那么这两个部分就叫做两个半平面。

然后，半平面怎么表示呢？二维坐标系下，直线可以表示为ax + by + c = 0，那么两个半平面则可以表示为ax + by + c >= 0 和ax + by + c < 0，这就是半平面的表示方法。

还有，半平面的交是神马玩意？其实就是一个方程组，让你画出满足若干个式子的坐标系上的区域（类似于线性规划的可行域），方程组就是由类似于上面的这些不等式组成的。

另外，半平面交可以干什么？半平面交虽然说是半平面的问题，但它其实就是关于直线的问题。一个一个的半平面其实就是一个一个有方向的直线而已。

半平面交的一个重要应用就是求多边形的核。多边形的核又是神马玩意？它是平面简单多边形的核是该多边形内部的一个点集，该点集中任意一点与多边形边界上一点的连线都处于这个多边形内部。就是一个在一个房子里面放一个摄像头，能将所有的地方监视到的放摄像头的地点的集合即为多边形的核。经常会遇到让你判定一个多边形是否有核的问题。

简易模板

int m;int cCnt,curCnt;

struct point

{

    double x,y;

    point(double x=,double y = ):x(x),y(y){}

};

point points[N],p[N],q[N];

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;

}

void initial()

{

    for(int i = ; i <= m; ++i)p[i] = points[i];

    p[m+] = p[]; p[] = p[m];

    cCnt = m;

}

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);

    point pt;

    pt.x=(x.x * v + y.x * u) / (u + v);

    pt.y=(x.y * v + y.y * u) / (u + v);

    return  pt;

}

void cut(double a,double b ,double c)

{

    curCnt = ;

    for(int i = ; i <= cCnt; ++i)

    {

        if(a*p[i].x + b*p[i].y + c >= )q[++curCnt] = p[i];

        else

        {

            if(a*p[i-].x + b*p[i-].y + c > )

            q[++curCnt] = intersect(p[i],p[i-],a,b,c);

            if(a*p[i+].x + b*p[i+].y + c > )

            q[++curCnt] = intersect(p[i],p[i+],a,b,c);

        }

    }

    for(int i = ; i <= curCnt; ++i)p[i] = q[i];

    p[curCnt+] = q[];

    p[] = p[curCnt];

    cCnt = curCnt;

}

void solve()

{

    initial();

    for(int i = ; i <= m; ++i)

    {

        double a,b,c;

        getline(points[i],points[i+],a,b,c);

        cut(a,b,c);

    }

    double area = ;

    for(int i = ; i <= cCnt; ++i)

        area += p[i].x * p[i + ].y - p[i + ].x * p[i].y;

    area = fabs(area / 2.0);

    printf("%.2f\n",area);

}

淘来一份模板

#include<iostream>

#include <stdio.h>

#include <math.h>

#define eps 1e-8

using namespace std;

const int MAXN=;

int m;int cCnt,curCnt;//此时cCnt为最终切割得到的多边形的顶点数、暂存顶点个数

struct point

{

    double x,y;

};

point points[MAXN],p[MAXN],q[MAXN];//读入的多边形的顶点（顺时针）、p为存放最终切割得到的多边形顶点的数组、暂存核的顶点

void getline(point x,point y,double &a,double &b,double   &c) //两点x、y确定一条直线a、b、c为其系数

{

    a = y.y - x.y;

    b = x.x - y.x;

    c = y.x * x.y - x.x * y.y;

}

void initial()

{

    for(int i = ; i <= m; ++i)p[i] = points[i];

    p[m+] = p[];

    p[] = p[m];

    cCnt = m;//cCnt为最终切割得到的多边形的顶点数，将其初始化为多边形的顶点的个数

}

point intersect(point x,point y,double a,double b,double c) //求x、y形成的直线与已知直线a、b、c、的交点

{

    double u = fabs(a * x.x + b * x.y + c);

    double v = fabs(a * y.x + b * y.y + c);

    point pt;

    pt.x=(x.x * v + y.x * u) / (u + v);

    pt.y=(x.y * v + y.y * u) / (u + v);

    return  pt;

}

void cut(double a,double b ,double c)

{

    curCnt = ;

    for(int i = ; i <= cCnt; ++i)

    {

        if(a*p[i].x + b*p[i].y + c >= )q[++curCnt] = p[i];// c由于精度问题，可能会偏小，所以有些点本应在右侧而没在，

        //故应该接着判断

        else

        {

            if(a*p[i-].x + b*p[i-].y + c > ) //如果p[i-1]在直线的右侧的话，

            {

                //则将p[i],p[i-1]形成的直线与已知直线的交点作为核的一个顶点(这样的话，由于精度的问题，核的面积可能会有所减少)

                q[++curCnt] = intersect(p[i],p[i-],a,b,c);

            }

            if(a*p[i+].x + b*p[i+].y + c > ) //原理同上

            {

                q[++curCnt] = intersect(p[i],p[i+],a,b,c);

            }

        }

    }

    for(int i = ; i <= curCnt; ++i)p[i] = q[i];//将q中暂存的核的顶点转移到p中

    p[curCnt+] = q[];

    p[] = p[curCnt];

    cCnt = curCnt;

}

void solve()

{

    //注意：默认点是顺时针，如果题目不是顺时针，规整化方向

    initial();

    for(int i = ; i <= m; ++i)

    {

        double a,b,c;

        getline(points[i],points[i+],a,b,c);

        cut(a,b,c);

    }

    /*

      如果要向内推进r，用该部分代替上个函数

      for(int i = 1; i <= m; ++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 = ;

    for(int i = ; i <= cCnt; ++i)

        area += p[i].x * p[i + ].y - p[i + ].x * p[i].y;

    area = fabs(area / 2.0);

    printf("%.2f\n",area);

}

/*void GuiZhengHua(){

     //规整化方向，逆时针变顺时针，顺时针变逆时针

    for(int i = 1; i < (m+1)/2; i ++)

      swap(points[i], points[m-i]);

}*/

int main()

{

    int t;

    cin>>t;

    while(t--)

    {

        cin>>m;

        int i;

        for(i=; i<=m; i++)

            cin>>points[i].x>>points[i].y;

        points[m+]=points[];

        solve();

    }

}

例题Poj1279

新增nlog(n)算法参考http://blog.csdn.net/acm_cxlove/article/details/7915167

#include<iostream>

#include<cstdio>

#include<algorithm>

#include<cmath>

#define eps 1e-10

#define N 50005

#define zero(a) (fabs(a)<eps)

using namespace std;

struct point

{

    double x,y;

    point() {}

    point(double tx,double ty)

    {

        x=tx;

        y=ty;

    }

} p[N],q[N];

int n,m;

struct Segment

{

    point s,e;

    double angle;

    void get_angle()

    {

        angle=atan2(e.y-s.y,e.x-s.x);

    }

} seg[N];

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);

}

double Get_Area(point pt[],int n)

{

    double area=;

    for(int i=; i<n-; i++)

        area+=xmul(pt[],pt[i],pt[i+]);

    return fabs(area)/;

}

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;

}

void HalfPlaneIntersect(Segment seg[],int n)

{

    int idx;

    for(int i=; i<n; i++)

        if(seg[i].angle+eps<seg[(i+)%n].angle&&seg[i].angle+eps<seg[(i-+n)%n].angle)

        {

            idx=i;

            break;

        }

    Segment deq[N];

    deq[]=seg[idx];

    deq[]=seg[(idx+)%n];

    int head=,tail=;

    idx=(idx+)%n;

    for(int i=; i<n; i++,idx=(idx+)%n)

    {

        while(head<tail&&xmul(seg[idx].s,seg[idx].e,Get_Intersect(deq[tail],deq[tail-]))<-eps) tail--;

        while(head<tail&&xmul(seg[idx].s,seg[idx].e,Get_Intersect(deq[head],deq[head+]))<-eps) head++;

        deq[++tail]=seg[idx];

    }

    while(head<tail&&xmul(deq[head].s,deq[head].e,Get_Intersect(deq[tail],deq[tail-]))<-eps) tail--;

    while(head<tail&&xmul(deq[tail].s,deq[tail].e,Get_Intersect(deq[head],deq[head+]))<-eps) head++;

    m=;

    if(tail==head) return;

    for(int i=head; i<tail; i++)

    {

        q[m++]=Get_Intersect(deq[i],deq[i+]);

    }

    if(tail>head+)

        q[m++]=Get_Intersect(deq[head],deq[tail]);

}

int slove(int mid)

{

    if(n-mid<=) return ;

    for(int i=; i<n; i++)

    {

        seg[i].s=p[i];

        seg[i].e=p[(i+mid+)%n];

        seg[i].get_angle();

    }

    HalfPlaneIntersect(seg,n);

    return zero(Get_Area(q,m));

}

int main()

{

    int t;

    scanf("%d",&t);

    while(t--)

    {

        scanf("%d",&n);

        for(int i=; i<n; i++)

            scanf("%lf%lf",&p[i].x,&p[i].y);

        for(int i=; i<=n/; i++) swap(p[i],p[n-i]);

        int ans,low=,high=n,mid;

        while(low<=high)

        {

            mid=(low+high)/;

            if(slove(mid))

            {

                ans=mid;

                high=mid-;

            }

            else low=mid+;

        }

        printf("%d\n",ans);

    }

    return ;

}

巴特西

半平面交模板（O（nn）&& O(nlog(n)）

最新文章

热门文章