105 auto outputBlock = std::make_unique<QgsRasterBlock>();
109 return outputBlock.release();
112 std::shared_ptr< QgsRasterBlock > inputBlock(
mInput->block( mBand,
extent, width, height, feedback ) );
114 if ( !inputBlock || inputBlock->isEmpty() )
117 return outputBlock.release();
120 std::shared_ptr< QgsRasterBlock > alphaBlock;
125 if ( !alphaBlock || alphaBlock->isEmpty() )
128 return outputBlock.release();
133 alphaBlock = inputBlock;
138 return outputBlock.release();
141 if ( width == 0 || height == 0 )
142 return outputBlock.release();
147 float cellXSize =
static_cast<float>(
extent.width() ) / width;
148 float cellYSize =
static_cast<float>(
extent.height() ) / height;
149 float zenithRad =
static_cast<float>( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 );
150 float azimuthRad =
static_cast<float>( -1 * mLightAzimuth * M_PI / 180.0 );
151 float cosZenithRad = std::cos( zenithRad );
152 float sinZenithRad = std::sin( zenithRad );
155 float cos_alt_mul_z = cosZenithRad *
static_cast<float>( mZFactor );
156 float cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z;
157 float sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z;
158 float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0f * cos_az_mul_cos_alt_mul_z;
159 float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0f * sin_az_mul_cos_alt_mul_z;
160 float square_z =
static_cast<float>( mZFactor * mZFactor );
161 float sin_altRadians_mul_254 = 254.0f * sinZenithRad;
164 float sin_altRadians_mul_127 = 127.0f * sinZenithRad;
166 float cos225_az_mul_cos_alt_mul_z_mul_127 = -32.87001872802012f * cos_alt_mul_z;
167 float cos_alt_mul_z_mul_127 = 127.0f * cos_alt_mul_z;
178 && inputBlock->dataTypeSize() <= 4 );
184 if ( source.isEmpty() )
192 std::chrono::time_point<std::chrono::system_clock> startTime( std::chrono::system_clock::now() );
200 std::size_t inputDataTypeSize = inputBlock->dataTypeSize();
201 std::size_t outputDataTypeSize = outputBlock->dataTypeSize();
204 switch ( inputBlock->dataType() )
207 typeName = QStringLiteral(
"unsigned char" );
210 typeName = QStringLiteral(
"unsigned int" );
213 typeName = QStringLiteral(
"short" );
216 typeName = QStringLiteral(
"unsigned int" );
219 typeName = QStringLiteral(
"int" );
222 typeName = QStringLiteral(
"float" );
225 throw QgsException( QStringLiteral(
"Unsupported data type for OpenCL processing." ) );
230 source.replace( QLatin1String(
"__global float *scanLine" ), QStringLiteral(
"__global %1 *scanLine" ).arg( typeName ) );
234 std::size_t scanLineWidth( inputBlock->width() + 2 );
235 std::size_t inputSize( inputDataTypeSize * inputBlock->width() );
238 std::size_t bufferWidth( width + 2 );
239 std::size_t bufferSize( inputDataTypeSize * bufferWidth );
244 auto scanLine = std::make_unique<QgsRasterBlock>( inputBlock->dataType(), scanLineWidth, 1 );
251 std::vector<float> rasterParams;
252 rasterParams.push_back( inputBlock->noDataValue() );
253 rasterParams.push_back( outputBlock->noDataValue() );
254 rasterParams.push_back( mZFactor );
255 rasterParams.push_back( cellXSize );
256 rasterParams.push_back( cellYSize );
257 rasterParams.push_back(
static_cast<float>(
mOpacity ) );
260 rasterParams.push_back( cos_az_mul_cos_alt_mul_z_mul_254 );
261 rasterParams.push_back( sin_az_mul_cos_alt_mul_z_mul_254 );
262 rasterParams.push_back( square_z );
263 rasterParams.push_back( sin_altRadians_mul_254 );
266 rasterParams.push_back( sin_altRadians_mul_127 );
267 rasterParams.push_back( cos225_az_mul_cos_alt_mul_z_mul_127 );
268 rasterParams.push_back( cos_alt_mul_z_mul_127 );
271 rasterParams.push_back(
static_cast<float>( qBlue( defaultNodataColor ) ) );
272 rasterParams.push_back(
static_cast<float>( qGreen( defaultNodataColor ) ) );
273 rasterParams.push_back(
static_cast<float>( qRed( defaultNodataColor ) ) );
274 rasterParams.push_back(
static_cast<float>( qAlpha( defaultNodataColor ) ) / 255.0f );
277 rasterParams.push_back(
static_cast<float>( mMultiDirectional ) );
279 cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(),
true,
false,
nullptr );
280 cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
281 cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
282 cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
283 cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer};
285 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, outputDataTypeSize * width,
nullptr,
nullptr );
287 static std::map<Qgis::DataType, cl::Program> programCache;
288 cl::Program program = programCache[inputBlock->dataType()];
289 if ( ! program.get() )
293 program = programCache[inputBlock->dataType()];
300 auto kernel = cl::KernelFunctor <
306 > ( program,
"processNineCellWindow" );
310 std::vector<int> rowIndex = {0, 1, 2};
312 for (
int i = 0; i < height; i++ )
321 feedback->
setProgress( 100.0 *
static_cast< double >( i ) / height );
327 scanLine->resetNoDataValue();
328 queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) );
330 memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i, 0 ), inputSize );
331 queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) );
333 memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i + 1, 0 ), inputSize );
334 queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) );
340 if ( i == inputBlock->height() - 1 )
342 scanLine->resetNoDataValue();
343 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine->bits( ) );
347 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, inputDataTypeSize * 1 , inputSize, inputBlock->bits( i + 1, 0 ) );
351 kernel( cl::EnqueueArgs(
355 *scanLineBuffer[rowIndex[0]],
356 *scanLineBuffer[rowIndex[1]],
357 *scanLineBuffer[rowIndex[2]],
362 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, outputDataTypeSize * outputBlock->width( ), outputBlock->bits( i, 0 ) );
363 std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
366 catch ( cl::Error &e )
381 double pixelValues[9] {0, 0, 0, 0, 0, 0, 0, 0, 0};
382 bool isNoData[9] {
false,
false,
false,
false,
false,
false,
false,
false,
false};
384 for (
int row = 0; row < height; row++ )
386 for (
int col = 0; col < width; col++ )
394 else if ( row == height - 1 )
402 pixelValues[ 0 ] = inputBlock->valueAndNoData( iUp, 0, isNoData[0] );
403 pixelValues[ 1 ] = pixelValues[0];
404 isNoData[1] = isNoData[0];
405 pixelValues[ 2 ] = pixelValues[0];
406 isNoData[2] = isNoData[0];
408 pixelValues[ 3 ] = inputBlock->valueAndNoData( row, 0, isNoData[3] );
409 pixelValues[ 4 ] = pixelValues[3];
410 isNoData[4] = isNoData[3];
411 pixelValues[ 5 ] = pixelValues[3];
412 isNoData[5] = isNoData[3];
414 pixelValues[ 6 ] = inputBlock->valueAndNoData( iDown, 0, isNoData[6] );
415 pixelValues[ 7 ] = pixelValues[6];
416 isNoData[7] = isNoData[6];
417 pixelValues[ 8 ] = pixelValues[6];
418 isNoData[8] = isNoData[6];
423 pixelValues[ 0 ] = pixelValues[1];
424 pixelValues[ 1 ] = pixelValues[2];
425 pixelValues[ 3 ] = pixelValues[4];
426 pixelValues[ 4 ] = pixelValues[5];
427 pixelValues[ 6 ] = pixelValues[7];
428 pixelValues[ 7 ] = pixelValues[8];
429 isNoData[ 0 ] = isNoData[1];
430 isNoData[ 1 ] = isNoData[2];
431 isNoData[ 3 ] = isNoData[4];
432 isNoData[ 4 ] = isNoData[5];
433 isNoData[ 6 ] = isNoData[7];
434 isNoData[ 7 ] = isNoData[8];
438 if ( col < width - 1 )
440 pixelValues[2] = inputBlock->valueAndNoData( iUp, col + 1, isNoData[2] );
441 pixelValues[5] = inputBlock->valueAndNoData( row, col + 1, isNoData[5] );
442 pixelValues[8] = inputBlock->valueAndNoData( iDown, col + 1, isNoData[8] );
447 outputBlock->setColor( row, col, defaultNodataColor );
452 const double x22 = pixelValues[4];
454 const double x11 = isNoData[0] ? x22 : pixelValues[0];
455 const double x21 = isNoData[3] ? x22 : pixelValues[3];
456 const double x31 = isNoData[6] ? x22 : pixelValues[6];
457 const double x12 = isNoData[1] ? x22 : pixelValues[1];
459 const double x32 = isNoData[7] ? x22 : pixelValues[7];
460 const double x13 = isNoData[2] ? x22 : pixelValues[2];
461 const double x23 = isNoData[5] ? x22 : pixelValues[5];
462 const double x33 = isNoData[8] ? x22 : pixelValues[8];
465 const double derX = ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellXSize );
466 const double derY = ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellYSize );
471 if ( !mMultiDirectional )
474 grayValue = std::clamp( ( sin_altRadians_mul_254 -
475 ( derY * cos_az_mul_cos_alt_mul_z_mul_254 -
476 derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) /
477 std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ),
484 const float xx = derX * derX;
485 const float yy = derY * derY;
486 const float xx_plus_yy = xx + yy;
488 if ( xx_plus_yy == 0.0 )
490 grayValue = std::clamp(
static_cast<float>( 1.0 + sin_altRadians_mul_254 ), 0.0f, 255.0f );
495 float val225_mul_127 = sin_altRadians_mul_127 +
496 ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
497 val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
498 float val270_mul_127 = sin_altRadians_mul_127 -
499 derX * cos_alt_mul_z_mul_127;
500 val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
501 float val315_mul_127 = sin_altRadians_mul_127 +
502 ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
503 val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
504 float val360_mul_127 = sin_altRadians_mul_127 -
505 derY * cos_alt_mul_z_mul_127;
506 val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
509 const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
510 const float weight_270 = xx;
511 const float weight_315 = xx_plus_yy - weight_225;
512 const float weight_360 = yy;
513 const float cang_mul_127 = (
514 ( weight_225 * val225_mul_127 +
515 weight_270 * val270_mul_127 +
516 weight_315 * val315_mul_127 +
517 weight_360 * val360_mul_127 ) / xx_plus_yy ) /
518 ( 1 + square_z * xx_plus_yy );
520 grayValue = std::clamp( 1.0f + cang_mul_127, 0.0f, 255.0f );
531 currentAlpha *= alphaBlock->value(
static_cast<qgssize>( row ) * width + col ) / 255.0;
536 outputBlock->setColor( row, col, qRgba( grayValue, grayValue, grayValue, 255 ) );
540 outputBlock->setColor( row, col, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
549 if (
QgsSettings().value( QStringLiteral(
"Map/logCanvasRefreshEvent" ),
false ).toBool() )
552 .arg( useOpenCL ? QStringLiteral(
"OpenCL" ) : QStringLiteral(
"CPU" ) )
555 .arg( std::chrono::duration_cast<std::chrono::milliseconds>( std::chrono::system_clock::now() - startTime ).count() ),
561 return outputBlock.release();
613 QDomNodeList elements = element.elementsByTagName( QStringLiteral(
"sld:RasterSymbolizer" ) );
614 if ( elements.size() == 0 )
618 QDomElement rasterSymbolizerElem = elements.at( 0 ).toElement();
625 QDomElement channelSelectionElem = doc.createElement( QStringLiteral(
"sld:ChannelSelection" ) );
626 elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral(
"sld:Opacity" ) );
627 if ( elements.size() != 0 )
629 rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
633 elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral(
"sld:Geometry" ) );
634 if ( elements.size() != 0 )
636 rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
640 rasterSymbolizerElem.insertBefore( channelSelectionElem, rasterSymbolizerElem.firstChild() );
645 QDomElement channelElem = doc.createElement( QStringLiteral(
"sld:GrayChannel" ) );
646 channelSelectionElem.appendChild( channelElem );
649 QDomElement sourceChannelNameElem = doc.createElement( QStringLiteral(
"sld:SourceChannelName" ) );
650 sourceChannelNameElem.appendChild( doc.createTextNode( QString::number( mBand ) ) );
651 channelElem.appendChild( sourceChannelNameElem );
655 QDomElement shadedReliefElem = doc.createElement( QStringLiteral(
"sld:ShadedRelief" ) );
656 rasterSymbolizerElem.appendChild( shadedReliefElem );
659 QDomElement brightnessOnlyElem = doc.createElement( QStringLiteral(
"sld:BrightnessOnly" ) );
660 brightnessOnlyElem.appendChild( doc.createTextNode( QStringLiteral(
"true" ) ) );
661 shadedReliefElem.appendChild( brightnessOnlyElem );
664 QDomElement reliefFactorElem = doc.createElement( QStringLiteral(
"sld:ReliefFactor" ) );
665 reliefFactorElem.appendChild( doc.createTextNode( QString::number(
zFactor() ) ) );
666 shadedReliefElem.appendChild( reliefFactorElem );
669 QDomElement altitudeVendorOptionElem = doc.createElement( QStringLiteral(
"sld:VendorOption" ) );
670 altitudeVendorOptionElem.setAttribute( QStringLiteral(
"name" ), QStringLiteral(
"altitude" ) );
671 altitudeVendorOptionElem.appendChild( doc.createTextNode( QString::number(
altitude() ) ) );
672 shadedReliefElem.appendChild( altitudeVendorOptionElem );
675 QDomElement azimutVendorOptionElem = doc.createElement( QStringLiteral(
"sld:VendorOption" ) );
676 azimutVendorOptionElem.setAttribute( QStringLiteral(
"name" ), QStringLiteral(
"azimuth" ) );
677 azimutVendorOptionElem.appendChild( doc.createTextNode( QString::number(
azimuth() ) ) );
678 shadedReliefElem.appendChild( azimutVendorOptionElem );
681 QDomElement multidirectionalVendorOptionElem = doc.createElement( QStringLiteral(
"sld:VendorOption" ) );
682 multidirectionalVendorOptionElem.setAttribute( QStringLiteral(
"name" ), QStringLiteral(
"multidirectional" ) );
683 multidirectionalVendorOptionElem.appendChild( doc.createTextNode( QString::number(
multiDirectional() ) ) );
684 shadedReliefElem.appendChild( multidirectionalVendorOptionElem );