二維幾何-半平面交
簡單的說,半平面交問題就是給出若干個半平面,求它們的公共部分。其中每個半平面用一條有向直線表示。它的左側就是它所代表的半平面。
在很多情況下,這個半平面交都是一個凸多邊形,但也有時候會得到一個無界多邊形,甚至是一條直線、線段或者是點。不管怎麼樣,結果移動是凸的(因為凸集的交是凸的)。當然,半平面交也可以為空。
計算半平面交的一個方法是增量法,即初始答案為整個平面,然後逐一加入各個半平面,維護當前的半平面交,為了程式設計方便,我們一般用一個很大的矩形(4個半平面的交)代替整個平面,計算出結果以後刪除這4個人工半平面,這樣,每加入一個平面都相當於用一條有向直線去切割多邊形。
切割的方法很簡單:按照逆時針順序考慮多邊形的所有頂點,保留在直線左側和直線上的點,而刪除直線右邊的點。如果有向直線和多邊形相交產生了新的點,這些點應加在新多邊形中。具體來說,沒考慮完一個頂點Pi,在考慮Pi+1之前先判斷PiPi+1是否和有向直線在PiPpi+1的內部(端點不算)相交。如果是,則還要把交點加入新多邊形中。
用有向直線A->B切割多邊形poly,返回左側。如果退化,可能會返回單點或者線段。
Polygon CutPolygon(const Polygon &poly, Point A, Point B) { int i, n; Polygon newpoly; Point C, D, ip; n = poly.size(); for(i = 0; i < n; i++) { C = poly[i]; D = poly[(i+1)%n]; if(dcmp(Cross(B-A, C-A)) >= 0) newpoly.push_back(C); if(dcmp(Cross(B-A, C-D)) != 0) { ip = GetLineIntersection(A, B-A, C, D-C); if(OnSegment(ip, C, D)) newpoly.push_back(ip); } } return newpoly; }
可惜每次切割需要O(n)時間,因此上述演算法的時間複雜度為O(n^2),不是很優秀。
和凸包類似,半平面交也可以通過排序、掃描的方法在O(nlogn)時間內解決,不同的是凸包用的棧,而半平面用的是雙端佇列。注意,按照極角排序後,每次新加入的半平面可能會讓隊尾的半平面變得無用,從而需要刪除。
注意,由於極角序是環形的,新加的半平面也可能饒了一圈以後讓隊首的半平面變得無用。
點p在直線L的左邊,線上不算。
bool OnLeft(Line L, Point p)
{
return Cross(L.v, p-L.p) > 0;
}
二直線交點,假定交點唯一存在
Point GetIntersection(Line a, Line b) { Vector u; double t; u = a.p - b.p; t = Cross(b.v, u) / Cross(a.v, b.v); return a.p + a.v*t; }
半平面交的主過程
int HalfplaneIntersection(Line *L, int n, Point *poly)
{
int i, m, first, last;
Point p[maxn]; //p[i]為q[i]和q[i+1]的交點
Line q[maxn]; //雙端佇列
sort(L, L+n); //按極角排序
first = 0; //雙端佇列的第一個元素的下標
last = 0; //雙端佇列的最後一個元素的下標
q[first] = L[0]; //雙端佇列初始化為只有一個半平面L[0]
for(i = 1; i < n; i++)
{
while(first < last && !OnLeft(L[i], p[last-1])) //新的半平面讓隊尾變得無用
last--;
while(first < last && !OnLeft(L[i], p[first])) //新的半平面讓隊首變得無用
first++;
q[++last] = L[i]; //將新的半平面加入佇列中
if(dcmp(Cross(q[last].v, q[last-1].v)) == 0)
{
//兩向量平行且同向,取內側的一個
last--;
if(OnLeft(q[last], L[i].p))
q[last] = L[i];
}
if(first < last) //獲得交點
p[last-1] = GetIntersection(q[last], q[last-1]);
}
//刪除無用平面 首部半平面可能讓隊尾元素無用
while(first < last && !OnLeft(q[first], p[last-1]))
last--;
if(last - first <= 1) //空集
return 0;
p[last] = GetIntersection(q[first], q[last]); // 計算首尾兩個半平面的交點
//從deque複製到輸出中
m = 0;
for(i = first; i <= last; i++)
poly[m++] = p[i];
return m;
}
在大多數情況下,若干半平面的交是一個凸多邊形,但也有例外,比如可能得到的是一個無界的區域,解決方法和前面一樣,在外面加一個座標很大(注意不要讓運算溢位)的包圍框(由4個半平面組成),最後再把框去掉。但對於保證不會產生無界區域的情況下,就不需要加這4個特殊的半平面了。