Skip to content

Commit df8dd11

Browse files
rouaultnyalldawson
authored andcommittedJun 19, 2020
Early resampling: replicate behaviour similar to QgsRasterResampleFilter when zooming out beyond the max oversampling factor
1 parent aaaf0f7 commit df8dd11

File tree

2 files changed

+119
-8
lines changed

2 files changed

+119
-8
lines changed
 

‎src/core/providers/gdal/qgsgdalprovider.cpp

Lines changed: 97 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -720,9 +720,6 @@ bool QgsGdalProvider::canDoResampling(
720720
int bufferWidthPix,
721721
int bufferHeightPix )
722722
{
723-
if ( !mProviderResamplingEnabled )
724-
return false;
725-
726723
QMutexLocker locker( mpMutex );
727724
if ( !initIfNeeded() )
728725
return false;
@@ -816,6 +813,7 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
816813
const double srcXRes = mGeoTransform[1];
817814
const double srcYRes = mGeoTransform[5]; // may be negative?
818815
QgsDebugMsgLevel( QStringLiteral( "reqXRes = %1 reqYRes = %2 srcXRes = %3 srcYRes = %4" ).arg( reqXRes ).arg( reqYRes ).arg( srcXRes ).arg( srcYRes ), 5 );
816+
const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / srcYRes );
819817

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

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

930-
const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / std::fabs( srcYRes ) );
931929
ResamplingMethod method;
932930
if ( resamplingFactor < 1 )
933931
{
@@ -972,6 +970,100 @@ bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int
972970
&sExtraArg ) == CE_None;
973971
}
974972
}
973+
// Provider resampling was asked but we cannot do it in a performant way
974+
// (too much downsampling compared to the allowed maximum resampling factor),
975+
// so fallback to something replicating QgsRasterResampleFilter behaviour
976+
else if ( mProviderResamplingEnabled &&
977+
mZoomedOutResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest &&
978+
resamplingFactor > 1 )
979+
{
980+
// Do the resampling in two steps:
981+
// - downsample with nearest neighbour down to tgtWidth * mMaxOversampling, tgtHeight * mMaxOversampling
982+
// - then downsample with mZoomedOutResamplingMethod down to tgtWidth, tgtHeight
983+
const int tgtWidth = tgtRightOri - tgtLeftOri + 1;
984+
const int tgtHeight = tgtBottomOri - tgtTopOri + 1;
985+
986+
const int tmpWidth = static_cast<int>( tgtWidth * mMaxOversampling + 0.5 );
987+
const int tmpHeight = static_cast<int>( tgtHeight * mMaxOversampling + 0.5 );
988+
989+
// Allocate temporary block
990+
size_t bufferSize = dataSize * static_cast<size_t>( tmpWidth ) * static_cast<size_t>( tmpHeight );
991+
#ifdef Q_PROCESSOR_X86_32
992+
// Safety check for 32 bit systems
993+
qint64 _buffer_size = dataSize * static_cast<qint64>( tmpWidth ) * static_cast<qint64>( tmpHeight );
994+
if ( _buffer_size != static_cast<qint64>( bufferSize ) )
995+
{
996+
QgsDebugMsg( QStringLiteral( "Integer overflow calculating buffer size on a 32 bit system." ) );
997+
return false;
998+
}
999+
#endif
1000+
char *tmpBlock = static_cast<char *>( qgsMalloc( bufferSize ) );
1001+
if ( ! tmpBlock )
1002+
{
1003+
QgsDebugMsgLevel( QStringLiteral( "Couldn't allocate temporary buffer of %1 bytes" ).arg( dataSize * tmpWidth * tmpHeight ), 5 );
1004+
return false;
1005+
}
1006+
CPLErrorReset();
1007+
1008+
1009+
CPLErr err = gdalRasterIO( gdalBand, GF_Read,
1010+
srcLeft, srcTop, srcWidth, srcHeight,
1011+
static_cast<void *>( tmpBlock ),
1012+
tmpWidth, tmpHeight, type,
1013+
0, 0, feedback );
1014+
1015+
if ( err != CPLE_None )
1016+
{
1017+
const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
1018+
if ( feedback )
1019+
feedback->appendError( lastError );
1020+
1021+
QgsLogger::warning( "RasterIO error: " + lastError );
1022+
qgsFree( tmpBlock );
1023+
return false;
1024+
}
1025+
1026+
GDALDriverH hDriverMem = GDALGetDriverByName( "MEM" );
1027+
if ( !hDriverMem )
1028+
{
1029+
qgsFree( tmpBlock );
1030+
return false;
1031+
}
1032+
gdal::dataset_unique_ptr hSrcDS( GDALCreate(
1033+
hDriverMem, "", tmpWidth, tmpHeight, 0, GDALGetRasterDataType( gdalBand ), nullptr ) );
1034+
1035+
char **papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
1036+
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( dataSize )
1037+
<< QStringLiteral( "LINEOFFSET=%1" ).arg( dataSize * tmpWidth )
1038+
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( tmpBlock ) ) );
1039+
GDALAddBand( hSrcDS.get(), GDALGetRasterDataType( gdalBand ), papszOptions );
1040+
CSLDestroy( papszOptions );
1041+
1042+
GDALRasterIOExtraArg sExtraArg;
1043+
INIT_RASTERIO_EXTRA_ARG( sExtraArg );
1044+
1045+
if ( mZoomedOutResamplingMethod == ResamplingMethod::Bilinear )
1046+
sExtraArg.eResampleAlg = GRIORA_Bilinear;
1047+
else if ( mZoomedOutResamplingMethod == ResamplingMethod::Cubic )
1048+
sExtraArg.eResampleAlg = GRIORA_Cubic;
1049+
else
1050+
sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
1051+
CPLErr eErr = GDALRasterIOEx( GDALGetRasterBand( hSrcDS.get(), 1 ),
1052+
GF_Read,
1053+
0, 0, tmpWidth, tmpHeight,
1054+
static_cast<char *>( data ) +
1055+
( tgtTopOri * bufferWidthPix + tgtLeftOri ) * dataSize,
1056+
tgtWidth,
1057+
tgtHeight,
1058+
type,
1059+
dataSize,
1060+
dataSize * bufferWidthPix,
1061+
&sExtraArg );
1062+
1063+
qgsFree( tmpBlock );
1064+
1065+
return eErr == CE_None;
1066+
}
9751067

9761068
const int tgtTop = tgtTopOri;
9771069
const int tgtBottom = tgtBottomOri;

‎tests/src/python/test_qgsrasterresampler.py

Lines changed: 22 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -429,16 +429,35 @@ def testGDALResampling_downsampling_cubic(self):
429429
block = provider.block(1, self._getExtentRequestInside(), 1, 1)
430430
self.checkRawBlockContents(block, [[86]])
431431

432-
def testGDALResampling_downsampling_bilinear_ignored(self):
432+
def testGDALResampling_downsampling_bilinear_beyond_max_oversampling_factor(self):
433433

434434
with self.setupGDALResampling() as provider:
435435
provider.setMaxOversampling(1.5)
436436
provider.setZoomedOutResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Bilinear)
437437
provider.setZoomedInResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Cubic) # ignored
438438
block = provider.block(1, self._getExtentRequestInside(), 1, 1)
439439
# as we request at a x2 oversampling factor and the limit is 1.5
440-
# fallback to nearest neighbour
441-
self.checkRawBlockContents(block, [[50]])
440+
# fallback to an alternate method using first nearest resampling
441+
# and then bilinear
442+
self.checkRawBlockContents(block, [[120]])
443+
444+
def testGDALResampling_downsampling_bilinear_beyond_max_oversampling_factor_containing_raster_extent(self):
445+
446+
xmin = 100
447+
ymax = 1000
448+
xres = 5
449+
yres = 5
450+
451+
extent = QgsRectangle(xmin - 10 * xres,
452+
ymax - (5 + 10) * yres,
453+
xmin + (5 + 10) * xres,
454+
ymax + 10 * yres)
455+
456+
with self.setupGDALResampling() as provider:
457+
provider.setZoomedInResamplingMethod(QgsRasterDataProvider.ResamplingMethod.Bilinear)
458+
provider.setMaxOversampling(1)
459+
block = provider.block(1, extent, 3, 3)
460+
self.checkRawBlockContents(block, [[0, 0, 0], [0, 50.0, 0], [0, 0, 0]])
442461

443462
def testGDALResampling_nominal_resolution_containing_raster_extent(self):
444463

0 commit comments

Comments
 (0)
Please sign in to comment.