1. 程式人生 > >二維幾何-半平面交

二維幾何-半平面交

簡單的說,半平面交問題就是給出若干個半平面,求它們的公共部分。其中每個半平面用一條有向直線表示。它的左側就是它所代表的半平面。

在很多情況下,這個半平面交都是一個凸多邊形,但也有時候會得到一個無界多邊形,甚至是一條直線、線段或者是點。不管怎麼樣,結果移動是凸的(因為凸集的交是凸的)。當然,半平面交也可以為空。

計算半平面交的一個方法是增量法,即初始答案為整個平面,然後逐一加入各個半平面,維護當前的半平面交,為了程式設計方便,我們一般用一個很大的矩形(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個特殊的半平面了。