侯凯宇，徐景，郑衡，李佑林，周钰笛，李智文

长沙市规划勘测设计研究院, 长沙 410000

1 引 言

在测绘地理信息、城市规划等行业，经常需要在某一个封闭区域构建面要素（伍梅等，2005；梁晨等，2019)。例如，在控制性详细规划编制过程中，需要在路网合围区域绘制用地地块，人工绘制费时费力，迫切需要一种自动识别封闭区域并绘制面要素的方法。

封闭区域识别本质上是计算机数据结构中图的闭合回路搜索问题（Du 等，2020；曹嘉琛和杨武东，2022)。图是由若干给定的点及连接两点的线所构成的图形，分为有向图（倪文辉等，2018；金恒旭等，2022；崔建明等，2023)和无向图（胡荣明等，2014；张傲，2021；王国良和任雪玉，2023)。有向图表示连接两点的线有方向；无向图表示连接两点的线没有方向，封闭区域识别不需要方向。传统的闭合回路搜索算法有深度优先搜索和广度优先搜索，这两种算法原本的目的是遍历全图和寻找最短路径。刘旭（2017)为了解决经典正方化树图布局算法无法获得平均长宽比最优结果的问题，结合深度优先搜索技术，提出了基于深度优先搜索（depth first search，DFS)的正方化树图布局算法。杨雅涵（2021)针对联锁系统办理进路功能中存在的通用性较差、搜索效率较低等问题，提出了以深度优先搜索算法为基础的进路搜索方法，减少程序调试的过程，进而提高搜索程序的安全性与可靠性。曹嘉琛和杨武东（2022)为了提高计算机联锁系统的效率，改进传统进路搜索算法，提出了一种动态有向图的站场数据结构和改进的深度优先搜索算法。这些研究改进了深度优先搜索算法，优化了图的结构或提高了图的搜索效率，但搜索图的节点时是无序的，搜索路线无法围绕指定坐标前进。

本文研究一种封闭区域识别算法，只需要使用鼠标点击合围区域即可识别包含鼠标点击位置的最小封闭区域。首先，以鼠标点击位置为基点建立缓冲区，将缓冲区内的线和面在相交处打断以构建无向图和邻接表（殷茗等，2019)；其次，利用本文提出的改进的深度优先搜索算法搜索最小封闭区域，本文算法拓展了图的邻接表结构，并使用旋转角度控制邻接点的访问顺序，实现了搜索路线始终围绕指定坐标前进。

2 基于ArcGIS 对象模型构建无向图

ArcGIS Engine 是一个完整的嵌入式的GIS组件库，可以用于构建定制应用（牟乃夏等，2015；兰红等，2019；黄立鑫等，2021)。封闭区域识别算法的输入是任意的ArcGIS 线和面，输出是封闭区域合围形成的面。邻接表中存储了节点之间的ArcGIS 路径信息，路径由段组成，段是ArcGIS 最小的线单元，段可以是任意曲线。如果依次将封闭区域的路径连接到一起，就可以在ArcGIS Engine 开发的软件中构建任意曲线组成的面。

2.1 无向图

假设存在图G=〈V,E〉，其中，V是非空集合，称为节点集；E是V中元素构成的二元组的集合，称为边集。如果边没有方向，则称G是无向图。图的节点与节点之间是“多对多”的关系。在无向图中，关联一对节点的无向边如果多于1 条，则称这些边为平行边。既不含平行边也不包含自环的图，称为简单图（黄竞伟等，2000)，本文讨论的是简单无向图（simple undirected graph，SUG)。如图1所示，是一个含有4 个节点和4 条边的简单无向图。节点与节点之间的连接关系可以存储在邻接矩阵或邻接表中，本文选用邻接表。邻接表是单向链表的集合，链表的头节点为图的节点，链表的其他节点表示该节点连接的其他节点。例如，图1 的邻接表表示法，如图2 所示。

图1 简单无向图Fig.1 Simple undirected graph

图2 邻接表Fig.2 Adjacency list

2.2 扩展邻接表

本文构建的图属于稀疏图范畴，邻接表表示稀疏图比邻接矩阵更加简洁紧凑。在算法中，将邻接表节点的结构进行扩展，每个节点增加3 个属性，P表示父节点，用于提取回路路径；Color 表示节点颜色，用于标识该节点是否已经被访问；M表示相遇点，用于判断是否构成封闭区域。

2.3 ArcGIS 线和面对象模型

以路径为简单无向图的边，以路径的端点为图的节点。路径是连续的段的集合，段是线的最小组成单元；ArcGIS 中段有直线、圆弧、椭圆弧、贝泽尔曲线共四种。路径是线和面的基本组成单元，线由路径组成，面由闭合的路径—— 环组成（图3)。

图3 ArcGIS 线和面的对象模型Fig.3 ArcGIS object model for lines and surfaces

构建无向图的步骤如下：首先，通过鼠标点击屏幕的位置获得地理坐标，以该坐标建立缓冲区并获取缓冲区内的线和面；其次，将所有的线和面在交点处分割成多个路径（不考虑路径自交叉的情况)；最后，以路径的起点和终点为节点建立邻接表。实现求交点功能可以使用ArcGIS Engine 中的ITopologicalOperator 接口（牟乃夏等，2015)，也可以利用已有研究（Wei 等，2019；成雯，2021)算法求解交点。

3 基于改进的深度优先搜索的封闭区域识别

3.1 深度优先搜索

深度优先搜索算法被广泛应用于计算机图形学，用于遍历图的节点及寻找两个节点之间的最短路径。深度优先搜索采用了一种“一直向前走，走不通就掉头”的思想。从图中某节点S出发，先访问S的一个邻接点V，然后访问V的所有邻接点直到经过V这条路径的节点全部被访问，再退回来访问S的另一个邻接点，直至图中和S有路径相通的节点都被访问。如图4 所示，若以节点A为起点，该图的节点访问顺序有四种可能，分别是ABEC FHGD、ABEDGHFC、ACFHGDBE和ADGHF CBE。

图4 深度优先搜索Fig.4 Depth-first search

3.2 改进的深度优先搜索

对深度优先搜索算法加以改进，使之能够快速寻找包含指定坐标的最小封闭区域。算法的基本思路，如图5（a)所示，在无向图中寻找包含坐标点C的最小封闭区域。

图5 改进的深度优先搜索算法Fig.5 Improved depth-first search algorithm

3.3 选择搜索起点

搜索起点的选择直接影响到算法的执行时间，因为要寻找包含指定坐标点C的最小封闭区域，所以搜索起点选择距离C最近的节点。如果以该节点为起点进行搜索未能遍历图的所有节点且未能得到包含C的封闭回路，则选择未访问的节点中距离C最近的节点为起点。

3.4 起点邻接点的访问顺序

图5（a)中，选择距离坐标点C最近的节点V8为起点，V8的邻接点有3 个，分别是V5、V7和V9。如果优先访问V7，则一定无法得到包含C的最小封闭区域；如果优先访问V5或V9，则有可能得到包含C的最小封闭区域，因此需要限制邻接点的访问顺序。

在传统的深度优先搜索算法中，邻接点的访问顺序是无序的，是一种盲目搜索法。本文利用旋转角度解决邻接点访问顺序问题，使起点优先访问邻接点V5或V9。先求出直线V8C顺时针或逆时针旋转到V5、V7和V9的角度；根据旋转角度将起点的3个邻接点排序，再依次访问这些节点。无论采用顺时针旋转还是逆时针旋转都能保证V8首先访问V5或V9，本文采用逆时针旋转（图5（b))。

3.5 非起点邻接点的访问顺序

与起点邻接点的访问顺序类似，非起点邻接点的访问顺序也采用旋转角度对邻接点排序。定义当前节点V，V的父节点为P，计算以V为顶点，VP逆时针旋转至其他邻接点的角度。如图5（c)所示，V5V8逆时针旋转至V5V6、V5V2和V5V4，按旋转角度从小到大访问邻接点，那么下一个被访问的节点是V6；但访问V6的时候发现当前节点V6的邻接点只有父节点V5，所以原路返回并继续访问V5的下一个邻接点V2。以此类推，V2之后依次访问V3、V9和V8。

3.6 封闭区域提取

V9的邻接点中最先被访问的是V8，当访问V8时发现它的颜色是黑色，说明该节点已经被访问过了，从而构成了封闭区域V8V5V2V3V9（图5（c))。算法结束前要在V9的相遇点信息中记录节点V8以便于后续的封闭区域提取。无向图中仅V9含有相遇点信息，因此，以V9为起点提取封闭区域，V9的父节点为V3，V3不是V9的相遇点，记录V9到V3的路径信息；V3的父节点为V2，V2也不是V9的相遇点，记录V3到V2的路径信息；同理，记录V2到V5的路径；直到发现V8是V9的相遇点，记录V5到V8和V8到V9的路径后封闭区域提取完成。此时所有记录的路径信息构成了环V8V5V2V3V9，ArcGIS 中环可以构建面。

3.7 非直线边的旋转角度

当2 个节点之间的边为曲线时，计算旋转角度不能直接使用2 个节点的连线进行旋转，应当使用过当前节点的切线。如图6 所示，假设指定坐标点位于最小封闭区域PACD内，节点A为当前点，节点P、B、C为A的邻接点，其中，P为A的父节点。按照本文定义的旋转角度计算方法，以A为顶点，边PA逆时针旋转至AC和AB的角度分别为α和β，角β小于角α，导致先访问B再访问C，最终得到一个非最小封闭区域PABCD。显然，PA逆时针旋转会先遇到边AC再遇到边AB，过节点A作边AC的切线AA′，PA旋转至AA′的角度γ才是正确的旋转角，角γ小于角β，则先访问C，最终得到正确的最小封闭区域PACD。如果路径A到C是一个混合了直线和各种曲线的复杂路径，那么取A到C路径上的第一个线段计算旋转角度。

图6 非直线边的旋转角度Fig.6 Rotation angle of non-straight edges

ArcGIS Engine 中提供了丰富的曲线接口，如圆弧接口（ICircularArc)、椭圆弧接口（IEllipticArc)、贝泽尔曲线接口（IBezierCurve)等，由此可以方便的求取曲线上某一点的切线。以圆弧为例，图7 是ArcGIS 中定义的圆弧属性，已知圆弧的起点、终点和圆心，易求得过点起点的切线，也可以使用曲线接口中的QueryTangent 方法求取切线。

图7 ArcGIS 圆弧属性Fig.7 ArcGIS circular arc properties

3.8 算法步骤

算法步骤如下所述。

（1)以指定坐标点C为基点建立缓冲区，缓冲区内的线、面要素的路径添加到集合List＜IPath＞。

（2)对集合List＜IPath＞中的所有元素求交点（利用ArcGIS Engine 中的ITopologicalOperator 接口)，并在交点处打断，得到简单无向图。

（3)遍历图中的边以建立邻接表，初始化表中的所有节点，颜色Color=白色，父节点P=空，相遇点M=空。

（4)计算图中所有节点到C点的距离，选择距离最小的节点作为深度优先搜索的起点S。

（5)开始访问当前节点S，修改S的颜色Color[S]=黑色。

（6)计算S的旋转角度：以S为顶点，直线SC逆时针旋转到S的其他邻接点的角度（如果S到邻接点的边是曲线，则计算SC逆时针旋转到S至邻接点路径上第一个段的切线，下同)。

（7)按照旋转角度从小到大对S的邻接点排序。

（8)然后拿到S的一个邻接点V。

（9)开始访问当前节点V，V的父节点P[V]=S，颜色 Color[V]=黑色。

（10)计算以V为顶点，VS逆时针旋转至V的邻接点的旋转角度。

（11)按照旋转角度从小到大对V的邻接点排序。

（12)然后拿到V的邻接点N。

（13)如果N不是V的父节点，且N未被访问过，即 Color[N]=白色，将N设为当前节点，跳到步骤（9)。

（14)如果N不是V的父节点，但N被访问过，即 Color[N]=黑色，记录相遇点M[V]=N，跳到步骤（19)。

（15)如果N是V的父节点，说明V的邻接点访问完毕（父节点的旋转角度为360°，排在其他所有邻接点之后)，将V的父节点设为当前节点，跳到步骤（9)。

（16)至此，与起点S相连通的所有节点都被访问了，但仍然没有找封闭区域，考虑无向图中是否存在其他的连通分量（谭屯子等，2018；廖小飞等，2019)。

（17)遍历图中所有节点，如果发现节点V的颜色 Color[V]=白色，将V设为当前节点，跳到步骤（5)。

（18)如果图中所有节点均已被访问，表示未寻找到包含C点的封闭区域，程序结束。

（19)开始提取封闭区域，拿到含有相遇点信息的节点E，E的相遇点为M[E]。

（20)获取当前节点V的父节点P[V]。

（21)根据邻接表获取节点V到节点P[V]的路径信息，并记录到路径集合List＜IPath＞。

（22)如果P[V] =M[E]，跳到步骤（24)。

（23)如果P[V] ≠M[E]，将P[V]设为当前节点，跳到步骤（20)。

（24)根据邻接表获取节点M[E]到节点E的路径信息，并记录到路径集合List＜IPath＞。

（25)根据路径集合List＜IPath＞构建ArcGIS 面几何（利用ArcGIS Engine 中的IGeometryCollection接口依次添加路径集合List＜IPath＞中的路径)。

（26)程序结束。

4 实 验

图8（a)是由12 个节点组成的简单无向图，其中P点表示搜索基点；图8（b)中圆形虚线表示缓冲区范围，与缓冲区相交或被缓冲区包含的路径参与简单无向图的构建，并在图中高亮显示这些路径。如果没有找到封闭区域则扩大缓冲区范围，重新构建无向图，直至找到封闭区域或缓冲区达到阈值，阈值的设定可采用如下三种方式：①设置构建缓冲区次数的上限；②设置缓冲区面积的上限；③设置缓冲区内路径数量的上限。

图8 实验Fig.8 Experiments

以P为基点建立缓冲区（图8（b))，则参与构建简单无向图的路径有AB、AD、AE、AF，其中，A点到P点的距离最近，则无向图的当前节点为A。根据前文所述知图中节点的访问顺序为ADEFB，节点D、E、F、B除父节点外，未找到其他被访问过的邻接点，所以封闭区域搜索失败。

扩大缓冲区半径（图8（c))，则参与构建简单无向图的路径有AB、AD、AE、AF、DG、DH、DC、ab、bc、cd、da。图中节点的访问顺序为ADCGHEFB，同理封闭区域搜索失败。此时，图中还有abcd节点未被访问，这表示abcd与其他节点不连通，其中，距离P点最近的节点为a，节点的访问顺序为adcb，节点b除父节点c外，还有邻接点a被访问过，所以得到封闭区域abcd，但由于该封闭区域未包含P点，所以封闭区域搜索失败。

继续扩大缓冲区半径（图8（d))，则参与构建简单无向图的路径有AB、BC、CD、DA、AE、AF、DG、DH、ab、bc、cd、da，图中节点的访问顺序为ADCB，可得到封闭区域ABCD且该封闭区域包含P点。注意，虽然封闭区域abcd未包含P点，但abcd被ABCD包含，所以ABCD应对abcd进行面内构岛，否则ABCD就不是最小封闭区域。图8（e)所示的黄色区域即为包含P点的最小封闭区域，在ArcGIS Engine 中利用IGeometry Collection 接口依次添加面ABCD和面abcd即可得到该最小封闭区域构成的面。

5 结 论

本文提出一种改进的深度优先搜索算法，用于搜索包含指定坐标点的最小封闭区域。在传统深度优先搜索算法的基础上，利用旋转角度限制邻接点的访问顺序，使得搜索路线始终围绕指定坐标前进，又因为邻接表存储了节点之间的ArcGIS 路径，可以方便地构建任意面要素。本文算法已应用于ArcGIS Engine 开发的控规软件（图9)，利用点选功能（图10)可以快速绘制规划地块，使用鼠标点击线与面构成的合围区域，即可搜索包含鼠标点击位置的最小封闭区域并创建规划地块面要素。

图9 软件菜单Fig.9 Software menu

图10 绘制用地Fig.10 Draw land

在实现算法的过程中，搜索算法本身有着严密的逻辑，而简单无向图的构建反而容易成为搜索失败的原因。理由有三点：一是地图中的线和面有着复杂的空间关系，如空间关系相等的要素要剔除；二是地图中的线和面不可避免地存在拓扑错误，如自相交需要提前判断；三是因为简单无向图是一种理想的图，它要求节点没有自环且节点之间有唯一的路径相连，这在实际情况中无法得到保证。综上，在执行封闭区域搜索算法前，需要对参与构建简单无向图的线和面进行合理的分析和预处理，才能保证搜索算法的顺利执行和面要素的成功构建。