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

Jump to bottom

Overview / RasterIO resampling: do not use nodata value as a valid output value #1097

Merged
merged 1 commit into from Nov 14, 2018
Merged

Overview / RasterIO resampling: do not use nodata value as a valid output value #1097

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.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
81 changes: 81 additions & 0 deletions autotest/gcore/rasterio.py
View file
Open in desktop
Expand Up @@ -1127,6 +1127,86 @@ def rasterio_lanczos_nodata():

return 'success'

###############################################################################


def rasterio_resampled_value_is_nodata():

gdal.FileFromMemBuffer('/vsimem/in.asc',
"""ncols 4
nrows 4
xllcorner 440720.000000000000
yllcorner 3750120.000000000000
cellsize 60.000000000000
nodata_value 0
-1.1 -1.1 1.1 1.1
-1.1 -1.1 1.1 1.1
-1.1 -1.1 1.1 1.1
-1.1 -1.1 1.1 1.1""")

ds = gdal.Open('/vsimem/in.asc')

data = ds.GetRasterBand(1).ReadRaster(buf_xsize=1,
buf_ysize=1,
resample_alg=gdal.GRIORA_Lanczos)
data_ar = struct.unpack('f' * 1, data)
expected_ar = (1.1754943508222875e-38, )
if data_ar != expected_ar:
gdaltest.post_reason('fail')
print(data_ar)
return 'fail'

data = ds.GetRasterBand(1).ReadRaster(buf_xsize=1,
buf_ysize=1,
resample_alg=gdal.GRIORA_Average)
data_ar = struct.unpack('f' * 1, data)
expected_ar = (1.1754943508222875e-38, )
if data_ar != expected_ar:
gdaltest.post_reason('fail')
print(data_ar)
return 'fail'

gdal.Unlink('/vsimem/in.asc')


gdal.FileFromMemBuffer('/vsimem/in.asc',
"""ncols 4
nrows 4
xllcorner 440720.000000000000
yllcorner 3750120.000000000000
cellsize 60.000000000000
nodata_value 0
-1 -1 1 1
-1 -1 1 1
-1 -1 1 1
-1 -1 1 1""")

ds = gdal.Open('/vsimem/in.asc')

data = ds.GetRasterBand(1).ReadRaster(buf_xsize=1,
buf_ysize=1,
resample_alg=gdal.GRIORA_Lanczos)
data_ar = struct.unpack('I' * 1, data)
expected_ar = (1, )
if data_ar != expected_ar:
gdaltest.post_reason('fail')
print(data_ar)
return 'fail'

data = ds.GetRasterBand(1).ReadRaster(buf_xsize=1,
buf_ysize=1,
resample_alg=gdal.GRIORA_Average)
data_ar = struct.unpack('I' * 1, data)
expected_ar = (1, )
if data_ar != expected_ar:
gdaltest.post_reason('fail')
print(data_ar)
return 'fail'

gdal.Unlink('/vsimem/in.asc')

return 'success'

gdaltest_list = [
rasterio_1,
rasterio_2,
Expand All @@ -1145,6 +1225,7 @@ def rasterio_lanczos_nodata():
rasterio_15,
rasterio_16,
rasterio_lanczos_nodata,
rasterio_resampled_value_is_nodata,
]

# gdaltest_list = [ rasterio_16 ]
Expand Down
4 changes: 2 additions & 2 deletions autotest/gdrivers/vrtpansharpen.py
View file
Open in desktop
Expand Up @@ -1681,7 +1681,7 @@ def vrtpansharpen_9():
gdaltest.post_reason('fail')
return 'fail'
cs = [vrt_ds.GetRasterBand(i + 1).Checksum() for i in range(vrt_ds.RasterCount)]
if cs not in([4179, 8767, 52257], [4175, 8758, 52249]):
if cs not in ([4640, 9158, 54450],):
gdaltest.post_reason('fail')
print(cs)
return 'fail'
Expand Down Expand Up @@ -1723,7 +1723,7 @@ def vrtpansharpen_9():
gdaltest.post_reason('fail')
return 'fail'
cs = [vrt_ds.GetRasterBand(i + 1).Checksum() for i in range(vrt_ds.RasterCount)]
if cs not in([4179, 8767, 52257], [4175, 8758, 52249]):
if cs not in ([4640, 9158, 54450],):
gdaltest.post_reason('fail')
print(cs)
return 'fail'
Expand Down
180 changes: 166 additions & 14 deletions gdal/gcore/overview.cpp
View file
Open in desktop
Expand Up @@ -249,6 +249,100 @@ static bool ReadColorTableAsArray( const GDALColorTable* poColorTable,
return true;
}

/************************************************************************/
/* GetReplacementValueIfNoData() */
/************************************************************************/

static float GetReplacementValueIfNoData(GDALDataType dt, int bHasNoData,
float fNoDataValue)
{
float fReplacementVal = 0.0f;
if( bHasNoData )
{
if( dt == GDT_Byte )
{
if( fNoDataValue == std::numeric_limits<unsigned char>::max() )
fReplacementVal = static_cast<float>(
std::numeric_limits<unsigned char>::max() - 1);
else
fReplacementVal = fNoDataValue + 1;
}
else if( dt == GDT_UInt16 )
{
if( fNoDataValue == std::numeric_limits<GUInt16>::max() )
fReplacementVal = static_cast<float>(
std::numeric_limits<GUInt16>::max() - 1);
else
fReplacementVal = fNoDataValue + 1;
}
else if( dt == GDT_Int16 )
{
if( fNoDataValue == std::numeric_limits<GInt16>::max() )
fReplacementVal = static_cast<float>(
std::numeric_limits<GInt16>::max() - 1);
else
fReplacementVal = fNoDataValue + 1;
}
else if( dt == GDT_UInt32 )
{
// Be careful to limited precision of float
fReplacementVal = fNoDataValue + 1;
double dfVal = fNoDataValue;
if( fReplacementVal >= static_cast<double>(std::numeric_limits<GUInt32>::max() - 128) )
{
while( fReplacementVal == fNoDataValue )
{
dfVal -= 1.0;
fReplacementVal = static_cast<float>(dfVal);
}
}
else
{
while( fReplacementVal == fNoDataValue )
{
dfVal += 1.0;
fReplacementVal = static_cast<float>(dfVal);
}
}
}
else if( dt == GDT_Int32 )
{
// Be careful to limited precision of float
fReplacementVal = fNoDataValue + 1;
double dfVal = fNoDataValue;
if( fReplacementVal >= static_cast<double>(std::numeric_limits<GInt32>::max() - 64) )
{
while( fReplacementVal == fNoDataValue )
{
dfVal -= 1.0;
fReplacementVal = static_cast<float>(dfVal);
}
}
else
{
while( fReplacementVal == fNoDataValue )
{
dfVal += 1.0;
fReplacementVal = static_cast<float>(dfVal);
}
}
}
else if( dt == GDT_Float32 || dt == GDT_Float64 )
{
if( fNoDataValue == 0 )
{
fReplacementVal = std::numeric_limits<float>::min();
}
else
{
fReplacementVal = static_cast<float>(
fNoDataValue + 1e-7 * fNoDataValue);
}
}
}
return fReplacementVal;
}

/************************************************************************/
/* GDALResampleChunk32R_Average() */
/************************************************************************/
Expand Down Expand Up @@ -284,6 +378,8 @@ GDALResampleChunk32R_AverageT( double dfXRatioDstToSrc,
tNoDataValue = 0;
else
tNoDataValue = static_cast<T>(fNoDataValue);
const T tReplacementVal = static_cast<T>(GetReplacementValueIfNoData(
poOverview->GetRasterDataType(), bHasNoData, fNoDataValue));

int nChunkRightXOff = nChunkXOff + nChunkXSize;
int nChunkBottomYOff = nChunkYOff + nChunkYSize;
Expand Down Expand Up @@ -441,7 +537,10 @@ GDALResampleChunk32R_AverageT( double dfXRatioDstToSrc,
+ pSrcScanlineShifted[nChunkXSize]
+ pSrcScanlineShifted[1+nChunkXSize];

pDstScanline[iDstPixel] = static_cast<T>((nTotal + 2) / 4);
auto nVal = static_cast<T>((nTotal + 2) / 4);
if( bHasNoData && nVal == tNoDataValue )
nVal = tReplacementVal;
pDstScanline[iDstPixel] = nVal;
pSrcScanlineShifted += 2;
}
}
Expand Down Expand Up @@ -480,11 +579,19 @@ GDALResampleChunk32R_AverageT( double dfXRatioDstToSrc,
}
else if( eWrkDataType == GDT_Byte ||
eWrkDataType == GDT_UInt16)
pDstScanline[iDstPixel] =
static_cast<T>((dfTotal + nCount / 2) / nCount);
{
auto nVal = static_cast<T>((dfTotal + nCount / 2) / nCount);
if( bHasNoData && nVal == tNoDataValue )
nVal = tReplacementVal;
pDstScanline[iDstPixel] = nVal;
}
else
pDstScanline[iDstPixel] =
static_cast<T>(dfTotal / nCount);
{
auto nVal = static_cast<T>(dfTotal / nCount);
if( bHasNoData && nVal == tNoDataValue )
nVal = tReplacementVal;
pDstScanline[iDstPixel] = nVal;
}
}
}
}
Expand Down Expand Up @@ -1772,6 +1879,33 @@ GDALResampleChunk32R_ConvolutionT( double dfXRatioDstToSrc,
{
if( !bHasNoData )
fNoDataValue = 0.0f;
const float fReplacementVal = GetReplacementValueIfNoData(
papoDstBands[0]->GetRasterDataType(), bHasNoData, fNoDataValue);
// cppcheck-suppress unreadVariable
const int isIntegerDT = GDALDataTypeIsInteger(papoDstBands[0]->GetRasterDataType());
const auto nNodataValueInt64 = static_cast<GInt64>(fNoDataValue);

auto replaceValIfNodata =
[bHasNoData, isIntegerDT, nNodataValueInt64, fNoDataValue, fReplacementVal](float fVal)
{
if( bHasNoData )
{
if( isIntegerDT )
{
if( nNodataValueInt64 == static_cast<GInt64>(fVal) )
{
// Do not use the nodata value
return fReplacementVal;
}
}
else if( fNoDataValue == fVal )
{
// Do not use the nodata value
return fReplacementVal;
}
}
return fVal;
};

/* -------------------------------------------------------------------- */
/* Allocate work buffers. */
Expand Down Expand Up @@ -2098,6 +2232,14 @@ GDALResampleChunk32R_ConvolutionT( double dfXRatioDstToSrc,
GDALResampleConvolutionVertical_16cols(
padfHorizontalFilteredBand + j, nDstXSize, padfWeights,
nSrcLineCount, pafDstScanline + iFilteredPixelOff );
if( bHasNoData )
{
for( int k = 0; k < 16; k++ )
{
pafDstScanline[iFilteredPixelOff + k] =
replaceValIfNodata(pafDstScanline[iFilteredPixelOff + k]);
}
}
}
#else
for( ;
Expand All @@ -2107,16 +2249,24 @@ GDALResampleChunk32R_ConvolutionT( double dfXRatioDstToSrc,
GDALResampleConvolutionVertical_8cols(
padfHorizontalFilteredBand + j, nDstXSize, padfWeights,
nSrcLineCount, pafDstScanline + iFilteredPixelOff );
if( bHasNoData )
{
for( int k = 0; k < 8; k++ )
{
pafDstScanline[iFilteredPixelOff + k] =
replaceValIfNodata(pafDstScanline[iFilteredPixelOff + k]);
}
}
}
#endif

for( ; iFilteredPixelOff < nDstXSize; iFilteredPixelOff++, j++ )
{
const double dfVal =
const float fVal = static_cast<float>(
GDALResampleConvolutionVertical(
padfHorizontalFilteredBand + j,
nDstXSize, padfWeights, nSrcLineCount );
pafDstScanline[iFilteredPixelOff] = static_cast<float>(dfVal);
nDstXSize, padfWeights, nSrcLineCount ));
pafDstScanline[iFilteredPixelOff] = replaceValIfNodata(fVal);
}
#else
for( ;
Expand All @@ -2128,17 +2278,19 @@ GDALResampleChunk32R_ConvolutionT( double dfXRatioDstToSrc,
GDALResampleConvolutionVertical_2cols(
padfHorizontalFilteredBand + j, nDstXSize, padfWeights,
nSrcLineCount, dfVal1, dfVal2 );
pafDstScanline[iFilteredPixelOff] = static_cast<float>(dfVal1);
pafDstScanline[iFilteredPixelOff+1] =
static_cast<float>(dfVal2);
pafDstScanline[iFilteredPixelOff] = replaceValIfNodata(
static_cast<float>(dfVal1));
pafDstScanline[iFilteredPixelOff+1] = replaceValIfNodata(
static_cast<float>(dfVal2));
}
if( iFilteredPixelOff < nDstXSize )
{
const double dfVal =
GDALResampleConvolutionVertical(
padfHorizontalFilteredBand + j,
nDstXSize, padfWeights, nSrcLineCount );
pafDstScanline[iFilteredPixelOff] = static_cast<float>(dfVal);
pafDstScanline[iFilteredPixelOff] = replaceValIfNodata(
static_cast<float>(dfVal));
}
#endif
}
Expand Down Expand Up @@ -2195,8 +2347,8 @@ GDALResampleChunk32R_ConvolutionT( double dfXRatioDstToSrc,
}
if( dfWeightSum > 0.0 )
{
pafDstScanline[iFilteredPixelOff] =
static_cast<float>(dfVal / dfWeightSum);
pafDstScanline[iFilteredPixelOff] = replaceValIfNodata(
static_cast<float>(dfVal / dfWeightSum));
}
else
{
Expand Down