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: avoid burning the same pixel twice in MERGE_ALG=ADD mode within a same geometry (fixes #8437) #8444

Merged
merged 1 commit into from
Sep 22, 2023
+211 −30
Merged

Rasterize: avoid burning the same pixel twice in MERGE_ALG=ADD mode within a same geometry (fixes #8437) #8444

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
9 changes: 7 additions & 2 deletions alg/gdal_alg_priv.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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