劉 京

(株洲市規(guī)劃設(shè)計(jì)院，湖南 株洲 412007)

自然界中可以用平面二維圖形表示的要素在GIS中大部分可以用多邊形來描述[1]，因此有關(guān)多邊形的操作也成為GIS專業(yè)軟件不可缺少的功能。多邊形的分割操作是與多邊形裁剪、合并、求交差密切相關(guān)，是關(guān)于多邊形運(yùn)算的一類算法[2]，目前在GIS運(yùn)用中越來越頻繁。典型的應(yīng)用有農(nóng)村土地承包經(jīng)營(yíng)權(quán)確權(quán)登記[3-5]、宗地分割、農(nóng)田整理規(guī)劃、礦山修復(fù)等。有關(guān)多邊形分割算法的研究主要集中在分割位置的計(jì)算和在GIS開發(fā)平臺(tái)下的實(shí)現(xiàn)。文獻(xiàn)[1]提出了一種基于二分法的簡(jiǎn)單平面多邊形快速分割算法。文獻(xiàn)[3]提出了基于面積差值的多邊形分割算法并應(yīng)用到了農(nóng)村土地確權(quán)登記中。文獻(xiàn)[4]則直接根據(jù)分割線和多邊形的關(guān)系推導(dǎo)了分割線交點(diǎn)位置的公式。文獻(xiàn)[6]提出了基于簡(jiǎn)單要素模型的算法，該算法采用掃描線及外包矩形的方法加快了分割相交線段的查找效率[6]。這些文獻(xiàn)最后都采用ArcGIS或者M(jìn)apGIS等GIS平臺(tái)的二次開發(fā)環(huán)境實(shí)現(xiàn)了算法。

由于AutoCAD具有友好的界面設(shè)計(jì)、便捷的交互模式和強(qiáng)大的制圖功能，國(guó)內(nèi)大部分測(cè)繪生產(chǎn)單位仍采用基于AutoCAD平臺(tái)開發(fā)的軟件作為測(cè)繪生產(chǎn)的主要工具，土地確權(quán)、宗地分割等項(xiàng)目也都以這些工具為平臺(tái)開展。因此，選擇一種合適的多邊形面積分割算法，并在AutoCAD的開發(fā)環(huán)境下實(shí)現(xiàn)該算法具有實(shí)際的意義。

1 多邊形面積分割算法選擇

1.1 多邊形面積分割算法選擇

文獻(xiàn)[3]提出的算法可以適用于任意的凸多邊形，但筆者在實(shí)現(xiàn)過程中發(fā)現(xiàn)了兩個(gè)問題：(1)如果初始分割子多邊形的面積和預(yù)期面積相差太大，分割方向線下一次掃描會(huì)超過分割多邊形最小外包矩形，并且不與分割多邊形的任何一條邊相交，計(jì)算出的目標(biāo)分割子多邊形面積恒為零，從而導(dǎo)致算法不收斂；(2)掃描線的步長(zhǎng)計(jì)算是基于梯形面積公式計(jì)算出來的，實(shí)際上當(dāng)前分割子多邊形與預(yù)期目標(biāo)分割子多邊形的差運(yùn)算得到的不一定是梯形，計(jì)算出來的步長(zhǎng)不一定是最優(yōu)的，從而減慢了算法收斂的速度。文獻(xiàn)[4]雖然提出了分割線交點(diǎn)位置的公式，可以直接實(shí)現(xiàn)簡(jiǎn)單多邊形按面積要求分割，但是該公式是基于CAD面積調(diào)整中調(diào)整一邊的方法設(shè)計(jì)的，且只能夠適用于分割目標(biāo)多邊形是梯形的情況；文獻(xiàn)[6]主要關(guān)注的是分割交點(diǎn)的計(jì)算及分割目標(biāo)子多邊形的構(gòu)造方法，用到了鏈表等復(fù)雜的數(shù)據(jù)結(jié)構(gòu)，不太適合于AutoCAD的VBA開發(fā)環(huán)境;文獻(xiàn)[1]的算法以嚴(yán)格的解析幾何為基礎(chǔ)，結(jié)合計(jì)算機(jī)領(lǐng)域廣泛應(yīng)用的二分查找算法，實(shí)現(xiàn)簡(jiǎn)單，應(yīng)用條件寬松，能夠滿足一般的項(xiàng)目要求。以上算法的原理、優(yōu)缺點(diǎn)及適用條件見表1。因此，本文以文獻(xiàn)[1]的算法為基礎(chǔ)，補(bǔ)充部分關(guān)鍵流程中需要用到的公式及計(jì)算方法，從底層的角度詳細(xì)論述算法在AutoCAD中實(shí)現(xiàn)的步驟。

表1 多邊形面積分割算法優(yōu)缺點(diǎn)對(duì)照表

1.2 多邊形面積分割二分法原理

算法的基本思想是首先輸入?yún)?shù)，包括要分割的多邊形、初始分割方向線、分割精度、預(yù)期分割面積，然后根據(jù)分割多邊形和初始分割方向線確定分割方向掃描線移動(dòng)的初始邊界條件，在邊界條件的范圍內(nèi)，逐步移動(dòng)分割方向掃描線，計(jì)算出掃描線與多邊形的交點(diǎn)，進(jìn)而計(jì)算出交點(diǎn)與被分割多邊形構(gòu)成的左邊(或右邊)分割子多邊形的面積，然后判斷該面積與預(yù)期分割面積的差值是否在分割精度范圍內(nèi)，如果是則取當(dāng)前分割掃描線為最終的分割方向線，否則就要修改分割掃描線移動(dòng)的邊界條件，如果分割面積比預(yù)期分割面積小，則修改下界為當(dāng)前分割掃描線，反之，則修改上界為當(dāng)前分割掃描線，如此不斷使用二分法的思想逐漸逼近零點(diǎn)，最終得到滿足分割精度的解。

2 AutoCAD中算法實(shí)現(xiàn)

Visual Basic for Application(VBA)是由微軟創(chuàng)建的用來自動(dòng)執(zhí)行任務(wù)的一個(gè)編程環(huán)境。它提供了一些用來創(chuàng)建圖形用戶界面的可拖拉工具和用來與AutoCAD對(duì)象交互的編程語言[7]。使用VBA with AutoCAD可以隨意定制專業(yè)級(jí)的CAD應(yīng)用程序。雖然目前CAD的主流開發(fā)方式逐漸由宿主式開發(fā)轉(zhuǎn)向獨(dú)立程序開發(fā)，但AutoCAD提供了一種可供開發(fā)者通過編程手段來操縱AutoCAD應(yīng)用程序的ActiveX機(jī)制，用戶可以自定義AutoCAD，與其他應(yīng)用程序共享圖形數(shù)據(jù)以及自動(dòng)完成任務(wù)[8]，這使得基于VBA開發(fā)的宿主式程序很容易向獨(dú)立的CAD程序遷移，而且AutoCAD允許VBA開發(fā)環(huán)境與AutoCAD同時(shí)運(yùn)行，使得算法運(yùn)行的每一步都能夠在CAD中以圖形化的方式可視化，便于直觀了解算法實(shí)際運(yùn)行效果。因此，本文選擇基于VBA的開發(fā)方式實(shí)現(xiàn)多邊形面積分割二分法算法。

2.1 基本解析方程

多邊形面積分割二分法是以解析幾何為基礎(chǔ)的的平行線掃描算法，不管是分割方向線掃描、多邊形面積計(jì)算還是掃描直線與多邊形邊求交過程都用到了直線方程，因此首先要定義各種直線方程，明確數(shù)學(xué)關(guān)系，才能實(shí)現(xiàn)算法各個(gè)關(guān)鍵步驟。掃描線與多邊形關(guān)系如圖1所示。

圖1 掃描線與多邊形關(guān)系示意圖

已知初始邊界點(diǎn)P1(x1,y1)、P4(x4,y4),直線P1P4與分割方向線AB的交點(diǎn)pt可表示：

pt(x)=(1-t)×x1+t×x4

pt(y)=(1-t)×y1+t×y4,t∈[0,1]

(1)

令a0=k,b0=-1,c0=pt(y)-k×pt(x),假設(shè)當(dāng)前分割線與邊P2P3相交，令a1=y2-y3，b1=x3-x2,c1=x2×y3-x3×y2，則分割線AB方程和邊P2P3方程可以可分別表示為：

a0×x+b0×y+c0=0

(2)

a1×x+b1×y+c1=0

(3)

在AutoCAD中對(duì)多邊形頂點(diǎn)和分割方向線的端點(diǎn)坐標(biāo)的提取可以采用文獻(xiàn)[9]所述的方法[9]。

2.2 初始邊界計(jì)算

二分法的核心思想是通過不斷迭代區(qū)間邊界逼近零點(diǎn)的過程來尋找問題的解。要進(jìn)行迭代求解，就必須要給出問題解的初始邊界區(qū)間。由文獻(xiàn)[1]可知初始邊界點(diǎn)與初始分割方向線在空間上距離最遠(yuǎn)。這里的最遠(yuǎn)是有方向性的，即分割線左邊的初始邊界點(diǎn)是所有位于分割線左邊的分割多邊形頂點(diǎn)中距離分割線最遠(yuǎn)的，反之亦然。由于初始分割線段的長(zhǎng)度是已知的，判斷距離的大小，可以直接轉(zhuǎn)化為判斷初始分割線段兩端點(diǎn)與分割多邊形頂點(diǎn)構(gòu)成的三角形面積的大小。AutoCAD中創(chuàng)建getBoundaryPoint函數(shù)計(jì)算初始邊界點(diǎn)部分代碼如下：

fp = xa * xb + xb * ny + nx * ya - xa * ny - xb * ya - nx * yb '三角形面積計(jì)算公式

If (fp > 0) and fp>max_s) Then '求位于分割線左邊的初始邊界點(diǎn)

max_s = fp:rst_pt(0).X = nx:rst_pt(0).Y = ny

End If

2.3 掃描直線與線段求交

由式(1)、式(2)可知，t每變化一次，分割線AB就會(huì)移動(dòng)掃描一次，掃描線與分割多邊形的邊相交會(huì)產(chǎn)生交點(diǎn)，要計(jì)算目標(biāo)分割子多邊形的面積首先就要先計(jì)算出分割交點(diǎn)的坐標(biāo)。由式(2)、式(3)可推導(dǎo)出分割掃描線AB與其中一條邊的交點(diǎn)坐標(biāo)公式如下：

(4)

由式(4)及數(shù)學(xué)知識(shí)可知，只要分割線AB與分割多邊形的邊不平行，則分割線在掃描移動(dòng)的過程中一定會(huì)跟所有的邊都有交點(diǎn)，只不過有的交點(diǎn)在邊上，有的在邊的延長(zhǎng)線上。而根據(jù)實(shí)際情況，只能選擇那些在邊上的交點(diǎn)作為當(dāng)前分割點(diǎn)。由幾何知識(shí)可知，在邊上的交點(diǎn)位于邊的兩個(gè)端點(diǎn)之間，其距離關(guān)系滿足如下的公式：

dik+djk-dij<ε

(5)

式(5)中,dik表示交點(diǎn)到邊上一個(gè)端點(diǎn)的距離，djk表示交點(diǎn)到邊上另外一個(gè)端點(diǎn)的距離，dij表示邊的長(zhǎng)度。ε是一個(gè)非常小的數(shù)字，在程序中可以設(shè)置為0.001。在AutoCAD中創(chuàng)建函數(shù)Intersect_Line_Segment和兩點(diǎn)間平面距離函數(shù)getDist，用于求掃描直線與分割多邊形邊上的交點(diǎn)。

2.4 構(gòu)造分割子多邊形

分割多邊形算法的最終目的是找到一個(gè)與預(yù)期面積逼近的子多邊形。得到分割子多邊形的基本思路是：判定交點(diǎn)所在的弧段，求出兩交點(diǎn)所載弧段的端點(diǎn)序列，然后構(gòu)造分割子多邊形[10]。由于分割線在每次掃描移動(dòng)中，會(huì)將多邊形虛擬的分割成左右兩個(gè)子多邊形，計(jì)算機(jī)是無法自動(dòng)選擇哪個(gè)方向的多邊形做為目標(biāo)多邊形的。因此需要人為的指定分割線左邊或者右邊的多邊形為目標(biāo)多邊形。不失一般性，假設(shè)選擇的目標(biāo)多邊形方向與初始邊界第一個(gè)頂點(diǎn)的方向一致，則在分割線掃描中，需要選擇多邊形上與初始邊界第一個(gè)頂點(diǎn)方向一致的點(diǎn)作為目標(biāo)多邊形的頂點(diǎn)。由此可以轉(zhuǎn)化為已知分割線的直線方程和位于直線一側(cè)的一個(gè)點(diǎn)，尋找多邊形上與這個(gè)點(diǎn)同側(cè)的頂點(diǎn)構(gòu)成分割子多邊形的頂點(diǎn)序列這一問題。已知直線的斜率k,直線上一點(diǎn)(xp,yp),則點(diǎn)(newx,newy)與直線的關(guān)系有：

a=k,b=-1,c=yp-k×xp

F=a×newx+b×newy+c

(6)

程序運(yùn)行時(shí)，先按式(6)計(jì)算初始邊界第一個(gè)頂點(diǎn)與當(dāng)前分割直線的值F1,然后逐個(gè)計(jì)算多邊形上頂點(diǎn)與當(dāng)前分割直線的值F2，如果F1與F2同號(hào)，則當(dāng)前頂點(diǎn)與初始邊界第一個(gè)頂點(diǎn)的方向相同。找到所有同方向的頂點(diǎn)，將其和分割線與多邊形邊的交點(diǎn)順序加入數(shù)組中，構(gòu)成目標(biāo)分割子多邊形。AutoCAD中創(chuàng)建getRetPolygon函數(shù)獲取目標(biāo)分割多邊形。部分關(guān)鍵代碼如下：

F1 = ComputeF(kab, xp, yp, xp1, yp1) '計(jì)算初始邊界第一個(gè)頂點(diǎn)與分割線的F值

For i = 0 To (UBound(polygon.Coordinates) - 1) / 2 - 1

F2 = ComputeF(kab, xp, yp, coords_x(i), coords_y(i))

If (F1 * F2 > 0) Then addVect left_x, left_y, coords_x(i), coords_y(i)

Next

2.5 二分法迭代求面積

由當(dāng)前分割子多邊形的頂點(diǎn)構(gòu)成的數(shù)組可以計(jì)算其面積，如果面積與預(yù)期面積的差值滿足精度要求，則停止計(jì)算，輸出當(dāng)前分割子多邊形；否則使用二分法修改邊界條件實(shí)現(xiàn)分割直線掃描移動(dòng)，再求交和計(jì)算新的面積直至面積差值達(dá)到精度。部分關(guān)鍵代碼如下：

Do '二分法逐步逼近要求面積

t = (sta_h + end_h) / 2 :xp = xp1 + t * (xp2 - xp1) :yp = yp1 + t * (yp2 - yp1)

s = getBRetPolygonArea(xp,yp,left_x, left_y) - rst '計(jì)算當(dāng)前面積與預(yù)期面積的差值

If Abs(s) < Sigma Then drawPoly left_x, left_y exit '達(dá)到精度要求則繪制分割子多邊形

If (s > 0) Then end_h = t '調(diào)整邊界值

If (s < 0) Then sta_h = t

j = j + 1

Loop While (Abs(s) > Sigma)

3 實(shí)驗(yàn)及結(jié)果

在AutoCAD環(huán)境下的VBA中設(shè)計(jì)了一個(gè)窗體用于進(jìn)行交互測(cè)試，如圖2所示，為了便于對(duì)比驗(yàn)證算法的可靠性，窗體除了實(shí)現(xiàn)了文獻(xiàn)[1]提出的算法還實(shí)現(xiàn)了文獻(xiàn)[3]的算法。實(shí)驗(yàn)中對(duì)一個(gè)多邊形以不同方向的分割線按照預(yù)期面積進(jìn)行分割，分割結(jié)果如表2 所示。從實(shí)驗(yàn)結(jié)果可以看出，對(duì)于一般情況下的簡(jiǎn)單多邊形，算法1和算法3都可以按精度要求實(shí)現(xiàn)面積的任意分割，但由于算法3對(duì)初始值

圖2 基于AutoCAD的簡(jiǎn)單多邊形面積分割算法實(shí)驗(yàn)

過于敏感，在某些分割方向下，當(dāng)初始分割面積與預(yù)期面積相差太大時(shí)，算法可能會(huì)出現(xiàn)不收斂的情況，從而導(dǎo)致分割失敗。

表2 面積分割實(shí)驗(yàn)結(jié)果表

4 結(jié) 論

有關(guān)文獻(xiàn)提出的多邊形按面積分割算法都是采用基于GIS平臺(tái)提供的二次開發(fā)接口實(shí)現(xiàn)的。本文以文獻(xiàn)[1]的算法為基礎(chǔ)，以嚴(yán)謹(jǐn)?shù)慕馕鰩缀芜\(yùn)算為依托，補(bǔ)充了部分關(guān)鍵步驟的計(jì)算方法，并采用AutoCAD的VBA開發(fā)環(huán)境，從底層的角度實(shí)現(xiàn)了多邊形面積分割的二分算法。通過對(duì)比試驗(yàn)，驗(yàn)證了算法的收斂性和可靠性。本文算法是采用VBA語言實(shí)現(xiàn)的，對(duì)于大量多邊形面積自動(dòng)分割操作而言，其效率會(huì)顯低下，以后可考慮將代碼移植到采用C++語言的ObjectARX環(huán)境下，以提升算法效率。