Skip to content

Commit 47b8a21

Browse files
author
ersts
committedJul 28, 2009
-Fix no data value from being included in histogram and stats
-Closes ticket #1380 git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@11189 c8812cc2-4d05-0410-92ff-de0c093fc19c

File tree

1 file changed

+33
-33
lines changed

1 file changed

+33
-33
lines changed
 

‎src/core/raster/qgsrasterlayer.cpp

Lines changed: 33 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,13 @@ email : tim at linfiniti.com
7070
# endif
7171
#endif
7272

73+
// Comparison value for equality; i.e., we shouldn't directly compare two
74+
// floats so it's better to take their difference and see if they're within
75+
// a certain range -- in this case twenty times the smallest value that
76+
// doubles can take for the current system. (Yes, 20 was arbitrary.)
77+
#define TINY_VALUE std::numeric_limits<double>::epsilon() * 20
78+
79+
7380
QgsRasterLayer::QgsRasterLayer(
7481
QString const & path,
7582
QString const & baseName,
@@ -106,7 +113,7 @@ QgsRasterLayer::QgsRasterLayer(
106113

107114
mBandCount = 0;
108115
mHasPyramids = false;
109-
mNoDataValue = -9999;
116+
mNoDataValue = -9999.0;
110117
mValidNoDataValue = false;
111118

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

762-
void *myData = CPLMalloc( myXBlockSize * myYBlockSize * GDALGetDataTypeSize( myDataType ) / 8 );
769+
void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( GDALGetDataTypeSize( myDataType ) / 8 ) );
763770

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

767-
// Comparison value for equality; i.e., we shouldn't directly compare two
768-
// floats so it's better to take their difference and see if they're within
769-
// a certain range -- in this case twenty times the smallest value that
770-
// doubles can take for the current system. (Yes, 20 was arbitrary.)
771-
double myPrecision = std::numeric_limits<double>::epsilon() * 20;
772-
Q_UNUSED( myPrecision );
773-
774774
//ifdefs below to remove compiler warning about unused vars
775775
#ifdef QGISDEBUG
776776
int success;
@@ -841,10 +841,9 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
841841
{
842842
for ( int iX = 0; iX < nXValid; iX++ )
843843
{
844-
double my = readValue( myData, myDataType, iX + iY * myXBlockSize );
844+
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
845845

846-
//if ( mValidNoDataValue && (fabs(my - mNoDataValue) < myPrecision || my == mNoDataValue || my != my))
847-
if ( mValidNoDataValue && ( my == mNoDataValue || my != my ) )
846+
if ( mValidNoDataValue && ( fabs( myValue - mNoDataValue ) <= TINY_VALUE || myValue != myValue ) )
848847
{
849848
continue; // NULL
850849
}
@@ -854,23 +853,23 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
854853
{
855854
//this is the first iteration so initialise vars
856855
myFirstIterationFlag = false;
857-
myRasterBandStats.minimumValue = my;
858-
myRasterBandStats.maximumValue = my;
856+
myRasterBandStats.minimumValue = myValue;
857+
myRasterBandStats.maximumValue = myValue;
859858
++myRasterBandStats.elementCount;
860859
} //end of true part for first iteration check
861860
else
862861
{
863862
//this is done for all subsequent iterations
864-
if ( my < myRasterBandStats.minimumValue )
863+
if ( myValue < myRasterBandStats.minimumValue )
865864
{
866-
myRasterBandStats.minimumValue = my;
865+
myRasterBandStats.minimumValue = myValue;
867866
}
868-
if ( my > myRasterBandStats.maximumValue )
867+
if ( myValue > myRasterBandStats.maximumValue )
869868
{
870-
myRasterBandStats.maximumValue = my;
869+
myRasterBandStats.maximumValue = myValue;
871870
}
872871

873-
myRasterBandStats.sum += my;
872+
myRasterBandStats.sum += myValue;
874873
++myRasterBandStats.elementCount;
875874
} //end of false part for first iteration check
876875
}
@@ -912,16 +911,15 @@ const QgsRasterBandStats QgsRasterLayer::bandStatistics( int theBandNo )
912911
{
913912
for ( int iX = 0; iX < nXValid; iX++ )
914913
{
915-
double my = readValue( myData, myDataType, iX + iY * myXBlockSize );
914+
double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) );
916915

917-
//if ( mValidNoDataValue && (fabs(my - mNoDataValue) < myPrecision || my == mNoDataValue || my != my))
918-
if ( mValidNoDataValue && ( my == mNoDataValue || my != my ) )
916+
if ( mValidNoDataValue && ( fabs( myValue - mNoDataValue ) <= TINY_VALUE || myValue != myValue ) )
919917
{
920918
continue; // NULL
921919
}
922920

923921
myRasterBandStats.sumOfSquares += static_cast < double >
924-
( pow( my - myRasterBandStats.mean, 2 ) );
922+
( pow( myValue - myRasterBandStats.mean, 2 ) );
925923
}
926924
}
927925
} //end of column wise loop
@@ -1893,7 +1891,7 @@ bool QgsRasterLayer::identify( const QgsPoint& thePoint, QMap<QString, QString>&
18931891
#endif
18941892
QString v;
18951893

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

30723070
void QgsRasterLayer::resetNoDataValue()
30733071
{
3074-
mNoDataValue = -9999;
3072+
mNoDataValue = std::numeric_limits<int>::max();
3073+
mValidNoDataValue = false;
30753074
if ( mGdalDataset != NULL && GDALGetRasterCount( mGdalDataset ) > 0 )
30763075
{
30773076
int myRequestValid;
@@ -3084,7 +3083,7 @@ void QgsRasterLayer::resetNoDataValue()
30843083
}
30853084
else
30863085
{
3087-
setNoDataValue( myValue );
3086+
setNoDataValue( -9999.0 );
30883087
mValidNoDataValue = false;
30893088
}
30903089
}
@@ -4340,7 +4339,7 @@ void QgsRasterLayer::drawMultiBandColor( QPainter * theQPainter, QgsRasterViewPo
43404339
myBlueValue = readValue( myGdalBlueData, ( GDALDataType )myBlueType,
43414340
myRow * theRasterViewPort->drawableAreaXDim + myColumn );
43424341

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

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

4557-
if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
4556+
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
45584557
{
45594558
myLineBuffer[ myColumn ] = myDefaultColor;
45604559
continue;
@@ -4664,7 +4663,7 @@ void QgsRasterLayer::drawPalettedSingleBandPseudoColor( QPainter * theQPainter,
46644663
myPixelValue = readValue( myGdalScanData, ( GDALDataType )myDataType,
46654664
myRow * theRasterViewPort->drawableAreaXDim + myColumn );
46664665

4667-
if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
4666+
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
46684667
{
46694668
myLineBuffer[ myColumn ] = myDefaultColor;
46704669
continue;
@@ -4773,7 +4772,8 @@ void QgsRasterLayer::drawSingleBandGray( QPainter * theQPainter, QgsRasterViewPo
47734772
// against myGrayVal will always fail ( nan==nan always
47744773
// returns false, by design), hence the slightly odd comparison
47754774
// of myGrayVal against itself.
4776-
if ( mValidNoDataValue && ( myGrayValue == mNoDataValue || myGrayValue != myGrayValue ) )
4775+
4776+
if ( mValidNoDataValue && ( fabs( myGrayValue - mNoDataValue ) <= TINY_VALUE || myGrayValue != myGrayValue ) )
47774777
{
47784778
myLineBuffer[ myColumn ] = myDefaultColor;
47794779
continue;
@@ -4876,7 +4876,7 @@ void QgsRasterLayer::drawSingleBandPseudoColor( QPainter * theQPainter,
48764876
myPixelValue = readValue( myGdalScanData, myDataType,
48774877
myRow * theRasterViewPort->drawableAreaXDim + myColumn );
48784878

4879-
if ( mValidNoDataValue && ( myPixelValue == mNoDataValue || myPixelValue != myPixelValue ) )
4879+
if ( mValidNoDataValue && ( fabs( myPixelValue - mNoDataValue ) <= TINY_VALUE || myPixelValue != myPixelValue ) )
48804880
{
48814881
myLineBuffer[ myColumn ] = myDefaultColor;
48824882
continue;
@@ -5198,7 +5198,7 @@ bool QgsRasterLayer::readFile( QString const &theFilename )
51985198
//
51995199
// Determin the nodatavalue
52005200
//
5201-
mNoDataValue = -9999; //Standard default?
5201+
mNoDataValue = -9999.0; //Standard default?
52025202
mValidNoDataValue = false;
52035203
int isValid = false;
52045204
double myNoDataValue = GDALGetRasterNoDataValue( GDALGetRasterBand( mGdalDataset, 1 ), &isValid );

0 commit comments

Comments
 (0)
Please sign in to comment.