1. 程式人生 > >初學計算幾何(五)——半平面交

初學計算幾何(五)——半平面交

前言

學了這麼久計算幾何,其實就是為了學半平面交

而學半平面交,其實就是為了做洛谷試煉場裡的這道題:【洛谷3222】[HNOI2012] 射箭


概念

  • 半平面:一條直線會把一個平面分成兩個半平面,通常我們把直線左邊的半平面視為我們所需要的半平面。
  • 半平面交:即若干半平面的交集,顯然它有以下幾種可能性:點/線段/直線/凸多邊形/無窮平面。為了避免產生無窮平面,我們一般會加入四個半平面作為限制。

那麼半平面交有什麼用?

比較常見的作用是求多邊形的核(如果多邊形中存在一個區域使得在該區域中可以看到多邊形中任意位置,則這個區域就是多邊形的核)。


大致思路

首先,我們將直線進行極角排序

,其中極角相同的直線我們讓較左邊的直線在前面,這一函式如下:

inline bool operator < (Line A,Line B) {return dcmp(A.Ang-B.Ang)?!~dcmp(A.Ang-B.Ang):!~dcmp(Cro(A.B-A.A,B.B-A.A));}

接下來,我們用一個雙向佇列\(ql\)來維護邊,同時開一個數組\(qp\),用\(qp_i\)儲存\(ql_i\)\(ql_{i+1}\)的交點。

則每當要加入一條新邊,我們就要維護一個凸多邊形,大致操作如下:

if(!dcmp(L[i].Ang-L[i-1].Ang)) continue;//如果出現極角相同的情況,將後出現的邊跳過
while(H<T&&IsOnRight(qp[T-1],L[i])) --T;while(H<T&&IsOnRight(qp[H],L[i])) ++H;//每次將首尾在當前邊右側的點彈出 
ql[++T]=L[i],qp[T-1]=LineIntersection(ql[T-1],ql[T]);//將邊加入佇列中,同時更新qp[T-1]

注意,其中的\(IsOnRight()\)函式利用了叉積性質,具體實現如下:

inline bool IsOnRight(Point P,Line L) {return !~dcmp(Cro(L.B-L.A,P-L.A));}//判斷一個點是否在一條直線右側

這樣操作一遍之後,前後其實還是有一些多餘的點,因此我們需要用首尾兩條邊再去判斷一次,將在這兩條邊右邊的點彈出:

while(H<T&&IsOnRight(qp[T-1],ql[H])) --T;while(H<T&&IsOnRight(qp[H],ql[T])) ++H;

最後判斷,如果佇列中只剩下小於等於\(2\)

個點,就說明面積為\(0\)

否則,可以把剩下點看成一個凸包,求凸包面積即可。


程式碼(板子題:【POJ2451】Uyuw's Concert

#include<cstdio>
#include<algorithm>
#include<cmath>
#define N 20000
using namespace std;
int n;
namespace ComputationGeometry//計算幾何
{
    #define eps 1e-10
    inline int dcmp(double x) {return fabs(x)<eps?0:(x>0?1:-1);}
    struct Point
    {
        double x,y;
        Point(double nx=0,double ny=0):x(nx),y(ny){}
    };
    typedef Point Vector;
    inline Vector operator + (Point A,Point B) {return Vector(A.x+B.x,A.y+B.y);}
    inline Vector operator - (Point A,Point B) {return Vector(A.x-B.x,A.y-B.y);}
    inline Vector operator * (Vector A,double x) {return Vector(A.x*x,A.y*x);}
    inline Vector operator / (Vector A,double x) {return Vector(A.x/x,A.y/x);} 
    inline bool operator < (Vector A,Vector B) {return fabs(A.x-B.x)>eps?A.x<B.x:A.y<B.y;}
    inline double Cro(Vector A,Vector B) {return A.x*B.y-A.y*B.x;}
    inline double Angle(Vector A) {return atan2(A.y,A.x);}
    struct Line
    {
        Point A,B;double Ang;
        Line(Point x=Point(0,0),Point y=Point(0,0)):A(x),B(y),Ang(Angle(y-x)){}
    };
    typedef Line Segment;
    inline bool operator < (Line A,Line B) {return dcmp(A.Ang-B.Ang)?!~dcmp(A.Ang-B.Ang):!~dcmp(Cro(A.B-A.A,B.B-A.A));}
    inline Point LineIntersection(Line L1,Line L2) {return L1.A+(L1.B-L1.A)*Cro(L2.B-L2.A,L1.A-L2.A)/Cro(L1.B-L1.A,L2.B-L2.A);}
    inline bool IsOnRight(Point P,Line L) {return !~dcmp(Cro(L.B-L.A,P-L.A));}
    struct Polygon
    {
        int n;Point p[N+5];
        Polygon() {n=0;}
    };
    typedef Polygon ConvexHull;
    inline double PolygonArea(Polygon P)
    {
        register int i;register double res=0;
        for(i=2;i<P.n;++i) res+=Cro(P.p[i]-P.p[1],P.p[i+1]-P.p[1])/2;
        return res;
    }
    inline double HalfPlaneIntersection(int tot,Line *L)//半平面交
    {
        register int i,H=1,T=1;static Line ql[N+5];static Point qp[N+5];sort(L+1,L+tot+1),ql[1]=L[1];
        for(i=2;i<=tot;++i)
        {
            if(!dcmp(L[i].Ang-L[i-1].Ang)) continue;
            while(H<T&&IsOnRight(qp[T-1],L[i])) --T;while(H<T&&IsOnRight(qp[H],L[i])) ++H; 
            ql[++T]=L[i],qp[T-1]=LineIntersection(ql[T-1],ql[T]);
        }
        while(H<T&&IsOnRight(qp[T-1],ql[H])) --T;while(H<T&&IsOnRight(qp[H],ql[T])) ++H;
        if(T-H<=1) return 0;
        register ConvexHull C;
        for(qp[T]=LineIntersection(ql[H],ql[T]),i=H;i<=T;++i) C.p[++C.n]=qp[i];
        return PolygonArea(C);
    }
}; 
using namespace ComputationGeometry;
Line L[N+5];
int main()
{
    register int i;register double X1,Y1,X2,Y2;
    for(scanf("%d",&n),i=1;i<=n;++i) scanf("%lf%lf%lf%lf",&X1,&Y1,&X2,&Y2),L[i]=Line(Point(X1,Y1),Point(X2,Y2));//讀入直線
    L[n+1]=Line(Point(0,0),Point(10000,0)),L[n+2]=Line(Point(10000,0),Point(10000,10000)),L[n+3]=Line(Point(10000,10000),Point(0,10000)),L[n+4]=Line(Point(0,10000),Point(0,0));//加入四個半平面作為限制
    return printf("%.1lf",HalfPlaneIntersection(n+4,L)),0;//求解答案
}