Skip to content

Commit

Permalink
-Fix no data value from being included in histogram and stats
Browse files Browse the repository at this point in the history 
-Closes ticket #1380

git-svn-id: http://svn.osgeo.org/qgis/trunk@11189 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
ersts committed Jul 28, 2009
1 parent da0c844 commit 1c3390e
Showing 1 changed file with 33 additions and 33 deletions.
66 changes: 33 additions & 33 deletions src/core/raster/qgsrasterlayer.cpp
View file
Expand Up @@ -70,6 +70,13 @@ email : tim at linfiniti.com
# endif
#endif

// Comparison value for equality; i.e., we shouldn't directly compare two
// floats so it's better to take their difference and see if they're within
// a certain range -- in this case twenty times the smallest value that
// doubles can take for the current system. (Yes, 20 was arbitrary.)
#define TINY_VALUE std::numeric_limits<double>::epsilon() * 20


QgsRasterLayer::QgsRasterLayer(
QString const & path,
QString const & baseName,
Expand Down Expand Up @@ -106,7 +113,7 @@ QgsRasterLayer::QgsRasterLayer(

mBandCount = 0;
mHasPyramids = false;
mNoDataValue = -9999;
mNoDataValue = -9999.0;
mValidNoDataValue = false;

mGdalBaseDataset = 0;
Expand Down Expand Up @@ -759,18 +766,11 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
myNXBlocks = ( GDALGetRasterXSize( myGdalBand ) + myXBlockSize - 1 ) / myXBlockSize;
myNYBlocks = ( GDALGetRasterYSize( myGdalBand ) + myYBlockSize - 1 ) / myYBlockSize;

void *myData = CPLMalloc( myXBlockSize * myYBlockSize * GDALGetDataTypeSize( myDataType ) / 8 );
void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( GDALGetDataTypeSize( myDataType ) / 8 ) );

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

// Comparison value for equality; i.e., we shouldn't directly compare two
// floats so it's better to take their difference and see if they're within
// a certain range -- in this case twenty times the smallest value that
// doubles can take for the current system. (Yes, 20 was arbitrary.)
double myPrecision = std::numeric_limits<double>::epsilon() * 20;
Q_UNUSED( myPrecision );

//ifdefs below to remove compiler warning about unused vars
#ifdef QGISDEBUG
int success;
Expand Down Expand Up @@ -841,10 +841,9 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
{
for ( int iX = 0; iX < nXValid; iX++ )
{
double my = readValue( myData, myDataType, iX + iY * myXBlockSize );
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );

//if ( mValidNoDataValue && (fabs(my - mNoDataValue) < myPrecision || my == mNoDataValue || my != my))
if ( mValidNoDataValue && ( my == mNoDataValue || my != my ) )
if ( mValidNoDataValue && ( fabs( myValue - mNoDataValue ) <= TINY_VALUE || myValue != myValue ) )
{
continue; // NULL
}
Expand All @@ -854,23 +853,23 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
{
//this is the first iteration so initialise vars
myFirstIterationFlag = false;
myRasterBandStats.minimumValue = my;
myRasterBandStats.maximumValue = my;
myRasterBandStats.minimumValue = myValue;
myRasterBandStats.maximumValue = myValue;
++myRasterBandStats.elementCount;
} //end of true part for first iteration check
else
{
//this is done for all subsequent iterations
if ( my < myRasterBandStats.minimumValue )
if ( myValue < myRasterBandStats.minimumValue )
{
myRasterBandStats.minimumValue = my;
myRasterBandStats.minimumValue = myValue;
}
if ( my > myRasterBandStats.maximumValue )
if ( myValue > myRasterBandStats.maximumValue )
{
myRasterBandStats.maximumValue = my;
myRasterBandStats.maximumValue = myValue;
}

myRasterBandStats.sum += my;
myRasterBandStats.sum += myValue;
++myRasterBandStats.elementCount;
} //end of false part for first iteration check
}
Expand Down Expand Up @@ -912,16 +911,15 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
{
for ( int iX = 0; iX < nXValid; iX++ )
{
double my = readValue( myData, myDataType, iX + iY * myXBlockSize );
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );

//if ( mValidNoDataValue && (fabs(my - mNoDataValue) < myPrecision || my == mNoDataValue || my != my))
if ( mValidNoDataValue && ( my == mNoDataValue || my != my ) )
if ( mValidNoDataValue && ( fabs( myValue - mNoDataValue ) <= TINY_VALUE || myValue != myValue ) )
{
continue; // NULL
}

myRasterBandStats.sumOfSquares += static_cast < double >
( pow( my - myRasterBandStats.mean, 2 ) );
( pow( myValue - myRasterBandStats.mean, 2 ) );
}
}
} //end of column wise loop
Expand Down Expand Up @@ -1893,7 +1891,7 @@ bool QgsRasterLayer::identify( const QgsPoint& thePoint, QMap<QString, QString>&
#endif
QString v;

if ( mValidNoDataValue && ( mNoDataValue == value || value != value ) )
if ( mValidNoDataValue && ( fabs( value - mNoDataValue ) <= TINY_VALUE || value != value ) )
{
v = tr( "null (no data)" );
}
Expand Down Expand Up @@ -3071,7 +3069,8 @@ bool QgsRasterLayer::readColorTable( int theBandNumber, QList<QgsColorRampShader

void QgsRasterLayer::resetNoDataValue()
{
mNoDataValue = -9999;
mNoDataValue = std::numeric_limits<int>::max();
mValidNoDataValue = false;
if ( mGdalDataset != NULL && GDALGetRasterCount( mGdalDataset ) > 0 )
{
int myRequestValid;
Expand All @@ -3084,7 +3083,7 @@ void QgsRasterLayer::resetNoDataValue()
}
else
{
setNoDataValue( myValue );
setNoDataValue( -9999.0 );
mValidNoDataValue = false;
}
}
Expand Down Expand Up @@ -4340,7 +4339,7 @@ void QgsRasterLayer::drawMultiBandColor( QPainter * theQPainter, QgsRasterViewPo
myBlueValue = readValue( myGdalBlueData, ( GDALDataType )myBlueType,
myRow * theRasterViewPort->drawableAreaXDim + myColumn );

if ( mValidNoDataValue && (( myRedValue == mNoDataValue || myRedValue != myRedValue ) || ( myGreenValue == mNoDataValue || myGreenValue != myGreenValue ) || ( myBlueValue == mNoDataValue || myBlueValue != myBlueValue ) ) )
if ( mValidNoDataValue && (( fabs( myRedValue - mNoDataValue ) <= TINY_VALUE || myRedValue != myRedValue ) || ( fabs( myGreenValue - mNoDataValue ) <= TINY_VALUE || myGreenValue != myGreenValue ) || ( fabs( myBlueValue - mNoDataValue ) <= TINY_VALUE || myBlueValue != myBlueValue ) ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -4466,7 +4465,7 @@ void QgsRasterLayer::drawPalettedSingleBandColor( QPainter * theQPainter, QgsRas
myPixelValue = readValue( myGdalScanData, ( GDALDataType )myDataType,
myRow * theRasterViewPort->drawableAreaXDim + myColumn );

if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -4554,7 +4553,7 @@ void QgsRasterLayer::drawPalettedSingleBandGray( QPainter * theQPainter, QgsRast
myPixelValue = readValue( myGdalScanData, ( GDALDataType )myDataType,
myRow * theRasterViewPort->drawableAreaXDim + myColumn );

if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -4664,7 +4663,7 @@ void QgsRasterLayer::drawPalettedSingleBandPseudoColor( QPainter * theQPainter,
myPixelValue = readValue( myGdalScanData, ( GDALDataType )myDataType,
myRow * theRasterViewPort->drawableAreaXDim + myColumn );

if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -4773,7 +4772,8 @@ void QgsRasterLayer::drawSingleBandGray( QPainter * theQPainter, QgsRasterViewPo
// against myGrayVal will always fail ( nan==nan always
// returns false, by design), hence the slightly odd comparison
// of myGrayVal against itself.
if ( mValidNoDataValue && ( myGrayValue == mNoDataValue || myGrayValue != myGrayValue ) )

if ( mValidNoDataValue && ( fabs( myGrayValue - mNoDataValue ) <= TINY_VALUE || myGrayValue != myGrayValue ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -4876,7 +4876,7 @@ void QgsRasterLayer::drawSingleBandPseudoColor( QPainter * theQPainter,
myPixelValue = readValue( myGdalScanData, myDataType,
myRow * theRasterViewPort->drawableAreaXDim + myColumn );

if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
{
myLineBuffer[ myColumn ] = myDefaultColor;
continue;
Expand Down Expand Up @@ -5198,7 +5198,7 @@ bool QgsRasterLayer::readFile( QString const &theFilename )
//
// Determin the nodatavalue
//
mNoDataValue = -9999; //Standard default?
mNoDataValue = -9999.0; //Standard default?
mValidNoDataValue = false;
int isValid = false;
double myNoDataValue = GDALGetRasterNoDataValue( GDALGetRasterBand( mGdalDataset, 1 ), &isValid );
Expand Down

0 comments on commit 1c3390e

Please sign in to comment.