Skip to content
Permalink
Browse files

Early resampling: replicate behaviour similar to QgsRasterResampleFil…

…ter when zooming out beyond the max oversampling factor
  • Loading branch information
rouault authored and nyalldawson committed Jun 12, 2020
1 parent aaaf0f7 commit df8dd11e40e61d542934b024c6ed677afdcd2667
Showing with 119 additions and 8 deletions.
  1. +97 −5 src/core/providers/gdal/qgsgdalprovider.cpp
  2. +22 −3 tests/src/python/test_qgsrasterresampler.py
@@ -720,9 +720,6 @@ bool QgsGdalProvider::canDoResampling(
int bufferWidthPix,
int bufferHeightPix )
{
if ( !mProviderResamplingEnabled )
return false;

QMutexLocker locker( mpMutex );
if ( !initIfNeeded() )
return false;
@@ -816,6 +813,7 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
const double srcXRes = mGeoTransform[1];
const double srcYRes = mGeoTransform[5]; // may be negative?
QgsDebugMsgLevel( QStringLiteral( "reqXRes = %1 reqYRes = %2 srcXRes = %3 srcYRes = %4" ).arg( reqXRes ).arg( reqYRes ).arg( srcXRes ).arg( srcYRes ), 5 );
const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / srcYRes );

GDALRasterBandH gdalBand = getBand( bandNo );
const GDALDataType type = static_cast<GDALDataType>( mGdalDataType.at( bandNo - 1 ) );
@@ -860,7 +858,8 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
const int srcHeight = srcBottom - srcTop + 1;

// Use GDAL resampling if asked and possible
if ( canDoResampling( bandNo, reqExtent, bufferWidthPix, bufferHeightPix ) )
if ( mProviderResamplingEnabled &&
canDoResampling( bandNo, reqExtent, bufferWidthPix, bufferHeightPix ) )
{
int tgtTop = tgtTopOri;
int tgtBottom = tgtBottomOri;
@@ -927,7 +926,6 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
GDALRasterIOExtraArg sExtraArg;
INIT_RASTERIO_EXTRA_ARG( sExtraArg );

const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / std::fabs( srcYRes ) );
ResamplingMethod method;
if ( resamplingFactor < 1 )
{
@@ -972,6 +970,100 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
&sExtraArg ) == CE_None;
}
}
// Provider resampling was asked but we cannot do it in a performant way
// (too much downsampling compared to the allowed maximum resampling factor),
// so fallback to something replicating QgsRasterResampleFilter behaviour
else if ( mProviderResamplingEnabled &&
mZoomedOutResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest &&
resamplingFactor > 1 )
{
// Do the resampling in two steps:
// - downsample with nearest neighbour down to tgtWidth * mMaxOversampling, tgtHeight * mMaxOversampling
// - then downsample with mZoomedOutResamplingMethod down to tgtWidth, tgtHeight
const int tgtWidth = tgtRightOri - tgtLeftOri + 1;
const int tgtHeight = tgtBottomOri - tgtTopOri + 1;

const int tmpWidth = static_cast<int>( tgtWidth * mMaxOversampling + 0.5 );
const int tmpHeight = static_cast<int>( tgtHeight * mMaxOversampling + 0.5 );

// Allocate temporary block
size_t bufferSize = dataSize * static_cast<size_t>( tmpWidth ) * static_cast<size_t>( tmpHeight );
#ifdef Q_PROCESSOR_X86_32
// Safety check for 32 bit systems
qint64 _buffer_size = dataSize * static_cast<qint64>( tmpWidth ) * static_cast<qint64>( tmpHeight );
if ( _buffer_size != static_cast<qint64>( bufferSize ) )
{
QgsDebugMsg( QStringLiteral( "Integer overflow calculating buffer size on a 32 bit system." ) );
return false;
}
#endif
char *tmpBlock = static_cast<char *>( qgsMalloc( bufferSize ) );
if ( ! tmpBlock )
{
QgsDebugMsgLevel( QStringLiteral( "Couldn't allocate temporary buffer of %1 bytes" ).arg( dataSize * tmpWidth * tmpHeight ), 5 );
return false;
}
CPLErrorReset();


CPLErr err = gdalRasterIO( gdalBand, GF_Read,
srcLeft, srcTop, srcWidth, srcHeight,
static_cast<void *>( tmpBlock ),
tmpWidth, tmpHeight, type,
0, 0, feedback );

if ( err != CPLE_None )
{
const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
if ( feedback )
feedback->appendError( lastError );

QgsLogger::warning( "RasterIO error: " + lastError );
qgsFree( tmpBlock );
return false;
}

GDALDriverH hDriverMem = GDALGetDriverByName( "MEM" );
if ( !hDriverMem )
{
qgsFree( tmpBlock );
return false;
}
gdal::dataset_unique_ptr hSrcDS( GDALCreate(
hDriverMem, "", tmpWidth, tmpHeight, 0, GDALGetRasterDataType( gdalBand ), nullptr ) );

char **papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( dataSize )
<< QStringLiteral( "LINEOFFSET=%1" ).arg( dataSize * tmpWidth )
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( tmpBlock ) ) );
GDALAddBand( hSrcDS.get(), GDALGetRasterDataType( gdalBand ), papszOptions );
CSLDestroy( papszOptions );

GDALRasterIOExtraArg sExtraArg;
INIT_RASTERIO_EXTRA_ARG( sExtraArg );

if ( mZoomedOutResamplingMethod == ResamplingMethod::Bilinear )
sExtraArg.eResampleAlg = GRIORA_Bilinear;
else if ( mZoomedOutResamplingMethod == ResamplingMethod::Cubic )
sExtraArg.eResampleAlg = GRIORA_Cubic;
else
sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
CPLErr eErr = GDALRasterIOEx( GDALGetRasterBand( hSrcDS.get(), 1 ),
GF_Read,
0, 0, tmpWidth, tmpHeight,
static_cast<char *>( data ) +
( tgtTopOri * bufferWidthPix + tgtLeftOri ) * dataSize,
tgtWidth,
tgtHeight,
type,
dataSize,
dataSize * bufferWidthPix,
&sExtraArg );

qgsFree( tmpBlock );

return eErr == CE_None;
}

const int tgtTop = tgtTopOri;
const int tgtBottom = tgtBottomOri;
@@ -429,16 +429,35 @@ def testGDALResampling_downsampling_cubic(self):
block = provider.block(1, self._getExtentRequestInside(), 1, 1)
self.checkRawBlockContents(block, [[86]])

def testGDALResampling_downsampling_bilinear_ignored(self):
def testGDALResampling_downsampling_bilinear_beyond_max_oversampling_factor(self):

with self.setupGDALResampling() as provider:
provider.setMaxOversampling(1.5)
provider.setZoomedOutResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Bilinear)
provider.setZoomedInResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Cubic) # ignored
block = provider.block(1, self._getExtentRequestInside(), 1, 1)
# as we request at a x2 oversampling factor and the limit is 1.5
# fallback to nearest neighbour
self.checkRawBlockContents(block, [[50]])
# fallback to an alternate method using first nearest resampling
# and then bilinear
self.checkRawBlockContents(block, [[120]])

def testGDALResampling_downsampling_bilinear_beyond_max_oversampling_factor_containing_raster_extent(self):

xmin = 100
ymax = 1000
xres = 5
yres = 5

extent = QgsRectangle(xmin - 10 * xres,
ymax - (5 + 10) * yres,
xmin + (5 + 10) * xres,
ymax + 10 * yres)

with self.setupGDALResampling() as provider:
provider.setZoomedInResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Bilinear)
provider.setMaxOversampling(1)
block = provider.block(1, extent, 3, 3)
self.checkRawBlockContents(block, [[0, 0, 0], [0, 50.0, 0], [0, 0, 0]])

def testGDALResampling_nominal_resolution_containing_raster_extent(self):

0 comments on commit df8dd11

Please sign in to comment.
You can’t perform that action at this time.