Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rasterize with MERGE_ALG=ADD: avoid burning several times intermediate points of linestrings (fixes #1307) #1309

Merged
merged 1 commit into from
Feb 22, 2019
Merged

Rasterize with MERGE_ALG=ADD: avoid burning several times intermediate points of linestrings (fixes #1307) #1309

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions autotest/alg/rasterize.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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