基于AutoCAD的多边形面积分割算法实现

刘 京

(株洲市规划设计院，湖南 株洲 412007)

自然界中可以用平面二维图形表示的要素在GIS中大部分可以用多边形来描述[1]，因此有关多边形的操作也成为GIS专业软件不可缺少的功能。多边形的分割操作是与多边形裁剪、合并、求交差密切相关，是关于多边形运算的一类算法[2]，目前在GIS运用中越来越频繁。典型的应用有农村土地承包经营权确权登记[3-5]、宗地分割、农田整理规划、矿山修复等。有关多边形分割算法的研究主要集中在分割位置的计算和在GIS开发平台下的实现。文献[1]提出了一种基于二分法的简单平面多边形快速分割算法。文献[3]提出了基于面积差值的多边形分割算法并应用到了农村土地确权登记中。文献[4]则直接根据分割线和多边形的关系推导了分割线交点位置的公式。文献[6]提出了基于简单要素模型的算法，该算法采用扫描线及外包矩形的方法加快了分割相交线段的查找效率[6]。这些文献最后都采用ArcGIS或者MapGIS等GIS平台的二次开发环境实现了算法。

由于AutoCAD具有友好的界面设计、便捷的交互模式和强大的制图功能，国内大部分测绘生产单位仍采用基于AutoCAD平台开发的软件作为测绘生产的主要工具，土地确权、宗地分割等项目也都以这些工具为平台开展。因此，选择一种合适的多边形面积分割算法，并在AutoCAD的开发环境下实现该算法具有实际的意义。

1 多边形面积分割算法选择

1.1 多边形面积分割算法选择

文献[3]提出的算法可以适用于任意的凸多边形，但笔者在实现过程中发现了两个问题：(1)如果初始分割子多边形的面积和预期面积相差太大，分割方向线下一次扫描会超过分割多边形最小外包矩形，并且不与分割多边形的任何一条边相交，计算出的目标分割子多边形面积恒为零，从而导致算法不收敛；(2)扫描线的步长计算是基于梯形面积公式计算出来的，实际上当前分割子多边形与预期目标分割子多边形的差运算得到的不一定是梯形，计算出来的步长不一定是最优的，从而减慢了算法收敛的速度。文献[4]虽然提出了分割线交点位置的公式，可以直接实现简单多边形按面积要求分割，但是该公式是基于CAD面积调整中调整一边的方法设计的，且只能够适用于分割目标多边形是梯形的情况；文献[6]主要关注的是分割交点的计算及分割目标子多边形的构造方法，用到了链表等复杂的数据结构，不太适合于AutoCAD的VBA开发环境;文献[1]的算法以严格的解析几何为基础，结合计算机领域广泛应用的二分查找算法，实现简单，应用条件宽松，能够满足一般的项目要求。以上算法的原理、优缺点及适用条件见表1。因此，本文以文献[1]的算法为基础，补充部分关键流程中需要用到的公式及计算方法，从底层的角度详细论述算法在AutoCAD中实现的步骤。

表1 多边形面积分割算法优缺点对照表

1.2 多边形面积分割二分法原理

算法的基本思想是首先输入参数，包括要分割的多边形、初始分割方向线、分割精度、预期分割面积，然后根据分割多边形和初始分割方向线确定分割方向扫描线移动的初始边界条件，在边界条件的范围内，逐步移动分割方向扫描线，计算出扫描线与多边形的交点，进而计算出交点与被分割多边形构成的左边(或右边)分割子多边形的面积，然后判断该面积与预期分割面积的差值是否在分割精度范围内，如果是则取当前分割扫描线为最终的分割方向线，否则就要修改分割扫描线移动的边界条件，如果分割面积比预期分割面积小，则修改下界为当前分割扫描线，反之，则修改上界为当前分割扫描线，如此不断使用二分法的思想逐渐逼近零点，最终得到满足分割精度的解。

2 AutoCAD中算法实现

Visual Basic for Application(VBA)是由微软创建的用来自动执行任务的一个编程环境。它提供了一些用来创建图形用户界面的可拖拉工具和用来与AutoCAD对象交互的编程语言[7]。使用VBA with AutoCAD可以随意定制专业级的CAD应用程序。虽然目前CAD的主流开发方式逐渐由宿主式开发转向独立程序开发，但AutoCAD提供了一种可供开发者通过编程手段来操纵AutoCAD应用程序的ActiveX机制，用户可以自定义AutoCAD，与其他应用程序共享图形数据以及自动完成任务[8]，这使得基于VBA开发的宿主式程序很容易向独立的CAD程序迁移，而且AutoCAD允许VBA开发环境与AutoCAD同时运行，使得算法运行的每一步都能够在CAD中以图形化的方式可视化，便于直观了解算法实际运行效果。因此，本文选择基于VBA的开发方式实现多边形面积分割二分法算法。

2.1 基本解析方程

多边形面积分割二分法是以解析几何为基础的的平行线扫描算法，不管是分割方向线扫描、多边形面积计算还是扫描直线与多边形边求交过程都用到了直线方程，因此首先要定义各种直线方程，明确数学关系，才能实现算法各个关键步骤。扫描线与多边形关系如图1所示。

图1 扫描线与多边形关系示意图

已知初始边界点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),假设当前分割线与边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中对多边形顶点和分割方向线的端点坐标的提取可以采用文献[9]所述的方法[9]。

2.2 初始边界计算

二分法的核心思想是通过不断迭代区间边界逼近零点的过程来寻找问题的解。要进行迭代求解，就必须要给出问题解的初始边界区间。由文献[1]可知初始边界点与初始分割方向线在空间上距离最远。这里的最远是有方向性的，即分割线左边的初始边界点是所有位于分割线左边的分割多边形顶点中距离分割线最远的，反之亦然。由于初始分割线段的长度是已知的，判断距离的大小，可以直接转化为判断初始分割线段两端点与分割多边形顶点构成的三角形面积的大小。AutoCAD中创建getBoundaryPoint函数计算初始边界点部分代码如下：

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就会移动扫描一次，扫描线与分割多边形的边相交会产生交点，要计算目标分割子多边形的面积首先就要先计算出分割交点的坐标。由式(2)、式(3)可推导出分割扫描线AB与其中一条边的交点坐标公式如下：

(4)

由式(4)及数学知识可知，只要分割线AB与分割多边形的边不平行，则分割线在扫描移动的过程中一定会跟所有的边都有交点，只不过有的交点在边上，有的在边的延长线上。而根据实际情况，只能选择那些在边上的交点作为当前分割点。由几何知识可知，在边上的交点位于边的两个端点之间，其距离关系满足如下的公式：

dik+djk-dij<ε

(5)

式(5)中,dik表示交点到边上一个端点的距离，djk表示交点到边上另外一个端点的距离，dij表示边的长度。ε是一个非常小的数字，在程序中可以设置为0.001。在AutoCAD中创建函数Intersect_Line_Segment和两点间平面距离函数getDist，用于求扫描直线与分割多边形边上的交点。

2.4 构造分割子多边形

分割多边形算法的最终目的是找到一个与预期面积逼近的子多边形。得到分割子多边形的基本思路是：判定交点所在的弧段，求出两交点所载弧段的端点序列，然后构造分割子多边形[10]。由于分割线在每次扫描移动中，会将多边形虚拟的分割成左右两个子多边形，计算机是无法自动选择哪个方向的多边形做为目标多边形的。因此需要人为的指定分割线左边或者右边的多边形为目标多边形。不失一般性，假设选择的目标多边形方向与初始边界第一个顶点的方向一致，则在分割线扫描中，需要选择多边形上与初始边界第一个顶点方向一致的点作为目标多边形的顶点。由此可以转化为已知分割线的直线方程和位于直线一侧的一个点，寻找多边形上与这个点同侧的顶点构成分割子多边形的顶点序列这一问题。已知直线的斜率k,直线上一点(xp,yp),则点(newx,newy)与直线的关系有：

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

F=a×newx+b×newy+c

(6)

程序运行时，先按式(6)计算初始边界第一个顶点与当前分割直线的值F1,然后逐个计算多边形上顶点与当前分割直线的值F2，如果F1与F2同号，则当前顶点与初始边界第一个顶点的方向相同。找到所有同方向的顶点，将其和分割线与多边形边的交点顺序加入数组中，构成目标分割子多边形。AutoCAD中创建getRetPolygon函数获取目标分割多边形。部分关键代码如下：

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 二分法迭代求面积

由当前分割子多边形的顶点构成的数组可以计算其面积，如果面积与预期面积的差值满足精度要求，则停止计算，输出当前分割子多边形；否则使用二分法修改边界条件实现分割直线扫描移动，再求交和计算新的面积直至面积差值达到精度。部分关键代码如下：

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 '计算当前面积与预期面积的差值

If Abs(s) < Sigma Then drawPoly left_x, left_y exit '达到精度要求则绘制分割子多边形

If (s > 0) Then end_h = t '调整边界值

If (s < 0) Then sta_h = t

j = j + 1

Loop While (Abs(s) > Sigma)

3 实验及结果

在AutoCAD环境下的VBA中设计了一个窗体用于进行交互测试，如图2所示，为了便于对比验证算法的可靠性，窗体除了实现了文献[1]提出的算法还实现了文献[3]的算法。实验中对一个多边形以不同方向的分割线按照预期面积进行分割，分割结果如表2 所示。从实验结果可以看出，对于一般情况下的简单多边形，算法1和算法3都可以按精度要求实现面积的任意分割，但由于算法3对初始值

图2 基于AutoCAD的简单多边形面积分割算法实验

过于敏感，在某些分割方向下，当初始分割面积与预期面积相差太大时，算法可能会出现不收敛的情况，从而导致分割失败。

表2 面积分割实验结果表

4 结 论

有关文献提出的多边形按面积分割算法都是采用基于GIS平台提供的二次开发接口实现的。本文以文献[1]的算法为基础，以严谨的解析几何运算为依托，补充了部分关键步骤的计算方法，并采用AutoCAD的VBA开发环境，从底层的角度实现了多边形面积分割的二分算法。通过对比试验，验证了算法的收敛性和可靠性。本文算法是采用VBA语言实现的，对于大量多边形面积自动分割操作而言，其效率会显低下，以后可考虑将代码移植到采用C++语言的ObjectARX环境下，以提升算法效率。