Skip to content

Commit

Permalink
Merge pull request #1309 from rouault/fix_1307
Browse files Browse the repository at this point in the history 
Rasterize with MERGE_ALG=ADD: avoid burning several times intermediate points of linestrings (fixes #1307)
  • Loading branch information
@rouault
rouault authored Feb 22, 2019
2 parents 2cdcdca + 932fc9a commit c5565eb
Show file tree
Hide file tree
Showing 4 changed files with 156 additions and 6 deletions.
55 changes: 55 additions & 0 deletions autotest/alg/rasterize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -329,4 +329,59 @@ def test_rasterize_6():
gdal.RasterizeLayer(mask_ds, [1], layer, burn_values=[1], options=["ALL_TOUCHED"])


###############################################################################
# Test rasterizing linestring with multiple segments and MERGE_ALG=ADD
# Tests https://github.com/OSGeo/gdal/issues/1307

def test_rasterize_merge_alg_add_multiple_segment_linestring():

# Setup working spatial reference
sr_wkt = 'LOCAL_CS["arbitrary"]'
sr = osr.SpatialReference(sr_wkt)

data_source = ogr.GetDriverByName('MEMORY').CreateDataSource('')
layer = data_source.CreateLayer('', sr, geom_type=ogr.wkbLineString)
feature = ogr.Feature(layer.GetLayerDefn())
# Diagonal segments
feature.SetGeometryDirectly(ogr.CreateGeometryFromWkt('LINESTRING(0.5 0.5,100.5 50.5,199.5 99.5)'))
layer.CreateFeature(feature)
feature = ogr.Feature(layer.GetLayerDefn())
# Vertical and horizontal segments
feature.SetGeometryDirectly(ogr.CreateGeometryFromWkt('LINESTRING(30.5 40.5,30.5 70.5,50.5 70.5)'))
layer.CreateFeature(feature)

ds = gdal.GetDriverByName('Mem').Create('', 10, 10, 1, gdal.GDT_Byte)
ds.SetGeoTransform([0, 20, 0, 100, 0, -10])
ds.SetProjection(sr_wkt)

ds.GetRasterBand(1).Fill(0)
gdal.RasterizeLayer(ds, [1], layer, burn_values=[1], options=["MERGE_ALG=ADD"])

got = struct.unpack('B' * 100, ds.ReadRaster())
expected = (0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 1, 1, 0, 0, 0, 0, 1, 0, 0,
0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 1, 0, 0, 0, 1, 0, 0, 0, 0,
0, 1, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
assert got == expected, '%s' % str(got)

ds.GetRasterBand(1).Fill(0)
gdal.RasterizeLayer(ds, [1], layer, burn_values=[1], options=["MERGE_ALG=ADD", "ALL_TOUCHED"])

got = struct.unpack('B' * 100, ds.ReadRaster())
expected = (0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
0, 1, 1, 0, 0, 0, 1, 1, 1, 0,
0, 1, 0, 0, 0, 1, 1, 0, 0, 0,
0, 1, 0, 0, 1, 1, 0, 0, 0, 0,
0, 1, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
assert got == expected, '%s' % str(got)
3 changes: 2 additions & 1 deletion gdal/alg/gdal_alg_priv.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,8 @@ void GDALdllImageLineAllTouched( int nRasterXSize, int nRasterYSize,
int nPartCount, int *panPartSize,
double *padfX, double *padfY,
double *padfVariant,
llPointFunc pfnPointFunc, void *pCBData );
llPointFunc pfnPointFunc, void *pCBData,
int bAvoidBurningSamePoints );

void GDALdllImageFilledPolygon( int nRasterXSize, int nRasterYSize,
int nPartCount, int *panPartSize,
Expand Down
9 changes: 6 additions & 3 deletions gdal/alg/gdalrasterize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -462,7 +462,8 @@ gv_rasterize_one_shape( unsigned char *pabyChunkBuf, int nXOff, int nYOff,
&(aPointX[0]), &(aPointY[0]),
(eBurnValueSrc == GBV_UserBurnValue)?
nullptr : &(aPointVariant[0]),
gvBurnPoint, &sInfo );
gvBurnPoint, &sInfo,
eMergeAlg == GRMA_Add );
else
GDALdllImageLine( sInfo.nXSize, nYSize,
static_cast<int>(aPartSize.size()),
Expand Down Expand Up @@ -496,7 +497,8 @@ gv_rasterize_one_shape( unsigned char *pabyChunkBuf, int nXOff, int nYOff,
static_cast<int>(aPartSize.size()), &(aPartSize[0]),
&(aPointX[0]), &(aPointY[0]),
nullptr,
gvBurnPoint, &sInfo );
gvBurnPoint, &sInfo,
eMergeAlg == GRMA_Add );
}
else
{
Expand All @@ -513,7 +515,8 @@ gv_rasterize_one_shape( unsigned char *pabyChunkBuf, int nXOff, int nYOff,
static_cast<int>(aPartSize.size()), &(aPartSize[0]),
&(aPointX[0]), &(aPointY[0]),
&(aPointVariant[0]),
gvBurnPoint, &sInfo );
gvBurnPoint, &sInfo,
eMergeAlg == GRMA_Add );
}
}
}
Expand Down
95 changes: 93 additions & 2 deletions gdal/alg/llrasterize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include <cstring>

#include <algorithm>
#include <set>
#include <utility>

#include "gdal_alg.h"
Expand Down Expand Up @@ -322,6 +323,14 @@ void GDALdllImageLine( int nRasterXSize, int nRasterYSize,
? 0.0
: (dfVariant1 - dfVariant) / nDeltaX;

// Do not burn the end point, unless we are in the last
// segment. This is to avoid burning twice intermediate points,
// which causes artifacts in Add mode
if( j != panPartSize[i] - 1 )
{
nDeltaX --;
}

while( nDeltaX-- >= 0 )
{
if( 0 <= iX && iX < nRasterXSize
Expand Down Expand Up @@ -354,6 +363,14 @@ void GDALdllImageLine( int nRasterXSize, int nRasterYSize,
? 0.0
: (dfVariant1 - dfVariant) / nDeltaY;

// Do not burn the end point, unless we are in the last
// segment. This is to avoid burning twice intermediate points,
// which causes artifacts in Add mode
if( j != panPartSize[i] - 1 )
{
nDeltaY --;
}

while( nDeltaY-- >= 0 )
{
if( 0 <= iX && iX < nRasterXSize
Expand Down Expand Up @@ -394,16 +411,23 @@ void
GDALdllImageLineAllTouched( int nRasterXSize, int nRasterYSize,
int nPartCount, int *panPartSize,
double *padfX, double *padfY, double *padfVariant,
llPointFunc pfnPointFunc, void *pCBData )
llPointFunc pfnPointFunc, void *pCBData,
int bAvoidBurningSamePoints )

{
if( !nPartCount )
return;

for( int i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
{
std::set<std::pair<int,int>> lastBurntPoints;
std::set<std::pair<int,int>> newBurntPoints;

for( int j = 1; j < panPartSize[i]; j++ )
{
lastBurntPoints = std::move(newBurntPoints);
newBurntPoints.clear();

double dfX = padfX[n + j - 1];
double dfY = padfY[n + j - 1];

Expand Down Expand Up @@ -465,11 +489,37 @@ GDALdllImageLineAllTouched( int nRasterXSize, int nRasterYSize,
dfVariant += dfDeltaVariant * (iY - dfY);

if( padfVariant == nullptr )
{
for( ; iY <= iYEnd; iY++ )
{
if( bAvoidBurningSamePoints )
{
auto yx = std::pair<int,int>(iY,iX);
if( lastBurntPoints.find( yx ) != lastBurntPoints.end() )
{
continue;
}
newBurntPoints.insert(yx);
}
pfnPointFunc( pCBData, iY, iX, 0.0 );
}
}
else
{
for( ; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant )
{
if( bAvoidBurningSamePoints )
{
auto yx = std::pair<int,int>(iY,iX);
if( lastBurntPoints.find( yx ) != lastBurntPoints.end() )
{
continue;
}
newBurntPoints.insert(yx);
}
pfnPointFunc( pCBData, iY, iX, dfVariant );
}
}

continue; // Next segment.
}
Expand Down Expand Up @@ -502,11 +552,37 @@ GDALdllImageLineAllTouched( int nRasterXSize, int nRasterYSize,
dfVariant += dfDeltaVariant * (iX - dfX);

if( padfVariant == nullptr )
{
for( ; iX <= iXEnd; iX++ )
{
if( bAvoidBurningSamePoints )
{
auto yx = std::pair<int,int>(iY,iX);
if( lastBurntPoints.find( yx ) != lastBurntPoints.end() )
{
continue;
}
newBurntPoints.insert(yx);
}
pfnPointFunc( pCBData, iY, iX, 0.0 );
}
}
else
{
for( ; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant )
{
if( bAvoidBurningSamePoints )
{
auto yx = std::pair<int,int>(iY,iX);
if( lastBurntPoints.find( yx ) != lastBurntPoints.end() )
{
continue;
}
newBurntPoints.insert(yx);
}
pfnPointFunc( pCBData, iY, iX, dfVariant );
}
}

continue; // Next segment.
}
Expand Down Expand Up @@ -576,7 +652,22 @@ GDALdllImageLineAllTouched( int nRasterXSize, int nRasterYSize,
// We should be able to drop the Y check because we clipped
// in Y, but there may be some error with all the small steps.
if( iY >= 0 && iY < nRasterYSize )
pfnPointFunc( pCBData, iY, iX, dfVariant );
{
if( bAvoidBurningSamePoints )
{
auto yx = std::pair<int,int>(iY,iX);
if( lastBurntPoints.find( yx ) == lastBurntPoints.end() &&
newBurntPoints.find( yx ) == newBurntPoints.end() )
{
newBurntPoints.insert(yx);
pfnPointFunc( pCBData, iY, iX, dfVariant );
}
}
else
{
pfnPointFunc( pCBData, iY, iX, dfVariant );
}
}

double dfStepX = floor(dfX+1.0) - dfX;
double dfStepY = dfStepX * dfSlope;
Expand Down

0 comments on commit c5565eb

Please sign in to comment.