Skip to content

Commit

Permalink
Merge pull request #8444 from rouault/fix_8437
Browse files Browse the repository at this point in the history 
Rasterize: avoid burning the same pixel twice in MERGE_ALG=ADD mode within a same geometry (fixes #8437)
  • Loading branch information
@rouault
rouault authored Sep 22, 2023
2 parents 68c32d3 + 03da2f2 commit 3ab1115
Show file tree
Hide file tree
Showing 4 changed files with 211 additions and 30 deletions.
9 changes: 7 additions & 2 deletions alg/gdal_alg_priv.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@

#include <cstdint>

#include <set>

#include "gdal_alg.h"
#include "ogr_spatialref.h"

Expand Down Expand Up @@ -73,6 +75,8 @@ typedef struct
} burnValues;
GDALBurnValueSrc eBurnValueSource;
GDALRasterMergeAlg eMergeAlg;
bool bFillSetVisitedPoints;
std::set<uint64_t> *poSetVisitedPoints;
} GDALRasterizeInfo;

typedef enum
Expand Down Expand Up @@ -104,14 +108,15 @@ void GDALdllImageLineAllTouched(int nRasterXSize, int nRasterYSize,
const double *padfX, const double *padfY,
const double *padfVariant,
llPointFunc pfnPointFunc, void *pCBData,
int bAvoidBurningSamePoints,
bool bAvoidBurningSamePoints,
bool bIntersectOnly);

void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
int nPartCount, const int *panPartSize,
const double *padfX, const double *padfY,
const double *padfVariant,
llScanlineFunc pfnScanlineFunc, void *pCBData);
llScanlineFunc pfnScanlineFunc, void *pCBData,
bool bAvoidBurningSamePoints);

CPL_C_END

Expand Down
123 changes: 101 additions & 22 deletions alg/gdalrasterize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,15 @@ template <typename T> static inline T SaturatedAddSigned(T a, T b)
}
}

/************************************************************************/
/* MakeKey() */
/************************************************************************/

inline uint64_t MakeKey(int y, int x)
{
return (static_cast<uint64_t>(y) << 32) | static_cast<uint64_t>(x);
}

/************************************************************************/
/* gvBurnScanlineBasic() */
/************************************************************************/
Expand All @@ -90,23 +99,43 @@ static inline void gvBurnScanlineBasic(GDALRasterizeInfo *psInfo, int nY,
unsigned char *pabyInsert =
psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
nY * psInfo->nLineSpace + nXStart * psInfo->nPixelSpace;
int nPixels = nXEnd - nXStart + 1;
if (psInfo->eMergeAlg == GRMA_Add)
{
while (nPixels-- > 0)
if (psInfo->poSetVisitedPoints)
{
double dfVal =
static_cast<double>(*reinterpret_cast<T *>(pabyInsert)) +
burnValue;
GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
pabyInsert += psInfo->nPixelSpace;
CPLAssert(!psInfo->bFillSetVisitedPoints);
uint64_t nKey = MakeKey(nY, nXStart);
auto &oSetVisitedPoints = *(psInfo->poSetVisitedPoints);
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
if (oSetVisitedPoints.find(nKey) == oSetVisitedPoints.end())
{
double dfVal = static_cast<double>(
*reinterpret_cast<T *>(pabyInsert)) +
burnValue;
GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
}
pabyInsert += psInfo->nPixelSpace;
++nKey;
}
}
else
{
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
double dfVal = static_cast<double>(
*reinterpret_cast<T *>(pabyInsert)) +
burnValue;
GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
pabyInsert += psInfo->nPixelSpace;
}
}
}
else
{
T nVal;
GDALCopyWord(burnValue, nVal);
while (nPixels-- > 0)
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
*reinterpret_cast<T *>(pabyInsert) = nVal;
pabyInsert += psInfo->nPixelSpace;
Expand All @@ -127,21 +156,41 @@ static inline void gvBurnScanlineInt64UserBurnValue(GDALRasterizeInfo *psInfo,
unsigned char *pabyInsert =
psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
nY * psInfo->nLineSpace + nXStart * psInfo->nPixelSpace;
int nPixels = nXEnd - nXStart + 1;
if (psInfo->eMergeAlg == GRMA_Add)
{
while (nPixels-- > 0)
if (psInfo->poSetVisitedPoints)
{
*reinterpret_cast<std::int64_t *>(pabyInsert) =
SaturatedAddSigned(
*reinterpret_cast<std::int64_t *>(pabyInsert),
burnValue);
pabyInsert += psInfo->nPixelSpace;
CPLAssert(!psInfo->bFillSetVisitedPoints);
uint64_t nKey = MakeKey(nY, nXStart);
auto &oSetVisitedPoints = *(psInfo->poSetVisitedPoints);
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
if (oSetVisitedPoints.find(nKey) == oSetVisitedPoints.end())
{
*reinterpret_cast<std::int64_t *>(pabyInsert) =
SaturatedAddSigned(
*reinterpret_cast<std::int64_t *>(pabyInsert),
burnValue);
}
pabyInsert += psInfo->nPixelSpace;
++nKey;
}
}
else
{
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
*reinterpret_cast<std::int64_t *>(pabyInsert) =
SaturatedAddSigned(
*reinterpret_cast<std::int64_t *>(pabyInsert),
burnValue);
pabyInsert += psInfo->nPixelSpace;
}
}
}
else
{
while (nPixels-- > 0)
for (int nX = nXStart; nX <= nXEnd; ++nX)
{
*reinterpret_cast<std::int64_t *>(pabyInsert) = burnValue;
pabyInsert += psInfo->nPixelSpace;
Expand Down Expand Up @@ -285,6 +334,21 @@ static void gvBurnPoint(void *pCBData, int nY, int nX, double dfVariant)
CPLAssert(nY >= 0 && nY < psInfo->nYSize);
CPLAssert(nX >= 0 && nX < psInfo->nXSize);

if (psInfo->poSetVisitedPoints)
{
const uint64_t nKey = MakeKey(nY, nX);
if (psInfo->poSetVisitedPoints->find(nKey) ==
psInfo->poSetVisitedPoints->end())
{
if (psInfo->bFillSetVisitedPoints)
psInfo->poSetVisitedPoints->insert(nKey);
}
else
{
return;
}
}

if (psInfo->eBurnValueType == GDT_Int64)
{
if (psInfo->eType == GDT_Int64 &&
Expand Down Expand Up @@ -552,6 +616,8 @@ static void gv_rasterize_one_shape(
}
sInfo.eBurnValueSource = eBurnValueSrc;
sInfo.eMergeAlg = eMergeAlg;
sInfo.bFillSetVisitedPoints = false;
sInfo.poSetVisitedPoints = nullptr;

/* -------------------------------------------------------------------- */
/* Transform polygon geometries into a set of rings and a part */
Expand Down Expand Up @@ -607,6 +673,11 @@ static void gv_rasterize_one_shape(
case wkbLineString:
case wkbMultiLineString:
{
if (eMergeAlg == GRMA_Add)
{
sInfo.bFillSetVisitedPoints = true;
sInfo.poSetVisitedPoints = new std::set<uint64_t>();
}
if (bAllTouched)
GDALdllImageLineAllTouched(
sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
Expand All @@ -626,12 +697,11 @@ static void gv_rasterize_one_shape(

default:
{
GDALdllImageFilledPolygon(
sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
aPartSize.data(), aPointX.data(), aPointY.data(),
(eBurnValueSrc == GBV_UserBurnValue) ? nullptr
: aPointVariant.data(),
gvBurnScanline, &sInfo);
if (eMergeAlg == GRMA_Add)
{
sInfo.bFillSetVisitedPoints = true;
sInfo.poSetVisitedPoints = new std::set<uint64_t>();
}
if (bAllTouched)
{
// Reverting the variants to the first value because the
Expand Down Expand Up @@ -662,9 +732,18 @@ static void gv_rasterize_one_shape(
gvBurnPoint, &sInfo, eMergeAlg == GRMA_Add, true);
}
}
sInfo.bFillSetVisitedPoints = false;
GDALdllImageFilledPolygon(
sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
aPartSize.data(), aPointX.data(), aPointY.data(),
(eBurnValueSrc == GBV_UserBurnValue) ? nullptr
: aPointVariant.data(),
gvBurnScanline, &sInfo, eMergeAlg == GRMA_Add);
}
break;
}

delete sInfo.poSetVisitedPoints;
}

/************************************************************************/
Expand Down
43 changes: 37 additions & 6 deletions alg/llrasterize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,8 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
int nPartCount, const int *panPartSize,
const double *padfX, const double *padfY,
const double *dfVariant,
llScanlineFunc pfnScanlineFunc, void *pCBData)
llScanlineFunc pfnScanlineFunc, void *pCBData,
bool bAvoidBurningSamePoints)
{
if (!nPartCount)
{
Expand All @@ -89,6 +90,9 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
n += panPartSize[part];

std::vector<int> polyInts(n);
std::vector<int> polyInts2;
if (bAvoidBurningSamePoints)
polyInts2.resize(n);

double dminy = padfY[0];
double dmaxy = padfY[0];
Expand Down Expand Up @@ -123,6 +127,7 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,

int part = 0;
int ints = 0;
int ints2 = 0;

for (int i = 0; i < n; i++)
{
Expand Down Expand Up @@ -169,7 +174,6 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,

// AE: DO NOT skip bottom horizontal segments
// -Fill them separately-
// They are not taken into account twice.
if (padfX[ind1] > padfX[ind2])
{
const int horizontal_x1 =
Expand All @@ -181,9 +185,17 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
continue;

// Fill the horizontal segment (separately from the rest).
pfnScanlineFunc(pCBData, y, horizontal_x1,
horizontal_x2 - 1,
(dfVariant == nullptr) ? 0 : dfVariant[0]);
if (bAvoidBurningSamePoints)
{
polyInts2[ints2++] = horizontal_x1;
polyInts2[ints2++] = horizontal_x2;
}
else
{
pfnScanlineFunc(
pCBData, y, horizontal_x1, horizontal_x2 - 1,
dfVariant == nullptr ? 0 : dfVariant[0]);
}
}
// else: Skip top horizontal segments.
// They are already filled in the regular loop.
Expand All @@ -200,6 +212,7 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
}

std::sort(polyInts.begin(), polyInts.begin() + ints);
std::sort(polyInts2.begin(), polyInts2.begin() + ints2);

for (int i = 0; i + 1 < ints; i += 2)
{
Expand All @@ -209,6 +222,24 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
dfVariant == nullptr ? 0 : dfVariant[0]);
}
}

for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
{
if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
{
// "synchronize" polyInts[i] with polyInts2[i2]
while (i + 1 < ints && polyInts[i] < polyInts2[i2])
i += 2;
// Only burn if we don't have a common segment between
// polyInts[] and polyInts2[]
if (i + 1 >= ints || polyInts[i] != polyInts2[i2])
{
pfnScanlineFunc(pCBData, y, polyInts2[i2],
polyInts2[i2 + 1] - 1,
dfVariant == nullptr ? 0 : dfVariant[0]);
}
}
}
}
}

Expand Down Expand Up @@ -372,7 +403,7 @@ void GDALdllImageLineAllTouched(int nRasterXSize, int nRasterYSize,
const double *padfX, const double *padfY,
const double *padfVariant,
llPointFunc pfnPointFunc, void *pCBData,
int bAvoidBurningSamePoints,
bool bAvoidBurningSamePoints,
bool bIntersectOnly)

{
Expand Down
Loading

0 comments on commit 3ab1115

Please sign in to comment.