基于AutoCAD的多邊形面積分割算法實現(xiàn)

劉 京

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

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

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

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

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

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

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

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

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

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

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

2.1 基本解析方程

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

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

已知初始邊界點P1(x1,y1)、P4(x4,y4),直線P1P4與分割方向線AB的交點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中對多邊形頂點和分割方向線的端點坐標(biāo)的提取可以采用文獻(xiàn)[9]所述的方法[9]。

2.2 初始邊界計算

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

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

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

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

End If

2.3 掃描直線與線段求交

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

(4)

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

dik+djk-dij<ε

(5)

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

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

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

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

F=a×newx+b×newy+c

(6)

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

F1 = ComputeF(kab, xp, yp, xp1, yp1) '計算初始邊界第一個頂點與分割線的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)前分割子多邊形的頂點構(gòu)成的數(shù)組可以計算其面積，如果面積與預(yù)期面積的差值滿足精度要求，則停止計算，輸出當(dāng)前分割子多邊形；否則使用二分法修改邊界條件實現(xiàn)分割直線掃描移動，再求交和計算新的面積直至面積差值達(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 '計算當(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 實驗及結(jié)果

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

圖2 基于AutoCAD的簡單多邊形面積分割算法實驗

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

表2 面積分割實驗結(jié)果表

4 結(jié) 論

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