sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-10-07 00:15:48 -04:00
Code Issues Packages Projects Releases Wiki Activity

Early resampling: replicate behaviour similar to QgsRasterResampleFilter when zooming out beyond the max oversampling factor

Browse Source
This commit is contained in:
Even Rouault 2020-06-12 18:50:50 +02:00 committed by Nyall Dawson
parent aaaf0f7dc4
commit df8dd11e40
2 changed files with 119 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

102
src/core/providers/gdal/qgsgdalprovider.cpp
View File

@ -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;

25
tests/src/python/test_qgsrasterresampler.py
View File

@ -429,7 +429,7 @@ class TestQgsRasterResampler(unittest.TestCase):
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)
@ -437,8 +437,27 @@ class TestQgsRasterResampler(unittest.TestCase):
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):