Skip to content

Commit

Permalink
fix histogram and stats compute (set bapprox=false) for gdal provider…
Browse files Browse the repository at this point in the history 
… ; fix histogram display for non-Byte data
  • Loading branch information
@etiennesky
etiennesky committed Jul 14, 2012
1 parent 3c13f3c commit e74a28b
Show file tree
Hide file tree
Showing 4 changed files with 77 additions and 208 deletions.
28 changes: 21 additions & 7 deletions src/app/qgsrasterlayerproperties.cpp
View file
Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
Expand All @@ -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 );
Expand All @@ -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();

Expand Down Expand Up @@ -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;
Expand Down
5 changes: 4 additions & 1 deletion src/core/raster/qgsrasterlayer.cpp
View file
Expand Up @@ -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;

Expand All @@ -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;
Expand Down Expand Up @@ -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 )
Expand Down
88 changes: 52 additions & 36 deletions src/providers/gdal/qgsgdalprovider.cpp
View file
Expand Up @@ -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 )
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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 );
Expand Down

0 comments on commit e74a28b

Please sign in to comment.