1. 程式人生 > >改進的point in polygon problem演算法介紹

改進的point in polygon problem演算法介紹

背景知識

點和多邊形的位置問題(point-in-polygon (PIP) problem), 一般指的是給定二維平面上的一個點Q以及一個多邊形P,怎樣判斷點Q是位於多邊形P內部還是外部。該演算法在計算機圖形學,地理空間資訊學等方面有廣泛的應用。目前有兩種通用的演算法實現: Ray casting algorithm(又稱even-odd演算法)和Winding number algorithm:

  • Ray casting algorithm
    從Q向無窮遠處引一條射線,計算該射線和多邊形的邊相交的次數,如果是奇數次,則點在多邊形內部,如果是偶數次,則點在多邊形外部。

    圖1

    圖1

    該演算法的理論依據是約當曲線定理 Jordan curve theorem: 一個簡單閉合多邊形一定會將平面分為內部和外部,連線內部任意一點與外部任意一點的線段一定經過閉合圖形的邊。雖然直觀上看上去理所當然,但是要給出嚴謹的數學證明卻很困難,數學家Jordan第一個給出了證明,後來雖然被找出該證明有不少瑕疵,這個定理仍然以他的名字來命名。
    根據約當曲線定理,可以對演算法有直觀的解釋:對於給定點Q和無窮遠點間S的連線QS,由於S必然位於多邊形P外部,而QS和多邊形的每個交點K都是一個外部和內部切換的臨界點:

    1. 如果是奇數個交點,從S出發,切換過程必然是:S(外部)–>
      k 1 (內部)–> k 2
      (外部)….–> k 2 n + 1 (內部)–>Q,這樣Q必然位於P內部。
    2. 如果是偶數個交點,從S出發,切換過程必然是:S(外部)–> k 1 (內部)–> k 2 (外部)….–> k 2 n (外部)–>Q,這樣Q必然位於P外部。
  • Winding number algorithm
    winding number algorithm 又稱為轉角累加法,知乎上有一個比較直觀的解釋如下:

    圖2

    圖2

    給多邊形的頂點依次編號(如上圖中為A,B,C,D,E),稱待判斷的點為O。依次考察角AOB、BOC、COD、DOE、EOA。這裡的角均取小於180度的那一側,並均為有向角,我們規定逆時針為正,順時針為負。上面左圖中,五個角都是正的;右圖中,角DOE是負的。把所有的角累加起來,我們發現左圖中角度之和為360度,右圖為零。
    形象一點兒理解,就是你站在O處,轉身用目光把整個多邊形掃描一遍。若你轉了奇數圈,則你在多邊形內;若你轉了偶數圈,則你在多邊形外。

對於簡單多邊形(邊不自交的多邊形),數學上說法是同構於單位圓的圖形,上面的兩種演算法通常能得到相同的結果。對於邊相交的非簡單多邊形,內部與外部通常有不同的定義(詳見page44):有的把重疊部分根據異或(XOR)原則,如下圖3中左側重疊部分定義為外部(Ray-casting演算法給出在此種定義下點與多邊形的關係),有的把重疊部分根據或(OR)原則定義為內部,如下圖3右側部分(Winding number演算法給出此種定義下點與多邊形的關係),

圖3

圖3 來源

由於一般的Ray-casting演算法不能夠處理如下圖4射線和頂點相交的情形(如點g,d)以及射線和多變形邊重合等特殊情形,本文的主要目的是介紹一種簡單有效的改進了的Ray-casting演算法。
圖4

圖4

演算法介紹

該演算法由Michael Galetzka和Patrick Glauner於2012年提出,2017年發表於Proceedings of the 12th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2017),在這裡可以下載到演算法全文。

由於點Q與多邊形P(頂點為 p 1 , p 2 , …, p n )位於二維平面上,可以考慮將Q和P通過座標平移,使Q平移到座標原點而其和P的相對位置不變,顯然平移不會改變點Q與多邊形P的位置關係(在多邊形內或外),下面演算法主要是針對平移後的情形進行說明,演算法步驟如下:

1. 判斷Q是否位於是P的頂點或者位於P的邊上,如果是,點Q在多邊形內部。
2.在頂點集合P中尋找不在x軸上的頂點 p s ,如果找不到,點Q在多邊形外部。
   (注:根據步驟1,Q不在P的頂點或者邊上,而找不到不在x軸上的點 p s
  明多邊形P的點都在X軸上,而Q位於原點,說明Q在P的外部。)
 3. 設 i = 1 ,從點 p s 開始通過重複下面步驟直到所有的頂點都被訪問到:
  a.判斷點 p s + i 是否位於x軸上,如果在,遞增i,如果 s + i > n ,那麼將i設定為 s ,從 p 0 開始繼續尋找,直到找到一個不位於x軸上的點 p s + i 為止。
  b. 根據步驟a中的查詢過程,採取下面的操作:
   i. 如果步驟a中找 p s + i 時沒有skip掉任何頂點,那麼判斷從 p s p s + i 的線段是否和x軸的正半軸相交,如相交,交點個數加1;
   ii. 如果步驟a中找 p s + i 時skip掉至少一個x軸座標為正的頂點,那麼判斷從 p s p s + i 的線段是否和整個x軸相交,如相交,交點個數加1;
   iii. 如果步驟a中找 p s + i 時skip掉至少一個x軸座標為負的頂點,不做任何操作;
   (注:不可能同時skip掉位於X軸正半軸和X軸負半軸的點,如果有,說明Q位於多邊形的邊上,步驟1就返回了)
  c. p s + i 為下一輪迭代的起點
4. 判斷交點個數是奇數還是偶數,如果是奇數,說明點Q位於多邊形P內部;如果是偶數,說明點Q位於多邊形P外部

這種方法實際是通過新增輔助線來判斷點是否位於多邊形的內部,這種新增輔助線的方法是否正確呢?我們來分析一下:

  • 情形i中無頂點skip掉則判斷線段是否和x正半軸相交的情形就是一般的Ray-casting演算法處理的情形,x軸的正半軸就是一般Ray-casting演算法中引出的射線;
  • 情形ii中如果skip掉的是x軸正向上的點,如下圖5所示,通過判斷線段 p 1 p 4 和X軸的交點,就可以替換線段 p 1 p 2 p 2 p 3 p 3 p 4 ,並且有效的繞開了射線和多邊形頂點相交以及與多邊形邊重合的情形。
    圖5

    圖5

    那麼為什麼是和整個X軸相交而不是像情形i中和x軸的正半軸相交那樣呢?這是為了處理下圖6中的情形:
    圖6

    圖6

    p 1 出發約過 p 2 ,我們實際上在 p 1 p 3 間添加了一條輔助線,如果這種情況下仍然只考慮線段 p 1 p 3 與x軸正半軸相交,將得不到任何交點,所以此種情況下,通過和x軸相交,得到一個交點。
  • 情形iii中為什麼越過的是x軸負半軸上的頂點時,不做任何處理呢?分析如下圖7可以得到答案:
    圖7

    圖7

    如圖7,從 p 1 開始,越過 p 2 ,我們找到 p 3 ,按照規則iii,計數為0,但是當從 p 3 開始到 p 1 ,我們沒有越過任何點,根據情形i,和x軸正半軸相交,計數為1,點在多邊形內部。當圖7三角向左移動變為如下圖8時,我們仍然可以結合情形iii和情形i,判斷點在多邊形外部。
    圖8

    圖8

那麼通過這種演算法能否覆蓋到所有情形呢?
根據Ray-casting演算法,我們不考慮越過的點在x負半軸的情形(因為Q位於原點,如果起點 p u 和終點