CN110008837B - 一种根据所属道路网络封闭多边形统计poi的快速实现方法 - Google Patents

Abstract

一种根据所属道路网络封闭多边形统计POI的快速实现方法，首先根据实际道路网络的拓扑关系按比例输出道路网络图像，明确实际平面坐标范围，利用图像处理将道路网络图像转化为二值像素点矩阵，其次基于4邻接的邻接关系识别图像中道路网络分割的封闭多边形，并为属于同一个封闭多边形的像素点赋予唯一标记数值，然后，对POI点的实际坐标进行像素坐标转化，得到POI在道路网络图像中的坐标值，根据该坐标像素点取值标记POI，最后，根据POI所属封闭多边形标记值统计汇总道路网络分割的封闭多边形内部POI。本发明基于封闭多边形内部像素点的连通性，将点与多边形位置关系的判断转化为像素点间连通关系的识别，在多边形和POI集较大情况下实现位置关系快速判断。

Description

一种根据所属道路网络封闭多边形统计POI的快速实现方法

技术领域

本发明涉及道路网络中几何元素空间位置关系识别领域，特别是涉及在大规模道路网络及POI情况下快速识别POI点与任意形状封闭多边形位置关系的技术，属于一种根据所属道路网络封闭多边形统计POI的快速实现方法。

背景技术

POI(Point of Interest，兴趣点)是将与人们生活密切相关的地理实体抽象为点的表达结果。在交通规划领域，由于具有数据可自动获取、更新速度快、覆盖面广等优势，POI广泛用于用地功能识别、交通区划、人口密度关联性等分析，具有较高的实际研究价值。

在上述提到的研究中，常常需要以人为分割或道路网络分割的封闭多边形为单元分析POI，为保证分析精度，对于一个大型研究区域，其分割形成的封闭多边形数量级达千级，POI数量级也达到十万级以上。根据实际分析需要，基于道路网络分割的封闭多边形常常不是规则的凸多边形，而是任意形状的封闭多边形，随着道路网络规模的增大和POI数量级的增加，POI所属封闭多边形的识别效率逐渐降低。

发明内容

为了解决以上问题，本发明提供一种根据所属道路网络封闭多边形统计POI的快速实现方法，为提高POI所属封闭多边形的识别速度，本发明将点与多边形的位置关系识别转化为像素点间连通性的识别，为达此目的，本发明提一种根据所属道路网络封闭多边形统计POI的快速实现方法，包括如下步骤，其特征在于：

步骤(1)根据城市道路网络拓扑关系按比例输出道路网络图像，确定图像覆盖的实际平面坐标范围：(x_min,y_min)、(x_max,y_max)，单位为m；图像分辨率为m×n，单位为像素，权衡分析精度和速度，建议比例范围：(x_max-x_min):m∈[10,25]；

步骤(2)依次对图像进行顺时针旋转90°、灰度化、二值化处理，得到二值像素点矩阵P₁，其中取值为0的像素点表示为非道路，取值为1的像素点为道路，便于邻接关系识别，将矩阵扩展为(m+2)×(n+2)的新像素点矩阵P₂，扩展公式如下：

步骤(3)遍历矩阵P₂中位于1至m+1行且位于1至n+1列的像素点，对于标记为0的非道路像素点基于4邻接的邻接关系识别连通区域，并为属于同一个连通区域的像素点赋予唯一标记值，将经过连通区域标记的像素矩阵记为P₃；

步骤(4)遍历矩阵P₃第2行、第m+1行、第2列、第n+1列，统计上述像素点中标记为连通区域的像素点取值，记为outer＝[poly₁,poly₂,…,poly_a]，记录位于道路网络覆盖范围外的连通区域标记值，其中a为位于道路网络覆盖范围外的连通区域数量；

步骤(5)遍历POI＝[POI₁,POI₂,…POI_s]，将POI_i实际平面坐标(x_i,y_i)转化为像素坐标(r_i,c_i)，转化公式如下：

其中，ceil表示向上取整；若r_i＜1、r_i＞m+1、c_i＜1、c_i＞n+1中任意表达式成立，或P₃[r_i,c_i]∈outer，该POI位于道路网络覆盖范围外，将该POI从集合中剔除，否则为POI_i构建所属连通区域及道路网络分割的封闭多边形的索引POI_i:P₃[r_i,c_i]；

步骤(6)根据POI与连通区域的索引关系，汇总各封闭多边形内部的POI，即将所属封闭多边形标记值相同的POI归为该封闭多边形的内部POI。

作为本发明进一步改进，所述步骤(1)中，为便于后续图像处理，根据实际道路网络拓扑关系输出的图像需要保证道路像素点与非道路像素点的RGB取值范围无交叉，并且横轴、纵轴按照相同的比例输出，即(x_max-x_min):m＝(y_max-y_min):n。

作为本发明进一步改进，所述步骤(3)中，识别全部连通区域包含以下步骤：

步骤(3.1)令r＝2，c＝2，count＝2，P₃＝P₂；

步骤(3.2)若r＝m+2，结束步骤(3)，全部连通区域均以识别，输出连通区域标记结果P₃；否则，执行步骤(3.3)；

步骤(3.3)若P₃[r,c]≠0，执行步骤(3.6)；否则，令temp1＝[(r,c)]，P₃[r,c]＝count，统计尚未标记的且取值为0的像素点坐标，执行步骤(3.4)；

步骤(3.4)令temp2＝[]，统计temp1中的像素点邻接的标记为0的像素点坐标，邻接关系选择4邻接，具体做法为：遍历temp1，对任意(r_j,c_j)∈temp1，提取像素点(r_j,c_j)上、下、左、右四个像素点，即(r_j-1,c_j)、(r_j+1,c_j)、(r_j,c_j-1)、(r_j,c_j+1)，若P₃[r_j-1,c_j]＝0成立，将(r_j-1,c_j)追加写入temp2，并令P₃[r_j-1,c_j]＝count，对其余邻接像素点执行相同的判断和操作；

步骤(3.5)令temp1＝temp2，若

count＝count+1，执行步骤(3.6)；否则，执行步骤(3.4)；

步骤(3.6)令c＝c+1，若c＝n+2，r＝r+1；执行步骤(3.2)。

作为本发明进一步改进，所述步骤(4)中，识别道路网络覆盖范围外的连通区域包含以下步骤：

步骤(4.1)令r＝2，c＝2，outer＝[]；

步骤(4.2)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.3)令c＝c+1，若c＝n+2，执行步骤(4.4)；否则，执行步骤(4.2)；

步骤(4.4)令r＝m+1，c＝2；

步骤(4.5)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.6)令c＝c+1，若c＝n+2，执行步骤(4.7)；否则，执行步骤(4.5)；

步骤(4.7)令r＝2，c＝2；

步骤(4.8)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.9)令r＝r+1，若r＝m+2，执行步骤(4.10)；否则，执行步骤(4.8)；

步骤(4.10)令r＝2，c＝n+1；

步骤(4.11)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.12)令r＝r+1，若r＝m+2，结束步骤(4)，道路网络覆盖范围外的连通区域均以识别，输出识别结果outer＝[poly₁,poly₂,…,poly_a]；否则，执行步骤(4.11)。

本发明提供一种根据所属道路网络封闭多边形统计POI的快速实现方法，首先根据实际道路网络的拓扑关系按比例输出道路网络图像，明确实际平面坐标范围，利用图像处理将道路网络图像转化为二值像素点矩阵，其次基于4邻接的邻接关系识别图像中道路网络分割的封闭多边形，并为属于同一个封闭多边形的像素点赋予唯一标记数值，然后，对POI点的实际坐标进行像素坐标转化，得到POI在道路网络图像中的坐标值，根据该坐标像素点取值标记POI，最后，根据POI所属封闭多边形标记值统计汇总道路网络分割的封闭多边形内部POI。

有益效果：在针对大型道路网络的分析背景下，该方法可以提高POI所属道路网络分割多边形的识别速度，由于该方法基于连通区域标记识别点与封闭多边形的位置关系，对任意方式划分和任意形状的封闭多边形均适用。

附图说明

图1为本发明根据所属道路网络封闭多边形统计POI的快速实现方法的流程图；

图2为本发明选用的城市道路网络及其边界示意图；

图3为本发明测试POI在实际地图中的分布示意图；

图4为本发明剔除道路网络覆盖范围外POI前后的POI分布对比示意图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

本发明提供一种根据所属道路网络封闭多边形统计POI的快速实现方法，，为提高POI所属封闭多边形的识别速度，本发明将点与多边形的位置关系识别转化为像素点间连通性的识别。

图1是实现一种根据所属道路网络封闭多边形统计POI的快速实现方法的流程图，下面将结合实例对本发明作进一步阐述：

步骤(1)根据城市道路网络拓扑关系按比例输出道路网络图像，确定图像覆盖的实际平面坐标范围：(x_min,y_min)、(x_max,y_max)，单位为m；图像分辨率为m×n，单位为像素，权衡分析精度和速度，建议比例范围：(x_max-x_min):m∈[10,25]；

本实施实例中，选择南京市主城区为测试对象，基于OSM获取道路网络拓扑关系，利用python的绘图功能输出道路网络(含边界)图像如图2所示。该图像覆盖的实际平面坐标范围为：(x_min,y_min)＝(649693.187,3517677.675)、(x_max,y_max)＝(690498.259,3566738.156)，单位为m；图像分辨率为1632×1962，单位为像素，实际距离与像素距离的比例为25，即一个像素点表示一个25×25m²的区域。

步骤(2)依次对图像进行顺时针旋转90°、灰度化、二值化处理，得到二值像素点矩阵P₁，其中取值为0的像素点表示为非道路，取值为1的像素点为道路，便于邻接关系识别，将矩阵扩展为(m+2)×(n+2)的新像素点矩阵P₂，扩展公式如下：

本实施实例中，使用python的PIL库进行图像的旋转、灰度化、二值化处理，其中二值化的全局阈值取150，得到1632×1962的二值像素点矩阵P₁，经公式(1)的扩展，得到1634×1964的二值像素点矩阵P₂，表1为从P₂中截取的部分像素矩阵。

表1P₂部分像素取值示例

步骤(3)遍历矩阵P₂中位于1至m+1行且位于1至n+1列的像素点，对于标记为0的非道路像素点基于4邻接的邻接关系识别连通区域，并为属于同一个连通区域的像素点赋予唯一标记值，将经过连通区域标记的像素矩阵记为P₃；

本实施实例中，以表1中标记为1的道路像素点围合的封闭多边形识别为例说明本步骤的实施过程。

a.已知当前识别的像素点是(734,1163)，count＝871，根据步骤(3.3)，P₃[734,1163]＝0，令temp2＝[(734,1163)]，P₃[734,1163]＝871；

b.根据步骤(3.4)，令temp2＝[]，遍历中temp1中像素点，标记其具有4邻接关系且取值为0的邻接像素点：像素点(734,1163)在矩阵中具有4邻接关系的且取值为0的像素点有(735,1163)、(734,1164)，将(735,1163)、(734,1164)追加写入temp2得到temp2＝[(735,1163),(734,1164)]，并令P₃[735,1163]＝871，P₃[734,1164]＝871；

c.根据步骤(3.5)，令temp1＝temp2＝[(735,1163),(734,1164)]，由于

继续执行步骤(3.4)；

d.根据步骤(3.4)，同理，令temp2＝[]，像素点(735,1163)、(734,1164)在矩阵中具有4邻接关系的且取值为0的像素点有(736,1163)、(735,1164)、(734,1165)，将上述像素点追加写入temp2得到temp2＝[(736,1163)、(735,1164)、(734,1165)]，并令P₃[736,1163]＝871，P₃[735,1164]＝871，P₃[734,1165]＝871；

e.根据步骤(3.5)令temp1＝temp2＝[(736,1163)、(735,1164)、(734,1165)]，由于

继续执行步骤(3.4)，后续过程以此类推，直至

完成一个连通区域的标记，该区域标记结果如表3所示。

表2像素点4邻接关系及标记示例

表3部分连通区域标记结果示例

步骤(4)遍历矩阵P₃第2行、第m+1行、第2列、第n+1列，统计上述像素点中标记为连通区域的像素点取值，记为outer＝[poly₁,poly₂,…,poly_a]，记录位于道路网络覆盖范围外的连通区域标记值，其中a为位于道路网络覆盖范围外的连通区域数量；

本实施实例中，遍历矩阵P₃第2行、第1633行、第2列、第1963列，统计大于等于2的取值，得到位于道路网络覆盖范围外的连通区域像素点取值outer＝[2,1008,3,959,1061]。

步骤(5)遍历POI＝[POI₁,POI₂,…POI_s]，将POI_i实际平面坐标(x_i,y_i)转化为像素坐标(r_i,c_i)，转化公式如下：

其中，ceil表示向上取整；若r_i＜1、r_i＞m+1、c_i＜1、c_i＞n+1中任意表达式成立，或P₃[r_i,c_i]∈outer，该POI位于道路网络覆盖范围外，将该POI从集合中剔除，否则为POI_i构建所属连通区域及道路网络分割的封闭多边形的索引POI_i:P₃[r_i,c_i]；

本实施实例中，所用测试POI如表4所示，图3展示了测试POI在地图中的实际分布情况。以POI“retrogallery”为例，实际平面坐标转为像素坐标的计算过程如下：

表4测试POI识别结果

名称	平面横坐标	平面纵坐标	像素横坐标	像素纵坐标	所属多边形标记值
						retrogallery	667993.549	3546775.399	733	1165	871
贝尔尼尼	668070.324	3546825.468	736	1167	871
						jorya	668065.494	3546819.358	736	1167	871
BORA AKSU店铺	668048.603	3546802.795	736	1166	871
						金鹰	668041.642	3547690.156	735	1202	864
Melvita	668018.877	3547656.396	734	1200	864

步骤(6)根据POI与连通区域的索引关系，汇总各封闭多边形内部的POI，即将所属封闭多边形标记值相同的POI归为该封闭多边形的内部POI。

本实施实例中，以连通区域871为例，经过程序测试，共计有499个POI所属多边形标记值为871，即该道路分割的封闭多边形内有499个POI点。图4展示了剔除属于道路网络覆盖范围外连通区域POI点前后的POI分布情况。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (4)

1.一种根据所属道路网络封闭多边形统计POI的快速实现方法，包括如下步骤，其特征在于：

步骤(1)根据城市道路网络拓扑关系按比例输出道路网络图像，确定图像覆盖的实际平面坐标范围：(x_min,y_min)、(x_max,y_max)，单位为m；图像分辨率为m×n，单位为像素，权衡分析精度和速度，比例范围：(x_max-x_min):m∈[10,25]；

步骤(2)依次对图像进行顺时针旋转90°、灰度化、二值化处理，得到二值像素点矩阵P₁，其中取值为0的像素点表示为非道路，取值为1的像素点为道路，便于邻接关系识别，将矩阵扩展为(m+2)×(n+2)的新像素点矩阵P₂，扩展公式如下：

步骤(3)遍历矩阵P₂中位于1至m+1行且位于1至n+1列的像素点，对于标记为0的非道路像素点基于4邻接的邻接关系识别连通区域，并为属于同一个连通区域的像素点赋予唯一标记值，将经过连通区域标记的像素矩阵记为P₃；

步骤(4)遍历矩阵P₃第2行、第m+1行、第2列、第n+1列，统计上述像素点中标记为连通区域的像素点取值，记为outer＝[poly₁,poly₂,…,poly_a]，记录位于道路网络覆盖范围外的连通区域标记值，其中a为位于道路网络覆盖范围外的连通区域数量；

步骤(5)遍历POI＝[POI₁,POI₂,…POI_s]，将POI_i实际平面坐标(x_i,y_i)转化为像素坐标(r_i,c_i)，转化公式如下：

其中，ceil表示向上取整；若r_i＜1、r_i＞m+1、c_i＜1、c_i＞n+1中任意表达式成立，或P₃[r_i,c_i]∈outer，该POI位于道路网络覆盖范围外，将该POI从集合中剔除，否则为POI_i构建所属连通区域及道路网络分割的封闭多边形的索引POI_i:P₃[r_i,c_i]；

步骤(6)根据POI与连通区域的索引关系，汇总各封闭多边形内部的POI，即将所属封闭多边形标记值相同的POI归为该封闭多边形的内部POI。

2.根据权利要求1所述的一种根据所属道路网络封闭多边形统计POI的快速实现方法，其特征在于：所述步骤(1)中，为便于后续图像处理，根据实际道路网络拓扑关系输出的图像需要保证道路像素点与非道路像素点的RGB取值范围无交叉，并且横轴、纵轴按照相同的比例输出，即(x_max-x_min):m＝(y_max-y_min):n。

3.根据权利要求1所述的一种根据所属道路网络封闭多边形统计POI的快速实现方法，其特征在于：所述步骤(3)中，识别全部连通区域包含以下步骤：

步骤(3.1)令r＝2，c＝2，count＝2，P₃＝P₂；

步骤(3.2)若r＝m+2，结束步骤(3)，全部连通区域均以识别，输出连通区域标记结果P₃；否则，执行步骤(3.3)；

步骤(3.3)若P₃[r,c]≠0，执行步骤(3.6)；否则，令temp1＝[(r,c)]，P₃[r,c]＝count，统计尚未标记的且取值为0的像素点坐标，执行步骤(3.4)；

步骤(3.4)令temp2＝[]，统计temp1中的像素点邻接的标记为0的像素点坐标，邻接关系选择4邻接，具体做法为：遍历temp1，对任意(r_j,c_j)∈temp1，提取像素点(r_j,c_j)上、下、左、右四个像素点，即(r_j-1,c_j)、(r_j+1,c_j)、(r_j,c_j-1)、(r_j,c_j+1)，若P₃[r_j-1,c_j]＝0成立，将(r_j-1,c_j)追加写入temp2，并令P₃[r_j-1,c_j]＝count，对其余邻接像素点执行相同的判断和操作；

步骤(3.5)令temp1＝temp2，若

count＝count+1，执行步骤(3.6)；否则，执行步骤(3.4)；

步骤(3.6)令c＝c+1，若c＝n+2，r＝r+1；执行步骤(3.2)。

4.根据权利要求1所述的一种根据所属道路网络封闭多边形统计POI的快速实现方法，其特征在于：所述步骤(4)中，识别道路网络覆盖范围外的连通区域包含以下步骤：

步骤(4.1)令r＝2，c＝2，outer＝[]；

步骤(4.2)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.3)令c＝c+1，若c＝n+2，执行步骤(4.4)；否则，执行步骤(4.2)；

步骤(4.4)令r＝m+1，c＝2；

步骤(4.5)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.6)令c＝c+1，若c＝n+2，执行步骤(4.7)；否则，执行步骤(4.5)；

步骤(4.7)令r＝2，c＝2；

步骤(4.8)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.9)令r＝r+1，若r＝m+2，执行步骤(4.10)；否则，执行步骤(4.8)；

步骤(4.10)令r＝2，c＝n+1；

步骤(4.11)若P₃[r,c]≥2，将P₃[r,c]追加写入outer；

步骤(4.12)令r＝r+1，若r＝m+2，结束步骤(4)，道路网络覆盖范围外的连通区域均以识别，输出识别结果outer＝[poly₁,poly₂,…,poly_a]；否则，执行步骤(4.11)。

Images