sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-03-03 00:02:25 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #9305 from elpaso/bugfix-21405-raster-calc-wrong-results

Browse Source 
[opencl] Fix raster calculator operator precedence
This commit is contained in:
Alessandro Pasotti 2019-03-01 14:30:28 +01:00 committed by GitHub
commit 83f286a9e5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 106 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/analysis/raster/qgsrastercalcnode.cpp
View File

@ -210,22 +210,17 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const
if ( mRight ) if ( mRight )
right = mRight->toString( cStyle ); right = mRight->toString( cStyle );
auto floatCast = [ ]( const QString s ) -> QString
{
return QStringLiteral( "(float) ( %1 )" ).arg( s );
};
switch ( mType ) switch ( mType )
{ {
case tOperator: case tOperator:
switch ( mOperator ) switch ( mOperator )
{ {
case opPLUS: case opPLUS:
result = QStringLiteral( "%1 + %2" ).arg( left ).arg( right ); result = QStringLiteral( "( %1 + %2 )" ).arg( left ).arg( right );
break; break;
case opMINUS: case opMINUS:
case opSIGN: case opSIGN:
result = QStringLiteral( "%1 - %2" ).arg( left ).arg( right ); result = QStringLiteral( "( %1 - %2 )" ).arg( left ).arg( right );
break; break;
case opMUL: case opMUL:
result = QStringLiteral( "%1 * %2" ).arg( left ).arg( right ); result = QStringLiteral( "%1 * %2" ).arg( left ).arg( right );
@ -235,7 +230,7 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const
break; break;
case opPOW: case opPOW:
if ( cStyle ) if ( cStyle )
result = QStringLiteral( "pow( %1, %2 )" ).arg( floatCast( left ) ).arg( floatCast( right ) ); result = QStringLiteral( "pow( %1, %2 )" ).arg( left ).arg( right );
else else
result = QStringLiteral( "%1^%2" ).arg( left ).arg( right ); result = QStringLiteral( "%1^%2" ).arg( left ).arg( right );
break; break;
@ -273,57 +268,30 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const
result = QStringLiteral( "%1 OR %2" ).arg( left ).arg( right ); result = QStringLiteral( "%1 OR %2" ).arg( left ).arg( right );
break; break;
case opSQRT: case opSQRT:
if ( cStyle )
result = QStringLiteral( "sqrt( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "sqrt( %1 )" ).arg( left ); result = QStringLiteral( "sqrt( %1 )" ).arg( left );
break; break;
case opSIN: case opSIN:
if ( cStyle )
result = QStringLiteral( "sin( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "sin( %1 )" ).arg( left ); result = QStringLiteral( "sin( %1 )" ).arg( left );
break; break;
case opCOS: case opCOS:
if ( cStyle )
result = QStringLiteral( "cos( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "cos( %1 )" ).arg( left ); result = QStringLiteral( "cos( %1 )" ).arg( left );
break; break;
case opTAN: case opTAN:
if ( cStyle )
result = QStringLiteral( "tan( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "tan( %1 )" ).arg( left ); result = QStringLiteral( "tan( %1 )" ).arg( left );
break; break;
case opASIN: case opASIN:
if ( cStyle )
result = QStringLiteral( "asin( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "asin( %1 )" ).arg( left ); result = QStringLiteral( "asin( %1 )" ).arg( left );
break; break;
case opACOS: case opACOS:
if ( cStyle )
result = QStringLiteral( "acos( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "acos( %1 )" ).arg( left ); result = QStringLiteral( "acos( %1 )" ).arg( left );
break; break;
case opATAN: case opATAN:
if ( cStyle )
result = QStringLiteral( "atan( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "atan( %1 )" ).arg( left ); result = QStringLiteral( "atan( %1 )" ).arg( left );
break; break;
case opLOG: case opLOG:
if ( cStyle )
result = QStringLiteral( "log( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "log( %1 )" ).arg( left ); result = QStringLiteral( "log( %1 )" ).arg( left );
break; break;
case opLOG10: case opLOG10:
if ( cStyle )
result = QStringLiteral( "log10( %1 )" ).arg( floatCast( left ) );
else
result = QStringLiteral( "log10( %1 )" ).arg( left ); result = QStringLiteral( "log10( %1 )" ).arg( left );
break; break;
case opNONE: case opNONE:
@ -335,6 +303,10 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const
break; break;
case tNumber: case tNumber:
result = QString::number( mNumber ); result = QString::number( mNumber );
if ( cStyle )
{
result = QStringLiteral( "(float) ( %1 )" ).arg( result );
}
break; break;
case tMatrix: case tMatrix:
break; break;

98
tests/src/analysis/testqgsrastercalculator.cpp
View File

@ -14,6 +14,10 @@ Email : nyall dot dawson at gmail dot com
***************************************************************************/ ***************************************************************************/
#include "qgstest.h" #include "qgstest.h"
#ifdef HAVE_OPENCL
#include "qgsopenclutils.h"
#endif
#include "qgsrastercalculator.h" #include "qgsrastercalculator.h"
#include "qgsrastercalcnode.h" #include "qgsrastercalcnode.h"
#include "qgsrasterdataprovider.h" #include "qgsrasterdataprovider.h"
@ -59,6 +63,7 @@ class TestQgsRasterCalculator : public QObject
void findNodes(); void findNodes();
void testRasterEntries(); void testRasterEntries();
void calcFormulasWithReprojectedLayers();
private: private:
@ -66,15 +71,21 @@ class TestQgsRasterCalculator : public QObject
QgsRasterLayer *mpLandsatRasterLayer4326 = nullptr; QgsRasterLayer *mpLandsatRasterLayer4326 = nullptr;
}; };
void TestQgsRasterCalculator::initTestCase() void TestQgsRasterCalculator::initTestCase()
{ {
// //
// Runs once before any tests are run // Runs once before any tests are run
// //
// init QGIS's paths - true means that all path will be inited from prefix // Set up the QgsSettings environment
QCoreApplication::setOrganizationName( QStringLiteral( "QGIS" ) );
QCoreApplication::setOrganizationDomain( QStringLiteral( "qgis.org" ) );
QCoreApplication::setApplicationName( QStringLiteral( "QGIS-TEST" ) );
QgsApplication::init(); QgsApplication::init();
QgsApplication::initQgis(); QgsApplication::initQgis();
QString testDataDir = QStringLiteral( TEST_DATA_DIR ) + '/'; //defined in CmakeLists.txt QString testDataDir = QStringLiteral( TEST_DATA_DIR ) + '/'; //defined in CmakeLists.txt
QString landsatFileName = testDataDir + "landsat.tif"; QString landsatFileName = testDataDir + "landsat.tif";
@ -99,8 +110,13 @@ void TestQgsRasterCalculator::cleanupTestCase()
void TestQgsRasterCalculator::init() void TestQgsRasterCalculator::init()
{ {
#ifdef HAVE_OPENCL
QgsOpenClUtils::setEnabled( false );
// Reset to default in case some tests mess it up
QgsOpenClUtils::setSourcePath( QDir( QgsApplication::pkgDataPath() ).absoluteFilePath( QStringLiteral( "resources/opencl_programs" ) ) );
#endif
} }
void TestQgsRasterCalculator::cleanup() void TestQgsRasterCalculator::cleanup()
{ {
@ -676,12 +692,78 @@ void TestQgsRasterCalculator::toString()
return error; return error;
return calcNode->toString( cStyle ); return calcNode->toString( cStyle );
}; };
QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), false ), QString( "\"raster@1\" + 2" ) ); QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), false ), QString( "( \"raster@1\" + 2 )" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), true ), QString( "\"raster@1\" + 2" ) ); QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), true ), QString( "( \"raster@1\" + (float) ( 2 ) )" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), false ), QString( "\"raster@1\"^3 + 2" ) ); QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), false ), QString( "( \"raster@1\"^3 + 2 )" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), true ), QString( "pow( (float) ( \"raster@1\" ), (float) ( 3 ) ) + 2" ) ); QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), true ), QString( "( pow( \"raster@1\", (float) ( 3 ) ) + (float) ( 2 ) )" ) );
QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), false ), QString( "atan( \"raster@1\" ) * cos( 3 + 2 )" ) ); QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), false ), QString( "atan( \"raster@1\" ) * cos( ( 3 + 2 ) )" ) );
QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), true ), QString( "atan( (float) ( \"raster@1\" ) ) * cos( (float) ( 3 + 2 ) )" ) ); QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), true ), QString( "atan( \"raster@1\" ) * cos( ( (float) ( 3 ) + (float) ( 2 ) ) )" ) );
QCOMPARE( _test( QStringLiteral( "0.5 * ( 1.4 * (\"raster@1\" + 2) )" ), false ), QString( "0.5 * 1.4 * ( \"raster@1\" + 2 )" ) );
QCOMPARE( _test( QStringLiteral( "0.5 * ( 1.4 * (\"raster@1\" + 2) )" ), true ), QString( "(float) ( 0.5 ) * (float) ( 1.4 ) * ( \"raster@1\" + (float) ( 2 ) )" ) );
}
void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()
{
QgsRasterCalculatorEntry entry1;
entry1.bandNumber = 1;
entry1.raster = mpLandsatRasterLayer;
entry1.ref = QStringLiteral( "landsat@1" );
QgsRasterCalculatorEntry entry2;
entry2.bandNumber = 2;
entry2.raster = mpLandsatRasterLayer4326;
entry2.ref = QStringLiteral( "landsat_4326@2" );
QVector<QgsRasterCalculatorEntry> entries;
entries << entry1 << entry2;
QgsCoordinateReferenceSystem crs;
crs.createFromId( 32633, QgsCoordinateReferenceSystem::EpsgCrsId );
QgsRectangle extent( 783235, 3348110, 783350, 3347960 );
auto _chk = [ = ]( const QString & formula, const std::vector<float> &values, bool useOpenCL )
{
#ifdef HAVE_OPENCL
if ( ! QgsOpenClUtils::available() )
return ;
QgsOpenClUtils::setEnabled( useOpenCL );
#endif
QTemporaryFile tmpFile;
tmpFile.open(); // fileName is not available until open
QString tmpName = tmpFile.fileName();
tmpFile.close();
QgsRasterCalculator rc( formula,
tmpName,
QStringLiteral( "GTiff" ),
extent, crs, 2, 3, entries );
QCOMPARE( static_cast< int >( rc.processCalculation() ), 0 );
//open output file and check results
QgsRasterLayer *result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) );
QCOMPARE( result->width(), 2 );
QCOMPARE( result->height(), 3 );
QgsRasterBlock *block = result->dataProvider()->block( 1, extent, 2, 3 );
qDebug() << block->value( 0, 0 ) << block->value( 0, 1 ) << block->value( 1, 0 ) << block->value( 1, 1 ) << block->value( 2, 0 ) << block->value( 2, 1 );
const float epsilon { 0.0001f };
QVERIFY( std::abs( ( static_cast<float>( block->value( 0, 0 ) ) - values[0] ) / values[0] ) < epsilon );
QVERIFY( std::abs( ( static_cast<float>( block->value( 0, 1 ) ) - values[1] ) / values[1] ) < epsilon );
QVERIFY( std::abs( ( static_cast<float>( block->value( 1, 0 ) ) - values[2] ) / values[2] ) < epsilon );
QVERIFY( std::abs( ( static_cast<float>( block->value( 1, 1 ) ) - values[3] ) / values[3] ) < epsilon );
QVERIFY( std::abs( ( static_cast<float>( block->value( 2, 0 ) ) - values[4] ) / values[4] ) < epsilon );
QVERIFY( std::abs( ( static_cast<float>( block->value( 2, 1 ) ) - values[5] ) / values[5] ) < epsilon );
delete result;
delete block;
};
_chk( QStringLiteral( "\"landsat@1\" + \"landsat_4326@2\"" ), {264.0, 263.0, 264.0, 264.0, 266.0, 261.0}, false );
_chk( QStringLiteral( "\"landsat@1\" + \"landsat_4326@2\"" ), {264.0, 263.0, 264.0, 264.0, 266.0, 261.0}, true );
_chk( QStringLiteral( "\"landsat@1\"^2 + 3 + \"landsat_4326@2\"" ), {15767, 15766, 15519, 15767, 15769, 15516}, false );
_chk( QStringLiteral( "\"landsat@1\"^2 + 3 + \"landsat_4326@2\"" ), {15767, 15766, 15519, 15767, 15769, 15516}, true );
_chk( QStringLiteral( "0.5*((2*\"landsat@1\"+1)-sqrt((2*\"landsat@1\"+1)^2-8*(\"landsat@1\"-\"landsat_4326@2\")))" ), {-0.111504, -0.103543, -0.128448, -0.111504, -0.127425, -0.104374}, false );
_chk( QStringLiteral( "0.5*((2*\"landsat@1\"+1)-sqrt((2*\"landsat@1\"+1)^2-8*(\"landsat@1\"-\"landsat_4326@2\")))" ), {-0.111504, -0.103543, -0.128448, -0.111504, -0.127425, -0.104374}, true );
} }

3
tests/src/core/testqgsopenclutils.cpp
View File

@ -29,9 +29,6 @@
class TestQgsOpenClUtils: public QObject class TestQgsOpenClUtils: public QObject
{ {
Q_OBJECT Q_OBJECT
public:
//void testRunMakeProgram();
private slots: private slots:
void initTestCase();// will be called before the first testfunction is executed. void initTestCase();// will be called before the first testfunction is executed.