From 03da2f29f6e2f3d792d15f3a3bdab7e06add5fcd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 21 Sep 2023 00:40:41 +0200 Subject: [PATCH] Rasterize: avoid burning the same pixel twice in MERGE_ALG=ADD mode within a same geometry (fixes #8437) --- alg/gdal_alg_priv.h | 9 ++- alg/gdalrasterize.cpp | 123 +++++++++++++++++++++++++++++++------- alg/llrasterize.cpp | 43 +++++++++++-- autotest/alg/rasterize.py | 66 ++++++++++++++++++++ 4 files changed, 211 insertions(+), 30 deletions(-) diff --git a/alg/gdal_alg_priv.h b/alg/gdal_alg_priv.h index 406d00365843..55dc351b4f96 100644 --- a/alg/gdal_alg_priv.h +++ b/alg/gdal_alg_priv.h @@ -36,6 +36,8 @@ #include +#include + #include "gdal_alg.h" #include "ogr_spatialref.h" @@ -73,6 +75,8 @@ typedef struct } burnValues; GDALBurnValueSrc eBurnValueSource; GDALRasterMergeAlg eMergeAlg; + bool bFillSetVisitedPoints; + std::set *poSetVisitedPoints; } GDALRasterizeInfo; typedef enum @@ -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 diff --git a/alg/gdalrasterize.cpp b/alg/gdalrasterize.cpp index 7060a5b06db8..6787184b8822 100644 --- a/alg/gdalrasterize.cpp +++ b/alg/gdalrasterize.cpp @@ -73,6 +73,15 @@ template static inline T SaturatedAddSigned(T a, T b) } } +/************************************************************************/ +/* MakeKey() */ +/************************************************************************/ + +inline uint64_t MakeKey(int y, int x) +{ + return (static_cast(y) << 32) | static_cast(x); +} + /************************************************************************/ /* gvBurnScanlineBasic() */ /************************************************************************/ @@ -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(*reinterpret_cast(pabyInsert)) + - burnValue; - GDALCopyWord(dfVal, *reinterpret_cast(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( + *reinterpret_cast(pabyInsert)) + + burnValue; + GDALCopyWord(dfVal, *reinterpret_cast(pabyInsert)); + } + pabyInsert += psInfo->nPixelSpace; + ++nKey; + } + } + else + { + for (int nX = nXStart; nX <= nXEnd; ++nX) + { + double dfVal = static_cast( + *reinterpret_cast(pabyInsert)) + + burnValue; + GDALCopyWord(dfVal, *reinterpret_cast(pabyInsert)); + pabyInsert += psInfo->nPixelSpace; + } } } else { T nVal; GDALCopyWord(burnValue, nVal); - while (nPixels-- > 0) + for (int nX = nXStart; nX <= nXEnd; ++nX) { *reinterpret_cast(pabyInsert) = nVal; pabyInsert += psInfo->nPixelSpace; @@ -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(pabyInsert) = - SaturatedAddSigned( - *reinterpret_cast(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(pabyInsert) = + SaturatedAddSigned( + *reinterpret_cast(pabyInsert), + burnValue); + } + pabyInsert += psInfo->nPixelSpace; + ++nKey; + } + } + else + { + for (int nX = nXStart; nX <= nXEnd; ++nX) + { + *reinterpret_cast(pabyInsert) = + SaturatedAddSigned( + *reinterpret_cast(pabyInsert), + burnValue); + pabyInsert += psInfo->nPixelSpace; + } } } else { - while (nPixels-- > 0) + for (int nX = nXStart; nX <= nXEnd; ++nX) { *reinterpret_cast(pabyInsert) = burnValue; pabyInsert += psInfo->nPixelSpace; @@ -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 && @@ -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 */ @@ -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(); + } if (bAllTouched) GDALdllImageLineAllTouched( sInfo.nXSize, nYSize, static_cast(aPartSize.size()), @@ -626,12 +697,11 @@ static void gv_rasterize_one_shape( default: { - GDALdllImageFilledPolygon( - sInfo.nXSize, nYSize, static_cast(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(); + } if (bAllTouched) { // Reverting the variants to the first value because the @@ -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(aPartSize.size()), + aPartSize.data(), aPointX.data(), aPointY.data(), + (eBurnValueSrc == GBV_UserBurnValue) ? nullptr + : aPointVariant.data(), + gvBurnScanline, &sInfo, eMergeAlg == GRMA_Add); } break; } + + delete sInfo.poSetVisitedPoints; } /************************************************************************/ diff --git a/alg/llrasterize.cpp b/alg/llrasterize.cpp index ebf2642d5a4f..6444a0b098ae 100644 --- a/alg/llrasterize.cpp +++ b/alg/llrasterize.cpp @@ -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) { @@ -89,6 +90,9 @@ void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize, n += panPartSize[part]; std::vector polyInts(n); + std::vector polyInts2; + if (bAvoidBurningSamePoints) + polyInts2.resize(n); double dminy = padfY[0]; double dmaxy = padfY[0]; @@ -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++) { @@ -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 = @@ -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. @@ -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) { @@ -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]); + } + } + } } } @@ -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) { diff --git a/autotest/alg/rasterize.py b/autotest/alg/rasterize.py index c618a87a7e04..59a172225c31 100755 --- a/autotest/alg/rasterize.py +++ b/autotest/alg/rasterize.py @@ -899,3 +899,69 @@ def test_rasterize_bugfix_gh6981(): ) == gdal.CE_None ) + + +############################################################################### + + +@pytest.mark.parametrize( + "wkt", + [ + "POLYGON((12.5 2.5, 12.5 12.5, 2.5 12.5, 2.5 2.5, 12.5 2.5))", + "POLYGON((12.5 2.5, 2.5 2.5, 2.5 12.5, 12.5 12.5, 12.5 2.5))", + "LINESTRING(12.5 2.5, 2.5 2.5, 2.5 12.5, 12.5 12.5, 12.5 2.5)", + ], + ids=["clockwise", "counterclockwise", "linestring"], +) +@pytest.mark.parametrize( + "options", + [ + ["MERGE_ALG=ADD", "ALL_TOUCHED=YES"], + ["ALL_TOUCHED=YES"], + ["MERGE_ALG=ADD"], + [], + ], +) +@pytest.mark.parametrize("nbands", [1, 2]) +def test_rasterize_bugfix_gh8437(wkt, options, nbands): + + # Setup working spatial reference + sr_wkt = 'LOCAL_CS["arbitrary"]' + sr = osr.SpatialReference(sr_wkt) + + # Create a memory raster to rasterize into. + target_ds = gdal.GetDriverByName("MEM").Create("", 15, 15, nbands, gdal.GDT_Byte) + target_ds.SetGeoTransform((15, -1, 0, 15, 0, -1)) + + target_ds.SetProjection(sr_wkt) + + # Create a memory layer to rasterize from. + rast_ogr_ds = ogr.GetDriverByName("Memory").CreateDataSource("wrk") + rast_mem_lyr = rast_ogr_ds.CreateLayer("poly", srs=sr) + + # Add a polygon. + feat = ogr.Feature(rast_mem_lyr.GetLayerDefn()) + feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt)) + + rast_mem_lyr.CreateFeature(feat) + + # Run the algorithm. + gdal.RasterizeLayer( + target_ds, + [i + 1 for i in range(nbands)], + rast_mem_lyr, + burn_values=[80] * nbands, + options=options, + ) + + expected_checksum = ( + 519 + if wkt.startswith("LINESTRING") + else 1727 + if "ALL_TOUCHED=YES" in options + else 1435 + ) + for i in range(nbands): + _, maxval = target_ds.GetRasterBand(i + 1).ComputeRasterMinMax() + assert maxval == 80 + assert target_ds.GetRasterBand(i + 1).Checksum() == expected_checksum