Skip to content

Commit

Permalink
Merge pull request #11713 from rouault/fix_11711
Browse files Browse the repository at this point in the history 
Warping: fix inconsistent replacement of valid value that collides with the dstnodata value
  • Loading branch information
@rouault
rouault authored Jan 26, 2025
2 parents 4622512 + 2b0b1c5 commit 6074ce4
Show file tree
Hide file tree
Showing 4 changed files with 179 additions and 88 deletions.
2 changes: 2 additions & 0 deletions alg/gdalwarper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -463,6 +463,8 @@ class CPL_DLL GDALWarpKernel

GWKTieStrategy eTieStrategy;

bool bWarnedAboutDstNoDataReplacement = false;

/*! @endcond */

GDALWarpKernel();
Expand Down
208 changes: 125 additions & 83 deletions alg/gdalwarpkernel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 *);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -1520,6 +1523,60 @@ template <> double GWKClampValueT<double>(double dfValue)
}
#endif

/************************************************************************/
/* AvoidNoData() */
/************************************************************************/

template <class T>
inline void AvoidNoData(const GDALWarpKernel *poWK, int iBand,
GPtrDiff_t iDstOffset)
{
GByte *pabyDst = poWK->papabyDstImage[iBand];
T *pDst = reinterpret_cast<T *>(pabyDst);

if (poWK->padfDstNoDataReal != nullptr &&
poWK->padfDstNoDataReal[iBand] == static_cast<double>(pDst[iDstOffset]))
{
if constexpr (std::numeric_limits<T>::is_integer)
{
if (pDst[iDstOffset] ==
static_cast<T>(std::numeric_limits<T>::lowest()))
{
pDst[iDstOffset] =
static_cast<T>(std::numeric_limits<T>::lowest() + 1);
}
else
pDst[iDstOffset]--;
}
else
{
if (pDst[iDstOffset] == std::numeric_limits<T>::max())
{
pDst[iDstOffset] =
std::nextafter(pDst[iDstOffset], static_cast<T>(0));
}
else
{
pDst[iDstOffset] = std::nextafter(
pDst[iDstOffset], std::numeric_limits<T>::max());
}
}

if (!poWK->bWarnedAboutDstNoDataReplacement)
{
const_cast<GDALWarpKernel *>(poWK)
->bWarnedAboutDstNoDataReplacement = true;
CPLError(CE_Warning, CPLE_AppDefined,
"Value %g in the source dataset has been changed to %g "
"in the destination dataset to avoid being treated as "
"NoData. To avoid this, select a different NoData value "
"for the destination dataset.",
poWK->padfDstNoDataReal[iBand],
static_cast<double>(pDst[iDstOffset]));
}
}
}

/************************************************************************/
/* GWKSetPixelValueRealT() */
/************************************************************************/
Expand Down Expand Up @@ -1580,16 +1637,39 @@ static bool GWKSetPixelValueRealT(const GDALWarpKernel *poWK, int iBand,
pDst[iDstOffset] = value;
}

if (poWK->padfDstNoDataReal != nullptr &&
poWK->padfDstNoDataReal[iBand] == static_cast<double>(pDst[iDstOffset]))
AvoidNoData<T>(poWK, iBand, iDstOffset);

return true;
}

/************************************************************************/
/* ClampRoundAndAvoidNoData() */
/************************************************************************/

template <class T>
inline void ClampRoundAndAvoidNoData(const GDALWarpKernel *poWK, int iBand,
GPtrDiff_t iDstOffset, double dfReal)
{
GByte *pabyDst = poWK->papabyDstImage[iBand];
T *pDst = reinterpret_cast<T *>(pabyDst);

if constexpr (std::numeric_limits<T>::is_integer)
{
if (pDst[iDstOffset] == std::numeric_limits<T>::min())
pDst[iDstOffset] = std::numeric_limits<T>::min() + 1;
if (dfReal < static_cast<double>(std::numeric_limits<T>::lowest()))
pDst[iDstOffset] = static_cast<T>(std::numeric_limits<T>::lowest());
else if (dfReal > static_cast<double>(std::numeric_limits<T>::max()))
pDst[iDstOffset] = static_cast<T>(std::numeric_limits<T>::max());
else
pDst[iDstOffset]--;
pDst[iDstOffset] = (std::numeric_limits<T>::is_signed)
? static_cast<T>(floor(dfReal + 0.5))
: static_cast<T>(dfReal + 0.5);
}
else
{
pDst[iDstOffset] = static_cast<T>(dfReal);
}

return true;
AvoidNoData<T>(poWK, iBand, iDstOffset);
}

/************************************************************************/
Expand Down Expand Up @@ -1724,83 +1804,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<type *>(pabyDst); \
if (dfReal < static_cast<double>(std::numeric_limits<type>::min())) \
_pDst[iDstOffset] = \
static_cast<type>(std::numeric_limits<type>::min()); \
else if (dfReal > \
static_cast<double>(std::numeric_limits<type>::max())) \
_pDst[iDstOffset] = \
static_cast<type>(std::numeric_limits<type>::max()); \
else \
_pDst[iDstOffset] = (std::numeric_limits<type>::is_signed) \
? static_cast<type>(floor(dfReal + 0.5)) \
: static_cast<type>(dfReal + 0.5); \
if (poWK->padfDstNoDataReal != nullptr && \
poWK->padfDstNoDataReal[iBand] == \
static_cast<double>(_pDst[iDstOffset])) \
{ \
if (_pDst[iDstOffset] == \
static_cast<type>(std::numeric_limits<type>::min())) \
_pDst[iDstOffset] = \
static_cast<type>(std::numeric_limits<type>::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<GByte>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int8:
CLAMP(GInt8);
ClampRoundAndAvoidNoData<GInt8>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int16:
CLAMP(GInt16);
ClampRoundAndAvoidNoData<GInt16>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt16:
CLAMP(GUInt16);
ClampRoundAndAvoidNoData<GUInt16>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt32:
CLAMP(GUInt32);
ClampRoundAndAvoidNoData<GUInt32>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int32:
CLAMP(GInt32);
ClampRoundAndAvoidNoData<GInt32>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt64:
CLAMP(std::uint64_t);
ClampRoundAndAvoidNoData<std::uint64_t>(poWK, iBand, iDstOffset,
dfReal);
break;

case GDT_Int64:
CLAMP(std::int64_t);
ClampRoundAndAvoidNoData<std::int64_t>(poWK, iBand, iDstOffset,
dfReal);
break;

case GDT_Float32:
reinterpret_cast<float *>(pabyDst)[iDstOffset] =
static_cast<float>(dfReal);
ClampRoundAndAvoidNoData<float>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Float64:
reinterpret_cast<double *>(pabyDst)[iDstOffset] = dfReal;
ClampRoundAndAvoidNoData<double>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_CInt16:
Expand Down Expand Up @@ -1983,44 +2035,45 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand,
switch (poWK->eWorkingDataType)
{
case GDT_Byte:
CLAMP(GByte);
ClampRoundAndAvoidNoData<GByte>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int8:
CLAMP(GInt8);
ClampRoundAndAvoidNoData<GInt8>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int16:
CLAMP(GInt16);
ClampRoundAndAvoidNoData<GInt16>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt16:
CLAMP(GUInt16);
ClampRoundAndAvoidNoData<GUInt16>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt32:
CLAMP(GUInt32);
ClampRoundAndAvoidNoData<GUInt32>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Int32:
CLAMP(GInt32);
ClampRoundAndAvoidNoData<GInt32>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_UInt64:
CLAMP(std::uint64_t);
ClampRoundAndAvoidNoData<std::uint64_t>(poWK, iBand, iDstOffset,
dfReal);
break;

case GDT_Int64:
CLAMP(std::int64_t);
ClampRoundAndAvoidNoData<std::int64_t>(poWK, iBand, iDstOffset,
dfReal);
break;

case GDT_Float32:
reinterpret_cast<float *>(pabyDst)[iDstOffset] =
static_cast<float>(dfReal);
ClampRoundAndAvoidNoData<float>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_Float64:
reinterpret_cast<double *>(pabyDst)[iDstOffset] = dfReal;
ClampRoundAndAvoidNoData<double>(poWK, iBand, iDstOffset, dfReal);
break;

case GDT_CInt16:
Expand Down Expand Up @@ -6678,24 +6731,8 @@ template <class T> 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<T *>(
poWK->papabyDstImage[iBand])[iDstOffset] = value;
}
GWKSetPixelValueRealT(poWK, iBand, iDstOffset,
dfBandDensity, value);
}
}

Expand Down Expand Up @@ -6810,6 +6847,11 @@ static CPLErr GWKNearestShort(GDALWarpKernel *poWK)
return GWKRun(poWK, "GWKNearestShort", GWKNearestThread<GInt16>);
}

static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK)
{
return GWKRun(poWK, "GWKNearestUnsignedShort", GWKNearestThread<GUInt16>);
}

static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK)
{
return GWKRun(
Expand Down
8 changes: 6 additions & 2 deletions autotest/alg/applyverticalshiftgrid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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

Expand Down
Loading

0 comments on commit 6074ce4

Please sign in to comment.