Add method to resample QImage using gdal resample to a gdal memory da…

…taset
qgis · Oct 31, 2019 · c1326ce · c1326ce
1 parent 2308c51
commit c1326ce
Show file tree

Hide file tree

Showing 3 changed files with 110 additions and 0 deletions.
diff --git a/src/core/qgsgdalutils.cpp b/src/core/qgsgdalutils.cpp
@@ -14,12 +14,14 @@
  ***************************************************************************/
 
 #include "qgsgdalutils.h"
+#include "qgslogger.h"
 
 #define CPL_SUPRESS_CPLUSPLUS  //#spellok
 #include "gdal.h"
 #include "cpl_string.h"
 
 #include <QString>
+#include <QImage>
 
 bool QgsGdalUtils::supportsRasterCreate( GDALDriverH driver )
 {
@@ -184,6 +186,51 @@ 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 )
+{
+  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;
+}
+
 QString QgsGdalUtils::helpCreationOptionsFormat( QString format )
 {
   QString message;

diff --git a/src/core/qgsgdalutils.h b/src/core/qgsgdalutils.h
@@ -74,6 +74,13 @@ class CORE_EXPORT QgsGdalUtils
      */
     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 );
+
     /**
      * Gets creation options metadata for a given format
      * \since QGIS 3.10

diff --git a/tests/src/core/testqgsgdalutils.cpp b/tests/src/core/testqgsgdalutils.cpp
@@ -40,6 +40,7 @@ class TestQgsGdalUtils: public QObject
     void testCreateSingleBandTiffDataset();
     void testResampleSingleBandRaster();
     void testImageToDataset();
+    void testResampleImage();
 
   private:
 
@@ -250,6 +251,61 @@ void TestQgsGdalUtils::testImageToDataset()
   QCOMPARE( identify( dstDS.get(), 4, 200, 200 ), 255.0 );
 }
 
+void TestQgsGdalUtils::testResampleImage()
+{
+  QString inputFilename = QString( TEST_DATA_DIR ) + "/rgb256x256.png";
+  QImage src = QImage( inputFilename );
+  src = src.convertToFormat( QImage::Format_ARGB32 );
+  QVERIFY( !src.isNull() );
+
+  gdal::dataset_unique_ptr dstDS = QgsGdalUtils::resampleImage( QImage(), QgsRectangle( 0, 0, 256, 256 ), QSize( 50, 50 ), GRA_NearestNeighbour );
+  QVERIFY( !dstDS );
+
+  dstDS = QgsGdalUtils::resampleImage( src, QgsRectangle( 0, 0, 256, 256 ), QSize( 50, 50 ), GRA_NearestNeighbour );
+  QVERIFY( dstDS );
+
+  QCOMPARE( GDALGetRasterCount( dstDS.get() ), 4 );
+  QCOMPARE( GDALGetRasterXSize( dstDS.get() ), 50 );
+  QCOMPARE( GDALGetRasterYSize( dstDS.get() ), 50 );
+
+  QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 1 ) ), GDT_Float32 );
+  QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 2 ) ), GDT_Float32 );
+  QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 3 ) ), GDT_Float32 );
+  QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 4 ) ), GDT_Float32 );
+
+  QCOMPARE( identify( dstDS.get(), 1, 10, 10 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 1, 40, 10 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 1, 10, 40 ), 0.0 );
+  QCOMPARE( identify( dstDS.get(), 1, 40, 40 ), 0.0 );
+
+  QCOMPARE( identify( dstDS.get(), 2, 10, 10 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 2, 40, 10 ), 0.0 );
+  QCOMPARE( identify( dstDS.get(), 2, 10, 40 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 2, 40, 40 ), 0.0 );
+
+  QCOMPARE( identify( dstDS.get(), 3, 10, 10 ), 0.0 );
+  QCOMPARE( identify( dstDS.get(), 3, 40, 10 ), 0.0 );
+  QCOMPARE( identify( dstDS.get(), 3, 10, 40 ), 0.0 );
+  QCOMPARE( identify( dstDS.get(), 3, 40, 40 ), 255.0 );
+
+  QCOMPARE( identify( dstDS.get(), 4, 10, 10 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 4, 40, 10 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 4, 10, 40 ), 255.0 );
+  QCOMPARE( identify( dstDS.get(), 4, 40, 40 ), 255.0 );
+
+  dstDS = QgsGdalUtils::resampleImage( src, QgsRectangle( 0, 0, 256, 256 ), QSize( 4, 4 ), GRA_Cubic );
+  QVERIFY( dstDS );
+
+  QCOMPARE( GDALGetRasterCount( dstDS.get() ), 4 );
+  QCOMPARE( GDALGetRasterXSize( dstDS.get() ), 4 );
+  QCOMPARE( GDALGetRasterYSize( dstDS.get() ), 4 );
+
+  QGSCOMPARENEAR( identify( dstDS.get(), 1, 0, 0 ), 257.901092529, 0.0001 );
+  QGSCOMPARENEAR( identify( dstDS.get(), 1, 0, 1 ), 238.678848267, 0.0001 );
+  QGSCOMPARENEAR( identify( dstDS.get(), 1, 0, 2 ), 62.8939933777, 0.00001 );
+  QGSCOMPARENEAR( identify( dstDS.get(), 1, 0, 3 ), 87.553146, 0.00001 );
+}
+
 double TestQgsGdalUtils::identify( GDALDatasetH dataset, int band, int px, int py )
 {
   GDALRasterBandH hBand = GDALGetRasterBand( dataset, band );