diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index c233d84ae107..911a620e629f 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -190,6 +190,7 @@ static CPLErr GWKBilinearNoMasksOrDstDensityOnlyDouble(GDALWarpKernel *poWK); static CPLErr GWKCubicNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK); static CPLErr GWKCubicSplineNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK); static CPLErr GWKNearestShort(GDALWarpKernel *poWK); +static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK); static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK); static CPLErr GWKNearestFloat(GDALWarpKernel *poWK); static CPLErr GWKAverageOrMode(GDALWarpKernel *); @@ -1307,10 +1308,12 @@ CPLErr GDALWarpKernel::PerformWarp() bNoMasksOrDstDensityOnly) return GWKBilinearNoMasksOrDstDensityOnlyUShort(this); - if ((eWorkingDataType == GDT_Int16 || eWorkingDataType == GDT_UInt16) && - eResample == GRA_NearestNeighbour) + if (eWorkingDataType == GDT_Int16 && eResample == GRA_NearestNeighbour) return GWKNearestShort(this); + if (eWorkingDataType == GDT_Int16 && eResample == GRA_NearestNeighbour) + return GWKNearestUnsignedShort(this); + if (eWorkingDataType == GDT_Float32 && eResample == GRA_NearestNeighbour && bNoMasksOrDstDensityOnly) return GWKNearestNoMasksOrDstDensityOnlyFloat(this); @@ -1520,6 +1523,47 @@ template <> double GWKClampValueT(double dfValue) } #endif +/************************************************************************/ +/* AvoidNoData() */ +/************************************************************************/ + +template +inline void AvoidNoData(const GDALWarpKernel *poWK, int iBand, + GPtrDiff_t iDstOffset) +{ + GByte *pabyDst = poWK->papabyDstImage[iBand]; + T *pDst = reinterpret_cast(pabyDst); + + if (poWK->padfDstNoDataReal != nullptr && + poWK->padfDstNoDataReal[iBand] == static_cast(pDst[iDstOffset])) + { + if constexpr (std::numeric_limits::is_integer) + { + if (pDst[iDstOffset] == + static_cast(std::numeric_limits::lowest())) + { + pDst[iDstOffset] = + static_cast(std::numeric_limits::lowest() + 1); + } + else + pDst[iDstOffset]--; + } + else + { + if (pDst[iDstOffset] == std::numeric_limits::max()) + { + pDst[iDstOffset] = + std::nextafter(pDst[iDstOffset], static_cast(0)); + } + else + { + pDst[iDstOffset] = std::nextafter( + pDst[iDstOffset], std::numeric_limits::max()); + } + } + } +} + /************************************************************************/ /* GWKSetPixelValueRealT() */ /************************************************************************/ @@ -1580,16 +1624,39 @@ static bool GWKSetPixelValueRealT(const GDALWarpKernel *poWK, int iBand, pDst[iDstOffset] = value; } - if (poWK->padfDstNoDataReal != nullptr && - poWK->padfDstNoDataReal[iBand] == static_cast(pDst[iDstOffset])) + AvoidNoData(poWK, iBand, iDstOffset); + + return true; +} + +/************************************************************************/ +/* ClampRoundAndAvoidNoData() */ +/************************************************************************/ + +template +inline void ClampRoundAndAvoidNoData(const GDALWarpKernel *poWK, int iBand, + GPtrDiff_t iDstOffset, double dfReal) +{ + GByte *pabyDst = poWK->papabyDstImage[iBand]; + T *pDst = reinterpret_cast(pabyDst); + + if constexpr (std::numeric_limits::is_integer) { - if (pDst[iDstOffset] == std::numeric_limits::min()) - pDst[iDstOffset] = std::numeric_limits::min() + 1; + if (dfReal < static_cast(std::numeric_limits::lowest())) + pDst[iDstOffset] = static_cast(std::numeric_limits::lowest()); + else if (dfReal > static_cast(std::numeric_limits::max())) + pDst[iDstOffset] = static_cast(std::numeric_limits::max()); else - pDst[iDstOffset]--; + pDst[iDstOffset] = (std::numeric_limits::is_signed) + ? static_cast(floor(dfReal + 0.5)) + : static_cast(dfReal + 0.5); + } + else + { + pDst[iDstOffset] = static_cast(dfReal); } - return true; + AvoidNoData(poWK, iBand, iDstOffset); } /************************************************************************/ @@ -1724,83 +1791,55 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand, (dfDensity + dfDstInfluence); } -/* -------------------------------------------------------------------- */ -/* Actually apply the destination value. */ -/* */ -/* Avoid using the destination nodata value for integer datatypes */ -/* if by chance it is equal to the computed pixel value. */ -/* -------------------------------------------------------------------- */ - -// TODO(schwehr): Can we make this a template? -#define CLAMP(type) \ - do \ - { \ - type *_pDst = reinterpret_cast(pabyDst); \ - if (dfReal < static_cast(std::numeric_limits::min())) \ - _pDst[iDstOffset] = \ - static_cast(std::numeric_limits::min()); \ - else if (dfReal > \ - static_cast(std::numeric_limits::max())) \ - _pDst[iDstOffset] = \ - static_cast(std::numeric_limits::max()); \ - else \ - _pDst[iDstOffset] = (std::numeric_limits::is_signed) \ - ? static_cast(floor(dfReal + 0.5)) \ - : static_cast(dfReal + 0.5); \ - if (poWK->padfDstNoDataReal != nullptr && \ - poWK->padfDstNoDataReal[iBand] == \ - static_cast(_pDst[iDstOffset])) \ - { \ - if (_pDst[iDstOffset] == \ - static_cast(std::numeric_limits::min())) \ - _pDst[iDstOffset] = \ - static_cast(std::numeric_limits::min() + 1); \ - else \ - _pDst[iDstOffset]--; \ - } \ - } while (false) + /* -------------------------------------------------------------------- */ + /* Actually apply the destination value. */ + /* */ + /* Avoid using the destination nodata value for integer datatypes */ + /* if by chance it is equal to the computed pixel value. */ + /* -------------------------------------------------------------------- */ switch (poWK->eWorkingDataType) { case GDT_Byte: - CLAMP(GByte); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int8: - CLAMP(GInt8); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int16: - CLAMP(GInt16); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt16: - CLAMP(GUInt16); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt32: - CLAMP(GUInt32); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int32: - CLAMP(GInt32); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt64: - CLAMP(std::uint64_t); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, + dfReal); break; case GDT_Int64: - CLAMP(std::int64_t); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, + dfReal); break; case GDT_Float32: - reinterpret_cast(pabyDst)[iDstOffset] = - static_cast(dfReal); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Float64: - reinterpret_cast(pabyDst)[iDstOffset] = dfReal; + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_CInt16: @@ -1983,44 +2022,45 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand, switch (poWK->eWorkingDataType) { case GDT_Byte: - CLAMP(GByte); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int8: - CLAMP(GInt8); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int16: - CLAMP(GInt16); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt16: - CLAMP(GUInt16); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt32: - CLAMP(GUInt32); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Int32: - CLAMP(GInt32); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_UInt64: - CLAMP(std::uint64_t); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, + dfReal); break; case GDT_Int64: - CLAMP(std::int64_t); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, + dfReal); break; case GDT_Float32: - reinterpret_cast(pabyDst)[iDstOffset] = - static_cast(dfReal); + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_Float64: - reinterpret_cast(pabyDst)[iDstOffset] = dfReal; + ClampRoundAndAvoidNoData(poWK, iBand, iDstOffset, dfReal); break; case GDT_CInt16: @@ -6678,24 +6718,8 @@ template static void GWKNearestThread(void *pData) padfZ[iDstX] * dfMultFactorVerticalShiftPipeline); } - if (dfBandDensity < 1.0) - { - if (dfBandDensity == 0.0) - { - // Do nothing. - } - else - { - // Let the general code take care of mixing. - GWKSetPixelValueRealT(poWK, iBand, iDstOffset, - dfBandDensity, value); - } - } - else - { - reinterpret_cast( - poWK->papabyDstImage[iBand])[iDstOffset] = value; - } + GWKSetPixelValueRealT(poWK, iBand, iDstOffset, + dfBandDensity, value); } } @@ -6810,6 +6834,11 @@ static CPLErr GWKNearestShort(GDALWarpKernel *poWK) return GWKRun(poWK, "GWKNearestShort", GWKNearestThread); } +static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK) +{ + return GWKRun(poWK, "GWKNearestUnsignedShort", GWKNearestThread); +} + static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK) { return GWKRun( diff --git a/autotest/alg/applyverticalshiftgrid.py b/autotest/alg/applyverticalshiftgrid.py index 95d8dc8d326a..407d6ae49791 100755 --- a/autotest/alg/applyverticalshiftgrid.py +++ b/autotest/alg/applyverticalshiftgrid.py @@ -265,7 +265,9 @@ def test_applyverticalshiftgrid_5(): outputType=gdal.GDT_Float32, scaleParams=[[0, 1, 0, 0.5]], ) - out_ds = gdal.ApplyVerticalShiftGrid(src_ds, grid_ds, srcUnitToMeter=2) + out_ds = gdal.ApplyVerticalShiftGrid( + src_ds, grid_ds, srcUnitToMeter=2, options=["ERROR_ON_MISSING_VERT_SHIFT=YES"] + ) cs = out_ds.GetRasterBand(1).Checksum() assert cs == 4672 @@ -279,7 +281,9 @@ def test_applyverticalshiftgrid_5(): outputType=gdal.GDT_Float32, scaleParams=[[0, 1, 0, 0.5]], ) - out_ds = gdal.ApplyVerticalShiftGrid(src_ds, grid_ds, dstUnitToMeter=0.5) + out_ds = gdal.ApplyVerticalShiftGrid( + src_ds, grid_ds, dstUnitToMeter=0.5, options=["ERROR_ON_MISSING_VERT_SHIFT=YES"] + ) cs = out_ds.GetRasterBand(1).Checksum() assert cs == 4672 diff --git a/autotest/alg/warp.py b/autotest/alg/warp.py index 24158ceaec6a..b060b327242f 100755 --- a/autotest/alg/warp.py +++ b/autotest/alg/warp.py @@ -694,7 +694,7 @@ def test_warp_24(): ds_ref = gdal.Open("data/test3658.tif") cs_ref = ds_ref.GetRasterBand(1).Checksum() - ds = gdal.Warp("", ds_ref, options="-of MEM -r bilinear") + ds = gdal.Warp("", ds_ref, options="-srcnodata none -of MEM -r bilinear") cs = ds.GetRasterBand(1).Checksum() assert cs == cs_ref, "did not get expected checksum" @@ -1084,7 +1084,11 @@ def test_warp_38(): ds.SetGCPs(gcp_list, gdaltest.user_srs_to_wkt("EPSG:32632")) ds = None - gdal.Warp(out_file, "data/test3658.tif", options="-to DST_METHOD=GCP_POLYNOMIAL") + gdal.Warp( + out_file, + "data/test3658.tif", + options="-srcnodata none -to DST_METHOD=GCP_POLYNOMIAL", + ) ds = gdal.Open(out_file) cs = ds.GetRasterBand(1).Checksum() @@ -1117,7 +1121,9 @@ def test_warp_39(): ds.SetGCPs(gcp_list, gdaltest.user_srs_to_wkt("EPSG:32632")) ds = None - gdal.Warp(out_file, "data/test3658.tif", options="-to DST_METHOD=GCP_TPS") + gdal.Warp( + out_file, "data/test3658.tif", options="-srcnodata none -to DST_METHOD=GCP_TPS" + ) ds = gdal.Open(out_file) cs = ds.GetRasterBand(1).Checksum() @@ -1469,6 +1475,7 @@ def test_warp_53(typestr, option, alg_name, expected_cs): src_ds.GetRasterBand(2).Fill(255) zero = struct.pack("B" * 1, 0) src_ds.GetRasterBand(2).WriteRaster(10, 10, 1, 1, zero, buf_type=gdal.GDT_Byte) + dst_ds = gdal.Translate( "", src_ds, options="-outsize 10 10 -of MEM -a_srs EPSG:32611" ) @@ -1917,3 +1924,39 @@ def test_warp_average_NODATA_VALUES_PCT_THRESHOLD(): options="-of MEM -ts 1 1 -r average -wo NODATA_VALUES_PCT_THRESHOLD=25", ) assert struct.unpack("B", out_ds.ReadRaster())[0] == 20 + + +############################################################################### +# + + +@pytest.mark.parametrize( + "dt,expected_val", + [ + (gdal.GDT_Byte, 1.0), + (gdal.GDT_Int8, -1.0), + (gdal.GDT_UInt16, 1.0), + (gdal.GDT_Int16, -1.0), + (gdal.GDT_UInt32, 1.0), + (gdal.GDT_Int32, -1.0), + (gdal.GDT_UInt64, 1.0), + (gdal.GDT_Int64, -1.0), + (gdal.GDT_Float32, 1.401298464324817e-45), + (gdal.GDT_Float64, 5e-324), + ], +) +@pytest.mark.parametrize("resampling", ["nearest", "bilinear"]) +def test_warp_nodata_substitution(dt, expected_val, resampling): + + src_ds = gdal.GetDriverByName("MEM").Create("", 4, 4, 1, dt) + src_ds.SetGeoTransform([1, 1, 0, 1, 0, 1]) + + out_ds = gdal.Warp( + "", + src_ds, + options=f"-of MEM -dstnodata 0 -r {resampling}", + ) + assert ( + struct.unpack("d", out_ds.ReadRaster(0, 0, 1, 1, buf_type=gdal.GDT_Float64))[0] + == expected_val + )