CN101719154B - 一种基于栅格结构的空间索引建立方法和系统 - Google Patents

Abstract

本发明公开了一种基于栅格结构的空间索引建立方法和系统。所述方法，包括下列步骤：读取矢量图层，获取组成该图层所有矢量要素的要素类型和矢量要素的要素索引ID；根据所述矢量要素的要素类型和要素索引ID，栅格化所述矢量要素得到该矢量要素的栅格图像，栅格化的过程包括：扫描所述矢量要素的轮廓；根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述栅格图像，输出矢量图层与图像像素相对应的索引底图。

Description

一种基于栅格结构的空间索引建立方法和系统

技术领域

本发明涉及地理信息系统(Geographic Information System，GIS)领域，特别是涉及一种基于栅格结构的空间索引建立方法和系统。

背景技术

地理信息系统(Geographic Information System或Geo-Informationsystem，GIS)有时又称为“地学信息系统”或“资源与环境信息系统”。用于采集、存储、管理、处理、检索、分析和表达地理空间数据的计算机系统，是分析和处理海量地理数据的通用技术。从GIS系统应用角度，它是由计算机系统、地理数据和用户组成，通过对地理数据的集成、存储、检索、操作和分析，生成并输出各种地理信息，从而为土地利用、资源评价与管理、环境监测、交通运输、经济建设、城市规划以及政府部门行政管理提供新的知识，为工程设计和规划、管理决策服务。

随着计算机技术的飞速发展、空间技术的日新月异及计算机图形学理论的日渐完善，GIS技术也日趋成熟，“数字城市”、“电子地图”等概念的提出为人们的生活提供了更大的便捷。基于电子地图的特殊性，GIS的主要任务之一是有效地在后台服务器检索空间数据及快速响应前台不同用户的在线查询需求，而检索的速度将直接影响到整个系统执行的效率。提高检索速度的主要方法就是建立高效的空间索引。空间索引是对存储在介质上的某一对象集合中所有数据位置信息的描述，通过对空间对象集合建立索引，使系统能够快速的检索出这个对象集合中的某个特定对象。

目前已经有很多学者提出了很多的基于栅格结构的空间索引方法，如四叉树、R-树、R*-树、Hilbert R-树等。其中R-树是一种经典的目前被广泛应用的空间索引结构，它是1984年由A.Guttman提出，采用了最小约束矩形(MinimumBounding Rectangles，简称MBR)法来近似的表示空间对象。R-树是B-树在多维空间的扩展，因此它具有B-树的很多优点，如自动平衡、空间利用率高、适合于外存存储等。然而R-树的主要问题是，由于空间对象在建立过程中未进行分割，即中间节点目录矩形允许且可能重叠，因此在查询的过程中往往存在多条路径，并且其中的某些路径不包含查找结果，这就影响了查找性能。研究表明，随着索引空间维数的增加，R-树中间结点目录矩形的重叠迅速增加，这种结果是我们在实际处理GIS空间信息数据所不希望的。另外，可以预见，随着索引数据量的增加，树的高度及R-树中间结点目录矩形的重叠均会增加，不仅查找性能会下降，而且这种下降的速度会高于索引数据量的增加速度。

发明内容

本发明的目的在于提供一种基于栅格结构的空间索引建立方法和系统。利用建立的索引底图进行检索，保证了没有搜索重叠区域，提高了检索效率。

为实现本发明的目的而提供的一种基于栅格结构的空间索引建立方法，所述方法，包括下列步骤：

步骤100.读取矢量图层，获取组成该图层所有矢量要素的要素类型和矢量要素的要素索引ID；

步骤200.根据所述矢量要素的要素类型和要素索引ID，扫描所述矢量要素的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

步骤300.根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述轮廓，输出矢量图层与图像像素相对应的索引底图。

所述步骤100，包括下列步骤：

步骤110.配置OGR/GDAL库，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源；

步骤120.根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

步骤130.循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID。

所述步骤200，包括下列步骤：

步骤210.配置AGG库，定义绘制缓存，设置填充规则；

步骤220.根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓。

所述步骤300，包括：

步骤310.计算填充单元跨段存入绘制缓存；

步骤320.根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

步骤330.根据所述填充单元跨段和其对应的RGB颜色分量，填充所述图像像素，完成索引颜色存入绘制缓存，输出索引底图。

所述步骤220，包括下列步骤：

步骤221.利用坐标转换通道对所述矢量要素进行坐标转换；

步骤222.根据转换后的坐标，扫描矢量要素轮廓，计算轮廓像素的权值和覆盖面积。

所述步骤221采用的方法是将原整数栅格再细分为256×256个小单位，通过将双精度型数据左移八位与255做“与”操作，将小数部分转化到0-255之间的数。

所述步骤222是采用亚像素精度的Bresenham生成直线算法。

所述步骤310，包括下列步骤：

将所有生成轮廓的单元格进行排序，然后按照行扫描顺序从最小行到最大行逐行扫描；

利用图形对每个扫描栅格点的覆盖面积area进行判断，area等于0说明该图形未压盖此像素，无需填充；Area不等于0说明该图形部分压盖或全部压盖了此像素，需要扫描存进绘制控制器中；

对于未完全填充的栅格利用add_cell()函数进行标记，对于完全填充的栅格利用add_span()函数进行标记。

所述步骤320中解析过程采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

所述步骤330中，根据输入的采用24位RGB颜色参数，不再计算其混合后的颜色(*p++＝(((c.r-r)*alpha)+(r＜＜16))＞＞16)，而是直接将颜色赋给指针p(*p++＝c.r)，G分量和B分量也做如上修改，对于要素边界上的栅格，只要该要素对其压盖，就认为该栅格的像素颜色为该要素索引。

根据所述矢量要素与图像象素相对应的索引底图，对于点查询，只需根据鼠标点选范围的坐标对应的颜色进行图像像素到矢量要素的反解析，便得到了与所述图像像素对应的矢量要素；对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行图像像素到矢量要素的反解析，便得到了与查询框相交的矢量要素。

当分别取出R、G、B分量后再采用对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

为实现本发明的目的还提供一种基于栅格结构的空间索引建立系统，所述系统，包括：

矢量图层读取模块，用于读取矢量图层，获取组成该矢量图层所有矢量要素的要素类型和矢量要素的要素索引ID；

矢量要素处理模块，用于根据所述矢量要素的要素类型和要素索引ID，扫描所述矢量要素的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

索引底图输出模块，用于根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述轮廓，输出矢量图层与图像像素相对应的索引底图。

所述矢量图层读取模块，包括：

OGR/GDAL库，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源；

坐标转换模块，用于根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

读取模块，用于循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID。

所述矢量要素处理模块，包括：

绘制缓存配置模块，用于定义AGG绘制缓存，设置填充规则；

扫描模块，用于根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓。

所述索引底图输出模块，包括：

跨段计算模块，用于计算填充单元跨段存入绘制缓存；

颜色解析模块，用于根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

绘制模块，用于根据所述填充单元跨段和其对应的RGB颜色分量，填充所述图像像素，完成索引颜色存入绘制缓存，输出索引底图。

所述扫描模块，包括：

坐标计算模块，用于利用坐标转换通道对所述矢量要素进行坐标转换；

计算模块，用于根据转换后的坐标，扫描矢量要素轮廓，计算轮廓线所覆盖像素的权值和覆盖面积。

所述坐标计算模块，采用的方法是将原整数栅格再细分为256×256个小单位，通过将双精度型数据左移八位与255做“与”操作，将小数部分转化到0-255之间的数。

所述计算模块采用亚像素精度的Bresenham生成直线算法。

所述跨段计算模块，执行下列操作：

将所有生成轮廓的单元格进行排序，然后按照行扫描顺序从最小行到最大行逐行扫描；

利用图形对每个扫描栅格点的覆盖面积area进行判断，area等于0说明该图形未压盖此像素，无需填充；Area不等于0说明该图形部分压盖或全部压盖了此像素，需要扫描存进绘制控制器中；

对于未完全填充的栅格利用add_cell()函数进行标记，对于完全填充的栅格利用add_span()函数进行标记。

所述颜色解析模块，采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

所述绘制模块，根据输入的采用24位RGB颜色参数，不再计算其混合后的颜色(*p++＝(((c.r-r)*alpha)+(r＜＜16))＞＞16)，而是直接将颜色赋给指针p(*p++＝c.r)，G分量和B分量也做如上修改，对于要素边界上的栅格，只要该要素对其压盖，就认为该栅格的像素颜色为该要素索引。

所述系统，还包括：

搜索模块，用于根据所述矢量要素与图像象素相对应的索引底图，对于点查询，只需根据鼠标点选范围的坐标对应的颜色进行图像像素到矢量要素的反解析，便得到了与所述图像像素对应的矢量要素；对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行图像像素到矢量要素的反解析，便得到了与查询框相交的矢量要素。

所述搜索模块，分别取出R、G、B分量后再采用对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

本发明的有益效果是：

1.本发明的基于栅格结构的空间索引建立方法和系统，提出了一种用于GIS中的栅格空间索引机制，此方法利用AGG绘制引擎不仅能够快速建立栅格索引，而且与经典的空间索引结构R树比较，在建立索引的内存使用方面具有一定优势；

2.本发明的基于栅格结构的空间索引建立方法和系统，提出了两素对应索引方法在查询精度上具有一定优势，能够达到90％以上，误差大小取决于栅格大小。并且查询效率上大部分情况要优于R数，该发明尤其适用于GIS中海量密集数据的检索上，其检索的时间复杂度为一维线性；

3.由于本发明提出的基于AGG引擎的索引机制在检索效率上的优势，为网络地图服务模式下的空间检索的性能提升提供了很好的探索方向。

附图说明

图1是本发明一种基于栅格结构的空间索引建立方法的步骤流程图；

图2是本发明中读取矢量图层的方法的步骤流程图；

图3是本发明中对栅格图像进行扫描的方法步骤流程图；

图4是本发明中扫描所述栅格图像的轮廓的步骤流程图；

图5是本发明中通过坐标转换通道对所述矢量要素进行坐标转换的示意图；

图6是本发明中增量计算的示意图；

图7是本发明中计算权值(cover)和覆盖面积(area)的示意图；

图8是本发明中折线细节计算方法的示意图；

图9是本发明中扫描效果的示意图；

图10是本发明中输出索引底图的方步骤流程图；

图11是本发明中绘制控制器的计算过程的示意图；

图12是本发明中索引解析的示意图；

图13是本发明中对AGG填充部分作出修改后的示意图；

图14是本发明中矢量地图的原图；

图15是本发明中根据矢量地图生成的索引底图；

图16是本发明中采用框查询的示例图；

图17是本发明一种基于栅格结构的空间索引建立系统的结构示意图；

图18(a)是本发明中试验数据为密集数据的示意图；

图18(b)是本发明中试验数据为稀疏数据的示意图；

图19(a)是采用R树过滤边界的示例图；

图19(b)是采用本发明的两素对应过滤边界的示例图；

图20是建立索引比对图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明的一种基于栅格结构的空间索引建立方法和系统进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明的一种基于栅格结构的空间索引建立方法和系统，是利用AGG(Anti-Grain Geometry)引擎将矢量图层中的每个矢量要素转换成栅格数据并对该栅格数据进行渲染，而达到矢量要素与图像像素相对应，生成索引底图，利用该索引底图进行检索，保证了没有搜索重叠区域，提高了检索效率。

本发明要解决的技术问题包括：

1.矢量要素的栅格渲染；

2.要素索引ID的解析及与图像像素的关联；

3.利用这种两素(矢量要素与图像像素)对应索引机制实现的检索过程。

下面结合上述目标详细介绍本发明一种基于栅格结构的空间索引建立方法，图1是本发明一种基于栅格结构的空间索引建立方法的步骤流程图，如图1所示，所述方法包括下列步骤：

步骤100.读取矢量图层，获取组成该图层所有矢量要素的要素类型和要素索引ID；

本发明是采用OGR库来读取.shp格式的矢量图层，获取组成该图层所有要素的信息。该过程中的关键点如下：

(1)多边形带洞的情况。

在处理面数据的时候，一定要判断其是否存在内环，由于AGG绘制引擎渲染过程就是栅格填充过程，如果不判断内环情况，当出现面要素带洞情况将无法识别而将洞内全部填充。反之，如果考虑内环，可采取利用背景色绘制内环部分，则新绘制的内环中的像素将覆盖原来部分。这样，就能够处理地理信息数据中面要素带洞的情况。

(2)逻辑坐标转换。

在处理真实地理数据的时候，为了使真实数据在屏幕上完整的显示，最关键要做的就是世界坐标系向屏幕坐标系的转换，该过程需要地理外包、屏幕外包、及比例尺等参数。

图2是本发明中读取矢量图层的方法的步骤流程图，如图2所示，所述步骤100，包括下列步骤：

步骤110.配置OGR/GDAL库，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息；

从官方网站上下载动态链接库gdal.dll放在C:\WINDOWS\system32目录下，静态链接库gdal.lib放在开发环境的lib安装目录下，如：D:\ProgramFiles\Microsoft Visual Studio 8\SDK\v2.0\Lib，头文件放在开发环境的include目录下，如：D:\Program Files\Microsoft Visual Studio8\SDK\v2.0\include。

通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源。其中利用gdal库建立数据源，属于现有技术，在此不再赘述。

步骤120.根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

步骤130.循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID。

以下代码为生成索引底图的过程，即将矢量要素栅格化的过程，包括绘制内容，在这段代码里，当要绘制时会调用别的函数，因此它应该算是生成索引底图的整体过程。

//读取.shp文件ReadShp(){利用GDAL库获取数据源OGRDataSource *poDS
	//获取该数据源的图层poLayer及该图层的外包EnvpoLayer＝poDS-＞GetLayer(0)poLayer-＞GetExtent(&Env，true)根据外包计算世界坐标向屏幕坐标转换公式//循环读取该图层的要素OGRFeature*poFeatureWhile(poFeature＝poLayer-＞GetNextFeature()！＝NULL){定义OGR几何体OGRGeometry*geom获取要素索引信息FeatureId//点数据if(geom-＞getGeometryType()＝＝wkbPoint){//通过geom定义OGRPoint*poPoint＝(OGRPoint*)geom；坐标系转换利用AGG引擎绘制点数据}//线数据else if(geom-＞getGeometryType()＝＝wkbLineString){//定义线集合OGRLineString*poLine＝(OGRLineString*)geom循环取遍该线上所有点存入数组坐标系转换利用AGG引擎绘制线数据}//面数据else if(geom-＞getGeometryType()＝＝wkbPolygon){//定义面OGRPolygon
*poPolygon＝(OGRPolygon*)geom；获取该面得内环、外环坐标系转换利用AGG引擎填充外环内的面数据利用AGG引擎填充内环外的面数据}//Multi点else if(geom-＞getGeometryType()＝＝wkbMultiPoint){//获取OGRGeometryCollection*PointCollection＝(OGRGeometryCollection*)geom；按照处理点数据的方式循环处理该集合}//Multi线elseif(geom-＞getGeometryType()＝＝wkbMultiLineString){//获取OGRGeometryCollection*PointCollection＝(OGRGeometryCollection*)geom；按照处理线数据的方式循环处理该集合}//Multi面else if(geom-＞getGeometryType()＝＝wkbMultiPolygon){//获取OGRGeometryCollection*PointCollection＝(OGRGeometryCollection*)geom；按照处理面数据的方式循环处理该集合}}}

步骤200，根据所述矢量要素的要素类型和要素索引ID，扫描所述栅格图像的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

图3是本发明中对栅格图像进行扫描的方法步骤流程图，如图3所示，所述步骤200，包括下列步骤：

步骤210.配置AGG库，定义绘制缓存，设置填充规则；

(1)配置AGG库：下载AGGLite，将其生成静态链接库agg.lib放入开发环境的lib安装目录下，如：D:\Program Files\Microsoft Visual Studio8\SDK\v2.0\Lib中，工程项目中加入agglite.h头文件。

(2)定义AGG绘制缓存：

//定义画布大小

void*buf＝new unsigned void[width*height*3]

//根据该缓存创建绘制缓存

agg::rendering_buffer rbuf((unsigned char*)buf，width，height，width*3)；

//定义AGG绘制控制器

agg::renderer<agg::span_rgb24>ren(rbuf)；

agg::rasterizer ras；

//边沿虚化程度

ras.gamma(1.3)；

//填充规则

ras.filling_rule(agg::fill_even_odd)；

//画布背景色

ren.clear(agg::rgba8(255，255，255))；

设置好规则后，所有的绘制部分都利用rasterizer的对象ras来完成，由于它是利用缓存buf创建的，所以，所有绘制完的图像像素都存储在buf中，如需将图形输出到屏幕，只需将其利用bitmap与前缓存进行交换。如需保存成ppm格式，只需按照ppm文件格式将buf写入即可。

步骤220.根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓；

此时栅格图像还没有生成，这里的轮廓是指得到的矢量要素的轮廓)

在本步骤中，扫描所述矢量要素的轮廓，并非是实际绘制出栅格图像的轮廓，而仅仅是扫描矢量要素的轮廓，将扫描得到后的数据存入缓存中，待计算填充单元及跨段时使用。

图4是本发明中扫描所述矢量要素的轮廓的步骤流程图，如图4所示，所述步骤220，包括下列步骤：

步骤221.利用坐标转换通道对所述矢量要素进行坐标转换；

只要是利用AGG引擎绘制，都必须进行这种坐标转换，当获得了某个要素的所有坐标点准备扫描轮廓前，都要将这些坐标点经过这种坐标转换。

图5是本发明中通过坐标转换通道对所述矢量要素进行坐标转换的示意图，如图5所示，由于在渲染的过程中要考虑到亚像素精度的反走样以便于绘制控制器记录跨段信息，因此需要考虑双精度(double)型的真实数据中小数部分对整数部分所做的贡献，该贡献是指小数部分的大小对整数部分的影响程度。采用的方法是将原整数栅格再细分为256×256个小单位。通过将double型数据左移八位与255做“与”操作，即可将小数部分转化到0-255之间的数。

步骤222.根据转换后的坐标，扫描矢量要素轮廓，进而计算权值(cover)和覆盖面积(area)；

图6是本发明中增量计算的示意图，图7是本发明中计算权值(cover)和覆盖面积(area)的示意图。

扫描轮廓采用Y方向上逐行增加，X方向上逐单元增加的方法，将矢量要素的轮廓线所经过的每一个栅格点计算权值(cover)和覆盖面积(area)。以便用加权采样方法进行反走样的实现。这里运用了亚像素精度(subpixelaccuracy)的Bresenham生成直线算法，Bresenham算法是一种基于误差判别式来生成直线的方法。与传统Bresenham所不同的是，该算法利用误差判别选择像素的过程是基于亚像素的，即将单位栅格平均分成了256*256个子像素。根据直线斜率及增量依次判断该直线与所有横纵栅格线的交点，而此交点的计算不再仅仅依赖于斜率将原坐标(x_i，y_i)进行增量(delta)运算，而是需要转换成亚像素精度，增量(delta)的计算公式描述如下：

$\mathrm{delta}=(256-\mathrm{fy})\&times;\frac{\mathrm{\&Delta;x}}{\mathrm{\&Delta;y}}---\left(1\right)$

其中fy为实数坐标分量Y部分经过坐标转换通道所得的小数部分，Δx＝x₂-x₁，Δy＝y₂-y₁。如图6所示。X方向的增量采用与Y方向相同的方法。这样，在每一步上，便可求出该直线与X轴栅格线或Y轴栅格线的交点。根据这些点，可方便计算出受该直线所影响的整数栅格点的area和cover，计算公式描述如下：

cov er＝fy2-ft1                (2)

area＝(fx2+fx1)×cov er        (3)

如图7所示，对整数坐标A点影响的点有直线顶点a，及直线与X方向栅格线的交点b，整数坐标B点影响的点有直线与X轴栅格线交点b及与Y轴栅格线交点c。因此，A、B两点的权值和覆盖面积如下：

cov er_A＝fy_b-fy_a              (4)

area_A＝(fx_a+fx_b)×cov er_A    (5)

cov er_B＝fy_c-fy_b              (6)

area_B＝(fx_c+fx_b)×cov er_B    (7)

图8是本发明中折线细节计算方法的示意图，图9是本发明中扫描效果的示意图。

由于AGG引擎是根据封闭多边形在闭合区域内处理的，因此对于点、线、面的处理要求将它们都转化为封闭多边形。下面分别对这三种类型进行详细说明。

①点的扫描：由于点要素的顶点序列中只有一个点，若想形成封闭区域可以以该点为圆心画实心圆，半径可以自己指定。具体为：先将ras移动到圆上的任一点，ras.move_to_d(x-r，y)，然后根据角度利用正弦余弦计算圆周上所有点的坐标，最后利用循环将角度从0到360循环一次，ras.line_to_d(x+cos(a)*r，y+sin(a)*r)，ras便将圆周所有点画了一遍。

②线的扫描：线要素的扫描比较复杂，线要素实际是由一系列点序列组成的折线，单纯的这条折线在AGG中是绘制不出来的，必须给它赋予一定的宽度，形成封闭多边形。由于地理数据特殊性，两条直线的夹角可能会非常小，由于直线斜率无限接近，利用宽度线求交点将导致交点在无限远处，地图整体看来会有毛刺，为避免上述情形，不能根据直线宽度求交点，采用保留原直线宽度一半的做法，上部分求出四个点，下部分直接利用折线上的点，然后在分别对斜率不同、折线点个数不同的情况进行分类。最后绘制封闭的多边形。如图8所示，对于线段AB与线段BC所组成的折线来说，组成其的端点序列为：D-E-F-G-C-B-A。则D点的坐标描述如下：

$x_D=x_A-\mathrm{width}\&times;\frac{(y_B-y_A)}{\sqrt{{(x_B-x_A)}^2+{(y_B-y_A)}^2}}---\left(8\right)$

$y_D=y_A+\mathrm{width}\&times;\frac{(x_B-x_A)}{\sqrt{{(x_B-x_A)}^2+{(y_B-y_A)}^2}}---\left(9\right)$

用同样的方法依次计算出E、F、G。将D、E、F、G存入容器vector_up中。当所有端点都循环计算一次后，该直线上部分的所有点便都获得。最后将它们形成封闭多边形，如图9所示。

③面的扫描：根据面要素的顶点序列，只需先将ras移动到面要素的点序列的第一个点ras.move_to(vec[0].first，vec[0].second)，然后对所有点序列中每个点进行循环，使ras绘制边界ras.line_to_d(vec[i].first，vec[i].second)，最后轮廓扫描完毕。

AGG引擎提供的接口最重要的就是ras.move_to_d(x，y)和ras.line_to_d(x，y)。这两个函数是用来生成栅格图形的轮廓的。

步骤300.根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述栅格图像，得到所述矢量图层与图像像素相对应的索引底图；

本步骤中，所述栅格图像是由多个图像象素组成，进而得到与所述矢量要素对应的图像像素。

图10是本发明中输出索引底图的方步骤流程图，如图10所示，所述步骤300，包括下列步骤：

步骤310.根据所述权值和覆盖面积，计算填充单元跨段；

将所有生成轮廓的单元格(cell)进行排序，然后按照行扫描顺序从最小行到最大行逐行扫描，利用图形对每个扫描栅格点的覆盖面积area进行判断，area等于0说明该图形未压盖此像素，无需填充。Area不等于0说明该图形部分压盖或全部压盖了此像素，需要扫描存进绘制控制器中。对于未完全填充的栅格利用add_cell()函数进行标记，对于完全填充的栅格利用add_span()函数进行标记。由于cover有正负，所以在扫描每一行的时候都将所有cell的cover相加，(因为多边形一定是闭合的，所以当该跨段填充完以后cover会为0)，通过这样的判断，就能自动找到哪部分需要填充，哪部分不需要填充。从而解决了带洞多边形填充的问题。图11是本发明中绘制控制器的计算过程的示意图，如图11所示。在1-12标记的cell中执行add_cell，在图形内部执行add_span。

步骤320.根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

两素对应就是矢量要素与图形像素对应，具体指在绘制矢量要素的同时，生成一张与原图匹配的栅格底图，该底图就是其对应的索引图。因此该索引底上的每个像素点都对应了两部分内容，一是对应了该真实图像的颜色，二是对应了该矢量要素索引号ID转换为RGB的颜色。这种索引结构的生成过程和查询过程由于都是基于栅格的运算，因此效率上具有很大的优势。

本发明中，颜色格式采用RGB24，也就是说使用24位来表示一个像素，R、G、B每个分量都用8位表示，取值范围为0-255，内存中RGB各分量的排列顺序为：RGB。解析过程采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

图12是本发明中索引解析的示意图，如图12所示，如下三个面要素的索引ID分别为：159，24202，5513482。其解析的过程是先将159转换为二进制数，为10011111，其扩充到24位后将前16位补0，因此后八位为10011111，即159，为B分量；中间八位为00000000，即0，为G分量；前八位为00000000，即0，为R分量。因此该要素索引解析为RGB颜色为(0，0，159)。同理，索引号为24202的要素ID转换为颜色是(0，94，138)，索引号为5513482的要素ID解析为颜色是(84，33，10)。

这样便实现了矢量要素索引ID向颜色空间的转换过程，巧妙的将矢量要素索引和图像像素结合到一起，完成了两素对应方法。

上述解析过程的输入参数为索引ID，输出参数为R、G、B分量。伪代码如下：

void convert(int id，int&r，int&g，int&b){定义颜色数组如果要素ID不为0将该十进制数转换为二进制数存入数组中
	for(循环数组的低8位){计算B分量；}for(循环数组的中间8位){计算G分量；}for(循环数组的高8位){计算R分量；}}

由于AGG引擎及基于亚像素精度和反走样技术的，因此在要素栅格化后，填充颜色时会有混合色的计算过程，这将使颜色深浅过渡更加平缓。但是，该发明要做的工作利用解析为RGB颜色信息的索引ID对栅格化的要素进行填充，如果采用反走样，在要素的边沿，对于没有完全占满整个栅格的cell来说，会和其它对该栅格做贡献的颜色产生混合效果，颜色的变化将会导致索引ID的变化，因此在要素边沿处，就无法保证要素ID的正确性。

较佳地，为了避免这种情况，可以对AGG颜色填充部分做修改。根据输入的24位RGB颜色参数c，不再计算其混合后的颜色(*p++＝(((c.r-r)*alpha)+(r＜＜16))＞＞16)，而是直接将颜色赋给指针p(*p++＝c.r)，G分量和B分量也做如上修改。图13是本发明中对AGG填充部分作出修改后的示意图，修改后产生的效果如图13所示：

此处修改，默认了对于要素边界上的栅格，只要该要素对其压盖，就认为该栅格的像素颜色为该要素索引。因此，通过这种方式建立起的栅格索引的误差均产生在要素边缘的像素点上。误差的大小将取决于栅格的大小。

按照这种生成索引的方式，我们就可以对已经栅格化的要素进行填充。根据如上计算的ID颜色，将每个像素的RGB分量存入绘制缓存中。

所有的绘制部分都利用rasterizer的对象ras来完成，因此可利用接口ras.render(ren，agg::rgba8(red，green，blue))填充颜色，其中red、green、blue就是进行解析后的要素ID。

步骤330.根据所述填充单元跨段和其对应的RGB颜色分量，填充所述图像像素，完成索引颜色信息存入绘制缓存，输出索引底图；

图14是本发明中根据矢量地图生成的索引底图，如图14所示，生成的索引底图，如果绘制出来，颜色会是从黑色到蓝色逐渐变化的，因为黑色的RGB是(0，0，0)，蓝色影响最后8位，因此会逐渐变蓝，如果超过255，就会影响绿色分量G。以此类推，再影响红色分量R。

针对这种矢量要素与图像象素相对应的索引生成方法，搜索运算将异常简单。对于点查询，只需根据鼠标点选范围的坐标，获取(X，Y)点，根据已生成的绘制缓存进行反解析。对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行反解析即可。解析的代码如下：

//获取渲染好的画布

unsigned char*buffer＝(unsigned char*)this-＞bufId；

//获取R、G、B分量

unsigned int red＝buffer[(height-point_y)*width*3+point_x*3]；

unsigned int green＝buffer[(height-point_y)*width*3+point_x*3+1]；

unsigned int blue＝buffer[(height-point_y)*width*3+point_x*3+2]；

当分别取出R、G、B分量后再采用和生成颜色索引相反的方式，对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

图15是本发明中采用框查询的示例图，如图15所示，如果采用框查询，给定查询点A，及查询框大小Width×Heigth，则以起始点A逐行像素扫描至B，取出与框相交的像素颜色，有(0，0，159)、(0，94，138)、(84，33，10)，将他们从R分量开始反解析，分别为159，24202，5513482，由此便得到了与查询框相交的矢量要素。

相应于本发明的一种基于栅格结构的空间索引建立方法，还提供一种基于栅格结构的空间索引建立系统。

图16是本发明一种基于栅格结构的空间索引建立系统的结构示意图，如图16所示，所述系统，包括：

矢量图层读取模块1，用于读取矢量图层，获取组成该矢量图层所有矢量要素的要素类型和矢量要素的要素索引ID；

所述矢量图层读取模块1，包括：

OGR/GDAL库11，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源；

坐标转换模块12，用于根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

读取模块13，用于循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID。

矢量要素处理模块2，用于根据所述矢量要素的要素类型和要素索引ID，扫描所述矢量要素的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

所述矢量要素处理模块2，包括：

绘制缓存配置模块21，用于定义AGG绘制缓存，设置填充规则；

扫描模块22，用于根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓。

所述扫描模块22，包括：

坐标计算模块221，用于利用坐标转换通道对所述矢量要素进行坐标转换；

计算模块222，用于根据转换后的坐标，扫描矢量要素轮廓，计算轮廓线所覆盖像素的权值和覆盖面积。

所述坐标计算模块，采用的方法是将原整数栅格再细分为256×256个小单位，通过将双精度型数据左移八位与255做“与”操作，将小数部分转化到0-255之间的数。

所述计算模块采用亚像素精度的Bresenham生成直线算法。

索引底图输出模块3，用于根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述轮廓，输出矢量图层与图像像素相对应的索引底图。

所述索引底图输出模块3，包括：

跨段计算模块31，用于计算填充单元跨段存入绘制缓存；

颜色解析模块32，用于根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

绘制模块33，用于根据所述填充单元跨段和其对应的RGB颜色分量，填充所述轮廓，完成索引颜色存入绘制缓存，输出索引底图。

所述跨段计算模块，执行下列操作：

将所有生成轮廓的单元格进行排序，然后按照行扫描顺序从最小行到最大行逐行扫描；

利用图形对每个扫描栅格点的覆盖面积area进行判断，area等于0说明该图形未压盖此像素，无需填充；Area不等于0说明该图形部分压盖或全部压盖了此像素，需要扫描存进绘制控制器中；

对于未完全填充的栅格利用add_cell()函数进行标记，对于完全填充的栅格利用add_span()函数进行标记。

所述颜色解析模块，采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

所述绘制模块，根据输入的采用24位RGB颜色参数，不再计算其混合后的颜色(*p++＝(((c.r-r)*alpha)+(r＜＜16))＞＞16)，而是直接将颜色赋给指针p(*p++＝c.r)，G分量和B分量也做如上修改，对于要素边界上的栅格，只要该要素对其压盖，就认为该栅格的像素颜色为该要素索引。

所述基于栅格结构的空间索引建立系统，还包括：

搜索模块4，用于根据所述矢量要素与图像象素相对应的索引底图，对于点查询，只需根据鼠标点选范围的坐标对应的颜色进行图像像素到矢量要素的反解析，便得到了与所述图像像素对应的矢量要素；对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行图像像素到矢量要素的反解析，便得到了与查询框相交的矢量要素。

所述搜索模块，分别取出R、G、B分量后再采用对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

较佳地，由于两素对应生成索引的方法是采用将ID解析为颜色，然后按照栅格扫描的方法完成的。整个过程是以渲染为基础的，图形要素的渲染效率将对整个索引生成过程产生很大的影响，选择恰当的图形渲染引擎很重要，因此有必要对本发明所使用的AGG渲染引擎与GDI+渲染引擎在渲染效率上做对比。

GDI+是微软的.Net Framework的一个重要组成部分，增强了Windows的图形图像开发接口，同时它也改变了传统的GDI编程模式，使开发工作大幅度简化。我们选择渲染多边形来进行测验，根据多边形的边数不同，即多边形的复杂情况不同，和渲染多边形的个数不同来测试渲染时间进行比对，图17是本发明中AGG与GDI+的渲染效率对比图，实验结果如图17所示。

由此看出，AGG绘图引擎的渲染效率整体要高于GDI+，有时会超过10倍，并且渲染的图形要素越复杂，AGG所具有的优势越明显。

建立索引的最终目的将是为了高效的查询，查询的性能将直接影响到电子地图响应用户需求的速度。因此，有必要对该发明提出的索引机制在过滤精度及查询速度方面与传统的空间索引结构R树作对比。

图18(a)是本发明中试验数据为密集数据的示意图，图18(b)是本发明中试验数据为稀疏数据的示意图。

为了对比R-树和两素对应的索引行为，本发明设计了三个实验针对以下几个指标来测试它们的性能：两种索引的过滤精度、建立索引所花费的时间及内存、针对密集数据和稀疏数据利用两种索引分别进行点查询和框查询的效率。实验的测试环境：操作系统为WindowsXP(Service Pack2)，CPU为IntelPentium 1.73GHZ，内存为1G。

(1)过滤精度对比

测试R树索引与两素对应索引的过滤精度可采用点查询和框查询，但对于两素对应法来讲框查询也将转换为点查询，因此框过滤精度无异于点过滤的精度。本实验以点过滤为例，为了进行对比还将选一种基于矢量的计算几何的方法来作为参照。这种计算几何方法采用射线算法能够精确判断所查询的随机点的索引号，即是否属于某一要素。实验采用的数据为表1测试数据集中的D2。该实验通过生成随机点利用三种方法分别判断出能够查询到索引的点的数目，试验结果见表2。

表1测试数据集

编号	名称	几何类型	大小(MB)	要素数目	对应图示
						D1	BOUND_poly.shp	MultiPolygon	17.3	3407	图18(a)
D2	CS_99_03c.shp	MultiPolygon	1.48	123	图18(b)

表2过滤精确率对比

图19(a)是采用R树过滤边界的示例图，图19(b)是采用本发明的两素对应过滤边界的示例图。

从实验结果可以看出两素对应的过滤集能达到95％以上，要比R树的过滤集精度更高一些。产生这种差异的原因是，R树空间目标采用最小约束矩形(MBR)来加以近似表达，如图19(a)中红线所示，而且允许兄弟结点之间的相互重叠，因此对于精确匹配查询，R树不能保证唯一的检索路径，和唯一的检索结果。而基于AGG的两素索引机制查询只需进行简单的内存地址计算就可以检索出与查询条件匹配的索引项，相对比较精确，由于采用的是栅格方法，它的误差只产生在轮廓线的边沿像素，如图19(b)红线所示。

(2)建立索引性能对比

将R-树与两素对应分别对建立索引的时间与内存占用作对比。实验数据为表1测试数据集中的D1。图20是建立索引比对图，实验结果如图20所示。

从该结果可以看出，R树建立索引占用的内存会随着数据量的增大而增加，而两素对应由于是在内存中生成Bitmap，只要长度与宽度指定，它建立索引占用的内存不会随着数据量的增大而增加，由此看来，两素对应更适合大数据量的索引建立。从在建立索引占用的时间来看，R树所用时间要比两素对应要少一些，但在实际电子地图服务器的应用中，索引是预生成的。由于两素对应方法在查询方面的效率更高，所以更适用于在电子地图服务器应用。

(3)查询性能对比

比较R树索引与两素对应索引在查询方面的效率。为了试验结果的客观性，将分别对密集数据和稀疏数据做测试。密集数据采用的是表1测试数据集中的D1，该面图层的重叠区域较多。而稀疏数据采用的是表1测试数据集中的D2，该面图层的重叠部分较少。

点查询的实验对比如表3所示。实验数据为空间目标随机点，结果表明，无论是对密集数据还是稀疏数据，两素对应的点查询效率更高一些，因为这种索引机制是基于AGG的栅格渲染过程，每一点的索引保存在bitmap，点查询时，只需根据点坐标计算该点在bitmap中的行列地址，便可根据颜色信息解析出该点索引，时间复杂度是一维线性的O(n)，而对于R树，无论是密集数据还是稀疏数据都要从根结点向下搜索，当索引数据量增加时，这种索引效率就会降低。

框查询的实验结果如表3所示。实验数据为空间目标随机点，框的大小为200×200。结果表明，针对密集数据的框查询，两素对应的效率要明显高于R树索引效率，因为数据越密集，MBR的重叠越多，这样在框查询时，R树搜索的路径不唯一，导致效率下降。针对稀疏数据的框查询，两素对应的效率却没有R树高，原因是：由于数据稀疏，MBR的重叠少，R树的效率将有所提升，同时由于在框查询过程中，两素对应的查询方法是以随机点为定点，框的大小为扫描的长度、宽度，循环扫描单位栅格点，从而找出与查询框相交的要素索引。因此在数据极稀疏的情况下，循环时间将影响两素对应索引的整体效率。由于网络地图应用场景中，尤其是城市电子地图应用中，绝大部分情况要素都较为密集，因此此种情形不会影响两素对应算法在实际应用中的性能表现。

表3查询对比表

本发明的有益效果在于：

1.本发明的基于栅格结构的空间索引建立方法和系统，提出了一种用于GIS中的栅格空间索引机制，此方法利用AGG绘制引擎不仅能够快速建立栅格索引，而且与经典的空间索引结构R树比较，在建立索引的内存使用方面具有一定优势；

2.本发明的基于栅格结构的空间索引建立方法和系统，提出了两素对应索引方法在查询精度上具有一定优势，能够达到90％以上，误差大小取决于栅格大小。并且查询效率上大部分情况要优于R数，该发明尤其适用于GIS中海量密集数据的检索上，其检索的时间复杂度为一维线性；

3.由于本发明提出的基于AGG引擎的索引机制在检索效率上的优势，为网络地图服务模式下的空间检索的性能提升提供了很好的探索方向。

通过结合附图对本发明具体实施例的描述，本发明的其它方面及特征对本领域的技术人员而言是显而易见的。

以上对本发明的具体实施例进行了描述和说明，这些实施例应被认为其只是示例性的，并不用于对本发明进行限制，本发明应根据所附的权利要求进行解释。

Claims (12)

1.一种基于栅格结构的空间索引建立方法，其特征在于，利用所述方法建立的索引底图进行检索，没有搜索重叠区域，所述方法，包括下列步骤：

步骤100.读取矢量图层，获取组成该图层所有矢量要素的要素类型和矢量要素的要素索引ID；

步骤200.根据所述矢量要素的要素类型和要素索引ID，扫描所述矢量要素的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

步骤300.根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述轮廓，输出矢量要素与图像像素相对应的索引底图；

所述步骤100，包括下列步骤110、120、130：

步骤110.配置OGR/GDAL库，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源；

步骤120.根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

步骤130.循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID；

所述步骤200，包括下列步骤210、220：

步骤210.配置AGG库，定义绘制缓存，设置填充规则；

步骤220.根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓；

所述步骤300，包括下列步骤310、320、330：

步骤310.计算填充单元跨段存入绘制缓存；

步骤320.根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

步骤330.根据所述填充单元跨段和其对应的RGB颜色分量，填充所述图像像素，完成索引颜色存入绘制缓存，输出索引底图；

所述步骤220，包括下列步骤221、222：

步骤221.利用坐标转换通道对所述矢量要素进行坐标转换；

步骤222.根据转换后的坐标，扫描矢量要素轮廓，计算轮廓线所覆盖像素的权值和覆盖面积。

2.根据权利要求1所述的基于栅格结构的空间索引建立方法，其特征在于，所述步骤221采用的方法是将原整数栅格再细分为256×256个小单位，通过将双精度型数据左移八位与255做“与”操作，将小数部分转化到0-255之间的数。

3.根据权利要求1所述的基于栅格结构的空间索引建立方法，其特征在于，所述步骤222是采用亚像素精度的Bresenham生成直线算法。

4.根据权利要求1所述的基于栅格结构的空间索引建立方法，其特征在于，所述步骤320中解析过程采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

5.根据权利要求1所述的基于栅格结构的空间索引建立方法，其特征在于，根据所述矢量要素与图像象素相对应的索引底图，对于点查询，只需根据鼠标点选范围的坐标对应的颜色进行图像像素到矢量要素的反解析，便得到了与所述图像像素对应的矢量要素；对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行图像像素到矢量要素的反解析，便得到了与查询框相交的矢量要素。

6.根据权利要求5所述的基于栅格结构的空间索引建立方法，其特征在于，当分别取出R、G、B分量后再采用对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

7.一种基于栅格结构的空间索引建立系统，其特征在于，利用所述系统建立的索引底图进行检索，没有搜索重叠区域，所述系统，包括：

矢量图层读取模块，用于读取矢量图层，获取组成该矢量图层所有矢量要素的要素类型和矢量要素的要素索引ID；

矢量要素处理模块，用于根据所述矢量要素的要素类型和要素索引ID，扫描所述矢量要素的轮廓，为栅格化所述矢量要素得到该矢量要素的栅格图像做准备；

索引底图输出模块，用于根据所述矢量要素的要素索引ID解析出的RGB颜色分量，填充所述轮廓，输出矢量要素与图像像素相对应的索引底图；

所述矢量图层读取模块，包括：

建立数据源模块，用于配置OGR/GDAL库，通过GDAL库获取.shp格式的矢量图层，得到该矢量图层的图层信息和地理外包信息以建立数据源；

坐标转换模块，用于根据所述地理外包信息计算世界坐标系向屏幕坐标系的转换公式，以便使真实地理数据能够全部显示在有限的屏幕内；

读取模块，用于循环读取该矢量图层中的所有矢量要素得到所述矢量要素的索引ID；

所述矢量要素处理模块，包括：

绘制缓存配置模块，用于配置AGG库，定义AGG绘制缓存，设置填充规则；

扫描模块，用于根据定义的AGG绘制缓存，扫描所述矢量要素的轮廓；

所述索引底图输出模块，包括：

跨段计算模块，用于计算填充单元跨段存入绘制缓存；

颜色解析模块，用于根据所述矢量要素的要素索引ID，解析出RGB颜色分量；

绘制模块，用于根据所述填充单元跨段和其对应的RGB颜色分量，填充所述图像像素，完成索引颜色存入绘制缓存，输出索引底图；

所述扫描模块，包括：

坐标计算模块，用于利用坐标转换通道对所述矢量要素进行坐标转换；

计算模块，用于根据转换后的坐标，扫描矢量要素轮廓，计算轮廓线所覆盖像素的权值和覆盖面积。

8.根据权利要求7所述的基于栅格结构的空间索引建立系统，其特征在于，所述坐标计算模块，采用的方法是将原整数栅格再细分为256×256个小单位，通过将双精度型数据左移八位与255做“与”操作，将小数部分转化到0-255之间的数。

9.根据权利要求7所述的基于栅格结构的空间索引建立系统，其特征在于，所述计算模块采用亚像素精度的Bresenham生成直线算法。

10.根据权利要求7所述的基于栅格结构的空间索引建立系统，其特征在于，所述颜色解析模块采用的做法是将十进制ID转换为二进制，先从B分量的低8位开始存，共能存256个ID，超过部分存入G分量，再超过部分存入R分量。

11.根据权利要求7所述的基于栅格结构的空间索引建立系统，其特征在于，所述系统，还包括：

搜索模块，用于根据所述矢量要素与图像象素相对应的索引底图，对于点查询，只需根据鼠标点选范围的坐标对应的颜色进行图像像素到矢量要素的反解析，便得到了与所述图像像素对应的矢量要素；对于框查询，只需根据框定点坐标及框的大小确定在缓存中的位置，利用行扫描循环取出相交栅格的颜色进行图像像素到矢量要素的反解析，便得到了与查询框相交的矢量要素。

12.根据权利要求11所述的基于栅格结构的空间索引建立系统，其特征在于，所述搜索模块分别取出R、G、B分量后再采用对颜色信息先从R分量的高8位解析，再解析G分量的中间8位，最后解析B分量的低8位。

Images