new raster histogram

qgis · Jul 26, 2012 · d794522 · d794522
1 parent 47cc4f2
commit d794522
Show file tree

Hide file tree

Showing 4 changed files with 271 additions and 1 deletion.
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
@@ -390,6 +390,7 @@ SET(QGIS_CORE_HDRS
   raster/qgspseudocolorshader.h
   raster/qgsrasterpyramid.h
   raster/qgsrasterbandstats.h
+  raster/qgsrasterhistogram.h
   raster/qgsrasterinterface.h
   raster/qgsrasterlayer.h
   raster/qgsrastertransparency.h

diff --git a/src/core/qgsrasterdataprovider.cpp b/src/core/qgsrasterdataprovider.cpp
@@ -23,6 +23,8 @@
 #include <QMap>
 #include <QByteArray>
 
+#include <cmath>
+
 void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle
                                        const & viewExtent, int width,
                                        int height,
@@ -270,7 +272,6 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
   return ba;
 }
 
-
 QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
 {
   double myNoDataValue = noDataValue();
@@ -426,6 +427,169 @@ QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
   return myRasterBandStats;
 }
 
+QgsRasterHistogram QgsRasterDataProvider::histogram( int theBandNo,
+    double theMinimum, double theMaximum,
+    int theBinCount,
+    const QgsRectangle & theExtent,
+    int theSampleSize,
+    bool theIncludeOutOfRange )
+{
+  QgsRasterHistogram myHistogram;
+  myHistogram.bandNumber = theBandNo;
+  myHistogram.minimum = theMinimum;
+  myHistogram.maximum = theMaximum;
+  myHistogram.includeOutOfRange = theIncludeOutOfRange;
+  //myHistogram.sampleSize = theSampleSize
+
+  // First calc defaults
+  QgsRectangle myExtent = theExtent.isEmpty() ? extent() : theExtent;
+  myHistogram.extent = myExtent;
+
+  int myWidth, myHeight;
+  if ( theSampleSize > 0 )
+  {
+    // Calc resolution from theSampleSize
+    double xRes, yRes;
+    xRes = yRes = sqrt(( myExtent.width() * myExtent.height() ) / theSampleSize );
+
+    // But limit by physical resolution
+    if ( capabilities() & Size )
+    {
+      double srcXRes = extent().width() / xSize();
+      double srcYRes = extent().height() / ySize();
+      if ( xRes < srcXRes ) xRes = srcXRes;
+      if ( yRes < srcYRes ) yRes = srcYRes;
+    }
+
+    myWidth = static_cast <int>( myExtent.width() / xRes );
+    myHeight = static_cast <int>( myExtent.height() / yRes );
+  }
+  else
+  {
+    if ( capabilities() & Size )
+    {
+      myWidth = xSize();
+      myHeight = ySize();
+    }
+    else
+    {
+      myWidth = 1000;
+      myHeight = 1000;
+    }
+  }
+  myHistogram.width = myWidth;
+  myHistogram.height = myHeight;
+  QgsDebugMsg( QString( "myWidth = %1 myHeight = %2" ).arg( myWidth ).arg( myHeight ) );
+
+  double myNoDataValue = noDataValue();
+  int myDataType = dataType( theBandNo );
+
+  int myBinCount = theBinCount;
+  if ( myBinCount == 0 )
+  {
+    if ( myDataType == QgsRasterDataProvider::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. In any any case, it does not make sense to use more than width*height;
+      myBinCount = myWidth * myHeight;
+      if ( myBinCount > 1000 )  myBinCount = 1000;
+    }
+  }
+  myHistogram.binCount = theBinCount;
+  QgsDebugMsg( QString( "myBinCount = %1" ).arg( myBinCount ) );
+
+  // Check if we have cached
+  foreach( QgsRasterHistogram histogram, mHistograms )
+  {
+    if ( histogram.bandNumber == theBandNo &&
+         histogram.minimum == theMinimum &&
+         histogram.maximum == theMaximum &&
+         histogram.binCount == myBinCount &&
+         histogram.extent == myExtent &&
+         histogram.width == myWidth &&
+         histogram.height == myHeight &&
+         histogram.includeOutOfRange == theIncludeOutOfRange )
+    {
+      return histogram;
+    }
+  }
+
+  myHistogram.histogramVector.resize( myBinCount );
+
+  int  myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
+  myXBlockSize = xBlockSize();
+  myYBlockSize = yBlockSize();
+
+  if ( myXBlockSize == 0 || myYBlockSize == 0 ) // should not happen
+  {
+    return myHistogram;
+  }
+
+  myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
+  myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
+
+  void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) );
+
+  int myBandXSize = xSize();
+  int myBandYSize = ySize();
+  double myXRes = myExtent.width() / myWidth;
+  double myYRes = myExtent.height() / myHeight;
+  double binSize = ( theMaximum - theMinimum ) / myBinCount;
+  for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
+  {
+    for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
+    {
+      int myPartWidth = qMin( myXBlockSize, myWidth - iXBlock * myXBlockSize );
+      int myPartHeight = qMin( myYBlockSize, myHeight - iYBlock * myYBlockSize );
+
+      double xmin = myExtent.xMinimum() + iXBlock * myXBlockSize * myXRes;
+      double xmax = xmin + myPartWidth * myXRes;
+      double ymin = myExtent.yMaximum() - iYBlock * myYBlockSize * myYRes;
+      double ymax = ymin - myPartHeight * myYRes;
+
+      QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
+
+      readBlock( theBandNo, myPartExtent, myPartWidth, myPartHeight, myData );
+
+      // Collect the histogram counts.
+      for ( int iY = 0; iY < myPartHeight; iY++ )
+      {
+        for ( int iX = 0; iX < myPartWidth; iX++ )
+        {
+          double myValue = readValue( myData, myDataType, iX + ( iY * myPartWidth ) );
+          //QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) );
+
+          if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
+          {
+            continue; // NULL
+          }
+
+          int myBinIndex = static_cast <int>( floor(( myValue - theMinimum ) /  binSize ) ) ;
+          if (( myBinIndex < 0 || myBinIndex > ( myBinCount - 1 ) ) && !theIncludeOutOfRange )
+          {
+            continue;
+          }
+          if ( myBinIndex < 0 ) myBinIndex = 0;
+          if ( myBinIndex > ( myBinCount - 1 ) ) myBinIndex = myBinCount - 1;
+
+          myHistogram.histogramVector[myBinIndex] += 1;
+          myHistogram.nonNullCount++;
+        }
+      }
+    } //end of column wise loop
+  } //end of row wise loop
+
+  CPLFree( myData );
+
+  myHistogram.valid = true;
+  mHistograms.append( myHistogram );
+
+  return myHistogram;
+}
+
 double QgsRasterDataProvider::readValue( void *data, int type, int index )
 {
   if ( !data )

diff --git a/src/core/qgsrasterdataprovider.h b/src/core/qgsrasterdataprovider.h
@@ -30,6 +30,7 @@
 #include "qgsrasterpyramid.h"
 #include "qgscoordinatereferencesystem.h"
 #include "qgsrasterbandstats.h"
+#include "qgsrasterhistogram.h"
 
 #include "cpl_conv.h"
 #include <cmath>
@@ -314,6 +315,22 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
       Q_UNUSED( theIgnoreOutOfRangeFlag ); Q_UNUSED( theThoroughBandScanFlag );
     }
 
+    /** \brief Get histogram. Histograms are cached in providers.
+     * @param theBandNo The band (number).
+     * @param theMinimum Minimum value.
+     * @param theMaximum Maximum value.
+     * @param theBinCount Number of bins (intervals,buckets). If 0, the number of bins is decided automaticaly according to data type, raster size etc.
+     * @param theExtent Extent used to calc histogram, if empty, whole raster extent is used.
+     * @param theSampleSize Approximate number of cells in sample. If 0, all cells (whole raster will be used). If raster does not have exact size (WCS without exact size for example), provider decides size of sample.
+     * @return Vector of non NULL cell counts for each bin.
+     */
+    virtual QgsRasterHistogram histogram( int theBandNo,
+                                          double theMinimum, double theMaximum,
+                                          int theBinCount = 0,
+                                          const QgsRectangle & theExtent = QgsRectangle(),
+                                          int theSampleSize = 0,
+                                          bool theIncludeOutOfRange = false );
+
     /** \brief Create pyramid overviews */
     virtual QString buildPyramids( const QList<QgsRasterPyramid>  & thePyramidList,
                                    const QString &  theResamplingMethod = "NEAREST",
@@ -503,6 +520,9 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
     bool mValidNoDataValue;
 
     QgsRectangle mExtent;
+
+    /** \brief List  of cached histograms, all bands mixed */
+    QList <QgsRasterHistogram> mHistograms;
 };
 
 #endif
diff --git a/src/core/raster/qgsrasterhistogram.h b/src/core/raster/qgsrasterhistogram.h
@@ -0,0 +1,85 @@
+/***************************************************************************
+                        qgsrasterhistogram.h
+                              -------------------
+ begin                : July 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSRASTERHISTOGRAM
+#define QGSRASTERHISTOGRAM
+
+#include <QString>
+#include <QVector>
+
+#include <limits>
+
+/** \ingroup core
+ * The QgsRasterHistogram is a container for histogram of a single raster band.
+ * It is used to cache computed histograms in raster providers.
+ */
+class CORE_EXPORT QgsRasterHistogram 
+{
+  public:
+    typedef QVector<int> HistogramVector;
+
+    QgsRasterHistogram()
+    {
+      bandNumber = 0;
+      binCount = 0;
+      nonNullCount = 0;
+      //sampleSize = 0;
+      includeOutOfRange = false;
+      maximum = 0;
+      minimum = 0;
+      width = 0;
+      height = 0;
+      valid = false;
+    }
+
+    /** \brief The gdal band number (starts at 1)*/
+    int bandNumber;
+
+    /** \brief Number of bins (intervals,buckets) in histogram. */
+
+    /** \brief The number of non NULL cells used to calculate histogram. */
+    int nonNullCount;
+
+    /** \brief Approximate number of cells used to calc histogram. Approximately 
+     *  width * height. */
+    //int sampleSize;
+
+    /** \brief Whether histogram includes out of range values (in first and last bin) */
+    bool includeOutOfRange;
+
+    /** \brief Store the histogram for a given layer */
+    HistogramVector histogramVector;
+
+    /** \brief The maximum histogram value. */
+    double maximum;
+
+    /** \brief The minimum histogram value. */
+    double minimum;
+
+    /** \brief Number of columns used to calc histogram */
+    int width;
+
+    /** \brief Number of rows used to calc histogram */
+    int height;
+
+    /** \brief Extent used to calc histogram */
+    QgsRectangle extent;
+
+    /** \brief Histogram is valid */
+    bool valid;
+};
+#endif