qgsrasterinterface.cpp

Go to the documentation of this file.

/***************************************************************************
    qgsrasterface.cpp - Internal raster processing modules interface
     --------------------------------------
    Date                 : Jun 21, 2012
    Copyright            : (C) 2012 by Radim Blazek
    email                : radim dot blazek at gmail dot com
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgsrasterinterface.h"
 
#include <limits>
#include <typeinfo>
 
#include "qgslogger.h"
#include "qgsrasterbandstats.h"
#include "qgsrasterhistogram.h"
#include "qgsrectangle.h"
 
#include <QByteArray>
#include <QString>
#include <QStringList>
#include <QTime>
 
#include "moc_qgsrasterinterface.cpp"
 
using namespace Qt::StringLiterals;
 
QgsRasterInterface::QgsRasterInterface( QgsRasterInterface *input )
  : mInput( input )
{}
 
Qgis::RasterInterfaceCapabilities QgsRasterInterface::capabilities() const
{
  return Qgis::RasterInterfaceCapability::NoCapabilities;
}
 
void QgsRasterInterface::initStatistics( QgsRasterBandStats &statistics, int bandNo, Qgis::RasterBandStatistics stats, const QgsRectangle &boundingBox, int sampleSize ) const
{
  QgsDebugMsgLevel( u"theBandNo = %1 sampleSize = %2"_s.arg( bandNo ).arg( sampleSize ), 4 );
 
  statistics.bandNumber = bandNo;
  statistics.statsGathered = stats;
 
  QgsRectangle finalExtent;
  if ( boundingBox.isEmpty() )
  {
    finalExtent = extent();
  }
  else
  {
    finalExtent = extent().intersect( boundingBox );
  }
  statistics.extent = finalExtent;
 
  if ( sampleSize > 0 )
  {
    // Calc resolution from theSampleSize
    double xRes, yRes;
    xRes = yRes = std::sqrt( ( finalExtent.width() * finalExtent.height() ) / sampleSize );
 
    // But limit by physical resolution
    if ( capabilities() & Qgis::RasterInterfaceCapability::Size )
    {
      const double srcXRes = extent().width() / xSize();
      const double srcYRes = extent().height() / ySize();
      if ( xRes < srcXRes )
        xRes = srcXRes;
      if ( yRes < srcYRes )
        yRes = srcYRes;
    }
    QgsDebugMsgLevel( u"xRes = %1 yRes = %2"_s.arg( xRes ).arg( yRes ), 4 );
 
    statistics.width = static_cast<int>( std::ceil( finalExtent.width() / xRes ) );
    statistics.height = static_cast<int>( std::ceil( finalExtent.height() / yRes ) );
  }
  else
  {
    if ( capabilities() & Qgis::RasterInterfaceCapability::Size )
    {
      statistics.width = xSize();
      statistics.height = ySize();
    }
    else
    {
      statistics.width = 1000;
      statistics.height = 1000;
    }
  }
  QgsDebugMsgLevel( u"theStatistics.width = %1 statistics.height = %2"_s.arg( statistics.width ).arg( statistics.height ), 4 );
}
 
bool QgsRasterInterface::hasStatistics( int bandNo, Qgis::RasterBandStatistics stats, const QgsRectangle &extent, int sampleSize )
{
  QgsDebugMsgLevel( u"theBandNo = %1 stats = %2 sampleSize = %3"_s.arg( bandNo ).arg( static_cast<int>( stats ) ).arg( sampleSize ), 4 );
  if ( mStatistics.isEmpty() )
    return false;
 
  QgsRasterBandStats myRasterBandStats;
  initStatistics( myRasterBandStats, bandNo, stats, extent, sampleSize );
 
  const auto constMStatistics = mStatistics;
  for ( const QgsRasterBandStats &stats : constMStatistics )
  {
    if ( stats.contains( myRasterBandStats ) )
    {
      QgsDebugMsgLevel( u"Has cached statistics."_s, 4 );
      return true;
    }
  }
  return false;
}
 
QgsRasterBandStats QgsRasterInterface::bandStatistics( int bandNo, Qgis::RasterBandStatistics stats, const QgsRectangle &extent, int sampleSize, QgsRasterBlockFeedback *feedback )
{
  QgsDebugMsgLevel( u"theBandNo = %1 stats = %2 sampleSize = %3"_s.arg( bandNo ).arg( static_cast<int>( stats ) ).arg( sampleSize ), 4 );
 
  // TODO: null values set on raster layer!!!
 
  QgsRasterBandStats myRasterBandStats;
  initStatistics( myRasterBandStats, bandNo, stats, extent, sampleSize );
 
  const auto constMStatistics = mStatistics;
  for ( const QgsRasterBandStats &stats : constMStatistics )
  {
    if ( stats.contains( myRasterBandStats ) )
    {
      QgsDebugMsgLevel( u"Using cached statistics."_s, 4 );
      return stats;
    }
  }
 
  const QgsRectangle myExtent = myRasterBandStats.extent;
  const int myWidth = myRasterBandStats.width;
  const int myHeight = myRasterBandStats.height;
 
  //int myDataType = dataType( bandNo );
 
  int myXBlockSize = xBlockSize();
  int myYBlockSize = yBlockSize();
  if ( myXBlockSize == 0 ) // should not happen, but happens
  {
    myXBlockSize = 500;
  }
  if ( myYBlockSize == 0 ) // should not happen, but happens
  {
    myYBlockSize = 500;
  }
 
  const int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
  const int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
 
  const double myXRes = myExtent.width() / myWidth;
  const double myYRes = myExtent.height() / myHeight;
  // TODO: progress signals
 
  // used by single pass stdev
  double myMean = 0;
  double mySumOfSquares = 0;
 
  bool myFirstIterationFlag = true;
  bool isNoData = false;
  for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
  {
    for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
    {
      if ( feedback && feedback->isCanceled() )
        return myRasterBandStats;
 
      QgsDebugMsgLevel( u"myYBlock = %1 myXBlock = %2"_s.arg( myYBlock ).arg( myXBlock ), 4 );
      const int myBlockWidth = std::min( myXBlockSize, myWidth - myXBlock * myXBlockSize );
      const int myBlockHeight = std::min( myYBlockSize, myHeight - myYBlock * myYBlockSize );
 
      const double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
      const double xmax = xmin + myBlockWidth * myXRes;
      const double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
      const double ymax = ymin - myBlockHeight * myYRes;
 
      const QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
 
      std::unique_ptr< QgsRasterBlock > blk( block( bandNo, myPartExtent, myBlockWidth, myBlockHeight, feedback ) );
 
      // Collect the histogram counts.
      for ( qgssize i = 0; i < ( static_cast< qgssize >( myBlockHeight ) ) * myBlockWidth; i++ )
      {
        const double myValue = blk->valueAndNoData( i, isNoData );
        if ( isNoData )
          continue; // NULL
 
        myRasterBandStats.sum += myValue;
        myRasterBandStats.elementCount++;
 
        if ( !std::isfinite( myValue ) )
          continue; // inf
 
        if ( myFirstIterationFlag )
        {
          myFirstIterationFlag = false;
          myRasterBandStats.minimumValue = myValue;
          myRasterBandStats.maximumValue = myValue;
        }
        else
        {
          if ( myValue < myRasterBandStats.minimumValue )
          {
            myRasterBandStats.minimumValue = myValue;
          }
          if ( myValue > myRasterBandStats.maximumValue )
          {
            myRasterBandStats.maximumValue = myValue;
          }
        }
 
        // Single pass stdev
        const double myDelta = myValue - myMean;
        myMean += myDelta / myRasterBandStats.elementCount;
        mySumOfSquares += myDelta * ( myValue - myMean );
      }
    }
  }
 
  myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
  myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;
 
  myRasterBandStats.sumOfSquares = mySumOfSquares; // OK with single pass?
 
  // stdDev may differ  from GDAL stats, because GDAL is using naive single pass
  // algorithm which is more error prone (because of rounding errors)
  // Divide result by sample size - 1 and get square root to get stdev
  myRasterBandStats.stdDev = std::sqrt( mySumOfSquares / ( static_cast< double >( myRasterBandStats.elementCount ) - 1 ) );
 
  QgsDebugMsgLevel( u"************ STATS **************"_s, 4 );
  QgsDebugMsgLevel( u"MIN %1"_s.arg( myRasterBandStats.minimumValue ), 4 );
  QgsDebugMsgLevel( u"MAX %1"_s.arg( myRasterBandStats.maximumValue ), 4 );
  QgsDebugMsgLevel( u"RANGE %1"_s.arg( myRasterBandStats.range ), 4 );
  QgsDebugMsgLevel( u"MEAN %1"_s.arg( myRasterBandStats.mean ), 4 );
  QgsDebugMsgLevel( u"STDDEV %1"_s.arg( myRasterBandStats.stdDev ), 4 );
 
  myRasterBandStats.statsGathered = Qgis::RasterBandStatistic::All;
  mStatistics.append( myRasterBandStats );
 
  return myRasterBandStats;
}
 
bool QgsRasterInterface::hasStatistics( int bandNo, int stats, const QgsRectangle &extent, int sampleSize )
{
  return hasStatistics( bandNo, static_cast< Qgis::RasterBandStatistics >( stats ), extent, sampleSize );
}
 
void QgsRasterInterface::initHistogram( QgsRasterHistogram &histogram, int bandNo, int binCount, double minimum, double maximum, const QgsRectangle &boundingBox, int sampleSize, bool includeOutOfRange )
{
  histogram.bandNumber = bandNo;
  histogram.minimum = minimum;
  histogram.maximum = maximum;
  histogram.includeOutOfRange = includeOutOfRange;
 
  const Qgis::DataType mySrcDataType = sourceDataType( bandNo );
 
  if ( std::isnan( histogram.minimum ) )
  {
    // TODO: this was OK when stats/histogram were calced in provider,
    // but what TODO in other interfaces? Check for mInput for now.
    if ( !mInput && mySrcDataType == Qgis::DataType::Byte )
    {
      histogram.minimum = 0; // see histogram() for shift for rounding
    }
    else
    {
      // We need statistics -> avoid histogramDefaults in hasHistogram if possible
      // TODO: use approximated statistics if approximated histogram is requested
      // (theSampleSize > 0)
      const QgsRasterBandStats stats = bandStatistics( bandNo, Qgis::RasterBandStatistic::Min, boundingBox, sampleSize );
      histogram.minimum = stats.minimumValue;
    }
  }
  if ( std::isnan( histogram.maximum ) )
  {
    if ( !mInput && mySrcDataType == Qgis::DataType::Byte )
    {
      histogram.maximum = 255;
    }
    else
    {
      const QgsRasterBandStats stats = bandStatistics( bandNo, Qgis::RasterBandStatistic::Max, boundingBox, sampleSize );
      histogram.maximum = stats.maximumValue;
    }
  }
 
  QgsRectangle finalExtent;
  if ( boundingBox.isEmpty() )
  {
    finalExtent = extent();
  }
  else
  {
    finalExtent = extent().intersect( boundingBox );
  }
  histogram.extent = finalExtent;
 
  if ( sampleSize > 0 )
  {
    // Calc resolution from theSampleSize
    double xRes, yRes;
    xRes = yRes = std::sqrt( ( static_cast<double>( finalExtent.width() ) * finalExtent.height() ) / sampleSize );
 
    // But limit by physical resolution
    if ( capabilities() & Qgis::RasterInterfaceCapability::Size )
    {
      const double srcXRes = extent().width() / xSize();
      const double srcYRes = extent().height() / ySize();
      if ( xRes < srcXRes )
        xRes = srcXRes;
      if ( yRes < srcYRes )
        yRes = srcYRes;
    }
    QgsDebugMsgLevel( u"xRes = %1 yRes = %2"_s.arg( xRes ).arg( yRes ), 4 );
 
    histogram.width = static_cast<int>( finalExtent.width() / xRes );
    histogram.height = static_cast<int>( finalExtent.height() / yRes );
  }
  else
  {
    if ( capabilities() & Qgis::RasterInterfaceCapability::Size )
    {
      histogram.width = xSize();
      histogram.height = ySize();
    }
    else
    {
      histogram.width = 1000;
      histogram.height = 1000;
    }
  }
  QgsDebugMsgLevel( u"theHistogram.width = %1 histogram.height = %2"_s.arg( histogram.width ).arg( histogram.height ), 4 );
 
  qint64 myBinCount = binCount;
  if ( myBinCount == 0 )
  {
    // TODO: this was OK when stats/histogram were calced in provider,
    // but what TODO in other interfaces? Check for mInput for now.
    if ( !mInput && mySrcDataType == Qgis::DataType::Byte )
    {
      myBinCount = 256; // Cannot store more values in byte
    }
    else
    {
      // There is no best default value, to display something reasonable in histogram chart,
      // binCount should be small, OTOH, to get precise data for cumulative cut, the number should be big.
      // Because it is easier to define fixed lower value for the chart, we calc optimum binCount
      // for higher resolution (to avoid calculating that where histogram() is used.
 
      // for Int16/Int32 make sure bin count <= actual range, because there is no sense in having
      // bins at fractional values
      if ( !mInput && ( mySrcDataType == Qgis::DataType::Int16 || mySrcDataType == Qgis::DataType::Int32 || mySrcDataType == Qgis::DataType::UInt16 || mySrcDataType == Qgis::DataType::UInt32 ) )
      {
        myBinCount = std::min( static_cast<qint64>( histogram.width ) * histogram.height, static_cast<qint64>( std::ceil( histogram.maximum - histogram.minimum + 1 ) ) );
      }
      else
      {
        // This is for not integer types
        const QgsRasterBandStats stats = bandStatistics( bandNo, Qgis::RasterBandStatistic::StdDev | Qgis::RasterBandStatistic::Range, boundingBox, sampleSize );
        // If stats are not available or stdDev is 0, use width*height as binCount
        // but this could lead to binCount being too large
        if ( stats.statsGathered.testFlags( Qgis::RasterBandStatistic::StdDev | Qgis::RasterBandStatistic::Range ) )
        {
          // Use Scott's Rule to determine bin count: binCount = 3.49 * stdDev / n^(1/3)
          // In standard C++ (prior to C++26), std::sqrt is not a constexpr function,
          // so we precompute the constant value  24 * sqrt( M_PI ) = 42.538892421732385
          constexpr double SCOTTS_COEFFICIENT = 42.538892421732385;
          myBinCount = static_cast<qint64>( stats.range / ( stats.stdDev * std::cbrt( SCOTTS_COEFFICIENT / ( sampleSize > 0 ? sampleSize : static_cast<double>( histogram.width ) * histogram.height ) ) ) );
        }
        else
        {
          // Fallback using the whole histogram size
          myBinCount = static_cast<qint64>( histogram.width ) * static_cast<qint64>( histogram.height );
        }
      }
    }
  }
  // Hard limit 10'000'000
  histogram.binCount = static_cast<int>( std::min( 10000000LL, myBinCount ) );
  QgsDebugMsgLevel( u"theHistogram.binCount = %1"_s.arg( histogram.binCount ), 4 );
}
 
void QgsRasterInterface::initStatistics( QgsRasterBandStats &statistics, int bandNo, int stats, const QgsRectangle &boundingBox, int binCount ) const
{
  initStatistics( statistics, bandNo, static_cast< Qgis::RasterBandStatistics>( stats ), boundingBox, binCount );
}
 
bool QgsRasterInterface::hasHistogram( int bandNo, int binCount, double minimum, double maximum, const QgsRectangle &extent, int sampleSize, bool includeOutOfRange )
{
  QgsDebugMsgLevel( u"theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5"_s.arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 4 );
  // histogramDefaults() needs statistics if minimum or maximum is NaN ->
  // do other checks which don't need statistics before histogramDefaults()
  if ( mHistograms.isEmpty() )
    return false;
 
  QgsRasterHistogram myHistogram;
  initHistogram( myHistogram, bandNo, binCount, minimum, maximum, extent, sampleSize, includeOutOfRange );
 
  const auto constMHistograms = mHistograms;
  for ( const QgsRasterHistogram &histogram : constMHistograms )
  {
    if ( histogram == myHistogram )
    {
      QgsDebugMsgLevel( u"Has cached histogram."_s, 4 );
      return true;
    }
  }
  return false;
}
 
QgsRasterHistogram QgsRasterInterface::histogram( int bandNo, int binCount, double minimum, double maximum, const QgsRectangle &extent, int sampleSize, bool includeOutOfRange, QgsRasterBlockFeedback *feedback )
{
  QgsDebugMsgLevel( u"theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5"_s.arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 4 );
 
  QgsRasterHistogram myHistogram;
  initHistogram( myHistogram, bandNo, binCount, minimum, maximum, extent, sampleSize, includeOutOfRange );
 
  // Find cached
  const auto constMHistograms = mHistograms;
  for ( const QgsRasterHistogram &histogram : constMHistograms )
  {
    if ( histogram == myHistogram )
    {
      QgsDebugMsgLevel( u"Using cached histogram."_s, 4 );
      return histogram;
    }
  }
 
  const int myBinCount = myHistogram.binCount;
  const int myWidth = myHistogram.width;
  const int myHeight = myHistogram.height;
  const QgsRectangle myExtent = myHistogram.extent;
  myHistogram.histogramVector.resize( myBinCount );
 
  int myXBlockSize = xBlockSize();
  int myYBlockSize = yBlockSize();
  if ( myXBlockSize == 0 ) // should not happen, but happens
  {
    myXBlockSize = 500;
  }
  if ( myYBlockSize == 0 ) // should not happen, but happens
  {
    myYBlockSize = 500;
  }
 
  const int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
  const int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
 
  const double myXRes = myExtent.width() / myWidth;
  const double myYRes = myExtent.height() / myHeight;
 
  double myMinimum = myHistogram.minimum;
  double myMaximum = myHistogram.maximum;
 
  // To avoid rounding errors
  // TODO: check this
  const double myerval = ( myMaximum - myMinimum ) / myHistogram.binCount;
  myMinimum -= 0.1 * myerval;
  myMaximum += 0.1 * myerval;
 
  QgsDebugMsgLevel( u"binCount = %1 myMinimum = %2 myMaximum = %3"_s.arg( myHistogram.binCount ).arg( myMinimum ).arg( myMaximum ), 4 );
 
  const double myBinSize = ( myMaximum - myMinimum ) / myBinCount;
 
  // TODO: progress signals
  bool isNoData = false;
  for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
  {
    for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
    {
      if ( feedback && feedback->isCanceled() )
        return myHistogram;
 
      const int myBlockWidth = std::min( myXBlockSize, myWidth - myXBlock * myXBlockSize );
      const int myBlockHeight = std::min( myYBlockSize, myHeight - myYBlock * myYBlockSize );
 
      const double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
      const double xmax = xmin + myBlockWidth * myXRes;
      const double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
      const double ymax = ymin - myBlockHeight * myYRes;
 
      const QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
 
      std::unique_ptr< QgsRasterBlock > blk( block( bandNo, myPartExtent, myBlockWidth, myBlockHeight, feedback ) );
 
      // Collect the histogram counts.
      for ( qgssize i = 0; i < ( static_cast< qgssize >( myBlockHeight ) ) * myBlockWidth; i++ )
      {
        const double myValue = blk->valueAndNoData( i, isNoData );
        if ( isNoData )
        {
          continue; // NULL
        }
 
        int myBinIndex = static_cast<int>( std::floor( ( myValue - myMinimum ) / myBinSize ) );
 
        if ( ( myBinIndex < 0 || myBinIndex > ( myBinCount - 1 ) ) && !includeOutOfRange )
        {
          continue;
        }
        if ( myBinIndex < 0 )
          myBinIndex = 0;
        if ( myBinIndex > ( myBinCount - 1 ) )
          myBinIndex = myBinCount - 1;
 
        myHistogram.histogramVector[myBinIndex] += 1;
        myHistogram.nonNullCount++;
      }
    }
  }
 
  myHistogram.valid = true;
  mHistograms.append( myHistogram );
 
#ifdef QGISDEBUG
  QString hist;
  for ( std::size_t i = 0; i < std::min< std::size_t >( myHistogram.histogramVector.size(), 500 ); i++ )
  {
    hist += QString::number( myHistogram.histogramVector.value( i ) ) + ' ';
  }
  QgsDebugMsgLevel( u"Histogram (max first 500 bins): "_s + hist, 4 );
#endif
 
  return myHistogram;
}
 
void QgsRasterInterface::cumulativeCut( int bandNo, double lowerCount, double upperCount, double &lowerValue, double &upperValue, const QgsRectangle &extent, int sampleSize )
{
  QgsDebugMsgLevel( u"theBandNo = %1 lowerCount = %2 upperCount = %3 sampleSize = %4"_s.arg( bandNo ).arg( lowerCount ).arg( upperCount ).arg( sampleSize ), 4 );
 
  const Qgis::DataType mySrcDataType = sourceDataType( bandNo );
 
  // Init to NaN is better than histogram min/max to catch errors
  lowerValue = std::numeric_limits<double>::quiet_NaN();
  upperValue = std::numeric_limits<double>::quiet_NaN();
 
  //get band stats to specify real histogram min/max (fix #9793 Byte bands)
  const QgsRasterBandStats stats = bandStatistics( bandNo, Qgis::RasterBandStatistic::Min, extent, sampleSize );
  if ( stats.maximumValue < stats.minimumValue )
    return;
 
  // for byte bands make sure bin count == actual range
  const int myBinCount = ( mySrcDataType == Qgis::DataType::Byte ) ? int( std::ceil( stats.maximumValue - stats.minimumValue + 1 ) ) : 0;
  const QgsRasterHistogram myHistogram = histogram( bandNo, myBinCount, stats.minimumValue, stats.maximumValue, extent, sampleSize );
  //QgsRasterHistogram myHistogram = histogram( bandNo, 0, std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), extent, sampleSize );
 
  const double myBinXStep = ( myHistogram.maximum - myHistogram.minimum ) / myHistogram.binCount;
  int myCount = 0;
  const int myMinCount = static_cast< int >( std::round( lowerCount * myHistogram.nonNullCount ) );
  const int myMaxCount = static_cast< int >( std::round( upperCount * myHistogram.nonNullCount ) );
  bool myLowerFound = false;
  QgsDebugMsgLevel( u"binCount = %1 minimum = %2 maximum = %3 myBinXStep = %4"_s.arg( myHistogram.binCount ).arg( myHistogram.minimum ).arg( myHistogram.maximum ).arg( myBinXStep ), 4 );
  QgsDebugMsgLevel( u"myMinCount = %1 myMaxCount = %2"_s.arg( myMinCount ).arg( myMaxCount ), 4 );
 
  for ( int myBin = 0; myBin < myHistogram.histogramVector.size(); myBin++ )
  {
    const int myBinValue = myHistogram.histogramVector.value( myBin );
    myCount += myBinValue;
    if ( !myLowerFound && myCount > myMinCount )
    {
      lowerValue = myHistogram.minimum + myBin * myBinXStep;
      myLowerFound = true;
      QgsDebugMsgLevel( u"found lowerValue %1 at bin %2"_s.arg( lowerValue ).arg( myBin ), 4 );
    }
    if ( myCount >= myMaxCount )
    {
      upperValue = myHistogram.minimum + myBin * myBinXStep;
      QgsDebugMsgLevel( u"found upperValue %1 at bin %2"_s.arg( upperValue ).arg( myBin ), 4 );
      break;
    }
  }
 
  // fix integer data - round down/up
  if ( mySrcDataType == Qgis::DataType::Byte
       || mySrcDataType == Qgis::DataType::Int16
       || mySrcDataType == Qgis::DataType::Int32
       || mySrcDataType == Qgis::DataType::UInt16
       || mySrcDataType == Qgis::DataType::UInt32 )
  {
    if ( !std::isnan( lowerValue ) )
      lowerValue = std::floor( lowerValue );
    if ( !std::isnan( upperValue ) )
      upperValue = std::ceil( upperValue );
  }
}
 
QString QgsRasterInterface::capabilitiesString() const
{
  QStringList abilitiesList;
 
  const Qgis::RasterInterfaceCapabilities abilities = capabilities();
 
  // Not all all capabilities are here (Size, IdentifyValue, IdentifyText,
  // IdentifyHtml, IdentifyFeature) because those are quite technical and probably
  // would be confusing for users
 
  if ( abilities & Qgis::RasterInterfaceCapability::Identify )
  {
    abilitiesList += tr( "Identify" );
  }
 
  if ( abilities & Qgis::RasterInterfaceCapability::BuildPyramids )
  {
    abilitiesList += tr( "Build Pyramids" );
  }
 
  QgsDebugMsgLevel( "Capability: " + abilitiesList.join( ", "_L1 ), 4 );
 
  return abilitiesList.join( ", "_L1 );
}
 
QString QgsRasterInterface::generateBandName( int bandNumber ) const
{
  if ( mInput )
    return mInput->generateBandName( bandNumber );
 
  // For bad layers bandCount is 0, no log!
  return tr( "Band" ) + u" %1"_s.arg( bandNumber, 1 + ( bandCount() > 0 ? static_cast< int >( std::log10( static_cast< double >( bandCount() ) ) ) : 0 ), 10, QChar( '0' ) );
}
 
QString QgsRasterInterface::colorInterpretationName( int bandNo ) const
{
  if ( mInput )
    return mInput->colorInterpretationName( bandNo );
 
  return QString();
}
 
QString QgsRasterInterface::displayBandName( int bandNumber ) const
{
  QString name = generateBandName( bandNumber );
  const QString colorInterp = colorInterpretationName( bandNumber );
  if ( colorInterp != tr( "Undefined" ) )
  {
    name.append( u" (%1)"_s.arg( colorInterp ) );
  }
  return name;
}
 
QgsRasterBandStats QgsRasterInterface::bandStatistics( int bandNo, int stats, const QgsRectangle &extent, int sampleSize, QgsRasterBlockFeedback *feedback )
{
  return bandStatistics( bandNo, static_cast< Qgis::RasterBandStatistics>( stats ), extent, sampleSize, feedback );
}
 
QgsRenderContext QgsRasterBlockFeedback::renderContext() const
{
  return mRenderContext;
}
 
void QgsRasterBlockFeedback::setRenderContext( const QgsRenderContext &renderContext )
{
  mRenderContext = renderContext;
}