Skip to content

Commit

Permalink
Heavily optimise GDAL based image resampling
Browse files Browse the repository at this point in the history
  • Loading branch information
@nyalldawson
nyalldawson committed Nov 4, 2019
1 parent 36072a6 commit 209130c
Show file tree
Hide file tree
Showing 5 changed files with 82 additions and 236 deletions.
194 changes: 52 additions & 142 deletions src/core/qgsgdalutils.cpp
View file
Expand Up @@ -97,61 +97,46 @@ gdal::dataset_unique_ptr QgsGdalUtils::createSingleBandTiffDataset( QString file
return hDstDS;
}

gdal::dataset_unique_ptr QgsGdalUtils::imageToMemoryDataset( const QImage &image, QgsRectangle extent, const QgsCoordinateReferenceSystem &crs )
gdal::dataset_unique_ptr QgsGdalUtils::imageToMemoryDataset( const QImage &image )
{
if ( image.isNull() )
return nullptr;
QByteArray red, green, blue, alpha;

const QRgb *rgb = reinterpret_cast<const QRgb *>( image.constBits() );

// potentially an overflow here -- that's on Qt!
const qgssize count = static_cast< qgssize >( image.width() ) * image.height();
red.resize( sizeof( float ) * count );
green.resize( sizeof( float ) * count );
blue.resize( sizeof( float ) * count );
alpha.resize( sizeof( float ) * count );

float *rData = reinterpret_cast<float *>( red.data() );
float *gData = reinterpret_cast<float *>( green.data() );
float *bData = reinterpret_cast<float *>( blue.data() );
float *aData = reinterpret_cast<float *>( alpha.data() );

gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createMultiBandMemoryDataset( GDT_Float32, 4, extent, image.width(), image.height(), crs ) );
for ( qgssize i = 0; i < count; ++i )
GDALDriverH hDriverMem = GDALGetDriverByName( "MEM" );
if ( !hDriverMem )
{
const QRgb c = *rgb++;
*rData++ = qRed( c );
*gData++ = qGreen( c );
*bData++ = qBlue( c );
*aData++ = qAlpha( c );
return nullptr;
}
gdal::dataset_unique_ptr hSrcDS( GDALCreate( hDriverMem, "", image.width(), image.height(), 0, GDT_Byte, nullptr ) );

char **papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( sizeof( QRgb ) )
<< QStringLiteral( "LINEOFFSET=%1" ).arg( image.bytesPerLine() )
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( rgb ) + 2 ) );
GDALAddBand( hSrcDS.get(), GDT_Byte, papszOptions );
CSLDestroy( papszOptions );

CPLErr err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 1 ), GF_Write, 0, 0, image.width(), image.height(), red.data(), image.width(), image.height(), GDT_Float32, 0, 0 );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to write red data to GDAL dataset" ) );
return nullptr;
}
err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 2 ), GF_Write, 0, 0, image.width(), image.height(), green.data(), image.width(), image.height(), GDT_Float32, 0, 0 );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to write green data to GDAL dataset" ) );
return nullptr;
}
err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 3 ), GF_Write, 0, 0, image.width(), image.height(), blue.data(), image.width(), image.height(), GDT_Float32, 0, 0 );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to write blue data to GDAL dataset" ) );
return nullptr;
}
err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 4 ), GF_Write, 0, 0, image.width(), image.height(), alpha.data(), image.width(), image.height(), GDT_Float32, 0, 0 );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to write alpha data to GDAL dataset" ) );
return nullptr;
}
papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( sizeof( QRgb ) )
<< QStringLiteral( "LINEOFFSET=%1" ).arg( image.bytesPerLine() )
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( rgb ) + 1 ) );
GDALAddBand( hSrcDS.get(), GDT_Byte, papszOptions );
CSLDestroy( papszOptions );

papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( sizeof( QRgb ) )
<< QStringLiteral( "LINEOFFSET=%1" ).arg( image.bytesPerLine() )
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( rgb ) ) );
GDALAddBand( hSrcDS.get(), GDT_Byte, papszOptions );
CSLDestroy( papszOptions );

papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
<< QStringLiteral( "PIXELOFFSET=%1" ).arg( sizeof( QRgb ) )
<< QStringLiteral( "LINEOFFSET=%1" ).arg( image.bytesPerLine() )
<< QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( rgb ) + 3 ) );
GDALAddBand( hSrcDS.get(), GDT_Byte, papszOptions );
CSLDestroy( papszOptions );

return hSrcDS;
}
Expand Down Expand Up @@ -186,133 +171,58 @@ void QgsGdalUtils::resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH h
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
}

gdal::dataset_unique_ptr QgsGdalUtils::resampleImage( const QImage &image, const QgsRectangle &extent, QSize outputSize, GDALResampleAlg resampleAlg )
QImage QgsGdalUtils::resampleImage( const QImage &image, QSize outputSize, GDALRIOResampleAlg resampleAlg )
{
gdal::dataset_unique_ptr srcDS = QgsGdalUtils::imageToMemoryDataset( image, extent, QgsCoordinateReferenceSystem() );
if ( !srcDS )
return nullptr;

gdal::dataset_unique_ptr destDS = QgsGdalUtils::createMultiBandMemoryDataset( GDT_Float32, 4, extent, outputSize.width(), outputSize.height(), QgsCoordinateReferenceSystem() );
if ( !destDS )
return nullptr;

gdal::warp_options_unique_ptr psWarpOptions( GDALCreateWarpOptions() );
psWarpOptions->hSrcDS = srcDS.get();
psWarpOptions->hDstDS = destDS.get();

psWarpOptions->nBandCount = 4;
psWarpOptions->panSrcBands = ( int * ) CPLMalloc( sizeof( int ) * 4 );
psWarpOptions->panDstBands = ( int * ) CPLMalloc( sizeof( int ) * 4 );
psWarpOptions->panSrcBands[0] = 1;
psWarpOptions->panSrcBands[1] = 2;
psWarpOptions->panSrcBands[2] = 3;
psWarpOptions->panSrcBands[3] = 4;
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->panDstBands[1] = 2;
psWarpOptions->panDstBands[2] = 3;
psWarpOptions->panDstBands[3] = 4;

psWarpOptions->eResampleAlg = resampleAlg;

// Establish reprojection transformer.
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( srcDS.get(), GDALGetProjectionRef( srcDS.get() ),
destDS.get(), GDALGetProjectionRef( destDS.get() ),
FALSE, 0.0, 1 );
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

// Initialize and execute the warp operation.
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions.get() );

oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( destDS.get() ), GDALGetRasterYSize( destDS.get() ) );

GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
return destDS;
}

QImage QgsGdalUtils::resampleImage( const QImage &image, QSize outputSize, GDALResampleAlg resampleAlg )
{
gdal::dataset_unique_ptr srcDS = resampleImage( image, QgsRectangle( 0, 0, image.width(), image.height() ), outputSize, resampleAlg );
gdal::dataset_unique_ptr srcDS = QgsGdalUtils::imageToMemoryDataset( image );
if ( !srcDS )
return QImage();

QImage res( outputSize, image.format() );
QRgb *rgb = reinterpret_cast<QRgb *>( res.bits() );
GDALRasterIOExtraArg extra;
INIT_RASTERIO_EXTRA_ARG( extra );
extra.eResampleAlg = resampleAlg;

GDALRasterBandH rBand = GDALGetRasterBand( srcDS.get(), 1 );
GDALRasterBandH gBand = GDALGetRasterBand( srcDS.get(), 2 );
GDALRasterBandH bBand = GDALGetRasterBand( srcDS.get(), 3 );
GDALRasterBandH aBand = GDALGetRasterBand( srcDS.get(), 4 );
QImage res( outputSize, QImage::Format_ARGB32 ); // premultiplied??
if ( res.isNull() )
return QImage();

const qgssize count = static_cast< qgssize >( outputSize.width() ) * outputSize.height();
GByte *rgb = reinterpret_cast<GByte *>( res.bits() );

float *rData = ( float * ) CPLMalloc( count * sizeof( float ) );
CPLErr err = GDALRasterIO( rBand, GF_Read, 0, 0, outputSize.width(), outputSize.height(),
rData, outputSize.width(), outputSize.height(), GDT_Float32, 0, 0 );
CPLErr err = GDALRasterIOEx( GDALGetRasterBand( srcDS.get(), 1 ), GF_Read, 0, 0, image.width(), image.height(), rgb + 2, outputSize.width(),
outputSize.height(), GDT_Byte, sizeof( QRgb ), res.bytesPerLine(), &extra );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to read red band" ) );
CPLFree( rData );
return QImage();
}
float *gData = ( float * ) CPLMalloc( count * sizeof( float ) );
err = GDALRasterIO( gBand, GF_Read, 0, 0, outputSize.width(), outputSize.height(),
gData, outputSize.width(), outputSize.height(), GDT_Float32, 0, 0 );

err = GDALRasterIOEx( GDALGetRasterBand( srcDS.get(), 2 ), GF_Read, 0, 0, image.width(), image.height(), rgb + 1, outputSize.width(),
outputSize.height(), GDT_Byte, sizeof( QRgb ), res.bytesPerLine(), &extra );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to read green band" ) );
CPLFree( rData );
CPLFree( gData );
return QImage();
}
float *bData = ( float * ) CPLMalloc( count * sizeof( float ) );
err = GDALRasterIO( bBand, GF_Read, 0, 0, outputSize.width(), outputSize.height(),
bData, outputSize.width(), outputSize.height(), GDT_Float32, 0, 0 );

err = GDALRasterIOEx( GDALGetRasterBand( srcDS.get(), 3 ), GF_Read, 0, 0, image.width(), image.height(), rgb, outputSize.width(),
outputSize.height(), GDT_Byte, sizeof( QRgb ), res.bytesPerLine(), &extra );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to read blue band" ) );
CPLFree( rData );
CPLFree( gData );
CPLFree( bData );
return QImage();
}
float *aData = ( float * ) CPLMalloc( count * sizeof( float ) );
err = GDALRasterIO( aBand, GF_Read, 0, 0, outputSize.width(), outputSize.height(),
aData, outputSize.width(), outputSize.height(), GDT_Float32, 0, 0 );

err = GDALRasterIOEx( GDALGetRasterBand( srcDS.get(), 4 ), GF_Read, 0, 0, image.width(), image.height(), rgb + 3, outputSize.width(),
outputSize.height(), GDT_Byte, sizeof( QRgb ), res.bytesPerLine(), &extra );
if ( err != CE_None )
{
QgsDebugMsg( QStringLiteral( "failed to read alpha band" ) );
CPLFree( rData );
CPLFree( gData );
CPLFree( bData );
CPLFree( aData );
return QImage();
}

float *rIn = rData;
float *gIn = gData;
float *bIn = bData;
float *aIn = aData;
for ( qgssize i = 0; i < count; ++i )
{
int red = std::max( 0, std::min( 255, static_cast< int >( std::round( *rIn++ ) ) ) );
int green = std::max( 0, std::min( 255, static_cast< int >( std::round( *gIn++ ) ) ) );
int blue = std::max( 0, std::min( 255, static_cast< int >( std::round( *bIn++ ) ) ) );
int alpha = std::max( 0, std::min( 255, static_cast< int >( std::round( *aIn++ ) ) ) );

*rgb++ = qRgba( red, green, blue, alpha );
}

CPLFree( rData );
CPLFree( gData );
CPLFree( bData );
CPLFree( aData );

return res;
}

QString QgsGdalUtils::helpCreationOptionsFormat( QString format )
QString QgsGdalUtils::helpCreationOptionsFormat( const QString &format )
{
QString message;
GDALDriverH myGdalDriver = GDALGetDriverByName( format.toLocal8Bit().constData() );
Expand Down
29 changes: 13 additions & 16 deletions src/core/qgsgdalutils.h
View file
Expand Up @@ -60,38 +60,24 @@ class CORE_EXPORT QgsGdalUtils
*/
static gdal::dataset_unique_ptr createSingleBandTiffDataset( QString filename, GDALDataType dataType, QgsRectangle extent, int width, int height, const QgsCoordinateReferenceSystem &crs );

/**
* Converts an \a image to a GDAL memory dataset.
*
* \since QGIS 3.12
*/
static gdal::dataset_unique_ptr imageToMemoryDataset( const QImage &image, QgsRectangle extent, const QgsCoordinateReferenceSystem &crs );

/**
* Resamples a single band raster to the destination dataset with different resolution (and possibly with different CRS).
* Ideally the source dataset should cover the whole area or the destination dataset.
* \since QGIS 3.8
*/
static void resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg );

/**
* Resamples a QImage \a image to a GDAL memory dataset with different resolution.
* \since QGIS 3.12
*/
static gdal::dataset_unique_ptr resampleImage( const QImage &image, const QgsRectangle &extent,
QSize outputSize, GDALResampleAlg resampleAlg );

/**
* Resamples a QImage \a image using GDAL resampler.
* \since QGIS 3.12
*/
static QImage resampleImage( const QImage &image, QSize outputSize, GDALResampleAlg resampleAlg );
static QImage resampleImage( const QImage &image, QSize outputSize, GDALRIOResampleAlg resampleAlg );

/**
* Gets creation options metadata for a given format
* \since QGIS 3.10
*/
static QString helpCreationOptionsFormat( QString format );
static QString helpCreationOptionsFormat( const QString &format );

/**
* Validates creation options for a given format, regardless of layer.
Expand All @@ -104,6 +90,17 @@ class CORE_EXPORT QgsGdalUtils
* \since QGIS 3.10
*/
static char **papszFromStringList( const QStringList &list );

private:

/**
* Converts an \a image to a GDAL memory dataset.
*
* \since QGIS 3.12
*/
static gdal::dataset_unique_ptr imageToMemoryDataset( const QImage &image );

friend class TestQgsGdalUtils;
};

#endif // QGSGDALUTILS_H
2 changes: 1 addition & 1 deletion src/core/raster/qgsbilinearrasterresampler.cpp
View file
Expand Up @@ -28,7 +28,7 @@ QgsBilinearRasterResampler *QgsBilinearRasterResampler::clone() const
Q_NOWARN_DEPRECATED_PUSH
void QgsBilinearRasterResampler::resample( const QImage &srcImage, QImage &dstImage )
{
dstImage = QgsGdalUtils::resampleImage( srcImage, dstImage.size(), GRA_Bilinear );
dstImage = QgsGdalUtils::resampleImage( srcImage, dstImage.size(), GRIORA_Bilinear );
}
Q_NOWARN_DEPRECATED_POP

Expand Down
10 changes: 8 additions & 2 deletions src/core/raster/qgscubicrasterresampler.cpp
View file
Expand Up @@ -27,12 +27,18 @@ QgsCubicRasterResampler *QgsCubicRasterResampler::clone() const

QImage QgsCubicRasterResampler::resampleV2( const QImage &source, const QSize &size )
{
return QgsGdalUtils::resampleImage( source, size, GRA_Cubic );
return QgsGdalUtils::resampleImage( source, size, GRIORA_Cubic );
}

Q_NOWARN_DEPRECATED_PUSH
void QgsCubicRasterResampler::resample( const QImage &srcImage, QImage &dstImage )
{
dstImage = QgsGdalUtils::resampleImage( srcImage, dstImage.size(), GRA_Cubic );
dstImage = QgsGdalUtils::resampleImage( srcImage, dstImage.size(), GRIORA_Cubic );
}
Q_NOWARN_DEPRECATED_POP

QString QgsCubicRasterResampler::type() const
{
return QStringLiteral( "cubic" );
}

0 comments on commit 209130c

Please sign in to comment.