Skip to content

Commit

Permalink
Backported raster speed up (using providers to compute stats natively…
Browse files Browse the repository at this point in the history
…) to release branch.
  • Loading branch information
timlinux committed Sep 8, 2011
1 parent 625f5fa commit 4db175e
Show file tree
Hide file tree
Showing 5 changed files with 385 additions and 101 deletions.
203 changes: 194 additions & 9 deletions src/core/qgsrasterdataprovider.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
* (at your option) any later version. *
* *
***************************************************************************/
/* $Id$ */

#include "qgsrasterdataprovider.h"
#include "qgsrasterprojector.h"
Expand Down Expand Up @@ -54,7 +53,8 @@ void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle const & viewExt
// TODO: init data by nulls

// If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 ) return;
if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 )
return;

// Allocate memory for not projected source data
int mySize = dataTypeSize( bandNo ) / 8;
Expand Down Expand Up @@ -174,6 +174,7 @@ QString QgsRasterDataProvider::capabilitiesString() const

bool QgsRasterDataProvider::identify( const QgsPoint& thePoint, QMap<QString, QString>& theResults )
{
Q_UNUSED( thePoint );
theResults.clear();
return false;
}
Expand All @@ -200,31 +201,31 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
double d;
switch ( type )
{
case QgsRasterDataProvider::Byte:
case Byte:
uc = ( unsigned char )noval;
memcpy( data, &uc, size );
break;
case QgsRasterDataProvider::UInt16:
case UInt16:
us = ( unsigned short )noval;
memcpy( data, &us, size );
break;
case QgsRasterDataProvider::Int16:
case Int16:
s = ( short )noval;
memcpy( data, &s, size );
break;
case QgsRasterDataProvider::UInt32:
case UInt32:
ui = ( unsigned int )noval;
memcpy( data, &ui, size );
break;
case QgsRasterDataProvider::Int32:
case Int32:
i = ( int )noval;
memcpy( data, &i, size );
break;
case QgsRasterDataProvider::Float32:
case Float32:
f = ( float )noval;
memcpy( data, &f, size );
break;
case QgsRasterDataProvider::Float64:
case Float64:
d = ( double )noval;
memcpy( data, &d, size );
break;
Expand All @@ -234,4 +235,188 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
return ba;
}


QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo )
{
double myNoDataValue = noDataValue();
QgsRasterBandStats myRasterBandStats;
myRasterBandStats.elementCount = 0; // because we'll be counting only VALID pixels later

int myDataType = dataType( theBandNo );

int myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize;
myXBlockSize = xBlockSize();
myYBlockSize = yBlockSize();

myNXBlocks = ( xSize() + myXBlockSize - 1 ) / myXBlockSize;
myNYBlocks = ( ySize() + myYBlockSize - 1 ) / myYBlockSize;

void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) );

// unfortunately we need to make two passes through the data to calculate stddev
bool myFirstIterationFlag = true;

int myBandXSize = xSize();
int myBandYSize = ySize();
for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
{
for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
{
int nXValid, nYValid;
readBlock( theBandNo, iXBlock, iYBlock, myData );

// Compute the portion of the block that is valid
// for partial edge blocks.
if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
nXValid = myBandXSize - iXBlock * myXBlockSize;
else
nXValid = myXBlockSize;

if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
nYValid = myBandYSize - iYBlock * myYBlockSize;
else
nYValid = myYBlockSize;

// Collect the histogram counts.
for ( int iY = 0; iY < nYValid; iY++ )
{
for ( int iX = 0; iX < nXValid; iX++ )
{
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
//QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) );

if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
{
continue; // NULL
}

myRasterBandStats.sum += myValue;
++myRasterBandStats.elementCount;
//only use this element if we have a non null element
if ( myFirstIterationFlag )
{
//this is the first iteration so initialise vars
myFirstIterationFlag = false;
myRasterBandStats.minimumValue = myValue;
myRasterBandStats.maximumValue = myValue;
} //end of true part for first iteration check
else
{
//this is done for all subsequent iterations
if ( myValue < myRasterBandStats.minimumValue )
{
myRasterBandStats.minimumValue = myValue;
}
if ( myValue > myRasterBandStats.maximumValue )
{
myRasterBandStats.maximumValue = myValue;
}
} //end of false part for first iteration check
}
}
} //end of column wise loop
} //end of row wise loop


//end of first pass through data now calculate the range
myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
//calculate the mean
myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;

//for the second pass we will get the sum of the squares / mean
for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ )
{
for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ )
{
int nXValid, nYValid;

readBlock( theBandNo, iXBlock, iYBlock, myData );

// Compute the portion of the block that is valid
// for partial edge blocks.
if (( iXBlock + 1 ) * myXBlockSize > myBandXSize )
nXValid = myBandXSize - iXBlock * myXBlockSize;
else
nXValid = myXBlockSize;

if (( iYBlock + 1 ) * myYBlockSize > myBandYSize )
nYValid = myBandYSize - iYBlock * myYBlockSize;
else
nYValid = myYBlockSize;

// Collect the histogram counts.
for ( int iY = 0; iY < nYValid; iY++ )
{
for ( int iX = 0; iX < nXValid; iX++ )
{
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
//QgsDebugMsg ( "myValue = " + QString::number(myValue) );

if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) )
{
continue; // NULL
}

myRasterBandStats.sumOfSquares += static_cast < double >
( pow( myValue - myRasterBandStats.mean, 2 ) );
}
}
} //end of column wise loop
} //end of row wise loop

//divide result by sample size - 1 and get square root to get stdev
myRasterBandStats.stdDev = static_cast < double >( sqrt( myRasterBandStats.sumOfSquares /
( myRasterBandStats.elementCount - 1 ) ) );

#ifdef QGISDEBUG
QgsLogger::debug( "************ STATS **************", 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "VALID NODATA", mValidNoDataValue, 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "NULL", noDataValue() , 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "MIN", myRasterBandStats.minimumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "MAX", myRasterBandStats.maximumValue, 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "RANGE", myRasterBandStats.range, 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "MEAN", myRasterBandStats.mean, 1, __FILE__, __FUNCTION__, __LINE__ );
QgsLogger::debug( "STDDEV", myRasterBandStats.stdDev, 1, __FILE__, __FUNCTION__, __LINE__ );
#endif

CPLFree( myData );
myRasterBandStats.statsGathered = true;
return myRasterBandStats;
}

double QgsRasterDataProvider::readValue( void *data, int type, int index )
{
if ( !data )
return mValidNoDataValue ? noDataValue() : 0.0;

switch ( type )
{
case Byte:
return ( double )(( GByte * )data )[index];
break;
case UInt16:
return ( double )(( GUInt16 * )data )[index];
break;
case Int16:
return ( double )(( GInt16 * )data )[index];
break;
case UInt32:
return ( double )(( GUInt32 * )data )[index];
break;
case Int32:
return ( double )(( GInt32 * )data )[index];
break;
case Float32:
return ( double )(( float * )data )[index];
break;
case Float64:
return ( double )(( double * )data )[index];
break;
default:
QgsLogger::warning( "GDAL data type is not supported" );
}

return mValidNoDataValue ? noDataValue() : 0.0;
}

// ENDS
41 changes: 27 additions & 14 deletions src/core/qgsrasterdataprovider.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
* (at your option) any later version. *
* *
***************************************************************************/
/* $Id$ */

/* Thank you to Marco Hugentobler for the original vector DataProvider */

Expand All @@ -29,12 +28,12 @@
#include "qgscolorrampshader.h"
#include "qgsrasterpyramid.h"
#include "qgscoordinatereferencesystem.h"

#include "qgsrasterbandstats.h"
#include "cpl_conv.h"
#include <cmath>

class QImage;
class QgsPoint;
class QgsRasterBandStats;
class QByteArray;

#define TINY_VALUE std::numeric_limits<double>::epsilon() * 20
Expand Down Expand Up @@ -110,10 +109,11 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider
};

// Progress types
enum Progress
enum RasterProgressType
{
ProgressHistogram = 0,
ProgressPyramids = 1
ProgressPyramids = 1,
ProgressStatistics = 2
};

QgsRasterDataProvider();
Expand Down Expand Up @@ -184,6 +184,7 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider
*/
virtual int srcDataType( int bandNo ) const
{
Q_UNUSED( bandNo );
return QgsRasterDataProvider::UnknownDataType;
}

Expand Down Expand Up @@ -234,6 +235,7 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider
/** Returns data type for the band specified by number */
virtual int colorInterpretation( int theBandNo ) const
{
Q_UNUSED( theBandNo );
return QgsRasterDataProvider::UndefinedColorInterpretation;
}

Expand Down Expand Up @@ -315,21 +317,27 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider

/** read block of data */
// TODO clarify what happens on the last block (the part outside raster)
virtual void readBlock( int bandNo, int xBlock, int yBlock, void *data ) {}
virtual void readBlock( int bandNo, int xBlock, int yBlock, void *data )
{ Q_UNUSED( bandNo ); Q_UNUSED( xBlock ); Q_UNUSED( yBlock ); Q_UNUSED( data ); }

/** read block of data using give extent and size */
virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data ) {};
virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data )
{ Q_UNUSED( bandNo ); Q_UNUSED( viewExtent ); Q_UNUSED( width ); Q_UNUSED( height ); Q_UNUSED( data ); }

/** read block of data using give extent and size */
virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, QgsCoordinateReferenceSystem theSrcCRS, QgsCoordinateReferenceSystem theDestCRS, void *data );

/* Read a value from a data block at a given index. */
virtual double readValue( void *data, int type, int index );

/** value representing null data */
virtual double noDataValue() const { return 0; }

virtual double minimumValue( int bandNo )const { return 0; }
virtual double maximumValue( int bandNo )const { return 0; }
virtual double minimumValue( int bandNo ) const { Q_UNUSED( bandNo ); return 0; }
virtual double maximumValue( int bandNo ) const { Q_UNUSED( bandNo ); return 0; }

virtual QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const { return QList<QgsColorRampShader::ColorRampItem>(); }
virtual QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo ) const
{ Q_UNUSED( bandNo ); return QList<QgsColorRampShader::ColorRampItem>(); }

// Defined in parent
/** \brief Returns the sublayers of this layer - Useful for providers that manage their own layers, such as WMS */
Expand All @@ -345,12 +353,14 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider
int theBinCountInt = 256,
bool theIgnoreOutOfRangeFlag = true,
bool theThoroughBandScanFlag = false
) {};
)
{ Q_UNUSED( theBandNoInt ); Q_UNUSED( theBandStats ); Q_UNUSED( theBinCountInt ); Q_UNUSED( theIgnoreOutOfRangeFlag ); Q_UNUSED( theThoroughBandScanFlag ); }

/** \brief Create pyramid overviews */
virtual QString buildPyramids( const QList<QgsRasterPyramid> & thePyramidList,
const QString & theResamplingMethod = "NEAREST",
bool theTryInternalFlag = false ) { return "FAILED_NOT_SUPPORTED"; };
bool theTryInternalFlag = false )
{ Q_UNUSED( thePyramidList ); Q_UNUSED( theResamplingMethod ); Q_UNUSED( theTryInternalFlag ); return "FAILED_NOT_SUPPORTED"; };

/** \brief Accessor for ths raster layers pyramid list. A pyramid list defines the
* POTENTIAL pyramids that can be in a raster. To know which of the pyramid layers
Expand All @@ -359,13 +369,16 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider
*/
virtual QList<QgsRasterPyramid> buildPyramidList() { return QList<QgsRasterPyramid>(); };

/** If the provider supports it, return band stats for the
given band. Default behaviour is to blockwise read the data
and generate the stats unless the provider overloads this function. */
virtual QgsRasterBandStats bandStatistics( int theBandNo );

/** \brief helper function to create zero padded band names */
QString generateBandName( int theBandNumber )
{
return tr( "Band" ) + QString( " %1" ) .arg( theBandNumber, 1 + ( int ) log10(( float ) bandCount() ), 10, QChar( '0' ) );
}

};

/**
* Get metadata in a format suitable for feeding directly
Expand Down
Loading

0 comments on commit 4db175e

Please sign in to comment.