fix histogram and stats compute (set bapprox=false) for gdal provider ; fix histogram display for non-Byte data

Author: etiennesky

src/app/qgsrasterlayerproperties.cpp

@@ -1300,8 +1300,7 @@ void QgsRasterLayerProperties::refreshHistogram()
 
   if ( ! computeHistogram( false ) )
   {
-    // TODO - check raster min/max and min/max of default histogram
-    QgsDebugMsg( QString( "raster does not have cached histo" ) );
+    QgsDebugMsg( QString( "raster does not have cached histogram" ) );
     stackedWidget2->setCurrentIndex( 2 );
     return;
   }
@@ -1419,6 +1418,8 @@ void QgsRasterLayerProperties::refreshHistogram()
   bool myFirstIteration = true;
   /* get selected band list, if mHistoShowBands != ShowAll */
   QList< int > mySelectedBands = histoSelectedBands();
+  double myBinXStep = 1;
+  double myBinX = 0;
 
   for ( int myIteratorInt = 1;
         myIteratorInt <= myBandCountInt;
@@ -1442,15 +1443,28 @@ void QgsRasterLayerProperties::refreshHistogram()
     QVector<double> myX2Data;
     QVector<double> myY2Data;
 #endif
+    // calculate first bin x value and bin step size if not Byte data
+    if ( mRasterLayer->dataProvider()->srcDataType( myIteratorInt ) != QgsRasterDataProvider::Byte )
+    {
+      myBinXStep = myRasterBandStats.range / BINCOUNT;
+      myBinX = myRasterBandStats.minimumValue + myBinXStep / 2.0;
+    }
+    else
+    {
+      myBinXStep = 1;
+      myBinX = 0;
+    }
+
     for ( int myBin = 0; myBin < BINCOUNT; myBin++ )
     {
       int myBinValue = myRasterBandStats.histogramVector->at( myBin );
 #if defined(QWT_VERSION) && QWT_VERSION>=0x060000
-      data << QPointF( myBin, myBinValue );
+      data << QPointF( myBinX, myBinValue );
 #else
-      myX2Data.append( double( myBin ) );
+      myX2Data.append( double( myBinX ) );
       myY2Data.append( double( myBinValue ) );
 #endif
+      myBinX += myBinXStep;
     }
 #if defined(QWT_VERSION) && QWT_VERSION>=0x060000
     mypCurve->setSamples( data );
@@ -1472,9 +1486,10 @@ void QgsRasterLayerProperties::refreshHistogram()
   // for x axis use band pixel values rather than gdal hist. bin values
   // subtract -0.5 to prevent rounding errors
   // see http://www.gdal.org/classGDALRasterBand.html#3f8889607d3b2294f7e0f11181c201c8
+  // fix x range for non-Byte data
   mpPlot->setAxisScale( QwtPlot::xBottom,
-                        mHistoMin - 0.5,
-                        mHistoMax + 0.5 );
+                        mHistoMin - myBinXStep / 2,
+                        mHistoMax + myBinXStep / 2 );
 
   mpPlot->replot();
 
@@ -1982,7 +1997,6 @@ void QgsRasterLayerProperties::histoActionTriggered( QAction* action )
   // Load actions
   else if ( actionName.left( 5 ) == "Load " )
   {
-    QgsRasterBandStats myRasterBandStats;
     QVector<int> myBands;
     double minMaxValues[2];
     bool ok = false;

src/core/raster/qgsrasterlayer.cpp

@@ -371,6 +371,8 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
     return myNullReturnStats;
   }
 
+  // TODO this is buggy - because the stats might have changed (e.g. theIgnoreOutOfRangeFlag in populateHistogram())
+  // should have a pointer to the stats instead
   QgsRasterBandStats myRasterBandStats = mRasterStatsList[theBandNo - 1];
   myRasterBandStats.bandNumber = theBandNo;
 
@@ -384,6 +386,7 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
   QgsDebugMsg( "adding stats to stats collection at position " + QString::number( theBandNo - 1 ) );
   //add this band to the class stats map
   mRasterStatsList[theBandNo - 1] = myRasterBandStats;
+
   emit drawingProgress( mHeight, mHeight ); //reset progress
   QgsDebugMsg( "Stats collection completed returning" );
   return myRasterBandStats;
@@ -1422,7 +1425,7 @@ QPixmap QgsRasterLayer::paletteAsPixmap( int theBandNumber )
 }
 
 /*
- * @param theBandNoInt - which band to compute the histogram for
+ * @param theBandNoInt - which band to find out if has a cached histogram
  * @param theBinCountInt - how many 'bins' to categorise the data into
  */
 bool QgsRasterLayer::hasCachedHistogram( int theBandNo, int theBinCount )

src/providers/gdal/qgsgdalprovider.cpp

@@ -1426,24 +1426,25 @@ bool QgsGdalProvider::hasCachedHistogram( int theBandNo, int theBinCount )
   if ( ! myGdalBand )
     return false;
 
-  // get default histo
+  // get default histogram with force=false to see if there is a cached histo
   double myMinVal, myMaxVal;
-  QgsRasterBandStats theBandStats = bandStatistics( theBandNo );
-  double myerval = ( theBandStats.maximumValue - theBandStats.minimumValue ) / theBinCount;
-  myMinVal = theBandStats.minimumValue - 0.1*myerval;
-  myMaxVal = theBandStats.maximumValue + 0.1*myerval;
-  int *myHistogramArray=0;
-  CPLErr myError = GDALGetDefaultHistogram( myGdalBand, &myMinVal, &myMaxVal, 
-                                            &theBinCount, &myHistogramArray, false,
-                                            NULL, NULL );
-  if( myHistogramArray )
+  int myBinCount;
+  int *myHistogramArray = 0;
+  CPLErr myError = GDALGetDefaultHistogram( myGdalBand, &myMinVal, &myMaxVal,
+                   &myBinCount, &myHistogramArray, false,
+                   NULL, NULL );
+  if ( myHistogramArray )
     VSIFree( myHistogramArray );
 
   // if there was any error/warning assume the histogram is not valid or non-existent
   if ( myError != CE_None )
     return false;
-  return true;
 
+  // make sure the cached histo has the same number of bins than requested
+  if ( myBinCount != theBinCount )
+    return false;
+
+  return true;
 }
 
 void QgsGdalProvider::populateHistogram( int theBandNo, QgsRasterBandStats & theBandStats, int theBinCount, bool theIgnoreOutOfRangeFlag, bool theHistogramEstimatedFlag )
@@ -1460,34 +1461,58 @@ void QgsGdalProvider::populateHistogram( int theBandNo, QgsRasterBandStats & the
        theIgnoreOutOfRangeFlag != theBandStats.isHistogramOutOfRange ||
        theHistogramEstimatedFlag != theBandStats.isHistogramEstimated )
   {
+    QgsDebugMsg( "Computing histogram" );
     theBandStats.histogramVector->clear();
     theBandStats.isHistogramEstimated = theHistogramEstimatedFlag;
     theBandStats.isHistogramOutOfRange = theIgnoreOutOfRangeFlag;
     int *myHistogramArray = new int[theBinCount];
 
-#if 0
-    CPLErr GDALRasterBand::GetHistogram(
-      double       dfMin,
-      double      dfMax,
-      int     nBuckets,
-      int *   panHistogram,
-      int     bIncludeOutOfRange,
-      int     bApproxOK,
-      GDALProgressFunc    pfnProgress,
-      void *      pProgressData
-    )
-#endif
-
     QgsGdalProgress myProg;
     myProg.type = ProgressHistogram;
     myProg.provider = this;
-    double myerval = ( theBandStats.maximumValue - theBandStats.minimumValue ) / theBinCount;
 
+#if 0 // this is the old method
+
+    double myerval = ( theBandStats.maximumValue - theBandStats.minimumValue ) / theBinCount;
     GDALGetRasterHistogram( myGdalBand, theBandStats.minimumValue - 0.1*myerval,
                             theBandStats.maximumValue + 0.1*myerval, theBinCount, myHistogramArray,
                             theIgnoreOutOfRangeFlag, theHistogramEstimatedFlag, progressCallback,
                             &myProg ); //this is the arg for our custom gdal progress callback
 
+#else // this is the new method, which gets a "Default" histogram
+
+    // calculate min/max like in GDALRasterBand::GetDefaultHistogram, but don't call it directly
+    // because there is no bApproxOK argument - that is lacking from the API
+    double myMinVal, myMaxVal;
+    const char* pszPixelType = GDALGetMetadataItem( myGdalBand, "PIXELTYPE", "IMAGE_STRUCTURE" );
+    int bSignedByte = ( pszPixelType != NULL && EQUAL( pszPixelType, "SIGNEDBYTE" ) );
+
+    if ( GDALGetRasterDataType( myGdalBand ) == GDT_Byte && !bSignedByte )
+    {
+      myMinVal = -0.5;
+      myMaxVal = 255.5;
+    }
+    else
+    {
+      CPLErr eErr = CE_Failure;
+      double dfHalfBucket = 0;
+      eErr = GDALGetRasterStatistics( myGdalBand, TRUE, TRUE, &myMinVal, &myMaxVal, NULL, NULL );
+      if ( eErr != CE_None )
+        return;
+      dfHalfBucket = ( myMaxVal - myMinVal ) / ( 2 * theBinCount );
+      myMinVal -= dfHalfBucket;
+      myMaxVal += dfHalfBucket;
+    }
+
+    CPLErr myError = GDALGetRasterHistogram( myGdalBand, myMinVal, myMaxVal,
+                     theBinCount, myHistogramArray,
+                     theIgnoreOutOfRangeFlag, theHistogramEstimatedFlag, progressCallback,
+                     &myProg ); //this is the arg for our custom gdal progress callback
+    if ( myError != CE_None )
+      return;
+
+#endif
+
     for ( int myBin = 0; myBin < theBinCount; myBin++ )
     {
       if ( myHistogramArray[myBin] < 0 ) //can't have less than 0 pixels of any value
@@ -2063,7 +2088,8 @@ QgsRasterBandStats QgsGdalProvider::bandStatistics( int theBandNo )
 {
   GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, theBandNo );
   QgsRasterBandStats myRasterBandStats;
-  int bApproxOK = true;
+  // int bApproxOK = true;
+  int bApproxOK = false; //as we asked for stats, don't get approx values
   double pdfMin;
   double pdfMax;
   double pdfMean;
@@ -2072,16 +2098,6 @@ QgsRasterBandStats QgsGdalProvider::bandStatistics( int theBandNo )
   myProg.type = ProgressHistogram;
   myProg.provider = this;
 
-  // double myerval =
-  //   GDALComputeRasterStatistics (
-  //      myGdalBand, bApproxOK, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev,
-  //      progressCallback, &myProg ) ;
-  // double myerval =
-  //   GDALGetRasterStatistics ( myGdalBand, bApproxOK, TRUE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev);
-  // double myerval =
-  //   GDALGetRasterStatisticsProgress ( myGdalBand, bApproxOK, TRUE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev,
-  //           progressCallback, &myProg );
-
   // try to fetch the cached stats (bForce=FALSE)
   CPLErr myerval =
     GDALGetRasterStatistics( myGdalBand, bApproxOK, FALSE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev );