GIS缓冲区算法对比研究（Buffer Algorithm）

来源：互联网发布：为什么会害怕别人知乎编辑：程序博客网时间：2024/05/02 05:02

----by wangsh

一．介绍算法：

缓冲区分析是邻近度分析的一种，缓冲区是为了识别某一地理实体或空间物体对其周围地物的影响度而在其周围建立具有一定宽度的带状区域。缓冲区作为独立的数据层进行叠加分析，可应用到道路、河流、环境污染源、居民点、辐射源等的空间分析，为某种应用目的提供科学依据，另外，结合不同的专业模型，可以在生活、军事、城乡规划等领域发挥重要的作用。

针对点、线、面不同的几何类型，建立缓冲区的方式相互有所不同。建立点缓冲区比较简单，即以某点要素为圆心，以缓冲半径 R作圆，得到点要素的缓冲区；线要素的缓冲区是以线为轴，以 R为距离作两侧的平行线，在线的两端构建两个半圆弧段，和平行线一起组成缓冲区；面缓冲区的建立，是以面要素的边界为基线向内外侧作平行线，平行线和基线里的区域就是面缓冲区。除此之外，还可以对栅格数据建立缓冲区，根据不同的模型方程建立动态缓冲区，不论对什么样的数据建立缓冲区，其基本方法都是相似的。(见参考文献4）

开源版本的的GRASS，GEOS和JTS里面都有Buffer算法模块，其中GEOS是JTS的C++版本。

使用GEOS的应用有PostGIS (C API),MapServer (C API),Quantum GIS (C API),OGR (C API),GRASS (C API),Shapely (C API),INGRES (C API),SpatiaLite (C API),MapGuide Open Source GeoDjango (C API),MapWindow GIS (C API), osm2pgsql (C++ API),osgEarth (C++ API),MonetDB (C API),rgeos (C API)，其中还有软件包的FME和Autodesk MapGuide Enterprise。

GRASS版本介绍（见http://grass.itc.it/gdp/html_grass63/v.buffer.html）（相关的工程介绍见http://download.osgeo.org/grass/grass6_progman/）

这里仅列出相关源码，待以后列出对比算法分析，算法性能等。

二．算法说明及伪代码演示

缓冲区实现算法有矢量方法和栅格方法两种。其中矢量方法数据量小，方法相对成熟，栅格图像需要进行栅格像元之间进行布尔运算，当缓冲区较大时会带来较重的运算负荷，实际运用中存在一定的局限性。矢量方法算法一般遵循以下步骤：

点：确定中心点——以中心点为圆心、 R为半径生成一个圆——得到缓冲区边界；

线、面：确定轴线——以距离 R生成中心轴线的平行曲线——处理转角弧段——对生成的弧段进行求交、合并运算——生成缓冲区边界；

常用的矢量数据中心线扩张算法：

角分线法

基本思想：即“简单平行线法”，在轴线的两边作出平行线，在转角处形成尖角，两端形成弧段，组成缓冲区。

缺陷：难以保证在尖角处缓冲区左右边线等宽；校正过程复杂，主要体现在轴线折角很大和很小时的情况；算法模型复杂，主要是因为几何生成过程中需要处理较多的异常。

凸角圆弧法

基本思想：顾名思义，即是在转角外侧用圆弧来代替尖角，内侧仍然使用尖角的方法，生成缓冲区。

实施步骤：

1. 直线性判断，判断相邻三点是否在同一直线上；

2. 折点凸凹性判断，确定转角的地方哪侧使用直线求交，哪侧使用圆弧连接；

3. 凸点圆弧的嵌入，即将转角外侧形成的圆弧和两边的线段相连；

4. 边线关系的判别与处理，岛屿多边形参与缓冲区边界的构成，重叠多边形不参与缓冲区边界的构成；

5. 缓冲区边界的形成，具体是将重叠区域进行合并，绘制外围的边线，包括岛屿多边形的轮廓，形成最终的缓冲区边界。

在缓冲区算法中，需要注意的一个问题是缓冲区多边形的重叠与合并，包括同一要素缓冲区的重叠和多个要素之间缓冲区的重叠。栅格数据缓冲区内的栅格具有一个与其影响度对应的一个值，如果重叠区域具有相同影响度则任取一值，如果不同则采取影响度大的代替影响度小的方法处理。对于矢量数据的处理算法有三种：数学运算法；矢量－栅格转换法；矢量－栅格混合法。(见参考文献4）

这里简单给出GEOS的源码(源码位置在GEOS/source/operation/buffer/BufferOp.cpp)：

doubleBufferOp::precisionScaleFactor(constGeometry *g,

doubledistance,

intmaxPrecisionDigits)

{

constEnvelope *env=g->getEnvelopeInternal();

doubleenvSize=std::max(env->getHeight(),env->getWidth());

doubleexpandByDistance=distance > 0.0 ?distance : 0.0;

doublebufEnvSize=envSize + 2 *expandByDistance;

// the smallest power of 10 greater than the buffer envelope

intbufEnvLog10=(int) (std::log(bufEnvSize) /std::log(10.0) + 1.0);

intminUnitLog10=bufEnvLog10 -maxPrecisionDigits;

// scale factor is inverse of min Unit size, so flip sign of exponent

doublescaleFactor=std::pow(10.0,-minUnitLog10);

returnscaleFactor;

}

/*public static*/

Geometry*BufferOp::bufferOp(constGeometry *g,doubledistance,

intquadrantSegments,

intnEndCapStyle)

{

BufferOpbufOp(g);

bufOp.setQuadrantSegments(quadrantSegments);

bufOp.setEndCapStyle(nEndCapStyle);

returnbufOp.getResultGeometry(distance);

}

/*public*/

Geometry*BufferOp::getResultGeometry(doublenDistance)

{

distance=nDistance;

computeGeometry();

returnresultGeometry;

}

/*public*/

Geometry*BufferOp::getResultGeometry(doublenDistance,int nQuadrantSegments)

{

distance=nDistance;

setQuadrantSegments(nQuadrantSegments);

computeGeometry();

returnresultGeometry;

}

/*private*/

VoidBufferOp::computeGeometry()

{

#ifGEOS_DEBUG

std::cerr<<"BufferOp::computeGeometry: trying with original precision"<<std::endl;

#endif

//bufferReducedPrecision(); return; // FIXME: remove this code

bufferOriginalPrecision();

if (resultGeometry!=NULL)return;

std::cerr <<"bufferOriginalPrecision failed (" << saveException.what() <<"), trying with reduced precision"

<< std::endl;

constPrecisionModel& argPM = *(argGeom->getFactory()->getPrecisionModel());

if (argPM.getType() ==PrecisionModel::FIXED )

bufferFixedPrecision(argPM);

else

bufferReducedPrecision();

}

/*private*/

VoidBufferOp::bufferReducedPrecision()

{

// try and compute with decreasing precision

for (intprecDigits=MAX_PRECISION_DIGITS;precDigits >= 0;precDigits--)

{

#ifGEOS_DEBUG

std::cerr<<"BufferOp::computeGeometry: trying with precDigits "<<precDigits<<std::endl;

#endif

try {

bufferReducedPrecision(precDigits);

}catch (constutil::TopologyException&ex) {

saveException=ex;

// don't propagate the exception - it will be detected by fact that resultGeometry is null

}

if (resultGeometry!=NULL) {

// debug

//if ( saveException ) std::cerr<<saveException->toString()<<std::endl;

return;

}

// tried everything - have to bail

throwsaveException;

}

/*private*/

VoidBufferOp::bufferOriginalPrecision()

{

BufferBuilderbufBuilder;

bufBuilder.setQuadrantSegments(quadrantSegments);

bufBuilder.setEndCapStyle(endCapStyle);

//std::cerr<<"computing with original precision"<<std::endl;

try

{

resultGeometry=bufBuilder.buffer(argGeom,distance);

}

catch (constutil::TopologyException&ex)

{

// don't propagate the exception - it will be detected by

// fact that resultGeometry is null

saveException=ex;

//std::cerr<<ex->toString()<<std::endl;

}

//std::cerr<<"done"<<std::endl;

}

VoidBufferOp::bufferReducedPrecision(intprecisionDigits)

{

doublesizeBasedScaleFactor=precisionScaleFactor(argGeom,distance,precisionDigits);

std::cerr <<"recomputing with precision scale factor = "

<<sizeBasedScaleFactor

<<std::endl;

assert(sizeBasedScaleFactor>0);

PrecisionModelfixedPM(sizeBasedScaleFactor);

bufferFixedPrecision(fixedPM);

}

/*private*/

VoidBufferOp::bufferFixedPrecision(constPrecisionModel&fixedPM)

{

PrecisionModelpm(1.0); // fixed as well

// FIXME: both MCIndexSnapRounder and SimpleSnapRounder

// makes buffer_snapround.xml test fail.

#if 0

snapround::MCIndexSnapRounder inoder(pm); // fail

snapround::SimpleSnapRounder inoder(pm); // fail

#else

algorithm::LineIntersectorli(&fixedPM);

IntersectionAdderia(li);

MCIndexNoderinoder(&ia);// This works fine (but does not snapround)

#endif

ScaledNodernoder(inoder,fixedPM.getScale());

BufferBuilderbufBuilder;

bufBuilder.setWorkingPrecisionModel(&fixedPM);

bufBuilder.setNoder(&noder);

bufBuilder.setQuadrantSegments(quadrantSegments);

bufBuilder.setEndCapStyle(endCapStyle);

// this may throw an exception, if robustness errors are encountered

resultGeometry=bufBuilder.buffer(argGeom,distance);

}

这里给出GRASS两个版本（源码位置grass/vector/v.buffer和grass-. src/vector/v.buffer，具体算法在grass/lib/vector/Vlib/buffer.c和grass/lib/vector/Vlib/buffer2.c）

第一个版本就是buffer.c版本（根据点，线和面类型分开处理：点和线归为一类处理；面分解为线对象处理，最后合并为面对象）：

/*!

/brief Create buffer around the line line.

This function is replaced by Vect_line_buffer().

Buffer is closed counter clockwise polygon. Warning: output line

may contain loops!

/param InPoints input line

/param distance create buffer in distance

/param tolerance maximum distance between theoretical arc and polygon segments

/param[out] OutPoints output line

VoidVect_line_buffer(conststructline_pnts *InPoints,double distance,

doubletolerance, struct line_pnts *OutPoints)

{

doubledangle;

intside, npoints;

staticstructline_pnts *Points =NULL;

staticstructline_pnts *PPoints =NULL;

distance =fabs(distance);

dangle = 2 *acos(1 - tolerance / fabs(distance)); /* angle step */

if (Points ==NULL)

Points =Vect_new_line_struct();

if (PPoints ==NULL)

PPoints =Vect_new_line_struct();

/* Copy and prune input */

Vect_reset_line(Points);

Vect_append_points(Points,InPoints,GV_FORWARD);

Vect_line_prune(Points);

Vect_reset_line(OutPoints);

npoints =Points->n_points;

if (npoints <= 0) {

return;

}

elseif (npoints == 1) { /* make a circle */

doubleangle, x, y;

for (angle = 0;angle < 2 * PI; angle += dangle) {

x =Points->x[0] +distance *cos(angle);

y =Points->y[0] +distance *sin(angle);

Vect_append_point(OutPoints,x,y, 0);

}

/* Close polygon */

Vect_append_point(OutPoints,OutPoints->x[0],OutPoints->y[0], 0);

}

else { /* 2 and more points */

for (side = 0;side < 2; side++) {

double angle, sangle;

double lx1, ly1, lx2, ly2;

double x, y, nx, ny, sx, sy, ex, ey;

/* Parallel on one side */

if (side == 0) {

Vect_line_parallel(Points,distance,tolerance, 0, PPoints);

Vect_append_points(OutPoints,PPoints,GV_FORWARD);

}

else {

Vect_line_parallel(Points, -distance,tolerance, 0,PPoints);

Vect_append_points(OutPoints,PPoints,GV_BACKWARD);

}

/* Arc at the end */

/* 2 points at theend of original line */

if (side == 0) {

lx1 =Points->x[npoints - 2];

ly1 =Points->y[npoints - 2];

lx2 =Points->x[npoints - 1];

ly2 =Points->y[npoints - 1];

}

else {

lx1 =Points->x[1];

ly1 =Points->y[1];

lx2 =Points->x[0];

ly2 =Points->y[0];

}

/* normalized vector */

vect(lx1,ly1,lx2, ly2, &nx, &ny);

/* starting point */

sangle =atan2(-nx,ny); /* starting angle */

sx =lx2 +ny * distance;

sy =ly2 -nx * distance;

/* end point */

ex =lx2 -ny * distance;

ey =ly2 +nx * distance;

Vect_append_point(OutPoints,sx,sy, 0);

/* arc */

for (angle =dangle;angle < PI; angle += dangle) {

x =lx2 + distance * cos(sangle +angle);

y =ly2 + distance * sin(sangle +angle);

Vect_append_point(OutPoints,x,y, 0);

}

Vect_append_point(OutPoints,ex,ey, 0);

}

/* Close polygon */

Vect_append_point(OutPoints,OutPoints->x[0],OutPoints->y[0], 0);

}

Vect_line_prune(OutPoints);

return;

}

第二个版本就是buffer2.c版本（这个是Google Summer of Code 2008的源码，将点，线和面的buffer分开处理的）：

/* area_outer and area_isles[i] must be closed non self-intersecting lines

side: 0 - auto, 1 - right, -1 left

staticvoidbuffer_lines(structline_pnts *area_outer,structline_pnts **area_isles,

intisles_count, int side, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints,struct line_pnts ***iPoints,int *inner_count)

{

structplanar_graph *pg2;

structline_pnts *sPoints, *cPoints;

structline_pnts **arrPoints;

inti, count = 0;

intres, winding;

intauto_side;

intmore = 8;

intallocated = 0;

doublepx, py;

G_debug(3,"buffer_lines()");

auto_side = (side == 0);

/* initializations */

sPoints =Vect_new_line_struct();

cPoints =Vect_new_line_struct();

arrPoints =NULL;

/* outer contour */

G_debug(3," processing outer contour");

*oPoints =Vect_new_line_struct();

if (auto_side)

side =

get_polygon_orientation(area_outer->x,area_outer->y,

area_outer->n_points -

1) ? LEFT_SIDE : RIGHT_SIDE;

convolution_line(area_outer,da,db, dalpha, side, round, caps, tol,

sPoints);

pg2 =pg_create(sPoints);

extract_outer_contour(pg2, 0, *oPoints);

res =extract_inner_contour(pg2, &winding,cPoints);

while (res != 0) {

if (winding == 0) {

if (!Vect_point_in_poly(cPoints->x[0],cPoints->y[0],area_outer)) {

if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)

G_fatal_error(_("Vect_get_point_in_poly() failed"));

if (!point_in_buf(area_outer,px,py, da, db, dalpha)) {

add_line_to_array(cPoints, &arrPoints, &count, &allocated,

more);

cPoints =Vect_new_line_struct();

}

res =extract_inner_contour(pg2, &winding,cPoints);

}

pg_destroy_struct(pg2);

/* inner contours */

G_debug(3," processing inner contours");

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

if (auto_side)

side =

get_polygon_orientation(area_isles[i]->x,area_isles[i]->y,

area_isles[i]->n_points -

1) ?RIGHT_SIDE : LEFT_SIDE;

convolution_line(area_isles[i],da,db, dalpha, side, round, caps,

tol,sPoints);

pg2 =pg_create(sPoints);

extract_outer_contour(pg2, 0,cPoints);

res =extract_inner_contour(pg2, &winding,cPoints);

while (res != 0) {

if (winding == -1) {

/* we need to check if the area is in the buffer.

I've simplfied convolution_line(), so that it runs faster,

however that leads to ocasional problems */

if (Vect_point_in_poly

(cPoints->x[0],cPoints->y[0],area_isles[i])) {

if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)

G_fatal_error(_("Vect_get_point_in_poly() failed"));

if (!point_in_buf(area_isles[i],px,py, da, db, dalpha)) {

add_line_to_array(cPoints, &arrPoints, &count,

&allocated,more);

cPoints =Vect_new_line_struct();

}

res =extract_inner_contour(pg2, &winding,cPoints);

}

pg_destroy_struct(pg2);

}

arrPoints =G_realloc(arrPoints,count * sizeof(structline_pnts *));

*inner_count =count;

*iPoints =arrPoints;

Vect_destroy_line_struct(sPoints);

Vect_destroy_line_struct(cPoints);

G_debug(3,"buffer_lines() ... done");

return;

}

/*!

/brief Creates buffer around line.