半平面交詳解

更好的閱讀體驗

定義：

半平面：

顧名思義，就是平面的一半。一條直線會把平面分成兩部分，就是兩個半平面。對於半平面，我們可以用直線方程式如：\(ax + by >= c\) 表示，更常用的是用直線表示。

半平面交：

顧名思義，就是多個半平面求交集。其結果可能是一個凸多邊形、無窮平面、直線、線段、點等。

多邊形的核：

如果多邊形中存在一個區域使得在區域中可以看到多邊形中任意位置（反之亦然），則這個區域就是多邊形的核。可以用半平面交來求解。

極點：

點 \((x,y)\) 與原點的連線與 \(x\) 軸的夾角，其范圍為 [0,360].

為了方便求解我們假設所有的直線的左側為我們所需要的半平面。

一般來說求解半平面交有兩種方法

① 分治法 \(O(n\log_2 n)\)

② 增量法 \(O(n\log_2 n)\)

但是在這里由於分治法常數較大，代碼實現較第二種復雜，所以我們着重介紹第二種方法。

算法流程

① 將所有直線極角排序,角度相同的保留下需要的一個

② 用一個雙端隊列存儲當前半平面交,每次通過判斷隊首與隊尾第一個交點是否滿足當前直線來更新

③ 先用隊尾判定隊首交點是否合法,再用隊首判斷隊尾交點是否合法

④ 最后求出來的半平面交是一個凸多邊形

極角排序

如下圖排序后

\((a,b)，(c,d)，(e,f)，(g,h)，(i,j)\)

所以我們可以更方便的逆時針依次構造，半平面交

由於我們規定了所有的直線的左側為我們所需要的半平面

所以極角相同的直線，我們保留最靠左的。

構造半平面交

我們可以依照以下流程來構造一個半平面交，並且構造完成的半平面交有多種情況

① 直線、線段、點不合法

② 凸多邊形，無窮平面（可以增加4個用於限制的半平面，使得平面變得有限）

我們維護兩個雙端隊列

一個儲存當前有用的直線（半平面），一個儲存半平面交的點。

我們依次加入每一條直線（半平面），在加入之前先將之前保存了的點，但不是最終半平面交中的點彈出隊列

如下圖我們首先讓 直線\(AB\) 進入隊列，再加入直線\(FG\)，並且求出 \(AB\) 與 \(FG\) 的交點 \(H\)，並且把它加入第二個隊列里

那么加入 直線\(CD\) 時我們會發現， \(H\) 這個點在半平面之外，那么就將它彈出隊列，同時將這條邊也彈出隊列。

所以只要一個點在加入的這條直線的右邊我們就將它彈出。

為什么我們要使用雙端隊列？

可以明顯的發現，雙端隊列的中的點是逆時針排列的

且滿足一定的單調性（這個需要自己畫圖思考，或者配合以下圖片思考）

當前加入一條直線，存放點的隊列中相鄰的兩個元素，所對應的直線的極角滿足 \(\angle a<=\angle b\)（\(a\) 對應點的在隊列中的位置位置小於 \(b\)）

那么也就是說，這條直線要么使隊首的點處於半平面外，要么使隊尾的點處於半平面外

但是因為構造半平面交是環形的構造，如下圖

我們順次處理 \(IJ\) ，\(AB\) ，\(CD\) ，\(EF\)，\(GH\) 所以當處理到 \(GH\) 的時候會發現，\(K\) 點是需要彈出的點但是，\(K\) 因為在處理 \(AB\) 時就已經加入了隊列而且處在隊列的首部，所以我們如果要彈出 \(K\) 就必須要維護一個雙端隊列來支持我們的彈出隊首的操作。

刪去不合法的點

最后隊首和隊尾都會剩下一些點，它們在半平面交之外

如下圖，由於半平面交構造時實際是一個上凸殼和下凸殼，然后是環形的，所以說，會發現某些情況下，有些點多余了

我們需要用隊首的直線，判斷一下多余的隊尾隊首的點。

POJ2451模板題，求半平面交面積

AC代碼，有注釋

#include<cmath>
#include<cstring>
#include<vector>
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
const double eps=1e-6;
const int maxn=2e5+10;
const double Pi=acos(-1.00);
inline int dcmp(double x)
{
	if(x>eps)return 1;
	return x<-eps?-1:0;
}
struct Vector
{
	double x,y;
	Vector(double X=0,double Y=0)
	{
		x=X,y=Y;
	}
	bool operator == (const Vector &b)const 
	{
		return dcmp(x-b.x)==0&&dcmp(y-b.y)==0;
	}
	double angle()
	{
		return atan2(y,x);//求出極角
	}
};
typedef Vector Point;
Vector operator + (Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator - (Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator * (Vector a,double b){return Vector(a.x*b,a.y*b);}
Vector operator / (Vector a,double b){return Vector(a.x/b,a.y/b);}
struct Line
{
	Point s,t;
	double ang;
	Line(Point X=Vector(),Point Y=Vector())
	{
		s=X,t=Y,ang=(Y-X).angle();
	}
};
typedef Line Segment;
double dot(Vector a,Vector b)
{
	return a.x*b.x+a.y*b.y;
}
double cross(Vector a,Vector b)
{
	return a.x*b.y-a.y*b.x;
}
bool is_parallel(Line a,Line b)//判斷a,b直線是否平行
{
	return dcmp(cross(a.t-a.s,b.t-b.s))==0;
}
Point intersection(Line a,Line b)//求出a,b的交點
{
	return a.s+(a.t-a.s)*(cross(b.t-b.s,a.s-b.s)/cross(a.t-a.s,b.t-b.s));
}
double area(Point *p,int n)//求出多邊形的面積
{
	double res=0;
	p[n+1]=p[1];
	for(int i=1;i<=n;i++)res+=cross(p[i],p[i+1]);
	return fabs(res/2);
}
bool operator < (const Line &a,const Line &b)//極角排序，如果極角相同則，選擇最靠左的直線
{
	double r=a.ang-b.ang;
	if(dcmp(r)!=0)return dcmp(r)==-1;
	return dcmp(cross(a.t-a.s,b.t-a.s))==-1;
}
bool OnRight(Line a,Point b)//檢查b是否在a直線的右邊
{
	return dcmp(cross(a.t-a.s,b-a.s))<0;
}
bool SI(Line *l,int n,Point *s,int &m)//增量法求半平面交
{
	static Line que[maxn];
	static Point que2[maxn];//兩個雙端隊列
	int head=0,tail=0;
	sort(l+1,l+1+n);
	que[0]=l[1];
	for(int i=2;i<=n;i++)
	if(dcmp(l[i].ang-l[i-1].ang)!=0)//極角相等的直線，取一個
	{
		if(head<tail&&(is_parallel(que[head],que[head+1])||is_parallel(que[tail],que[tail-1])))return false;//如果兩個直線共線，但是極角不同，則沒有半平面交
		while(head<tail&&OnRight(l[i],que2[tail-1]))tail--;//如果在直線右邊，刪除點
		while(head<tail&&OnRight(l[i],que2[head]))head++;
		que[++tail]=l[i];
		if(head<tail)que2[tail-1]=intersection(que[tail],que[tail-1]);//加入新點
	}
	while(head<tail&&OnRight(que[head],que2[tail-1]))tail--;//刪去多余點
	while(head<tail&&OnRight(que[tail],que2[head]))head++;
	if(tail-head<=1)return false;//只有一個點或零個點，沒有半平面交
	que2[tail]=intersection(que[head],que[tail]);//加入最后一條邊，和第一條邊的交點
	m=0;
	for(int i=head;i<=tail;i++)s[++m]=que2[i];
	return true;
}
const double lim=10000;
int n,m;
Point p[maxn];
Line l[maxn];
double solve()
{
	Point a=Point(0,0);//加入最大限制，防止半平面交無限大
	Point b=Point(lim,0);
	Point c=Point(lim,lim);
	Point d=Point(0,lim);
	l[++n]=Line(a,b);
	l[++n]=Line(b,c);
	l[++n]=Line(c,d);
	l[++n]=Line(d,a);
	if(!SI(l,n,p,m))return 0;
	return area(p,m);
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		Point a,b;
		scanf("%lf%lf%lf%lf",&a.x,&a.y,&b.x,&b.y);
		l[i]=Line(a,b);
	}
	printf("%.1f\n",solve());
}