From 4263637ba637857a9488dc583480591062daf453 Mon Sep 17 00:00:00 2001 From: Alessandro Pasotti Date: Tue, 29 Oct 2019 17:54:29 +0100 Subject: [PATCH] Fix raster calc OpenCL < operator Fixes #32477 also: - catch build exceptions - remove parenthesis after casting numbers - add cast to raster ref (for fabs overload) with tests --- src/analysis/raster/qgsrastercalcnode.cpp | 9 +- src/analysis/raster/qgsrastercalculator.cpp | 314 +++++++++--------- .../src/analysis/testqgsrastercalculator.cpp | 21 +- 3 files changed, 180 insertions(+), 164 deletions(-) diff --git a/src/analysis/raster/qgsrastercalcnode.cpp b/src/analysis/raster/qgsrastercalcnode.cpp index bd06d3ed1f0..678a26a7249 100644 --- a/src/analysis/raster/qgsrastercalcnode.cpp +++ b/src/analysis/raster/qgsrastercalcnode.cpp @@ -265,7 +265,7 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const break; case opLT: if ( cStyle ) - result = QStringLiteral( "( float ) ( %1 < %2" ).arg( left ).arg( right ); + result = QStringLiteral( "( float ) ( %1 < %2 )" ).arg( left ).arg( right ); else result = QStringLiteral( "%1 < %2" ).arg( left ).arg( right ); break; @@ -344,13 +344,16 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const } break; case tRasterRef: - result = QStringLiteral( "\"%1\"" ).arg( mRasterName ); + if ( cStyle ) + result = QStringLiteral( "( float ) \"%1\"" ).arg( mRasterName ); + else + result = QStringLiteral( "\"%1\"" ).arg( mRasterName ); break; case tNumber: result = QString::number( mNumber ); if ( cStyle ) { - result = QStringLiteral( "( float ) ( %1 )" ).arg( result ); + result = QStringLiteral( "( float ) %1" ).arg( result ); } break; case tMatrix: diff --git a/src/analysis/raster/qgsrastercalculator.cpp b/src/analysis/raster/qgsrastercalculator.cpp index fc8a778cf1c..8bb05ae0a59 100644 --- a/src/analysis/raster/qgsrastercalculator.cpp +++ b/src/analysis/raster/qgsrastercalculator.cpp @@ -427,172 +427,182 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni inputRefs.push_back( entry ); } - // Prepare context and queue - cl::Context ctx( QgsOpenClUtils::context() ); - cl::CommandQueue queue( QgsOpenClUtils::commandQueue() ); - - // Create the C expression - std::vector inputBuffers; - inputBuffers.reserve( inputRefs.size() ); - QStringList inputArgs; - for ( const auto &ref : inputRefs ) + // May throw an openCL exception + try { - cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) ); - inputArgs.append( QStringLiteral( "__global %1 *%2" ) - .arg( ref.typeName ) - .arg( ref.varName ) ); - inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) ); - } + // Prepare context and queue + cl::Context ctx( QgsOpenClUtils::context() ); + cl::CommandQueue queue( QgsOpenClUtils::commandQueue() ); - //qDebug() << cExpression; - - // Create the program - QString programTemplate( R"CL( - // Inputs: - ##INPUT_DESC## - // Expression: ##EXPRESSION_ORIGINAL## - __kernel void rasterCalculator( ##INPUT## - __global float *resultLine - ) - { - // Get the index of the current element - const int i = get_global_id(0); - // Expression - resultLine[i] = ##EXPRESSION##; - } - )CL" ); - - QStringList inputDesc; - for ( const auto &ref : inputRefs ) - { - inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) ); - } - programTemplate = programTemplate.replace( QStringLiteral( "##INPUT_DESC##" ), inputDesc.join( '\n' ) ); - programTemplate = programTemplate.replace( QStringLiteral( "##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join( ',' ).append( ',' ) ) : QChar( ' ' ) ); - programTemplate = programTemplate.replace( QStringLiteral( "##EXPRESSION##" ), cExpression ); - programTemplate = programTemplate.replace( QStringLiteral( "##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) ); - - // qDebug() << programTemplate; - - // Create a program from the kernel source - cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) ); - - // Create the buffers, output is float32 (4 bytes) - // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754 - Q_ASSERT( sizeof( float ) == 4 ); - std::size_t resultBufferSize( 4 * static_cast( mNumOutputColumns ) ); - cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, - resultBufferSize, nullptr, nullptr ); - - auto kernel = cl::Kernel( program, "rasterCalculator" ); - - for ( unsigned int i = 0; i < inputBuffers.size() ; i++ ) - { - kernel.setArg( i, inputBuffers.at( i ) ); - } - kernel.setArg( static_cast( inputBuffers.size() ), resultLineBuffer ); - - QgsOpenClUtils::CPLAllocator resultLine( static_cast( mNumOutputColumns ) ); - - //open output dataset for writing - GDALDriverH outputDriver = openOutputDriver(); - if ( !outputDriver ) - { - mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat ); - return CreateOutputError; - } - - gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) ); - if ( !outputDataset ) - { - mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile ); - return CreateOutputError; - } - - GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() ); - - GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 ); - if ( !outputRasterBand ) - return BandError; - - // Input block (buffer) - std::unique_ptr block; - - // Run kernel on all scanlines - auto rowHeight = mOutputRectangle.height() / mNumOutputRows; - for ( int line = 0; line < mNumOutputRows; line++ ) - { - if ( feedback && feedback->isCanceled() ) - { - break; - } - - if ( feedback ) - { - feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows ); - } - - // Read lines from rasters into the buffers + // Create the C expression + std::vector inputBuffers; + inputBuffers.reserve( inputRefs.size() ); + QStringList inputArgs; for ( const auto &ref : inputRefs ) { - // Read one row - QgsRectangle rect( mOutputRectangle ); - rect.setYMaximum( rect.yMaximum() - rowHeight * line ); - rect.setYMinimum( rect.yMaximum() - rowHeight ); - - // TODO: check if this is too slow - // if crs transform needed - if ( ref.layer->crs() != mOutputCrs ) - { - QgsRasterProjector proj; - proj.setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() ); - proj.setInput( ref.layer->dataProvider() ); - proj.setPrecision( QgsRasterProjector::Exact ); - block.reset( proj.block( ref.band, rect, mNumOutputColumns, 1 ) ); - } - else - { - block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) ); - } - - //for ( int i = 0; i < mNumOutputColumns; i++ ) - // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i ); - //qDebug() << "Writing buffer " << ref.index; - - Q_ASSERT( ref.bufferSize == static_cast( block->data().size( ) ) ); - queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0, - ref.bufferSize, block->bits() ); - + cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) ); + inputArgs.append( QStringLiteral( "__global %1 *%2" ) + .arg( ref.typeName ) + .arg( ref.varName ) ); + inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) ); } - // Run the kernel - queue.enqueueNDRangeKernel( - kernel, - 0, - cl::NDRange( mNumOutputColumns ) - ); - // Write the result - queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, - resultBufferSize, resultLine.get() ); + //qDebug() << cExpression; - //for ( int i = 0; i < mNumOutputColumns; i++ ) - // qDebug() << "Output: " << line << i << " = " << resultLine[i]; - - if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None ) + // Create the program + QString programTemplate( R"CL( + // Inputs: + ##INPUT_DESC## + // Expression: ##EXPRESSION_ORIGINAL## + __kernel void rasterCalculator( ##INPUT## + __global float *resultLine + ) { + // Get the index of the current element + const int i = get_global_id(0); + // Expression + resultLine[i] = ##EXPRESSION##; + } + )CL" ); + + QStringList inputDesc; + for ( const auto &ref : inputRefs ) + { + inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) ); + } + programTemplate = programTemplate.replace( QStringLiteral( "##INPUT_DESC##" ), inputDesc.join( '\n' ) ); + programTemplate = programTemplate.replace( QStringLiteral( "##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join( ',' ).append( ',' ) ) : QChar( ' ' ) ); + programTemplate = programTemplate.replace( QStringLiteral( "##EXPRESSION##" ), cExpression ); + programTemplate = programTemplate.replace( QStringLiteral( "##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) ); + + // qDebug() << programTemplate; + + // Create a program from the kernel source + cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) ); + + // Create the buffers, output is float32 (4 bytes) + // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754 + Q_ASSERT( sizeof( float ) == 4 ); + std::size_t resultBufferSize( 4 * static_cast( mNumOutputColumns ) ); + cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, + resultBufferSize, nullptr, nullptr ); + + auto kernel = cl::Kernel( program, "rasterCalculator" ); + + for ( unsigned int i = 0; i < inputBuffers.size() ; i++ ) + { + kernel.setArg( i, inputBuffers.at( i ) ); + } + kernel.setArg( static_cast( inputBuffers.size() ), resultLineBuffer ); + + QgsOpenClUtils::CPLAllocator resultLine( static_cast( mNumOutputColumns ) ); + + //open output dataset for writing + GDALDriverH outputDriver = openOutputDriver(); + if ( !outputDriver ) + { + mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat ); return CreateOutputError; } - } - if ( feedback && feedback->isCanceled() ) + gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) ); + if ( !outputDataset ) + { + mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile ); + return CreateOutputError; + } + + GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() ); + + GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 ); + if ( !outputRasterBand ) + return BandError; + + // Input block (buffer) + std::unique_ptr block; + + // Run kernel on all scanlines + auto rowHeight = mOutputRectangle.height() / mNumOutputRows; + for ( int line = 0; line < mNumOutputRows; line++ ) + { + if ( feedback && feedback->isCanceled() ) + { + break; + } + + if ( feedback ) + { + feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows ); + } + + // Read lines from rasters into the buffers + for ( const auto &ref : inputRefs ) + { + // Read one row + QgsRectangle rect( mOutputRectangle ); + rect.setYMaximum( rect.yMaximum() - rowHeight * line ); + rect.setYMinimum( rect.yMaximum() - rowHeight ); + + // TODO: check if this is too slow + // if crs transform needed + if ( ref.layer->crs() != mOutputCrs ) + { + QgsRasterProjector proj; + proj.setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() ); + proj.setInput( ref.layer->dataProvider() ); + proj.setPrecision( QgsRasterProjector::Exact ); + block.reset( proj.block( ref.band, rect, mNumOutputColumns, 1 ) ); + } + else + { + block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) ); + } + + //for ( int i = 0; i < mNumOutputColumns; i++ ) + // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i ); + //qDebug() << "Writing buffer " << ref.index; + + Q_ASSERT( ref.bufferSize == static_cast( block->data().size( ) ) ); + queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0, + ref.bufferSize, block->bits() ); + + } + // Run the kernel + queue.enqueueNDRangeKernel( + kernel, + 0, + cl::NDRange( mNumOutputColumns ) + ); + + // Write the result + queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, + resultBufferSize, resultLine.get() ); + + //for ( int i = 0; i < mNumOutputColumns; i++ ) + // qDebug() << "Output: " << line << i << " = " << resultLine[i]; + + if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + return CreateOutputError; + } + } + + if ( feedback && feedback->isCanceled() ) + { + //delete the dataset without closing (because it is faster) + gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile ); + return Canceled; + } + + inputBuffers.clear(); + + } + catch ( cl::Error &e ) { - //delete the dataset without closing (because it is faster) - gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile ); - return Canceled; + mLastError = e.what(); + return CreateOutputError; } - inputBuffers.clear(); - return Success; } #endif diff --git a/tests/src/analysis/testqgsrastercalculator.cpp b/tests/src/analysis/testqgsrastercalculator.cpp index 6897a6e0cf2..d2fe6c251d3 100644 --- a/tests/src/analysis/testqgsrastercalculator.cpp +++ b/tests/src/analysis/testqgsrastercalculator.cpp @@ -724,25 +724,28 @@ void TestQgsRasterCalculator::toString() return calcNode->toString( cStyle ); }; QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), false ), QString( "( \"raster@1\" + 2 )" ) ); - QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), true ), QString( "( \"raster@1\" + ( float ) ( 2 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), true ), QString( "( ( float ) \"raster@1\" + ( float ) 2 )" ) ); QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), false ), QString( "( \"raster@1\"^3 + 2 )" ) ); - QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), true ), QString( "( pow( \"raster@1\", ( float ) ( 3 ) ) + ( float ) ( 2 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), true ), QString( "( pow( ( float ) \"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 )" ), true ), QString( "atan( \"raster@1\" ) * cos( ( ( float ) ( 3 ) + ( float ) ( 2 ) ) )" ) ); + QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), true ), QString( "atan( ( float ) \"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 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * ( 1.4 * (\"raster@1\" + 2) )" ), true ), QString( "( float ) 0.5 * ( float ) 1.4 * ( ( float ) \"raster@1\" + ( float ) 2 )" ) ); QCOMPARE( _test( QStringLiteral( "0.5 * ( 1 > 0 )" ), false ), QString( "0.5 * 1 > 0" ) ); - QCOMPARE( _test( QStringLiteral( "0.5 * ( 1 > 0 )" ), true ), QString( "( float ) ( 0.5 ) * ( float ) ( ( float ) ( 1 ) > ( float ) ( 0 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * ( 1 > 0 )" ), true ), QString( "( float ) 0.5 * ( float ) ( ( float ) 1 > ( float ) 0 )" ) ); // Test negative numbers QCOMPARE( _test( QStringLiteral( "0.5 * ( -1 > 0 )" ), false ), QString( "0.5 * -1 > 0" ) ); - QCOMPARE( _test( QStringLiteral( "0.5 * ( -1 > 0 )" ), true ), QString( "( float ) ( 0.5 ) * ( float ) ( -( float ) ( 1 ) > ( float ) ( 0 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * ( -1 > 0 )" ), true ), QString( "( float ) 0.5 * ( float ) ( -( float ) 1 > ( float ) 0 )" ) ); // Test new functions QCOMPARE( _test( QStringLiteral( "0.5 * abs( -1 )" ), false ), QString( "0.5 * abs( -1 )" ) ); - QCOMPARE( _test( QStringLiteral( "0.5 * abs( -1 )" ), true ), QString( "( float ) ( 0.5 ) * fabs( -( float ) ( 1 ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * abs( -1 )" ), true ), QString( "( float ) 0.5 * fabs( -( float ) 1 )" ) ); QCOMPARE( _test( QStringLiteral( "0.5 * min( -1, 1 )" ), false ), QString( "0.5 * min( -1, 1 )" ) ); - QCOMPARE( _test( QStringLiteral( "0.5 * min( -1, 1 )" ), true ), QString( "( float ) ( 0.5 ) * min( ( float ) ( -( float ) ( 1 ) ), ( float ) ( ( float ) ( 1 ) ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * min( -1, 1 )" ), true ), QString( "( float ) 0.5 * min( ( float ) ( -( float ) 1 ), ( float ) ( ( float ) 1 ) )" ) ); QCOMPARE( _test( QStringLiteral( "0.5 * max( -1, 1 )" ), false ), QString( "0.5 * max( -1, 1 )" ) ); - QCOMPARE( _test( QStringLiteral( "0.5 * max( -1, 1 )" ), true ), QString( "( float ) ( 0.5 ) * max( ( float ) ( -( float ) ( 1 ) ), ( float ) ( ( float ) ( 1 ) ) )" ) ); + QCOMPARE( _test( QStringLiteral( "0.5 * max( -1, 1 )" ), true ), QString( "( float ) 0.5 * max( ( float ) ( -( float ) 1 ), ( float ) ( ( float ) 1 ) )" ) ); + // Test regression #32477 + QCOMPARE( _test( QStringLiteral( R"raw(("r@1"<100.09)*0.1)raw" ), true ), + QString( R"raw(( float ) ( ( float ) "r@1" < ( float ) 100.09 ) * ( float ) 0.1)raw" ) ); } void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()