CN102542035A

CN102542035A - 基于扫描线法的多边形栅格化并行转换方法

Info

Publication number: CN102542035A
Application number: CN2011104423514A
Authority: CN
Inventors: 陈振杰; 张帅; 李飞雪; 王亚飞; 李满春; 蒲英霞; 王加胜; 程亮
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2011-12-20
Filing date: 2011-12-27
Publication date: 2012-07-04
Anticipated expiration: 2031-12-27
Also published as: CN102542035B

Abstract

本发明公开了基于扫描线法的多边形栅格化并行转换方法，属于地理信息系统领域。其包括输入命令行参数；MPI并行初始化，获取总的进程数和当前进程数，采用对等式并行模式，各进程分别解析命令行参数，分别收集引导符后的参数值，利用OGROpen方法读取矢量数据源，判断是否为0号进程；采用数据并行策略，划分栅格数据集合矢量多边形，然后分发各个进程，每个进程同时进行多边形的栅格化；写栅格数据，各进程分别更新栅格分块，并输出转换后的栅格数据。利用本发明进行大数据量的多边形栅格化的操作，可以得到较高的效率和满意的转换结果，充分提高了高性能服务器的多核/多处理器对多边形栅格化的转换处理速度，极大地缩小了多边形栅格化的转换时间。

Description

基于扫描线法的多边形栅格化并行转换方法

技术领域

本发明涉及一种矢量数据的栅格化方法，特别是涉及基于扫描线法的多边形栅格化并行转换方法。

背景技术

地理信息系统( GIS )是以空间数据为基础，获取、表达、处理、管理、分析和显示空间数据并为地理研究和地理决策服务的计算机服务系统。空间数据通常有矢量数据(Vector Data)和栅格数据(Raster Data)两种形式。矢量数据是通过记录坐标的方式，表示点、线、多边形等地理实体，自然地理实体的位置是用其在坐标参考系中的空间位置来定义的，坐标空间设为连续，其特点是定位明显，属性隐含。而栅格数据又称为网格数据(grid cell)，即将平面划分为m×n个像元（正方形小方格），每个像元由行列号唯一地确定其所在平面位置，给像元赋予属性以表达该覆盖的自然地理实体的类型，其最明显的特点是属性明显，定位隐含。

在GIS空间分析时，由于栅格形式的GIS 数据非常适合诸如空间叠加、空间相关和空间模拟等空间分析，因而通常需要把矢量数据转化成栅格数据。矢量数据栅格化被广泛认为是地理信息系统中的基础问题。矢量数据栅格化包括点的栅格化、线的栅格化以及多边形的栅格化。点和线的栅格化方法目前已经比较成熟，方法也趋于固定。多边形的栅格化就是对矢量数据的面状图斑根据给定的栅格化像元的大小离散化为像元的集合，像元值为矢量面状图斑所具有的某种属性值。长期以来以矢量多边形栅格化的研究最为热点。多边形的栅格化已有很多算法，传统的串行算法如内部点扩散法、复数积分算法、射线法、扫描法和边界代数法等，这些方法各有优缺点，目前还没有一种标准统一的最优算法。随着计算机的快速发展，又产生了许多新方法，比如：2004年，王建等在《地理与地理信息科学》20卷第3期中发表“矢量数据向栅格数据转换的一种改进算法”一文，总结和分析了多边形栅格化的传统算法与新方法，提出了一种改进的折线边界跟踪方法，保证了多边形填充的精度；2005年，章孝灿等在《计算机辅助设计与图形学学报》17卷第6期中发表“面状矢量拓扑数据快速栅格化算法”一文，提出了一种快速栅格化算法—差分边界标志与累加扫描算法，2009年，武广臣等在《测绘科学》43卷第1期中发表“矢量数据栅格化的一种有效方法—环绕数法”一文，提出了一种基于计算几何转角理论的环绕数法，着重处理了自相交多边形的栅格化问题；2010年，李青元等在《武汉大学学报信息科学版》35卷第8期中发表“基于绘制—检出的矢量数据栅格化方法研究”一文，探讨了基于绘制—检出的矢量数据栅格化方法。然而研究的重点都是围绕改进串行算法展开的，对于海量多边形的栅格化效率的提升相当有限。

随着对地观测技术的长足发展，海量栅格数据需求迅速激增，数据量为T级的栅格数据普遍存在（1TB=1024GB）。海量矢量数据栅格化呈现出计算高度密集的特点，耗时巨大。现有的矢量数据栅格化串行算法模式和传统的硬件平台，已经无法满足海量地理数据处理的需求。基于并行计算集群与多核处理器的新型硬件架构的逐渐普及，为受制于计算性能而难以展开的地理数据转换提供了契机。本发明充分利用现有的高性能计算机和并行处理技术，基于数据并行策略采用对等式的并行程序设计模式，提出了一种基于矢量多边形扫描线的数据并行方法，有效地解决了海量的矢量数据栅格化的问题。

发明内容

1.发明要解决的技术问题

针对如上所述，从数据需求方面说，矢量数据向栅格数据的转换是GIS一直研究的基础问题；从软硬件上来说，逐渐普及的并行计算集群与多核处理器的新型硬件架构需要得到有效利用；最重要的从效率上来说，海量的矢量数据的栅格化运行时间过长、效率过低的问题，本发明提供了基于扫描线法的多边形栅格化并行转换方法，该方法采用数据并行策略，即将待处理的矢量多边形按进程数进行划分，然后分发给各个进程，每个进程同时进行多边形的栅格化。这样数据划分策略是基于目标栅格数据的逻辑划分，可以有效完成大数据量的矢量多边形的栅格化，且不必考虑生成的栅格数据的拼接问题，取得了良好的栅格化效果和较高的效率，满足海量的矢量数据的栅格化要求。

2.技术方案

发明原理：一般来说，数据并行的实现过程是主进程将待处理的数据分派到其他若干子进程分别处理，再由主进程负责收集不同子进程的数据处理结果并进行组合，达到多处理器共同完成某一个任务的目的。本发明中利用栅格数据行列规整的特征，先由一个进程生成一个像素值初始化为0的栅格数据集，接着按给定的进程数划分该栅格数据，得到与进程数相等的栅格数据分块。然后查询栅格分块范围内的多边形（包括与该栅格分块相交的多边形）并提取出来分发给各个进程，每个进程进行相同的多边形栅格化的操作。

1. 基于扫描线法的多边形栅格化并行转换方法，包括以下步骤：

步骤1：输入命令行参数：mpirun -np 8 hpgc_rasterize -a GHDLDM -l hpgc_data -of HFA -tr 20 20 ~\data\hpgc_data.shp ~\data\test_result.img；

步骤2：

(1) MPI并行初始化，获取总的进程数和当前进程数，并注册GDAL/OGR格式驱动；

(2) 采用对等式并行模式，各进程分别解析命令行参数，分别收集引导符后的参数值；

(3) 以GDAL为矢量数据读写工具，利用OGROpen方法读取矢量数据源；

(4) 判断是否为0号进程，若是0号进程，进行以下操作：首先，判断栅格文件是否存在，若存在，以GDAL中的Update方式打开，若不存在就获取HFA格式驱动，创建目标栅格数据集；目标数据集是根据格式驱动类型、输出的栅格文件名、像元的长和宽、波段数、数据类型参数利用GDAL中GDALCreate方法创建；

步骤3：采用数据并行策略，划分栅格数据集和矢量多边形：

（1）划分生成的栅格数据集，根据对等式并行模式按总的进程数划分生成的栅格数据集，每个进程分别负责处理栅格分块范围内的多边形填充，即通过计算栅格数据的总行数RasterYSize，根据总的进程数np划分生成的栅格数据，每个进程处理的栅格分块行数为nYChunkSize =ceil[RasterYSize/np]， ceil[n]表示不小于n的最小整数；每个栅格块的起始行坐标为iY=cp*nYChunkSize，对第i个进程的起始行坐标iY+nYChunkSize进行判断，iY+nYChunkSize > RasterYSize时，第i+1号进程处理的栅格分块的行数为nYChunkSize = RasterYSize – iY；

（2）按总的进程数划分矢量多边形，即基于划分的栅格数据分块的范围对所有矢量多边形按进程数进行划分，实现多边形划分的具体操作为：各进程分别进行空间查询，根据各栅格分块的左上角、右上角、左下角、右下角四个角点组成的矩形为查询范围，获取与该矩形区域相交的多边形，包括位于该矩形区域内的多边形和与该矩形区域的边界相交的多边形，同时获取用于填充像元的图层的属性字段GHDLDM；

步骤4：读取栅格分块数据，即指定栅格分块数据的数据类型，计算每行栅格数据所需要的存储量，继而计算出每个栅格分块数据的总数据量，以GDAL为影像数据读写工具，利用GDALDataset.RasterIO方法读取生成的栅格数据；

步骤5：

(1) 定义数组获取所有多边形的顶点坐标并按单个多边形的顶点顺序一一存储到数组里；

（2）将多边形的顶点坐标一一进行坐标转换，从地理坐标转换为行列坐标；

（3）根据各栅格分块的矢量多边形最小外接矩形确定扫描线的起止行坐标，扫描线是像素中心点所在的直线；

（4）然后逐行扫描，循环处理每个多边形的每一条边与扫描线的交点，然后逐行根据列坐标的大小对交点进行排序并两两组合，填充同一行中两交点间的栅格单元，依次扫描，直至所有扫描线扫描完毕；

步骤6：以GDAL为影像数据读写工具，利用GDALDataset.RasterIO方法写栅格数据，各进程分别更新栅格分块，并输出转换后的栅格数据；

步骤7：在IBM System x3500-M3X系列服务器Openmpi并行环境下，编译并使用实验数据测试。

步骤2中采用对等式并行模式，即首先指定一个进程创建一个初始像元值为0的栅格数据集，然后采用数据并行策略，划分栅格数据集和矢量多边形，最后每个进程分别完成各栅格分块范围内的多边形的填充，各进程间相对独立，避免了数据并行时的通信，同时也避免了栅格化后结果的拼接。

步骤3中数据并行策略，采用基于划分栅格数据集进而对矢量数据进行划分的数据划分方式，即采用空间求交查询过滤的方式查询各栅格分块里的矢量多边形，完成对矢量多边形的划分。

所述步骤5中为了保证栅格化的精度，多边形的所有顶点均以双精度字符类型进行存储，判读边缘像元采用中心点的方法，若该像素的中心点位于多边形内，该像素属于该多边形。

3.有益效果

相对于现有的矢量栅格化方法，本发明将扫描线栅格化方法与并行处理技术结合，基于数据并行策略采用对等式的并行程序设计模式，实现了矢量多边形栅格化的并行处理，在一定程度上解决了海量矢量数据的栅格化问题。具体有益效果如下：

第一，改进了传统扫描线法的矢量栅格化方法。本发明考虑已有算法的并行性，针对栅格数据行列规整的特点，采用行扫描线法，对多边形的每条边界与扫描线进行相交判断，计算交点坐标。本方法先计算矢量多边形的最小包容窗，缩小了扫描范围；通过比较多边形顶点坐标与当前扫描线的坐标的大小，能够快速的完成对所有多边形边界的遍历。

第二，改善了栅格化的精度控制。本发明以像素中心点所在的射线作为扫描线；对边界像素的归属问题，采用像素中心点的判断方法，即当像素中心点落在多边形内，则认为该像素在多边形内。对于矢量多边形的外边界（不与任何多边形相邻的边），本发明另外提供一个参数选项设置，选用该参数即表示与多边形边界相交的像元单元都被认为处在该多边形内，这种处理方法在很多空间分析中是需要的。

第三，提升多边形栅格化的效率，缩短了多边形栅格化的时间，在海量数据处理方面尤其明显。本发明充分利用现有的多核处理器和计算机集群，并行处理多边形的栅格化问题。利用一个进程创建目标栅格数据，然后采用数据并行策略，在逻辑层面上对生成的栅格数据和矢量数据进行划分，然后多进程分别同时进行处理。该并行策略是在逻辑上进行的，避免了进程间的通信，同时有效地避免了栅格分块的接边处理。

综上，本发明充分利用多核处理器等新型计算机设备，能够快速地实现海量多边形的栅格化。实践证明，该方法具有较高的并行性和良好的并行效率，达到了海量多边形快速栅格化的目的。本发明可以应用于矢量多边形向栅格数据的大批量的转换。

附图说明

图1实施例1中实验数据土地利用矢量多边形图；

图2基于扫描线的多边形栅格化并行转换流程图；

图3基于数据并行的栅格数据的数据划分图；

图4每条行扫描线扫描填充多边形的流程图；

图5多边形栅格化结果图；

图6串行方法与并行方法运行结果的对比图；

图7栅格化并行加速比结果图。

具体实施方式

实施例1

本实施例指定源数据为中国长沙市土地现状调查地类图斑数据，区域范围为西经111.877度，东经114.256度，南纬27.836度，北纬28.666度；区域总面积11819.46平方公里。数据格式为ESRI shapefile格式，图斑总数为692177，数据量为938MB.。数据空间参考系为1980西安坐标系。局部图如图1所示。本实施例具体实施按照图2所示的技术路线，采用如图3所示的数据并行策略，采用标准C++编程语言在Microsoft Visual Studio 2008开发平台下开发，并在MPI并行环境下实现，矢量和栅格数据的读写操作通过开源地理数据格式转换类库GDAL实现。

程序运行选择IBM System x3500-M3X系列服务器环境。服务器硬件配置为CPU 2颗，规格为Intel Xeon Quad Core E5620(主频2.4GHz，12MB Cache，四核）；内存为8GB（2根4GB内存条，规格为DDR3 1333MHz LP RDIMM）；硬盘为2TB（4个500GB硬盘，规格为7.2K 6Gbps NL SAS 2.5-inch SFF Slim-HS HD），网络为集成的双口千兆以太网。软件配置：操作系统为Centos Linux 6.0，其中MPI的实现产品选择Openmpi 1.4.1。GDAL版本为1.8.1。

具体实施步骤如下：

步骤1：输入命令行参数：mpirun -np 8 hpgc_rasterize -a GHDLDM -l hpgc_data -of HFA -tr 20 20 ~\data\hpgc_data.shp ~\data\test_result.img

其中，‘-a’ ‘-l’ ‘-of’ ‘-tr’ 均为命令行的引导符，‘mpirun’表示调用MPI应用程序的声明参数；‘-np 8’表示调用8个进程数；‘ hpgc_rasterize’为本发明编译后生成的exe运行程序文件名；‘-a GHDLDM’用于像元填充的多边形图层的属性字段，属性字段名为GHDLDM；‘-l hpgc_data’用于待栅格化的矢量图层，矢量图层名称为hpgc_data；‘ -of HFA’表示转换后的栅格格式，HFA是常用栅格格式Erdas Imagine .img的格式代码；‘-tr 20 20’表示生成的栅格单元长和宽分别为20m*20m，~\data\hpgc_data.shp为源数据，~\data\test_result.img为栅格化后的输出文件。

步骤2：(1) MPI并行初始化，获取总的进程数和当前进程数，并注册GDAL/OGR格式驱动；(2) 采用对等式并行模式，即各进程相对独立，分别解析命令行参数，分别收集引导符后的参数值，采用对等式并行模式，首先指定一个进程创建一个初始像元值为0的栅格数据集，然后采用数据并行策略，划分栅格数据集和矢量多边形，最后每个进程分别完成各栅格分块范围内的多边形的填充。各进程间相对独立，避免了数据并行时的通信，同时也避免了栅格化后结果的拼接；(3) 以GDAL为矢量数据读写工具，利用OGROpen方法读取矢量数据源；(4) 判断是否为0号进程，若是0号进程，进行以下操作：首先，判断栅格文件是否存在，若存在，以GDAL中的Update方式打开，若不存在就获取HFA格式驱动，创建目标栅格数据集；目标数据集是根据格式驱动类型、输出的栅格文件名、像元的长和宽、波段数、数据类型等参数利用GDAL中GDALCreate方法创建。利用进程号为0的进程生成栅格数据的核心代码为：

GDALDatasetH hDstDS = NULL; //定义目标栅格数据集

int nLayerCount = CSLCount(papszLayers); //定义图层个数

OGREnvelope sEnvelop; //定义栅格数据外围的矩形窗

double adfProjection[6]; //定义仿射变换参数数组

if(cp==0) //判断是否为0号进程

{

hDstDS = GDALOpen(pszDstFilename, GA_Update);//以GA_Update方式打开栅格数据

if (hDstDS == NULL) //若栅格数据不存在

{

hDriver = GDALGetDriverByName(pszFormat); //获取格式驱动

std::vector<OGRLayerH> ahLayers; //定义矢量图层对象

for (i = 0; i < nLayerCount; i++)

{

//获取数据源文件中的图层

OGRLayerH hLayer = OGR_DS_GetLayerByName(hSrcDS, papszLayers[i]);

if (hLayer == NULL) //若图层不存在

continue;

ahLayers.push_back(hLayer); //将矢量图层入栈

}

//设置栅格数据的范围

sEnvelop.MinX = sLayerEnvelop.MinX;

sEnvelop.MinY = sLayerEnvelop.MinY;

sEnvelop.MaxX = sLayerEnvelop.MaxX;

sEnvelop.MaxY = sLayerEnvelop.MaxY;

//设置仿射变换参数

adfProjection[0] = sEnvelop.MinX;

adfProjection[1] = dfXRes;

adfProjection[2] = 0;

adfProjection[3] = sEnvelop.MaxY;

adfProjection[4] = 0;

adfProjection[5] = -dfYRes;

//利用GDALCreate方法生成栅格数据集

hDstDS = GDALCreate(hDriver, pszDstFilename, nXSize, nYSize, nBandCount,

eOutputType, papszCreateOptions);

GDALSetGeoTransform(hDstDS, adfProjection); //对栅格数据进行投影变换

}

MPI_Barrier(MPI_COMM_WORLD); //其他进程均等待

步骤3：采用数据并行策略，划分栅格数据集和矢量多边形，所述数据并行策略，采用基于划分栅格数据集进而对矢量数据进行划分的数据划分方式，即采用空间求交查询过滤的方式查询各栅格分块里的矢量多边形，完成对矢量多边形的划分，这种划分在逻辑上的，因此不需要考虑跨分块多边形的特殊情况的考虑。划分生成的栅格数据集，根据对等式并行模式按总的进程数划分生成的栅格数据集，每个进程分别负责处理栅格分块范围内的多边形填充，即先计算栅格数据集的总行数RasterYSize=4510，总进程数为np=12，栅格分块的行数为Ychunksize，取不小于[RasterYSize/np]的最小整数376，即进程号0-10进程处理的栅格分块行数为376，第11号进程处理的栅格分块行数为374。按总的进程数划分矢量多边形，即基于划分的栅格数据分块的范围对所有矢量多边形进行划分，实现多边形划分的具体操作为：各进程分别进行空间查询，根据各栅格分块的左上角、右上角、左下角、右下角四个角点组成的矩形为查询范围，获取与该矩形区域相交的多边形（包括位于该矩形区域内的多边形和与该矩形区域的边界相交的多边形），同时获取用于填充像元的图层的属性字段GHDLDM。划分栅格数据集和矢量多边形的核心代码为：

int nYChunkSize, iY; //定义栅格分块的行数和每个进程处理的栅格分块起始行号

int pulx, puly, plrx, plry; //定义栅格分块矩形的左上角和右下角的行列坐标

double dulx,duly,dlrx,dlry; //仿射变换后栅格分块矩形的左上角和右下角的地理坐标

double gT[6]; //定义仿射变换参数数组

if(cp==0) //0进程号输出栅格数据的多边形个数，列数和行数

{

printf("\nThe total feature size of this layer is %d.\n"

"The X size of the destination raster is %d.\n"

"The Y size of the destination raster is %d.\n",

OGR_L_GetFeatureCount(hSrcLayer,true),

((GDALDataset*)hDstDS)->GetRasterXSize(),

((GDALDataset*)hDstDS)->GetRasterYSize());

}

//获取仿射变换参数

GDALGetGeoTransform(hDstDS,gT);

//计算栅格分块的行数

nYChunkSize = ceil(((GDALDataset*)hDstDS)->GetRasterYSize()/(double)np);

//计算每个进程处理的栅格分块的起始行数

iY=cp*nYChunkSize;

//计算最后一个进程处理的栅格分块的行数

if( nYChunkSize + iY >( (GDALDataset*)hDstDS)->GetRasterYSize() )

nYChunkSize = ((GDALDataset*)hDstDS)->GetRasterYSize() - iY;

//计算栅格分块矩形的左上角和右下角的行列坐标

pulx=0;

puly=iY;

plrx=((GDALDataset*)hDstDS)->GetRasterXSize()-1;

plry=iY+nYChunkSize-1;

//通过仿射变换公式计算栅格分块矩形的左上角和右下角的地理坐标

dulx=gT[0] + gT[1] * pulx + gT[2]*puly;

duly=gT[3] + gT[4] * pulx + gT[5]*puly;

dlrx=gT[0] + gT[1] * plrx + gT[2]*plry;

dlry=gT[3] + gT[4] * plrx + gT[5]*plry;

OGRGeometry *poSpatialFilter=NULL;

//由栅格分块的矩形的四个角点坐标定义该环

OGRLinearRing oRing;

oRing.addPoint( dulx, duly);

oRing.addPoint( dulx, dlry);

oRing.addPoint( dlrx, dlry);

oRing.addPoint( dlrx, duly);

oRing.addPoint( dulx, duly);

//创建复合多边形的几何类型

poSpatialFilter = OGRGeometryFactory::createGeometry(wkbPolygon);

((OGRPolygon *) poSpatialFilter)->addRing( &oRing );

//利用空间过滤器进行空间查询

OGR_L_SetSpatialFilter (hSrcLayer,(OGRGeometryH )poSpatialFilter);

步骤4：指定栅格分块数据的数据类型，计算每行栅格数据所需要的存储量，继而计算出每个栅格分块数据的总数据量。以GDAL为影像数据读写工具，利用GDALDataset.RasterIO方法读取生成的栅格数据。读取栅格分块数据的核心代码为：

int nYChunkSize, iY; //定义栅格分块行数和每个分块的起始行坐标

unsigned char *pabyChunkBuf; //定义栅格分块的数据量

int nScanlineBytes; //定义每行栅格所占的数据量

//指定栅格数据的数据类型

if( poBand->GetRasterDataType() == GDT_Byte )

eType = GDT_Byte;

else

eType = GDT_Float32;

//计算每行栅格数据所占的数据量

nScanlineBytes = nBandCount * poDS->GetRasterXSize()* (GDALGetDataTypeSize(eType)/8);

//每个进程处理的栅格数据块

nYChunkSize = ceil(poDS->GetRasterYSize()/(double)np);

//每个进程处理的栅格块的起始行

iY=cp*nYChunkSize;

//当前每个进程处理的栅格分块的行数

int nThisYChunkSize;

int iShape;

nThisYChunkSize = nYChunkSize;

if( nThisYChunkSize + iY > poDS->GetRasterYSize() )

nThisYChunkSize = poDS->GetRasterYSize() - iY;

pabyChunkBuf = (unsigned char *) VSIMalloc(nThisYChunkSize * nScanlineBytes);

//初始化错误类型

CPLErr eErr = CE_None;

//判断读取栅格数据是否成功

eErr =poDS->RasterIO(GF_Read,

0, iY, poDS->GetRasterXSize(), nThisYChunkSize,

pabyChunkBuf,poDS->GetRasterXSize(),nThisYChunkSize,

eType, nBandCount, panBandList, 0, 0, 0 );

步骤5：

(2) 将多边形的顶点坐标一一进行坐标转换，从地理坐标转换为行列坐标；

(3) 根据各栅格分块的矢量多边形接矩形确定扫描线的起止行坐标；以像素中心点所在的直线作为扫描线。

(4) 然后逐行扫描，循环处理每个多边形的每一条边与扫描线的交点。然后逐行根据列坐标的大小对交点进行排序并两两组合，填充同一行中两交点间的栅格单元。为了保证栅格化的精度，多边形的所有顶点均以双精度字符类型进行存储，判读边缘像元采用中心点的方法，若该像素的中心点位于多边形内，则将该像素属于该多边形。依次扫描，直至所有扫描线扫描完毕。每条扫描线进行扫描填充多边形的流程图如图4所示。提取多边形顶点坐标的核心代码为：

//获取几何体的类型

OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());

int i;

//处理几何体类型为线环的情况

if ( EQUAL(poShape->getGeometryName(),"LINEARRING") )

{

OGRLinearRing *poRing = (OGRLinearRing *) poShape; //指定几何体类型为线环

int nCount = poRing->getNumPoints(); //获取环的顶点数

int nNewCount = aPointX.size() + nCount; //存储顶点坐标的数组的当前序号

//确定存储X、Y坐标数组的最小长度

aPointX.reserve( nNewCount );

aPointY.reserve( nNewCount );

//将线环上的所有节点依次存储到X、Y

for ( i = nCount - 1; i >= 0; i-- )

{

aPointX.push_back( poRing->getX(i) );

aPointY.push_back( poRing->getY(i) );

}

//将多边形的顶点个数入栈

aPartSize.push_back( nCount );

}

//若几何体类型为单个多边形

else if( eFlatType == wkbPolygon )

{

OGRPolygon *poPolygon = (OGRPolygon *) poShape;

//依次收集多边形外边界的顶点坐标

GDALCollectRingsFromGeometry( poPolygon->getExteriorRing(),

aPointX, aPointY, aPointVariant,

aPartSize, eBurnValueSrc );

//依次处理多边形内环的顶点坐标

for( i = 0; i < poPolygon->getNumInteriorRings(); i++ )

GDALCollectRingsFromGeometry( poPolygon->getInteriorRing(i),

aPointX, aPointY, aPointVariant,

aPartSize, eBurnValueSrc );

}

//若几何体类型为复合多边形

else if(eFlatType == wkbGeometryCollection )

{

OGRGeometryCollection *poGC = (OGRGeometryCollection *) poShape;

//依次获取每个多边形的顶点坐标

for( i = 0; i < poGC->getNumGeometries(); i++ )

GDALCollectRingsFromGeometry( poGC->getGeometryRef(i),

aPointX, aPointY, aPointVariant,

aPartSize, eBurnValueSrc );

}

扫描线法填充多边形的核心代码为：

int i; //顶点个数

int y; //扫描线的行坐标

int miny, maxy,minx,maxx; //最小矩形范围

double dminy, dmaxy; //行坐标最小值，最大值

//相邻的两坐标点

double dx1, dy1;

double dx2, dy2;

double dy; //扫描线的行坐标

double intersect; //扫描线与多边形边界交点的列坐标

int ind1, ind2; //多边形两顶点

int ints, part; //交点的个数，多边形的序号

int *polyInts, polyAllocated;

int horizontal_x1, horizontal_x2; //位于同一行的两个点的列坐标

//计算所有顶点的个数

int n = 0;

for( part = 0; part < nPartCount; part++ )

n += panPartSize[part];

polyInts = (int *) malloc(sizeof(int) * n); //分配n个整型空间指针数组

polyAllocated = n;

//初始化行坐标的最小值和最大值

dminy = padfY[0];

dmaxy = padfY[0];

for (i=1; (i < n); i++) {

if (padfY[i] < dminy) {

dminy = padfY[i];

}

if (padfY[i] > dmaxy) {

dmaxy = padfY[i];

}

//确定最小外接矩形的范围

miny = (int) dminy;

maxy = (int) dmaxy;

if( miny < 0 )

miny = 0;

if( maxy >= nRasterYSize )

maxy = nRasterYSize-1;

minx = 0;

maxx = nRasterXSize - 1;

//从最小行开始，逐行扫描，循环处理至最大行

for (y=miny; y <= maxy; y++) {

int partoffset = 0;

dy = y +0.5; //扫描线取栅格中心点所在的直线

part = 0; //多边形的序号

ints = 0;

memset(polyInts, -1, sizeof(int) * n); //初始化指针数组

//循环处理每个多边形的边界

for (i=0; (i < n); i++) {

//一个多边形的边界遍历结束，取下一个多边形顶点

if( i == partoffset + panPartSize[part] ) {

partoffset += panPartSize[part];

part++;

}

//顶点遍历到多边形最后一点

if( i == partoffset ) {

ind1 = partoffset + panPartSize[part] - 1;

ind2 = partoffset;

}

//取多边形边界的相邻两点

else {

ind1 = i-1;

ind2 = i;

}

//确定相邻两点的行坐标

dy1 = padfY[ind1];

dy2 = padfY[ind2];

//当两点所在直线不与扫描线相交，结束本次循环，取下一顶点

if( (dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy) )

continue;

//所取的两顶点按行坐标大小排序

if (dy1 < dy2) {

dx1 = padfX[ind1];

dx2 = padfX[ind2];

} else if (dy1 > dy2) {

dy2 = padfY[ind1];

dy1 = padfY[ind2];

dx2 = padfX[ind1];

dx1 = padfX[ind2];

} else

//当两顶点位于同一行

if (padfX[ind1] > padfX[ind2])

{

horizontal_x1 = (int) floor(padfX[ind2]+0.5);

horizontal_x2 = (int) floor(padfX[ind1]+0.5);

if ( (horizontal_x1 > maxx) || (horizontal_x2 <= minx) )

continue;

//填充该行两点间的所有像元

pfnScanlineFunc( pCBData, y, horizontal_x1, horizontal_x2 - 1, (dfVariant == NULL)?0:dfVariant[0] );

continue;

}

else

continue;

}

//计算扫描线与两顶点间线段的交点

if(( dy < dy2 ) && (dy >= dy1))

{

intersect = (dy-dy1) * (dx2-dx1) / (dy2-dy1) + dx1;

polyInts[ints++] = (int) floor(intersect+0.5);

}

//将交点坐标按列坐标大小排序

qsort(polyInts, ints, sizeof(int), llCompareInt);

//取偶数点组成的线段，填充两点间的像元

for (i=0; (i < (ints)); i+=2)

{

if( polyInts[i] <= maxx && polyInts[i+1] > minx )

{

pfnScanlineFunc( pCBData, y, polyInts[i], polyInts[i+1] - 1, (dfVariant == NULL)?0:dfVariant[0] );

}

步骤6：以GDAL为影像数据读写工具，利用GDALDataset.RasterIO方法写栅格数据。各进程分别更新栅格分块。并输出转换后的栅格数据。生成的栅格数据局部图如图5所示。

步骤7：在IBM System x3500-M3X系列服务器Openmpi并行环境下，编译并使用实验数据测试，依次测试总进程数为1—24个进程的并行加速比。加速比为最优的串行算法的运行时间与稳定的并行算法运行时间之比。每种进程测试10次，取10次运行时间的平均值作为该进程执行并行转换方法的运行时间。当进程数等于1时，执行多边形栅格化的串行转换方法，此时由一个进程负责整个数据的栅格化任务。当进程数大于1时，执行多边形栅格化的并行转换方法。最终取得的串行时间与并行时间对比图如图6所示，图7表示加速比。从图6和图7可以发现，在串行环境下，处理相同的实验数据运行时间为42.755s，而并行转换方法采用12个进程时运行的时间最少，为9.508s，即该并行转换方法取得的最大加速比为4.5，并行的效率较高。通过本方法充分提高了多核/多处理器的高性能服务器对多边形栅格化的转换处理速度，极大地缩小了多边形栅格化的转换时间。

Claims

1.基于扫描线法的多边形栅格化并行转换方法，包括以下步骤：

步骤2：

步骤3：采用数据并行策略，划分栅格数据集和矢量多边形：

步骤5：

2.根据权利要求1所述的基于扫描线法的多边形栅格化并行转换方法，其特征在于：所述步骤2中采用对等式并行模式，即首先指定一个进程创建一个初始像元值为0的栅格数据集，然后采用数据并行策略，划分栅格数据集和矢量多边形，最后每个进程分别完成各栅格分块范围内的多边形的填充，各进程间相对独立，避免了数据并行时的通信，同时也避免了栅格化后结果的拼接。

3.根据权利要求2所述的基于扫描线法的多边形栅格化并行转换方法，其特征在于：所述步骤3中数据并行策略，采用基于划分栅格数据集进而对矢量数据进行划分的数据划分方式，即采用空间求交查询过滤的方式查询各栅格分块里的矢量多边形，完成对矢量多边形的划分。

4.根据权利要求1或3所述的基于扫描线法的多边形栅格化并行转换方法，其特征在于：所述步骤5中为了保证栅格化的精度，多边形的所有顶点均以双精度字符类型进行存储，判读边缘像元采用中心点的方法，若该像素的中心点位于多边形内，该像素属于该多边形。