加強局部簡便計算的點在多邊形內的高效判定

王盛春，王文成，譚雪晗，李 靜

?

加強局部簡便計算的點在多邊形內的高效判定

王盛春1,2，王文成1,2，譚雪晗1,2，李 靜3

(1. 中國科學院軟件研究所計算機科學國家重點實驗室，北京 100190； 2. 中國科學院大學計算機科學與技術學院，北京 100190； 3. 中國科學院動物研究所動物進化與系統學院重點實驗室，北京 100101)

多邊形；網格；簡便計算；點在多邊形內的檢測

點在多邊形內的判定檢測(point-in-polygon test)是計算幾何中一個基本問題，在計算機圖形學、模式識別、衛星定位系統及地理信息系統等方面有著廣泛的應用[1]。該問題的基本描述為：給定一個多邊形與一個任意的點，判斷點是否位于多邊形內。

點在多邊形內的判定檢測算法可以分為2類：①算法不需要預處理，直接對多邊形進行某些參數的計算得到判定結果，如常見的射線法[2-3]、轉角法[4]、面積和法[5]等。其中，經典的射線法(ray-crossing)通過從測試點向某方向作一條射線，根據該射線與多邊形邊的交點數目的奇偶性進行判定。若交點數目為奇數，則該點位于多邊形內；否則位于多邊形外。這類算法簡單直觀，可處理任意多邊形，通常被用作與其他算法對比的基準算法。但這類算法需要逐邊操作，計算時間復雜度為()。②算法通過預處理，將多邊形分解為簡單幾何形狀并進行有序管理，如梯形法[6-7]、凸剖分法[8]、層次三角形法[9]等，或建立高效的空間劃分結構，如均勻網格法[1,10-11]、分層法[12]、BSP樹法[8,13]等。通過預處理的組織管理，這類算法能夠降低判定計算的復雜度，如凸剖分法[8]復雜度為(log2)，其中為多邊形的邊數。文獻[1]對這些算法進行了對比分析，并提出一種基于均勻網格的簡便快速算法。近年來，學界雖然提出了一些新的相關算法，如文獻[14]提出的一種基于sign()函數的點在多邊形內外判別算法，將二維平面內的點轉化為三維空間在平面上的點，借用符號函數判斷待測點與多邊形的位置關系，但在計算效率上，仍與文獻[1]有所差距。

本文在文獻[1]算法的基礎上進行如下改進，以進一步提高檢測速度：

(1) 本文將網格中心點屬性的計算，轉換為網格線交點屬性的計算，并記錄各網格交點一側(本文為其右側)的網格邊與多邊形邊的相交情況。這樣，在計算多邊形的邊與網格單元相交的情況時，即可同時進行這些計算，避免了另行計算網格單元中心點屬性的計算。

(2) 不同于文獻[1]僅記錄與其相交的多邊形的邊，本文方法記錄位于該網格中的多邊形的邊片段。由于邊片段的縮小，在測試線與多邊形的邊進行求交計算時，可基于簡單的坐標對比大幅剔除不相交的邊，提高求交測試效率。

(3) 將測試點與附近已知屬性的網格交點的連線，轉換為與坐標軸平行的2條相連的直線段。這依然遵循了射線法的基本原則，但與坐標軸平行的直線段與多邊形的邊的求交計算，可簡化以加速。

1 改進算法

本文提出的新方法仍基于均勻網格方法，包括2個階段：預處理和判定計算。

(1) 預處理。首先建立均勻網格；然后記錄網格線與多邊形邊的交點坐標以及各個網格單元中所包含的邊的片段；最后逐個計算網格交點位于多邊形內/外的位置屬性。在該屬性計算時，位于包圍盒最外圍的網格交點一定是位于多邊形外的(所取的包圍盒比多邊形的緊致包圍盒稍大一點)，而其他網格交點的位置屬性可依據其鄰近已知屬性的網格點、其之間連線與多邊形的邊的相交次數，即可獲得。

(2) 判定計算。首先確定待測點所在的網格單元；然后沿軸方向向鄰近網格邊做垂線，基于該垂線在網格邊上的交點得到最近網格交點，形成一條平行軸的連線；最后根據以上這兩條平行于坐標軸的線段與多邊形邊的交點個數的奇偶性與網格交點的位置屬性，即可判斷待測點的位置屬性。

1.1 算法概述

1.1.1 預處理

預處理過程包括：創建均勻網格、記錄多邊形與網格邊的交點坐標、記錄各個網格單元包含的多邊形邊的片段和計算網格頂點的位置屬性4部分。

與其他均勻網格法類似，由于網格分辨率的大小決定網格創建的效率和后期測試點屬性判斷的速度，因此，創建網格時需確定網格的分辨率。通常，網格的分辨率越高，則判定的時間越短。然而，當網格的分辨率越高，預處理時間與存儲空間也將增加。因此，在實際應用中，需要同時考慮到預處理時間、判定時間和存儲空間來決定優化的網格分辨率。本文采用與文獻[1]相同的方法確定網格的分辨率，即設網格的寬長比為，多邊形的邊數為，則和方向上的單元數分別為

文獻[1]在確定網格分辨率之后，采用數字微分分析器(digital differential analyzer，DDA)[15]算法將通過每個網格單元的整條邊記錄在相關單元中。由于整條邊與多個網格單元相關，求交測試時難以基于簡單的坐標對比剔除不相交的邊。對此，本文通過計算多邊形邊與網格邊的交點，得到各網格單元中多邊形邊的片段。這樣，基于簡單的坐標對比，可提高剔除不相交邊的概率，有利于加速。

各網格線與多邊形邊的交點，可形成順序排列的數組。由于各條網格線的某一坐標值是相同的，因此本文方法僅存儲網格線上變化的坐標值，從而降低存儲空間的消耗量。由此，各網格交點，可基于其之間連線上的交點個數，得到其位置屬性。

圖1為計算網格頂點位置屬性的一個實例。最外層網格邊上的網格點均位于多邊形外，如藍色點所示。與文獻[1]類似，根據已知位置屬性的網格交點，考察相鄰網格點間連線與多邊形邊的相交次數，可推知相鄰網格點的位置屬性。以未知鄰近網格點(3,1)為例，已知(3,0)點位于多邊形外，網格邊(3,0)(3,1)與1條多邊形的邊相交，因此(3,1)點位于多邊形內。同理，由于網格邊(3,1)(3,2)與多邊形的邊沒有相交，可推知(3,2)點位于多邊形內。

1.1.2 判定計算

在確定一個待測點的位置屬性時，本文方法首先確定待測點所在的網格單元；然后沿軸方向向鄰近的網格邊做垂線，基于該垂線與網格邊的交點找到最近網格頂點，形成一條平行軸的連線，最后依據兩條平行于坐標軸的線段與多邊形邊交點個數的奇偶性與網格交點的位置屬性，判斷待測點的位置屬性。

圖1 計算網格點位置屬性的實例圖

以圖2中的待測點1為例。首先，定位1點位于網格單元[1,1]中，并從待測點1出發向鄰近的網格邊做垂線交于點1。連接點1與其左側的網格點(1,1)。最后，遍歷網格單元[1,1]中的所有邊，計數與線段11有交點的邊。同時依據網格邊1對應的交點存儲數組，統計線段1(1,1)間的交點個數。將2條測試邊與多邊形邊相交點的個數相加，交點總個數為1。因此，待測點1與網格點(1,1)的位置屬性相反，位于多邊形外。

圖2 點包容性判定實例圖

在本文算法的判定計算中，最主要的工作是檢測測試邊的垂線是否與網格單元內的邊相交。對此，可在該垂線向上與向下兩個方向的線段中，選取較短者進行處理，以利于減少相交邊的計算。如圖3所示，對于垂線段1：1，考察其與一條多邊形邊的片段1：(設的坐標為(1,1)，的坐標為(2,2))的相交情況，可按下式進行

當進一步處理上述結果時，根據邊片段1的方程計算垂線1與邊片段1交點的軸坐標1；然后，與垂線段1的軸坐標范圍(設1的軸坐標范圍為(,))比較，判定垂線段1是否與邊片段1相交。具體比較方法如下

由此，對于不相交的邊，本文能夠基于簡單的坐標比較實現快速剔除；對于相交的邊，與傳統斜邊求交中需要4次叉積運算和4次比較不同，本文方法只需1次乘法、1次加法和2次比較即可完成判定。

圖3 判斷垂線與多邊形邊是否相交實例圖

1.2 奇異情況的處理

本文方法的本質是射線法的局部化實現，當多邊形的邊或頂點位于射線上時(奇異情況)，將提高相交次數統計的復雜性。

奇異情況主要分為2大類：測試線與多邊形頂點相交的奇異點情況和測試線與多邊形的邊共線的奇異邊情況。經分析，該奇異情況可分為以下4種情況處理，圖4和圖5分別展示了奇異點與奇異邊的4種奇異情況。

(1) 若多邊形頂點位于網格邊上(或與網格交點重合)且連接該頂點的2條邊位于該網格邊異側，如圖4A區和圖4C區，則該網格邊(如1)關于這個多邊形頂點的交點個數記為1；

(2) 若多邊形頂點位于網格邊上(或與網格頂點重合)且連接該頂點的2條邊位于網格邊的同側，如圖4B區和圖4D區，則交點個數記為2；

(3) 若多邊形的邊與網格邊重合且連接該邊的2條邊位于網格邊的異側，如圖5B區和圖5D區關于網格線1和3的情況，則交點個數記為1；

(4) 若多邊形的邊與網格邊重合且連接該邊的2條邊位于網格邊的同側，如圖5A區和圖5C區關于網格線0和2的情況，則交點個數記為2。

圖4 奇異點的相關情況

圖5 奇異邊的相關情況

在預處理過程中，只需對奇異點/邊進行標注即可。在判定計算中，測試線遇到奇異點/邊時，則沿著其方向繼續前行，直至遇到非奇異情況，然后進行相關處理即可。文獻[16]將文獻[1]中的方法擴展到三維空間中，并對奇異情況進行了詳細的分析。一般情況下，在奇異邊繼續向前延伸中，最多延伸兩次即可得到非奇異情況。

1.2.1 預處理

預處理階段的基本操作就是統計2個網格點之間的交點個數，然后根據交點個數的奇偶性判定網格點的位置屬性。就奇異點而言，可以圖4中網格點(2,2)和(3,2)進行分析，假設網格點(2,0)和(3,0)的位置屬性為多邊形外。在確定網格點(2,2)位置屬性時，首先需要向上搜索鄰近的網格點，由于網格點(2,1)為奇異點，因此向上搜索至網格點(2,0)。統計網格點(2,0)和(2,2)之間的交點個數，并依據上述的方法，連接頂點3的2條邊位于網格邊的異側，交點個數為1，因此網格點(2,2)和(2,0)的位置屬性相反，位于多邊形內，結論正確。同理，確定網格點(3,2)的位置屬性，依據上述的方法，網格點(3,2)和(3,0)之間的交點個數為2個，因此網格點(3,2)與網格點(3,0)位置屬性相同，位于多邊形外，結論正確。就奇異邊而言，只需將奇異邊看作成一個奇異點，采用與奇異點相同的方法進行計算。

在奇異情況的處理過程中，最主要的操作是判斷奇異點/邊兩端的邊是否處于奇異點/邊的同側或者異側。對于該問題，可通過比較2條邊的頂點坐標解決。以圖5A區為例，該情況的奇異邊垂直軸，可設方程為，因此，通過對比奇異邊兩端邊的坐標值就能夠判斷出相應邊的分布情況。若1邊2個端點坐標的值均不小于(或均不大于)值且1邊兩個端點坐標的值均不小于(或均不大于)值，則邊1和1處于奇異邊1的同側；反之，邊1和1處于奇異邊1的異側。同理，對于平行軸的奇異邊，只需通過對應的軸坐標值的比較即可得出結果。

1.2.2 判定計算

判定階段的奇異情況主要包括：測試邊與網格邊重合、測試邊與多邊形的邊重合和測試邊經過多邊形的頂點。對于該類奇異情況，可以采用類似奇異邊的方式進行計算。

(1) 對于測試邊與網格邊重合的奇異情況。其處理方法為：直接根據待測點與鄰近網格點間的交點個數的奇偶性確定待測點的位置屬性。如圖2中的待測點4，其測試邊與網格線重合。因此可以通過連接待測點4與網格點(5,3)，然后根據該線段之間的交點個數確定待測點4的位置屬性。由于該測試邊通過多邊形的頂點6，且與該頂點相連的兩條邊位于測試邊的異側，交點個數為1。因此待測點4的位置屬性與網格點(5,3)位置屬性相反，位于多邊形外。

(2) 對于測試邊與多邊形的邊重合的奇異情況。其處理方法為：延伸測試邊，直至穿出多邊形的邊，并交于鄰近的網格邊，然后連接到與其交點最近的網格點，最后統計測試邊與多邊形的交點個數，以此確定待測點的位置屬性。如圖2中的待測點2，其測試邊與多邊形邊1211重合。依據上述規定，延長測試邊使其與網格線0交于點2，同時連接點2與網格點(0,0)。統計測試邊與多邊形邊的交點個數，由于邊1211兩端連接的邊處于測試邊的同側，其交點個數記為2。另外，測試邊與多邊形邊10相交，交點總數為3，因此，待測點2的位置屬性與網格點(0,0)位置屬性相反，位于多邊形內。

(3) 對于測試邊經過多邊形的頂點的奇異情況。其處理方法與預處理階段的奇異點情況相似，可利用相同的處理方法進行計算。如圖2中的待測點3和5，其測試邊分別經過頂點8和10。由于連接頂點8的兩條邊處于3測試邊的兩側，且連接頂點10的兩條邊處于5測試邊的同側，測試邊33，55與多邊形的交點個數分別為1和4。因此3的位置屬性與網格點(3,1)的位置屬性相反，位于多邊形外；5的位置屬性與網格點(1,2)的位置屬性相同，位于多邊形內。

2 復雜度分析

根據文獻[1]對網格中心點算法的分析，網格單元數的復雜度為()，空間復雜度為()。由于本文方法采用與其相同的公式確定網格分辨率，而在預處理過程中僅為每個網格交點增加了一個bit位記錄其位置屬性，且增加了()個字節存儲多邊形邊與網格邊的交點，因此本文方法的空間復雜度仍為()。

本文方法的預處理過程如下：

(1) 創建均勻網格。本文方法與文獻[1]采用相同的公式確定網格分辨率，時間復雜度為()。

(2) 計算網格線與多邊形邊的交點坐標，并將其存儲到相應的數組中。由于一條邊一般只與數個網格單元相交，檢測條邊，其時間復雜度為()。

(3) 將所有的多邊形邊片段加入到對應的網格單元中，類似文獻[10]，時間復雜度為()。

(4) 確定網格交點的位置屬性。對于每個網格交點，只需檢測其臨近網格點的位置屬性和查詢其之間網格邊對應的交點數目即可，時間復雜度為()。

本文方法的判定計算過程如下：

(1) 被測點位于無邊的網格內。被測點與臨近網格頂點具有相同的屬性，時間復雜度為(1)。

雖然本文方法與文獻[1]具有相同的時間復雜度，但本文方法通過形成平行坐標軸的測試線、記錄多邊形位于各個網格單元中的邊片段，由此可基于一些簡單計算，大幅提高測試線與多邊形的邊的求交計算效率。

3 實驗與討論

本文在一臺微機上進行了算法實現，該微機的硬件環境為Intel(R) Core(TM) i7-8700 CPU @ 3.20 GHz，16 GB ARM，軟件開發環境為MS Windows 10操作系統，visual studio 2017。近年來，點在多邊形內判定算法的研究不多且進展緩慢，其中較好的算法有凸剖分法[8]、CBCA[10]和QCPM[11]等，這些方法在文獻[1]中均有對比分析。由于本文方法是在文獻[1]的基礎上進行改進，且上述方法在性能方面均次于文獻[1]方法，因此，本文方法僅與文獻[1]方法進行對比實驗。測試用例來自文獻[10]，如圖6所示。測試用點1000×1000個，其均勻分布在比多邊形包圍盒稍大的范圍內。

Pol10Pol100 Pol1249Pol28000

實驗統計數據見表1，其中，判定時間為檢測所有測試點的總時間。由表1可知，由于本文方法需要計算多邊形邊與網格邊的交點，并對這些交點進行排序和存儲，因此在預處理階段本文方法耗費了更多時間。但在判定計算過程中，相較于網格中心點算法[1]，本文方法的效率提升了2倍多，且隨著多邊形的邊數增多、復雜性的增大，本文方法的加速更多。

表1 新算法與網格中心點算法的對比實驗數據

多邊形邊數網格分辨率算法預處理時間(加速率rp)# (ms)判定時間(加速率rd)^ (ms) Pol10108×6網格中心點算法0.33331.965 本文新算法0.339 (–0.017 70)15.199 (1.103 10) Pol10010038×18網格中心點算法0.40929.036 本文新算法1.422 (–0.712 38)12.747 (1.277 87) Pol12491 249100×46網格中心點算法0.90026.064 本文新算法9.889 (–0.908 99)12.918 (1.017 65) Pol2800028 000100×100網格中心點算法4.73484.590 本文新算法195.272 (–0.975 76)24.913 (2.395 42)

(注：# r=(T–T)/T，T為網格中心點算法的預處理時間；T為新算法的預處理時間；^ r=(T?T)/T，T為網格中心點算法的判定計算時間；T為新算法的判定計算時間)

表2 新算法與網格中心點算法的數值計算次數的對比實驗數據

多邊形邊數網格分辨率網格中心點算法本文新算法 乘法次數加法次數比較次數乘法次數加法次數比較次數 Pol10108×67 296 05614 592 1123 648 028164 778164 778329 556 Pol10010038×185 309 34410 618 6882 654 672968 8496 884193 768 Pol12491 249100×465 251 41610 502 8322 625 70888 93788 937177 874 Pol2800028 000100×10027 896 33655 792 67213 948 16897 41397 413194 826

為了進一步分析本文方法加速的原因，本文針對網格中心點算法和本文方法的判定計算過程的數值計算次數進行了統計分析，即測試邊與多邊形邊的求交過程的計算量，實驗結果見表2，其中，乘法次數和加法次數分別代表測試邊與多邊形邊在求交過程中需要進行乘法運算和加法運算的次數。對于比較次數，在網格中心點算法中，主要是比較叉積值與0值的次數，而在本文算法中，主要是比較測試邊和多邊形邊交點坐標值與測試邊坐標值的次數。如第1.1.2節中所述，網格中心點算法中計算兩條斜邊是否相交需要4次叉積運算和4次比較運算，即需要進行16次加法和8次乘法，而本文方法在一次求交過程中只需要1次乘法、1次加法和2次比較。從表2可看出，在確定網格分辨率和測試點的情況下，本文方法所需的乘法計算和加法計算次數明顯低于網格中心點算法所需的計算次數，同時，本文方法所需的比較次數也明顯低于網格中心點算法所需的比較次數。

4 結束語

本文提出一種新的基于均勻網格的點在多邊形內的判定算法。通過加強簡便計算的處理，很好地減少測試線與多邊形邊的求交計算開銷，提高了射線法局部化實現的效率。實驗表明，相比目前最快的類似算法，可將測試計算的效率提高2倍多，且多邊形的邊越多越復雜，加速性能越高。

[1] 李靜, 王文成. 基于網格中心點的點在多邊形內的高效判定[J]. 軟件學報, 2012, 23(9): 2481-2488.

[2] HUGHES J F, VAN DAM A, FOLEY J D, et al. Computer graphics: Principles and practice [M]. New Jersey: Pearson Education, 2014: 81-98.

[3] FEITO F, TORRES J C, URE?A A. Orientation, simplicity, and inclusion test for planar polygons [J]. Mathematical Gazette, 1995, 19(4): 595-600.

[4] FEITO F R, TORRES J C. Inclusion test for general polyhedral [J]. Computers and Graphics, 1997, 21(1): 23-30.

[5] 張宏, 溫永寧, 劉愛利. 地理信息系統算法基礎[M]. 北京: 科學出版社, 2006: 25-28.

[6] BERG M, CHEONG O, KREVELD M, et al. Computational geometry: Algorithms and applications [M]. Berlin: Springer-Verlag, 2008: 121-144.

[7] ZALIK B, CLAPWORTHY G J. A universal trapezoidation algorithm for planar polygons [J]. Computers and Graphics, 1999, 23(3): 353-363.

[8] LI J, WANG W C, WU E H. Point-in-polygon tests by convex decomposition [J]. Computers and Graphics, 2007, 31(4): 636-648.

[9] JIMéNEZ J J, FEITO F R, SEGURA R J. A new hierarchical triangle-based point-in-polygon data structure [J]. Computers and Geoscienes, 2009, 35(9): 1843-1853.

[10] ZALIK B, KOLINGEROVA I. A cell-based point-in-polygon algorithm suitable for large sets of points [J]. Computers and Geosciences, 2001, 27(10): 1135-1145.

[11] YANG S, YONG J H, SUN J, et al. A point-in-polygon method based on a quasi-closest point [J]. Computers and Geosciences, 2010, 36(2): 205-213.

[12] WANG W C, LI J, WU E H. 2D point-in-polygon test by classifying edges into layers [J]. Computers and Graphics, 2005, 29(3): 427-439.

[13] 李楠, 肖克炎. 一種改進的點在多邊形內外判斷算法[J]. 計算機工程, 2012, 38(5): 30-34.

[14] 孫愛玲, 趙光華, 趙敏華, 等. 基于sign(x)函數的點在多邊形內外判別算法及應用[J]. 計算機工程與科學, 2017, 39(4): 785-790.

[15] CLEARY J G, GEOFF W. Analysis of an algorithm for fast ray tracing using uniform space subdivision [J]. Visual Computer, 1988, 4(2): 65-83.

[16] LI J, WANG W C. Fast and robust GPU-based point-in-polyhedron determination [J]. Computer-Aided Design, 2017, 87: 20-28.

Efficient Point-in-Polygon Tests with Local and Simple Computation

WANG Sheng-chun1,2, WANG Wen-cheng1,2, TAN Xue-han1,2, LI Jing3

(1. State Key Laboratory of Computer Science, Institute of Software, The Chinese Academy of Sciences, Beijing 100190, China; 2. School of Computer Science and Technology, University of Chinese Academy of Sciences, Beijing 100190, China; 3. Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences, Beijing 100101, China)

polygon; grid; simple computation; point-in-polygon tests

TP 391

10.11996/JG.j.2095-302X.2019020267

A

2095-302X(2019)02-0267-07

2018-09-01；

2018-09-18

國家重點研發計劃項目(2017YFB1002700)；國家自然科學基金項目(61661146002，61872348)

王盛春(1991-)，男，陜西延安人，博士研究生。主要研究方向為計算機圖形學、紋理分析等。E-mail：wangsc@ios.ac.cn

王文成(1967-)，男，湖南雙峰人，研究員，博士，博士生導師。主要研究方向為計算機圖形學、虛擬現實及可視分析等。 E-mail：whn@ios.ac.cn