Skip to content

Commit 6074ce4

Browse files
rouaultrouault
authored
Merge pull request OSGeo#11713 from rouault/fix_11711
Warping: fix inconsistent replacement of valid value that collides with the dstnodata value
2 parents 4622512 + 2b0b1c5 commit 6074ce4

File tree

4 files changed

+179
-88
lines changed

4 files changed

+179
-88
lines changed

‎alg/gdalwarper.h

+2
Original file line numberDiff line numberDiff line change
@@ -463,6 +463,8 @@ class CPL_DLL GDALWarpKernel
463463

464464
GWKTieStrategy eTieStrategy;
465465

466+
bool bWarnedAboutDstNoDataReplacement = false;
467+
466468
/*! @endcond */
467469

468470
GDALWarpKernel();

‎alg/gdalwarpkernel.cpp

+125-83
Original file line numberDiff line numberDiff line change
@@ -190,6 +190,7 @@ static CPLErr GWKBilinearNoMasksOrDstDensityOnlyDouble(GDALWarpKernel *poWK);
190190
static CPLErr GWKCubicNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);
191191
static CPLErr GWKCubicSplineNoMasksOrDstDensityOnlyShort(GDALWarpKernel *poWK);
192192
static CPLErr GWKNearestShort(GDALWarpKernel *poWK);
193+
static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK);
193194
static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK);
194195
static CPLErr GWKNearestFloat(GDALWarpKernel *poWK);
195196
static CPLErr GWKAverageOrMode(GDALWarpKernel *);
@@ -1307,10 +1308,12 @@ CPLErr GDALWarpKernel::PerformWarp()
13071308
bNoMasksOrDstDensityOnly)
13081309
return GWKBilinearNoMasksOrDstDensityOnlyUShort(this);
13091310

1310-
if ((eWorkingDataType == GDT_Int16 || eWorkingDataType == GDT_UInt16) &&
1311-
eResample == GRA_NearestNeighbour)
1311+
if (eWorkingDataType == GDT_Int16 && eResample == GRA_NearestNeighbour)
13121312
return GWKNearestShort(this);
13131313

1314+
if (eWorkingDataType == GDT_Int16 && eResample == GRA_NearestNeighbour)
1315+
return GWKNearestUnsignedShort(this);
1316+
13141317
if (eWorkingDataType == GDT_Float32 && eResample == GRA_NearestNeighbour &&
13151318
bNoMasksOrDstDensityOnly)
13161319
return GWKNearestNoMasksOrDstDensityOnlyFloat(this);
@@ -1520,6 +1523,60 @@ template <> double GWKClampValueT<double>(double dfValue)
15201523
}
15211524
#endif
15221525

1526+
/************************************************************************/
1527+
/* AvoidNoData() */
1528+
/************************************************************************/
1529+
1530+
template <class T>
1531+
inline void AvoidNoData(const GDALWarpKernel *poWK, int iBand,
1532+
GPtrDiff_t iDstOffset)
1533+
{
1534+
GByte *pabyDst = poWK->papabyDstImage[iBand];
1535+
T *pDst = reinterpret_cast<T *>(pabyDst);
1536+
1537+
if (poWK->padfDstNoDataReal != nullptr &&
1538+
poWK->padfDstNoDataReal[iBand] == static_cast<double>(pDst[iDstOffset]))
1539+
{
1540+
if constexpr (std::numeric_limits<T>::is_integer)
1541+
{
1542+
if (pDst[iDstOffset] ==
1543+
static_cast<T>(std::numeric_limits<T>::lowest()))
1544+
{
1545+
pDst[iDstOffset] =
1546+
static_cast<T>(std::numeric_limits<T>::lowest() + 1);
1547+
}
1548+
else
1549+
pDst[iDstOffset]--;
1550+
}
1551+
else
1552+
{
1553+
if (pDst[iDstOffset] == std::numeric_limits<T>::max())
1554+
{
1555+
pDst[iDstOffset] =
1556+
std::nextafter(pDst[iDstOffset], static_cast<T>(0));
1557+
}
1558+
else
1559+
{
1560+
pDst[iDstOffset] = std::nextafter(
1561+
pDst[iDstOffset], std::numeric_limits<T>::max());
1562+
}
1563+
}
1564+
1565+
if (!poWK->bWarnedAboutDstNoDataReplacement)
1566+
{
1567+
const_cast<GDALWarpKernel *>(poWK)
1568+
->bWarnedAboutDstNoDataReplacement = true;
1569+
CPLError(CE_Warning, CPLE_AppDefined,
1570+
"Value %g in the source dataset has been changed to %g "
1571+
"in the destination dataset to avoid being treated as "
1572+
"NoData. To avoid this, select a different NoData value "
1573+
"for the destination dataset.",
1574+
poWK->padfDstNoDataReal[iBand],
1575+
static_cast<double>(pDst[iDstOffset]));
1576+
}
1577+
}
1578+
}
1579+
15231580
/************************************************************************/
15241581
/* GWKSetPixelValueRealT() */
15251582
/************************************************************************/
@@ -1580,16 +1637,39 @@ static bool GWKSetPixelValueRealT(const GDALWarpKernel *poWK, int iBand,
15801637
pDst[iDstOffset] = value;
15811638
}
15821639

1583-
if (poWK->padfDstNoDataReal != nullptr &&
1584-
poWK->padfDstNoDataReal[iBand] == static_cast<double>(pDst[iDstOffset]))
1640+
AvoidNoData<T>(poWK, iBand, iDstOffset);
1641+
1642+
return true;
1643+
}
1644+
1645+
/************************************************************************/
1646+
/* ClampRoundAndAvoidNoData() */
1647+
/************************************************************************/
1648+
1649+
template <class T>
1650+
inline void ClampRoundAndAvoidNoData(const GDALWarpKernel *poWK, int iBand,
1651+
GPtrDiff_t iDstOffset, double dfReal)
1652+
{
1653+
GByte *pabyDst = poWK->papabyDstImage[iBand];
1654+
T *pDst = reinterpret_cast<T *>(pabyDst);
1655+
1656+
if constexpr (std::numeric_limits<T>::is_integer)
15851657
{
1586-
if (pDst[iDstOffset] == std::numeric_limits<T>::min())
1587-
pDst[iDstOffset] = std::numeric_limits<T>::min() + 1;
1658+
if (dfReal < static_cast<double>(std::numeric_limits<T>::lowest()))
1659+
pDst[iDstOffset] = static_cast<T>(std::numeric_limits<T>::lowest());
1660+
else if (dfReal > static_cast<double>(std::numeric_limits<T>::max()))
1661+
pDst[iDstOffset] = static_cast<T>(std::numeric_limits<T>::max());
15881662
else
1589-
pDst[iDstOffset]--;
1663+
pDst[iDstOffset] = (std::numeric_limits<T>::is_signed)
1664+
? static_cast<T>(floor(dfReal + 0.5))
1665+
: static_cast<T>(dfReal + 0.5);
1666+
}
1667+
else
1668+
{
1669+
pDst[iDstOffset] = static_cast<T>(dfReal);
15901670
}
15911671

1592-
return true;
1672+
AvoidNoData<T>(poWK, iBand, iDstOffset);
15931673
}
15941674

15951675
/************************************************************************/
@@ -1724,83 +1804,55 @@ static bool GWKSetPixelValue(const GDALWarpKernel *poWK, int iBand,
17241804
(dfDensity + dfDstInfluence);
17251805
}
17261806

1727-
/* -------------------------------------------------------------------- */
1728-
/* Actually apply the destination value. */
1729-
/* */
1730-
/* Avoid using the destination nodata value for integer datatypes */
1731-
/* if by chance it is equal to the computed pixel value. */
1732-
/* -------------------------------------------------------------------- */
1733-
1734-
// TODO(schwehr): Can we make this a template?
1735-
#define CLAMP(type) \
1736-
do \
1737-
{ \
1738-
type *_pDst = reinterpret_cast<type *>(pabyDst); \
1739-
if (dfReal < static_cast<double>(std::numeric_limits<type>::min())) \
1740-
_pDst[iDstOffset] = \
1741-
static_cast<type>(std::numeric_limits<type>::min()); \
1742-
else if (dfReal > \
1743-
static_cast<double>(std::numeric_limits<type>::max())) \
1744-
_pDst[iDstOffset] = \
1745-
static_cast<type>(std::numeric_limits<type>::max()); \
1746-
else \
1747-
_pDst[iDstOffset] = (std::numeric_limits<type>::is_signed) \
1748-
? static_cast<type>(floor(dfReal + 0.5)) \
1749-
: static_cast<type>(dfReal + 0.5); \
1750-
if (poWK->padfDstNoDataReal != nullptr && \
1751-
poWK->padfDstNoDataReal[iBand] == \
1752-
static_cast<double>(_pDst[iDstOffset])) \
1753-
{ \
1754-
if (_pDst[iDstOffset] == \
1755-
static_cast<type>(std::numeric_limits<type>::min())) \
1756-
_pDst[iDstOffset] = \
1757-
static_cast<type>(std::numeric_limits<type>::min() + 1); \
1758-
else \
1759-
_pDst[iDstOffset]--; \
1760-
} \
1761-
} while (false)
1807+
/* -------------------------------------------------------------------- */
1808+
/* Actually apply the destination value. */
1809+
/* */
1810+
/* Avoid using the destination nodata value for integer datatypes */
1811+
/* if by chance it is equal to the computed pixel value. */
1812+
/* -------------------------------------------------------------------- */
17621813

17631814
switch (poWK->eWorkingDataType)
17641815
{
17651816
case GDT_Byte:
1766-
CLAMP(GByte);
1817+
ClampRoundAndAvoidNoData<GByte>(poWK, iBand, iDstOffset, dfReal);
17671818
break;
17681819

17691820
case GDT_Int8:
1770-
CLAMP(GInt8);
1821+
ClampRoundAndAvoidNoData<GInt8>(poWK, iBand, iDstOffset, dfReal);
17711822
break;
17721823

17731824
case GDT_Int16:
1774-
CLAMP(GInt16);
1825+
ClampRoundAndAvoidNoData<GInt16>(poWK, iBand, iDstOffset, dfReal);
17751826
break;
17761827

17771828
case GDT_UInt16:
1778-
CLAMP(GUInt16);
1829+
ClampRoundAndAvoidNoData<GUInt16>(poWK, iBand, iDstOffset, dfReal);
17791830
break;
17801831

17811832
case GDT_UInt32:
1782-
CLAMP(GUInt32);
1833+
ClampRoundAndAvoidNoData<GUInt32>(poWK, iBand, iDstOffset, dfReal);
17831834
break;
17841835

17851836
case GDT_Int32:
1786-
CLAMP(GInt32);
1837+
ClampRoundAndAvoidNoData<GInt32>(poWK, iBand, iDstOffset, dfReal);
17871838
break;
17881839

17891840
case GDT_UInt64:
1790-
CLAMP(std::uint64_t);
1841+
ClampRoundAndAvoidNoData<std::uint64_t>(poWK, iBand, iDstOffset,
1842+
dfReal);
17911843
break;
17921844

17931845
case GDT_Int64:
1794-
CLAMP(std::int64_t);
1846+
ClampRoundAndAvoidNoData<std::int64_t>(poWK, iBand, iDstOffset,
1847+
dfReal);
17951848
break;
17961849

17971850
case GDT_Float32:
1798-
reinterpret_cast<float *>(pabyDst)[iDstOffset] =
1799-
static_cast<float>(dfReal);
1851+
ClampRoundAndAvoidNoData<float>(poWK, iBand, iDstOffset, dfReal);
18001852
break;
18011853

18021854
case GDT_Float64:
1803-
reinterpret_cast<double *>(pabyDst)[iDstOffset] = dfReal;
1855+
ClampRoundAndAvoidNoData<double>(poWK, iBand, iDstOffset, dfReal);
18041856
break;
18051857

18061858
case GDT_CInt16:
@@ -1983,44 +2035,45 @@ static bool GWKSetPixelValueReal(const GDALWarpKernel *poWK, int iBand,
19832035
switch (poWK->eWorkingDataType)
19842036
{
19852037
case GDT_Byte:
1986-
CLAMP(GByte);
2038+
ClampRoundAndAvoidNoData<GByte>(poWK, iBand, iDstOffset, dfReal);
19872039
break;
19882040

19892041
case GDT_Int8:
1990-
CLAMP(GInt8);
2042+
ClampRoundAndAvoidNoData<GInt8>(poWK, iBand, iDstOffset, dfReal);
19912043
break;
19922044

19932045
case GDT_Int16:
1994-
CLAMP(GInt16);
2046+
ClampRoundAndAvoidNoData<GInt16>(poWK, iBand, iDstOffset, dfReal);
19952047
break;
19962048

19972049
case GDT_UInt16:
1998-
CLAMP(GUInt16);
2050+
ClampRoundAndAvoidNoData<GUInt16>(poWK, iBand, iDstOffset, dfReal);
19992051
break;
20002052

20012053
case GDT_UInt32:
2002-
CLAMP(GUInt32);
2054+
ClampRoundAndAvoidNoData<GUInt32>(poWK, iBand, iDstOffset, dfReal);
20032055
break;
20042056

20052057
case GDT_Int32:
2006-
CLAMP(GInt32);
2058+
ClampRoundAndAvoidNoData<GInt32>(poWK, iBand, iDstOffset, dfReal);
20072059
break;
20082060

20092061
case GDT_UInt64:
2010-
CLAMP(std::uint64_t);
2062+
ClampRoundAndAvoidNoData<std::uint64_t>(poWK, iBand, iDstOffset,
2063+
dfReal);
20112064
break;
20122065

20132066
case GDT_Int64:
2014-
CLAMP(std::int64_t);
2067+
ClampRoundAndAvoidNoData<std::int64_t>(poWK, iBand, iDstOffset,
2068+
dfReal);
20152069
break;
20162070

20172071
case GDT_Float32:
2018-
reinterpret_cast<float *>(pabyDst)[iDstOffset] =
2019-
static_cast<float>(dfReal);
2072+
ClampRoundAndAvoidNoData<float>(poWK, iBand, iDstOffset, dfReal);
20202073
break;
20212074

20222075
case GDT_Float64:
2023-
reinterpret_cast<double *>(pabyDst)[iDstOffset] = dfReal;
2076+
ClampRoundAndAvoidNoData<double>(poWK, iBand, iDstOffset, dfReal);
20242077
break;
20252078

20262079
case GDT_CInt16:
@@ -6678,24 +6731,8 @@ template <class T> static void GWKNearestThread(void *pData)
66786731
padfZ[iDstX] * dfMultFactorVerticalShiftPipeline);
66796732
}
66806733

6681-
if (dfBandDensity < 1.0)
6682-
{
6683-
if (dfBandDensity == 0.0)
6684-
{
6685-
// Do nothing.
6686-
}
6687-
else
6688-
{
6689-
// Let the general code take care of mixing.
6690-
GWKSetPixelValueRealT(poWK, iBand, iDstOffset,
6691-
dfBandDensity, value);
6692-
}
6693-
}
6694-
else
6695-
{
6696-
reinterpret_cast<T *>(
6697-
poWK->papabyDstImage[iBand])[iDstOffset] = value;
6698-
}
6734+
GWKSetPixelValueRealT(poWK, iBand, iDstOffset,
6735+
dfBandDensity, value);
66996736
}
67006737
}
67016738

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

6850+
static CPLErr GWKNearestUnsignedShort(GDALWarpKernel *poWK)
6851+
{
6852+
return GWKRun(poWK, "GWKNearestUnsignedShort", GWKNearestThread<GUInt16>);
6853+
}
6854+
68136855
static CPLErr GWKNearestNoMasksOrDstDensityOnlyFloat(GDALWarpKernel *poWK)
68146856
{
68156857
return GWKRun(

‎autotest/alg/applyverticalshiftgrid.py

+6-2
Original file line numberDiff line numberDiff line change
@@ -265,7 +265,9 @@ def test_applyverticalshiftgrid_5():
265265
outputType=gdal.GDT_Float32,
266266
scaleParams=[[0, 1, 0, 0.5]],
267267
)
268-
out_ds = gdal.ApplyVerticalShiftGrid(src_ds, grid_ds, srcUnitToMeter=2)
268+
out_ds = gdal.ApplyVerticalShiftGrid(
269+
src_ds, grid_ds, srcUnitToMeter=2, options=["ERROR_ON_MISSING_VERT_SHIFT=YES"]
270+
)
269271
cs = out_ds.GetRasterBand(1).Checksum()
270272
assert cs == 4672
271273

@@ -279,7 +281,9 @@ def test_applyverticalshiftgrid_5():
279281
outputType=gdal.GDT_Float32,
280282
scaleParams=[[0, 1, 0, 0.5]],
281283
)
282-
out_ds = gdal.ApplyVerticalShiftGrid(src_ds, grid_ds, dstUnitToMeter=0.5)
284+
out_ds = gdal.ApplyVerticalShiftGrid(
285+
src_ds, grid_ds, dstUnitToMeter=0.5, options=["ERROR_ON_MISSING_VERT_SHIFT=YES"]
286+
)
283287
cs = out_ds.GetRasterBand(1).Checksum()
284288
assert cs == 4672
285289

0 commit comments

Comments
 (0)