Skip to content

Commit d794522

Browse files
blazekblazek
committedJul 26, 2012
new raster histogram
·
ltr-3_40Before-merge-new_vector_api
1 parent 47cc4f2 commit d794522

File tree

4 files changed

+271
-1
lines changed

4 files changed

+271
-1
lines changed
 

‎src/core/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -390,6 +390,7 @@ SET(QGIS_CORE_HDRS
390390
raster/qgspseudocolorshader.h
391391
raster/qgsrasterpyramid.h
392392
raster/qgsrasterbandstats.h
393+
raster/qgsrasterhistogram.h
393394
raster/qgsrasterinterface.h
394395
raster/qgsrasterlayer.h
395396
raster/qgsrastertransparency.h

‎src/core/qgsrasterdataprovider.cpp

Lines changed: 165 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@
2323
#include <QMap>
2424
#include <QByteArray>
2525

26+
#include <cmath>
27+
2628
void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle
2729
const & viewExtent, int width,
2830
int height,
@@ -270,7 +272,6 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
270272
return ba;
271273
}
272274

273-
274275
QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
275276
{
276277
double myNoDataValue = noDataValue();
@@ -426,6 +427,169 @@ QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
426427
return myRasterBandStats;
427428
}
428429

430+
QgsRasterHistogram QgsRasterDataProvider::histogram( int theBandNo,
431+
double theMinimum, double theMaximum,
432+
int theBinCount,
433+
const QgsRectangle & theExtent,
434+
int theSampleSize,
435+
bool theIncludeOutOfRange )
436+
{
437+
QgsRasterHistogram myHistogram;
438+
myHistogram.bandNumber = theBandNo;
439+
myHistogram.minimum = theMinimum;
440+
myHistogram.maximum = theMaximum;
441+
myHistogram.includeOutOfRange = theIncludeOutOfRange;
442+
//myHistogram.sampleSize = theSampleSize
443+
444+
// First calc defaults
445+
QgsRectangle myExtent = theExtent.isEmpty() ? extent() : theExtent;
446+
myHistogram.extent = myExtent;
447+
448+
int myWidth, myHeight;
449+
if ( theSampleSize > 0 )
450+
{
451+
// Calc resolution from theSampleSize
452+
double xRes, yRes;
453+
xRes = yRes = sqrt(( myExtent.width() * myExtent.height() ) / theSampleSize );
454+
455+
// But limit by physical resolution
456+
if ( capabilities() & Size )
457+
{
458+
double srcXRes = extent().width() / xSize();
459+
double srcYRes = extent().height() / ySize();
460+
if ( xRes < srcXRes ) xRes = srcXRes;
461+
if ( yRes < srcYRes ) yRes = srcYRes;
462+
}
463+
464+
myWidth = static_cast <int>( myExtent.width() / xRes );
465+
myHeight = static_cast <int>( myExtent.height() / yRes );
466+
}
467+
else
468+
{
469+
if ( capabilities() & Size )
470+
{
471+
myWidth = xSize();
472+
myHeight = ySize();
473+
}
474+
else
475+
{
476+
myWidth = 1000;
477+
myHeight = 1000;
478+
}
479+
}
480+
myHistogram.width = myWidth;
481+
myHistogram.height = myHeight;
482+
QgsDebugMsg( QString( "myWidth = %1 myHeight = %2" ).arg( myWidth ).arg( myHeight ) );
483+
484+
double myNoDataValue = noDataValue();
485+
int myDataType = dataType( theBandNo );
486+
487+
int myBinCount = theBinCount;
488+
if ( myBinCount == 0 )
489+
{
490+
if ( myDataType == QgsRasterDataProvider::Byte )
491+
{
492+
myBinCount = 256; // Cannot store more values in byte
493+
}
494+
else
495+
{
496+
// 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;
497+
myBinCount = myWidth * myHeight;
498+
if ( myBinCount > 1000 ) myBinCount = 1000;
499+
}
500+
}
501+
myHistogram.binCount = theBinCount;
502+
QgsDebugMsg( QString( "myBinCount = %1" ).arg( myBinCount ) );
503+
504+
// Check if we have cached
505+
foreach( QgsRasterHistogram histogram, mHistograms )
506+
{
507+
if ( histogram.bandNumber == theBandNo &&
508+
histogram.minimum == theMinimum &&
509+
histogram.maximum == theMaximum &&
510+
histogram.binCount == myBinCount &&
511+
histogram.extent == myExtent &&
512+
histogram.width == myWidth &&
513+
histogram.height == myHeight &&
514+
histogram.includeOutOfRange == theIncludeOutOfRange )
515+
{
516+
return histogram;
517+
}
518+
}
519+
520+
myHistogram.histogramVector.resize( myBinCount );
521+
522+
int myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
523+
myXBlockSize = xBlockSize();
524+
myYBlockSize = yBlockSize();
525+
526+
if ( myXBlockSize == 0 || myYBlockSize == 0 ) // should not happen
527+
{
528+
return myHistogram;
529+
}
530+
531+
myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
532+
myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
533+
534+
void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) );
535+
536+
int myBandXSize = xSize();
537+
int myBandYSize = ySize();
538+
double myXRes = myExtent.width() / myWidth;
539+
double myYRes = myExtent.height() / myHeight;
540+
double binSize = ( theMaximum - theMinimum ) / myBinCount;
541+
for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
542+
{
543+
for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
544+
{
545+
int myPartWidth = qMin( myXBlockSize, myWidth - iXBlock * myXBlockSize );
546+
int myPartHeight = qMin( myYBlockSize, myHeight - iYBlock * myYBlockSize );
547+
548+
double xmin = myExtent.xMinimum() + iXBlock * myXBlockSize * myXRes;
549+
double xmax = xmin + myPartWidth * myXRes;
550+
double ymin = myExtent.yMaximum() - iYBlock * myYBlockSize * myYRes;
551+
double ymax = ymin - myPartHeight * myYRes;
552+
553+
QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
554+
555+
readBlock( theBandNo, myPartExtent, myPartWidth, myPartHeight, myData );
556+
557+
// Collect the histogram counts.
558+
for ( int iY = 0; iY < myPartHeight; iY++ )
559+
{
560+
for ( int iX = 0; iX < myPartWidth; iX++ )
561+
{
562+
double myValue = readValue( myData, myDataType, iX + ( iY * myPartWidth ) );
563+
//QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) );
564+
565+
if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
566+
{
567+
continue; // NULL
568+
}
569+
570+
int myBinIndex = static_cast <int>( floor(( myValue - theMinimum ) / binSize ) ) ;
571+
if (( myBinIndex < 0 || myBinIndex > ( myBinCount - 1 ) ) && !theIncludeOutOfRange )
572+
{
573+
continue;
574+
}
575+
if ( myBinIndex < 0 ) myBinIndex = 0;
576+
if ( myBinIndex > ( myBinCount - 1 ) ) myBinIndex = myBinCount - 1;
577+
578+
myHistogram.histogramVector[myBinIndex] += 1;
579+
myHistogram.nonNullCount++;
580+
}
581+
}
582+
} //end of column wise loop
583+
} //end of row wise loop
584+
585+
CPLFree( myData );
586+
587+
myHistogram.valid = true;
588+
mHistograms.append( myHistogram );
589+
590+
return myHistogram;
591+
}
592+
429593
double QgsRasterDataProvider::readValue( void *data, int type, int index )
430594
{
431595
if ( !data )

‎src/core/qgsrasterdataprovider.h

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030
#include "qgsrasterpyramid.h"
3131
#include "qgscoordinatereferencesystem.h"
3232
#include "qgsrasterbandstats.h"
33+
#include "qgsrasterhistogram.h"
3334

3435
#include "cpl_conv.h"
3536
#include <cmath>
@@ -314,6 +315,22 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
314315
Q_UNUSED( theIgnoreOutOfRangeFlag ); Q_UNUSED( theThoroughBandScanFlag );
315316
}
316317

318+
/** \brief Get histogram. Histograms are cached in providers.
319+
* @param theBandNo The band (number).
320+
* @param theMinimum Minimum value.
321+
* @param theMaximum Maximum value.
322+
* @param theBinCount Number of bins (intervals,buckets). If 0, the number of bins is decided automaticaly according to data type, raster size etc.
323+
* @param theExtent Extent used to calc histogram, if empty, whole raster extent is used.
324+
* @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.
325+
* @return Vector of non NULL cell counts for each bin.
326+
*/
327+
virtual QgsRasterHistogram histogram( int theBandNo,
328+
double theMinimum, double theMaximum,
329+
int theBinCount = 0,
330+
const QgsRectangle & theExtent = QgsRectangle(),
331+
int theSampleSize = 0,
332+
bool theIncludeOutOfRange = false );
333+
317334
/** \brief Create pyramid overviews */
318335
virtual QString buildPyramids( const QList<QgsRasterPyramid> & thePyramidList,
319336
const QString & theResamplingMethod = "NEAREST",
@@ -503,6 +520,9 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
503520
bool mValidNoDataValue;
504521

505522
QgsRectangle mExtent;
523+
524+
/** \brief List of cached histograms, all bands mixed */
525+
QList <QgsRasterHistogram> mHistograms;
506526
};
507527

508528
#endif

‎src/core/raster/qgsrasterhistogram.h

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
/***************************************************************************
2+
qgsrasterhistogram.h
3+
-------------------
4+
begin : July 2012
5+
copyright : (C) 2012 by Radim Blazek
6+
email : radim dot blazek at gmail dot com
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#ifndef QGSRASTERHISTOGRAM
19+
#define QGSRASTERHISTOGRAM
20+
21+
#include <QString>
22+
#include <QVector>
23+
24+
#include <limits>
25+
26+
/** \ingroup core
27+
* The QgsRasterHistogram is a container for histogram of a single raster band.
28+
* It is used to cache computed histograms in raster providers.
29+
*/
30+
class CORE_EXPORT QgsRasterHistogram
31+
{
32+
public:
33+
typedef QVector<int> HistogramVector;
34+
35+
QgsRasterHistogram()
36+
{
37+
bandNumber = 0;
38+
binCount = 0;
39+
nonNullCount = 0;
40+
//sampleSize = 0;
41+
includeOutOfRange = false;
42+
maximum = 0;
43+
minimum = 0;
44+
width = 0;
45+
height = 0;
46+
valid = false;
47+
}
48+
49+
/** \brief The gdal band number (starts at 1)*/
50+
int bandNumber;
51+
52+
/** \brief Number of bins (intervals,buckets) in histogram. */
53+
54+
/** \brief The number of non NULL cells used to calculate histogram. */
55+
int nonNullCount;
56+
57+
/** \brief Approximate number of cells used to calc histogram. Approximately
58+
* width * height. */
59+
//int sampleSize;
60+
61+
/** \brief Whether histogram includes out of range values (in first and last bin) */
62+
bool includeOutOfRange;
63+
64+
/** \brief Store the histogram for a given layer */
65+
HistogramVector histogramVector;
66+
67+
/** \brief The maximum histogram value. */
68+
double maximum;
69+
70+
/** \brief The minimum histogram value. */
71+
double minimum;
72+
73+
/** \brief Number of columns used to calc histogram */
74+
int width;
75+
76+
/** \brief Number of rows used to calc histogram */
77+
int height;
78+
79+
/** \brief Extent used to calc histogram */
80+
QgsRectangle extent;
81+
82+
/** \brief Histogram is valid */
83+
bool valid;
84+
};
85+
#endif

0 commit comments

Comments
 (0)
Please sign in to comment.