Fix rendering artifacts on the edges of resampled raster tiles

We now buffer the request for the input to a resample filter by a variable number of pixels (depending on the resampling type) in order to fetch the neighbouring pixels to the edges of individual raster blocks. This allows the resampling to utilise these neighbouring pixels when resampling the edges of tiles, avoiding discrepancies and rendering artifacts over the borders of raster tiles. Fixes #30152
qgis · Nov 15, 2019 · b7ae145 · b7ae145
1 parent 73023e8
commit b7ae145
Show file tree

Hide file tree

Showing 11 changed files with 240 additions and 15 deletions.
diff --git a/python/core/auto_generated/raster/qgsbilinearrasterresampler.sip.in b/python/core/auto_generated/raster/qgsbilinearrasterresampler.sip.in
@@ -35,6 +35,8 @@ Constructor for QgsBilinearRasterResampler.
 
     virtual QgsBilinearRasterResampler *clone() const /Factory/;
 
+    virtual int tileBufferPixels() const;
+
 };
 
 /************************************************************************

diff --git a/python/core/auto_generated/raster/qgscubicrasterresampler.sip.in b/python/core/auto_generated/raster/qgscubicrasterresampler.sip.in
@@ -34,6 +34,8 @@ Constructor for QgsCubicRasterResampler.
 
     virtual QString type() const;
 
+    virtual int tileBufferPixels() const;
+
 };
 
 /************************************************************************

diff --git a/python/core/auto_generated/raster/qgsrasterresampler.sip.in b/python/core/auto_generated/raster/qgsrasterresampler.sip.in
@@ -54,6 +54,15 @@ Needs to be implemented by subclasses.
 Gets a deep copy of this object.
 Needs to be reimplemented by subclasses.
 Ownership is transferred to the caller.
+%End
+
+    virtual int tileBufferPixels() const;
+%Docstring
+Returns the optional tile buffer size in pixels. This represents
+the size to buffer individual resampled tile requests prior to resampling,
+in order to avoid rendering artifacts at the edges of raster tile boundaries.
+
+.. versionadded:: 3.10.1
 %End
 };
 

diff --git a/src/core/raster/qgsbilinearrasterresampler.cpp b/src/core/raster/qgsbilinearrasterresampler.cpp
@@ -25,6 +25,11 @@ QgsBilinearRasterResampler *QgsBilinearRasterResampler::clone() const
   return new QgsBilinearRasterResampler();
 }
 
+int QgsBilinearRasterResampler::tileBufferPixels() const
+{
+  return 1;
+}
+
 Q_NOWARN_DEPRECATED_PUSH
 void QgsBilinearRasterResampler::resample( const QImage &srcImage, QImage &dstImage )
 {

diff --git a/src/core/raster/qgsbilinearrasterresampler.h b/src/core/raster/qgsbilinearrasterresampler.h
@@ -43,6 +43,7 @@ class CORE_EXPORT QgsBilinearRasterResampler: public QgsRasterResamplerV2
     QImage resampleV2( const QImage &source, const QSize &size ) override;
     QString type() const override;
     QgsBilinearRasterResampler *clone() const override SIP_FACTORY;
+    int tileBufferPixels() const override;
 };
 
 #endif // QGSBILINEARRASTERRESAMPLER_H
diff --git a/src/core/raster/qgscubicrasterresampler.cpp b/src/core/raster/qgscubicrasterresampler.cpp
@@ -42,3 +42,8 @@ QString QgsCubicRasterResampler::type() const
   return QStringLiteral( "cubic" );
 }
 
+int QgsCubicRasterResampler::tileBufferPixels() const
+{
+  return 2;
+}
+
diff --git a/src/core/raster/qgscubicrasterresampler.h b/src/core/raster/qgscubicrasterresampler.h
@@ -42,6 +42,7 @@ class CORE_EXPORT QgsCubicRasterResampler: public QgsRasterResamplerV2
     QImage resampleV2( const QImage &source, const QSize &size ) override;
     Q_DECL_DEPRECATED void resample( const QImage &srcImage, QImage &dstImage ) override SIP_DEPRECATED;
     QString type() const override;
+    int tileBufferPixels() const override;
 };
 
 #endif // QGSCUBICRASTERRESAMPLER_H
diff --git a/src/core/raster/qgsrasterresamplefilter.cpp b/src/core/raster/qgsrasterresamplefilter.cpp
@@ -127,17 +127,18 @@ QgsRasterBlock *QgsRasterResampleFilter::block( int bandNo, QgsRectangle  const
     return outputBlock.release();
 
   double oversampling = 1.0; // approximate global oversampling factor
-
+  double outputXRes;
+  double providerXRes = 0;
   if ( mZoomedInResampler || mZoomedOutResampler )
   {
     QgsRasterDataProvider *provider = dynamic_cast<QgsRasterDataProvider *>( mInput->sourceInput() );
     if ( provider && ( provider->capabilities() & QgsRasterDataProvider::Size ) )
     {
-      double xRes = extent.width() / width;
-      double providerXRes = provider->extent().width() / provider->xSize();
-      double pixelRatio = xRes / providerXRes;
+      outputXRes = extent.width() / width;
+      providerXRes = provider->extent().width() / provider->xSize();
+      double pixelRatio = outputXRes / providerXRes;
       oversampling = ( pixelRatio > mMaxOversampling ) ? mMaxOversampling : pixelRatio;
-      QgsDebugMsgLevel( QStringLiteral( "xRes = %1 providerXRes = %2 pixelRatio = %3 oversampling = %4" ).arg( xRes ).arg( providerXRes ).arg( pixelRatio ).arg( oversampling ), 4 );
+      QgsDebugMsgLevel( QStringLiteral( "xRes = %1 providerXRes = %2 pixelRatio = %3 oversampling = %4" ).arg( outputXRes ).arg( providerXRes ).arg( pixelRatio ).arg( oversampling ), 4 );
     }
     else
     {
@@ -165,12 +166,31 @@ QgsRasterBlock *QgsRasterResampleFilter::block( int bandNo, QgsRectangle  const
   double oversamplingX = ( static_cast< double >( width ) * oversampling ) / width;
   double oversamplingY = ( static_cast< double >( height ) * oversampling ) / height;
 
-  // TODO: we must also increase the extent to get correct result on borders of parts
+  // we must also increase the extent to get correct result on borders of parts
+  int tileBufferPixels = 0;
+  if ( providerXRes != 0 )
+  {
+    if ( mZoomedInResampler && ( oversamplingX < 1.0 || qgsDoubleNear( oversampling, 1.0 ) ) )
+    {
+      tileBufferPixels = mZoomedInResampler->tileBufferPixels();
+    }
+    else if ( mZoomedOutResampler && oversamplingX > 1.0 )
+    {
+      tileBufferPixels = mZoomedOutResampler->tileBufferPixels();
+    }
+  }
+  const double sourceTileBufferSize = providerXRes * tileBufferPixels;
+
+  QgsRectangle bufferedExtent( extent.xMinimum() - sourceTileBufferSize,
+                               extent.yMinimum() - sourceTileBufferSize,
+                               extent.xMaximum() + sourceTileBufferSize,
+                               extent.yMaximum() + sourceTileBufferSize
+                             );
 
-  int resWidth = width * oversamplingX;
-  int resHeight = height * oversamplingY;
+  int resWidth = static_cast< int >( std::round( width * oversamplingX ) ) + 2 * tileBufferPixels;
+  int resHeight = static_cast< int >( std::round( height * oversamplingY ) ) + 2 * tileBufferPixels;
 
-  std::unique_ptr< QgsRasterBlock > inputBlock( mInput->block( bandNumber, extent, resWidth, resHeight, feedback ) );
+  std::unique_ptr< QgsRasterBlock > inputBlock( mInput->block( bandNumber, bufferedExtent, resWidth, resHeight, feedback ) );
   if ( !inputBlock || inputBlock->isEmpty() )
   {
     QgsDebugMsg( QStringLiteral( "No raster data!" ) );
@@ -184,21 +204,24 @@ QgsRasterBlock *QgsRasterResampleFilter::block( int bandNo, QgsRectangle  const
 
   //resample image
   QImage img = inputBlock->image();
-  QImage dstImg;
 
+  int resampleWidth = static_cast< int >( std::round( width * ( bufferedExtent.width() / extent.width() ) ) );
+  int resampleHeight = static_cast< int >( std::round( height * ( bufferedExtent.height() / extent.height() ) ) );
+
+  QImage dstImg;
   if ( mZoomedInResampler && ( oversamplingX < 1.0 || qgsDoubleNear( oversampling, 1.0 ) ) )
   {
     QgsDebugMsgLevel( QStringLiteral( "zoomed in resampling" ), 4 );
 
     if ( QgsRasterResamplerV2 *resamplerV2 = dynamic_cast< QgsRasterResamplerV2 * >( mZoomedInResampler.get( ) ) )
     {
-      dstImg = resamplerV2->resampleV2( img, QSize( width, height ) );
+      dstImg = resamplerV2->resampleV2( img, QSize( resampleWidth, resampleHeight ) );
     }
     else
     {
       // old inefficient interface
       Q_NOWARN_DEPRECATED_PUSH
-      QImage dstImg = QImage( width, height, QImage::Format_ARGB32_Premultiplied );
+      QImage dstImg = QImage( resampleWidth, resampleHeight, QImage::Format_ARGB32_Premultiplied );
       mZoomedInResampler->resample( img, dstImg );
       Q_NOWARN_DEPRECATED_POP
     }
@@ -209,13 +232,13 @@ QgsRasterBlock *QgsRasterResampleFilter::block( int bandNo, QgsRectangle  const
 
     if ( QgsRasterResamplerV2 *resamplerV2 = dynamic_cast< QgsRasterResamplerV2 * >( mZoomedOutResampler.get( ) ) )
     {
-      dstImg = resamplerV2->resampleV2( img, QSize( width, height ) );
+      dstImg = resamplerV2->resampleV2( img, QSize( resampleWidth, resampleHeight ) );
     }
     else
     {
       // old inefficient interface
       Q_NOWARN_DEPRECATED_PUSH
-      QImage dstImg = QImage( width, height, QImage::Format_ARGB32_Premultiplied );
+      QImage dstImg = QImage( resampleWidth, resampleHeight, QImage::Format_ARGB32_Premultiplied );
       mZoomedOutResampler->resample( img, dstImg );
       Q_NOWARN_DEPRECATED_POP
     }
@@ -227,7 +250,10 @@ QgsRasterBlock *QgsRasterResampleFilter::block( int bandNo, QgsRectangle  const
     dstImg = img.scaled( width, height );
   }
 
-  outputBlock->setImage( &dstImg );
+  // extract desired part of dstImage
+  QImage cropped = tileBufferPixels > 0 ? dstImg.copy( ( resampleWidth - width ) / 2, ( resampleHeight - height ) / 2, width, height )
+                   : dstImg; // otherwise implicity copy, nice and cheap
+  outputBlock->setImage( &cropped );
 
   return outputBlock.release(); // No resampling
 }

diff --git a/src/core/raster/qgsrasterresampler.h b/src/core/raster/qgsrasterresampler.h
@@ -74,6 +74,15 @@ class CORE_EXPORT QgsRasterResampler
      * Ownership is transferred to the caller.
      */
     virtual QgsRasterResampler *clone() const = 0 SIP_FACTORY;
+
+    /**
+     * Returns the optional tile buffer size in pixels. This represents
+     * the size to buffer individual resampled tile requests prior to resampling,
+     * in order to avoid rendering artifacts at the edges of raster tile boundaries.
+     *
+     * \since QGIS 3.10.1
+     */
+    virtual int tileBufferPixels() const { return 0; }
 };
 
 

diff --git a/tests/src/python/CMakeLists.txt b/tests/src/python/CMakeLists.txt
@@ -185,6 +185,7 @@ ADD_PYTHON_TEST(PyQgsRasterFileWriterTask test_qgsrasterfilewritertask.py)
 ADD_PYTHON_TEST(PyQgsRasterLayer test_qgsrasterlayer.py)
 ADD_PYTHON_TEST(PyQgsRasterColorRampShader test_qgsrastercolorrampshader.py)
 ADD_PYTHON_TEST(PyQgsRasterRange test_qgsrasterrange.py)
+ADD_PYTHON_TEST(PyQgsRasterResampler test_qgsrasterresampler.py)
 ADD_PYTHON_TEST(PyQgsRatioLockButton test_qgsratiolockbutton.py)
 ADD_PYTHON_TEST(PyQgsRectangle test_qgsrectangle.py)
 ADD_PYTHON_TEST(PyQgsReferencedGeometry test_qgsreferencedgeometry.py)

diff --git a/tests/src/python/test_qgsrasterresampler.py b/tests/src/python/test_qgsrasterresampler.py
@@ -0,0 +1,164 @@
+# -*- coding: utf-8 -*-
+"""QGIS Unit tests for QgsRasterResampler.
+
+.. note:: This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+"""
+from builtins import str
+
+__author__ = 'Nyall Dawson'
+__date__ = '14/11/2019'
+__copyright__ = 'Copyright 2019, The QGIS Project'
+
+import qgis  # NOQA
+
+import os
+
+from qgis.PyQt.QtGui import qRed
+
+from qgis.core import (QgsRasterLayer,
+                       QgsRectangle,
+                       QgsRasterResampleFilter,
+                       QgsSingleBandGrayRenderer,
+                       QgsCubicRasterResampler,
+                       QgsBilinearRasterResampler
+                       )
+from utilities import unitTestDataPath
+from qgis.testing import start_app, unittest
+
+# Convenience instances in case you may need them
+# not used in this test
+start_app()
+
+
+class TestQgsRasterResampler(unittest.TestCase):
+
+    def checkBlockContents(self, block, expected):
+        res = []
+        for r in range(block.height()):
+            res.append([qRed(block.color(r, c)) for c in range(block.width())])
+        self.assertEqual(res, expected)
+
+    def testBilinearResample(self):
+        path = os.path.join(unitTestDataPath(), 'landsat.tif')
+        raster_layer = QgsRasterLayer(path, 'test')
+        self.assertTrue(raster_layer.isValid())
+
+        extent = QgsRectangle(785994.37761193525511771,
+                              3346249.2209800467826426,
+                              786108.49096253968309611,
+                              3346362.94137834152206779)
+
+        renderer = QgsSingleBandGrayRenderer(raster_layer.dataProvider(), 1)
+        filter = QgsRasterResampleFilter(renderer)
+
+        # default (nearest neighbour) resampling
+        block = filter.block(1, extent, 2, 2)
+        self.checkBlockContents(block, [[124, 127], [125, 126]])
+
+        block = filter.block(1, extent, 4, 4)
+        self.checkBlockContents(block, [[124, 124, 127, 127],
+                                        [124, 124, 127, 127],
+                                        [125, 125, 126, 126],
+                                        [125, 125, 126, 126]]
+                                )
+
+        block = filter.block(1, extent, 8, 8)
+        self.checkBlockContents(block, [[124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126]])
+
+        # with resampling
+        filter.setZoomedInResampler(QgsBilinearRasterResampler())
+        block = filter.block(1, extent, 2, 2)
+        self.checkBlockContents(block, [[124, 127], [125, 126]])
+
+        block = filter.block(1, extent, 4, 4)
+        self.checkBlockContents(block,
+                                [[124, 124, 126, 126],
+                                 [124, 124, 125, 126],
+                                 [124, 124, 125, 126],
+                                 [125, 125, 125, 126]]
+                                )
+
+        block = filter.block(1, extent, 8, 8)
+        self.checkBlockContents(block,
+                                [[124, 124, 124, 125, 125, 126, 126, 126],
+                                 [124, 124, 124, 125, 125, 126, 126, 126],
+                                 [124, 124, 124, 124, 125, 125, 126, 126],
+                                 [124, 124, 124, 124, 125, 125, 126, 126],
+                                 [124, 124, 124, 124, 125, 125, 126, 126],
+                                 [124, 124, 124, 124, 125, 125, 126, 126],
+                                 [125, 125, 125, 125, 125, 125, 126, 126],
+                                 [125, 125, 125, 125, 125, 125, 126, 126]]
+                                )
+
+    def testCubicResample(self):
+        path = os.path.join(unitTestDataPath(), 'landsat.tif')
+        raster_layer = QgsRasterLayer(path, 'test')
+        self.assertTrue(raster_layer.isValid())
+
+        extent = QgsRectangle(785994.37761193525511771,
+                              3346249.2209800467826426,
+                              786108.49096253968309611,
+                              3346362.94137834152206779)
+
+        renderer = QgsSingleBandGrayRenderer(raster_layer.dataProvider(), 1)
+        filter = QgsRasterResampleFilter(renderer)
+
+        # default (nearest neighbour) resampling
+        block = filter.block(1, extent, 2, 2)
+        self.checkBlockContents(block, [[124, 127], [125, 126]])
+
+        block = filter.block(1, extent, 4, 4)
+        self.checkBlockContents(block, [[124, 124, 127, 127],
+                                        [124, 124, 127, 127],
+                                        [125, 125, 126, 126],
+                                        [125, 125, 126, 126]]
+                                )
+
+        block = filter.block(1, extent, 8, 8)
+        self.checkBlockContents(block, [[124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [124, 124, 124, 124, 127, 127, 127, 127],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126],
+                                        [125, 125, 125, 125, 126, 126, 126, 126]])
+
+        # with resampling
+        filter.setZoomedInResampler(QgsCubicRasterResampler())
+        block = filter.block(1, extent, 2, 2)
+        self.checkBlockContents(block, [[124, 127], [125, 126]])
+
+        block = filter.block(1, extent, 4, 4)
+        self.checkBlockContents(block,
+                                [[124, 125, 127, 127],
+                                 [124, 125, 126, 127],
+                                 [125, 125, 126, 126],
+                                 [125, 125, 126, 126]]
+                                )
+
+        block = filter.block(1, extent, 8, 8)
+        self.checkBlockContents(block,
+                                [[124, 124, 124, 125, 126, 127, 127, 127],
+                                 [124, 124, 124, 125, 126, 127, 127, 127],
+                                 [124, 124, 124, 125, 126, 127, 127, 127],
+                                 [125, 124, 124, 125, 126, 126, 127, 126],
+                                 [125, 125, 125, 125, 126, 126, 126, 126],
+                                 [125, 125, 125, 125, 126, 126, 126, 126],
+                                 [125, 125, 125, 125, 126, 126, 126, 126],
+                                 [125, 125, 125, 125, 126, 126, 126, 126]]
+                                )
+
+
+if __name__ == '__main__':
+    unittest.main()