Skip to content
Permalink
Browse files

new raster histogram

  • Loading branch information
blazek committed Jul 26, 2012
1 parent 47cc4f2 commit d79452289320ccf9d852d422d75082d2c8c29e4f
@@ -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
@@ -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 )
@@ -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
@@ -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

0 comments on commit d794522

Please sign in to comment.
You can’t perform that action at this time.